safety checkin
[physik/posic.git] / mdrun.c
diff --git a/mdrun.c b/mdrun.c
index 4a2822e..6f7a2f0 100644 (file)
--- 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;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]);
                }
@@ -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(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;
 }
 
@@ -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;
 }