From 25e7289ac8fd66c6e5089df58e72885917158c07 Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Wed, 5 May 2021 06:41:53 +0000 Subject: [PATCH] Remove mdatoms from coupling code and use more ArrayRef --- src/gromacs/mdlib/coupling.cpp | 257 ++++++++++++++++---------------- src/gromacs/mdlib/coupling.h | 149 ++++++++++-------- src/gromacs/mdlib/tgroup.cpp | 2 +- src/gromacs/mdlib/update_vv.cpp | 42 +++++- src/gromacs/mdrun/md.cpp | 50 ++++++- 5 files changed, 295 insertions(+), 205 deletions(-) diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp index 4a5bd1e606..e6c3613f49 100644 --- a/src/gromacs/mdlib/coupling.cpp +++ b/src/gromacs/mdlib/coupling.cpp @@ -61,7 +61,6 @@ #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" -#include "gromacs/mdtypes/mdatom.h" #include "gromacs/mdtypes/state.h" #include "gromacs/pbcutil/boxutilities.h" #include "gromacs/pbcutil/pbc.h" @@ -97,12 +96,13 @@ static constexpr std::array sy_const = { nullptr, sy_const_ sy_const_3, nullptr, sy_const_5 }; -void update_tcouple(int64_t step, - const t_inputrec* inputrec, - t_state* state, - gmx_ekindata_t* ekind, - const t_extmass* MassQ, - const t_mdatoms* md) +void update_tcouple(int64_t step, + const t_inputrec* inputrec, + t_state* state, + gmx_ekindata_t* ekind, + const t_extmass* MassQ, + int homenr, + gmx::ArrayRef cTC) { // This condition was explicitly checked in previous version, but should have never been satisfied @@ -144,22 +144,18 @@ void update_tcouple(int64_t step, berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral); break; case TemperatureCoupling::NoseHoover: - nosehoover_tcoupl(&(inputrec->opts), - ekind, - dttc, - state->nosehoover_xi.data(), - state->nosehoover_vxi.data(), - MassQ); + nosehoover_tcoupl( + &(inputrec->opts), ekind, dttc, state->nosehoover_xi, state->nosehoover_vxi, MassQ); break; case TemperatureCoupling::VRescale: - vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data()); + vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral); break; default: gmx_fatal(FARGS, "Unknown temperature coupling algorithm"); } /* rescale in place here */ if (EI_VV(inputrec->eI)) { - rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array()); + rescale_velocities(ekind, cTC, 0, homenr, state->v); } } else @@ -179,7 +175,7 @@ void update_pcouple_before_coordinates(FILE* fplog, t_state* state, matrix parrinellorahmanMu, matrix M, - gmx_bool bInitStep) + bool bInitStep) { /* Berendsen P-coupling is completely handled after the coordinate update. * Trotter P-coupling is handled by separate calls to trotter_update(). @@ -194,21 +190,21 @@ void update_pcouple_before_coordinates(FILE* fplog, } } -void update_pcouple_after_coordinates(FILE* fplog, - int64_t step, - const t_inputrec* inputrec, - const t_mdatoms* md, - const matrix pressure, - const matrix forceVirial, - const matrix constraintVirial, - matrix pressureCouplingMu, - t_state* state, - t_nrnb* nrnb, - gmx::BoxDeformation* boxDeformation, - const bool scaleCoordinates) +void update_pcouple_after_coordinates(FILE* fplog, + int64_t step, + const t_inputrec* inputrec, + const int homenr, + gmx::ArrayRef cFREEZE, + const matrix pressure, + const matrix forceVirial, + const matrix constraintVirial, + matrix pressureCouplingMu, + t_state* state, + t_nrnb* nrnb, + gmx::BoxDeformation* boxDeformation, + const bool scaleCoordinates) { - int start = 0; - int homenr = md->homenr; + int start = 0; /* Cast to real for faster code, no loss in precision (see comment above) */ real dt = inputrec->delta_t; @@ -232,17 +228,18 @@ void update_pcouple_after_coordinates(FILE* fplog, constraintVirial, pressureCouplingMu, &state->baros_integral); - pressureCouplingScaleBoxAndCoordinates(inputrec, - pressureCouplingMu, - state->box, - state->box_rel, - start, - homenr, - state->x.rvec_array(), - nullptr, - md->cFREEZE, - nrnb, - scaleCoordinates); + pressureCouplingScaleBoxAndCoordinates( + inputrec, + pressureCouplingMu, + state->box, + state->box_rel, + start, + homenr, + state->x, + gmx::ArrayRef(), + cFREEZE, + nrnb, + scaleCoordinates); } break; case (PressureCoupling::CRescale): @@ -259,18 +256,17 @@ void update_pcouple_after_coordinates(FILE* fplog, constraintVirial, pressureCouplingMu, &state->baros_integral); - pressureCouplingScaleBoxAndCoordinates( - inputrec, - pressureCouplingMu, - state->box, - state->box_rel, - start, - homenr, - state->x.rvec_array(), - state->v.rvec_array(), - md->cFREEZE, - nrnb, - scaleCoordinates); + pressureCouplingScaleBoxAndCoordinates(inputrec, + pressureCouplingMu, + state->box, + state->box_rel, + start, + homenr, + state->x, + state->v, + cFREEZE, + nrnb, + scaleCoordinates); } break; case (PressureCoupling::ParrinelloRahman): @@ -335,13 +331,15 @@ void update_pcouple_after_coordinates(FILE* fplog, } } -extern gmx_bool update_randomize_velocities(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - const t_mdatoms* md, - gmx::ArrayRef v, - const gmx::Update* upd, - const gmx::Constraints* constr) +extern bool update_randomize_velocities(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + int homenr, + gmx::ArrayRef cTC, + gmx::ArrayRef invMass, + gmx::ArrayRef v, + const gmx::Update* upd, + const gmx::Constraints* constr) { real rate = (ir->delta_t) / ir->opts.tau_t[0]; @@ -364,22 +362,12 @@ extern gmx_bool update_randomize_velocities(const t_inputrec* ir, if ((ir->etc == TemperatureCoupling::Andersen) || do_per_step(step, gmx::roundToInt(1.0 / rate))) { andersen_tcoupl( - ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(), upd->getBoltzmanFactor()); + ir, step, cr, homenr, cTC, invMass, v, rate, upd->getAndersenRandomizeGroup(), upd->getBoltzmanFactor()); return TRUE; } return FALSE; } -/* - static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = { - {}, - {1}, - {}, - {0.828981543588751,-0.657963087177502,0.828981543588751}, - {}, - {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065} - };*/ - /* these integration routines are only referenced inside this file */ static void NHC_trotter(const t_grpopts* opts, int nvar, @@ -390,20 +378,20 @@ static void NHC_trotter(const t_grpopts* opts, double scalefac[], real* veta, const t_extmass* MassQ, - gmx_bool bEkinAveVel) + bool bEkinAveVel) { /* general routine for both barostat and thermostat nose hoover chains */ - int i, j, mi, mj; - double Ekin, Efac, reft, kT, nd; - double dt; - double * ivxi, *ixi; - double* GQ; - gmx_bool bBarostat; - int mstepsi, mstepsj; - int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */ - int nh = opts->nhchainlength; + int i, j, mi, mj; + double Ekin, Efac, reft, kT, nd; + double dt; + double *ivxi, *ixi; + double* GQ; + bool bBarostat; + int mstepsi, mstepsj; + int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */ + int nh = opts->nhchainlength; snew(GQ, nh); mstepsi = mstepsj = ns; @@ -663,7 +651,7 @@ void parrinellorahman_pcoupl(FILE* fplog, tensor boxv, tensor M, matrix mu, - gmx_bool bFirstStep) + bool bFirstStep) { /* This doesn't do any coordinate updating. It just * integrates the box vector equations from the calculated @@ -1106,17 +1094,17 @@ void pressureCouplingCalculateScalingMatrix(FILE* fplog, } template -void pressureCouplingScaleBoxAndCoordinates(const t_inputrec* ir, - const matrix mu, - matrix box, - matrix box_rel, - int start, - int nr_atoms, - rvec x[], - rvec v[], - const unsigned short cFREEZE[], - t_nrnb* nrnb, - const bool scaleCoordinates) +void pressureCouplingScaleBoxAndCoordinates(const t_inputrec* ir, + const matrix mu, + matrix box, + matrix box_rel, + int start, + int nr_atoms, + gmx::ArrayRef x, + gmx::ArrayRef v, + gmx::ArrayRef cFREEZE, + t_nrnb* nrnb, + const bool scaleCoordinates) { static_assert(pressureCouplingType == PressureCoupling::Berendsen || pressureCouplingType == PressureCoupling::CRescale, @@ -1139,7 +1127,7 @@ void pressureCouplingScaleBoxAndCoordinates(const t_inputrec* ir, { // Trivial OpenMP region that does not throw int g = 0; - if (cFREEZE != nullptr) + if (!cFREEZE.empty()) { g = cFREEZE[n]; } @@ -1227,14 +1215,16 @@ void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std: } } -void andersen_tcoupl(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - const t_mdatoms* md, - gmx::ArrayRef v, - real rate, - const std::vector& randomize, - gmx::ArrayRef boltzfac) +void andersen_tcoupl(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + const int homenr, + gmx::ArrayRef cTC, + gmx::ArrayRef invMass, + gmx::ArrayRef v, + real rate, + const std::vector& randomize, + gmx::ArrayRef boltzfac) { const int* gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr); int i; @@ -1245,16 +1235,16 @@ void andersen_tcoupl(const t_inputrec* ir, /* randomize the velocities of the selected particles */ - for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */ + for (i = 0; i < homenr; i++) /* now loop over the list of atoms */ { - int ng = gatindex ? gatindex[i] : i; - gmx_bool bRandomize; + int ng = gatindex ? gatindex[i] : i; + bool bRandomize; rng.restart(step, ng); - if (md->cTC) + if (!cTC.empty()) { - gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */ + gc = cTC[i]; /* assign the atom to a temperature group if there are more than one */ } if (randomize[gc]) { @@ -1274,7 +1264,7 @@ void andersen_tcoupl(const t_inputrec* ir, real scal; int d; - scal = std::sqrt(boltzfac[gc] * md->invmass[i]); + scal = std::sqrt(boltzfac[gc] * invMass[i]); normalDist.reset(); @@ -1291,8 +1281,8 @@ void andersen_tcoupl(const t_inputrec* ir, void nosehoover_tcoupl(const t_grpopts* opts, const gmx_ekindata_t* ekind, real dt, - double xi[], - double vxi[], + gmx::ArrayRef xi, + gmx::ArrayRef vxi, const t_extmass* MassQ) { int i; @@ -1309,16 +1299,18 @@ void nosehoover_tcoupl(const t_grpopts* opts, } } -void trotter_update(const t_inputrec* ir, - int64_t step, - gmx_ekindata_t* ekind, - const gmx_enerdata_t* enerd, - t_state* state, - const tensor vir, - const t_mdatoms* md, - const t_extmass* MassQ, - gmx::ArrayRef> trotter_seqlist, - int trotter_seqno) +void trotter_update(const t_inputrec* ir, + int64_t step, + gmx_ekindata_t* ekind, + const gmx_enerdata_t* enerd, + t_state* state, + const tensor vir, + int homenr, + gmx::ArrayRef cTC, + gmx::ArrayRef invMass, + const t_extmass* MassQ, + gmx::ArrayRef> trotter_seqlist, + int trotter_seqno) { int n, i, d, ngtc, gc = 0, t; @@ -1328,7 +1320,7 @@ void trotter_update(const t_inputrec* ir, real dt; double * scalefac, dtc; rvec sumv = { 0, 0, 0 }; - gmx_bool bCouple; + bool bCouple; if (trotter_seqno <= ettTSEQ2) { @@ -1420,11 +1412,11 @@ void trotter_update(const t_inputrec* ir, /* but do we actually need the total? */ /* modify the velocities as well */ - for (n = 0; n < md->homenr; n++) + for (n = 0; n < homenr; n++) { - if (md->cTC) /* does this conditional need to be here? is this always true?*/ + if (!cTC.empty()) /* does this conditional need to be here? is this always true?*/ { - gc = md->cTC[n]; + gc = cTC[n]; } for (d = 0; d < DIM; d++) { @@ -1435,7 +1427,7 @@ void trotter_update(const t_inputrec* ir, { for (d = 0; d < DIM; d++) { - sumv[d] += (v[n][d]) / md->invmass[n]; + sumv[d] += (v[n][d]) / invMass[n]; } } } @@ -1448,7 +1440,7 @@ void trotter_update(const t_inputrec* ir, } -extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bInit) +extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, bool bInit) { int n, i, j, d, ngtc, nh; const t_grpopts* opts; @@ -1546,7 +1538,7 @@ extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* Mas } std::array, ettTSEQMAX> -init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bTrotter) +init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, bool bTrotter) { int i, j, nnhpres, nh; const t_grpopts* opts; @@ -1992,7 +1984,7 @@ real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t ste return ekin_new; } -void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[]) +void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, gmx::ArrayRef therm_integral) { const t_grpopts* opts; int i; @@ -2048,15 +2040,18 @@ void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, } } -void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[]) +void rescale_velocities(const gmx_ekindata_t* ekind, + gmx::ArrayRef cTC, + int start, + int end, + gmx::ArrayRef v) { - const unsigned short* cTC = mdatoms->cTC; gmx::ArrayRef tcstat = ekind->tcstat; for (int n = start; n < end; n++) { int gt = 0; - if (cTC) + if (!cTC.empty()) { gt = cTC[n]; } diff --git a/src/gromacs/mdlib/coupling.h b/src/gromacs/mdlib/coupling.h index affd45bd44..92333c5cca 100644 --- a/src/gromacs/mdlib/coupling.h +++ b/src/gromacs/mdlib/coupling.h @@ -38,11 +38,13 @@ #ifndef GMX_MDLIB_COUPLING_H #define GMX_MDLIB_COUPLING_H -#include "gromacs/math/paddedvector.h" +#include + +#include +#include + #include "gromacs/math/vectypes.h" #include "gromacs/mdtypes/md_enums.h" -#include "gromacs/utility/arrayref.h" -#include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/real.h" class gmx_ekindata_t; @@ -51,7 +53,6 @@ struct t_commrec; struct t_extmass; struct t_grpopts; struct t_inputrec; -struct t_mdatoms; struct t_nrnb; class t_state; @@ -62,17 +63,20 @@ namespace gmx class BoxDeformation; class Constraints; class Update; +template +class ArrayRef; }; // namespace gmx /* Update the size of per-atom arrays (e.g. after DD re-partitioning, which might increase the number of home atoms). */ -void update_tcouple(int64_t step, - const t_inputrec* inputrec, - t_state* state, - gmx_ekindata_t* ekind, - const t_extmass* MassQ, - const t_mdatoms* md); +void update_tcouple(int64_t step, + const t_inputrec* inputrec, + t_state* state, + gmx_ekindata_t* ekind, + const t_extmass* MassQ, + int homenr, + gmx::ArrayRef cTC); /* Update Parrinello-Rahman, to be called before the coordinate update */ void update_pcouple_before_coordinates(FILE* fplog, @@ -81,78 +85,93 @@ void update_pcouple_before_coordinates(FILE* fplog, t_state* state, matrix parrinellorahmanMu, matrix M, - gmx_bool bInitStep); + bool bInitStep); /* Update the box, to be called after the coordinate update. * For Berendsen P-coupling, also calculates the scaling factor * and scales the coordinates. * When the deform option is used, scales coordinates and box here. */ -void update_pcouple_after_coordinates(FILE* fplog, - int64_t step, - const t_inputrec* inputrec, - const t_mdatoms* md, - const matrix pressure, - const matrix forceVirial, - const matrix constraintVirial, - matrix pressureCouplingMu, - t_state* state, - t_nrnb* nrnb, - gmx::BoxDeformation* boxDeformation, - bool scaleCoordinates); +void update_pcouple_after_coordinates(FILE* fplog, + int64_t step, + const t_inputrec* inputrec, + int homenr, + gmx::ArrayRef cFREEZE, + const matrix pressure, + const matrix forceVirial, + const matrix constraintVirial, + matrix pressureCouplingMu, + t_state* state, + t_nrnb* nrnb, + gmx::BoxDeformation* boxDeformation, + bool scaleCoordinates); /* Return TRUE if OK, FALSE in case of Shake Error */ -extern gmx_bool update_randomize_velocities(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - const t_mdatoms* md, - gmx::ArrayRef v, - const gmx::Update* upd, - const gmx::Constraints* constr); +extern bool update_randomize_velocities(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + int homenr, + gmx::ArrayRef cTC, + gmx::ArrayRef invMass, + gmx::ArrayRef v, + const gmx::Update* upd, + const gmx::Constraints* constr); void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector& therm_integral); //NOLINT(google-runtime-references) -void andersen_tcoupl(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - const t_mdatoms* md, - gmx::ArrayRef v, - real rate, - const std::vector& randomize, - gmx::ArrayRef boltzfac); +void andersen_tcoupl(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + int homenr, + gmx::ArrayRef cTC, + gmx::ArrayRef invMass, + gmx::ArrayRef v, + real rate, + const std::vector& randomize, + gmx::ArrayRef boltzfac); void nosehoover_tcoupl(const t_grpopts* opts, const gmx_ekindata_t* ekind, real dt, - double xi[], - double vxi[], + gmx::ArrayRef xi, + gmx::ArrayRef vxi, const t_extmass* MassQ); -void trotter_update(const t_inputrec* ir, - int64_t step, - gmx_ekindata_t* ekind, - const gmx_enerdata_t* enerd, - t_state* state, - const tensor vir, - const t_mdatoms* md, - const t_extmass* MassQ, - gmx::ArrayRef> trotter_seqlist, - int trotter_seqno); +void trotter_update(const t_inputrec* ir, + int64_t step, + gmx_ekindata_t* ekind, + const gmx_enerdata_t* enerd, + t_state* state, + const tensor vir, + int homenr, + gmx::ArrayRef cTC, + gmx::ArrayRef invMass, + const t_extmass* MassQ, + gmx::ArrayRef> trotter_seqlist, + int trotter_seqno); std::array, ettTSEQMAX> -init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* Mass, gmx_bool bTrotter); +init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* Mass, bool bTrotter); real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ); /* computes all the pressure/tempertature control energy terms to get a conserved energy */ -void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[]); +void vrescale_tcoupl(const t_inputrec* ir, + int64_t step, + gmx_ekindata_t* ekind, + real dt, + gmx::ArrayRef therm_integral); /* Compute temperature scaling. For V-rescale it is done in update. */ -void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[]); +void rescale_velocities(const gmx_ekindata_t* ekind, + gmx::ArrayRef cTC, + int start, + int end, + gmx::ArrayRef v); /* Rescale the velocities with the scaling factor in ekind */ //! Check whether we do simulated annealing. @@ -183,7 +202,7 @@ void parrinellorahman_pcoupl(FILE* fplog, tensor boxv, tensor M, matrix mu, - gmx_bool bFirstStep); + bool bFirstStep); /*! \brief Calculate the pressure coupling scaling matrix * @@ -211,17 +230,17 @@ void pressureCouplingCalculateScalingMatrix(FILE* fplog, * coupling algorithm. */ template -void pressureCouplingScaleBoxAndCoordinates(const t_inputrec* ir, - const matrix mu, - matrix box, - matrix box_rel, - int start, - int nr_atoms, - rvec x[], - rvec v[], - const unsigned short cFREEZE[], - t_nrnb* nrnb, - bool scaleCoordinates); +void pressureCouplingScaleBoxAndCoordinates(const t_inputrec* ir, + const matrix mu, + matrix box, + matrix box_rel, + int start, + int nr_atoms, + gmx::ArrayRef x, + gmx::ArrayRef v, + gmx::ArrayRef cFREEZE, + t_nrnb* nrnb, + bool scaleCoordinates); void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir); diff --git a/src/gromacs/mdlib/tgroup.cpp b/src/gromacs/mdlib/tgroup.cpp index f86fc068b1..86f9d5b905 100644 --- a/src/gromacs/mdlib/tgroup.cpp +++ b/src/gromacs/mdlib/tgroup.cpp @@ -46,7 +46,7 @@ #include "gromacs/mdtypes/inputrec.h" -real sum_ekin(const t_grpopts* opts, gmx_ekindata_t* ekind, real* dekindlambda, gmx_bool bEkinAveVel, gmx_bool bScaleEkin) +real sum_ekin(const t_grpopts* opts, gmx_ekindata_t* ekind, real* dekindlambda, bool bEkinAveVel, bool bScaleEkin) { int i, j, m, ngtc; real T; diff --git a/src/gromacs/mdlib/update_vv.cpp b/src/gromacs/mdlib/update_vv.cpp index fe874ca89b..57d1a14f96 100644 --- a/src/gromacs/mdlib/update_vv.cpp +++ b/src/gromacs/mdlib/update_vv.cpp @@ -134,7 +134,19 @@ void integrateVVFirstStep(int64_t step, 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); + trotter_update(ir, + step, + ekind, + enerd, + state, + total_vir, + mdatoms->homenr, + mdatoms->cTC ? gmx::arrayRefFromArray(mdatoms->cTC, mdatoms->nr) + : gmx::ArrayRef(), + gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr), + MassQ, + trotter_seq, + ettTSEQ1); } upd->update_coords(*ir, @@ -233,7 +245,19 @@ void integrateVVFirstStep(int64_t step, 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); + trotter_update(ir, + step, + ekind, + enerd, + state, + total_vir, + mdatoms->homenr, + mdatoms->cTC ? gmx::arrayRefFromArray(mdatoms->cTC, mdatoms->nr) + : gmx::ArrayRef(), + gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr), + MassQ, + trotter_seq, + ettTSEQ2); /* TODO This is only needed when we're about to write * a checkpoint, because we use it after the restart @@ -442,7 +466,19 @@ void integrateVVSecondStep(int64_t step, bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE); wallcycle_start(wcycle, WallCycleCounter::Update); - trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ, trotter_seq, ettTSEQ4); + trotter_update(ir, + step, + ekind, + enerd, + state, + total_vir, + mdatoms->homenr, + mdatoms->cTC ? gmx::arrayRefFromArray(mdatoms->cTC, mdatoms->nr) + : gmx::ArrayRef(), + gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr), + MassQ, + trotter_seq, + ettTSEQ4); /* now we know the scaling, we can compute the positions again */ std::copy(cbuf->begin(), cbuf->end(), state->x.begin()); diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index a2226a8c85..789b892cba 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -1396,7 +1396,16 @@ void gmx::LegacySimulator::do_md() 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, md, state->v, &upd, constr); + bIfRandomize = update_randomize_velocities(ir, + step, + cr, + md->homenr, + md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr) + : gmx::ArrayRef(), + gmx::arrayRefFromArray(md->invmass, md->nr), + state->v, + &upd, + constr); /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */ if (constr && bIfRandomize) { @@ -1419,7 +1428,19 @@ void gmx::LegacySimulator::do_md() /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */ if (bTrotter) { - trotter_update(ir, step, ekind, enerd, state, total_vir, md, &MassQ, trotter_seq, ettTSEQ3); + trotter_update(ir, + step, + ekind, + enerd, + state, + total_vir, + md->homenr, + md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr) + : gmx::ArrayRef(), + gmx::arrayRefFromArray(md->invmass, md->nr), + &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 @@ -1428,7 +1449,14 @@ void gmx::LegacySimulator::do_md() } else { - update_tcouple(step, ir, state, ekind, &MassQ, md); + update_tcouple(step, + ir, + state, + ekind, + &MassQ, + md->homenr, + md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr) + : gmx::ArrayRef()); update_pcouple_before_coordinates(fplog, step, ir, state, pressureCouplingMu, M, bInitStep); } @@ -1728,8 +1756,20 @@ void gmx::LegacySimulator::do_md() accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals); } - update_pcouple_after_coordinates( - fplog, step, ir, md, pres, force_vir, shake_vir, pressureCouplingMu, state, nrnb, upd.deform(), !useGpuForUpdate); + update_pcouple_after_coordinates(fplog, + step, + ir, + md->homenr, + md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr) + : gmx::ArrayRef(), + pres, + force_vir, + shake_vir, + pressureCouplingMu, + state, + nrnb, + upd.deform(), + !useGpuForUpdate); const bool doBerendsenPressureCoupling = (inputrec->epc == PressureCoupling::Berendsen && do_per_step(step, inputrec->nstpcouple)); -- 2.22.0