gmx_edsam_t ed, t_forcerec *fr,
int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t membed,
real cpt_period, real max_hours,
- const char gmx_unused *deviceOptions,
int imdport,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting)
t_graph *graph = NULL;
gmx_signalling_t gs;
gmx_groups_t *groups;
- gmx_ekindata_t *ekind, *ekind_save;
+ gmx_ekindata_t *ekind;
gmx_shellfc_t shellfc;
int count, nconverged = 0;
double tcount = 0;
nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
bGStatEveryStep = (nstglobalcomm == 1);
- if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
- {
- fprintf(fplog,
- "To reduce the energy communication with nstlist = -1\n"
- "the neighbor list validity should not be checked at every step,\n"
- "this means that exact integration is not guaranteed.\n"
- "The neighbor list validity is checked after:\n"
- " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
- "In most cases this will result in exact integration.\n"
- "This reduces the energy communication by a factor of 2 to 3.\n"
- "If you want less energy communication, set nstlist > 3.\n\n");
- }
-
if (bRerunMD)
{
ir->nstxout_compressed = 0;
/* Kinetic energy data */
snew(ekind, 1);
init_ekindata(fplog, top_global, &(ir->opts), ekind);
- /* needed for iteration of constraints */
- snew(ekind_save, 1);
- init_ekindata(fplog, top_global, &(ir->opts), ekind_save);
/* Copy the cos acceleration to the groups struct */
ekind->cosacc.cos_accel = ir->cos_accel;
nrnb, wcycle,
do_verbose && !bPMETuneRunning);
wallcycle_stop(wcycle, ewcDOMDEC);
- /* If using an iterative integrator, reallocate space to match the decomposition */
}
}
ekind, M, upd, bInitStep, etrtVELOCITY1,
cr, nrnb, constr, &top->idef);
- /* TODO remove the brace below, once iteration is removed */
+ if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
{
- if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
- {
- wallcycle_stop(wcycle, ewcUPDATE);
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
- &top->idef, shake_vir,
- cr, nrnb, wcycle, upd, constr,
- TRUE, bCalcVir);
- wallcycle_start(wcycle, ewcUPDATE);
- if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
- {
- /* Correct the virial for multiple time stepping */
- m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
- }
- }
- else if (graph)
- {
- /* Need to unshift here if a do_force has been
- called in the previous step */
- unshift_self(graph, state->box, state->x);
- }
- /* if VV, compute the pressure and constraints */
- /* For VV2, we strictly only need this if using pressure
- * control, but we really would like to have accurate pressures
- * printed out.
- * Think about ways around this in the future?
- * For now, keep this choice in comments.
- */
- /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
- /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
- bPres = TRUE;
- bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
- if (bCalcEner && ir->eI == eiVVAK)
+ wallcycle_stop(wcycle, ewcUPDATE);
+ update_constraints(fplog, step, NULL, ir, mdatoms,
+ state, fr->bMolPBC, graph, f,
+ &top->idef, shake_vir,
+ cr, nrnb, wcycle, upd, constr,
+ TRUE, bCalcVir);
+ wallcycle_start(wcycle, ewcUPDATE);
+ if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
{
- bSumEkinhOld = TRUE;
+ /* Correct the virial for multiple time stepping */
+ m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
}
- /* for vv, the first half of the integration actually corresponds to the previous step.
- So we need information from the last step in the first half of the integration */
- if (bGStat || do_per_step(step-1, nstglobalcomm))
+ }
+ else if (graph)
+ {
+ /* Need to unshift here if a do_force has been
+ called in the previous step */
+ unshift_self(graph, state->box, state->x);
+ }
+ /* if VV, compute the pressure and constraints */
+ /* For VV2, we strictly only need this if using pressure
+ * control, but we really would like to have accurate pressures
+ * printed out.
+ * Think about ways around this in the future?
+ * For now, keep this choice in comments.
+ */
+ /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
+ /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
+ bPres = TRUE;
+ bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
+ if (bCalcEner && ir->eI == eiVVAK)
+ {
+ bSumEkinhOld = TRUE;
+ }
+ /* for vv, the first half of the integration actually corresponds to the previous step.
+ So we need information from the last step in the first half of the integration */
+ if (bGStat || do_per_step(step-1, nstglobalcomm))
+ {
+ wallcycle_stop(wcycle, ewcUPDATE);
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, NULL, FALSE, state->box,
+ top_global, &bSumEkinhOld,
+ cglo_flags
+ | CGLO_ENERGY
+ | (bTemp ? CGLO_TEMPERATURE : 0)
+ | (bPres ? CGLO_PRESSURE : 0)
+ | (bPres ? CGLO_CONSTRAINT : 0)
+ | CGLO_SCALEEKIN
+ );
+ /* explanation of above:
+ a) We compute Ekin at the full time step
+ if 1) we are using the AveVel Ekin, and it's not the
+ initial step, or 2) if we are using AveEkin, but need the full
+ time step kinetic energy for the pressure (always true now, since we want accurate statistics).
+ b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
+ EkinAveVel because it's needed for the pressure */
+ wallcycle_start(wcycle, ewcUPDATE);
+ }
+ /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
+ if (!bInitStep)
+ {
+ if (bTrotter)
{
- wallcycle_stop(wcycle, ewcUPDATE);
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
- wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, NULL, FALSE, state->box,
- top_global, &bSumEkinhOld,
- cglo_flags
- | CGLO_ENERGY
- | (bTemp ? CGLO_TEMPERATURE : 0)
- | (bPres ? CGLO_PRESSURE : 0)
- | (bPres ? CGLO_CONSTRAINT : 0)
- | CGLO_SCALEEKIN
- );
- /* explanation of above:
- a) We compute Ekin at the full time step
- if 1) we are using the AveVel Ekin, and it's not the
- initial step, or 2) if we are using AveEkin, but need the full
- time step kinetic energy for the pressure (always true now, since we want accurate statistics).
- b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
- EkinAveVel because it's needed for the pressure */
- wallcycle_start(wcycle, ewcUPDATE);
+ m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
}
- /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
- if (!bInitStep)
+ else
{
- if (bTrotter)
- {
- m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
- }
- else
+ if (bExchanged)
{
- if (bExchanged)
- {
- wallcycle_stop(wcycle, ewcUPDATE);
- /* We need the kinetic energy at minus the half step for determining
- * the full step kinetic energy and possibly for T-coupling.*/
- /* This may not be quite working correctly yet . . . . */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
- wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
- constr, NULL, FALSE, state->box,
- top_global, &bSumEkinhOld,
- CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
- wallcycle_start(wcycle, ewcUPDATE);
- }
+ wallcycle_stop(wcycle, ewcUPDATE);
+ /* We need the kinetic energy at minus the half step for determining
+ * the full step kinetic energy and possibly for T-coupling.*/
+ /* This may not be quite working correctly yet . . . . */
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ constr, NULL, FALSE, state->box,
+ top_global, &bSumEkinhOld,
+ CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
+ wallcycle_start(wcycle, ewcUPDATE);
}
}
}
- /* TODO remove the brace above, once iteration is removed */
if (bTrotter && !bInitStep)
{
copy_mat(shake_vir, state->svir_prev);
}
}
}
- /* TODO remove the brace below, once iteration is removed */
+ /* ######### START SECOND UPDATE STEP ################# */
+ /* Box is changed in update() when we do pressure coupling,
+ * but we should still use the old box for energy corrections and when
+ * writing it to the energy file, so it matches the trajectory files for
+ * the same timestep above. Make a copy in a separate array.
+ */
+ copy_mat(state->box, lastbox);
+
+ dvdl_constr = 0;
+
+ if (!bRerunMD || rerun_fr.bV || bForceUpdate)
{
- /* ######### START SECOND UPDATE STEP ################# */
- /* Box is changed in update() when we do pressure coupling,
- * but we should still use the old box for energy corrections and when
- * writing it to the energy file, so it matches the trajectory files for
- * the same timestep above. Make a copy in a separate array.
- */
- copy_mat(state->box, lastbox);
+ wallcycle_start(wcycle, ewcUPDATE);
+ /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
+ if (bTrotter)
+ {
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
+ /* We can only do Berendsen coupling after we have summed
+ * the kinetic energy or virial. Since the happens
+ * in global_state after update, we should only do it at
+ * step % nstlist = 1 with bGStatEveryStep=FALSE.
+ */
+ }
+ else
+ {
+ update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
+ update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
+ }
- dvdl_constr = 0;
+ if (bVV)
+ {
+ bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
+
+ /* velocity half-step update */
+ update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
+ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
+ ekind, M, upd, FALSE, etrtVELOCITY2,
+ cr, nrnb, constr, &top->idef);
+ }
+
+ /* Above, initialize just copies ekinh into ekin,
+ * it doesn't copy position (for VV),
+ * and entire integrator for MD.
+ */
- if (!bRerunMD || rerun_fr.bV || bForceUpdate)
+ if (ir->eI == eiVVAK)
{
- wallcycle_start(wcycle, ewcUPDATE);
- /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
- if (bTrotter)
+ /* We probably only need md->homenr, not state->natoms */
+ if (state->natoms > cbuf_nalloc)
{
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
- /* We can only do Berendsen coupling after we have summed
- * the kinetic energy or virial. Since the happens
- * in global_state after update, we should only do it at
- * step % nstlist = 1 with bGStatEveryStep=FALSE.
- */
- }
- else
- {
- update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
- update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
+ cbuf_nalloc = state->natoms;
+ srenew(cbuf, cbuf_nalloc);
}
+ copy_rvecn(state->x, cbuf, 0, state->natoms);
+ }
+ bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
- if (bVV)
- {
- bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
+ update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
+ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
+ ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
+ wallcycle_stop(wcycle, ewcUPDATE);
- /* velocity half-step update */
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
- ekind, M, upd, FALSE, etrtVELOCITY2,
- cr, nrnb, constr, &top->idef);
- }
+ update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
+ fr->bMolPBC, graph, f,
+ &top->idef, shake_vir,
+ cr, nrnb, wcycle, upd, constr,
+ FALSE, bCalcVir);
- /* Above, initialize just copies ekinh into ekin,
- * it doesn't copy position (for VV),
- * and entire integrator for MD.
- */
+ if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
+ {
+ /* Correct the virial for multiple time stepping */
+ m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
+ }
+
+ if (ir->eI == eiVVAK)
+ {
+ /* erase F_EKIN and F_TEMP here? */
+ /* just compute the kinetic energy at the half step to perform a trotter step */
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, NULL, FALSE, lastbox,
+ top_global, &bSumEkinhOld,
+ cglo_flags | CGLO_TEMPERATURE
+ );
+ wallcycle_start(wcycle, ewcUPDATE);
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
+ /* now we know the scaling, we can compute the positions again again */
+ copy_rvecn(cbuf, state->x, 0, state->natoms);
- if (ir->eI == eiVVAK)
- {
- /* We probably only need md->homenr, not state->natoms */
- if (state->natoms > cbuf_nalloc)
- {
- cbuf_nalloc = state->natoms;
- srenew(cbuf, cbuf_nalloc);
- }
- copy_rvecn(state->x, cbuf, 0, state->natoms);
- }
bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
wallcycle_stop(wcycle, ewcUPDATE);
- update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
- fr->bMolPBC, graph, f,
- &top->idef, shake_vir,
- cr, nrnb, wcycle, upd, constr,
+ /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
+ /* are the small terms in the shake_vir here due
+ * to numerical errors, or are they important
+ * physically? I'm thinking they are just errors, but not completely sure.
+ * For now, will call without actually constraining, constr=NULL*/
+ update_constraints(fplog, step, NULL, ir, mdatoms,
+ state, fr->bMolPBC, graph, f,
+ &top->idef, tmp_vir,
+ cr, nrnb, wcycle, upd, NULL,
FALSE, bCalcVir);
-
- if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
- {
- /* Correct the virial for multiple time stepping */
- m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
- }
-
- if (ir->eI == eiVVAK)
- {
- /* erase F_EKIN and F_TEMP here? */
- /* just compute the kinetic energy at the half step to perform a trotter step */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
- wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, NULL, FALSE, lastbox,
- top_global, &bSumEkinhOld,
- cglo_flags | CGLO_TEMPERATURE
- );
- wallcycle_start(wcycle, ewcUPDATE);
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
- /* now we know the scaling, we can compute the positions again again */
- copy_rvecn(cbuf, state->x, 0, state->natoms);
-
- bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
-
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
- ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
- wallcycle_stop(wcycle, ewcUPDATE);
-
- /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
- /* are the small terms in the shake_vir here due
- * to numerical errors, or are they important
- * physically? I'm thinking they are just errors, but not completely sure.
- * For now, will call without actually constraining, constr=NULL*/
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
- &top->idef, tmp_vir,
- cr, nrnb, wcycle, upd, NULL,
- FALSE, bCalcVir);
- }
- if (bVV)
- {
- /* this factor or 2 correction is necessary
- because half of the constraint force is removed
- in the vv step, so we have to double it. See
- the Redmine issue #1255. It is not yet clear
- if the factor of 2 is exact, or just a very
- good approximation, and this will be
- investigated. The next step is to see if this
- can be done adding a dhdl contribution from the
- rattle step, but this is somewhat more
- complicated with the current code. Will be
- investigated, hopefully for 4.6.3. However,
- this current solution is much better than
- having it completely wrong.
- */
- enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
- }
- else
- {
- enerd->term[F_DVDL_CONSTR] += dvdl_constr;
- }
}
- else if (graph)
+ if (bVV)
{
- /* Need to unshift here */
- unshift_self(graph, state->box, state->x);
+ /* this factor or 2 correction is necessary
+ because half of the constraint force is removed
+ in the vv step, so we have to double it. See
+ the Redmine issue #1255. It is not yet clear
+ if the factor of 2 is exact, or just a very
+ good approximation, and this will be
+ investigated. The next step is to see if this
+ can be done adding a dhdl contribution from the
+ rattle step, but this is somewhat more
+ complicated with the current code. Will be
+ investigated, hopefully for 4.6.3. However,
+ this current solution is much better than
+ having it completely wrong.
+ */
+ enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
}
-
- if (vsite != NULL)
+ else
{
- wallcycle_start(wcycle, ewcVSITECONSTR);
- if (graph != NULL)
- {
- shift_self(graph, state->box, state->x);
- }
- construct_vsites(vsite, state->x, ir->delta_t, state->v,
- top->idef.iparams, top->idef.il,
- fr->ePBC, fr->bMolPBC, cr, state->box);
-
- if (graph != NULL)
- {
- unshift_self(graph, state->box, state->x);
- }
- wallcycle_stop(wcycle, ewcVSITECONSTR);
+ enerd->term[F_DVDL_CONSTR] += dvdl_constr;
}
+ }
+ else if (graph)
+ {
+ /* Need to unshift here */
+ unshift_self(graph, state->box, state->x);
+ }
- /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
- /* With Leap-Frog we can skip compute_globals at
- * non-communication steps, but we need to calculate
- * the kinetic energy one step before communication.
- */
- if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
+ if (vsite != NULL)
+ {
+ wallcycle_start(wcycle, ewcVSITECONSTR);
+ if (graph != NULL)
{
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
- wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, &gs,
- (step_rel % gs.nstms == 0) &&
- (multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
- lastbox,
- top_global, &bSumEkinhOld,
- cglo_flags
- | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
- | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
- | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
- | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
- | CGLO_CONSTRAINT
- );
+ shift_self(graph, state->box, state->x);
}
+ construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ top->idef.iparams, top->idef.il,
+ fr->ePBC, fr->bMolPBC, cr, state->box);
- /* ############# END CALC EKIN AND PRESSURE ################# */
+ if (graph != NULL)
+ {
+ unshift_self(graph, state->box, state->x);
+ }
+ wallcycle_stop(wcycle, ewcVSITECONSTR);
+ }
- /* Note: this is OK, but there are some numerical precision issues with using the convergence of
- the virial that should probably be addressed eventually. state->veta has better properies,
- but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
- generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
+ /* ############## IF NOT VV, Calculate globals HERE ############ */
+ /* With Leap-Frog we can skip compute_globals at
+ * non-communication steps, but we need to calculate
+ * the kinetic energy one step before communication.
+ */
+ if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
+ {
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, &gs,
+ (step_rel % gs.nstms == 0) &&
+ (multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
+ lastbox,
+ top_global, &bSumEkinhOld,
+ cglo_flags
+ | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
+ | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
+ | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
+ | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
+ | CGLO_CONSTRAINT
+ );
}
- /* TODO remove the brace above, once iteration is removed */
+
+ /* ############# END CALC EKIN AND PRESSURE ################# */
+
+ /* Note: this is OK, but there are some numerical precision issues with using the convergence of
+ the virial that should probably be addressed eventually. state->veta has better properies,
+ but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
+ generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
if (ir->efep != efepNO && (!bVV || bRerunMD))
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff