2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, 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.
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.
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.
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.
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.
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.
58 #include "md_support.h"
59 #include "md_logging.h"
74 #include "mpelogging.h"
76 #include "domdec_network.h"
82 #include "compute_io.h"
84 #include "checkpoint.h"
85 #include "mtop_util.h"
86 #include "sighandler.h"
89 #include "pme_loadbal.h"
92 #include "types/nlistheuristics.h"
93 #include "types/iteratedconstraints.h"
94 #include "nbnxn_cuda_data_mgmt.h"
104 #include "corewrap.h"
107 static void reset_all_counters(FILE *fplog
, t_commrec
*cr
,
108 gmx_large_int_t step
,
109 gmx_large_int_t
*step_rel
, t_inputrec
*ir
,
110 gmx_wallcycle_t wcycle
, t_nrnb
*nrnb
,
111 gmx_runtime_t
*runtime
,
112 nbnxn_cuda_ptr_t cu_nbv
)
114 char sbuf
[STEPSTRSIZE
];
116 /* Reset all the counters related to performance over the run */
117 md_print_warn(cr
, fplog
, "step %s: resetting all time and cycle counters\n",
118 gmx_step_str(step
, sbuf
));
122 nbnxn_cuda_reset_timings(cu_nbv
);
125 wallcycle_stop(wcycle
, ewcRUN
);
126 wallcycle_reset_all(wcycle
);
127 if (DOMAINDECOMP(cr
))
129 reset_dd_statistics_counters(cr
->dd
);
132 ir
->init_step
+= *step_rel
;
133 ir
->nsteps
-= *step_rel
;
135 wallcycle_start(wcycle
, ewcRUN
);
136 runtime_start(runtime
);
137 print_date_and_time(fplog
, cr
->nodeid
, "Restarted time", runtime
);
140 double do_md(FILE *fplog
, t_commrec
*cr
, int nfile
, const t_filenm fnm
[],
141 const output_env_t oenv
, gmx_bool bVerbose
, gmx_bool bCompact
,
143 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
144 int stepout
, t_inputrec
*ir
,
145 gmx_mtop_t
*top_global
,
147 t_state
*state_global
,
149 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
150 gmx_edsam_t ed
, t_forcerec
*fr
,
151 int repl_ex_nst
, int repl_ex_nex
, int repl_ex_seed
, gmx_membed_t membed
,
152 real cpt_period
, real max_hours
,
153 const char *deviceOptions
,
155 gmx_runtime_t
*runtime
)
158 gmx_large_int_t step
, step_rel
;
160 double t
, t0
, lam0
[efptNR
];
161 gmx_bool bGStatEveryStep
, bGStat
, bCalcVir
, bCalcEner
;
162 gmx_bool bNS
, bNStList
, bSimAnn
, bStopCM
, bRerunMD
, bNotLastFrame
= FALSE
,
163 bFirstStep
, bStateFromCP
, bStateFromTPX
, bInitStep
, bLastStep
,
164 bBornRadii
, bStartingFromCpt
;
165 gmx_bool bDoDHDL
= FALSE
, bDoFEP
= FALSE
, bDoExpanded
= FALSE
;
166 gmx_bool do_ene
, do_log
, do_verbose
, bRerunWarnNoV
= TRUE
,
167 bForceUpdate
= FALSE
, bCPT
;
169 gmx_bool bMasterState
;
170 int force_flags
, cglo_flags
;
171 tensor force_vir
, shake_vir
, total_vir
, tmp_vir
, pres
;
176 t_state
*bufstate
= NULL
;
177 matrix
*scale_tot
, pcoupl_mu
, M
, ebox
;
180 gmx_repl_ex_t repl_ex
= NULL
;
183 t_mdebin
*mdebin
= NULL
;
184 t_state
*state
= NULL
;
185 rvec
*f_global
= NULL
;
188 gmx_enerdata_t
*enerd
;
190 gmx_global_stat_t gstat
;
191 gmx_update_t upd
= NULL
;
192 t_graph
*graph
= NULL
;
194 gmx_rng_t mcrng
= NULL
;
196 gmx_groups_t
*groups
;
197 gmx_ekindata_t
*ekind
, *ekind_save
;
198 gmx_shellfc_t shellfc
;
199 int count
, nconverged
= 0;
202 gmx_bool bIonize
= FALSE
;
203 gmx_bool bTCR
= FALSE
, bConverged
= TRUE
, bOK
, bSumEkinhOld
, bExchanged
;
205 gmx_bool bResetCountersHalfMaxH
= FALSE
;
206 gmx_bool bVV
, bIterativeCase
, bFirstIterate
, bTemp
, bPres
, bTrotter
;
207 gmx_bool bUpdateDoLR
;
208 real mu_aver
= 0, dvdl_constr
;
209 int a0
, a1
, gnx
= 0, ii
;
210 atom_id
*grpindex
= NULL
;
212 t_coupl_rec
*tcr
= NULL
;
213 rvec
*xcopy
= NULL
, *vcopy
= NULL
, *cbuf
= NULL
;
214 matrix boxcopy
= {{0}}, lastbox
;
216 real fom
, oldfom
, veta_save
, pcurr
, scalevir
, tracevir
;
223 real saved_conserved_quantity
= 0;
228 char sbuf
[STEPSTRSIZE
], sbuf2
[STEPSTRSIZE
];
229 int handled_stop_condition
= gmx_stop_cond_none
; /* compare to get_stop_condition*/
230 gmx_iterate_t iterate
;
231 gmx_large_int_t multisim_nsteps
= -1; /* number of steps to do before first multisim
232 simulation stops. If equal to zero, don't
233 communicate any more between multisims.*/
234 /* PME load balancing data for GPU kernels */
235 pme_load_balancing_t pme_loadbal
= NULL
;
237 gmx_bool bPMETuneTry
= FALSE
, bPMETuneRunning
= FALSE
;
240 /* Temporary addition for FAHCORE checkpointing */
244 /* Check for special mdrun options */
245 bRerunMD
= (Flags
& MD_RERUN
);
246 bIonize
= (Flags
& MD_IONIZE
);
247 bFFscan
= (Flags
& MD_FFSCAN
);
248 bAppend
= (Flags
& MD_APPENDFILES
);
249 if (Flags
& MD_RESETCOUNTERSHALFWAY
)
253 /* Signal to reset the counters half the simulation steps. */
254 wcycle_set_reset_counters(wcycle
, ir
->nsteps
/2);
256 /* Signal to reset the counters halfway the simulation time. */
257 bResetCountersHalfMaxH
= (max_hours
> 0);
260 /* md-vv uses averaged full step velocities for T-control
261 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
262 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
264 if (bVV
) /* to store the initial velocities while computing virial */
266 snew(cbuf
, top_global
->natoms
);
268 /* all the iteratative cases - only if there are constraints */
269 bIterativeCase
= ((IR_NPH_TROTTER(ir
) || IR_NPT_TROTTER(ir
)) && (constr
) && (!bRerunMD
));
270 gmx_iterate_init(&iterate
, FALSE
); /* The default value of iterate->bIterationActive is set to
271 false in this step. The correct value, true or false,
272 is set at each step, as it depends on the frequency of temperature
273 and pressure control.*/
274 bTrotter
= (bVV
&& (IR_NPT_TROTTER(ir
) || IR_NPH_TROTTER(ir
) || IR_NVT_TROTTER(ir
)));
278 /* Since we don't know if the frames read are related in any way,
279 * rebuild the neighborlist at every step.
282 ir
->nstcalcenergy
= 1;
286 check_ir_old_tpx_versions(cr
, fplog
, ir
, top_global
);
288 nstglobalcomm
= check_nstglobalcomm(fplog
, cr
, nstglobalcomm
, ir
);
289 bGStatEveryStep
= (nstglobalcomm
== 1);
291 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
294 "To reduce the energy communication with nstlist = -1\n"
295 "the neighbor list validity should not be checked at every step,\n"
296 "this means that exact integration is not guaranteed.\n"
297 "The neighbor list validity is checked after:\n"
298 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
299 "In most cases this will result in exact integration.\n"
300 "This reduces the energy communication by a factor of 2 to 3.\n"
301 "If you want less energy communication, set nstlist > 3.\n\n");
304 if (bRerunMD
|| bFFscan
)
308 groups
= &top_global
->groups
;
311 init_md(fplog
, cr
, ir
, oenv
, &t
, &t0
, state_global
->lambda
,
312 &(state_global
->fep_state
), lam0
,
313 nrnb
, top_global
, &upd
,
314 nfile
, fnm
, &outf
, &mdebin
,
315 force_vir
, shake_vir
, mu_tot
, &bSimAnn
, &vcm
, state_global
, Flags
);
317 clear_mat(total_vir
);
319 /* Energy terms and groups */
321 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
, ir
->fepvals
->n_lambda
,
323 if (DOMAINDECOMP(cr
))
329 snew(f
, top_global
->natoms
);
332 /* Kinetic energy data */
334 init_ekindata(fplog
, top_global
, &(ir
->opts
), ekind
);
335 /* needed for iteration of constraints */
337 init_ekindata(fplog
, top_global
, &(ir
->opts
), ekind_save
);
338 /* Copy the cos acceleration to the groups struct */
339 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
341 gstat
= global_stat_init(ir
);
344 /* Check for polarizable models and flexible constraints */
345 shellfc
= init_shell_flexcon(fplog
,
346 top_global
, n_flexible_constraints(constr
),
347 (ir
->bContinuation
||
348 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
349 NULL
: state_global
->x
);
353 #ifdef GMX_THREAD_MPI
354 tMPI_Thread_mutex_lock(&deform_init_box_mutex
);
356 set_deform_reference_box(upd
,
357 deform_init_init_step_tpx
,
358 deform_init_box_tpx
);
359 #ifdef GMX_THREAD_MPI
360 tMPI_Thread_mutex_unlock(&deform_init_box_mutex
);
365 double io
= compute_io(ir
, top_global
->natoms
, groups
, mdebin
->ebin
->nener
, 1);
366 if ((io
> 2000) && MASTER(cr
))
369 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
374 if (DOMAINDECOMP(cr
))
376 top
= dd_init_local_top(top_global
);
379 dd_init_local_state(cr
->dd
, state_global
, state
);
381 if (DDMASTER(cr
->dd
) && ir
->nstfout
)
383 snew(f_global
, state_global
->natoms
);
390 /* Initialize the particle decomposition and split the topology */
391 top
= split_system(fplog
, top_global
, ir
, cr
);
393 pd_cg_range(cr
, &fr
->cg0
, &fr
->hcg
);
394 pd_at_range(cr
, &a0
, &a1
);
398 top
= gmx_mtop_generate_local_top(top_global
, ir
);
401 a1
= top_global
->natoms
;
404 forcerec_set_excl_load(fr
, top
, cr
);
406 state
= partdec_init_local_state(cr
, state_global
);
409 atoms2md(top_global
, ir
, 0, NULL
, a0
, a1
-a0
, mdatoms
);
413 set_vsite_top(vsite
, top
, mdatoms
, cr
);
416 if (ir
->ePBC
!= epbcNONE
&& !fr
->bMolPBC
)
418 graph
= mk_graph(fplog
, &(top
->idef
), 0, top_global
->natoms
, FALSE
, FALSE
);
423 make_local_shells(cr
, mdatoms
, shellfc
);
426 setup_bonded_threading(fr
, &top
->idef
);
428 if (ir
->pull
&& PAR(cr
))
430 dd_make_local_pull_groups(NULL
, ir
->pull
, mdatoms
);
434 if (DOMAINDECOMP(cr
))
436 /* Distribute the charge groups over the nodes from the master node */
437 dd_partition_system(fplog
, ir
->init_step
, cr
, TRUE
, 1,
438 state_global
, top_global
, ir
,
439 state
, &f
, mdatoms
, top
, fr
,
440 vsite
, shellfc
, constr
,
441 nrnb
, wcycle
, FALSE
);
445 update_mdatoms(mdatoms
, state
->lambda
[efptMASS
]);
447 if (opt2bSet("-cpi", nfile
, fnm
))
449 bStateFromCP
= gmx_fexist_master(opt2fn_master("-cpi", nfile
, fnm
, cr
), cr
);
453 bStateFromCP
= FALSE
;
458 init_expanded_ensemble(bStateFromCP
,ir
,&mcrng
,&state
->dfhist
);
465 /* Update mdebin with energy history if appending to output files */
466 if (Flags
& MD_APPENDFILES
)
468 restore_energyhistory_from_state(mdebin
, &state_global
->enerhist
);
472 /* We might have read an energy history from checkpoint,
473 * free the allocated memory and reset the counts.
475 done_energyhistory(&state_global
->enerhist
);
476 init_energyhistory(&state_global
->enerhist
);
479 /* Set the initial energy history in state by updating once */
480 update_energyhistory(&state_global
->enerhist
, mdebin
);
483 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
))
485 /* Set the random state if we read a checkpoint file */
486 set_stochd_state(upd
, state
);
489 if (state
->flags
& (1<<estMC_RNG
))
491 set_mc_state(mcrng
, state
);
494 /* Initialize constraints */
497 if (!DOMAINDECOMP(cr
))
499 set_constraints(constr
, top
, ir
, mdatoms
, cr
);
503 /* Check whether we have to GCT stuff */
504 bTCR
= ftp2bSet(efGCT
, nfile
, fnm
);
509 fprintf(stderr
, "Will do General Coupling Theory!\n");
511 gnx
= top_global
->mols
.nr
;
513 for (i
= 0; (i
< gnx
); i
++)
521 /* We need to be sure replica exchange can only occur
522 * when the energies are current */
523 check_nst_param(fplog
, cr
, "nstcalcenergy", ir
->nstcalcenergy
,
524 "repl_ex_nst", &repl_ex_nst
);
525 /* This check needs to happen before inter-simulation
526 * signals are initialized, too */
528 if (repl_ex_nst
> 0 && MASTER(cr
))
530 repl_ex
= init_replica_exchange(fplog
, cr
->ms
, state_global
, ir
,
531 repl_ex_nst
, repl_ex_nex
, repl_ex_seed
);
534 /* PME tuning is only supported with GPUs or PME nodes and not with rerun.
535 * With perturbed charges with soft-core we should not change the cut-off.
537 if ((Flags
& MD_TUNEPME
) &&
538 EEL_PME(fr
->eeltype
) &&
539 ( (fr
->cutoff_scheme
== ecutsVERLET
&& fr
->nbv
->bUseGPU
) || !(cr
->duty
& DUTY_PME
)) &&
540 !(ir
->efep
!= efepNO
&& mdatoms
->nChargePerturbed
> 0 && ir
->fepvals
->bScCoul
) &&
543 pme_loadbal_init(&pme_loadbal
, ir
, state
->box
, fr
->ic
, fr
->pmedata
);
545 if (cr
->duty
& DUTY_PME
)
547 /* Start tuning right away, as we can't measure the load */
548 bPMETuneRunning
= TRUE
;
552 /* Separate PME nodes, we can measure the PP/PME load balance */
557 if (!ir
->bContinuation
&& !bRerunMD
)
559 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
561 /* Set the velocities of frozen particles to zero */
562 for (i
= mdatoms
->start
; i
< mdatoms
->start
+mdatoms
->homenr
; i
++)
564 for (m
= 0; m
< DIM
; m
++)
566 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
576 /* Constrain the initial coordinates and velocities */
577 do_constrain_first(fplog
, constr
, ir
, mdatoms
, state
, f
,
578 graph
, cr
, nrnb
, fr
, top
, shake_vir
);
582 /* Construct the virtual sites for the initial configuration */
583 construct_vsites(fplog
, vsite
, state
->x
, nrnb
, ir
->delta_t
, NULL
,
584 top
->idef
.iparams
, top
->idef
.il
,
585 fr
->ePBC
, fr
->bMolPBC
, graph
, cr
, state
->box
);
591 /* set free energy calculation frequency as the minimum
592 greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/
593 nstfep
= ir
->fepvals
->nstdhdl
;
596 nstfep
= gmx_greatest_common_divisor(ir
->fepvals
->nstdhdl
,nstfep
);
600 nstfep
= gmx_greatest_common_divisor(repl_ex_nst
,nstfep
);
603 /* I'm assuming we need global communication the first time! MRS */
604 cglo_flags
= (CGLO_TEMPERATURE
| CGLO_GSTAT
605 | ((ir
->comm_mode
!= ecmNO
) ? CGLO_STOPCM
: 0)
606 | (bVV
? CGLO_PRESSURE
: 0)
607 | (bVV
? CGLO_CONSTRAINT
: 0)
608 | (bRerunMD
? CGLO_RERUNMD
: 0)
609 | ((Flags
& MD_READ_EKIN
) ? CGLO_READEKIN
: 0));
611 bSumEkinhOld
= FALSE
;
612 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, state_global
, mdatoms
, nrnb
, vcm
,
613 NULL
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
614 constr
, NULL
, FALSE
, state
->box
,
615 top_global
, &pcurr
, top_global
->natoms
, &bSumEkinhOld
, cglo_flags
);
616 if (ir
->eI
== eiVVAK
)
618 /* a second call to get the half step temperature initialized as well */
619 /* we do the same call as above, but turn the pressure off -- internally to
620 compute_globals, this is recognized as a velocity verlet half-step
621 kinetic energy calculation. This minimized excess variables, but
622 perhaps loses some logic?*/
624 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, state_global
, mdatoms
, nrnb
, vcm
,
625 NULL
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
626 constr
, NULL
, FALSE
, state
->box
,
627 top_global
, &pcurr
, top_global
->natoms
, &bSumEkinhOld
,
628 cglo_flags
&~(CGLO_STOPCM
| CGLO_PRESSURE
));
631 /* Calculate the initial half step temperature, and save the ekinh_old */
632 if (!(Flags
& MD_STARTFROMCPT
))
634 for (i
= 0; (i
< ir
->opts
.ngtc
); i
++)
636 copy_mat(ekind
->tcstat
[i
].ekinh
, ekind
->tcstat
[i
].ekinh_old
);
641 enerd
->term
[F_TEMP
] *= 2; /* result of averages being done over previous and current step,
642 and there is no previous step */
645 /* if using an iterative algorithm, we need to create a working directory for the state. */
648 bufstate
= init_bufstate(state
);
652 snew(xcopy
, state
->natoms
);
653 snew(vcopy
, state
->natoms
);
654 copy_rvecn(state
->x
, xcopy
, 0, state
->natoms
);
655 copy_rvecn(state
->v
, vcopy
, 0, state
->natoms
);
656 copy_mat(state
->box
, boxcopy
);
659 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
660 temperature control */
661 trotter_seq
= init_npt_vars(ir
, state
, &MassQ
, bTrotter
);
665 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
668 "RMS relative constraint deviation after constraining: %.2e\n",
669 constr_rmsd(constr
, FALSE
));
671 if (EI_STATE_VELOCITY(ir
->eI
))
673 fprintf(fplog
, "Initial temperature: %g K\n", enerd
->term
[F_TEMP
]);
677 fprintf(stderr
, "starting md rerun '%s', reading coordinates from"
678 " input trajectory '%s'\n\n",
679 *(top_global
->name
), opt2fn("-rerun", nfile
, fnm
));
682 fprintf(stderr
, "Calculated time to finish depends on nsteps from "
683 "run input file,\nwhich may not correspond to the time "
684 "needed to process input trajectory.\n\n");
690 fprintf(stderr
, "starting mdrun '%s'\n",
691 *(top_global
->name
));
694 sprintf(tbuf
, "%8.1f", (ir
->init_step
+ir
->nsteps
)*ir
->delta_t
);
698 sprintf(tbuf
, "%s", "infinite");
700 if (ir
->init_step
> 0)
702 fprintf(stderr
, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
703 gmx_step_str(ir
->init_step
+ir
->nsteps
, sbuf
), tbuf
,
704 gmx_step_str(ir
->init_step
, sbuf2
),
705 ir
->init_step
*ir
->delta_t
);
709 fprintf(stderr
, "%s steps, %s ps.\n",
710 gmx_step_str(ir
->nsteps
, sbuf
), tbuf
);
713 fprintf(fplog
, "\n");
716 /* Set and write start time */
717 runtime_start(runtime
);
718 print_date_and_time(fplog
, cr
->nodeid
, "Started mdrun", runtime
);
719 wallcycle_start(wcycle
, ewcRUN
);
722 fprintf(fplog
, "\n");
725 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
727 chkpt_ret
= fcCheckPointParallel( cr
->nodeid
,
731 gmx_fatal( 3, __FILE__
, __LINE__
, "Checkpoint error on step %d\n", 0 );
736 /***********************************************************
740 ************************************************************/
742 /* if rerunMD then read coordinates and velocities from input trajectory */
745 if (getenv("GMX_FORCE_UPDATE"))
753 bNotLastFrame
= read_first_frame(oenv
, &status
,
754 opt2fn("-rerun", nfile
, fnm
),
755 &rerun_fr
, TRX_NEED_X
| TRX_READ_V
);
756 if (rerun_fr
.natoms
!= top_global
->natoms
)
759 "Number of atoms in trajectory (%d) does not match the "
760 "run input file (%d)\n",
761 rerun_fr
.natoms
, top_global
->natoms
);
763 if (ir
->ePBC
!= epbcNONE
)
767 gmx_fatal(FARGS
, "Rerun trajectory frame step %d time %f does not contain a box, while pbc is used", rerun_fr
.step
, rerun_fr
.time
);
769 if (max_cutoff2(ir
->ePBC
, rerun_fr
.box
) < sqr(fr
->rlistlong
))
771 gmx_fatal(FARGS
, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr
.step
, rerun_fr
.time
);
778 rerun_parallel_comm(cr
, &rerun_fr
, &bNotLastFrame
);
781 if (ir
->ePBC
!= epbcNONE
)
783 /* Set the shift vectors.
784 * Necessary here when have a static box different from the tpr box.
786 calc_shifts(rerun_fr
.box
, fr
->shift_vec
);
790 /* loop over MD steps or if rerunMD to end of input trajectory */
792 /* Skip the first Nose-Hoover integration when we get the state from tpx */
793 bStateFromTPX
= !bStateFromCP
;
794 bInitStep
= bFirstStep
&& (bStateFromTPX
|| bVV
);
795 bStartingFromCpt
= (Flags
& MD_STARTFROMCPT
) && bInitStep
;
797 bSumEkinhOld
= FALSE
;
800 init_global_signals(&gs
, cr
, ir
, repl_ex_nst
);
802 step
= ir
->init_step
;
805 if (ir
->nstlist
== -1)
807 init_nlistheuristics(&nlh
, bGStatEveryStep
, step
);
810 if (MULTISIM(cr
) && (repl_ex_nst
<= 0 ))
812 /* check how many steps are left in other sims */
813 multisim_nsteps
= get_multisim_nsteps(cr
, ir
->nsteps
);
817 /* and stop now if we should */
818 bLastStep
= (bRerunMD
|| (ir
->nsteps
>= 0 && step_rel
> ir
->nsteps
) ||
819 ((multisim_nsteps
>= 0) && (step_rel
>= multisim_nsteps
)));
820 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
))
823 wallcycle_start(wcycle
, ewcSTEP
);
825 GMX_MPE_LOG(ev_timestep1
);
831 step
= rerun_fr
.step
;
832 step_rel
= step
- ir
->init_step
;
845 bLastStep
= (step_rel
== ir
->nsteps
);
846 t
= t0
+ step
*ir
->delta_t
;
849 if (ir
->efep
!= efepNO
|| ir
->bSimTemp
)
851 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
852 requiring different logic. */
854 set_current_lambdas(step
, ir
->fepvals
, bRerunMD
, &rerun_fr
, state_global
, state
, lam0
);
855 bDoDHDL
= do_per_step(step
, ir
->fepvals
->nstdhdl
);
856 bDoFEP
= (do_per_step(step
, nstfep
) && (ir
->efep
!= efepNO
));
857 bDoExpanded
= (do_per_step(step
, ir
->expandedvals
->nstexpanded
)
858 && (ir
->bExpanded
) && (step
> 0) && (!bStartingFromCpt
));
863 update_annealing_target_temp(&(ir
->opts
), t
);
868 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
870 for (i
= 0; i
< state_global
->natoms
; i
++)
872 copy_rvec(rerun_fr
.x
[i
], state_global
->x
[i
]);
876 for (i
= 0; i
< state_global
->natoms
; i
++)
878 copy_rvec(rerun_fr
.v
[i
], state_global
->v
[i
]);
883 for (i
= 0; i
< state_global
->natoms
; i
++)
885 clear_rvec(state_global
->v
[i
]);
889 fprintf(stderr
, "\nWARNING: Some frames do not contain velocities.\n"
890 " Ekin, temperature and pressure are incorrect,\n"
891 " the virial will be incorrect when constraints are present.\n"
893 bRerunWarnNoV
= FALSE
;
897 copy_mat(rerun_fr
.box
, state_global
->box
);
898 copy_mat(state_global
->box
, state
->box
);
900 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
902 if (DOMAINDECOMP(cr
))
904 gmx_fatal(FARGS
, "Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
908 /* Following is necessary because the graph may get out of sync
909 * with the coordinates if we only have every N'th coordinate set
911 mk_mshift(fplog
, graph
, fr
->ePBC
, state
->box
, state
->x
);
912 shift_self(graph
, state
->box
, state
->x
);
914 construct_vsites(fplog
, vsite
, state
->x
, nrnb
, ir
->delta_t
, state
->v
,
915 top
->idef
.iparams
, top
->idef
.il
,
916 fr
->ePBC
, fr
->bMolPBC
, graph
, cr
, state
->box
);
919 unshift_self(graph
, state
->box
, state
->x
);
924 /* Stop Center of Mass motion */
925 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
, ir
->nstcomm
));
927 /* Copy back starting coordinates in case we're doing a forcefield scan */
930 for (ii
= 0; (ii
< state
->natoms
); ii
++)
932 copy_rvec(xcopy
[ii
], state
->x
[ii
]);
933 copy_rvec(vcopy
[ii
], state
->v
[ii
]);
935 copy_mat(boxcopy
, state
->box
);
940 /* for rerun MD always do Neighbour Searching */
941 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
946 /* Determine whether or not to do Neighbour Searching and LR */
947 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
949 bNS
= (bFirstStep
|| bExchanged
|| bNStList
|| bDoFEP
||
950 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
952 if (bNS
&& ir
->nstlist
== -1)
954 set_nlistheuristics(&nlh
, bFirstStep
|| bExchanged
|| bDoFEP
, step
);
958 /* check whether we should stop because another simulation has
962 if ( (multisim_nsteps
>= 0) && (step_rel
>= multisim_nsteps
) &&
963 (multisim_nsteps
!= ir
->nsteps
) )
970 "Stopping simulation %d because another one has finished\n",
974 gs
.sig
[eglsCHKPT
] = 1;
979 /* < 0 means stop at next step, > 0 means stop at next NS step */
980 if ( (gs
.set
[eglsSTOPCOND
] < 0) ||
981 ( (gs
.set
[eglsSTOPCOND
] > 0) && (bNStList
|| ir
->nstlist
== 0) ) )
986 /* Determine whether or not to update the Born radii if doing GB */
987 bBornRadii
= bFirstStep
;
988 if (ir
->implicit_solvent
&& (step
% ir
->nstgbradii
== 0))
993 do_log
= do_per_step(step
, ir
->nstlog
) || bFirstStep
|| bLastStep
;
994 do_verbose
= bVerbose
&&
995 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
997 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
1001 bMasterState
= TRUE
;
1005 bMasterState
= FALSE
;
1006 /* Correct the new box if it is too skewed */
1007 if (DYNAMIC_BOX(*ir
))
1009 if (correct_box(fplog
, step
, state
->box
, graph
))
1011 bMasterState
= TRUE
;
1014 if (DOMAINDECOMP(cr
) && bMasterState
)
1016 dd_collect_state(cr
->dd
, state
, state_global
);
1020 if (DOMAINDECOMP(cr
))
1022 /* Repartition the domain decomposition */
1023 wallcycle_start(wcycle
, ewcDOMDEC
);
1024 dd_partition_system(fplog
, step
, cr
,
1025 bMasterState
, nstglobalcomm
,
1026 state_global
, top_global
, ir
,
1027 state
, &f
, mdatoms
, top
, fr
,
1028 vsite
, shellfc
, constr
,
1030 do_verbose
&& !bPMETuneRunning
);
1031 wallcycle_stop(wcycle
, ewcDOMDEC
);
1032 /* If using an iterative integrator, reallocate space to match the decomposition */
1036 if (MASTER(cr
) && do_log
&& !bFFscan
)
1038 print_ebin_header(fplog
, step
, t
, state
->lambda
[efptFEP
]); /* can we improve the information printed here? */
1041 if (ir
->efep
!= efepNO
)
1043 update_mdatoms(mdatoms
, state
->lambda
[efptMASS
]);
1046 if ((bRerunMD
&& rerun_fr
.bV
) || bExchanged
)
1049 /* We need the kinetic energy at minus the half step for determining
1050 * the full step kinetic energy and possibly for T-coupling.*/
1051 /* This may not be quite working correctly yet . . . . */
1052 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, state_global
, mdatoms
, nrnb
, vcm
,
1053 wcycle
, enerd
, NULL
, NULL
, NULL
, NULL
, mu_tot
,
1054 constr
, NULL
, FALSE
, state
->box
,
1055 top_global
, &pcurr
, top_global
->natoms
, &bSumEkinhOld
,
1056 CGLO_RERUNMD
| CGLO_GSTAT
| CGLO_TEMPERATURE
);
1058 clear_mat(force_vir
);
1060 /* Ionize the atoms if necessary */
1063 ionize(fplog
, oenv
, mdatoms
, top_global
, t
, ir
, state
->x
, state
->v
,
1064 mdatoms
->start
, mdatoms
->start
+mdatoms
->homenr
, state
->box
, cr
);
1067 /* Update force field in ffscan program */
1070 if (update_forcefield(fplog
,
1072 mdatoms
->nr
, state
->x
, state
->box
))
1080 GMX_MPE_LOG(ev_timestep2
);
1082 /* We write a checkpoint at this MD step when:
1083 * either at an NS step when we signalled through gs,
1084 * or at the last step (but not when we do not want confout),
1085 * but never at the first step or with rerun.
1087 bCPT
= (((gs
.set
[eglsCHKPT
] && (bNS
|| ir
->nstlist
== 0)) ||
1088 (bLastStep
&& (Flags
& MD_CONFOUT
))) &&
1089 step
> ir
->init_step
&& !bRerunMD
);
1092 gs
.set
[eglsCHKPT
] = 0;
1095 /* Determine the energy and pressure:
1096 * at nstcalcenergy steps and at energy output steps (set below).
1098 if (EI_VV(ir
->eI
) && (!bInitStep
))
1100 /* for vv, the first half of the integration actually corresponds
1101 to the previous step. bCalcEner is only required to be evaluated on the 'next' step,
1102 but the virial needs to be calculated on both the current step and the 'next' step. Future
1103 reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
1105 bCalcEner
= do_per_step(step
-1, ir
->nstcalcenergy
);
1106 bCalcVir
= bCalcEner
||
1107 (ir
->epc
!= epcNO
&& (do_per_step(step
, ir
->nstpcouple
) || do_per_step(step
-1, ir
->nstpcouple
)));
1111 bCalcEner
= do_per_step(step
, ir
->nstcalcenergy
);
1112 bCalcVir
= bCalcEner
||
1113 (ir
->epc
!= epcNO
&& do_per_step(step
, ir
->nstpcouple
));
1116 /* Do we need global communication ? */
1117 bGStat
= (bCalcVir
|| bCalcEner
|| bStopCM
||
1118 do_per_step(step
, nstglobalcomm
) || (bVV
&& IR_NVT_TROTTER(ir
) && do_per_step(step
-1, nstglobalcomm
)) ||
1119 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
1121 do_ene
= (do_per_step(step
, ir
->nstenergy
) || bLastStep
);
1123 if (do_ene
|| do_log
)
1130 /* these CGLO_ options remain the same throughout the iteration */
1131 cglo_flags
= ((bRerunMD
? CGLO_RERUNMD
: 0) |
1132 (bGStat
? CGLO_GSTAT
: 0)
1135 force_flags
= (GMX_FORCE_STATECHANGED
|
1136 ((DYNAMIC_BOX(*ir
) || bRerunMD
) ? GMX_FORCE_DYNAMICBOX
: 0) |
1137 GMX_FORCE_ALLFORCES
|
1139 (bCalcVir
? GMX_FORCE_VIRIAL
: 0) |
1140 (bCalcEner
? GMX_FORCE_ENERGY
: 0) |
1141 (bDoFEP
? GMX_FORCE_DHDL
: 0)
1146 if (do_per_step(step
, ir
->nstcalclr
))
1148 force_flags
|= GMX_FORCE_DO_LR
;
1154 /* Now is the time to relax the shells */
1155 count
= relax_shell_flexcon(fplog
, cr
, bVerbose
, bFFscan
? step
+1 : step
,
1156 ir
, bNS
, force_flags
,
1157 bStopCM
, top
, top_global
,
1159 state
, f
, force_vir
, mdatoms
,
1160 nrnb
, wcycle
, graph
, groups
,
1161 shellfc
, fr
, bBornRadii
, t
, mu_tot
,
1162 state
->natoms
, &bConverged
, vsite
,
1173 /* The coordinates (x) are shifted (to get whole molecules)
1175 * This is parallellized as well, and does communication too.
1176 * Check comments in sim_util.c
1178 do_force(fplog
, cr
, ir
, step
, nrnb
, wcycle
, top
, top_global
, groups
,
1179 state
->box
, state
->x
, &state
->hist
,
1180 f
, force_vir
, mdatoms
, enerd
, fcd
,
1181 state
->lambda
, graph
,
1182 fr
, vsite
, mu_tot
, t
, outf
->fp_field
, ed
, bBornRadii
,
1183 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
1186 GMX_BARRIER(cr
->mpi_comm_mygroup
);
1190 mu_aver
= calc_mu_aver(cr
, state
->x
, mdatoms
->chargeA
,
1191 mu_tot
, &top_global
->mols
, mdatoms
, gnx
, grpindex
);
1194 if (bTCR
&& bFirstStep
)
1196 tcr
= init_coupling(fplog
, nfile
, fnm
, cr
, fr
, mdatoms
, &(top
->idef
));
1197 fprintf(fplog
, "Done init_coupling\n");
1201 if (bVV
&& !bStartingFromCpt
&& !bRerunMD
)
1202 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1204 if (ir
->eI
== eiVV
&& bInitStep
)
1206 /* if using velocity verlet with full time step Ekin,
1207 * take the first half step only to compute the
1208 * virial for the first step. From there,
1209 * revert back to the initial coordinates
1210 * so that the input is actually the initial step.
1212 copy_rvecn(state
->v
, cbuf
, 0, state
->natoms
); /* should make this better for parallelizing? */
1216 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1217 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ1
);
1220 /* If we are using twin-range interactions where the long-range component
1221 * is only evaluated every nstcalclr>1 steps, we should do a special update
1222 * step to combine the long-range forces on these steps.
1223 * For nstcalclr=1 this is not done, since the forces would have been added
1224 * directly to the short-range forces already.
1226 bUpdateDoLR
= (fr
->bTwinRange
&& do_per_step(step
, ir
->nstcalclr
));
1228 update_coords(fplog
, step
, ir
, mdatoms
, state
, fr
->bMolPBC
,
1229 f
, bUpdateDoLR
, fr
->f_twin
, fcd
,
1230 ekind
, M
, wcycle
, upd
, bInitStep
, etrtVELOCITY1
,
1231 cr
, nrnb
, constr
, &top
->idef
);
1233 if (bIterativeCase
&& do_per_step(step
-1, ir
->nstpcouple
) && !bInitStep
)
1235 gmx_iterate_init(&iterate
, TRUE
);
1237 /* for iterations, we save these vectors, as we will be self-consistently iterating
1240 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1242 /* save the state */
1243 if (iterate
.bIterationActive
)
1245 copy_coupling_state(state
, bufstate
, ekind
, ekind_save
, &(ir
->opts
));
1248 bFirstIterate
= TRUE
;
1249 while (bFirstIterate
|| iterate
.bIterationActive
)
1251 if (iterate
.bIterationActive
)
1253 copy_coupling_state(bufstate
, state
, ekind_save
, ekind
, &(ir
->opts
));
1254 if (bFirstIterate
&& bTrotter
)
1256 /* The first time through, we need a decent first estimate
1257 of veta(t+dt) to compute the constraints. Do
1258 this by computing the box volume part of the
1259 trotter integration at this time. Nothing else
1260 should be changed by this routine here. If
1261 !(first time), we start with the previous value
1264 veta_save
= state
->veta
;
1265 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ0
);
1266 vetanew
= state
->veta
;
1267 state
->veta
= veta_save
;
1272 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
) /* Why is rerun_fr.bV here? Unclear. */
1274 update_constraints(fplog
, step
, NULL
, ir
, ekind
, mdatoms
,
1275 state
, fr
->bMolPBC
, graph
, f
,
1276 &top
->idef
, shake_vir
, NULL
,
1277 cr
, nrnb
, wcycle
, upd
, constr
,
1278 bInitStep
, TRUE
, bCalcVir
, vetanew
);
1280 if (!bOK
&& !bFFscan
)
1282 gmx_fatal(FARGS
, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
1288 /* Need to unshift here if a do_force has been
1289 called in the previous step */
1290 unshift_self(graph
, state
->box
, state
->x
);
1293 /* if VV, compute the pressure and constraints */
1294 /* For VV2, we strictly only need this if using pressure
1295 * control, but we really would like to have accurate pressures
1297 * Think about ways around this in the future?
1298 * For now, keep this choice in comments.
1300 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1301 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1303 bTemp
= ((ir
->eI
== eiVV
&& (!bInitStep
)) || (ir
->eI
== eiVVAK
));
1304 if (bCalcEner
&& ir
->eI
== eiVVAK
) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
1306 bSumEkinhOld
= TRUE
;
1308 /* for vv, the first half of the integration actually corresponds to the previous step.
1309 So we need information from the last step in the first half of the integration */
1310 if (bGStat
|| do_per_step(step
-1, nstglobalcomm
))
1312 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, state_global
, mdatoms
, nrnb
, vcm
,
1313 wcycle
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
1314 constr
, NULL
, FALSE
, state
->box
,
1315 top_global
, &pcurr
, top_global
->natoms
, &bSumEkinhOld
,
1318 | (bTemp
? CGLO_TEMPERATURE
: 0)
1319 | (bPres
? CGLO_PRESSURE
: 0)
1320 | (bPres
? CGLO_CONSTRAINT
: 0)
1321 | ((iterate
.bIterationActive
) ? CGLO_ITERATE
: 0)
1322 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
1325 /* explanation of above:
1326 a) We compute Ekin at the full time step
1327 if 1) we are using the AveVel Ekin, and it's not the
1328 initial step, or 2) if we are using AveEkin, but need the full
1329 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1330 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1331 EkinAveVel because it's needed for the pressure */
1333 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1338 m_add(force_vir
, shake_vir
, total_vir
); /* we need the un-dispersion corrected total vir here */
1339 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ2
);
1346 /* We need the kinetic energy at minus the half step for determining
1347 * the full step kinetic energy and possibly for T-coupling.*/
1348 /* This may not be quite working correctly yet . . . . */
1349 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, state_global
, mdatoms
, nrnb
, vcm
,
1350 wcycle
, enerd
, NULL
, NULL
, NULL
, NULL
, mu_tot
,
1351 constr
, NULL
, FALSE
, state
->box
,
1352 top_global
, &pcurr
, top_global
->natoms
, &bSumEkinhOld
,
1353 CGLO_RERUNMD
| CGLO_GSTAT
| CGLO_TEMPERATURE
);
1358 if (iterate
.bIterationActive
&&
1359 done_iterating(cr
, fplog
, step
, &iterate
, bFirstIterate
,
1360 state
->veta
, &vetanew
))
1364 bFirstIterate
= FALSE
;
1367 if (bTrotter
&& !bInitStep
)
1369 copy_mat(shake_vir
, state
->svir_prev
);
1370 copy_mat(force_vir
, state
->fvir_prev
);
1371 if (IR_NVT_TROTTER(ir
) && ir
->eI
== eiVV
)
1373 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1374 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
), ekind
, NULL
, (ir
->eI
== eiVV
), FALSE
, FALSE
);
1375 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
1378 /* if it's the initial step, we performed this first step just to get the constraint virial */
1379 if (bInitStep
&& ir
->eI
== eiVV
)
1381 copy_rvecn(cbuf
, state
->v
, 0, state
->natoms
);
1384 GMX_MPE_LOG(ev_timestep1
);
1387 /* MRS -- now done iterating -- compute the conserved quantity */
1390 saved_conserved_quantity
= compute_conserved_from_auxiliary(ir
, state
, &MassQ
);
1393 last_ekin
= enerd
->term
[F_EKIN
];
1395 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
1397 saved_conserved_quantity
-= enerd
->term
[F_DISPCORR
];
1399 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1402 sum_dhdl(enerd
, state
->lambda
, ir
->fepvals
);
1406 /* ######## END FIRST UPDATE STEP ############## */
1407 /* ######## If doing VV, we now have v(dt) ###### */
1410 /* perform extended ensemble sampling in lambda - we don't
1411 actually move to the new state before outputting
1412 statistics, but if performing simulated tempering, we
1413 do update the velocities and the tau_t. */
1415 lamnew
= ExpandedEnsembleDynamics(fplog
, ir
, enerd
, state
, &MassQ
, state
->fep_state
, &state
->dfhist
, step
, mcrng
, state
->v
, mdatoms
);
1416 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1417 copy_df_history(&state_global
->dfhist
,&state
->dfhist
);
1419 /* ################## START TRAJECTORY OUTPUT ################# */
1421 /* Now we have the energies and forces corresponding to the
1422 * coordinates at time t. We must output all of this before
1424 * for RerunMD t is read from input trajectory
1426 GMX_MPE_LOG(ev_output_start
);
1429 if (do_per_step(step
, ir
->nstxout
))
1431 mdof_flags
|= MDOF_X
;
1433 if (do_per_step(step
, ir
->nstvout
))
1435 mdof_flags
|= MDOF_V
;
1437 if (do_per_step(step
, ir
->nstfout
))
1439 mdof_flags
|= MDOF_F
;
1441 if (do_per_step(step
, ir
->nstxtcout
))
1443 mdof_flags
|= MDOF_XTC
;
1447 mdof_flags
|= MDOF_CPT
;
1451 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1454 /* Enforce writing positions and velocities at end of run */
1455 mdof_flags
|= (MDOF_X
| MDOF_V
);
1461 fcReportProgress( ir
->nsteps
, step
);
1464 #if defined(__native_client__)
1465 fcCheckin(MASTER(cr
));
1468 /* sync bCPT and fc record-keeping */
1469 if (bCPT
&& MASTER(cr
))
1471 fcRequestCheckPoint();
1475 if (mdof_flags
!= 0)
1477 wallcycle_start(wcycle
, ewcTRAJ
);
1480 if (state
->flags
& (1<<estLD_RNG
))
1482 get_stochd_state(upd
, state
);
1484 if (state
->flags
& (1<<estMC_RNG
))
1486 get_mc_state(mcrng
, state
);
1492 state_global
->ekinstate
.bUpToDate
= FALSE
;
1496 update_ekinstate(&state_global
->ekinstate
, ekind
);
1497 state_global
->ekinstate
.bUpToDate
= TRUE
;
1499 update_energyhistory(&state_global
->enerhist
, mdebin
);
1502 write_traj(fplog
, cr
, outf
, mdof_flags
, top_global
,
1503 step
, t
, state
, state_global
, f
, f_global
, &n_xtc
, &x_xtc
);
1510 if (bLastStep
&& step_rel
== ir
->nsteps
&&
1511 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
1512 !bRerunMD
&& !bFFscan
)
1514 /* x and v have been collected in write_traj,
1515 * because a checkpoint file will always be written
1518 fprintf(stderr
, "\nWriting final coordinates.\n");
1521 /* Make molecules whole only for confout writing */
1522 do_pbc_mtop(fplog
, ir
->ePBC
, state
->box
, top_global
, state_global
->x
);
1524 write_sto_conf_mtop(ftp2fn(efSTO
, nfile
, fnm
),
1525 *top_global
->name
, top_global
,
1526 state_global
->x
, state_global
->v
,
1527 ir
->ePBC
, state
->box
);
1530 wallcycle_stop(wcycle
, ewcTRAJ
);
1532 GMX_MPE_LOG(ev_output_finish
);
1534 /* kludge -- virial is lost with restart for NPT control. Must restart */
1535 if (bStartingFromCpt
&& bVV
)
1537 copy_mat(state
->svir_prev
, shake_vir
);
1538 copy_mat(state
->fvir_prev
, force_vir
);
1540 /* ################## END TRAJECTORY OUTPUT ################ */
1542 /* Determine the wallclock run time up till now */
1543 run_time
= gmx_gettime() - (double)runtime
->real
;
1545 /* Check whether everything is still allright */
1546 if (((int)gmx_get_stop_condition() > handled_stop_condition
)
1547 #ifdef GMX_THREAD_MPI
1552 /* this is just make gs.sig compatible with the hack
1553 of sending signals around by MPI_Reduce with together with
1555 if (gmx_get_stop_condition() == gmx_stop_cond_next_ns
)
1557 gs
.sig
[eglsSTOPCOND
] = 1;
1559 if (gmx_get_stop_condition() == gmx_stop_cond_next
)
1561 gs
.sig
[eglsSTOPCOND
] = -1;
1563 /* < 0 means stop at next step, > 0 means stop at next NS step */
1567 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1568 gmx_get_signal_name(),
1569 gs
.sig
[eglsSTOPCOND
] == 1 ? "NS " : "");
1573 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1574 gmx_get_signal_name(),
1575 gs
.sig
[eglsSTOPCOND
] == 1 ? "NS " : "");
1577 handled_stop_condition
= (int)gmx_get_stop_condition();
1579 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
1580 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
1581 gs
.sig
[eglsSTOPCOND
] == 0 && gs
.set
[eglsSTOPCOND
] == 0)
1583 /* Signal to terminate the run */
1584 gs
.sig
[eglsSTOPCOND
] = 1;
1587 fprintf(fplog
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step
, sbuf
), max_hours
*0.99);
1589 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step
, sbuf
), max_hours
*0.99);
1592 if (bResetCountersHalfMaxH
&& MASTER(cr
) &&
1593 run_time
> max_hours
*60.0*60.0*0.495)
1595 gs
.sig
[eglsRESETCOUNTERS
] = 1;
1598 if (ir
->nstlist
== -1 && !bRerunMD
)
1600 /* When bGStatEveryStep=FALSE, global_stat is only called
1601 * when we check the atom displacements, not at NS steps.
1602 * This means that also the bonded interaction count check is not
1603 * performed immediately after NS. Therefore a few MD steps could
1604 * be performed with missing interactions.
1605 * But wrong energies are never written to file,
1606 * since energies are only written after global_stat
1609 if (step
>= nlh
.step_nscheck
)
1611 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
, fr
, &top
->cgs
,
1612 nlh
.scale_tot
, state
->x
);
1616 /* This is not necessarily true,
1617 * but step_nscheck is determined quite conservatively.
1623 /* In parallel we only have to check for checkpointing in steps
1624 * where we do global communication,
1625 * otherwise the other nodes don't know.
1627 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
1630 run_time
>= nchkpt
*cpt_period
*60.0)) &&
1631 gs
.set
[eglsCHKPT
] == 0)
1633 gs
.sig
[eglsCHKPT
] = 1;
1636 /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */
1641 update_tcouple(fplog
, step
, ir
, state
, ekind
, wcycle
, upd
, &MassQ
, mdatoms
);
1643 if (ETC_ANDERSEN(ir
->etc
)) /* keep this outside of update_tcouple because of the extra info required to pass */
1645 gmx_bool bIfRandomize
;
1646 bIfRandomize
= update_randomize_velocities(ir
, step
, mdatoms
, state
, upd
, &top
->idef
, constr
);
1647 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1648 if (constr
&& bIfRandomize
)
1650 update_constraints(fplog
, step
, NULL
, ir
, ekind
, mdatoms
,
1651 state
, fr
->bMolPBC
, graph
, f
,
1652 &top
->idef
, tmp_vir
, NULL
,
1653 cr
, nrnb
, wcycle
, upd
, constr
,
1654 bInitStep
, TRUE
, bCalcVir
, vetanew
);
1659 if (bIterativeCase
&& do_per_step(step
, ir
->nstpcouple
))
1661 gmx_iterate_init(&iterate
, TRUE
);
1662 /* for iterations, we save these vectors, as we will be redoing the calculations */
1663 copy_coupling_state(state
, bufstate
, ekind
, ekind_save
, &(ir
->opts
));
1666 bFirstIterate
= TRUE
;
1667 while (bFirstIterate
|| iterate
.bIterationActive
)
1669 /* We now restore these vectors to redo the calculation with improved extended variables */
1670 if (iterate
.bIterationActive
)
1672 copy_coupling_state(bufstate
, state
, ekind_save
, ekind
, &(ir
->opts
));
1675 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1676 so scroll down for that logic */
1678 /* ######### START SECOND UPDATE STEP ################# */
1679 GMX_MPE_LOG(ev_update_start
);
1680 /* Box is changed in update() when we do pressure coupling,
1681 * but we should still use the old box for energy corrections and when
1682 * writing it to the energy file, so it matches the trajectory files for
1683 * the same timestep above. Make a copy in a separate array.
1685 copy_mat(state
->box
, lastbox
);
1690 if (!(bRerunMD
&& !rerun_fr
.bV
&& !bForceUpdate
))
1692 wallcycle_start(wcycle
, ewcUPDATE
);
1693 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1696 if (iterate
.bIterationActive
)
1704 /* we use a new value of scalevir to converge the iterations faster */
1705 scalevir
= tracevir
/trace(shake_vir
);
1707 msmul(shake_vir
, scalevir
, shake_vir
);
1708 m_add(force_vir
, shake_vir
, total_vir
);
1709 clear_mat(shake_vir
);
1711 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ3
);
1712 /* We can only do Berendsen coupling after we have summed
1713 * the kinetic energy or virial. Since the happens
1714 * in global_state after update, we should only do it at
1715 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1720 update_tcouple(fplog
, step
, ir
, state
, ekind
, wcycle
, upd
, &MassQ
, mdatoms
);
1721 update_pcouple(fplog
, step
, ir
, state
, pcoupl_mu
, M
, wcycle
,
1727 bUpdateDoLR
= (fr
->bTwinRange
&& do_per_step(step
, ir
->nstcalclr
));
1729 /* velocity half-step update */
1730 update_coords(fplog
, step
, ir
, mdatoms
, state
, fr
->bMolPBC
, f
,
1731 bUpdateDoLR
, fr
->f_twin
, fcd
,
1732 ekind
, M
, wcycle
, upd
, FALSE
, etrtVELOCITY2
,
1733 cr
, nrnb
, constr
, &top
->idef
);
1736 /* Above, initialize just copies ekinh into ekin,
1737 * it doesn't copy position (for VV),
1738 * and entire integrator for MD.
1741 if (ir
->eI
== eiVVAK
)
1743 copy_rvecn(state
->x
, cbuf
, 0, state
->natoms
);
1745 bUpdateDoLR
= (fr
->bTwinRange
&& do_per_step(step
, ir
->nstcalclr
));
1747 update_coords(fplog
, step
, ir
, mdatoms
, state
, fr
->bMolPBC
, f
,
1748 bUpdateDoLR
, fr
->f_twin
, fcd
,
1749 ekind
, M
, wcycle
, upd
, bInitStep
, etrtPOSITION
, cr
, nrnb
, constr
, &top
->idef
);
1750 wallcycle_stop(wcycle
, ewcUPDATE
);
1752 update_constraints(fplog
, step
, &dvdl_constr
, ir
, ekind
, mdatoms
, state
,
1753 fr
->bMolPBC
, graph
, f
,
1754 &top
->idef
, shake_vir
, force_vir
,
1755 cr
, nrnb
, wcycle
, upd
, constr
,
1756 bInitStep
, FALSE
, bCalcVir
, state
->veta
);
1758 if (ir
->eI
== eiVVAK
)
1760 /* erase F_EKIN and F_TEMP here? */
1761 /* just compute the kinetic energy at the half step to perform a trotter step */
1762 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, state_global
, mdatoms
, nrnb
, vcm
,
1763 wcycle
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
1764 constr
, NULL
, FALSE
, lastbox
,
1765 top_global
, &pcurr
, top_global
->natoms
, &bSumEkinhOld
,
1766 cglo_flags
| CGLO_TEMPERATURE
1768 wallcycle_start(wcycle
, ewcUPDATE
);
1769 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ4
);
1770 /* now we know the scaling, we can compute the positions again again */
1771 copy_rvecn(cbuf
, state
->x
, 0, state
->natoms
);
1773 bUpdateDoLR
= (fr
->bTwinRange
&& do_per_step(step
, ir
->nstcalclr
));
1775 update_coords(fplog
, step
, ir
, mdatoms
, state
, fr
->bMolPBC
, f
,
1776 bUpdateDoLR
, fr
->f_twin
, fcd
,
1777 ekind
, M
, wcycle
, upd
, bInitStep
, etrtPOSITION
, cr
, nrnb
, constr
, &top
->idef
);
1778 wallcycle_stop(wcycle
, ewcUPDATE
);
1780 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1781 /* are the small terms in the shake_vir here due
1782 * to numerical errors, or are they important
1783 * physically? I'm thinking they are just errors, but not completely sure.
1784 * For now, will call without actually constraining, constr=NULL*/
1785 update_constraints(fplog
, step
, NULL
, ir
, ekind
, mdatoms
,
1786 state
, fr
->bMolPBC
, graph
, f
,
1787 &top
->idef
, tmp_vir
, force_vir
,
1788 cr
, nrnb
, wcycle
, upd
, NULL
,
1789 bInitStep
, FALSE
, bCalcVir
,
1792 if (!bOK
&& !bFFscan
)
1794 gmx_fatal(FARGS
, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
1797 if (fr
->bSepDVDL
&& fplog
&& do_log
)
1799 fprintf(fplog
, sepdvdlformat
, "Constraint dV/dl", 0.0, dvdl_constr
);
1803 /* this factor or 2 correction is necessary
1804 because half of the constraint force is removed
1805 in the vv step, so we have to double it. See
1806 the Redmine issue #1255. It is not yet clear
1807 if the factor of 2 is exact, or just a very
1808 good approximation, and this will be
1809 investigated. The next step is to see if this
1810 can be done adding a dhdl contribution from the
1811 rattle step, but this is somewhat more
1812 complicated with the current code. Will be
1813 investigated, hopefully for 4.6.3. However,
1814 this current solution is much better than
1815 having it completely wrong.
1817 enerd
->term
[F_DVDL_CONSTR
] += 2*dvdl_constr
;
1821 enerd
->term
[F_DVDL_CONSTR
] += dvdl_constr
;
1826 /* Need to unshift here */
1827 unshift_self(graph
, state
->box
, state
->x
);
1830 GMX_BARRIER(cr
->mpi_comm_mygroup
);
1831 GMX_MPE_LOG(ev_update_finish
);
1835 wallcycle_start(wcycle
, ewcVSITECONSTR
);
1838 shift_self(graph
, state
->box
, state
->x
);
1840 construct_vsites(fplog
, vsite
, state
->x
, nrnb
, ir
->delta_t
, state
->v
,
1841 top
->idef
.iparams
, top
->idef
.il
,
1842 fr
->ePBC
, fr
->bMolPBC
, graph
, cr
, state
->box
);
1846 unshift_self(graph
, state
->box
, state
->x
);
1848 wallcycle_stop(wcycle
, ewcVSITECONSTR
);
1851 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1852 /* With Leap-Frog we can skip compute_globals at
1853 * non-communication steps, but we need to calculate
1854 * the kinetic energy one step before communication.
1856 if (bGStat
|| (!EI_VV(ir
->eI
) && do_per_step(step
+1, nstglobalcomm
)))
1858 if (ir
->nstlist
== -1 && bFirstIterate
)
1860 gs
.sig
[eglsNABNSB
] = nlh
.nabnsb
;
1862 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, state_global
, mdatoms
, nrnb
, vcm
,
1863 wcycle
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
1865 bFirstIterate
? &gs
: NULL
,
1866 (step_rel
% gs
.nstms
== 0) &&
1867 (multisim_nsteps
< 0 || (step_rel
< multisim_nsteps
)),
1869 top_global
, &pcurr
, top_global
->natoms
, &bSumEkinhOld
,
1871 | (!EI_VV(ir
->eI
) || bRerunMD
? CGLO_ENERGY
: 0)
1872 | (!EI_VV(ir
->eI
) && bStopCM
? CGLO_STOPCM
: 0)
1873 | (!EI_VV(ir
->eI
) ? CGLO_TEMPERATURE
: 0)
1874 | (!EI_VV(ir
->eI
) || bRerunMD
? CGLO_PRESSURE
: 0)
1875 | (iterate
.bIterationActive
? CGLO_ITERATE
: 0)
1876 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
1879 if (ir
->nstlist
== -1 && bFirstIterate
)
1881 nlh
.nabnsb
= gs
.set
[eglsNABNSB
];
1882 gs
.set
[eglsNABNSB
] = 0;
1885 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1886 /* ############# END CALC EKIN AND PRESSURE ################# */
1888 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1889 the virial that should probably be addressed eventually. state->veta has better properies,
1890 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1891 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1893 if (iterate
.bIterationActive
&&
1894 done_iterating(cr
, fplog
, step
, &iterate
, bFirstIterate
,
1895 trace(shake_vir
), &tracevir
))
1899 bFirstIterate
= FALSE
;
1902 if (!bVV
|| bRerunMD
)
1904 /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
1905 sum_dhdl(enerd
, state
->lambda
, ir
->fepvals
);
1907 update_box(fplog
, step
, ir
, mdatoms
, state
, graph
, f
,
1908 ir
->nstlist
== -1 ? &nlh
.scale_tot
: NULL
, pcoupl_mu
, nrnb
, wcycle
, upd
, bInitStep
, FALSE
);
1910 /* ################# END UPDATE STEP 2 ################# */
1911 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1913 /* The coordinates (x) were unshifted in update */
1914 if (bFFscan
&& (shellfc
== NULL
|| bConverged
))
1916 if (print_forcefield(fplog
, enerd
->term
, mdatoms
->homenr
,
1918 &(top_global
->mols
), mdatoms
->massT
, pres
))
1922 fprintf(stderr
, "\n");
1928 /* We will not sum ekinh_old,
1929 * so signal that we still have to do it.
1931 bSumEkinhOld
= TRUE
;
1936 /* Only do GCT when the relaxation of shells (minimization) has converged,
1937 * otherwise we might be coupling to bogus energies.
1938 * In parallel we must always do this, because the other sims might
1942 /* Since this is called with the new coordinates state->x, I assume
1943 * we want the new box state->box too. / EL 20040121
1945 do_coupling(fplog
, oenv
, nfile
, fnm
, tcr
, t
, step
, enerd
->term
, fr
,
1947 mdatoms
, &(top
->idef
), mu_aver
,
1948 top_global
->mols
.nr
, cr
,
1949 state
->box
, total_vir
, pres
,
1950 mu_tot
, state
->x
, f
, bConverged
);
1954 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1956 /* use the directly determined last velocity, not actually the averaged half steps */
1957 if (bTrotter
&& ir
->eI
== eiVV
)
1959 enerd
->term
[F_EKIN
] = last_ekin
;
1961 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
1965 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + saved_conserved_quantity
;
1969 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + compute_conserved_from_auxiliary(ir
, state
, &MassQ
);
1971 /* Check for excessively large energies */
1975 real etot_max
= 1e200
;
1977 real etot_max
= 1e30
;
1979 if (fabs(enerd
->term
[F_ETOT
]) > etot_max
)
1981 fprintf(stderr
, "Energy too large (%g), giving up\n",
1982 enerd
->term
[F_ETOT
]);
1985 /* ######### END PREPARING EDR OUTPUT ########### */
1987 /* Time for performance */
1988 if (((step
% stepout
) == 0) || bLastStep
)
1990 runtime_upd_proc(runtime
);
1996 gmx_bool do_dr
, do_or
;
1998 if (fplog
&& do_log
&& bDoExpanded
)
2000 /* only needed if doing expanded ensemble */
2001 PrintFreeEnergyInfoToFile(fplog
, ir
->fepvals
, ir
->expandedvals
, ir
->bSimTemp
? ir
->simtempvals
: NULL
,
2002 &state_global
->dfhist
, state
->fep_state
, ir
->nstlog
, step
);
2004 if (!(bStartingFromCpt
&& (EI_VV(ir
->eI
))))
2008 upd_mdebin(mdebin
, bDoDHDL
, TRUE
,
2009 t
, mdatoms
->tmass
, enerd
, state
,
2010 ir
->fepvals
, ir
->expandedvals
, lastbox
,
2011 shake_vir
, force_vir
, total_vir
, pres
,
2012 ekind
, mu_tot
, constr
);
2016 upd_mdebin_step(mdebin
);
2019 do_dr
= do_per_step(step
, ir
->nstdisreout
);
2020 do_or
= do_per_step(step
, ir
->nstorireout
);
2022 print_ebin(outf
->fp_ene
, do_ene
, do_dr
, do_or
, do_log
? fplog
: NULL
,
2024 eprNORMAL
, bCompact
, mdebin
, fcd
, groups
, &(ir
->opts
));
2026 if (ir
->ePull
!= epullNO
)
2028 pull_print_output(ir
->pull
, step
, t
);
2031 if (do_per_step(step
, ir
->nstlog
))
2033 if (fflush(fplog
) != 0)
2035 gmx_fatal(FARGS
, "Cannot flush logfile - maybe you are out of disk space?");
2041 /* Have to do this part _after_ outputting the logfile and the edr file */
2042 /* Gets written into the state at the beginning of next loop*/
2043 state
->fep_state
= lamnew
;
2046 /* Remaining runtime */
2047 if (MULTIMASTER(cr
) && (do_verbose
|| gmx_got_usr_signal()) && !bPMETuneRunning
)
2051 fprintf(stderr
, "\n");
2053 print_time(stderr
, runtime
, step
, ir
, cr
);
2056 /* Replica exchange */
2058 if ((repl_ex_nst
> 0) && (step
> 0) && !bLastStep
&&
2059 do_per_step(step
, repl_ex_nst
))
2061 bExchanged
= replica_exchange(fplog
, cr
, repl_ex
,
2062 state_global
, enerd
,
2065 if (bExchanged
&& DOMAINDECOMP(cr
))
2067 dd_partition_system(fplog
, step
, cr
, TRUE
, 1,
2068 state_global
, top_global
, ir
,
2069 state
, &f
, mdatoms
, top
, fr
,
2070 vsite
, shellfc
, constr
,
2071 nrnb
, wcycle
, FALSE
);
2077 bStartingFromCpt
= FALSE
;
2079 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2080 /* With all integrators, except VV, we need to retain the pressure
2081 * at the current step for coupling at the next step.
2083 if ((state
->flags
& (1<<estPRES_PREV
)) &&
2085 (ir
->nstpcouple
> 0 && step
% ir
->nstpcouple
== 0)))
2087 /* Store the pressure in t_state for pressure coupling
2088 * at the next MD step.
2090 copy_mat(pres
, state
->pres_prev
);
2093 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2095 if ( (membed
!= NULL
) && (!bLastStep
) )
2097 rescale_membed(step_rel
, membed
, state_global
->x
);
2104 /* read next frame from input trajectory */
2105 bNotLastFrame
= read_next_frame(oenv
, status
, &rerun_fr
);
2110 rerun_parallel_comm(cr
, &rerun_fr
, &bNotLastFrame
);
2114 if (!bRerunMD
|| !rerun_fr
.bStep
)
2116 /* increase the MD step number */
2121 cycles
= wallcycle_stop(wcycle
, ewcSTEP
);
2122 if (DOMAINDECOMP(cr
) && wcycle
)
2124 dd_cycles_add(cr
->dd
, cycles
, ddCyclStep
);
2127 if (bPMETuneRunning
|| bPMETuneTry
)
2129 /* PME grid + cut-off optimization with GPUs or PME nodes */
2131 /* Count the total cycles over the last steps */
2132 cycles_pmes
+= cycles
;
2134 /* We can only switch cut-off at NS steps */
2135 if (step
% ir
->nstlist
== 0)
2137 /* PME grid + cut-off optimization with GPUs or PME nodes */
2140 if (DDMASTER(cr
->dd
))
2142 /* PME node load is too high, start tuning */
2143 bPMETuneRunning
= (dd_pme_f_ratio(cr
->dd
) >= 1.05);
2145 dd_bcast(cr
->dd
, sizeof(gmx_bool
), &bPMETuneRunning
);
2147 if (bPMETuneRunning
|| step_rel
> ir
->nstlist
*50)
2149 bPMETuneTry
= FALSE
;
2152 if (bPMETuneRunning
)
2154 /* init_step might not be a multiple of nstlist,
2155 * but the first cycle is always skipped anyhow.
2158 pme_load_balance(pme_loadbal
, cr
,
2159 (bVerbose
&& MASTER(cr
)) ? stderr
: NULL
,
2161 ir
, state
, cycles_pmes
,
2162 fr
->ic
, fr
->nbv
, &fr
->pmedata
,
2165 /* Update constants in forcerec/inputrec to keep them in sync with fr->ic */
2166 fr
->ewaldcoeff
= fr
->ic
->ewaldcoeff
;
2167 fr
->rlist
= fr
->ic
->rlist
;
2168 fr
->rlistlong
= fr
->ic
->rlistlong
;
2169 fr
->rcoulomb
= fr
->ic
->rcoulomb
;
2170 fr
->rvdw
= fr
->ic
->rvdw
;
2176 if (step_rel
== wcycle_get_reset_counters(wcycle
) ||
2177 gs
.set
[eglsRESETCOUNTERS
] != 0)
2179 /* Reset all the counters related to performance over the run */
2180 reset_all_counters(fplog
, cr
, step
, &step_rel
, ir
, wcycle
, nrnb
, runtime
,
2181 fr
->nbv
!= NULL
&& fr
->nbv
->bUseGPU
? fr
->nbv
->cu_nbv
: NULL
);
2182 wcycle_set_reset_counters(wcycle
, -1);
2183 if (!(cr
->duty
& DUTY_PME
))
2185 /* Tell our PME node to reset its counters */
2186 gmx_pme_send_resetcounters(cr
, step
);
2188 /* Correct max_hours for the elapsed time */
2189 max_hours
-= run_time
/(60.0*60.0);
2190 bResetCountersHalfMaxH
= FALSE
;
2191 gs
.set
[eglsRESETCOUNTERS
] = 0;
2195 /* End of main MD loop */
2199 runtime_end(runtime
);
2201 if (bRerunMD
&& MASTER(cr
))
2206 if (!(cr
->duty
& DUTY_PME
))
2208 /* Tell the PME only node to finish */
2209 gmx_pme_send_finish(cr
);
2214 if (ir
->nstcalcenergy
> 0 && !bRerunMD
)
2216 print_ebin(outf
->fp_ene
, FALSE
, FALSE
, FALSE
, fplog
, step
, t
,
2217 eprAVER
, FALSE
, mdebin
, fcd
, groups
, &(ir
->opts
));
2225 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
2227 fprintf(fplog
, "Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n", nlh
.s1
/nlh
.nns
, sqrt(nlh
.s2
/nlh
.nns
- sqr(nlh
.s1
/nlh
.nns
)));
2228 fprintf(fplog
, "Average number of atoms that crossed the half buffer length: %.1f\n\n", nlh
.ab
/nlh
.nns
);
2231 if (pme_loadbal
!= NULL
)
2233 pme_loadbal_done(pme_loadbal
, cr
, fplog
,
2234 fr
->nbv
!= NULL
&& fr
->nbv
->bUseGPU
);
2237 if (shellfc
&& fplog
)
2239 fprintf(fplog
, "Fraction of iterations that converged: %.2f %%\n",
2240 (nconverged
*100.0)/step_rel
);
2241 fprintf(fplog
, "Average number of force evaluations per MD step: %.2f\n\n",
2245 if (repl_ex_nst
> 0 && MASTER(cr
))
2247 print_replica_exchange_statistics(fplog
, repl_ex
);
2250 runtime
->nsteps_done
= step_rel
;