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