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         if (opt2bSet("-cpi",nfile,fnm))
2081         {
2082             /* Update mdebin with energy history if appending to output files */
2083             if ( Flags & MD_APPENDFILES )
2084             {
2085                 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
2086             }
2087             else
2088             {
2089                 /* We might have read an energy history from checkpoint,
2090                  * free the allocated memory and reset the counts.
2091                  */
2092                 done_energyhistory(&state_global->enerhist);
2093                 init_energyhistory(&state_global->enerhist);
2094             }
2095         }
2096         /* Set the initial energy history in state by updating once */
2097         update_energyhistory(&state_global->enerhist,mdebin);
2098     }
2099
2100     if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
2101         /* Set the random state if we read a checkpoint file */
2102         set_stochd_state(upd,state);
2103     }
2104
2105     /* Initialize constraints */
2106     if (constr) {
2107         if (!DOMAINDECOMP(cr))
2108             set_constraints(constr,top,ir,mdatoms,cr);
2109     }
2110
2111     /* Check whether we have to GCT stuff */
2112  /*   bTCR = ftp2bSet(efGCT,nfile,fnm);
2113     if (bTCR) {
2114         if (MASTER(cr)) {
2115             fprintf(stderr,"Will do General Coupling Theory!\n");
2116         }
2117         gnx = top_global->mols.nr;
2118         snew(grpindex,gnx);
2119         for(i=0; (i<gnx); i++) {
2120             grpindex[i] = i;
2121         }
2122     }*/
2123
2124 /*    if (repl_ex_nst > 0 && MASTER(cr))
2125         repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
2126                                         repl_ex_nst,repl_ex_seed);*/
2127
2128     if (!ir->bContinuation && !bRerunMD)
2129     {
2130         if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
2131         {
2132             /* Set the velocities of frozen particles to zero */
2133             for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
2134             {
2135                 for(m=0; m<DIM; m++)
2136                 {
2137                     if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
2138                     {
2139                         state->v[i][m] = 0;
2140                     }
2141                 }
2142             }
2143         }
2144
2145         if (constr)
2146         {
2147             /* Constrain the initial coordinates and velocities */
2148             do_constrain_first(fplog,constr,ir,mdatoms,state,f,
2149                                graph,cr,nrnb,fr,top,shake_vir);
2150         }
2151         if (vsite)
2152         {
2153             /* Construct the virtual sites for the initial configuration */
2154             construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
2155                              top->idef.iparams,top->idef.il,
2156                              fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2157         }
2158     }
2159
2160     debug_gmx();
2161
2162     /* I'm assuming we need global communication the first time! MRS */
2163     cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
2164                   | (bVV ? CGLO_PRESSURE:0)
2165                   | (bVV ? CGLO_CONSTRAINT:0)
2166                   | (bRerunMD ? CGLO_RERUNMD:0)
2167                   | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
2168
2169     bSumEkinhOld = FALSE;
2170     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2171                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2172                     constr,NULL,FALSE,state->box,
2173                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
2174     if (ir->eI == eiVVAK) {
2175         /* a second call to get the half step temperature initialized as well */
2176         /* we do the same call as above, but turn the pressure off -- internally, this
2177            is recognized as a velocity verlet half-step kinetic energy calculation.
2178            This minimized excess variables, but perhaps loses some logic?*/
2179
2180         compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2181                         wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2182                         constr,NULL,FALSE,state->box,
2183                         top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2184                         cglo_flags &~ CGLO_PRESSURE);
2185     }
2186
2187     /* Calculate the initial half step temperature, and save the ekinh_old */
2188     if (!(Flags & MD_STARTFROMCPT))
2189     {
2190         for(i=0; (i<ir->opts.ngtc); i++)
2191         {
2192             copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
2193         }
2194     }
2195     if (ir->eI != eiVV)
2196     {
2197         enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
2198                                      and there is no previous step */
2199     }
2200     temp0 = enerd->term[F_TEMP];
2201
2202     /* if using an iterative algorithm, we need to create a working directory for the state. */
2203     if (bIterations)
2204     {
2205             bufstate = init_bufstate(state);
2206     }
2207     if (bFFscan)
2208     {
2209         snew(xcopy,state->natoms);
2210         snew(vcopy,state->natoms);
2211         copy_rvecn(state->x,xcopy,0,state->natoms);
2212         copy_rvecn(state->v,vcopy,0,state->natoms);
2213         copy_mat(state->box,boxcopy);
2214     }
2215
2216     /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
2217        temperature control */
2218     trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
2219
2220     if (MASTER(cr))
2221     {
2222         if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
2223         {
2224             fprintf(fplog,
2225                     "RMS relative constraint deviation after constraining: %.2e\n",
2226                     constr_rmsd(constr,FALSE));
2227         }
2228         fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
2229         if (bRerunMD)
2230         {
2231             fprintf(stderr,"starting md rerun '%s', reading coordinates from"
2232                     " input trajectory '%s'\n\n",
2233                     *(top_global->name),opt2fn("-rerun",nfile,fnm));
2234             if (bVerbose)
2235             {
2236                 fprintf(stderr,"Calculated time to finish depends on nsteps from "
2237                         "run input file,\nwhich may not correspond to the time "
2238                         "needed to process input trajectory.\n\n");
2239             }
2240         }
2241         else
2242         {
2243             char tbuf[20];
2244             fprintf(stderr,"starting mdrun '%s'\n",
2245                     *(top_global->name));
2246             if (ir->nsteps >= 0)
2247             {
2248                 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
2249             }
2250             else
2251             {
2252                 sprintf(tbuf,"%s","infinite");
2253             }
2254             if (ir->init_step > 0)
2255             {
2256                 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
2257                         gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
2258                         gmx_step_str(ir->init_step,sbuf2),
2259                         ir->init_step*ir->delta_t);
2260             }
2261             else
2262             {
2263                 fprintf(stderr,"%s steps, %s ps.\n",
2264                         gmx_step_str(ir->nsteps,sbuf),tbuf);
2265             }
2266         }
2267         fprintf(fplog,"\n");
2268     }
2269
2270     /* Set and write start time */
2271     runtime_start(runtime);
2272     print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
2273     wallcycle_start(wcycle,ewcRUN);
2274     if (fplog)
2275         fprintf(fplog,"\n");
2276
2277     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
2278 /*#ifdef GMX_FAHCORE
2279     chkpt_ret=fcCheckPointParallel( cr->nodeid,
2280                                     NULL,0);
2281     if ( chkpt_ret == 0 )
2282         gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
2283 #endif*/
2284
2285     debug_gmx();
2286     /***********************************************************
2287      *
2288      *             Loop over MD steps
2289      *
2290      ************************************************************/
2291
2292     /* if rerunMD then read coordinates and velocities from input trajectory */
2293     if (bRerunMD)
2294     {
2295         if (getenv("GMX_FORCE_UPDATE"))
2296         {
2297             bForceUpdate = TRUE;
2298         }
2299
2300         bNotLastFrame = read_first_frame(oenv,&status,
2301                                          opt2fn("-rerun",nfile,fnm),
2302                                          &rerun_fr,TRX_NEED_X | TRX_READ_V);
2303         if (rerun_fr.natoms != top_global->natoms)
2304         {
2305             gmx_fatal(FARGS,
2306                       "Number of atoms in trajectory (%d) does not match the "
2307                       "run input file (%d)\n",
2308                       rerun_fr.natoms,top_global->natoms);
2309         }
2310         if (ir->ePBC != epbcNONE)
2311         {
2312             if (!rerun_fr.bBox)
2313             {
2314                 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);
2315             }
2316             if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
2317             {
2318                 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
2319             }
2320
2321             /* Set the shift vectors.
2322              * Necessary here when have a static box different from the tpr box.
2323              */
2324             calc_shifts(rerun_fr.box,fr->shift_vec);
2325         }
2326     }
2327
2328     /* loop over MD steps or if rerunMD to end of input trajectory */
2329     bFirstStep = TRUE;
2330     /* Skip the first Nose-Hoover integration when we get the state from tpx */
2331     bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
2332     bInitStep = bFirstStep && (bStateFromTPX || bVV);
2333     bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
2334     bLastStep    = FALSE;
2335     bSumEkinhOld = FALSE;
2336     bExchanged   = FALSE;
2337
2338     init_global_signals(&gs,cr,ir,repl_ex_nst);
2339
2340     step = ir->init_step;
2341     step_rel = 0;
2342
2343     if (ir->nstlist == -1)
2344     {
2345         init_nlistheuristics(&nlh,bGStatEveryStep,step);
2346     }
2347
2348     bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
2349     while (!bLastStep || (bRerunMD && bNotLastFrame)) {
2350
2351         wallcycle_start(wcycle,ewcSTEP);
2352
2353         GMX_MPE_LOG(ev_timestep1);
2354
2355         if (bRerunMD) {
2356             if (rerun_fr.bStep) {
2357                 step = rerun_fr.step;
2358                 step_rel = step - ir->init_step;
2359             }
2360             if (rerun_fr.bTime) {
2361                 t = rerun_fr.time;
2362             }
2363             else
2364             {
2365                 t = step;
2366             }
2367         }
2368         else
2369         {
2370             bLastStep = (step_rel == ir->nsteps);
2371             t = t0 + step*ir->delta_t;
2372         }
2373
2374         if (ir->efep != efepNO)
2375         {
2376             if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
2377             {
2378                 state_global->lambda = rerun_fr.lambda;
2379             }
2380             else
2381             {
2382                 state_global->lambda = lam0 + step*ir->delta_lambda;
2383             }
2384             state->lambda = state_global->lambda;
2385             bDoDHDL = do_per_step(step,ir->nstdhdl);
2386         }
2387
2388         if (bSimAnn)
2389         {
2390             update_annealing_target_temp(&(ir->opts),t);
2391         }
2392
2393         if (bRerunMD)
2394         {
2395             if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
2396             {
2397                 for(i=0; i<state_global->natoms; i++)
2398                 {
2399                     copy_rvec(rerun_fr.x[i],state_global->x[i]);
2400                 }
2401                 if (rerun_fr.bV)
2402                 {
2403                     for(i=0; i<state_global->natoms; i++)
2404                     {
2405                         copy_rvec(rerun_fr.v[i],state_global->v[i]);
2406                     }
2407                 }
2408                 else
2409                 {
2410                     for(i=0; i<state_global->natoms; i++)
2411                     {
2412                         clear_rvec(state_global->v[i]);
2413                     }
2414                     if (bRerunWarnNoV)
2415                     {
2416                         fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
2417                                 "         Ekin, temperature and pressure are incorrect,\n"
2418                                 "         the virial will be incorrect when constraints are present.\n"
2419                                 "\n");
2420                         bRerunWarnNoV = FALSE;
2421                     }
2422                 }
2423             }
2424             copy_mat(rerun_fr.box,state_global->box);
2425             copy_mat(state_global->box,state->box);
2426
2427             if (vsite && (Flags & MD_RERUN_VSITE))
2428             {
2429                 if (DOMAINDECOMP(cr))
2430                 {
2431                     gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
2432                 }
2433                 if (graph)
2434                 {
2435                     /* Following is necessary because the graph may get out of sync
2436                      * with the coordinates if we only have every N'th coordinate set
2437                      */
2438                     mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
2439                     shift_self(graph,state->box,state->x);
2440                 }
2441                 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2442                                  top->idef.iparams,top->idef.il,
2443                                  fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2444                 if (graph)
2445                 {
2446                     unshift_self(graph,state->box,state->x);
2447                 }
2448             }
2449         }
2450
2451         /* Stop Center of Mass motion */
2452         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
2453
2454         /* Copy back starting coordinates in case we're doing a forcefield scan */
2455         if (bFFscan)
2456         {
2457             for(ii=0; (ii<state->natoms); ii++)
2458             {
2459                 copy_rvec(xcopy[ii],state->x[ii]);
2460                 copy_rvec(vcopy[ii],state->v[ii]);
2461             }
2462             copy_mat(boxcopy,state->box);
2463         }
2464
2465         if (bRerunMD)
2466         {
2467             /* for rerun MD always do Neighbour Searching */
2468             bNS = (bFirstStep || ir->nstlist != 0);
2469             bNStList = bNS;
2470         }
2471         else
2472         {
2473             /* Determine whether or not to do Neighbour Searching and LR */
2474             bNStList = (ir->nstlist > 0  && step % ir->nstlist == 0);
2475
2476             bNS = (bFirstStep || bExchanged || bNStList ||
2477                    (ir->nstlist == -1 && nlh.nabnsb > 0));
2478
2479             if (bNS && ir->nstlist == -1)
2480             {
2481                 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
2482             }
2483         }
2484
2485         /* < 0 means stop at next step, > 0 means stop at next NS step */
2486         if ( (gs.set[eglsSTOPCOND] < 0 ) ||
2487              ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
2488         {
2489             bLastStep = TRUE;
2490         }
2491
2492         /* Determine whether or not to update the Born radii if doing GB */
2493         bBornRadii=bFirstStep;
2494         if (ir->implicit_solvent && (step % ir->nstgbradii==0))
2495         {
2496             bBornRadii=TRUE;
2497         }
2498
2499         do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
2500         do_verbose = bVerbose &&
2501                   (step % stepout == 0 || bFirstStep || bLastStep);
2502
2503         if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
2504         {
2505             if (bRerunMD)
2506             {
2507                 bMasterState = TRUE;
2508             }
2509             else
2510             {
2511                 bMasterState = FALSE;
2512                 /* Correct the new box if it is too skewed */
2513                 if (DYNAMIC_BOX(*ir))
2514                 {
2515                     if (correct_box(fplog,step,state->box,graph))
2516                     {
2517                         bMasterState = TRUE;
2518                     }
2519                 }
2520                 if (DOMAINDECOMP(cr) && bMasterState)
2521                 {
2522                     dd_collect_state(cr->dd,state,state_global);
2523                 }
2524             }
2525
2526             if (DOMAINDECOMP(cr))
2527             {
2528                 /* Repartition the domain decomposition */
2529                 wallcycle_start(wcycle,ewcDOMDEC);
2530                 dd_partition_system(fplog,step,cr,
2531                                     bMasterState,nstglobalcomm,
2532                                     state_global,top_global,ir,
2533                                     state,&f,mdatoms,top,fr,
2534                                     vsite,shellfc,constr,
2535                                     nrnb,wcycle,do_verbose);
2536                 wallcycle_stop(wcycle,ewcDOMDEC);
2537                 /* If using an iterative integrator, reallocate space to match the decomposition */
2538             }
2539         }
2540
2541         if (MASTER(cr) && do_log && !bFFscan)
2542         {
2543             print_ebin_header(fplog,step,t,state->lambda);
2544         }
2545
2546         if (ir->efep != efepNO)
2547         {
2548             update_mdatoms(mdatoms,state->lambda);
2549         }
2550
2551         if (bRerunMD && rerun_fr.bV)
2552         {
2553
2554             /* We need the kinetic energy at minus the half step for determining
2555              * the full step kinetic energy and possibly for T-coupling.*/
2556             /* This may not be quite working correctly yet . . . . */
2557             compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2558                             wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
2559                             constr,NULL,FALSE,state->box,
2560                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2561                             CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
2562         }
2563         clear_mat(force_vir);
2564
2565         /* Ionize the atoms if necessary */
2566 /*        if (bIonize)
2567         {
2568             ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
2569                    mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
2570         }*/
2571
2572         /* Update force field in ffscan program */
2573 /*        if (bFFscan)
2574         {
2575             if (update_forcefield(fplog,
2576                                   nfile,fnm,fr,
2577                                   mdatoms->nr,state->x,state->box)) {
2578                 if (gmx_parallel_env_initialized())
2579                 {
2580                     gmx_finalize();
2581                 }
2582                 exit(0);
2583             }
2584         }*/
2585
2586         GMX_MPE_LOG(ev_timestep2);
2587
2588         /* We write a checkpoint at this MD step when:
2589          * either at an NS step when we signalled through gs,
2590          * or at the last step (but not when we do not want confout),
2591          * but never at the first step or with rerun.
2592          */
2593 /*        bCPT = (((gs.set[eglsCHKPT] && bNS) ||
2594                  (bLastStep && (Flags & MD_CONFOUT))) &&
2595                 step > ir->init_step && !bRerunMD);
2596         if (bCPT)
2597         {
2598             gs.set[eglsCHKPT] = 0;
2599         }*/
2600
2601         /* Determine the energy and pressure:
2602          * at nstcalcenergy steps and at energy output steps (set below).
2603          */
2604         bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2605         bCalcEnerPres = bNstEner;
2606
2607         /* Do we need global communication ? */
2608         bGStat = (bCalcEnerPres || bStopCM ||
2609                   (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2610
2611         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2612
2613         if (do_ene || do_log)
2614         {
2615             bCalcEnerPres = TRUE;
2616             bGStat    = TRUE;
2617         }
2618
2619         /* these CGLO_ options remain the same throughout the iteration */
2620         cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
2621                       (bStopCM ? CGLO_STOPCM : 0) |
2622                       (bGStat ? CGLO_GSTAT : 0)
2623             );
2624
2625         force_flags = (GMX_FORCE_STATECHANGED |
2626                        ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
2627                        GMX_FORCE_ALLFORCES |
2628                        (bNStList ? GMX_FORCE_DOLR : 0) |
2629                        GMX_FORCE_SEPLRF |
2630                        (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
2631                        (bDoDHDL ? GMX_FORCE_DHDL : 0)
2632             );
2633
2634         if (shellfc)
2635         {
2636             /* Now is the time to relax the shells */
2637             count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
2638                                       ir,bNS,force_flags,
2639                                       bStopCM,top,top_global,
2640                                       constr,enerd,fcd,
2641                                       state,f,force_vir,mdatoms,
2642                                       nrnb,wcycle,graph,groups,
2643                                       shellfc,fr,bBornRadii,t,mu_tot,
2644                                       state->natoms,&bConverged,vsite,
2645                                       outf->fp_field);
2646             tcount+=count;
2647
2648             if (bConverged)
2649             {
2650                 nconverged++;
2651             }
2652         }
2653         else
2654         {
2655             /* The coordinates (x) are shifted (to get whole molecules)
2656              * in do_force.
2657              * This is parallellized as well, and does communication too.
2658              * Check comments in sim_util.c
2659              */
2660
2661             do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
2662                      state->box,state->x,&state->hist,
2663                      f,force_vir,mdatoms,enerd,fcd,
2664                      state->lambda,graph,
2665                      fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
2666                      (bNS ? GMX_FORCE_NS : 0) | force_flags);
2667         }
2668
2669         GMX_BARRIER(cr->mpi_comm_mygroup);
2670
2671  /*       if (bTCR)
2672         {
2673             mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
2674                                    mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
2675         }
2676
2677         if (bTCR && bFirstStep)
2678         {
2679             tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
2680             fprintf(fplog,"Done init_coupling\n");
2681             fflush(fplog);
2682         }*/
2683
2684         /*  ############### START FIRST UPDATE HALF-STEP ############### */
2685
2686         if (bVV && !bStartingFromCpt && !bRerunMD)
2687         {
2688             if (ir->eI == eiVV)
2689             {
2690                 if (bInitStep)
2691                 {
2692                     /* if using velocity verlet with full time step Ekin,
2693                      * take the first half step only to compute the
2694                      * virial for the first step. From there,
2695                      * revert back to the initial coordinates
2696                      * so that the input is actually the initial step.
2697                      */
2698                     copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
2699                 }
2700
2701                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
2702                 if (!bInitStep)
2703                 {
2704                   trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
2705                 }
2706
2707                 if (ir->eI == eiVVAK)
2708                 {
2709                   update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2710                 }
2711
2712                 update_coords(fplog,step,ir,mdatoms,state,
2713                               f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2714                               ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
2715                               cr,nrnb,constr,&top->idef);
2716
2717                 if (bIterations)
2718                 {
2719                     gmx_iterate_init(&iterate,bIterations && !bInitStep);
2720                 }
2721                 /* for iterations, we save these vectors, as we will be self-consistently iterating
2722                    the calculations */
2723                 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
2724
2725                 /* save the state */
2726                 if (bIterations && iterate.bIterate) {
2727                     copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2728                 }
2729             }
2730
2731             bFirstIterate = TRUE;
2732             while (bFirstIterate || (bIterations && iterate.bIterate))
2733             {
2734                 if (bIterations && iterate.bIterate)
2735                 {
2736                     copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2737                     if (bFirstIterate && bTrotter)
2738                     {
2739                         /* The first time through, we need a decent first estimate
2740                            of veta(t+dt) to compute the constraints.  Do
2741                            this by computing the box volume part of the
2742                            trotter integration at this time. Nothing else
2743                            should be changed by this routine here.  If
2744                            !(first time), we start with the previous value
2745                            of veta.  */
2746
2747                         veta_save = state->veta;
2748                         trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
2749                         vetanew = state->veta;
2750                         state->veta = veta_save;
2751                     }
2752                 }
2753
2754                 bOK = TRUE;
2755                 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) {  /* Why is rerun_fr.bV here?  Unclear. */
2756                     dvdl = 0;
2757
2758                     update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2759                                        &top->idef,shake_vir,NULL,
2760                                        cr,nrnb,wcycle,upd,constr,
2761                                        bInitStep,TRUE,bCalcEnerPres,vetanew);
2762
2763                     if (!bOK && !bFFscan)
2764                     {
2765                         gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2766                     }
2767
2768                 }
2769                 else if (graph)
2770                 { /* Need to unshift here if a do_force has been
2771                      called in the previous step */
2772                     unshift_self(graph,state->box,state->x);
2773                 }
2774
2775
2776                 if (bVV) {
2777                     /* if VV, compute the pressure and constraints */
2778                     /* if VV2, the pressure and constraints only if using pressure control.*/
2779                     bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir));
2780                     bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));
2781                     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2782                                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2783                                     constr,NULL,FALSE,state->box,
2784                                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2785                                     cglo_flags
2786                                     | CGLO_ENERGY
2787                                     | (bTemp ? CGLO_TEMPERATURE:0)
2788                                     | (bPres ? CGLO_PRESSURE : 0)
2789                                     | (bPres ? CGLO_CONSTRAINT : 0)
2790                                     | (iterate.bIterate ? CGLO_ITERATE : 0)
2791                                     | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2792                                     | CGLO_SCALEEKIN
2793                         );
2794                 }
2795                 /* explanation of above:
2796                    a) We compute Ekin at the full time step
2797                    if 1) we are using the AveVel Ekin, and it's not the
2798                    initial step, or 2) if we are using AveEkin, but need the full
2799                    time step kinetic energy for the pressure.
2800                    b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2801                    EkinAveVel because it's needed for the pressure */
2802
2803                 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2804                 if (bVV && !bInitStep)
2805                 {
2806                   trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ2);
2807                 }
2808
2809                 if (bIterations &&
2810                     done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2811                                    state->veta,&vetanew))
2812                 {
2813                     break;
2814                 }
2815                 bFirstIterate = FALSE;
2816             }
2817
2818             if (bTrotter && !bInitStep) {
2819                 copy_mat(shake_vir,state->svir_prev);
2820                 copy_mat(force_vir,state->fvir_prev);
2821                 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
2822                     /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2823                     enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
2824                     enerd->term[F_EKIN] = trace(ekind->ekin);
2825                 }
2826             }
2827             /* if it's the initial step, we performed this first step just to get the constraint virial */
2828             if (bInitStep && ir->eI==eiVV) {
2829                 copy_rvecn(cbuf,state->v,0,state->natoms);
2830             }
2831
2832             if (fr->bSepDVDL && fplog && do_log)
2833             {
2834                 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2835             }
2836             enerd->term[F_DHDL_CON] += dvdl;
2837
2838             GMX_MPE_LOG(ev_timestep1);
2839
2840         }
2841
2842         /* MRS -- now done iterating -- compute the conserved quantity */
2843         if (bVV) {
2844             last_conserved = 0;
2845             if (IR_NVT_TROTTER(ir) || IR_NPT_TROTTER(ir))
2846             {
2847                 last_conserved =
2848                     NPT_energy(ir,state,&MassQ);
2849                 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
2850                 {
2851                     last_conserved -= enerd->term[F_DISPCORR];
2852                 }
2853             }
2854             if (ir->eI==eiVV) {
2855                 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
2856             }
2857         }
2858
2859         /* ########  END FIRST UPDATE STEP  ############## */
2860         /* ########  If doing VV, we now have v(dt) ###### */
2861
2862         /* ################## START TRAJECTORY OUTPUT ################# */
2863
2864         /* Now we have the energies and forces corresponding to the
2865          * coordinates at time t. We must output all of this before
2866          * the update.
2867          * for RerunMD t is read from input trajectory
2868          */
2869         GMX_MPE_LOG(ev_output_start);
2870
2871         mdof_flags = 0;
2872         if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
2873         if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
2874         if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
2875         if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
2876 /*        if (bCPT) { mdof_flags |= MDOF_CPT; };*/
2877
2878 #ifdef GMX_FAHCORE
2879         if (MASTER(cr))
2880             fcReportProgress( ir->nsteps, step );
2881
2882         if (bLastStep)
2883         {
2884             /* Enforce writing positions and velocities at end of run */
2885             mdof_flags |= (MDOF_X | MDOF_V);
2886         }
2887             /* sync bCPT and fc record-keeping */
2888 /*            if (bCPT && MASTER(cr))
2889                 fcRequestCheckPoint();*/
2890 #endif
2891
2892         if (mdof_flags != 0)
2893         {
2894             wallcycle_start(wcycle,ewcTRAJ);
2895 /*            if (bCPT)
2896             {
2897                 if (state->flags & (1<<estLD_RNG))
2898                 {
2899                     get_stochd_state(upd,state);
2900                 }
2901                 if (MASTER(cr))
2902                 {
2903                     if (bSumEkinhOld)
2904                     {
2905                         state_global->ekinstate.bUpToDate = FALSE;
2906                     }
2907                     else
2908                     {
2909                         update_ekinstate(&state_global->ekinstate,ekind);
2910                         state_global->ekinstate.bUpToDate = TRUE;
2911                     }
2912                     update_energyhistory(&state_global->enerhist,mdebin);
2913                 }
2914             }*/
2915             write_traj(fplog,cr,outf,mdof_flags,top_global,
2916                        step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2917 /*            if (bCPT)
2918             {
2919                 nchkpt++;
2920                 bCPT = FALSE;
2921             }*/
2922             debug_gmx();
2923             if (bLastStep && step_rel == ir->nsteps &&
2924                 (Flags & MD_CONFOUT) && MASTER(cr) &&
2925                 !bRerunMD && !bFFscan)
2926             {
2927                 /* x and v have been collected in write_traj,
2928                  * because a checkpoint file will always be written
2929                  * at the last step.
2930                  */
2931                 fprintf(stderr,"\nWriting final coordinates.\n");
2932                 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2933                     DOMAINDECOMP(cr))
2934                 {
2935                     /* Make molecules whole only for confout writing */
2936                     do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2937                 }
2938 /*                write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2939                                     *top_global->name,top_global,
2940                                     state_global->x,state_global->v,
2941                                     ir->ePBC,state->box);*/
2942                 debug_gmx();
2943             }
2944             wallcycle_stop(wcycle,ewcTRAJ);
2945         }
2946         GMX_MPE_LOG(ev_output_finish);
2947
2948         /* kludge -- virial is lost with restart for NPT control. Must restart */
2949         if (bStartingFromCpt && bVV)
2950         {
2951             copy_mat(state->svir_prev,shake_vir);
2952             copy_mat(state->fvir_prev,force_vir);
2953         }
2954         /*  ################## END TRAJECTORY OUTPUT ################ */
2955
2956         /* Determine the pressure:
2957          * always when we want exact averages in the energy file,
2958          * at ns steps when we have pressure coupling,
2959          * otherwise only at energy output steps (set below).
2960          */
2961
2962         bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2963         bCalcEnerPres = bNstEner;
2964
2965         /* Do we need global communication ? */
2966         bGStat = (bGStatEveryStep || bStopCM || bNS ||
2967                   (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2968
2969         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2970
2971         if (do_ene || do_log)
2972         {
2973             bCalcEnerPres = TRUE;
2974             bGStat        = TRUE;
2975         }
2976
2977         /* Determine the wallclock run time up till now */
2978         run_time = gmx_gettime() - (double)runtime->real;
2979
2980         /* Check whether everything is still allright */
2981         if (((int)gmx_get_stop_condition() > handled_stop_condition)
2982 #ifdef GMX_THREADS
2983             && MASTER(cr)
2984 #endif
2985             )
2986         {
2987             /* this is just make gs.sig compatible with the hack
2988                of sending signals around by MPI_Reduce with together with
2989                other floats */
2990             if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2991                 gs.sig[eglsSTOPCOND]=1;
2992             if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2993                 gs.sig[eglsSTOPCOND]=-1;
2994             /* < 0 means stop at next step, > 0 means stop at next NS step */
2995             if (fplog)
2996             {
2997                 fprintf(fplog,
2998                         "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2999                         gmx_get_signal_name(),
3000                         gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
3001                 fflush(fplog);
3002             }
3003             fprintf(stderr,
3004                     "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
3005                     gmx_get_signal_name(),
3006                     gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
3007             fflush(stderr);
3008             handled_stop_condition=(int)gmx_get_stop_condition();
3009         }
3010         else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
3011                  (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
3012                  gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
3013         {
3014             /* Signal to terminate the run */
3015             gs.sig[eglsSTOPCOND] = 1;
3016             if (fplog)
3017             {
3018                 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
3019             }
3020             fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
3021         }
3022
3023         if (bResetCountersHalfMaxH && MASTER(cr) &&
3024             run_time > max_hours*60.0*60.0*0.495)
3025         {
3026             gs.sig[eglsRESETCOUNTERS] = 1;
3027         }
3028
3029         if (ir->nstlist == -1 && !bRerunMD)
3030         {
3031             /* When bGStatEveryStep=FALSE, global_stat is only called
3032              * when we check the atom displacements, not at NS steps.
3033              * This means that also the bonded interaction count check is not
3034              * performed immediately after NS. Therefore a few MD steps could
3035              * be performed with missing interactions.
3036              * But wrong energies are never written to file,
3037              * since energies are only written after global_stat
3038              * has been called.
3039              */
3040             if (step >= nlh.step_nscheck)
3041             {
3042                 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
3043                                                      nlh.scale_tot,state->x);
3044             }
3045             else
3046             {
3047                 /* This is not necessarily true,
3048                  * but step_nscheck is determined quite conservatively.
3049                  */
3050                 nlh.nabnsb = 0;
3051             }
3052         }
3053
3054         /* In parallel we only have to check for checkpointing in steps
3055          * where we do global communication,
3056          *  otherwise the other nodes don't know.
3057          */
3058         if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
3059                            cpt_period >= 0 &&
3060                            (cpt_period == 0 ||
3061                             run_time >= nchkpt*cpt_period*60.0)) &&
3062             gs.set[eglsCHKPT] == 0)
3063         {
3064             gs.sig[eglsCHKPT] = 1;
3065         }
3066
3067         if (bIterations)
3068         {
3069             gmx_iterate_init(&iterate,bIterations);
3070         }
3071
3072         /* for iterations, we save these vectors, as we will be redoing the calculations */
3073         if (bIterations && iterate.bIterate)
3074         {
3075             copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
3076         }
3077         bFirstIterate = TRUE;
3078         while (bFirstIterate || (bIterations && iterate.bIterate))
3079         {
3080             /* We now restore these vectors to redo the calculation with improved extended variables */
3081             if (bIterations)
3082             {
3083                 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
3084             }
3085
3086             /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
3087                so scroll down for that logic */
3088
3089             /* #########   START SECOND UPDATE STEP ################# */
3090             GMX_MPE_LOG(ev_update_start);
3091             bOK = TRUE;
3092             if (!bRerunMD || rerun_fr.bV || bForceUpdate)
3093             {
3094                 wallcycle_start(wcycle,ewcUPDATE);
3095                 dvdl = 0;
3096                 /* Box is changed in update() when we do pressure coupling,
3097                  * but we should still use the old box for energy corrections and when
3098                  * writing it to the energy file, so it matches the trajectory files for
3099                  * the same timestep above. Make a copy in a separate array.
3100                  */
3101                 copy_mat(state->box,lastbox);
3102                 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
3103                 if (bTrotter)
3104                 {
3105                     if (bIterations && iterate.bIterate)
3106                     {
3107                         if (bFirstIterate)
3108                         {
3109                             scalevir = 1;
3110                         }
3111                         else
3112                         {
3113                             /* we use a new value of scalevir to converge the iterations faster */
3114                             scalevir = tracevir/trace(shake_vir);
3115                         }
3116                         msmul(shake_vir,scalevir,shake_vir);
3117                         m_add(force_vir,shake_vir,total_vir);
3118                         clear_mat(shake_vir);
3119                     }
3120                     trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ3);
3121                 }
3122                 /* We can only do Berendsen coupling after we have summed
3123                  * the kinetic energy or virial. Since the happens
3124                  * in global_state after update, we should only do it at
3125                  * step % nstlist = 1 with bGStatEveryStep=FALSE.
3126                  */
3127
3128                 if (ir->eI != eiVVAK)
3129                 {
3130                   update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
3131                 }
3132                 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
3133                                 upd,bInitStep);
3134
3135                 if (bVV)
3136                 {
3137                     /* velocity half-step update */
3138                     update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3139                                   ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,cr,nrnb,constr,&top->idef);
3140                 }
3141
3142                 /* Above, initialize just copies ekinh into ekin,
3143                  * it doesn't copy position (for VV),
3144                  * and entire integrator for MD.
3145                  */
3146
3147                 if (ir->eI==eiVVAK)
3148                 {
3149                     copy_rvecn(state->x,cbuf,0,state->natoms);
3150                 }
3151
3152                 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3153                               ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
3154                 wallcycle_stop(wcycle,ewcUPDATE);
3155
3156                 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
3157                                    &top->idef,shake_vir,force_vir,
3158                                    cr,nrnb,wcycle,upd,constr,
3159                                    bInitStep,FALSE,bCalcEnerPres,state->veta);
3160
3161                 if (ir->eI==eiVVAK)
3162                 {
3163                     /* erase F_EKIN and F_TEMP here? */
3164                     /* just compute the kinetic energy at the half step to perform a trotter step */
3165                     compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
3166                                     wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
3167                                     constr,NULL,FALSE,lastbox,
3168                                     top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
3169                                     cglo_flags | CGLO_TEMPERATURE | CGLO_CONSTRAINT
3170                         );
3171                     wallcycle_start(wcycle,ewcUPDATE);
3172                     trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ4);
3173                     /* now we know the scaling, we can compute the positions again again */
3174                     copy_rvecn(cbuf,state->x,0,state->natoms);
3175
3176                     update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3177                                   ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
3178                     wallcycle_stop(wcycle,ewcUPDATE);
3179
3180                     /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
3181                     /* are the small terms in the shake_vir here due
3182                      * to numerical errors, or are they important
3183                      * physically? I'm thinking they are just errors, but not completely sure.
3184                      * For now, will call without actually constraining, constr=NULL*/
3185                     update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
3186                                        &top->idef,tmp_vir,force_vir,
3187                                        cr,nrnb,wcycle,upd,NULL,
3188                                        bInitStep,FALSE,bCalcEnerPres,state->veta);
3189                 }
3190                 if (!bOK && !bFFscan)
3191                 {
3192                     gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
3193                 }
3194
3195                 if (fr->bSepDVDL && fplog && do_log)
3196                 {
3197                     fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
3198                 }
3199                 enerd->term[F_DHDL_CON] += dvdl;
3200             }
3201             else if (graph)
3202             {
3203                 /* Need to unshift here */
3204                 unshift_self(graph,state->box,state->x);
3205             }
3206
3207             GMX_BARRIER(cr->mpi_comm_mygroup);
3208             GMX_MPE_LOG(ev_update_finish);
3209
3210             if (vsite != NULL)
3211             {
3212                 wallcycle_start(wcycle,ewcVSITECONSTR);
3213                 if (graph != NULL)
3214                 {
3215                     shift_self(graph,state->box,state->x);
3216                 }
3217                 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
3218                                  top->idef.iparams,top->idef.il,
3219                                  fr->ePBC,fr->bMolPBC,graph,cr,state->box);
3220
3221                 if (graph != NULL)
3222                 {
3223                     unshift_self(graph,state->box,state->x);
3224                 }
3225                 wallcycle_stop(wcycle,ewcVSITECONSTR);
3226             }
3227
3228             /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
3229             if (ir->nstlist == -1 && bFirstIterate)
3230             {
3231                 gs.sig[eglsNABNSB] = nlh.nabnsb;
3232             }
3233             compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
3234                             wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
3235                             constr,
3236                             bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
3237                             lastbox,
3238                             top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
3239                             cglo_flags
3240                             | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
3241                             | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
3242                             | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
3243                             | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
3244                             | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
3245                             | CGLO_CONSTRAINT
3246                 );
3247             if (ir->nstlist == -1 && bFirstIterate)
3248             {
3249                 nlh.nabnsb = gs.set[eglsNABNSB];
3250                 gs.set[eglsNABNSB] = 0;
3251             }
3252             /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
3253             /* #############  END CALC EKIN AND PRESSURE ################# */
3254
3255             /* Note: this is OK, but there are some numerical precision issues with using the convergence of
3256                the virial that should probably be addressed eventually. state->veta has better properies,
3257                but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
3258                generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
3259
3260             if (bIterations &&
3261                 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
3262                                trace(shake_vir),&tracevir))
3263             {
3264                 break;
3265             }
3266             bFirstIterate = FALSE;
3267         }
3268
3269         update_box(fplog,step,ir,mdatoms,state,graph,f,
3270                    ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
3271
3272         /* ################# END UPDATE STEP 2 ################# */
3273         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
3274
3275         /* The coordinates (x) were unshifted in update */
3276 /*        if (bFFscan && (shellfc==NULL || bConverged))
3277         {
3278             if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
3279                                  f,NULL,xcopy,
3280                                  &(top_global->mols),mdatoms->massT,pres))
3281             {
3282                 if (gmx_parallel_env_initialized())
3283                 {
3284                     gmx_finalize();
3285                 }
3286                 fprintf(stderr,"\n");
3287                 exit(0);
3288             }
3289         }*/
3290         if (!bGStat)
3291         {
3292             /* We will not sum ekinh_old,
3293              * so signal that we still have to do it.
3294              */
3295             bSumEkinhOld = TRUE;
3296         }
3297
3298 /*        if (bTCR)
3299         {*/
3300             /* Only do GCT when the relaxation of shells (minimization) has converged,
3301              * otherwise we might be coupling to bogus energies.
3302              * In parallel we must always do this, because the other sims might
3303              * update the FF.
3304              */
3305
3306             /* Since this is called with the new coordinates state->x, I assume
3307              * we want the new box state->box too. / EL 20040121
3308              */
3309 /*            do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
3310                         ir,MASTER(cr),
3311                         mdatoms,&(top->idef),mu_aver,
3312                         top_global->mols.nr,cr,
3313                         state->box,total_vir,pres,
3314                         mu_tot,state->x,f,bConverged);
3315             debug_gmx();
3316         }*/
3317
3318         /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
3319
3320         sum_dhdl(enerd,state->lambda,ir);
3321         /* use the directly determined last velocity, not actually the averaged half steps */
3322         if (bTrotter && ir->eI==eiVV)
3323         {
3324             enerd->term[F_EKIN] = last_ekin;
3325         }
3326         enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
3327
3328         switch (ir->etc)
3329         {
3330         case etcNO:
3331             break;
3332         case etcBERENDSEN:
3333             break;
3334         case etcNOSEHOOVER:
3335             if (IR_NVT_TROTTER(ir)) {
3336                 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + last_conserved;
3337             } else {
3338                 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] +
3339                     NPT_energy(ir,state,&MassQ);
3340             }
3341             break;
3342         case etcVRESCALE:
3343             enerd->term[F_ECONSERVED] =
3344                 enerd->term[F_ETOT] + vrescale_energy(&(ir->opts),
3345                                                       state->therm_integral);
3346             break;
3347         default:
3348             break;
3349         }
3350
3351         /* Check for excessively large energies */
3352 /*        if (bIonize)
3353         {
3354 #ifdef GMX_DOUBLE
3355             real etot_max = 1e200;
3356 #else
3357             real etot_max = 1e30;
3358 #endif
3359             if (fabs(enerd->term[F_ETOT]) > etot_max)
3360             {
3361                 fprintf(stderr,"Energy too large (%g), giving up\n",
3362                         enerd->term[F_ETOT]);
3363             }
3364         }*/
3365         /* #########  END PREPARING EDR OUTPUT  ###########  */
3366
3367         /* Time for performance */
3368         if (((step % stepout) == 0) || bLastStep)
3369         {
3370             runtime_upd_proc(runtime);
3371         }
3372
3373         /* Output stuff */
3374         if (MASTER(cr))
3375         {
3376             gmx_bool do_dr,do_or;
3377
3378             if (!(bStartingFromCpt && (EI_VV(ir->eI))))
3379             {
3380                 if (bNstEner)
3381                 {
3382                     upd_mdebin(mdebin,bDoDHDL,TRUE,
3383                                t,mdatoms->tmass,enerd,state,lastbox,
3384                                shake_vir,force_vir,total_vir,pres,
3385                                ekind,mu_tot,constr);
3386                 }
3387                 else
3388                 {
3389                     upd_mdebin_step(mdebin);
3390                 }
3391
3392                 do_dr  = do_per_step(step,ir->nstdisreout);
3393                 do_or  = do_per_step(step,ir->nstorireout);
3394
3395                 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
3396                            step,t,
3397                            eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
3398             }
3399             if (ir->ePull != epullNO)
3400             {
3401                 pull_print_output(ir->pull,step,t);
3402             }
3403
3404             if (do_per_step(step,ir->nstlog))
3405             {
3406                 if(fflush(fplog) != 0)
3407                 {
3408                     gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
3409                 }
3410             }
3411         }
3412
3413
3414         /* Remaining runtime */
3415         if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
3416         {
3417             if (shellfc)
3418             {
3419                 fprintf(stderr,"\n");
3420             }
3421             print_time(stderr,runtime,step,ir,cr);
3422         }
3423
3424                 /* Set new positions for the group to embed */
3425                 if(!bLastStep){
3426                         if(step_rel<=it_xy)
3427                         {
3428                                 fac[0]+=xy_step;
3429                                 fac[1]+=xy_step;
3430                         } else if (step_rel<=(it_xy+it_z))
3431                         {
3432                                 fac[2]+=z_step;
3433                         }
3434                         resize(ins_at,r_ins,state_global->x,pos_ins,fac);
3435                 }
3436
3437         /* Replica exchange */
3438 /*        bExchanged = FALSE;
3439         if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
3440             do_per_step(step,repl_ex_nst))
3441         {
3442             bExchanged = replica_exchange(fplog,cr,repl_ex,
3443                                           state_global,enerd->term,
3444                                           state,step,t);
3445         }
3446         if (bExchanged && PAR(cr))
3447         {
3448             if (DOMAINDECOMP(cr))
3449             {
3450                 dd_partition_system(fplog,step,cr,TRUE,1,
3451                                     state_global,top_global,ir,
3452                                     state,&f,mdatoms,top,fr,
3453                                     vsite,shellfc,constr,
3454                                     nrnb,wcycle,FALSE);
3455             }
3456             else
3457             {
3458                 bcast_state(cr,state,FALSE);
3459             }
3460         }*/
3461
3462         bFirstStep = FALSE;
3463         bInitStep = FALSE;
3464         bStartingFromCpt = FALSE;
3465
3466         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
3467         /* With all integrators, except VV, we need to retain the pressure
3468          * at the current step for coupling at the next step.
3469          */
3470         if ((state->flags & (1<<estPRES_PREV)) &&
3471             (bGStatEveryStep ||
3472              (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
3473         {
3474             /* Store the pressure in t_state for pressure coupling
3475              * at the next MD step.
3476              */
3477             copy_mat(pres,state->pres_prev);
3478         }
3479
3480         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
3481
3482         if (bRerunMD)
3483         {
3484             /* read next frame from input trajectory */
3485             bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
3486         }
3487
3488         if (!bRerunMD || !rerun_fr.bStep)
3489         {
3490             /* increase the MD step number */
3491             step++;
3492             step_rel++;
3493         }
3494
3495         cycles = wallcycle_stop(wcycle,ewcSTEP);
3496         if (DOMAINDECOMP(cr) && wcycle)
3497         {
3498             dd_cycles_add(cr->dd,cycles,ddCyclStep);
3499         }
3500
3501         if (step_rel == wcycle_get_reset_counters(wcycle) ||
3502             gs.set[eglsRESETCOUNTERS] != 0)
3503         {
3504             /* Reset all the counters related to performance over the run */
3505             reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
3506             wcycle_set_reset_counters(wcycle,-1);
3507             bResetCountersHalfMaxH = FALSE;
3508             gs.set[eglsRESETCOUNTERS] = 0;
3509         }
3510     }
3511     /* End of main MD loop */
3512     debug_gmx();
3513     write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
3514                                         *top_global->name,top_global,
3515                                         state_global->x,state_global->v,
3516                                         ir->ePBC,state->box);
3517
3518     /* Stop the time */
3519     runtime_end(runtime);
3520
3521     if (bRerunMD)
3522     {
3523         close_trj(status);
3524     }
3525
3526     if (!(cr->duty & DUTY_PME))
3527     {
3528         /* Tell the PME only node to finish */
3529         gmx_pme_finish(cr);
3530     }
3531
3532     if (MASTER(cr))
3533     {
3534         if (ir->nstcalcenergy > 0 && !bRerunMD)
3535         {
3536             print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
3537                        eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
3538         }
3539     }
3540
3541     done_mdoutf(outf);
3542
3543     debug_gmx();
3544
3545     if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
3546     {
3547         fprintf(fplog,"Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n",nlh.s1/nlh.nns,sqrt(nlh.s2/nlh.nns - sqr(nlh.s1/nlh.nns)));
3548         fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
3549     }
3550
3551     if (shellfc && fplog)
3552     {
3553         fprintf(fplog,"Fraction of iterations that converged:           %.2f %%\n",
3554                 (nconverged*100.0)/step_rel);
3555         fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
3556                 tcount/step_rel);
3557     }
3558
3559 /*    if (repl_ex_nst > 0 && MASTER(cr))
3560     {
3561         print_replica_exchange_statistics(fplog,repl_ex);
3562     }*/
3563
3564     runtime->nsteps_done = step_rel;
3565
3566     return 0;
3567 }
3568
3569
3570 int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
3571              const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
3572              int nstglobalcomm,
3573              ivec ddxyz,int dd_node_order,real rdd,real rconstr,
3574              const char *dddlb_opt,real dlb_scale,
3575              const char *ddcsx,const char *ddcsy,const char *ddcsz,
3576              int nstepout,int resetstep,int nmultisim,int repl_ex_nst,int repl_ex_seed,
3577              real pforce,real cpt_period,real max_hours,
3578              const char *deviceOptions,
3579              unsigned long Flags,
3580              real xy_fac, real xy_max, real z_fac, real z_max,
3581              int it_xy, int it_z, real probe_rad, int low_up_rm,
3582              int pieces, gmx_bool bALLOW_ASYMMETRY, int maxwarn)
3583 {
3584     double     nodetime=0,realtime;
3585     t_inputrec *inputrec;
3586     t_state    *state=NULL;
3587     matrix     box;
3588     gmx_ddbox_t ddbox;
3589     int        npme_major,npme_minor;
3590     real       tmpr1,tmpr2;
3591     t_nrnb     *nrnb;
3592     gmx_mtop_t *mtop=NULL;
3593     t_mdatoms  *mdatoms=NULL;
3594     t_forcerec *fr=NULL;
3595     t_fcdata   *fcd=NULL;
3596     real       ewaldcoeff=0;
3597     gmx_pme_t  *pmedata=NULL;
3598     gmx_vsite_t *vsite=NULL;
3599     gmx_constr_t constr;
3600     int        i,m,nChargePerturbed=-1,status,nalloc;
3601     char       *gro;
3602     gmx_wallcycle_t wcycle;
3603     gmx_bool       bReadRNG,bReadEkin;
3604     int        list;
3605     gmx_runtime_t runtime;
3606     int        rc;
3607     gmx_large_int_t reset_counters;
3608     gmx_edsam_t ed=NULL;
3609     t_commrec   *cr_old=cr;
3610     int        nthreads=1,nthreads_requested=1;
3611
3612
3613         char                    *ins;
3614         int                     rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
3615         int                     ng,j,max_lip_rm,ins_grp_id,ins_nat,mem_nat,ntype,lip_rm,tpr_version;
3616         real                    xy_step=0,z_step=0;
3617         real                    prot_area;
3618         rvec                    *r_ins=NULL,fac;
3619         t_block                 *ins_at,*rest_at;
3620         pos_ins_t               *pos_ins;
3621         mem_t                   *mem_p;
3622         rm_t                    *rm_p;
3623         gmx_groups_t            *groups;
3624         gmx_bool                        bExcl=FALSE;
3625         t_atoms                 atoms;
3626         t_pbc                   *pbc;
3627         char                    **piecename=NULL;
3628
3629     /* CAUTION: threads may be started later on in this function, so
3630        cr doesn't reflect the final parallel state right now */
3631     snew(inputrec,1);
3632     snew(mtop,1);
3633
3634     if (bVerbose && SIMMASTER(cr))
3635     {
3636         fprintf(stderr,"Getting Loaded...\n");
3637     }
3638
3639     if (Flags & MD_APPENDFILES)
3640     {
3641         fplog = NULL;
3642     }
3643
3644     snew(state,1);
3645     if (MASTER(cr))
3646     {
3647         /* Read (nearly) all data required for the simulation */
3648         read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
3649
3650         /* NOW the threads will be started: */
3651 #ifdef GMX_THREADS
3652 #endif
3653     }
3654     /* END OF CAUTION: cr is now reliable */
3655
3656     if (PAR(cr))
3657     {
3658         /* now broadcast everything to the non-master nodes/threads: */
3659         init_parallel(fplog, cr, inputrec, mtop);
3660     }
3661     /* now make sure the state is initialized and propagated */
3662     set_state_entries(state,inputrec,cr->nnodes);
3663
3664     if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
3665     {
3666         /* All-vs-all loops do not work with domain decomposition */
3667         Flags |= MD_PARTDEC;
3668     }
3669
3670     if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
3671     {
3672         cr->npmenodes = 0;
3673     }
3674
3675         snew(ins_at,1);
3676         snew(pos_ins,1);
3677         if(MASTER(cr))
3678         {
3679                 tpr_version = get_tpr_version(ftp2fn(efTPX,nfile,fnm));
3680                 if (tpr_version<58)
3681                         gmx_fatal(FARGS,"Version of *.tpr file to old (%d). Rerun grompp with gromacs VERSION 4.0.3 or newer.\n",tpr_version);
3682
3683                 if( inputrec->eI != eiMD )
3684                         gmx_input("Change integrator to md in mdp file.");
3685
3686                 if(PAR(cr))
3687                         gmx_input("Sorry, parallel g_membed is not yet fully functrional.");
3688
3689                 groups=&(mtop->groups);
3690
3691                 atoms=gmx_mtop_global_atoms(mtop);
3692                 snew(mem_p,1);
3693                 fprintf(stderr,"\nSelect a group to embed in the membrane:\n");
3694                 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(ins_at->nr),&(ins_at->index),&ins);
3695                 ins_grp_id = search_string(ins,groups->ngrpname,(groups->grpname));
3696                 fprintf(stderr,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins);
3697                 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(mem_p->mem_at.nr),&(mem_p->mem_at.index),&(mem_p->name));
3698
3699                 pos_ins->pieces=pieces;
3700                 snew(pos_ins->nidx,pieces);
3701                 snew(pos_ins->subindex,pieces);
3702                 snew(piecename,pieces);
3703                 if (pieces>1)
3704                 {
3705                         fprintf(stderr,"\nSelect pieces to embed:\n");
3706                         get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),pieces,pos_ins->nidx,pos_ins->subindex,piecename);
3707                 }
3708                 else
3709                 {
3710                         /*use whole embedded group*/
3711                         snew(pos_ins->nidx,1);
3712                         snew(pos_ins->subindex,1);
3713                         pos_ins->nidx[0]=ins_at->nr;
3714                         pos_ins->subindex[0]=ins_at->index;
3715                 }
3716
3717                 if(probe_rad<0.2199999)
3718                 {
3719                         warn++;
3720                         fprintf(stderr,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 can result in overlap between waters "
3721                                         "and the group to embed, which will result in Lincs errors etc.\nIf you are sure, you can increase maxwarn.\n\n",warn);
3722                 }
3723
3724                 if(xy_fac<0.09999999)
3725                 {
3726                         warn++;
3727                         fprintf(stderr,"\nWarning %d:\nThe initial size of %s is probably too smal.\n"
3728                                         "If you are sure, you can increase maxwarn.\n\n",warn,ins);
3729                 }
3730
3731                 if(it_xy<1000)
3732                 {
3733                         warn++;
3734                         fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d) is probably too small.\n"
3735                                         "Increase -nxy or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_xy);
3736                 }
3737
3738                 if( (it_z<100) && ( z_fac<0.99999999 || z_fac>1.0000001) )
3739                 {
3740                         warn++;
3741                         fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d) is probably too small.\n"
3742                                        "Increase -nz or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_z);
3743                 }
3744
3745                 if(it_xy+it_z>inputrec->nsteps)
3746                 {
3747                         warn++;
3748                         fprintf(stderr,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the number of steps in the tpr.\n"
3749                                         "If you are sure, you can increase maxwarn.\n\n",warn);
3750                 }
3751
3752                 fr_id=-1;
3753                 if( inputrec->opts.ngfrz==1)
3754                         gmx_fatal(FARGS,"You did not specify \"%s\" as a freezegroup.",ins);
3755                 for(i=0;i<inputrec->opts.ngfrz;i++)
3756                 {
3757                         tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
3758                         if(ins_grp_id==tmp_id)
3759                         {
3760                                 fr_id=tmp_id;
3761                                 fr_i=i;
3762                         }
3763                 }
3764                 if (fr_id == -1 )
3765                         gmx_fatal(FARGS,"\"%s\" not as freezegroup defined in the mdp-file.",ins);
3766
3767                 for(i=0;i<DIM;i++)
3768                         if( inputrec->opts.nFreeze[fr_i][i] != 1)
3769                                 gmx_fatal(FARGS,"freeze dimensions for %s are not Y Y Y\n",ins);
3770
3771                 ng = groups->grps[egcENER].nr;
3772                 if (ng == 1)
3773                         gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
3774
3775                 for(i=0;i<ng;i++)
3776                 {
3777                         for(j=0;j<ng;j++)
3778                         {
3779                                 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
3780                                 {
3781                                         bExcl = TRUE;
3782                                         if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
3783                                                 gmx_fatal(FARGS,"Energy exclusions \"%s\" and  \"%s\" do not match the group to embed \"%s\"",
3784                                                                 *groups->grpname[groups->grps[egcENER].nm_ind[i]],
3785                                                                 *groups->grpname[groups->grps[egcENER].nm_ind[j]],ins);
3786                                 }
3787                         }
3788                 }
3789                 if (!bExcl)
3790                         gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in the freeze group");
3791
3792                 /* Set all atoms in box*/
3793                 /*set_inbox(state->natoms,state->x);*/
3794
3795                 /* Guess the area the protein will occupy in the membrane plane  Calculate area per lipid*/
3796                 snew(rest_at,1);
3797                 ins_nat = init_ins_at(ins_at,rest_at,state,pos_ins,groups,ins_grp_id,xy_max);
3798                 /* Check moleculetypes in insertion group */
3799                 check_types(ins_at,rest_at,mtop);
3800
3801                 mem_nat = init_mem_at(mem_p,mtop,state->x,state->box,pos_ins);
3802
3803                 prot_area = est_prot_area(pos_ins,state->x,ins_at,mem_p);
3804                 if ( (prot_area>7.5) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX])<50) )
3805                 {
3806                         warn++;
3807                         fprintf(stderr,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
3808                                         "This might cause pressure problems during the growth phase. Just try with\n"
3809                                         "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
3810                                         "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn);
3811                 }
3812                 if(warn>maxwarn)
3813                                         gmx_fatal(FARGS,"Too many warnings.\n");
3814
3815                 printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area);
3816                 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\nThe area per lipid is %.4f nm^2.\n",mem_p->nmol,mem_p->lip_area);
3817
3818                 /* Maximum number of lipids to be removed*/
3819                 max_lip_rm=(int)(2*prot_area/mem_p->lip_area);
3820                 printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm);
3821
3822                 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
3823                                 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
3824                                 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",xy_fac,z_fac,mem_p->zmin,mem_p->zmax);
3825
3826                 /* resize the protein by xy and by z if necessary*/
3827                 snew(r_ins,ins_at->nr);
3828                 init_resize(ins_at,r_ins,pos_ins,mem_p,state->x,bALLOW_ASYMMETRY);
3829                 fac[0]=fac[1]=xy_fac;
3830                 fac[2]=z_fac;
3831
3832                 xy_step =(xy_max-xy_fac)/(double)(it_xy);
3833                 z_step  =(z_max-z_fac)/(double)(it_z-1);
3834
3835                 resize(ins_at,r_ins,state->x,pos_ins,fac);
3836
3837                 /* remove overlapping lipids and water from the membrane box*/
3838                 /*mark molecules to be removed*/
3839                 snew(pbc,1);
3840                 set_pbc(pbc,inputrec->ePBC,state->box);
3841
3842                 snew(rm_p,1);
3843                 lip_rm = gen_rm_list(rm_p,ins_at,rest_at,pbc,mtop,state->x, r_ins, mem_p,pos_ins,probe_rad,low_up_rm,bALLOW_ASYMMETRY);
3844         lip_rm -= low_up_rm;
3845
3846                 if(fplog)
3847                         for(i=0;i<rm_p->nr;i++)
3848                                 fprintf(fplog,"rm mol %d\n",rm_p->mol[i]);
3849
3850                 for(i=0;i<mtop->nmolblock;i++)
3851                 {
3852                         ntype=0;
3853                         for(j=0;j<rm_p->nr;j++)
3854                                 if(rm_p->block[j]==i)
3855                                         ntype++;
3856                         printf("Will remove %d %s molecules\n",ntype,*(mtop->moltype[mtop->molblock[i].type].name));
3857                 }
3858
3859                 if(lip_rm>max_lip_rm)
3860                 {
3861                         warn++;
3862                         fprintf(stderr,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated protein area\n"
3863                                         "Try making the -xyinit resize factor smaller. If you are sure about this increase maxwarn.\n\n",warn);
3864                 }
3865
3866                 /*remove all lipids and waters overlapping and update all important structures*/
3867                 rm_group(inputrec,groups,mtop,rm_p,state,ins_at,pos_ins);
3868
3869                 rm_bonded_at = rm_bonded(ins_at,mtop);
3870                 if (rm_bonded_at != ins_at->nr)
3871                 {
3872                         fprintf(stderr,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
3873                                         "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
3874                                         "molecule type as atoms that are not to be embedded.\n",rm_bonded_at,ins_at->nr);
3875                 }
3876
3877                 if(warn>maxwarn)
3878                         gmx_fatal(FARGS,"Too many warnings.\nIf you are sure these warnings are harmless, you can increase -maxwarn");
3879
3880                 if (MASTER(cr))
3881                 {
3882                         if (ftp2bSet(efTOP,nfile,fnm))
3883                                 top_update(opt2fn("-p",nfile,fnm),ins,rm_p,mtop);
3884                 }
3885
3886                 sfree(pbc);
3887                 sfree(rest_at);
3888         }
3889
3890 #ifdef GMX_FAHCORE
3891     fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
3892 #endif
3893
3894     /* NMR restraints must be initialized before load_checkpoint,
3895      * since with time averaging the history is added to t_state.
3896      * For proper consistency check we therefore need to extend
3897      * t_state here.
3898      * So the PME-only nodes (if present) will also initialize
3899      * the distance restraints.
3900      */
3901     snew(fcd,1);
3902
3903     /* This needs to be called before read_checkpoint to extend the state */
3904     init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
3905
3906     if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
3907     {
3908         if (PAR(cr) && !(Flags & MD_PARTDEC))
3909         {
3910             gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
3911         }
3912         /* Orientation restraints */
3913         if (MASTER(cr))
3914         {
3915             init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
3916                         state);
3917         }
3918     }
3919
3920     if (DEFORM(*inputrec))
3921     {
3922         /* Store the deform reference box before reading the checkpoint */
3923         if (SIMMASTER(cr))
3924         {
3925             copy_mat(state->box,box);
3926         }
3927         if (PAR(cr))
3928         {
3929             gmx_bcast(sizeof(box),box,cr);
3930         }
3931         /* Because we do not have the update struct available yet
3932          * in which the reference values should be stored,
3933          * we store them temporarily in static variables.
3934          * This should be thread safe, since they are only written once
3935          * and with identical values.
3936          */
3937 /*        deform_init_init_step_tpx = inputrec->init_step;*/
3938 /*        copy_mat(box,deform_init_box_tpx);*/
3939     }
3940
3941     if (opt2bSet("-cpi",nfile,fnm))
3942     {
3943         /* Check if checkpoint file exists before doing continuation.
3944          * This way we can use identical input options for the first and subsequent runs...
3945          */
3946         if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
3947         {
3948             load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
3949                             cr,Flags & MD_PARTDEC,ddxyz,
3950                             inputrec,state,&bReadRNG,&bReadEkin,
3951                             (Flags & MD_APPENDFILES));
3952
3953             if (bReadRNG)
3954             {
3955                 Flags |= MD_READ_RNG;
3956             }
3957             if (bReadEkin)
3958             {
3959                 Flags |= MD_READ_EKIN;
3960             }
3961         }
3962     }
3963
3964     if ((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
3965     {
3966         gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
3967                              Flags,&fplog);
3968     }
3969
3970     if (SIMMASTER(cr))
3971     {
3972         copy_mat(state->box,box);
3973     }
3974
3975     if (PAR(cr))
3976     {
3977         gmx_bcast(sizeof(box),box,cr);
3978     }
3979
3980     if (bVerbose && SIMMASTER(cr))
3981     {
3982         fprintf(stderr,"Loaded with Money\n\n");
3983     }
3984
3985     if (PAR(cr) && !((Flags & MD_PARTDEC) || EI_TPI(inputrec->eI)))
3986     {
3987         cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
3988                                            dddlb_opt,dlb_scale,
3989                                            ddcsx,ddcsy,ddcsz,
3990                                            mtop,inputrec,
3991                                            box,state->x,
3992                                            &ddbox,&npme_major,&npme_minor);
3993
3994         make_dd_communicators(fplog,cr,dd_node_order);
3995
3996         /* Set overallocation to avoid frequent reallocation of arrays */
3997         set_over_alloc_dd(TRUE);
3998     }
3999     else
4000     {
4001         /* PME, if used, is done on all nodes with 1D decomposition */
4002         cr->npmenodes = 0;
4003         cr->duty = (DUTY_PP | DUTY_PME);
4004         npme_major = cr->nnodes;
4005         npme_minor = 1;
4006
4007         if (inputrec->ePBC == epbcSCREW)
4008         {
4009             gmx_fatal(FARGS,
4010                       "pbc=%s is only implemented with domain decomposition",
4011                       epbc_names[inputrec->ePBC]);
4012         }
4013     }
4014
4015     if (PAR(cr))
4016     {
4017         /* After possible communicator splitting in make_dd_communicators.
4018          * we can set up the intra/inter node communication.
4019          */
4020         gmx_setup_nodecomm(fplog,cr);
4021     }
4022
4023     wcycle = wallcycle_init(fplog,resetstep,cr);
4024     if (PAR(cr))
4025     {
4026         /* Master synchronizes its value of reset_counters with all nodes
4027          * including PME only nodes */
4028         reset_counters = wcycle_get_reset_counters(wcycle);
4029         gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
4030         wcycle_set_reset_counters(wcycle, reset_counters);
4031     }
4032
4033
4034     snew(nrnb,1);
4035     if (cr->duty & DUTY_PP)
4036     {
4037         /* For domain decomposition we allocate dynamically
4038          * in dd_partition_system.
4039          */
4040         if (DOMAINDECOMP(cr))
4041         {
4042             bcast_state_setup(cr,state);
4043         }
4044         else
4045         {
4046             if (PAR(cr))
4047             {
4048                 if (!MASTER(cr))
4049                 {
4050                     snew(state,1);
4051                 }
4052                 bcast_state(cr,state,TRUE);
4053             }
4054         }
4055
4056         /* Dihedral Restraints */
4057         if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
4058         {
4059             init_dihres(fplog,mtop,inputrec,fcd);
4060         }
4061
4062         /* Initiate forcerecord */
4063         fr = mk_forcerec();
4064         init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
4065                       opt2fn("-table",nfile,fnm),
4066                       opt2fn("-tablep",nfile,fnm),
4067                       opt2fn("-tableb",nfile,fnm),FALSE,pforce);
4068
4069         /* version for PCA_NOT_READ_NODE (see md.c) */
4070         /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
4071           "nofile","nofile","nofile",FALSE,pforce);
4072           */
4073         fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
4074
4075         /* Initialize QM-MM */
4076         if(fr->bQMMM)
4077         {
4078             init_QMMMrec(cr,box,mtop,inputrec,fr);
4079         }
4080
4081         /* Initialize the mdatoms structure.
4082          * mdatoms is not filled with atom data,
4083          * as this can not be done now with domain decomposition.
4084          */
4085         mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
4086
4087         /* Initialize the virtual site communication */
4088         vsite = init_vsite(mtop,cr);
4089
4090         calc_shifts(box,fr->shift_vec);
4091
4092         /* With periodic molecules the charge groups should be whole at start up
4093          * and the virtual sites should not be far from their proper positions.
4094          */
4095         if (!inputrec->bContinuation && MASTER(cr) &&
4096             !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
4097         {
4098             /* Make molecules whole at start of run */
4099             if (fr->ePBC != epbcNONE)
4100             {
4101                 do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
4102             }
4103             if (vsite)
4104             {
4105                 /* Correct initial vsite positions are required
4106                  * for the initial distribution in the domain decomposition
4107                  * and for the initial shell prediction.
4108                  */
4109                 construct_vsites_mtop(fplog,vsite,mtop,state->x);
4110             }
4111         }
4112
4113         /* Initiate PPPM if necessary */
4114         if (fr->eeltype == eelPPPM)
4115         {
4116             if (mdatoms->nChargePerturbed)
4117             {
4118                 gmx_fatal(FARGS,"Free energy with %s is not implemented",
4119                           eel_names[fr->eeltype]);
4120             }
4121             status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
4122                                    getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
4123             if (status != 0)
4124             {
4125                 gmx_fatal(FARGS,"Error %d initializing PPPM",status);
4126             }
4127         }
4128
4129         if (EEL_PME(fr->eeltype))
4130         {
4131             ewaldcoeff = fr->ewaldcoeff;
4132             pmedata = &fr->pmedata;
4133         }
4134         else
4135         {
4136             pmedata = NULL;
4137         }
4138     }
4139     else
4140     {
4141         /* This is a PME only node */
4142
4143         /* We don't need the state */
4144         done_state(state);
4145
4146         ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
4147         snew(pmedata,1);
4148     }
4149
4150     /* Initiate PME if necessary,
4151      * either on all nodes or on dedicated PME nodes only. */
4152     if (EEL_PME(inputrec->coulombtype))
4153     {
4154         if (mdatoms)
4155         {
4156             nChargePerturbed = mdatoms->nChargePerturbed;
4157         }
4158         if (cr->npmenodes > 0)
4159         {
4160             /* The PME only nodes need to know nChargePerturbed */
4161             gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
4162         }
4163         if (cr->duty & DUTY_PME)
4164         {
4165             status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
4166                                   mtop ? mtop->natoms : 0,nChargePerturbed,
4167                                   (Flags & MD_REPRODUCIBLE));
4168             if (status != 0)
4169             {
4170                 gmx_fatal(FARGS,"Error %d initializing PME",status);
4171             }
4172         }
4173     }
4174
4175
4176 /*    if (integrator[inputrec->eI].func == do_md
4177 #ifdef GMX_OPENMM
4178         ||
4179         integrator[inputrec->eI].func == do_md_openmm
4180 #endif
4181         )
4182     {*/
4183         /* Turn on signal handling on all nodes */
4184         /*
4185          * (A user signal from the PME nodes (if any)
4186          * is communicated to the PP nodes.
4187          */
4188         signal_handler_install();
4189 /*    }*/
4190
4191     if (cr->duty & DUTY_PP)
4192     {
4193         if (inputrec->ePull != epullNO)
4194         {
4195             /* Initialize pull code */
4196             init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
4197                       EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
4198         }
4199
4200         constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
4201
4202         if (DOMAINDECOMP(cr))
4203         {
4204             dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
4205                             Flags & MD_DDBONDCHECK,fr->cginfo_mb);
4206
4207             set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
4208
4209             setup_dd_grid(fplog,cr->dd);
4210         }
4211
4212         /* Now do whatever the user wants us to do (how flexible...) */
4213         do_md_membed(fplog,cr,nfile,fnm,
4214                                       oenv,bVerbose,bCompact,
4215                                       nstglobalcomm,
4216                                       vsite,constr,
4217                                       nstepout,inputrec,mtop,
4218                                       fcd,state,
4219                                       mdatoms,nrnb,wcycle,ed,fr,
4220                                       repl_ex_nst,repl_ex_seed,
4221                                       cpt_period,max_hours,
4222                                       deviceOptions,
4223                                       Flags,
4224                                       &runtime,
4225                                       fac, r_ins, pos_ins, ins_at,
4226                                       xy_step, z_step, it_xy, it_z);
4227
4228         if (inputrec->ePull != epullNO)
4229         {
4230             finish_pull(fplog,inputrec->pull);
4231         }
4232     }
4233     else
4234     {
4235         /* do PME only */
4236         gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
4237     }
4238
4239     if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
4240     {
4241         /* Some timing stats */
4242         if (MASTER(cr))
4243         {
4244             if (runtime.proc == 0)
4245             {
4246                 runtime.proc = runtime.real;
4247             }
4248         }
4249         else
4250         {
4251             runtime.real = 0;
4252         }
4253     }
4254
4255     wallcycle_stop(wcycle,ewcRUN);
4256
4257     /* Finish up, write some stuff
4258      * if rerunMD, don't write last frame again
4259      */
4260     finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
4261                inputrec,nrnb,wcycle,&runtime,
4262                EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
4263
4264     /* Does what it says */
4265     print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
4266
4267     /* Close logfile already here if we were appending to it */
4268     if (MASTER(cr) && (Flags & MD_APPENDFILES))
4269     {
4270         gmx_log_close(fplog);
4271     }
4272
4273     if (pieces>1)
4274     {
4275         sfree(piecename);
4276     }
4277
4278     rc=(int)gmx_get_stop_condition();
4279
4280     return rc;
4281 }
4282
4283 int gmx_membed(int argc,char *argv[])
4284 {
4285         const char *desc[] = {
4286                         "g_membed embeds a membrane protein into an equilibrated lipid bilayer at the position",
4287                         "and orientation specified by the user.\n",
4288                         "\n",
4289                         "SHORT MANUAL\n------------\n",
4290                         "The user should merge the structure files of the protein and membrane (+solvent), creating a",
4291                         "single structure file with the protein overlapping the membrane at the desired position and",
4292                         "orientation. Box size should be taken from the membrane structure file. The corresponding topology",
4293                         "files should also be merged. Consecutively, create a tpr file (input for g_membed) from these files,"
4294                         "with the following options included in the mdp file.\n",
4295                         " - integrator      = md\n",
4296                         " - energygrp       = Protein (or other group that you want to insert)\n",
4297                         " - freezegrps      = Protein\n",
4298                         " - freezedim       = Y Y Y\n",
4299                         " - energygrp_excl  = Protein Protein\n",
4300                         "The output is a structure file containing the protein embedded in the membrane. If a topology",
4301                         "file is provided, the number of lipid and ",
4302                         "solvent molecules will be updated to match the new structure file.\n",
4303                         "For a more extensive manual see Wolf et al, J Comp Chem 31 (2010) 2169-2174, Appendix.\n",
4304                         "\n",
4305                         "SHORT METHOD DESCRIPTION\n",
4306                         "------------------------\n",
4307                         "1. The protein is resized around its center of mass by a factor -xy in the xy-plane",
4308                         "(the membrane plane) and a factor -z in the z-direction (if the size of the",
4309                         "protein in the z-direction is the same or smaller than the width of the membrane, a",
4310                         "-z value larger than 1 can prevent that the protein will be enveloped by the lipids).\n",
4311                         "2. All lipid and solvent molecules overlapping with the resized protein are removed. All",
4312                         "intraprotein interactions are turned off to prevent numerical issues for small values of -xy",
4313                         " or -z\n",
4314                         "3. One md step is performed.\n",
4315                         "4. The resize factor (-xy or -z) is incremented by a small amount ((1-xy)/nxy or (1-z)/nz) and the",
4316                         "protein is resized again around its center of mass. The resize factor for the xy-plane",
4317                         "is incremented first. The resize factor for the z-direction is not changed until the -xy factor",
4318                         "is 1 (thus after -nxy iteration).\n",
4319                         "5. Repeat step 3 and 4 until the protein reaches its original size (-nxy + -nz iterations).\n",
4320                         "For a more extensive method descrition see Wolf et al, J Comp Chem, 31 (2010) 2169-2174.\n",
4321                         "\n",
4322                         "NOTE\n----\n",
4323                         " - Protein can be any molecule you want to insert in the membrane.\n",
4324                         " - It is recommended to perform a short equilibration run after the embedding",
4325                         "(see Wolf et al, J Comp Chem 31 (2010) 2169-2174, to re-equilibrate the membrane. Clearly",
4326                         "protein equilibration might require longer.\n",
4327                         "\n"
4328         };
4329         t_commrec    *cr;
4330         t_filenm fnm[] = {
4331                         { efTPX, "-f",      "into_mem", ffREAD },
4332                         { efNDX, "-n",      "index",    ffOPTRD },
4333                         { efTOP, "-p",      "topol",    ffOPTRW },
4334                         { efTRN, "-o",      NULL,       ffWRITE },
4335                         { efXTC, "-x",      NULL,       ffOPTWR },
4336                         { efCPT, "-cpi",    NULL,       ffOPTRD },
4337                         { efCPT, "-cpo",    NULL,       ffOPTWR },
4338                         { efSTO, "-c",      "membedded",  ffWRITE },
4339                         { efEDR, "-e",      "ener",     ffWRITE },
4340                         { efLOG, "-g",      "md",       ffWRITE },
4341                         { efEDI, "-ei",     "sam",      ffOPTRD },
4342                         { efTRX, "-rerun",  "rerun",    ffOPTRD },
4343                         { efXVG, "-table",  "table",    ffOPTRD },
4344                         { efXVG, "-tablep", "tablep",   ffOPTRD },
4345                         { efXVG, "-tableb", "table",    ffOPTRD },
4346                         { efXVG, "-dhdl",   "dhdl",     ffOPTWR },
4347                         { efXVG, "-field",  "field",    ffOPTWR },
4348                         { efXVG, "-table",  "table",    ffOPTRD },
4349                         { efXVG, "-tablep", "tablep",   ffOPTRD },
4350                         { efXVG, "-tableb", "table",    ffOPTRD },
4351                         { efTRX, "-rerun",  "rerun",    ffOPTRD },
4352                         { efXVG, "-tpi",    "tpi",      ffOPTWR },
4353                         { efXVG, "-tpid",   "tpidist",  ffOPTWR },
4354                         { efEDI, "-ei",     "sam",      ffOPTRD },
4355                         { efEDO, "-eo",     "sam",      ffOPTWR },
4356                         { efGCT, "-j",      "wham",     ffOPTRD },
4357                         { efGCT, "-jo",     "bam",      ffOPTWR },
4358                         { efXVG, "-ffout",  "gct",      ffOPTWR },
4359                         { efXVG, "-devout", "deviatie", ffOPTWR },
4360                         { efXVG, "-runav",  "runaver",  ffOPTWR },
4361                         { efXVG, "-px",     "pullx",    ffOPTWR },
4362                         { efXVG, "-pf",     "pullf",    ffOPTWR },
4363                         { efMTX, "-mtx",    "nm",       ffOPTWR },
4364                         { efNDX, "-dn",     "dipole",   ffOPTWR }
4365         };
4366 #define NFILE asize(fnm)
4367
4368         /* Command line options ! */
4369         gmx_bool bCart        = FALSE;
4370         gmx_bool bPPPME       = FALSE;
4371         gmx_bool bPartDec     = FALSE;
4372         gmx_bool bDDBondCheck = TRUE;
4373         gmx_bool bDDBondComm  = TRUE;
4374         gmx_bool bVerbose     = FALSE;
4375         gmx_bool bCompact     = TRUE;
4376         gmx_bool bSepPot      = FALSE;
4377         gmx_bool bRerunVSite  = FALSE;
4378         gmx_bool bIonize      = FALSE;
4379         gmx_bool bConfout     = TRUE;
4380         gmx_bool bReproducible = FALSE;
4381
4382         int  npme=-1;
4383         int  nmultisim=0;
4384         int  nstglobalcomm=-1;
4385         int  repl_ex_nst=0;
4386         int  repl_ex_seed=-1;
4387         int  nstepout=100;
4388         int  nthreads=0; /* set to determine # of threads automatically */
4389         int  resetstep=-1;
4390
4391         rvec realddxyz={0,0,0};
4392         const char *ddno_opt[ddnoNR+1] =
4393         { NULL, "interleave", "pp_pme", "cartesian", NULL };
4394         const char *dddlb_opt[] =
4395         { NULL, "auto", "no", "yes", NULL };
4396         real rdd=0.0,rconstr=0.0,dlb_scale=0.8,pforce=-1;
4397         char *ddcsx=NULL,*ddcsy=NULL,*ddcsz=NULL;
4398         real cpt_period=15.0,max_hours=-1;
4399         gmx_bool bAppendFiles=TRUE,bAddPart=TRUE;
4400         gmx_bool bResetCountersHalfWay=FALSE;
4401         output_env_t oenv=NULL;
4402         const char *deviceOptions = "";
4403
4404         real xy_fac = 0.5;
4405         real xy_max = 1.0;
4406         real z_fac = 1.0;
4407         real z_max = 1.0;
4408         int it_xy = 1000;
4409         int it_z = 0;
4410         real probe_rad = 0.22;
4411         int low_up_rm = 0;
4412         int maxwarn=0;
4413         int pieces=1;
4414     gmx_bool bALLOW_ASYMMETRY=FALSE;
4415
4416
4417 /* arguments relevant to OPENMM only*/
4418 #ifdef GMX_OPENMM
4419     gmx_input("g_membed not functional in openmm");
4420 #endif
4421
4422         t_pargs pa[] = {
4423                         { "-xyinit",   FALSE, etREAL,  {&xy_fac},       "Resize factor for the protein in the xy dimension before starting embedding" },
4424                         { "-xyend",   FALSE, etREAL,  {&xy_max},                "Final resize factor in the xy dimension" },
4425                         { "-zinit",    FALSE, etREAL,  {&z_fac},                "Resize factor for the protein in the z dimension before starting embedding" },
4426                         { "-zend",    FALSE, etREAL,  {&z_max},                 "Final resize faction in the z dimension" },
4427                         { "-nxy",     FALSE,  etINT,  {&it_xy},         "Number of iteration for the xy dimension" },
4428                         { "-nz",      FALSE,  etINT,  {&it_z},          "Number of iterations for the z dimension" },
4429                         { "-rad",     FALSE, etREAL,  {&probe_rad},     "Probe radius to check for overlap between the group to embed and the membrane"},
4430                         { "-pieces",  FALSE,  etINT,  {&pieces},        "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
4431             { "-asymmetry",FALSE, etBOOL,{&bALLOW_ASYMMETRY}, "Allow asymmetric insertion, i.e. the number of lipids removed from the upper and lower leaflet will not be checked." },
4432             { "-ndiff" ,  FALSE, etINT, {&low_up_rm},       "Number of lipids that will additionally be removed from the lower (negative number) or upper (positive number) membrane leaflet." },
4433                         { "-maxwarn", FALSE, etINT, {&maxwarn},                 "Maximum number of warning allowed" },
4434   { "-pd",      FALSE, etBOOL,{&bPartDec},
4435     "HIDDENUse particle decompostion" },
4436   { "-dd",      FALSE, etRVEC,{&realddxyz},
4437     "HIDDENDomain decomposition grid, 0 is optimize" },
4438   { "-nt",      FALSE, etINT, {&nthreads},
4439     "HIDDENNumber of threads to start (0 is guess)" },
4440   { "-npme",    FALSE, etINT, {&npme},
4441     "HIDDENNumber of separate nodes to be used for PME, -1 is guess" },
4442   { "-ddorder", FALSE, etENUM, {ddno_opt},
4443     "HIDDENDD node order" },
4444   { "-ddcheck", FALSE, etBOOL, {&bDDBondCheck},
4445     "HIDDENCheck for all bonded interactions with DD" },
4446   { "-ddbondcomm", FALSE, etBOOL, {&bDDBondComm},
4447     "HIDDENUse special bonded atom communication when -rdd > cut-off" },
4448   { "-rdd",     FALSE, etREAL, {&rdd},
4449     "HIDDENThe maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates" },
4450   { "-rcon",    FALSE, etREAL, {&rconstr},
4451     "HIDDENMaximum distance for P-LINCS (nm), 0 is estimate" },
4452   { "-dlb",     FALSE, etENUM, {dddlb_opt},
4453     "HIDDENDynamic load balancing (with DD)" },
4454   { "-dds",     FALSE, etREAL, {&dlb_scale},
4455     "HIDDENMinimum allowed dlb scaling of the DD cell size" },
4456   { "-ddcsx",   FALSE, etSTR, {&ddcsx},
4457     "HIDDENThe DD cell sizes in x" },
4458   { "-ddcsy",   FALSE, etSTR, {&ddcsy},
4459     "HIDDENThe DD cell sizes in y" },
4460   { "-ddcsz",   FALSE, etSTR, {&ddcsz},
4461     "HIDDENThe DD cell sizes in z" },
4462   { "-gcom",    FALSE, etINT,{&nstglobalcomm},
4463     "HIDDENGlobal communication frequency" },
4464   { "-compact", FALSE, etBOOL,{&bCompact},
4465     "Write a compact log file" },
4466   { "-seppot",  FALSE, etBOOL, {&bSepPot},
4467     "HIDDENWrite separate V and dVdl terms for each interaction type and node to the log file(s)" },
4468   { "-pforce",  FALSE, etREAL, {&pforce},
4469     "HIDDENPrint all forces larger than this (kJ/mol nm)" },
4470   { "-reprod",  FALSE, etBOOL,{&bReproducible},
4471     "HIDDENTry to avoid optimizations that affect binary reproducibility" },
4472   { "-multi",   FALSE, etINT,{&nmultisim},
4473     "HIDDENDo multiple simulations in parallel" },
4474   { "-replex",  FALSE, etINT, {&repl_ex_nst},
4475     "HIDDENAttempt replica exchange every # steps" },
4476   { "-reseed",  FALSE, etINT, {&repl_ex_seed},
4477     "HIDDENSeed for replica exchange, -1 is generate a seed" },
4478   { "-rerunvsite", FALSE, etBOOL, {&bRerunVSite},
4479     "HIDDENRecalculate virtual site coordinates with -rerun" },
4480   { "-ionize",  FALSE, etBOOL,{&bIonize},
4481     "HIDDENDo a simulation including the effect of an X-Ray bombardment on your system" },
4482   { "-confout", TRUE, etBOOL, {&bConfout},
4483     "HIDDENWrite the last configuration with -c and force checkpointing at the last step" },
4484   { "-stepout", FALSE, etINT, {&nstepout},
4485     "HIDDENFrequency of writing the remaining runtime" },
4486   { "-resetstep", FALSE, etINT, {&resetstep},
4487     "HIDDENReset cycle counters after these many time steps" },
4488   { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
4489     "HIDDENReset the cycle counters after half the number of steps or halfway -maxh" },
4490   { "-v",       FALSE, etBOOL,{&bVerbose},
4491     "Be loud and noisy" },
4492   { "-maxh",   FALSE, etREAL, {&max_hours},
4493     "HIDDENTerminate after 0.99 times this time (hours)" },
4494   { "-cpt",     FALSE, etREAL, {&cpt_period},
4495     "HIDDENCheckpoint interval (minutes)" },
4496   { "-append",  FALSE, etBOOL, {&bAppendFiles},
4497     "HIDDENAppend to previous output files when continuing from checkpoint" },
4498   { "-addpart",  FALSE, etBOOL, {&bAddPart},
4499     "HIDDENAdd the simulation part number to all output files when continuing from checkpoint" },
4500         };
4501         gmx_edsam_t  ed;
4502         unsigned long Flags, PCA_Flags;
4503         ivec     ddxyz;
4504         int      dd_node_order;
4505         gmx_bool     HaveCheckpoint;
4506         FILE     *fplog,*fptest;
4507         int      sim_part,sim_part_fn;
4508         const char *part_suffix=".part";
4509         char     suffix[STRLEN];
4510         int      rc;
4511
4512
4513         cr = init_par(&argc,&argv);
4514
4515         PCA_Flags = (PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_CAN_SET_DEFFNM
4516                         | (MASTER(cr) ? 0 : PCA_QUIET));
4517
4518
4519         /* Comment this in to do fexist calls only on master
4520          * works not with rerun or tables at the moment
4521          * also comment out the version of init_forcerec in md.c
4522          * with NULL instead of opt2fn
4523          */
4524         /*
4525    if (!MASTER(cr))
4526    {
4527    PCA_Flags |= PCA_NOT_READ_NODE;
4528    }
4529          */
4530
4531         parse_common_args(&argc,argv,PCA_Flags, NFILE,fnm,asize(pa),pa,
4532                         asize(desc),desc,0,NULL, &oenv);
4533
4534         /* we set these early because they might be used in init_multisystem()
4535    Note that there is the potential for npme>nnodes until the number of
4536    threads is set later on, if there's thread parallelization. That shouldn't
4537    lead to problems. */
4538         dd_node_order = nenum(ddno_opt);
4539         cr->npmenodes = npme;
4540
4541 #ifdef GMX_THREADS
4542         /* now determine the number of threads automatically. The threads are
4543    only started at mdrunner_threads, though. */
4544         if (nthreads<1)
4545         {
4546                 nthreads=tMPI_Get_recommended_nthreads();
4547         }
4548 #else
4549         nthreads=1;
4550 #endif
4551
4552
4553         if (repl_ex_nst != 0 && nmultisim < 2)
4554                 gmx_fatal(FARGS,"Need at least two replicas for replica exchange (option -multi)");
4555
4556         if (nmultisim > 1) {
4557 #ifndef GMX_THREADS
4558                 init_multisystem(cr,nmultisim,NFILE,fnm,TRUE);
4559 #else
4560                 gmx_fatal(FARGS,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
4561 #endif
4562         }
4563
4564         /* Check if there is ANY checkpoint file available */
4565         sim_part    = 1;
4566         sim_part_fn = sim_part;
4567         if (opt2bSet("-cpi",NFILE,fnm))
4568         {
4569                 bAppendFiles =
4570                         read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE,fnm,cr),
4571                                                         &sim_part_fn,NULL,cr,
4572                                                         bAppendFiles,NFILE,fnm,
4573                                                         part_suffix,&bAddPart);
4574                 if (sim_part_fn==0 && MASTER(cr))
4575                 {
4576                         fprintf(stdout,"No previous checkpoint file present, assuming this is a new run.\n");
4577                 }
4578                 else
4579                 {
4580                         sim_part = sim_part_fn + 1;
4581                 }
4582         }
4583         else
4584         {
4585                 bAppendFiles = FALSE;
4586         }
4587
4588         if (!bAppendFiles)
4589         {
4590                 sim_part_fn = sim_part;
4591         }
4592
4593         if (bAddPart && sim_part_fn > 1)
4594         {
4595                 /* This is a continuation run, rename trajectory output files
4596        (except checkpoint files) */
4597                 /* create new part name first (zero-filled) */
4598                 sprintf(suffix,"%s%04d",part_suffix,sim_part_fn);
4599
4600                 add_suffix_to_output_names(fnm,NFILE,suffix);
4601                 fprintf(stdout,"Checkpoint file is from part %d, new output files will be suffixed '%s'.\n",sim_part-1,suffix);
4602         }
4603
4604         Flags = opt2bSet("-rerun",NFILE,fnm) ? MD_RERUN : 0;
4605         Flags = Flags | (bSepPot       ? MD_SEPPOT       : 0);
4606         Flags = Flags | (bIonize       ? MD_IONIZE       : 0);
4607         Flags = Flags | (bPartDec      ? MD_PARTDEC      : 0);
4608         Flags = Flags | (bDDBondCheck  ? MD_DDBONDCHECK  : 0);
4609         Flags = Flags | (bDDBondComm   ? MD_DDBONDCOMM   : 0);
4610         Flags = Flags | (bConfout      ? MD_CONFOUT      : 0);
4611         Flags = Flags | (bRerunVSite   ? MD_RERUN_VSITE  : 0);
4612         Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
4613         Flags = Flags | (bAppendFiles  ? MD_APPENDFILES  : 0);
4614         Flags = Flags | (sim_part>1    ? MD_STARTFROMCPT : 0);
4615         Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
4616
4617
4618         /* We postpone opening the log file if we are appending, so we can
4619    first truncate the old log file and append to the correct position
4620    there instead.  */
4621         if ((MASTER(cr) || bSepPot) && !bAppendFiles)
4622         {
4623                 gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,!bSepPot,Flags,&fplog);
4624                 CopyRight(fplog,argv[0]);
4625                 please_cite(fplog,"Hess2008b");
4626                 please_cite(fplog,"Spoel2005a");
4627                 please_cite(fplog,"Lindahl2001a");
4628                 please_cite(fplog,"Berendsen95a");
4629         }
4630         else
4631         {
4632                 fplog = NULL;
4633         }
4634
4635         ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
4636         ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
4637         ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
4638
4639         /* even if nthreads = 1, we still call this one */
4640
4641         rc = mdrunner_membed(fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
4642                         nstglobalcomm,
4643                         ddxyz, dd_node_order, rdd, rconstr, dddlb_opt[0], dlb_scale,
4644                         ddcsx, ddcsy, ddcsz, nstepout, resetstep, nmultisim, repl_ex_nst,
4645                         repl_ex_seed, pforce, cpt_period, max_hours, deviceOptions, Flags,
4646                         xy_fac,xy_max,z_fac,z_max,
4647                         it_xy,it_z,probe_rad,low_up_rm,
4648                         pieces,bALLOW_ASYMMETRY,maxwarn);
4649
4650         if (gmx_parallel_env_initialized())
4651                 gmx_finalize();
4652
4653         if (MULTIMASTER(cr)) {
4654                 thanx(stderr);
4655         }
4656
4657         /* Log file has to be closed in mdrunner if we are appending to it
4658    (fplog not set here) */
4659         fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
4660
4661         if (MASTER(cr) && !bAppendFiles)
4662         {
4663                 gmx_log_close(fplog);
4664         }
4665
4666         return rc;
4667 }