From fd323f81addb515cf68de5bebfbe8388abcf0b0b Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Thu, 25 Mar 2021 06:11:17 +0000 Subject: [PATCH] Use ArrayRef in signatures of constraints and some pulling In order to use ArrayRef in constr, settle, shake, and lincs, some pull code also needed to be refactored. Part of preparation for refactoring of t_mdatoms. --- src/gromacs/applied_forces/awh/awh.cpp | 2 +- src/gromacs/applied_forces/awh/awh.h | 2 +- src/gromacs/domdec/mdsetup.cpp | 10 +- src/gromacs/gmxpreprocess/grompp.cpp | 8 +- src/gromacs/gmxpreprocess/readir.h | 7 +- src/gromacs/gmxpreprocess/readpull.cpp | 26 ++-- src/gromacs/mdlib/constr.cpp | 64 +++++---- src/gromacs/mdlib/constr.h | 16 +-- src/gromacs/mdlib/lincs.cpp | 28 ++-- src/gromacs/mdlib/lincs.h | 9 +- src/gromacs/mdlib/settle.cpp | 12 +- src/gromacs/mdlib/settle.h | 10 +- src/gromacs/mdlib/shake.cpp | 13 +- src/gromacs/mdlib/shake.h | 26 ++-- src/gromacs/mdlib/sim_util.cpp | 25 ++-- src/gromacs/mdlib/tests/constrtestrunners.cpp | 8 +- src/gromacs/mdlib/tests/settletestrunners.cpp | 8 +- src/gromacs/mdlib/tests/shake.cpp | 4 +- src/gromacs/mdrun/md.cpp | 9 +- src/gromacs/pulling/pull.cpp | 134 ++++++++++-------- src/gromacs/pulling/pull.h | 87 ++++++------ src/gromacs/pulling/pullutil.cpp | 109 ++++++++------ 22 files changed, 329 insertions(+), 288 deletions(-) diff --git a/src/gromacs/applied_forces/awh/awh.cpp b/src/gromacs/applied_forces/awh/awh.cpp index d987b09f2c..4a5e5aeeb4 100644 --- a/src/gromacs/applied_forces/awh/awh.cpp +++ b/src/gromacs/applied_forces/awh/awh.cpp @@ -289,7 +289,7 @@ bool Awh::isOutputStep(int64_t step) const } real Awh::applyBiasForcesAndUpdateBias(PbcType pbcType, - const real* masses, + ArrayRef masses, ArrayRef neighborLambdaEnergies, ArrayRef neighborLambdaDhdl, const matrix box, diff --git a/src/gromacs/applied_forces/awh/awh.h b/src/gromacs/applied_forces/awh/awh.h index c1d3632b02..7529537394 100644 --- a/src/gromacs/applied_forces/awh/awh.h +++ b/src/gromacs/applied_forces/awh/awh.h @@ -186,7 +186,7 @@ public: * \returns the potential energy for the bias. */ real applyBiasForcesAndUpdateBias(PbcType pbcType, - const real* masses, + ArrayRef masses, ArrayRef neighborLambdaEnergies, ArrayRef neighborLambdaDhdl, const matrix box, diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp index ccf793b0bd..0baae937e7 100644 --- a/src/gromacs/domdec/mdsetup.cpp +++ b/src/gromacs/domdec/mdsetup.cpp @@ -52,6 +52,7 @@ #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/topology/topology.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/gmxassert.h" #include "gromacs/utility/smalloc.h" @@ -105,7 +106,7 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, numHomeAtoms, mdAtoms); - auto mdatoms = mdAtoms->mdatoms(); + t_mdatoms* mdatoms = mdAtoms->mdatoms(); if (usingDomDec) { dd_sort_local_top(*cr->dd, mdatoms, top); @@ -150,11 +151,12 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, constr->setConstraints(top, mdatoms->nr, mdatoms->homenr, - mdatoms->massT, - mdatoms->invmass, + gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr), + gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr), mdatoms->nMassPerturbed != 0, mdatoms->lambda, - mdatoms->cFREEZE); + mdatoms->cFREEZE ? gmx::arrayRefFromArray(mdatoms->cFREEZE, mdatoms->nr) + : gmx::ArrayRef()); } } diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index a6d17f12b0..85f6e376cc 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -2391,12 +2391,8 @@ int gmx_grompp(int argc, char* argv[]) if (ir->bPull) { - pull = set_pull_init(ir, - sys, - state.x.rvec_array(), - state.box, - state.lambda[FreeEnergyPerturbationCouplingType::Mass], - wi); + pull = set_pull_init( + ir, sys, state.x, state.box, state.lambda[FreeEnergyPerturbationCouplingType::Mass], wi); } /* Modules that supply external potential for pull coordinates diff --git a/src/gromacs/gmxpreprocess/readir.h b/src/gromacs/gmxpreprocess/readir.h index dfec8859a6..da326072f5 100644 --- a/src/gromacs/gmxpreprocess/readir.h +++ b/src/gromacs/gmxpreprocess/readir.h @@ -165,7 +165,12 @@ void checkPullCoords(gmx::ArrayRef pullGroups, gmx::ArrayRef pullCoords); /* Process the pull coordinates after reading the pull groups */ -pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix box, real lambda, warninp_t wi); +pull_t* set_pull_init(t_inputrec* ir, + const gmx_mtop_t& mtop, + gmx::ArrayRef x, + matrix box, + real lambda, + warninp_t wi); /* Prints the initial pull group distances in x. * If requested, adds the current distance to the initial reference location. * Returns the pull_t pull work struct. This should be passed to finish_pull() diff --git a/src/gromacs/gmxpreprocess/readpull.cpp b/src/gromacs/gmxpreprocess/readpull.cpp index 5af3b4b9e6..2f8e417c8a 100644 --- a/src/gromacs/gmxpreprocess/readpull.cpp +++ b/src/gromacs/gmxpreprocess/readpull.cpp @@ -549,7 +549,12 @@ void checkPullCoords(gmx::ArrayRef pullGroups, gmx::ArrayRef } } -pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix box, real lambda, warninp_t wi) +pull_t* set_pull_init(t_inputrec* ir, + const gmx_mtop_t& mtop, + gmx::ArrayRef x, + matrix box, + real lambda, + warninp_t wi) { pull_t* pull_work; t_pbc pbc; @@ -573,18 +578,14 @@ pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix bo if (pull->bSetPbcRefToPrevStepCOM) { - initPullComFromPrevStep(nullptr, pull_work, md->massT, &pbc, x); + initPullComFromPrevStep(nullptr, pull_work, gmx::arrayRefFromArray(md->massT, md->nr), &pbc, x); } - pull_calc_coms(nullptr, pull_work, md->massT, &pbc, t_start, x, nullptr); + pull_calc_coms(nullptr, pull_work, gmx::arrayRefFromArray(md->massT, md->nr), &pbc, t_start, x, {}); for (int g = 0; g < pull->ngroup; g++) { - bool groupObeysPbc = pullCheckPbcWithinGroup( - *pull_work, - gmx::arrayRefFromArray(reinterpret_cast(x), mtop.natoms), - pbc, - g, - c_pullGroupSmallGroupThreshold); + bool groupObeysPbc = + pullCheckPbcWithinGroup(*pull_work, x, pbc, g, c_pullGroupSmallGroupThreshold); if (!groupObeysPbc) { char buf[STRLEN]; @@ -616,12 +617,7 @@ pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t& mtop, rvec* x, matrix bo } if (groupObeysPbc) { - groupObeysPbc = pullCheckPbcWithinGroup( - *pull_work, - gmx::arrayRefFromArray(reinterpret_cast(x), mtop.natoms), - pbc, - g, - c_pullGroupPbcMargin); + groupObeysPbc = pullCheckPbcWithinGroup(*pull_work, x, pbc, g, c_pullGroupPbcMargin); if (!groupObeysPbc) { char buf[STRLEN]; diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 1cb575ad6a..d8cddba26e 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -59,6 +59,7 @@ #include "gromacs/fileio/gmxfio.h" #include "gromacs/fileio/pdbio.h" #include "gromacs/gmxlib/nrnb.h" +#include "gromacs/math/arrayrefwithpadding.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" @@ -75,6 +76,7 @@ #include "gromacs/topology/ifunc.h" #include "gromacs/topology/mtop_lookup.h" #include "gromacs/topology/mtop_util.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/gmxassert.h" @@ -109,14 +111,14 @@ public: int numConstraints, int numSettles); ~Impl(); - void setConstraints(gmx_localtop_t* top, - int numAtoms, - int numHomeAtoms, - real* masses, - real* inverseMasses, - bool hasMassPerturbedAtoms, - real lambda, - unsigned short* cFREEZE); + void setConstraints(gmx_localtop_t* top, + int numAtoms, + int numHomeAtoms, + gmx::ArrayRef masses, + gmx::ArrayRef inverseMasses, + bool hasMassPerturbedAtoms, + real lambda, + gmx::ArrayRef cFREEZE); bool apply(bool bLog, bool bEner, int64_t step, @@ -169,15 +171,15 @@ public: //! Number of local atoms. int numHomeAtoms_ = 0; //! Atoms masses. - real* masses_; + gmx::ArrayRef masses_; //! Inverse masses. - real* inverseMasses_; + gmx::ArrayRef inverseMasses_; //! If there are atoms with perturbed mass (for FEP). bool hasMassPerturbedAtoms_; //! FEP lambda value. real lambda_; //! Freeze groups data - unsigned short* cFREEZE_; + gmx::ArrayRef cFREEZE_; //! Whether we need to do pbc for handling bonds. bool pbcHandlingRequired_ = false; @@ -768,9 +770,9 @@ bool Constraints::Impl::apply(bool bLog, cr, ir.delta_t, t, - as_rvec_array(x.unpaddedArrayRef().data()), - as_rvec_array(xprime.unpaddedArrayRef().data()), - as_rvec_array(v.unpaddedArrayRef().data()), + x.unpaddedArrayRef(), + xprime.unpaddedArrayRef(), + v.unpaddedArrayRef(), constraintsVirial); } if (ed && delta_step > 0) @@ -788,7 +790,7 @@ bool Constraints::Impl::apply(bool bLog, wallcycle_stop(wcycle, ewcCONSTR); const bool haveVelocities = (!v.empty() || econq == ConstraintVariable::Velocities); - if (haveVelocities && cFREEZE_) + if (haveVelocities && !cFREEZE_.empty()) { /* Set the velocities of frozen dimensions to zero */ ArrayRef vRef; @@ -993,14 +995,14 @@ static std::vector make_at2settle(int natoms, const InteractionList& ilist) return at2s; } -void Constraints::Impl::setConstraints(gmx_localtop_t* top, - int numAtoms, - int numHomeAtoms, - real* masses, - real* inverseMasses, - const bool hasMassPerturbedAtoms, - const real lambda, - unsigned short* cFREEZE) +void Constraints::Impl::setConstraints(gmx_localtop_t* top, + int numAtoms, + int numHomeAtoms, + gmx::ArrayRef masses, + gmx::ArrayRef inverseMasses, + const bool hasMassPerturbedAtoms, + const real lambda, + gmx::ArrayRef cFREEZE) { numAtoms_ = numAtoms; numHomeAtoms_ = numHomeAtoms; @@ -1050,14 +1052,14 @@ void Constraints::Impl::setConstraints(gmx_localtop_t* top, } } -void Constraints::setConstraints(gmx_localtop_t* top, - const int numAtoms, - const int numHomeAtoms, - real* masses, - real* inverseMasses, - const bool hasMassPerturbedAtoms, - const real lambda, - unsigned short* cFREEZE) +void Constraints::setConstraints(gmx_localtop_t* top, + const int numAtoms, + const int numHomeAtoms, + gmx::ArrayRef masses, + gmx::ArrayRef inverseMasses, + const bool hasMassPerturbedAtoms, + const real lambda, + gmx::ArrayRef cFREEZE) { impl_->setConstraints( top, numAtoms, numHomeAtoms, masses, inverseMasses, hasMassPerturbedAtoms, lambda, cFREEZE); diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h index eb51ffb602..03eb4772ac 100644 --- a/src/gromacs/mdlib/constr.h +++ b/src/gromacs/mdlib/constr.h @@ -132,14 +132,14 @@ public: * * \todo Make this a callback that is called automatically * once a new domain has been made. */ - void setConstraints(gmx_localtop_t* top, - int numAtoms, - int numHomeAtoms, - real* masses, - real* inverseMasses, - bool hasMassPerturbedAtoms, - real lambda, - unsigned short* cFREEZE); + void setConstraints(gmx_localtop_t* top, + int numAtoms, + int numHomeAtoms, + gmx::ArrayRef masses, + gmx::ArrayRef inverseMasses, + bool hasMassPerturbedAtoms, + real lambda, + gmx::ArrayRef cFREEZE); /*! \brief Applies constraints to coordinates. * diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp index b5fe54be92..1bbfd8330e 100644 --- a/src/gromacs/mdlib/lincs.cpp +++ b/src/gromacs/mdlib/lincs.cpp @@ -349,10 +349,10 @@ static void lincs_update_atoms_noind(int ncons, real preFactor, gmx::ArrayRef fac, gmx::ArrayRef r, - const real* invmass, + gmx::ArrayRef invmass, rvec* x) { - if (invmass != nullptr) + if (!invmass.empty()) { for (int b = 0; b < ncons; b++) { @@ -398,10 +398,10 @@ static void lincs_update_atoms_ind(gmx::ArrayRef ind, real preFactor, gmx::ArrayRef fac, gmx::ArrayRef r, - const real* invmass, + gmx::ArrayRef invmass, rvec* x) { - if (invmass != nullptr) + if (!invmass.empty()) { for (int b : ind) { @@ -447,7 +447,7 @@ static void lincs_update_atoms(Lincs* li, real preFactor, gmx::ArrayRef fac, gmx::ArrayRef r, - const real* invmass, + gmx::ArrayRef invmass, rvec* x) { if (li->ntask == 1) @@ -574,7 +574,7 @@ static void do_lincsp(ArrayRefWithPadding xPadded, t_pbc* pbc, Lincs* lincsd, int th, - const real* invmass, + ArrayRef invmass, ConstraintVariable econq, bool bCalcDHDL, bool bCalcVir, @@ -713,7 +713,7 @@ static void do_lincsp(ArrayRefWithPadding xPadded, 1.0, sol, r, - (econq != ConstraintVariable::Force) ? invmass : nullptr, + (econq != ConstraintVariable::Force) ? invmass : gmx::ArrayRef(), as_rvec_array(fp.data())); if (bCalcDHDL) @@ -958,7 +958,7 @@ static void do_lincs(ArrayRefWithPadding xPadded, t_pbc* pbc, Lincs* lincsd, int th, - const real* invmass, + ArrayRef invmass, const t_commrec* cr, bool bCalcDHDL, real wangle, @@ -1217,7 +1217,11 @@ static void do_lincs(ArrayRefWithPadding xPadded, } /*! \brief Sets the elements in the LINCS matrix for task task. */ -static void set_lincs_matrix_task(Lincs* li, Task* li_task, const real* invmass, int* ncc_triangle, int* nCrossTaskTriangles) +static void set_lincs_matrix_task(Lincs* li, + Task* li_task, + ArrayRef invmass, + int* ncc_triangle, + int* nCrossTaskTriangles) { /* Construct the coupling coefficient matrix blmf */ li_task->ntriangle = 0; @@ -1305,7 +1309,7 @@ static void set_lincs_matrix_task(Lincs* li, Task* li_task, const real* invmass, } /*! \brief Sets the elements in the LINCS matrix. */ -static void set_lincs_matrix(Lincs* li, const real* invmass, real lambda) +static void set_lincs_matrix(Lincs* li, ArrayRef invmass, real lambda) { const real invsqrt2 = 0.7071067811865475244; @@ -1854,7 +1858,7 @@ static void set_matrix_indices(Lincs* li, const Task& li_task, const ListOfLists void set_lincs(const InteractionDefinitions& idef, const int numAtoms, - const real* invmass, + ArrayRef invmass, const real lambda, bool bDynamics, const t_commrec* cr, @@ -2268,7 +2272,7 @@ bool constrain_lincs(bool computeRmsd, const t_inputrec& ir, int64_t step, Lincs* lincsd, - const real* invmass, + ArrayRef invmass, const t_commrec* cr, const gmx_multisim_t* ms, ArrayRefWithPadding xPadded, diff --git a/src/gromacs/mdlib/lincs.h b/src/gromacs/mdlib/lincs.h index 0b00370c01..fcaf91a985 100644 --- a/src/gromacs/mdlib/lincs.h +++ b/src/gromacs/mdlib/lincs.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, 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. @@ -61,9 +61,10 @@ struct t_pbc; namespace gmx { - template class ArrayRefWithPadding; +template +class ArrayRef; enum class ConstraintVariable : int; class Lincs; template @@ -91,7 +92,7 @@ void done_lincs(Lincs* li); /*! \brief Initialize lincs stuff */ void set_lincs(const InteractionDefinitions& idef, int numAtoms, - const real* invmass, + ArrayRef invmass, real lambda, bool bDynamics, const t_commrec* cr, @@ -104,7 +105,7 @@ bool constrain_lincs(bool computeRmsd, const t_inputrec& ir, int64_t step, Lincs* lincsd, - const real* invmass, + ArrayRef invmass, const t_commrec* cr, const gmx_multisim_t* ms, ArrayRefWithPadding x, diff --git a/src/gromacs/mdlib/settle.cpp b/src/gromacs/mdlib/settle.cpp index 43e21eb486..0ccadf81e9 100644 --- a/src/gromacs/mdlib/settle.cpp +++ b/src/gromacs/mdlib/settle.cpp @@ -4,7 +4,7 @@ * 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 + * Copyright (c) 2018,2019,2020,2021, 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. @@ -57,7 +57,6 @@ #include "gromacs/math/invertmatrix.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/constr.h" -#include "gromacs/mdtypes/mdatom.h" #include "gromacs/pbcutil/ishift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pbcutil/pbc_simd.h" @@ -66,6 +65,7 @@ #include "gromacs/topology/idef.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/mtop_util.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/gmxassert.h" @@ -196,10 +196,10 @@ SettleData::SettleData(const gmx_mtop_t& mtop) : parametersAllMasses1_ = settleParameters(1.0, 1.0, 1.0, 1.0, dOH, dHH); } -void SettleData::setConstraints(const InteractionList& il_settle, - const int numHomeAtoms, - const real* masses, - const real* inverseMasses) +void SettleData::setConstraints(const InteractionList& il_settle, + const int numHomeAtoms, + gmx::ArrayRef masses, + gmx::ArrayRef inverseMasses) { #if GMX_SIMD_HAVE_REAL const int pack_size = GMX_SIMD_REAL_WIDTH; diff --git a/src/gromacs/mdlib/settle.h b/src/gromacs/mdlib/settle.h index 4714bd669f..86bf81fbf1 100644 --- a/src/gromacs/mdlib/settle.h +++ b/src/gromacs/mdlib/settle.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, 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. @@ -121,10 +121,10 @@ public: SettleData(const gmx_mtop_t& mtop); //! Sets the constraints from the interaction list and the masses - void setConstraints(const InteractionList& il_settle, - int numHomeAtoms, - const real* masses, - const real* inverseMasses); + void setConstraints(const InteractionList& il_settle, + int numHomeAtoms, + gmx::ArrayRef masses, + gmx::ArrayRef inverseMasses); //! Returns settle parameters for constraining coordinates and forces const SettleParameters& parametersMassWeighted() const { return parametersMassWeighted_; } diff --git a/src/gromacs/mdlib/shake.cpp b/src/gromacs/mdlib/shake.cpp index 80e2deee7e..47ef4167f2 100644 --- a/src/gromacs/mdlib/shake.cpp +++ b/src/gromacs/mdlib/shake.cpp @@ -61,6 +61,7 @@ #include "gromacs/mdtypes/md_enums.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/invblock.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" @@ -282,7 +283,7 @@ void cshake(const int iatom[], ArrayRef initial_displacements, ArrayRef half_of_reduced_mass, real omega, - const real invmass[], + ArrayRef invmass, ArrayRef distance_squared_tolerance, ArrayRef scaled_lagrange_multiplier, int* nerror) @@ -371,7 +372,7 @@ static void crattle(const int iatom[], ArrayRef rij, ArrayRef m2, real omega, - const real invmass[], + ArrayRef invmass, ArrayRef distance_squared_tolerance, ArrayRef scaled_lagrange_multiplier, int* nerror, @@ -440,7 +441,7 @@ static void crattle(const int iatom[], //! Applies SHAKE static int vec_shakef(FILE* fplog, shakedata* shaked, - const real invmass[], + ArrayRef invmass, int ncon, ArrayRef ip, const int* iatom, @@ -636,7 +637,7 @@ static void check_cons(FILE* log, const t_pbc* pbc, ArrayRef ip, const int* iatom, - const real invmass[], + ArrayRef invmass, ConstraintVariable econq) { int ai, aj; @@ -699,7 +700,7 @@ static void check_cons(FILE* log, //! Applies SHAKE. static bool bshakef(FILE* log, shakedata* shaked, - const real invmass[], + ArrayRef invmass, const InteractionDefinitions& idef, const t_inputrec& ir, ArrayRef x_s, @@ -821,7 +822,7 @@ static bool bshakef(FILE* log, bool constrain_shake(FILE* log, shakedata* shaked, - const real invmass[], + ArrayRef invmass, const InteractionDefinitions& idef, const t_inputrec& ir, ArrayRef x_s, diff --git a/src/gromacs/mdlib/shake.h b/src/gromacs/mdlib/shake.h index 8559c57553..f0723818bc 100644 --- a/src/gromacs/mdlib/shake.h +++ b/src/gromacs/mdlib/shake.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, 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. @@ -108,17 +108,17 @@ void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon); * sblock[n] to sblock[n+1]. Array sblock should be large enough. * Return TRUE when OK, FALSE when shake-error */ -bool constrain_shake(FILE* log, /* Log file */ - shakedata* shaked, /* Total number of atoms */ - const real invmass[], /* Atomic masses */ - const InteractionDefinitions& idef, /* The interaction def */ - const t_inputrec& ir, /* Input record */ - ArrayRef x_s, /* Coords before update */ - ArrayRef xprime, /* Output coords when constraining x */ - ArrayRef vprime, /* Output coords when constraining v */ - const t_pbc* pbc, /* PBC information */ - t_nrnb* nrnb, /* Performance measure */ - real lambda, /* FEP lambda */ +bool constrain_shake(FILE* log, /* Log file */ + shakedata* shaked, /* Total number of atoms */ + gmx::ArrayRef invmass, /* Atomic masses */ + const InteractionDefinitions& idef, /* The interaction def */ + const t_inputrec& ir, /* Input record */ + ArrayRef x_s, /* Coords before update */ + ArrayRef xprime, /* Output coords when constraining x */ + ArrayRef vprime, /* Output coords when constraining v */ + const t_pbc* pbc, /* PBC information */ + t_nrnb* nrnb, /* Performance measure */ + real lambda, /* FEP lambda */ real* dvdlambda, /* FEP force */ real invdt, /* 1/delta_t */ ArrayRef v, /* Also constrain v if not empty */ @@ -138,7 +138,7 @@ void cshake(const int iatom[], ArrayRef rij, ArrayRef half_of_reduced_mass, real omega, - const real invmass[], + ArrayRef invmass, ArrayRef distance_squared_tolerance, ArrayRef scaled_lagrange_multiplier, int* nerror); diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 6f42594ed0..424e9844fc 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -205,12 +205,12 @@ static void pull_potential_wrapper(const t_commrec* cr, dvdl = 0; enerd->term[F_COM_PULL] += pull_potential(pull_work, - mdatoms->massT, + gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr), &pbc, cr, t, lambda[static_cast(FreeEnergyPerturbationCouplingType::Restraint)], - as_rvec_array(x.data()), + x, force, &dvdl); enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Restraint] += dvdl; @@ -680,16 +680,17 @@ static void computeSpecialForces(FILE* fplog, } auto& forceWithVirial = (mtsLevel == 0) ? forceWithVirialMtsLevel0 : forceWithVirialMtsLevel1; - enerd->term[F_COM_PULL] += awh->applyBiasForcesAndUpdateBias(inputrec.pbcType, - mdatoms->massT, - foreignLambdaDeltaH, - foreignLambdaDhDl, - box, - forceWithVirial, - t, - step, - wcycle, - fplog); + enerd->term[F_COM_PULL] += awh->applyBiasForcesAndUpdateBias( + inputrec.pbcType, + gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr), + foreignLambdaDeltaH, + foreignLambdaDhDl, + box, + forceWithVirial, + t, + step, + wcycle, + fplog); } } diff --git a/src/gromacs/mdlib/tests/constrtestrunners.cpp b/src/gromacs/mdlib/tests/constrtestrunners.cpp index d8b0bfc53c..32ebfbeb55 100644 --- a/src/gromacs/mdlib/tests/constrtestrunners.cpp +++ b/src/gromacs/mdlib/tests/constrtestrunners.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, 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. @@ -86,7 +86,7 @@ void ShakeConstraintsRunner::applyConstraints(ConstraintsTestData* testData, t_p make_shake_sblock_serial(&shaked, testData->idef_.get(), testData->numAtoms_); bool success = constrain_shake(nullptr, &shaked, - testData->invmass_.data(), + testData->invmass_, *testData->idef_, testData->ir_, testData->x_, @@ -141,7 +141,7 @@ void LincsConstraintsRunner::applyConstraints(ConstraintsTestData* testData, t_p testData->ir_.nProjOrder); set_lincs(*testData->idef_, testData->numAtoms_, - testData->invmass_.data(), + testData->invmass_, testData->lambda_, EI_DYNAMICS(testData->ir_.eI), &cr, @@ -152,7 +152,7 @@ void LincsConstraintsRunner::applyConstraints(ConstraintsTestData* testData, t_p testData->ir_, 0, lincsd, - testData->invmass_.data(), + testData->invmass_, &cr, &ms, testData->x_.arrayRefWithPadding(), diff --git a/src/gromacs/mdlib/tests/settletestrunners.cpp b/src/gromacs/mdlib/tests/settletestrunners.cpp index 6dc80bf66e..c494341cc2 100644 --- a/src/gromacs/mdlib/tests/settletestrunners.cpp +++ b/src/gromacs/mdlib/tests/settletestrunners.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020,2021, 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. @@ -65,10 +65,8 @@ void SettleHostTestRunner::applySettle(SettleTestData* testData, { SettleData settled(testData->mtop_); - settled.setConstraints(testData->idef_->il[F_SETTLE], - testData->numAtoms_, - testData->masses_.data(), - testData->inverseMasses_.data()); + settled.setConstraints( + testData->idef_->il[F_SETTLE], testData->numAtoms_, testData->masses_, testData->inverseMasses_); bool errorOccured; int numThreads = 1; diff --git a/src/gromacs/mdlib/tests/shake.cpp b/src/gromacs/mdlib/tests/shake.cpp index 50ef3e0fb2..7f58294382 100644 --- a/src/gromacs/mdlib/tests/shake.cpp +++ b/src/gromacs/mdlib/tests/shake.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2014,2015,2017,2018,2019,2020, by the GROMACS development team, led by + * Copyright (c) 2014,2015,2017,2018,2019,2020,2021, 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. @@ -205,7 +205,7 @@ public: initialDisplacements, halfOfReducedMasses, omega_, - inverseMasses.data(), + inverseMasses, distanceSquaredTolerances, lagrangianValues, &numErrors); diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 77d3641a6d..5584b60232 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -515,8 +515,13 @@ void gmx::LegacySimulator::do_md() EnergyData::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput); } - preparePrevStepPullCom( - ir, pull_work, mdatoms->massT, state, state_global, cr, startingBehavior != StartingBehavior::NewSimulation); + preparePrevStepPullCom(ir, + pull_work, + gmx::arrayRefFromArray(mdatoms->massT, mdatoms->nr), + state, + state_global, + cr, + startingBehavior != StartingBehavior::NewSimulation); // TODO: Remove this by converting AWH into a ForceProvider auto awh = prepareAwhModule(fplog, diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp index ae29d66801..5a5538ccb9 100644 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@ -37,7 +37,6 @@ */ #include "gmxpre.h" -#include "gromacs/utility/stringutil.h" #include "pull.h" #include "config.h" @@ -71,6 +70,7 @@ #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/mtop_lookup.h" #include "gromacs/topology/topology.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/cstringutil.h" #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" @@ -80,6 +80,7 @@ #include "gromacs/utility/real.h" #include "gromacs/utility/smalloc.h" #include "gromacs/utility/strconvert.h" +#include "gromacs/utility/stringutil.h" #include "pull_internal.h" @@ -175,13 +176,13 @@ double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd) } /* Apply forces in a mass weighted fashion for part of the pull group */ -static void apply_forces_grp_part(const pull_group_work_t* pgrp, - int ind_start, - int ind_end, - const real* masses, - const dvec f_pull, - int sign, - rvec* f) +static void apply_forces_grp_part(const pull_group_work_t* pgrp, + int ind_start, + int ind_end, + gmx::ArrayRef masses, + const dvec f_pull, + int sign, + rvec* f) { double inv_wm = pgrp->mwscale; @@ -203,12 +204,12 @@ static void apply_forces_grp_part(const pull_group_work_t* pgrp, } /* Apply forces in a mass weighted fashion */ -static void apply_forces_grp(const pull_group_work_t* pgrp, - const real* masses, - const dvec f_pull, - int sign, - rvec* f, - int nthreads) +static void apply_forces_grp(const pull_group_work_t* pgrp, + gmx::ArrayRef masses, + const dvec f_pull, + int sign, + rvec* f, + int nthreads) { auto localAtomIndices = pgrp->atomSet.localIndex(); @@ -243,13 +244,13 @@ static void apply_forces_grp(const pull_group_work_t* pgrp, } /* Apply forces in a mass weighted fashion to a cylinder group */ -static void apply_forces_cyl_grp(const pull_group_work_t* pgrp, - const double dv_corr, - const real* masses, - const dvec f_pull, - double f_scal, - int sign, - rvec* f, +static void apply_forces_cyl_grp(const pull_group_work_t* pgrp, + const double dv_corr, + gmx::ArrayRef masses, + const dvec f_pull, + double f_scal, + int sign, + rvec* f, int gmx_unused nthreads) { double inv_wm = pgrp->mwscale; @@ -290,10 +291,10 @@ static void apply_forces_cyl_grp(const pull_group_work_t* pgrp, /* Apply torque forces in a mass weighted fashion to the groups that define * the pull vector direction for pull coordinate pcrd. */ -static void apply_forces_vec_torque(const struct pull_t* pull, - const pull_coord_work_t* pcrd, - const real* masses, - rvec* f) +static void apply_forces_vec_torque(const struct pull_t* pull, + const pull_coord_work_t* pcrd, + gmx::ArrayRef masses, + rvec* f) { const PullCoordSpatialData& spatialData = pcrd->spatialData; @@ -326,7 +327,7 @@ static void apply_forces_vec_torque(const struct pull_t* pull, static void apply_forces_coord(struct pull_t* pull, int coord, const PullCoordVectorForces& forces, - const real* masses, + gmx::ArrayRef masses, rvec* f) { /* Here it would be more efficient to use one large thread-parallel @@ -817,8 +818,14 @@ void clear_pull_forces(pull_t* pull) } /* Apply constraint using SHAKE */ -static void -do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t) +static void do_constraint(struct pull_t* pull, + t_pbc* pbc, + gmx::ArrayRef x, + gmx::ArrayRef v, + gmx_bool bMaster, + tensor vir, + double dt, + double t) { dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */ @@ -1111,7 +1118,7 @@ do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaste /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */ - if (v) + if (!v.empty()) { invdt = 1 / dt; } @@ -1147,7 +1154,7 @@ do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaste { x[ii][m] += tmp[m]; } - if (v) + if (!v.empty()) { for (m = 0; m < DIM; m++) { @@ -1532,11 +1539,11 @@ static void check_external_potential_registration(const struct pull_t* pull) * potential energy is added either to the pull term or to a term * specific to the external module. */ -void apply_external_pull_coord_force(struct pull_t* pull, - int coord_index, - double coord_force, - const real* masses, - gmx::ForceWithVirial* forceWithVirial) +void apply_external_pull_coord_force(struct pull_t* pull, + int coord_index, + double coord_force, + gmx::ArrayRef masses, + gmx::ForceWithVirial* forceWithVirial) { pull_coord_work_t* pcrd; @@ -1602,15 +1609,15 @@ static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull, return pullCoordForces; } -real pull_potential(struct pull_t* pull, - const real* masses, - t_pbc* pbc, - const t_commrec* cr, - double t, - real lambda, - const rvec* x, - gmx::ForceWithVirial* force, - real* dvdlambda) +real pull_potential(struct pull_t* pull, + gmx::ArrayRef masses, + t_pbc* pbc, + const t_commrec* cr, + double t, + real lambda, + gmx::ArrayRef x, + gmx::ForceWithVirial* force, + real* dvdlambda) { real V = 0; @@ -1626,7 +1633,7 @@ real pull_potential(struct pull_t* pull, { real dVdl = 0; - pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr); + pull_calc_coms(cr, pull, masses, pbc, t, x, {}); rvec* f = as_rvec_array(force->force_.data()); matrix virial = { { 0 } }; @@ -1666,16 +1673,16 @@ real pull_potential(struct pull_t* pull, return (MASTER(cr) ? V : 0.0); } -void pull_constraint(struct pull_t* pull, - const real* masses, - t_pbc* pbc, - const t_commrec* cr, - double dt, - double t, - rvec* x, - rvec* xp, - rvec* v, - tensor vir) +void pull_constraint(struct pull_t* pull, + gmx::ArrayRef masses, + t_pbc* pbc, + const t_commrec* cr, + double dt, + double t, + gmx::ArrayRef x, + gmx::ArrayRef xp, + gmx::ArrayRef v, + tensor vir) { assert(pull != nullptr); @@ -2400,13 +2407,13 @@ static void destroy_pull(struct pull_t* pull) delete pull; } -void preparePrevStepPullCom(const t_inputrec* ir, - pull_t* pull_work, - const real* masses, - t_state* state, - const t_state* state_global, - const t_commrec* cr, - bool startingFromCheckpoint) +void preparePrevStepPullCom(const t_inputrec* ir, + pull_t* pull_work, + gmx::ArrayRef masses, + t_state* state, + const t_state* state_global, + const t_commrec* cr, + bool startingFromCheckpoint) { if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM) { @@ -2432,7 +2439,8 @@ void preparePrevStepPullCom(const t_inputrec* ir, { t_pbc pbc; set_pbc(&pbc, ir->pbcType, state->box); - initPullComFromPrevStep(cr, pull_work, masses, &pbc, state->x.rvec_array()); + initPullComFromPrevStep( + cr, pull_work, masses, &pbc, state->x.arrayRefWithPadding().unpaddedArrayRef()); updatePrevStepPullCom(pull_work, state); } } diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h index a91784cf83..4016bed9fa 100644 --- a/src/gromacs/pulling/pull.h +++ b/src/gromacs/pulling/pull.h @@ -72,6 +72,8 @@ class t_state; namespace gmx { +template +class ArrayRef; class ForceWithVirial; class LocalAtomSetManager; } // namespace gmx @@ -155,12 +157,11 @@ void register_external_pull_potential(struct pull_t* pull, int coord_index, cons * \param[in] masses Atoms masses. * \param[in,out] forceWithVirial Force and virial buffers. */ -void apply_external_pull_coord_force(struct pull_t* pull, - int coord_index, - double coord_force, - const real* masses, - gmx::ForceWithVirial* forceWithVirial); - +void apply_external_pull_coord_force(struct pull_t* pull, + int coord_index, + double coord_force, + gmx::ArrayRef masses, + gmx::ForceWithVirial* forceWithVirial); /*! \brief Set the all the pull forces to zero. * @@ -183,15 +184,15 @@ void clear_pull_forces(struct pull_t* pull); * * \returns The pull potential energy. */ -real pull_potential(struct pull_t* pull, - const real* masses, - struct t_pbc* pbc, - const t_commrec* cr, - double t, - real lambda, - const rvec* x, - gmx::ForceWithVirial* force, - real* dvdlambda); +real pull_potential(struct pull_t* pull, + gmx::ArrayRef masses, + struct t_pbc* pbc, + const t_commrec* cr, + double t, + real lambda, + gmx::ArrayRef x, + gmx::ForceWithVirial* force, + real* dvdlambda); /*! \brief Constrain the coordinates xp in the directions in x @@ -208,16 +209,16 @@ real pull_potential(struct pull_t* pull, * \param[in,out] v Velocities, which may get a pull correction. * \param[in,out] vir The virial, which, if != NULL, gets a pull correction. */ -void pull_constraint(struct pull_t* pull, - const real* masses, - struct t_pbc* pbc, - const t_commrec* cr, - double dt, - double t, - rvec* x, - rvec* xp, - rvec* v, - tensor vir); +void pull_constraint(struct pull_t* pull, + gmx::ArrayRef masses, + struct t_pbc* pbc, + const t_commrec* cr, + double dt, + double t, + gmx::ArrayRef x, + gmx::ArrayRef xp, + gmx::ArrayRef v, + tensor vir); /*! \brief Make a selection of the home atoms for all pull groups. @@ -266,13 +267,13 @@ void finish_pull(struct pull_t* pull); * \param[in,out] xp Updated x, can be NULL. * */ -void pull_calc_coms(const t_commrec* cr, - pull_t* pull, - const real* masses, - t_pbc* pbc, - double t, - const rvec x[], - rvec* xp); +void pull_calc_coms(const t_commrec* cr, + pull_t* pull, + gmx::ArrayRef masses, + t_pbc* pbc, + double t, + gmx::ArrayRef x, + gmx::ArrayRef xp); /*! \brief Margin for checking pull group PBC distances compared to half the box size */ static constexpr real c_pullGroupPbcMargin = 0.9; @@ -297,7 +298,7 @@ static constexpr real c_pullGroupSmallGroupThreshold = 0.5; * \param[in] pbcMargin The minimum margin (as a fraction) to half the box size * \returns -1 when all groups obey PBC or the first group index that fails PBC */ -int pullCheckPbcWithinGroups(const pull_t& pull, const rvec* x, const t_pbc& pbc, real pbcMargin); +int pullCheckPbcWithinGroups(const pull_t& pull, gmx::ArrayRef x, const t_pbc& pbc, real pbcMargin); /*! \brief Checks whether a specific group that uses a reference atom is within PBC restrictions * @@ -376,13 +377,13 @@ void updatePrevStepPullCom(struct pull_t* pull, t_state* state); * \param[in] cr Struct for communication info. * \param[in] startingFromCheckpoint Is the simulation starting from a checkpoint? */ -void preparePrevStepPullCom(const t_inputrec* ir, - pull_t* pull_work, - const real* masses, - t_state* state, - const t_state* state_global, - const t_commrec* cr, - bool startingFromCheckpoint); +void preparePrevStepPullCom(const t_inputrec* ir, + pull_t* pull_work, + gmx::ArrayRef masses, + t_state* state, + const t_state* state_global, + const t_commrec* cr, + bool startingFromCheckpoint); /*! \brief Initializes the COM of the previous step (set to initial COM) * @@ -392,6 +393,10 @@ void preparePrevStepPullCom(const t_inputrec* ir, * \param[in] pbc Information struct about periodicity. * \param[in] x The local positions. */ -void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, const rvec x[]); +void initPullComFromPrevStep(const t_commrec* cr, + pull_t* pull, + gmx::ArrayRef masses, + t_pbc* pbc, + gmx::ArrayRef x); #endif diff --git a/src/gromacs/pulling/pullutil.cpp b/src/gromacs/pulling/pullutil.cpp index c4985a8f80..b9dbb2d68e 100644 --- a/src/gromacs/pulling/pullutil.cpp +++ b/src/gromacs/pulling/pullutil.cpp @@ -48,13 +48,14 @@ #include "gromacs/math/units.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" +#include "gromacs/math/vectypes.h" #include "gromacs/mdtypes/commrec.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/pbc.h" #include "gromacs/pulling/pull.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/futil.h" #include "gromacs/utility/gmxassert.h" @@ -119,7 +120,7 @@ static void pullAllReduce(const t_commrec* cr, pull_comm_t* comm, int n, T* data /* Copies the coordinates of the PBC atom of pgrp to x_pbc. * When those coordinates are not available on this rank, clears x_pbc. */ -static void setPbcAtomCoords(const pull_group_work_t& pgrp, const rvec* x, rvec x_pbc) +static void setPbcAtomCoords(const pull_group_work_t& pgrp, gmx::ArrayRef x, rvec x_pbc) { if (pgrp.pbcAtomSet != nullptr) { @@ -140,7 +141,10 @@ static void setPbcAtomCoords(const pull_group_work_t& pgrp, const rvec* x, rvec } } -static void pull_set_pbcatoms(const t_commrec* cr, struct pull_t* pull, const rvec* x, gmx::ArrayRef x_pbc) +static void pull_set_pbcatoms(const t_commrec* cr, + struct pull_t* pull, + gmx::ArrayRef x, + gmx::ArrayRef x_pbc) { int numPbcAtoms = 0; for (size_t g = 0; g < pull->group.size(); g++) @@ -167,8 +171,12 @@ static void pull_set_pbcatoms(const t_commrec* cr, struct pull_t* pull, const rv } } -static void -make_cyl_refgrps(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, double t, const rvec* x) +static void make_cyl_refgrps(const t_commrec* cr, + pull_t* pull, + gmx::ArrayRef masses, + t_pbc* pbc, + double t, + gmx::ArrayRef x) { pull_comm_t* comm = &pull->comm; @@ -380,15 +388,15 @@ static double atan2_0_2pi(double y, double x) return a; } -static void sum_com_part(const pull_group_work_t* pgrp, - int ind_start, - int ind_end, - const rvec* x, - const rvec* xp, - const real* mass, - const t_pbc* pbc, - const rvec x_pbc, - ComSums* sum_com) +static void sum_com_part(const pull_group_work_t* pgrp, + int ind_start, + int ind_end, + gmx::ArrayRef x, + gmx::ArrayRef xp, + gmx::ArrayRef mass, + const t_pbc* pbc, + const rvec x_pbc, + ComSums* sum_com) { double sum_wm = 0; double sum_wwm = 0; @@ -421,7 +429,7 @@ static void sum_com_part(const pull_group_work_t* pgrp, { sum_wmx[d] += wm * x[ii][d]; } - if (xp) + if (!xp.empty()) { for (int d = 0; d < DIM; d++) { @@ -439,7 +447,7 @@ static void sum_com_part(const pull_group_work_t* pgrp, { sum_wmx[d] += wm * dx[d]; } - if (xp) + if (!xp.empty()) { /* For xp add the difference between xp and x to dx, * such that we use the same periodic image, @@ -456,21 +464,21 @@ static void sum_com_part(const pull_group_work_t* pgrp, sum_com->sum_wm = sum_wm; sum_com->sum_wwm = sum_wwm; copy_dvec(sum_wmx, sum_com->sum_wmx); - if (xp) + if (!xp.empty()) { copy_dvec(sum_wmxp, sum_com->sum_wmxp); } } -static void sum_com_part_cosweight(const pull_group_work_t* pgrp, - int ind_start, - int ind_end, - int cosdim, - real twopi_box, - const rvec* x, - const rvec* xp, - const real* mass, - ComSums* sum_com) +static void sum_com_part_cosweight(const pull_group_work_t* pgrp, + int ind_start, + int ind_end, + int cosdim, + real twopi_box, + gmx::ArrayRef x, + gmx::ArrayRef xp, + gmx::ArrayRef mass, + ComSums* sum_com) { /* Cosine weighting geometry */ double sum_cm = 0; @@ -496,7 +504,7 @@ static void sum_com_part_cosweight(const pull_group_work_t* pgrp, sum_csm += static_cast(cw * sw * m); sum_ssm += static_cast(sw * sw * m); - if (xp != nullptr) + if (!xp.empty()) { real cw = std::cos(xp[ii][cosdim] * twopi_box); real sw = std::sin(xp[ii][cosdim] * twopi_box); @@ -515,7 +523,13 @@ static void sum_com_part_cosweight(const pull_group_work_t* pgrp, } /* calculates center of mass of selection index from all coordinates x */ -void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, double t, const rvec x[], rvec* xp) +void pull_calc_coms(const t_commrec* cr, + pull_t* pull, + gmx::ArrayRef masses, + t_pbc* pbc, + double t, + gmx::ArrayRef x, + gmx::ArrayRef xp) { real twopi_box = 0; pull_comm_t* comm; @@ -605,9 +619,9 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1 && masses[pgrp->atomSet.localIndex()[0]] == 0) { - GMX_ASSERT(xp == nullptr, + GMX_ASSERT(xp.empty(), "We should not have groups with zero mass with constraints, i.e. " - "xp!=NULL"); + "xp not empty"); /* Copy the single atom coordinate */ for (int d = 0; d < DIM; d++) @@ -743,14 +757,14 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc for (m = 0; m < DIM; m++) { pgrp->x[m] = comBuffer[0][m] * pgrp->mwscale; - if (xp) + if (!xp.empty()) { pgrp->xp[m] = comBuffer[1][m] * pgrp->mwscale; } if (pgrp->epgrppbc == epgrppbcREFAT || pgrp->epgrppbc == epgrppbcPREVSTEPCOM) { pgrp->x[m] += comm->pbcAtomBuffer[g][m]; - if (xp) + if (!xp.empty()) { pgrp->xp[m] += comm->pbcAtomBuffer[g][m]; } @@ -784,7 +798,7 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc pgrp->localWeights[i] = csw * std::cos(twopi_box * x[ii][pull->cosdim]) + snw * std::sin(twopi_box * x[ii][pull->cosdim]); } - if (xp) + if (!xp.empty()) { csw = comBuffer[2][0]; snw = comBuffer[2][1]; @@ -808,12 +822,12 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc using BoolVec = gmx::BasicVector; /* Returns whether the pull group obeys the PBC restrictions */ -static bool pullGroupObeysPbcRestrictions(const pull_group_work_t& group, - const BoolVec& dimUsed, - const rvec* x, - const t_pbc& pbc, - const gmx::RVec& x_pbc, - const real pbcMargin) +static bool pullGroupObeysPbcRestrictions(const pull_group_work_t& group, + const BoolVec& dimUsed, + gmx::ArrayRef x, + const t_pbc& pbc, + const gmx::RVec& x_pbc, + const real pbcMargin) { /* Determine which dimensions are relevant for PBC */ BoolVec dimUsesPbc = { false, false, false }; @@ -896,7 +910,7 @@ static bool pullGroupObeysPbcRestrictions(const pull_group_work_t& group, return true; } -int pullCheckPbcWithinGroups(const pull_t& pull, const rvec* x, const t_pbc& pbc, real pbcMargin) +int pullCheckPbcWithinGroups(const pull_t& pull, gmx::ArrayRef x, const t_pbc& pbc, real pbcMargin) { if (pbc.pbcType == PbcType::No) { @@ -976,7 +990,7 @@ bool pullCheckPbcWithinGroup(const pull_t& pull, } return (pullGroupObeysPbcRestrictions( - group, dimUsed, as_rvec_array(x.data()), pbc, pull.comm.pbcAtomBuffer[groupNr], pbcMargin)); + group, dimUsed, x, pbc, pull.comm.pbcAtomBuffer[groupNr], pbcMargin)); } void setPrevStepPullComFromState(struct pull_t* pull, const t_state* state) @@ -1019,7 +1033,11 @@ void allocStatePrevStepPullCom(t_state* state, const pull_t* pull) } } -void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc* pbc, const rvec x[]) +void initPullComFromPrevStep(const t_commrec* cr, + pull_t* pull, + gmx::ArrayRef masses, + t_pbc* pbc, + gmx::ArrayRef x) { pull_comm_t* comm = &pull->comm; size_t ngroup = pull->group.size(); @@ -1048,7 +1066,7 @@ void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* mass "Groups with no atoms, or only one atom, should not " "use the COM from the previous step as reference."); - rvec x_pbc = { 0, 0, 0 }; + gmx::RVec x_pbc = { 0, 0, 0 }; copy_rvec(comm->pbcAtomBuffer[g], x_pbc); if (debug) @@ -1068,7 +1086,7 @@ void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* mass if (pgrp->atomSet.numAtomsLocal() <= c_pullMaxNumLocalAtomsSingleThreaded) { - sum_com_part(pgrp, 0, pgrp->atomSet.numAtomsLocal(), x, nullptr, masses, pbc, x_pbc, &comSumsTotal); + sum_com_part(pgrp, 0, pgrp->atomSet.numAtomsLocal(), x, {}, masses, pbc, x_pbc, &comSumsTotal); } else { @@ -1077,8 +1095,7 @@ void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* mass { int ind_start = (pgrp->atomSet.numAtomsLocal() * (t + 0)) / pull->nthreads; int ind_end = (pgrp->atomSet.numAtomsLocal() * (t + 1)) / pull->nthreads; - sum_com_part( - pgrp, ind_start, ind_end, x, nullptr, masses, pbc, x_pbc, &pull->comSums[t]); + sum_com_part(pgrp, ind_start, ind_end, x, {}, masses, pbc, x_pbc, &pull->comSums[t]); } /* Reduce the thread contributions to sum_com[0] */ -- 2.22.0