From 78e7637ad0cf41cc62131b13103c09e8a0f48722 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Fri, 17 Jan 2020 10:59:18 +0100 Subject: [PATCH] Use ArrayRef(WithPadding) in constraint code This change is only refactoring, except for adding padding to two buffers used for flexible constraints. This avoid potential illegal memory access one element beyond the buffer with certain SIMD implementations. Note the this memory was not used nor modified. Change-Id: I385b3007a27888d15741e737ccfcf3e3a4369d1e --- src/gromacs/domdec/domdec.h | 6 +- src/gromacs/domdec/domdec_constraints.cpp | 9 +- src/gromacs/mdlib/constr.cpp | 191 +++++++++--------- src/gromacs/mdlib/constr.h | 29 +-- src/gromacs/mdlib/lincs.cpp | 129 ++++++------ src/gromacs/mdlib/lincs.h | 47 ++--- src/gromacs/mdlib/settle.cpp | 56 ++--- src/gromacs/mdlib/settle.h | 45 +++-- src/gromacs/mdlib/tests/constrtestrunners.cpp | 11 +- src/gromacs/mdlib/tests/settletestrunners.cpp | 9 +- src/gromacs/mdlib/update.cpp | 22 +- src/gromacs/mdrun/minimize.cpp | 24 ++- src/gromacs/mdrun/shellfc.cpp | 24 ++- .../modularsimulator/constraintelement.cpp | 21 +- 14 files changed, 330 insertions(+), 293 deletions(-) diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h index b7d4fb56dd..439e86b5e7 100644 --- a/src/gromacs/domdec/domdec.h +++ b/src/gromacs/domdec/domdec.h @@ -237,7 +237,11 @@ void dd_move_f_vsites(struct gmx_domdec_t* dd, rvec* f, rvec* fshift); void dd_clear_f_vsites(struct gmx_domdec_t* dd, rvec* f); /*! \brief Move x0 and also x1 if x1!=NULL. bX1IsCoord tells if to do PBC on x1 */ -void dd_move_x_constraints(struct gmx_domdec_t* dd, const matrix box, rvec* x0, rvec* x1, gmx_bool bX1IsCoord); +void dd_move_x_constraints(struct gmx_domdec_t* dd, + const matrix box, + gmx::ArrayRef x0, + gmx::ArrayRef x1, + gmx_bool bX1IsCoord); /*! \brief Communicates the coordinates involved in virtual sites */ void dd_move_x_vsites(struct gmx_domdec_t* dd, const matrix box, rvec* x); diff --git a/src/gromacs/domdec/domdec_constraints.cpp b/src/gromacs/domdec/domdec_constraints.cpp index 5bbe432cc2..ec39ab1494 100644 --- a/src/gromacs/domdec/domdec_constraints.cpp +++ b/src/gromacs/domdec/domdec_constraints.cpp @@ -103,11 +103,16 @@ struct gmx_domdec_constraints_t //! @endcond }; -void dd_move_x_constraints(gmx_domdec_t* dd, const matrix box, rvec* x0, rvec* x1, gmx_bool bX1IsCoord) +void dd_move_x_constraints(gmx_domdec_t* dd, + const matrix box, + gmx::ArrayRef x0, + gmx::ArrayRef x1, + gmx_bool bX1IsCoord) { if (dd->constraint_comm) { - dd_move_x_specat(dd, dd->constraint_comm, box, x0, x1, bX1IsCoord); + dd_move_x_specat(dd, dd->constraint_comm, box, as_rvec_array(x0.data()), + as_rvec_array(x1.data()), bX1IsCoord); ddReopenBalanceRegionCpu(dd); } diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index e17685d74c..14defa9ab0 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -82,7 +82,6 @@ #include "gromacs/utility/gmxassert.h" #include "gromacs/utility/listoflists.h" #include "gromacs/utility/pleasecite.h" -#include "gromacs/utility/smalloc.h" #include "gromacs/utility/txtdump.h" namespace gmx @@ -114,20 +113,20 @@ public: int numSettles); ~Impl(); void setConstraints(const gmx_localtop_t& top, const t_mdatoms& md); - bool apply(bool bLog, - bool bEner, - int64_t step, - int delta_step, - real step_scaling, - rvec* x, - rvec* xprime, - rvec* min_proj, - const matrix box, - real lambda, - real* dvdlambda, - rvec* v, - tensor* vir, - ConstraintVariable econq); + bool apply(bool bLog, + bool bEner, + int64_t step, + int delta_step, + real step_scaling, + ArrayRefWithPadding x, + ArrayRefWithPadding xprime, + ArrayRef min_proj, + const matrix box, + real lambda, + real* dvdlambda, + ArrayRefWithPadding v, + tensor* vir, + ConstraintVariable econq); //! The total number of constraints. int ncon_tot = 0; //! The number of flexible constraints. @@ -207,7 +206,7 @@ bool Constraints::havePerturbedConstraints() const } //! Clears constraint quantities for atoms in nonlocal region. -static void clear_constraint_quantity_nonlocal(gmx_domdec_t* dd, rvec* q) +static void clear_constraint_quantity_nonlocal(gmx_domdec_t* dd, ArrayRef q) { int nonlocal_at_start, nonlocal_at_end, at; @@ -232,14 +231,14 @@ void too_many_constraint_warnings(int eConstrAlg, int warncount) } //! Writes out coordinates. -static void write_constr_pdb(const char* fn, - const char* title, - const gmx_mtop_t& mtop, - int start, - int homenr, - const t_commrec* cr, - const rvec x[], - const matrix box) +static void write_constr_pdb(const char* fn, + const char* title, + const gmx_mtop_t& mtop, + int start, + int homenr, + const t_commrec* cr, + ArrayRef x, + const matrix box) { char fname[STRLEN]; FILE* out; @@ -294,15 +293,15 @@ static void write_constr_pdb(const char* fn, } //! Writes out domain contents to help diagnose crashes. -static void dump_confs(FILE* log, - int64_t step, - const gmx_mtop_t& mtop, - int start, - int homenr, - const t_commrec* cr, - const rvec x[], - rvec xprime[], - const matrix box) +static void dump_confs(FILE* log, + int64_t step, + const gmx_mtop_t& mtop, + int start, + int homenr, + const t_commrec* cr, + ArrayRef x, + ArrayRef xprime, + const matrix box) { char buf[STRLEN], buf2[22]; @@ -323,39 +322,39 @@ static void dump_confs(FILE* log, fprintf(stderr, "Wrote pdb files with previous and current coordinates\n"); } -bool Constraints::apply(bool bLog, - bool bEner, - int64_t step, - int delta_step, - real step_scaling, - rvec* x, - rvec* xprime, - rvec* min_proj, - const matrix box, - real lambda, - real* dvdlambda, - rvec* v, - tensor* vir, - ConstraintVariable econq) +bool Constraints::apply(bool bLog, + bool bEner, + int64_t step, + int delta_step, + real step_scaling, + ArrayRefWithPadding x, + ArrayRefWithPadding xprime, + ArrayRef min_proj, + const matrix box, + real lambda, + real* dvdlambda, + ArrayRefWithPadding v, + tensor* vir, + ConstraintVariable econq) { - return impl_->apply(bLog, bEner, step, delta_step, step_scaling, x, xprime, min_proj, box, - lambda, dvdlambda, v, vir, econq); + return impl_->apply(bLog, bEner, step, delta_step, step_scaling, std::move(x), std::move(xprime), + min_proj, box, lambda, dvdlambda, std::move(v), vir, econq); } -bool Constraints::Impl::apply(bool bLog, - bool bEner, - int64_t step, - int delta_step, - real step_scaling, - rvec* x, - rvec* xprime, - rvec* min_proj, - const matrix box, - real lambda, - real* dvdlambda, - rvec* v, - tensor* vir, - ConstraintVariable econq) +bool Constraints::Impl::apply(bool bLog, + bool bEner, + int64_t step, + int delta_step, + real step_scaling, + ArrayRefWithPadding x, + ArrayRefWithPadding xprime, + ArrayRef min_proj, + const matrix box, + real lambda, + real* dvdlambda, + ArrayRefWithPadding v, + tensor* vir, + ConstraintVariable econq) { bool bOK, bDump; int start, homenr; @@ -444,23 +443,24 @@ bool Constraints::Impl::apply(bool bLog, */ if (cr->dd) { - dd_move_x_constraints(cr->dd, box, x, xprime, econq == ConstraintVariable::Positions); + dd_move_x_constraints(cr->dd, box, x.unpaddedArrayRef(), xprime.unpaddedArrayRef(), + econq == ConstraintVariable::Positions); - if (v != nullptr) + if (!v.empty()) { /* We need to initialize the non-local components of v. * We never actually use these values, but we do increment them, * so we should avoid uninitialized variables and overflows. */ - clear_constraint_quantity_nonlocal(cr->dd, v); + clear_constraint_quantity_nonlocal(cr->dd, v.unpaddedArrayRef()); } } if (lincsd != nullptr) { bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md, cr, ms, x, xprime, min_proj, box, - pbc_null, lambda, dvdlambda, invdt, v, vir != nullptr, vir_r_m_dr, - econq, nrnb, maxwarn, &warncount_lincs); + pbc_null, lambda, dvdlambda, invdt, v.unpaddedArrayRef(), + vir != nullptr, vir_r_m_dr, econq, nrnb, maxwarn, &warncount_lincs); if (!bOK && maxwarn < INT_MAX) { if (log != nullptr) @@ -474,8 +474,11 @@ bool Constraints::Impl::apply(bool bLog, if (shaked != nullptr) { - bOK = constrain_shake(log, shaked, md.invmass, *idef, ir, x, xprime, min_proj, nrnb, lambda, - dvdlambda, invdt, v, vir != nullptr, vir_r_m_dr, maxwarn < INT_MAX, econq); + bOK = constrain_shake( + log, shaked, md.invmass, *idef, ir, as_rvec_array(x.unpaddedArrayRef().data()), + as_rvec_array(xprime.unpaddedArrayRef().data()), as_rvec_array(min_proj.data()), + nrnb, lambda, dvdlambda, invdt, as_rvec_array(v.unpaddedArrayRef().data()), + vir != nullptr, vir_r_m_dr, maxwarn < INT_MAX, econq); if (!bOK && maxwarn < INT_MAX) { @@ -505,14 +508,14 @@ bool Constraints::Impl::apply(bool bLog, clear_mat(vir_r_m_dr_th[th]); } - csettle(settled, nth, th, pbc_null, x[0], xprime[0], invdt, v ? v[0] : nullptr, - vir != nullptr, th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th], + csettle(settled, nth, th, pbc_null, x, xprime, invdt, v, vir != nullptr, + th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th], th == 0 ? &bSettleErrorHasOccurred0 : &bSettleErrorHasOccurred[th]); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } inc_nrnb(nrnb, eNR_SETTLE, nsettle); - if (v != nullptr) + if (!v.empty()) { inc_nrnb(nrnb, eNR_CONSTR_V, nsettle * 3); } @@ -553,8 +556,8 @@ bool Constraints::Impl::apply(bool bLog, { settle_proj(settled, econq, end_th - start_th, settle->iatoms + start_th * (1 + NRAL(F_SETTLE)), pbc_null, - x, xprime, min_proj, calcvir_atom_end, - th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]); + x.unpaddedArrayRef(), xprime.unpaddedArrayRef(), min_proj, + calcvir_atom_end, th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]); } } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR @@ -644,7 +647,8 @@ bool Constraints::Impl::apply(bool bLog, if (bDump) { - dump_confs(log, step, mtop, start, homenr, cr, x, xprime, box); + dump_confs(log, step, mtop, start, homenr, cr, x.unpaddedArrayRef(), + xprime.unpaddedArrayRef(), box); } if (econq == ConstraintVariable::Positions) @@ -660,19 +664,24 @@ bool Constraints::Impl::apply(bool bLog, t = ir.init_t; } set_pbc(&pbc, ir.pbcType, box); - pull_constraint(pull_work, &md, &pbc, cr, ir.delta_t, t, x, xprime, v, *vir); + pull_constraint(pull_work, &md, &pbc, cr, ir.delta_t, t, + as_rvec_array(x.unpaddedArrayRef().data()), + as_rvec_array(xprime.unpaddedArrayRef().data()), + as_rvec_array(v.unpaddedArrayRef().data()), *vir); } if (ed && delta_step > 0) { /* apply the essential dynamics constraints here */ - do_edsam(&ir, step, cr, xprime, v, box, ed); + do_edsam(&ir, step, cr, as_rvec_array(xprime.unpaddedArrayRef().data()), + as_rvec_array(v.unpaddedArrayRef().data()), box, ed); } } wallcycle_stop(wcycle, ewcCONSTR); - if (v != nullptr && md.cFREEZE) + if (!v.empty() && md.cFREEZE) { /* Set the velocities of frozen dimensions to zero */ + ArrayRef vRef = v.unpaddedArrayRef(); int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate); @@ -685,7 +694,7 @@ bool Constraints::Impl::apply(bool bLog, { if (ir.opts.nFreeze[freezeGroup][d]) { - v[i][d] = 0; + vRef[i][d] = 0; } } } @@ -1136,15 +1145,8 @@ void do_constrain_first(FILE* fplog, int64_t step; real dt = ir->delta_t; real dvdl_dum; - rvec* savex; - auto xRvec = as_rvec_array(x.paddedArrayRef().data()); - auto vRvec = as_rvec_array(v.paddedArrayRef().data()); - - /* We need to allocate one element extra, since we might use - * (unaligned) 4-wide SIMD loads to access rvec entries. - */ - snew(savex, natoms + 1); + PaddedVector savex(natoms); start = 0; end = md->homenr; @@ -1163,14 +1165,14 @@ void do_constrain_first(FILE* fplog, dvdl_dum = 0; /* constrain the current position */ - constr->apply(TRUE, FALSE, step, 0, 1.0, xRvec, xRvec, nullptr, box, lambda, &dvdl_dum, nullptr, - nullptr, gmx::ConstraintVariable::Positions); + constr->apply(TRUE, FALSE, step, 0, 1.0, x, x, ArrayRef(), box, lambda, &dvdl_dum, + ArrayRefWithPadding(), nullptr, gmx::ConstraintVariable::Positions); if (EI_VV(ir->eI)) { /* constrain the inital velocity, and save it */ /* also may be useful if we need the ekin from the halfstep for velocity verlet */ - constr->apply(TRUE, FALSE, step, 0, 1.0, xRvec, vRvec, vRvec, box, lambda, &dvdl_dum, - nullptr, nullptr, gmx::ConstraintVariable::Velocities); + constr->apply(TRUE, FALSE, step, 0, 1.0, x, v, v.unpaddedArrayRef(), box, lambda, &dvdl_dum, + ArrayRefWithPadding(), nullptr, gmx::ConstraintVariable::Velocities); } /* constrain the inital velocities at t-dt/2 */ if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV) @@ -1196,8 +1198,8 @@ void do_constrain_first(FILE* fplog, fprintf(fplog, "\nConstraining the coordinates at t0-dt (step %s)\n", gmx_step_str(step, buf)); } dvdl_dum = 0; - constr->apply(TRUE, FALSE, step, -1, 1.0, xRvec, savex, nullptr, box, lambda, &dvdl_dum, - vRvec, nullptr, gmx::ConstraintVariable::Positions); + constr->apply(TRUE, FALSE, step, -1, 1.0, x, savex.arrayRefWithPadding(), ArrayRef(), + box, lambda, &dvdl_dum, v, nullptr, gmx::ConstraintVariable::Positions); for (i = start; i < end; i++) { @@ -1208,7 +1210,6 @@ void do_constrain_first(FILE* fplog, } } } - sfree(savex); } } // namespace gmx diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h index 9824e43a97..c06d7d3576 100644 --- a/src/gromacs/mdlib/constr.h +++ b/src/gromacs/mdlib/constr.h @@ -49,6 +49,7 @@ #include +#include "gromacs/math/arrayrefwithpadding.h" #include "gromacs/math/vectypes.h" #include "gromacs/topology/idef.h" #include "gromacs/utility/arrayref.h" @@ -170,20 +171,20 @@ public: * * /note x is non-const, because non-local atoms need to be communicated. */ - bool apply(bool bLog, - bool bEner, - int64_t step, - int delta_step, - real step_scaling, - rvec* x, - rvec* xprime, - rvec* min_proj, - const matrix box, - real lambda, - real* dvdlambda, - rvec* v, - tensor* vir, - ConstraintVariable econq); + bool apply(bool bLog, + bool bEner, + int64_t step, + int delta_step, + real step_scaling, + ArrayRefWithPadding x, + ArrayRefWithPadding xprime, + ArrayRef min_proj, + const matrix box, + real lambda, + real* dvdlambda, + ArrayRefWithPadding v, + tensor* vir, + ConstraintVariable econq); //! Links the essentialdynamics and constraint code. void saveEdsamPointer(gmx_edsam* ed); //! Getter for use by domain decomposition. diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp index 98fdc2e03c..b425d12eaf 100644 --- a/src/gromacs/mdlib/lincs.cpp +++ b/src/gromacs/mdlib/lincs.cpp @@ -569,17 +569,17 @@ static void gmx_simdcall calc_dr_x_f_simd(int b0, #endif // GMX_SIMD_HAVE_REAL /*! \brief LINCS projection, works on derivatives of the coordinates. */ -static void do_lincsp(const rvec* x, - rvec* f, - rvec* fp, - t_pbc* pbc, - Lincs* lincsd, - int th, - real* invmass, - ConstraintVariable econq, - bool bCalcDHDL, - bool bCalcVir, - tensor rmdf) +static void do_lincsp(ArrayRefWithPadding xPadded, + ArrayRefWithPadding fPadded, + ArrayRef fp, + t_pbc* pbc, + Lincs* lincsd, + int th, + real* invmass, + ConstraintVariable econq, + bool bCalcDHDL, + bool bCalcVir, + tensor rmdf) { const int b0 = lincsd->task[th].b0; const int b1 = lincsd->task[th].b1; @@ -608,6 +608,9 @@ static void do_lincsp(const rvec* x, gmx::ArrayRef rhs2 = lincsd->tmp2; gmx::ArrayRef sol = lincsd->tmp3; + const rvec* x = as_rvec_array(xPadded.paddedArrayRef().data()); + rvec* f = as_rvec_array(fPadded.paddedArrayRef().data()); + #if GMX_SIMD_HAVE_REAL /* This SIMD code does the same as the plain-C code after the #else. * The only difference is that we always call pbc code, as with SIMD @@ -706,8 +709,8 @@ static void do_lincsp(const rvec* x, /* When constraining forces, we should not use mass weighting, * so we pass invmass=NULL, which results in the use of 1 for all atoms. */ - lincs_update_atoms(lincsd, th, 1.0, sol, r, - (econq != ConstraintVariable::Force) ? invmass : nullptr, fp); + lincs_update_atoms(lincsd, th, 1.0, sol, r, (econq != ConstraintVariable::Force) ? invmass : nullptr, + as_rvec_array(fp.data())); if (bCalcDHDL) { @@ -945,22 +948,26 @@ static void gmx_simdcall calc_dist_iter_simd(int b0, #endif // GMX_SIMD_HAVE_REAL //! Implements LINCS constraining. -static void do_lincs(const rvec* x, - rvec* xp, - const matrix box, - t_pbc* pbc, - Lincs* lincsd, - int th, - const real* invmass, - const t_commrec* cr, - bool bCalcDHDL, - real wangle, - bool* bWarn, - real invdt, - rvec* gmx_restrict v, - bool bCalcVir, - tensor vir_r_m_dr) +static void do_lincs(ArrayRefWithPadding xPadded, + ArrayRefWithPadding xpPadded, + const matrix box, + t_pbc* pbc, + Lincs* lincsd, + int th, + const real* invmass, + const t_commrec* cr, + bool bCalcDHDL, + real wangle, + bool* bWarn, + real invdt, + ArrayRef vRef, + bool bCalcVir, + tensor vir_r_m_dr) { + const rvec* x = as_rvec_array(xPadded.paddedArrayRef().data()); + rvec* xp = as_rvec_array(xpPadded.paddedArrayRef().data()); + rvec* gmx_restrict v = as_rvec_array(vRef.data()); + const int b0 = lincsd->task[th].b0; const int b1 = lincsd->task[th].b1; @@ -1098,7 +1105,7 @@ static void do_lincs(const rvec* x, /* Communicate the corrected non-local coordinates */ if (DOMAINDECOMP(cr)) { - dd_move_x_constraints(cr->dd, box, xp, nullptr, FALSE); + dd_move_x_constraints(cr->dd, box, xpPadded.unpaddedArrayRef(), ArrayRef(), FALSE); } } #pragma omp barrier @@ -2124,8 +2131,8 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ //! Issues a warning when LINCS constraints cannot be satisfied. static void lincs_warning(gmx_domdec_t* dd, - const rvec* x, - rvec* xprime, + ArrayRef x, + ArrayRef xprime, t_pbc* pbc, int ncons, gmx::ArrayRef atoms, @@ -2192,7 +2199,7 @@ struct LincsDeviations }; //! Determine how well the constraints have been satisfied. -static LincsDeviations makeLincsDeviations(const Lincs& lincsd, const rvec* x, const t_pbc* pbc) +static LincsDeviations makeLincsDeviations(const Lincs& lincsd, ArrayRef x, const t_pbc* pbc) { LincsDeviations result; const ArrayRef atoms = lincsd.atoms; @@ -2241,28 +2248,28 @@ static LincsDeviations makeLincsDeviations(const Lincs& lincsd, const rvec* x, c return result; } -bool constrain_lincs(bool computeRmsd, - const t_inputrec& ir, - int64_t step, - Lincs* lincsd, - const t_mdatoms& md, - const t_commrec* cr, - const gmx_multisim_t* ms, - const rvec* x, - rvec* xprime, - rvec* min_proj, - const matrix box, - t_pbc* pbc, - real lambda, - real* dvdlambda, - real invdt, - rvec* v, - bool bCalcVir, - tensor vir_r_m_dr, - ConstraintVariable econq, - t_nrnb* nrnb, - int maxwarn, - int* warncount) +bool constrain_lincs(bool computeRmsd, + const t_inputrec& ir, + int64_t step, + Lincs* lincsd, + const t_mdatoms& md, + const t_commrec* cr, + const gmx_multisim_t* ms, + ArrayRefWithPadding xPadded, + ArrayRefWithPadding xprimePadded, + ArrayRef min_proj, + const matrix box, + t_pbc* pbc, + real lambda, + real* dvdlambda, + real invdt, + ArrayRef v, + bool bCalcVir, + tensor vir_r_m_dr, + ConstraintVariable econq, + t_nrnb* nrnb, + int maxwarn, + int* warncount) { bool bOK = TRUE; @@ -2282,6 +2289,9 @@ bool constrain_lincs(bool computeRmsd, return bOK; } + ArrayRef x = xPadded.unpaddedArrayRef(); + ArrayRef xprime = xprimePadded.unpaddedArrayRef(); + if (econq == ConstraintVariable::Positions) { /* We can't use bCalcDHDL here, since NULL can be passed for dvdlambda @@ -2355,8 +2365,9 @@ bool constrain_lincs(bool computeRmsd, clear_mat(lincsd->task[th].vir_r_m_dr); - do_lincs(x, xprime, box, pbc, lincsd, th, md.invmass, cr, bCalcDHDL, ir.LincsWarnAngle, - &bWarn, invdt, v, bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr); + do_lincs(xPadded, xprimePadded, box, pbc, lincsd, th, md.invmass, cr, bCalcDHDL, + ir.LincsWarnAngle, &bWarn, invdt, v, bCalcVir, + th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } @@ -2433,8 +2444,8 @@ bool constrain_lincs(bool computeRmsd, { int th = gmx_omp_get_thread_num(); - do_lincsp(x, xprime, min_proj, pbc, lincsd, th, md.invmass, econq, bCalcDHDL, - bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr); + do_lincsp(xPadded, xprimePadded, min_proj, pbc, lincsd, th, md.invmass, econq, + bCalcDHDL, bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } @@ -2475,7 +2486,7 @@ bool constrain_lincs(bool computeRmsd, { inc_nrnb(nrnb, eNR_LINCSMAT, lincsd->nOrder * lincsd->ncc_triangle); } - if (v) + if (!v.empty()) { inc_nrnb(nrnb, eNR_CONSTR_V, lincsd->nc_real * 2); } diff --git a/src/gromacs/mdlib/lincs.h b/src/gromacs/mdlib/lincs.h index ef43c3a586..649c58d661 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, by the GROMACS development team, led by + * 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. @@ -46,6 +46,7 @@ #include +#include "gromacs/math/arrayrefwithpadding.h" #include "gromacs/math/vectypes.h" #include "gromacs/utility/arrayref.h" #include "gromacs/utility/basedefinitions.h" @@ -93,28 +94,28 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ /*! \brief Applies LINCS constraints. * * \returns true if the constraining succeeded. */ -bool constrain_lincs(bool computeRmsd, - const t_inputrec& ir, - int64_t step, - Lincs* lincsd, - const t_mdatoms& md, - const t_commrec* cr, - const gmx_multisim_t* ms, - const rvec* x, - rvec* xprime, - rvec* min_proj, - const matrix box, - t_pbc* pbc, - real lambda, - real* dvdlambda, - real invdt, - rvec* v, - bool bCalcVir, - tensor vir_r_m_dr, - ConstraintVariable econq, - t_nrnb* nrnb, - int maxwarn, - int* warncount); +bool constrain_lincs(bool computeRmsd, + const t_inputrec& ir, + int64_t step, + Lincs* lincsd, + const t_mdatoms& md, + const t_commrec* cr, + const gmx_multisim_t* ms, + ArrayRefWithPadding x, + ArrayRefWithPadding xprime, + ArrayRef min_proj, + const matrix box, + t_pbc* pbc, + real lambda, + real* dvdlambda, + real invdt, + ArrayRef v, + bool bCalcVir, + tensor vir_r_m_dr, + ConstraintVariable econq, + t_nrnb* nrnb, + int maxwarn, + int* warncount); } // namespace gmx diff --git a/src/gromacs/mdlib/settle.cpp b/src/gromacs/mdlib/settle.cpp index f8f183bd46..2fa116a051 100644 --- a/src/gromacs/mdlib/settle.cpp +++ b/src/gromacs/mdlib/settle.cpp @@ -306,16 +306,16 @@ void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const } } -void settle_proj(settledata* settled, - ConstraintVariable econq, - int nsettle, - const t_iatom iatoms[], - const t_pbc* pbc, - const rvec x[], - rvec* der, - rvec* derp, - int calcvir_atom_end, - tensor vir_r_m_dder) +void settle_proj(settledata* settled, + ConstraintVariable econq, + int nsettle, + const t_iatom iatoms[], + const t_pbc* pbc, + ArrayRef x, + ArrayRef der, + ArrayRef derp, + int calcvir_atom_end, + tensor vir_r_m_dder) { /* Settle for projection out constraint components * of derivatives of the coordinates. @@ -812,18 +812,22 @@ static void settleTemplateWrapper(settledata* settled, } } -void csettle(settledata* settled, - int nthread, - int thread, - const t_pbc* pbc, - const real x[], - real xprime[], - real invdt, - real* v, - bool bCalcVirial, - tensor vir_r_m_dr, - bool* bErrorHasOccurred) +void csettle(settledata* settled, + int nthread, + int thread, + const t_pbc* pbc, + ArrayRefWithPadding x, + ArrayRefWithPadding xprime, + real invdt, + ArrayRefWithPadding v, + bool bCalcVirial, + tensor vir_r_m_dr, + bool* bErrorHasOccurred) { + const real* xPtr = as_rvec_array(x.paddedArrayRef().data())[0]; + real* xprimePtr = as_rvec_array(xprime.paddedArrayRef().data())[0]; + real* vPtr = as_rvec_array(v.paddedArrayRef().data())[0]; + #if GMX_SIMD_HAVE_REAL if (settled->bUseSimd) { @@ -832,8 +836,8 @@ void csettle(settledata* settled, set_pbc_simd(pbc, pbcSimd); settleTemplateWrapper( - settled, nthread, thread, pbcSimd, x, xprime, invdt, v, bCalcVirial, vir_r_m_dr, - bErrorHasOccurred); + settled, nthread, thread, pbcSimd, xPtr, xprimePtr, invdt, vPtr, bCalcVirial, + vir_r_m_dr, bErrorHasOccurred); } else #endif @@ -852,9 +856,9 @@ void csettle(settledata* settled, pbcNonNull = &pbcNo; } - settleTemplateWrapper(settled, nthread, thread, pbcNonNull, x, - xprime, invdt, v, bCalcVirial, - vir_r_m_dr, bErrorHasOccurred); + settleTemplateWrapper(settled, nthread, thread, pbcNonNull, + &xPtr[0], &xprimePtr[0], invdt, &vPtr[0], + bCalcVirial, vir_r_m_dr, bErrorHasOccurred); } } diff --git a/src/gromacs/mdlib/settle.h b/src/gromacs/mdlib/settle.h index f503ac8481..7cc0a89bea 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, by the GROMACS development team, led by + * 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. @@ -44,6 +44,7 @@ #ifndef GMX_MDLIB_SETTLE_H #define GMX_MDLIB_SETTLE_H +#include "gromacs/math/arrayrefwithpadding.h" #include "gromacs/topology/idef.h" struct gmx_cmap_t; @@ -72,32 +73,32 @@ void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const /*! \brief Constrain coordinates using SETTLE. * Can be called on any number of threads. */ -void csettle(settledata* settled, /* The SETTLE structure */ - int nthread, /* The number of threads used */ - int thread, /* Our thread index */ - const t_pbc* pbc, /* PBC data pointer, can be NULL */ - const real x[], /* Reference coordinates */ - real xprime[], /* New coords, to be settled */ - real invdt, /* 1/delta_t */ - real* v, /* Also constrain v if v!=NULL */ - bool bCalcVirial, /* Calculate the virial contribution */ - tensor vir_r_m_dr, /* sum r x m delta_r */ - bool* bErrorHasOccurred /* True if a settle error occurred */ +void csettle(settledata* settled, /* The SETTLE structure */ + int nthread, /* The number of threads used */ + int thread, /* Our thread index */ + const t_pbc* pbc, /* PBC data pointer, can be NULL */ + ArrayRefWithPadding x, /* Reference coordinates */ + ArrayRefWithPadding xprime, /* New coords, to be settled */ + real invdt, /* 1/delta_t */ + ArrayRefWithPadding v, /* Also constrain v if v!=NULL */ + bool bCalcVirial, /* Calculate the virial contribution */ + tensor vir_r_m_dr, /* sum r x m delta_r */ + bool* bErrorHasOccurred /* True if a settle error occurred */ ); /*! \brief Analytical algorithm to subtract the components of derivatives * of coordinates working on settle type constraint. */ -void settle_proj(settledata* settled, - ConstraintVariable econq, - int nsettle, - const t_iatom iatoms[], - const t_pbc* pbc, /* PBC data pointer, can be NULL */ - const rvec x[], - rvec* der, - rvec* derp, - int CalcVirAtomEnd, - tensor vir_r_m_dder); +void settle_proj(settledata* settled, + ConstraintVariable econq, + int nsettle, + const t_iatom iatoms[], + const t_pbc* pbc, /* PBC data pointer, can be NULL */ + ArrayRef x, + ArrayRef der, + ArrayRef derp, + int CalcVirAtomEnd, + tensor vir_r_m_dder); } // namespace gmx diff --git a/src/gromacs/mdlib/tests/constrtestrunners.cpp b/src/gromacs/mdlib/tests/constrtestrunners.cpp index b1f3356092..5f974aaf8b 100644 --- a/src/gromacs/mdlib/tests/constrtestrunners.cpp +++ b/src/gromacs/mdlib/tests/constrtestrunners.cpp @@ -131,11 +131,12 @@ void applyLincs(ConstraintsTestData* testData, t_pbc pbc) // Evaluate constraints bool success = constrain_lincs( false, testData->ir_, 0, lincsd, testData->md_, &testData->cr_, &testData->ms_, - as_rvec_array(testData->x_.data()), as_rvec_array(testData->xPrime_.data()), - as_rvec_array(testData->xPrime2_.data()), pbc.box, &pbc, testData->md_.lambda, - &testData->dHdLambda_, testData->invdt_, as_rvec_array(testData->v_.data()), - testData->computeVirial_, testData->virialScaled_, gmx::ConstraintVariable::Positions, - &testData->nrnb_, maxwarn, &warncount_lincs); + testData->x_.arrayRefWithPadding(), testData->xPrime_.arrayRefWithPadding(), + testData->xPrime2_.arrayRefWithPadding().unpaddedArrayRef(), pbc.box, &pbc, + testData->md_.lambda, &testData->dHdLambda_, testData->invdt_, + testData->v_.arrayRefWithPadding().unpaddedArrayRef(), testData->computeVirial_, + testData->virialScaled_, gmx::ConstraintVariable::Positions, &testData->nrnb_, maxwarn, + &warncount_lincs); EXPECT_TRUE(success) << "Test failed with a false return value in LINCS."; EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS."; done_lincs(lincsd); diff --git a/src/gromacs/mdlib/tests/settletestrunners.cpp b/src/gromacs/mdlib/tests/settletestrunners.cpp index 88a4c1774e..59cc352b87 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, by the GROMACS development team, led by + * 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. @@ -70,10 +70,9 @@ void applySettle(SettleTestData* testData, bool errorOccured; int numThreads = 1; int threadIndex = 0; - csettle(settled, numThreads, threadIndex, &pbc, - reinterpret_cast(as_rvec_array(testData->x_.data())), - reinterpret_cast(as_rvec_array(testData->xPrime_.data())), testData->reciprocalTimeStep_, - updateVelocities ? reinterpret_cast(as_rvec_array(testData->v_.data())) : nullptr, + csettle(settled, numThreads, threadIndex, &pbc, testData->x_.arrayRefWithPadding(), + testData->xPrime_.arrayRefWithPadding(), testData->reciprocalTimeStep_, + updateVelocities ? testData->v_.arrayRefWithPadding() : ArrayRefWithPadding(), calcVirial, testData->virial_, &errorOccured); settle_free(settled); EXPECT_FALSE(errorOccured) << testDescription; diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp index 1afa13fecf..275bdf68c1 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -1465,9 +1465,10 @@ void constrain_velocities(int64_t step, clear_mat(vir_part); /* Constrain the coordinates upd->xp */ - constr->apply(do_log, do_ene, step, 1, 1.0, state->x.rvec_array(), state->v.rvec_array(), - state->v.rvec_array(), state->box, state->lambda[efptBONDED], dvdlambda, - nullptr, bCalcVir ? &vir_con : nullptr, ConstraintVariable::Velocities); + constr->apply(do_log, do_ene, step, 1, 1.0, state->x.arrayRefWithPadding(), + state->v.arrayRefWithPadding(), state->v.arrayRefWithPadding().unpaddedArrayRef(), + state->box, state->lambda[efptBONDED], dvdlambda, ArrayRefWithPadding(), + bCalcVir ? &vir_con : nullptr, ConstraintVariable::Velocities); if (bCalcVir) { @@ -1498,10 +1499,10 @@ void constrain_coordinates(int64_t step, clear_mat(vir_part); /* Constrain the coordinates upd->xp */ - constr->apply(do_log, do_ene, step, 1, 1.0, state->x.rvec_array(), upd->xp()->rvec_array(), - nullptr, state->box, state->lambda[efptBONDED], dvdlambda, - as_rvec_array(state->v.data()), bCalcVir ? &vir_con : nullptr, - ConstraintVariable::Positions); + constr->apply(do_log, do_ene, step, 1, 1.0, state->x.arrayRefWithPadding(), + upd->xp()->arrayRefWithPadding(), ArrayRef(), state->box, + state->lambda[efptBONDED], dvdlambda, state->v.arrayRefWithPadding(), + bCalcVir ? &vir_con : nullptr, ConstraintVariable::Positions); if (bCalcVir) { @@ -1564,9 +1565,10 @@ void update_sd_second_half(int64_t step, wallcycle_stop(wcycle, ewcUPDATE); /* Constrain the coordinates upd->xp for half a time step */ - constr->apply(do_log, do_ene, step, 1, 0.5, state->x.rvec_array(), upd->xp()->rvec_array(), - nullptr, state->box, state->lambda[efptBONDED], dvdlambda, - as_rvec_array(state->v.data()), nullptr, ConstraintVariable::Positions); + constr->apply(do_log, do_ene, step, 1, 0.5, state->x.arrayRefWithPadding(), + upd->xp()->arrayRefWithPadding(), ArrayRef(), state->box, + state->lambda[efptBONDED], dvdlambda, state->v.arrayRefWithPadding(), nullptr, + ConstraintVariable::Positions); } } diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 593d432b37..64da82baac 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -108,7 +108,9 @@ #include "legacysimulator.h" #include "shellfc.h" +using gmx::ArrayRef; using gmx::MdrunScheduleWorkload; +using gmx::RVec; //! Utility structure for manipulating states during EM typedef struct @@ -455,8 +457,9 @@ static void init_em(FILE* fplog, { /* Constrain the starting coordinates */ dvdl_constr = 0; - constr->apply(TRUE, TRUE, -1, 0, 1.0, ems->s.x.rvec_array(), ems->s.x.rvec_array(), - nullptr, ems->s.box, ems->s.lambda[efptFEP], &dvdl_constr, nullptr, + constr->apply(TRUE, TRUE, -1, 0, 1.0, ems->s.x.arrayRefWithPadding(), + ems->s.x.arrayRefWithPadding(), ArrayRef(), ems->s.box, + ems->s.lambda[efptFEP], &dvdl_constr, gmx::ArrayRefWithPadding(), nullptr, gmx::ConstraintVariable::Positions); } } @@ -680,9 +683,10 @@ static bool do_em_step(const t_commrec* cr, if (constr) { dvdl_constr = 0; - validStep = constr->apply(TRUE, TRUE, count, 0, 1.0, s1->x.rvec_array(), s2->x.rvec_array(), - nullptr, s2->box, s2->lambda[efptBONDED], &dvdl_constr, nullptr, - nullptr, gmx::ConstraintVariable::Positions); + validStep = constr->apply( + TRUE, TRUE, count, 0, 1.0, s1->x.arrayRefWithPadding(), s2->x.arrayRefWithPadding(), + ArrayRef(), s2->box, s2->lambda[efptBONDED], &dvdl_constr, + gmx::ArrayRefWithPadding(), nullptr, gmx::ConstraintVariable::Positions); if (cr->nnodes > 1) { @@ -887,11 +891,11 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, if (constr) { /* Project out the constraint components of the force */ - dvdl_constr = 0; - rvec* f_rvec = ems->f.rvec_array(); - constr->apply(FALSE, FALSE, count, 0, 1.0, ems->s.x.rvec_array(), f_rvec, f_rvec, - ems->s.box, ems->s.lambda[efptBONDED], &dvdl_constr, nullptr, &shake_vir, - gmx::ConstraintVariable::ForceDispl); + dvdl_constr = 0; + auto f = ems->f.arrayRefWithPadding(); + constr->apply(FALSE, FALSE, count, 0, 1.0, ems->s.x.arrayRefWithPadding(), f, + f.unpaddedArrayRef(), ems->s.box, ems->s.lambda[efptBONDED], &dvdl_constr, + gmx::ArrayRefWithPadding(), &shake_vir, gmx::ConstraintVariable::ForceDispl); enerd->term[F_DVDL_CONSTR] += dvdl_constr; m_add(force_vir, shake_vir, vir); } diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index a864ef0d48..094674d983 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -115,8 +115,8 @@ struct gmx_shellfc_t /* Flexible constraint working data */ std::vector acc_dir; /* Acceleration direction for flexcon */ gmx::PaddedVector x_old; /* Old coordinates for flexcon */ - std::vector adir_xnold; /* Work space for init_adir */ - std::vector adir_xnew; /* Work space for init_adir */ + gmx::PaddedVector adir_xnold; /* Work space for init_adir */ + gmx::PaddedVector adir_xnew; /* Work space for init_adir */ std::int64_t numForceEvaluations; /* Total number of force evaluations */ int numConvergedIterations; /* Total number of iterations that converged */ }; @@ -872,8 +872,8 @@ static void init_adir(gmx_shellfc_t* shfc, { n = end; } - shfc->adir_xnold.resize(n); - shfc->adir_xnew.resize(n); + shfc->adir_xnold.resizeWithPadding(n); + shfc->adir_xnew.resizeWithPadding(n); rvec* xnold = as_rvec_array(shfc->adir_xnold.data()); rvec* xnew = as_rvec_array(shfc->adir_xnew.data()); rvec* x_old = as_rvec_array(xOld.paddedArrayRef().data()); @@ -902,10 +902,12 @@ static void init_adir(gmx_shellfc_t* shfc, } } } - constr->apply(FALSE, FALSE, step, 0, 1.0, x, xnold, nullptr, box, lambda[efptBONDED], - &(dvdlambda[efptBONDED]), nullptr, nullptr, gmx::ConstraintVariable::Positions); - constr->apply(FALSE, FALSE, step, 0, 1.0, x, xnew, nullptr, box, lambda[efptBONDED], - &(dvdlambda[efptBONDED]), nullptr, nullptr, gmx::ConstraintVariable::Positions); + constr->apply(FALSE, FALSE, step, 0, 1.0, xCurrent, shfc->adir_xnold.arrayRefWithPadding(), + ArrayRef(), box, lambda[efptBONDED], &(dvdlambda[efptBONDED]), + ArrayRefWithPadding(), nullptr, gmx::ConstraintVariable::Positions); + constr->apply(FALSE, FALSE, step, 0, 1.0, xCurrent, shfc->adir_xnew.arrayRefWithPadding(), + ArrayRef(), box, lambda[efptBONDED], &(dvdlambda[efptBONDED]), + ArrayRefWithPadding(), nullptr, gmx::ConstraintVariable::Positions); for (n = 0; n < end; n++) { @@ -918,9 +920,9 @@ static void init_adir(gmx_shellfc_t* shfc, } /* Project the acceleration on the old bond directions */ - constr->apply(FALSE, FALSE, step, 0, 1.0, x_old, xnew, as_rvec_array(acc_dir.data()), box, - lambda[efptBONDED], &(dvdlambda[efptBONDED]), nullptr, nullptr, - gmx::ConstraintVariable::Deriv_FlexCon); + constr->apply(FALSE, FALSE, step, 0, 1.0, xOld, shfc->adir_xnew.arrayRefWithPadding(), acc_dir, + box, lambda[efptBONDED], &(dvdlambda[efptBONDED]), ArrayRefWithPadding(), + nullptr, gmx::ConstraintVariable::Deriv_FlexCon); } void relax_shell_flexcon(FILE* fplog, diff --git a/src/gromacs/modularsimulator/constraintelement.cpp b/src/gromacs/modularsimulator/constraintelement.cpp index 90f993f005..4a0b7fe545 100644 --- a/src/gromacs/modularsimulator/constraintelement.cpp +++ b/src/gromacs/modularsimulator/constraintelement.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2019, by the GROMACS development team, led by + * Copyright (c) 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. @@ -125,7 +125,10 @@ void ConstraintsElement::apply(Step step, bool calculateVirial, bool w { tensor vir_con; - rvec *x, *xprime, *min_proj, *v; + ArrayRefWithPadding x; + ArrayRefWithPadding xprime; + ArrayRef min_proj; + ArrayRefWithPadding v; // disabled functionality real lambda = 0; @@ -134,16 +137,14 @@ void ConstraintsElement::apply(Step step, bool calculateVirial, bool w switch (variable) { case ConstraintVariable::Positions: - x = as_rvec_array(statePropagatorData_->previousPositionsView().paddedArrayRef().data()); - xprime = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data()); - min_proj = nullptr; - v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data()); + x = statePropagatorData_->previousPositionsView(); + xprime = statePropagatorData_->positionsView(); + v = statePropagatorData_->velocitiesView(); break; case ConstraintVariable::Velocities: - x = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data()); - xprime = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data()); - min_proj = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data()); - v = nullptr; + x = statePropagatorData_->positionsView(); + xprime = statePropagatorData_->velocitiesView(); + min_proj = statePropagatorData_->velocitiesView().unpaddedArrayRef(); break; default: gmx_fatal(FARGS, "Constraint algorithm not implemented for modular simulator."); } -- 2.22.0