From 0787bf7e66a9a1c432578ca0bd93d27432d795a7 Mon Sep 17 00:00:00 2001 From: Artem Zhmurov Date: Fri, 6 Nov 2020 10:14:13 +0300 Subject: [PATCH] Separate VV first and second steps from MD loop This combines code blocks that are only active for VV integrators, and moves them into a separate source file. Makes logic in the do_md(..) more clear and prepares for the removal of VV from the legacy integrator. --- src/gromacs/mdlib/update_vv.cpp | 367 ++++++++++++++++++++++++++++++++ src/gromacs/mdlib/update_vv.h | 238 +++++++++++++++++++++ src/gromacs/mdrun/md.cpp | 364 ++++++++----------------------- 3 files changed, 696 insertions(+), 273 deletions(-) create mode 100644 src/gromacs/mdlib/update_vv.cpp create mode 100644 src/gromacs/mdlib/update_vv.h diff --git a/src/gromacs/mdlib/update_vv.cpp b/src/gromacs/mdlib/update_vv.cpp new file mode 100644 index 0000000000..afc6b4d069 --- /dev/null +++ b/src/gromacs/mdlib/update_vv.cpp @@ -0,0 +1,367 @@ +/* + * 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 +#include + +#include +#include + +#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, 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* cbuf, + gmx::Update* upd, + gmx::Constraints* constr, + gmx::SimulationSignaller* nullSignaller, + std::array, 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; +} diff --git a/src/gromacs/mdlib/update_vv.h b/src/gromacs/mdlib/update_vv.h new file mode 100644 index 0000000000..8699f6b9ef --- /dev/null +++ b/src/gromacs/mdlib/update_vv.h @@ -0,0 +1,238 @@ +/* + * 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 + +#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, 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* cbuf, + gmx::Update* upd, + gmx::Constraints* constr, + gmx::SimulationSignaller* nullSignaller, + std::array, ettTSEQMAX> trotter_seq, + t_nrnb* nrnb, + gmx_wallcycle* wcycle); + + +#endif // GMX_MDLIB_UPDATE_VV_H diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index f5c0fd393e..4c274e8611 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -101,6 +101,7 @@ #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" @@ -178,7 +179,7 @@ void gmx::LegacySimulator::do_md() 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 cbuf; matrix lastbox; @@ -857,7 +858,6 @@ void gmx::LegacySimulator::do_md() 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 . . . . */ @@ -952,151 +952,15 @@ void gmx::LegacySimulator::do_md() // 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 ############## */ @@ -1222,24 +1086,6 @@ void gmx::LegacySimulator::do_md() 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 @@ -1253,133 +1099,105 @@ void gmx::LegacySimulator::do_md() 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 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 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; } -- 2.22.0