From e2c5d8b0eb4ce6faeb48830634eef522dbdb52b0 Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 10 Jul 2007 17:27:01 +0000 Subject: [PATCH] safety checkin, gergo here! --- fluctuation_calc.c | 1 - moldyn.c | 109 ++++++++++++++++++++++----------------------- moldyn.h | 24 +++++----- potentials/albe.c | 2 +- sic.c | 10 +---- 5 files changed, 67 insertions(+), 79 deletions(-) diff --git a/fluctuation_calc.c b/fluctuation_calc.c index 0ef0920..469371d 100644 --- a/fluctuation_calc.c +++ b/fluctuation_calc.c @@ -100,7 +100,6 @@ int main(int argc,char **argv) { } } -printf(" count = %d\n",count); de2=e2/count-e*e/(count*count); printf("--> fluctuation [eV/atom]: %.20f\n",de2); if(argc==8) { diff --git a/moldyn.c b/moldyn.c index 16c811e..e0881de 100644 --- a/moldyn.c +++ b/moldyn.c @@ -218,10 +218,10 @@ int set_potential_params(t_moldyn *moldyn,void *params) { return 0; } -int set_mean_skip(t_moldyn *moldyn,int skip) { +int set_avg_skip(t_moldyn *moldyn,int skip) { printf("[moldyn] skip %d steps before starting average calc\n",skip); - moldyn->mean_skip=skip; + moldyn->avg_skip=skip; return 0; } @@ -737,12 +737,6 @@ double temperature_calc(t_moldyn *moldyn) { moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); - if(moldyn->total_stepsmean_skip) - return 0; - - moldyn->t_sum+=moldyn->t; - moldyn->mean_t=moldyn->t_sum/(moldyn->total_steps+1-moldyn->mean_skip); - return moldyn->t; } @@ -824,7 +818,7 @@ double pressure_calc(t_moldyn *moldyn) { /* * PV = NkT + - * W = 1/3 sum_i f_i r_i + * with W = 1/3 sum_i f_i r_i (- skipped!) * virial = sum_i f_i r_i * * => P = (2 Ekin + virial) / (3V) @@ -836,19 +830,16 @@ double pressure_calc(t_moldyn *moldyn) { v+=(virial->xx+virial->yy+virial->zz); } - /* virial sum and mean virial */ - moldyn->virial_sum+=v; - if(moldyn->total_steps>=moldyn->mean_skip) - moldyn->mean_v=moldyn->virial_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); - - /* assume up to date kinetic energy */ - moldyn->p=2.0*moldyn->ekin+moldyn->mean_v; - moldyn->p/=(3.0*moldyn->volume); - if(moldyn->total_steps>=moldyn->mean_skip) { + /* virial sum and average virial */ + if(moldyn->total_steps>=moldyn->avg_skip) { + moldyn->virial_sum+=v; + moldyn->virial_avg=moldyn->virial_sum/ + (moldyn->total_steps+1-moldyn->avg_skip); + moldyn->p=2.0*moldyn->k_avg+moldyn->virial_avg; + moldyn->p/=(3.0*moldyn->volume); moldyn->p_sum+=moldyn->p; - moldyn->mean_p=moldyn->p_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); + moldyn->p_avg=moldyn->p_sum/ + (moldyn->total_steps+1-moldyn->avg_skip); } /* pressure from 'absolute coordinates' virial */ @@ -856,37 +847,45 @@ double pressure_calc(t_moldyn *moldyn) { v=virial->xx+virial->yy+virial->zz; moldyn->gp=2.0*moldyn->ekin+v; moldyn->gp/=(3.0*moldyn->volume); - if(moldyn->total_steps>=moldyn->mean_skip) { + if(moldyn->total_steps>=moldyn->avg_skip) { moldyn->gp_sum+=moldyn->gp; - moldyn->mean_gp=moldyn->gp_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); + moldyn->gp_avg=moldyn->gp_sum/ + (moldyn->total_steps+1-moldyn->avg_skip); } return moldyn->p; } -int energy_fluctuation_calc(t_moldyn *moldyn) { +int average_and_fluctuation_calc(t_moldyn *moldyn) { - if(moldyn->total_stepsmean_skip) + if(moldyn->total_stepsavg_skip) return 0; - /* assume up to date energies */ + /* assume up to date energies, temperature, pressure etc */ /* kinetic energy */ moldyn->k_sum+=moldyn->ekin; moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin); - moldyn->k_mean=moldyn->k_sum/(moldyn->total_steps+1-moldyn->mean_skip); - moldyn->k2_mean=moldyn->k2_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); - moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean); + moldyn->k_avg=moldyn->k_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->k2_avg=moldyn->k2_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg); /* potential energy */ moldyn->v_sum+=moldyn->energy; moldyn->v2_sum+=(moldyn->energy*moldyn->energy); - moldyn->v_mean=moldyn->v_sum/(moldyn->total_steps+1-moldyn->mean_skip); - moldyn->v2_mean=moldyn->v2_sum/ - (moldyn->total_steps+1-moldyn->mean_skip); - moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean); + moldyn->v_avg=moldyn->v_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->v2_avg=moldyn->v2_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg); + + /* temperature */ + moldyn->t_sum+=moldyn->t; + moldyn->t_avg=moldyn->t_sum/(moldyn->total_steps+1-moldyn->avg_skip); + + /* virial */ + + + /* pressure */ + return 0; } @@ -896,13 +895,13 @@ int get_heat_capacity(t_moldyn *moldyn) { double temp2,ighc; /* averages needed for heat capacity calc */ - if(moldyn->total_stepsmean_skip) + if(moldyn->total_stepsavg_skip) return 0; /* (temperature average)^2 */ - temp2=moldyn->mean_t*moldyn->mean_t; + temp2=moldyn->t_avg*moldyn->t_avg; printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n", - moldyn->mean_t); + moldyn->t_avg); /* ideal gas contribution */ ighc=3.0*moldyn->count*K_BOLTZMANN/2.0; @@ -910,16 +909,16 @@ int get_heat_capacity(t_moldyn *moldyn) { 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->dv2_avg/(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=ighc/(1.0-(moldyn->dv2_avg/(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); -printf(" --> sim: %f experimental: %f\n",moldyn->dv2_mean,1.5*moldyn->count*K_B2*moldyn->mean_t*moldyn->mean_t*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM))); +printf(" --> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM))); return 0; } @@ -1404,9 +1403,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); + average_and_fluctuation_calc(moldyn); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1436,15 +1433,15 @@ return 0; if(!(i%p)) { dprintf(moldyn->pfd, "%f %f %f %f %f\n",moldyn->time, - moldyn->p/BAR,moldyn->mean_p/BAR, - moldyn->gp/BAR,moldyn->mean_gp/BAR); + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->gp/BAR,moldyn->gp_avg/BAR); } } if(t) { if(!(i%t)) { dprintf(moldyn->tfd, "%f %f %f\n", - moldyn->time,moldyn->t,moldyn->mean_t); + moldyn->time,moldyn->t,moldyn->t_avg); } } if(s) { @@ -1472,13 +1469,12 @@ return 0; if(!(i%10)) { printf("\rsched: %d, steps: %d, T: %f, P: %f %f V: %f", sched->count,i, - moldyn->mean_t, - moldyn->mean_p/BAR, - moldyn->mean_gp/BAR, + moldyn->t_avg, + //moldyn->p_avg/BAR, + moldyn->p/BAR, + moldyn->gp_avg/BAR, moldyn->volume); fflush(stdout); -printf("\n"); -get_heat_capacity(moldyn); } /* increase absolute time */ @@ -1489,11 +1485,12 @@ get_heat_capacity(moldyn); /* check for hooks */ if(sched->count+1total_sched) - if(sched->hook) + if(sched->hook) { + printf("\n ## schedule hook %d/%d start ##\n", + sched->count+1,sched->total_sched); sched->hook(moldyn,sched->hook_params); - - /* get a new info line */ - printf("\n"); + printf(" ## schedule hook end ##\n"); + } } diff --git a/moldyn.h b/moldyn.h index 868db62..212db14 100644 --- a/moldyn.h +++ b/moldyn.h @@ -102,25 +102,25 @@ typedef struct s_moldyn { t_linkcell lc; /* linked cell list interface */ - int mean_skip; /* amount of steps without average calc */ + int avg_skip; /* amount of steps without average calc */ double t_ref; /* reference temperature */ double t; /* actual temperature */ double t_sum; /* sum over all t */ - double mean_t; /* mean value of t */ + double t_avg; /* average value of t */ t_virial virial; /* global virial (absolute coordinates) */ double gp; /* pressure computed from global virial */ double gp_sum; /* sum over all gp */ - double mean_gp; /* mean value of gp */ + double gp_avg; /* average value of gp */ - double mean_v; /* mean of virial */ + double virial_avg; /* average of virial */ double virial_sum; /* sum over all calculated virials */ double p_ref; /* reference pressure */ double p; /* actual pressure (computed by virial) */ double p_sum; /* sum over all p */ - double mean_p; /* mean value of p */ + double p_avg; /* average value of p */ t_3dvec tp; /* thermodynamic pressure dU/dV */ double dv; /* dV for thermodynamic pressure calc */ @@ -149,14 +149,14 @@ typedef struct s_moldyn { /* energy averages & fluctuations */ 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 k_avg; /* average of kinetic energy */ + double v_avg; /* 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 */ + double k2_avg; /* average of kinetic energy squared */ + double v2_avg; /* average of potential energy squared */ + double dk2_avg; /* mean square kinetic energy fluctuations */ + double dv2_avg; /* mean square potential energy fluctuations */ /* response functions */ double c_v_nve; /* constant volume heat capacity (nve) */ @@ -395,7 +395,7 @@ int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func); int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func); int set_potential_params(t_moldyn *moldyn,void *params); -int set_mean_skip(t_moldyn *moldyn,int skip); +int set_avg_skip(t_moldyn *moldyn,int skip); int moldyn_set_log_dir(t_moldyn *moldyn,char *dir); int moldyn_set_report(t_moldyn *moldyn,char *author,char *title); diff --git a/potentials/albe.c b/potentials/albe.c index a79f58d..338fb03 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -450,7 +450,7 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, /* virial */ //v3_scale(&force,&force,-1.0); - //virial_calc(ai,&force,&dist_ik); + virial_calc(ai,&force,&dist_ik); /* increase k counter */ exchange->kcount++; diff --git a/sic.c b/sic.c index 5bdfeea..f8e6ce1 100644 --- a/sic.c +++ b/sic.c @@ -35,14 +35,6 @@ int hook(void *moldyn,void *hook_params) { md=moldyn; -// vortrag -set_temperature(moldyn,(4-md->schedule.count)*1000.0); -set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); - -return 0; - - printf("\nschedule hook: "); - if(!(md->schedule.count%2)) { /* add carbon at random place, and enable t scaling */ for(j=0;j