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