]> www.hackdaworld.org Git - physik/posic.git/blob - mdrun.c
e5feda3aa4ceed4c52e53f6327d79738ebd97518
[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])*BAR;
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                                         iap.cr=atof(word[9]);
408                                 }
409                                 else {
410                                         iap.type=INS_REGION;
411                                         iap.x0=atof(word[8]);
412                                         iap.y0=atof(word[9]);
413                                         iap.z0=atof(word[10]);
414                                         iap.x1=atof(word[11]);
415                                         iap.y1=atof(word[12]);
416                                         iap.z1=atof(word[13]);
417                                         iap.cr=atof(word[14]);
418                                 }
419                                 add_stage(mdrun,STAGE_INSERT_ATOMS,&iap);
420                         }
421                         else if(!strncmp(word[1],"continue",8)) {
422                                 cp.runs=atoi(word[2]);
423                                 add_stage(mdrun,STAGE_CONTINUE,&cp);
424                         }
425                         else if(!strncmp(word[1],"anneal",6)) {
426                                 ap.count=0;
427                                 ap.runs=atoi(word[2]);
428                                 ap.dt=atof(word[3]);
429                                 add_stage(mdrun,STAGE_ANNEAL,&ap);
430                         }
431                         else {
432                                 printf("%s unknown stage type: %s\n",
433                                        ME,word[1]);
434                                 return -1;
435                         }
436                 }
437                 else {
438                         printf("%s unknown keyword '%s', skipped!\n",
439                                ME,word[0]);
440                         continue;
441                 }
442         }
443
444         /* close file */
445         close(fd);
446
447         return 0;
448 }
449
450 int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
451
452         double delta;
453
454         if(!(mdrun->sattr&SATTR_PRELAX)) {
455                 return TRUE;
456         }
457
458         delta=moldyn->p_avg-moldyn->p_ref;
459
460         if(delta<0)
461                 delta=-delta;
462
463         if(delta>mdrun->dp)
464                 return FALSE;
465
466         return TRUE;
467 }
468
469 int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
470
471         double delta;
472
473         if(!(mdrun->sattr&SATTR_TRELAX))
474                 return TRUE;
475
476         delta=moldyn->t_avg-moldyn->t_ref;
477
478         if(delta<0)
479                 delta=-delta;
480
481         if(delta>mdrun->dt)
482                 return FALSE;
483
484         return TRUE;
485 }
486
487 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
488
489         t_insert_atoms_params *iap;
490         t_stage *stage;
491         t_atom *atom;
492         t_3dvec r,v,dist;
493         double d,dmin;
494         int cnt,i;
495         double x,y,z;
496         double x0,y0,z0;
497         u8 cr_check,run;
498         
499         stage=mdrun->stage.current->data;
500         iap=stage->params;
501
502         cr_check=FALSE;
503
504         v.x=0.0; v.y=0.0; v.z=0.0;
505
506         switch(iap->type) {
507                 case INS_TOTAL:
508                         x=moldyn->dim.x;
509                         x0=0.0;
510                         y=moldyn->dim.y;
511                         y0=0.0;
512                         z=moldyn->dim.z;
513                         z0=0.0;
514                         cr_check=TRUE;
515                         break;
516                 case INS_REGION:
517                         x=iap->x1-iap->x0;
518                         x0=iap->x0;
519                         y=iap->y1-iap->y0;
520                         y0=iap->y0;
521                         z=iap->z1-iap->z0;
522                         z0=iap->z0;
523                         cr_check=TRUE;
524                         break;
525                 default:
526                         printf("%s unknown insertion mode: %02x\n",
527                                ME,iap->type);
528                         return -1;
529         }
530
531         cnt=0;
532         while(cnt<iap->ins_atoms) {
533                 run=1;
534                 while(run) {
535                         r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0;
536                         r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0;
537                         r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0;
538                         run=0;
539                         if(cr_check==TRUE) {
540                                 dmin=1000;      // for sure too high!
541                                 for(i=0;i<moldyn->count;i++) {
542                                         atom=&(moldyn->atom[i]);
543                                         v3_sub(&dist,&(atom->r),&r);
544                                         check_per_bound(moldyn,&dist);
545                                         d=v3_absolute_square(&dist);
546                                         if(d<(iap->cr*iap->cr)) {
547                                                 run=1;
548                                                 break;
549                                         }
550                                         if(d<dmin)
551                                                 dmin=d;
552                                 }
553                         }
554                 }
555                 add_atom(moldyn,iap->element,pse_mass[iap->element],
556                          iap->brand,iap->attr,&r,&v);
557                 printf("%s atom inserted (%d/%d): %f %f %f\n",
558                        ME,(iap->cnt_steps+1)*iap->ins_atoms,
559                        iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z);
560                 printf("  -> d2 = %f/%f\n",dmin,iap->cr*iap->cr);
561                 cnt+=1;
562         }
563
564         return 0;
565 }
566
567 int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) {
568
569         t_stage *stage;
570         t_chaattr_params *cap;
571         t_atom *atom;
572         int i;
573
574         stage=mdrun->stage.current->data;
575         cap=stage->params;
576
577         for(i=0;i<moldyn->count;i++) {
578                 atom=&(moldyn->atom[i]);
579                 if(cap->type&CHAATTR_ELEMENT) {
580                         if(cap->element!=atom->element)
581                                 continue;
582                 }
583                 if(cap->type&CHAATTR_REGION) {
584                         if(cap->x0<atom->r.x)
585                                 continue;
586                         if(cap->y0<atom->r.y)
587                                 continue;
588                         if(cap->z0<atom->r.z)
589                                 continue;
590                         if(cap->x1>atom->r.x)
591                                 continue;
592                         if(cap->y1>atom->r.y)
593                                 continue;
594                         if(cap->z1>atom->r.z)
595                                 continue;
596                 }
597                 atom->attr=cap->attr;
598         }
599
600         return 0;
601 }
602
603 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
604
605         t_stage *stage;
606         t_chsattr_params *csp;
607
608         stage=mdrun->stage.current->data;
609         csp=stage->params;
610
611         if(csp->type&CHSATTR_PCTRL) {
612                 if(csp->ptau>0)
613                         set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
614                 else
615                         set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
616         }
617         if(csp->type&CHSATTR_TCTRL) {
618                 if(csp->ttau>0)
619                         set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
620                 else
621                         set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
622         }
623         if(csp->type&CHSATTR_PRELAX) {
624                 if(csp->dp<0)
625                         mdrun->sattr&=(~(SATTR_PRELAX));
626                 else
627                         mdrun->sattr|=SATTR_PRELAX;
628                 mdrun->dp=csp->dp;
629         }
630         if(csp->type&CHSATTR_TRELAX) {
631                 if(csp->dt<0)
632                         mdrun->sattr&=(~(SATTR_TRELAX));
633                 else
634                         mdrun->sattr|=SATTR_TRELAX;
635                 mdrun->dt=csp->dt;
636         }
637         if(csp->type&CHSATTR_AVGRST) {
638                 if(csp->avgrst)
639                         mdrun->sattr|=SATTR_AVGRST;
640                 else
641                         mdrun->sattr&=(~(SATTR_AVGRST));
642         }
643         if(csp->type&CHSATTR_RSTEPS) {
644                 mdrun->relax_steps=csp->rsteps;
645         }
646
647         return 0;
648 }
649
650 #define stage_print(m)  if(!(stage->executed)) \
651                                 printf("%s",m)
652
653 int mdrun_hook(void *ptr1,void *ptr2) {
654
655         t_moldyn *moldyn;
656         t_mdrun *mdrun;
657         t_stage *stage;
658         t_list *sl;
659         int steps;
660         double tau;
661         u8 change_stage;
662
663         t_insert_atoms_params *iap;
664         t_continue_params *cp;
665         t_anneal_params *ap;
666
667         moldyn=ptr1;
668         mdrun=ptr2;
669
670         sl=&(mdrun->stage);
671
672         change_stage=FALSE;
673
674         /* return immediately if there are no more stage */
675         if(sl->current==NULL)
676                 return 0;
677
678         /* get stage description */
679         stage=sl->current->data;
680
681         /* default steps and tau values */
682         steps=mdrun->relax_steps;
683         tau=mdrun->timestep;
684
685         /* check whether relaxation steps are necessary */
686         if((check_pressure(moldyn,mdrun)==TRUE)&\
687            (check_temperature(moldyn,mdrun)==TRUE)) {
688         
689                 /* be verbose */
690                 stage_print("\n###########################\n");
691                 stage_print("# [mdrun] executing stage #\n");
692                 stage_print("###########################\n\n");
693                 
694                 /* stage specific stuff */
695                 switch(stage->type) {
696                         case STAGE_INSERT_ATOMS:
697                                 stage_print("  -> insert atoms\n\n");
698                                 iap=stage->params;
699                                 if(iap->cnt_steps==iap->ins_steps) {
700                                         change_stage=TRUE;
701                                         break;
702                                 }
703                                 insert_atoms(moldyn,mdrun);
704                                 iap->cnt_steps+=1;
705                                         break;
706                         case STAGE_CONTINUE:
707                                 stage_print("  -> continue\n\n");
708                                 if(stage->executed==TRUE) {
709                                         change_stage=TRUE;
710                                         break;
711                                 }
712                                 cp=stage->params;
713                                 steps=cp->runs;
714                                         break;
715                         case STAGE_ANNEAL:
716                                 stage_print("  -> anneal\n\n");
717                                 ap=stage->params;
718                                 if(ap->count==ap->runs) {
719                                         change_stage=TRUE;
720                                         break;
721                                 }
722                                 set_temperature(moldyn,moldyn->t_ref+ap->dt);
723                                 ap->count+=1;
724                                 break;
725                         case STAGE_CHAATTR:
726                                 stage_print("  -> chaattr\n\n");
727                                 chaatr(moldyn,mdrun);
728                                 change_stage=TRUE;
729                                 break;
730                         case STAGE_CHSATTR:
731                                 stage_print("  -> chsattr\n\n");
732                                 chsattr(moldyn,mdrun);
733                                 change_stage=TRUE;
734                                 break;
735                         default:
736                                 printf("%s unknwon stage type\n",ME);
737                                 break;
738                 }
739         
740                 /* mark as executed */
741                 stage->executed=TRUE;
742         
743                 /* change state */
744                 if(change_stage==TRUE) {
745                         printf("%s finished stage\n",ME);
746                         if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
747                                 printf("%s no more stages\n",ME);
748                                 return 0;
749                         }
750                         steps=0;
751                         tau=mdrun->timestep;
752                 }
753
754         }
755         else {
756
757                 /* averages */
758                 if(mdrun->sattr&SATTR_AVGRST)
759                         average_reset(moldyn);
760
761         }
762
763         /* continue simulation */
764         moldyn_add_schedule(moldyn,steps,tau);
765
766         return 0;
767 }
768
769 int main(int argc,char **argv) {
770
771         t_mdrun mdrun;
772         t_moldyn moldyn;
773         t_3dvec o;
774
775         /* clear structs */
776         memset(&mdrun,0,sizeof(t_mdrun));
777         memset(&moldyn,0,sizeof(t_moldyn));
778
779         /* parse arguments */
780         if(mdrun_parse_argv(&mdrun,argc,argv)<0)
781                 return -1;
782
783         /* initialize list system */
784         list_init_f(&(mdrun.stage));
785
786         /* parse config file */
787         mdrun_parse_config(&mdrun);
788
789         /* reset the stage list */
790         list_reset_f(&(mdrun.stage));
791
792         /* sanity check! */
793
794         /* prepare simulation */
795         moldyn_init(&moldyn,argc,argv);
796         if(set_int_alg(&moldyn,mdrun.intalgo)<0)
797                 return -1;
798         set_cutoff(&moldyn,mdrun.cutoff);
799         if(set_potential(&moldyn,mdrun.potential)<0)
800                 return -1;
801         switch(mdrun.potential) {
802                 case MOLDYN_POTENTIAL_AM:
803                         albe_mult_set_params(&moldyn,
804                                              mdrun.element1,
805                                              mdrun.element2);
806                         break;
807                 case MOLDYN_POTENTIAL_TM:
808                         tersoff_mult_set_params(&moldyn,
809                                                 mdrun.element1,
810                                                 mdrun.element2);
811                         break;
812                 case MOLDYN_POTENTIAL_HO:
813                         harmonic_oscillator_set_params(&moldyn,mdrun.element1);
814                         break;
815                 case MOLDYN_POTENTIAL_LJ:
816                         lennard_jones_set_params(&moldyn,mdrun.element1);
817                         break;
818                 default:
819                         printf("%s unknown potential: %02x\n",
820                                ME,mdrun.potential);
821                         return -1;
822         }
823         set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
824         set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
825         switch(mdrun.lattice) {
826                 case FCC:
827                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
828                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
829                                        mdrun.ly,mdrun.lz,NULL);
830                         break;
831                 case DIAMOND:
832                         create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1,
833                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
834                                        mdrun.ly,mdrun.lz,NULL);
835                         break;
836                 case ZINCBLENDE:
837                         o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
838                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
839                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
840                                        mdrun.ly,mdrun.lz,&o);
841                         o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
842                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
843                                        mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx,
844                                        mdrun.ly,mdrun.lz,&o);
845                         break;
846                 default:
847                         printf("%s unknown lattice type: %02x\n",
848                                ME,mdrun.lattice);
849                         return -1;
850         }
851         moldyn_bc_check(&moldyn);
852         set_temperature(&moldyn,mdrun.temperature);
853         set_pressure(&moldyn,mdrun.pressure);
854         thermal_init(&moldyn,TRUE);
855         moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
856         moldyn_set_log_dir(&moldyn,mdrun.sdir);
857         moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
858         if(mdrun.elog)
859                 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
860         if(mdrun.tlog)
861                 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
862         if(mdrun.plog)
863                 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
864         if(mdrun.vlog)
865                 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
866         if(mdrun.visualize)
867                 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
868         if(mdrun.save)
869                 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
870         moldyn_set_log(&moldyn,CREATE_REPORT,0);
871         set_avg_skip(&moldyn,mdrun.avgskip);
872
873         /* prepare the hook function */
874         moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
875
876         /* run the simulation */
877         moldyn_integrate(&moldyn);
878
879         /* shutdown */
880         moldyn_shutdown(&moldyn);
881         del_stages(&mdrun);
882         list_destroy_f(&(mdrun.stage));
883
884         return 0;
885 }