gs.sig[eglsCHKPT] = 1;
}
- /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */
- if (EI_VV(ir->eI))
+ /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
+ in preprocessing */
+
+ if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
{
- if (!bInitStep)
+ gmx_bool bIfRandomize;
+ bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr);
+ /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
+ if (constr && bIfRandomize)
{
- update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
- }
- if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
- {
- gmx_bool bIfRandomize;
- bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr);
- /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
- if (constr && bIfRandomize)
- {
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
- &top->idef, tmp_vir,
- cr, nrnb, wcycle, upd, constr,
- TRUE, bCalcVir);
- }
- update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
++ update_constraints(fplog, step, NULL, ir, mdatoms,
+ state, fr->bMolPBC, graph, f,
+ &top->idef, tmp_vir,
+ cr, nrnb, wcycle, upd, constr,
- TRUE, bCalcVir, vetanew);
++ TRUE, bCalcVir);
}
}
+ /* ######### 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);
- if (bIterativeCase && do_per_step(step, ir->nstpcouple))
- {
- gmx_iterate_init(&iterate, TRUE);
- /* for iterations, we save these vectors, as we will be redoing the calculations */
- copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
- }
+ dvdl_constr = 0;
- bFirstIterate = TRUE;
- while (bFirstIterate || iterate.bIterationActive)
+ if (!bRerunMD || rerun_fr.bV || bForceUpdate)
{
- /* We now restore these vectors to redo the calculation with improved extended variables */
- if (iterate.bIterationActive)
+ wallcycle_start(wcycle, ewcUPDATE);
+ /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
+ if (bTrotter)
{
- copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
+ 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);
}
- /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
- so scroll down for that logic */
+ if (bVV)
+ {
+ bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
- /* ######### 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);
+ /* 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);
+ }
- bOK = TRUE;
- dvdl_constr = 0;
+ /* 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)
{
- if (iterate.bIterationActive)
- {
- if (bFirstIterate)
- {
- scalevir = 1;
- }
- else
- {
- /* we use a new value of scalevir to converge the iterations faster */
- scalevir = tracevir/trace(shake_vir);
- }
- msmul(shake_vir, scalevir, shake_vir);
- m_add(force_vir, shake_vir, total_vir);
- clear_mat(shake_vir);
- }
- 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,
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff