From fb951c04e522e4637618bf622fc67194c2a7b15f Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 17 Jan 2007 17:33:16 +0000 Subject: [PATCH] ... --- clean | 3 +- moldyn.c | 190 ++++++++++++++++++++++++++++++++++++++++++++++--------- moldyn.h | 7 +- sic.c | 29 +++++---- 4 files changed, 185 insertions(+), 44 deletions(-) diff --git a/clean b/clean index ca57788..1f359db 100755 --- a/clean +++ b/clean @@ -1 +1,2 @@ -rm -f $1/* > /dev/null 2>&1 +rm -rf $1 +mkdir -p $1 diff --git a/moldyn.c b/moldyn.c index 9db5cd9..94d596e 100644 --- a/moldyn.c +++ b/moldyn.c @@ -121,12 +121,15 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->vis.dim.z=z; } + moldyn->dv=0.0001*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; } @@ -279,6 +282,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, count=moldyn->count; /* how many atoms do we expect */ + if(type==CUBIC) new*=1; if(type==FCC) new*=4; if(type==DIAMOND) new*=8; @@ -294,6 +298,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, v3_zero(&origin); switch(type) { + case CUBIC: + ret=cubic_init(a,b,c,atom,&origin); + break; case FCC: ret=fcc_init(a,b,c,lc,atom,&origin); break; @@ -329,6 +336,15 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, return ret; } +/* cubic init */ +int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { + + int count; + t_3dvec r + +HIER WEITER ! +} + /* fcc lattice init */ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { @@ -507,6 +523,7 @@ double temperature_calc(t_moldyn *moldyn) { atom=moldyn->atom; + double_ekin=0; for(i=0;icount;i++) double_ekin+=atom[i].mass*v3_absolute_square(&(atom[i].v)); @@ -577,34 +594,153 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { return 0; } +double ideal_gas_law_pressure(t_moldyn *moldyn) { + + double p; + + p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume; +printf("temp = %f => ideal gas law pressure: %f\n",moldyn->t,p/ATM); + + return p; +} + double pressure_calc(t_moldyn *moldyn) { int i; - t_atom *atom; - double p1,p2,p=0; - - for(i=0;icount;i++) { - + double p1,p; + double v,h; + t_virial *virial; + v=0.0; + for(i=0;icount;i++) { + virial=&(moldyn->atom[i].virial); + v+=(virial->xx+virial->yy+virial->zz); } + h=moldyn->count*K_BOLTZMANN*moldyn->t; + + p=(h-ONE_THIRD*v); + p/=moldyn->volume; + p1=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt1); p1/=moldyn->volume; - p2=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt2); - p2/=moldyn->volume; + printf("debug: vt1=%f v=%f nkt=%f\n",moldyn->vt1,v,h); - printf("compare pressures: %f %f\n",p1/ATM,p2/ATM); + printf("compare pressures: %f %f %f\n",p1/ATM,p/ATM,h/moldyn->volume/ATM); return moldyn->p; } +double thermodynamic_pressure_calc(t_moldyn *moldyn) { + + t_3dvec dim,*tp; + double u,p; + double scale; + t_atom *store; + + tp=&(moldyn->tp); + 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 */ + u=moldyn->energy; + memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom)); + dim=moldyn->dim; + + /* derivative with respect to x direction */ + scale=1.0+moldyn->dv/(moldyn->dim.y*moldyn->dim.z); + scale_dim(moldyn,scale,TRUE,0,0); + scale_atoms(moldyn,scale,TRUE,0,0); + potential_force_calc(moldyn); + tp->x=(moldyn->energy-u)/moldyn->dv; + p=tp->x*tp->x; + + /* restore atomic configuration + dim */ + memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); + moldyn->dim=dim; + + /* derivative with respect to y direction */ + scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.z); + scale_dim(moldyn,scale,0,TRUE,0); + scale_atoms(moldyn,scale,0,TRUE,0); + potential_force_calc(moldyn); + tp->y=(moldyn->energy-u)/moldyn->dv; + p+=tp->y*tp->y; + + /* restore atomic configuration + dim */ + memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); + moldyn->dim=dim; + + /* derivative with respect to z direction */ + scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.y); + scale_dim(moldyn,scale,0,0,TRUE); + scale_atoms(moldyn,scale,0,0,TRUE); + potential_force_calc(moldyn); + tp->z=(moldyn->energy-u)/moldyn->dv; + p+=tp->z*tp->z; + + /* restore atomic configuration + dim */ + memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); + moldyn->dim=dim; + + printf("dU/dV komp addiert = %f\n",(tp->x+tp->y+tp->z)/ATM); + + scale=1.0+pow(moldyn->dv/moldyn->volume,ONE_THIRD); + + scale_dim(moldyn,scale,1,1,1); + scale_dim(moldyn,scale,1,1,1); + potential_force_calc(moldyn); + + printf("dU/dV einfach = %f\n",(moldyn->energy-u)/moldyn->dv/ATM); + + /* restore atomic configuration + dim */ + memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); + moldyn->dim=dim; + + /* restore energy */ + moldyn->energy=u; + + return sqrt(p); +} + double get_pressure(t_moldyn *moldyn) { return moldyn->p; } +int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { + + t_3dvec *dim; + + dim=&(moldyn->dim); + + if(x) dim->x*=scale; + if(y) dim->y*=scale; + if(z) dim->z*=scale; + + return 0; +} + +int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { + + int i; + t_3dvec *r; + + for(i=0;icount;i++) { + r=&(moldyn->atom[i].r); + if(x) r->x*=scale; + if(y) r->y*=scale; + if(z) r->z*=scale; + } + + return 0; +} + int scale_volume(t_moldyn *moldyn) { t_atom *atom; @@ -979,22 +1115,17 @@ int moldyn_integrate(t_moldyn *moldyn) { if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) scale_volume(moldyn); +printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM); + /* check for log & visualization */ -//double ax; -//double ao; -//double av; if(e) { if(!(i%e)) -//ao=sqrt(0.1/M_SI); -//ax=((0.28-0.25)*sqrt(3)*LC_SI/2)*cos(ao*i); -//av=ao*(0.28-0.25)*sqrt(3)*LC_SI/2*sin(ao*i); update_e_kin(moldyn); dprintf(moldyn->efd, "%f %f %f %f\n", moldyn->time,moldyn->ekin, moldyn->energy, get_total_energy(moldyn)); -//moldyn->atom[0].r.x,ax,av*av*M_SI,0.1*ax*ax,av*av*M_SI+0.1*ax*ax); } if(m) { if(!(i%m)) { @@ -1115,27 +1246,32 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset energy */ moldyn->energy=0.0; - moldyn->vt2=0.0; - - /* get energy and force of every atom */ + /* reset virial */ + moldyn->vt1=0.0; + + /* reset force, site energy and virial of every atom */ for(i=0;ixx=0.0; virial->yy=0.0; virial->zz=0.0; virial->xy=0.0; virial->xz=0.0; virial->yz=0.0; - moldyn->vt1=0.0; - + /* reset site energy */ itom[i].e=0.0; + } + + /* get energy,force and virial of every atom */ + for(i=0;ifunc1b(moldyn,&(itom[i])); @@ -1241,13 +1377,9 @@ printf("\n\n"); printf("\n\n"); #endif - moldyn->vt2=0.0; - for(i=0;ivt2-=v3_scalar_product(&(itom[i].r),&(itom[i].f)); - -//printf("compare: vt1: %f vt2: %f\n",moldyn->vt1,moldyn->vt2); - -//pressure_calc(moldyn); +temperature_calc(moldyn); +pressure_calc(moldyn); +ideal_gas_law_pressure(moldyn); return 0; } diff --git a/moldyn.h b/moldyn.h index 8254568..c910d7f 100644 --- a/moldyn.h +++ b/moldyn.h @@ -105,7 +105,9 @@ typedef struct s_moldyn { double t; /* actual temperature */ double p_ref; /* reference pressure */ - double p; /* actual pressure */ + double p; /* actual pressure (computed by virial) */ + t_3dvec tp; /* thermodynamic pressure dU/dV */ + double dv; /* dV for thermodynamic pressure calc */ /* pressure and temperature control (velocity/volume scaling) */ /* (t_tc in units of tau, p_tc in units of tau * isoth. compressib.) */ @@ -407,8 +409,11 @@ 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); +double thermodynamic_pressure_calc(t_moldyn *moldyn); double get_pressure(t_moldyn *moldyn); int scale_volume(t_moldyn *moldyn); +int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z); +int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z); double get_e_kin(t_moldyn *moldyn); double update_e_kin(t_moldyn *moldyn); diff --git a/sic.c b/sic.c index a890270..5cfbb03 100644 --- a/sic.c +++ b/sic.c @@ -17,8 +17,16 @@ int hook(void *moldyn,void *hook_params) { md=moldyn; - /* decrease temperature in every hook */ - set_temperature(md,md->t_ref-100.0); + /* switch to direct scaling in first hook */ + if(md->schedule.count==0) + set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); + /* switch off temp scaling in second hook */ + if(md->schedule.count==1) + set_pt_scale(md,0,0,0,0); + + + //set_temperature(md,md->t_ref-100.0); + return 0; } @@ -153,29 +161,24 @@ int main(int argc,char **argv) { /* set p/t scaling */ //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0, // T_SCALE_BERENDSEN,100.0); - set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); + set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0); //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,0,0); /* initial thermal fluctuations of particles (in equilibrium) */ thermal_init(&md,TRUE); /* create the simulation schedule */ - moldyn_add_schedule(&md,100,1.0); - moldyn_add_schedule(&md,100,1.0); - moldyn_add_schedule(&md,100,1.0); - moldyn_add_schedule(&md,100,1.0); - moldyn_add_schedule(&md,100,1.0); - moldyn_add_schedule(&md,100,1.0); - moldyn_add_schedule(&md,100,1.0); - moldyn_add_schedule(&md,100,1.0); + moldyn_add_schedule(&md,1001,3.0); + moldyn_add_schedule(&md,501,1.0); + moldyn_add_schedule(&md,501,1.0); /* schedule hook function */ moldyn_set_schedule_hook(&md,&hook,NULL); /* activate logging */ moldyn_set_log_dir(&md,argv[1]); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,1); - moldyn_set_log(&md,VISUAL_STEP,10); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,5); + moldyn_set_log(&md,VISUAL_STEP,50); /* * let's do the actual md algorithm now -- 2.20.1