still playing around with heat capacity ...

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,

diff --git a/sic.c b/sic.c
index d50c35e..07f32bc 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -209,19 +209,25 @@ int main(int argc,char **argv) {
        albe_mult_complete_params(&ap);
 
        /* set (initial) dimensions of simulation volume */
+#ifdef ALBLE
        set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE);
-       //set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
        //set_dim(&md,6*LC_C_ALBE,6*LC_C_ALBE,6*LC_C_ALBE,TRUE);
-       //set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE);
        //set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE);
+#else
+       set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+       //set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE);
+#endif
 
        /* set periodic boundary conditions in all directions */
        set_pbc(&md,TRUE,TRUE,TRUE);
 
        /* create the lattice / place atoms */
-       //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-       //create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C,
+#ifdef ALBLE
        create_lattice(&md,DIAMOND,LC_SI_ALBE,SI,M_SI,
+#else
+       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+       //create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C,
+#endif
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
        //               ATOM_ATTR_2BP|ATOM_ATTR_HB,
                       0,6,6,6,NULL);
@@ -278,6 +284,7 @@ int main(int argc,char **argv) {
        set_pressure(&md,BAR);
 
        /* set p/t scaling */
+       //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
        //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,
        //                 T_SCALE_BERENDSEN,100.0);
        //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
@@ -288,7 +295,7 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        /* initial configuration */
-       moldyn_add_schedule(&md,10000,1.0);
+       moldyn_add_schedule(&md,3000,1.0);
        //moldyn_add_schedule(&md,1000,1.0);
        //moldyn_add_schedule(&md,1000,1.0);
        //moldyn_add_schedule(&md,1000,1.0);
@@ -331,7 +338,7 @@ return 0;
         */
 
        /* response functions expressed by energy fluctuations */
-       calc_fluctuations(1000.0,9999.0,&md);
+       calc_fluctuations(1000.0,2999.0,&md);
        get_heat_capacity(&md);
 
        /* close */