still playing around with heat capacity ...

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index c130ef4..4c43f2e 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -1838,6 +1838,7 @@ int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
        int count,ret;
        double time,pot,kin,tot;
        double p_sum,k_sum,t_sum;
+       double p2_sum,k2_sum,t2_sum;
        char buf[64];
        char file[128+7];
 
@@ -1850,7 +1851,7 @@ int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
                return fd;
        }
 
-       /* first calc the averages */
+       /* calc the averages of A and A^2 */
        p_sum=0.0;
        k_sum=0.0;
        t_sum=0.0;
@@ -1865,38 +1866,21 @@ int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
                p_sum+=pot;
                k_sum+=kin;
                t_sum+=tot;
+               p2_sum+=(pot*pot);
+               k2_sum+=(kin*kin);
+               t2_sum+=(tot*tot);
                count+=1;
        }
 
-       moldyn->p_m=p_sum/count;
+       /* averages */
        moldyn->k_m=k_sum/count;
+       moldyn->p_m=p_sum/count;
        moldyn->t_m=t_sum/count;
 
-       /* mean square fluctuations */
-       if(lseek(fd,SEEK_SET,0)<0) {
-               perror("[moldyn] lseek");
-               return -1;
-       }
-       count=0;
-       p_sum=0.0;
-       k_sum=0.0;
-       t_sum=0.0;
-       while(1) {
-               ret=get_line(fd,buf,63);
-               if(ret<=0) break;
-               if(buf[0]=='#') continue;
-               sscanf(buf,"%lf %lf %lf %lf",&time,&kin,&pot,&tot);
-               if(time<start) continue;
-               if(time>end) break;
-               k_sum+=((kin-moldyn->k_m)*(kin-moldyn->k_m));
-               p_sum+=((pot-moldyn->p_m)*(pot-moldyn->p_m));
-               t_sum+=((tot-moldyn->t_m)*(tot-moldyn->t_m));
-               count+=1;
-       }
-
-       moldyn->dp2_m=p_sum/count;
-       moldyn->dk2_m=k_sum/count;
-       moldyn->dt2_m=t_sum/count;
+       /* rms */
+       moldyn->dk2_m=k2_sum/count-moldyn->k_m*moldyn->k_m;
+       moldyn->dp2_m=p2_sum/count-moldyn->p_m*moldyn->p_m;
+       moldyn->dt2_m=t2_sum/count-moldyn->t_m*moldyn->t_m;
 
        printf("  averages   : %f %f %f\n",moldyn->k_m,
                                        moldyn->p_m,