src/tools/gmx_membed.c

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