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