start of diffusion calc function

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index e2c84a5..93c72da 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -1738,6 +1738,9 @@ int velocity_verlet(t_moldyn *moldyn) {
        tau_square=moldyn->tau_square;
 
        for(i=0;i<count;i++) {
+               /* check whether fixed atom */
+               if(atom[i].attr&ATOM_ATTR_FP)
+                       continue;
                /* new positions */
                h=0.5/atom[i].mass;
                v3_scale(&delta,&(atom[i].v),tau);
@@ -1761,6 +1764,9 @@ int velocity_verlet(t_moldyn *moldyn) {
        potential_force_calc(moldyn);
 
        for(i=0;i<count;i++) {
+               /* check whether fixed atom */
+               if(atom[i].attr&ATOM_ATTR_FP)
+                       continue;
                /* again velocities [actually v(t+tau)] */
                v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
                v3_add(&(atom[i].v),&(atom[i].v),&delta);
@@ -2088,9 +2094,9 @@ int potential_force_calc(t_moldyn *moldyn) {
        //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
        if(moldyn->time>DSTART&&moldyn->time<DEND) {
                printf("force:\n");
-               printf("  x: %0.40f\n",moldyn->atom[5832].f.x);
-               printf("  y: %0.40f\n",moldyn->atom[5832].f.y);
-               printf("  z: %0.40f\n",moldyn->atom[5832].f.z);
+               printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
+               printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
+               printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
        }
 #endif
 
@@ -2316,6 +2322,42 @@ int pair_correlation_init(t_moldyn *moldyn,double dr) {
        return 0;
 }
 
+int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
+
+       int i;
+       t_atom *atom;
+       t_3dvec dist;
+       double d2;
+       int a_cnt;
+       int b_cnt;
+
+       atom=moldyn->atom;
+       dc[0]=0;
+       dc[1]=0;
+       dc[2]=0;
+       a_cnt=0;
+       b_cnt=0;
+
+       for(i=0;i<moldyn->count;i++) {
+
+               v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
+               check_per_bound(moldyn,&dist);
+               d2=v3_absolute_square(&dist);
+
+               if(atom[i].brand) {
+                       b_cnt+=1;
+                       dc[1]+=d2;
+               }
+               else {
+                       a_cnt+=1;
+               }
+
+               dc[2]+=d2;
+       }       
+               
+       return 0;
+}
+
 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 
        int slots;