src/kernel/md_openmm.c

   1 #ifdef HAVE_CONFIG_H
   2 #include <config.h>
   3 #endif
   4
   5 #include <signal.h>
   6 #include <stdlib.h>
   7
   8 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
   9 /* _isnan() */
  10 #include <float.h>
  11 #endif
  12
  13 #include "typedefs.h"
  14 #include "smalloc.h"
  15 #include "sysstuff.h"
  16 #include "vec.h"
  17 #include "statutil.h"
  18 #include "vcm.h"
  19 #include "mdebin.h"
  20 #include "nrnb.h"
  21 #include "calcmu.h"
  22 #include "index.h"
  23 #include "vsite.h"
  24 #include "update.h"
  25 #include "ns.h"
  26 #include "trnio.h"
  27 #include "xtcio.h"
  28 #include "mdrun.h"
  29 #include "confio.h"
  30 #include "network.h"
  31 #include "pull.h"
  32 #include "xvgr.h"
  33 #include "physics.h"
  34 #include "names.h"
  35 #include "xmdrun.h"
  36 #include "ionize.h"
  37 #include "disre.h"
  38 #include "orires.h"
  39 #include "dihre.h"
  40 #include "pppm.h"
  41 #include "pme.h"
  42 #include "mdatoms.h"
  43 #include "repl_ex.h"
  44 #include "qmmm.h"
  45 #include "mpelogging.h"
  46 #include "domdec.h"
  47 #include "partdec.h"
  48 #include "topsort.h"
  49 #include "coulomb.h"
  50 #include "constr.h"
  51 #include "shellfc.h"
  52 #include "compute_io.h"
  53 #include "mvdata.h"
  54 #include "checkpoint.h"
  55 #include "mtop_util.h"
  56 #include "genborn.h"
  57 #include "string2.h"
  58
  59 #ifdef GMX_THREADS
  60 #include "tmpi.h"
  61 #endif
  62
  63 #ifdef GMX_OPENMM
  64 #include "openmm_wrapper.h"
  65 #endif
  66
  67 double do_md_openmm(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
  68              const output_env_t oenv, bool bVerbose,bool bCompact,
  69              int nstglobalcomm,
  70              gmx_vsite_t *vsite,gmx_constr_t constr,
  71              int stepout,t_inputrec *ir,
  72              gmx_mtop_t *top_global,
  73              t_fcdata *fcd,
  74              t_state *state_global,
  75              t_mdatoms *mdatoms,
  76              t_nrnb *nrnb,gmx_wallcycle_t wcycle,
  77              gmx_edsam_t ed,t_forcerec *fr,
  78              int repl_ex_nst,int repl_ex_seed,
  79              real cpt_period,real max_hours,
  80              const char *deviceOptions,
  81              unsigned long Flags,
  82              gmx_runtime_t *runtime)
  83 {
  84     int        fp_trn=0,fp_xtc=0;
  85     ener_file_t fp_ene=NULL;
  86     gmx_large_int_t step,step_rel;
  87     const char *fn_cpt;
  88     FILE       *fp_dhdl=NULL,*fp_field=NULL;
  89     double     run_time;
  90     double     t,t0,lam0;
  91     bool       bGStatEveryStep,bGStat,bNstEner,bCalcPres,bCalcEner;
  92     bool       bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
  93                bFirstStep,bStateFromTPX,bInitStep,bLastStep,
  94                bBornRadii,bStartingFromCpt;
  95     bool       bDoDHDL=FALSE;
  96     bool       bNEMD,do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
  97                bForceUpdate=FALSE,bX,bV,bF,bXTC,bCPT=FALSE;
  98     bool       bMasterState;
  99     int        force_flags,cglo_flags;
 100     tensor     force_vir,shake_vir,total_vir,tmp_vir,pres;
 101     int        i,m,status;
 102     rvec       mu_tot;
 103     t_vcm      *vcm;
 104     t_state    *bufstate=NULL;
 105     matrix     *scale_tot,pcoupl_mu,M,ebox;
 106 //    gmx_nlheur_t nlh;
 107     t_trxframe rerun_fr;
 108     gmx_repl_ex_t repl_ex=NULL;
 109     int        nchkpt=1;
 110     /* Booleans (disguised as a reals) to checkpoint and terminate mdrun */
 111     real       chkpt=0,terminate=0,terminate_now=0;
 112
 113     gmx_localtop_t *top;
 114     t_mdebin *mdebin=NULL;
 115     t_state    *state=NULL;
 116     rvec       *f_global=NULL;
 117     int        n_xtc=-1;
 118     rvec       *x_xtc=NULL;
 119     gmx_enerdata_t *enerd;
 120     rvec       *f=NULL;
 121     gmx_global_stat_t gstat;
 122     gmx_update_t upd=NULL;
 123     t_graph    *graph=NULL;
 124
 125     bool        bFFscan;
 126     gmx_groups_t *groups;
 127     gmx_ekindata_t *ekind, *ekind_save;
 128     gmx_shellfc_t shellfc;
 129     int         count,nconverged=0;
 130     real        timestep=0;
 131     double      tcount=0;
 132     bool        bIonize=FALSE;
 133     bool        bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
 134     bool        bAppend;
 135     bool        bResetCountersHalfMaxH=FALSE;
 136     bool        bVV,bIterations,bIterate,bFirstIterate,bTemp,bPres,bTrotter;
 137     real        temp0,mu_aver=0,dvdl;
 138     int         a0,a1,gnx=0,ii;
 139     atom_id     *grpindex=NULL;
 140     char        *grpname;
 141     t_coupl_rec *tcr=NULL;
 142     rvec        *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
 143     matrix      boxcopy={{0}},lastbox;
 144         tensor      tmpvir;
 145         real        fom,oldfom,veta_save,pcurr,scalevir,tracevir;
 146         real        vetanew = 0;
 147     double      cycles;
 148     real        reset_counters=0,reset_counters_now=0;
 149         real        last_conserved = 0;
 150     real        last_ekin = 0;
 151         int         iter_i;
 152         t_extmass   MassQ;
 153     int         **trotter_seq;
 154     char        sbuf[22],sbuf2[22];
 155     bool        bHandledSignal=FALSE;
 156 //    gmx_iterate_t iterate;
 157
 158     const char *ommOptions = NULL;
 159     void   *openmmData;
 160
 161     /* Check for special mdrun options */
 162     bRerunMD = (Flags & MD_RERUN);
 163     bIonize  = (Flags & MD_IONIZE);
 164     bFFscan  = (Flags & MD_FFSCAN);
 165     bAppend  = (Flags & MD_APPENDFILES);
 166     if (Flags & MD_RESETCOUNTERSHALFWAY)
 167     {
 168         if (ir->nsteps > 0)
 169         {
 170             /* Signal to reset the counters half the simulation steps. */
 171             wcycle_set_reset_counters(wcycle,ir->nsteps/2);
 172         }
 173         /* Signal to reset the counters halfway the simulation time. */
 174         bResetCountersHalfMaxH = (max_hours > 0);
 175     }
 176
 177     /* md-vv uses averaged full step velocities for T-control
 178        md-vv2 uses averaged half step velocities for T-control (but full step ekin for P control)
 179        md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
 180     bVV = EI_VV(ir->eI);
 181     if (bVV) /* to store the initial velocities while computing virial */
 182     {
 183         snew(cbuf,top_global->natoms);
 184     }
 185     /* all the iteratative cases - only if there are constraints */
 186     bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
 187     bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
 188
 189     if (bRerunMD)
 190     {
 191         /* Since we don't know if the frames read are related in any way,
 192          * rebuild the neighborlist at every step.
 193          */
 194         ir->nstlist       = 1;
 195         ir->nstcalcenergy = 1;
 196         nstglobalcomm     = 1;
 197     }
 198
 199     check_ir_old_tpx_versions(cr,fplog,ir,top_global);
 200
 201 //****    nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
 202     bGStatEveryStep = (nstglobalcomm == 1);
 203
 204     if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
 205     {
 206         fprintf(fplog,
 207                 "To reduce the energy communication with nstlist = -1\n"
 208                 "the neighbor list validity should not be checked at every step,\n"
 209                 "this means that exact integration is not guaranteed.\n"
 210                 "The neighbor list validity is checked after:\n"
 211                 "  <n.list life time> - 2*std.dev.(n.list life time)  steps.\n"
 212                 "In most cases this will result in exact integration.\n"
 213                 "This reduces the energy communication by a factor of 2 to 3.\n"
 214                 "If you want less energy communication, set nstlist > 3.\n\n");
 215     }
 216
 217     if (bRerunMD || bFFscan)
 218     {
 219         ir->nstxtcout = 0;
 220     }
 221     groups = &top_global->groups;
 222
 223     /* Initial values */
 224     init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
 225             nrnb,top_global,&upd,
 226             nfile,fnm,&fp_trn,&fp_xtc,&fp_ene,&fn_cpt,
 227             &fp_dhdl,&fp_field,&mdebin,
 228             force_vir,shake_vir,mu_tot,&bNEMD,&bSimAnn,&vcm,state_global,Flags);
 229
 230     clear_mat(total_vir);
 231     clear_mat(pres);
 232     /* Energy terms and groups */
 233     snew(enerd,1);
 234     init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
 235     if (DOMAINDECOMP(cr))
 236     {
 237         f = NULL;
 238     }
 239     else
 240     {
 241         snew(f,top_global->natoms);
 242     }
 243
 244     /* Kinetic energy data */
 245     snew(ekind,1);
 246     init_ekindata(fplog,top_global,&(ir->opts),ekind);
 247     /* needed for iteration of constraints */
 248     snew(ekind_save,1);
 249     init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
 250     /* Copy the cos acceleration to the groups struct */
 251     ekind->cosacc.cos_accel = ir->cos_accel;
 252
 253     gstat = global_stat_init(ir);
 254     debug_gmx();
 255
 256     /* Check for polarizable models and flexible constraints */
 257     shellfc = init_shell_flexcon(fplog,
 258                                  top_global,n_flexible_constraints(constr),
 259                                  (ir->bContinuation ||
 260                                   (DOMAINDECOMP(cr) && !MASTER(cr))) ?
 261                                  NULL : state_global->x);
 262
 263     if (DEFORM(*ir))
 264     {
 265 #ifdef GMX_THREADS
 266         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
 267 #endif
 268         set_deform_reference_box(upd,
 269                                  deform_init_init_step_tpx,
 270                                  deform_init_box_tpx);
 271 #ifdef GMX_THREADS
 272         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
 273 #endif
 274     }
 275
 276     {
 277         double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
 278         if ((io > 2000) && MASTER(cr))
 279             fprintf(stderr,
 280                     "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
 281                     io);
 282     }
 283
 284     if (DOMAINDECOMP(cr)) {
 285         top = dd_init_local_top(top_global);
 286
 287         snew(state,1);
 288         dd_init_local_state(cr->dd,state_global,state);
 289
 290         if (DDMASTER(cr->dd) && ir->nstfout) {
 291             snew(f_global,state_global->natoms);
 292         }
 293     } else {
 294         if (PAR(cr)) {
 295             /* Initialize the particle decomposition and split the topology */
 296             top = split_system(fplog,top_global,ir,cr);
 297
 298             pd_cg_range(cr,&fr->cg0,&fr->hcg);
 299             pd_at_range(cr,&a0,&a1);
 300         } else {
 301             top = gmx_mtop_generate_local_top(top_global,ir);
 302
 303             a0 = 0;
 304             a1 = top_global->natoms;
 305         }
 306
 307         state = partdec_init_local_state(cr,state_global);
 308         f_global = f;
 309
 310         atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
 311
 312         if (vsite) {
 313             set_vsite_top(vsite,top,mdatoms,cr);
 314         }
 315
 316         if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
 317             graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
 318         }
 319
 320         if (shellfc) {
 321             make_local_shells(cr,mdatoms,shellfc);
 322         }
 323
 324         if (ir->pull && PAR(cr)) {
 325             dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
 326         }
 327     }
 328
 329     if (DOMAINDECOMP(cr))
 330     {
 331         /* Distribute the charge groups over the nodes from the master node */
 332         dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
 333                             state_global,top_global,ir,
 334                             state,&f,mdatoms,top,fr,
 335                             vsite,shellfc,constr,
 336                             nrnb,wcycle,FALSE);
 337     }
 338
 339     /* If not DD, copy gb data */
 340     if(ir->implicit_solvent && !DOMAINDECOMP(cr))
 341     {
 342         make_local_gb(cr,fr->born,ir->gb_algorithm);
 343     }
 344
 345     update_mdatoms(mdatoms,state->lambda);
 346 //#ifdef GMX_OPENMM
 347
 348       if(deviceOptions[0]=='\0')
 349       {
 350           /* empty options, which should default to OpenMM in this build */
 351           ommOptions=deviceOptions;
 352       }
 353       else
 354       {
 355           if(gmx_strncasecmp(deviceOptions,"OpenMM",6)!=0)
 356           {
 357               gmx_fatal(FARGS, "This Gromacs version currently only works with OpenMM. Use -device \"OpenMM:<options>\"");
 358           }
 359           else
 360           {
 361               ommOptions=strchr(deviceOptions,':');
 362               if(NULL!=ommOptions)
 363               {
 364                   /* Increase the pointer to skip the colon */
 365                   ommOptions++;
 366             }
 367           }
 368       }
 369       openmmData = openmm_init(fplog, ommOptions, cr, ir, top_global, top, mdatoms, fr, state);
 370 //#endif
 371
 372     if (MASTER(cr))
 373     {
 374         /* Update mdebin with energy history if appending to output files */
 375         if ( Flags & MD_APPENDFILES )
 376         {
 377             restore_energyhistory_from_state(mdebin,&state_global->enerhist);
 378         }
 379         /* Set the initial energy history in state to zero by updating once */
 380         update_energyhistory(&state_global->enerhist,mdebin);
 381     }
 382
 383     if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
 384         /* Set the random state if we read a checkpoint file */
 385         set_stochd_state(upd,state);
 386     }
 387
 388     /* Initialize constraints */
 389     if (constr) {
 390         if (!DOMAINDECOMP(cr))
 391             set_constraints(constr,top,ir,mdatoms,cr);
 392     }
 393
 394     /* Check whether we have to GCT stuff */
 395     bTCR = ftp2bSet(efGCT,nfile,fnm);
 396     if (bTCR) {
 397         if (MASTER(cr)) {
 398             fprintf(stderr,"Will do General Coupling Theory!\n");
 399         }
 400         gnx = top_global->mols.nr;
 401         snew(grpindex,gnx);
 402         for(i=0; (i<gnx); i++) {
 403             grpindex[i] = i;
 404         }
 405     }
 406
 407     if (repl_ex_nst > 0 && MASTER(cr))
 408         repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
 409                                         repl_ex_nst,repl_ex_seed);
 410
 411     if (!ir->bContinuation && !bRerunMD)
 412     {
 413         if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
 414         {
 415             /* Set the velocities of frozen particles to zero */
 416             for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
 417             {
 418                 for(m=0; m<DIM; m++)
 419                 {
 420                     if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 421                     {
 422                         state->v[i][m] = 0;
 423                     }
 424                 }
 425             }
 426         }
 427
 428         if (constr)
 429         {
 430             /* Constrain the initial coordinates and velocities */
 431             do_constrain_first(fplog,constr,ir,mdatoms,state,f,
 432                                graph,cr,nrnb,fr,top,shake_vir);
 433         }
 434         if (vsite)
 435         {
 436             /* Construct the virtual sites for the initial configuration */
 437             construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
 438                              top->idef.iparams,top->idef.il,
 439                              fr->ePBC,fr->bMolPBC,graph,cr,state->box);
 440         }
 441     }
 442
 443     debug_gmx();
 444
 445     /* I'm assuming we need global communication the first time! MRS */
 446 //*************
 447 //  cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
 448 //                  | (bVV ? CGLO_PRESSURE:0)
 449 //                  | (bVV ? CGLO_CONSTRAINT:0)
 450 //                  | (bRerunMD ? CGLO_RERUNMD:0)
 451 //                  | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
 452 //
 453 //    bSumEkinhOld = FALSE;
 454 //    compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
 455 //                    wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
 456 //                    constr,NULL,NULL,NULL,NULL,NULL,
 457 //                    NULL,state->box,
 458 //                    top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
 459 //    if (ir->eI == eiVVAK) {
 460 //        /* a second call to get the half step temperature initialized as well */
 461 //        /* we do the same call as above, but turn the pressure off -- internally, this
 462 //           is recognized as a velocity verlet half-step kinetic energy calculation.
 463 //           This minimized excess variables, but perhaps loses some logic?*/
 464 //
 465 //        compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
 466 //                        wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
 467 //                        constr,NULL,NULL,NULL,NULL,NULL,NULL,state->box,
 468 //                        top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
 469 //                        cglo_flags &~ CGLO_PRESSURE);
 470 //    }
 471
 472     /* Calculate the initial half step temperature, and save the ekinh_old */
 473     if (!(Flags & MD_STARTFROMCPT))
 474     {
 475         for(i=0; (i<ir->opts.ngtc); i++)
 476         {
 477             copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
 478         }
 479     }
 480     temp0 = enerd->term[F_TEMP];
 481
 482     /* if using an iterative algorithm, we need to create a working directory for the state. */
 483     if (bIterations)
 484     {
 485             bufstate = init_bufstate(state);
 486     }
 487     if (bFFscan)
 488     {
 489         snew(xcopy,state->natoms);
 490         snew(vcopy,state->natoms);
 491         copy_rvecn(state->x,xcopy,0,state->natoms);
 492         copy_rvecn(state->v,vcopy,0,state->natoms);
 493         copy_mat(state->box,boxcopy);
 494     }
 495
 496     /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
 497        temperature control */
 498     trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
 499
 500     if (MASTER(cr))
 501     {
 502         if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
 503         {
 504             fprintf(fplog,
 505                     "RMS relative constraint deviation after constraining: %.2e\n",
 506                     constr_rmsd(constr,FALSE));
 507         }
 508         if (!EI_VV(ir->eI))
 509         {
 510             enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
 511                                          and there is no previous step */
 512         }
 513         fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
 514         if (bRerunMD)
 515         {
 516             fprintf(stderr,"starting md rerun '%s', reading coordinates from"
 517                     " input trajectory '%s'\n\n",
 518                     *(top_global->name),opt2fn("-rerun",nfile,fnm));
 519             if (bVerbose)
 520             {
 521                 fprintf(stderr,"Calculated time to finish depends on nsteps from "
 522                         "run input file,\nwhich may not correspond to the time "
 523                         "needed to process input trajectory.\n\n");
 524             }
 525         }
 526         else
 527         {
 528             char tbuf[20];
 529             fprintf(stderr,"starting mdrun '%s'\n",
 530                     *(top_global->name));
 531             if (ir->nsteps >= 0)
 532             {
 533                 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
 534             }
 535             else
 536             {
 537                 sprintf(tbuf,"%s","infinite");
 538             }
 539             if (ir->init_step > 0)
 540             {
 541                 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 542                         gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
 543                         gmx_step_str(ir->init_step,sbuf2),
 544                         ir->init_step*ir->delta_t);
 545             }
 546             else
 547             {
 548                 fprintf(stderr,"%s steps, %s ps.\n",
 549                         gmx_step_str(ir->nsteps,sbuf),tbuf);
 550             }
 551         }
 552         fprintf(fplog,"\n");
 553     }
 554
 555     /* Set and write start time */
 556     runtime_start(runtime);
 557     print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
 558     wallcycle_start(wcycle,ewcRUN);
 559     if (fplog)
 560         fprintf(fplog,"\n");
 561
 562     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
 563
 564     debug_gmx();
 565     /***********************************************************
 566      *
 567      *             Loop over MD steps
 568      *
 569      ************************************************************/
 570
 571     /* if rerunMD then read coordinates and velocities from input trajectory */
 572     if (bRerunMD)
 573     {
 574         if (getenv("GMX_FORCE_UPDATE"))
 575         {
 576             bForceUpdate = TRUE;
 577         }
 578
 579         bNotLastFrame = read_first_frame(oenv,&status,
 580                                          opt2fn("-rerun",nfile,fnm),
 581                                          &rerun_fr,TRX_NEED_X | TRX_READ_V);
 582         if (rerun_fr.natoms != top_global->natoms)
 583         {
 584             gmx_fatal(FARGS,
 585                       "Number of atoms in trajectory (%d) does not match the "
 586                       "run input file (%d)\n",
 587                       rerun_fr.natoms,top_global->natoms);
 588         }
 589         if (ir->ePBC != epbcNONE)
 590         {
 591             if (!rerun_fr.bBox)
 592             {
 593                 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);
 594             }
 595             if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
 596             {
 597                 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
 598             }
 599
 600             /* Set the shift vectors.
 601              * Necessary here when have a static box different from the tpr box.
 602              */
 603             calc_shifts(rerun_fr.box,fr->shift_vec);
 604         }
 605     }
 606
 607     /* loop over MD steps or if rerunMD to end of input trajectory */
 608     bFirstStep = TRUE;
 609     /* Skip the first Nose-Hoover integration when we get the state from tpx */
 610     bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
 611     bInitStep = bFirstStep && (bStateFromTPX || bVV);
 612     bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
 613     bLastStep = FALSE;
 614     bSumEkinhOld = FALSE;
 615     bExchanged = FALSE;
 616
 617     step = ir->init_step;
 618     step_rel = 0;
 619
 620 //    if (ir->nstlist == -1)
 621 //    {
 622 //        init_nlistheuristics(&nlh,bGStatEveryStep,step);
 623 //    }
 624
 625     bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
 626     while (!bLastStep || (bRerunMD && bNotLastFrame)) {
 627
 628         wallcycle_start(wcycle,ewcSTEP);
 629
 630         GMX_MPE_LOG(ev_timestep1);
 631
 632         /* Now we have the energies and forces corresponding to the
 633          * coordinates at time t. We must output all of this before
 634          * the update.
 635          * for RerunMD t is read from input trajectory
 636          */
 637         GMX_MPE_LOG(ev_output_start);
 638
 639         bX   = do_per_step(step,ir->nstxout) || (bLastStep && ir->nstxout);
 640         bV   = do_per_step(step,ir->nstvout) || (bLastStep && ir->nstvout);
 641         bF   = do_per_step(step,ir->nstfout) || (bLastStep && ir->nstfout);
 642         bXTC = do_per_step(step,ir->nstxtcout) || (bLastStep && ir->nstxtcout);
 643         do_ene = do_per_step(step,ir->nstenergy) || (bLastStep && ir->nstenergy);
 644
 645 //#ifdef GMX_OPENMM
 646
 647       if( bX || bXTC || bV ){
 648         wallcycle_start(wcycle,ewcTRAJ);
 649         openmm_copy_state(openmmData, state, &t, f, enerd, bX||bXTC, bV, 0, 0);
 650         wallcycle_stop(wcycle,ewcTRAJ);
 651       }
 652
 653       openmm_take_one_step(openmmData);
 654       bLastStep = (step_rel == ir->nsteps);
 655       if (bX || bV || bF || bXTC || do_ene) {
 656         wallcycle_start(wcycle,ewcTRAJ);
 657         if( bF || do_ene ){
 658            openmm_copy_state(openmmData, state, &t, f, enerd, 0, 0, bF, do_ene);
 659       }
 660       upd_mdebin(mdebin,fp_dhdl,bGStatEveryStep && !bRerunMD,
 661                  t,mdatoms->tmass,enerd,state,lastbox,
 662                  shake_vir,force_vir,total_vir,pres,
 663                  ekind,mu_tot,constr);
 664       print_ebin(fp_ene,do_ene,FALSE,FALSE,do_log?fplog:NULL,step,t,
 665                  eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts)); // XXX check this do_log stuff!
 666 //#else
 667
 668             write_traj(fplog,cr,fp_trn,bX,bV,bF,fp_xtc,bXTC,ir->xtcprec,
 669                        fn_cpt,bCPT,top_global,ir->eI,ir->simulation_part,
 670                        step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
 671             debug_gmx();
 672             if (bLastStep && step_rel == ir->nsteps &&
 673                 (Flags & MD_CONFOUT) && MASTER(cr) &&
 674                 !bRerunMD && !bFFscan)
 675             {
 676                 /* x and v have been collected in write_traj */
 677                 fprintf(stderr,"\nWriting final coordinates.\n");
 678                 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
 679                     DOMAINDECOMP(cr))
 680                 {
 681                     /* Make molecules whole only for confout writing */
 682                     do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
 683                 }
 684                 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
 685                                     *top_global->name,top_global,
 686                                     state_global->x,state_global->v,
 687                                     ir->ePBC,state->box);
 688                 debug_gmx();
 689             }
 690             wallcycle_stop(wcycle,ewcTRAJ);
 691         }
 692         GMX_MPE_LOG(ev_output_finish);
 693
 694         bFirstStep = FALSE;
 695         bInitStep = FALSE;
 696         step++;
 697         step_rel++;
 698     }
 699     /* End of main MD loop */
 700     debug_gmx();
 701
 702 //#ifdef GMX_OPENMM
 703     openmm_cleanup(fplog, openmmData);
 704 //#endif
 705
 706     /* Stop the time */
 707     runtime_end(runtime);
 708
 709     if (bRerunMD)
 710     {
 711         close_trj(status);
 712     }
 713
 714     if (!(cr->duty & DUTY_PME))
 715     {
 716         /* Tell the PME only node to finish */
 717         gmx_pme_finish(cr);
 718     }
 719
 720     if (MASTER(cr))
 721     {
 722         if (bGStatEveryStep && !bRerunMD)
 723         {
 724             print_ebin(fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
 725                        eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
 726         }
 727         close_enx(fp_ene);
 728         if (ir->nstxtcout)
 729         {
 730             close_xtc(fp_xtc);
 731         }
 732         close_trn(fp_trn);
 733         if (fp_dhdl)
 734         {
 735             gmx_fio_fclose(fp_dhdl);
 736         }
 737         if (fp_field)
 738         {
 739             gmx_fio_fclose(fp_field);
 740         }
 741     }
 742     debug_gmx();
 743
 744 //    if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
 745 //    {
 746 //        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)));
 747 //        fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
 748 //    }
 749
 750     if (shellfc && fplog)
 751     {
 752         fprintf(fplog,"Fraction of iterations that converged:           %.2f %%\n",
 753                 (nconverged*100.0)/step_rel);
 754         fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
 755                 tcount/step_rel);
 756     }
 757
 758     if (repl_ex_nst > 0 && MASTER(cr))
 759     {
 760         print_replica_exchange_statistics(fplog,repl_ex);
 761     }
 762
 763     runtime->nsteps_done = step_rel;
 764
 765     return 0;
 766 }
 767