src/kernel/md_openmm.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2010, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <signal.h>
  44 #include <stdlib.h>
  45
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "sysstuff.h"
  49 #include "vec.h"
  50 #include "statutil.h"
  51 #include "vcm.h"
  52 #include "mdebin.h"
  53 #include "nrnb.h"
  54 #include "calcmu.h"
  55 #include "index.h"
  56 #include "vsite.h"
  57 #include "update.h"
  58 #include "ns.h"
  59 #include "trnio.h"
  60 #include "xtcio.h"
  61 #include "mdrun.h"
  62 #include "md_support.h"
  63 #include "confio.h"
  64 #include "network.h"
  65 #include "pull.h"
  66 #include "xvgr.h"
  67 #include "physics.h"
  68 #include "names.h"
  69 #include "xmdrun.h"
  70 #include "ionize.h"
  71 #include "disre.h"
  72 #include "orires.h"
  73 #include "pme.h"
  74 #include "mdatoms.h"
  75 #include "qmmm.h"
  76 #include "mpelogging.h"
  77 #include "domdec.h"
  78 #include "partdec.h"
  79 #include "topsort.h"
  80 #include "coulomb.h"
  81 #include "constr.h"
  82 #include "compute_io.h"
  83 #include "mvdata.h"
  84 #include "checkpoint.h"
  85 #include "mtop_util.h"
  86 #include "sighandler.h"
  87 #include "genborn.h"
  88 #include "string2.h"
  89 #include "copyrite.h"
  90 #include "membed.h"
  91
  92 #ifdef GMX_THREAD_MPI
  93 #include "tmpi.h"
  94 #endif
  95
  96 /* include even when OpenMM not used to force compilation of do_md_openmm */
  97 #include "openmm_wrapper.h"
  98
  99 double do_md_openmm(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 100                     const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
 101                     int nstglobalcomm,
 102                     gmx_vsite_t *vsite,gmx_constr_t constr,
 103                     int stepout,t_inputrec *ir,
 104                     gmx_mtop_t *top_global,
 105                     t_fcdata *fcd,
 106                     t_state *state_global,
 107                     t_mdatoms *mdatoms,
 108                     t_nrnb *nrnb,gmx_wallcycle_t wcycle,
 109                     gmx_edsam_t ed,t_forcerec *fr,
 110                     int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
 111                     gmx_membed_t membed,
 112                     real cpt_period,real max_hours,
 113                     const char *deviceOptions,
 114                     unsigned long Flags,
 115                     gmx_runtime_t *runtime)
 116 {
 117     gmx_mdoutf_t *outf;
 118     gmx_large_int_t step,step_rel;
 119     double     run_time;
 120     double     t,t0,lam0;
 121     gmx_bool       bSimAnn,
 122     bFirstStep,bStateFromTPX,bLastStep,bStartingFromCpt;
 123     gmx_bool       bInitStep=TRUE;
 124     gmx_bool       do_ene,do_log, do_verbose,
 125     bX,bV,bF,bCPT;
 126     tensor     force_vir,shake_vir,total_vir,pres;
 127     int        i,m;
 128     int        mdof_flags;
 129     rvec       mu_tot;
 130     t_vcm      *vcm;
 131     int        nchkpt=1;
 132     gmx_localtop_t *top;
 133     t_mdebin *mdebin;
 134     t_state    *state=NULL;
 135     rvec       *f_global=NULL;
 136     int        n_xtc=-1;
 137     rvec       *x_xtc=NULL;
 138     gmx_enerdata_t *enerd;
 139     rvec       *f=NULL;
 140     gmx_global_stat_t gstat;
 141     gmx_update_t upd=NULL;
 142     t_graph    *graph=NULL;
 143     globsig_t   gs;
 144
 145     gmx_groups_t *groups;
 146     gmx_ekindata_t *ekind, *ekind_save;
 147     gmx_bool        bAppend;
 148     int         a0,a1;
 149     matrix      lastbox;
 150     real        reset_counters=0,reset_counters_now=0;
 151     char        sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
 152     int         handled_stop_condition=gmx_stop_cond_none;
 153
 154     const char *ommOptions = NULL;
 155     void   *openmmData;
 156
 157 #ifdef GMX_DOUBLE
 158     /* Checks in cmake should prevent the compilation in double precision
 159      * with OpenMM, but just to be sure we check here.
 160      */
 161     gmx_fatal(FARGS,"Compilation was performed in double precision, but OpenMM only supports single precision. If you want to use to OpenMM, compile in single precision.");
 162 #endif
 163
 164     bAppend  = (Flags & MD_APPENDFILES);
 165     check_ir_old_tpx_versions(cr,fplog,ir,top_global);
 166
 167     groups = &top_global->groups;
 168
 169     /* Initial values */
 170     init_md(fplog,cr,ir,oenv,&t,&t0,state_global->lambda,
 171             &(state_global->fep_state),&lam0,
 172             nrnb,top_global,&upd,
 173             nfile,fnm,&outf,&mdebin,
 174             force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
 175
 176     clear_mat(total_vir);
 177     clear_mat(pres);
 178     /* Energy terms and groups */
 179     snew(enerd,1);
 180     init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,
 181                   enerd);
 182     snew(f,top_global->natoms);
 183
 184     /* Kinetic energy data */
 185     snew(ekind,1);
 186     init_ekindata(fplog,top_global,&(ir->opts),ekind);
 187     /* needed for iteration of constraints */
 188     snew(ekind_save,1);
 189     init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
 190     /* Copy the cos acceleration to the groups struct */
 191     ekind->cosacc.cos_accel = ir->cos_accel;
 192
 193     gstat = global_stat_init(ir);
 194     debug_gmx();
 195
 196     {
 197         double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
 198         if ((io > 2000) && MASTER(cr))
 199             fprintf(stderr,
 200                     "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
 201                     io);
 202     }
 203
 204     top = gmx_mtop_generate_local_top(top_global,ir);
 205
 206     a0 = 0;
 207     a1 = top_global->natoms;
 208
 209     state = partdec_init_local_state(cr,state_global);
 210     f_global = f;
 211
 212     atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
 213
 214     if (vsite)
 215     {
 216         set_vsite_top(vsite,top,mdatoms,cr);
 217     }
 218
 219     if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
 220     {
 221         graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
 222     }
 223
 224     update_mdatoms(mdatoms,state->lambda[efptMASS]);
 225
 226     if (deviceOptions[0]=='\0')
 227     {
 228         /* empty options, which should default to OpenMM in this build */
 229         ommOptions=deviceOptions;
 230     }
 231     else
 232     {
 233         if (gmx_strncasecmp(deviceOptions,"OpenMM",6)!=0)
 234         {
 235             gmx_fatal(FARGS, "This Gromacs version currently only works with OpenMM. Use -device \"OpenMM:<options>\"");
 236         }
 237         else
 238         {
 239             ommOptions=strchr(deviceOptions,':');
 240             if (NULL!=ommOptions)
 241             {
 242                 /* Increase the pointer to skip the colon */
 243                 ommOptions++;
 244             }
 245         }
 246     }
 247
 248     openmmData = openmm_init(fplog, ommOptions, ir, top_global, top, mdatoms, fr, state);
 249     please_cite(fplog,"Friedrichs2009");
 250
 251     if (MASTER(cr))
 252     {
 253         /* Update mdebin with energy history if appending to output files */
 254         if ( Flags & MD_APPENDFILES )
 255         {
 256             restore_energyhistory_from_state(mdebin,&state_global->enerhist);
 257         }
 258         /* Set the initial energy history in state to zero by updating once */
 259         update_energyhistory(&state_global->enerhist,mdebin);
 260     }
 261
 262     if (constr)
 263     {
 264         set_constraints(constr,top,ir,mdatoms,cr);
 265     }
 266
 267     if (!ir->bContinuation)
 268     {
 269         if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
 270         {
 271             /* Set the velocities of frozen particles to zero */
 272             for (i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
 273             {
 274                 for (m=0; m<DIM; m++)
 275                 {
 276                     if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 277                     {
 278                         state->v[i][m] = 0;
 279                     }
 280                 }
 281             }
 282         }
 283
 284         if (constr)
 285         {
 286             /* Constrain the initial coordinates and velocities */
 287             do_constrain_first(fplog,constr,ir,mdatoms,state,f,
 288                                graph,cr,nrnb,fr,top,shake_vir);
 289         }
 290         if (vsite)
 291         {
 292             /* Construct the virtual sites for the initial configuration */
 293             construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
 294                              top->idef.iparams,top->idef.il,
 295                              fr->ePBC,fr->bMolPBC,graph,cr,state->box);
 296         }
 297     }
 298
 299     debug_gmx();
 300
 301     if (MASTER(cr))
 302     {
 303         char tbuf[20];
 304         fprintf(stderr,"starting mdrun '%s'\n",
 305                 *(top_global->name));
 306         if (ir->nsteps >= 0)
 307         {
 308             sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
 309         }
 310         else
 311         {
 312             sprintf(tbuf,"%s","infinite");
 313         }
 314         if (ir->init_step > 0)
 315         {
 316             fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 317                     gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
 318                     gmx_step_str(ir->init_step,sbuf2),
 319                     ir->init_step*ir->delta_t);
 320         }
 321         else
 322         {
 323             fprintf(stderr,"%s steps, %s ps.\n",
 324                     gmx_step_str(ir->nsteps,sbuf),tbuf);
 325         }
 326     }
 327
 328     fprintf(fplog,"\n");
 329
 330     /* Set and write start time */
 331     runtime_start(runtime);
 332     print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
 333     wallcycle_start(wcycle,ewcRUN);
 334     if (fplog)
 335         fprintf(fplog,"\n");
 336
 337     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
 338
 339     debug_gmx();
 340     /***********************************************************
 341      *
 342      *             Loop over MD steps
 343      *
 344      ************************************************************/
 345
 346     /* loop over MD steps or if rerunMD to end of input trajectory */
 347     bFirstStep = TRUE;
 348     /* Skip the first Nose-Hoover integration when we get the state from tpx */
 349     bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
 350     bInitStep = bFirstStep && bStateFromTPX;
 351     bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
 352     bLastStep = FALSE;
 353
 354     init_global_signals(&gs,cr,ir,repl_ex_nst);
 355
 356     step = ir->init_step;
 357     step_rel = 0;
 358
 359     while (!bLastStep)
 360     {
 361         wallcycle_start(wcycle,ewcSTEP);
 362
 363         GMX_MPE_LOG(ev_timestep1);
 364
 365         bLastStep = (step_rel == ir->nsteps);
 366         t = t0 + step*ir->delta_t;
 367
 368         if (gs.set[eglsSTOPCOND] != 0)
 369         {
 370             bLastStep = TRUE;
 371         }
 372
 373         do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
 374         do_verbose = bVerbose &&
 375                      (step % stepout == 0 || bFirstStep || bLastStep);
 376
 377         if (MASTER(cr) && do_log)
 378         {
 379             print_ebin_header(fplog,step,t,state->lambda[efptFEP]);
 380         }
 381
 382         clear_mat(force_vir);
 383         GMX_MPE_LOG(ev_timestep2);
 384
 385         /* We write a checkpoint at this MD step when:
 386          * either when we signalled through gs (in OpenMM NS works different),
 387          * or at the last step (but not when we do not want confout),
 388          * but never at the first step.
 389          */
 390         bCPT = ((gs.set[eglsCHKPT] ||
 391                  (bLastStep && (Flags & MD_CONFOUT))) &&
 392                 step > ir->init_step );
 393         if (bCPT)
 394         {
 395             gs.set[eglsCHKPT] = 0;
 396         }
 397
 398         /* Now we have the energies and forces corresponding to the
 399          * coordinates at time t. We must output all of this before
 400          * the update.
 401          * for RerunMD t is read from input trajectory
 402          */
 403         GMX_MPE_LOG(ev_output_start);
 404
 405         mdof_flags = 0;
 406         if (do_per_step(step,ir->nstxout))
 407         {
 408             mdof_flags |= MDOF_X;
 409         }
 410         if (do_per_step(step,ir->nstvout))
 411         {
 412             mdof_flags |= MDOF_V;
 413         }
 414         if (do_per_step(step,ir->nstfout))
 415         {
 416             mdof_flags |= MDOF_F;
 417         }
 418         if (do_per_step(step,ir->nstxtcout))
 419         {
 420             mdof_flags |= MDOF_XTC;
 421         }
 422         if (bCPT)
 423         {
 424             mdof_flags |= MDOF_CPT;
 425         };
 426         do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
 427
 428         if (mdof_flags != 0 || do_ene || do_log)
 429         {
 430             wallcycle_start(wcycle,ewcTRAJ);
 431             bF = (mdof_flags & MDOF_F);
 432             bX = (mdof_flags & (MDOF_X | MDOF_XTC | MDOF_CPT));
 433             bV = (mdof_flags & (MDOF_V | MDOF_CPT));
 434
 435             openmm_copy_state(openmmData, state, &t, f, enerd, bX, bV, bF, do_ene);
 436
 437             upd_mdebin(mdebin,FALSE,TRUE,
 438                        t,mdatoms->tmass,enerd,state,ir->fepvals,ir->expandedvals,lastbox,
 439                        shake_vir,force_vir,total_vir,pres,
 440                        ekind,mu_tot,constr);
 441             print_ebin(outf->fp_ene,do_ene,FALSE,FALSE,do_log?fplog:NULL,
 442                        step,t,
 443                        eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
 444             write_traj(fplog,cr,outf,mdof_flags,top_global,
 445                        step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
 446             if (bCPT)
 447             {
 448                 nchkpt++;
 449                 bCPT = FALSE;
 450             }
 451             debug_gmx();
 452             if (bLastStep && step_rel == ir->nsteps &&
 453                     (Flags & MD_CONFOUT) && MASTER(cr))
 454             {
 455                 /* x and v have been collected in write_traj,
 456                  * because a checkpoint file will always be written
 457                  * at the last step.
 458                  */
 459                 fprintf(stderr,"\nWriting final coordinates.\n");
 460                 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
 461                 {
 462                     /* Make molecules whole only for confout writing */
 463                     do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
 464                 }
 465                 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
 466                                     *top_global->name,top_global,
 467                                     state_global->x,state_global->v,
 468                                     ir->ePBC,state->box);
 469                 debug_gmx();
 470             }
 471             wallcycle_stop(wcycle,ewcTRAJ);
 472         }
 473         GMX_MPE_LOG(ev_output_finish);
 474
 475
 476         /* Determine the wallclock run time up till now */
 477         run_time = gmx_gettime() - (double)runtime->real;
 478
 479         /* Check whether everything is still allright */
 480         if (((int)gmx_get_stop_condition() > handled_stop_condition)
 481 #ifdef GMX_THREAD_MPI
 482             && MASTER(cr)
 483 #endif
 484             )
 485         {
 486            /* this is just make gs.sig compatible with the hack
 487                of sending signals around by MPI_Reduce with together with
 488                other floats */
 489             /* NOTE: this only works for serial code. For code that allows
 490                MPI nodes to propagate their condition, see kernel/md.c*/
 491             if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
 492                 gs.set[eglsSTOPCOND]=1;
 493             if ( gmx_get_stop_condition() == gmx_stop_cond_next )
 494                 gs.set[eglsSTOPCOND]=1;
 495             /* < 0 means stop at next step, > 0 means stop at next NS step */
 496             if (fplog)
 497             {
 498                 fprintf(fplog,
 499                         "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
 500                         gmx_get_signal_name(),
 501                         gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
 502                 fflush(fplog);
 503             }
 504             fprintf(stderr,
 505                     "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
 506                     gmx_get_signal_name(),
 507                     gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
 508             fflush(stderr);
 509             handled_stop_condition=(int)gmx_get_stop_condition();
 510         }
 511         else if (MASTER(cr) &&
 512                  (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
 513                  gs.set[eglsSTOPCOND] == 0)
 514         {
 515             /* Signal to terminate the run */
 516             gs.set[eglsSTOPCOND] = 1;
 517             if (fplog)
 518             {
 519                 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
 520             }
 521             fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
 522         }
 523
 524         /* checkpoints */
 525         if (MASTER(cr) && (cpt_period >= 0 &&
 526                            (cpt_period == 0 ||
 527                             run_time >= nchkpt*cpt_period*60.0)) &&
 528                 gs.set[eglsCHKPT] == 0)
 529         {
 530             gs.set[eglsCHKPT] = 1;
 531         }
 532
 533         /* Time for performance */
 534         if (((step % stepout) == 0) || bLastStep)
 535         {
 536             runtime_upd_proc(runtime);
 537         }
 538
 539         if (do_per_step(step,ir->nstlog))
 540         {
 541             if (fflush(fplog) != 0)
 542             {
 543                 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
 544             }
 545         }
 546
 547         /* Remaining runtime */
 548         if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
 549         {
 550             print_time(stderr,runtime,step,ir,cr);
 551         }
 552
 553         bFirstStep = FALSE;
 554         bInitStep = FALSE;
 555         bStartingFromCpt = FALSE;
 556         step++;
 557         step_rel++;
 558
 559         openmm_take_one_step(openmmData);
 560     }
 561     /* End of main MD loop */
 562     debug_gmx();
 563
 564     /* Stop the time */
 565     runtime_end(runtime);
 566
 567     if (MASTER(cr))
 568     {
 569         if (ir->nstcalcenergy > 0)
 570         {
 571             print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
 572                        eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
 573         }
 574     }
 575
 576     openmm_cleanup(fplog, openmmData);
 577
 578     done_mdoutf(outf);
 579
 580     debug_gmx();
 581
 582     runtime->nsteps_done = step_rel;
 583
 584     return 0;
 585 }