src/tools/gmx_membed.c

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