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;
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);
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;i<strlen(wptr)) {
+ switch(word[2][i]) {
+ case 'b':
+ cap.attr|=ATOM_ATTR_VB;
+ break;
+ case 'h':
+ cap.attr|=ATOM_ATTR_HB;
+ break;
+ case 'v':
+ cap.attr|=ATOM_ATTR_VA;
+ break;
+ case 'f':
+ cap.attr|=ATOM_ATTR_FP;
+ break;
+ case '1':
+ cap.attr|=ATOM_ATTR_1BP;
+ break;
+ case '2':
+ cap.attr|=ATOM_ATTR_2BP;
+ break;
+ case '3':
+ cap.attr|=ATOM_ATTR_3BP;
+ break;
+ default:
+ break;
+ }
+ }
+ add_stage(mdrun,STAGE_CHAATTR,&cap);
}
- else if(!strncmp(word[0],"eattrib",7)) {
+ else if(!strncmp(word[0],"ectrl",5)) {
mdrun->p_tau=atof(word[1]);
mdrun->t_tau=atof(word[2]);
}
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(cnt<iap->ins_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;i<moldyn->count;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(d<dmin)
+ dmin=d;
+ }
+ }
+ }
+ add_atom(moldyn,iap->element,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;i<moldyn->count;i++) {
+ atom=&(moldyn->atom[i]);
+ if(cap->type&CHAATTR_ELEMENT) {
+ if(cap->element!=atom->element)
+ continue;
+ }
+ if(cap->type&CHAATTR_REGION) {
+ if(cap->x0<atom->r.x)
+ continue;
+ if(cap->y0<atom->r.y)
+ continue;
+ if(cap->z0<atom->r.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;
}
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 */
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;
}
#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
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;
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
*/
#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 */
#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
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);
projects / physik / posic.git / commitdiff