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