--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "update_vv.h"
+
+#include <cmath>
+#include <cstdio>
+
+#include <algorithm>
+#include <memory>
+
+#include "gromacs/domdec/partition.h"
+#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/coupling.h"
+#include "gromacs/mdlib/enerdata_utils.h"
+#include "gromacs/mdlib/mdatoms.h"
+#include "gromacs/mdlib/md_support.h"
+#include "gromacs/mdlib/stat.h"
+#include "gromacs/mdlib/tgroup.h"
+#include "gromacs/mdlib/update.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/enerdata.h"
+#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcebuffers.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/group.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdtypes/state.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/topology/topology.h"
+
+void integrateVVFirstStep(int64_t step,
+ bool bFirstStep,
+ bool bInitStep,
+ gmx::StartingBehavior startingBehavior,
+ int nstglobalcomm,
+ t_inputrec* ir,
+ t_forcerec* fr,
+ t_commrec* cr,
+ t_state* state,
+ t_mdatoms* mdatoms,
+ const t_fcdata& fcdata,
+ t_extmass* MassQ,
+ t_vcm* vcm,
+ const gmx_mtop_t* top_global,
+ const gmx_localtop_t& top,
+ gmx_enerdata_t* enerd,
+ gmx_ekindata_t* ekind,
+ gmx_global_stat* gstat,
+ real* last_ekin,
+ bool bCalcVir,
+ tensor total_vir,
+ tensor shake_vir,
+ tensor force_vir,
+ tensor pres,
+ matrix M,
+ bool do_log,
+ bool do_ene,
+ bool bCalcEner,
+ bool bGStat,
+ bool bStopCM,
+ bool bTrotter,
+ bool bExchanged,
+ bool* bSumEkinhOld,
+ bool* shouldCheckNumberOfBondedInteractions,
+ real* saved_conserved_quantity,
+ gmx::ForceBuffers* f,
+ gmx::Update* upd,
+ gmx::Constraints* constr,
+ gmx::SimulationSignaller* nullSignaller,
+ std::array<std::vector<int>, ettTSEQMAX> trotter_seq,
+ t_nrnb* nrnb,
+ const gmx::MDLogger& mdlog,
+ FILE* fplog,
+ gmx_wallcycle* wcycle)
+{
+ if (!bFirstStep || startingBehavior == gmx::StartingBehavior::NewSimulation)
+ {
+ /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
+ rvec* vbuf = nullptr;
+
+ wallcycle_start(wcycle, ewcUPDATE);
+ if (ir->eI == eiVV && bInitStep)
+ {
+ /* if using velocity verlet with full time step Ekin,
+ * take the first half step only to compute the
+ * virial for the first step. From there,
+ * revert back to the initial coordinates
+ * so that the input is actually the initial step.
+ */
+ snew(vbuf, state->natoms);
+ copy_rvecn(state->v.rvec_array(), vbuf, 0,
+ state->natoms); /* should make this better for parallelizing? */
+ }
+ else
+ {
+ /* this is for NHC in the Ekin(t+dt/2) version of vv */
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ, trotter_seq, ettTSEQ1);
+ }
+
+ upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind,
+ M, etrtVELOCITY1, cr, constr != nullptr);
+
+ wallcycle_stop(wcycle, ewcUPDATE);
+ constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
+ wallcycle_start(wcycle, ewcUPDATE);
+ /* 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 || inputrecNptTrotter(ir)); */
+ /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
+ bool bPres = TRUE;
+ bool 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);
+ int totalNumberOfBondedInteractions = -1;
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
+ enerd, force_vir, shake_vir, total_vir, pres, constr, nullSignaller,
+ state->box, &totalNumberOfBondedInteractions, bSumEkinhOld,
+ (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
+ | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
+ | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
+ | (*shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+ : 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 */
+ checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global,
+ &top, makeConstArrayRef(state->x), state->box,
+ shouldCheckNumberOfBondedInteractions);
+ if (bStopCM)
+ {
+ process_and_stopcm_grp(fplog, vcm, *mdatoms, makeArrayRef(state->x),
+ makeArrayRef(state->v));
+ inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
+ }
+ wallcycle_start(wcycle, ewcUPDATE);
+ }
+ /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
+ if (!bInitStep)
+ {
+ 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);
+
+ /* TODO This is only needed when we're about to write
+ * a checkpoint, because we use it after the restart
+ * (in a kludge?). But what should we be doing if
+ * the startingBehavior is NewSimulation or bInitStep are true? */
+ if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
+ {
+ copy_mat(shake_vir, state->svir_prev);
+ copy_mat(force_vir, state->fvir_prev);
+ }
+ if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
+ {
+ /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
+ enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
+ enerd->term[F_EKIN] = trace(ekind->ekin);
+ }
+ }
+ else 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(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
+ enerd, nullptr, nullptr, nullptr, nullptr, constr, nullSignaller,
+ state->box, nullptr, bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
+ wallcycle_start(wcycle, ewcUPDATE);
+ }
+ }
+ /* if it's the initial step, we performed this first step just to get the constraint virial */
+ if (ir->eI == eiVV && bInitStep)
+ {
+ copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
+ sfree(vbuf);
+ }
+ wallcycle_stop(wcycle, ewcUPDATE);
+ }
+
+ /* compute the conserved quantity */
+ *saved_conserved_quantity = NPT_energy(ir, state, MassQ);
+ if (ir->eI == eiVV)
+ {
+ *last_ekin = enerd->term[F_EKIN];
+ }
+ if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
+ {
+ *saved_conserved_quantity -= enerd->term[F_DISPCORR];
+ }
+ /* sum up the foreign kinetic energy and dK/dl terms for vv. currently done every step so that dhdl is correct in the .edr */
+ if (ir->efep != efepNO)
+ {
+ accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
+ }
+}
+
+void integrateVVSecondStep(int64_t step,
+ t_inputrec* ir,
+ t_forcerec* fr,
+ t_commrec* cr,
+ t_state* state,
+ t_mdatoms* mdatoms,
+ const t_fcdata& fcdata,
+ t_extmass* MassQ,
+ t_vcm* vcm,
+ pull_t* pull_work,
+ gmx_enerdata_t* enerd,
+ gmx_ekindata_t* ekind,
+ gmx_global_stat* gstat,
+ real* dvdl_constr,
+ bool bCalcVir,
+ tensor total_vir,
+ tensor shake_vir,
+ tensor force_vir,
+ tensor pres,
+ matrix M,
+ matrix lastbox,
+ bool do_log,
+ bool do_ene,
+ bool bGStat,
+ bool* bSumEkinhOld,
+ gmx::ForceBuffers* f,
+ std::vector<gmx::RVec>* cbuf,
+ gmx::Update* upd,
+ gmx::Constraints* constr,
+ gmx::SimulationSignaller* nullSignaller,
+ std::array<std::vector<int>, ettTSEQMAX> trotter_seq,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle)
+{
+ /* velocity half-step update */
+ upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M,
+ etrtVELOCITY2, cr, constr != nullptr);
+
+
+ /* Above, initialize just copies ekinh into ekin,
+ * it doesn't copy position (for VV),
+ * and entire integrator for MD.
+ */
+
+ if (ir->eI == eiVVAK)
+ {
+ cbuf->resize(state->x.size());
+ std::copy(state->x.begin(), state->x.end(), cbuf->begin());
+ }
+
+ if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
+ {
+ updatePrevStepPullCom(pull_work, state);
+ }
+
+ upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M,
+ etrtPOSITION, cr, constr != nullptr);
+
+ wallcycle_stop(wcycle, ewcUPDATE);
+
+ constrain_coordinates(constr, do_log, do_ene, step, state, upd->xp()->arrayRefWithPadding(),
+ dvdl_constr, bCalcVir, shake_vir);
+
+ upd->update_sd_second_half(*ir, step, dvdl_constr, mdatoms, state, cr, nrnb, wcycle, constr,
+ do_log, do_ene);
+ upd->finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
+
+ 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(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle, enerd,
+ force_vir, shake_vir, total_vir, pres, constr, nullSignaller, lastbox,
+ nullptr, bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | 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 */
+ std::copy(cbuf->begin(), cbuf->end(), state->x.begin());
+
+ upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind,
+ M, etrtPOSITION, cr, constr != nullptr);
+ wallcycle_stop(wcycle, ewcUPDATE);
+
+ /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) 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*/
+ upd->finish_update(*ir, mdatoms, state, wcycle, false);
+ }
+ /* 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 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;
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifndef GMX_MDLIB_UPDATE_VV_H
+#define GMX_MDLIB_UPDATE_VV_H
+
+#include <vector>
+
+#include "gromacs/math/vectypes.h"
+#include "gromacs/mdrunutility/handlerestart.h"
+#include "gromacs/mdtypes/md_enums.h"
+
+struct gmx_ekindata_t;
+struct gmx_enerdata_t;
+struct gmx_global_stat;
+struct gmx_localtop_t;
+struct gmx_mtop_t;
+struct gmx_wallcycle;
+struct pull_t;
+struct t_commrec;
+struct t_extmass;
+struct t_fcdata;
+struct t_forcerec;
+struct t_inputrec;
+struct t_mdatoms;
+struct t_nrnb;
+class t_state;
+struct t_vcm;
+
+namespace gmx
+{
+class Constraints;
+class ForceBuffers;
+class MDLogger;
+class SimulationSignaller;
+class Update;
+} // namespace gmx
+
+/*! \brief Make the first step of Velocity Verlet integration
+ *
+ * \param[in] step Current timestep.
+ * \param[in] bFirstStep Is it a first step.
+ * \param[in] bInitStep Is it an initialization step.
+ * \param[in] startingBehavior Describes whether this is a restart appending to output files.
+ * \param[in] nstglobalcomm Will globals be computed on this step.
+ * \param[in] ir Input record.
+ * \param[in] fr Force record.
+ * \param[in] cr Comunication record.
+ * \param[in] state Simulation state.
+ * \param[in] mdatoms MD atoms data.
+ * \param[in] fcdata Force calculation data.
+ * \param[in] MassQ Mass/pressure data.
+ * \param[in] vcm Center of mass motion removal.
+ * \param[in] top_global Global topology.
+ * \param[in] top Local topology.
+ * \param[in] enerd Energy data.
+ * \param[in] ekind Kinetic energy data.
+ * \param[in] gstat Storage of thermodynamic parameters data.
+ * \param[out] last_ekin Kinetic energies of the last step.
+ * \param[in] bCalcVir If the virial is computed on this step.
+ * \param[in] total_vir Total virial tensor.
+ * \param[in] shake_vir Constraints virial.
+ * \param[in] force_vir Force virial.
+ * \param[in] pres Pressure tensor.
+ * \param[in] M Parrinello-Rahman velocity scaling matrix.
+ * \param[in] do_log Do logging on this step.
+ * \param[in] do_ene Print energies on this step.
+ * \param[in] bCalcEner Compute energies on this step.
+ * \param[in] bGStat Collect globals this step.
+ * \param[in] bStopCM Stop the center of mass motion on this step.
+ * \param[in] bTrotter Do trotter routines this step.
+ * \param[in] bExchanged If this is a replica exchange step.
+ * \param[out] bSumEkinhOld Old kinetic energies will need to be summed up.
+ * \param[out] shouldCheckNumberOfBondedInteractions If checks for missing bonded
+ * interactions will be needed.
+ * \param[out] saved_conserved_quantity Place to store the conserved energy.
+ * \param[in] f Force buffers.
+ * \param[in] upd Update object.
+ * \param[in] constr Constraints object.
+ * \param[in] nullSignaller Simulation signaller.
+ * \param[in] trotter_seq NPT variables.
+ * \param[in] nrnb Cycle counters.
+ * \param[in] mdlog Logger.
+ * \param[in] fplog Another logger.
+ * \param[in] wcycle Wall-clock cycle counter.
+ */
+void integrateVVFirstStep(int64_t step,
+ bool bFirstStep,
+ bool bInitStep,
+ gmx::StartingBehavior startingBehavior,
+ int nstglobalcomm,
+ t_inputrec* ir,
+ t_forcerec* fr,
+ t_commrec* cr,
+ t_state* state,
+ t_mdatoms* mdatoms,
+ const t_fcdata& fcdata,
+ t_extmass* MassQ,
+ t_vcm* vcm,
+ const gmx_mtop_t* top_global,
+ const gmx_localtop_t& top,
+ gmx_enerdata_t* enerd,
+ gmx_ekindata_t* ekind,
+ gmx_global_stat* gstat,
+ real* last_ekin,
+ bool bCalcVir,
+ tensor total_vir,
+ tensor shake_vir,
+ tensor force_vir,
+ tensor pres,
+ matrix M,
+ bool do_log,
+ bool do_ene,
+ bool bCalcEner,
+ bool bGStat,
+ bool bStopCM,
+ bool bTrotter,
+ bool bExchanged,
+ bool* bSumEkinhOld,
+ bool* shouldCheckNumberOfBondedInteractions,
+ real* saved_conserved_quantity,
+ gmx::ForceBuffers* f,
+ gmx::Update* upd,
+ gmx::Constraints* constr,
+ gmx::SimulationSignaller* nullSignaller,
+ std::array<std::vector<int>, ettTSEQMAX> trotter_seq,
+ t_nrnb* nrnb,
+ const gmx::MDLogger& mdlog,
+ FILE* fplog,
+ gmx_wallcycle* wcycle);
+
+
+/*! \brief Make the second step of Velocity Verlet integration
+ *
+ * \param[in] step Current timestep.
+ * \param[in] ir Input record.
+ * \param[in] fr Force record.
+ * \param[in] cr Comunication record.
+ * \param[in] state Simulation state.
+ * \param[in] mdatoms MD atoms data.
+ * \param[in] fcdata Force calculation data.
+ * \param[in] MassQ Mass/pressure data.
+ * \param[in] vcm Center of mass motion removal.
+ * \param[in] pull_work Pulling data.
+ * \param[in] enerd Energy data.
+ * \param[in] ekind Kinetic energy data.
+ * \param[in] gstat Storage of thermodynamic parameters data.
+ * \param[out] dvdl_constr FEP data for constraints.
+ * \param[in] bCalcVir If the virial is computed on this step.
+ * \param[in] total_vir Total virial tensor.
+ * \param[in] shake_vir Constraints virial.
+ * \param[in] force_vir Force virial.
+ * \param[in] pres Pressure tensor.
+ * \param[in] M Parrinello-Rahman velocity scaling matrix.
+ * \param[in] lastbox Last recorded PBC box.
+ * \param[in] do_log Do logging on this step.
+ * \param[in] do_ene Print energies on this step.
+ * \param[in] bGStat Collect globals this step.
+ * \param[out] bSumEkinhOld Old kinetic energies need to be summed up.
+ * \param[in] f Force buffers.
+ * \param[in] cbuf Buffer to store intermediate coordinates
+ * \param[in] upd Update object.
+ * \param[in] constr Constraints object.
+ * \param[in] nullSignaller Simulation signaller.
+ * \param[in] trotter_seq NPT variables.
+ * \param[in] nrnb Cycle counters.
+ * \param[in] wcycle Wall-clock cycle counter.
+ */
+void integrateVVSecondStep(int64_t step,
+ t_inputrec* ir,
+ t_forcerec* fr,
+ t_commrec* cr,
+ t_state* state,
+ t_mdatoms* mdatoms,
+ const t_fcdata& fcdata,
+ t_extmass* MassQ,
+ t_vcm* vcm,
+ pull_t* pull_work,
+ gmx_enerdata_t* enerd,
+ gmx_ekindata_t* ekind,
+ gmx_global_stat* gstat,
+ real* dvdl_constr,
+ bool bCalcVir,
+ tensor total_vir,
+ tensor shake_vir,
+ tensor force_vir,
+ tensor pres,
+ matrix M,
+ matrix lastbox,
+ bool do_log,
+ bool do_ene,
+ bool bGStat,
+ bool* bSumEkinhOld,
+ gmx::ForceBuffers* f,
+ std::vector<gmx::RVec>* cbuf,
+ gmx::Update* upd,
+ gmx::Constraints* constr,
+ gmx::SimulationSignaller* nullSignaller,
+ std::array<std::vector<int>, ettTSEQMAX> trotter_seq,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle);
+
+
+#endif // GMX_MDLIB_UPDATE_VV_H
#include "gromacs/mdlib/trajectory_writing.h"
#include "gromacs/mdlib/update.h"
#include "gromacs/mdlib/update_constrain_gpu.h"
+#include "gromacs/mdlib/update_vv.h"
#include "gromacs/mdlib/vcm.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdrunutility/handlerestart.h"
gmx_global_stat_t gstat;
gmx_shellfc_t* shellfc;
gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
- gmx_bool bTemp, bPres, bTrotter;
+ gmx_bool bTrotter;
real dvdl_constr;
std::vector<RVec> cbuf;
matrix lastbox;
if (bExchanged)
{
-
/* 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 . . . . */
// if it is the first step after starting from a checkpoint.
// That is, the half step is needed on all other steps, and
// also the first step when starting from a .tpr file.
- if (EI_VV(ir->eI) && (!bFirstStep || startingBehavior == StartingBehavior::NewSimulation))
- /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
- {
- rvec* vbuf = nullptr;
-
- wallcycle_start(wcycle, ewcUPDATE);
- if (ir->eI == eiVV && bInitStep)
- {
- /* if using velocity verlet with full time step Ekin,
- * take the first half step only to compute the
- * virial for the first step. From there,
- * revert back to the initial coordinates
- * so that the input is actually the initial step.
- */
- snew(vbuf, state->natoms);
- copy_rvecn(state->v.rvec_array(), vbuf, 0,
- state->natoms); /* should make this better for parallelizing? */
- }
- else
- {
- /* this is for NHC in the Ekin(t+dt/2) version of vv */
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ,
- trotter_seq, ettTSEQ1);
- }
-
- upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
- M, etrtVELOCITY1, cr, constr != nullptr);
-
- wallcycle_stop(wcycle, ewcUPDATE);
- constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
- wallcycle_start(wcycle, ewcUPDATE);
- /* 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 || inputrecNptTrotter(ir)); */
- /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(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(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
- makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
- enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
- state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
- (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
- | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
- | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
- | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
- : 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 */
- checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
- top_global, &top, makeConstArrayRef(state->x),
- state->box, &shouldCheckNumberOfBondedInteractions);
- if (bStopCM)
- {
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, makeArrayRef(state->x),
- makeArrayRef(state->v));
- inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
- }
- wallcycle_start(wcycle, ewcUPDATE);
- }
- /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
- if (!bInitStep)
- {
- 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);
-
- /* TODO This is only needed when we're about to write
- * a checkpoint, because we use it after the restart
- * (in a kludge?). But what should we be doing if
- * the startingBehavior is NewSimulation or bInitStep are true? */
- if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
- {
- copy_mat(shake_vir, state->svir_prev);
- copy_mat(force_vir, state->fvir_prev);
- }
- if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
- {
- /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
- enerd->term[F_TEMP] =
- sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
- enerd->term[F_EKIN] = trace(ekind->ekin);
- }
- }
- else 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(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
- makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
- enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
- state->box, nullptr, &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
- wallcycle_start(wcycle, ewcUPDATE);
- }
- }
- /* if it's the initial step, we performed this first step just to get the constraint virial */
- if (ir->eI == eiVV && bInitStep)
- {
- copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
- sfree(vbuf);
- }
- wallcycle_stop(wcycle, ewcUPDATE);
- }
-
- /* compute the conserved quantity */
if (EI_VV(ir->eI))
{
- saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
- if (ir->eI == eiVV)
- {
- last_ekin = enerd->term[F_EKIN];
- }
- if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
- {
- saved_conserved_quantity -= enerd->term[F_DISPCORR];
- }
- /* sum up the foreign kinetic energy and dK/dl terms for vv. currently done every step so that dhdl is correct in the .edr */
- if (ir->efep != efepNO)
- {
- accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
- }
+ integrateVVFirstStep(step, bFirstStep, bInitStep, startingBehavior, nstglobalcomm, ir,
+ fr, cr, state, mdatoms, fcdata, &MassQ, &vcm, top_global, top, enerd,
+ ekind, gstat, &last_ekin, bCalcVir, total_vir, shake_vir, force_vir,
+ pres, M, do_log, do_ene, bCalcEner, bGStat, bStopCM, bTrotter,
+ bExchanged, &bSumEkinhOld, &shouldCheckNumberOfBondedInteractions,
+ &saved_conserved_quantity, &f, &upd, constr, &nullSignaller,
+ trotter_seq, nrnb, mdlog, fplog, wcycle);
}
/* ######## END FIRST UPDATE STEP ############## */
update_pcouple_before_coordinates(fplog, step, ir, state, pressureCouplingMu, M, bInitStep);
}
- if (EI_VV(ir->eI))
- {
- /* velocity half-step update */
- upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
- M, etrtVELOCITY2, cr, constr != nullptr);
- }
-
- /* Above, initialize just copies ekinh into ekin,
- * it doesn't copy position (for VV),
- * and entire integrator for MD.
- */
-
- if (ir->eI == eiVVAK)
- {
- cbuf.resize(state->x.size());
- std::copy(state->x.begin(), state->x.end(), cbuf.begin());
- }
-
/* With leap-frog type integrators we compute the kinetic energy
* at a whole time step as the average of the half-time step kinetic
* energies of two subsequent steps. Therefore we need to compute the
const bool doParrinelloRahman = (ir->epc == epcPARRINELLORAHMAN
&& do_per_step(step + ir->nstpcouple - 1, ir->nstpcouple));
- if (useGpuForUpdate)
+ if (EI_VV(ir->eI))
{
- wallcycle_stop(wcycle, ewcUPDATE);
+ GMX_ASSERT(!useGpuForUpdate, "GPU update is not supported with VVAK integrator.");
- if (bNS && (bFirstStep || DOMAINDECOMP(cr)))
+ integrateVVSecondStep(step, ir, fr, cr, state, mdatoms, fcdata, &MassQ, &vcm, pull_work,
+ enerd, ekind, gstat, &dvdl_constr, bCalcVir, total_vir, shake_vir,
+ force_vir, pres, M, lastbox, do_log, do_ene, bGStat, &bSumEkinhOld,
+ &f, &cbuf, &upd, constr, &nullSignaller, trotter_seq, nrnb, wcycle);
+ }
+ else
+ {
+ if (useGpuForUpdate)
{
- integrator->set(stateGpu->getCoordinates(), stateGpu->getVelocities(),
- stateGpu->getForces(), top.idef, *mdatoms, ekind->ngtc);
- // Copy data to the GPU after buffers might have being reinitialized
- stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
- stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
- }
+ wallcycle_stop(wcycle, ewcUPDATE);
- if (simulationWork.useGpuPme && !runScheduleWork->simulationWork.useGpuPmePpCommunication
- && !thisRankHasDuty(cr, DUTY_PME))
- {
- // The PME forces were recieved to the host, so have to be copied
- stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::All);
- }
- else if (!runScheduleWork->stepWork.useGpuFBufferOps)
- {
- // The buffer ops were not offloaded this step, so the forces are on the
- // host and have to be copied
- stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local);
- }
+ if (bNS && (bFirstStep || DOMAINDECOMP(cr)))
+ {
+ integrator->set(stateGpu->getCoordinates(), stateGpu->getVelocities(),
+ stateGpu->getForces(), top.idef, *mdatoms, ekind->ngtc);
- const bool doTemperatureScaling =
- (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
+ // Copy data to the GPU after buffers might have being reinitialized
+ stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
+ stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
+ }
- // This applies Leap-Frog, LINCS and SETTLE in succession
- integrator->integrate(stateGpu->getForcesReadyOnDeviceEvent(
- AtomLocality::Local, runScheduleWork->stepWork.useGpuFBufferOps),
- ir->delta_t, true, bCalcVir, shake_vir, doTemperatureScaling,
- ekind->tcstat, doParrinelloRahman, ir->nstpcouple * ir->delta_t, M);
+ if (simulationWork.useGpuPme && !runScheduleWork->simulationWork.useGpuPmePpCommunication
+ && !thisRankHasDuty(cr, DUTY_PME))
+ {
+ // The PME forces were recieved to the host, so have to be copied
+ stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::All);
+ }
+ else if (!runScheduleWork->stepWork.useGpuFBufferOps)
+ {
+ // The buffer ops were not offloaded this step, so the forces are on the
+ // host and have to be copied
+ stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local);
+ }
- // Copy velocities D2H after update if:
- // - Globals are computed this step (includes the energy output steps).
- // - Temperature is needed for the next step.
- if (bGStat || needHalfStepKineticEnergy)
- {
- stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
- stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
+ const bool doTemperatureScaling =
+ (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
+
+ // This applies Leap-Frog, LINCS and SETTLE in succession
+ integrator->integrate(
+ stateGpu->getForcesReadyOnDeviceEvent(
+ AtomLocality::Local, runScheduleWork->stepWork.useGpuFBufferOps),
+ ir->delta_t, true, bCalcVir, shake_vir, doTemperatureScaling, ekind->tcstat,
+ doParrinelloRahman, ir->nstpcouple * ir->delta_t, M);
+
+ // Copy velocities D2H after update if:
+ // - Globals are computed this step (includes the energy output steps).
+ // - Temperature is needed for the next step.
+ if (bGStat || needHalfStepKineticEnergy)
+ {
+ stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
+ stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
+ }
}
- }
- else
- {
- /* With multiple time stepping we need to do an additional normal
- * update step to obtain the virial, as the actual MTS integration
- * using an acceleration where the slow forces are multiplied by mtsFactor.
- * Using that acceleration would result in a virial with the slow
- * force contribution would be a factor mtsFactor too large.
- */
- if (fr->useMts && bCalcVir && constr != nullptr)
+ else
{
- upd.update_for_constraint_virial(*ir, *mdatoms, *state, f.view().forceWithPadding(), *ekind);
+ /* With multiple time stepping we need to do an additional normal
+ * update step to obtain the virial, as the actual MTS integration
+ * using an acceleration where the slow forces are multiplied by mtsFactor.
+ * Using that acceleration would result in a virial with the slow
+ * force contribution would be a factor mtsFactor too large.
+ */
+ if (fr->useMts && bCalcVir && constr != nullptr)
+ {
+ upd.update_for_constraint_virial(*ir, *mdatoms, *state,
+ f.view().forceWithPadding(), *ekind);
- constrain_coordinates(constr, do_log, do_ene, step, state,
- upd.xp()->arrayRefWithPadding(), &dvdl_constr, bCalcVir, shake_vir);
- }
+ constrain_coordinates(constr, do_log, do_ene, step, state,
+ upd.xp()->arrayRefWithPadding(), &dvdl_constr, bCalcVir,
+ shake_vir);
+ }
- ArrayRefWithPadding<const RVec> forceCombined =
- (fr->useMts && step % ir->mtsLevels[1].stepFactor == 0)
- ? f.view().forceMtsCombinedWithPadding()
- : f.view().forceWithPadding();
- upd.update_coords(*ir, step, mdatoms, state, forceCombined, fcdata, ekind, M,
- etrtPOSITION, cr, constr != nullptr);
+ ArrayRefWithPadding<const RVec> forceCombined =
+ (fr->useMts && step % ir->mtsLevels[1].stepFactor == 0)
+ ? f.view().forceMtsCombinedWithPadding()
+ : f.view().forceWithPadding();
+ upd.update_coords(*ir, step, mdatoms, state, forceCombined, fcdata, ekind, M,
+ etrtPOSITION, cr, constr != nullptr);
- wallcycle_stop(wcycle, ewcUPDATE);
+ wallcycle_stop(wcycle, ewcUPDATE);
- constrain_coordinates(constr, do_log, do_ene, step, state, upd.xp()->arrayRefWithPadding(),
- &dvdl_constr, bCalcVir && !fr->useMts, shake_vir);
+ constrain_coordinates(constr, do_log, do_ene, step, state,
+ upd.xp()->arrayRefWithPadding(), &dvdl_constr,
+ bCalcVir && !fr->useMts, shake_vir);
- upd.update_sd_second_half(*ir, step, &dvdl_constr, mdatoms, state, cr, nrnb, wcycle,
- constr, do_log, do_ene);
- upd.finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
- }
+ upd.update_sd_second_half(*ir, step, &dvdl_constr, mdatoms, state, cr, nrnb, wcycle,
+ constr, do_log, do_ene);
+ upd.finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
+ }
- if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
- {
- updatePrevStepPullCom(pull_work, state);
- }
+ if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
+ {
+ updatePrevStepPullCom(pull_work, state);
+ }
- 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(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
- makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle, enerd,
- force_vir, shake_vir, total_vir, pres, constr, &nullSignaller, lastbox,
- nullptr, &bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | 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 */
- std::copy(cbuf.begin(), cbuf.end(), state->x.begin());
-
- upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
- M, etrtPOSITION, cr, constr != nullptr);
- wallcycle_stop(wcycle, ewcUPDATE);
-
- /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) 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*/
- upd.finish_update(*ir, mdatoms, state, wcycle, false);
- }
- if (EI_VV(ir->eI))
- {
- /* 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 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;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff