src/tools/gmx_membed.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <signal.h>
  40 #include <stdlib.h>
  41 #include "typedefs.h"
  42 #include "sysstuff.h"
  43 #include "statutil.h"
  44 #include "macros.h"
  45 #include "copyrite.h"
  46 #include "main.h"
  47 #include "gmx_ana.h"
  48
  49 #ifdef GMX_LIB_MPI
  50 #include <mpi.h>
  51 #endif
  52 #ifdef GMX_THREADS
  53 #include "tmpi.h"
  54 #endif
  55
  56 /* afm stuf */
  57 #include "pull.h"
  58
  59 /* We use the same defines as in mvdata.c here */
  60 #define  block_bc(cr,   d) gmx_bcast(     sizeof(d),     &(d),(cr))
  61 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
  62 #define   snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
  63
  64 /* The following two variables and the signal_handler function
  65  * is used from pme.c as well
  66  */
  67
  68 typedef struct {
  69         t_state s;
  70         rvec    *f;
  71         real    epot;
  72         real    fnorm;
  73         real    fmax;
  74         int     a_fmax;
  75 } em_state_t;
  76
  77 typedef struct {
  78         int    it_xy;
  79         int    it_z;
  80         int    xy_step;
  81         int    z_step;
  82         rvec    xmin;
  83         rvec    xmax;
  84         rvec    *geom_cent;
  85         int    pieces;
  86         int    *nidx;
  87         atom_id **subindex;
  88 } pos_ins_t;
  89
  90 typedef struct {
  91         int             id;
  92         char    *name;
  93         int     nr;
  94         int     natoms;     /*nr of atoms per lipid*/
  95         int     mol1;       /*id of the first lipid molecule*/
  96         real    area;
  97 } lip_t;
  98
  99 typedef struct {
 100         char    *name;
 101         t_block mem_at;
 102         int             *mol_id;
 103         int             nmol;
 104         real    lip_area;
 105         real    zmin;
 106         real    zmax;
 107         real    zmed;
 108 } mem_t;
 109
 110 typedef struct {
 111         int             *mol;
 112         int             *block;
 113         int     nr;
 114 } rmm_t;
 115
 116 int search_string(char *s,int ng,char ***gn)
 117 {
 118         int i;
 119
 120         for(i=0; (i<ng); i++)
 121                 if (gmx_strcasecmp(s,*gn[i]) == 0)
 122                         return i;
 123
 124         gmx_fatal(FARGS,"Group %s not found in indexfile.\nMaybe you have non-default groups in your mdp file, while not using the '-n' option of grompp.\nIn that case use the '-n' option.\n",s);
 125
 126         return -1;
 127 }
 128
 129 int get_mol_id(int at,int nmblock,gmx_molblock_t *mblock, int *type, int *block)
 130 {
 131         int mol_id=0;
 132         int i;
 133
 134         for(i=0;i<nmblock;i++)
 135         {
 136                 if(at<(mblock[i].nmol*mblock[i].natoms_mol))
 137                 {
 138                         mol_id+=at/mblock[i].natoms_mol;
 139                         *type = mblock[i].type;
 140                         *block = i;
 141                         return mol_id;
 142                 } else {
 143                         at-= mblock[i].nmol*mblock[i].natoms_mol;
 144                         mol_id+=mblock[i].nmol;
 145                 }
 146         }
 147
 148         gmx_fatal(FARGS,"Something is wrong in mol ids, at %d, mol_id %d",at,mol_id);
 149
 150         return -1;
 151 }
 152
 153 int get_block(int mol_id,int nmblock,gmx_molblock_t *mblock)
 154 {
 155         int i;
 156         int nmol=0;
 157
 158         for(i=0;i<nmblock;i++)
 159         {
 160                 nmol+=mblock[i].nmol;
 161                 if(mol_id<nmol)
 162                         return i;
 163         }
 164
 165         gmx_fatal(FARGS,"mol_id %d larger than total number of molecules %d.\n",mol_id,nmol);
 166
 167         return -1;
 168 }
 169
 170 int get_tpr_version(const char *infile)
 171 {
 172         char    buf[STRLEN];
 173         gmx_bool        bDouble;
 174         int     precision,fver;
 175         t_fileio *fio;
 176
 177         fio = open_tpx(infile,"r");
 178         gmx_fio_checktype(fio);
 179
 180         precision = sizeof(real);
 181
 182         gmx_fio_do_string(fio,buf);
 183         if (strncmp(buf,"VERSION",7))
 184                 gmx_fatal(FARGS,"Can not read file %s,\n"
 185                                 "             this file is from a Gromacs version which is older than 2.0\n"
 186                                 "             Make a new one with grompp or use a gro or pdb file, if possible",
 187                                 gmx_fio_getname(fio));
 188         gmx_fio_do_int(fio,precision);
 189         bDouble = (precision == sizeof(double));
 190         if ((precision != sizeof(float)) && !bDouble)
 191                 gmx_fatal(FARGS,"Unknown precision in file %s: real is %d bytes "
 192                                 "instead of %d or %d",
 193                                 gmx_fio_getname(fio),precision,sizeof(float),sizeof(double));
 194         gmx_fio_setprecision(fio,bDouble);
 195         fprintf(stderr,"Reading file %s, %s (%s precision)\n",
 196                         gmx_fio_getname(fio),buf,bDouble ? "double" : "single");
 197
 198         gmx_fio_do_int(fio,fver);
 199
 200         close_tpx(fio);
 201
 202         return fver;
 203 }
 204
 205 void set_inbox(int natom, rvec *x)
 206 {
 207         rvec tmp;
 208         int  i;
 209
 210         tmp[XX]=tmp[YY]=tmp[ZZ]=0.0;
 211         for(i=0;i<natom;i++)
 212         {
 213                 if(x[i][XX]<tmp[XX])            tmp[XX]=x[i][XX];
 214                 if(x[i][YY]<tmp[YY])            tmp[YY]=x[i][YY];
 215                 if(x[i][ZZ]<tmp[ZZ])            tmp[ZZ]=x[i][ZZ];
 216         }
 217
 218         for(i=0;i<natom;i++)
 219                         rvec_inc(x[i],tmp);
 220 }
 221
 222 int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
 223 {
 224         int i,j,nr,mol_id;
 225         int type=0,block=0;
 226         gmx_bool bNEW;
 227
 228         nr=0;
 229         snew(tlist->index,at->nr);
 230         for (i=0;i<at->nr;i++)
 231         {
 232                 bNEW=TRUE;
 233                 mol_id = get_mol_id(at->index[i],mtop->nmolblock,mtop->molblock,&type,&block);
 234                 for(j=0;j<nr;j++)
 235                 {
 236                         if(tlist->index[j]==type)
 237                                                 bNEW=FALSE;
 238                 }
 239                 if(bNEW==TRUE)
 240                 {
 241                         tlist->index[nr]=type;
 242                         nr++;
 243                 }
 244         }
 245
 246         srenew(tlist->index,nr);
 247         return nr;
 248 }
 249
 250 void check_types(t_block *ins_at,t_block *rest_at,gmx_mtop_t *mtop)
 251 {
 252         t_block         *ins_mtype,*rest_mtype;
 253         int                     i,j;
 254
 255         snew(ins_mtype,1);
 256         snew(rest_mtype,1);
 257     ins_mtype->nr  = get_mtype_list(ins_at , mtop, ins_mtype );
 258     rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
 259
 260     for(i=0;i<ins_mtype->nr;i++)
 261     {
 262         for(j=0;j<rest_mtype->nr;j++)
 263         {
 264                 if(ins_mtype->index[i]==rest_mtype->index[j])
 265                         gmx_fatal(FARGS,"Moleculetype %s is found both in the group to insert and the rest of the system.\n"
 266                                         "Because we need to exclude all interactions between the atoms in the group to\n"
 267                                         "insert, the same moleculetype can not be used in both groups. Change the\n"
 268                                         "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
 269                                         "an appropriate *.itp file",*(mtop->moltype[rest_mtype->index[j]].name),
 270                                         *(mtop->moltype[rest_mtype->index[j]].name));
 271         }
 272     }
 273
 274     sfree(ins_mtype->index);
 275     sfree(rest_mtype->index);
 276     sfree(ins_mtype);
 277     sfree(rest_mtype);
 278 }
 279
 280 int init_ins_at(t_block *ins_at,t_block *rest_at,t_state *state, pos_ins_t *pos_ins,gmx_groups_t *groups,int ins_grp_id, real xy_max)
 281 {
 282         int i,gid,c=0;
 283         real x,xmin,xmax,y,ymin,ymax,z,zmin,zmax;
 284
 285         snew(rest_at->index,state->natoms);
 286
 287         xmin=xmax=state->x[ins_at->index[0]][XX];
 288         ymin=ymax=state->x[ins_at->index[0]][YY];
 289         zmin=zmax=state->x[ins_at->index[0]][ZZ];
 290
 291         for(i=0;i<state->natoms;i++)
 292         {
 293                 gid = groups->grpnr[egcFREEZE][i];
 294                 if(groups->grps[egcFREEZE].nm_ind[gid]==ins_grp_id)
 295                 {
 296                         x=state->x[i][XX];
 297                         if (x<xmin)                     xmin=x;
 298                         if (x>xmax)                     xmax=x;
 299                         y=state->x[i][YY];
 300                         if (y<ymin)                             ymin=y;
 301                         if (y>ymax)                             ymax=y;
 302                         z=state->x[i][ZZ];
 303                         if (z<zmin)                             zmin=z;
 304                         if (z>zmax)                             zmax=z;
 305                 } else {
 306                         rest_at->index[c]=i;
 307                         c++;
 308                 }
 309         }
 310
 311         rest_at->nr=c;
 312         srenew(rest_at->index,c);
 313
 314         if(xy_max>1.000001)
 315         {
 316                 pos_ins->xmin[XX]=xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
 317                 pos_ins->xmin[YY]=ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
 318
 319                 pos_ins->xmax[XX]=xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
 320                 pos_ins->xmax[YY]=ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
 321         } else {
 322                 pos_ins->xmin[XX]=xmin;
 323                 pos_ins->xmin[YY]=ymin;
 324
 325                 pos_ins->xmax[XX]=xmax;
 326                 pos_ins->xmax[YY]=ymax;
 327         }
 328
 329         /* 6.0 is estimated thickness of bilayer */
 330         if( (zmax-zmin) < 6.0 )
 331         {
 332                 pos_ins->xmin[ZZ]=zmin+(zmax-zmin)/2.0-3.0;
 333                 pos_ins->xmax[ZZ]=zmin+(zmax-zmin)/2.0+3.0;
 334         } else {
 335                 pos_ins->xmin[ZZ]=zmin;
 336                 pos_ins->xmax[ZZ]=zmax;
 337         }
 338
 339         return c;
 340 }
 341
 342 real est_prot_area(pos_ins_t *pos_ins,rvec *r,t_block *ins_at, mem_t *mem_p)
 343 {
 344         real x,y,dx=0.15,dy=0.15,area=0.0;
 345         real add;
 346         int c,at;
 347
 348         for(x=pos_ins->xmin[XX];x<pos_ins->xmax[XX];x+=dx)
 349         {
 350                 for(y=pos_ins->xmin[YY];y<pos_ins->xmax[YY];y+=dy)
 351                 {
 352                         c=0;
 353                         add=0.0;
 354                         do
 355                         {
 356                                 at=ins_at->index[c];
 357                                 if ( (r[at][XX]>=x) && (r[at][XX]<x+dx) &&
 358                                                 (r[at][YY]>=y) && (r[at][YY]<y+dy) &&
 359                                                 (r[at][ZZ]>mem_p->zmin+1.0) && (r[at][ZZ]<mem_p->zmax-1.0) )
 360                                         add=1.0;
 361                                 c++;
 362                         } while ( (c<ins_at->nr) && (add<0.5) );
 363                         area+=add;
 364                 }
 365         }
 366         area=area*dx*dy;
 367
 368         return area;
 369 }
 370
 371 void init_lip(matrix box, gmx_mtop_t *mtop, lip_t *lip)
 372 {
 373         int i;
 374         real mem_area;
 375         int mol1=0;
 376
 377         mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
 378         for(i=0;i<mtop->nmolblock;i++)
 379         {
 380                 if(mtop->molblock[i].type == lip->id)
 381                 {
 382                         lip->nr=mtop->molblock[i].nmol;
 383                         lip->natoms=mtop->molblock[i].natoms_mol;
 384                 }
 385         }
 386         lip->area=2.0*mem_area/(double)lip->nr;
 387
 388         for (i=0;i<lip->id;i++)
 389                 mol1+=mtop->molblock[i].nmol;
 390         lip->mol1=mol1;
 391 }
 392
 393 int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
 394 {
 395         int i,j,at,mol,nmol,nmolbox,count;
 396         t_block *mem_a;
 397         real z,zmin,zmax,mem_area;
 398         gmx_bool bNew;
 399         atom_id *mol_id;
 400         int type=0,block=0;
 401
 402         nmol=count=0;
 403         mem_a=&(mem_p->mem_at);
 404         snew(mol_id,mem_a->nr);
 405 /*      snew(index,mem_a->nr); */
 406         zmin=pos_ins->xmax[ZZ];
 407         zmax=pos_ins->xmin[ZZ];
 408         for(i=0;i<mem_a->nr;i++)
 409         {
 410                 at=mem_a->index[i];
 411                 if(     (r[at][XX]>pos_ins->xmin[XX]) && (r[at][XX]<pos_ins->xmax[XX]) &&
 412                         (r[at][YY]>pos_ins->xmin[YY]) && (r[at][YY]<pos_ins->xmax[YY]) &&
 413                         (r[at][ZZ]>pos_ins->xmin[ZZ]) && (r[at][ZZ]<pos_ins->xmax[ZZ]) )
 414                 {
 415                         mol = get_mol_id(at,mtop->nmolblock,mtop->molblock,&type,&block);
 416
 417                         bNew=TRUE;
 418                         for(j=0;j<nmol;j++)
 419                                 if(mol == mol_id[j])
 420                                         bNew=FALSE;
 421
 422                         if(bNew)
 423                         {
 424                                 mol_id[nmol]=mol;
 425                                 nmol++;
 426                         }
 427
 428                         z=r[at][ZZ];
 429                         if(z<zmin)                                      zmin=z;
 430                         if(z>zmax)                                      zmax=z;
 431
 432 /*                      index[count]=at;*/
 433                         count++;
 434                 }
 435         }
 436
 437         mem_p->nmol=nmol;
 438         srenew(mol_id,nmol);
 439         mem_p->mol_id=mol_id;
 440 /*      srenew(index,count);*/
 441 /*      mem_p->mem_at.nr=count;*/
 442 /*      sfree(mem_p->mem_at.index);*/
 443 /*      mem_p->mem_at.index=index;*/
 444
 445         if((zmax-zmin)>(box[ZZ][ZZ]-0.5))
 446                 gmx_fatal(FARGS,"Something is wrong with your membrane. Max and min z values are %f and %f.\n"
 447                                 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
 448                                 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
 449                                 "your water layer is not thick enough.\n",zmax,zmin);
 450         mem_p->zmin=zmin;
 451         mem_p->zmax=zmax;
 452         mem_p->zmed=(zmax-zmin)/2+zmin;
 453
 454         /*number of membrane molecules in protein box*/
 455         nmolbox = count/mtop->molblock[block].natoms_mol;
 456         /*mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
 457         mem_p->lip_area = 2.0*mem_area/(double)mem_p->nmol;*/
 458         mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
 459         mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
 460
 461         return mem_p->mem_at.nr;
 462 }
 463
 464 void init_resize(t_block *ins_at,rvec *r_ins,pos_ins_t *pos_ins,mem_t *mem_p,rvec *r, gmx_bool bALLOW_ASYMMETRY)
 465 {
 466         int i,j,at,c,outsidesum,gctr=0;
 467     int idxsum=0;
 468
 469     /*sanity check*/
 470     for (i=0;i<pos_ins->pieces;i++)
 471           idxsum+=pos_ins->nidx[i];
 472     if (idxsum!=ins_at->nr)
 473           gmx_fatal(FARGS,"Piecewise sum of inserted atoms not same as size of group selected to insert.");
 474
 475     snew(pos_ins->geom_cent,pos_ins->pieces);
 476     for (i=0;i<pos_ins->pieces;i++)
 477     {
 478         c=0;
 479         outsidesum=0;
 480         for(j=0;j<DIM;j++)
 481                 pos_ins->geom_cent[i][j]=0;
 482
 483         for(j=0;j<DIM;j++)
 484                 pos_ins->geom_cent[i][j]=0;
 485         for (j=0;j<pos_ins->nidx[i];j++)
 486         {
 487                 at=pos_ins->subindex[i][j];
 488                 copy_rvec(r[at],r_ins[gctr]);
 489                 if( (r_ins[gctr][ZZ]<mem_p->zmax) && (r_ins[gctr][ZZ]>mem_p->zmin) )
 490                 {
 491                         rvec_inc(pos_ins->geom_cent[i],r_ins[gctr]);
 492                         c++;
 493                 }
 494                 else
 495                         outsidesum++;
 496                 gctr++;
 497         }
 498         if (c>0)
 499                 svmul(1/(double)c,pos_ins->geom_cent[i],pos_ins->geom_cent[i]);
 500         if (!bALLOW_ASYMMETRY)
 501                 pos_ins->geom_cent[i][ZZ]=mem_p->zmed;
 502
 503         fprintf(stderr,"Embedding piece %d with center of geometry: %f %f %f\n",i,pos_ins->geom_cent[i][XX],pos_ins->geom_cent[i][YY],pos_ins->geom_cent[i][ZZ]);
 504     }
 505     fprintf(stderr,"\n");
 506 }
 507
 508 void resize(t_block *ins_at, rvec *r_ins, rvec *r, pos_ins_t *pos_ins,rvec fac)
 509 {
 510         int i,j,k,at,c=0;
 511         for (k=0;k<pos_ins->pieces;k++)
 512                 for(i=0;i<pos_ins->nidx[k];i++)
 513                 {
 514                         at=pos_ins->subindex[k][i];
 515                         for(j=0;j<DIM;j++)
 516                                 r[at][j]=pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
 517                         c++;
 518                 }
 519 }
 520
 521 int gen_rm_list(rmm_t *rm_p,t_block *ins_at,t_block *rest_at,t_pbc *pbc, gmx_mtop_t *mtop,
 522                 rvec *r, rvec *r_ins, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad, int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
 523 {
 524         int i,j,k,l,at,at2,mol_id;
 525         int type=0,block=0;
 526         int nrm,nupper,nlower;
 527         real r_min_rad,z_lip,min_norm;
 528         gmx_bool bRM;
 529         rvec dr,dr_tmp;
 530         real *dist;
 531         int *order;
 532
 533         r_min_rad=probe_rad*probe_rad;
 534         snew(rm_p->mol,mtop->mols.nr);
 535         snew(rm_p->block,mtop->mols.nr);
 536         nrm=nupper=0;
 537         nlower=low_up_rm;
 538         for(i=0;i<ins_at->nr;i++)
 539         {
 540                 at=ins_at->index[i];
 541                 for(j=0;j<rest_at->nr;j++)
 542                 {
 543                         at2=rest_at->index[j];
 544                         pbc_dx(pbc,r[at],r[at2],dr);
 545
 546                         if(norm2(dr)<r_min_rad)
 547                         {
 548                                 mol_id = get_mol_id(at2,mtop->nmolblock,mtop->molblock,&type,&block);
 549                                 bRM=TRUE;
 550                                 for(l=0;l<nrm;l++)
 551                                         if(rm_p->mol[l]==mol_id)
 552                                                 bRM=FALSE;
 553                                 if(bRM)
 554                                 {
 555                                         /*fprintf(stderr,"%d wordt toegevoegd\n",mol_id);*/
 556                                         rm_p->mol[nrm]=mol_id;
 557                                         rm_p->block[nrm]=block;
 558                                         nrm++;
 559                                         z_lip=0.0;
 560                                         for(l=0;l<mem_p->nmol;l++)
 561                                         {
 562                                                 if(mol_id==mem_p->mol_id[l])
 563                                                 {
 564                                                         for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
 565                                                                 z_lip+=r[k][ZZ];
 566                                                         z_lip/=mtop->molblock[block].natoms_mol;
 567                                                         if(z_lip<mem_p->zmed)
 568                                                                 nlower++;
 569                                                         else
 570                                                                 nupper++;
 571                                                 }
 572                                         }
 573                                 }
 574                         }
 575                 }
 576         }
 577
 578         /*make sure equal number of lipids from upper and lower layer are removed */
 579         if( (nupper!=nlower) && (!bALLOW_ASYMMETRY) )
 580         {
 581                 snew(dist,mem_p->nmol);
 582                 snew(order,mem_p->nmol);
 583                 for(i=0;i<mem_p->nmol;i++)
 584                 {
 585                         at = mtop->mols.index[mem_p->mol_id[i]];
 586                         pbc_dx(pbc,r[at],pos_ins->geom_cent[0],dr);
 587                         if (pos_ins->pieces>1)
 588                         {
 589                                 /*minimum dr value*/
 590                                 min_norm=norm2(dr);
 591                                 for (k=1;k<pos_ins->pieces;k++)
 592                                 {
 593                                         pbc_dx(pbc,r[at],pos_ins->geom_cent[k],dr_tmp);
 594                                         if (norm2(dr_tmp) < min_norm)
 595                                         {
 596                                                 min_norm=norm2(dr_tmp);
 597                                                 copy_rvec(dr_tmp,dr);
 598                                         }
 599                                 }
 600                         }
 601                         dist[i]=dr[XX]*dr[XX]+dr[YY]*dr[YY];
 602                         j=i-1;
 603                         while (j>=0 && dist[i]<dist[order[j]])
 604                         {
 605                                 order[j+1]=order[j];
 606                                 j--;
 607                         }
 608                         order[j+1]=i;
 609                 }
 610
 611                 i=0;
 612                 while(nupper!=nlower)
 613                 {
 614                         mol_id=mem_p->mol_id[order[i]];
 615                         block=get_block(mol_id,mtop->nmolblock,mtop->molblock);
 616
 617                         bRM=TRUE;
 618                         for(l=0;l<nrm;l++)
 619                                 if(rm_p->mol[l]==mol_id)
 620                                         bRM=FALSE;
 621                         if(bRM)
 622                         {
 623                                 z_lip=0;
 624                                 for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
 625                                         z_lip+=r[k][ZZ];
 626                                 z_lip/=mtop->molblock[block].natoms_mol;
 627                                 if(nupper>nlower && z_lip<mem_p->zmed)
 628                                 {
 629                                         rm_p->mol[nrm]=mol_id;
 630                                         rm_p->block[nrm]=block;
 631                                         nrm++;
 632                                         nlower++;
 633                                 }
 634                                 else if (nupper<nlower && z_lip>mem_p->zmed)
 635                                 {
 636                                         rm_p->mol[nrm]=mol_id;
 637                                         rm_p->block[nrm]=block;
 638                                         nrm++;
 639                                         nupper++;
 640                                 }
 641                         }
 642                         i++;
 643
 644                         if(i>mem_p->nmol)
 645                                 gmx_fatal(FARGS,"Trying to remove more lipid molecules than there are in the membrane");
 646                 }
 647                 sfree(dist);
 648                 sfree(order);
 649         }
 650
 651         rm_p->nr=nrm;
 652         srenew(rm_p->mol,nrm);
 653         srenew(rm_p->block,nrm);
 654
 655         return nupper+nlower;
 656 }
 657
 658 void rm_group(t_inputrec *ir, gmx_groups_t *groups, gmx_mtop_t *mtop, rmm_t *rm_p, t_state *state, t_block *ins_at, pos_ins_t *pos_ins)
 659 {
 660         int i,j,k,n,rm,mol_id,at,block;
 661         rvec *x_tmp,*v_tmp;
 662         atom_id *list,*new_mols;
 663         unsigned char  *new_egrp[egcNR];
 664         gmx_bool bRM;
 665
 666         snew(list,state->natoms);
 667         n=0;
 668         for(i=0;i<rm_p->nr;i++)
 669         {
 670                 mol_id=rm_p->mol[i];
 671                 at=mtop->mols.index[mol_id];
 672                 block =rm_p->block[i];
 673                 mtop->molblock[block].nmol--;
 674                 for(j=0;j<mtop->molblock[block].natoms_mol;j++)
 675                 {
 676                         list[n]=at+j;
 677                         n++;
 678                 }
 679
 680                 mtop->mols.index[mol_id]=-1;
 681         }
 682
 683         mtop->mols.nr-=rm_p->nr;
 684         mtop->mols.nalloc_index-=rm_p->nr;
 685         snew(new_mols,mtop->mols.nr);
 686         for(i=0;i<mtop->mols.nr+rm_p->nr;i++)
 687         {
 688                 j=0;
 689                 if(mtop->mols.index[i]!=-1)
 690                 {
 691                         new_mols[j]=mtop->mols.index[i];
 692                         j++;
 693                 }
 694         }
 695         sfree(mtop->mols.index);
 696         mtop->mols.index=new_mols;
 697
 698
 699         mtop->natoms-=n;
 700         state->natoms-=n;
 701         state->nalloc=state->natoms;
 702         snew(x_tmp,state->nalloc);
 703         snew(v_tmp,state->nalloc);
 704
 705         for(i=0;i<egcNR;i++)
 706         {
 707                 if(groups->grpnr[i]!=NULL)
 708                 {
 709                         groups->ngrpnr[i]=state->natoms;
 710                         snew(new_egrp[i],state->natoms);
 711                 }
 712         }
 713
 714         rm=0;
 715         for (i=0;i<state->natoms+n;i++)
 716         {
 717                 bRM=FALSE;
 718                 for(j=0;j<n;j++)
 719                 {
 720                         if(i==list[j])
 721                         {
 722                                 bRM=TRUE;
 723                                 rm++;
 724                         }
 725                 }
 726
 727                 if(!bRM)
 728                 {
 729                         for(j=0;j<egcNR;j++)
 730                         {
 731                                 if(groups->grpnr[j]!=NULL)
 732                                 {
 733                                         new_egrp[j][i-rm]=groups->grpnr[j][i];
 734                                 }
 735                         }
 736                         copy_rvec(state->x[i],x_tmp[i-rm]);
 737                         copy_rvec(state->v[i],v_tmp[i-rm]);
 738                         for(j=0;j<ins_at->nr;j++)
 739                         {
 740                                 if (i==ins_at->index[j])
 741                                         ins_at->index[j]=i-rm;
 742                         }
 743                         for(j=0;j<pos_ins->pieces;j++)
 744                         {
 745                                 for(k=0;k<pos_ins->nidx[j];k++)
 746                                 {
 747                                         if (i==pos_ins->subindex[j][k])
 748                                                 pos_ins->subindex[j][k]=i-rm;
 749                                 }
 750                         }
 751                 }
 752         }
 753         sfree(state->x);
 754         state->x=x_tmp;
 755         sfree(state->v);
 756         state->v=v_tmp;
 757
 758         for(i=0;i<egcNR;i++)
 759         {
 760                 if(groups->grpnr[i]!=NULL)
 761                 {
 762                         sfree(groups->grpnr[i]);
 763                         groups->grpnr[i]=new_egrp[i];
 764                 }
 765         }
 766 }
 767
 768 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
 769 {
 770         int i,j,m;
 771         int type,natom,nmol,at,atom1=0,rm_at=0;
 772         gmx_bool *bRM,bINS;
 773         /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
 774         /*this routine does not live as dangerously as it seems. There is namely a check in mdrunner_membed to make
 775          *sure that g_membed exits with a warning when there are molecules of the same type not in the
 776          *ins_at index group. MGWolf 050710 */
 777
 778
 779         snew(bRM,mtop->nmoltype);
 780         for (i=0;i<mtop->nmoltype;i++)
 781         {
 782                 bRM[i]=TRUE;
 783         }
 784
 785         for (i=0;i<mtop->nmolblock;i++)
 786         {
 787             /*loop over molecule blocks*/
 788                 type        =mtop->molblock[i].type;
 789                 natom       =mtop->molblock[i].natoms_mol;
 790                 nmol            =mtop->molblock[i].nmol;
 791
 792                 for(j=0;j<natom*nmol && bRM[type]==TRUE;j++)
 793                 {
 794                     /*loop over atoms in the block*/
 795                         at=j+atom1; /*atom index = block index + offset*/
 796                         bINS=FALSE;
 797
 798                         for (m=0;(m<ins_at->nr) && (bINS==FALSE);m++)
 799                         {
 800                             /*loop over atoms in insertion index group to determine if we're inserting one*/
 801                                 if(at==ins_at->index[m])
 802                                 {
 803                                         bINS=TRUE;
 804                                 }
 805                         }
 806                         bRM[type]=bINS;
 807                 }
 808                 atom1+=natom*nmol; /*update offset*/
 809                 if(bRM[type])
 810                 {
 811                         rm_at+=natom*nmol; /*increment bonded removal counter by # atoms in block*/
 812                 }
 813         }
 814
 815         for(i=0;i<mtop->nmoltype;i++)
 816         {
 817                 if(bRM[i])
 818                 {
 819                         for(j=0;j<F_LJ;j++)
 820                         {
 821                                 mtop->moltype[i].ilist[j].nr=0;
 822                         }
 823                         for(j=F_POSRES;j<=F_VSITEN;j++)
 824                         {
 825                                 mtop->moltype[i].ilist[j].nr=0;
 826                         }
 827                 }
 828         }
 829         sfree(bRM);
 830
 831         return rm_at;
 832 }
 833
 834 void top_update(const char *topfile, char *ins, rmm_t *rm_p, gmx_mtop_t *mtop)
 835 {
 836 #define TEMP_FILENM "temp.top"
 837         int     bMolecules=0;
 838         FILE    *fpin,*fpout;
 839         char    buf[STRLEN],buf2[STRLEN],*temp;
 840         int             i,*nmol_rm,nmol,line;
 841
 842         fpin  = ffopen(topfile,"r");
 843         fpout = ffopen(TEMP_FILENM,"w");
 844
 845         snew(nmol_rm,mtop->nmoltype);
 846         for(i=0;i<rm_p->nr;i++)
 847                 nmol_rm[rm_p->block[i]]++;
 848
 849         line=0;
 850         while(fgets(buf,STRLEN,fpin))
 851         {
 852                 line++;
 853                 if(buf[0]!=';')
 854                 {
 855                         strcpy(buf2,buf);
 856                         if ((temp=strchr(buf2,'\n')) != NULL)
 857                                 temp[0]='\0';
 858                         ltrim(buf2);
 859
 860                         if (buf2[0]=='[')
 861                         {
 862                                 buf2[0]=' ';
 863                                 if ((temp=strchr(buf2,'\n')) != NULL)
 864                                         temp[0]='\0';
 865                                 rtrim(buf2);
 866                                 if (buf2[strlen(buf2)-1]==']')
 867                                 {
 868                                         buf2[strlen(buf2)-1]='\0';
 869                                         ltrim(buf2);
 870                                         rtrim(buf2);
 871                                         if (gmx_strcasecmp(buf2,"molecules")==0)
 872                                                 bMolecules=1;
 873                                 }
 874                                 fprintf(fpout,"%s",buf);
 875                         } else if (bMolecules==1)
 876                         {
 877                                 for(i=0;i<mtop->nmolblock;i++)
 878                                 {
 879                                         nmol=mtop->molblock[i].nmol;
 880                                         sprintf(buf,"%-15s %5d\n",*(mtop->moltype[mtop->molblock[i].type].name),nmol);
 881                                         fprintf(fpout,"%s",buf);
 882                                 }
 883                                 bMolecules=2;
 884                         } else if (bMolecules==2)
 885                         {
 886                                 /* print nothing */
 887                         } else
 888                         {
 889                                 fprintf(fpout,"%s",buf);
 890                         }
 891                 } else
 892                 {
 893                         fprintf(fpout,"%s",buf);
 894                 }
 895         }
 896
 897         fclose(fpout);
 898         /* use ffopen to generate backup of topinout */
 899         fpout=ffopen(topfile,"w");
 900         fclose(fpout);
 901         rename(TEMP_FILENM,topfile);
 902 #undef TEMP_FILENM
 903 }
 904
 905 double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 906              const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
 907              int nstglobalcomm,
 908              gmx_vsite_t *vsite,gmx_constr_t constr,
 909              int stepout,t_inputrec *ir,
 910              gmx_mtop_t *top_global,
 911              t_fcdata *fcd,
 912              t_state *state_global,
 913              t_mdatoms *mdatoms,
 914              t_nrnb *nrnb,gmx_wallcycle_t wcycle,
 915              gmx_edsam_t ed,t_forcerec *fr,
 916              int repl_ex_nst,int repl_ex_seed,
 917              real cpt_period,real max_hours,
 918              const char *deviceOptions,
 919              unsigned long Flags,
 920              gmx_runtime_t *runtime,
 921              rvec fac, rvec *r_ins, pos_ins_t *pos_ins, t_block *ins_at,
 922              real xy_step, real z_step, int it_xy, int it_z)
 923 {
 924     gmx_mdoutf_t *outf;
 925     gmx_large_int_t step,step_rel;
 926     double     run_time;
 927     double     t,t0,lam0;
 928     gmx_bool       bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
 929     gmx_bool       bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
 930                bFirstStep,bStateFromTPX,bInitStep,bLastStep,
 931                bBornRadii,bStartingFromCpt;
 932     gmx_bool       bDoDHDL=FALSE;
 933     gmx_bool       do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
 934                bForceUpdate=FALSE,bCPT;
 935     int        mdof_flags;
 936     gmx_bool       bMasterState;
 937     int        force_flags,cglo_flags;
 938     tensor     force_vir,shake_vir,total_vir,tmp_vir,pres;
 939     int        i,m;
 940     t_trxstatus *status;
 941     rvec       mu_tot;
 942     t_vcm      *vcm;
 943     t_state    *bufstate=NULL;
 944     matrix     *scale_tot,pcoupl_mu,M,ebox;
 945     gmx_nlheur_t nlh;
 946     t_trxframe rerun_fr;
 947 /*    gmx_repl_ex_t repl_ex=NULL;*/
 948     int        nchkpt=1;
 949
 950     gmx_localtop_t *top;
 951     t_mdebin *mdebin=NULL;
 952     t_state    *state=NULL;
 953     rvec       *f_global=NULL;
 954     int        n_xtc=-1;
 955     rvec       *x_xtc=NULL;
 956     gmx_enerdata_t *enerd;
 957     rvec       *f=NULL;
 958     gmx_global_stat_t gstat;
 959     gmx_update_t upd=NULL;
 960     t_graph    *graph=NULL;
 961     globsig_t   gs;
 962
 963     gmx_bool        bFFscan;
 964     gmx_groups_t *groups;
 965     gmx_ekindata_t *ekind, *ekind_save;
 966     gmx_shellfc_t shellfc;
 967     int         count,nconverged=0;
 968     real        timestep=0;
 969     double      tcount=0;
 970     gmx_bool        bIonize=FALSE;
 971     gmx_bool        bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
 972     gmx_bool        bAppend;
 973     gmx_bool        bResetCountersHalfMaxH=FALSE;
 974     gmx_bool        bVV,bIterations,bIterate,bFirstIterate,bTemp,bPres,bTrotter;
 975     real        temp0,dvdl;
 976     int         a0,a1,ii;
 977     rvec        *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
 978     matrix      boxcopy={{0}},lastbox;
 979         real        veta_save,pcurr,scalevir,tracevir;
 980         real        vetanew = 0;
 981     double      cycles;
 982         real        last_conserved = 0;
 983     real        last_ekin = 0;
 984         t_extmass   MassQ;
 985     int         **trotter_seq;
 986     char        sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
 987     int         handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
 988     gmx_iterate_t iterate;
 989 #ifdef GMX_FAHCORE
 990     /* Temporary addition for FAHCORE checkpointing */
 991     int chkpt_ret;
 992 #endif
 993
 994     /* Check for special mdrun options */
 995     bRerunMD = (Flags & MD_RERUN);
 996     bIonize  = (Flags & MD_IONIZE);
 997     bFFscan  = (Flags & MD_FFSCAN);
 998     bAppend  = (Flags & MD_APPENDFILES);
 999     bGStatEveryStep = FALSE;
1000     if (Flags & MD_RESETCOUNTERSHALFWAY)
1001     {
1002         if (ir->nsteps > 0)
1003         {
1004             /* Signal to reset the counters half the simulation steps. */
1005             wcycle_set_reset_counters(wcycle,ir->nsteps/2);
1006         }
1007         /* Signal to reset the counters halfway the simulation time. */
1008         bResetCountersHalfMaxH = (max_hours > 0);
1009     }
1010
1011     /* md-vv uses averaged full step velocities for T-control
1012        md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1013        md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1014     bVV = EI_VV(ir->eI);
1015     if (bVV) /* to store the initial velocities while computing virial */
1016     {
1017         snew(cbuf,top_global->natoms);
1018     }
1019     /* all the iteratative cases - only if there are constraints */
1020     bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
1021     bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
1022
1023     if (bRerunMD)
1024     {
1025         /* Since we don't know if the frames read are related in any way,
1026          * rebuild the neighborlist at every step.
1027          */
1028         ir->nstlist       = 1;
1029         ir->nstcalcenergy = 1;
1030         nstglobalcomm     = 1;
1031     }
1032
1033     check_ir_old_tpx_versions(cr,fplog,ir,top_global);
1034
1035     nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
1036     /*bGStatEveryStep = (nstglobalcomm == 1);*/
1037     bGStatEveryStep = FALSE;
1038
1039     if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
1040     {
1041         fprintf(fplog,
1042                 "To reduce the energy communication with nstlist = -1\n"
1043                 "the neighbor list validity should not be checked at every step,\n"
1044                 "this means that exact integration is not guaranteed.\n"
1045                 "The neighbor list validity is checked after:\n"
1046                 "  <n.list life time> - 2*std.dev.(n.list life time)  steps.\n"
1047                 "In most cases this will result in exact integration.\n"
1048                 "This reduces the energy communication by a factor of 2 to 3.\n"
1049                 "If you want less energy communication, set nstlist > 3.\n\n");
1050     }
1051
1052     if (bRerunMD || bFFscan)
1053     {
1054         ir->nstxtcout = 0;
1055     }
1056     groups = &top_global->groups;
1057
1058     /* Initial values */
1059     init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
1060             nrnb,top_global,&upd,
1061             nfile,fnm,&outf,&mdebin,
1062             force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
1063
1064     clear_mat(total_vir);
1065     clear_mat(pres);
1066     /* Energy terms and groups */
1067     snew(enerd,1);
1068     init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
1069     if (DOMAINDECOMP(cr))
1070     {
1071         f = NULL;
1072     }
1073     else
1074     {
1075         snew(f,top_global->natoms);
1076     }
1077
1078     /* Kinetic energy data */
1079     snew(ekind,1);
1080     init_ekindata(fplog,top_global,&(ir->opts),ekind);
1081     /* needed for iteration of constraints */
1082     snew(ekind_save,1);
1083     init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
1084     /* Copy the cos acceleration to the groups struct */
1085     ekind->cosacc.cos_accel = ir->cos_accel;
1086
1087     gstat = global_stat_init(ir);
1088     debug_gmx();
1089
1090     /* Check for polarizable models and flexible constraints */
1091     shellfc = init_shell_flexcon(fplog,
1092                                  top_global,n_flexible_constraints(constr),
1093                                  (ir->bContinuation ||
1094                                   (DOMAINDECOMP(cr) && !MASTER(cr))) ?
1095                                  NULL : state_global->x);
1096
1097 /*    if (DEFORM(*ir))
1098     {
1099 #ifdef GMX_THREADS
1100         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
1101 #endif
1102         set_deform_reference_box(upd,
1103                                  deform_init_init_step_tpx,
1104                                  deform_init_box_tpx);
1105 #ifdef GMX_THREADS
1106         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
1107 #endif
1108     }*/
1109
1110 /*    {
1111         double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
1112         if ((io > 2000) && MASTER(cr))
1113             fprintf(stderr,
1114                     "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1115                     io);
1116     }*/
1117
1118     if (DOMAINDECOMP(cr)) {
1119         top = dd_init_local_top(top_global);
1120
1121         snew(state,1);
1122         dd_init_local_state(cr->dd,state_global,state);
1123
1124         if (DDMASTER(cr->dd) && ir->nstfout) {
1125             snew(f_global,state_global->natoms);
1126         }
1127     } else {
1128         if (PAR(cr)) {
1129             /* Initialize the particle decomposition and split the topology */
1130             top = split_system(fplog,top_global,ir,cr);
1131
1132             pd_cg_range(cr,&fr->cg0,&fr->hcg);
1133             pd_at_range(cr,&a0,&a1);
1134         } else {
1135             top = gmx_mtop_generate_local_top(top_global,ir);
1136
1137             a0 = 0;
1138             a1 = top_global->natoms;
1139         }
1140
1141         state = partdec_init_local_state(cr,state_global);
1142         f_global = f;
1143
1144         atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
1145
1146         if (vsite) {
1147             set_vsite_top(vsite,top,mdatoms,cr);
1148         }
1149
1150         if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
1151             graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
1152         }
1153
1154         if (shellfc) {
1155             make_local_shells(cr,mdatoms,shellfc);
1156         }
1157
1158         if (ir->pull && PAR(cr)) {
1159             dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
1160         }
1161     }
1162
1163     if (DOMAINDECOMP(cr))
1164     {
1165         /* Distribute the charge groups over the nodes from the master node */
1166         dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
1167                             state_global,top_global,ir,
1168                             state,&f,mdatoms,top,fr,
1169                             vsite,shellfc,constr,
1170                             nrnb,wcycle,FALSE);
1171     }
1172
1173     update_mdatoms(mdatoms,state->lambda);
1174
1175     if (MASTER(cr))
1176     {
1177         if (opt2bSet("-cpi",nfile,fnm))
1178         {
1179             /* Update mdebin with energy history if appending to output files */
1180             if ( Flags & MD_APPENDFILES )
1181             {
1182                 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
1183             }
1184             else
1185             {
1186                 /* We might have read an energy history from checkpoint,
1187                  * free the allocated memory and reset the counts.
1188                  */
1189                 done_energyhistory(&state_global->enerhist);
1190                 init_energyhistory(&state_global->enerhist);
1191             }
1192         }
1193         /* Set the initial energy history in state by updating once */
1194         update_energyhistory(&state_global->enerhist,mdebin);
1195     }
1196
1197     if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
1198         /* Set the random state if we read a checkpoint file */
1199         set_stochd_state(upd,state);
1200     }
1201
1202     /* Initialize constraints */
1203     if (constr) {
1204         if (!DOMAINDECOMP(cr))
1205             set_constraints(constr,top,ir,mdatoms,cr);
1206     }
1207
1208     /* Check whether we have to GCT stuff */
1209  /*   bTCR = ftp2bSet(efGCT,nfile,fnm);
1210     if (bTCR) {
1211         if (MASTER(cr)) {
1212             fprintf(stderr,"Will do General Coupling Theory!\n");
1213         }
1214         gnx = top_global->mols.nr;
1215         snew(grpindex,gnx);
1216         for(i=0; (i<gnx); i++) {
1217             grpindex[i] = i;
1218         }
1219     }*/
1220
1221 /*    if (repl_ex_nst > 0 && MASTER(cr))
1222         repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
1223                                         repl_ex_nst,repl_ex_seed);*/
1224
1225     if (!ir->bContinuation && !bRerunMD)
1226     {
1227         if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
1228         {
1229             /* Set the velocities of frozen particles to zero */
1230             for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
1231             {
1232                 for(m=0; m<DIM; m++)
1233                 {
1234                     if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
1235                     {
1236                         state->v[i][m] = 0;
1237                     }
1238                 }
1239             }
1240         }
1241
1242         if (constr)
1243         {
1244             /* Constrain the initial coordinates and velocities */
1245             do_constrain_first(fplog,constr,ir,mdatoms,state,f,
1246                                graph,cr,nrnb,fr,top,shake_vir);
1247         }
1248         if (vsite)
1249         {
1250             /* Construct the virtual sites for the initial configuration */
1251             construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
1252                              top->idef.iparams,top->idef.il,
1253                              fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1254         }
1255     }
1256
1257     debug_gmx();
1258
1259     /* I'm assuming we need global communication the first time! MRS */
1260     cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
1261                   | (bVV ? CGLO_PRESSURE:0)
1262                   | (bVV ? CGLO_CONSTRAINT:0)
1263                   | (bRerunMD ? CGLO_RERUNMD:0)
1264                   | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
1265
1266     bSumEkinhOld = FALSE;
1267     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1268                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1269                     constr,NULL,FALSE,state->box,
1270                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
1271     if (ir->eI == eiVVAK) {
1272         /* a second call to get the half step temperature initialized as well */
1273         /* we do the same call as above, but turn the pressure off -- internally, this
1274            is recognized as a velocity verlet half-step kinetic energy calculation.
1275            This minimized excess variables, but perhaps loses some logic?*/
1276
1277         compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1278                         wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1279                         constr,NULL,FALSE,state->box,
1280                         top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1281                         cglo_flags &~ CGLO_PRESSURE);
1282     }
1283
1284     /* Calculate the initial half step temperature, and save the ekinh_old */
1285     if (!(Flags & MD_STARTFROMCPT))
1286     {
1287         for(i=0; (i<ir->opts.ngtc); i++)
1288         {
1289             copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
1290         }
1291     }
1292     if (ir->eI != eiVV)
1293     {
1294         enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
1295                                      and there is no previous step */
1296     }
1297     temp0 = enerd->term[F_TEMP];
1298
1299     /* if using an iterative algorithm, we need to create a working directory for the state. */
1300     if (bIterations)
1301     {
1302             bufstate = init_bufstate(state);
1303     }
1304     if (bFFscan)
1305     {
1306         snew(xcopy,state->natoms);
1307         snew(vcopy,state->natoms);
1308         copy_rvecn(state->x,xcopy,0,state->natoms);
1309         copy_rvecn(state->v,vcopy,0,state->natoms);
1310         copy_mat(state->box,boxcopy);
1311     }
1312
1313     /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
1314        temperature control */
1315     trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
1316
1317     if (MASTER(cr))
1318     {
1319         if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
1320         {
1321             fprintf(fplog,
1322                     "RMS relative constraint deviation after constraining: %.2e\n",
1323                     constr_rmsd(constr,FALSE));
1324         }
1325         fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
1326         if (bRerunMD)
1327         {
1328             fprintf(stderr,"starting md rerun '%s', reading coordinates from"
1329                     " input trajectory '%s'\n\n",
1330                     *(top_global->name),opt2fn("-rerun",nfile,fnm));
1331             if (bVerbose)
1332             {
1333                 fprintf(stderr,"Calculated time to finish depends on nsteps from "
1334                         "run input file,\nwhich may not correspond to the time "
1335                         "needed to process input trajectory.\n\n");
1336             }
1337         }
1338         else
1339         {
1340             char tbuf[20];
1341             fprintf(stderr,"starting mdrun '%s'\n",
1342                     *(top_global->name));
1343             if (ir->nsteps >= 0)
1344             {
1345                 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
1346             }
1347             else
1348             {
1349                 sprintf(tbuf,"%s","infinite");
1350             }
1351             if (ir->init_step > 0)
1352             {
1353                 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
1354                         gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
1355                         gmx_step_str(ir->init_step,sbuf2),
1356                         ir->init_step*ir->delta_t);
1357             }
1358             else
1359             {
1360                 fprintf(stderr,"%s steps, %s ps.\n",
1361                         gmx_step_str(ir->nsteps,sbuf),tbuf);
1362             }
1363         }
1364         fprintf(fplog,"\n");
1365     }
1366
1367     /* Set and write start time */
1368     runtime_start(runtime);
1369     print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
1370     wallcycle_start(wcycle,ewcRUN);
1371     if (fplog)
1372         fprintf(fplog,"\n");
1373
1374     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
1375 /*#ifdef GMX_FAHCORE
1376     chkpt_ret=fcCheckPointParallel( cr->nodeid,
1377                                     NULL,0);
1378     if ( chkpt_ret == 0 )
1379         gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
1380 #endif*/
1381
1382     debug_gmx();
1383     /***********************************************************
1384      *
1385      *             Loop over MD steps
1386      *
1387      ************************************************************/
1388
1389     /* if rerunMD then read coordinates and velocities from input trajectory */
1390     if (bRerunMD)
1391     {
1392         if (getenv("GMX_FORCE_UPDATE"))
1393         {
1394             bForceUpdate = TRUE;
1395         }
1396
1397         bNotLastFrame = read_first_frame(oenv,&status,
1398                                          opt2fn("-rerun",nfile,fnm),
1399                                          &rerun_fr,TRX_NEED_X | TRX_READ_V);
1400         if (rerun_fr.natoms != top_global->natoms)
1401         {
1402             gmx_fatal(FARGS,
1403                       "Number of atoms in trajectory (%d) does not match the "
1404                       "run input file (%d)\n",
1405                       rerun_fr.natoms,top_global->natoms);
1406         }
1407         if (ir->ePBC != epbcNONE)
1408         {
1409             if (!rerun_fr.bBox)
1410             {
1411                 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f does not contain a box, while pbc is used",rerun_fr.step,rerun_fr.time);
1412             }
1413             if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
1414             {
1415                 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
1416             }
1417
1418             /* Set the shift vectors.
1419              * Necessary here when have a static box different from the tpr box.
1420              */
1421             calc_shifts(rerun_fr.box,fr->shift_vec);
1422         }
1423     }
1424
1425     /* loop over MD steps or if rerunMD to end of input trajectory */
1426     bFirstStep = TRUE;
1427     /* Skip the first Nose-Hoover integration when we get the state from tpx */
1428     bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
1429     bInitStep = bFirstStep && (bStateFromTPX || bVV);
1430     bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
1431     bLastStep    = FALSE;
1432     bSumEkinhOld = FALSE;
1433     bExchanged   = FALSE;
1434
1435     init_global_signals(&gs,cr,ir,repl_ex_nst);
1436
1437     step = ir->init_step;
1438     step_rel = 0;
1439
1440     if (ir->nstlist == -1)
1441     {
1442         init_nlistheuristics(&nlh,bGStatEveryStep,step);
1443     }
1444
1445     bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
1446     while (!bLastStep || (bRerunMD && bNotLastFrame)) {
1447
1448         wallcycle_start(wcycle,ewcSTEP);
1449
1450         GMX_MPE_LOG(ev_timestep1);
1451
1452         if (bRerunMD) {
1453             if (rerun_fr.bStep) {
1454                 step = rerun_fr.step;
1455                 step_rel = step - ir->init_step;
1456             }
1457             if (rerun_fr.bTime) {
1458                 t = rerun_fr.time;
1459             }
1460             else
1461             {
1462                 t = step;
1463             }
1464         }
1465         else
1466         {
1467             bLastStep = (step_rel == ir->nsteps);
1468             t = t0 + step*ir->delta_t;
1469         }
1470
1471         if (ir->efep != efepNO)
1472         {
1473             if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
1474             {
1475                 state_global->lambda = rerun_fr.lambda;
1476             }
1477             else
1478             {
1479                 state_global->lambda = lam0 + step*ir->delta_lambda;
1480             }
1481             state->lambda = state_global->lambda;
1482             bDoDHDL = do_per_step(step,ir->nstdhdl);
1483         }
1484
1485         if (bSimAnn)
1486         {
1487             update_annealing_target_temp(&(ir->opts),t);
1488         }
1489
1490         if (bRerunMD)
1491         {
1492             if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
1493             {
1494                 for(i=0; i<state_global->natoms; i++)
1495                 {
1496                     copy_rvec(rerun_fr.x[i],state_global->x[i]);
1497                 }
1498                 if (rerun_fr.bV)
1499                 {
1500                     for(i=0; i<state_global->natoms; i++)
1501                     {
1502                         copy_rvec(rerun_fr.v[i],state_global->v[i]);
1503                     }
1504                 }
1505                 else
1506                 {
1507                     for(i=0; i<state_global->natoms; i++)
1508                     {
1509                         clear_rvec(state_global->v[i]);
1510                     }
1511                     if (bRerunWarnNoV)
1512                     {
1513                         fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
1514                                 "         Ekin, temperature and pressure are incorrect,\n"
1515                                 "         the virial will be incorrect when constraints are present.\n"
1516                                 "\n");
1517                         bRerunWarnNoV = FALSE;
1518                     }
1519                 }
1520             }
1521             copy_mat(rerun_fr.box,state_global->box);
1522             copy_mat(state_global->box,state->box);
1523
1524             if (vsite && (Flags & MD_RERUN_VSITE))
1525             {
1526                 if (DOMAINDECOMP(cr))
1527                 {
1528                     gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1529                 }
1530                 if (graph)
1531                 {
1532                     /* Following is necessary because the graph may get out of sync
1533                      * with the coordinates if we only have every N'th coordinate set
1534                      */
1535                     mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
1536                     shift_self(graph,state->box,state->x);
1537                 }
1538                 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1539                                  top->idef.iparams,top->idef.il,
1540                                  fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1541                 if (graph)
1542                 {
1543                     unshift_self(graph,state->box,state->x);
1544                 }
1545             }
1546         }
1547
1548         /* Stop Center of Mass motion */
1549         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
1550
1551         /* Copy back starting coordinates in case we're doing a forcefield scan */
1552         if (bFFscan)
1553         {
1554             for(ii=0; (ii<state->natoms); ii++)
1555             {
1556                 copy_rvec(xcopy[ii],state->x[ii]);
1557                 copy_rvec(vcopy[ii],state->v[ii]);
1558             }
1559             copy_mat(boxcopy,state->box);
1560         }
1561
1562         if (bRerunMD)
1563         {
1564             /* for rerun MD always do Neighbour Searching */
1565             bNS = (bFirstStep || ir->nstlist != 0);
1566             bNStList = bNS;
1567         }
1568         else
1569         {
1570             /* Determine whether or not to do Neighbour Searching and LR */
1571             bNStList = (ir->nstlist > 0  && step % ir->nstlist == 0);
1572
1573             bNS = (bFirstStep || bExchanged || bNStList ||
1574                    (ir->nstlist == -1 && nlh.nabnsb > 0));
1575
1576             if (bNS && ir->nstlist == -1)
1577             {
1578                 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
1579             }
1580         }
1581
1582         /* < 0 means stop at next step, > 0 means stop at next NS step */
1583         if ( (gs.set[eglsSTOPCOND] < 0 ) ||
1584              ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
1585         {
1586             bLastStep = TRUE;
1587         }
1588
1589         /* Determine whether or not to update the Born radii if doing GB */
1590         bBornRadii=bFirstStep;
1591         if (ir->implicit_solvent && (step % ir->nstgbradii==0))
1592         {
1593             bBornRadii=TRUE;
1594         }
1595
1596         do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
1597         do_verbose = bVerbose &&
1598                   (step % stepout == 0 || bFirstStep || bLastStep);
1599
1600         if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
1601         {
1602             if (bRerunMD)
1603             {
1604                 bMasterState = TRUE;
1605             }
1606             else
1607             {
1608                 bMasterState = FALSE;
1609                 /* Correct the new box if it is too skewed */
1610                 if (DYNAMIC_BOX(*ir))
1611                 {
1612                     if (correct_box(fplog,step,state->box,graph))
1613                     {
1614                         bMasterState = TRUE;
1615                     }
1616                 }
1617                 if (DOMAINDECOMP(cr) && bMasterState)
1618                 {
1619                     dd_collect_state(cr->dd,state,state_global);
1620                 }
1621             }
1622
1623             if (DOMAINDECOMP(cr))
1624             {
1625                 /* Repartition the domain decomposition */
1626                 wallcycle_start(wcycle,ewcDOMDEC);
1627                 dd_partition_system(fplog,step,cr,
1628                                     bMasterState,nstglobalcomm,
1629                                     state_global,top_global,ir,
1630                                     state,&f,mdatoms,top,fr,
1631                                     vsite,shellfc,constr,
1632                                     nrnb,wcycle,do_verbose);
1633                 wallcycle_stop(wcycle,ewcDOMDEC);
1634                 /* If using an iterative integrator, reallocate space to match the decomposition */
1635             }
1636         }
1637
1638         if (MASTER(cr) && do_log && !bFFscan)
1639         {
1640             print_ebin_header(fplog,step,t,state->lambda);
1641         }
1642
1643         if (ir->efep != efepNO)
1644         {
1645             update_mdatoms(mdatoms,state->lambda);
1646         }
1647
1648         if (bRerunMD && rerun_fr.bV)
1649         {
1650
1651             /* We need the kinetic energy at minus the half step for determining
1652              * the full step kinetic energy and possibly for T-coupling.*/
1653             /* This may not be quite working correctly yet . . . . */
1654             compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1655                             wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
1656                             constr,NULL,FALSE,state->box,
1657                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1658                             CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1659         }
1660         clear_mat(force_vir);
1661
1662         /* Ionize the atoms if necessary */
1663 /*        if (bIonize)
1664         {
1665             ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
1666                    mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
1667         }*/
1668
1669         /* Update force field in ffscan program */
1670 /*        if (bFFscan)
1671         {
1672             if (update_forcefield(fplog,
1673                                   nfile,fnm,fr,
1674                                   mdatoms->nr,state->x,state->box)) {
1675                 if (gmx_parallel_env_initialized())
1676                 {
1677                     gmx_finalize();
1678                 }
1679                 exit(0);
1680             }
1681         }*/
1682
1683         GMX_MPE_LOG(ev_timestep2);
1684
1685         /* We write a checkpoint at this MD step when:
1686          * either at an NS step when we signalled through gs,
1687          * or at the last step (but not when we do not want confout),
1688          * but never at the first step or with rerun.
1689          */
1690 /*        bCPT = (((gs.set[eglsCHKPT] && bNS) ||
1691                  (bLastStep && (Flags & MD_CONFOUT))) &&
1692                 step > ir->init_step && !bRerunMD);
1693         if (bCPT)
1694         {
1695             gs.set[eglsCHKPT] = 0;
1696         }*/
1697
1698         /* Determine the energy and pressure:
1699          * at nstcalcenergy steps and at energy output steps (set below).
1700          */
1701         bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
1702         bCalcEnerPres = bNstEner;
1703
1704         /* Do we need global communication ? */
1705         bGStat = (bCalcEnerPres || bStopCM ||
1706                   (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1707
1708         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1709
1710         if (do_ene || do_log)
1711         {
1712             bCalcEnerPres = TRUE;
1713             bGStat    = TRUE;
1714         }
1715
1716         /* these CGLO_ options remain the same throughout the iteration */
1717         cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1718                       (bStopCM ? CGLO_STOPCM : 0) |
1719                       (bGStat ? CGLO_GSTAT : 0)
1720             );
1721
1722         force_flags = (GMX_FORCE_STATECHANGED |
1723                        ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1724                        GMX_FORCE_ALLFORCES |
1725                        (bNStList ? GMX_FORCE_DOLR : 0) |
1726                        GMX_FORCE_SEPLRF |
1727                        (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
1728                        (bDoDHDL ? GMX_FORCE_DHDL : 0)
1729             );
1730
1731         if (shellfc)
1732         {
1733             /* Now is the time to relax the shells */
1734             count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1735                                       ir,bNS,force_flags,
1736                                       bStopCM,top,top_global,
1737                                       constr,enerd,fcd,
1738                                       state,f,force_vir,mdatoms,
1739                                       nrnb,wcycle,graph,groups,
1740                                       shellfc,fr,bBornRadii,t,mu_tot,
1741                                       state->natoms,&bConverged,vsite,
1742                                       outf->fp_field);
1743             tcount+=count;
1744
1745             if (bConverged)
1746             {
1747                 nconverged++;
1748             }
1749         }
1750         else
1751         {
1752             /* The coordinates (x) are shifted (to get whole molecules)
1753              * in do_force.
1754              * This is parallellized as well, and does communication too.
1755              * Check comments in sim_util.c
1756              */
1757
1758             do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1759                      state->box,state->x,&state->hist,
1760                      f,force_vir,mdatoms,enerd,fcd,
1761                      state->lambda,graph,
1762                      fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1763                      (bNS ? GMX_FORCE_NS : 0) | force_flags);
1764         }
1765
1766         GMX_BARRIER(cr->mpi_comm_mygroup);
1767
1768  /*       if (bTCR)
1769         {
1770             mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1771                                    mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1772         }
1773
1774         if (bTCR && bFirstStep)
1775         {
1776             tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1777             fprintf(fplog,"Done init_coupling\n");
1778             fflush(fplog);
1779         }*/
1780
1781         /*  ############### START FIRST UPDATE HALF-STEP ############### */
1782
1783         if (bVV && !bStartingFromCpt && !bRerunMD)
1784         {
1785             if (ir->eI == eiVV)
1786             {
1787                 if (bInitStep)
1788                 {
1789                     /* if using velocity verlet with full time step Ekin,
1790                      * take the first half step only to compute the
1791                      * virial for the first step. From there,
1792                      * revert back to the initial coordinates
1793                      * so that the input is actually the initial step.
1794                      */
1795                     copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1796                 }
1797
1798                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1799                 if (!bInitStep)
1800                 {
1801                   trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1802                 }
1803
1804                 if (ir->eI == eiVVAK)
1805                 {
1806                   update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1807                 }
1808
1809                 update_coords(fplog,step,ir,mdatoms,state,
1810                               f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1811                               ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1812                               cr,nrnb,constr,&top->idef);
1813
1814                 if (bIterations)
1815                 {
1816                     gmx_iterate_init(&iterate,bIterations && !bInitStep);
1817                 }
1818                 /* for iterations, we save these vectors, as we will be self-consistently iterating
1819                    the calculations */
1820                 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1821
1822                 /* save the state */
1823                 if (bIterations && iterate.bIterate) {
1824                     copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1825                 }
1826             }
1827
1828             bFirstIterate = TRUE;
1829             while (bFirstIterate || (bIterations && iterate.bIterate))
1830             {
1831                 if (bIterations && iterate.bIterate)
1832                 {
1833                     copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1834                     if (bFirstIterate && bTrotter)
1835                     {
1836                         /* The first time through, we need a decent first estimate
1837                            of veta(t+dt) to compute the constraints.  Do
1838                            this by computing the box volume part of the
1839                            trotter integration at this time. Nothing else
1840                            should be changed by this routine here.  If
1841                            !(first time), we start with the previous value
1842                            of veta.  */
1843
1844                         veta_save = state->veta;
1845                         trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1846                         vetanew = state->veta;
1847                         state->veta = veta_save;
1848                     }
1849                 }
1850
1851                 bOK = TRUE;
1852                 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) {  /* Why is rerun_fr.bV here?  Unclear. */
1853                     dvdl = 0;
1854
1855                     update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1856                                        &top->idef,shake_vir,NULL,
1857                                        cr,nrnb,wcycle,upd,constr,
1858                                        bInitStep,TRUE,bCalcEnerPres,vetanew);
1859
1860                     if (!bOK && !bFFscan)
1861                     {
1862                         gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1863                     }
1864
1865                 }
1866                 else if (graph)
1867                 { /* Need to unshift here if a do_force has been
1868                      called in the previous step */
1869                     unshift_self(graph,state->box,state->x);
1870                 }
1871
1872
1873                 if (bVV) {
1874                     /* if VV, compute the pressure and constraints */
1875                     /* if VV2, the pressure and constraints only if using pressure control.*/
1876                     bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir));
1877                     bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));
1878                     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1879                                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1880                                     constr,NULL,FALSE,state->box,
1881                                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1882                                     cglo_flags
1883                                     | CGLO_ENERGY
1884                                     | (bTemp ? CGLO_TEMPERATURE:0)
1885                                     | (bPres ? CGLO_PRESSURE : 0)
1886                                     | (bPres ? CGLO_CONSTRAINT : 0)
1887                                     | (iterate.bIterate ? CGLO_ITERATE : 0)
1888                                     | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1889                                     | CGLO_SCALEEKIN
1890                         );
1891                 }
1892                 /* explanation of above:
1893                    a) We compute Ekin at the full time step
1894                    if 1) we are using the AveVel Ekin, and it's not the
1895                    initial step, or 2) if we are using AveEkin, but need the full
1896                    time step kinetic energy for the pressure.
1897                    b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1898                    EkinAveVel because it's needed for the pressure */
1899
1900                 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1901                 if (bVV && !bInitStep)
1902                 {
1903                   trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ2);
1904                 }
1905
1906                 if (bIterations &&
1907                     done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1908                                    state->veta,&vetanew))
1909                 {
1910                     break;
1911                 }
1912                 bFirstIterate = FALSE;
1913             }
1914
1915             if (bTrotter && !bInitStep) {
1916                 copy_mat(shake_vir,state->svir_prev);
1917                 copy_mat(force_vir,state->fvir_prev);
1918                 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1919                     /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1920                     enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1921                     enerd->term[F_EKIN] = trace(ekind->ekin);
1922                 }
1923             }
1924             /* if it's the initial step, we performed this first step just to get the constraint virial */
1925             if (bInitStep && ir->eI==eiVV) {
1926                 copy_rvecn(cbuf,state->v,0,state->natoms);
1927             }
1928
1929             if (fr->bSepDVDL && fplog && do_log)
1930             {
1931                 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1932             }
1933             enerd->term[F_DHDL_CON] += dvdl;
1934
1935             GMX_MPE_LOG(ev_timestep1);
1936
1937         }
1938
1939         /* MRS -- now done iterating -- compute the conserved quantity */
1940         if (bVV) {
1941             last_conserved = 0;
1942             if (IR_NVT_TROTTER(ir) || IR_NPT_TROTTER(ir))
1943             {
1944                 last_conserved =
1945                     NPT_energy(ir,state,&MassQ);
1946                 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1947                 {
1948                     last_conserved -= enerd->term[F_DISPCORR];
1949                 }
1950             }
1951             if (ir->eI==eiVV) {
1952                 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
1953             }
1954         }
1955
1956         /* ########  END FIRST UPDATE STEP  ############## */
1957         /* ########  If doing VV, we now have v(dt) ###### */
1958
1959         /* ################## START TRAJECTORY OUTPUT ################# */
1960
1961         /* Now we have the energies and forces corresponding to the
1962          * coordinates at time t. We must output all of this before
1963          * the update.
1964          * for RerunMD t is read from input trajectory
1965          */
1966         GMX_MPE_LOG(ev_output_start);
1967
1968         mdof_flags = 0;
1969         if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
1970         if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
1971         if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
1972         if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
1973 /*        if (bCPT) { mdof_flags |= MDOF_CPT; };*/
1974
1975 #ifdef GMX_FAHCORE
1976         if (MASTER(cr))
1977             fcReportProgress( ir->nsteps, step );
1978
1979         if (bLastStep)
1980         {
1981             /* Enforce writing positions and velocities at end of run */
1982             mdof_flags |= (MDOF_X | MDOF_V);
1983         }
1984             /* sync bCPT and fc record-keeping */
1985 /*            if (bCPT && MASTER(cr))
1986                 fcRequestCheckPoint();*/
1987 #endif
1988
1989         if (mdof_flags != 0)
1990         {
1991             wallcycle_start(wcycle,ewcTRAJ);
1992 /*            if (bCPT)
1993             {
1994                 if (state->flags & (1<<estLD_RNG))
1995                 {
1996                     get_stochd_state(upd,state);
1997                 }
1998                 if (MASTER(cr))
1999                 {
2000                     if (bSumEkinhOld)
2001                     {
2002                         state_global->ekinstate.bUpToDate = FALSE;
2003                     }
2004                     else
2005                     {
2006                         update_ekinstate(&state_global->ekinstate,ekind);
2007                         state_global->ekinstate.bUpToDate = TRUE;
2008                     }
2009                     update_energyhistory(&state_global->enerhist,mdebin);
2010                 }
2011             }*/
2012             write_traj(fplog,cr,outf,mdof_flags,top_global,
2013                        step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2014 /*            if (bCPT)
2015             {
2016                 nchkpt++;
2017                 bCPT = FALSE;
2018             }*/
2019             debug_gmx();
2020             if (bLastStep && step_rel == ir->nsteps &&
2021                 (Flags & MD_CONFOUT) && MASTER(cr) &&
2022                 !bRerunMD && !bFFscan)
2023             {
2024                 /* x and v have been collected in write_traj,
2025                  * because a checkpoint file will always be written
2026                  * at the last step.
2027                  */
2028                 fprintf(stderr,"\nWriting final coordinates.\n");
2029                 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2030                     DOMAINDECOMP(cr))
2031                 {
2032                     /* Make molecules whole only for confout writing */
2033                     do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2034                 }
2035 /*                write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2036                                     *top_global->name,top_global,
2037                                     state_global->x,state_global->v,
2038                                     ir->ePBC,state->box);*/
2039                 debug_gmx();
2040             }
2041             wallcycle_stop(wcycle,ewcTRAJ);
2042         }
2043         GMX_MPE_LOG(ev_output_finish);
2044
2045         /* kludge -- virial is lost with restart for NPT control. Must restart */
2046         if (bStartingFromCpt && bVV)
2047         {
2048             copy_mat(state->svir_prev,shake_vir);
2049             copy_mat(state->fvir_prev,force_vir);
2050         }
2051         /*  ################## END TRAJECTORY OUTPUT ################ */
2052
2053         /* Determine the pressure:
2054          * always when we want exact averages in the energy file,
2055          * at ns steps when we have pressure coupling,
2056          * otherwise only at energy output steps (set below).
2057          */
2058
2059         bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2060         bCalcEnerPres = bNstEner;
2061
2062         /* Do we need global communication ? */
2063         bGStat = (bGStatEveryStep || bStopCM || bNS ||
2064                   (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2065
2066         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2067
2068         if (do_ene || do_log)
2069         {
2070             bCalcEnerPres = TRUE;
2071             bGStat        = TRUE;
2072         }
2073
2074         /* Determine the wallclock run time up till now */
2075         run_time = gmx_gettime() - (double)runtime->real;
2076
2077         /* Check whether everything is still allright */
2078         if (((int)gmx_get_stop_condition() > handled_stop_condition)
2079 #ifdef GMX_THREADS
2080             && MASTER(cr)
2081 #endif
2082             )
2083         {
2084             /* this is just make gs.sig compatible with the hack
2085                of sending signals around by MPI_Reduce with together with
2086                other floats */
2087             if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2088                 gs.sig[eglsSTOPCOND]=1;
2089             if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2090                 gs.sig[eglsSTOPCOND]=-1;
2091             /* < 0 means stop at next step, > 0 means stop at next NS step */
2092             if (fplog)
2093             {
2094                 fprintf(fplog,
2095                         "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2096                         gmx_get_signal_name(),
2097                         gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2098                 fflush(fplog);
2099             }
2100             fprintf(stderr,
2101                     "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2102                     gmx_get_signal_name(),
2103                     gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2104             fflush(stderr);
2105             handled_stop_condition=(int)gmx_get_stop_condition();
2106         }
2107         else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
2108                  (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
2109                  gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
2110         {
2111             /* Signal to terminate the run */
2112             gs.sig[eglsSTOPCOND] = 1;
2113             if (fplog)
2114             {
2115                 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2116             }
2117             fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2118         }
2119
2120         if (bResetCountersHalfMaxH && MASTER(cr) &&
2121             run_time > max_hours*60.0*60.0*0.495)
2122         {
2123             gs.sig[eglsRESETCOUNTERS] = 1;
2124         }
2125
2126         if (ir->nstlist == -1 && !bRerunMD)
2127         {
2128             /* When bGStatEveryStep=FALSE, global_stat is only called
2129              * when we check the atom displacements, not at NS steps.
2130              * This means that also the bonded interaction count check is not
2131              * performed immediately after NS. Therefore a few MD steps could
2132              * be performed with missing interactions.
2133              * But wrong energies are never written to file,
2134              * since energies are only written after global_stat
2135              * has been called.
2136              */
2137             if (step >= nlh.step_nscheck)
2138             {
2139                 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
2140                                                      nlh.scale_tot,state->x);
2141             }
2142             else
2143             {
2144                 /* This is not necessarily true,
2145                  * but step_nscheck is determined quite conservatively.
2146                  */
2147                 nlh.nabnsb = 0;
2148             }
2149         }
2150
2151         /* In parallel we only have to check for checkpointing in steps
2152          * where we do global communication,
2153          *  otherwise the other nodes don't know.
2154          */
2155         if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
2156                            cpt_period >= 0 &&
2157                            (cpt_period == 0 ||
2158                             run_time >= nchkpt*cpt_period*60.0)) &&
2159             gs.set[eglsCHKPT] == 0)
2160         {
2161             gs.sig[eglsCHKPT] = 1;
2162         }
2163
2164         if (bIterations)
2165         {
2166             gmx_iterate_init(&iterate,bIterations);
2167         }
2168
2169         /* for iterations, we save these vectors, as we will be redoing the calculations */
2170         if (bIterations && iterate.bIterate)
2171         {
2172             copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2173         }
2174         bFirstIterate = TRUE;
2175         while (bFirstIterate || (bIterations && iterate.bIterate))
2176         {
2177             /* We now restore these vectors to redo the calculation with improved extended variables */
2178             if (bIterations)
2179             {
2180                 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2181             }
2182
2183             /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2184                so scroll down for that logic */
2185
2186             /* #########   START SECOND UPDATE STEP ################# */
2187             GMX_MPE_LOG(ev_update_start);
2188             bOK = TRUE;
2189             if (!bRerunMD || rerun_fr.bV || bForceUpdate)
2190             {
2191                 wallcycle_start(wcycle,ewcUPDATE);
2192                 dvdl = 0;
2193                 /* Box is changed in update() when we do pressure coupling,
2194                  * but we should still use the old box for energy corrections and when
2195                  * writing it to the energy file, so it matches the trajectory files for
2196                  * the same timestep above. Make a copy in a separate array.
2197                  */
2198                 copy_mat(state->box,lastbox);
2199                 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
2200                 if (bTrotter)
2201                 {
2202                     if (bIterations && iterate.bIterate)
2203                     {
2204                         if (bFirstIterate)
2205                         {
2206                             scalevir = 1;
2207                         }
2208                         else
2209                         {
2210                             /* we use a new value of scalevir to converge the iterations faster */
2211                             scalevir = tracevir/trace(shake_vir);
2212                         }
2213                         msmul(shake_vir,scalevir,shake_vir);
2214                         m_add(force_vir,shake_vir,total_vir);
2215                         clear_mat(shake_vir);
2216                     }
2217                     trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ3);
2218                 }
2219                 /* We can only do Berendsen coupling after we have summed
2220                  * the kinetic energy or virial. Since the happens
2221                  * in global_state after update, we should only do it at
2222                  * step % nstlist = 1 with bGStatEveryStep=FALSE.
2223                  */
2224
2225                 if (ir->eI != eiVVAK)
2226                 {
2227                   update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2228                 }
2229                 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
2230                                 upd,bInitStep);
2231
2232                 if (bVV)
2233                 {
2234                     /* velocity half-step update */
2235                     update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2236                                   ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,cr,nrnb,constr,&top->idef);
2237                 }
2238
2239                 /* Above, initialize just copies ekinh into ekin,
2240                  * it doesn't copy position (for VV),
2241                  * and entire integrator for MD.
2242                  */
2243
2244                 if (ir->eI==eiVVAK)
2245                 {
2246                     copy_rvecn(state->x,cbuf,0,state->natoms);
2247                 }
2248
2249                 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2250                               ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2251                 wallcycle_stop(wcycle,ewcUPDATE);
2252
2253                 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2254                                    &top->idef,shake_vir,force_vir,
2255                                    cr,nrnb,wcycle,upd,constr,
2256                                    bInitStep,FALSE,bCalcEnerPres,state->veta);
2257
2258                 if (ir->eI==eiVVAK)
2259                 {
2260                     /* erase F_EKIN and F_TEMP here? */
2261                     /* just compute the kinetic energy at the half step to perform a trotter step */
2262                     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2263                                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2264                                     constr,NULL,FALSE,lastbox,
2265                                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2266                                     cglo_flags | CGLO_TEMPERATURE | CGLO_CONSTRAINT
2267                         );
2268                     wallcycle_start(wcycle,ewcUPDATE);
2269                     trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ4);
2270                     /* now we know the scaling, we can compute the positions again again */
2271                     copy_rvecn(cbuf,state->x,0,state->natoms);
2272
2273                     update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2274                                   ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2275                     wallcycle_stop(wcycle,ewcUPDATE);
2276
2277                     /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2278                     /* are the small terms in the shake_vir here due
2279                      * to numerical errors, or are they important
2280                      * physically? I'm thinking they are just errors, but not completely sure.
2281                      * For now, will call without actually constraining, constr=NULL*/
2282                     update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2283                                        &top->idef,tmp_vir,force_vir,
2284                                        cr,nrnb,wcycle,upd,NULL,
2285                                        bInitStep,FALSE,bCalcEnerPres,state->veta);
2286                 }
2287                 if (!bOK && !bFFscan)
2288                 {
2289                     gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2290                 }
2291
2292                 if (fr->bSepDVDL && fplog && do_log)
2293                 {
2294                     fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2295                 }
2296                 enerd->term[F_DHDL_CON] += dvdl;
2297             }
2298             else if (graph)
2299             {
2300                 /* Need to unshift here */
2301                 unshift_self(graph,state->box,state->x);
2302             }
2303
2304             GMX_BARRIER(cr->mpi_comm_mygroup);
2305             GMX_MPE_LOG(ev_update_finish);
2306
2307             if (vsite != NULL)
2308             {
2309                 wallcycle_start(wcycle,ewcVSITECONSTR);
2310                 if (graph != NULL)
2311                 {
2312                     shift_self(graph,state->box,state->x);
2313                 }
2314                 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2315                                  top->idef.iparams,top->idef.il,
2316                                  fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2317
2318                 if (graph != NULL)
2319                 {
2320                     unshift_self(graph,state->box,state->x);
2321                 }
2322                 wallcycle_stop(wcycle,ewcVSITECONSTR);
2323             }
2324
2325             /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
2326             if (ir->nstlist == -1 && bFirstIterate)
2327             {
2328                 gs.sig[eglsNABNSB] = nlh.nabnsb;
2329             }
2330             compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2331                             wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2332                             constr,
2333                             bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
2334                             lastbox,
2335                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2336                             cglo_flags
2337                             | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
2338                             | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
2339                             | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
2340                             | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
2341                             | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2342                             | CGLO_CONSTRAINT
2343                 );
2344             if (ir->nstlist == -1 && bFirstIterate)
2345             {
2346                 nlh.nabnsb = gs.set[eglsNABNSB];
2347                 gs.set[eglsNABNSB] = 0;
2348             }
2349             /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
2350             /* #############  END CALC EKIN AND PRESSURE ################# */
2351
2352             /* Note: this is OK, but there are some numerical precision issues with using the convergence of
2353                the virial that should probably be addressed eventually. state->veta has better properies,
2354                but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
2355                generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
2356
2357             if (bIterations &&
2358                 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2359                                trace(shake_vir),&tracevir))
2360             {
2361                 break;
2362             }
2363             bFirstIterate = FALSE;
2364         }
2365
2366         update_box(fplog,step,ir,mdatoms,state,graph,f,
2367                    ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
2368
2369         /* ################# END UPDATE STEP 2 ################# */
2370         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
2371
2372         /* The coordinates (x) were unshifted in update */
2373 /*        if (bFFscan && (shellfc==NULL || bConverged))
2374         {
2375             if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
2376                                  f,NULL,xcopy,
2377                                  &(top_global->mols),mdatoms->massT,pres))
2378             {
2379                 if (gmx_parallel_env_initialized())
2380                 {
2381                     gmx_finalize();
2382                 }
2383                 fprintf(stderr,"\n");
2384                 exit(0);
2385             }
2386         }*/
2387         if (!bGStat)
2388         {
2389             /* We will not sum ekinh_old,
2390              * so signal that we still have to do it.
2391              */
2392             bSumEkinhOld = TRUE;
2393         }
2394
2395 /*        if (bTCR)
2396         {*/
2397             /* Only do GCT when the relaxation of shells (minimization) has converged,
2398              * otherwise we might be coupling to bogus energies.
2399              * In parallel we must always do this, because the other sims might
2400              * update the FF.
2401              */
2402
2403             /* Since this is called with the new coordinates state->x, I assume
2404              * we want the new box state->box too. / EL 20040121
2405              */
2406 /*            do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
2407                         ir,MASTER(cr),
2408                         mdatoms,&(top->idef),mu_aver,
2409                         top_global->mols.nr,cr,
2410                         state->box,total_vir,pres,
2411                         mu_tot,state->x,f,bConverged);
2412             debug_gmx();
2413         }*/
2414
2415         /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
2416
2417         sum_dhdl(enerd,state->lambda,ir);
2418         /* use the directly determined last velocity, not actually the averaged half steps */
2419         if (bTrotter && ir->eI==eiVV)
2420         {
2421             enerd->term[F_EKIN] = last_ekin;
2422         }
2423         enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
2424
2425         switch (ir->etc)
2426         {
2427         case etcNO:
2428             break;
2429         case etcBERENDSEN:
2430             break;
2431         case etcNOSEHOOVER:
2432             if (IR_NVT_TROTTER(ir)) {
2433                 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + last_conserved;
2434             } else {
2435                 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] +
2436                     NPT_energy(ir,state,&MassQ);
2437             }
2438             break;
2439         case etcVRESCALE:
2440             enerd->term[F_ECONSERVED] =
2441                 enerd->term[F_ETOT] + vrescale_energy(&(ir->opts),
2442                                                       state->therm_integral);
2443             break;
2444         default:
2445             break;
2446         }
2447
2448         /* Check for excessively large energies */
2449 /*        if (bIonize)
2450         {
2451 #ifdef GMX_DOUBLE
2452             real etot_max = 1e200;
2453 #else
2454             real etot_max = 1e30;
2455 #endif
2456             if (fabs(enerd->term[F_ETOT]) > etot_max)
2457             {
2458                 fprintf(stderr,"Energy too large (%g), giving up\n",
2459                         enerd->term[F_ETOT]);
2460             }
2461         }*/
2462         /* #########  END PREPARING EDR OUTPUT  ###########  */
2463
2464         /* Time for performance */
2465         if (((step % stepout) == 0) || bLastStep)
2466         {
2467             runtime_upd_proc(runtime);
2468         }
2469
2470         /* Output stuff */
2471         if (MASTER(cr))
2472         {
2473             gmx_bool do_dr,do_or;
2474
2475             if (!(bStartingFromCpt && (EI_VV(ir->eI))))
2476             {
2477                 if (bNstEner)
2478                 {
2479                     upd_mdebin(mdebin,bDoDHDL,TRUE,
2480                                t,mdatoms->tmass,enerd,state,lastbox,
2481                                shake_vir,force_vir,total_vir,pres,
2482                                ekind,mu_tot,constr);
2483                 }
2484                 else
2485                 {
2486                     upd_mdebin_step(mdebin);
2487                 }
2488
2489                 do_dr  = do_per_step(step,ir->nstdisreout);
2490                 do_or  = do_per_step(step,ir->nstorireout);
2491
2492                 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
2493                            step,t,
2494                            eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
2495             }
2496             if (ir->ePull != epullNO)
2497             {
2498                 pull_print_output(ir->pull,step,t);
2499             }
2500
2501             if (do_per_step(step,ir->nstlog))
2502             {
2503                 if(fflush(fplog) != 0)
2504                 {
2505                     gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
2506                 }
2507             }
2508         }
2509
2510
2511         /* Remaining runtime */
2512         if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
2513         {
2514             if (shellfc)
2515             {
2516                 fprintf(stderr,"\n");
2517             }
2518             print_time(stderr,runtime,step,ir,cr);
2519         }
2520
2521                 /* Set new positions for the group to embed */
2522                 if(!bLastStep){
2523                         if(step_rel<=it_xy)
2524                         {
2525                                 fac[0]+=xy_step;
2526                                 fac[1]+=xy_step;
2527                         } else if (step_rel<=(it_xy+it_z))
2528                         {
2529                                 fac[2]+=z_step;
2530                         }
2531                         resize(ins_at,r_ins,state_global->x,pos_ins,fac);
2532                 }
2533
2534         /* Replica exchange */
2535 /*        bExchanged = FALSE;
2536         if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
2537             do_per_step(step,repl_ex_nst))
2538         {
2539             bExchanged = replica_exchange(fplog,cr,repl_ex,
2540                                           state_global,enerd->term,
2541                                           state,step,t);
2542         }
2543         if (bExchanged && PAR(cr))
2544         {
2545             if (DOMAINDECOMP(cr))
2546             {
2547                 dd_partition_system(fplog,step,cr,TRUE,1,
2548                                     state_global,top_global,ir,
2549                                     state,&f,mdatoms,top,fr,
2550                                     vsite,shellfc,constr,
2551                                     nrnb,wcycle,FALSE);
2552             }
2553             else
2554             {
2555                 bcast_state(cr,state,FALSE);
2556             }
2557         }*/
2558
2559         bFirstStep = FALSE;
2560         bInitStep = FALSE;
2561         bStartingFromCpt = FALSE;
2562
2563         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2564         /* With all integrators, except VV, we need to retain the pressure
2565          * at the current step for coupling at the next step.
2566          */
2567         if ((state->flags & (1<<estPRES_PREV)) &&
2568             (bGStatEveryStep ||
2569              (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
2570         {
2571             /* Store the pressure in t_state for pressure coupling
2572              * at the next MD step.
2573              */
2574             copy_mat(pres,state->pres_prev);
2575         }
2576
2577         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
2578
2579         if (bRerunMD)
2580         {
2581             /* read next frame from input trajectory */
2582             bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
2583         }
2584
2585         if (!bRerunMD || !rerun_fr.bStep)
2586         {
2587             /* increase the MD step number */
2588             step++;
2589             step_rel++;
2590         }
2591
2592         cycles = wallcycle_stop(wcycle,ewcSTEP);
2593         if (DOMAINDECOMP(cr) && wcycle)
2594         {
2595             dd_cycles_add(cr->dd,cycles,ddCyclStep);
2596         }
2597
2598         if (step_rel == wcycle_get_reset_counters(wcycle) ||
2599             gs.set[eglsRESETCOUNTERS] != 0)
2600         {
2601             /* Reset all the counters related to performance over the run */
2602             reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
2603             wcycle_set_reset_counters(wcycle,-1);
2604             bResetCountersHalfMaxH = FALSE;
2605             gs.set[eglsRESETCOUNTERS] = 0;
2606         }
2607     }
2608     /* End of main MD loop */
2609     debug_gmx();
2610     write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2611                                         *top_global->name,top_global,
2612                                         state_global->x,state_global->v,
2613                                         ir->ePBC,state->box);
2614
2615     /* Stop the time */
2616     runtime_end(runtime);
2617
2618     if (bRerunMD)
2619     {
2620         close_trj(status);
2621     }
2622
2623     if (!(cr->duty & DUTY_PME))
2624     {
2625         /* Tell the PME only node to finish */
2626         gmx_pme_finish(cr);
2627     }
2628
2629     if (MASTER(cr))
2630     {
2631         if (ir->nstcalcenergy > 0 && !bRerunMD)
2632         {
2633             print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
2634                        eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
2635         }
2636     }
2637
2638     done_mdoutf(outf);
2639
2640     debug_gmx();
2641
2642     if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2643     {
2644         fprintf(fplog,"Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n",nlh.s1/nlh.nns,sqrt(nlh.s2/nlh.nns - sqr(nlh.s1/nlh.nns)));
2645         fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
2646     }
2647
2648     if (shellfc && fplog)
2649     {
2650         fprintf(fplog,"Fraction of iterations that converged:           %.2f %%\n",
2651                 (nconverged*100.0)/step_rel);
2652         fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
2653                 tcount/step_rel);
2654     }
2655
2656 /*    if (repl_ex_nst > 0 && MASTER(cr))
2657     {
2658         print_replica_exchange_statistics(fplog,repl_ex);
2659     }*/
2660
2661     runtime->nsteps_done = step_rel;
2662
2663     return 0;
2664 }
2665
2666
2667 int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
2668              const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
2669              int nstglobalcomm,
2670              ivec ddxyz,int dd_node_order,real rdd,real rconstr,
2671              const char *dddlb_opt,real dlb_scale,
2672              const char *ddcsx,const char *ddcsy,const char *ddcsz,
2673              int nstepout,int resetstep,int nmultisim,int repl_ex_nst,int repl_ex_seed,
2674              real pforce,real cpt_period,real max_hours,
2675              const char *deviceOptions,
2676              unsigned long Flags,
2677              real xy_fac, real xy_max, real z_fac, real z_max,
2678              int it_xy, int it_z, real probe_rad, int low_up_rm,
2679              int pieces, gmx_bool bALLOW_ASYMMETRY, int maxwarn)
2680 {
2681     double     nodetime=0,realtime;
2682     t_inputrec *inputrec;
2683     t_state    *state=NULL;
2684     matrix     box;
2685     gmx_ddbox_t ddbox;
2686     int        npme_major,npme_minor;
2687     real       tmpr1,tmpr2;
2688     t_nrnb     *nrnb;
2689     gmx_mtop_t *mtop=NULL;
2690     t_mdatoms  *mdatoms=NULL;
2691     t_forcerec *fr=NULL;
2692     t_fcdata   *fcd=NULL;
2693     real       ewaldcoeff=0;
2694     gmx_pme_t  *pmedata=NULL;
2695     gmx_vsite_t *vsite=NULL;
2696     gmx_constr_t constr;
2697     int        i,m,nChargePerturbed=-1,status,nalloc;
2698     char       *gro;
2699     gmx_wallcycle_t wcycle;
2700     gmx_bool       bReadRNG,bReadEkin;
2701     int        list;
2702     gmx_runtime_t runtime;
2703     int        rc;
2704     gmx_large_int_t reset_counters;
2705     gmx_edsam_t ed=NULL;
2706     t_commrec   *cr_old=cr;
2707     int        nthreads=1,nthreads_requested=1;
2708
2709
2710         char                    *ins;
2711         int                     rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
2712         int                     ng,j,max_lip_rm,ins_grp_id,ins_nat,mem_nat,ntype,lip_rm,tpr_version;
2713         real                    xy_step=0,z_step=0;
2714         real                    prot_area;
2715         rvec                    *r_ins=NULL,fac;
2716         t_block                 *ins_at,*rest_at;
2717         pos_ins_t               *pos_ins;
2718         mem_t                   *mem_p;
2719         rmm_t                   *rm_p;
2720         gmx_groups_t            *groups;
2721         gmx_bool                        bExcl=FALSE;
2722         t_atoms                 atoms;
2723         t_pbc                   *pbc;
2724         char                    **piecename=NULL;
2725
2726     /* CAUTION: threads may be started later on in this function, so
2727        cr doesn't reflect the final parallel state right now */
2728     snew(inputrec,1);
2729     snew(mtop,1);
2730
2731     if (bVerbose && SIMMASTER(cr))
2732     {
2733         fprintf(stderr,"Getting Loaded...\n");
2734     }
2735
2736     if (Flags & MD_APPENDFILES)
2737     {
2738         fplog = NULL;
2739     }
2740
2741     snew(state,1);
2742     if (MASTER(cr))
2743     {
2744         /* Read (nearly) all data required for the simulation */
2745         read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
2746
2747         /* NOW the threads will be started: */
2748 #ifdef GMX_THREADS
2749 #endif
2750     }
2751     /* END OF CAUTION: cr is now reliable */
2752
2753     if (PAR(cr))
2754     {
2755         /* now broadcast everything to the non-master nodes/threads: */
2756         init_parallel(fplog, cr, inputrec, mtop);
2757     }
2758     /* now make sure the state is initialized and propagated */
2759     set_state_entries(state,inputrec,cr->nnodes);
2760
2761     if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
2762     {
2763         /* All-vs-all loops do not work with domain decomposition */
2764         Flags |= MD_PARTDEC;
2765     }
2766
2767     if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
2768     {
2769         cr->npmenodes = 0;
2770     }
2771
2772         snew(ins_at,1);
2773         snew(pos_ins,1);
2774         if(MASTER(cr))
2775         {
2776                 tpr_version = get_tpr_version(ftp2fn(efTPX,nfile,fnm));
2777                 if (tpr_version<58)
2778                         gmx_fatal(FARGS,"Version of *.tpr file to old (%d). Rerun grompp with gromacs VERSION 4.0.3 or newer.\n",tpr_version);
2779
2780                 if( inputrec->eI != eiMD )
2781                         gmx_input("Change integrator to md in mdp file.");
2782
2783                 if(PAR(cr))
2784                         gmx_input("Sorry, parallel g_membed is not yet fully functrional.");
2785
2786                 groups=&(mtop->groups);
2787
2788                 atoms=gmx_mtop_global_atoms(mtop);
2789                 snew(mem_p,1);
2790                 fprintf(stderr,"\nSelect a group to embed in the membrane:\n");
2791                 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(ins_at->nr),&(ins_at->index),&ins);
2792                 ins_grp_id = search_string(ins,groups->ngrpname,(groups->grpname));
2793                 fprintf(stderr,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins);
2794                 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(mem_p->mem_at.nr),&(mem_p->mem_at.index),&(mem_p->name));
2795
2796                 pos_ins->pieces=pieces;
2797                 snew(pos_ins->nidx,pieces);
2798                 snew(pos_ins->subindex,pieces);
2799                 snew(piecename,pieces);
2800                 if (pieces>1)
2801                 {
2802                         fprintf(stderr,"\nSelect pieces to embed:\n");
2803                         get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),pieces,pos_ins->nidx,pos_ins->subindex,piecename);
2804                 }
2805                 else
2806                 {
2807                         /*use whole embedded group*/
2808                         snew(pos_ins->nidx,1);
2809                         snew(pos_ins->subindex,1);
2810                         pos_ins->nidx[0]=ins_at->nr;
2811                         pos_ins->subindex[0]=ins_at->index;
2812                 }
2813
2814                 if(probe_rad<0.2199999)
2815                 {
2816                         warn++;
2817                         fprintf(stderr,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 can result in overlap between waters "
2818                                         "and the group to embed, which will result in Lincs errors etc.\nIf you are sure, you can increase maxwarn.\n\n",warn);
2819                 }
2820
2821                 if(xy_fac<0.09999999)
2822                 {
2823                         warn++;
2824                         fprintf(stderr,"\nWarning %d:\nThe initial size of %s is probably too smal.\n"
2825                                         "If you are sure, you can increase maxwarn.\n\n",warn,ins);
2826                 }
2827
2828                 if(it_xy<1000)
2829                 {
2830                         warn++;
2831                         fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d) is probably too small.\n"
2832                                         "Increase -nxy or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_xy);
2833                 }
2834
2835                 if( (it_z<100) && ( z_fac<0.99999999 || z_fac>1.0000001) )
2836                 {
2837                         warn++;
2838                         fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d) is probably too small.\n"
2839                                        "Increase -nz or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_z);
2840                 }
2841
2842                 if(it_xy+it_z>inputrec->nsteps)
2843                 {
2844                         warn++;
2845                         fprintf(stderr,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the number of steps in the tpr.\n"
2846                                         "If you are sure, you can increase maxwarn.\n\n",warn);
2847                 }
2848
2849                 fr_id=-1;
2850                 if( inputrec->opts.ngfrz==1)
2851                         gmx_fatal(FARGS,"You did not specify \"%s\" as a freezegroup.",ins);
2852                 for(i=0;i<inputrec->opts.ngfrz;i++)
2853                 {
2854                         tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
2855                         if(ins_grp_id==tmp_id)
2856                         {
2857                                 fr_id=tmp_id;
2858                                 fr_i=i;
2859                         }
2860                 }
2861                 if (fr_id == -1 )
2862                         gmx_fatal(FARGS,"\"%s\" not as freezegroup defined in the mdp-file.",ins);
2863
2864                 for(i=0;i<DIM;i++)
2865                         if( inputrec->opts.nFreeze[fr_i][i] != 1)
2866                                 gmx_fatal(FARGS,"freeze dimensions for %s are not Y Y Y\n",ins);
2867
2868                 ng = groups->grps[egcENER].nr;
2869                 if (ng == 1)
2870                         gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
2871
2872                 for(i=0;i<ng;i++)
2873                 {
2874                         for(j=0;j<ng;j++)
2875                         {
2876                                 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
2877                                 {
2878                                         bExcl = TRUE;
2879                                         if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
2880                                                 gmx_fatal(FARGS,"Energy exclusions \"%s\" and  \"%s\" do not match the group to embed \"%s\"",
2881                                                                 *groups->grpname[groups->grps[egcENER].nm_ind[i]],
2882                                                                 *groups->grpname[groups->grps[egcENER].nm_ind[j]],ins);
2883                                 }
2884                         }
2885                 }
2886                 if (!bExcl)
2887                         gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in the freeze group");
2888
2889                 /* Set all atoms in box*/
2890                 /*set_inbox(state->natoms,state->x);*/
2891
2892                 /* Guess the area the protein will occupy in the membrane plane  Calculate area per lipid*/
2893                 snew(rest_at,1);
2894                 ins_nat = init_ins_at(ins_at,rest_at,state,pos_ins,groups,ins_grp_id,xy_max);
2895                 /* Check moleculetypes in insertion group */
2896                 check_types(ins_at,rest_at,mtop);
2897
2898                 mem_nat = init_mem_at(mem_p,mtop,state->x,state->box,pos_ins);
2899
2900                 prot_area = est_prot_area(pos_ins,state->x,ins_at,mem_p);
2901                 if ( (prot_area>7.5) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX])<50) )
2902                 {
2903                         warn++;
2904                         fprintf(stderr,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
2905                                         "This might cause pressure problems during the growth phase. Just try with\n"
2906                                         "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
2907                                         "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn);
2908                 }
2909                 if(warn>maxwarn)
2910                                         gmx_fatal(FARGS,"Too many warnings.\n");
2911
2912                 printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area);
2913                 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\nThe area per lipid is %.4f nm^2.\n",mem_p->nmol,mem_p->lip_area);
2914
2915                 /* Maximum number of lipids to be removed*/
2916                 max_lip_rm=(int)(2*prot_area/mem_p->lip_area);
2917                 printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm);
2918
2919                 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
2920                                 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
2921                                 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",xy_fac,z_fac,mem_p->zmin,mem_p->zmax);
2922
2923                 /* resize the protein by xy and by z if necessary*/
2924                 snew(r_ins,ins_at->nr);
2925                 init_resize(ins_at,r_ins,pos_ins,mem_p,state->x,bALLOW_ASYMMETRY);
2926                 fac[0]=fac[1]=xy_fac;
2927                 fac[2]=z_fac;
2928
2929                 xy_step =(xy_max-xy_fac)/(double)(it_xy);
2930                 z_step  =(z_max-z_fac)/(double)(it_z-1);
2931
2932                 resize(ins_at,r_ins,state->x,pos_ins,fac);
2933
2934                 /* remove overlapping lipids and water from the membrane box*/
2935                 /*mark molecules to be removed*/
2936                 snew(pbc,1);
2937                 set_pbc(pbc,inputrec->ePBC,state->box);
2938
2939                 snew(rm_p,1);
2940                 lip_rm = gen_rm_list(rm_p,ins_at,rest_at,pbc,mtop,state->x, r_ins, mem_p,pos_ins,probe_rad,low_up_rm,bALLOW_ASYMMETRY);
2941         lip_rm -= low_up_rm;
2942
2943                 if(fplog)
2944                         for(i=0;i<rm_p->nr;i++)
2945                                 fprintf(fplog,"rm mol %d\n",rm_p->mol[i]);
2946
2947                 for(i=0;i<mtop->nmolblock;i++)
2948                 {
2949                         ntype=0;
2950                         for(j=0;j<rm_p->nr;j++)
2951                                 if(rm_p->block[j]==i)
2952                                         ntype++;
2953                         printf("Will remove %d %s molecules\n",ntype,*(mtop->moltype[mtop->molblock[i].type].name));
2954                 }
2955
2956                 if(lip_rm>max_lip_rm)
2957                 {
2958                         warn++;
2959                         fprintf(stderr,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated protein area\n"
2960                                         "Try making the -xyinit resize factor smaller. If you are sure about this increase maxwarn.\n\n",warn);
2961                 }
2962
2963                 /*remove all lipids and waters overlapping and update all important structures*/
2964                 rm_group(inputrec,groups,mtop,rm_p,state,ins_at,pos_ins);
2965
2966                 rm_bonded_at = rm_bonded(ins_at,mtop);
2967                 if (rm_bonded_at != ins_at->nr)
2968                 {
2969                         fprintf(stderr,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
2970                                         "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
2971                                         "molecule type as atoms that are not to be embedded.\n",rm_bonded_at,ins_at->nr);
2972                 }
2973
2974                 if(warn>maxwarn)
2975                         gmx_fatal(FARGS,"Too many warnings.\nIf you are sure these warnings are harmless, you can increase -maxwarn");
2976
2977                 if (MASTER(cr))
2978                 {
2979                         if (ftp2bSet(efTOP,nfile,fnm))
2980                                 top_update(opt2fn("-p",nfile,fnm),ins,rm_p,mtop);
2981                 }
2982
2983                 sfree(pbc);
2984                 sfree(rest_at);
2985         }
2986
2987 #ifdef GMX_FAHCORE
2988     fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
2989 #endif
2990
2991     /* NMR restraints must be initialized before load_checkpoint,
2992      * since with time averaging the history is added to t_state.
2993      * For proper consistency check we therefore need to extend
2994      * t_state here.
2995      * So the PME-only nodes (if present) will also initialize
2996      * the distance restraints.
2997      */
2998     snew(fcd,1);
2999
3000     /* This needs to be called before read_checkpoint to extend the state */
3001     init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
3002
3003     if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
3004     {
3005         if (PAR(cr) && !(Flags & MD_PARTDEC))
3006         {
3007             gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
3008         }
3009         /* Orientation restraints */
3010         if (MASTER(cr))
3011         {
3012             init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
3013                         state);
3014         }
3015     }
3016
3017     if (DEFORM(*inputrec))
3018     {
3019         /* Store the deform reference box before reading the checkpoint */
3020         if (SIMMASTER(cr))
3021         {
3022             copy_mat(state->box,box);
3023         }
3024         if (PAR(cr))
3025         {
3026             gmx_bcast(sizeof(box),box,cr);
3027         }
3028         /* Because we do not have the update struct available yet
3029          * in which the reference values should be stored,
3030          * we store them temporarily in static variables.
3031          * This should be thread safe, since they are only written once
3032          * and with identical values.
3033          */
3034 /*        deform_init_init_step_tpx = inputrec->init_step;*/
3035 /*        copy_mat(box,deform_init_box_tpx);*/
3036     }
3037
3038     if (opt2bSet("-cpi",nfile,fnm))
3039     {
3040         /* Check if checkpoint file exists before doing continuation.
3041          * This way we can use identical input options for the first and subsequent runs...
3042          */
3043         if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
3044         {
3045             load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
3046                             cr,Flags & MD_PARTDEC,ddxyz,
3047                             inputrec,state,&bReadRNG,&bReadEkin,
3048                             (Flags & MD_APPENDFILES));
3049
3050             if (bReadRNG)
3051             {
3052                 Flags |= MD_READ_RNG;
3053             }
3054             if (bReadEkin)
3055             {
3056                 Flags |= MD_READ_EKIN;
3057             }
3058         }
3059     }
3060
3061     if ((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
3062     {
3063         gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
3064                              Flags,&fplog);
3065     }
3066
3067     if (SIMMASTER(cr))
3068     {
3069         copy_mat(state->box,box);
3070     }
3071
3072     if (PAR(cr))
3073     {
3074         gmx_bcast(sizeof(box),box,cr);
3075     }
3076
3077     if (bVerbose && SIMMASTER(cr))
3078     {
3079         fprintf(stderr,"Loaded with Money\n\n");
3080     }
3081
3082     if (PAR(cr) && !((Flags & MD_PARTDEC) || EI_TPI(inputrec->eI)))
3083     {
3084         cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
3085                                            dddlb_opt,dlb_scale,
3086                                            ddcsx,ddcsy,ddcsz,
3087                                            mtop,inputrec,
3088                                            box,state->x,
3089                                            &ddbox,&npme_major,&npme_minor);
3090
3091         make_dd_communicators(fplog,cr,dd_node_order);
3092
3093         /* Set overallocation to avoid frequent reallocation of arrays */
3094         set_over_alloc_dd(TRUE);
3095     }
3096     else
3097     {
3098         /* PME, if used, is done on all nodes with 1D decomposition */
3099         cr->npmenodes = 0;
3100         cr->duty = (DUTY_PP | DUTY_PME);
3101         npme_major = cr->nnodes;
3102         npme_minor = 1;
3103
3104         if (inputrec->ePBC == epbcSCREW)
3105         {
3106             gmx_fatal(FARGS,
3107                       "pbc=%s is only implemented with domain decomposition",
3108                       epbc_names[inputrec->ePBC]);
3109         }
3110     }
3111
3112     if (PAR(cr))
3113     {
3114         /* After possible communicator splitting in make_dd_communicators.
3115          * we can set up the intra/inter node communication.
3116          */
3117         gmx_setup_nodecomm(fplog,cr);
3118     }
3119
3120     wcycle = wallcycle_init(fplog,resetstep,cr);
3121     if (PAR(cr))
3122     {
3123         /* Master synchronizes its value of reset_counters with all nodes
3124          * including PME only nodes */
3125         reset_counters = wcycle_get_reset_counters(wcycle);
3126         gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
3127         wcycle_set_reset_counters(wcycle, reset_counters);
3128     }
3129
3130
3131     snew(nrnb,1);
3132     if (cr->duty & DUTY_PP)
3133     {
3134         /* For domain decomposition we allocate dynamically
3135          * in dd_partition_system.
3136          */
3137         if (DOMAINDECOMP(cr))
3138         {
3139             bcast_state_setup(cr,state);
3140         }
3141         else
3142         {
3143             if (PAR(cr))
3144             {
3145                 if (!MASTER(cr))
3146                 {
3147                     snew(state,1);
3148                 }
3149                 bcast_state(cr,state,TRUE);
3150             }
3151         }
3152
3153         /* Dihedral Restraints */
3154         if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
3155         {
3156             init_dihres(fplog,mtop,inputrec,fcd);
3157         }
3158
3159         /* Initiate forcerecord */
3160         fr = mk_forcerec();
3161         init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
3162                       opt2fn("-table",nfile,fnm),
3163                       opt2fn("-tablep",nfile,fnm),
3164                       opt2fn("-tableb",nfile,fnm),FALSE,pforce);
3165
3166         /* version for PCA_NOT_READ_NODE (see md.c) */
3167         /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
3168           "nofile","nofile","nofile",FALSE,pforce);
3169           */
3170         fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
3171
3172         /* Initialize QM-MM */
3173         if(fr->bQMMM)
3174         {
3175             init_QMMMrec(cr,box,mtop,inputrec,fr);
3176         }
3177
3178         /* Initialize the mdatoms structure.
3179          * mdatoms is not filled with atom data,
3180          * as this can not be done now with domain decomposition.
3181          */
3182         mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
3183
3184         /* Initialize the virtual site communication */
3185         vsite = init_vsite(mtop,cr);
3186
3187         calc_shifts(box,fr->shift_vec);
3188
3189         /* With periodic molecules the charge groups should be whole at start up
3190          * and the virtual sites should not be far from their proper positions.
3191          */
3192         if (!inputrec->bContinuation && MASTER(cr) &&
3193             !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
3194         {
3195             /* Make molecules whole at start of run */
3196             if (fr->ePBC != epbcNONE)
3197             {
3198                 do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
3199             }
3200             if (vsite)
3201             {
3202                 /* Correct initial vsite positions are required
3203                  * for the initial distribution in the domain decomposition
3204                  * and for the initial shell prediction.
3205                  */
3206                 construct_vsites_mtop(fplog,vsite,mtop,state->x);
3207             }
3208         }
3209
3210         /* Initiate PPPM if necessary */
3211         if (fr->eeltype == eelPPPM)
3212         {
3213             if (mdatoms->nChargePerturbed)
3214             {
3215                 gmx_fatal(FARGS,"Free energy with %s is not implemented",
3216                           eel_names[fr->eeltype]);
3217             }
3218             status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
3219                                    getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
3220             if (status != 0)
3221             {
3222                 gmx_fatal(FARGS,"Error %d initializing PPPM",status);
3223             }
3224         }
3225
3226         if (EEL_PME(fr->eeltype))
3227         {
3228             ewaldcoeff = fr->ewaldcoeff;
3229             pmedata = &fr->pmedata;
3230         }
3231         else
3232         {
3233             pmedata = NULL;
3234         }
3235     }
3236     else
3237     {
3238         /* This is a PME only node */
3239
3240         /* We don't need the state */
3241         done_state(state);
3242
3243         ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
3244         snew(pmedata,1);
3245     }
3246
3247     /* Initiate PME if necessary,
3248      * either on all nodes or on dedicated PME nodes only. */
3249     if (EEL_PME(inputrec->coulombtype))
3250     {
3251         if (mdatoms)
3252         {
3253             nChargePerturbed = mdatoms->nChargePerturbed;
3254         }
3255         if (cr->npmenodes > 0)
3256         {
3257             /* The PME only nodes need to know nChargePerturbed */
3258             gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
3259         }
3260         if (cr->duty & DUTY_PME)
3261         {
3262             status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
3263                                   mtop ? mtop->natoms : 0,nChargePerturbed,
3264                                   (Flags & MD_REPRODUCIBLE));
3265             if (status != 0)
3266             {
3267                 gmx_fatal(FARGS,"Error %d initializing PME",status);
3268             }
3269         }
3270     }
3271
3272
3273 /*    if (integrator[inputrec->eI].func == do_md
3274 #ifdef GMX_OPENMM
3275         ||
3276         integrator[inputrec->eI].func == do_md_openmm
3277 #endif
3278         )
3279     {*/
3280         /* Turn on signal handling on all nodes */
3281         /*
3282          * (A user signal from the PME nodes (if any)
3283          * is communicated to the PP nodes.
3284          */
3285         signal_handler_install();
3286 /*    }*/
3287
3288     if (cr->duty & DUTY_PP)
3289     {
3290         if (inputrec->ePull != epullNO)
3291         {
3292             /* Initialize pull code */
3293             init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
3294                       EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
3295         }
3296
3297         constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
3298
3299         if (DOMAINDECOMP(cr))
3300         {
3301             dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
3302                             Flags & MD_DDBONDCHECK,fr->cginfo_mb);
3303
3304             set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
3305
3306             setup_dd_grid(fplog,cr->dd);
3307         }
3308
3309         /* Now do whatever the user wants us to do (how flexible...) */
3310         do_md_membed(fplog,cr,nfile,fnm,
3311                                       oenv,bVerbose,bCompact,
3312                                       nstglobalcomm,
3313                                       vsite,constr,
3314                                       nstepout,inputrec,mtop,
3315                                       fcd,state,
3316                                       mdatoms,nrnb,wcycle,ed,fr,
3317                                       repl_ex_nst,repl_ex_seed,
3318                                       cpt_period,max_hours,
3319                                       deviceOptions,
3320                                       Flags,
3321                                       &runtime,
3322                                       fac, r_ins, pos_ins, ins_at,
3323                                       xy_step, z_step, it_xy, it_z);
3324
3325         if (inputrec->ePull != epullNO)
3326         {
3327             finish_pull(fplog,inputrec->pull);
3328         }
3329     }
3330     else
3331     {
3332         /* do PME only */
3333         gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
3334     }
3335
3336     if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
3337     {
3338         /* Some timing stats */
3339         if (MASTER(cr))
3340         {
3341             if (runtime.proc == 0)
3342             {
3343                 runtime.proc = runtime.real;
3344             }
3345         }
3346         else
3347         {
3348             runtime.real = 0;
3349         }
3350     }
3351
3352     wallcycle_stop(wcycle,ewcRUN);
3353
3354     /* Finish up, write some stuff
3355      * if rerunMD, don't write last frame again
3356      */
3357     finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
3358                inputrec,nrnb,wcycle,&runtime,
3359                EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
3360
3361     /* Does what it says */
3362     print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
3363
3364     /* Close logfile already here if we were appending to it */
3365     if (MASTER(cr) && (Flags & MD_APPENDFILES))
3366     {
3367         gmx_log_close(fplog);
3368     }
3369
3370     if (pieces>1)
3371     {
3372         sfree(piecename);
3373     }
3374
3375     rc=(int)gmx_get_stop_condition();
3376
3377     return rc;
3378 }
3379
3380 int gmx_membed(int argc,char *argv[])
3381 {
3382         const char *desc[] = {
3383                         "[TT]g_membed[tt] embeds a membrane protein into an equilibrated lipid bilayer at the position",
3384                         "and orientation specified by the user.[PAR]",
3385                         "SHORT MANUAL[BR]------------[BR]",
3386                         "The user should merge the structure files of the protein and membrane (+solvent), creating a",
3387                         "single structure file with the protein overlapping the membrane at the desired position and",
3388                         "orientation. The box size is taken from the membrane structure file. The corresponding topology",
3389                         "files should also be merged. Consecutively, create a [TT].tpr[tt] file (input for [TT]g_membed[tt]) from these files,"
3390                         "with the following options included in the [TT].mdp[tt] file.[BR]",
3391                         " - [TT]integrator      = md[tt][BR]",
3392                         " - [TT]energygrp       = Protein[tt] (or other group that you want to insert)[BR]",
3393                         " - [TT]freezegrps      = Protein[tt][BR]",
3394                         " - [TT]freezedim       = Y Y Y[tt][BR]",
3395                         " - [TT]energygrp_excl  = Protein Protein[tt][BR]",
3396                         "The output is a structure file containing the protein embedded in the membrane. If a topology",
3397                         "file is provided, the number of lipid and ",
3398                         "solvent molecules will be updated to match the new structure file.[BR]",
3399                         "For a more extensive manual see Wolf et al, J Comp Chem 31 (2010) 2169-2174, Appendix.[PAR]",
3400                         "SHORT METHOD DESCRIPTION[BR]",
3401                         "------------------------[BR]",
3402                         "1. The protein is resized around its center of mass by a factor [TT]-xy[tt] in the xy-plane",
3403                         "(the membrane plane) and a factor [TT]-z[tt] in the [IT]z[it]-direction (if the size of the",
3404                         "protein in the z-direction is the same or smaller than the width of the membrane, a",
3405                         "[TT]-z[tt] value larger than 1 can prevent that the protein will be enveloped by the lipids).[BR]",
3406                         "2. All lipid and solvent molecules overlapping with the resized protein are removed. All",
3407                         "intraprotein interactions are turned off to prevent numerical issues for small values of [TT]-xy[tt]",
3408                         " or [TT]-z[tt][BR]",
3409                         "3. One md step is performed.[BR]",
3410                         "4. The resize factor ([TT]-xy[tt] or [TT]-z[tt]) is incremented by a small amount ((1-xy)/nxy or (1-z)/nz) and the",
3411                         "protein is resized again around its center of mass. The resize factor for the xy-plane",
3412                         "is incremented first. The resize factor for the z-direction is not changed until the [TT]-xy[tt] factor",
3413                         "is 1 (thus after [TT]-nxy[tt] iterations).[BR]",
3414                         "5. Repeat step 3 and 4 until the protein reaches its original size ([TT]-nxy[tt] + [TT]-nz[tt] iterations).[BR]",
3415                         "For a more extensive method description see Wolf et al, J Comp Chem, 31 (2010) 2169-2174.[PAR]",
3416                         "NOTE[BR]----[BR]",
3417                         " - Protein can be any molecule you want to insert in the membrane.[BR]",
3418                         " - It is recommended to perform a short equilibration run after the embedding",
3419                         "(see Wolf et al, J Comp Chem 31 (2010) 2169-2174, to re-equilibrate the membrane. Clearly",
3420                         "protein equilibration might require longer.\n",
3421                         " - It is now also possible to use the g_membed functionality with mdrun. You should than pass",
3422                         "a data file containing the command line options of g_membed following the -membed option, for",
3423                         "example mdrun -s into_mem.tpr -membed membed.dat.",
3424                         "\n"
3425         };
3426         t_filenm fnm[] = {
3427                         { efTPX, "-f",      "into_mem", ffREAD },
3428                         { efNDX, "-n",      "index",    ffOPTRD },
3429                         { efTOP, "-p",      "topol",    ffOPTRW },
3430                         { efTRN, "-o",      NULL,       ffWRITE },
3431                         { efXTC, "-x",      NULL,       ffOPTWR },
3432                         { efSTO, "-c",      "membedded",  ffWRITE },
3433                         { efEDR, "-e",      "ener",     ffWRITE },
3434                         { efDAT, "-dat",    "membed",   ffWRITE }
3435                         { efLOG, "-g",      "md",       ffWRITE },
3436                         { efEDI, "-ei",     "sam",      ffOPTRD },
3437                         { efTRX, "-rerun",  "rerun",    ffOPTRD },
3438                         { efXVG, "-table",  "table",    ffOPTRD },
3439                         { efXVG, "-tablep", "tablep",   ffOPTRD },
3440                         { efXVG, "-tableb", "table",    ffOPTRD },
3441                         { efXVG, "-dhdl",   "dhdl",     ffOPTWR },
3442                         { efXVG, "-field",  "field",    ffOPTWR },
3443                         { efXVG, "-table",  "table",    ffOPTRD },
3444                         { efXVG, "-tablep", "tablep",   ffOPTRD },
3445                         { efXVG, "-tableb", "table",    ffOPTRD },
3446                         { efTRX, "-rerun",  "rerun",    ffOPTRD },
3447                         { efXVG, "-tpi",    "tpi",      ffOPTWR },
3448                         { efXVG, "-tpid",   "tpidist",  ffOPTWR },
3449                         { efEDI, "-ei",     "sam",      ffOPTRD },
3450                         { efEDO, "-eo",     "sam",      ffOPTWR },
3451                         { efGCT, "-j",      "wham",     ffOPTRD },
3452                         { efGCT, "-jo",     "bam",      ffOPTWR },
3453                         { efXVG, "-ffout",  "gct",      ffOPTWR },
3454                         { efXVG, "-devout", "deviatie", ffOPTWR },
3455                         { efXVG, "-runav",  "runaver",  ffOPTWR },
3456                         { efXVG, "-px",     "pullx",    ffOPTWR },
3457                         { efXVG, "-pf",     "pullf",    ffOPTWR },
3458                         { efMTX, "-mtx",    "nm",       ffOPTWR },
3459                         { efNDX, "-dn",     "dipole",   ffOPTWR },
3460                         { efRND, "-multidir",NULL,      ffOPTRDMULT}
3461         };
3462 #define NFILE asize(fnm)
3463
3464         /* Command line options ! */
3465         real xy_fac = 0.5;
3466         real xy_max = 1.0;
3467         real z_fac = 1.0;
3468         real z_max = 1.0;
3469         int it_xy = 1000;
3470         int it_z = 0;
3471         real probe_rad = 0.22;
3472         int low_up_rm = 0;
3473         int maxwarn=0;
3474         int pieces=1;
3475         gmx_bool bALLOW_ASYMMETRY=FALSE;
3476         gmx_bool bStart=FALSE;
3477         int nstepout=100;
3478         gmx_bool bVerbose=FALSE;
3479         char *mdrun_path=NULL;
3480
3481 /* arguments relevant to OPENMM only*/
3482 #ifdef GMX_OPENMM
3483     gmx_input("g_membed not functional in openmm");
3484 #endif
3485
3486         t_pargs pa[] = {
3487                         { "-xyinit",   FALSE, etREAL,  {&xy_fac},
3488                                 "Resize factor for the protein in the xy dimension before starting embedding" },
3489                         { "-xyend",   FALSE, etREAL,  {&xy_max},
3490                                 "Final resize factor in the xy dimension" },
3491                         { "-zinit",    FALSE, etREAL,  {&z_fac},
3492                                 "Resize factor for the protein in the z dimension before starting embedding" },
3493                         { "-zend",    FALSE, etREAL,  {&z_max},
3494                                 "Final resize faction in the z dimension" },
3495                         { "-nxy",     FALSE,  etINT,  {&it_xy},
3496                                 "Number of iteration for the xy dimension" },
3497                         { "-nz",      FALSE,  etINT,  {&it_z},
3498                                 "Number of iterations for the z dimension" },
3499                         { "-rad",     FALSE, etREAL,  {&probe_rad},
3500                                 "Probe radius to check for overlap between the group to embed and the membrane"},
3501                         { "-pieces",  FALSE,  etINT,  {&pieces},
3502                                 "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
3503                         { "-asymmetry",FALSE, etBOOL,{&bALLOW_ASYMMETRY},
3504                                 "Allow asymmetric insertion, i.e. the number of lipids removed from the upper and lower leaflet will not be checked." },
3505                         { "-ndiff" ,  FALSE, etINT, {&low_up_rm},
3506                                 "Number of lipids that will additionally be removed from the lower (negative number) or upper (positive number) membrane leaflet." },
3507                         { "-maxwarn", FALSE, etINT, {&maxwarn},
3508                                 "Maximum number of warning allowed" },
3509                         { "-start",   FALSE, etBOOL, {&bStart},
3510                                 "Call mdrun with membed options" },
3511                         { "-stepout", FALSE, etINT, {&nstepout},
3512                                 "HIDDENFrequency of writing the remaining runtime" },
3513                         { "-v",       FALSE, etBOOL,{&bVerbose},
3514                                 "Be loud and noisy" },
3515                         { "-mdrun_path", FALSE, etSTR, {&mdrun_path},
3516                                 "Path to the mdrun executable compiled with this g_membed version" }
3517         };
3518
3519         FILE *data_out;
3520         output_env_t oenv;
3521         char buf[256],buf2[64];
3522         gmx_edsam_t  ed;
3523         unsigned long Flags, PCA_Flags;
3524         ivec     ddxyz;
3525         int      dd_node_order;
3526         gmx_bool     HaveCheckpoint;
3527         FILE     *fplog,*fptest;
3528         int      sim_part,sim_part_fn;
3529         const char *part_suffix=".part";
3530         char     suffix[STRLEN];
3531         int      rc;
3532         char **multidir=NULL;
3533
3534         cr = init_par(&argc,&argv);
3535
3536         PCA_Flags = (PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_CAN_SET_DEFFNM
3537                         | (MASTER(cr) ? 0 : PCA_QUIET));
3538
3539
3540         /* Comment this in to do fexist calls only on master
3541          * works not with rerun or tables at the moment
3542          * also comment out the version of init_forcerec in md.c
3543          * with NULL instead of opt2fn
3544          */
3545         /*
3546    if (!MASTER(cr))
3547    {
3548    PCA_Flags |= PCA_NOT_READ_NODE;
3549    }
3550          */
3551
3552         parse_common_args(&argc,argv,PCA_Flags, NFILE,fnm,asize(pa),pa,
3553                         asize(desc),desc,0,NULL, &oenv);
3554
3555         /* we set these early because they might be used in init_multisystem()
3556    Note that there is the potential for npme>nnodes until the number of
3557    threads is set later on, if there's thread parallelization. That shouldn't
3558    lead to problems. */
3559         dd_node_order = nenum(ddno_opt);
3560         cr->npmenodes = npme;
3561
3562 #ifdef GMX_THREADS
3563         /* now determine the number of threads automatically. The threads are
3564    only started at mdrunner_threads, though. */
3565         if (nthreads<1)
3566         {
3567                 nthreads=tMPI_Thread_get_hw_number();
3568         }
3569 #else
3570         nthreads=1;
3571 #endif
3572
3573         /* now check the -multi and -multidir option */
3574         if (opt2bSet("-multidir", NFILE, fnm))
3575         {
3576             int i;
3577             if (nmultisim > 0)
3578             {
3579                 gmx_fatal(FARGS, "mdrun -multi and -multidir options are mutually     exclusive.");
3580             }
3581             nmultisim = opt2fns(&multidir, "-multidir", NFILE, fnm);
3582         }
3583
3584
3585
3586         parse_common_args(&argc,argv,0, NFILE,fnm,asize(pa),pa,
3587                     asize(desc),desc,0,NULL, &oenv);
3588         if (nmultisim > 1) {
3589 #ifndef GMX_THREADS
3590                 gmx_bool bParFn = (multidir == NULL);
3591                 init_multisystem(cr,nmultisim,multidir,NFILE,fnm,TRUE);
3592 #else
3593                 gmx_fatal(FARGS,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
3594 #endif
3595         }
3596
3597         /* Check if there is ANY checkpoint file available */
3598         sim_part    = 1;
3599         sim_part_fn = sim_part;
3600         if (opt2bSet("-cpi",NFILE,fnm))
3601         {
3602                 bAppendFiles =
3603                         read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE,fnm,cr),
3604                                                         &sim_part_fn,NULL,cr,
3605                                                         bAppendFiles,NFILE,fnm,
3606                                                         part_suffix,&bAddPart);
3607                 if (sim_part_fn==0 && MASTER(cr))
3608                 {
3609                         fprintf(stdout,"No previous checkpoint file present, assuming this is a new run.\n");
3610                 }
3611                 else
3612                 {
3613                         sim_part = sim_part_fn + 1;
3614                 }
3615         }
3616         else
3617         {
3618                 bAppendFiles = FALSE;
3619         }
3620
3621         data_out = ffopen(opt2fn("-dat",NFILE,fnm),"w");
3622         fprintf(data_out,"nxy = %d\nnz = %d\nxyinit = %f\nxyend = %f\nzinit = %f\nzend = %f\n"
3623                         "rad = %f\npieces = %d\nasymmetry = %s\nndiff = %d\nmaxwarn = %d\n",
3624                         it_xy,it_z,xy_fac,xy_max,z_fac,z_max,probe_rad,pieces,
3625                         bALLOW_ASYMMETRY ? "yes" : "no",low_up_rm,maxwarn);
3626         fclose(data_out);
3627
3628         sprintf(buf,"%s -s %s -membed %s -o %s -c %s -e %s -nt 1 -cpt -1",
3629                     (mdrun_path==NULL) ? "mdrun" : mdrun_path,
3630                     opt2fn("-f",NFILE,fnm),opt2fn("-dat",NFILE,fnm),opt2fn("-o",NFILE,fnm),
3631                     opt2fn("-c",NFILE,fnm),opt2fn("-e",NFILE,fnm));
3632         if (opt2bSet("-n",NFILE,fnm))
3633         {
3634                 sprintf(buf2," -mn %s",opt2fn("-n",NFILE,fnm));
3635                 strcat(buf,buf2);
3636         }
3637         if (opt2bSet("-x",NFILE,fnm))
3638         {
3639                 sprintf(buf2," -x %s",opt2fn("-x",NFILE,fnm));
3640                 strcat(buf,buf2);
3641         }
3642         if (opt2bSet("-p",NFILE,fnm))
3643         {
3644                 sprintf(buf2," -mp %s",opt2fn("-p",NFILE,fnm));
3645                 strcat(buf,buf2);
3646         }
3647         if (bVerbose)
3648         {
3649                 sprintf(buf2," -v -stepout %d",nstepout);
3650                 strcat(buf,buf2);
3651         }
3652
3653         printf("%s\n",buf);
3654         if (bStart)
3655         {
3656                 system(buf);
3657         } else {
3658                 printf("You can membed your protein now by:\n%s\n",buf);
3659         }
3660
3661         fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
3662
3663         return 0;
3664 }