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                                         case 'n':
 426                                                 cap.type|=CHAATTR_NUMBER;
 427                                         default:
 428                                                 break;
 429                                 }
 430                         }
 431                         i=2;
 432                         if(cap.type&CHAATTR_REGION) {
 433                                 cap.x0=atof(word[2]);
 434                                 cap.y0=atof(word[3]);
 435                                 cap.z0=atof(word[4]);
 436                                 cap.x1=atof(word[5]);
 437                                 cap.y1=atof(word[6]);
 438                                 cap.z1=atof(word[7]);
 439                                 i+=6;
 440                         }
 441                         if(cap.type&CHAATTR_ELEMENT) {
 442                                 cap.element=atoi(word[i]);
 443                                 i+=1;
 444                         }
 445                         if(cap.type&CHAATTR_NUMBER) {
 446                                 cap.element=atoi(word[i]);
 447                                 i+=1;
 448                         }
 449                         for(o=0;o<strlen(word[i]);o++) {
 450                                 switch(word[i][o]) {
 451                                         case 'b':
 452                                                 cap.attr|=ATOM_ATTR_VB;
 453                                                 break;
 454                                         case 'h':
 455                                                 cap.attr|=ATOM_ATTR_HB;
 456                                                 break;
 457                                         case 'v':
 458                                                 cap.attr|=ATOM_ATTR_VA;
 459                                                 break;
 460                                         case 'f':
 461                                                 cap.attr|=ATOM_ATTR_FP;
 462                                                 break;
 463                                         case '1':
 464                                                 cap.attr|=ATOM_ATTR_1BP;
 465                                                 break;
 466                                         case '2':
 467                                                 cap.attr|=ATOM_ATTR_2BP;
 468                                                 break;
 469                                         case '3':
 470                                                 cap.attr|=ATOM_ATTR_3BP;
 471                                                 break;
 472                                         default:
 473                                                 break;
 474                                 }
 475                         }
 476                         add_stage(mdrun,STAGE_CHAATTR,&cap);
 477                 }
 478                 else if(!strncmp(word[0],"sattr",5)) {
 479                         // for satrib line we need a special stage
 480                         // containing one schedule of 0 loops ...
 481                         csp.type=0;
 482                         for(i=1;i<wcnt;i++) {
 483                                 if(!strncmp(word[i],"pctrl",5)) {
 484                                         csp.ptau=atof(word[++i]);
 485                                         if(csp.ptau>0)
 486                                                 csp.ptau=0.01/(csp.ptau*GPA);
 487                                         csp.type|=CHSATTR_PCTRL;
 488                                 }
 489                                 if(!strncmp(word[i],"tctrl",5)) {
 490                                         csp.ttau=atof(word[++i]);
 491                                         csp.type|=CHSATTR_TCTRL;
 492                                 }
 493                                 if(!strncmp(word[i],"prelax",6)) {
 494                                         csp.dp=atof(word[++i])*BAR;
 495                                         csp.type|=CHSATTR_PRELAX;
 496                                 }
 497                                 if(!strncmp(word[i],"trelax",6)) {
 498                                         csp.dt=atof(word[++i]);
 499                                         csp.type|=CHSATTR_TRELAX;
 500                                 }
 501                                 if(!strncmp(word[i],"rsteps",6)) {
 502                                         csp.rsteps=atoi(word[++i]);
 503                                         csp.type|=CHSATTR_RSTEPS;
 504                                 }
 505                                 if(!strncmp(word[i],"avgrst",6)) {
 506                                         csp.avgrst=atoi(word[++i]);
 507                                         csp.type|=CHSATTR_AVGRST;
 508                                 }
 509                         }
 510                         add_stage(mdrun,STAGE_CHSATTR,&csp);
 511                 }
 512                 else if(!strncmp(word[0],"prerun",6))
 513                         mdrun->prerun=atoi(word[1]);
 514                 else if(!strncmp(word[0],"avgskip",7))
 515                         mdrun->avgskip=atoi(word[1]);
 516                 else if(!strncmp(word[0],"elog",4))
 517                         mdrun->elog=atoi(word[1]);
 518                 else if(!strncmp(word[0],"tlog",4))
 519                         mdrun->tlog=atoi(word[1]);
 520                 else if(!strncmp(word[0],"plog",4))
 521                         mdrun->plog=atoi(word[1]);
 522                 else if(!strncmp(word[0],"vlog",4))
 523                         mdrun->vlog=atoi(word[1]);
 524                 else if(!strncmp(word[0],"save",4))
 525                         mdrun->save=atoi(word[1]);
 526                 else if(!strncmp(word[0],"visualize",9))
 527                         mdrun->visualize=atoi(word[1]);
 528                 else if(!strncmp(word[0],"stage",5)) {
 529                         // for every stage line, add a stage
 530                         if(!strncmp(word[1],"displace",8)) {
 531                                 dap.nr=atoi(word[2]);
 532                                 dap.dx=atof(word[3]);
 533                                 dap.dy=atof(word[4]);
 534                                 dap.dz=atof(word[5]);
 535                                 add_stage(mdrun,STAGE_DISPLACE_ATOM,&dap);
 536                         }
 537                         else if(!strncmp(word[1],"del_atoms",9)) {
 538                                 delp.o.x=atof(word[2]);
 539                                 delp.o.y=atof(word[3]);
 540                                 delp.o.z=atof(word[4]);
 541                                 delp.r=atof(word[5]);
 542                                 add_stage(mdrun,STAGE_DEL_ATOMS,&delp);
 543                         }
 544                         else if(!strncmp(word[1],"mod_atoms",8)) {
 545                                 i=2;
 546                                 while(i<wcnt) {
 547                                         if(!strncmp(word[i],"t",1)) {
 548                                                 map.tag=atoi(word[++i]);
 549                                                 i+=1;
 550                                         }
 551                                         if(!strncmp(word[i],"ekin",5)) {
 552                                                 map.ekin.x=atof(word[++i])*EV;
 553                                                 map.ekin.y=atof(word[++i])*EV;
 554                                                 map.ekin.z=atof(word[++i])*EV;
 555                                                 i+=1;
 556                                         }
 557                                 }
 558                                 add_stage(mdrun,STAGE_MODIFY_ATOMS,&map);
 559                         }
 560                         else if(!strncmp(word[1],"ins_atoms",9)) {
 561                                 iap.ins_steps=atoi(word[2]);
 562                                 iap.ins_atoms=atoi(word[3]);
 563                                 iap.element=atoi(word[4]);
 564                                 iap.brand=atoi(word[5]);
 565                                 for(i=0;i<strlen(word[6]);i++) {
 566                                         switch(word[6][i]) {
 567                                                 case 'b':
 568                                                         iap.attr|=ATOM_ATTR_VB;
 569                                                         break;
 570                                                 case 'h':
 571                                                         iap.attr|=ATOM_ATTR_HB;
 572                                                         break;
 573                                                 case 'v':
 574                                                         iap.attr|=ATOM_ATTR_VA;
 575                                                         break;
 576                                                 case 'f':
 577                                                         iap.attr|=ATOM_ATTR_FP;
 578                                                         break;
 579                                                 case '1':
 580                                                         iap.attr|=ATOM_ATTR_1BP;
 581                                                         break;
 582                                                 case '2':
 583                                                         iap.attr|=ATOM_ATTR_2BP;
 584                                                         break;
 585                                                 case '3':
 586                                                         iap.attr|=ATOM_ATTR_3BP;
 587                                                         break;
 588                                                 default:
 589                                                         break;
 590                                         }
 591                                 }
 592                                 switch(word[7][0]) {
 593                                         // insertion types
 594                                         case 'p':
 595                                                 iap.type=INS_POS;
 596                                                 iap.x0=atof(word[8]);
 597                                                 iap.y0=atof(word[9]);
 598                                                 iap.z0=atof(word[10]);
 599                                                 break;
 600                                         case 'P':
 601                                                 iap.type=INS_RELPOS;
 602                                                 iap.x0=atof(word[8]);
 603                                                 iap.y0=atof(word[9]);
 604                                                 iap.z0=atof(word[10]);
 605                                                 break;
 606                                         case 'r':
 607                                                 switch(word[8][0]) {
 608
 609                                                 case 't':
 610                                                         iap.type=INS_TOTAL;
 611                                                         iap.cr=atof(word[9]);
 612                                                         break;
 613                                                 case 'r':
 614                                                         iap.type=INS_RECT;
 615                                                         iap.x0=atof(word[9]);
 616                                                         iap.y0=atof(word[10]);
 617                                                         iap.z0=atof(word[11]);
 618                                                         iap.x1=atof(word[12]);
 619                                                         iap.y1=atof(word[13]);
 620                                                         iap.z1=atof(word[14]);
 621                                                         iap.cr=atof(word[15]);
 622                                                         break;
 623                                                 case 's':
 624                                                         iap.type=INS_SPHERE;
 625                                                         iap.x0=atof(word[9]);
 626                                                         iap.y0=atof(word[10]);
 627                                                         iap.z0=atof(word[11]);
 628                                                         iap.x1=atof(word[12]);
 629                                                         iap.cr=atof(word[13]);
 630                                                         break;
 631                                                 default:
 632                                                         break;
 633                                                 }
 634                                         default:
 635                                                 break;
 636                                 }
 637                                 add_stage(mdrun,STAGE_INSERT_ATOMS,&iap);
 638                         }
 639
 640
 641                         // HERE WE GO ...
 642
 643                         else if(!strncmp(word[1],"ins_m_atoms",11)) {
 644                                 imp.element1=atoi(word[2]);
 645                                 imp.element2=atoi(word[3]);
 646                                 imp.amount1=atoi(word[4]);
 647                                 imp.amount2=atoi(word[5]);
 648                                 imp.brand1=atoi(word[6]);
 649                                 imp.brand2=atoi(word[7]);
 650                                 imp.crmin=atof(word[8]);
 651                                 imp.crmax=atof(word[9]);
 652                                 /* do this later ...
 653                                 for(i=0;i<strlen(word[8]);i++) {
 654                                         switch(word[8][i]) {
 655                                                 case 'b':
 656                                                         imp.attr|=ATOM_ATTR_VB;
 657                                                         break;
 658                                                 case 'h':
 659                                                         imp.attr|=ATOM_ATTR_HB;
 660                                                         break;
 661                                                 case 'v':
 662                                                         imp.attr|=ATOM_ATTR_VA;
 663                                                         break;
 664                                                 case 'f':
 665                                                         imp.attr|=ATOM_ATTR_FP;
 666                                                         break;
 667                                                 case '1':
 668                                                         imp.attr|=ATOM_ATTR_1BP;
 669                                                         break;
 670                                                 case '2':
 671                                                         imp.attr|=ATOM_ATTR_2BP;
 672                                                         break;
 673                                                 case '3':
 674                                                         imp.attr|=ATOM_ATTR_3BP;
 675                                                         break;
 676                                                 default:
 677                                                         break;
 678                                         }
 679                                 }
 680                                 */
 681                                 imp.attr1=ATOM_ATTR_HB|ATOM_ATTR_VA|ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_FP;
 682                                 imp.attr2=ATOM_ATTR_HB|ATOM_ATTR_VA|ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_FP;
 683                                 add_stage(mdrun,STAGE_INSERT_MIXED_ATOMS,&imp);
 684                         }
 685
 686
 687
 688
 689
 690                         else if(!strncmp(word[1],"continue",8)) {
 691                                 cp.runs=atoi(word[2]);
 692                                 add_stage(mdrun,STAGE_CONTINUE,&cp);
 693                         }
 694                         else if(!strncmp(word[1],"anneal",6)) {
 695                                 ap.count=0;
 696                                 ap.runs=atoi(word[2]);
 697                                 ap.dt=atof(word[3]);
 698                                 ap.interval=atoi(word[4]);
 699                                 add_stage(mdrun,STAGE_ANNEAL,&ap);
 700                         }
 701                         else if(!strncmp(word[1],"set_temp",8)) {
 702                                 if(word[2][0]=='c') {
 703                                         stp.type=SET_TEMP_CURRENT;
 704                                         stp.val=0.0;
 705                                 }
 706                                 else {
 707                                         stp.type=SET_TEMP_VALUE;
 708                                         stp.val=atof(word[2]);
 709                                 }
 710                                 add_stage(mdrun,STAGE_SET_TEMP,&stp);
 711                         }
 712                         else if(!strncmp(word[1],"set_timestep",12)) {
 713                                 stsp.tau=atof(word[2]);
 714                                 add_stage(mdrun,STAGE_SET_TIMESTEP,&stsp);
 715                         }
 716                         else {
 717                                 printf("%s unknown stage type: %s\n",
 718                                        ME,word[1]);
 719                                 return -1;
 720                         }
 721                 }
 722                 else {
 723                         printf("%s unknown keyword '%s', skipped!\n",
 724                                ME,word[0]);
 725                         continue;
 726                 }
 727         }
 728
 729         /* close file */
 730         close(fd);
 731
 732         return 0;
 733 }
 734
 735 int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
 736
 737         double delta;
 738
 739         if(!(mdrun->sattr&SATTR_PRELAX)) {
 740                 return TRUE;
 741         }
 742
 743         delta=moldyn->p_avg-moldyn->p_ref;
 744
 745         if(delta<0)
 746                 delta=-delta;
 747
 748         if(delta>mdrun->dp)
 749                 return FALSE;
 750
 751         return TRUE;
 752 }
 753
 754 int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
 755
 756         double delta;
 757
 758         if(!(mdrun->sattr&SATTR_TRELAX))
 759                 return TRUE;
 760
 761         delta=moldyn->t_avg-moldyn->t_ref;
 762
 763         if(delta<0)
 764                 delta=-delta;
 765
 766         if(delta>mdrun->dt)
 767                 return FALSE;
 768
 769         return TRUE;
 770 }
 771
 772 int displace_atom(t_moldyn *moldyn,t_mdrun *mdrun) {
 773
 774         t_displace_atom_params *dap;
 775         t_stage *stage;
 776         t_atom *atom;
 777
 778         stage=mdrun->stage.current->data;
 779         dap=stage->params;
 780
 781         atom=&(moldyn->atom[dap->nr]);
 782         atom->r.x+=dap->dx;
 783         atom->r.y+=dap->dy;
 784         atom->r.z+=dap->dz;
 785
 786         return 0;
 787 }
 788
 789 int del_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 790
 791         t_stage *stage;
 792         t_del_atoms_params *delp;
 793         int i;
 794         t_3dvec dist;
 795
 796         stage=mdrun->stage.current->data;
 797         delp=stage->params;
 798
 799         for(i=0;i<moldyn->count;i++) {
 800                 v3_sub(&dist,&(delp->o),&(moldyn->atom[i].r));
 801 //printf("%d ----> %f %f %f = %f | %f\n",i,dist.x,dist.y,dist.z,v3_absolute_square(&dist),delp->r*delp->r);
 802                 if(v3_absolute_square(&dist)<=(delp->r*delp->r)) {
 803                         del_atom(moldyn,moldyn->atom[i].tag);
 804                         printf("%s atom deleted: %d %d %d\n",ME,
 805                                moldyn->atom[i].tag,moldyn->atom[i].element,
 806                                moldyn->atom[i].brand);
 807                 }
 808         }
 809
 810         return 0;
 811 }
 812
 813 int modify_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 814
 815         t_modify_atoms_params *map;
 816         t_stage *stage;
 817         t_atom *atom;
 818         t_3dvec v;
 819         int i;
 820
 821         atom=moldyn->atom;
 822         stage=mdrun->stage.current->data;
 823         map=stage->params;
 824         v.x=0.0; v.y=0.0; v.z=0.0;
 825
 826         for(i=0;i<moldyn->count;i++) {
 827                 if(atom[i].tag==map->tag) {
 828                         v.x=sqrt(2.0*fabs(map->ekin.x)/atom[i].mass);
 829                         if(map->ekin.x<0.0)
 830                                 v.x=-v.x;
 831                         v.y=sqrt(2.0*fabs(map->ekin.y)/atom[i].mass);
 832                         if(map->ekin.y<0.0)
 833                                 v.y=-v.y;
 834                         v.z=sqrt(2.0*fabs(map->ekin.z)/atom[i].mass);
 835                         if(map->ekin.z<0.0)
 836                                 v.z=-v.z;
 837                         v3_copy(&(atom[i].v),&v);
 838                         printf("%s atom modified: v = (%f %f %f)\n",
 839                                ME,v.x,v.y,v.z);
 840                 }
 841         }
 842
 843         return 0;
 844 }
 845
 846 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 847
 848         t_insert_atoms_params *iap;
 849         t_stage *stage;
 850         t_atom *atom;
 851         t_3dvec r,v,dist;
 852         double d,dmin,o;
 853         int cnt,i;
 854         double x,y,z;
 855         double x0,y0,z0;
 856         u8 cr_check,run;
 857
 858         stage=mdrun->stage.current->data;
 859         iap=stage->params;
 860
 861         cr_check=FALSE;
 862
 863         v.x=0.0; v.y=0.0; v.z=0.0;
 864         x=0; y=0; z=0;
 865         o=0; dmin=0;
 866
 867         switch(mdrun->lattice) {
 868                 case CUBIC:
 869                         o=0.5*mdrun->lc;
 870                         break;
 871                 case FCC:
 872                         o=0.25*mdrun->lc;
 873                         break;
 874                 case DIAMOND:
 875                 case ZINCBLENDE:
 876                         o=0.125*mdrun->lc;
 877                         break;
 878                 default:
 879                         break;
 880         }
 881
 882         switch(iap->type) {
 883                 case INS_TOTAL:
 884                         x=moldyn->dim.x;
 885                         x0=-x/2.0;
 886                         y=moldyn->dim.y;
 887                         y0=-y/2.0;
 888                         z=moldyn->dim.z;
 889                         z0=-z/2.0;
 890                         cr_check=TRUE;
 891                         break;
 892                 case INS_RECT:
 893                         x=iap->x1-iap->x0;
 894                         x0=iap->x0;
 895                         y=iap->y1-iap->y0;
 896                         y0=iap->y0;
 897                         z=iap->z1-iap->z0;
 898                         z0=iap->z0;
 899                         cr_check=TRUE;
 900                         break;
 901                 case INS_SPHERE:
 902                         x=2.0*iap->x1;
 903                         x0=iap->x0-iap->x1;
 904                         y=x;
 905                         y0=iap->y0-iap->x1;
 906                         z=x;
 907                         z0=iap->z0-iap->x1;
 908                         cr_check=TRUE;
 909                         break;
 910                 case INS_POS:
 911                 case INS_RELPOS:
 912                         x0=iap->x0;
 913                         y0=iap->y0;
 914                         z0=iap->z0;
 915                         cr_check=FALSE;
 916                         break;
 917                 default:
 918                         printf("%s unknown insertion mode: %02x\n",
 919                                ME,iap->type);
 920                         return -1;
 921         }
 922
 923         cnt=0;
 924         while(cnt<iap->ins_atoms) {
 925                 run=1;
 926                 while(run) {
 927                         if((iap->type!=INS_POS)&&(iap->type!=INS_RELPOS)) {
 928                                 r.x=rand_get_double(&(moldyn->random))*x;
 929                                 r.y=rand_get_double(&(moldyn->random))*y;
 930                                 r.z=rand_get_double(&(moldyn->random))*z;
 931                         }
 932                         else {
 933                                 r.x=0.0;
 934                                 r.y=0.0;
 935                                 r.z=0.0;
 936                         }
 937                         if(iap->type==INS_RELPOS) {
 938                                 r.x+=x0*mdrun->lc;
 939                                 r.y+=y0*mdrun->lc;
 940                                 r.z+=z0*mdrun->lc;
 941                         }
 942                         else {
 943                                 r.x+=x0;
 944                                 r.y+=y0;
 945                                 r.z+=z0;
 946                         }
 947                         // offset
 948                         if(iap->type!=INS_TOTAL) {
 949                                 r.x+=o;
 950                                 r.y+=o;
 951                                 r.z+=o;
 952                         }
 953                         run=0;
 954                         if(cr_check==TRUE) {
 955                                 dmin=1000;      // for sure too high!
 956                                 for(i=0;i<moldyn->count;i++) {
 957                                         atom=&(moldyn->atom[i]);
 958                                         v3_sub(&dist,&(atom->r),&r);
 959                                         check_per_bound(moldyn,&dist);
 960                                         d=v3_absolute_square(&dist);
 961                                         if(d<(iap->cr*iap->cr)) {
 962                                                 run=1;
 963                                                 break;
 964                                         }
 965                                         if(d<dmin)
 966                                                 dmin=d;
 967                                 }
 968                         }
 969                         if(iap->type==INS_SPHERE) {
 970                                 if((r.x-iap->x0)*(r.x-iap->x0)+
 971                                    (r.y-iap->y0)*(r.y-iap->y0)+
 972                                    (r.z-iap->z0)*(r.z-iap->z0)>
 973                                    (iap->x1*iap->x1)) {
 974                                         run=1;
 975                                 }
 976                         }
 977                 }
 978                 add_atom(moldyn,iap->element,
 979                          iap->brand,iap->attr,&r,&v);
 980                 printf("%s atom inserted (%d/%d): %f %f %f\n",
 981                        ME,(iap->cnt_steps+1)*iap->ins_atoms,
 982                        iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z);
 983                 printf("  attributes: ");
 984                 if(iap->attr&ATOM_ATTR_VB)
 985                         printf("b ");
 986                 if(iap->attr&ATOM_ATTR_HB)
 987                         printf("h ");
 988                 if(iap->attr&ATOM_ATTR_VA)
 989                         printf("v ");
 990                 if(iap->attr&ATOM_ATTR_FP)
 991                         printf("f ");
 992                 if(iap->attr&ATOM_ATTR_1BP)
 993                         printf("1 ");
 994                 if(iap->attr&ATOM_ATTR_2BP)
 995                         printf("2 ");
 996                 if(iap->attr&ATOM_ATTR_3BP)
 997                         printf("3 ");
 998                 printf("\n");
 999                 printf("  d2 = %f/%f\n",dmin,iap->cr*iap->cr);
1000                 cnt+=1;
1001         }
1002
1003         return 0;
1004 }
1005
1006 int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
1007
1008         t_stage *stage;
1009         t_insert_mixed_atoms_params *imp;
1010         t_atom *atom;
1011         double x,x0,y,y0,z,z0;
1012         double dmin,d,cmin,cmax;
1013         u8 retry;
1014         t_3dvec r,v,dist;
1015         int i;
1016
1017
1018         stage=mdrun->stage.current->data;
1019         imp=stage->params;
1020
1021         x=moldyn->dim.x;
1022         x0=-x/2.0;
1023         y=moldyn->dim.y;
1024         y0=-y/2.0;
1025         z=moldyn->dim.z;
1026         z0=-z/2.0;
1027
1028         v.x=0.0;
1029         v.y=0.0;
1030         v.z=0.0;
1031
1032         cmin=imp->crmin*imp->crmin;
1033         cmax=imp->crmax*imp->crmax;
1034
1035         while(imp->amount1|imp->amount2) {
1036                 if(imp->amount1) {
1037                         retry=1;
1038                         while(retry) {
1039                                 retry=0;
1040                                 r.x=rand_get_double(&(moldyn->random))*x;
1041                                 r.y=rand_get_double(&(moldyn->random))*y;
1042                                 r.z=rand_get_double(&(moldyn->random))*z;
1043                                 r.x+=x0;
1044                                 r.y+=y0;
1045                                 r.z+=z0;
1046                                 dmin=1000.0;    // for sure too high!
1047                                 for(i=0;i<moldyn->count;i++) {
1048                                         atom=&(moldyn->atom[i]);
1049                                         v3_sub(&dist,&(atom->r),&r);
1050                                         check_per_bound(moldyn,&dist);
1051                                         d=v3_absolute_square(&dist);
1052                                         if(d<cmin) {
1053                                                 retry=1;
1054                                                 break;
1055                                         }
1056                                         if(d<dmin)
1057                                                 dmin=d;
1058                                 }
1059                                 if(dmin!=1000.0)
1060                                         if(dmin>cmax)
1061                                                 retry=1;
1062                         }
1063                         add_atom(moldyn,imp->element1,
1064                                  imp->brand1,imp->attr1,&r,&v);
1065                         printf("%s (mixed) atom inserted (%d): %f %f %f\n",
1066                                ME,imp->amount1,r.x,r.y,r.z);
1067                         printf("  -> d2 = %f/%f/%f\n",dmin,cmin,cmax);
1068                         imp->amount1-=1;
1069                 }
1070                 if(imp->amount2) {
1071                         retry=1;
1072                         while(retry) {
1073                                 retry=0;
1074                                 r.x=rand_get_double(&(moldyn->random))*x;
1075                                 r.y=rand_get_double(&(moldyn->random))*y;
1076                                 r.z=rand_get_double(&(moldyn->random))*z;
1077                                 r.x+=x0;
1078                                 r.y+=y0;
1079                                 r.z+=z0;
1080                                 dmin=1000.0;    // for sure too high!
1081                                 for(i=0;i<moldyn->count;i++) {
1082                                         atom=&(moldyn->atom[i]);
1083                                         v3_sub(&dist,&(atom->r),&r);
1084                                         check_per_bound(moldyn,&dist);
1085                                         d=v3_absolute_square(&dist);
1086                                         if(d<cmin) {
1087                                                 retry=1;
1088                                                 break;
1089                                         }
1090                                         if(d<dmin)
1091                                                 dmin=d;
1092                                 }
1093                                 if(dmin!=1000.0)
1094                                         if(dmin>cmax)
1095                                                 retry=1;
1096                         }
1097                         add_atom(moldyn,imp->element2,
1098                                  imp->brand2,imp->attr2,&r,&v);
1099                         printf("%s (mixed) atom inserted (%d): %f %f %f\n",
1100                                ME,imp->amount2,r.x,r.y,r.z);
1101                         printf("  -> d2 = %f/%f/%f\n",dmin,cmin,cmax);
1102                         imp->amount2-=1;
1103                 }
1104         }
1105
1106         return 0;
1107 }
1108
1109 int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) {
1110
1111         t_stage *stage;
1112         t_chaattr_params *cap;
1113         t_atom *atom;
1114         int i;
1115
1116         stage=mdrun->stage.current->data;
1117         cap=stage->params;
1118
1119         for(i=0;i<moldyn->count;i++) {
1120                 atom=&(moldyn->atom[i]);
1121                 if(cap->type&CHAATTR_ELEMENT) {
1122                         if(cap->element!=atom->element)
1123                                 continue;
1124                 }
1125                 if(cap->type&CHAATTR_NUMBER) {
1126                         if(cap->element!=atom->tag)
1127                                 continue;
1128                 }
1129                 if(cap->type&CHAATTR_REGION) {
1130                         if(cap->x0>atom->r.x)
1131                                 continue;
1132                         if(cap->y0>atom->r.y)
1133                                 continue;
1134                         if(cap->z0>atom->r.z)
1135                                 continue;
1136                         if(cap->x1<atom->r.x)
1137                                 continue;
1138                         if(cap->y1<atom->r.y)
1139                                 continue;
1140                         if(cap->z1<atom->r.z)
1141                                 continue;
1142                 }
1143                 if(!(cap->type&CHAATTR_TOTALV))
1144                         printf("  changing attributes of atom %d (0x%x)\n",
1145                                i,cap->attr);
1146                 atom->attr=cap->attr;
1147         }
1148
1149         printf("\n\n");
1150
1151         return 0;
1152 }
1153
1154 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
1155
1156         t_stage *stage;
1157         t_chsattr_params *csp;
1158
1159         stage=mdrun->stage.current->data;
1160         csp=stage->params;
1161
1162         if(csp->type&CHSATTR_PCTRL) {
1163                 if(csp->ptau>0)
1164                         set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
1165                 else
1166                         set_p_scale(moldyn,P_SCALE_NONE,1.0);
1167         }
1168         if(csp->type&CHSATTR_TCTRL) {
1169                 if(csp->ttau>0)
1170                         set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
1171                 else
1172                         set_t_scale(moldyn,T_SCALE_NONE,1.0);
1173         }
1174         if(csp->type&CHSATTR_PRELAX) {
1175                 if(csp->dp<0)
1176                         mdrun->sattr&=(~(SATTR_PRELAX));
1177                 else
1178                         mdrun->sattr|=SATTR_PRELAX;
1179                 mdrun->dp=csp->dp;
1180         }
1181         if(csp->type&CHSATTR_TRELAX) {
1182                 if(csp->dt<0)
1183                         mdrun->sattr&=(~(SATTR_TRELAX));
1184                 else
1185                         mdrun->sattr|=SATTR_TRELAX;
1186                 mdrun->dt=csp->dt;
1187         }
1188         if(csp->type&CHSATTR_AVGRST) {
1189                 if(csp->avgrst)
1190                         mdrun->sattr|=SATTR_AVGRST;
1191                 else
1192                         mdrun->sattr&=(~(SATTR_AVGRST));
1193         }
1194         if(csp->type&CHSATTR_RSTEPS) {
1195                 mdrun->relax_steps=csp->rsteps;
1196         }
1197
1198         return 0;
1199 }
1200
1201 #define stage_print(m)  if(!(stage->executed)) \
1202                                 printf("%s",m)
1203
1204 int mdrun_hook(void *ptr1,void *ptr2) {
1205
1206         t_moldyn *moldyn;
1207         t_mdrun *mdrun;
1208         t_stage *stage;
1209         t_list *sl;
1210         int steps;
1211         u8 change_stage;
1212         t_3dvec o;
1213
1214         t_insert_atoms_params *iap;
1215         t_insert_mixed_atoms_params *imp;
1216         t_continue_params *cp;
1217         t_anneal_params *ap;
1218         t_set_temp_params *stp;
1219         t_set_timestep_params *stsp;
1220         t_fill_params *fp;
1221
1222         moldyn=ptr1;
1223         mdrun=ptr2;
1224
1225         sl=&(mdrun->stage);
1226
1227         change_stage=FALSE;
1228
1229         /* return immediately if there are no more stage */
1230         if(sl->current==NULL)
1231                 return 0;
1232
1233         /* get stage description */
1234         stage=sl->current->data;
1235
1236         /* steps in next schedule */
1237         steps=mdrun->relax_steps;
1238
1239         /* check whether relaxation steps are necessary */
1240         if((check_pressure(moldyn,mdrun)==TRUE)&\
1241            (check_temperature(moldyn,mdrun)==TRUE)) {
1242
1243                 /* be verbose */
1244                 stage_print("\n###########################\n");
1245                 stage_print("# [mdrun] executing stage #\n");
1246                 stage_print("###########################\n\n");
1247
1248                 /* stage specific stuff */
1249                 switch(stage->type) {
1250                         case STAGE_DISPLACE_ATOM:
1251                                 stage_print("  -> displace atom\n\n");
1252                                 displace_atom(moldyn,mdrun);
1253                                 change_stage=TRUE;
1254                                 break;
1255                         case STAGE_DEL_ATOMS:
1256                                 stage_print(" -> del atoms\n\n");
1257                                 del_atoms(moldyn,mdrun);
1258                                 change_stage=TRUE;
1259                                 break;
1260                         case STAGE_MODIFY_ATOMS:
1261                                 stage_print(" -> modify atoms\n\n");
1262                                 modify_atoms(moldyn,mdrun);
1263                                 change_stage=TRUE;
1264                                 break;
1265                         case STAGE_INSERT_ATOMS:
1266                                 stage_print("  -> insert atoms\n\n");
1267                                 iap=stage->params;
1268                                 if(iap->cnt_steps==iap->ins_steps) {
1269                                         change_stage=TRUE;
1270                                         break;
1271                                 }
1272                                 insert_atoms(moldyn,mdrun);
1273                                 iap->cnt_steps+=1;
1274                                 break;
1275
1276
1277
1278                         case STAGE_INSERT_MIXED_ATOMS:
1279                                 stage_print("  -> insert mixed atoms\n\n");
1280                                 imp=stage->params;
1281                                 insert_mixed_atoms(moldyn,mdrun);
1282                                 change_stage=TRUE;
1283                                 break;
1284
1285
1286
1287                         case STAGE_CONTINUE:
1288                                 stage_print("  -> continue\n\n");
1289                                 if(stage->executed==TRUE) {
1290                                         change_stage=TRUE;
1291                                         break;
1292                                 }
1293                                 cp=stage->params;
1294                                 steps=cp->runs;
1295                                 break;
1296                         case STAGE_ANNEAL:
1297                                 stage_print("  -> anneal\n\n");
1298                                 ap=stage->params;
1299                                 if(ap->count==ap->runs) {
1300                                         change_stage=TRUE;
1301                                         break;
1302                                 }
1303                                 if(moldyn->t_ref+ap->dt>=0.0)
1304                                         set_temperature(moldyn,
1305                                                         moldyn->t_ref+ap->dt);
1306                                 ap->count+=1;
1307                                 steps=ap->interval;
1308                                 break;
1309                         case STAGE_CHAATTR:
1310                                 stage_print("  -> change atom attributes\n\n");
1311                                 chaatr(moldyn,mdrun);
1312                                 change_stage=TRUE;
1313                                 break;
1314                         case STAGE_CHSATTR:
1315                                 stage_print("  -> change sys attributes\n\n");
1316                                 chsattr(moldyn,mdrun);
1317                                 change_stage=TRUE;
1318                                 break;
1319                         case STAGE_SET_TEMP:
1320                                 stage_print("  -> set temperature\n\n");
1321                                 stp=stage->params;
1322                                 if(stp->type&SET_TEMP_CURRENT) {
1323                                         set_temperature(moldyn,moldyn->t_avg);
1324                                 }
1325                                 else {
1326                                         set_temperature(moldyn,stp->val);
1327                                 }
1328                                 change_stage=TRUE;
1329                                 break;
1330                         case STAGE_SET_TIMESTEP:
1331                                 stage_print("  -> set timestep\n\n");
1332                                 stsp=stage->params;
1333                                 mdrun->timestep=stsp->tau;
1334                                 change_stage=TRUE;
1335                                 break;
1336                         case STAGE_FILL:
1337                                 stage_print("  -> fill lattice\n\n");
1338                                 fp=stage->params;
1339                                 switch(fp->lattice) {
1340                                         case ZINCBLENDE:
1341
1342                                         o.x=0.5*0.25*fp->lc;
1343                                         o.y=o.x;
1344                                         o.z=o.x;
1345                                         create_lattice(moldyn,
1346                                                        FCC,fp->lc,
1347                                                        mdrun->element1,
1348                                                        DEFAULT_ATOM_ATTR,0,
1349                                                        fp->lx,fp->ly,fp->lz,
1350                                                        &o,
1351                                                        &(fp->p_params),
1352                                                        &(fp->d_params),
1353                                                        &(fp->o_params));
1354                                         o.x+=0.25*fp->lc;
1355                                         o.y=o.x;
1356                                         o.z=o.x;
1357                                         create_lattice(moldyn,
1358                                                        FCC,fp->lc,
1359                                                        mdrun->element2,
1360                                                        DEFAULT_ATOM_ATTR,1,
1361                                                        fp->lx,fp->ly,fp->lz,
1362                                                        &o,
1363                                                        &(fp->p_params),
1364                                                        &(fp->d_params),
1365                                                        &(fp->o_params));
1366                                         break;
1367
1368                                         default:
1369
1370                                         create_lattice(moldyn,
1371                                                        fp->lattice,fp->lc,
1372                                                        fp->fill_element,
1373                                                        DEFAULT_ATOM_ATTR,
1374                                                        fp->fill_brand,
1375                                                        fp->lx,fp->ly,fp->lz,
1376                                                        NULL,
1377                                                        &(fp->p_params),
1378                                                        &(fp->d_params),
1379                                                        &(fp->o_params));
1380                                         break;
1381                                 }
1382                                 moldyn_bc_check(moldyn);
1383                                 change_stage=TRUE;
1384                                 break;
1385                         case STAGE_THERMAL_INIT:
1386                                 stage_print("  -> thermal init\n\n");
1387                                 thermal_init(moldyn,TRUE);
1388                                 change_stage=TRUE;
1389                                 break;
1390                         default:
1391                                 printf("%s unknwon stage type\n",ME);
1392                                 break;
1393                 }
1394
1395                 /* mark as executed */
1396                 stage->executed=TRUE;
1397
1398                 /* change state */
1399                 if(change_stage==TRUE) {
1400                         printf("%s finished stage\n",ME);
1401                         if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
1402                                 printf("%s no more stages\n",ME);
1403                                 return 0;
1404                         }
1405                         steps=0;
1406                 }
1407
1408         }
1409         else {
1410
1411                 /* averages */
1412                 if(mdrun->sattr&SATTR_AVGRST)
1413                         average_reset(moldyn);
1414
1415         }
1416
1417         /* continue simulation */
1418         moldyn_add_schedule(moldyn,steps,mdrun->timestep);
1419
1420         return 0;
1421 }
1422
1423 int main(int argc,char **argv) {
1424
1425         t_mdrun mdrun;
1426         t_moldyn moldyn;
1427         //t_3dvec o;
1428
1429         /* clear structs */
1430         memset(&mdrun,0,sizeof(t_mdrun));
1431         memset(&moldyn,0,sizeof(t_moldyn));
1432
1433         /* parse arguments */
1434         if(mdrun_parse_argv(&mdrun,argc,argv)<0)
1435                 return -1;
1436
1437         /* initialize list system */
1438         list_init_f(&(mdrun.stage));
1439
1440         /* parse config file */
1441         mdrun_parse_config(&mdrun);
1442
1443         /* reset the stage list */
1444         list_reset_f(&(mdrun.stage));
1445
1446         /* sanity check! */
1447
1448         /* prepare simulation */
1449
1450         if(mdrun.continue_file[0]) {
1451                 // read the save file
1452                 moldyn_read_save_file(&moldyn,mdrun.continue_file);
1453                 // manualaadjusting some stuff
1454                 moldyn.argc=argc;
1455                 moldyn.args=argv;
1456                 rand_init(&(moldyn.random),NULL,1);
1457                 moldyn.random.status|=RAND_STAT_VERBOSE;
1458         }
1459         else {
1460                 moldyn_init(&moldyn,argc,argv);
1461         }
1462
1463         if(set_int_alg(&moldyn,mdrun.intalgo)<0)
1464                 return -1;
1465
1466         /* potential */
1467         set_cutoff(&moldyn,mdrun.cutoff);
1468         if(set_potential(&moldyn,mdrun.potential)<0)
1469                 return -1;
1470         switch(mdrun.potential) {
1471                 case MOLDYN_POTENTIAL_AM:
1472                         albe_mult_set_params(&moldyn,
1473                                              mdrun.element1,
1474                                              mdrun.element2);
1475                         break;
1476                 case MOLDYN_POTENTIAL_TM:
1477                         tersoff_mult_set_params(&moldyn,
1478                                                 mdrun.element1,
1479                                                 mdrun.element2);
1480                         break;
1481                 case MOLDYN_POTENTIAL_HO:
1482                         harmonic_oscillator_set_params(&moldyn,mdrun.element1);
1483                         break;
1484                 case MOLDYN_POTENTIAL_LJ:
1485                         lennard_jones_set_params(&moldyn,mdrun.element1);
1486                         break;
1487                 default:
1488                         printf("%s unknown potential: %02x\n",
1489                                ME,mdrun.potential);
1490                         return -1;
1491         }
1492
1493         /* if it is a continue run, reset schedule and skip lattice init */
1494         if(mdrun.continue_file[0]) {
1495                 memset(&(moldyn.schedule),0,sizeof(t_moldyn_schedule));
1496                 goto addsched;
1497         }
1498
1499         /* initial lattice and dimensions */
1500         set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
1501         set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
1502         /* replaced by fill stage !!
1503         switch(mdrun.lattice) {
1504                 case FCC:
1505                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.fill_element,
1506                                        mdrun.m1,DEFAULT_ATOM_ATTR,
1507                                        mdrun.fill_brand,mdrun.lx,
1508                                        mdrun.ly,mdrun.lz,NULL,0,NULL);
1509                         break;
1510                 case DIAMOND:
1511                         create_lattice(&moldyn,DIAMOND,mdrun.lc,
1512                                        mdrun.fill_element,
1513                                        mdrun.m1,DEFAULT_ATOM_ATTR,
1514                                        mdrun.fill_brand,mdrun.lx,
1515                                        mdrun.ly,mdrun.lz,NULL,0,NULL);
1516                         break;
1517                 case ZINCBLENDE:
1518                         o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
1519                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
1520                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
1521                                        mdrun.ly,mdrun.lz,&o,0,NULL);
1522                         o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
1523                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
1524                                        mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx,
1525                                        mdrun.ly,mdrun.lz,&o,0,NULL);
1526                         break;
1527                 case NONE:
1528                         set_nn_dist(&moldyn,mdrun.nnd);
1529                         break;
1530                 default:
1531                         printf("%s unknown lattice type: %02x\n",
1532                                ME,mdrun.lattice);
1533                         return -1;
1534         }
1535         moldyn_bc_check(&moldyn);
1536         replaced by fill stage !! */
1537
1538         /* temperature and pressure */
1539         set_temperature(&moldyn,mdrun.temperature);
1540         set_pressure(&moldyn,mdrun.pressure);
1541         /* replaced thermal_init stage
1542         thermal_init(&moldyn,TRUE);
1543         */
1544
1545 addsched:
1546         /* first schedule */
1547         moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
1548
1549         /* log */
1550         moldyn_set_log_dir(&moldyn,mdrun.sdir);
1551         moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
1552         if(mdrun.elog)
1553                 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
1554         if(mdrun.tlog)
1555                 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
1556         if(mdrun.plog)
1557                 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
1558         if(mdrun.vlog)
1559                 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
1560         if(mdrun.visualize)
1561                 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
1562         if(mdrun.save)
1563                 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
1564         moldyn_set_log(&moldyn,CREATE_REPORT,0);
1565         set_avg_skip(&moldyn,mdrun.avgskip);
1566
1567         /* prepare the hook function */
1568         moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
1569
1570         /* run the simulation */
1571         moldyn_integrate(&moldyn);
1572
1573         /* shutdown */
1574         moldyn_shutdown(&moldyn);
1575         del_stages(&mdrun);
1576         list_destroy_f(&(mdrun.stage));
1577
1578         return 0;
1579 }