mdrun.c

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