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