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