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 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 add_stage(t_mdrun *mdrun,u8 type,void *params) {
  87
  88         int psize;
  89
  90         t_stage *stage;
  91
  92         switch(type) {
  93                 case STAGE_INSERT_ATOMS:
  94                         psize=sizeof(t_insert_atoms_params);
  95                         break;
  96                 case STAGE_CONTINUE:
  97                         psize=sizeof(t_continue_params);
  98                         break;
  99                 case STAGE_ANNEAL:
 100                         psize=sizeof(t_anneal_params);
 101                         break;
 102                 case STAGE_CHAATTR:
 103                         psize=sizeof(t_chaattr_params);
 104                         break;
 105                 case STAGE_CHSATTR:
 106                         psize=sizeof(t_chsattr_params);
 107                         break;
 108                 default:
 109
 110         }
 111
 112         stage=malloc(sizeof(t_stage));
 113         if(stage==NULL) {
 114                 perror("[mdrun] malloc (add stage)");
 115                 return -1;
 116         }
 117
 118         stage->type=type;
 119         stage->executed=FALSE;
 120
 121         stage->params=malloc(psize);
 122         if(stage->params==NULL) {
 123                 perror("[mdrun] malloc (stage params)");
 124                 return -1;
 125         }
 126
 127         memcpy(stage->params,params,psize);
 128
 129         list_add_immediate_f(mdrun->stage,stage);
 130
 131         return 0;
 132 }
 133
 134 int mdrun_parse_config(t_mdrun *mdrun) {
 135
 136         int fd,ret;
 137         char error[128];
 138         char line[128];
 139         char *wptr;
 140         char word[16][32];
 141         int wcnt;
 142         int i;
 143
 144         t_insert_atoms_params iap;
 145         t_continue_params cp;
 146         t_anneal_params ap;
 147         t_chaattr_params cap;
 148         t_chsattr_params csp;
 149
 150         /* open config file */
 151         fd=open(mdrun->cfile,O_RDONLY);
 152         if(fd<0) {
 153                 snprintf(error,128,"%s open cfile %s",ME,mdrun->cfile);
 154                 perror(error);
 155                 return fd;
 156         }
 157
 158         /* read, parse, set */
 159         ret=1;
 160         while(ret>0) {
 161
 162                 /* read */
 163                 ret=get_line(fd,line,128);
 164
 165                 /* parse */
 166
 167                 // ignore # lines and \n
 168                 if((line[0]=='#')|(ret==1))
 169                         continue;
 170
 171                 // get command + args
 172                 wcnt=0;
 173                 while(1) {
 174                         if(wcnt)
 175                                 wptr=strtok(NULL," \t");
 176                         else
 177                                 wptr=strtok(line," \t");
 178                         if(wptr==NULL)
 179                                 break;
 180                         strncpy(word[wcnt],wptr,32);
 181                         wcnt+=1;
 182                 }
 183
 184                 if(!strncmp(word[0],"potential",9)) {
 185                         if(!strncmp(word[1],"albe",4))
 186                                 mdrun->potential=MOLDYN_POTENTIAL_AM;
 187                         if(!strncmp(word[1],"tersoff",7))
 188                                 mdrun->potential=MOLDYN_POTENTIAL_TM;
 189                         if(!strncmp(word[1],"ho",2))
 190                                 mdrun->potential=MOLDYN_POTENTIAL_HO;
 191                         if(!strncmp(word[1],"lj",2))
 192                                 mdrun->potential=MOLDYN_POTENTIAL_LJ;
 193                 }
 194                 else if(!strncmp(word[0],"cutoff",6))
 195                         mdrun->cutoff=atof(word[1]);
 196                 else if(!strncmp(word[0],"intalgo",7)) {
 197                         if(!strncmp(word[1],"verlet",5))
 198                                 mdrun->intalgo=MOLDYN_INTEGRATE_VERLET;
 199                 }
 200                 else if(!strncmp(word[0],"timestep",8))
 201                         mdrun->timestep=atof(word[1]);
 202                 else if(!strncmp(word[0],"volume",6)) {
 203                         mdrun->dim.x=atof(word[1]);
 204                         mdrun->dim.y=atof(word[2]);
 205                         mdrun->dim.z=atof(word[3]);
 206                         if(strncmp(word[4],"0",1))
 207                                 mdrun->vis=TRUE;
 208                 }
 209                 else if(!strncmp(word[0],"pbc",3)) {
 210                         if(strncmp(word[1],"0",1))
 211                                 mdrun->pbcx=TRUE;
 212                         else
 213                                 mdrun->pbcx=FALSE;
 214                         if(strncmp(word[2],"0",1))
 215                                 mdrun->pbcy=TRUE;
 216                         else
 217                                 mdrun->pbcy=FALSE;
 218                         if(strncmp(word[3],"0",1))
 219                                 mdrun->pbcz=TRUE;
 220                         else
 221                                 mdrun->pbcz=FALSE;
 222                 }
 223                 else if(!strncmp(word[0],"temperature",11))
 224                         mdrun->temperature=atof(word[1]);
 225                 else if(!strncmp(word[0],"pressure",8))
 226                         mdrun->pressure=atof(word[1]);
 227                 else if(!strncmp(word[0],"lattice",7)) {
 228                         if(!strncmp(word[1],"zincblende",10))
 229                                 mdrun->lattice=ZINCBLENDE;
 230                         if(!strncmp(word[1],"fcc",3))
 231                                 mdrun->lattice=FCC;
 232                         if(!strncmp(word[1],"diamond",7))
 233                                 mdrun->lattice=DIAMOND;
 234                 }
 235                 else if(!strncmp(word[0],"element1",8)) {
 236                         mdrun->element1=atoi(word[1]);
 237                         switch(mdrun->element1) {
 238                                 case SI:
 239                                         mdrun->m1=M_SI;
 240                                         break;
 241                                 case C:
 242                                         mdrun->m1=M_C;
 243                                         break;
 244                                 default:
 245                                         printf("%s unknown element1: %s|%d\n",
 246                                                ME,word[1],mdrun->element1);
 247                                         return -1;
 248                         }
 249                 }
 250                 else if(!strncmp(word[0],"element2",8)) {
 251                         mdrun->element2=atoi(word[1]);
 252                         switch(mdrun->element2) {
 253                                 case SI:
 254                                         mdrun->m2=M_SI;
 255                                         break;
 256                                 case C:
 257                                         mdrun->m2=M_C;
 258                                         break;
 259                                 default:
 260                                         printf("%s unknown element2: %s|%d\n",
 261                                                ME,word[1],mdrun->element2);
 262                                         return -1;
 263                         }
 264                 }
 265                 else if(!strncmp(word[0],"fill",6)) {
 266                         // only lc mode by now
 267                         mdrun->lx=atoi(word[2]);
 268                         mdrun->ly=atoi(word[3]);
 269                         mdrun->lz=atoi(word[4]);
 270                 }
 271                 else if(!strncmp(word[0],"aattr",5)) {
 272                         // for aatrib line we need a special stage
 273                         // containing one schedule of 0 loops ...
 274                         if(!strncmp(word[1],"all",3)) {
 275                                 cap.type=
 276
 277
 278                                 HIER WEITER
 279                         }
 280                         for(i=0;i<strlen(wptr)) {
 281                                 switch(word[2][i]) {
 282                                         case 'b':
 283                                                 cap.attr|=ATOM_ATTR_VB;
 284                                                 break;
 285                                         case 'h':
 286                                                 cap.attr|=ATOM_ATTR_HB;
 287                                                 break;
 288                                         case 'v':
 289                                                 cap.attr|=ATOM_ATTR_VA;
 290                                                 break;
 291                                         case 'f':
 292                                                 cap.attr|=ATOM_ATTR_FP;
 293                                                 break;
 294                                         case '1':
 295                                                 cap.attr|=ATOM_ATTR_1BP;
 296                                                 break;
 297                                         case '2':
 298                                                 cap.attr|=ATOM_ATTR_2BP;
 299                                                 break;
 300                                         case '3':
 301                                                 cap.attr|=ATOM_ATTR_3BP;
 302                                                 break;
 303                                         default:
 304                                                 break;
 305                                 }
 306                         }
 307                         add_stage(mdrun,STAGE_CHAATTR,&cap);
 308                 }
 309                 else if(!strncmp(word[0],"ectrl",5)) {
 310                         mdrun->p_tau=atof(word[1]);
 311                         mdrun->t_tau=atof(word[2]);
 312                 }
 313                 else if(!strncmp(word[0],"prerun",6))
 314                         mdrun->prerun=atoi(word[1]);
 315                 else if(!strncmp(word[0],"avgskip",7))
 316                         mdrun->avgskip=atoi(word[1]);
 317                 else if(!strncmp(word[0],"elog",4))
 318                         mdrun->elog=atoi(word[1]);
 319                 else if(!strncmp(word[0],"tlog",4))
 320                         mdrun->tlog=atoi(word[1]);
 321                 else if(!strncmp(word[0],"plog",4))
 322                         mdrun->plog=atoi(word[1]);
 323                 else if(!strncmp(word[0],"vlog",4))
 324                         mdrun->vlog=atoi(word[1]);
 325                 else if(!strncmp(word[0],"save",4))
 326                         mdrun->save=atoi(word[1]);
 327                 else if(!strncmp(word[0],"visualize",9))
 328                         mdrun->visualize=atoi(word[1]);
 329                 else if(!strncmp(word[0],"stage",5)) {
 330                         // for every stage line, add a stage
 331                         printf("%s stage: %s\n",ME,word[1]);
 332                 }
 333                 else {
 334                         printf("%s unknown command %s, skipped!\n",ME,word[0]);
 335                         continue;
 336                 }
 337         }
 338
 339         /* close file */
 340         close(fd);
 341
 342         return 0;
 343 }
 344
 345 int check_pressure(t_moldyn *moldyn,t_mdrun *mdrun) {
 346
 347         double delta;
 348
 349         if(!(mdrun->sattr&SATTR_PRELAX))
 350                 return TRUE;
 351
 352         delta=moldyn->p_avg-moldyn->p_ref;
 353
 354         if(delta<0)
 355                 delta=-delta;
 356
 357         if(delta>mdrun->dp)
 358                 return FALSE;
 359
 360         return TRUE;
 361 }
 362
 363 int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) {
 364
 365         double delta;
 366
 367         if(!(mdrun->sattr&SATTR_TRELAX))
 368                 return TRUE;
 369
 370         delta=moldyn->t_avg-moldyn->t_ref;
 371
 372         if(delta<0)
 373                 delta=-delta;
 374
 375         if(delta>mdrun->dt)
 376                 return FALSE;
 377
 378         return TRUE;
 379 }
 380
 381 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 382
 383         insert_atoms_params *iap;
 384         t_stage *stage;
 385         t_atom *atom;
 386         t_3dvec r,v,dist;
 387         double d,dmin;
 388         int cnt,i;
 389         double x,y,z;
 390         double x0,y0,z0;
 391         u8 cr_check,run;
 392
 393         stage=mdrun->stage->current->data;
 394         iap=stage->params;
 395
 396         cr_check=FALSE;
 397
 398         v.x=0.0; v.y=0.0; v.z=0.0;
 399
 400         switch(iap->type) {
 401                 case INS_TOTAL:
 402                         x=moldyn->dim.x;
 403                         x0=0.0;
 404                         y=moldyn->dim.y;
 405                         y0=0.0;
 406                         z=moldyn->dim.z;
 407                         z0=0.0;
 408                         cr_check=TRUE;
 409                         break;
 410                 case INS_REGION:
 411                         x=iap->x1-iap->x0;
 412                         x0=iap->x0;
 413                         y=iap->y1-iap->y0;
 414                         y0=iap->y0;
 415                         z=iap->z1-iap->z0;
 416                         z0=iap->z0;
 417                         cr_check=TRUE;
 418                         break;
 419                 default:
 420                         printf("%s unknown insertion mode\n");
 421                         return -1;
 422         }
 423
 424         cnt=0;
 425         while(cnt<iap->ins_atoms) {
 426                 run=1;
 427                 while(run) {
 428                         r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0;
 429                         r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0;
 430                         r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0;
 431                         run=0;
 432                         if(cr_check==TRUE) {
 433                                 dmin=1000;      // for sure too high!
 434                                 for(i=0;i<moldyn->count;i++) {
 435                                         atom=&(moldyn->atom[i]);
 436                                         v3_sub(&dist,&(atom->r),&r);
 437                                         check_per_bound(moldyn,&dist);
 438                                         d=v3_absolute_square(&dist);
 439                                         if(d<(iap->cr*iap->cr)) {
 440                                                 run=1;
 441                                                 break;
 442                                         }
 443                                         if(d<dmin)
 444                                                 dmin=d;
 445                                 }
 446                         }
 447                 }
 448                 add_atom(moldyn,iap->element,pse_mass[iap->element],
 449                          iap->brand,iap->attr,&r,&v);
 450                 printf("%s atom inserted: %f %f %f | d squared = %f\n",
 451                        ME,r.x,r.y,r.z,dmin);
 452                 cnt+=1;
 453         }
 454
 455         return 0;
 456 }
 457
 458 int chaatrib(t_moldyn *moldyn,t_mdrun *mdrun) {
 459
 460         t_stage *stage;
 461         t_chaattr_params *cap;
 462         t_atom *atom;
 463         u8 assigne;
 464
 465         stage=mdrun->stage->current->data;
 466         cap=stage->params;
 467
 468         for(i=0;i<moldyn->count;i++) {
 469                 atom=&(moldyn->atom[i]);
 470                 if(cap->type&CHAATTR_ELEMENT) {
 471                         if(cap->element!=atom->element)
 472                                 continue;
 473                 }
 474                 if(cap->type&CHAATTR_REGION) {
 475                         if(cap->x0<atom->r.x)
 476                                 continue;
 477                         if(cap->y0<atom->r.y)
 478                                 continue;
 479                         if(cap->z0<atom->r.z)
 480                                 continue;
 481                         if(cap->x1>atom->r.x)
 482                                 continue;
 483                         if(cap->y1>atom->r.y)
 484                                 continue;
 485                         if(cap->z1>atom->r.z)
 486                                 continue;
 487                 }
 488                 atom->attr=cap->attr;
 489         }
 490
 491         return 0;
 492 }
 493
 494 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
 495
 496         t_stage *stage;
 497         t_chsattr_params *csp;
 498
 499         stage=mdrun->stage->current->data;
 500         csp=stage->params;
 501
 502         switch(csp->type) {
 503                 case CHSATTR_PCTRL:
 504                         set_p_scale(moldyn,csp->ctrl_type,csp->tau);
 505                         break;
 506                 case CHSATTR_TCTRL:
 507                         set_t_scale(moldyn,csp->ctrl_type,csp->tau);
 508                         break;
 509                 case CHSATTR_PRELAX:
 510                         mdrun->sattr&=(^(SATTR_PRELAX));
 511                         if(csp->ctrl==TRUE)
 512                                 mdrun->sattr|=SATTR_PRELAX;
 513                         break;
 514                 case CHSATTR_TRELAX:
 515                         mdrun->sattr&=(^(SATTR_TRELAX));
 516                         if(csp->ctrl==TRUE)
 517                                 mdrun->sattr|=SATTR_TRELAX;
 518                         break;
 519                 case CHSATTR_AVGRST:
 520                         mdrun->sattr&=(^(SATTR_AVGRST));
 521                         if(csp->ctrl==TRUE)
 522                                 mdrun->sattr|=SATTR_AVGRST;
 523                         break;
 524                 default:
 525                         printf("%s unknown system attribute change\n",ME);
 526                         return -1;
 527         }
 528
 529         return 0;
 530 }
 531
 532 int mdrun_hook(t_moldyn *moldyn,void *ptr) {
 533
 534         t_mdrun *mdrun;
 535         t_stage *stage;
 536         t_list *sl;
 537         int steps;
 538         double tau;
 539         u8 change_stage;
 540
 541         t_insert_atoms_params *iap;
 542         t_continue_params *cp;
 543         t_anneal_params *ap;
 544
 545         mdrun=ptr;
 546         sl=mdrun->stage;
 547
 548         change_stage=FALSE;
 549
 550         /* return immediately if there are no more stage */
 551         if(sl->current==NULL)
 552                 return 0;
 553
 554         /* get stage description */
 555         stage=mdrun->sl->current->data;
 556
 557         /* default steps and tau values */
 558         steps=mdrun->relax_steps;
 559         tau=mdrun->timestep;
 560
 561         /* check whether relaxation steps are necessary */
 562         if(!((check_pressure==FALSE)|(check_temperature==FALSE))) {
 563
 564                 /* stage specific stuff */
 565                 switch(stage->type) {
 566                         case STAGE_INSERT_ATOMS:
 567                                 iap=stage->params;
 568                                 if(iap->cnt_steps==iap->ins_steps) {
 569                                         change_stage=TRUE;
 570                                         break;
 571                                 }
 572                                 insert_atoms(moldyn,mdrun);
 573                                 iap->cnt_steps+=1;
 574                                         break;
 575                         case STAGE_CONTINUE:
 576                                 if(stage->executed==TRUE) {
 577                                         change_stage=TRUE;
 578                                         break;
 579                                 }
 580                                 cp=stage->params;
 581                                 steps=cp->runs;
 582                                         break;
 583                         case STAGE_ANNEAL:
 584                                 ap=stage->params;
 585                                 if(ap->count==ap->runs) {
 586                                         change_stage=TRUE;
 587                                         break;
 588                                 }
 589                                 set_temperature(moldyn,moldyn->t_ref+ap->dt);
 590                                 ap->count+=1;
 591                                 break;
 592                         case STAGE_CHAATTR:
 593                                 chaatrib(moldyn,mdrun);
 594                                 change_stage=TRUE;
 595                                 break;
 596                         case STAGE_CHSATTR:
 597                                 chsatrib(moldyn,mdrun);
 598                                 change_stage=TRUE;
 599                                 break;
 600                         default:
 601                                 printf("%s unknwon stage type\n",ME);
 602                                 break;
 603                 }
 604
 605                 /* mark as executed */
 606                 stage->executed=TRUE;
 607
 608                 /* change state */
 609                 if(change_stage==TRUE) {
 610                         printf("%s finished stage\n",ME);
 611                         if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
 612                                 printf("%s no more stages\n",ME);
 613                                 return 0;
 614                         }
 615                         steps=0;
 616                         tau=mdrun->timestep;
 617                 }
 618
 619         }
 620         else {
 621
 622                 /* averages */
 623                 if(mdrun->sattr&SATTR_AVGRST)
 624                         average_reset(moldyn);
 625
 626         }
 627
 628         /* continue simulation */
 629         moldyn_add_schedule(moldyn,steps,tau);
 630
 631         return 0;
 632 }
 633
 634 int main(int argc,char **argv) {
 635
 636         t_mdrun mdrun;
 637         t_moldyn moldyn;
 638         t_3dvec o;
 639
 640         /* clear mdrun struct */
 641         memset(&mdrun,0,sizeof(t_mdrun));
 642
 643         /* parse arguments */
 644         if(mdrun_parse_argv(&mdrun,argc,argv)<0)
 645                 return -1;
 646
 647         /* initialize list system */
 648         list_init_f(&(mdrun.stage));
 649
 650         /* parse config file */
 651         mdrun_parse_config(&mdrun);
 652
 653         /* reset the stage list */
 654         list_reset_f(&(mdrun.stage));
 655
 656         /* sanity check! */
 657
 658         /* prepare simulation */
 659         moldyn_init(&moldyn,argc,argv);
 660         if(set_int_alg(&moldyn,mdrun.intalgo)<0)
 661                 return -1;
 662         set_cutoff(&moldyn,mdrun.cutoff);
 663         if(set_potential(&moldyn,mdrun.potential)<0)
 664                 return -1;
 665         switch(mdrun.potential) {
 666                 case MOLDYN_POTENTIAL_AM:
 667                         albe_mult_set_params(&moldyn,
 668                                              mdrun.element1,
 669                                              mdrun.element2);
 670                         break;
 671                 case MOLDYN_POTENTIAL_TM:
 672                         tersoff_mult_set_params(&moldyn,
 673                                                 mdrun.element1,
 674                                                 mdrun.element2);
 675                         break;
 676                 case MOLDYN_POTENTIAL_HO:
 677                         harmonic_oscillator_set_params(&moldyn,mdrun.element1);
 678                         break;
 679                 case MOLDYN_POTENTIAL_LJ:
 680                         lennard_jones_set_params(&moldyn,mdrun.element1);
 681                         break;
 682                 default:
 683                         printf("%s unknown potential: %02x\n",
 684                                ME,mdrun.potential);
 685                         return -1;
 686         }
 687         set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
 688         set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
 689         switch(mdrun.lattice) {
 690                 case FCC:
 691                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
 692                                        mdrun.m1,0,0,mdrun.lx,
 693                                        mdrun.ly,mdrun.lz,NULL);
 694                         break;
 695                 case DIAMOND:
 696                         create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1,
 697                                        mdrun.m1,0,0,mdrun.lx,
 698                                        mdrun.ly,mdrun.lz,NULL);
 699                         break;
 700                 case ZINCBLENDE:
 701                         o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
 702                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
 703                                        mdrun.m1,0,0,mdrun.lx,
 704                                        mdrun.ly,mdrun.lz,&o);
 705                         o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
 706                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
 707                                        mdrun.m2,0,0,mdrun.lx,
 708                                        mdrun.ly,mdrun.lz,&o);
 709                         break;
 710                 default:
 711                         printf("%s unknown lattice type: %02x\n",
 712                                ME,mdrun.lattice);
 713                         return -1;
 714         }
 715         moldyn_bc_check(&moldyn);
 716         set_temperature(&moldyn,mdrun.temperature);
 717         set_pressure(&moldyn,mdrun.pressure);
 718         thermal_init(&moldyn,TRUE);
 719         moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
 720         moldyn_set_log_dir(&moldyn,mdrun.sdir);
 721         moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
 722         if(mdrun.elog)
 723                 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
 724         if(mdrun.tlog)
 725                 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
 726         if(mdrun.plog)
 727                 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
 728         if(mdrun.vlog)
 729                 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
 730         if(mdrun.visualize)
 731                 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
 732         if(mdrun.save)
 733                 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
 734         moldyn_set_log(&moldyn,CREATE_REPORT,0);
 735         set_avg_skip(&moldyn,mdrun.avgskip);
 736
 737         /* pt scaling */
 738         if(mdrun.p_tau!=0)
 739                 set_p_scale(&moldyn,P_SCALE_BERENDSEN,mdrun.p_tau);
 740         if(mdrun.t_tau!=0)
 741                 set_t_scale(&moldyn,T_SCALE_BERENDSEN,mdrun.t_tau);
 742
 743         /* prepare the hook function */
 744         moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
 745
 746         /* run the simulation */
 747         moldyn_integrate(&moldyn);
 748
 749         /* shutdown */
 750         moldyn_shutdown(&moldyn);
 751         del_stages(&mdrun);
 752         list_destroy_f(&(mdrun.stage));
 753
 754         return 0;
 755 }