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
 807         stage=mdrun->stage.current->data;
 808         delp=stage->params;
 809
 810         for(i=0;i<moldyn->count;i++) {
 811                 v3_sub(&dist,&(delp->o),&(moldyn->atom[i].r));
 812 //printf("%d ----> %f %f %f = %f | %f\n",i,dist.x,dist.y,dist.z,v3_absolute_square(&dist),delp->r*delp->r);
 813                 if(v3_absolute_square(&dist)<=(delp->r*delp->r)) {
 814                         del_atom(moldyn,moldyn->atom[i].tag);
 815                         printf("%s atom deleted: %d %d %d\n",ME,
 816                                moldyn->atom[i].tag,moldyn->atom[i].element,
 817                                moldyn->atom[i].brand);
 818                 }
 819         }
 820
 821         return 0;
 822 }
 823
 824 int modify_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 825
 826         t_modify_atoms_params *map;
 827         t_stage *stage;
 828         t_atom *atom;
 829         t_3dvec v;
 830         int i;
 831
 832         atom=moldyn->atom;
 833         stage=mdrun->stage.current->data;
 834         map=stage->params;
 835         v.x=0.0; v.y=0.0; v.z=0.0;
 836
 837         for(i=0;i<moldyn->count;i++) {
 838                 if(atom[i].tag==map->tag) {
 839                         v.x=sqrt(2.0*fabs(map->ekin.x)/atom[i].mass);
 840                         if(map->ekin.x<0.0)
 841                                 v.x=-v.x;
 842                         v.y=sqrt(2.0*fabs(map->ekin.y)/atom[i].mass);
 843                         if(map->ekin.y<0.0)
 844                                 v.y=-v.y;
 845                         v.z=sqrt(2.0*fabs(map->ekin.z)/atom[i].mass);
 846                         if(map->ekin.z<0.0)
 847                                 v.z=-v.z;
 848                         v3_copy(&(atom[i].v),&v);
 849                         printf("%s atom modified: v = (%f %f %f)\n",
 850                                ME,v.x,v.y,v.z);
 851                 }
 852         }
 853
 854         return 0;
 855 }
 856
 857 int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 858
 859         t_insert_atoms_params *iap;
 860         t_stage *stage;
 861         t_atom *atom;
 862         t_3dvec r,v,dist;
 863         double d,dmin,o;
 864         int cnt,i;
 865         double x,y,z;
 866         double x0,y0,z0;
 867         u8 cr_check,run;
 868
 869         stage=mdrun->stage.current->data;
 870         iap=stage->params;
 871
 872         cr_check=FALSE;
 873
 874         v.x=0.0; v.y=0.0; v.z=0.0;
 875         x=0; y=0; z=0;
 876         o=0; dmin=0;
 877
 878         switch(mdrun->lattice) {
 879                 case CUBIC:
 880                         o=0.5*mdrun->lc;
 881                         break;
 882                 case FCC:
 883                         o=0.25*mdrun->lc;
 884                         break;
 885                 case DIAMOND:
 886                 case ZINCBLENDE:
 887                         o=0.125*mdrun->lc;
 888                         break;
 889                 default:
 890                         break;
 891         }
 892
 893         switch(iap->type) {
 894                 case INS_TOTAL:
 895                         x=moldyn->dim.x;
 896                         x0=-x/2.0;
 897                         y=moldyn->dim.y;
 898                         y0=-y/2.0;
 899                         z=moldyn->dim.z;
 900                         z0=-z/2.0;
 901                         cr_check=TRUE;
 902                         break;
 903                 case INS_RECT:
 904                         x=iap->x1-iap->x0;
 905                         x0=iap->x0;
 906                         y=iap->y1-iap->y0;
 907                         y0=iap->y0;
 908                         z=iap->z1-iap->z0;
 909                         z0=iap->z0;
 910                         cr_check=TRUE;
 911                         break;
 912                 case INS_SPHERE:
 913                         x=2.0*iap->x1;
 914                         x0=iap->x0-iap->x1;
 915                         y=x;
 916                         y0=iap->y0-iap->x1;
 917                         z=x;
 918                         z0=iap->z0-iap->x1;
 919                         cr_check=TRUE;
 920                         break;
 921                 case INS_POS:
 922                 case INS_RELPOS:
 923                         x0=iap->x0;
 924                         y0=iap->y0;
 925                         z0=iap->z0;
 926                         cr_check=FALSE;
 927                         break;
 928                 default:
 929                         printf("%s unknown insertion mode: %02x\n",
 930                                ME,iap->type);
 931                         return -1;
 932         }
 933
 934         cnt=0;
 935         while(cnt<iap->ins_atoms) {
 936                 run=1;
 937                 while(run) {
 938                         if((iap->type!=INS_POS)&&(iap->type!=INS_RELPOS)) {
 939                                 r.x=rand_get_double(&(moldyn->random))*x;
 940                                 r.y=rand_get_double(&(moldyn->random))*y;
 941                                 r.z=rand_get_double(&(moldyn->random))*z;
 942                         }
 943                         else {
 944                                 r.x=0.0;
 945                                 r.y=0.0;
 946                                 r.z=0.0;
 947                         }
 948                         if(iap->type==INS_RELPOS) {
 949                                 r.x+=x0*mdrun->lc;
 950                                 r.y+=y0*mdrun->lc;
 951                                 r.z+=z0*mdrun->lc;
 952                         }
 953                         else {
 954                                 r.x+=x0;
 955                                 r.y+=y0;
 956                                 r.z+=z0;
 957                         }
 958                         // offset
 959                         if(iap->type!=INS_TOTAL) {
 960                                 r.x+=o;
 961                                 r.y+=o;
 962                                 r.z+=o;
 963                         }
 964                         run=0;
 965                         if(cr_check==TRUE) {
 966                                 dmin=1000;      // for sure too high!
 967                                 for(i=0;i<moldyn->count;i++) {
 968                                         atom=&(moldyn->atom[i]);
 969                                         v3_sub(&dist,&(atom->r),&r);
 970                                         check_per_bound(moldyn,&dist);
 971                                         d=v3_absolute_square(&dist);
 972                                         if(d<(iap->cr*iap->cr)) {
 973                                                 run=1;
 974                                                 break;
 975                                         }
 976                                         if(d<dmin)
 977                                                 dmin=d;
 978                                 }
 979                         }
 980                         if(iap->type==INS_SPHERE) {
 981                                 if((r.x-iap->x0)*(r.x-iap->x0)+
 982                                    (r.y-iap->y0)*(r.y-iap->y0)+
 983                                    (r.z-iap->z0)*(r.z-iap->z0)>
 984                                    (iap->x1*iap->x1)) {
 985                                         run=1;
 986                                 }
 987                         }
 988                 }
 989                 add_atom(moldyn,iap->element,
 990                          iap->brand,iap->attr,&r,&v);
 991                 printf("%s atom inserted (%d/%d): %f %f %f\n",
 992                        ME,(iap->cnt_steps+1)*iap->ins_atoms,
 993                        iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z);
 994                 printf("  attributes: ");
 995                 if(iap->attr&ATOM_ATTR_VB)
 996                         printf("b ");
 997                 if(iap->attr&ATOM_ATTR_HB)
 998                         printf("h ");
 999                 if(iap->attr&ATOM_ATTR_VA)
1000                         printf("v ");
1001                 if(iap->attr&ATOM_ATTR_FP)
1002                         printf("f ");
1003                 if(iap->attr&ATOM_ATTR_1BP)
1004                         printf("1 ");
1005                 if(iap->attr&ATOM_ATTR_2BP)
1006                         printf("2 ");
1007                 if(iap->attr&ATOM_ATTR_3BP)
1008                         printf("3 ");
1009                 printf("\n");
1010                 printf("  d2 = %f/%f\n",dmin,iap->cr*iap->cr);
1011                 cnt+=1;
1012         }
1013
1014         return 0;
1015 }
1016
1017 int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
1018
1019         t_stage *stage;
1020         t_insert_mixed_atoms_params *imp;
1021         t_atom *atom;
1022         double x,x0,y,y0,z,z0;
1023         double dmin,d,cmin,cmax;
1024         u8 retry;
1025         t_3dvec r,v,dist;
1026         int i;
1027
1028
1029         stage=mdrun->stage.current->data;
1030         imp=stage->params;
1031
1032         x=moldyn->dim.x;
1033         x0=-x/2.0;
1034         y=moldyn->dim.y;
1035         y0=-y/2.0;
1036         z=moldyn->dim.z;
1037         z0=-z/2.0;
1038
1039         v.x=0.0;
1040         v.y=0.0;
1041         v.z=0.0;
1042
1043         cmin=imp->crmin*imp->crmin;
1044         cmax=imp->crmax*imp->crmax;
1045
1046         while(imp->amount1|imp->amount2) {
1047                 if(imp->amount1) {
1048                         retry=1;
1049                         while(retry) {
1050                                 retry=0;
1051                                 r.x=rand_get_double(&(moldyn->random))*x;
1052                                 r.y=rand_get_double(&(moldyn->random))*y;
1053                                 r.z=rand_get_double(&(moldyn->random))*z;
1054                                 r.x+=x0;
1055                                 r.y+=y0;
1056                                 r.z+=z0;
1057                                 dmin=1000.0;    // for sure too high!
1058                                 for(i=0;i<moldyn->count;i++) {
1059                                         atom=&(moldyn->atom[i]);
1060                                         v3_sub(&dist,&(atom->r),&r);
1061                                         check_per_bound(moldyn,&dist);
1062                                         d=v3_absolute_square(&dist);
1063                                         if(d<cmin) {
1064                                                 retry=1;
1065                                                 break;
1066                                         }
1067                                         if(d<dmin)
1068                                                 dmin=d;
1069                                 }
1070                                 if(dmin!=1000.0)
1071                                         if(dmin>cmax)
1072                                                 retry=1;
1073                         }
1074                         add_atom(moldyn,imp->element1,
1075                                  imp->brand1,imp->attr1,&r,&v);
1076                         printf("%s (mixed) atom inserted (%d): %f %f %f\n",
1077                                ME,imp->amount1,r.x,r.y,r.z);
1078                         printf("  -> d2 = %f/%f/%f\n",dmin,cmin,cmax);
1079                         imp->amount1-=1;
1080                 }
1081                 if(imp->amount2) {
1082                         retry=1;
1083                         while(retry) {
1084                                 retry=0;
1085                                 r.x=rand_get_double(&(moldyn->random))*x;
1086                                 r.y=rand_get_double(&(moldyn->random))*y;
1087                                 r.z=rand_get_double(&(moldyn->random))*z;
1088                                 r.x+=x0;
1089                                 r.y+=y0;
1090                                 r.z+=z0;
1091                                 dmin=1000.0;    // for sure too high!
1092                                 for(i=0;i<moldyn->count;i++) {
1093                                         atom=&(moldyn->atom[i]);
1094                                         v3_sub(&dist,&(atom->r),&r);
1095                                         check_per_bound(moldyn,&dist);
1096                                         d=v3_absolute_square(&dist);
1097                                         if(d<cmin) {
1098                                                 retry=1;
1099                                                 break;
1100                                         }
1101                                         if(d<dmin)
1102                                                 dmin=d;
1103                                 }
1104                                 if(dmin!=1000.0)
1105                                         if(dmin>cmax)
1106                                                 retry=1;
1107                         }
1108                         add_atom(moldyn,imp->element2,
1109                                  imp->brand2,imp->attr2,&r,&v);
1110                         printf("%s (mixed) atom inserted (%d): %f %f %f\n",
1111                                ME,imp->amount2,r.x,r.y,r.z);
1112                         printf("  -> d2 = %f/%f/%f\n",dmin,cmin,cmax);
1113                         imp->amount2-=1;
1114                 }
1115         }
1116
1117         return 0;
1118 }
1119
1120 int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) {
1121
1122         t_stage *stage;
1123         t_chaattr_params *cap;
1124         t_atom *atom;
1125         int i;
1126
1127         stage=mdrun->stage.current->data;
1128         cap=stage->params;
1129
1130         for(i=0;i<moldyn->count;i++) {
1131                 atom=&(moldyn->atom[i]);
1132                 if(cap->type&CHAATTR_ELEMENT) {
1133                         if(cap->element!=atom->element)
1134                                 continue;
1135                 }
1136                 if(cap->type&CHAATTR_NUMBER) {
1137                         if(cap->element!=atom->tag)
1138                                 continue;
1139                 }
1140                 if(cap->type&CHAATTR_REGION) {
1141                         if(cap->x0>atom->r.x)
1142                                 continue;
1143                         if(cap->y0>atom->r.y)
1144                                 continue;
1145                         if(cap->z0>atom->r.z)
1146                                 continue;
1147                         if(cap->x1<atom->r.x)
1148                                 continue;
1149                         if(cap->y1<atom->r.y)
1150                                 continue;
1151                         if(cap->z1<atom->r.z)
1152                                 continue;
1153                 }
1154                 if(!(cap->type&CHAATTR_TOTALV))
1155                         printf("  changing attributes of atom %d (0x%x)\n",
1156                                i,cap->attr);
1157                 atom->attr=cap->attr;
1158         }
1159
1160         printf("\n\n");
1161
1162         return 0;
1163 }
1164
1165 int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) {
1166
1167         t_stage *stage;
1168         t_chsattr_params *csp;
1169
1170         stage=mdrun->stage.current->data;
1171         csp=stage->params;
1172
1173         if(csp->type&CHSATTR_PCTRL) {
1174                 if(csp->ptau>0)
1175                         set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau);
1176                 else
1177                         set_p_scale(moldyn,P_SCALE_NONE,1.0);
1178         }
1179         if(csp->type&CHSATTR_TCTRL) {
1180                 if(csp->ttau>0)
1181                         set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau);
1182                 else
1183                         set_t_scale(moldyn,T_SCALE_NONE,1.0);
1184         }
1185         if(csp->type&CHSATTR_PRELAX) {
1186                 if(csp->dp<0)
1187                         mdrun->sattr&=(~(SATTR_PRELAX));
1188                 else
1189                         mdrun->sattr|=SATTR_PRELAX;
1190                 mdrun->dp=csp->dp;
1191         }
1192         if(csp->type&CHSATTR_TRELAX) {
1193                 if(csp->dt<0)
1194                         mdrun->sattr&=(~(SATTR_TRELAX));
1195                 else
1196                         mdrun->sattr|=SATTR_TRELAX;
1197                 mdrun->dt=csp->dt;
1198         }
1199         if(csp->type&CHSATTR_AVGRST) {
1200                 if(csp->avgrst)
1201                         mdrun->sattr|=SATTR_AVGRST;
1202                 else
1203                         mdrun->sattr&=(~(SATTR_AVGRST));
1204         }
1205         if(csp->type&CHSATTR_RSTEPS) {
1206                 mdrun->relax_steps=csp->rsteps;
1207         }
1208
1209         return 0;
1210 }
1211
1212 int crt(t_moldyn *moldyn,t_mdrun *mdrun) {
1213
1214         t_stage *stage;
1215         t_crt_params *crtp;
1216
1217         int fd;
1218         char line[128];
1219         char *wptr;
1220         int acount;
1221         int ret;
1222         void *ptr;
1223
1224         t_atom *atom;
1225         t_3dvec disp;
1226         double frac;
1227         int i;
1228
1229         stage=mdrun->stage.current->data;
1230         crtp=stage->params;
1231
1232         acount=0;
1233
1234         /* initial stuff */
1235
1236         if(crtp->count==0) {
1237                 printf("  crt init\n");
1238                 // read final positions, constraints and do the alloc
1239                 fd=open(crtp->file,O_RDONLY);
1240                 if(fd<0) {
1241                         perror("[mdrun] FATAL reading constraints file");
1242                         return fd;
1243                 }
1244                 while(1) {
1245                         ret=get_line(fd,line,128);
1246                         // check for end of file
1247                         if(ret<=0) {
1248                                 printf("  read %d atom positions\n",acount);
1249                                 if(acount!=moldyn->count)
1250                                         printf("  atom count mismatch!!!\n");
1251                                 printf("\n");
1252                                 break;
1253                         }
1254                         // ignore # lines and \n
1255                         if((line[0]=='#')|(ret==1))
1256                                 continue;
1257                         // allocate new memory
1258                         ptr=realloc(crtp->r_fin,(acount+1)*sizeof(t_3dvec));
1259                         if(ptr==NULL) {
1260                                 perror("[mdrun] FATAL realloc crt positions");
1261                                 return -1;
1262                         }
1263                         crtp->r_fin=ptr;
1264                         ptr=realloc(constraints,(acount+1)*3*sizeof(u8));
1265                         if(ptr==NULL) {
1266                                 perror("[mdrun] FATAL realloc crt constraints");
1267                                 return -1;
1268                         }
1269                         constraints=ptr;
1270                         // ignore type
1271                         wptr=strtok(line," \t");
1272                         // read x y z
1273                         wptr=strtok(NULL," \t");
1274                         crtp->r_fin[acount].x=atof(wptr);
1275                         wptr=strtok(NULL," \t");
1276                         crtp->r_fin[acount].y=atof(wptr);
1277                         wptr=strtok(NULL," \t");
1278                         crtp->r_fin[acount].z=atof(wptr);
1279                         // read constraints
1280                         wptr=strtok(NULL," \t");
1281                         constraints[3*acount]=atoi(wptr);
1282                         wptr=strtok(NULL," \t");
1283                         constraints[3*acount+1]=atoi(wptr);
1284                         wptr=strtok(NULL," \t");
1285                         constraints[3*acount+2]=atoi(wptr);
1286                         // done reading
1287                         acount+=1;
1288                 }
1289                 close(fd);
1290                 // allocate trafo angles
1291                 trafo_angle=malloc(acount*2*sizeof(double));
1292                 if(trafo_angle==NULL) {
1293                         perror("[mdrun] FATAL alloc trafo angles");
1294                         return -1;
1295                 }
1296                 // set crt mode
1297                 crtt=crtp->type;
1298         }
1299
1300         /* write a save file s-crt_xofy.save */
1301         snprintf(line,128,"%s/s-crt_%03dof%03d.save",
1302                  moldyn->vlsdir,crtp->count,crtp->steps);
1303         fd=open(line,O_WRONLY|O_TRUNC|O_CREAT,S_IRUSR|S_IWUSR);
1304         if(fd<0) perror("[mdrun] crt save fd open");
1305         else {
1306                 write(fd,moldyn,sizeof(t_moldyn));
1307                 write(fd,moldyn->atom,
1308                       moldyn->count*sizeof(t_atom));
1309         }
1310         close(fd);
1311         /* visualize atoms */
1312         visual_atoms(moldyn);
1313
1314         /* output energy */
1315         printf("  crt energy: %d - %f\n\n",
1316                crtp->count,(moldyn->ekin+moldyn->energy)/EV);
1317
1318         /* crt routines: calculate displacement + set individual constraints */
1319
1320         printf("  crt step %d of %d in total\n\n",crtp->count+1,crtp->steps);
1321
1322         if((crtp->type==1)|(crtp->count==0))
1323                 printf("  crt angle update\n\n");
1324
1325         for(i=0;i<moldyn->count;i++) {
1326                 // calc displacements
1327                 atom=moldyn->atom;
1328                 v3_sub(&disp,&(crtp->r_fin[i]),&(atom[i].r));
1329                 // angles
1330                 if((crtp->type==1)|(crtp->count==0)) {
1331                         trafo_angle[2*i]=atan2(disp.x,disp.y);
1332                         trafo_angle[2*i+1]=-atan2(disp.z,
1333                                            sqrt(disp.x*disp.x+disp.y*disp.y));
1334                 }
1335                 // move atoms
1336                 frac=1.0/(crtp->steps-crtp->count);
1337                 v3_scale(&disp,&disp,frac);
1338                 v3_add(&(atom[i].r),&(atom[i].r),&disp);
1339         }
1340
1341         return 0;
1342 }
1343
1344 #define stage_print(m)  if(!(stage->executed)) \
1345                                 printf("%s",m)
1346
1347 int mdrun_hook(void *ptr1,void *ptr2) {
1348
1349         t_moldyn *moldyn;
1350         t_mdrun *mdrun;
1351         t_stage *stage;
1352         t_list *sl;
1353         int steps;
1354         u8 change_stage;
1355         t_3dvec o;
1356
1357         t_insert_atoms_params *iap;
1358         t_insert_mixed_atoms_params *imp;
1359         t_continue_params *cp;
1360         t_anneal_params *ap;
1361         t_set_temp_params *stp;
1362         t_set_timestep_params *stsp;
1363         t_fill_params *fp;
1364         t_crt_params *crtp;
1365
1366         moldyn=ptr1;
1367         mdrun=ptr2;
1368
1369         sl=&(mdrun->stage);
1370
1371         change_stage=FALSE;
1372
1373         /* return immediately if there are no more stage */
1374         if(sl->current==NULL)
1375                 return 0;
1376
1377         /* get stage description */
1378         stage=sl->current->data;
1379
1380         /* steps in next schedule */
1381         steps=mdrun->relax_steps;
1382
1383         /* check whether relaxation steps are necessary */
1384         if((check_pressure(moldyn,mdrun)==TRUE)&\
1385            (check_temperature(moldyn,mdrun)==TRUE)) {
1386
1387                 /* be verbose */
1388                 stage_print("\n###########################\n");
1389                 stage_print("# [mdrun] executing stage #\n");
1390                 stage_print("###########################\n\n");
1391
1392                 /* stage specific stuff */
1393                 switch(stage->type) {
1394                         case STAGE_DISPLACE_ATOM:
1395                                 stage_print("  -> displace atom\n\n");
1396                                 displace_atom(moldyn,mdrun);
1397                                 change_stage=TRUE;
1398                                 break;
1399                         case STAGE_DEL_ATOMS:
1400                                 stage_print(" -> del atoms\n\n");
1401                                 del_atoms(moldyn,mdrun);
1402                                 change_stage=TRUE;
1403                                 break;
1404                         case STAGE_MODIFY_ATOMS:
1405                                 stage_print(" -> modify atoms\n\n");
1406                                 modify_atoms(moldyn,mdrun);
1407                                 change_stage=TRUE;
1408                                 break;
1409                         case STAGE_INSERT_ATOMS:
1410                                 stage_print("  -> insert atoms\n\n");
1411                                 iap=stage->params;
1412                                 if(iap->cnt_steps==iap->ins_steps) {
1413                                         change_stage=TRUE;
1414                                         break;
1415                                 }
1416                                 insert_atoms(moldyn,mdrun);
1417                                 iap->cnt_steps+=1;
1418                                 break;
1419
1420
1421
1422                         case STAGE_INSERT_MIXED_ATOMS:
1423                                 stage_print("  -> insert mixed atoms\n\n");
1424                                 imp=stage->params;
1425                                 insert_mixed_atoms(moldyn,mdrun);
1426                                 change_stage=TRUE;
1427                                 break;
1428
1429
1430
1431                         case STAGE_CONTINUE:
1432                                 stage_print("  -> continue\n\n");
1433                                 if(stage->executed==TRUE) {
1434                                         change_stage=TRUE;
1435                                         break;
1436                                 }
1437                                 cp=stage->params;
1438                                 steps=cp->runs;
1439                                 break;
1440                         case STAGE_ANNEAL:
1441                                 stage_print("  -> anneal\n\n");
1442                                 ap=stage->params;
1443                                 if(ap->count==ap->runs) {
1444                                         change_stage=TRUE;
1445                                         break;
1446                                 }
1447                                 if(moldyn->t_ref+ap->dt>=0.0)
1448                                         set_temperature(moldyn,
1449                                                         moldyn->t_ref+ap->dt);
1450                                 ap->count+=1;
1451                                 steps=ap->interval;
1452                                 break;
1453                         case STAGE_CHAATTR:
1454                                 stage_print("  -> change atom attributes\n\n");
1455                                 chaatr(moldyn,mdrun);
1456                                 change_stage=TRUE;
1457                                 break;
1458                         case STAGE_CHSATTR:
1459                                 stage_print("  -> change sys attributes\n\n");
1460                                 chsattr(moldyn,mdrun);
1461                                 change_stage=TRUE;
1462                                 break;
1463                         case STAGE_SET_TEMP:
1464                                 stage_print("  -> set temperature\n\n");
1465                                 stp=stage->params;
1466                                 if(stp->type&SET_TEMP_CURRENT) {
1467                                         set_temperature(moldyn,moldyn->t_avg);
1468                                 }
1469                                 else {
1470                                         set_temperature(moldyn,stp->val);
1471                                 }
1472                                 change_stage=TRUE;
1473                                 break;
1474                         case STAGE_SET_TIMESTEP:
1475                                 stage_print("  -> set timestep\n\n");
1476                                 stsp=stage->params;
1477                                 mdrun->timestep=stsp->tau;
1478                                 change_stage=TRUE;
1479                                 break;
1480                         case STAGE_FILL:
1481                                 stage_print("  -> fill lattice\n\n");
1482                                 fp=stage->params;
1483                                 switch(fp->lattice) {
1484                                         case ZINCBLENDE:
1485
1486                                         o.x=0.5*0.25*fp->lc;
1487                                         o.y=o.x;
1488                                         o.z=o.x;
1489                                         create_lattice(moldyn,
1490                                                        FCC,fp->lc,
1491                                                        mdrun->element1,
1492                                                        DEFAULT_ATOM_ATTR,0,
1493                                                        fp->lx,fp->ly,fp->lz,
1494                                                        &o,
1495                                                        &(fp->p_params),
1496                                                        &(fp->d_params),
1497                                                        &(fp->o_params));
1498                                         o.x+=0.25*fp->lc;
1499                                         o.y=o.x;
1500                                         o.z=o.x;
1501                                         create_lattice(moldyn,
1502                                                        FCC,fp->lc,
1503                                                        mdrun->element2,
1504                                                        DEFAULT_ATOM_ATTR,1,
1505                                                        fp->lx,fp->ly,fp->lz,
1506                                                        &o,
1507                                                        &(fp->p_params),
1508                                                        &(fp->d_params),
1509                                                        &(fp->o_params));
1510                                         break;
1511
1512                                         default:
1513
1514                                         create_lattice(moldyn,
1515                                                        fp->lattice,fp->lc,
1516                                                        fp->fill_element,
1517                                                        DEFAULT_ATOM_ATTR,
1518                                                        fp->fill_brand,
1519                                                        fp->lx,fp->ly,fp->lz,
1520                                                        NULL,
1521                                                        &(fp->p_params),
1522                                                        &(fp->d_params),
1523                                                        &(fp->o_params));
1524                                         break;
1525                                 }
1526                                 moldyn_bc_check(moldyn);
1527                                 change_stage=TRUE;
1528                                 break;
1529                         case STAGE_THERMAL_INIT:
1530                                 stage_print("  -> thermal init\n\n");
1531                                 thermal_init(moldyn,TRUE);
1532                                 change_stage=TRUE;
1533                                 break;
1534                         case STAGE_CRT:
1535                                 stage_print("  -> constraint relaxation");
1536                                 stage_print(" technique\n\n");
1537                                 crtp=stage->params;
1538                                 if(crtp->count==crtp->steps) {
1539                                         free(constraints);
1540                                         free(trafo_angle);
1541                                         free(crtp->r_fin);
1542                                         change_stage=TRUE;
1543                                         break;
1544                                 }
1545                                 crt(moldyn,mdrun);
1546                                 crtp->count+=1;
1547                                 break;
1548                         default:
1549                                 printf("%s unknwon stage type\n",ME);
1550                                 break;
1551                 }
1552
1553                 /* mark as executed */
1554                 stage->executed=TRUE;
1555
1556                 /* change state */
1557                 if(change_stage==TRUE) {
1558                         printf("%s finished stage\n",ME);
1559                         if(list_next_f(sl)==L_NO_NEXT_ELEMENT) {
1560                                 printf("%s no more stages\n",ME);
1561                                 return 0;
1562                         }
1563                         steps=0;
1564                 }
1565
1566         }
1567         else {
1568
1569                 /* averages */
1570                 if(mdrun->sattr&SATTR_AVGRST)
1571                         average_reset(moldyn);
1572
1573         }
1574
1575         /* continue simulation */
1576         moldyn_add_schedule(moldyn,steps,mdrun->timestep);
1577
1578         return 0;
1579 }
1580
1581 int main(int argc,char **argv) {
1582
1583         t_mdrun mdrun;
1584         t_moldyn moldyn;
1585         //t_3dvec o;
1586
1587         /* clear structs */
1588         memset(&mdrun,0,sizeof(t_mdrun));
1589         memset(&moldyn,0,sizeof(t_moldyn));
1590
1591         /* init crt variables */
1592         crtt=0;
1593         constraints=NULL;
1594         trafo_angle=NULL;
1595
1596         /* parse arguments */
1597         if(mdrun_parse_argv(&mdrun,argc,argv)<0)
1598                 return -1;
1599
1600         /* initialize list system */
1601         list_init_f(&(mdrun.stage));
1602
1603         /* parse config file */
1604         mdrun_parse_config(&mdrun);
1605
1606         /* reset the stage list */
1607         list_reset_f(&(mdrun.stage));
1608
1609         /* sanity check! */
1610
1611         /* prepare simulation */
1612
1613         if(mdrun.continue_file[0]) {
1614                 // read the save file
1615                 moldyn_read_save_file(&moldyn,mdrun.continue_file);
1616                 // manualaadjusting some stuff
1617                 moldyn.argc=argc;
1618                 moldyn.args=argv;
1619                 rand_init(&(moldyn.random),NULL,1);
1620                 moldyn.random.status|=RAND_STAT_VERBOSE;
1621         }
1622         else {
1623                 moldyn_init(&moldyn,argc,argv);
1624         }
1625
1626         if(set_int_alg(&moldyn,mdrun.intalgo)<0)
1627                 return -1;
1628
1629         /* potential */
1630         set_cutoff(&moldyn,mdrun.cutoff);
1631         if(set_potential(&moldyn,mdrun.potential)<0)
1632                 return -1;
1633         switch(mdrun.potential) {
1634                 case MOLDYN_POTENTIAL_AM:
1635                         albe_mult_set_params(&moldyn,
1636                                              mdrun.element1,
1637                                              mdrun.element2);
1638                         break;
1639                 case MOLDYN_POTENTIAL_TM:
1640                         tersoff_mult_set_params(&moldyn,
1641                                                 mdrun.element1,
1642                                                 mdrun.element2);
1643                         break;
1644                 case MOLDYN_POTENTIAL_HO:
1645                         harmonic_oscillator_set_params(&moldyn,mdrun.element1);
1646                         break;
1647                 case MOLDYN_POTENTIAL_LJ:
1648                         lennard_jones_set_params(&moldyn,mdrun.element1);
1649                         break;
1650                 default:
1651                         printf("%s unknown potential: %02x\n",
1652                                ME,mdrun.potential);
1653                         return -1;
1654         }
1655
1656         /* if it is a continue run, reset schedule and skip lattice init */
1657         if(mdrun.continue_file[0]) {
1658                 memset(&(moldyn.schedule),0,sizeof(t_moldyn_schedule));
1659                 goto addsched;
1660         }
1661
1662         /* initial lattice and dimensions */
1663         set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
1664         set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
1665         /* replaced by fill stage !!
1666         switch(mdrun.lattice) {
1667                 case FCC:
1668                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.fill_element,
1669                                        mdrun.m1,DEFAULT_ATOM_ATTR,
1670                                        mdrun.fill_brand,mdrun.lx,
1671                                        mdrun.ly,mdrun.lz,NULL,0,NULL);
1672                         break;
1673                 case DIAMOND:
1674                         create_lattice(&moldyn,DIAMOND,mdrun.lc,
1675                                        mdrun.fill_element,
1676                                        mdrun.m1,DEFAULT_ATOM_ATTR,
1677                                        mdrun.fill_brand,mdrun.lx,
1678                                        mdrun.ly,mdrun.lz,NULL,0,NULL);
1679                         break;
1680                 case ZINCBLENDE:
1681                         o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
1682                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
1683                                        mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx,
1684                                        mdrun.ly,mdrun.lz,&o,0,NULL);
1685                         o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
1686                         create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
1687                                        mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx,
1688                                        mdrun.ly,mdrun.lz,&o,0,NULL);
1689                         break;
1690                 case NONE:
1691                         set_nn_dist(&moldyn,mdrun.nnd);
1692                         break;
1693                 default:
1694                         printf("%s unknown lattice type: %02x\n",
1695                                ME,mdrun.lattice);
1696                         return -1;
1697         }
1698         moldyn_bc_check(&moldyn);
1699         replaced by fill stage !! */
1700
1701         /* temperature and pressure */
1702         set_temperature(&moldyn,mdrun.temperature);
1703         set_pressure(&moldyn,mdrun.pressure);
1704         /* replaced thermal_init stage
1705         thermal_init(&moldyn,TRUE);
1706         */
1707
1708 addsched:
1709         /* first schedule */
1710         moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
1711
1712         /* log */
1713         moldyn_set_log_dir(&moldyn,mdrun.sdir);
1714         moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
1715         if(mdrun.elog)
1716                 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
1717         if(mdrun.tlog)
1718                 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
1719         if(mdrun.plog)
1720                 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
1721         if(mdrun.vlog)
1722                 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
1723         if(mdrun.visualize)
1724                 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
1725         if(mdrun.save)
1726                 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
1727         moldyn_set_log(&moldyn,CREATE_REPORT,0);
1728         set_avg_skip(&moldyn,mdrun.avgskip);
1729
1730         /* prepare the hook function */
1731         moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);
1732
1733         /* run the simulation */
1734         moldyn_integrate(&moldyn);
1735
1736         /* shutdown */
1737         moldyn_shutdown(&moldyn);
1738         del_stages(&mdrun);
1739         list_destroy_f(&(mdrun.stage));
1740
1741         return 0;
1742 }