]> hackdaworld.org Git - physik/posic.git/commitdiff
safety checkin
authorhackbard <hackbard@sage.physik.uni-augsburg.de>
Mon, 28 Apr 2008 17:17:06 +0000 (19:17 +0200)
committerhackbard <hackbard@sage.physik.uni-augsburg.de>
Mon, 28 Apr 2008 17:17:06 +0000 (19:17 +0200)
config.default
mdrun.c
mdrun.h
moldyn.c
moldyn.h

index f221555303db5d97fd4260caf092f5da8a2b9902..0bd26d244cf09ac26ac8ce20b52ff7d8a8cdc68c 100644 (file)
@@ -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 4a2822e86d10c5d967bcc55ce00de6964c26ba8b..6f7a2f02a0704491b5d8bf50762031e0149711ba 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;
 }
diff --git a/mdrun.h b/mdrun.h
index ed9d88e0ac410057ec87ad8b28a3774a8d6b4dbc..be8e2a40df51db55552df1366182de9d6eae38b1 100644 (file)
--- a/mdrun.h
+++ b/mdrun.h
 #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
  */
index b52d51d4b790728b9e14f1fd38f48ebd41c8a7b2..988b4671780d6fb68a8826c92390c2cedd2ff663 100644 (file)
--- 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);
index e2b2ee40d0cbac9286c61b96a50a355afe637b1e..9d49dfa809468f8d8b77fbe59fb08cc6bbc2ee2b 100644 (file)
--- 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);