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