From a32468230b319b32819f1b20fd28aa9659574d45 Mon Sep 17 00:00:00 2001 From: hackbard Date: Sat, 5 Apr 2008 20:19:00 +0200 Subject: [PATCH] implemented basic p ctrl stuff + video with 13 fps --- config.h | 33 ++++++++------ moldyn.c | 112 ++++++++++++++++++++++++++++++---------------- moldyn.h | 21 ++++++--- potentials/albe.c | 13 +++--- ppm2avi | 1 + sic.c | 36 +++++++++++---- 6 files changed, 144 insertions(+), 72 deletions(-) diff --git a/config.h b/config.h index 9f04dbd..b69f41a 100644 --- a/config.h +++ b/config.h @@ -3,11 +3,13 @@ * */ +#include "moldyn.h" + // simulation volume -#define LCNTX 9 -#define LCNTY 9 -#define LCNTZ 9 +#define LCNTX 31 +#define LCNTY 31 +#define LCNTZ 31 // initial lattice @@ -19,16 +21,17 @@ //#define T_SCALE_TAU 10.0 #define T_SCALE_TAU 100.0 +#define P_SCALE_TAU (0.01/(100*GPA)) // prerun -#define PRERUN 0 +#define PRERUN 500 #define PRE_TAU 1.0 // insertrun -#define INS_RUNS 1 -#define INS_ATOMS 1 +#define INS_RUNS 600 +#define INS_ATOMS 10 #define INS_CARBON /* comment for silicon */ @@ -42,11 +45,11 @@ #define INS_BRAND 0 #endif -//#define INS_RAND // random nsert +#define INS_RAND // random nsert //#define INS_HEXA // hexagonal interstitial position //#define INS_TETRA // tetrahedral interstitial position //#define INS_110DB // 110 dumbbell interstitial position -#define INS_001DB // 001 dumbbell interstitial position +//#define INS_001DB // 001 dumbbell interstitial position //#define INS_USER // user defined insertion location #define INS_UX -0.25 #define INS_UY -0.25 @@ -63,17 +66,20 @@ #define INS_LENY (9*ALBE_LC_SI) #define INS_LENZ (9*ALBE_LC_SI) #define INS_OFFSET (ALBE_LC_SI/8.0) +#define INS_DYNAMIC_LEN // with p-ctrl we have a varying volume #define INS_R_C 1.5 #define INS_DELTA_TC 1.0 +#define INS_DELTA_PC (1.0*BAR) #define INS_RELAX 100 #define INS_TAU 1.0 // postrun -#define POST_RUNS 1 -#define POST_DELTA_TC 0.01 -#define POST_DT -0.0 +#define POST_RUNS 430 +#define POST_DELTA_TC 1.0 +#define POST_DELTA_PC (1.0*BAR) +#define POST_DT -1.0 #define POST_RELAX 100 #define POST_TAU 1.0 @@ -83,7 +89,8 @@ #define LOG_E 10 #define LOG_T 10 #define LOG_P 10 -#define LOG_S 20 -#define LOG_V 20 +#define LOG_V 10 +#define LOG_S 2000 +#define LOG_A 2000 #define AVG_SKIP 0 diff --git a/moldyn.c b/moldyn.c index 7496223..91fead7 100644 --- a/moldyn.c +++ b/moldyn.c @@ -191,15 +191,12 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->vis.dim.z=z; } - moldyn->dv=0.000001*moldyn->volume; - printf("[moldyn] dimensions in A and A^3 respectively:\n"); printf(" x: %f\n",moldyn->dim.x); printf(" y: %f\n",moldyn->dim.y); printf(" z: %f\n",moldyn->dim.z); printf(" volume: %f\n",moldyn->volume); printf(" visualize simulation box: %s\n",visualize?"yes":"no"); - printf(" delta volume (pressure calc): %f\n",moldyn->dv); return 0; } @@ -370,12 +367,25 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { dprintf(moldyn->tfd,"# temperature log file\n"); printf("temperature (%d)\n",timer); break; + case LOG_VOLUME: + moldyn->vwrite=timer; + snprintf(filename,127,"%s/volume",moldyn->vlsdir); + moldyn->vfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->vfd<0) { + perror("[moldyn] volume log file\n"); + return moldyn->vfd; + } + dprintf(moldyn->vfd,"# volume log file\n"); + printf("volume (%d)\n",timer); + break; case SAVE_STEP: moldyn->swrite=timer; printf("save file (%d)\n",timer); break; case VISUAL_STEP: - moldyn->vwrite=timer; + moldyn->awrite=timer; ret=visual_init(moldyn,moldyn->vlsdir); if(ret<0) { printf("[moldyn] visual init failure\n"); @@ -916,16 +926,26 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) { double virial_sum(t_moldyn *moldyn) { int i; - double v; t_virial *virial; /* virial (sum over atom virials) */ - v=0.0; + moldyn->virial=0.0; + moldyn->vir.xx=0.0; + moldyn->vir.yy=0.0; + moldyn->vir.zz=0.0; + moldyn->vir.xy=0.0; + moldyn->vir.xz=0.0; + moldyn->vir.yz=0.0; for(i=0;icount;i++) { virial=&(moldyn->atom[i].virial); - v+=(virial->xx+virial->yy+virial->zz); + moldyn->virial+=(virial->xx+virial->yy+virial->zz); + moldyn->vir.xx+=virial->xx; + moldyn->vir.yy+=virial->yy; + moldyn->vir.zz+=virial->zz; + moldyn->vir.xy+=virial->xy; + moldyn->vir.xz+=virial->xz; + moldyn->vir.yz+=virial->yz; } - moldyn->virial=v; /* global virial (absolute coordinates) */ virial=&(moldyn->gvir); @@ -982,6 +1002,7 @@ int average_reset(t_moldyn *moldyn) { /* pressure */ moldyn->p_sum=0.0; moldyn->gp_sum=0.0; + moldyn->tp_sum=0.0; return 0; } @@ -1026,6 +1047,8 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) { moldyn->p_avg=moldyn->p_sum/denom; moldyn->gp_sum+=moldyn->gp; moldyn->gp_avg=moldyn->gp_sum/denom; + moldyn->tp_sum+=moldyn->tp; + moldyn->tp_avg=moldyn->tp_sum/denom; return 0; } @@ -1067,8 +1090,9 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { t_3dvec dim; //t_3dvec *tp; - double u_up,u_down,dv; - double scale,p; + double h,dv; + double y0,y1; + double su,sd; t_atom *store; /* @@ -1078,54 +1102,56 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { * */ - scale=0.00001; - dv=8*scale*scale*scale*moldyn->volume; - + /* store atomic configuration + dimension */ store=malloc(moldyn->count*sizeof(t_atom)); if(store==NULL) { printf("[moldyn] allocating store mem failed\n"); return -1; } - - /* save unscaled potential energy + atom/dim configuration */ memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom)); dim=moldyn->dim; + /* x1, y1 */ + sd=0.00001; + h=(1.0-sd)*(1.0-sd)*(1.0-sd); + su=pow(2.0-h,ONE_THIRD)-1.0; + dv=(1.0-h)*moldyn->volume; + /* scale up dimension and atom positions */ - scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); - scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - u_up=moldyn->energy; + y1=moldyn->energy; /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; /* scale down dimension and atom positions */ - scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); - scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE); + scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE); link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - u_down=moldyn->energy; + y0=moldyn->energy; /* calculate pressure */ - p=-(u_up-u_down)/dv; -printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR); + moldyn->tp=-(y1-y0)/(2.0*dv); - /* restore atomic configuration + dim */ + /* restore atomic configuration */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - - /* restore energy */ - potential_force_calc(moldyn); - link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); + //potential_force_calc(moldyn); - return p; + /* free store buffer */ + if(store) + free(store); + + return moldyn->tp; } double get_pressure(t_moldyn *moldyn) { @@ -1186,13 +1212,12 @@ int scale_volume(t_moldyn *moldyn) { /* scaling factor */ if(moldyn->pt_scale&P_SCALE_BERENDSEN) { - scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc; + scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc; scale=pow(scale,ONE_THIRD); } else { scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD); } -moldyn->debug=scale; /* scale the atoms and dimensions */ scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); @@ -1538,7 +1563,7 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) { int moldyn_integrate(t_moldyn *moldyn) { int i; - unsigned int e,m,s,v,p,t; + unsigned int e,m,s,v,p,t,a; t_3dvec momentum; t_moldyn_schedule *sched; t_atom *atom; @@ -1560,6 +1585,7 @@ int moldyn_integrate(t_moldyn *moldyn) { m=moldyn->mwrite; s=moldyn->swrite; v=moldyn->vwrite; + a=moldyn->awrite; p=moldyn->pwrite; t=moldyn->twrite; @@ -1621,6 +1647,9 @@ int moldyn_integrate(t_moldyn *moldyn) { temperature_calc(moldyn); virial_sum(moldyn); pressure_calc(moldyn); + //thermodynamic_pressure_calc(moldyn); + + /* calculate fluctuations + averages */ average_and_fluctuation_calc(moldyn); /* p/t scaling */ @@ -1650,9 +1679,10 @@ int moldyn_integrate(t_moldyn *moldyn) { if(p) { if(!(moldyn->total_steps%p)) { dprintf(moldyn->pfd, - "%f %f %f %f %f\n",moldyn->time, + "%f %f %f %f %f %f %f\n",moldyn->time, moldyn->p/BAR,moldyn->p_avg/BAR, - moldyn->gp/BAR,moldyn->gp_avg/BAR); + moldyn->gp/BAR,moldyn->gp_avg/BAR, + moldyn->tp/BAR,moldyn->tp_avg/BAR); } } if(t) { @@ -1662,6 +1692,12 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->time,moldyn->t,moldyn->t_avg); } } + if(v) { + if(!(moldyn->total_steps%v)) { + dprintf(moldyn->vfd, + "%f %f\n",moldyn->time,moldyn->volume); + } + } if(s) { if(!(moldyn->total_steps%s)) { snprintf(dir,128,"%s/s-%07.f.save", @@ -1677,8 +1713,8 @@ int moldyn_integrate(t_moldyn *moldyn) { close(fd); } } - if(v) { - if(!(moldyn->total_steps%v)) { + if(a) { + if(!(moldyn->total_steps%a)) { visual_atoms(moldyn); } } @@ -1688,10 +1724,10 @@ int moldyn_integrate(t_moldyn *moldyn) { /* get current time */ gettimeofday(&t2,NULL); -printf("\rsched:%d, steps:%d/%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)", +printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)", sched->count,i,moldyn->total_steps, moldyn->t,moldyn->t_avg, - moldyn->p_avg/BAR,moldyn->gp_avg/BAR, + moldyn->p/BAR,moldyn->p_avg/BAR, moldyn->volume, (int)(t2.tv_sec-t1.tv_sec)); diff --git a/moldyn.h b/moldyn.h index 3867fc0..83b3a86 100644 --- a/moldyn.h +++ b/moldyn.h @@ -137,7 +137,8 @@ typedef struct s_moldyn { double gp_sum; /* sum over all gp */ double gp_avg; /* average value of gp */ - double virial; /* actual virial */ + t_virial vir; /* actual virial */ + double virial; double virial_sum; /* sum over all calculated virials */ double virial_avg; /* average of virial */ @@ -146,8 +147,10 @@ typedef struct s_moldyn { double p_sum; /* sum over all p */ double p_avg; /* average value of p */ - t_3dvec tp; /* thermodynamic pressure dU/dV */ - double dv; /* dV for thermodynamic pressure calc */ + double tp; /* thermodynamic pressure dU/dV */ + double tp_sum; /* sum over dU/dV pressure */ + double tp_avg; /* average value of dU/dV pressure */ + int tp_cnt; /* how often to do thermodynamic p calc */ /* pressure and temperature control (velocity/volume scaling) */ /* (t_tc in units of tau, p_tc in units of tau * isoth. compressib.) */ @@ -198,7 +201,9 @@ typedef struct s_moldyn { int pfd; /* fd for pressure log */ unsigned int twrite; /* how often to log temperature */ int tfd; /* fd for temperature log */ - unsigned int vwrite; /* how often to visualize atom information */ + unsigned int vwrite; /* how often to log volume */ + int vfd; /* fd for volume log */ + unsigned int awrite; /* how often to visualize atom information */ unsigned int swrite; /* how often to create a save file */ int rfd; /* report file descriptor */ char rtitle[64]; /* report title */ @@ -252,6 +257,7 @@ typedef struct s_moldyn { #define KILOGRAM (1.0/AMU) /* amu */ #define NEWTON (METER*KILOGRAM/(SECOND*SECOND)) /* A amu / fs^2 */ #define PASCAL (NEWTON/(METER*METER)) /* N / A^2 */ +#define GPA (1e9*PASCAL) /* N / A^2 */ #define BAR ((1.0e5*PASCAL)) /* N / A^2 */ #define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) /* NA/K */ #define K_B2 (K_BOLTZMANN*K_BOLTZMANN) /* (NA)^2/K^2 */ @@ -274,9 +280,10 @@ typedef struct s_moldyn { #define LOG_TOTAL_MOMENTUM 0x02 #define LOG_PRESSURE 0x04 #define LOG_TEMPERATURE 0x08 -#define SAVE_STEP 0x10 -#define VISUAL_STEP 0x20 -#define CREATE_REPORT 0x40 +#define LOG_VOLUME 0x10 +#define SAVE_STEP 0x20 +#define VISUAL_STEP 0x40 +#define CREATE_REPORT 0x80 #define TRUE 1 #define FALSE 0 diff --git a/potentials/albe.c b/potentials/albe.c index 677d6e3..d157bae 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -275,11 +275,17 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { db=-0.5*b/(1.0+exchange->zeta_ij); } - /* force contribution */ + /* force contribution for atom i */ scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism v3_scale(&force,&(exchange->dist_ij),scale); v3_add(&(ai->f),&(ai->f),&force); - v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij + + /* force contribution for atom j */ + v3_scale(&force,&force,-1.0); // dri rij = - drj rij + v3_add(&(aj->f),&(aj->f),&force); + + /* virial */ + virial_calc(aj,&force,&(exchange->dist_ij)); #ifdef DEBUG if(moldyn->time>DSTART&&moldyn->timetime>DSTART&&moldyn->timedist_ij)); - /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ exchange->pre_dzeta=0.5*f_a*f_c*db; diff --git a/ppm2avi b/ppm2avi index 6c5d9e5..74c8370 100755 --- a/ppm2avi +++ b/ppm2avi @@ -3,6 +3,7 @@ #done MFOPTS="-mf fps=13:type=png" +#MFOPTS="-mf fps=25:type=png" OUTFILE=md.avi [ -n "$2" ] && OUTFILE=$2 diff --git a/sic.c b/sic.c index 8901104..dea4d4b 100644 --- a/sic.c +++ b/sic.c @@ -75,8 +75,8 @@ int insert_atoms(t_moldyn *moldyn) { // 001 dumbbell r.x=(-0.5+0.25)*ALBE_LC_SI; r.y=(-0.5+0.25)*ALBE_LC_SI; - r.z=(-0.5+0.25+0.125)*ALBE_LC_SI; - moldyn->atom[4372].r.z=(-0.5+0.25-0.125)*ALBE_LC_SI; + r.z=(-0.1)*ALBE_LC_SI; + moldyn->atom[4372].r.z=(-0.4)*ALBE_LC_SI; #endif #ifdef INS_USER // 001 dumbbell @@ -86,9 +86,18 @@ int insert_atoms(t_moldyn *moldyn) { #endif #ifdef INS_RAND // random +#ifdef INS_DYNAMIC_LEN + r.x=(rand_get_double(&(moldyn->random))-0.5)*\ + moldyn->dim.x; + r.y=(rand_get_double(&(moldyn->random))-0.5)*\ + moldyn->dim.y; + r.z=(rand_get_double(&(moldyn->random))-0.5)*\ + moldyn->dim.z; +#else r.x=(rand_get_double(&(moldyn->random))-0.5)*INS_LENX; r.y=(rand_get_double(&(moldyn->random))-0.5)*INS_LENY; r.z=(rand_get_double(&(moldyn->random))-0.5)*INS_LENZ; +#endif #endif // offset r.x+=INS_OFFSET; @@ -130,6 +139,7 @@ int sic_hook(void *moldyn,void *hook_params) { int steps; double tau; double dt; + double dp; hp=hook_params; md=moldyn; @@ -139,7 +149,8 @@ int sic_hook(void *moldyn,void *hook_params) { /* switch on t scaling */ if(md->schedule.count==0) - set_pt_scale(md,0,0,T_SCALE_BERENDSEN,T_SCALE_TAU); + set_pt_scale(md,P_SCALE_BERENDSEN,P_SCALE_TAU, + T_SCALE_BERENDSEN,T_SCALE_TAU); /* my lousy state machine ! */ @@ -169,9 +180,12 @@ insert: /* check temperature */ dt=md->t_avg-md->t_ref; + dp=md->p_avg-md->p_ref; if(dt<0) dt=-dt; - if(dt>INS_DELTA_TC) + if(dp<0) + dp=-dp; + if((dt>INS_DELTA_TC)|(dp>INS_DELTA_PC)) goto addsched; /* immediately go on if no job is to be done */ @@ -196,9 +210,12 @@ postrun: /* check temperature */ dt=md->t_avg-md->t_ref; + dp=md->p_avg-md->p_ref; if(dt<0) dt=-dt; - if(dt>POST_DELTA_TC) + if(dp<0) + dp=-dp; + if((dt>POST_DELTA_TC)|(dp>POST_DELTA_PC)) goto addsched; /* immediately return if no job is to be done */ @@ -416,7 +433,7 @@ int main(int argc,char **argv) { 0,LCNTX,LCNTY,LCNTZ,&r); r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x; create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C, - ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB|ATOM_ATTR_VB, 1,LCNTX,LCNTY,LCNTZ,&r); #else r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x; @@ -467,14 +484,14 @@ int main(int argc,char **argv) { /* set temperature & pressure */ set_temperature(&md,atof(argv[2])+273.0); - set_pressure(&md,BAR); + set_pressure(&md,0.0); /* set amount of steps to skip before average calc */ set_avg_skip(&md,AVG_SKIP); /* set p/t scaling */ //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); - //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001, + //set_pt_scale(&md,P_SCALE_BERENDSEN,0.01/(100*GPA), // T_SCALE_BERENDSEN,100.0); //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0); //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0); @@ -499,7 +516,8 @@ int main(int argc,char **argv) { moldyn_set_log(&md,LOG_TOTAL_ENERGY,LOG_E); moldyn_set_log(&md,LOG_TEMPERATURE,LOG_T); moldyn_set_log(&md,LOG_PRESSURE,LOG_P); - moldyn_set_log(&md,VISUAL_STEP,LOG_V); + moldyn_set_log(&md,LOG_VOLUME,LOG_V); + moldyn_set_log(&md,VISUAL_STEP,LOG_A); moldyn_set_log(&md,SAVE_STEP,LOG_S); moldyn_set_log(&md,CREATE_REPORT,0); -- 2.20.1