added normalization to pair correlation calc

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 80cf572..f2923c3 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -2322,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;
@@ -2339,13 +2375,13 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
        unsigned char bc;
        t_3dvec dist;
        double d;
-       //double norm;
+       double norm;
        int o,s;
        unsigned char ibrand;
 
        lc=&(moldyn->lc);
 
-       slots=moldyn->cutoff/dr;
+       slots=2.0*moldyn->cutoff/dr;
        o=2*slots;
 
        if(slots*dr<=moldyn->cutoff)
@@ -2422,7 +2458,7 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
                                d=v3_absolute_square(&dist);
 
                                /* ignore if greater or equal cutoff */
-                               if(d>=moldyn->cutoff_square)
+                               if(d>=4.0*moldyn->cutoff_square)
                                        continue;
 
                                /* fill the slots */
@@ -2458,16 +2494,17 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
                }
        }
 
-       /* normalization 
+       /* normalization */
        for(i=1;i<slots;i++) {
-                // normalization: 4 pi r r dr
+                // normalization: 4 pi r^2 dr
                 // here: not double counting pairs -> 2 pi r r dr
-               norm=2*M_PI*moldyn->count*(i*dr*i*dr)*dr;
+                // ... and actually it's a constant times r^2
+               norm=i*i*dr*dr;
                stat[i]/=norm;
                stat[slots+i]/=norm;
                stat[o+i]/=norm;
        }
-       */
+       /* */
 
        if(ptr==NULL) {
                /* todo: store/print pair correlation function */