]> www.hackdaworld.org Git - physik/posic.git/blob - mdrun.c
save in src dir + app(endix)
[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_DISPLACE_ATOM:
94                         psize=sizeof(t_displace_atom_params);
95                         break;
96                 case STAGE_DEL_ATOMS:
97                         psize=sizeof(t_del_atoms_params);
98                         break;
99                 case STAGE_MODIFY_ATOMS:
100                         psize=sizeof(t_modify_atoms_params);
101                         break;
102                 case STAGE_INSERT_ATOMS:
103                         psize=sizeof(t_insert_atoms_params);
104                         break;
105                 case STAGE_INSERT_MIXED_ATOMS:
106                         psize=sizeof(t_insert_mixed_atoms_params);
107                         break;
108                 case STAGE_CONTINUE:
109                         psize=sizeof(t_continue_params);
110                         break;
111                 case STAGE_ANNEAL:
112                         psize=sizeof(t_anneal_params);
113                         break;
114                 case STAGE_CHAATTR:
115                         psize=sizeof(t_chaattr_params);
116                         break;
117                 case STAGE_CHSATTR:
118                         psize=sizeof(t_chsattr_params);
119                         break;
120                 case STAGE_SET_TEMP:
121                         psize=sizeof(t_set_temp_params);
122                         break;
123                 case STAGE_SET_TIMESTEP:
124                         psize=sizeof(t_set_timestep_params);
125                         break;
126                 case STAGE_FILL:
127                         psize=sizeof(t_fill_params);
128                         break;
129                 case STAGE_THERMAL_INIT:
130                         psize=0;
131                         break;
132                 default:
133                         printf("%s unknown stage type: %02x\n",ME,type);
134                         return -1;
135         }
136
137         stage=malloc(sizeof(t_stage));
138         if(stage==NULL) {
139                 perror("[mdrun] malloc (add stage)");
140                 return -1;
141         }
142
143         stage->type=type;
144         stage->executed=FALSE;
145
146         stage->params=malloc(psize);
147         if(stage->params==NULL) {
148                 perror("[mdrun] malloc (stage params)");
149                 return -1;
150         }
151
152         memcpy(stage->params,params,psize);
153
154         list_add_immediate_f(&(mdrun->stage),stage);
155
156         return 0;
157 }
158
159 int mdrun_parse_config(t_mdrun *mdrun) {
160
161         int fd,ret;
162         char error[128];
163         char line[128];
164         char *wptr;
165         char word[32][64];
166         int wcnt;
167         int i,o;
168
169         t_displace_atom_params dap;
170         t_modify_atoms_params map;
171         t_insert_atoms_params iap;
172         t_insert_mixed_atoms_params imp;
173         t_continue_params cp;
174         t_anneal_params ap;
175         t_chaattr_params cap;
176         t_chsattr_params csp;
177         t_set_temp_params stp;
178         t_set_timestep_params stsp;
179         t_fill_params fp;
180         t_del_atoms_params delp;
181
182         /* open config file */
183         fd=open(mdrun->cfile,O_RDONLY);
184         if(fd<0) {
185                 snprintf(error,128,"%s open cfile %s",ME,mdrun->cfile);
186                 perror(error);
187                 return fd;
188         }
189
190         /* read, parse, set */
191         ret=1;
192         while(ret>0) {
193
194                 /* read */
195                 ret=get_line(fd,line,128);
196
197                 /* parse */
198
199                 // ignore # lines and \n
200                 if((line[0]=='#')|(ret==1))
201                         continue;
202
203                 // reset
204                 memset(&iap,0,sizeof(t_insert_atoms_params));
205                 memset(&map,0,sizeof(t_modify_atoms_params));
206                 memset(&imp,0,sizeof(t_insert_mixed_atoms_params));
207                 memset(&cp,0,sizeof(t_continue_params));
208                 memset(&ap,0,sizeof(t_anneal_params));
209                 memset(&cap,0,sizeof(t_chaattr_params));
210                 memset(&csp,0,sizeof(t_chsattr_params));
211                 memset(&stp,0,sizeof(t_set_temp_params));
212                 memset(&stsp,0,sizeof(t_set_timestep_params));
213                 memset(&fp,0,sizeof(t_fill_params));
214                 memset(&delp,0,sizeof(t_del_atoms_params));
215
216                 // get command + args
217                 wcnt=0;
218                 while(1) {
219                         if(wcnt)
220                                 wptr=strtok(NULL," \t");
221                         else
222                                 wptr=strtok(line," \t");
223                         if(wptr==NULL)
224                                 break;
225                         strncpy(word[wcnt],wptr,64);
226                         wcnt+=1;
227                 }
228
229                 if(!strncmp(word[0],"potential",9)) {
230                         if(!strncmp(word[1],"albe",4))
231                                 mdrun->potential=MOLDYN_POTENTIAL_AM;
232                         if(!strncmp(word[1],"tersoff",7))
233                                 mdrun->potential=MOLDYN_POTENTIAL_TM;
234                         if(!strncmp(word[1],"ho",2))
235                                 mdrun->potential=MOLDYN_POTENTIAL_HO;
236                         if(!strncmp(word[1],"lj",2))
237                                 mdrun->potential=MOLDYN_POTENTIAL_LJ;
238                 }
239                 else if(!strncmp(word[0],"continue",8))
240                         strncpy(mdrun->continue_file,word[1],128);
241                 else if(!strncmp(word[0],"cutoff",6))
242                         mdrun->cutoff=atof(word[1]);
243                 else if(!strncmp(word[0],"nnd",3))
244                         mdrun->nnd=atof(word[1]);
245                 else if(!strncmp(word[0],"intalgo",7)) {
246                         if(!strncmp(word[1],"verlet",5))
247                                 mdrun->intalgo=MOLDYN_INTEGRATE_VERLET;
248                 }
249                 else if(!strncmp(word[0],"timestep",8))
250                         mdrun->timestep=atof(word[1]);
251                 else if(!strncmp(word[0],"volume",6)) {
252                         mdrun->dim.x=atof(word[1]);
253                         mdrun->dim.y=atof(word[2]);
254                         mdrun->dim.z=atof(word[3]);
255                         if(strncmp(word[4],"0",1))
256                                 mdrun->vis=TRUE;
257                 }
258                 else if(!strncmp(word[0],"pbc",3)) {
259                         if(strncmp(word[1],"0",1))
260                                 mdrun->pbcx=TRUE;
261                         else
262                                 mdrun->pbcx=FALSE;
263                         if(strncmp(word[2],"0",1))
264                                 mdrun->pbcy=TRUE;
265                         else
266                                 mdrun->pbcy=FALSE;
267                         if(strncmp(word[3],"0",1))
268                                 mdrun->pbcz=TRUE;
269                         else
270                                 mdrun->pbcz=FALSE;
271                 }
272                 else if(!strncmp(word[0],"temperature",11))
273                         mdrun->temperature=atof(word[1]);
274                 else if(!strncmp(word[0],"pressure",8))
275                         mdrun->pressure=atof(word[1]);
276                 else if(!strncmp(word[0],"lattice",7)) {
277                         if(!strncmp(word[1],"zincblende",10))
278                                 mdrun->lattice=ZINCBLENDE;
279                         if(!strncmp(word[1],"fcc",3))
280                                 mdrun->lattice=FCC;
281                         if(!strncmp(word[1],"diamond",7))
282                                 mdrun->lattice=DIAMOND;
283                         if(!strncmp(word[1],"none",4))
284                                 mdrun->lattice=NONE;
285                         if(wcnt==3)
286                                 mdrun->lc=atof(word[2]);
287                 }
288                 else if(!strncmp(word[0],"element1",8)) {
289                         mdrun->element1=atoi(word[1]);
290                 }
291                 else if(!strncmp(word[0],"element2",8)) {
292                         mdrun->element2=atoi(word[1]);
293                 }
294                 else if(!strncmp(word[0],"fill",6)) {
295                         // default values
296                         fp.fill_element=mdrun->element1;
297                         fp.fill_brand=0;
298                         fp.lattice=mdrun->lattice;
299                         fp.p_params.type=0;
300                         fp.d_params.type=0;
301                         fp.d_params.stype=0;
302                         // parse fill command
303                         i=1;
304                         while(i<wcnt) {
305                                 if(!strncmp(word[i],"lc",2)) {
306                                         fp.lx=atoi(word[++i]);
307                                         fp.ly=atoi(word[++i]);
308                                         fp.lz=atoi(word[++i]);
309                                         fp.lc=atof(word[++i]);
310                                         mdrun->lc=fp.lc;
311                                         continue;
312                                 }
313                                 if(!strncmp(word[i],"eb",2)) {
314                                         fp.fill_element=atoi(word[++i]);
315                                         fp.fill_brand=atoi(word[++i]);
316                                         continue;
317                                 }
318                                 if(word[i][0]=='p') {
319                                         i+=1;
320                                         switch(word[i][0]) {
321                                         case 'i':
322                                                 if(word[i][1]=='r')
323                                                 fp.p_params.type=PART_INSIDE_R;
324                                                 else
325                                                 fp.p_params.type=PART_INSIDE_D;
326                                                 break;
327                                         case 'o':
328                                                 if(word[i][1]=='r')
329                                                 fp.p_params.type=PART_OUTSIDE_R;
330                                                 else
331                                                 fp.p_params.type=PART_OUTSIDE_D;
332                                                 break;
333                                         default:
334                                                 break;
335                                         }
336                                         if((fp.p_params.type==PART_INSIDE_R)||
337                                           (fp.p_params.type==PART_OUTSIDE_R)) {
338                                                 fp.p_params.r=atof(word[++i]);  
339                                                 fp.p_params.p.x=atof(word[++i]);
340                                                 fp.p_params.p.y=atof(word[++i]);
341                                                 fp.p_params.p.z=atof(word[++i]);
342                                         }
343                                         if((fp.p_params.type==PART_INSIDE_D)||
344                                            (fp.p_params.type==PART_OUTSIDE_D)) {
345                                                 fp.p_params.p.x=atof(word[++i]);
346                                                 fp.p_params.p.y=atof(word[++i]);
347                                                 fp.p_params.p.z=atof(word[++i]);
348                                                 fp.p_params.d.x=atof(word[++i]);
349                                                 fp.p_params.d.y=atof(word[++i]);
350                                                 fp.p_params.d.z=atof(word[++i]);
351                                         }
352                                         continue;
353                                 }
354                                 if(word[i][0]=='d') {
355                                         switch(word[++i][0]) {
356                                                 case '0':
357
358                                 fp.d_params.type=DEFECT_TYPE_0D;
359                                 if(!strncmp(word[i+1],"dbx",3)) {
360                                         fp.d_params.stype=DEFECT_STYPE_DB_X;
361                                 }
362                                 if(!strncmp(word[i+1],"dby",3)) {
363                                         fp.d_params.stype=DEFECT_STYPE_DB_Y;
364                                 }
365                                 if(!strncmp(word[i+1],"dbz",3)) {
366                                         fp.d_params.stype=DEFECT_STYPE_DB_Z;
367                                 }
368                                 if(!strncmp(word[i+1],"dbr",3)) {
369                                         fp.d_params.stype=DEFECT_STYPE_DB_R;
370                                 }
371                                 i+=1;
372                                 fp.d_params.od=atof(word[++i]);
373                                 fp.d_params.dd=atof(word[++i]);
374                                 fp.d_params.element=atoi(word[++i]);
375                                 fp.d_params.brand=atoi(word[++i]);
376                                 // parsed in future
377                 fp.d_params.attr=ATOM_ATTR_HB|ATOM_ATTR_VA;
378                 fp.d_params.attr|=ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP;
379                                 break;
380
381                                                 case '1':
382                                 fp.d_params.type=DEFECT_TYPE_1D;
383                                 break;
384                                                 case '2':
385                                 fp.d_params.type=DEFECT_TYPE_2D;
386                                 break;
387                                                 case '3':
388                                 fp.d_params.type=DEFECT_TYPE_3D;
389                                 break;
390                                                 default:
391                                                         break;
392                                         }
393                                         continue;
394
395                                 }
396                                 // offset
397                                 if(word[i][0]=='o') {
398                                         fp.o_params.o.x=atof(word[++i])*fp.lc;
399                                         fp.o_params.o.y=atof(word[++i])*fp.lc;
400                                         fp.o_params.o.z=atof(word[++i])*fp.lc;
401                                         fp.o_params.use=1;
402                                         continue;
403                                 }
404                                 i+=1;
405                         }
406                         add_stage(mdrun,STAGE_FILL,&fp);
407                 }
408                 else if(!strncmp(word[0],"thermal_init",12)) {
409                         add_stage(mdrun,STAGE_THERMAL_INIT,NULL);
410                 }
411                 else if(!strncmp(word[0],"aattr",5)) {
412                         // for aatrib line we need a special stage
413                         // containing one schedule of 0 loops ...
414                         for(i=0;i<strlen(word[1]);i++) {
415                                 switch(word[1][i]) {
416                                         case 't':
417                                                 cap.type|=CHAATTR_TOTALV;
418                                                 break;
419                                         case 'r':
420                                                 cap.type|=CHAATTR_REGION;
421                                                 break;
422                                         case 'e':
423                                                 cap.type|=CHAATTR_ELEMENT;
424                                                 break;
425                                         default:
426                                                 break;
427                                 }
428                         }
429                         i=2;
430                         if(cap.type&CHAATTR_REGION) {
431                                 cap.x0=atof(word[2]);   
432                                 cap.y0=atof(word[3]);   
433                                 cap.z0=atof(word[4]);   
434                                 cap.x1=atof(word[5]);   
435                                 cap.y1=atof(word[6]);   
436                                 cap.z1=atof(word[7]);
437                                 i+=6;
438                         }
439                         if(cap.type&CHAATTR_ELEMENT) {
440                                 cap.element=atoi(word[i]);
441                                 i+=1;
442                         }
443                         for(o=0;o<strlen(word[i]);o++) {
444                                 switch(word[i][o]) {
445                                         case 'b':
446                                                 cap.attr|=ATOM_ATTR_VB;
447                                                 break;
448                                         case 'h':
449                                                 cap.attr|=ATOM_ATTR_HB;
450                                                 break;
451                                         case 'v':
452                                                 cap.attr|=ATOM_ATTR_VA;
453                                                 break;
454                                         case 'f':
455                                                 cap.attr|=ATOM_ATTR_FP;
456                                                 break;
457                                         case '1':
458                                                 cap.attr|=ATOM_ATTR_1BP;
459                                                 break;
460                                         case '2':
461                                                 cap.attr|=ATOM_ATTR_2BP;
462                                                 break;
463                                         case '3':
464                                                 cap.attr|=ATOM_ATTR_3BP;
465                                                 break;
466                                         default:
467                                                 break;
468                                 }
469                         }
470                         add_stage(mdrun,STAGE_CHAATTR,&cap);
471                 }
472                 else if(!strncmp(word[0],"sattr",5)) {
473                         // for satrib line we need a special stage
474                         // containing one schedule of 0 loops ...
475                         csp.type=0;
476                         for(i=1;i<wcnt;i++) {
477                                 if(!strncmp(word[i],"pctrl",5)) {
478                                         csp.ptau=atof(word[++i]);
479                                         if(csp.ptau>0)
480                                                 csp.ptau=0.01/(csp.ptau*GPA);
481                                         csp.type|=CHSATTR_PCTRL;
482                                 }
483                                 if(!strncmp(word[i],"tctrl",5)) {
484                                         csp.ttau=atof(word[++i]);
485                                         csp.type|=CHSATTR_TCTRL;
486                                 }
487                                 if(!strncmp(word[i],"prelax",6)) {
488                                         csp.dp=atof(word[++i])*BAR;
489                                         csp.type|=CHSATTR_PRELAX;
490                                 }
491                                 if(!strncmp(word[i],"trelax",6)) {
492                                         csp.dt=atof(word[++i]);
493                                         csp.type|=CHSATTR_TRELAX;
494                                 }
495                                 if(!strncmp(word[i],"rsteps",6)) {
496                                         csp.rsteps=atoi(word[++i]);
497                                         csp.type|=CHSATTR_RSTEPS;
498                                 }
499                                 if(!strncmp(word[i],"avgrst",6)) {
500                                         csp.avgrst=atoi(word[++i]);
501                                         csp.type|=CHSATTR_AVGRST;
502                                 }
503                         }
504                         add_stage(mdrun,STAGE_CHSATTR,&csp);
505                 }
506                 else if(!strncmp(word[0],"prerun",6))
507                         mdrun->prerun=atoi(word[1]);
508                 else if(!strncmp(word[0],"avgskip",7))
509                         mdrun->avgskip=atoi(word[1]);
510                 else if(!strncmp(word[0],"elog",4))
511                         mdrun->elog=atoi(word[1]);
512                 else if(!strncmp(word[0],"tlog",4))
513                         mdrun->tlog=atoi(word[1]);
514                 else if(!strncmp(word[0],"plog",4))
515                         mdrun->plog=atoi(word[1]);
516                 else if(!strncmp(word[0],"vlog",4))
517                         mdrun->vlog=atoi(word[1]);
518                 else if(!strncmp(word[0],"save",4))
519                         mdrun->save=atoi(word[1]);
520                 else if(!strncmp(word[0],"visualize",9))
521                         mdrun->visualize=atoi(word[1]);
522                 else if(!strncmp(word[0],"stage",5)) {
523                         // for every stage line, add a stage
524                         if(!strncmp(word[1],"displace",8)) {
525                                 dap.nr=atoi(word[2]);   
526                                 dap.dx=atof(word[3]);
527                                 dap.dy=atof(word[4]);
528                                 dap.dz=atof(word[5]);
529                                 add_stage(mdrun,STAGE_DISPLACE_ATOM,&dap);
530                         }
531                         else if(!strncmp(word[1],"del_atoms",9)) {
532                                 delp.o.x=atof(word[2]);
533                                 delp.o.y=atof(word[3]);
534                                 delp.o.z=atof(word[4]);
535                                 delp.r=atof(word[5]);
536                                 add_stage(mdrun,STAGE_DEL_ATOMS,&delp);
537                         }
538                         else if(!strncmp(word[1],"mod_atoms",8)) {
539                                 i=2;
540                                 while(i<wcnt) {
541                                         if(!strncmp(word[i],"t",1)) {
542                                                 map.tag=atoi(word[++i]);
543                                                 i+=1;
544                                         }
545                                         if(!strncmp(word[i],"ekin",5)) {
546                                                 map.ekin.x=atof(word[++i])*EV;
547                                                 map.ekin.y=atof(word[++i])*EV;
548                                                 map.ekin.z=atof(word[++i])*EV;
549                                                 i+=1;
550                                         }
551                                 }
552                                 add_stage(mdrun,STAGE_MODIFY_ATOMS,&map);
553                         }
554                         else if(!strncmp(word[1],"ins_atoms",9)) {
555                                 iap.ins_steps=atoi(word[2]);
556                                 iap.ins_atoms=atoi(word[3]);
557                                 iap.element=atoi(word[4]);
558                                 iap.brand=atoi(word[5]);
559                                 for(i=0;i<strlen(word[6]);i++) {
560                                         switch(word[6][i]) {
561                                                 case 'b':
562                                                         iap.attr|=ATOM_ATTR_VB;
563                                                         break;
564                                                 case 'h':
565                                                         iap.attr|=ATOM_ATTR_HB;
566                                                         break;
567                                                 case 'v':
568                                                         iap.attr|=ATOM_ATTR_VA;
569                                                         break;
570                                                 case 'f':
571                                                         iap.attr|=ATOM_ATTR_FP;
572                                                         break;
573                                                 case '1':
574                                                         iap.attr|=ATOM_ATTR_1BP;
575                                                         break;
576                                                 case '2':
577                                                         iap.attr|=ATOM_ATTR_2BP;
578                                                         break;
579                                                 case '3':
580                                                         iap.attr|=ATOM_ATTR_3BP;
581                                                         break;
582                                                 default:
583                                                         break;
584                                         }
585                                 }
586                                 switch(word[7][0]) {
587                                         // insertion types
588                                         case 'p':
589                                                 iap.type=INS_POS;
590                                                 iap.x0=atof(word[8]);
591                                                 iap.y0=atof(word[9]);
592                                                 iap.z0=atof(word[10]);
593                                                 break;
594                                         case 'r':
595                                                 switch(word[8][0]) {
596
597                                                 case 't':
598                                                         iap.type=INS_TOTAL;
599                                                         iap.cr=atof(word[9]);
600                                                         break;
601                                                 case 'r':
602                                                         iap.type=INS_RECT;
603                                                         iap.x0=atof(word[9]);
604                                                         iap.y0=atof(word[10]);
605                                                         iap.z0=atof(word[11]);
606                                                         iap.x1=atof(word[12]);
607                                                         iap.y1=atof(word[13]);
608                                                         iap.z1=atof(word[14]);
609                                                         iap.cr=atof(word[15]);
610                                                         break;
611                                                 case 's':
612                                                         iap.type=INS_SPHERE;
613                                                         iap.x0=atof(word[9]);
614                                                         iap.y0=atof(word[10]);
615                                                         iap.z0=atof(word[11]);
616                                                         iap.x1=atof(word[12]);
617                                                         iap.cr=atof(word[13]);
618                                                         break;
619                                                 default:
620                                                         break;
621                                                 }
622                                         default:
623                                                 break;
624                                 }
625                                 add_stage(mdrun,STAGE_INSERT_ATOMS,&iap);
626                         }
627
628
629                         // HERE WE GO ...
630
631                         else if(!strncmp(word[1],"ins_m_atoms",11)) {
632                                 imp.element1=atoi(word[2]);
633                                 imp.element2=atoi(word[3]);
634                                 imp.amount1=atoi(word[4]);
635                                 imp.amount2=atoi(word[5]);
636                                 imp.brand1=atoi(word[6]);
637                                 imp.brand2=atoi(word[7]);
638                                 imp.crmin=atof(word[8]);
639                                 imp.crmax=atof(word[9]);
640                                 /* do this later ...
641                                 for(i=0;i<strlen(word[8]);i++) {
642                                         switch(word[8][i]) {
643                                                 case 'b':
644                                                         imp.attr|=ATOM_ATTR_VB;
645                                                         break;
646                                                 case 'h':
647                                                         imp.attr|=ATOM_ATTR_HB;
648                                                         break;
649                                                 case 'v':
650                                                         imp.attr|=ATOM_ATTR_VA;
651                                                         break;
652                                                 case 'f':
653                                                         imp.attr|=ATOM_ATTR_FP;
654                                                         break;
655                                                 case '1':
656                                                         imp.attr|=ATOM_ATTR_1BP;
657                                                         break;
658                                                 case '2':
659                                                         imp.attr|=ATOM_ATTR_2BP;
660                                                         break;
661                                                 case '3':
662                                                         imp.attr|=ATOM_ATTR_3BP;
663                                                         break;
664                                                 default:
665                                                         break;
666                                         }
667                                 }
668                                 */
669                                 imp.attr1=ATOM_ATTR_HB|ATOM_ATTR_VA|ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_FP;
670                                 imp.attr2=ATOM_ATTR_HB|ATOM_ATTR_VA|ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_FP;
671                                 add_stage(mdrun,STAGE_INSERT_MIXED_ATOMS,&imp);
672                         }
673
674
675
676
677
678                         else if(!strncmp(word[1],"continue",8)) {
679                                 cp.runs=atoi(word[2]);
680                                 add_stage(mdrun,STAGE_CONTINUE,&cp);
681                         }
682                         else if(!strncmp(word[1],"anneal",6)) {
683                                 ap.count=0;
684                                 ap.runs=atoi(word[2]);
685                                 ap.dt=atof(word[3]);
686                                 ap.interval=atoi(word[4]);
687                                 add_stage(mdrun,STAGE_ANNEAL,&ap);
688                         }
689                         else if(!strncmp(word[1],"set_temp",8)) {
690                                 if(word[2][0]=='c') {
691                                         stp.type=SET_TEMP_CURRENT;
692                                         stp.val=0.0;
693                                 }
694                                 else {
695                                         stp.type=SET_TEMP_VALUE;
696                                         stp.val=atof(word[2]);
697                                 }
698                                 add_stage(mdrun,STAGE_SET_TEMP,&stp);
699                         }
700                         else if(!strncmp(word[1],"set_timestep",12)) {
701                                 stsp.tau=atof(word[2]);
702                                 add_stage(mdrun,STAGE_SET_TIMESTEP,&stsp);
703                         }
704                         else {
705                                 printf("%s unknown stage type: %s\n",
706                                        ME,word[1]);
707                                 return -1;
708                         }
709                 }
710                 else {
711                         printf("%s unknown keyword '%s', skipped!\n",
712                                ME,word[0]);
713                         continue;
714                 }
715         }
716
717         /* close file */
718         close(fd);
719
720         return 0;
721 }
722
723 int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
724
725         double delta;
726
727         if(!(mdrun->sattr&SATTR_PRELAX)) {
728                 return TRUE;
729         }
730
731         delta=moldyn->p_avg-moldyn->p_ref;
732
733         if(delta<0)
734                 delta=-delta;
735
736         if(delta>mdrun->dp)
737                 return FALSE;
738
739         return TRUE;
740 }
741
742 int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
743
744         double delta;
745
746         if(!(mdrun->sattr&SATTR_TRELAX))
747                 return TRUE;
748
749         delta=moldyn->t_avg-moldyn->t_ref;
750
751         if(delta<0)
752                 delta=-delta;
753
754         if(delta>mdrun->dt)
755                 return FALSE;
756
757         return TRUE;
758 }
759
760 int displace_atom(t_moldyn *moldyn,t_mdrun *mdrun) {
761
762         t_displace_atom_params *dap;
763         t_stage *stage;
764         t_atom *atom;
765
766         stage=mdrun->stage.current->data;
767         dap=stage->params;
768
769         atom=&(moldyn->atom[dap->nr]);
770         atom->r.x+=dap->dx;
771         atom->r.y+=dap->dy;
772         atom->r.z+=dap->dz;
773
774         return 0;
775 }
776
777 int del_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
778
779         t_stage *stage;
780         t_del_atoms_params *delp;
781         int i;
782         t_3dvec dist;
783         
784         stage=mdrun->stage.current->data;
785         delp=stage->params;
786
787         for(i=0;i<moldyn->count;i++) {
788                 v3_sub(&dist,&(delp->o),&(moldyn->atom[i].r));
789 //printf("%d ----> %f %f %f = %f | %f\n",i,dist.x,dist.y,dist.z,v3_absolute_square(&dist),delp->r*delp->r);
790                 if(v3_absolute_square(&dist)<=(delp->r*delp->r)) {
791                         del_atom(moldyn,moldyn->atom[i].tag);
792                         printf("%s atom deleted: %d %d %d\n",ME,
793                                moldyn->atom[i].tag,moldyn->atom[i].element,
794                                moldyn->atom[i].brand);
795                 }
796         }
797
798         return 0;
799 }
800
801 int modify_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
802
803         t_modify_atoms_params *map;
804         t_stage *stage;
805         t_atom *atom;
806         t_3dvec v;
807         int i;
808
809         atom=moldyn->atom;
810         stage=mdrun->stage.current->data;
811         map=stage->params;
812         v.x=0.0; v.y=0.0; v.z=0.0;
813
814         for(i=0;i<moldyn->count;i++) {
815                 if(atom[i].tag==map->tag) {
816                         v.x=sqrt(2.0*fabs(map->ekin.x)/atom[i].mass);
817                         if(map->ekin.x<0.0)
818                                 v.x=-v.x;
819                         v.y=sqrt(2.0*fabs(map->ekin.y)/atom[i].mass);
820                         if(map->ekin.y<0.0)
821                                 v.y=-v.y;
822                         v.z=sqrt(2.0*fabs(map->ekin.z)/atom[i].mass);
823                         if(map->ekin.z<0.0)
824                                 v.z=-v.z;
825                         v3_copy(&(atom[i].v),&v);
826                         printf("%s atom modified: v = (%f %f %f)\n",
827                                ME,v.x,v.y,v.z);
828                 }
829         }       
830
831         return 0;
832 }
833
834 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
835
836         t_insert_atoms_params *iap;
837         t_stage *stage;
838         t_atom *atom;
839         t_3dvec r,v,dist;
840         double d,dmin,o;
841         int cnt,i;
842         double x,y,z;
843         double x0,y0,z0;
844         u8 cr_check,run;
845         
846         stage=mdrun->stage.current->data;
847         iap=stage->params;
848
849         cr_check=FALSE;
850
851         v.x=0.0; v.y=0.0; v.z=0.0;
852         x=0; y=0; z=0;
853         o=0; dmin=0;
854
855         switch(mdrun->lattice) {
856                 case CUBIC:
857                         o=0.5*mdrun->lc;
858                         break;
859                 case FCC:
860                         o=0.25*mdrun->lc;
861                         break;
862                 case DIAMOND:
863                 case ZINCBLENDE:
864                         o=0.125*mdrun->lc;
865                         break;
866                 default:
867                         break;
868         }
869
870         switch(iap->type) {
871                 case INS_TOTAL:
872                         x=moldyn->dim.x;
873                         x0=-x/2.0;
874                         y=moldyn->dim.y;
875                         y0=-y/2.0;
876                         z=moldyn->dim.z;
877                         z0=-z/2.0;
878                         cr_check=TRUE;
879                         break;
880                 case INS_RECT:
881                         x=iap->x1-iap->x0;
882                         x0=iap->x0;
883                         y=iap->y1-iap->y0;
884                         y0=iap->y0;
885                         z=iap->z1-iap->z0;
886                         z0=iap->z0;
887                         cr_check=TRUE;
888                         break;
889                 case INS_SPHERE:
890                         x=2.0*iap->x1;
891                         x0=iap->x0-iap->x1;
892                         y=x;
893                         y0=iap->y0-iap->x1;
894                         z=x;
895                         z0=iap->z0-iap->x1;
896                         cr_check=TRUE;
897                         break;
898                 case INS_POS:
899                         x0=iap->x0;
900                         y0=iap->y0;
901                         z0=iap->z0;
902                         cr_check=FALSE;
903                         break;
904                 default:
905                         printf("%s unknown insertion mode: %02x\n",
906                                ME,iap->type);
907                         return -1;
908         }
909
910         cnt=0;
911         while(cnt<iap->ins_atoms) {
912                 run=1;
913                 while(run) {
914                         if(iap->type!=INS_POS) {
915                                 r.x=rand_get_double(&(moldyn->random))*x;
916                                 r.y=rand_get_double(&(moldyn->random))*y;
917                                 r.z=rand_get_double(&(moldyn->random))*z;
918                         }
919                         else {
920                                 r.x=0.0;
921                                 r.y=0.0;
922                                 r.z=0.0;
923                         }
924                         r.x+=x0;
925                         r.y+=y0;
926                         r.z+=z0;
927                         // offset
928                         if(iap->type!=INS_TOTAL) {
929                                 r.x+=o;
930                                 r.y+=o;
931                                 r.z+=o;
932                         }
933                         run=0;
934                         if(cr_check==TRUE) {
935                                 dmin=1000;      // for sure too high!
936                                 for(i=0;i<moldyn->count;i++) {
937                                         atom=&(moldyn->atom[i]);
938                                         v3_sub(&dist,&(atom->r),&r);
939                                         check_per_bound(moldyn,&dist);
940                                         d=v3_absolute_square(&dist);
941                                         if(d<(iap->cr*iap->cr)) {
942                                                 run=1;
943                                                 break;
944                                         }
945                                         if(d<dmin)
946                                                 dmin=d;
947                                 }
948                         }
949                         if(iap->type==INS_SPHERE) {
950                                 if((r.x-iap->x0)*(r.x-iap->x0)+
951                                    (r.y-iap->y0)*(r.y-iap->y0)+
952                                    (r.z-iap->z0)*(r.z-iap->z0)>
953                                    (iap->x1*iap->x1)) {
954                                         run=1;
955                                 }
956                         }
957                 }
958                 add_atom(moldyn,iap->element,
959                          iap->brand,iap->attr,&r,&v);
960                 printf("%s atom inserted (%d/%d): %f %f %f\n",
961                        ME,(iap->cnt_steps+1)*iap->ins_atoms,
962                        iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z);
963                 printf("  attributes: ");
964                 if(iap->attr&ATOM_ATTR_VB)
965                         printf("b ");
966                 if(iap->attr&ATOM_ATTR_HB)
967                         printf("h ");
968                 if(iap->attr&ATOM_ATTR_VA)
969                         printf("v ");
970                 if(iap->attr&ATOM_ATTR_FP)
971                         printf("f ");
972                 if(iap->attr&ATOM_ATTR_1BP)
973                         printf("1 ");
974                 if(iap->attr&ATOM_ATTR_2BP)
975                         printf("2 ");
976                 if(iap->attr&ATOM_ATTR_3BP)
977                         printf("3 ");
978                 printf("\n");
979                 printf("  d2 = %f/%f\n",dmin,iap->cr*iap->cr);
980                 cnt+=1;
981         }
982
983         return 0;
984 }
985
986 int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
987
988         t_stage *stage;
989         t_insert_mixed_atoms_params *imp;
990         t_atom *atom;
991         double x,x0,y,y0,z,z0;
992         double dmin,d,cmin,cmax;
993         u8 retry;
994         t_3dvec r,v,dist;
995         int i;
996         
997
998         stage=mdrun->stage.current->data;
999         imp=stage->params;
1000
1001         x=moldyn->dim.x;
1002         x0=-x/2.0;
1003         y=moldyn->dim.y;
1004         y0=-y/2.0;
1005         z=moldyn->dim.z;
1006         z0=-z/2.0;
1007
1008         v.x=0.0;
1009         v.y=0.0;
1010         v.z=0.0;
1011
1012         cmin=imp->crmin*imp->crmin;
1013         cmax=imp->crmax*imp->crmax;
1014
1015         while(imp->amount1|imp->amount2) {
1016                 if(imp->amount1) {
1017                         retry=1;
1018                         while(retry) {
1019                                 retry=0;
1020                                 r.x=rand_get_double(&(moldyn->random))*x;
1021                                 r.y=rand_get_double(&(moldyn->random))*y;
1022                                 r.z=rand_get_double(&(moldyn->random))*z;
1023                                 r.x+=x0;
1024                                 r.y+=y0;
1025                                 r.z+=z0;
1026                                 dmin=1000.0;    // for sure too high!
1027                                 for(i=0;i<moldyn->count;i++) {
1028                                         atom=&(moldyn->atom[i]);
1029                                         v3_sub(&dist,&(atom->r),&r);
1030                                         check_per_bound(moldyn,&dist);
1031                                         d=v3_absolute_square(&dist);
1032                                         if(d<cmin) {
1033                                                 retry=1;
1034                                                 break;
1035                                         }
1036                                         if(d<dmin)
1037                                                 dmin=d;
1038                                 }
1039                                 if(dmin!=1000.0)
1040                                         if(dmin>cmax)
1041                                                 retry=1;
1042                         }
1043                         add_atom(moldyn,imp->element1,
1044                                  imp->brand1,imp->attr1,&r,&v);
1045                         printf("%s (mixed) atom inserted (%d): %f %f %f\n",
1046                                ME,imp->amount1,r.x,r.y,r.z);
1047                         printf("  -> d2 = %f/%f/%f\n",dmin,cmin,cmax);
1048                         imp->amount1-=1;
1049                 }
1050                 if(imp->amount2) {
1051                         retry=1;
1052                         while(retry) {
1053                                 retry=0;
1054                                 r.x=rand_get_double(&(moldyn->random))*x;
1055                                 r.y=rand_get_double(&(moldyn->random))*y;
1056                                 r.z=rand_get_double(&(moldyn->random))*z;
1057                                 r.x+=x0;
1058                                 r.y+=y0;
1059                                 r.z+=z0;
1060                                 dmin=1000.0;    // for sure too high!
1061                                 for(i=0;i<moldyn->count;i++) {
1062                                         atom=&(moldyn->atom[i]);
1063                                         v3_sub(&dist,&(atom->r),&r);
1064                                         check_per_bound(moldyn,&dist);
1065                                         d=v3_absolute_square(&dist);
1066                                         if(d<cmin) {
1067                                                 retry=1;
1068                                                 break;
1069                                         }
1070                                         if(d<dmin)
1071                                                 dmin=d;
1072                                 }
1073                                 if(dmin!=1000.0)
1074                                         if(dmin>cmax)
1075                                                 retry=1;
1076                         }
1077                         add_atom(moldyn,imp->element2,
1078                                  imp->brand2,imp->attr2,&r,&v);
1079                         printf("%s (mixed) atom inserted (%d): %f %f %f\n",
1080                                ME,imp->amount2,r.x,r.y,r.z);
1081                         printf("  -> d2 = %f/%f/%f\n",dmin,cmin,cmax);
1082                         imp->amount2-=1;
1083                 }
1084         }
1085
1086         return 0;
1087 }
1088
1089 int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) {
1090
1091         t_stage *stage;
1092         t_chaattr_params *cap;
1093         t_atom *atom;
1094         int i;
1095
1096         stage=mdrun->stage.current->data;
1097         cap=stage->params;
1098
1099         for(i=0;i<moldyn->count;i++) {
1100                 atom=&(moldyn->atom[i]);
1101                 if(cap->type&CHAATTR_ELEMENT) {
1102                         if(cap->element!=atom->element)
1103                                 continue;
1104                 }
1105                 if(cap->type&CHAATTR_REGION) {
1106                         if(cap->x0>atom->r.x)
1107                                 continue;
1108                         if(cap->y0>atom->r.y)
1109                                 continue;
1110                         if(cap->z0>atom->r.z)
1111                                 continue;
1112                         if(cap->x1<atom->r.x)
1113                                 continue;
1114                         if(cap->y1<atom->r.y)
1115                                 continue;
1116                         if(cap->z1<atom->r.z)
1117                                 continue;
1118                 }
1119                 if(!(cap->type&CHAATTR_TOTALV))
1120                         printf("  changing attributes of atom %d (0x%x)\n",
1121                                i,cap->attr);
1122                 atom->attr=cap->attr;
1123         }
1124
1125         printf("\n\n");
1126
1127         return 0;
1128 }
1129
1130 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
1131
1132         t_stage *stage;
1133         t_chsattr_params *csp;
1134
1135         stage=mdrun->stage.current->data;
1136         csp=stage->params;
1137
1138         if(csp->type&CHSATTR_PCTRL) {
1139                 if(csp->ptau>0)
1140                         set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
1141                 else
1142                         set_p_scale(moldyn,P_SCALE_NONE,1.0);
1143         }
1144         if(csp->type&CHSATTR_TCTRL) {
1145                 if(csp->ttau>0)
1146                         set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
1147                 else
1148                         set_t_scale(moldyn,T_SCALE_NONE,1.0);
1149         }
1150         if(csp->type&CHSATTR_PRELAX) {
1151                 if(csp->dp<0)
1152                         mdrun->sattr&=(~(SATTR_PRELAX));
1153                 else
1154                         mdrun->sattr|=SATTR_PRELAX;
1155                 mdrun->dp=csp->dp;
1156         }
1157         if(csp->type&CHSATTR_TRELAX) {
1158                 if(csp->dt<0)
1159                         mdrun->sattr&=(~(SATTR_TRELAX));
1160                 else
1161                         mdrun->sattr|=SATTR_TRELAX;
1162                 mdrun->dt=csp->dt;
1163         }
1164         if(csp->type&CHSATTR_AVGRST) {
1165                 if(csp->avgrst)
1166                         mdrun->sattr|=SATTR_AVGRST;
1167                 else
1168                         mdrun->sattr&=(~(SATTR_AVGRST));
1169         }
1170         if(csp->type&CHSATTR_RSTEPS) {
1171                 mdrun->relax_steps=csp->rsteps;
1172         }
1173
1174         return 0;
1175 }
1176
1177 #define stage_print(m)  if(!(stage->executed)) \
1178                                 printf("%s",m)
1179
1180 int mdrun_hook(void *ptr1,void *ptr2) {
1181
1182         t_moldyn *moldyn;
1183         t_mdrun *mdrun;
1184         t_stage *stage;
1185         t_list *sl;
1186         int steps;
1187         u8 change_stage;
1188         t_3dvec o;
1189
1190         t_insert_atoms_params *iap;
1191         t_insert_mixed_atoms_params *imp;
1192         t_continue_params *cp;
1193         t_anneal_params *ap;
1194         t_set_temp_params *stp;
1195         t_set_timestep_params *stsp;
1196         t_fill_params *fp;
1197
1198         moldyn=ptr1;
1199         mdrun=ptr2;
1200
1201         sl=&(mdrun->stage);
1202
1203         change_stage=FALSE;
1204
1205         /* return immediately if there are no more stage */
1206         if(sl->current==NULL)
1207                 return 0;
1208
1209         /* get stage description */
1210         stage=sl->current->data;
1211
1212         /* steps in next schedule */
1213         steps=mdrun->relax_steps;
1214
1215         /* check whether relaxation steps are necessary */
1216         if((check_pressure(moldyn,mdrun)==TRUE)&\
1217            (check_temperature(moldyn,mdrun)==TRUE)) {
1218         
1219                 /* be verbose */
1220                 stage_print("\n###########################\n");
1221                 stage_print("# [mdrun] executing stage #\n");
1222                 stage_print("###########################\n\n");
1223                 
1224                 /* stage specific stuff */
1225                 switch(stage->type) {
1226                         case STAGE_DISPLACE_ATOM:
1227                                 stage_print("  -> displace atom\n\n");
1228                                 displace_atom(moldyn,mdrun);
1229                                 change_stage=TRUE;
1230                                 break;
1231                         case STAGE_DEL_ATOMS:
1232                                 stage_print(" -> del atoms\n\n");
1233                                 del_atoms(moldyn,mdrun);
1234                                 change_stage=TRUE;
1235                                 break;
1236                         case STAGE_MODIFY_ATOMS:
1237                                 stage_print(" -> modify atoms\n\n");
1238                                 modify_atoms(moldyn,mdrun);
1239                                 change_stage=TRUE;
1240                                 break;
1241                         case STAGE_INSERT_ATOMS:
1242                                 stage_print("  -> insert atoms\n\n");
1243                                 iap=stage->params;
1244                                 if(iap->cnt_steps==iap->ins_steps) {
1245                                         change_stage=TRUE;
1246                                         break;
1247                                 }
1248                                 insert_atoms(moldyn,mdrun);
1249                                 iap->cnt_steps+=1;
1250                                 break;
1251
1252
1253
1254                         case STAGE_INSERT_MIXED_ATOMS:
1255                                 stage_print("  -> insert mixed atoms\n\n");
1256                                 imp=stage->params;
1257                                 insert_mixed_atoms(moldyn,mdrun);
1258                                 change_stage=TRUE;
1259                                 break;
1260
1261
1262
1263                         case STAGE_CONTINUE:
1264                                 stage_print("  -> continue\n\n");
1265                                 if(stage->executed==TRUE) {
1266                                         change_stage=TRUE;
1267                                         break;
1268                                 }
1269                                 cp=stage->params;
1270                                 steps=cp->runs;
1271                                 break;
1272                         case STAGE_ANNEAL:
1273                                 stage_print("  -> anneal\n\n");
1274                                 ap=stage->params;
1275                                 if(ap->count==ap->runs) {
1276                                         change_stage=TRUE;
1277                                         break;
1278                                 }
1279                                 if(moldyn->t_ref+ap->dt>=0.0)
1280                                         set_temperature(moldyn,
1281                                                         moldyn->t_ref+ap->dt);
1282                                 ap->count+=1;
1283                                 steps=ap->interval;
1284                                 break;
1285                         case STAGE_CHAATTR:
1286                                 stage_print("  -> change atom attributes\n\n");
1287                                 chaatr(moldyn,mdrun);
1288                                 change_stage=TRUE;
1289                                 break;
1290                         case STAGE_CHSATTR:
1291                                 stage_print("  -> change sys attributes\n\n");
1292                                 chsattr(moldyn,mdrun);
1293                                 change_stage=TRUE;
1294                                 break;
1295                         case STAGE_SET_TEMP:
1296                                 stage_print("  -> set temperature\n\n");
1297                                 stp=stage->params;
1298                                 if(stp->type&SET_TEMP_CURRENT) {
1299                                         set_temperature(moldyn,moldyn->t_avg);
1300                                 }
1301                                 else {
1302                                         set_temperature(moldyn,stp->val);
1303                                 }
1304                                 change_stage=TRUE;
1305                                 break;
1306                         case STAGE_SET_TIMESTEP:
1307                                 stage_print("  -> set timestep\n\n");
1308                                 stsp=stage->params;
1309                                 mdrun->timestep=stsp->tau;
1310                                 change_stage=TRUE;
1311                                 break;
1312                         case STAGE_FILL:
1313                                 stage_print("  -> fill lattice\n\n");
1314                                 fp=stage->params;
1315                                 switch(fp->lattice) {
1316                                         case ZINCBLENDE:
1317
1318                                         o.x=0.5*0.25*fp->lc;
1319                                         o.y=o.x;
1320                                         o.z=o.x;
1321                                         create_lattice(moldyn,
1322                                                        FCC,fp->lc,
1323                                                        mdrun->element1,
1324                                                        DEFAULT_ATOM_ATTR,0,
1325                                                        fp->lx,fp->ly,fp->lz,
1326                                                        &o,
1327                                                        &(fp->p_params),
1328                                                        &(fp->d_params),
1329                                                        &(fp->o_params));
1330                                         o.x+=0.25*fp->lc;
1331                                         o.y=o.x;
1332                                         o.z=o.x;
1333                                         create_lattice(moldyn,
1334                                                        FCC,fp->lc,
1335                                                        mdrun->element2,
1336                                                        DEFAULT_ATOM_ATTR,1,
1337                                                        fp->lx,fp->ly,fp->lz,
1338                                                        &o,
1339                                                        &(fp->p_params),
1340                                                        &(fp->d_params),
1341                                                        &(fp->o_params));
1342                                         break;
1343
1344                                         default:
1345
1346                                         create_lattice(moldyn,
1347                                                        fp->lattice,fp->lc,
1348                                                        fp->fill_element,
1349                                                        DEFAULT_ATOM_ATTR,
1350                                                        fp->fill_brand,
1351                                                        fp->lx,fp->ly,fp->lz,
1352                                                        NULL,
1353                                                        &(fp->p_params),
1354                                                        &(fp->d_params),
1355                                                        &(fp->o_params));
1356                                         break;
1357                                 }
1358                                 moldyn_bc_check(moldyn);
1359                                 change_stage=TRUE;
1360                                 break;
1361                         case STAGE_THERMAL_INIT:
1362                                 stage_print("  -> thermal init\n\n");
1363                                 thermal_init(moldyn,TRUE);
1364                                 change_stage=TRUE;
1365                                 break;
1366                         default:
1367                                 printf("%s unknwon stage type\n",ME);
1368                                 break;
1369                 }
1370         
1371                 /* mark as executed */
1372                 stage->executed=TRUE;
1373         
1374                 /* change state */
1375                 if(change_stage==TRUE) {
1376                         printf("%s finished stage\n",ME);
1377                         if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
1378                                 printf("%s no more stages\n",ME);
1379                                 return 0;
1380                         }
1381                         steps=0;
1382                 }
1383
1384         }
1385         else {
1386
1387                 /* averages */
1388                 if(mdrun->sattr&SATTR_AVGRST)
1389                         average_reset(moldyn);
1390
1391         }
1392
1393         /* continue simulation */
1394         moldyn_add_schedule(moldyn,steps,mdrun->timestep);
1395
1396         return 0;
1397 }
1398
1399 int main(int argc,char **argv) {
1400
1401         t_mdrun mdrun;
1402         t_moldyn moldyn;
1403         //t_3dvec o;
1404
1405         /* clear structs */
1406         memset(&mdrun,0,sizeof(t_mdrun));
1407         memset(&moldyn,0,sizeof(t_moldyn));
1408
1409         /* parse arguments */
1410         if(mdrun_parse_argv(&mdrun,argc,argv)<0)
1411                 return -1;
1412
1413         /* initialize list system */
1414         list_init_f(&(mdrun.stage));
1415
1416         /* parse config file */
1417         mdrun_parse_config(&mdrun);
1418
1419         /* reset the stage list */
1420         list_reset_f(&(mdrun.stage));
1421
1422         /* sanity check! */
1423
1424         /* prepare simulation */
1425
1426         if(mdrun.continue_file[0]) {
1427                 // read the save file
1428                 moldyn_read_save_file(&moldyn,mdrun.continue_file);
1429                 // manualaadjusting some stuff
1430                 moldyn.argc=argc;
1431                 moldyn.args=argv;
1432                 rand_init(&(moldyn.random),NULL,1);
1433                 moldyn.random.status|=RAND_STAT_VERBOSE;
1434         }
1435         else {
1436                 moldyn_init(&moldyn,argc,argv);
1437         }
1438         
1439         if(set_int_alg(&moldyn,mdrun.intalgo)<0)
1440                 return -1;
1441
1442         /* potential */
1443         set_cutoff(&moldyn,mdrun.cutoff);
1444         if(set_potential(&moldyn,mdrun.potential)<0)
1445                 return -1;
1446         switch(mdrun.potential) {
1447                 case MOLDYN_POTENTIAL_AM:
1448                         albe_mult_set_params(&moldyn,
1449                                              mdrun.element1,
1450                                              mdrun.element2);
1451                         break;
1452                 case MOLDYN_POTENTIAL_TM:
1453                         tersoff_mult_set_params(&moldyn,
1454                                                 mdrun.element1,
1455                                                 mdrun.element2);
1456                         break;
1457                 case MOLDYN_POTENTIAL_HO:
1458                         harmonic_oscillator_set_params(&moldyn,mdrun.element1);
1459                         break;
1460                 case MOLDYN_POTENTIAL_LJ:
1461                         lennard_jones_set_params(&moldyn,mdrun.element1);
1462                         break;
1463                 default:
1464                         printf("%s unknown potential: %02x\n",
1465                                ME,mdrun.potential);
1466                         return -1;
1467         }
1468
1469         /* if it is a continue run, reset schedule and skip lattice init */
1470         if(mdrun.continue_file[0]) {
1471                 memset(&(moldyn.schedule),0,sizeof(t_moldyn_schedule));
1472                 goto addsched;
1473         }
1474
1475         /* initial lattice and dimensions */
1476         set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
1477         set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
1478         /* replaced by fill stage !! 
1479         switch(mdrun.lattice) {
1480                 case FCC:
1481                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.fill_element,
1482                                        mdrun.m1,DEFAULT_ATOM_ATTR,
1483                                        mdrun.fill_brand,mdrun.lx,
1484                                        mdrun.ly,mdrun.lz,NULL,0,NULL);
1485                         break;
1486                 case DIAMOND:
1487                         create_lattice(&moldyn,DIAMOND,mdrun.lc,
1488                                        mdrun.fill_element,
1489                                        mdrun.m1,DEFAULT_ATOM_ATTR,
1490                                        mdrun.fill_brand,mdrun.lx,
1491                                        mdrun.ly,mdrun.lz,NULL,0,NULL);
1492                         break;
1493                 case ZINCBLENDE:
1494                         o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
1495                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
1496                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
1497                                        mdrun.ly,mdrun.lz,&o,0,NULL);
1498                         o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
1499                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
1500                                        mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx,
1501                                        mdrun.ly,mdrun.lz,&o,0,NULL);
1502                         break;
1503                 case NONE:
1504                         set_nn_dist(&moldyn,mdrun.nnd);
1505                         break;
1506                 default:
1507                         printf("%s unknown lattice type: %02x\n",
1508                                ME,mdrun.lattice);
1509                         return -1;
1510         }
1511         moldyn_bc_check(&moldyn);
1512         replaced by fill stage !! */
1513
1514         /* temperature and pressure */
1515         set_temperature(&moldyn,mdrun.temperature);
1516         set_pressure(&moldyn,mdrun.pressure);
1517         /* replaced thermal_init stage
1518         thermal_init(&moldyn,TRUE);
1519         */
1520
1521 addsched:
1522         /* first schedule */
1523         moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
1524
1525         /* log */
1526         moldyn_set_log_dir(&moldyn,mdrun.sdir);
1527         moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
1528         if(mdrun.elog)
1529                 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
1530         if(mdrun.tlog)
1531                 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
1532         if(mdrun.plog)
1533                 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
1534         if(mdrun.vlog)
1535                 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
1536         if(mdrun.visualize)
1537                 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
1538         if(mdrun.save)
1539                 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
1540         moldyn_set_log(&moldyn,CREATE_REPORT,0);
1541         set_avg_skip(&moldyn,mdrun.avgskip);
1542
1543         /* prepare the hook function */
1544         moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
1545
1546         /* run the simulation */
1547         moldyn_integrate(&moldyn);
1548
1549         /* shutdown */
1550         moldyn_shutdown(&moldyn);
1551         del_stages(&mdrun);
1552         list_destroy_f(&(mdrun.stage));
1553
1554         return 0;
1555 }