src/kernel/md_openmm.c

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