src/kernel/md_openmm.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  * 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-2010, 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
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <signal.h>
  41 #include <stdlib.h>
  42
  43 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
  44 /* _isnan() */
  45 #include <float.h>
  46 #endif
  47
  48 #include "typedefs.h"
  49 #include "smalloc.h"
  50 #include "sysstuff.h"
  51 #include "vec.h"
  52 #include "statutil.h"
  53 #include "vcm.h"
  54 #include "mdebin.h"
  55 #include "nrnb.h"
  56 #include "calcmu.h"
  57 #include "index.h"
  58 #include "vsite.h"
  59 #include "update.h"
  60 #include "ns.h"
  61 #include "trnio.h"
  62 #include "xtcio.h"
  63 #include "mdrun.h"
  64 #include "confio.h"
  65 #include "network.h"
  66 #include "pull.h"
  67 #include "xvgr.h"
  68 #include "physics.h"
  69 #include "names.h"
  70 #include "xmdrun.h"
  71 #include "ionize.h"
  72 #include "disre.h"
  73 #include "orires.h"
  74 #include "dihre.h"
  75 #include "pppm.h"
  76 #include "pme.h"
  77 #include "mdatoms.h"
  78 #include "qmmm.h"
  79 #include "mpelogging.h"
  80 #include "domdec.h"
  81 #include "partdec.h"
  82 #include "topsort.h"
  83 #include "coulomb.h"
  84 #include "constr.h"
  85 #include "compute_io.h"
  86 #include "mvdata.h"
  87 #include "checkpoint.h"
  88 #include "mtop_util.h"
  89 #include "sighandler.h"
  90 #include "genborn.h"
  91 #include "string2.h"
  92
  93 #ifdef GMX_THREADS
  94 #include "tmpi.h"
  95 #endif
  96
  97 /* include even when OpenMM not used to force compilation of do_md_openmm */
  98 #include "openmm_wrapper.h"
  99
 100 /* simulation conditions to transmit */
 101 enum { eglsNABNSB, eglsCHKPT, eglsSTOPCOND, eglsRESETCOUNTERS, eglsNR };
 102
 103 typedef struct
 104 {
 105     int nstms;       /* The frequency for intersimulation communication */
 106     int sig[eglsNR]; /* The signal set by one process in do_md */
 107     int set[eglsNR]; /* The communicated signal, equal for all processes */
 108 } globsig_t;
 109
 110
 111 static int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
 112 {
 113     int  *buf;
 114     bool bPos,bEqual;
 115     int  s,d;
 116
 117     snew(buf,ms->nsim);
 118     buf[ms->sim] = n;
 119     gmx_sumi_sim(ms->nsim,buf,ms);
 120     bPos   = TRUE;
 121     bEqual = TRUE;
 122     for (s=0; s<ms->nsim; s++)
 123     {
 124         bPos   = bPos   && (buf[s] > 0);
 125         bEqual = bEqual && (buf[s] == buf[0]);
 126     }
 127     if (bPos)
 128     {
 129         if (bEqual)
 130         {
 131             nmin = min(nmin,buf[0]);
 132         }
 133         else
 134         {
 135             /* Find the least common multiple */
 136             for (d=2; d<nmin; d++)
 137             {
 138                 s = 0;
 139                 while (s < ms->nsim && d % buf[s] == 0)
 140                 {
 141                     s++;
 142                 }
 143                 if (s == ms->nsim)
 144                 {
 145                     /* We found the LCM and it is less than nmin */
 146                     nmin = d;
 147                     break;
 148                 }
 149             }
 150         }
 151     }
 152     sfree(buf);
 153
 154     return nmin;
 155 }
 156
 157 static int multisim_nstsimsync(const t_commrec *cr,
 158                                const t_inputrec *ir,int repl_ex_nst)
 159 {
 160     int nmin;
 161
 162     if (MASTER(cr))
 163     {
 164         nmin = INT_MAX;
 165         nmin = multisim_min(cr->ms,nmin,ir->nstlist);
 166         nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
 167         nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
 168         if (nmin == INT_MAX)
 169         {
 170             gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
 171         }
 172         /* Avoid inter-simulation communication at every (second) step */
 173         if (nmin <= 2)
 174         {
 175             nmin = 10;
 176         }
 177     }
 178
 179     gmx_bcast(sizeof(int),&nmin,cr);
 180
 181     return nmin;
 182 }
 183
 184 static void init_global_signals(globsig_t *gs,const t_commrec *cr,
 185                                 const t_inputrec *ir,int repl_ex_nst)
 186 {
 187     int i;
 188
 189     gs->nstms = 1;
 190
 191     for (i=0; i<eglsNR; i++)
 192     {
 193         gs->sig[i] = 0;
 194         gs->set[i] = 0;
 195     }
 196 }
 197
 198
 199 double do_md_openmm(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 200                     const output_env_t oenv, bool bVerbose,bool bCompact,
 201                     int nstglobalcomm,
 202                     gmx_vsite_t *vsite,gmx_constr_t constr,
 203                     int stepout,t_inputrec *ir,
 204                     gmx_mtop_t *top_global,
 205                     t_fcdata *fcd,
 206                     t_state *state_global,
 207                     t_mdatoms *mdatoms,
 208                     t_nrnb *nrnb,gmx_wallcycle_t wcycle,
 209                     gmx_edsam_t ed,t_forcerec *fr,
 210                     int repl_ex_nst,int repl_ex_seed,
 211                     real cpt_period,real max_hours,
 212                     const char *deviceOptions,
 213                     unsigned long Flags,
 214                     gmx_runtime_t *runtime)
 215 {
 216     gmx_mdoutf_t *outf;
 217     gmx_large_int_t step,step_rel;
 218     double     run_time;
 219     double     t,t0,lam0;
 220     bool       bSimAnn,
 221     bFirstStep,bStateFromTPX,bLastStep,bStartingFromCpt;
 222     bool       bInitStep=TRUE;
 223     bool       do_ene,do_log, do_verbose,
 224     bX,bV,bF,bCPT;
 225     tensor     force_vir,shake_vir,total_vir,pres;
 226     int        i,m;
 227     int        mdof_flags;
 228     rvec       mu_tot;
 229     t_vcm      *vcm;
 230     int        nchkpt=1;
 231     gmx_localtop_t *top;
 232     t_mdebin *mdebin=NULL;
 233     t_state    *state=NULL;
 234     rvec       *f_global=NULL;
 235     int        n_xtc=-1;
 236     rvec       *x_xtc=NULL;
 237     gmx_enerdata_t *enerd;
 238     rvec       *f=NULL;
 239     gmx_global_stat_t gstat;
 240     gmx_update_t upd=NULL;
 241     t_graph    *graph=NULL;
 242     globsig_t   gs;
 243
 244     gmx_groups_t *groups;
 245     gmx_ekindata_t *ekind, *ekind_save;
 246     bool        bAppend;
 247     int         a0,a1;
 248     matrix      lastbox;
 249     real        reset_counters=0,reset_counters_now=0;
 250     char        sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
 251     int         handled_stop_condition=gmx_stop_cond_none;
 252
 253     const char *ommOptions = NULL;
 254     void   *openmmData;
 255
 256     bAppend  = (Flags & MD_APPENDFILES);
 257     check_ir_old_tpx_versions(cr,fplog,ir,top_global);
 258
 259     groups = &top_global->groups;
 260
 261     /* Initial values */
 262     init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
 263             nrnb,top_global,&upd,
 264             nfile,fnm,&outf,&mdebin,
 265             force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
 266
 267     clear_mat(total_vir);
 268     clear_mat(pres);
 269     /* Energy terms and groups */
 270     snew(enerd,1);
 271     init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
 272     snew(f,top_global->natoms);
 273
 274     /* Kinetic energy data */
 275     snew(ekind,1);
 276     init_ekindata(fplog,top_global,&(ir->opts),ekind);
 277     /* needed for iteration of constraints */
 278     snew(ekind_save,1);
 279     init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
 280     /* Copy the cos acceleration to the groups struct */
 281     ekind->cosacc.cos_accel = ir->cos_accel;
 282
 283     gstat = global_stat_init(ir);
 284     debug_gmx();
 285
 286     {
 287         double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
 288         if ((io > 2000) && MASTER(cr))
 289             fprintf(stderr,
 290                     "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
 291                     io);
 292     }
 293
 294     top = gmx_mtop_generate_local_top(top_global,ir);
 295
 296     a0 = 0;
 297     a1 = top_global->natoms;
 298
 299     state = partdec_init_local_state(cr,state_global);
 300     f_global = f;
 301
 302     atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
 303
 304     if (vsite)
 305     {
 306         set_vsite_top(vsite,top,mdatoms,cr);
 307     }
 308
 309     if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
 310     {
 311         graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
 312     }
 313
 314     update_mdatoms(mdatoms,state->lambda);
 315
 316     if (deviceOptions[0]=='\0')
 317     {
 318         /* empty options, which should default to OpenMM in this build */
 319         ommOptions=deviceOptions;
 320     }
 321     else
 322     {
 323         if (gmx_strncasecmp(deviceOptions,"OpenMM",6)!=0)
 324         {
 325             gmx_fatal(FARGS, "This Gromacs version currently only works with OpenMM. Use -device \"OpenMM:<options>\"");
 326         }
 327         else
 328         {
 329             ommOptions=strchr(deviceOptions,':');
 330             if (NULL!=ommOptions)
 331             {
 332                 /* Increase the pointer to skip the colon */
 333                 ommOptions++;
 334             }
 335         }
 336     }
 337
 338     openmmData = openmm_init(fplog, ommOptions, ir, top_global, top, mdatoms, fr, state);
 339
 340     if (MASTER(cr))
 341     {
 342         /* Update mdebin with energy history if appending to output files */
 343         if ( Flags & MD_APPENDFILES )
 344         {
 345             restore_energyhistory_from_state(mdebin,&state_global->enerhist);
 346         }
 347         /* Set the initial energy history in state to zero by updating once */
 348         update_energyhistory(&state_global->enerhist,mdebin);
 349     }
 350
 351     if (constr)
 352     {
 353         set_constraints(constr,top,ir,mdatoms,cr);
 354     }
 355
 356     if (!ir->bContinuation)
 357     {
 358         if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
 359         {
 360             /* Set the velocities of frozen particles to zero */
 361             for (i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
 362             {
 363                 for (m=0; m<DIM; m++)
 364                 {
 365                     if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 366                     {
 367                         state->v[i][m] = 0;
 368                     }
 369                 }
 370             }
 371         }
 372
 373         if (constr)
 374         {
 375             /* Constrain the initial coordinates and velocities */
 376             do_constrain_first(fplog,constr,ir,mdatoms,state,f,
 377                                graph,cr,nrnb,fr,top,shake_vir);
 378         }
 379         if (vsite)
 380         {
 381             /* Construct the virtual sites for the initial configuration */
 382             construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
 383                              top->idef.iparams,top->idef.il,
 384                              fr->ePBC,fr->bMolPBC,graph,cr,state->box);
 385         }
 386     }
 387
 388     debug_gmx();
 389
 390     if (MASTER(cr))
 391     {
 392         char tbuf[20];
 393         fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
 394         fprintf(stderr,"starting mdrun '%s'\n",
 395                 *(top_global->name));
 396         if (ir->nsteps >= 0)
 397         {
 398             sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
 399         }
 400         else
 401         {
 402             sprintf(tbuf,"%s","infinite");
 403         }
 404         if (ir->init_step > 0)
 405         {
 406             fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 407                     gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
 408                     gmx_step_str(ir->init_step,sbuf2),
 409                     ir->init_step*ir->delta_t);
 410         }
 411         else
 412         {
 413             fprintf(stderr,"%s steps, %s ps.\n",
 414                     gmx_step_str(ir->nsteps,sbuf),tbuf);
 415         }
 416     }
 417
 418     fprintf(fplog,"\n");
 419
 420     /* Set and write start time */
 421     runtime_start(runtime);
 422     print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
 423     wallcycle_start(wcycle,ewcRUN);
 424     if (fplog)
 425         fprintf(fplog,"\n");
 426
 427     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
 428
 429     debug_gmx();
 430     /***********************************************************
 431      *
 432      *             Loop over MD steps
 433      *
 434      ************************************************************/
 435
 436     /* loop over MD steps or if rerunMD to end of input trajectory */
 437     bFirstStep = TRUE;
 438     /* Skip the first Nose-Hoover integration when we get the state from tpx */
 439     bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
 440     bInitStep = bFirstStep && bStateFromTPX;
 441     bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
 442     bLastStep = FALSE;
 443
 444     init_global_signals(&gs,cr,ir,repl_ex_nst);
 445
 446     step = ir->init_step;
 447     step_rel = 0;
 448
 449     while (!bLastStep)
 450     {
 451         wallcycle_start(wcycle,ewcSTEP);
 452
 453         GMX_MPE_LOG(ev_timestep1);
 454
 455         bLastStep = (step_rel == ir->nsteps);
 456         t = t0 + step*ir->delta_t;
 457
 458         if (gs.set[eglsSTOPCOND] != 0)
 459         {
 460             bLastStep = TRUE;
 461         }
 462
 463         do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
 464         do_verbose = bVerbose &&
 465                      (step % stepout == 0 || bFirstStep || bLastStep);
 466
 467         if (MASTER(cr) && do_log)
 468         {
 469             print_ebin_header(fplog,step,t,state->lambda);
 470         }
 471
 472         clear_mat(force_vir);
 473         GMX_MPE_LOG(ev_timestep2);
 474
 475         /* We write a checkpoint at this MD step when:
 476          * either when we signalled through gs (in OpenMM NS works different),
 477          * or at the last step (but not when we do not want confout),
 478          * but never at the first step.
 479          */
 480         bCPT = ((gs.set[eglsCHKPT] ||
 481                  (bLastStep && (Flags & MD_CONFOUT))) &&
 482                 step > ir->init_step );
 483         if (bCPT)
 484         {
 485             gs.set[eglsCHKPT] = 0;
 486         }
 487
 488         /* Now we have the energies and forces corresponding to the
 489          * coordinates at time t. We must output all of this before
 490          * the update.
 491          * for RerunMD t is read from input trajectory
 492          */
 493         GMX_MPE_LOG(ev_output_start);
 494
 495         mdof_flags = 0;
 496         if (do_per_step(step,ir->nstxout))
 497         {
 498             mdof_flags |= MDOF_X;
 499         }
 500         if (do_per_step(step,ir->nstvout))
 501         {
 502             mdof_flags |= MDOF_V;
 503         }
 504         if (do_per_step(step,ir->nstfout))
 505         {
 506             mdof_flags |= MDOF_F;
 507         }
 508         if (do_per_step(step,ir->nstxtcout))
 509         {
 510             mdof_flags |= MDOF_XTC;
 511         }
 512         if (bCPT)
 513         {
 514             mdof_flags |= MDOF_CPT;
 515         };
 516         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
 517
 518         if (mdof_flags != 0 || do_ene || do_log)
 519         {
 520             wallcycle_start(wcycle,ewcTRAJ);
 521             bF = (mdof_flags & MDOF_F);
 522             bX = (mdof_flags & (MDOF_X | MDOF_XTC | MDOF_CPT));
 523             bV = (mdof_flags & (MDOF_V | MDOF_CPT));
 524
 525             openmm_copy_state(openmmData, state, &t, f, enerd, bX, bV, bF, do_ene);
 526
 527             upd_mdebin(mdebin, NULL,TRUE,
 528                        t,mdatoms->tmass,enerd,state,lastbox,
 529                        shake_vir,force_vir,total_vir,pres,
 530                        ekind,mu_tot,constr);
 531             print_ebin(outf->fp_ene,do_ene,FALSE,FALSE,do_log?fplog:NULL,
 532                        step,t,
 533                        eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
 534             write_traj(fplog,cr,outf,mdof_flags,top_global,
 535                        step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
 536             if (bCPT)
 537             {
 538                 nchkpt++;
 539                 bCPT = FALSE;
 540             }
 541             debug_gmx();
 542             if (bLastStep && step_rel == ir->nsteps &&
 543                     (Flags & MD_CONFOUT) && MASTER(cr))
 544             {
 545                 /* x and v have been collected in write_traj,
 546                  * because a checkpoint file will always be written
 547                  * at the last step.
 548                  */
 549                 fprintf(stderr,"\nWriting final coordinates.\n");
 550                 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
 551                 {
 552                     /* Make molecules whole only for confout writing */
 553                     do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
 554                 }
 555                 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
 556                                     *top_global->name,top_global,
 557                                     state_global->x,state_global->v,
 558                                     ir->ePBC,state->box);
 559                 debug_gmx();
 560             }
 561             wallcycle_stop(wcycle,ewcTRAJ);
 562         }
 563         GMX_MPE_LOG(ev_output_finish);
 564
 565
 566         /* Determine the wallclock run time up till now */
 567         run_time = gmx_gettime() - (double)runtime->real;
 568
 569         /* Check whether everything is still allright */
 570         if (((int)gmx_get_stop_condition() > handled_stop_condition)
 571 #ifdef GMX_THREADS
 572             && MASTER(cr)
 573 #endif
 574             )
 575         {
 576            /* this is just make gs.sig compatible with the hack
 577                of sending signals around by MPI_Reduce with together with
 578                other floats */
 579             /* NOTE: this only works for serial code. For code that allows
 580                MPI nodes to propagate their condition, see kernel/md.c*/
 581             if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
 582                 gs.set[eglsSTOPCOND]=1;
 583             if ( gmx_get_stop_condition() == gmx_stop_cond_next )
 584                 gs.set[eglsSTOPCOND]=1;
 585             /* < 0 means stop at next step, > 0 means stop at next NS step */
 586             if (fplog)
 587             {
 588                 fprintf(fplog,
 589                         "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
 590                         gmx_get_signal_name(),
 591                         gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
 592                 fflush(fplog);
 593             }
 594             fprintf(stderr,
 595                     "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
 596                     gmx_get_signal_name(),
 597                     gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
 598             fflush(stderr);
 599             handled_stop_condition=(int)gmx_get_stop_condition();
 600         }
 601         else if (MASTER(cr) &&
 602                  (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
 603                  gs.set[eglsSTOPCOND] == 0)
 604         {
 605             /* Signal to terminate the run */
 606             gs.set[eglsSTOPCOND] = 1;
 607             if (fplog)
 608             {
 609                 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
 610             }
 611             fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
 612         }
 613
 614         /* checkpoints */
 615         if (MASTER(cr) && (cpt_period >= 0 &&
 616                            (cpt_period == 0 ||
 617                             run_time >= nchkpt*cpt_period*60.0)) &&
 618                 gs.set[eglsCHKPT] == 0)
 619         {
 620             gs.set[eglsCHKPT] = 1;
 621         }
 622
 623         /* Time for performance */
 624         if (((step % stepout) == 0) || bLastStep)
 625         {
 626             runtime_upd_proc(runtime);
 627         }
 628
 629         if (do_per_step(step,ir->nstlog))
 630         {
 631             if (fflush(fplog) != 0)
 632             {
 633                 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
 634             }
 635         }
 636
 637         /* Remaining runtime */
 638         if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
 639         {
 640             print_time(stderr,runtime,step,ir,cr);
 641         }
 642
 643         bFirstStep = FALSE;
 644         bInitStep = FALSE;
 645         bStartingFromCpt = FALSE;
 646         step++;
 647         step_rel++;
 648
 649         openmm_take_one_step(openmmData);
 650     }
 651     /* End of main MD loop */
 652     debug_gmx();
 653
 654     /* Stop the time */
 655     runtime_end(runtime);
 656
 657     if (MASTER(cr))
 658     {
 659         if (ir->nstcalcenergy > 0)
 660         {
 661             print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
 662                        eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
 663         }
 664     }
 665
 666     openmm_cleanup(fplog, openmmData);
 667
 668     done_mdoutf(outf);
 669
 670     debug_gmx();
 671
 672     runtime->nsteps_done = step_rel;
 673
 674     return 0;
 675 }