2 * mdrun.c - main code to run a md simulation
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
11 #include "potentials/harmonic_oscillator.h"
12 #include "potentials/lennard_jones.h"
13 #include "potentials/albe.h"
15 #include "potentials/tersoff_orig.h"
17 #include "potentials/tersoff.h"
26 int mdrun_usage(void) {
28 printf("%s usage:\n",ME);
33 int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) {
40 printf("%s unknown switch: %s\n",ME,argv[i]);
46 strncpy(mdrun->cfile,argv[++i],128);
49 strncpy(mdrun->sdir,argv[++i],128);
55 printf("%s unknown option: %s\n",ME,argv[i]);
65 del_stages(t_mdrun *mdrun) {
78 stage=sl->current->data;
81 } while(list_next_f(sl)!=L_NO_NEXT_ELEMENT);
86 add_stage(t_mdrun *mdrun,u8 type,void *params) {
93 case STAGE_INSERT_ATOMS:
94 psize=sizeof(t_insert_atoms_params);
97 psize=sizeof(t_continue_params);
100 psize=sizeof(t_anneal_params);
103 psize=sizeof(t_chaattr_params);
106 psize=sizeof(t_chsattr_params);
112 stage=malloc(sizeof(t_stage));
114 perror("[mdrun] malloc (add stage)");
119 stage->executed=FALSE;
121 stage->params=malloc(psize);
122 if(stage->params==NULL) {
123 perror("[mdrun] malloc (stage params)");
127 memcpy(stage->params,params,psize);
129 list_add_immediate_f(mdrun->stage,stage);
134 int mdrun_parse_config(t_mdrun *mdrun) {
144 t_insert_atoms_params iap;
145 t_continue_params cp;
147 t_chaattr_params cap;
148 t_chsattr_params csp;
150 /* open config file */
151 fd=open(mdrun->cfile,O_RDONLY);
153 snprintf(error,128,"%s open cfile %s",ME,mdrun->cfile);
158 /* read, parse, set */
163 ret=get_line(fd,line,128);
167 // ignore # lines and \n
168 if((line[0]=='#')|(ret==1))
171 // get command + args
175 wptr=strtok(NULL," \t");
177 wptr=strtok(line," \t");
180 strncpy(word[wcnt],wptr,32);
184 if(!strncmp(word[0],"potential",9)) {
185 if(!strncmp(word[1],"albe",4))
186 mdrun->potential=MOLDYN_POTENTIAL_AM;
187 if(!strncmp(word[1],"tersoff",7))
188 mdrun->potential=MOLDYN_POTENTIAL_TM;
189 if(!strncmp(word[1],"ho",2))
190 mdrun->potential=MOLDYN_POTENTIAL_HO;
191 if(!strncmp(word[1],"lj",2))
192 mdrun->potential=MOLDYN_POTENTIAL_LJ;
194 else if(!strncmp(word[0],"cutoff",6))
195 mdrun->cutoff=atof(word[1]);
196 else if(!strncmp(word[0],"intalgo",7)) {
197 if(!strncmp(word[1],"verlet",5))
198 mdrun->intalgo=MOLDYN_INTEGRATE_VERLET;
200 else if(!strncmp(word[0],"timestep",8))
201 mdrun->timestep=atof(word[1]);
202 else if(!strncmp(word[0],"volume",6)) {
203 mdrun->dim.x=atof(word[1]);
204 mdrun->dim.y=atof(word[2]);
205 mdrun->dim.z=atof(word[3]);
206 if(strncmp(word[4],"0",1))
209 else if(!strncmp(word[0],"pbc",3)) {
210 if(strncmp(word[1],"0",1))
214 if(strncmp(word[2],"0",1))
218 if(strncmp(word[3],"0",1))
223 else if(!strncmp(word[0],"temperature",11))
224 mdrun->temperature=atof(word[1]);
225 else if(!strncmp(word[0],"pressure",8))
226 mdrun->pressure=atof(word[1]);
227 else if(!strncmp(word[0],"lattice",7)) {
228 if(!strncmp(word[1],"zincblende",10))
229 mdrun->lattice=ZINCBLENDE;
230 if(!strncmp(word[1],"fcc",3))
232 if(!strncmp(word[1],"diamond",7))
233 mdrun->lattice=DIAMOND;
235 else if(!strncmp(word[0],"element1",8)) {
236 mdrun->element1=atoi(word[1]);
237 switch(mdrun->element1) {
245 printf("%s unknown element1: %s|%d\n",
246 ME,word[1],mdrun->element1);
250 else if(!strncmp(word[0],"element2",8)) {
251 mdrun->element2=atoi(word[1]);
252 switch(mdrun->element2) {
260 printf("%s unknown element2: %s|%d\n",
261 ME,word[1],mdrun->element2);
265 else if(!strncmp(word[0],"fill",6)) {
266 // only lc mode by now
267 mdrun->lx=atoi(word[2]);
268 mdrun->ly=atoi(word[3]);
269 mdrun->lz=atoi(word[4]);
271 else if(!strncmp(word[0],"aattr",5)) {
272 // for aatrib line we need a special stage
273 // containing one schedule of 0 loops ...
274 if(!strncmp(word[1],"all",3)) {
280 for(i=0;i<strlen(wptr)) {
283 cap.attr|=ATOM_ATTR_VB;
286 cap.attr|=ATOM_ATTR_HB;
289 cap.attr|=ATOM_ATTR_VA;
292 cap.attr|=ATOM_ATTR_FP;
295 cap.attr|=ATOM_ATTR_1BP;
298 cap.attr|=ATOM_ATTR_2BP;
301 cap.attr|=ATOM_ATTR_3BP;
307 add_stage(mdrun,STAGE_CHAATTR,&cap);
309 else if(!strncmp(word[0],"ectrl",5)) {
310 mdrun->p_tau=atof(word[1]);
311 mdrun->t_tau=atof(word[2]);
313 else if(!strncmp(word[0],"prerun",6))
314 mdrun->prerun=atoi(word[1]);
315 else if(!strncmp(word[0],"avgskip",7))
316 mdrun->avgskip=atoi(word[1]);
317 else if(!strncmp(word[0],"elog",4))
318 mdrun->elog=atoi(word[1]);
319 else if(!strncmp(word[0],"tlog",4))
320 mdrun->tlog=atoi(word[1]);
321 else if(!strncmp(word[0],"plog",4))
322 mdrun->plog=atoi(word[1]);
323 else if(!strncmp(word[0],"vlog",4))
324 mdrun->vlog=atoi(word[1]);
325 else if(!strncmp(word[0],"save",4))
326 mdrun->save=atoi(word[1]);
327 else if(!strncmp(word[0],"visualize",9))
328 mdrun->visualize=atoi(word[1]);
329 else if(!strncmp(word[0],"stage",5)) {
330 // for every stage line, add a stage
331 printf("%s stage: %s\n",ME,word[1]);
334 printf("%s unknown command %s, skipped!\n",ME,word[0]);
345 int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
349 if(!(mdrun->sattr&SATTR_PRELAX))
352 delta=moldyn->p_avg-moldyn->p_ref;
363 int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
367 if(!(mdrun->sattr&SATTR_TRELAX))
370 delta=moldyn->t_avg-moldyn->t_ref;
381 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
383 insert_atoms_params *iap;
393 stage=mdrun->stage->current->data;
398 v.x=0.0; v.y=0.0; v.z=0.0;
420 printf("%s unknown insertion mode\n");
425 while(cnt<iap->ins_atoms) {
428 r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0;
429 r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0;
430 r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0;
433 dmin=1000; // for sure too high!
434 for(i=0;i<moldyn->count;i++) {
435 atom=&(moldyn->atom[i]);
436 v3_sub(&dist,&(atom->r),&r);
437 check_per_bound(moldyn,&dist);
438 d=v3_absolute_square(&dist);
439 if(d<(iap->cr*iap->cr)) {
448 add_atom(moldyn,iap->element,pse_mass[iap->element],
449 iap->brand,iap->attr,&r,&v);
450 printf("%s atom inserted: %f %f %f | d squared = %f\n",
451 ME,r.x,r.y,r.z,dmin);
458 int chaatrib(t_moldyn *moldyn,t_mdrun *mdrun) {
461 t_chaattr_params *cap;
465 stage=mdrun->stage->current->data;
468 for(i=0;i<moldyn->count;i++) {
469 atom=&(moldyn->atom[i]);
470 if(cap->type&CHAATTR_ELEMENT) {
471 if(cap->element!=atom->element)
474 if(cap->type&CHAATTR_REGION) {
475 if(cap->x0<atom->r.x)
477 if(cap->y0<atom->r.y)
479 if(cap->z0<atom->r.z)
481 if(cap->x1>atom->r.x)
483 if(cap->y1>atom->r.y)
485 if(cap->z1>atom->r.z)
488 atom->attr=cap->attr;
494 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
497 t_chsattr_params *csp;
499 stage=mdrun->stage->current->data;
504 set_p_scale(moldyn,csp->ctrl_type,csp->tau);
507 set_t_scale(moldyn,csp->ctrl_type,csp->tau);
510 mdrun->sattr&=(^(SATTR_PRELAX));
512 mdrun->sattr|=SATTR_PRELAX;
515 mdrun->sattr&=(^(SATTR_TRELAX));
517 mdrun->sattr|=SATTR_TRELAX;
520 mdrun->sattr&=(^(SATTR_AVGRST));
522 mdrun->sattr|=SATTR_AVGRST;
525 printf("%s unknown system attribute change\n",ME);
532 int mdrun_hook(t_moldyn *moldyn,void *ptr) {
541 t_insert_atoms_params *iap;
542 t_continue_params *cp;
550 /* return immediately if there are no more stage */
551 if(sl->current==NULL)
554 /* get stage description */
555 stage=mdrun->sl->current->data;
557 /* default steps and tau values */
558 steps=mdrun->relax_steps;
561 /* check whether relaxation steps are necessary */
562 if(!((check_pressure==FALSE)|(check_temperature==FALSE))) {
564 /* stage specific stuff */
565 switch(stage->type) {
566 case STAGE_INSERT_ATOMS:
568 if(iap->cnt_steps==iap->ins_steps) {
572 insert_atoms(moldyn,mdrun);
576 if(stage->executed==TRUE) {
585 if(ap->count==ap->runs) {
589 set_temperature(moldyn,moldyn->t_ref+ap->dt);
593 chaatrib(moldyn,mdrun);
597 chsatrib(moldyn,mdrun);
601 printf("%s unknwon stage type\n",ME);
605 /* mark as executed */
606 stage->executed=TRUE;
609 if(change_stage==TRUE) {
610 printf("%s finished stage\n",ME);
611 if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
612 printf("%s no more stages\n",ME);
623 if(mdrun->sattr&SATTR_AVGRST)
624 average_reset(moldyn);
628 /* continue simulation */
629 moldyn_add_schedule(moldyn,steps,tau);
634 int main(int argc,char **argv) {
640 /* clear mdrun struct */
641 memset(&mdrun,0,sizeof(t_mdrun));
643 /* parse arguments */
644 if(mdrun_parse_argv(&mdrun,argc,argv)<0)
647 /* initialize list system */
648 list_init_f(&(mdrun.stage));
650 /* parse config file */
651 mdrun_parse_config(&mdrun);
653 /* reset the stage list */
654 list_reset_f(&(mdrun.stage));
658 /* prepare simulation */
659 moldyn_init(&moldyn,argc,argv);
660 if(set_int_alg(&moldyn,mdrun.intalgo)<0)
662 set_cutoff(&moldyn,mdrun.cutoff);
663 if(set_potential(&moldyn,mdrun.potential)<0)
665 switch(mdrun.potential) {
666 case MOLDYN_POTENTIAL_AM:
667 albe_mult_set_params(&moldyn,
671 case MOLDYN_POTENTIAL_TM:
672 tersoff_mult_set_params(&moldyn,
676 case MOLDYN_POTENTIAL_HO:
677 harmonic_oscillator_set_params(&moldyn,mdrun.element1);
679 case MOLDYN_POTENTIAL_LJ:
680 lennard_jones_set_params(&moldyn,mdrun.element1);
683 printf("%s unknown potential: %02x\n",
687 set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
688 set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
689 switch(mdrun.lattice) {
691 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
692 mdrun.m1,0,0,mdrun.lx,
693 mdrun.ly,mdrun.lz,NULL);
696 create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1,
697 mdrun.m1,0,0,mdrun.lx,
698 mdrun.ly,mdrun.lz,NULL);
701 o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
702 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
703 mdrun.m1,0,0,mdrun.lx,
704 mdrun.ly,mdrun.lz,&o);
705 o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
706 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
707 mdrun.m2,0,0,mdrun.lx,
708 mdrun.ly,mdrun.lz,&o);
711 printf("%s unknown lattice type: %02x\n",
715 moldyn_bc_check(&moldyn);
716 set_temperature(&moldyn,mdrun.temperature);
717 set_pressure(&moldyn,mdrun.pressure);
718 thermal_init(&moldyn,TRUE);
719 moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
720 moldyn_set_log_dir(&moldyn,mdrun.sdir);
721 moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
723 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
725 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
727 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
729 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
731 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
733 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
734 moldyn_set_log(&moldyn,CREATE_REPORT,0);
735 set_avg_skip(&moldyn,mdrun.avgskip);
739 set_p_scale(&moldyn,P_SCALE_BERENDSEN,mdrun.p_tau);
741 set_t_scale(&moldyn,T_SCALE_BERENDSEN,mdrun.t_tau);
743 /* prepare the hook function */
744 moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
746 /* run the simulation */
747 moldyn_integrate(&moldyn);
750 moldyn_shutdown(&moldyn);
752 list_destroy_f(&(mdrun.stage));