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"
14 #include "potentials/albe_orig.h"
16 #include "potentials/tersoff_orig.h"
18 #include "potentials/tersoff.h"
27 int mdrun_usage(void) {
29 printf("%s usage:\n",ME);
34 int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) {
41 printf("%s unknown switch: %s\n",ME,argv[i]);
47 strncpy(mdrun->cfile,argv[++i],128);
50 strncpy(mdrun->sdir,argv[++i],128);
56 printf("%s unknown option: %s\n",ME,argv[i]);
66 int del_stages(t_mdrun *mdrun) {
79 stage=sl->current->data;
82 } while(list_next_f(sl)!=L_NO_NEXT_ELEMENT);
87 int add_stage(t_mdrun *mdrun,u8 type,void *params) {
94 case STAGE_DISPLACE_ATOM:
95 psize=sizeof(t_displace_atom_params);
97 case STAGE_INSERT_ATOMS:
98 psize=sizeof(t_insert_atoms_params);
101 psize=sizeof(t_continue_params);
104 psize=sizeof(t_anneal_params);
107 psize=sizeof(t_chaattr_params);
110 psize=sizeof(t_chsattr_params);
113 printf("%s unknown stage type: %02x\n",ME,type);
117 stage=malloc(sizeof(t_stage));
119 perror("[mdrun] malloc (add stage)");
124 stage->executed=FALSE;
126 stage->params=malloc(psize);
127 if(stage->params==NULL) {
128 perror("[mdrun] malloc (stage params)");
132 memcpy(stage->params,params,psize);
134 list_add_immediate_f(&(mdrun->stage),stage);
139 int mdrun_parse_config(t_mdrun *mdrun) {
149 t_displace_atom_params dap;
150 t_insert_atoms_params iap;
151 t_continue_params cp;
153 t_chaattr_params cap;
154 t_chsattr_params csp;
156 /* open config file */
157 fd=open(mdrun->cfile,O_RDONLY);
159 snprintf(error,128,"%s open cfile %s",ME,mdrun->cfile);
164 /* read, parse, set */
169 ret=get_line(fd,line,128);
173 // ignore # lines and \n
174 if((line[0]=='#')|(ret==1))
178 memset(&iap,0,sizeof(t_insert_atoms_params));
179 memset(&cp,0,sizeof(t_continue_params));
180 memset(&ap,0,sizeof(t_anneal_params));
181 memset(&cap,0,sizeof(t_chaattr_params));
182 memset(&csp,0,sizeof(t_chsattr_params));
184 // get command + args
188 wptr=strtok(NULL," \t");
190 wptr=strtok(line," \t");
193 strncpy(word[wcnt],wptr,64);
197 if(!strncmp(word[0],"potential",9)) {
198 if(!strncmp(word[1],"albe",4))
199 mdrun->potential=MOLDYN_POTENTIAL_AM;
200 if(!strncmp(word[1],"albe_o",6))
201 mdrun->potential=MOLDYN_POTENTIAL_AO;
202 if(!strncmp(word[1],"tersoff",7))
203 mdrun->potential=MOLDYN_POTENTIAL_TM;
204 if(!strncmp(word[1],"ho",2))
205 mdrun->potential=MOLDYN_POTENTIAL_HO;
206 if(!strncmp(word[1],"lj",2))
207 mdrun->potential=MOLDYN_POTENTIAL_LJ;
209 else if(!strncmp(word[0],"continue",8))
210 strncpy(mdrun->continue_file,word[1],128);
211 else if(!strncmp(word[0],"cutoff",6))
212 mdrun->cutoff=atof(word[1]);
213 else if(!strncmp(word[0],"nnd",3))
214 mdrun->nnd=atof(word[1]);
215 else if(!strncmp(word[0],"intalgo",7)) {
216 if(!strncmp(word[1],"verlet",5))
217 mdrun->intalgo=MOLDYN_INTEGRATE_VERLET;
219 else if(!strncmp(word[0],"timestep",8))
220 mdrun->timestep=atof(word[1]);
221 else if(!strncmp(word[0],"volume",6)) {
222 mdrun->dim.x=atof(word[1]);
223 mdrun->dim.y=atof(word[2]);
224 mdrun->dim.z=atof(word[3]);
225 if(strncmp(word[4],"0",1))
228 else if(!strncmp(word[0],"pbc",3)) {
229 if(strncmp(word[1],"0",1))
233 if(strncmp(word[2],"0",1))
237 if(strncmp(word[3],"0",1))
242 else if(!strncmp(word[0],"temperature",11))
243 mdrun->temperature=atof(word[1]);
244 else if(!strncmp(word[0],"pressure",8))
245 mdrun->pressure=atof(word[1]);
246 else if(!strncmp(word[0],"lattice",7)) {
247 if(!strncmp(word[1],"zincblende",10))
248 mdrun->lattice=ZINCBLENDE;
249 if(!strncmp(word[1],"fcc",3))
251 if(!strncmp(word[1],"diamond",7))
252 mdrun->lattice=DIAMOND;
254 else if(!strncmp(word[0],"element1",8)) {
255 mdrun->element1=atoi(word[1]);
256 mdrun->m1=pse_mass[mdrun->element1];
258 else if(!strncmp(word[0],"element2",8)) {
259 mdrun->element2=atoi(word[1]);
260 mdrun->m2=pse_mass[mdrun->element2];
262 else if(!strncmp(word[0],"fill",6)) {
263 // only lc mode by now
264 mdrun->lx=atoi(word[2]);
265 mdrun->ly=atoi(word[3]);
266 mdrun->lz=atoi(word[4]);
267 mdrun->lc=atof(word[5]);
269 else if(!strncmp(word[0],"aattr",5)) {
270 // for aatrib line we need a special stage
271 // containing one schedule of 0 loops ...
272 for(i=0;i<strlen(word[1]);i++) {
275 cap.type|=CHAATTR_TOTALV;
278 cap.type|=CHAATTR_REGION;
281 cap.type|=CHAATTR_ELEMENT;
288 if(cap.type&CHAATTR_REGION) {
289 cap.x0=atof(word[1]);
290 cap.y0=atof(word[2]);
291 cap.z0=atof(word[3]);
292 cap.x1=atof(word[4]);
293 cap.y1=atof(word[5]);
294 cap.z1=atof(word[6]);
297 if(cap.type&CHAATTR_ELEMENT) {
298 cap.element=atoi(word[i]);
301 for(o=0;o<strlen(word[i]);o++) {
304 cap.attr|=ATOM_ATTR_VB;
307 cap.attr|=ATOM_ATTR_HB;
310 cap.attr|=ATOM_ATTR_VA;
313 cap.attr|=ATOM_ATTR_FP;
316 cap.attr|=ATOM_ATTR_1BP;
319 cap.attr|=ATOM_ATTR_2BP;
322 cap.attr|=ATOM_ATTR_3BP;
328 add_stage(mdrun,STAGE_CHAATTR,&cap);
330 else if(!strncmp(word[0],"sattr",5)) {
331 // for satrib line we need a special stage
332 // containing one schedule of 0 loops ...
334 for(i=1;i<wcnt;i++) {
335 if(!strncmp(word[i],"pctrl",5)) {
336 csp.ptau=0.01/(atof(word[++i])*GPA);
337 csp.type|=CHSATTR_PCTRL;
339 if(!strncmp(word[i],"tctrl",5)) {
340 csp.ttau=atof(word[++i]);
341 csp.type|=CHSATTR_TCTRL;
343 if(!strncmp(word[i],"prelax",6)) {
344 csp.dp=atof(word[++i])*BAR;
345 csp.type|=CHSATTR_PRELAX;
347 if(!strncmp(word[i],"trelax",6)) {
348 csp.dt=atof(word[++i]);
349 csp.type|=CHSATTR_TRELAX;
351 if(!strncmp(word[i],"rsteps",6)) {
352 csp.rsteps=atoi(word[++i]);
353 csp.type|=CHSATTR_RSTEPS;
355 if(!strncmp(word[i],"avgrst",6)) {
356 csp.avgrst=atoi(word[++i]);
357 csp.type|=CHSATTR_AVGRST;
360 add_stage(mdrun,STAGE_CHSATTR,&csp);
362 else if(!strncmp(word[0],"prerun",6))
363 mdrun->prerun=atoi(word[1]);
364 else if(!strncmp(word[0],"avgskip",7))
365 mdrun->avgskip=atoi(word[1]);
366 else if(!strncmp(word[0],"elog",4))
367 mdrun->elog=atoi(word[1]);
368 else if(!strncmp(word[0],"tlog",4))
369 mdrun->tlog=atoi(word[1]);
370 else if(!strncmp(word[0],"plog",4))
371 mdrun->plog=atoi(word[1]);
372 else if(!strncmp(word[0],"vlog",4))
373 mdrun->vlog=atoi(word[1]);
374 else if(!strncmp(word[0],"save",4))
375 mdrun->save=atoi(word[1]);
376 else if(!strncmp(word[0],"visualize",9))
377 mdrun->visualize=atoi(word[1]);
378 else if(!strncmp(word[0],"stage",5)) {
379 // for every stage line, add a stage
380 if(!strncmp(word[1],"displace",8)) {
381 dap.nr=atoi(word[2]);
382 dap.dx=atof(word[3]);
383 dap.dy=atof(word[4]);
384 dap.dz=atof(word[5]);
385 add_stage(mdrun,STAGE_DISPLACE_ATOM,&dap);
387 else if(!strncmp(word[1],"ins_atoms",9)) {
388 iap.ins_steps=atoi(word[2]);
389 iap.ins_atoms=atoi(word[3]);
390 iap.element=atoi(word[4]);
391 iap.brand=atoi(word[5]);
392 for(i=0;i<strlen(word[6]);i++) {
395 iap.attr|=ATOM_ATTR_VB;
398 iap.attr|=ATOM_ATTR_HB;
401 iap.attr|=ATOM_ATTR_VA;
404 iap.attr|=ATOM_ATTR_FP;
407 iap.attr|=ATOM_ATTR_1BP;
410 iap.attr|=ATOM_ATTR_2BP;
413 iap.attr|=ATOM_ATTR_3BP;
423 iap.x0=atof(word[8]);
424 iap.y0=atof(word[9]);
425 iap.z0=atof(word[10]);
428 if(word[8][0]=='t') {
430 iap.cr=atof(word[9]);
434 iap.x0=atof(word[8]);
435 iap.y0=atof(word[9]);
436 iap.z0=atof(word[10]);
437 iap.x1=atof(word[11]);
438 iap.y1=atof(word[12]);
439 iap.z1=atof(word[13]);
440 iap.cr=atof(word[14]);
445 add_stage(mdrun,STAGE_INSERT_ATOMS,&iap);
447 else if(!strncmp(word[1],"continue",8)) {
448 cp.runs=atoi(word[2]);
449 add_stage(mdrun,STAGE_CONTINUE,&cp);
451 else if(!strncmp(word[1],"anneal",6)) {
453 ap.runs=atoi(word[2]);
455 add_stage(mdrun,STAGE_ANNEAL,&ap);
458 printf("%s unknown stage type: %s\n",
464 printf("%s unknown keyword '%s', skipped!\n",
476 int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
480 if(!(mdrun->sattr&SATTR_PRELAX)) {
484 delta=moldyn->p_avg-moldyn->p_ref;
495 int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
499 if(!(mdrun->sattr&SATTR_TRELAX))
502 delta=moldyn->t_avg-moldyn->t_ref;
513 int displace_atom(t_moldyn *moldyn,t_mdrun *mdrun) {
515 t_displace_atom_params *dap;
519 stage=mdrun->stage.current->data;
522 atom=&(moldyn->atom[dap->nr]);
530 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
532 t_insert_atoms_params *iap;
542 stage=mdrun->stage.current->data;
547 v.x=0.0; v.y=0.0; v.z=0.0;
551 switch(mdrun->lattice) {
592 printf("%s unknown insertion mode: %02x\n",
598 while(cnt<iap->ins_atoms) {
601 if(iap->type!=INS_POS) {
602 r.x=rand_get_double(&(moldyn->random))*x;
603 r.y=rand_get_double(&(moldyn->random))*y;
604 r.z=rand_get_double(&(moldyn->random))*z;
615 if(iap->type!=INS_TOTAL) {
622 dmin=1000; // for sure too high!
623 for(i=0;i<moldyn->count;i++) {
624 atom=&(moldyn->atom[i]);
625 v3_sub(&dist,&(atom->r),&r);
626 check_per_bound(moldyn,&dist);
627 d=v3_absolute_square(&dist);
628 if(d<(iap->cr*iap->cr)) {
637 add_atom(moldyn,iap->element,pse_mass[iap->element],
638 iap->brand,iap->attr,&r,&v);
639 printf("%s atom inserted (%d/%d): %f %f %f\n",
640 ME,(iap->cnt_steps+1)*iap->ins_atoms,
641 iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z);
642 printf(" -> d2 = %f/%f\n",dmin,iap->cr*iap->cr);
649 int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) {
652 t_chaattr_params *cap;
656 stage=mdrun->stage.current->data;
659 for(i=0;i<moldyn->count;i++) {
660 atom=&(moldyn->atom[i]);
661 if(cap->type&CHAATTR_ELEMENT) {
662 if(cap->element!=atom->element)
665 if(cap->type&CHAATTR_REGION) {
666 if(cap->x0<atom->r.x)
668 if(cap->y0<atom->r.y)
670 if(cap->z0<atom->r.z)
672 if(cap->x1>atom->r.x)
674 if(cap->y1>atom->r.y)
676 if(cap->z1>atom->r.z)
679 atom->attr=cap->attr;
685 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
688 t_chsattr_params *csp;
690 stage=mdrun->stage.current->data;
693 if(csp->type&CHSATTR_PCTRL) {
695 set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
697 set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
699 if(csp->type&CHSATTR_TCTRL) {
701 set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
703 set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
705 if(csp->type&CHSATTR_PRELAX) {
707 mdrun->sattr&=(~(SATTR_PRELAX));
709 mdrun->sattr|=SATTR_PRELAX;
712 if(csp->type&CHSATTR_TRELAX) {
714 mdrun->sattr&=(~(SATTR_TRELAX));
716 mdrun->sattr|=SATTR_TRELAX;
719 if(csp->type&CHSATTR_AVGRST) {
721 mdrun->sattr|=SATTR_AVGRST;
723 mdrun->sattr&=(~(SATTR_AVGRST));
725 if(csp->type&CHSATTR_RSTEPS) {
726 mdrun->relax_steps=csp->rsteps;
732 #define stage_print(m) if(!(stage->executed)) \
735 int mdrun_hook(void *ptr1,void *ptr2) {
745 t_insert_atoms_params *iap;
746 t_continue_params *cp;
756 /* return immediately if there are no more stage */
757 if(sl->current==NULL)
760 /* get stage description */
761 stage=sl->current->data;
763 /* default steps and tau values */
764 steps=mdrun->relax_steps;
767 /* check whether relaxation steps are necessary */
768 if((check_pressure(moldyn,mdrun)==TRUE)&\
769 (check_temperature(moldyn,mdrun)==TRUE)) {
772 stage_print("\n###########################\n");
773 stage_print("# [mdrun] executing stage #\n");
774 stage_print("###########################\n\n");
776 /* stage specific stuff */
777 switch(stage->type) {
778 case STAGE_DISPLACE_ATOM:
779 stage_print(" -> displace atom\n\n");
780 displace_atom(moldyn,mdrun);
783 case STAGE_INSERT_ATOMS:
784 stage_print(" -> insert atoms\n\n");
786 if(iap->cnt_steps==iap->ins_steps) {
790 insert_atoms(moldyn,mdrun);
794 stage_print(" -> continue\n\n");
795 if(stage->executed==TRUE) {
803 stage_print(" -> anneal\n\n");
805 if(ap->count==ap->runs) {
809 if(moldyn->t_ref+ap->dt>=0.0)
810 set_temperature(moldyn,
811 moldyn->t_ref+ap->dt);
815 stage_print(" -> chaattr\n\n");
816 chaatr(moldyn,mdrun);
820 stage_print(" -> chsattr\n\n");
821 chsattr(moldyn,mdrun);
825 printf("%s unknwon stage type\n",ME);
829 /* mark as executed */
830 stage->executed=TRUE;
833 if(change_stage==TRUE) {
834 printf("%s finished stage\n",ME);
835 if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
836 printf("%s no more stages\n",ME);
847 if(mdrun->sattr&SATTR_AVGRST)
848 average_reset(moldyn);
852 /* continue simulation */
853 moldyn_add_schedule(moldyn,steps,tau);
858 int main(int argc,char **argv) {
865 memset(&mdrun,0,sizeof(t_mdrun));
866 memset(&moldyn,0,sizeof(t_moldyn));
868 /* parse arguments */
869 if(mdrun_parse_argv(&mdrun,argc,argv)<0)
872 /* initialize list system */
873 list_init_f(&(mdrun.stage));
875 /* parse config file */
876 mdrun_parse_config(&mdrun);
878 /* reset the stage list */
879 list_reset_f(&(mdrun.stage));
883 /* prepare simulation */
885 if(mdrun.continue_file[0]) {
886 // read the save file
887 moldyn_read_save_file(&moldyn,mdrun.continue_file);
888 // manualaadjusting some stuff
891 rand_init(&(moldyn.random),NULL,1);
892 moldyn.random.status|=RAND_STAT_VERBOSE;
895 moldyn_init(&moldyn,argc,argv);
898 if(set_int_alg(&moldyn,mdrun.intalgo)<0)
902 set_cutoff(&moldyn,mdrun.cutoff);
903 if(set_potential(&moldyn,mdrun.potential)<0)
905 switch(mdrun.potential) {
906 case MOLDYN_POTENTIAL_AM:
907 albe_mult_set_params(&moldyn,
911 case MOLDYN_POTENTIAL_AO:
912 albe_orig_mult_set_params(&moldyn,
916 case MOLDYN_POTENTIAL_TM:
917 tersoff_mult_set_params(&moldyn,
921 case MOLDYN_POTENTIAL_HO:
922 harmonic_oscillator_set_params(&moldyn,mdrun.element1);
924 case MOLDYN_POTENTIAL_LJ:
925 lennard_jones_set_params(&moldyn,mdrun.element1);
928 printf("%s unknown potential: %02x\n",
933 /* if it is a continue run, reset schedule and skip lattice init */
934 if(mdrun.continue_file[0]) {
935 memset(&(moldyn.schedule),0,sizeof(t_moldyn_schedule));
939 /* initial lattice and dimensions */
940 set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
941 set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
942 switch(mdrun.lattice) {
944 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
945 mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
946 mdrun.ly,mdrun.lz,NULL);
949 create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1,
950 mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
951 mdrun.ly,mdrun.lz,NULL);
954 o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
955 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
956 mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
957 mdrun.ly,mdrun.lz,&o);
958 o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
959 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
960 mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx,
961 mdrun.ly,mdrun.lz,&o);
964 printf("%s unknown lattice type: %02x\n",
968 moldyn_bc_check(&moldyn);
970 /* temperature and pressure */
971 set_temperature(&moldyn,mdrun.temperature);
972 set_pressure(&moldyn,mdrun.pressure);
973 thermal_init(&moldyn,TRUE);
977 moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
980 moldyn_set_log_dir(&moldyn,mdrun.sdir);
981 moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
983 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
985 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
987 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
989 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
991 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
993 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
994 moldyn_set_log(&moldyn,CREATE_REPORT,0);
995 set_avg_skip(&moldyn,mdrun.avgskip);
997 /* prepare the hook function */
998 moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
1000 /* run the simulation */
1001 moldyn_integrate(&moldyn);
1004 moldyn_shutdown(&moldyn);
1006 list_destroy_f(&(mdrun.stage));