src/tools/gmx_membed.c

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