X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=cd2a80309c6cf9de56e4b23132debbe9e46b02cd;hb=0fbd532d8ddce9848df592ed586ffba439369284;hp=c130ef442b838175272e62f375daa4fbaa8e199b;hpb=58fd691b276fbe87593036714f26dbfe7486cbeb;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index c130ef4..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,110 +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; - 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; - } - - /* first calc the averages */ - 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; - count+=1; - } - - moldyn->p_m=p_sum/count; - moldyn->k_m=k_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(timeend) 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; - - 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; -}