mdrun.c

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