1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
43 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
79 #include "mpelogging.h"
85 #include "compute_io.h"
87 #include "checkpoint.h"
88 #include "mtop_util.h"
89 #include "sighandler.h"
97 /* include even when OpenMM not used to force compilation of do_md_openmm */
98 #include "openmm_wrapper.h"
100 /* simulation conditions to transmit */
101 enum { eglsNABNSB, eglsCHKPT, eglsSTOPCOND, eglsRESETCOUNTERS, eglsNR };
105 int nstms; /* The frequency for intersimulation communication */
106 int sig[eglsNR]; /* The signal set by one process in do_md */
107 int set[eglsNR]; /* The communicated signal, equal for all processes */
111 static int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
119 gmx_sumi_sim(ms->nsim,buf,ms);
122 for (s=0; s<ms->nsim; s++)
124 bPos = bPos && (buf[s] > 0);
125 bEqual = bEqual && (buf[s] == buf[0]);
131 nmin = min(nmin,buf[0]);
135 /* Find the least common multiple */
136 for (d=2; d<nmin; d++)
139 while (s < ms->nsim && d % buf[s] == 0)
145 /* We found the LCM and it is less than nmin */
157 static int multisim_nstsimsync(const t_commrec *cr,
158 const t_inputrec *ir,int repl_ex_nst)
165 nmin = multisim_min(cr->ms,nmin,ir->nstlist);
166 nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
167 nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
170 gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
172 /* Avoid inter-simulation communication at every (second) step */
179 gmx_bcast(sizeof(int),&nmin,cr);
184 static void init_global_signals(globsig_t *gs,const t_commrec *cr,
185 const t_inputrec *ir,int repl_ex_nst)
191 for (i=0; i<eglsNR; i++)
199 double do_md_openmm(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
200 const output_env_t oenv, bool bVerbose,bool bCompact,
202 gmx_vsite_t *vsite,gmx_constr_t constr,
203 int stepout,t_inputrec *ir,
204 gmx_mtop_t *top_global,
206 t_state *state_global,
208 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
209 gmx_edsam_t ed,t_forcerec *fr,
210 int repl_ex_nst,int repl_ex_seed,
211 real cpt_period,real max_hours,
212 const char *deviceOptions,
214 gmx_runtime_t *runtime)
217 gmx_large_int_t step,step_rel;
221 bFirstStep,bStateFromTPX,bLastStep,bStartingFromCpt;
223 bool do_ene,do_log, do_verbose,
225 tensor force_vir,shake_vir,total_vir,pres;
232 t_mdebin *mdebin=NULL;
237 gmx_enerdata_t *enerd;
239 gmx_global_stat_t gstat;
240 gmx_update_t upd=NULL;
244 gmx_groups_t *groups;
245 gmx_ekindata_t *ekind, *ekind_save;
249 real reset_counters=0,reset_counters_now=0;
250 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
251 int handled_stop_condition=gmx_stop_cond_none;
253 const char *ommOptions = NULL;
256 bAppend = (Flags & MD_APPENDFILES);
257 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
259 groups = &top_global->groups;
262 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
263 nrnb,top_global,&upd,
264 nfile,fnm,&outf,&mdebin,
265 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
267 clear_mat(total_vir);
269 /* Energy terms and groups */
271 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
272 snew(f,top_global->natoms);
274 /* Kinetic energy data */
276 init_ekindata(fplog,top_global,&(ir->opts),ekind);
277 /* needed for iteration of constraints */
279 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
280 /* Copy the cos acceleration to the groups struct */
281 ekind->cosacc.cos_accel = ir->cos_accel;
283 gstat = global_stat_init(ir);
287 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
288 if ((io > 2000) && MASTER(cr))
290 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
294 top = gmx_mtop_generate_local_top(top_global,ir);
297 a1 = top_global->natoms;
299 state = partdec_init_local_state(cr,state_global);
302 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
306 set_vsite_top(vsite,top,mdatoms,cr);
309 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
311 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
314 update_mdatoms(mdatoms,state->lambda);
316 if (deviceOptions[0]=='\0')
318 /* empty options, which should default to OpenMM in this build */
319 ommOptions=deviceOptions;
323 if (gmx_strncasecmp(deviceOptions,"OpenMM",6)!=0)
325 gmx_fatal(FARGS, "This Gromacs version currently only works with OpenMM. Use -device \"OpenMM:<options>\"");
329 ommOptions=strchr(deviceOptions,':');
330 if (NULL!=ommOptions)
332 /* Increase the pointer to skip the colon */
338 openmmData = openmm_init(fplog, ommOptions, ir, top_global, top, mdatoms, fr, state);
342 /* Update mdebin with energy history if appending to output files */
343 if ( Flags & MD_APPENDFILES )
345 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
347 /* Set the initial energy history in state to zero by updating once */
348 update_energyhistory(&state_global->enerhist,mdebin);
353 set_constraints(constr,top,ir,mdatoms,cr);
356 if (!ir->bContinuation)
358 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
360 /* Set the velocities of frozen particles to zero */
361 for (i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
363 for (m=0; m<DIM; m++)
365 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
375 /* Constrain the initial coordinates and velocities */
376 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
377 graph,cr,nrnb,fr,top,shake_vir);
381 /* Construct the virtual sites for the initial configuration */
382 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
383 top->idef.iparams,top->idef.il,
384 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
393 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
394 fprintf(stderr,"starting mdrun '%s'\n",
395 *(top_global->name));
398 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
402 sprintf(tbuf,"%s","infinite");
404 if (ir->init_step > 0)
406 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
407 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
408 gmx_step_str(ir->init_step,sbuf2),
409 ir->init_step*ir->delta_t);
413 fprintf(stderr,"%s steps, %s ps.\n",
414 gmx_step_str(ir->nsteps,sbuf),tbuf);
420 /* Set and write start time */
421 runtime_start(runtime);
422 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
423 wallcycle_start(wcycle,ewcRUN);
427 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
430 /***********************************************************
434 ************************************************************/
436 /* loop over MD steps or if rerunMD to end of input trajectory */
438 /* Skip the first Nose-Hoover integration when we get the state from tpx */
439 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
440 bInitStep = bFirstStep && bStateFromTPX;
441 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
444 init_global_signals(&gs,cr,ir,repl_ex_nst);
446 step = ir->init_step;
451 wallcycle_start(wcycle,ewcSTEP);
453 GMX_MPE_LOG(ev_timestep1);
455 bLastStep = (step_rel == ir->nsteps);
456 t = t0 + step*ir->delta_t;
458 if (gs.set[eglsSTOPCOND] != 0)
463 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
464 do_verbose = bVerbose &&
465 (step % stepout == 0 || bFirstStep || bLastStep);
467 if (MASTER(cr) && do_log)
469 print_ebin_header(fplog,step,t,state->lambda);
472 clear_mat(force_vir);
473 GMX_MPE_LOG(ev_timestep2);
475 /* We write a checkpoint at this MD step when:
476 * either when we signalled through gs (in OpenMM NS works different),
477 * or at the last step (but not when we do not want confout),
478 * but never at the first step.
480 bCPT = ((gs.set[eglsCHKPT] ||
481 (bLastStep && (Flags & MD_CONFOUT))) &&
482 step > ir->init_step );
485 gs.set[eglsCHKPT] = 0;
488 /* Now we have the energies and forces corresponding to the
489 * coordinates at time t. We must output all of this before
491 * for RerunMD t is read from input trajectory
493 GMX_MPE_LOG(ev_output_start);
496 if (do_per_step(step,ir->nstxout))
498 mdof_flags |= MDOF_X;
500 if (do_per_step(step,ir->nstvout))
502 mdof_flags |= MDOF_V;
504 if (do_per_step(step,ir->nstfout))
506 mdof_flags |= MDOF_F;
508 if (do_per_step(step,ir->nstxtcout))
510 mdof_flags |= MDOF_XTC;
514 mdof_flags |= MDOF_CPT;
516 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
520 bX = (mdof_flags & (MDOF_X | MDOF_XTC | MDOF_CPT));
521 bV = (mdof_flags & (MDOF_V | MDOF_CPT));
523 wallcycle_start(wcycle,ewcTRAJ);
524 openmm_copy_state(openmmData, state, &t, f, enerd, bX, bV, 0, 0);
525 wallcycle_stop(wcycle,ewcTRAJ);
528 openmm_take_one_step(openmmData);
530 if (mdof_flags != 0 || do_ene || do_log)
532 wallcycle_start(wcycle,ewcTRAJ);
533 bF = (mdof_flags & MDOF_F);
534 if (bF || do_ene || do_log)
536 openmm_copy_state(openmmData, state, &t, f, enerd, 0, 0, bF, do_ene);
538 upd_mdebin(mdebin, NULL,TRUE,
539 t,mdatoms->tmass,enerd,state,lastbox,
540 shake_vir,force_vir,total_vir,pres,
541 ekind,mu_tot,constr);
542 print_ebin(outf->fp_ene,do_ene,FALSE,FALSE,do_log?fplog:NULL,
544 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
545 write_traj(fplog,cr,outf,mdof_flags,top_global,
546 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
553 if (bLastStep && step_rel == ir->nsteps &&
554 (Flags & MD_CONFOUT) && MASTER(cr))
556 /* x and v have been collected in write_traj,
557 * because a checkpoint file will always be written
560 fprintf(stderr,"\nWriting final coordinates.\n");
561 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
563 /* Make molecules whole only for confout writing */
564 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
566 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
567 *top_global->name,top_global,
568 state_global->x,state_global->v,
569 ir->ePBC,state->box);
572 wallcycle_stop(wcycle,ewcTRAJ);
574 GMX_MPE_LOG(ev_output_finish);
577 /* Determine the wallclock run time up till now */
578 run_time = gmx_gettime() - (double)runtime->real;
580 /* Check whether everything is still allright */
581 if (((int)gmx_get_stop_condition() > handled_stop_condition)
587 /* this is just make gs.sig compatible with the hack
588 of sending signals around by MPI_Reduce with together with
590 /* NOTE: this only works for serial code. For code that allows
591 MPI nodes to propagate their condition, see kernel/md.c*/
592 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
593 gs.set[eglsSTOPCOND]=1;
594 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
595 gs.set[eglsSTOPCOND]=1;
596 /* < 0 means stop at next step, > 0 means stop at next NS step */
600 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
601 gmx_get_signal_name(),
602 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
606 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
607 gmx_get_signal_name(),
608 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
610 handled_stop_condition=(int)gmx_get_stop_condition();
612 else if (MASTER(cr) &&
613 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
614 gs.set[eglsSTOPCOND] == 0)
616 /* Signal to terminate the run */
617 gs.set[eglsSTOPCOND] = 1;
620 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
622 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
626 if (MASTER(cr) && (cpt_period >= 0 &&
628 run_time >= nchkpt*cpt_period*60.0)) &&
629 gs.set[eglsCHKPT] == 0)
631 gs.set[eglsCHKPT] = 1;
634 /* Time for performance */
635 if (((step % stepout) == 0) || bLastStep)
637 runtime_upd_proc(runtime);
640 if (do_per_step(step,ir->nstlog))
642 if (fflush(fplog) != 0)
644 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
648 /* Remaining runtime */
649 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
651 print_time(stderr,runtime,step,ir,cr);
656 bStartingFromCpt = FALSE;
660 /* End of main MD loop */
664 runtime_end(runtime);
668 if (ir->nstcalcenergy > 0)
670 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
671 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
675 openmm_cleanup(fplog, openmmData);
681 runtime->nsteps_done = step_rel;