changed energy fluctuation calc and heat cap calc

diff --git a/moldyn.c b/moldyn.c
index 4c43f2e..cd2a803 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -504,6 +504,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
                check_per_bound(moldyn,&(atom[ret].r));
        }
 
+       /* update total system mass */
+       total_mass_calc(moldyn);
+
        return ret;
 }
 
@@ -645,6 +648,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        atom[count].tag=count;
        atom[count].attr=attr;
 
+       /* update total system mass */
+       total_mass_calc(moldyn);
+
        return 0;
 }
 
@@ -705,6 +711,18 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) {
        return 0;
 }
 
+double total_mass_calc(t_moldyn *moldyn) {
+
+       int i;
+
+       moldyn->mass=0.0;
+
+       for(i=0;i<moldyn->count;i++)
+               moldyn->mass+=moldyn->atom[i].mass;
+
+       return moldyn->mass;
+}
+
 double temperature_calc(t_moldyn *moldyn) {
 
        /* assume up to date kinetic energy, which is 3/2 N k_B T */
@@ -825,7 +843,56 @@ double pressure_calc(t_moldyn *moldyn) {
        moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps;
 
        return moldyn->p;
-}      
+}
+
+int energy_fluctuation_calc(t_moldyn *moldyn) {
+
+       /* assume up to date energies */
+
+       /* kinetic energy */
+       moldyn->k_sum+=moldyn->ekin;
+       moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
+       moldyn->k_mean=moldyn->k_sum/moldyn->total_steps;
+       moldyn->k2_mean=moldyn->k2_sum/moldyn->total_steps;
+       moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean);
+
+       /* potential energy */
+       moldyn->v_sum+=moldyn->energy;
+       moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
+       moldyn->v_mean=moldyn->v_sum/moldyn->total_steps;
+       moldyn->v2_mean=moldyn->v2_sum/moldyn->total_steps;
+       moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean);
+
+       return 0;
+}
+
+int get_heat_capacity(t_moldyn *moldyn) {
+
+       double temp2,ighc;
+
+       /* (temperature average)^2 */
+       temp2=moldyn->mean_t*moldyn->mean_t;
+       printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
+              moldyn->mean_t);
+
+       /* ideal gas contribution */
+       ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
+       printf("  ideal gas contribution: %f\n",
+              ighc/moldyn->mass*KILOGRAM/JOULE);
+
+       /* specific heat for nvt ensemble */
+       moldyn->c_v_nvt=moldyn->dv2_mean/(K_BOLTZMANN*temp2)+ighc;
+       moldyn->c_v_nvt/=moldyn->mass;
+
+       /* specific heat for nve ensemble */
+       moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_mean/(ighc*K_BOLTZMANN*temp2)));
+       moldyn->c_v_nve/=moldyn->mass;
+
+       printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
+       printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+
+       return 0;
+}
 
 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
@@ -1307,6 +1374,7 @@ return 0;
                e_kin_calc(moldyn);
                temperature_calc(moldyn);
                pressure_calc(moldyn);
+               energy_fluctuation_calc(moldyn);
                //tp=thermodynamic_pressure_calc(moldyn);
 //printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR);
 
@@ -1379,6 +1447,8 @@ return 0;
                               moldyn->mean_gp/BAR,
                               moldyn->volume);
                        fflush(stdout);
+printf("\n");
+get_heat_capacity(moldyn);
                }
 
                /* increase absolute time */
@@ -1811,7 +1881,7 @@ int moldyn_bc_check(t_moldyn *moldyn) {
 }
 
 /*
- * postprocessing functions
+ * post processing functions
  */
 
 int get_line(int fd,char *line,int max) {
@@ -1832,94 +1902,3 @@ int get_line(int fd,char *line,int max) {
        }
 }
 
