]> www.hackdaworld.org Git - physik/posic.git/blob - mdrun.c
31abd1c069a0cbe6dd5660cbc82fae22990ef672
[physik/posic.git] / mdrun.c
1 /*
2  * mdrun.c - main code to run a md simulation
3  *
4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
5  *
6  */
7
8 #include "mdrun.h"
9
10 /* potential */
11 #include "potentials/harmonic_oscillator.h"
12 #include "potentials/lennard_jones.h"
13 #include "potentials/albe.h"
14 #ifdef TERSOFF_ORIG
15 #include "potentials/tersoff_orig.h"
16 #else
17 #include "potentials/tersoff.h"
18 #endif
19
20 #define ME      "[mdrun]"
21
22 /*
23  * parse function
24  */
25
26 int mdrun_usage(void) {
27
28         printf("%s usage:\n",ME);
29
30         return 0;
31 }
32
33 int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) {
34
35         int i;
36
37         for(i=1;i<argc;i++) {
38
39                 if(argv[i][0]!='-') {
40                         printf("%s unknown switch: %s\n",ME,argv[i]);
41                         return -1;
42                 }
43
44                 switch(argv[i][1]) {
45                         case 'c':
46                                 strncpy(mdrun->cfile,argv[++i],128);
47                                 break;
48                         case 's':
49                                 strncpy(mdrun->sdir,argv[++i],128);
50                                 break;
51                         case 'h':
52                                 mdrun_usage();
53                                 break;
54                         default:
55                                 printf("%s unknown option: %s\n",ME,argv[i]);
56                                 mdrun_usage();
57                                 return -1;
58                 }
59
60         }
61
62         return 0;
63 }
64
65 int del_stages(t_mdrun *mdrun) {
66
67         t_list *sl;
68         t_stage *stage;
69
70         sl=&(mdrun->stage);
71
72         list_reset_f(sl);
73
74         if(sl->start==NULL)
75                 return 0;
76
77         do {
78                 stage=sl->current->data;
79                 free(stage->params);
80                 free(stage);
81         } while(list_next_f(sl)!=L_NO_NEXT_ELEMENT);
82
83         return 0;
84 }
85
86 int add_stage(t_mdrun *mdrun,u8 type,void *params) {
87
88         int psize;
89
90         t_stage *stage;
91
92         switch(type) {
93                 case STAGE_INSERT_ATOMS:
94                         psize=sizeof(t_insert_atoms_params);
95                         break;
96                 case STAGE_CONTINUE:
97                         psize=sizeof(t_continue_params);
98                         break;
99                 case STAGE_ANNEAL:
100                         psize=sizeof(t_anneal_params);
101                         break;
102                 case STAGE_CHAATTR:
103                         psize=sizeof(t_chaattr_params);
104                         break;
105                 case STAGE_CHSATTR:
106                         psize=sizeof(t_chsattr_params);
107                         break;
108                 default:
109                         printf("%s unknown stage type: %02x\n",ME,type);
110                         return -1;
111         }
112
113         stage=malloc(sizeof(t_stage));
114         if(stage==NULL) {
115                 perror("[mdrun] malloc (add stage)");
116                 return -1;
117         }
118
119         stage->type=type;
120         stage->executed=FALSE;
121
122         stage->params=malloc(psize);
123         if(stage->params==NULL) {
124                 perror("[mdrun] malloc (stage params)");
125                 return -1;
126         }
127
128         memcpy(stage->params,params,psize);
129
130         list_add_immediate_f(&(mdrun->stage),stage);
131
132         return 0;
133 }
134
135 int mdrun_parse_config(t_mdrun *mdrun) {
136
137         int fd,ret;
138         char error[128];
139         char line[128];
140         char *wptr;
141         char word[16][32];
142         int wcnt;
143         int i,o;
144
145         t_insert_atoms_params iap;
146         t_continue_params cp;
147         t_anneal_params ap;
148         t_chaattr_params cap;
149         t_chsattr_params csp;
150
151         /* open config file */
152         fd=open(mdrun->cfile,O_RDONLY);
153         if(fd<0) {
154                 snprintf(error,128,"%s open cfile %s",ME,mdrun->cfile);
155                 perror(error);
156                 return fd;
157         }
158
159         /* read, parse, set */
160         ret=1;
161         while(ret>0) {
162
163                 /* read */
164                 ret=get_line(fd,line,128);
165
166                 /* parse */
167
168                 // ignore # lines and \n
169                 if((line[0]=='#')|(ret==1))
170                         continue;
171
172                 // reset
173                 memset(&iap,0,sizeof(t_insert_atoms_params));
174                 memset(&cp,0,sizeof(t_continue_params));
175                 memset(&ap,0,sizeof(t_anneal_params));
176                 memset(&cap,0,sizeof(t_chaattr_params));
177                 memset(&csp,0,sizeof(t_chsattr_params));
178
179                 // get command + args
180                 wcnt=0;
181                 while(1) {
182                         if(wcnt)
183                                 wptr=strtok(NULL," \t");
184                         else
185                                 wptr=strtok(line," \t");
186                         if(wptr==NULL)
187                                 break;
188                         strncpy(word[wcnt],wptr,32);
189                         wcnt+=1;
190                 }
191
192                 if(!strncmp(word[0],"potential",9)) {
193                         if(!strncmp(word[1],"albe",4))
194                                 mdrun->potential=MOLDYN_POTENTIAL_AM;
195                         if(!strncmp(word[1],"tersoff",7))
196                                 mdrun->potential=MOLDYN_POTENTIAL_TM;
197                         if(!strncmp(word[1],"ho",2))
198                                 mdrun->potential=MOLDYN_POTENTIAL_HO;
199                         if(!strncmp(word[1],"lj",2))
200                                 mdrun->potential=MOLDYN_POTENTIAL_LJ;
201                 }
202                 else if(!strncmp(word[0],"cutoff",6))
203                         mdrun->cutoff=atof(word[1]);
204                 else if(!strncmp(word[0],"nnd",3))
205                         mdrun->nnd=atof(word[1]);
206                 else if(!strncmp(word[0],"intalgo",7)) {
207                         if(!strncmp(word[1],"verlet",5))
208                                 mdrun->intalgo=MOLDYN_INTEGRATE_VERLET;
209                 }
210                 else if(!strncmp(word[0],"timestep",8))
211                         mdrun->timestep=atof(word[1]);
212                 else if(!strncmp(word[0],"volume",6)) {
213                         mdrun->dim.x=atof(word[1]);
214                         mdrun->dim.y=atof(word[2]);
215                         mdrun->dim.z=atof(word[3]);
216                         if(strncmp(word[4],"0",1))
217                                 mdrun->vis=TRUE;
218                 }
219                 else if(!strncmp(word[0],"pbc",3)) {
220                         if(strncmp(word[1],"0",1))
221                                 mdrun->pbcx=TRUE;
222                         else
223                                 mdrun->pbcx=FALSE;
224                         if(strncmp(word[2],"0",1))
225                                 mdrun->pbcy=TRUE;
226                         else
227                                 mdrun->pbcy=FALSE;
228                         if(strncmp(word[3],"0",1))
229                                 mdrun->pbcz=TRUE;
230                         else
231                                 mdrun->pbcz=FALSE;
232                 }
233                 else if(!strncmp(word[0],"temperature",11))
234                         mdrun->temperature=atof(word[1]);
235                 else if(!strncmp(word[0],"pressure",8))
236                         mdrun->pressure=atof(word[1]);
237                 else if(!strncmp(word[0],"lattice",7)) {
238                         if(!strncmp(word[1],"zincblende",10))
239                                 mdrun->lattice=ZINCBLENDE;
240                         if(!strncmp(word[1],"fcc",3))
241                                 mdrun->lattice=FCC;
242                         if(!strncmp(word[1],"diamond",7))
243                                 mdrun->lattice=DIAMOND;
244                 }
245                 else if(!strncmp(word[0],"element1",8)) {
246                         mdrun->element1=atoi(word[1]);
247                         mdrun->m1=pse_mass[mdrun->element1];
248                 }
249                 else if(!strncmp(word[0],"element2",8)) {
250                         mdrun->element2=atoi(word[1]);
251                         mdrun->m2=pse_mass[mdrun->element2];
252                 }
253                 else if(!strncmp(word[0],"fill",6)) {
254                         // only lc mode by now
255                         mdrun->lx=atoi(word[2]);
256                         mdrun->ly=atoi(word[3]);
257                         mdrun->lz=atoi(word[4]);
258                         mdrun->lc=atof(word[5]);
259                 }
260                 else if(!strncmp(word[0],"aattr",5)) {
261                         // for aatrib line we need a special stage
262                         // containing one schedule of 0 loops ...
263                         for(i=0;i<strlen(word[1]);i++) {
264                                 switch(word[1][i]) {
265                                         case 't':
266                                                 cap.type|=CHAATTR_TOTALV;
267                                                 break;
268                                         case 'r':
269                                                 cap.type|=CHAATTR_REGION;
270                                                 break;
271                                         case 'e':
272                                                 cap.type|=CHAATTR_ELEMENT;
273                                                 break;
274                                         default:
275                                                 break;
276                                 }
277                         }
278                         i=2;
279                         if(cap.type&CHAATTR_REGION) {
280                                 cap.x0=atof(word[1]);   
281                                 cap.y0=atof(word[2]);   
282                                 cap.z0=atof(word[3]);   
283                                 cap.x1=atof(word[4]);   
284                                 cap.y1=atof(word[5]);   
285                                 cap.z1=atof(word[6]);
286                                 i+=6;
287                         }
288                         if(cap.type&CHAATTR_ELEMENT) {
289                                 cap.element=atoi(word[i]);
290                                 i+=1;
291                         }
292                         for(o=0;o<strlen(word[i]);o++) {
293                                 switch(word[i][o]) {
294                                         case 'b':
295                                                 cap.attr|=ATOM_ATTR_VB;
296                                                 break;
297                                         case 'h':
298                                                 cap.attr|=ATOM_ATTR_HB;
299                                                 break;
300                                         case 'v':
301                                                 cap.attr|=ATOM_ATTR_VA;
302                                                 break;
303                                         case 'f':
304                                                 cap.attr|=ATOM_ATTR_FP;
305                                                 break;
306                                         case '1':
307                                                 cap.attr|=ATOM_ATTR_1BP;
308                                                 break;
309                                         case '2':
310                                                 cap.attr|=ATOM_ATTR_2BP;
311                                                 break;
312                                         case '3':
313                                                 cap.attr|=ATOM_ATTR_3BP;
314                                                 break;
315                                         default:
316                                                 break;
317                                 }
318                         }
319                         add_stage(mdrun,STAGE_CHAATTR,&cap);
320                 }
321                 else if(!strncmp(word[0],"sattr",5)) {
322                         // for satrib line we need a special stage
323                         // containing one schedule of 0 loops ...
324                         csp.type=0;
325                         for(i=1;i<wcnt;i++) {
326                                 if(!strncmp(word[i],"pctrl",5)) {
327                                         csp.ptau=0.01/(atof(word[++i])*GPA);
328                                         csp.type|=CHSATTR_PCTRL;
329                                 }
330                                 if(!strncmp(word[i],"tctrl",5)) {
331                                         csp.ttau=atof(word[++i]);
332                                         csp.type|=CHSATTR_TCTRL;
333                                 }
334                                 if(!strncmp(word[i],"prelax",6)) {
335                                         csp.dp=atof(word[++i]);
336                                         csp.type|=CHSATTR_PRELAX;
337                                 }
338                                 if(!strncmp(word[i],"trelax",6)) {
339                                         csp.dt=atof(word[++i]);
340                                         csp.type|=CHSATTR_TRELAX;
341                                 }
342                                 if(!strncmp(word[i],"rsteps",6)) {
343                                         csp.rsteps=atoi(word[++i]);
344                                         csp.type|=CHSATTR_RSTEPS;
345                                 }
346                                 if(!strncmp(word[i],"avgrst",6)) {
347                                         csp.avgrst=atoi(word[++i]);
348                                         csp.type|=CHSATTR_AVGRST;
349                                 }
350                         }
351                         add_stage(mdrun,STAGE_CHSATTR,&csp);
352                 }
353                 else if(!strncmp(word[0],"prerun",6))
354                         mdrun->prerun=atoi(word[1]);
355                 else if(!strncmp(word[0],"avgskip",7))
356                         mdrun->avgskip=atoi(word[1]);
357                 else if(!strncmp(word[0],"elog",4))
358                         mdrun->elog=atoi(word[1]);
359                 else if(!strncmp(word[0],"tlog",4))
360                         mdrun->tlog=atoi(word[1]);
361                 else if(!strncmp(word[0],"plog",4))
362                         mdrun->plog=atoi(word[1]);
363                 else if(!strncmp(word[0],"vlog",4))
364                         mdrun->vlog=atoi(word[1]);
365                 else if(!strncmp(word[0],"save",4))
366                         mdrun->save=atoi(word[1]);
367                 else if(!strncmp(word[0],"visualize",9))
368                         mdrun->visualize=atoi(word[1]);
369                 else if(!strncmp(word[0],"stage",5)) {
370                         // for every stage line, add a stage
371                         if(!strncmp(word[1],"ins_atoms",9)) {
372                                 iap.ins_steps=atoi(word[2]);
373                                 iap.ins_atoms=atoi(word[3]);
374                                 iap.element=atoi(word[4]);
375                                 iap.element=atoi(word[4]);
376                                 iap.brand=atoi(word[5]);
377                                 for(i=0;i<strlen(word[6]);i++) {
378                                         switch(word[6][i]) {
379                                                 case 'b':
380                                                         iap.attr|=ATOM_ATTR_VB;
381                                                         break;
382                                                 case 'h':
383                                                         iap.attr|=ATOM_ATTR_HB;
384                                                         break;
385                                                 case 'v':
386                                                         iap.attr|=ATOM_ATTR_VA;
387                                                         break;
388                                                 case 'f':
389                                                         iap.attr|=ATOM_ATTR_FP;
390                                                         break;
391                                                 case '1':
392                                                         iap.attr|=ATOM_ATTR_1BP;
393                                                         break;
394                                                 case '2':
395                                                         iap.attr|=ATOM_ATTR_2BP;
396                                                         break;
397                                                 case '3':
398                                                         iap.attr|=ATOM_ATTR_3BP;
399                                                         break;
400                                                 default:
401                                                         break;
402                                         }
403                                 }
404                                 // only rand mode by now
405                                 if(word[8][0]=='t')
406                                         iap.type=INS_TOTAL;
407                                 else {
408                                         iap.type=INS_REGION;
409                                         iap.x0=atof(word[8]);
410                                         iap.y0=atof(word[9]);
411                                         iap.z0=atof(word[10]);
412                                         iap.x1=atof(word[11]);
413                                         iap.y1=atof(word[12]);
414                                         iap.z1=atof(word[13]);
415                                 }
416                                 add_stage(mdrun,STAGE_INSERT_ATOMS,&iap);
417                         }
418                         else if(!strncmp(word[1],"continue",8)) {
419                                 cp.runs=atoi(word[2]);
420                                 add_stage(mdrun,STAGE_CONTINUE,&cp);
421                         }
422                         else if(!strncmp(word[1],"anneal",6)) {
423                                 ap.count=0;
424                                 ap.runs=atoi(word[2]);
425                                 ap.dt=atof(word[3]);
426                                 add_stage(mdrun,STAGE_ANNEAL,&ap);
427                         }
428                         else {
429                                 printf("%s unknown stage type: %s\n",
430                                        ME,word[1]);
431                                 return -1;
432                         }
433                 }
434                 else {
435                         printf("%s unknown keyword '%s', skipped!\n",
436                                ME,word[0]);
437                         continue;
438                 }
439         }
440
441         /* close file */
442         close(fd);
443
444         return 0;
445 }
446
447 int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
448
449         double delta;
450
451         if(!(mdrun->sattr&SATTR_PRELAX))
452                 return TRUE;
453
454         delta=moldyn->p_avg-moldyn->p_ref;
455
456         if(delta<0)
457                 delta=-delta;
458
459         if(delta>mdrun->dp)
460                 return FALSE;
461
462         return TRUE;
463 }
464
465 int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
466
467         double delta;
468
469         if(!(mdrun->sattr&SATTR_TRELAX))
470                 return TRUE;
471
472         delta=moldyn->t_avg-moldyn->t_ref;
473
474         if(delta<0)
475                 delta=-delta;
476
477         if(delta>mdrun->dt)
478                 return FALSE;
479
480         return TRUE;
481 }
482
483 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
484
485         t_insert_atoms_params *iap;
486         t_stage *stage;
487         t_atom *atom;
488         t_3dvec r,v,dist;
489         double d,dmin;
490         int cnt,i;
491         double x,y,z;
492         double x0,y0,z0;
493         u8 cr_check,run;
494         
495         stage=mdrun->stage.current->data;
496         iap=stage->params;
497
498         cr_check=FALSE;
499
500         v.x=0.0; v.y=0.0; v.z=0.0;
501
502         switch(iap->type) {
503                 case INS_TOTAL:
504                         x=moldyn->dim.x;
505                         x0=0.0;
506                         y=moldyn->dim.y;
507                         y0=0.0;
508                         z=moldyn->dim.z;
509                         z0=0.0;
510                         cr_check=TRUE;
511                         break;
512                 case INS_REGION:
513                         x=iap->x1-iap->x0;
514                         x0=iap->x0;
515                         y=iap->y1-iap->y0;
516                         y0=iap->y0;
517                         z=iap->z1-iap->z0;
518                         z0=iap->z0;
519                         cr_check=TRUE;
520                         break;
521                 default:
522                         printf("%s unknown insertion mode: %02x\n",
523                                ME,iap->type);
524                         return -1;
525         }
526
527         cnt=0;
528         while(cnt<iap->ins_atoms) {
529                 run=1;
530                 while(run) {
531                         r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0;
532                         r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0;
533                         r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0;
534                         run=0;
535                         if(cr_check==TRUE) {
536                                 dmin=1000;      // for sure too high!
537                                 for(i=0;i<moldyn->count;i++) {
538                                         atom=&(moldyn->atom[i]);
539                                         v3_sub(&dist,&(atom->r),&r);
540                                         check_per_bound(moldyn,&dist);
541                                         d=v3_absolute_square(&dist);
542                                         if(d<(iap->cr*iap->cr)) {
543                                                 run=1;
544                                                 break;
545                                         }
546                                         if(d<dmin)
547                                                 dmin=d;
548                                 }
549                         }
550                 }
551                 add_atom(moldyn,iap->element,pse_mass[iap->element],
552                          iap->brand,iap->attr,&r,&v);
553                 printf("%s atom inserted: %f %f %f | d squared = %f\n",
554                        ME,r.x,r.y,r.z,dmin);
555                 cnt+=1;
556         }
557
558         return 0;
559 }
560
561 int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) {
562
563         t_stage *stage;
564         t_chaattr_params *cap;
565         t_atom *atom;
566         int i;
567
568         stage=mdrun->stage.current->data;
569         cap=stage->params;
570
571         for(i=0;i<moldyn->count;i++) {
572                 atom=&(moldyn->atom[i]);
573                 if(cap->type&CHAATTR_ELEMENT) {
574                         if(cap->element!=atom->element)
575                                 continue;
576                 }
577                 if(cap->type&CHAATTR_REGION) {
578                         if(cap->x0<atom->r.x)
579                                 continue;
580                         if(cap->y0<atom->r.y)
581                                 continue;
582                         if(cap->z0<atom->r.z)
583                                 continue;
584                         if(cap->x1>atom->r.x)
585                                 continue;
586                         if(cap->y1>atom->r.y)
587                                 continue;
588                         if(cap->z1>atom->r.z)
589                                 continue;
590                 }
591                 atom->attr=cap->attr;
592         }
593
594         return 0;
595 }
596
597 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
598
599         t_stage *stage;
600         t_chsattr_params *csp;
601
602         stage=mdrun->stage.current->data;
603         csp=stage->params;
604
605         if(csp->type&CHSATTR_PCTRL) {
606                 if(csp->ptau>0)
607                         set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
608                 else
609                         set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
610         }
611         if(csp->type&CHSATTR_TCTRL) {
612                 if(csp->ttau>0)
613                         set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
614                 else
615                         set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
616         }
617         if(csp->type&CHSATTR_PRELAX) {
618                 if(csp->dp<0)
619                         mdrun->sattr&=(~(SATTR_PRELAX));
620                 else
621                         mdrun->sattr|=SATTR_PRELAX;
622                 mdrun->dp=csp->dp;
623         }
624         if(csp->type&CHSATTR_TRELAX) {
625                 if(csp->dt<0)
626                         mdrun->sattr&=(~(SATTR_TRELAX));
627                 else
628                         mdrun->sattr|=SATTR_TRELAX;
629                 mdrun->dt=csp->dt;
630         }
631         if(csp->type&CHSATTR_AVGRST) {
632                 if(csp->avgrst)
633                         mdrun->sattr|=SATTR_AVGRST;
634                 else
635                         mdrun->sattr&=(~(SATTR_AVGRST));
636         }
637         if(csp->type&CHSATTR_RSTEPS) {
638                 mdrun->relax_steps=csp->rsteps;
639         }
640
641         return 0;
642 }
643
644 #define stage_print(m)  if(!(stage->executed)) \
645                                 printf("%s",m)
646
647 int mdrun_hook(void *ptr1,void *ptr2) {
648
649         t_moldyn *moldyn;
650         t_mdrun *mdrun;
651         t_stage *stage;
652         t_list *sl;
653         int steps;
654         double tau;
655         u8 change_stage;
656
657         t_insert_atoms_params *iap;
658         t_continue_params *cp;
659         t_anneal_params *ap;
660
661         moldyn=ptr1;
662         mdrun=ptr2;
663
664         sl=&(mdrun->stage);
665
666         change_stage=FALSE;
667
668         /* return immediately if there are no more stage */
669         if(sl->current==NULL)
670                 return 0;
671
672         /* get stage description */
673         stage=sl->current->data;
674
675         /* default steps and tau values */
676         steps=mdrun->relax_steps;
677         tau=mdrun->timestep;
678
679         /* check whether relaxation steps are necessary */
680         if(!((check_pressure(moldyn,mdrun)==FALSE)|\
681              (check_temperature(moldyn,mdrun)==FALSE))) {
682         
683                 /* be verbose */
684                 stage_print("\n###########################\n");
685                 stage_print("# [mdrun] executing stage #\n");
686                 stage_print("###########################\n\n");
687                 
688                 /* stage specific stuff */
689                 switch(stage->type) {
690                         case STAGE_INSERT_ATOMS:
691                                 stage_print("  -> insert atoms\n\n");
692                                 iap=stage->params;
693                                 if(iap->cnt_steps==iap->ins_steps) {
694                                         change_stage=TRUE;
695                                         break;
696                                 }
697                                 insert_atoms(moldyn,mdrun);
698                                 iap->cnt_steps+=1;
699                                         break;
700                         case STAGE_CONTINUE:
701                                 stage_print("  -> continue\n\n");
702                                 if(stage->executed==TRUE) {
703                                         change_stage=TRUE;
704                                         break;
705                                 }
706                                 cp=stage->params;
707                                 steps=cp->runs;
708                                         break;
709                         case STAGE_ANNEAL:
710                                 stage_print("  -> anneal\n\n");
711                                 ap=stage->params;
712                                 if(ap->count==ap->runs) {
713                                         change_stage=TRUE;
714                                         break;
715                                 }
716                                 set_temperature(moldyn,moldyn->t_ref+ap->dt);
717                                 ap->count+=1;
718                                 break;
719                         case STAGE_CHAATTR:
720                                 stage_print("  -> chaattr\n\n");
721                                 chaatr(moldyn,mdrun);
722                                 change_stage=TRUE;
723                                 break;
724                         case STAGE_CHSATTR:
725                                 stage_print("  -> chsattr\n\n");
726                                 chsattr(moldyn,mdrun);
727                                 change_stage=TRUE;
728                                 break;
729                         default:
730                                 printf("%s unknwon stage type\n",ME);
731                                 break;
732                 }
733         
734                 /* mark as executed */
735                 stage->executed=TRUE;
736         
737                 /* change state */
738                 if(change_stage==TRUE) {
739                         printf("%s finished stage\n",ME);
740                         if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
741                                 printf("%s no more stages\n",ME);
742                                 return 0;
743                         }
744                         steps=0;
745                         tau=mdrun->timestep;
746                 }
747
748         }
749         else {
750
751                 /* averages */
752                 if(mdrun->sattr&SATTR_AVGRST)
753                         average_reset(moldyn);
754
755         }
756
757         /* continue simulation */
758         moldyn_add_schedule(moldyn,steps,tau);
759
760         return 0;
761 }
762
763 int main(int argc,char **argv) {
764
765         t_mdrun mdrun;
766         t_moldyn moldyn;
767         t_3dvec o;
768
769         /* clear structs */
770         memset(&mdrun,0,sizeof(t_mdrun));
771         memset(&moldyn,0,sizeof(t_moldyn));
772
773         /* parse arguments */
774         if(mdrun_parse_argv(&mdrun,argc,argv)<0)
775                 return -1;
776
777         /* initialize list system */
778         list_init_f(&(mdrun.stage));
779
780         /* parse config file */
781         mdrun_parse_config(&mdrun);
782
783         /* reset the stage list */
784         list_reset_f(&(mdrun.stage));
785
786         /* sanity check! */
787
788         /* prepare simulation */
789         moldyn_init(&moldyn,argc,argv);
790         if(set_int_alg(&moldyn,mdrun.intalgo)<0)
791                 return -1;
792         set_cutoff(&moldyn,mdrun.cutoff);
793         if(set_potential(&moldyn,mdrun.potential)<0)
794                 return -1;
795         switch(mdrun.potential) {
796                 case MOLDYN_POTENTIAL_AM:
797                         albe_mult_set_params(&moldyn,
798                                              mdrun.element1,
799                                              mdrun.element2);
800                         break;
801                 case MOLDYN_POTENTIAL_TM:
802                         tersoff_mult_set_params(&moldyn,
803                                                 mdrun.element1,
804                                                 mdrun.element2);
805                         break;
806                 case MOLDYN_POTENTIAL_HO:
807                         harmonic_oscillator_set_params(&moldyn,mdrun.element1);
808                         break;
809                 case MOLDYN_POTENTIAL_LJ:
810                         lennard_jones_set_params(&moldyn,mdrun.element1);
811                         break;
812                 default:
813                         printf("%s unknown potential: %02x\n",
814                                ME,mdrun.potential);
815                         return -1;
816         }
817         set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
818         set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
819         switch(mdrun.lattice) {
820                 case FCC:
821                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
822                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
823                                        mdrun.ly,mdrun.lz,NULL);
824                         break;
825                 case DIAMOND:
826                         create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1,
827                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
828                                        mdrun.ly,mdrun.lz,NULL);
829                         break;
830                 case ZINCBLENDE:
831                         o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
832                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
833                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
834                                        mdrun.ly,mdrun.lz,&o);
835                         o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
836                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
837                                        mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx,
838                                        mdrun.ly,mdrun.lz,&o);
839                         break;
840                 default:
841                         printf("%s unknown lattice type: %02x\n",
842                                ME,mdrun.lattice);
843                         return -1;
844         }
845         moldyn_bc_check(&moldyn);
846         set_temperature(&moldyn,mdrun.temperature);
847         set_pressure(&moldyn,mdrun.pressure);
848         thermal_init(&moldyn,TRUE);
849         moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
850         moldyn_set_log_dir(&moldyn,mdrun.sdir);
851         moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
852         if(mdrun.elog)
853                 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
854         if(mdrun.tlog)
855                 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
856         if(mdrun.plog)
857                 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
858         if(mdrun.vlog)
859                 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
860         if(mdrun.visualize)
861                 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
862         if(mdrun.save)
863                 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
864         moldyn_set_log(&moldyn,CREATE_REPORT,0);
865         set_avg_skip(&moldyn,mdrun.avgskip);
866
867         /* prepare the hook function */
868         moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
869
870         /* run the simulation */
871         moldyn_integrate(&moldyn);
872
873         /* shutdown */
874         moldyn_shutdown(&moldyn);
875         del_stages(&mdrun);
876         list_destroy_f(&(mdrun.stage));
877
878         return 0;
879 }