src/kernel/do_md_openmm.h

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