From 0fbd532d8ddce9848df592ed586ffba439369284 Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 3 Jul 2007 16:04:05 +0000 Subject: [PATCH] changed energy fluctuation calc and heat cap calc --- moldyn.c | 165 ++++++++++++++++++++++++------------------------------- moldyn.h | 22 +++++--- sic.c | 4 -- 3 files changed, 86 insertions(+), 105 deletions(-) diff --git a/moldyn.c b/moldyn.c index 4c43f2e..cd2a803 100644 --- 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;icount;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(timeend) 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;icount;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 --- 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 --- 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); -- 2.20.1