debug_gmx();
- /* set free energy calculation frequency as the minimum
- greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/
- nstfep = ir->fepvals->nstdhdl;
- if (ir->bExpanded)
+ if (IR_TWINRANGE(*ir) && repl_ex_nst % ir->nstcalclr != 0)
{
- nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
+ /* We should exchange at nstcalclr steps to get correct integration */
+ gmx_fatal(FARGS, "The replica exchange period (%d) is not divisible by nstcalclr (%d)", repl_ex_nst, ir->nstcalclr);
}
- if (repl_ex_nst > 0)
+
+ if (ir->efep != efepNO)
{
- nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
+ /* Set free energy calculation frequency as the greatest common
+ * denominator of nstdhdl and repl_ex_nst.
+ * Check for nstcalclr with twin-range, since we need the long-range
+ * contribution to the free-energy at the correct (nstcalclr) steps.
+ */
+ nstfep = ir->fepvals->nstdhdl;
+ if (ir->bExpanded)
+ {
+ if (IR_TWINRANGE(*ir) &&
+ ir->expandedvals->nstexpanded % ir->nstcalclr != 0)
+ {
+ gmx_fatal(FARGS, "nstexpanded should be divisible by nstcalclr");
+ }
++ nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
+ }
+ if (repl_ex_nst > 0)
+ {
+ nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
+ }
+ /* We checked divisibility of repl_ex_nst and nstcalclr above */
+ if (IR_TWINRANGE(*ir) && nstfep % ir->nstcalclr != 0)
+ {
+ gmx_incons("nstfep not divisible by nstcalclr");
+ }
}
/* Be REALLY careful about what flags you set here. You CANNOT assume
ekind, M, upd, bInitStep, etrtVELOCITY1,
cr, nrnb, constr, &top->idef);
- if (bIterativeCase && do_per_step(step-1, ir->nstpcouple) && !bInitStep)
+ if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
{
- gmx_iterate_init(&iterate, TRUE);
+ 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);
+ }
}
- /* for iterations, we save these vectors, as we will be self-consistently iterating
- the calculations */
-
- /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
-
- /* save the state */
- if (iterate.bIterationActive)
+ else if (graph)
{
- copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
+ /* Need to unshift here if a do_force has been
+ called in the previous step */
+ unshift_self(graph, state->box, state->x);
}
-
- bFirstIterate = TRUE;
- while (bFirstIterate || iterate.bIterationActive)
+ /* 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)
{
- if (iterate.bIterationActive)
- {
- copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
- if (bFirstIterate && bTrotter)
- {
- /* The first time through, we need a decent first estimate
- of veta(t+dt) to compute the constraints. Do
- this by computing the box volume part of the
- trotter integration at this time. Nothing else
- should be changed by this routine here. If
- !(first time), we start with the previous value
- of veta. */
-
- veta_save = state->veta;
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ0);
- vetanew = state->veta;
- state->veta = veta_save;
- }
- }
-
- bOK = TRUE;
- if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
- {
- wallcycle_stop(wcycle, ewcUPDATE);
- update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
- state, fr->bMolPBC, graph, f,
- &top->idef, shake_vir,
- cr, nrnb, wcycle, upd, constr,
- TRUE, bCalcVir, vetanew);
- 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);
- }
-
- if (!bOK)
- {
- gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
- }
-
- }
- 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) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
- {
- 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))
+ 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)
++ | (bStopCM ? CGLO_STOPCM : 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)
- | (bStopCM ? CGLO_STOPCM : 0)
- | ((iterate.bIterationActive) ? CGLO_ITERATE : 0)
- | (bFirstIterate ? CGLO_FIRSTITERATE : 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);
}
}
-
- if (iterate.bIterationActive &&
- done_iterating(cr, fplog, step, &iterate, bFirstIterate,
- state->veta, &vetanew))
- {
- break;
- }
- bFirstIterate = FALSE;
}
-
if (bTrotter && !bInitStep)
{
copy_mat(shake_vir, state->svir_prev);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff