From d7f67c88195ab155f2737e57cc5e81973d3feb0c Mon Sep 17 00:00:00 2001 From: hackbard Date: Mon, 28 Apr 2008 19:17:06 +0200 Subject: [PATCH] safety checkin --- config.default | 4 +- mdrun.c | 433 ++++++++++++++++++++++++++++++++++++++++++++++++- mdrun.h | 82 +++++++++- moldyn.c | 86 ++++++++++ moldyn.h | 9 +- 5 files changed, 598 insertions(+), 16 deletions(-) diff --git a/config.default b/config.default index f221555..0bd26d2 100644 --- a/config.default +++ b/config.default @@ -24,7 +24,7 @@ pressure 0.0 ## ensemble control ## -eattrib 100.0 100.0 +ectrl 100.0 100.0 ## equi ctrl ## @@ -38,7 +38,7 @@ element1 silicon #element2 carbon fill lc 9 9 9 -aattrib all h +aattr all h ## initial simulation run ## diff --git a/mdrun.c b/mdrun.c index 4a2822e..6f7a2f0 100644 --- a/mdrun.c +++ b/mdrun.c @@ -62,6 +62,75 @@ int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) { return 0; } +del_stages(t_mdrun *mdrun) { + + t_list *sl; + t_stage *stage; + + sl=mdrun->stage; + + list_reset_f(sl); + + if(sl->start==NULL) + return 0; + + do { + stage=sl->current->data; + free(stage->params); + free(stage); + } while(list_next_f(sl)!=L_NO_NEXT_ELEMENT); + + return 0; +} + +add_stage(t_mdrun *mdrun,u8 type,void *params) { + + int psize; + + t_stage *stage; + + switch(type) { + case STAGE_INSERT_ATOMS: + psize=sizeof(t_insert_atoms_params); + break; + case STAGE_CONTINUE: + psize=sizeof(t_continue_params); + break; + case STAGE_ANNEAL: + psize=sizeof(t_anneal_params); + break; + case STAGE_CHAATTR: + psize=sizeof(t_chaattr_params); + break; + case STAGE_CHSATTR: + psize=sizeof(t_chsattr_params); + break; + default: + + } + + stage=malloc(sizeof(t_stage)); + if(stage==NULL) { + perror("[mdrun] malloc (add stage)"); + return -1; + } + + stage->type=type; + stage->executed=FALSE; + + stage->params=malloc(psize); + if(stage->params==NULL) { + perror("[mdrun] malloc (stage params)"); + return -1; + } + + memcpy(stage->params,params,psize); + + list_add_immediate_f(mdrun->stage,stage); + + return 0; +} + int mdrun_parse_config(t_mdrun *mdrun) { int fd,ret; @@ -70,6 +139,13 @@ int mdrun_parse_config(t_mdrun *mdrun) { char *wptr; char word[16][32]; int wcnt; + int i; + + t_insert_atoms_params iap; + t_continue_params cp; + t_anneal_params ap; + t_chaattr_params cap; + t_chsattr_params csp; /* open config file */ fd=open(mdrun->cfile,O_RDONLY); @@ -186,17 +262,51 @@ int mdrun_parse_config(t_mdrun *mdrun) { return -1; } } - else if(!strncmp(word[0],"filllc",6)) { - mdrun->lx=atoi(word[1]); - mdrun->ly=atoi(word[2]); - mdrun->lz=atoi(word[3]); + else if(!strncmp(word[0],"fill",6)) { + // only lc mode by now + mdrun->lx=atoi(word[2]); + mdrun->ly=atoi(word[3]); + mdrun->lz=atoi(word[4]); } - else if(!strncmp(word[0],"aattrib",7)) { + else if(!strncmp(word[0],"aattr",5)) { // for aatrib line we need a special stage // containing one schedule of 0 loops ... - printf("%s aattrib: %s\n",ME,word[1]); + if(!strncmp(word[1],"all",3)) { + cap.type= + + + HIER WEITER + } + for(i=0;ip_tau=atof(word[1]); mdrun->t_tau=atof(word[2]); } @@ -221,11 +331,303 @@ int mdrun_parse_config(t_mdrun *mdrun) { printf("%s stage: %s\n",ME,word[1]); } else { - printf("%s unknown command %s, skipped!\n",ME,wptr); + printf("%s unknown command %s, skipped!\n",ME,word[0]); continue; } } + /* close file */ + close(fd); + + return 0; +} + +int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) { + + double delta; + + if(!(mdrun->sattr&SATTR_PRELAX)) + return TRUE; + + delta=moldyn->p_avg-moldyn->p_ref; + + if(delta<0) + delta=-delta; + + if(delta>mdrun->dp) + return FALSE; + + return TRUE; +} + +int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) { + + double delta; + + if(!(mdrun->sattr&SATTR_TRELAX)) + return TRUE; + + delta=moldyn->t_avg-moldyn->t_ref; + + if(delta<0) + delta=-delta; + + if(delta>mdrun->dt) + return FALSE; + + return TRUE; +} + +int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { + + insert_atoms_params *iap; + t_stage *stage; + t_atom *atom; + t_3dvec r,v,dist; + double d,dmin; + int cnt,i; + double x,y,z; + double x0,y0,z0; + u8 cr_check,run; + + stage=mdrun->stage->current->data; + iap=stage->params; + + cr_check=FALSE; + + v.x=0.0; v.y=0.0; v.z=0.0; + + switch(iap->type) { + case INS_TOTAL: + x=moldyn->dim.x; + x0=0.0; + y=moldyn->dim.y; + y0=0.0; + z=moldyn->dim.z; + z0=0.0; + cr_check=TRUE; + break; + case INS_REGION: + x=iap->x1-iap->x0; + x0=iap->x0; + y=iap->y1-iap->y0; + y0=iap->y0; + z=iap->z1-iap->z0; + z0=iap->z0; + cr_check=TRUE; + break; + default: + printf("%s unknown insertion mode\n"); + return -1; + } + + cnt=0; + while(cntins_atoms) { + run=1; + while(run) { + r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0; + r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0; + r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0; + run=0; + if(cr_check==TRUE) { + dmin=1000; // for sure too high! + for(i=0;icount;i++) { + atom=&(moldyn->atom[i]); + v3_sub(&dist,&(atom->r),&r); + check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + if(d<(iap->cr*iap->cr)) { + run=1; + break; + } + if(delement,pse_mass[iap->element], + iap->brand,iap->attr,&r,&v); + printf("%s atom inserted: %f %f %f | d squared = %f\n", + ME,r.x,r.y,r.z,dmin); + cnt+=1; + } + + return 0; +} + +int chaatrib(t_moldyn *moldyn,t_mdrun *mdrun) { + + t_stage *stage; + t_chaattr_params *cap; + t_atom *atom; + u8 assigne; + + stage=mdrun->stage->current->data; + cap=stage->params; + + for(i=0;icount;i++) { + atom=&(moldyn->atom[i]); + if(cap->type&CHAATTR_ELEMENT) { + if(cap->element!=atom->element) + continue; + } + if(cap->type&CHAATTR_REGION) { + if(cap->x0r.x) + continue; + if(cap->y0r.y) + continue; + if(cap->z0r.z) + continue; + if(cap->x1>atom->r.x) + continue; + if(cap->y1>atom->r.y) + continue; + if(cap->z1>atom->r.z) + continue; + } + atom->attr=cap->attr; + } + + return 0; +} + +int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) { + + t_stage *stage; + t_chsattr_params *csp; + + stage=mdrun->stage->current->data; + csp=stage->params; + + switch(csp->type) { + case CHSATTR_PCTRL: + set_p_scale(moldyn,csp->ctrl_type,csp->tau); + break; + case CHSATTR_TCTRL: + set_t_scale(moldyn,csp->ctrl_type,csp->tau); + break; + case CHSATTR_PRELAX: + mdrun->sattr&=(^(SATTR_PRELAX)); + if(csp->ctrl==TRUE) + mdrun->sattr|=SATTR_PRELAX; + break; + case CHSATTR_TRELAX: + mdrun->sattr&=(^(SATTR_TRELAX)); + if(csp->ctrl==TRUE) + mdrun->sattr|=SATTR_TRELAX; + break; + case CHSATTR_AVGRST: + mdrun->sattr&=(^(SATTR_AVGRST)); + if(csp->ctrl==TRUE) + mdrun->sattr|=SATTR_AVGRST; + break; + default: + printf("%s unknown system attribute change\n",ME); + return -1; + } + + return 0; +} + +int mdrun_hook(t_moldyn *moldyn,void *ptr) { + + t_mdrun *mdrun; + t_stage *stage; + t_list *sl; + int steps; + double tau; + u8 change_stage; + + t_insert_atoms_params *iap; + t_continue_params *cp; + t_anneal_params *ap; + + mdrun=ptr; + sl=mdrun->stage; + + change_stage=FALSE; + + /* return immediately if there are no more stage */ + if(sl->current==NULL) + return 0; + + /* get stage description */ + stage=mdrun->sl->current->data; + + /* default steps and tau values */ + steps=mdrun->relax_steps; + tau=mdrun->timestep; + + /* check whether relaxation steps are necessary */ + if(!((check_pressure==FALSE)|(check_temperature==FALSE))) { + + /* stage specific stuff */ + switch(stage->type) { + case STAGE_INSERT_ATOMS: + iap=stage->params; + if(iap->cnt_steps==iap->ins_steps) { + change_stage=TRUE; + break; + } + insert_atoms(moldyn,mdrun); + iap->cnt_steps+=1; + break; + case STAGE_CONTINUE: + if(stage->executed==TRUE) { + change_stage=TRUE; + break; + } + cp=stage->params; + steps=cp->runs; + break; + case STAGE_ANNEAL: + ap=stage->params; + if(ap->count==ap->runs) { + change_stage=TRUE; + break; + } + set_temperature(moldyn,moldyn->t_ref+ap->dt); + ap->count+=1; + break; + case STAGE_CHAATTR: + chaatrib(moldyn,mdrun); + change_stage=TRUE; + break; + case STAGE_CHSATTR: + chsatrib(moldyn,mdrun); + change_stage=TRUE; + break; + default: + printf("%s unknwon stage type\n",ME); + break; + } + + /* mark as executed */ + stage->executed=TRUE; + + /* change state */ + if(change_stage==TRUE) { + printf("%s finished stage\n",ME); + if(list_next_f(sl)==L_NO_NEXT_ELEMENT) { + printf("%s no more stages\n",ME); + return 0; + } + steps=0; + tau=mdrun->timestep; + } + + } + else { + + /* averages */ + if(mdrun->sattr&SATTR_AVGRST) + average_reset(moldyn); + + } + + /* continue simulation */ + moldyn_add_schedule(moldyn,steps,tau); + return 0; } @@ -242,9 +644,15 @@ int main(int argc,char **argv) { if(mdrun_parse_argv(&mdrun,argc,argv)<0) return -1; + /* initialize list system */ + list_init_f(&(mdrun.stage)); + /* parse config file */ mdrun_parse_config(&mdrun); + /* reset the stage list */ + list_reset_f(&(mdrun.stage)); + /* sanity check! */ /* prepare simulation */ @@ -326,13 +734,22 @@ int main(int argc,char **argv) { moldyn_set_log(&moldyn,CREATE_REPORT,0); set_avg_skip(&moldyn,mdrun.avgskip); + /* pt scaling */ + if(mdrun.p_tau!=0) + set_p_scale(&moldyn,P_SCALE_BERENDSEN,mdrun.p_tau); + if(mdrun.t_tau!=0) + set_t_scale(&moldyn,T_SCALE_BERENDSEN,mdrun.t_tau); + /* prepare the hook function */ + moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun); /* run the simulation */ moldyn_integrate(&moldyn); /* shutdown */ moldyn_shutdown(&moldyn); + del_stages(&mdrun); + list_destroy_f(&(mdrun.stage)); return 0; } diff --git a/mdrun.h b/mdrun.h index ed9d88e..be8e2a4 100644 --- a/mdrun.h +++ b/mdrun.h @@ -31,25 +31,35 @@ #endif /* - * datatypes + * datatypes & definitions */ typedef struct s_stage { u8 type; - u8 attr; - int runs; + void *params; + u8 executed; } t_stage; +#define STAGE_INSERT_ATOMS 0x01 +#define STAGE_CONTINUE 0x02 +#define STAGE_ANNEAL 0x03 +#define STAGE_CHAATTR 0x04 +#define STAGE_CHSATTR 0x05 + typedef struct s_mdrun { char cfile[128]; // config file + u8 intalgo; // integration algorithm double timestep; // timestep + u8 potential; // potential + double cutoff; // cutoff radius t_3dvec dim; // simulation volume u8 pbcx; // periodic boundary conditions u8 pbcy; u8 pbcz; + int element1; // element 1 double m1; int element2; // element 2 @@ -58,13 +68,19 @@ typedef struct s_mdrun { int lx; // amount of lc units int ly; int lz; - u8 aattrib; // atom attributes u8 lattice; // type of lattice + + u8 sattr; // system attributes double temperature; // temperature double pressure; // pressure double p_tau; // pressure tau double t_tau; // temperature tau + double dp; // delta p fpr pctrl + double dt; // delta t for tctrl + int relax_steps; // amount of relaxation steps + int prerun; // amount of loops in first run + int elog; // logging int tlog; int plog; @@ -74,9 +90,65 @@ typedef struct s_mdrun { u8 vis; int avgskip; // average skip char sdir[128]; // save root - t_list stage; // list of stages + + t_list *stage; // stages + int s_cnt; // stage counter } t_mdrun; +#define SATTR_PRELAX 0x01 +#define SATTR_TRELAX 0x02 +#define SATTR_AVGRST 0x04 + +typedef struct s_insert_atoms_params { + u8 type; + double x0,y0,z0,x1,y1,z1; + double cr; + int ins_steps; + int cnt_steps; + int ins_atoms; + int element; + u8 brand; + u8 aattr; +} t_insert_atoms_params; + +#define INS_TOTAL 0x01 +#define INS_REGION 0x02 + +typedef struct s_continue_params { + int runs; +} t_continue_params; + +typedef struct s_anneal_params { + int runs; + int count; + double dt; +} t_anneal_params; + +typedef struct s_chaattr_params { + u8 type; + double x0,y0,z0; + double x1,y1,z1; + int element; + u8 attr; +} t_chaattr_params; + +#define CHAATTR_TOTALV 0x01 +#define CHAATTR_REGION 0x02 +#define CHAATTR_ELEMENT 0x04 + +typedef struct s_chsattr_params { + u8 type; + double tau; + u8 ctrl; + double delta; +} t_chsattr_params; + +#define CHSATTR_PCTRL 0x01 +#define CHSATTR_TCTRL 0x02 +#define CHSATTR_PRELAX 0x04 +#define CHSATTR_TRELAX 0x08 +#define CHSATTR_AVGRST 0x10 + /* * function prototypes */ diff --git a/moldyn.c b/moldyn.c index b52d51d..988b467 100644 --- a/moldyn.c +++ b/moldyn.c @@ -79,6 +79,50 @@ static char *pse_col[]={ "Ar", }; +static double *pse_mass[]={ + 0, + 0, + 0, + 0, + 0, + 0, + M_C, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + M_SI, + 0, + 0, + 0, + 0, +}; + +static double *pse_lc[]={ + 0, + 0, + 0, + 0, + 0, + 0, + LC_C, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + LC_SI, + 0, + 0, + 0, + 0, +}; + /* * the moldyn functions */ @@ -156,6 +200,48 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { return 0; } +int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) { + + moldyn->pt_scalei&=(^(P_SCALE_MASK)); + moldyn->pt_scale|=ptype; + moldyn->p_tc=ptc; + + printf("[moldyn] p/t scaling:\n"); + + printf(" p: %s",ptype?"yes":"no "); + if(ptype) + printf(" | type: %02x | factor: %f",ptype,ptc); + printf("\n"); + + printf(" t: %s",ttype?"yes":"no "); + if(ttype) + printf(" | type: %02x | factor: %f",ttype,ttc); + printf("\n"); + + return 0; +} + +int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) { + + moldyn->pt_scalei&=(^(T_SCALE_MASK)); + moldyn->pt_scale|=ttype; + moldyn->t_tc=ttc; + + printf("[moldyn] p/t scaling:\n"); + + printf(" p: %s",ptype?"yes":"no "); + if(ptype) + printf(" | type: %02x | factor: %f",ptype,ptc); + printf("\n"); + + printf(" t: %s",ttype?"yes":"no "); + if(ttype) + printf(" | type: %02x | factor: %f",ttype,ttc); + printf("\n"); + + return 0; +} + int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) { moldyn->pt_scale=(ptype|ttype); diff --git a/moldyn.h b/moldyn.h index e2b2ee4..9d49dfa 100644 --- a/moldyn.h +++ b/moldyn.h @@ -49,7 +49,7 @@ typedef struct s_atom { #define ATOM_ATTR_FP 0x01 /* fixed position (bulk material) */ #define ATOM_ATTR_HB 0x02 /* coupled to heat bath (velocity scaling) */ -#define ATOM_ATTR_VA 0x04 /* visualize this atom */ +#define ATOM_ATTR_VA 0x04 /* visualize this atom */ // TODO #define ATOM_ATTR_VB 0x08 /* visualize the bond of this atom */ #define ATOM_ATTR_1BP 0x10 /* single paricle potential */ @@ -255,10 +255,15 @@ typedef struct s_vb { #define MOLDYN_2BP 0x20 /* 2 body */ #define MOLDYN_3BP 0x40 /* and 3 body particle pots */ +#define T_SCALE_NONE 0x00 #define T_SCALE_BERENDSEN 0x01 /* berendsen t control */ #define T_SCALE_DIRECT 0x02 /* direct t control */ +#define T_SCALE_MASK 0x03 + +#define P_SCALE_NONE 0x00 #define P_SCALE_BERENDSEN 0x04 /* berendsen p control */ #define P_SCALE_DIRECT 0x08 /* direct p control */ +#define P_SCALE_MASK 0x0c /* * default values & units @@ -358,6 +363,8 @@ int set_int_alg(t_moldyn *moldyn,u8 algo); int set_cutoff(t_moldyn *moldyn,double cutoff); int set_temperature(t_moldyn *moldyn,double t_ref); int set_pressure(t_moldyn *moldyn,double p_ref); +int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc); +int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc); int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc); int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize); int set_nn_dist(t_moldyn *moldyn,double dist); -- 2.20.1