-int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
-
-       int fd;
-       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];
-
-       printf("[moldyn] calculating energy fluctuations [eV]:\n");
-
-       snprintf(file,128+7,"%s/energy",moldyn->vlsdir);
-       fd=open(file,O_RDONLY);
-       if(fd<0) {
-               perror("[moldyn] post proc energy open");
-               return fd;
-       }
-
-       /* calc the averages of A and A^2 */
-       p_sum=0.0;
-       k_sum=0.0;
-       t_sum=0.0;
-       count=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;
-               p_sum+=pot;
-               k_sum+=kin;
-               t_sum+=tot;
-               p2_sum+=(pot*pot);
-               k2_sum+=(kin*kin);
-               t2_sum+=(tot*tot);
-               count+=1;
-       }
-
-       /* averages */
-       moldyn->k_m=k_sum/count;
-       moldyn->p_m=p_sum/count;
-       moldyn->t_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,
-                                       moldyn->t_m);
-       printf("  mean square: %f %f %f\n",moldyn->dk2_m,
-                                          moldyn->dp2_m,
-                                          moldyn->dt2_m);
-
-       close(fd);
-
-       return 0;
-}
-
-int get_heat_capacity(t_moldyn *moldyn) {
-
-       double temp2,mass,ighc;
-       int i;
-
-       /* (temperature average)^2 */
-       temp2=2.0*moldyn->k_m*EV/(3.0*K_BOLTZMANN);
-       printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",temp2);
-       temp2*=temp2;
-
-       /* total mass */
-       mass=0.0;
-       for(i=0;i<moldyn->count;i++)
-               mass+=moldyn->atom[i].mass;
-
-       /* ideal gas contribution */
-       ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
-       printf("  ideal gas contribution: %f\n",ighc/mass*KILOGRAM/JOULE);
-
-       moldyn->c_v_nvt=moldyn->dp2_m*moldyn->count*moldyn->count*EV/(K_BOLTZMANN*temp2)+ighc;
-       moldyn->c_v_nvt/=mass;
-       moldyn->c_v_nve=ighc/(1.0-(moldyn->dp2_m*moldyn->count*moldyn->count*EV/(ighc*K_BOLTZMANN*temp2)));
-       moldyn->c_v_nve/=mass;
-
-       printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
-       printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
-
-       return 0;
-}

diff --git a/moldyn.h b/moldyn.h
index 484c456..2d1a438 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -77,6 +77,7 @@ typedef struct s_moldyn_schedule {
 /* moldyn main structure */
 typedef struct s_moldyn {
        int count;              /* total amount of atoms */
+       double mass;            /* total system mass */
        t_atom *atom;           /* pointer to the atoms */
 
        t_3dvec dim;            /* dimensions of the simulation volume */
@@ -144,12 +145,16 @@ typedef struct s_moldyn {
        double ekin;            /* kinetic energy */
 
        /* energy averages & fluctuations */
-       double k_m;
-       double p_m;
-       double t_m;
-       double dk2_m;           /* mean square kinetic energy fluctuations */
-       double dp2_m;           /* mean square potential energy fluctuations */
-       double dt2_m;           /* mean square fluctuations in total energy */
+       double k_sum;           /* sum of kinetic energy */
+       double v_sum;           /* sum of potential energy */
+       double k_mean;          /* average of kinetic energy */
+       double v_mean;          /* average of potential energy */
+       double k2_sum;          /* sum of kinetic energy squared */
+       double v2_sum;          /* sum of potential energy squared */
+       double k2_mean;         /* average of kinetic energy squared */
+       double v2_mean;         /* average of potential energy squared */
+       double dk2_mean;        /* mean square kinetic energy fluctuations */
+       double dv2_mean;        /* mean square potential energy fluctuations */
        
        /* response functions */
        double c_v_nve;         /* constant volume heat capacity (nve) */
@@ -403,10 +408,13 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 int destroy_atoms(t_moldyn *moldyn);
 
 int thermal_init(t_moldyn *moldyn,u8 equi_init);
+double total_mass_calc(t_moldyn *moldyn);
 double temperature_calc(t_moldyn *moldyn);
 double get_temperature(t_moldyn *moldyn);
 int scale_velocity(t_moldyn *moldyn,u8 equi_init);
 double pressure_calc(t_moldyn *moldyn);
+int energy_fluctuation_calc(t_moldyn *moldyn);
+int get_heat_capacity(t_moldyn *moldyn);
 double thermodynamic_pressure_calc(t_moldyn *moldyn);
 double get_pressure(t_moldyn *moldyn);
 int scale_volume(t_moldyn *moldyn);
@@ -443,7 +451,5 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a);
 int moldyn_bc_check(t_moldyn *moldyn);
 
 int get_line(int fd,char *line,int max);
-int calc_fluctuations(double start,double end,t_moldyn *moldyn);
-int get_heat_capacity(t_moldyn *moldyn);
 
 #endif

diff --git a/sic.c b/sic.c
index 07f32bc..426998a 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -337,10 +337,6 @@ return 0;
         * post processing the data
         */
 
-       /* response functions expressed by energy fluctuations */
-       calc_fluctuations(1000.0,2999.0,&md);
-       get_heat_capacity(&md);
-
        /* close */
        moldyn_shutdown(&md);