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     bReadEkin = (flags & CGLO_READEKIN);
1160     bScaleEkin = (flags & CGLO_SCALEEKIN);
1161     bEner = flags & CGLO_ENERGY;
1162     bTemp = flags & CGLO_TEMPERATURE;
1163     bPres  = (flags & CGLO_PRESSURE);
1164     bConstrain = (flags & CGLO_CONSTRAINT);
1165     bIterate = (flags & CGLO_ITERATE);
1166     bFirstIterate = (flags & CGLO_FIRSTITERATE);
1167
1168     /* we calculate a full state kinetic energy either with full-step velocity verlet
1169        or half step where we need the pressure */
1170     bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir) && bPres) || bReadEkin);
1171
1172     /* in initalization, it sums the shake virial in vv, and to
1173        sums ekinh_old in leapfrog (or if we are calculating ekinh_old for other reasons */
1174
1175     /* ########## Kinetic energy  ############## */
1176
1177     if (bTemp)
1178     {
1179         /* Non-equilibrium MD: this is parallellized, but only does communication
1180          * when there really is NEMD.
1181          */
1182
1183         if (PAR(cr) && (ekind->bNEMD))
1184         {
1185             accumulate_u(cr,&(ir->opts),ekind);
1186         }
1187         debug_gmx();
1188         if (bReadEkin)
1189         {
1190             restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
1191         }
1192         else
1193         {
1194
1195             calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
1196         }
1197
1198         debug_gmx();
1199
1200         /* Calculate center of mass velocity if necessary, also parallellized */
1201         if (bStopCM && !bRerunMD && bEner)
1202         {
1203             calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
1204                          state->x,state->v,vcm);
1205         }
1206     }
1207
1208     if (bTemp || bPres || bEner || bConstrain)
1209     {
1210         if (!bGStat)
1211         {
1212             /* We will not sum ekinh_old,
1213              * so signal that we still have to do it.
1214              */
1215             *bSumEkinhOld = TRUE;
1216
1217         }
1218         else
1219         {
1220             if (gs != NULL)
1221             {
1222                 for(i=0; i<eglsNR; i++)
1223                 {
1224                     gs_buf[i] = gs->sig[i];
1225                 }
1226             }
1227             if (PAR(cr))
1228             {
1229                 wallcycle_start(wcycle,ewcMoveE);
1230                 GMX_MPE_LOG(ev_global_stat_start);
1231                 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
1232                             ir,ekind,constr,vcm,
1233                             gs != NULL ? eglsNR : 0,gs_buf,
1234                             top_global,state,
1235                             *bSumEkinhOld,flags);
1236                 GMX_MPE_LOG(ev_global_stat_finish);
1237                 wallcycle_stop(wcycle,ewcMoveE);
1238             }
1239             if (gs != NULL)
1240             {
1241                 if (MULTISIM(cr) && bInterSimGS)
1242                 {
1243                     if (MASTER(cr))
1244                     {
1245                         /* Communicate the signals between the simulations */
1246                         gmx_sum_sim(eglsNR,gs_buf,cr->ms);
1247                     }
1248                     /* Communicate the signals form the master to the others */
1249                     gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
1250                 }
1251                 for(i=0; i<eglsNR; i++)
1252                 {
1253                     if (bInterSimGS || gs_simlocal[i])
1254                     {
1255                         /* Set the communicated signal only when it is non-zero,
1256                          * since signals might not be processed at each MD step.
1257                          */
1258                         gsi = (gs_buf[i] >= 0 ?
1259                                (int)(gs_buf[i] + 0.5) :
1260                                (int)(gs_buf[i] - 0.5));
1261                         if (gsi != 0)
1262                         {
1263                             gs->set[i] = gsi;
1264                         }
1265                         /* Turn off the local signal */
1266                         gs->sig[i] = 0;
1267                     }
1268                 }
1269             }
1270             *bSumEkinhOld = FALSE;
1271         }
1272     }
1273
1274     if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
1275     {
1276         correct_ekin(debug,
1277                      mdatoms->start,mdatoms->start+mdatoms->homenr,
1278                      state->v,vcm->group_p[0],
1279                      mdatoms->massT,mdatoms->tmass,ekind->ekin);
1280     }
1281
1282     if (bEner) {
1283         /* Do center of mass motion removal */
1284         if (bStopCM && !bRerunMD) /* is this correct?  Does it get called too often with this logic? */
1285         {
1286             check_cm_grp(fplog,vcm,ir,1);
1287             do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
1288                           state->x,state->v,vcm);
1289             inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
1290         }
1291     }
1292
1293     if (bTemp)
1294     {
1295         /* Sum the kinetic energies of the groups & calc temp */
1296         /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
1297         /* three maincase:  VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
1298            Leap with AveVel is also an option for the future but not supported now.
1299            bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
1300            If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
1301            bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
1302            If FALSE, we go ahead and erase over it.
1303         */
1304         enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
1305                                        bEkinAveVel,bIterate,bScaleEkin);
1306
1307         enerd->term[F_EKIN] = trace(ekind->ekin);
1308     }
1309
1310     /* ##########  Long range energy information ###### */
1311
1312     if (bEner || bPres || bConstrain)
1313     {
1314         calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda,
1315                       corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
1316     }
1317
1318     if (bEner && bFirstIterate)
1319     {
1320         enerd->term[F_DISPCORR] = enercorr;
1321         enerd->term[F_EPOT] += enercorr;
1322         enerd->term[F_DVDL] += dvdlcorr;
1323         if (fr->efep != efepNO) {
1324             enerd->dvdl_lin += dvdlcorr;
1325         }
1326     }
1327
1328     /* ########## Now pressure ############## */
1329     if (bPres || bConstrain)
1330     {
1331
1332         m_add(force_vir,shake_vir,total_vir);
1333
1334         /* Calculate pressure and apply LR correction if PPPM is used.
1335          * Use the box from last timestep since we already called update().
1336          */
1337
1338         enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
1339                                         (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
1340
1341         /* Calculate long range corrections to pressure and energy */
1342         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
1343            and computes enerd->term[F_DISPCORR].  Also modifies the
1344            total_vir and pres tesors */
1345
1346         m_add(total_vir,corr_vir,total_vir);
1347         m_add(pres,corr_pres,pres);
1348         enerd->term[F_PDISPCORR] = prescorr;
1349         enerd->term[F_PRES] += prescorr;
1350         *pcurr = enerd->term[F_PRES];
1351         /* calculate temperature using virial */
1352         enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
1353
1354     }
1355 }
1356
1357
1358 /* Definitions for convergence of iterated constraints */
1359
1360 /* iterate constraints up to 50 times  */
1361 #define MAXITERCONST       50
1362
1363 /* data type */
1364 typedef struct
1365 {
1366     real f,fprev,x,xprev;
1367     int iter_i;
1368     gmx_bool bIterate;
1369     real allrelerr[MAXITERCONST+2];
1370     int num_close; /* number of "close" violations, caused by limited precision. */
1371 } gmx_iterate_t;
1372
1373 #ifdef GMX_DOUBLE
1374 #define CONVERGEITER  0.000000001
1375 #define CLOSE_ENOUGH  0.000001000
1376 #else
1377 #define CONVERGEITER  0.0001
1378 #define CLOSE_ENOUGH  0.0050
1379 #endif
1380
1381 /* we want to keep track of the close calls.  If there are too many, there might be some other issues.
1382    so we make sure that it's either less than some predetermined number, or if more than that number,
1383    only some small fraction of the total. */
1384 #define MAX_NUMBER_CLOSE        50
1385 #define FRACTION_CLOSE       0.001
1386
1387 /* maximum length of cyclic traps to check, emerging from limited numerical precision  */
1388 #define CYCLEMAX            20
1389
1390 static void gmx_iterate_init(gmx_iterate_t *iterate,gmx_bool bIterate)
1391 {
1392     int i;
1393
1394     iterate->iter_i = 0;
1395     iterate->bIterate = bIterate;
1396     iterate->num_close = 0;
1397     for (i=0;i<MAXITERCONST+2;i++)
1398     {
1399         iterate->allrelerr[i] = 0;
1400     }
1401 }
1402
1403 static gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate_t *iterate, gmx_bool bFirstIterate, real fom, real *newf)
1404 {
1405     /* monitor convergence, and use a secant search to propose new
1406        values.
1407                                                                   x_{i} - x_{i-1}
1408        The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
1409                                                                 f(x_{i}) - f(x_{i-1})
1410
1411        The function we are trying to zero is fom-x, where fom is the
1412        "figure of merit" which is the pressure (or the veta value) we
1413        would get by putting in an old value of the pressure or veta into
1414        the incrementor function for the step or half step.  I have
1415        verified that this gives the same answer as self consistent
1416        iteration, usually in many fewer steps, especially for small tau_p.
1417
1418        We could possibly eliminate an iteration with proper use
1419        of the value from the previous step, but that would take a bit
1420        more bookkeeping, especially for veta, since tests indicate the
1421        function of veta on the last step is not sufficiently close to
1422        guarantee convergence this step. This is
1423        good enough for now.  On my tests, I could use tau_p down to
1424        0.02, which is smaller that would ever be necessary in
1425        practice. Generally, 3-5 iterations will be sufficient */
1426
1427     real relerr,err;
1428     char buf[256];
1429     int i;
1430     gmx_bool incycle;
1431
1432     if (bFirstIterate)
1433     {
1434         iterate->x = fom;
1435         iterate->f = fom-iterate->x;
1436         iterate->xprev = 0;
1437         iterate->fprev = 0;
1438         *newf = fom;
1439     }
1440     else
1441     {
1442         iterate->f = fom-iterate->x; /* we want to zero this difference */
1443         if ((iterate->iter_i > 1) && (iterate->iter_i < MAXITERCONST))
1444         {
1445             if (iterate->f==iterate->fprev)
1446             {
1447                 *newf = iterate->f;
1448             }
1449             else
1450             {
1451                 *newf = iterate->x - (iterate->x-iterate->xprev)*(iterate->f)/(iterate->f-iterate->fprev);
1452             }
1453         }
1454         else
1455         {
1456             /* just use self-consistent iteration the first step to initialize, or
1457                if it's not converging (which happens occasionally -- need to investigate why) */
1458             *newf = fom;
1459         }
1460     }
1461     /* Consider a slight shortcut allowing us to exit one sooner -- we check the
1462        difference between the closest of x and xprev to the new
1463        value. To be 100% certain, we should check the difference between
1464        the last result, and the previous result, or
1465
1466        relerr = (fabs((x-xprev)/fom));
1467
1468        but this is pretty much never necessary under typical conditions.
1469        Checking numerically, it seems to lead to almost exactly the same
1470        trajectories, but there are small differences out a few decimal
1471        places in the pressure, and eventually in the v_eta, but it could
1472        save an interation.
1473
1474        if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
1475        relerr = (fabs((*newf-xmin) / *newf));
1476     */
1477
1478     err = fabs((iterate->f-iterate->fprev));
1479     relerr = fabs(err/fom);
1480
1481     iterate->allrelerr[iterate->iter_i] = relerr;
1482
1483     if (iterate->iter_i > 0)
1484     {
1485         if (debug)
1486         {
1487             fprintf(debug,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
1488                     iterate->iter_i,fom,relerr,*newf);
1489         }
1490
1491         if ((relerr < CONVERGEITER) || (err < CONVERGEITER) || (fom==0) || ((iterate->x == iterate->xprev) && iterate->iter_i > 1))
1492         {
1493             iterate->bIterate = FALSE;
1494             if (debug)
1495             {
1496                 fprintf(debug,"Iterating NPT constraints: CONVERGED\n");
1497             }
1498             return TRUE;
1499         }
1500         if (iterate->iter_i > MAXITERCONST)
1501         {
1502             if (relerr < CLOSE_ENOUGH)
1503             {
1504                 incycle = FALSE;
1505                 for (i=1;i<CYCLEMAX;i++) {
1506                     if ((iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-1]) &&
1507                         (iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-(1+2*i)])) {
1508                         incycle = TRUE;
1509                         if (debug)
1510                         {
1511                             fprintf(debug,"Exiting from an NPT iterating cycle of length %d\n",i);
1512                         }
1513                         break;
1514                     }
1515                 }
1516
1517                 if (incycle) {
1518                     /* step 1: trapped in a numerical attractor */
1519                     /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
1520                        Better to give up convergence here than have the simulation die.
1521                     */
1522                     iterate->num_close++;
1523                     return TRUE;
1524                 }
1525                 else
1526                 {
1527                     /* Step #2: test if we are reasonably close for other reasons, then monitor the number.  If not, die */
1528
1529                     /* how many close calls have we had?  If less than a few, we're OK */
1530                     if (iterate->num_close < MAX_NUMBER_CLOSE)
1531                     {
1532                         sprintf(buf,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps,relerr);
1533                         md_print_warning(cr,fplog,buf);
1534                         iterate->num_close++;
1535                         return TRUE;
1536                         /* if more than a few, check the total fraction.  If too high, die. */
1537                     } else if (iterate->num_close/(double)nsteps > FRACTION_CLOSE) {
1538                         gmx_fatal(FARGS,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate->num_close/(double)nsteps);
1539                     }
1540                 }
1541             }
1542             else
1543             {
1544                 gmx_fatal(FARGS,"Could not converge NPT constraints\n");
1545             }
1546         }
1547     }
1548
1549     iterate->xprev = iterate->x;
1550     iterate->x = *newf;
1551     iterate->fprev = iterate->f;
1552     iterate->iter_i++;
1553
1554     return FALSE;
1555 }
1556
1557 static void check_nst_param(FILE *fplog,t_commrec *cr,
1558                             const char *desc_nst,int nst,
1559                             const char *desc_p,int *p)
1560 {
1561     char buf[STRLEN];
1562
1563     if (*p > 0 && *p % nst != 0)
1564     {
1565         /* Round up to the next multiple of nst */
1566         *p = ((*p)/nst + 1)*nst;
1567         sprintf(buf,"NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
1568         md_print_warning(cr,fplog,buf);
1569     }
1570 }
1571
1572 static void reset_all_counters(FILE *fplog,t_commrec *cr,
1573                                gmx_large_int_t step,
1574                                gmx_large_int_t *step_rel,t_inputrec *ir,
1575                                gmx_wallcycle_t wcycle,t_nrnb *nrnb,
1576                                gmx_runtime_t *runtime)
1577 {
1578     char buf[STRLEN],sbuf[STEPSTRSIZE];
1579
1580     /* Reset all the counters related to performance over the run */
1581     sprintf(buf,"Step %s: resetting all time and cycle counters\n",
1582             gmx_step_str(step,sbuf));
1583     md_print_warning(cr,fplog,buf);
1584
1585     wallcycle_stop(wcycle,ewcRUN);
1586     wallcycle_reset_all(wcycle);
1587     if (DOMAINDECOMP(cr))
1588     {
1589         reset_dd_statistics_counters(cr->dd);
1590     }
1591     init_nrnb(nrnb);
1592     ir->init_step += *step_rel;
1593     ir->nsteps    -= *step_rel;
1594     *step_rel = 0;
1595     wallcycle_start(wcycle,ewcRUN);
1596     runtime_start(runtime);
1597     print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
1598 }
1599
1600 static int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
1601                                int nstglobalcomm,t_inputrec *ir)
1602 {
1603     char buf[STRLEN];
1604
1605     if (!EI_DYNAMICS(ir->eI))
1606     {
1607         nstglobalcomm = 1;
1608     }
1609
1610     if (nstglobalcomm == -1)
1611     {
1612         if (ir->nstcalcenergy == 0 && ir->nstlist == 0)
1613         {
1614             nstglobalcomm = 10;
1615             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
1616             {
1617                 nstglobalcomm = ir->nstenergy;
1618             }
1619         }
1620         else
1621         {
1622             /* We assume that if nstcalcenergy > nstlist,
1623              * nstcalcenergy is a multiple of nstlist.
1624              */
1625             if (ir->nstcalcenergy == 0 ||
1626                 (ir->nstlist > 0 && ir->nstlist < ir->nstcalcenergy))
1627             {
1628                 nstglobalcomm = ir->nstlist;
1629             }
1630             else
1631             {
1632                 nstglobalcomm = ir->nstcalcenergy;
1633             }
1634         }
1635     }
1636     else
1637     {
1638         if (ir->nstlist > 0 &&
1639             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
1640         {
1641             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
1642             sprintf(buf,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
1643             md_print_warning(cr,fplog,buf);
1644         }
1645         if (nstglobalcomm > ir->nstcalcenergy)
1646         {
1647             check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
1648                             "nstcalcenergy",&ir->nstcalcenergy);
1649         }
1650
1651         check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
1652                         "nstenergy",&ir->nstenergy);
1653
1654         check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
1655                         "nstlog",&ir->nstlog);
1656     }
1657
1658     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
1659     {
1660         sprintf(buf,"WARNING: Changing nstcomm from %d to %d\n",
1661                 ir->nstcomm,nstglobalcomm);
1662         md_print_warning(cr,fplog,buf);
1663         ir->nstcomm = nstglobalcomm;
1664     }
1665
1666     return nstglobalcomm;
1667 }
1668
1669 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
1670                                t_inputrec *ir,gmx_mtop_t *mtop)
1671 {
1672     /* Check required for old tpx files */
1673     if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
1674         ir->nstcalcenergy % ir->nstlist != 0)
1675     {
1676         md_print_warning(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
1677
1678         if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
1679             gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
1680             ir->eConstrAlg == econtSHAKE)
1681         {
1682             md_print_warning(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
1683             if (ir->epc != epcNO)
1684             {
1685                 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
1686             }
1687         }
1688         check_nst_param(fplog,cr,"nstlist",ir->nstlist,
1689                         "nstcalcenergy",&ir->nstcalcenergy);
1690         check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
1691                         "nstenergy",&ir->nstenergy);
1692         check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
1693                         "nstlog",&ir->nstlog);
1694         if (ir->efep != efepNO)
1695         {
1696             check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
1697                             "nstdhdl",&ir->nstdhdl);
1698         }
1699     }
1700 }
1701
1702 typedef struct {
1703     gmx_bool       bGStatEveryStep;
1704     gmx_large_int_t step_ns;
1705     gmx_large_int_t step_nscheck;
1706     gmx_large_int_t nns;
1707     matrix     scale_tot;
1708     int        nabnsb;
1709     double     s1;
1710     double     s2;
1711     double     ab;
1712     double     lt_runav;
1713     double     lt_runav2;
1714 } gmx_nlheur_t;
1715
1716 static void reset_nlistheuristics(gmx_nlheur_t *nlh,gmx_large_int_t step)
1717 {
1718     nlh->lt_runav  = 0;
1719     nlh->lt_runav2 = 0;
1720     nlh->step_nscheck = step;
1721 }
1722
1723 static void init_nlistheuristics(gmx_nlheur_t *nlh,
1724                                  gmx_bool bGStatEveryStep,gmx_large_int_t step)
1725 {
1726     nlh->bGStatEveryStep = bGStatEveryStep;
1727     nlh->nns       = 0;
1728     nlh->nabnsb    = 0;
1729     nlh->s1        = 0;
1730     nlh->s2        = 0;
1731     nlh->ab        = 0;
1732
1733     reset_nlistheuristics(nlh,step);
1734 }
1735
1736 static void update_nliststatistics(gmx_nlheur_t *nlh,gmx_large_int_t step)
1737 {
1738     gmx_large_int_t nl_lt;
1739     char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1740
1741     /* Determine the neighbor list life time */
1742     nl_lt = step - nlh->step_ns;
1743     if (debug)
1744     {
1745         fprintf(debug,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh->nabnsb,gmx_step_str(nl_lt,sbuf));
1746     }
1747     nlh->nns++;
1748     nlh->s1 += nl_lt;
1749     nlh->s2 += nl_lt*nl_lt;
1750     nlh->ab += nlh->nabnsb;
1751     if (nlh->lt_runav == 0)
1752     {
1753         nlh->lt_runav  = nl_lt;
1754         /* Initialize the fluctuation average
1755          * such that at startup we check after 0 steps.
1756          */
1757         nlh->lt_runav2 = sqr(nl_lt/2.0);
1758     }
1759     /* Running average with 0.9 gives an exp. history of 9.5 */
1760     nlh->lt_runav2 = 0.9*nlh->lt_runav2 + 0.1*sqr(nlh->lt_runav - nl_lt);
1761     nlh->lt_runav  = 0.9*nlh->lt_runav  + 0.1*nl_lt;
1762     if (nlh->bGStatEveryStep)
1763     {
1764         /* Always check the nlist validity */
1765         nlh->step_nscheck = step;
1766     }
1767     else
1768     {
1769         /* We check after:  <life time> - 2*sigma
1770          * The factor 2 is quite conservative,
1771          * but we assume that with nstlist=-1 the user
1772          * prefers exact integration over performance.
1773          */
1774         nlh->step_nscheck = step
1775                   + (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)) - 1;
1776     }
1777     if (debug)
1778     {
1779         fprintf(debug,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
1780                 gmx_step_str(nl_lt,sbuf),nlh->lt_runav,sqrt(nlh->lt_runav2),
1781                 gmx_step_str(nlh->step_nscheck-step+1,sbuf2),
1782                 (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)));
1783     }
1784 }
1785
1786 static void set_nlistheuristics(gmx_nlheur_t *nlh,gmx_bool bReset,gmx_large_int_t step)
1787 {
1788     int d;
1789
1790     if (bReset)
1791     {
1792         reset_nlistheuristics(nlh,step);
1793     }
1794     else
1795     {
1796         update_nliststatistics(nlh,step);
1797     }
1798
1799     nlh->step_ns = step;
1800     /* Initialize the cumulative coordinate scaling matrix */
1801     clear_mat(nlh->scale_tot);
1802     for(d=0; d<DIM; d++)
1803     {
1804         nlh->scale_tot[d][d] = 1.0;
1805     }
1806 }
1807
1808 double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
1809              const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
1810              int nstglobalcomm,
1811              gmx_vsite_t *vsite,gmx_constr_t constr,
1812              int stepout,t_inputrec *ir,
1813              gmx_mtop_t *top_global,
1814              t_fcdata *fcd,
1815              t_state *state_global,
1816              t_mdatoms *mdatoms,
1817              t_nrnb *nrnb,gmx_wallcycle_t wcycle,
1818              gmx_edsam_t ed,t_forcerec *fr,
1819              int repl_ex_nst,int repl_ex_seed,
1820              real cpt_period,real max_hours,
1821              const char *deviceOptions,
1822              unsigned long Flags,
1823              gmx_runtime_t *runtime,
1824              rvec fac, rvec *r_ins, pos_ins_t *pos_ins, t_block *ins_at,
1825              real xy_step, real z_step, int it_xy, int it_z)
1826 {
1827     gmx_mdoutf_t *outf;
1828     gmx_large_int_t step,step_rel;
1829     double     run_time;
1830     double     t,t0,lam0;
1831     gmx_bool       bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
1832     gmx_bool       bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
1833                bFirstStep,bStateFromTPX,bInitStep,bLastStep,
1834                bBornRadii,bStartingFromCpt;
1835     gmx_bool       bDoDHDL=FALSE;
1836     gmx_bool       do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
1837                bForceUpdate=FALSE,bCPT;
1838     int        mdof_flags;
1839     gmx_bool       bMasterState;
1840     int        force_flags,cglo_flags;
1841     tensor     force_vir,shake_vir,total_vir,tmp_vir,pres;
1842     int        i,m;
1843     t_trxstatus *status;
1844     rvec       mu_tot;
1845     t_vcm      *vcm;
1846     t_state    *bufstate=NULL;
1847     matrix     *scale_tot,pcoupl_mu,M,ebox;
1848     gmx_nlheur_t nlh;
1849     t_trxframe rerun_fr;
1850 /*    gmx_repl_ex_t repl_ex=NULL;*/
1851     int        nchkpt=1;
1852
1853     gmx_localtop_t *top;
1854     t_mdebin *mdebin=NULL;
1855     t_state    *state=NULL;
1856     rvec       *f_global=NULL;
1857     int        n_xtc=-1;
1858     rvec       *x_xtc=NULL;
1859     gmx_enerdata_t *enerd;
1860     rvec       *f=NULL;
1861     gmx_global_stat_t gstat;
1862     gmx_update_t upd=NULL;
1863     t_graph    *graph=NULL;
1864     globsig_t   gs;
1865
1866     gmx_bool        bFFscan;
1867     gmx_groups_t *groups;
1868     gmx_ekindata_t *ekind, *ekind_save;
1869     gmx_shellfc_t shellfc;
1870     int         count,nconverged=0;
1871     real        timestep=0;
1872     double      tcount=0;
1873     gmx_bool        bIonize=FALSE;
1874     gmx_bool        bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
1875     gmx_bool        bAppend;
1876     gmx_bool        bResetCountersHalfMaxH=FALSE;
1877     gmx_bool        bVV,bIterations,bIterate,bFirstIterate,bTemp,bPres,bTrotter;
1878     real        temp0,dvdl;
1879     int         a0,a1,ii;
1880     rvec        *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
1881     matrix      boxcopy={{0}},lastbox;
1882         real        veta_save,pcurr,scalevir,tracevir;
1883         real        vetanew = 0;
1884     double      cycles;
1885         real        last_conserved = 0;
1886     real        last_ekin = 0;
1887         t_extmass   MassQ;
1888     int         **trotter_seq;
1889     char        sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1890     int         handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
1891     gmx_iterate_t iterate;
1892 #ifdef GMX_FAHCORE
1893     /* Temporary addition for FAHCORE checkpointing */
1894     int chkpt_ret;
1895 #endif
1896
1897     /* Check for special mdrun options */
1898     bRerunMD = (Flags & MD_RERUN);
1899     bIonize  = (Flags & MD_IONIZE);
1900     bFFscan  = (Flags & MD_FFSCAN);
1901     bAppend  = (Flags & MD_APPENDFILES);
1902     bGStatEveryStep = FALSE;
1903     if (Flags & MD_RESETCOUNTERSHALFWAY)
1904     {
1905         if (ir->nsteps > 0)
1906         {
1907             /* Signal to reset the counters half the simulation steps. */
1908             wcycle_set_reset_counters(wcycle,ir->nsteps/2);
1909         }
1910         /* Signal to reset the counters halfway the simulation time. */
1911         bResetCountersHalfMaxH = (max_hours > 0);
1912     }
1913
1914     /* md-vv uses averaged full step velocities for T-control
1915        md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1916        md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1917     bVV = EI_VV(ir->eI);
1918     if (bVV) /* to store the initial velocities while computing virial */
1919     {
1920         snew(cbuf,top_global->natoms);
1921     }
1922     /* all the iteratative cases - only if there are constraints */
1923     bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
1924     bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
1925
1926     if (bRerunMD)
1927     {
1928         /* Since we don't know if the frames read are related in any way,
1929          * rebuild the neighborlist at every step.
1930          */
1931         ir->nstlist       = 1;
1932         ir->nstcalcenergy = 1;
1933         nstglobalcomm     = 1;
1934     }
1935
1936     check_ir_old_tpx_versions(cr,fplog,ir,top_global);
1937
1938     nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
1939     /*bGStatEveryStep = (nstglobalcomm == 1);*/
1940     bGStatEveryStep = FALSE;
1941
1942     if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
1943     {
1944         fprintf(fplog,
1945                 "To reduce the energy communication with nstlist = -1\n"
1946                 "the neighbor list validity should not be checked at every step,\n"
1947                 "this means that exact integration is not guaranteed.\n"
1948                 "The neighbor list validity is checked after:\n"
1949                 "  <n.list life time> - 2*std.dev.(n.list life time)  steps.\n"
1950                 "In most cases this will result in exact integration.\n"
1951                 "This reduces the energy communication by a factor of 2 to 3.\n"
1952                 "If you want less energy communication, set nstlist > 3.\n\n");
1953     }
1954
1955     if (bRerunMD || bFFscan)
1956     {
1957         ir->nstxtcout = 0;
1958     }
1959     groups = &top_global->groups;
1960
1961     /* Initial values */
1962     init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
1963             nrnb,top_global,&upd,
1964             nfile,fnm,&outf,&mdebin,
1965             force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
1966
1967     clear_mat(total_vir);
1968     clear_mat(pres);
1969     /* Energy terms and groups */
1970     snew(enerd,1);
1971     init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
1972     if (DOMAINDECOMP(cr))
1973     {
1974         f = NULL;
1975     }
1976     else
1977     {
1978         snew(f,top_global->natoms);
1979     }
1980
1981     /* Kinetic energy data */
1982     snew(ekind,1);
1983     init_ekindata(fplog,top_global,&(ir->opts),ekind);
1984     /* needed for iteration of constraints */
1985     snew(ekind_save,1);
1986     init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
1987     /* Copy the cos acceleration to the groups struct */
1988     ekind->cosacc.cos_accel = ir->cos_accel;
1989
1990     gstat = global_stat_init(ir);
1991     debug_gmx();
1992
1993     /* Check for polarizable models and flexible constraints */
1994     shellfc = init_shell_flexcon(fplog,
1995                                  top_global,n_flexible_constraints(constr),
1996                                  (ir->bContinuation ||
1997                                   (DOMAINDECOMP(cr) && !MASTER(cr))) ?
1998                                  NULL : state_global->x);
1999
2000 /*    if (DEFORM(*ir))
2001     {
2002 #ifdef GMX_THREADS
2003         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
2004 #endif
2005         set_deform_reference_box(upd,
2006                                  deform_init_init_step_tpx,
2007                                  deform_init_box_tpx);
2008 #ifdef GMX_THREADS
2009         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
2010 #endif
2011     }*/
2012
2013 /*    {
2014         double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
2015         if ((io > 2000) && MASTER(cr))
2016             fprintf(stderr,
2017                     "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
2018                     io);
2019     }*/
2020
2021     if (DOMAINDECOMP(cr)) {
2022         top = dd_init_local_top(top_global);
2023
2024         snew(state,1);
2025         dd_init_local_state(cr->dd,state_global,state);
2026
2027         if (DDMASTER(cr->dd) && ir->nstfout) {
2028             snew(f_global,state_global->natoms);
2029         }
2030     } else {
2031         if (PAR(cr)) {
2032             /* Initialize the particle decomposition and split the topology */
2033             top = split_system(fplog,top_global,ir,cr);
2034
2035             pd_cg_range(cr,&fr->cg0,&fr->hcg);
2036             pd_at_range(cr,&a0,&a1);
2037         } else {
2038             top = gmx_mtop_generate_local_top(top_global,ir);
2039
2040             a0 = 0;
2041             a1 = top_global->natoms;
2042         }
2043
2044         state = partdec_init_local_state(cr,state_global);
2045         f_global = f;
2046
2047         atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
2048
2049         if (vsite) {
2050             set_vsite_top(vsite,top,mdatoms,cr);
2051         }
2052
2053         if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
2054             graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
2055         }
2056
2057         if (shellfc) {
2058             make_local_shells(cr,mdatoms,shellfc);
2059         }
2060
2061         if (ir->pull && PAR(cr)) {
2062             dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
2063         }
2064     }
2065
2066     if (DOMAINDECOMP(cr))
2067     {
2068         /* Distribute the charge groups over the nodes from the master node */
2069         dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
2070                             state_global,top_global,ir,
2071                             state,&f,mdatoms,top,fr,
2072                             vsite,shellfc,constr,
2073                             nrnb,wcycle,FALSE);
2074     }
2075
2076     update_mdatoms(mdatoms,state->lambda);
2077
2078     if (MASTER(cr))
2079     {
2080         /* Update mdebin with energy history if appending to output files */
2081         if ( Flags & MD_APPENDFILES )
2082         {
2083             restore_energyhistory_from_state(mdebin,&state_global->enerhist);
2084         }
2085         /* Set the initial energy history in state to zero by updating once */
2086         update_energyhistory(&state_global->enerhist,mdebin);
2087     }
2088
2089     if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
2090         /* Set the random state if we read a checkpoint file */
2091         set_stochd_state(upd,state);
2092     }
2093
2094     /* Initialize constraints */
2095     if (constr) {
2096         if (!DOMAINDECOMP(cr))
2097             set_constraints(constr,top,ir,mdatoms,cr);
2098     }
2099
2100     /* Check whether we have to GCT stuff */
2101  /*   bTCR = ftp2bSet(efGCT,nfile,fnm);
2102     if (bTCR) {
2103         if (MASTER(cr)) {
2104             fprintf(stderr,"Will do General Coupling Theory!\n");
2105         }
2106         gnx = top_global->mols.nr;
2107         snew(grpindex,gnx);
2108         for(i=0; (i<gnx); i++) {
2109             grpindex[i] = i;
2110         }
2111     }*/
2112
2113 /*    if (repl_ex_nst > 0 && MASTER(cr))
2114         repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
2115                                         repl_ex_nst,repl_ex_seed);*/
2116
2117     if (!ir->bContinuation && !bRerunMD)
2118     {
2119         if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
2120         {
2121             /* Set the velocities of frozen particles to zero */
2122             for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
2123             {
2124                 for(m=0; m<DIM; m++)
2125                 {
2126                     if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
2127                     {
2128                         state->v[i][m] = 0;
2129                     }
2130                 }
2131             }
2132         }
2133
2134         if (constr)
2135         {
2136             /* Constrain the initial coordinates and velocities */
2137             do_constrain_first(fplog,constr,ir,mdatoms,state,f,
2138                                graph,cr,nrnb,fr,top,shake_vir);
2139         }
2140         if (vsite)
2141         {
2142             /* Construct the virtual sites for the initial configuration */
2143             construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
2144                              top->idef.iparams,top->idef.il,
2145                              fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2146         }
2147     }
2148
2149     debug_gmx();
2150
2151     /* I'm assuming we need global communication the first time! MRS */
2152     cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
2153                   | (bVV ? CGLO_PRESSURE:0)
2154                   | (bVV ? CGLO_CONSTRAINT:0)
2155                   | (bRerunMD ? CGLO_RERUNMD:0)
2156                   | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
2157
2158     bSumEkinhOld = FALSE;
2159     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2160                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2161                     constr,NULL,FALSE,state->box,
2162                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
2163     if (ir->eI == eiVVAK) {
2164         /* a second call to get the half step temperature initialized as well */
2165         /* we do the same call as above, but turn the pressure off -- internally, this
2166            is recognized as a velocity verlet half-step kinetic energy calculation.
2167            This minimized excess variables, but perhaps loses some logic?*/
2168
2169         compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2170                         wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2171                         constr,NULL,FALSE,state->box,
2172                         top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2173                         cglo_flags &~ CGLO_PRESSURE);
2174     }
2175
2176     /* Calculate the initial half step temperature, and save the ekinh_old */
2177     if (!(Flags & MD_STARTFROMCPT))
2178     {
2179         for(i=0; (i<ir->opts.ngtc); i++)
2180         {
2181             copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
2182         }
2183     }
2184     if (ir->eI != eiVV)
2185     {
2186         enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
2187                                      and there is no previous step */
2188     }
2189     temp0 = enerd->term[F_TEMP];
2190
2191     /* if using an iterative algorithm, we need to create a working directory for the state. */
2192     if (bIterations)
2193     {
2194             bufstate = init_bufstate(state);
2195     }
2196     if (bFFscan)
2197     {
2198         snew(xcopy,state->natoms);
2199         snew(vcopy,state->natoms);
2200         copy_rvecn(state->x,xcopy,0,state->natoms);
2201         copy_rvecn(state->v,vcopy,0,state->natoms);
2202         copy_mat(state->box,boxcopy);
2203     }
2204
2205     /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
2206        temperature control */
2207     trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
2208
2209     if (MASTER(cr))
2210     {
2211         if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
2212         {
2213             fprintf(fplog,
2214                     "RMS relative constraint deviation after constraining: %.2e\n",
2215                     constr_rmsd(constr,FALSE));
2216         }
2217         fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
2218         if (bRerunMD)
2219         {
2220             fprintf(stderr,"starting md rerun '%s', reading coordinates from"
2221                     " input trajectory '%s'\n\n",
2222                     *(top_global->name),opt2fn("-rerun",nfile,fnm));
2223             if (bVerbose)
2224             {
2225                 fprintf(stderr,"Calculated time to finish depends on nsteps from "
2226                         "run input file,\nwhich may not correspond to the time "
2227                         "needed to process input trajectory.\n\n");
2228             }
2229         }
2230         else
2231         {
2232             char tbuf[20];
2233             fprintf(stderr,"starting mdrun '%s'\n",
2234                     *(top_global->name));
2235             if (ir->nsteps >= 0)
2236             {
2237                 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
2238             }
2239             else
2240             {
2241                 sprintf(tbuf,"%s","infinite");
2242             }
2243             if (ir->init_step > 0)
2244             {
2245                 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
2246                         gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
2247                         gmx_step_str(ir->init_step,sbuf2),
2248                         ir->init_step*ir->delta_t);
2249             }
2250             else
2251             {
2252                 fprintf(stderr,"%s steps, %s ps.\n",
2253                         gmx_step_str(ir->nsteps,sbuf),tbuf);
2254             }
2255         }
2256         fprintf(fplog,"\n");
2257     }
2258
2259     /* Set and write start time */
2260     runtime_start(runtime);
2261     print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
2262     wallcycle_start(wcycle,ewcRUN);
2263     if (fplog)
2264         fprintf(fplog,"\n");
2265
2266     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
2267 /*#ifdef GMX_FAHCORE
2268     chkpt_ret=fcCheckPointParallel( cr->nodeid,
2269                                     NULL,0);
2270     if ( chkpt_ret == 0 )
2271         gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
2272 #endif*/
2273
2274     debug_gmx();
2275     /***********************************************************
2276      *
2277      *             Loop over MD steps
2278      *
2279      ************************************************************/
2280
2281     /* if rerunMD then read coordinates and velocities from input trajectory */
2282     if (bRerunMD)
2283     {
2284         if (getenv("GMX_FORCE_UPDATE"))
2285         {
2286             bForceUpdate = TRUE;
2287         }
2288
2289         bNotLastFrame = read_first_frame(oenv,&status,
2290                                          opt2fn("-rerun",nfile,fnm),
2291                                          &rerun_fr,TRX_NEED_X | TRX_READ_V);
2292         if (rerun_fr.natoms != top_global->natoms)
2293         {
2294             gmx_fatal(FARGS,
2295                       "Number of atoms in trajectory (%d) does not match the "
2296                       "run input file (%d)\n",
2297                       rerun_fr.natoms,top_global->natoms);
2298         }
2299         if (ir->ePBC != epbcNONE)
2300         {
2301             if (!rerun_fr.bBox)
2302             {
2303                 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);
2304             }
2305             if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
2306             {
2307                 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
2308             }
2309
2310             /* Set the shift vectors.
2311              * Necessary here when have a static box different from the tpr box.
2312              */
2313             calc_shifts(rerun_fr.box,fr->shift_vec);
2314         }
2315     }
2316
2317     /* loop over MD steps or if rerunMD to end of input trajectory */
2318     bFirstStep = TRUE;
2319     /* Skip the first Nose-Hoover integration when we get the state from tpx */
2320     bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
2321     bInitStep = bFirstStep && (bStateFromTPX || bVV);
2322     bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
2323     bLastStep    = FALSE;
2324     bSumEkinhOld = FALSE;
2325     bExchanged   = FALSE;
2326
2327     init_global_signals(&gs,cr,ir,repl_ex_nst);
2328
2329     step = ir->init_step;
2330     step_rel = 0;
2331
2332     if (ir->nstlist == -1)
2333     {
2334         init_nlistheuristics(&nlh,bGStatEveryStep,step);
2335     }
2336
2337     bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
2338     while (!bLastStep || (bRerunMD && bNotLastFrame)) {
2339
2340         wallcycle_start(wcycle,ewcSTEP);
2341
2342         GMX_MPE_LOG(ev_timestep1);
2343
2344         if (bRerunMD) {
2345             if (rerun_fr.bStep) {
2346                 step = rerun_fr.step;
2347                 step_rel = step - ir->init_step;
2348             }
2349             if (rerun_fr.bTime) {
2350                 t = rerun_fr.time;
2351             }
2352             else
2353             {
2354                 t = step;
2355             }
2356         }
2357         else
2358         {
2359             bLastStep = (step_rel == ir->nsteps);
2360             t = t0 + step*ir->delta_t;
2361         }
2362
2363         if (ir->efep != efepNO)
2364         {
2365             if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
2366             {
2367                 state_global->lambda = rerun_fr.lambda;
2368             }
2369             else
2370             {
2371                 state_global->lambda = lam0 + step*ir->delta_lambda;
2372             }
2373             state->lambda = state_global->lambda;
2374             bDoDHDL = do_per_step(step,ir->nstdhdl);
2375         }
2376
2377         if (bSimAnn)
2378         {
2379             update_annealing_target_temp(&(ir->opts),t);
2380         }
2381
2382         if (bRerunMD)
2383         {
2384             if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
2385             {
2386                 for(i=0; i<state_global->natoms; i++)
2387                 {
2388                     copy_rvec(rerun_fr.x[i],state_global->x[i]);
2389                 }
2390                 if (rerun_fr.bV)
2391                 {
2392                     for(i=0; i<state_global->natoms; i++)
2393                     {
2394                         copy_rvec(rerun_fr.v[i],state_global->v[i]);
2395                     }
2396                 }
2397                 else
2398                 {
2399                     for(i=0; i<state_global->natoms; i++)
2400                     {
2401                         clear_rvec(state_global->v[i]);
2402                     }
2403                     if (bRerunWarnNoV)
2404                     {
2405                         fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
2406                                 "         Ekin, temperature and pressure are incorrect,\n"
2407                                 "         the virial will be incorrect when constraints are present.\n"
2408                                 "\n");
2409                         bRerunWarnNoV = FALSE;
2410                     }
2411                 }
2412             }
2413             copy_mat(rerun_fr.box,state_global->box);
2414             copy_mat(state_global->box,state->box);
2415
2416             if (vsite && (Flags & MD_RERUN_VSITE))
2417             {
2418                 if (DOMAINDECOMP(cr))
2419                 {
2420                     gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
2421                 }
2422                 if (graph)
2423                 {
2424                     /* Following is necessary because the graph may get out of sync
2425                      * with the coordinates if we only have every N'th coordinate set
2426                      */
2427                     mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
2428                     shift_self(graph,state->box,state->x);
2429                 }
2430                 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2431                                  top->idef.iparams,top->idef.il,
2432                                  fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2433                 if (graph)
2434                 {
2435                     unshift_self(graph,state->box,state->x);
2436                 }
2437             }
2438         }
2439
2440         /* Stop Center of Mass motion */
2441         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
2442
2443         /* Copy back starting coordinates in case we're doing a forcefield scan */
2444         if (bFFscan)
2445         {
2446             for(ii=0; (ii<state->natoms); ii++)
2447             {
2448                 copy_rvec(xcopy[ii],state->x[ii]);
2449                 copy_rvec(vcopy[ii],state->v[ii]);
2450             }
2451             copy_mat(boxcopy,state->box);
2452         }
2453
2454         if (bRerunMD)
2455         {
2456             /* for rerun MD always do Neighbour Searching */
2457             bNS = (bFirstStep || ir->nstlist != 0);
2458             bNStList = bNS;
2459         }
2460         else
2461         {
2462             /* Determine whether or not to do Neighbour Searching and LR */
2463             bNStList = (ir->nstlist > 0  && step % ir->nstlist == 0);
2464
2465             bNS = (bFirstStep || bExchanged || bNStList ||
2466                    (ir->nstlist == -1 && nlh.nabnsb > 0));
2467
2468             if (bNS && ir->nstlist == -1)
2469             {
2470                 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
2471             }
2472         }
2473
2474         /* < 0 means stop at next step, > 0 means stop at next NS step */
2475         if ( (gs.set[eglsSTOPCOND] < 0 ) ||
2476              ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
2477         {
2478             bLastStep = TRUE;
2479         }
2480
2481         /* Determine whether or not to update the Born radii if doing GB */
2482         bBornRadii=bFirstStep;
2483         if (ir->implicit_solvent && (step % ir->nstgbradii==0))
2484         {
2485             bBornRadii=TRUE;
2486         }
2487
2488         do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
2489         do_verbose = bVerbose &&
2490                   (step % stepout == 0 || bFirstStep || bLastStep);
2491
2492         if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
2493         {
2494             if (bRerunMD)
2495             {
2496                 bMasterState = TRUE;
2497             }
2498             else
2499             {
2500                 bMasterState = FALSE;
2501                 /* Correct the new box if it is too skewed */
2502                 if (DYNAMIC_BOX(*ir))
2503                 {
2504                     if (correct_box(fplog,step,state->box,graph))
2505                     {
2506                         bMasterState = TRUE;
2507                     }
2508                 }
2509                 if (DOMAINDECOMP(cr) && bMasterState)
2510                 {
2511                     dd_collect_state(cr->dd,state,state_global);
2512                 }
2513             }
2514
2515             if (DOMAINDECOMP(cr))
2516             {
2517                 /* Repartition the domain decomposition */
2518                 wallcycle_start(wcycle,ewcDOMDEC);
2519                 dd_partition_system(fplog,step,cr,
2520                                     bMasterState,nstglobalcomm,
2521                                     state_global,top_global,ir,
2522                                     state,&f,mdatoms,top,fr,
2523                                     vsite,shellfc,constr,
2524                                     nrnb,wcycle,do_verbose);
2525                 wallcycle_stop(wcycle,ewcDOMDEC);
2526                 /* If using an iterative integrator, reallocate space to match the decomposition */
2527             }
2528         }
2529
2530         if (MASTER(cr) && do_log && !bFFscan)
2531         {
2532             print_ebin_header(fplog,step,t,state->lambda);
2533         }
2534
2535         if (ir->efep != efepNO)
2536         {
2537             update_mdatoms(mdatoms,state->lambda);
2538         }
2539
2540         if (bRerunMD && rerun_fr.bV)
2541         {
2542
2543             /* We need the kinetic energy at minus the half step for determining
2544              * the full step kinetic energy and possibly for T-coupling.*/
2545             /* This may not be quite working correctly yet . . . . */
2546             compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2547                             wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
2548                             constr,NULL,FALSE,state->box,
2549                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2550                             CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
2551         }
2552         clear_mat(force_vir);
2553
2554         /* Ionize the atoms if necessary */
2555 /*        if (bIonize)
2556         {
2557             ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
2558                    mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
2559         }*/
2560
2561         /* Update force field in ffscan program */
2562 /*        if (bFFscan)
2563         {
2564             if (update_forcefield(fplog,
2565                                   nfile,fnm,fr,
2566                                   mdatoms->nr,state->x,state->box)) {
2567                 if (gmx_parallel_env_initialized())
2568                 {
2569                     gmx_finalize();
2570                 }
2571                 exit(0);
2572             }
2573         }*/
2574
2575         GMX_MPE_LOG(ev_timestep2);
2576
2577         /* We write a checkpoint at this MD step when:
2578          * either at an NS step when we signalled through gs,
2579          * or at the last step (but not when we do not want confout),
2580          * but never at the first step or with rerun.
2581          */
2582 /*        bCPT = (((gs.set[eglsCHKPT] && bNS) ||
2583                  (bLastStep && (Flags & MD_CONFOUT))) &&
2584                 step > ir->init_step && !bRerunMD);
2585         if (bCPT)
2586         {
2587             gs.set[eglsCHKPT] = 0;
2588         }*/
2589
2590         /* Determine the energy and pressure:
2591          * at nstcalcenergy steps and at energy output steps (set below).
2592          */
2593         bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2594         bCalcEnerPres = bNstEner;
2595
2596         /* Do we need global communication ? */
2597         bGStat = (bCalcEnerPres || bStopCM ||
2598                   (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2599
2600         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2601
2602         if (do_ene || do_log)
2603         {
2604             bCalcEnerPres = TRUE;
2605             bGStat    = TRUE;
2606         }
2607
2608         /* these CGLO_ options remain the same throughout the iteration */
2609         cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
2610                       (bStopCM ? CGLO_STOPCM : 0) |
2611                       (bGStat ? CGLO_GSTAT : 0)
2612             );
2613
2614         force_flags = (GMX_FORCE_STATECHANGED |
2615                        ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
2616                        GMX_FORCE_ALLFORCES |
2617                        (bNStList ? GMX_FORCE_DOLR : 0) |
2618                        GMX_FORCE_SEPLRF |
2619                        (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
2620                        (bDoDHDL ? GMX_FORCE_DHDL : 0)
2621             );
2622
2623         if (shellfc)
2624         {
2625             /* Now is the time to relax the shells */
2626             count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
2627                                       ir,bNS,force_flags,
2628                                       bStopCM,top,top_global,
2629                                       constr,enerd,fcd,
2630                                       state,f,force_vir,mdatoms,
2631                                       nrnb,wcycle,graph,groups,
2632                                       shellfc,fr,bBornRadii,t,mu_tot,
2633                                       state->natoms,&bConverged,vsite,
2634                                       outf->fp_field);
2635             tcount+=count;
2636
2637             if (bConverged)
2638             {
2639                 nconverged++;
2640             }
2641         }
2642         else
2643         {
2644             /* The coordinates (x) are shifted (to get whole molecules)
2645              * in do_force.
2646              * This is parallellized as well, and does communication too.
2647              * Check comments in sim_util.c
2648              */
2649
2650             do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
2651                      state->box,state->x,&state->hist,
2652                      f,force_vir,mdatoms,enerd,fcd,
2653                      state->lambda,graph,
2654                      fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
2655                      (bNS ? GMX_FORCE_NS : 0) | force_flags);
2656         }
2657
2658         GMX_BARRIER(cr->mpi_comm_mygroup);
2659
2660  /*       if (bTCR)
2661         {
2662             mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
2663                                    mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
2664         }
2665
2666         if (bTCR && bFirstStep)
2667         {
2668             tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
2669             fprintf(fplog,"Done init_coupling\n");
2670             fflush(fplog);
2671         }*/
2672
2673         /*  ############### START FIRST UPDATE HALF-STEP ############### */
2674
2675         if (bVV && !bStartingFromCpt && !bRerunMD)
2676         {
2677             if (ir->eI == eiVV)
2678             {
2679                 if (bInitStep)
2680                 {
2681                     /* if using velocity verlet with full time step Ekin,
2682                      * take the first half step only to compute the
2683                      * virial for the first step. From there,
2684                      * revert back to the initial coordinates
2685                      * so that the input is actually the initial step.
2686                      */
2687                     copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
2688                 }
2689
2690                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
2691                 if (!bInitStep)
2692                 {
2693                   trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
2694                 }
2695
2696                 if (ir->eI == eiVVAK)
2697                 {
2698                   update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2699                 }
2700
2701                 update_coords(fplog,step,ir,mdatoms,state,
2702                               f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2703                               ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
2704                               cr,nrnb,constr,&top->idef);
2705
2706                 if (bIterations)
2707                 {
2708                     gmx_iterate_init(&iterate,bIterations && !bInitStep);
2709                 }
2710                 /* for iterations, we save these vectors, as we will be self-consistently iterating
2711                    the calculations */
2712                 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
2713
2714                 /* save the state */
2715                 if (bIterations && iterate.bIterate) {
2716                     copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2717                 }
2718             }
2719
2720             bFirstIterate = TRUE;
2721             while (bFirstIterate || (bIterations && iterate.bIterate))
2722             {
2723                 if (bIterations && iterate.bIterate)
2724                 {
2725                     copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2726                     if (bFirstIterate && bTrotter)
2727                     {
2728                         /* The first time through, we need a decent first estimate
2729                            of veta(t+dt) to compute the constraints.  Do
2730                            this by computing the box volume part of the
2731                            trotter integration at this time. Nothing else
2732                            should be changed by this routine here.  If
2733                            !(first time), we start with the previous value
2734                            of veta.  */
2735
2736                         veta_save = state->veta;
2737                         trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
2738                         vetanew = state->veta;
2739                         state->veta = veta_save;
2740                     }
2741                 }
2742
2743                 bOK = TRUE;
2744                 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) {  /* Why is rerun_fr.bV here?  Unclear. */
2745                     dvdl = 0;
2746
2747                     update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2748                                        &top->idef,shake_vir,NULL,
2749                                        cr,nrnb,wcycle,upd,constr,
2750                                        bInitStep,TRUE,bCalcEnerPres,vetanew);
2751
2752                     if (!bOK && !bFFscan)
2753                     {
2754                         gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2755                     }
2756
2757                 }
2758                 else if (graph)
2759                 { /* Need to unshift here if a do_force has been
2760                      called in the previous step */
2761                     unshift_self(graph,state->box,state->x);
2762                 }
2763
2764
2765                 if (bVV) {
2766                     /* if VV, compute the pressure and constraints */
2767                     /* if VV2, the pressure and constraints only if using pressure control.*/
2768                     bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir));
2769                     bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));
2770                     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2771                                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2772                                     constr,NULL,FALSE,state->box,
2773                                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2774                                     cglo_flags
2775                                     | CGLO_ENERGY
2776                                     | (bTemp ? CGLO_TEMPERATURE:0)
2777                                     | (bPres ? CGLO_PRESSURE : 0)
2778                                     | (bPres ? CGLO_CONSTRAINT : 0)
2779                                     | (iterate.bIterate ? CGLO_ITERATE : 0)
2780                                     | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2781                                     | CGLO_SCALEEKIN
2782                         );
2783                 }
2784                 /* explanation of above:
2785                    a) We compute Ekin at the full time step
2786                    if 1) we are using the AveVel Ekin, and it's not the
2787                    initial step, or 2) if we are using AveEkin, but need the full
2788                    time step kinetic energy for the pressure.
2789                    b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2790                    EkinAveVel because it's needed for the pressure */
2791
2792                 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2793                 if (bVV && !bInitStep)
2794                 {
2795                   trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ2);
2796                 }
2797
2798                 if (bIterations &&
2799                     done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2800                                    state->veta,&vetanew))
2801                 {
2802                     break;
2803                 }
2804                 bFirstIterate = FALSE;
2805             }
2806
2807             if (bTrotter && !bInitStep) {
2808                 copy_mat(shake_vir,state->svir_prev);
2809                 copy_mat(force_vir,state->fvir_prev);
2810                 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
2811                     /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2812                     enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
2813                     enerd->term[F_EKIN] = trace(ekind->ekin);
2814                 }
2815             }
2816             /* if it's the initial step, we performed this first step just to get the constraint virial */
2817             if (bInitStep && ir->eI==eiVV) {
2818                 copy_rvecn(cbuf,state->v,0,state->natoms);
2819             }
2820
2821             if (fr->bSepDVDL && fplog && do_log)
2822             {
2823                 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2824             }
2825             enerd->term[F_DHDL_CON] += dvdl;
2826
2827             GMX_MPE_LOG(ev_timestep1);
2828
2829         }
2830
2831         /* MRS -- now done iterating -- compute the conserved quantity */
2832         if (bVV) {
2833             last_conserved = 0;
2834             if (IR_NVT_TROTTER(ir) || IR_NPT_TROTTER(ir))
2835             {
2836                 last_conserved =
2837                     NPT_energy(ir,state,&MassQ);
2838                 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
2839                 {
2840                     last_conserved -= enerd->term[F_DISPCORR];
2841                 }
2842             }
2843             if (ir->eI==eiVV) {
2844                 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
2845             }
2846         }
2847
2848         /* ########  END FIRST UPDATE STEP  ############## */
2849         /* ########  If doing VV, we now have v(dt) ###### */
2850
2851         /* ################## START TRAJECTORY OUTPUT ################# */
2852
2853         /* Now we have the energies and forces corresponding to the
2854          * coordinates at time t. We must output all of this before
2855          * the update.
2856          * for RerunMD t is read from input trajectory
2857          */
2858         GMX_MPE_LOG(ev_output_start);
2859
2860         mdof_flags = 0;
2861         if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
2862         if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
2863         if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
2864         if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
2865 /*        if (bCPT) { mdof_flags |= MDOF_CPT; };*/
2866
2867 #ifdef GMX_FAHCORE
2868         if (MASTER(cr))
2869             fcReportProgress( ir->nsteps, step );
2870
2871         if (bLastStep)
2872         {
2873             /* Enforce writing positions and velocities at end of run */
2874             mdof_flags |= (MDOF_X | MDOF_V);
2875         }
2876             /* sync bCPT and fc record-keeping */
2877 /*            if (bCPT && MASTER(cr))
2878                 fcRequestCheckPoint();*/
2879 #endif
2880
2881         if (mdof_flags != 0)
2882         {
2883             wallcycle_start(wcycle,ewcTRAJ);
2884 /*            if (bCPT)
2885             {
2886                 if (state->flags & (1<<estLD_RNG))
2887                 {
2888                     get_stochd_state(upd,state);
2889                 }
2890                 if (MASTER(cr))
2891                 {
2892                     if (bSumEkinhOld)
2893                     {
2894                         state_global->ekinstate.bUpToDate = FALSE;
2895                     }
2896                     else
2897                     {
2898                         update_ekinstate(&state_global->ekinstate,ekind);
2899                         state_global->ekinstate.bUpToDate = TRUE;
2900                     }
2901                     update_energyhistory(&state_global->enerhist,mdebin);
2902                 }
2903             }*/
2904             write_traj(fplog,cr,outf,mdof_flags,top_global,
2905                        step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2906 /*            if (bCPT)
2907             {
2908                 nchkpt++;
2909                 bCPT = FALSE;
2910             }*/
2911             debug_gmx();
2912             if (bLastStep && step_rel == ir->nsteps &&
2913                 (Flags & MD_CONFOUT) && MASTER(cr) &&
2914                 !bRerunMD && !bFFscan)
2915             {
2916                 /* x and v have been collected in write_traj,
2917                  * because a checkpoint file will always be written
2918                  * at the last step.
2919                  */
2920                 fprintf(stderr,"\nWriting final coordinates.\n");
2921                 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2922                     DOMAINDECOMP(cr))
2923                 {
2924                     /* Make molecules whole only for confout writing */
2925                     do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2926                 }
2927 /*                write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2928                                     *top_global->name,top_global,
2929                                     state_global->x,state_global->v,
2930                                     ir->ePBC,state->box);*/
2931                 debug_gmx();
2932             }
2933             wallcycle_stop(wcycle,ewcTRAJ);
2934         }
2935         GMX_MPE_LOG(ev_output_finish);
2936
2937         /* kludge -- virial is lost with restart for NPT control. Must restart */
2938         if (bStartingFromCpt && bVV)
2939         {
2940             copy_mat(state->svir_prev,shake_vir);
2941             copy_mat(state->fvir_prev,force_vir);
2942         }
2943         /*  ################## END TRAJECTORY OUTPUT ################ */
2944
2945         /* Determine the pressure:
2946          * always when we want exact averages in the energy file,
2947          * at ns steps when we have pressure coupling,
2948          * otherwise only at energy output steps (set below).
2949          */
2950
2951         bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2952         bCalcEnerPres = bNstEner;
2953
2954         /* Do we need global communication ? */
2955         bGStat = (bGStatEveryStep || bStopCM || bNS ||
2956                   (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2957
2958         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2959
2960         if (do_ene || do_log)
2961         {
2962             bCalcEnerPres = TRUE;
2963             bGStat        = TRUE;
2964         }
2965
2966         /* Determine the wallclock run time up till now */
2967         run_time = gmx_gettime() - (double)runtime->real;
2968
2969         /* Check whether everything is still allright */
2970         if (((int)gmx_get_stop_condition() > handled_stop_condition)
2971 #ifdef GMX_THREADS
2972             && MASTER(cr)
2973 #endif
2974             )
2975         {
2976             /* this is just make gs.sig compatible with the hack
2977                of sending signals around by MPI_Reduce with together with
2978                other floats */
2979             if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2980                 gs.sig[eglsSTOPCOND]=1;
2981             if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2982                 gs.sig[eglsSTOPCOND]=-1;
2983             /* < 0 means stop at next step, > 0 means stop at next NS step */
2984             if (fplog)
2985             {
2986                 fprintf(fplog,
2987                         "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2988                         gmx_get_signal_name(),
2989                         gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2990                 fflush(fplog);
2991             }
2992             fprintf(stderr,
2993                     "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2994                     gmx_get_signal_name(),
2995                     gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2996             fflush(stderr);
2997             handled_stop_condition=(int)gmx_get_stop_condition();
2998         }
2999         else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
3000                  (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
3001                  gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
3002         {
3003             /* Signal to terminate the run */
3004             gs.sig[eglsSTOPCOND] = 1;
3005             if (fplog)
3006             {
3007                 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
3008             }
3009             fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
3010         }
3011
3012         if (bResetCountersHalfMaxH && MASTER(cr) &&
3013             run_time > max_hours*60.0*60.0*0.495)
3014         {
3015             gs.sig[eglsRESETCOUNTERS] = 1;
3016         }
3017
3018         if (ir->nstlist == -1 && !bRerunMD)
3019         {
3020             /* When bGStatEveryStep=FALSE, global_stat is only called
3021              * when we check the atom displacements, not at NS steps.
3022              * This means that also the bonded interaction count check is not
3023              * performed immediately after NS. Therefore a few MD steps could
3024              * be performed with missing interactions.
3025              * But wrong energies are never written to file,
3026              * since energies are only written after global_stat
3027              * has been called.
3028              */
3029             if (step >= nlh.step_nscheck)
3030             {
3031                 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
3032                                                      nlh.scale_tot,state->x);
3033             }
3034             else
3035             {
3036                 /* This is not necessarily true,
3037                  * but step_nscheck is determined quite conservatively.
3038                  */
3039                 nlh.nabnsb = 0;
3040             }
3041         }
3042
3043         /* In parallel we only have to check for checkpointing in steps
3044          * where we do global communication,
3045          *  otherwise the other nodes don't know.
3046          */
3047         if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
3048                            cpt_period >= 0 &&
3049                            (cpt_period == 0 ||
3050                             run_time >= nchkpt*cpt_period*60.0)) &&
3051             gs.set[eglsCHKPT] == 0)
3052         {
3053             gs.sig[eglsCHKPT] = 1;
3054         }
3055
3056         if (bIterations)
3057         {
3058             gmx_iterate_init(&iterate,bIterations);
3059         }
3060
3061         /* for iterations, we save these vectors, as we will be redoing the calculations */
3062         if (bIterations && iterate.bIterate)
3063         {
3064             copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
3065         }
3066         bFirstIterate = TRUE;
3067         while (bFirstIterate || (bIterations && iterate.bIterate))
3068         {
3069             /* We now restore these vectors to redo the calculation with improved extended variables */
3070             if (bIterations)
3071             {
3072                 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
3073             }
3074
3075             /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
3076                so scroll down for that logic */
3077
3078             /* #########   START SECOND UPDATE STEP ################# */
3079             GMX_MPE_LOG(ev_update_start);
3080             bOK = TRUE;
3081             if (!bRerunMD || rerun_fr.bV || bForceUpdate)
3082             {
3083                 wallcycle_start(wcycle,ewcUPDATE);
3084                 dvdl = 0;
3085                 /* Box is changed in update() when we do pressure coupling,
3086                  * but we should still use the old box for energy corrections and when
3087                  * writing it to the energy file, so it matches the trajectory files for
3088                  * the same timestep above. Make a copy in a separate array.
3089                  */
3090                 copy_mat(state->box,lastbox);
3091                 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
3092                 if (bTrotter)
3093                 {
3094                     if (bIterations && iterate.bIterate)
3095                     {
3096                         if (bFirstIterate)
3097                         {
3098                             scalevir = 1;
3099                         }
3100                         else
3101                         {
3102                             /* we use a new value of scalevir to converge the iterations faster */
3103                             scalevir = tracevir/trace(shake_vir);
3104                         }
3105                         msmul(shake_vir,scalevir,shake_vir);
3106                         m_add(force_vir,shake_vir,total_vir);
3107                         clear_mat(shake_vir);
3108                     }
3109                     trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ3);
3110                 }
3111                 /* We can only do Berendsen coupling after we have summed
3112                  * the kinetic energy or virial. Since the happens
3113                  * in global_state after update, we should only do it at
3114                  * step % nstlist = 1 with bGStatEveryStep=FALSE.
3115                  */
3116
3117                 if (ir->eI != eiVVAK)
3118                 {
3119                   update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
3120                 }
3121                 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
3122                                 upd,bInitStep);
3123
3124                 if (bVV)
3125                 {
3126                     /* velocity half-step update */
3127                     update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3128                                   ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,cr,nrnb,constr,&top->idef);
3129                 }
3130
3131                 /* Above, initialize just copies ekinh into ekin,
3132                  * it doesn't copy position (for VV),
3133                  * and entire integrator for MD.
3134                  */
3135
3136                 if (ir->eI==eiVVAK)
3137                 {
3138                     copy_rvecn(state->x,cbuf,0,state->natoms);
3139                 }
3140
3141                 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3142                               ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
3143                 wallcycle_stop(wcycle,ewcUPDATE);
3144
3145                 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
3146                                    &top->idef,shake_vir,force_vir,
3147                                    cr,nrnb,wcycle,upd,constr,
3148                                    bInitStep,FALSE,bCalcEnerPres,state->veta);
3149
3150                 if (ir->eI==eiVVAK)
3151                 {
3152                     /* erase F_EKIN and F_TEMP here? */
3153                     /* just compute the kinetic energy at the half step to perform a trotter step */
3154                     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
3155                                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
3156                                     constr,NULL,FALSE,lastbox,
3157                                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
3158                                     cglo_flags | CGLO_TEMPERATURE | CGLO_CONSTRAINT
3159                         );
3160                     wallcycle_start(wcycle,ewcUPDATE);
3161                     trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ4);
3162                     /* now we know the scaling, we can compute the positions again again */
3163                     copy_rvecn(cbuf,state->x,0,state->natoms);
3164
3165                     update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3166                                   ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
3167                     wallcycle_stop(wcycle,ewcUPDATE);
3168
3169                     /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
3170                     /* are the small terms in the shake_vir here due
3171                      * to numerical errors, or are they important
3172                      * physically? I'm thinking they are just errors, but not completely sure.
3173                      * For now, will call without actually constraining, constr=NULL*/
3174                     update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
3175                                        &top->idef,tmp_vir,force_vir,
3176                                        cr,nrnb,wcycle,upd,NULL,
3177                                        bInitStep,FALSE,bCalcEnerPres,state->veta);
3178                 }
3179                 if (!bOK && !bFFscan)
3180                 {
3181                     gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
3182                 }
3183
3184                 if (fr->bSepDVDL && fplog && do_log)
3185                 {
3186                     fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
3187                 }
3188                 enerd->term[F_DHDL_CON] += dvdl;
3189             }
3190             else if (graph)
3191             {
3192                 /* Need to unshift here */
3193                 unshift_self(graph,state->box,state->x);
3194             }
3195
3196             GMX_BARRIER(cr->mpi_comm_mygroup);
3197             GMX_MPE_LOG(ev_update_finish);
3198
3199             if (vsite != NULL)
3200             {
3201                 wallcycle_start(wcycle,ewcVSITECONSTR);
3202                 if (graph != NULL)
3203                 {
3204                     shift_self(graph,state->box,state->x);
3205                 }
3206                 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
3207                                  top->idef.iparams,top->idef.il,
3208                                  fr->ePBC,fr->bMolPBC,graph,cr,state->box);
3209
3210                 if (graph != NULL)
3211                 {
3212                     unshift_self(graph,state->box,state->x);
3213                 }
3214                 wallcycle_stop(wcycle,ewcVSITECONSTR);
3215             }
3216
3217             /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
3218             if (ir->nstlist == -1 && bFirstIterate)
3219             {
3220                 gs.sig[eglsNABNSB] = nlh.nabnsb;
3221             }
3222             compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
3223                             wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
3224                             constr,
3225                             bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
3226                             lastbox,
3227                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
3228                             cglo_flags
3229                             | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
3230                             | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
3231                             | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
3232                             | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
3233                             | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
3234                             | CGLO_CONSTRAINT
3235                 );
3236             if (ir->nstlist == -1 && bFirstIterate)
3237             {
3238                 nlh.nabnsb = gs.set[eglsNABNSB];
3239                 gs.set[eglsNABNSB] = 0;
3240             }
3241             /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
3242             /* #############  END CALC EKIN AND PRESSURE ################# */
3243
3244             /* Note: this is OK, but there are some numerical precision issues with using the convergence of
3245                the virial that should probably be addressed eventually. state->veta has better properies,
3246                but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
3247                generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
3248
3249             if (bIterations &&
3250                 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
3251                                trace(shake_vir),&tracevir))
3252             {
3253                 break;
3254             }
3255             bFirstIterate = FALSE;
3256         }
3257
3258         update_box(fplog,step,ir,mdatoms,state,graph,f,
3259                    ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
3260
3261         /* ################# END UPDATE STEP 2 ################# */
3262         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
3263
3264         /* The coordinates (x) were unshifted in update */
3265 /*        if (bFFscan && (shellfc==NULL || bConverged))
3266         {
3267             if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
3268                                  f,NULL,xcopy,
3269                                  &(top_global->mols),mdatoms->massT,pres))
3270             {
3271                 if (gmx_parallel_env_initialized())
3272                 {
3273                     gmx_finalize();
3274                 }
3275                 fprintf(stderr,"\n");
3276                 exit(0);
3277             }
3278         }*/
3279         if (!bGStat)
3280         {
3281             /* We will not sum ekinh_old,
3282              * so signal that we still have to do it.
3283              */
3284             bSumEkinhOld = TRUE;
3285         }
3286
3287 /*        if (bTCR)
3288         {*/
3289             /* Only do GCT when the relaxation of shells (minimization) has converged,
3290              * otherwise we might be coupling to bogus energies.
3291              * In parallel we must always do this, because the other sims might
3292              * update the FF.
3293              */
3294
3295             /* Since this is called with the new coordinates state->x, I assume
3296              * we want the new box state->box too. / EL 20040121
3297              */
3298 /*            do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
3299                         ir,MASTER(cr),
3300                         mdatoms,&(top->idef),mu_aver,
3301                         top_global->mols.nr,cr,
3302                         state->box,total_vir,pres,
3303                         mu_tot,state->x,f,bConverged);
3304             debug_gmx();
3305         }*/
3306
3307         /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
3308
3309         sum_dhdl(enerd,state->lambda,ir);
3310         /* use the directly determined last velocity, not actually the averaged half steps */
3311         if (bTrotter && ir->eI==eiVV)
3312         {
3313             enerd->term[F_EKIN] = last_ekin;
3314         }
3315         enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
3316
3317         switch (ir->etc)
3318         {
3319         case etcNO:
3320             break;
3321         case etcBERENDSEN:
3322             break;
3323         case etcNOSEHOOVER:
3324             if (IR_NVT_TROTTER(ir)) {
3325                 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + last_conserved;
3326             } else {
3327                 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] +
3328                     NPT_energy(ir,state,&MassQ);
3329             }
3330             break;
3331         case etcVRESCALE:
3332             enerd->term[F_ECONSERVED] =
3333                 enerd->term[F_ETOT] + vrescale_energy(&(ir->opts),
3334                                                       state->therm_integral);
3335             break;
3336         default:
3337             break;
3338         }
3339
3340         /* Check for excessively large energies */
3341 /*        if (bIonize)
3342         {
3343 #ifdef GMX_DOUBLE
3344             real etot_max = 1e200;
3345 #else
3346             real etot_max = 1e30;
3347 #endif
3348             if (fabs(enerd->term[F_ETOT]) > etot_max)
3349             {
3350                 fprintf(stderr,"Energy too large (%g), giving up\n",
3351                         enerd->term[F_ETOT]);
3352             }
3353         }*/
3354         /* #########  END PREPARING EDR OUTPUT  ###########  */
3355
3356         /* Time for performance */
3357         if (((step % stepout) == 0) || bLastStep)
3358         {
3359             runtime_upd_proc(runtime);
3360         }
3361
3362         /* Output stuff */
3363         if (MASTER(cr))
3364         {
3365             gmx_bool do_dr,do_or;
3366
3367             if (!(bStartingFromCpt && (EI_VV(ir->eI))))
3368             {
3369                 if (bNstEner)
3370                 {
3371                     upd_mdebin(mdebin,bDoDHDL,TRUE,
3372                                t,mdatoms->tmass,enerd,state,lastbox,
3373                                shake_vir,force_vir,total_vir,pres,
3374                                ekind,mu_tot,constr);
3375                 }
3376                 else
3377                 {
3378                     upd_mdebin_step(mdebin);
3379                 }
3380
3381                 do_dr  = do_per_step(step,ir->nstdisreout);
3382                 do_or  = do_per_step(step,ir->nstorireout);
3383
3384                 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
3385                            step,t,
3386                            eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
3387             }
3388             if (ir->ePull != epullNO)
3389             {
3390                 pull_print_output(ir->pull,step,t);
3391             }
3392
3393             if (do_per_step(step,ir->nstlog))
3394             {
3395                 if(fflush(fplog) != 0)
3396                 {
3397                     gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
3398                 }
3399             }
3400         }
3401
3402
3403         /* Remaining runtime */
3404         if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
3405         {
3406             if (shellfc)
3407             {
3408                 fprintf(stderr,"\n");
3409             }
3410             print_time(stderr,runtime,step,ir,cr);
3411         }
3412
3413                 /* Set new positions for the group to embed */
3414                 if(!bLastStep){
3415                         if(step_rel<=it_xy)
3416                         {
3417                                 fac[0]+=xy_step;
3418                                 fac[1]+=xy_step;
3419                         } else if (step_rel<=(it_xy+it_z))
3420                         {
3421                                 fac[2]+=z_step;
3422                         }
3423                         resize(ins_at,r_ins,state_global->x,pos_ins,fac);
3424                 }
3425
3426         /* Replica exchange */
3427 /*        bExchanged = FALSE;
3428         if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
3429             do_per_step(step,repl_ex_nst))
3430         {
3431             bExchanged = replica_exchange(fplog,cr,repl_ex,
3432                                           state_global,enerd->term,
3433                                           state,step,t);
3434         }
3435         if (bExchanged && PAR(cr))
3436         {
3437             if (DOMAINDECOMP(cr))
3438             {
3439                 dd_partition_system(fplog,step,cr,TRUE,1,
3440                                     state_global,top_global,ir,
3441                                     state,&f,mdatoms,top,fr,
3442                                     vsite,shellfc,constr,
3443                                     nrnb,wcycle,FALSE);
3444             }
3445             else
3446             {
3447                 bcast_state(cr,state,FALSE);
3448             }
3449         }*/
3450
3451         bFirstStep = FALSE;
3452         bInitStep = FALSE;
3453         bStartingFromCpt = FALSE;
3454
3455         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
3456         /* Complicated conditional when bGStatEveryStep=FALSE.
3457          * We can not just use bGStat, since then the simulation results
3458          * would depend on nstenergy and nstlog or step_nscheck.
3459          */
3460         if (((state->flags & (1<<estPRES_PREV)) ||
3461              (state->flags & (1<<estSVIR_PREV)) ||
3462              (state->flags & (1<<estFVIR_PREV))) &&
3463             (bGStatEveryStep ||
3464              (ir->nstlist > 0 && step % ir->nstlist == 0) ||
3465              (ir->nstlist < 0 && nlh.nabnsb > 0) ||
3466              (ir->nstlist == 0 && bGStat)))
3467         {
3468             /* Store the pressure in t_state for pressure coupling
3469              * at the next MD step.
3470              */
3471             if (state->flags & (1<<estPRES_PREV))
3472             {
3473                 copy_mat(pres,state->pres_prev);
3474             }
3475         }
3476
3477         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
3478
3479         if (bRerunMD)
3480         {
3481             /* read next frame from input trajectory */
3482             bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
3483         }
3484
3485         if (!bRerunMD || !rerun_fr.bStep)
3486         {
3487             /* increase the MD step number */
3488             step++;
3489             step_rel++;
3490         }
3491
3492         cycles = wallcycle_stop(wcycle,ewcSTEP);
3493         if (DOMAINDECOMP(cr) && wcycle)
3494         {
3495             dd_cycles_add(cr->dd,cycles,ddCyclStep);
3496         }
3497
3498         if (step_rel == wcycle_get_reset_counters(wcycle) ||
3499             gs.set[eglsRESETCOUNTERS] != 0)
3500         {
3501             /* Reset all the counters related to performance over the run */
3502             reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
3503             wcycle_set_reset_counters(wcycle,-1);
3504             bResetCountersHalfMaxH = FALSE;
3505             gs.set[eglsRESETCOUNTERS] = 0;
3506         }
3507     }
3508     /* End of main MD loop */
3509     debug_gmx();
3510     write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
3511                                         *top_global->name,top_global,
3512                                         state_global->x,state_global->v,
3513                                         ir->ePBC,state->box);
3514
3515     /* Stop the time */
3516     runtime_end(runtime);
3517
3518     if (bRerunMD)
3519     {
3520         close_trj(status);
3521     }
3522
3523     if (!(cr->duty & DUTY_PME))
3524     {
3525         /* Tell the PME only node to finish */
3526         gmx_pme_finish(cr);
3527     }
3528
3529     if (MASTER(cr))
3530     {
3531         if (ir->nstcalcenergy > 0 && !bRerunMD)
3532         {
3533             print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
3534                        eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
3535         }
3536     }
3537
3538     done_mdoutf(outf);
3539
3540     debug_gmx();
3541
3542     if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
3543     {
3544         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)));
3545         fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
3546     }
3547
3548     if (shellfc && fplog)
3549     {
3550         fprintf(fplog,"Fraction of iterations that converged:           %.2f %%\n",
3551                 (nconverged*100.0)/step_rel);
3552         fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
3553                 tcount/step_rel);
3554     }
3555
3556 /*    if (repl_ex_nst > 0 && MASTER(cr))
3557     {
3558         print_replica_exchange_statistics(fplog,repl_ex);
3559     }*/
3560
3561     runtime->nsteps_done = step_rel;
3562
3563     return 0;
3564 }
3565
3566
3567 int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
3568              const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
3569              int nstglobalcomm,
3570              ivec ddxyz,int dd_node_order,real rdd,real rconstr,
3571              const char *dddlb_opt,real dlb_scale,
3572              const char *ddcsx,const char *ddcsy,const char *ddcsz,
3573              int nstepout,int resetstep,int nmultisim,int repl_ex_nst,int repl_ex_seed,
3574              real pforce,real cpt_period,real max_hours,
3575              const char *deviceOptions,
3576              unsigned long Flags,
3577              real xy_fac, real xy_max, real z_fac, real z_max,
3578              int it_xy, int it_z, real probe_rad, int low_up_rm,
3579              int pieces, gmx_bool bALLOW_ASYMMETRY, int maxwarn)
3580 {
3581     double     nodetime=0,realtime;
3582     t_inputrec *inputrec;
3583     t_state    *state=NULL;
3584     matrix     box;
3585     gmx_ddbox_t ddbox;
3586     int        npme_major,npme_minor;
3587     real       tmpr1,tmpr2;
3588     t_nrnb     *nrnb;
3589     gmx_mtop_t *mtop=NULL;
3590     t_mdatoms  *mdatoms=NULL;
3591     t_forcerec *fr=NULL;
3592     t_fcdata   *fcd=NULL;
3593     real       ewaldcoeff=0;
3594     gmx_pme_t  *pmedata=NULL;
3595     gmx_vsite_t *vsite=NULL;
3596     gmx_constr_t constr;
3597     int        i,m,nChargePerturbed=-1,status,nalloc;
3598     char       *gro;
3599     gmx_wallcycle_t wcycle;
3600     gmx_bool       bReadRNG,bReadEkin;
3601     int        list;
3602     gmx_runtime_t runtime;
3603     int        rc;
3604     gmx_large_int_t reset_counters;
3605     gmx_edsam_t ed=NULL;
3606     t_commrec   *cr_old=cr;
3607     int        nthreads=1,nthreads_requested=1;
3608
3609
3610         char                    *ins;
3611         int                     rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
3612         int                     ng,j,max_lip_rm,ins_grp_id,ins_nat,mem_nat,ntype,lip_rm,tpr_version;
3613         real                    xy_step=0,z_step=0;
3614         real                    prot_area;
3615         rvec                    *r_ins=NULL,fac;
3616         t_block                 *ins_at,*rest_at;
3617         pos_ins_t               *pos_ins;
3618         mem_t                   *mem_p;
3619         rm_t                    *rm_p;
3620         gmx_groups_t            *groups;
3621         gmx_bool                        bExcl=FALSE;
3622         t_atoms                 atoms;
3623         t_pbc                   *pbc;
3624         char                    **piecename=NULL;
3625
3626     /* CAUTION: threads may be started later on in this function, so
3627        cr doesn't reflect the final parallel state right now */
3628     snew(inputrec,1);
3629     snew(mtop,1);
3630
3631     if (bVerbose && SIMMASTER(cr))
3632     {
3633         fprintf(stderr,"Getting Loaded...\n");
3634     }
3635
3636     if (Flags & MD_APPENDFILES)
3637     {
3638         fplog = NULL;
3639     }
3640
3641     snew(state,1);
3642     if (MASTER(cr))
3643     {
3644         /* Read (nearly) all data required for the simulation */
3645         read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
3646
3647         /* NOW the threads will be started: */
3648 #ifdef GMX_THREADS
3649 #endif
3650     }
3651     /* END OF CAUTION: cr is now reliable */
3652
3653     if (PAR(cr))
3654     {
3655         /* now broadcast everything to the non-master nodes/threads: */
3656         init_parallel(fplog, cr, inputrec, mtop);
3657     }
3658     /* now make sure the state is initialized and propagated */
3659     set_state_entries(state,inputrec,cr->nnodes);
3660
3661     if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
3662     {
3663         /* All-vs-all loops do not work with domain decomposition */
3664         Flags |= MD_PARTDEC;
3665     }
3666
3667     if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
3668     {
3669         cr->npmenodes = 0;
3670     }
3671
3672         snew(ins_at,1);
3673         snew(pos_ins,1);
3674         if(MASTER(cr))
3675         {
3676                 tpr_version = get_tpr_version(ftp2fn(efTPX,nfile,fnm));
3677                 if (tpr_version<58)
3678                         gmx_fatal(FARGS,"Version of *.tpr file to old (%d). Rerun grompp with gromacs VERSION 4.0.3 or newer.\n",tpr_version);
3679
3680                 if( inputrec->eI != eiMD )
3681                         gmx_input("Change integrator to md in mdp file.");
3682
3683                 if(PAR(cr))
3684                         gmx_input("Sorry, parallel g_membed is not yet fully functrional.");
3685
3686                 groups=&(mtop->groups);
3687
3688                 atoms=gmx_mtop_global_atoms(mtop);
3689                 snew(mem_p,1);
3690                 fprintf(stderr,"\nSelect a group to embed in the membrane:\n");
3691                 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(ins_at->nr),&(ins_at->index),&ins);
3692                 ins_grp_id = search_string(ins,groups->ngrpname,(groups->grpname));
3693                 fprintf(stderr,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins);
3694                 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(mem_p->mem_at.nr),&(mem_p->mem_at.index),&(mem_p->name));
3695
3696                 pos_ins->pieces=pieces;
3697                 snew(pos_ins->nidx,pieces);
3698                 snew(pos_ins->subindex,pieces);
3699                 snew(piecename,pieces);
3700                 if (pieces>1)
3701                 {
3702                         fprintf(stderr,"\nSelect pieces to embed:\n");
3703                         get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),pieces,pos_ins->nidx,pos_ins->subindex,piecename);
3704                 }
3705                 else
3706                 {
3707                         /*use whole embedded group*/
3708                         snew(pos_ins->nidx,1);
3709                         snew(pos_ins->subindex,1);
3710                         pos_ins->nidx[0]=ins_at->nr;
3711                         pos_ins->subindex[0]=ins_at->index;
3712                 }
3713
3714                 if(probe_rad<0.2199999)
3715                 {
3716                         warn++;
3717                         fprintf(stderr,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 can result in overlap between waters "
3718                                         "and the group to embed, which will result in Lincs errors etc.\nIf you are sure, you can increase maxwarn.\n\n",warn);
3719                 }
3720
3721                 if(xy_fac<0.09999999)
3722                 {
3723                         warn++;
3724                         fprintf(stderr,"\nWarning %d:\nThe initial size of %s is probably too smal.\n"
3725                                         "If you are sure, you can increase maxwarn.\n\n",warn,ins);
3726                 }
3727
3728                 if(it_xy<1000)
3729                 {
3730                         warn++;
3731                         fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d) is probably too small.\n"
3732                                         "Increase -nxy or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_xy);
3733                 }
3734
3735                 if( (it_z<100) && ( z_fac<0.99999999 || z_fac>1.0000001) )
3736                 {
3737                         warn++;
3738                         fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d) is probably too small.\n"
3739                                        "Increase -nz or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_z);
3740                 }
3741
3742                 if(it_xy+it_z>inputrec->nsteps)
3743                 {
3744                         warn++;
3745                         fprintf(stderr,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the number of steps in the tpr.\n"
3746                                         "If you are sure, you can increase maxwarn.\n\n",warn);
3747                 }
3748
3749                 fr_id=-1;
3750                 if( inputrec->opts.ngfrz==1)
3751                         gmx_fatal(FARGS,"You did not specify \"%s\" as a freezegroup.",ins);
3752                 for(i=0;i<inputrec->opts.ngfrz;i++)
3753                 {
3754                         tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
3755                         if(ins_grp_id==tmp_id)
3756                         {
3757                                 fr_id=tmp_id;
3758                                 fr_i=i;
3759                         }
3760                 }
3761                 if (fr_id == -1 )
3762                         gmx_fatal(FARGS,"\"%s\" not as freezegroup defined in the mdp-file.",ins);
3763
3764                 for(i=0;i<DIM;i++)
3765                         if( inputrec->opts.nFreeze[fr_i][i] != 1)
3766                                 gmx_fatal(FARGS,"freeze dimensions for %s are not Y Y Y\n",ins);
3767
3768                 ng = groups->grps[egcENER].nr;
3769                 if (ng == 1)
3770                         gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
3771
3772                 for(i=0;i<ng;i++)
3773                 {
3774                         for(j=0;j<ng;j++)
3775                         {
3776                                 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
3777                                 {
3778                                         bExcl = TRUE;
3779                                         if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
3780                                                 gmx_fatal(FARGS,"Energy exclusions \"%s\" and  \"%s\" do not match the group to embed \"%s\"",
3781                                                                 *groups->grpname[groups->grps[egcENER].nm_ind[i]],
3782                                                                 *groups->grpname[groups->grps[egcENER].nm_ind[j]],ins);
3783                                 }
3784                         }
3785                 }
3786                 if (!bExcl)
3787                         gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in the freeze group");
3788
3789                 /* Set all atoms in box*/
3790                 /*set_inbox(state->natoms,state->x);*/
3791
3792                 /* Guess the area the protein will occupy in the membrane plane  Calculate area per lipid*/
3793                 snew(rest_at,1);
3794                 ins_nat = init_ins_at(ins_at,rest_at,state,pos_ins,groups,ins_grp_id,xy_max);
3795                 /* Check moleculetypes in insertion group */
3796                 check_types(ins_at,rest_at,mtop);
3797
3798                 mem_nat = init_mem_at(mem_p,mtop,state->x,state->box,pos_ins);
3799
3800                 prot_area = est_prot_area(pos_ins,state->x,ins_at,mem_p);
3801                 if ( (prot_area>7.5) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX])<50) )
3802                 {
3803                         warn++;
3804                         fprintf(stderr,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
3805                                         "This might cause pressure problems during the growth phase. Just try with\n"
3806                                         "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
3807                                         "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn);
3808                 }
3809                 if(warn>maxwarn)
3810                                         gmx_fatal(FARGS,"Too many warnings.\n");
3811
3812                 printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area);
3813                 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);
3814
3815                 /* Maximum number of lipids to be removed*/
3816                 max_lip_rm=(int)(2*prot_area/mem_p->lip_area);
3817                 printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm);
3818
3819                 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
3820                                 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
3821                                 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",xy_fac,z_fac,mem_p->zmin,mem_p->zmax);
3822
3823                 /* resize the protein by xy and by z if necessary*/
3824                 snew(r_ins,ins_at->nr);
3825                 init_resize(ins_at,r_ins,pos_ins,mem_p,state->x,bALLOW_ASYMMETRY);
3826                 fac[0]=fac[1]=xy_fac;
3827                 fac[2]=z_fac;
3828
3829                 xy_step =(xy_max-xy_fac)/(double)(it_xy);
3830                 z_step  =(z_max-z_fac)/(double)(it_z-1);
3831
3832                 resize(ins_at,r_ins,state->x,pos_ins,fac);
3833
3834                 /* remove overlapping lipids and water from the membrane box*/
3835                 /*mark molecules to be removed*/
3836                 snew(pbc,1);
3837                 set_pbc(pbc,inputrec->ePBC,state->box);
3838
3839                 snew(rm_p,1);
3840                 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);
3841         lip_rm -= low_up_rm;
3842
3843                 if(fplog)
3844                         for(i=0;i<rm_p->nr;i++)
3845                                 fprintf(fplog,"rm mol %d\n",rm_p->mol[i]);
3846
3847                 for(i=0;i<mtop->nmolblock;i++)
3848                 {
3849                         ntype=0;
3850                         for(j=0;j<rm_p->nr;j++)
3851                                 if(rm_p->block[j]==i)
3852                                         ntype++;
3853                         printf("Will remove %d %s molecules\n",ntype,*(mtop->moltype[mtop->molblock[i].type].name));
3854                 }
3855
3856                 if(lip_rm>max_lip_rm)
3857                 {
3858                         warn++;
3859                         fprintf(stderr,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated protein area\n"
3860                                         "Try making the -xyinit resize factor smaller. If you are sure about this increase maxwarn.\n\n",warn);
3861                 }
3862
3863                 /*remove all lipids and waters overlapping and update all important structures*/
3864                 rm_group(inputrec,groups,mtop,rm_p,state,ins_at,pos_ins);
3865
3866                 rm_bonded_at = rm_bonded(ins_at,mtop);
3867                 if (rm_bonded_at != ins_at->nr)
3868                 {
3869                         fprintf(stderr,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
3870                                         "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
3871                                         "molecule type as atoms that are not to be embedded.\n",rm_bonded_at,ins_at->nr);
3872                 }
3873
3874                 if(warn>maxwarn)
3875                         gmx_fatal(FARGS,"Too many warnings.\nIf you are sure these warnings are harmless, you can increase -maxwarn");
3876
3877                 if (MASTER(cr))
3878                 {
3879                         if (ftp2bSet(efTOP,nfile,fnm))
3880                                 top_update(opt2fn("-p",nfile,fnm),ins,rm_p,mtop);
3881                 }
3882
3883                 sfree(pbc);
3884                 sfree(rest_at);
3885         }
3886
3887 #ifdef GMX_FAHCORE
3888     fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
3889 #endif
3890
3891     /* NMR restraints must be initialized before load_checkpoint,
3892      * since with time averaging the history is added to t_state.
3893      * For proper consistency check we therefore need to extend
3894      * t_state here.
3895      * So the PME-only nodes (if present) will also initialize
3896      * the distance restraints.
3897      */
3898     snew(fcd,1);
3899
3900     /* This needs to be called before read_checkpoint to extend the state */
3901     init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
3902
3903     if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
3904     {
3905         if (PAR(cr) && !(Flags & MD_PARTDEC))
3906         {
3907             gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
3908         }
3909         /* Orientation restraints */
3910         if (MASTER(cr))
3911         {
3912             init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
3913                         state);
3914         }
3915     }
3916
3917     if (DEFORM(*inputrec))
3918     {
3919         /* Store the deform reference box before reading the checkpoint */
3920         if (SIMMASTER(cr))
3921         {
3922             copy_mat(state->box,box);
3923         }
3924         if (PAR(cr))
3925         {
3926             gmx_bcast(sizeof(box),box,cr);
3927         }
3928         /* Because we do not have the update struct available yet
3929          * in which the reference values should be stored,
3930          * we store them temporarily in static variables.
3931          * This should be thread safe, since they are only written once
3932          * and with identical values.
3933          */
3934 /*        deform_init_init_step_tpx = inputrec->init_step;*/
3935 /*        copy_mat(box,deform_init_box_tpx);*/
3936     }
3937
3938     if (opt2bSet("-cpi",nfile,fnm))
3939     {
3940         /* Check if checkpoint file exists before doing continuation.
3941          * This way we can use identical input options for the first and subsequent runs...
3942          */
3943         if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
3944         {
3945             load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
3946                             cr,Flags & MD_PARTDEC,ddxyz,
3947                             inputrec,state,&bReadRNG,&bReadEkin,
3948                             (Flags & MD_APPENDFILES));
3949
3950             if (bReadRNG)
3951             {
3952                 Flags |= MD_READ_RNG;
3953             }
3954             if (bReadEkin)
3955             {
3956                 Flags |= MD_READ_EKIN;
3957             }
3958         }
3959     }
3960
3961     if ((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
3962     {
3963         gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
3964                              Flags,&fplog);
3965     }
3966
3967     if (SIMMASTER(cr))
3968     {
3969         copy_mat(state->box,box);
3970     }
3971
3972     if (PAR(cr))
3973     {
3974         gmx_bcast(sizeof(box),box,cr);
3975     }
3976
3977     if (bVerbose && SIMMASTER(cr))
3978     {
3979         fprintf(stderr,"Loaded with Money\n\n");
3980     }
3981
3982     if (PAR(cr) && !((Flags & MD_PARTDEC) || EI_TPI(inputrec->eI)))
3983     {
3984         cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
3985                                            dddlb_opt,dlb_scale,
3986                                            ddcsx,ddcsy,ddcsz,
3987                                            mtop,inputrec,
3988                                            box,state->x,
3989                                            &ddbox,&npme_major,&npme_minor);
3990
3991         make_dd_communicators(fplog,cr,dd_node_order);
3992
3993         /* Set overallocation to avoid frequent reallocation of arrays */
3994         set_over_alloc_dd(TRUE);
3995     }
3996     else
3997     {
3998         /* PME, if used, is done on all nodes with 1D decomposition */
3999         cr->npmenodes = 0;
4000         cr->duty = (DUTY_PP | DUTY_PME);
4001         npme_major = cr->nnodes;
4002         npme_minor = 1;
4003
4004         if (inputrec->ePBC == epbcSCREW)
4005         {
4006             gmx_fatal(FARGS,
4007                       "pbc=%s is only implemented with domain decomposition",
4008                       epbc_names[inputrec->ePBC]);
4009         }
4010     }
4011
4012     if (PAR(cr))
4013     {
4014         /* After possible communicator splitting in make_dd_communicators.
4015          * we can set up the intra/inter node communication.
4016          */
4017         gmx_setup_nodecomm(fplog,cr);
4018     }
4019
4020     wcycle = wallcycle_init(fplog,resetstep,cr);
4021     if (PAR(cr))
4022     {
4023         /* Master synchronizes its value of reset_counters with all nodes
4024          * including PME only nodes */
4025         reset_counters = wcycle_get_reset_counters(wcycle);
4026         gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
4027         wcycle_set_reset_counters(wcycle, reset_counters);
4028     }
4029
4030
4031     snew(nrnb,1);
4032     if (cr->duty & DUTY_PP)
4033     {
4034         /* For domain decomposition we allocate dynamically
4035          * in dd_partition_system.
4036          */
4037         if (DOMAINDECOMP(cr))
4038         {
4039             bcast_state_setup(cr,state);
4040         }
4041         else
4042         {
4043             if (PAR(cr))
4044             {
4045                 if (!MASTER(cr))
4046                 {
4047                     snew(state,1);
4048                 }
4049                 bcast_state(cr,state,TRUE);
4050             }
4051         }
4052
4053         /* Dihedral Restraints */
4054         if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
4055         {
4056             init_dihres(fplog,mtop,inputrec,fcd);
4057         }
4058
4059         /* Initiate forcerecord */
4060         fr = mk_forcerec();
4061         init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
4062                       opt2fn("-table",nfile,fnm),
4063                       opt2fn("-tablep",nfile,fnm),
4064                       opt2fn("-tableb",nfile,fnm),FALSE,pforce);
4065
4066         /* version for PCA_NOT_READ_NODE (see md.c) */
4067         /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
4068           "nofile","nofile","nofile",FALSE,pforce);
4069           */
4070         fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
4071
4072         /* Initialize QM-MM */
4073         if(fr->bQMMM)
4074         {
4075             init_QMMMrec(cr,box,mtop,inputrec,fr);
4076         }
4077
4078         /* Initialize the mdatoms structure.
4079          * mdatoms is not filled with atom data,
4080          * as this can not be done now with domain decomposition.
4081          */
4082         mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
4083
4084         /* Initialize the virtual site communication */
4085         vsite = init_vsite(mtop,cr);
4086
4087         calc_shifts(box,fr->shift_vec);
4088
4089         /* With periodic molecules the charge groups should be whole at start up
4090          * and the virtual sites should not be far from their proper positions.
4091          */
4092         if (!inputrec->bContinuation && MASTER(cr) &&
4093             !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
4094         {
4095             /* Make molecules whole at start of run */
4096             if (fr->ePBC != epbcNONE)
4097             {
4098                 do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
4099             }
4100             if (vsite)
4101             {
4102                 /* Correct initial vsite positions are required
4103                  * for the initial distribution in the domain decomposition
4104                  * and for the initial shell prediction.
4105                  */
4106                 construct_vsites_mtop(fplog,vsite,mtop,state->x);
4107             }
4108         }
4109
4110         /* Initiate PPPM if necessary */
4111         if (fr->eeltype == eelPPPM)
4112         {
4113             if (mdatoms->nChargePerturbed)
4114             {
4115                 gmx_fatal(FARGS,"Free energy with %s is not implemented",
4116                           eel_names[fr->eeltype]);
4117             }
4118             status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
4119                                    getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
4120             if (status != 0)
4121             {
4122                 gmx_fatal(FARGS,"Error %d initializing PPPM",status);
4123             }
4124         }
4125
4126         if (EEL_PME(fr->eeltype))
4127         {
4128             ewaldcoeff = fr->ewaldcoeff;
4129             pmedata = &fr->pmedata;
4130         }
4131         else
4132         {
4133             pmedata = NULL;
4134         }
4135     }
4136     else
4137     {
4138         /* This is a PME only node */
4139
4140         /* We don't need the state */
4141         done_state(state);
4142
4143         ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
4144         snew(pmedata,1);
4145     }
4146
4147     /* Initiate PME if necessary,
4148      * either on all nodes or on dedicated PME nodes only. */
4149     if (EEL_PME(inputrec->coulombtype))
4150     {
4151         if (mdatoms)
4152         {
4153             nChargePerturbed = mdatoms->nChargePerturbed;
4154         }
4155         if (cr->npmenodes > 0)
4156         {
4157             /* The PME only nodes need to know nChargePerturbed */
4158             gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
4159         }
4160         if (cr->duty & DUTY_PME)
4161         {
4162             status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
4163                                   mtop ? mtop->natoms : 0,nChargePerturbed,
4164                                   (Flags & MD_REPRODUCIBLE));
4165             if (status != 0)
4166             {
4167                 gmx_fatal(FARGS,"Error %d initializing PME",status);
4168             }
4169         }
4170     }
4171
4172
4173 /*    if (integrator[inputrec->eI].func == do_md
4174 #ifdef GMX_OPENMM
4175         ||
4176         integrator[inputrec->eI].func == do_md_openmm
4177 #endif
4178         )
4179     {*/
4180         /* Turn on signal handling on all nodes */
4181         /*
4182          * (A user signal from the PME nodes (if any)
4183          * is communicated to the PP nodes.
4184          */
4185         signal_handler_install();
4186 /*    }*/
4187
4188     if (cr->duty & DUTY_PP)
4189     {
4190         if (inputrec->ePull != epullNO)
4191         {
4192             /* Initialize pull code */
4193             init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
4194                       EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
4195         }
4196
4197         constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
4198
4199         if (DOMAINDECOMP(cr))
4200         {
4201             dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
4202                             Flags & MD_DDBONDCHECK,fr->cginfo_mb);
4203
4204             set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
4205
4206             setup_dd_grid(fplog,cr->dd);
4207         }
4208
4209         /* Now do whatever the user wants us to do (how flexible...) */
4210         do_md_membed(fplog,cr,nfile,fnm,
4211                                       oenv,bVerbose,bCompact,
4212                                       nstglobalcomm,
4213                                       vsite,constr,
4214                                       nstepout,inputrec,mtop,
4215                                       fcd,state,
4216                                       mdatoms,nrnb,wcycle,ed,fr,
4217                                       repl_ex_nst,repl_ex_seed,
4218                                       cpt_period,max_hours,
4219                                       deviceOptions,
4220                                       Flags,
4221                                       &runtime,
4222                                       fac, r_ins, pos_ins, ins_at,
4223                                       xy_step, z_step, it_xy, it_z);
4224
4225         if (inputrec->ePull != epullNO)
4226         {
4227             finish_pull(fplog,inputrec->pull);
4228         }
4229     }
4230     else
4231     {
4232         /* do PME only */
4233         gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
4234     }
4235
4236     if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
4237     {
4238         /* Some timing stats */
4239         if (MASTER(cr))
4240         {
4241             if (runtime.proc == 0)
4242             {
4243                 runtime.proc = runtime.real;
4244             }
4245         }
4246         else
4247         {
4248             runtime.real = 0;
4249         }
4250     }
4251
4252     wallcycle_stop(wcycle,ewcRUN);
4253
4254     /* Finish up, write some stuff
4255      * if rerunMD, don't write last frame again
4256      */
4257     finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
4258                inputrec,nrnb,wcycle,&runtime,
4259                EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
4260
4261     /* Does what it says */
4262     print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
4263
4264     /* Close logfile already here if we were appending to it */
4265     if (MASTER(cr) && (Flags & MD_APPENDFILES))
4266     {
4267         gmx_log_close(fplog);
4268     }
4269
4270     if (pieces>1)
4271     {
4272         sfree(piecename);
4273     }
4274
4275     rc=(int)gmx_get_stop_condition();
4276
4277     return rc;
4278 }
4279
4280 int gmx_membed(int argc,char *argv[])
4281 {
4282         const char *desc[] = {
4283                         "g_membed embeds a membrane protein into an equilibrated lipid bilayer at the position",
4284                         "and orientation specified by the user.\n",
4285                         "\n",
4286                         "SHORT MANUAL\n------------\n",
4287                         "The user should merge the structure files of the protein and membrane (+solvent), creating a",
4288                         "single structure file with the protein overlapping the membrane at the desired position and",
4289                         "orientation. Box size should be taken from the membrane structure file. The corresponding topology",
4290                         "files should also be merged. Consecutively, create a tpr file (input for g_membed) from these files,"
4291                         "with the following options included in the mdp file.\n",
4292                         " - integrator      = md\n",
4293                         " - energygrp       = Protein (or other group that you want to insert)\n",
4294                         " - freezegrps      = Protein\n",
4295                         " - freezedim       = Y Y Y\n",
4296                         " - energygrp_excl  = Protein Protein\n",
4297                         "The output is a structure file containing the protein embedded in the membrane. If a topology",
4298                         "file is provided, the number of lipid and ",
4299                         "solvent molecules will be updated to match the new structure file.\n",
4300                         "For a more extensive manual see Wolf et al, J Comp Chem 31 (2010) 2169-2174, Appendix.\n",
4301                         "\n",
4302                         "SHORT METHOD DESCRIPTION\n",
4303                         "------------------------\n",
4304                         "1. The protein is resized around its center of mass by a factor -xy in the xy-plane",
4305                         "(the membrane plane) and a factor -z in the z-direction (if the size of the",
4306                         "protein in the z-direction is the same or smaller than the width of the membrane, a",
4307                         "-z value larger than 1 can prevent that the protein will be enveloped by the lipids).\n",
4308                         "2. All lipid and solvent molecules overlapping with the resized protein are removed. All",
4309                         "intraprotein interactions are turned off to prevent numerical issues for small values of -xy",
4310                         " or -z\n",
4311                         "3. One md step is performed.\n",
4312                         "4. The resize factor (-xy or -z) is incremented by a small amount ((1-xy)/nxy or (1-z)/nz) and the",
4313                         "protein is resized again around its center of mass. The resize factor for the xy-plane",
4314                         "is incremented first. The resize factor for the z-direction is not changed until the -xy factor",
4315                         "is 1 (thus after -nxy iteration).\n",
4316                         "5. Repeat step 3 and 4 until the protein reaches its original size (-nxy + -nz iterations).\n",
4317                         "For a more extensive method descrition see Wolf et al, J Comp Chem, 31 (2010) 2169-2174.\n",
4318                         "\n",
4319                         "NOTE\n----\n",
4320                         " - Protein can be any molecule you want to insert in the membrane.\n",
4321                         " - It is recommended to perform a short equilibration run after the embedding",
4322                         "(see Wolf et al, J Comp Chem 31 (2010) 2169-2174, to re-equilibrate the membrane. Clearly",
4323                         "protein equilibration might require longer.\n",
4324                         "\n"
4325         };
4326         t_commrec    *cr;
4327         t_filenm fnm[] = {
4328                         { efTPX, "-f",      "into_mem", ffREAD },
4329                         { efNDX, "-n",      "index",    ffOPTRD },
4330                         { efTOP, "-p",      "topol",    ffOPTRW },
4331                         { efTRN, "-o",      NULL,       ffWRITE },
4332                         { efXTC, "-x",      NULL,       ffOPTWR },
4333                         { efCPT, "-cpi",    NULL,       ffOPTRD },
4334                         { efCPT, "-cpo",    NULL,       ffOPTWR },
4335                         { efSTO, "-c",      "membedded",  ffWRITE },
4336                         { efEDR, "-e",      "ener",     ffWRITE },
4337                         { efLOG, "-g",      "md",       ffWRITE },
4338                         { efEDI, "-ei",     "sam",      ffOPTRD },
4339                         { efTRX, "-rerun",  "rerun",    ffOPTRD },
4340                         { efXVG, "-table",  "table",    ffOPTRD },
4341                         { efXVG, "-tablep", "tablep",   ffOPTRD },
4342                         { efXVG, "-tableb", "table",    ffOPTRD },
4343                         { efXVG, "-dhdl",   "dhdl",     ffOPTWR },
4344                         { efXVG, "-field",  "field",    ffOPTWR },
4345                         { efXVG, "-table",  "table",    ffOPTRD },
4346                         { efXVG, "-tablep", "tablep",   ffOPTRD },
4347                         { efXVG, "-tableb", "table",    ffOPTRD },
4348                         { efTRX, "-rerun",  "rerun",    ffOPTRD },
4349                         { efXVG, "-tpi",    "tpi",      ffOPTWR },
4350                         { efXVG, "-tpid",   "tpidist",  ffOPTWR },
4351                         { efEDI, "-ei",     "sam",      ffOPTRD },
4352                         { efEDO, "-eo",     "sam",      ffOPTWR },
4353                         { efGCT, "-j",      "wham",     ffOPTRD },
4354                         { efGCT, "-jo",     "bam",      ffOPTWR },
4355                         { efXVG, "-ffout",  "gct",      ffOPTWR },
4356                         { efXVG, "-devout", "deviatie", ffOPTWR },
4357                         { efXVG, "-runav",  "runaver",  ffOPTWR },
4358                         { efXVG, "-px",     "pullx",    ffOPTWR },
4359                         { efXVG, "-pf",     "pullf",    ffOPTWR },
4360                         { efMTX, "-mtx",    "nm",       ffOPTWR },
4361                         { efNDX, "-dn",     "dipole",   ffOPTWR }
4362         };
4363 #define NFILE asize(fnm)
4364
4365         /* Command line options ! */
4366         gmx_bool bCart        = FALSE;
4367         gmx_bool bPPPME       = FALSE;
4368         gmx_bool bPartDec     = FALSE;
4369         gmx_bool bDDBondCheck = TRUE;
4370         gmx_bool bDDBondComm  = TRUE;
4371         gmx_bool bVerbose     = FALSE;
4372         gmx_bool bCompact     = TRUE;
4373         gmx_bool bSepPot      = FALSE;
4374         gmx_bool bRerunVSite  = FALSE;
4375         gmx_bool bIonize      = FALSE;
4376         gmx_bool bConfout     = TRUE;
4377         gmx_bool bReproducible = FALSE;
4378
4379         int  npme=-1;
4380         int  nmultisim=0;
4381         int  nstglobalcomm=-1;
4382         int  repl_ex_nst=0;
4383         int  repl_ex_seed=-1;
4384         int  nstepout=100;
4385         int  nthreads=0; /* set to determine # of threads automatically */
4386         int  resetstep=-1;
4387
4388         rvec realddxyz={0,0,0};
4389         const char *ddno_opt[ddnoNR+1] =
4390         { NULL, "interleave", "pp_pme", "cartesian", NULL };
4391         const char *dddlb_opt[] =
4392         { NULL, "auto", "no", "yes", NULL };
4393         real rdd=0.0,rconstr=0.0,dlb_scale=0.8,pforce=-1;
4394         char *ddcsx=NULL,*ddcsy=NULL,*ddcsz=NULL;
4395         real cpt_period=15.0,max_hours=-1;
4396         gmx_bool bAppendFiles=TRUE,bAddPart=TRUE;
4397         gmx_bool bResetCountersHalfWay=FALSE;
4398         output_env_t oenv=NULL;
4399         const char *deviceOptions = "";
4400
4401         real xy_fac = 0.5;
4402         real xy_max = 1.0;
4403         real z_fac = 1.0;
4404         real z_max = 1.0;
4405         int it_xy = 1000;
4406         int it_z = 0;
4407         real probe_rad = 0.22;
4408         int low_up_rm = 0;
4409         int maxwarn=0;
4410         int pieces=1;
4411     gmx_bool bALLOW_ASYMMETRY=FALSE;
4412
4413
4414 /* arguments relevant to OPENMM only*/
4415 #ifdef GMX_OPENMM
4416     gmx_input("g_membed not functional in openmm");
4417 #endif
4418
4419         t_pargs pa[] = {
4420                         { "-xyinit",   FALSE, etREAL,  {&xy_fac},       "Resize factor for the protein in the xy dimension before starting embedding" },
4421                         { "-xyend",   FALSE, etREAL,  {&xy_max},                "Final resize factor in the xy dimension" },
4422                         { "-zinit",    FALSE, etREAL,  {&z_fac},                "Resize factor for the protein in the z dimension before starting embedding" },
4423                         { "-zend",    FALSE, etREAL,  {&z_max},                 "Final resize faction in the z dimension" },
4424                         { "-nxy",     FALSE,  etINT,  {&it_xy},         "Number of iteration for the xy dimension" },
4425                         { "-nz",      FALSE,  etINT,  {&it_z},          "Number of iterations for the z dimension" },
4426                         { "-rad",     FALSE, etREAL,  {&probe_rad},     "Probe radius to check for overlap between the group to embed and the membrane"},
4427                         { "-pieces",  FALSE,  etINT,  {&pieces},        "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
4428             { "-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." },
4429             { "-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." },
4430                         { "-maxwarn", FALSE, etINT, {&maxwarn},                 "Maximum number of warning allowed" },
4431   { "-pd",      FALSE, etBOOL,{&bPartDec},
4432     "HIDDENUse particle decompostion" },
4433   { "-dd",      FALSE, etRVEC,{&realddxyz},
4434     "HIDDENDomain decomposition grid, 0 is optimize" },
4435   { "-nt",      FALSE, etINT, {&nthreads},
4436     "HIDDENNumber of threads to start (0 is guess)" },
4437   { "-npme",    FALSE, etINT, {&npme},
4438     "HIDDENNumber of separate nodes to be used for PME, -1 is guess" },
4439   { "-ddorder", FALSE, etENUM, {ddno_opt},
4440     "HIDDENDD node order" },
4441   { "-ddcheck", FALSE, etBOOL, {&bDDBondCheck},
4442     "HIDDENCheck for all bonded interactions with DD" },
4443   { "-ddbondcomm", FALSE, etBOOL, {&bDDBondComm},
4444     "HIDDENUse special bonded atom communication when -rdd > cut-off" },
4445   { "-rdd",     FALSE, etREAL, {&rdd},
4446     "HIDDENThe maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates" },
4447   { "-rcon",    FALSE, etREAL, {&rconstr},
4448     "HIDDENMaximum distance for P-LINCS (nm), 0 is estimate" },
4449   { "-dlb",     FALSE, etENUM, {dddlb_opt},
4450     "HIDDENDynamic load balancing (with DD)" },
4451   { "-dds",     FALSE, etREAL, {&dlb_scale},
4452     "HIDDENMinimum allowed dlb scaling of the DD cell size" },
4453   { "-ddcsx",   FALSE, etSTR, {&ddcsx},
4454     "HIDDENThe DD cell sizes in x" },
4455   { "-ddcsy",   FALSE, etSTR, {&ddcsy},
4456     "HIDDENThe DD cell sizes in y" },
4457   { "-ddcsz",   FALSE, etSTR, {&ddcsz},
4458     "HIDDENThe DD cell sizes in z" },
4459   { "-gcom",    FALSE, etINT,{&nstglobalcomm},
4460     "HIDDENGlobal communication frequency" },
4461   { "-compact", FALSE, etBOOL,{&bCompact},
4462     "Write a compact log file" },
4463   { "-seppot",  FALSE, etBOOL, {&bSepPot},
4464     "HIDDENWrite separate V and dVdl terms for each interaction type and node to the log file(s)" },
4465   { "-pforce",  FALSE, etREAL, {&pforce},
4466     "HIDDENPrint all forces larger than this (kJ/mol nm)" },
4467   { "-reprod",  FALSE, etBOOL,{&bReproducible},
4468     "HIDDENTry to avoid optimizations that affect binary reproducibility" },
4469   { "-multi",   FALSE, etINT,{&nmultisim},
4470     "HIDDENDo multiple simulations in parallel" },
4471   { "-replex",  FALSE, etINT, {&repl_ex_nst},
4472     "HIDDENAttempt replica exchange every # steps" },
4473   { "-reseed",  FALSE, etINT, {&repl_ex_seed},
4474     "HIDDENSeed for replica exchange, -1 is generate a seed" },
4475   { "-rerunvsite", FALSE, etBOOL, {&bRerunVSite},
4476     "HIDDENRecalculate virtual site coordinates with -rerun" },
4477   { "-ionize",  FALSE, etBOOL,{&bIonize},
4478     "HIDDENDo a simulation including the effect of an X-Ray bombardment on your system" },
4479   { "-confout", TRUE, etBOOL, {&bConfout},
4480     "HIDDENWrite the last configuration with -c and force checkpointing at the last step" },
4481   { "-stepout", FALSE, etINT, {&nstepout},
4482     "HIDDENFrequency of writing the remaining runtime" },
4483   { "-resetstep", FALSE, etINT, {&resetstep},
4484     "HIDDENReset cycle counters after these many time steps" },
4485   { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
4486     "HIDDENReset the cycle counters after half the number of steps or halfway -maxh" },
4487   { "-v",       FALSE, etBOOL,{&bVerbose},
4488     "Be loud and noisy" },
4489   { "-maxh",   FALSE, etREAL, {&max_hours},
4490     "HIDDENTerminate after 0.99 times this time (hours)" },
4491   { "-cpt",     FALSE, etREAL, {&cpt_period},
4492     "HIDDENCheckpoint interval (minutes)" },
4493   { "-append",  FALSE, etBOOL, {&bAppendFiles},
4494     "HIDDENAppend to previous output files when continuing from checkpoint" },
4495   { "-addpart",  FALSE, etBOOL, {&bAddPart},
4496     "HIDDENAdd the simulation part number to all output files when continuing from checkpoint" },
4497         };
4498         gmx_edsam_t  ed;
4499         unsigned long Flags, PCA_Flags;
4500         ivec     ddxyz;
4501         int      dd_node_order;
4502         gmx_bool     HaveCheckpoint;
4503         FILE     *fplog,*fptest;
4504         int      sim_part,sim_part_fn;
4505         const char *part_suffix=".part";
4506         char     suffix[STRLEN];
4507         int      rc;
4508
4509
4510         cr = init_par(&argc,&argv);
4511
4512         PCA_Flags = (PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_CAN_SET_DEFFNM
4513                         | (MASTER(cr) ? 0 : PCA_QUIET));
4514
4515
4516         /* Comment this in to do fexist calls only on master
4517          * works not with rerun or tables at the moment
4518          * also comment out the version of init_forcerec in md.c
4519          * with NULL instead of opt2fn
4520          */
4521         /*
4522    if (!MASTER(cr))
4523    {
4524    PCA_Flags |= PCA_NOT_READ_NODE;
4525    }
4526          */
4527
4528         parse_common_args(&argc,argv,PCA_Flags, NFILE,fnm,asize(pa),pa,
4529                         asize(desc),desc,0,NULL, &oenv);
4530
4531         /* we set these early because they might be used in init_multisystem()
4532    Note that there is the potential for npme>nnodes until the number of
4533    threads is set later on, if there's thread parallelization. That shouldn't
4534    lead to problems. */
4535         dd_node_order = nenum(ddno_opt);
4536         cr->npmenodes = npme;
4537
4538 #ifdef GMX_THREADS
4539         /* now determine the number of threads automatically. The threads are
4540    only started at mdrunner_threads, though. */
4541         if (nthreads<1)
4542         {
4543                 nthreads=tMPI_Get_recommended_nthreads();
4544         }
4545 #else
4546         nthreads=1;
4547 #endif
4548
4549
4550         if (repl_ex_nst != 0 && nmultisim < 2)
4551                 gmx_fatal(FARGS,"Need at least two replicas for replica exchange (option -multi)");
4552
4553         if (nmultisim > 1) {
4554 #ifndef GMX_THREADS
4555                 init_multisystem(cr,nmultisim,NFILE,fnm,TRUE);
4556 #else
4557                 gmx_fatal(FARGS,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
4558 #endif
4559         }
4560
4561         /* Check if there is ANY checkpoint file available */
4562         sim_part    = 1;
4563         sim_part_fn = sim_part;
4564         if (opt2bSet("-cpi",NFILE,fnm))
4565         {
4566                 bAppendFiles =
4567                         read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE,fnm,cr),
4568                                                         &sim_part_fn,NULL,cr,
4569                                                         bAppendFiles,NFILE,fnm,
4570                                                         part_suffix,&bAddPart);
4571                 if (sim_part_fn==0 && MASTER(cr))
4572                 {
4573                         fprintf(stdout,"No previous checkpoint file present, assuming this is a new run.\n");
4574                 }
4575                 else
4576                 {
4577                         sim_part = sim_part_fn + 1;
4578                 }
4579         }
4580         else
4581         {
4582                 bAppendFiles = FALSE;
4583         }
4584
4585         if (!bAppendFiles)
4586         {
4587                 sim_part_fn = sim_part;
4588         }
4589
4590         if (bAddPart && sim_part_fn > 1)
4591         {
4592                 /* This is a continuation run, rename trajectory output files
4593        (except checkpoint files) */
4594                 /* create new part name first (zero-filled) */
4595                 sprintf(suffix,"%s%04d",part_suffix,sim_part_fn);
4596
4597                 add_suffix_to_output_names(fnm,NFILE,suffix);
4598                 fprintf(stdout,"Checkpoint file is from part %d, new output files will be suffixed '%s'.\n",sim_part-1,suffix);
4599         }
4600
4601         Flags = opt2bSet("-rerun",NFILE,fnm) ? MD_RERUN : 0;
4602         Flags = Flags | (bSepPot       ? MD_SEPPOT       : 0);
4603         Flags = Flags | (bIonize       ? MD_IONIZE       : 0);
4604         Flags = Flags | (bPartDec      ? MD_PARTDEC      : 0);
4605         Flags = Flags | (bDDBondCheck  ? MD_DDBONDCHECK  : 0);
4606         Flags = Flags | (bDDBondComm   ? MD_DDBONDCOMM   : 0);
4607         Flags = Flags | (bConfout      ? MD_CONFOUT      : 0);
4608         Flags = Flags | (bRerunVSite   ? MD_RERUN_VSITE  : 0);
4609         Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
4610         Flags = Flags | (bAppendFiles  ? MD_APPENDFILES  : 0);
4611         Flags = Flags | (sim_part>1    ? MD_STARTFROMCPT : 0);
4612         Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
4613
4614
4615         /* We postpone opening the log file if we are appending, so we can
4616    first truncate the old log file and append to the correct position
4617    there instead.  */
4618         if ((MASTER(cr) || bSepPot) && !bAppendFiles)
4619         {
4620                 gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,!bSepPot,Flags,&fplog);
4621                 CopyRight(fplog,argv[0]);
4622                 please_cite(fplog,"Hess2008b");
4623                 please_cite(fplog,"Spoel2005a");
4624                 please_cite(fplog,"Lindahl2001a");
4625                 please_cite(fplog,"Berendsen95a");
4626         }
4627         else
4628         {
4629                 fplog = NULL;
4630         }
4631
4632         ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
4633         ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
4634         ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
4635
4636         /* even if nthreads = 1, we still call this one */
4637
4638         rc = mdrunner_membed(fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
4639                         nstglobalcomm,
4640                         ddxyz, dd_node_order, rdd, rconstr, dddlb_opt[0], dlb_scale,
4641                         ddcsx, ddcsy, ddcsz, nstepout, resetstep, nmultisim, repl_ex_nst,
4642                         repl_ex_seed, pforce, cpt_period, max_hours, deviceOptions, Flags,
4643                         xy_fac,xy_max,z_fac,z_max,
4644                         it_xy,it_z,probe_rad,low_up_rm,
4645                         pieces,bALLOW_ASYMMETRY,maxwarn);
4646
4647         if (gmx_parallel_env_initialized())
4648                 gmx_finalize();
4649
4650         if (MULTIMASTER(cr)) {
4651                 thanx(stderr);
4652         }
4653
4654         /* Log file has to be closed in mdrunner if we are appending to it
4655    (fplog not set here) */
4656         fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
4657
4658         if (MASTER(cr) && !bAppendFiles)
4659         {
4660                 gmx_log_close(fplog);
4661         }
4662
4663         return rc;
4664 }