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