From 2eda57c96ee1189401f8f2484952e29bc149f4fc Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Thu, 25 Mar 2021 10:48:38 +0000 Subject: [PATCH] Use ArrayRef in special forces Pass ArrayRef instead of rvec* in edsam, pull_rotation, and imd. Also removed unneeded header basedefinitions. A followup change can refactor communicate_group_positions to take ArrayRefs. --- src/gromacs/essentialdynamics/edsam.cpp | 58 ++++++++-------- src/gromacs/essentialdynamics/edsam.h | 40 +++++------ src/gromacs/imd/imd.cpp | 48 +++++++------- src/gromacs/imd/imd.h | 65 +++++++++--------- src/gromacs/mdlib/constr.cpp | 8 +-- src/gromacs/mdlib/sim_util.cpp | 12 ++-- src/gromacs/mdrun/md.cpp | 2 +- src/gromacs/mdrun/minimize.cpp | 6 +- src/gromacs/mdrun/runner.cpp | 2 +- src/gromacs/pulling/pull_rotation.cpp | 88 +++++++++++++------------ src/gromacs/pulling/pull_rotation.h | 17 +++-- 11 files changed, 180 insertions(+), 166 deletions(-) diff --git a/src/gromacs/essentialdynamics/edsam.cpp b/src/gromacs/essentialdynamics/edsam.cpp index fa612d238d..7f33f49cbb 100644 --- a/src/gromacs/essentialdynamics/edsam.cpp +++ b/src/gromacs/essentialdynamics/edsam.cpp @@ -904,14 +904,14 @@ static void update_adaption(t_edpar* edi) } -static void do_single_flood(FILE* edo, - const rvec x[], - rvec force[], - t_edpar* edi, - int64_t step, - const matrix box, - const t_commrec* cr, - gmx_bool bNS) /* Are we in a neighbor searching step? */ +static void do_single_flood(FILE* edo, + gmx::ArrayRef coords, + gmx::ArrayRef force, + t_edpar* edi, + int64_t step, + const matrix box, + const t_commrec* cr, + gmx_bool bNS) /* Are we in a neighbor searching step? */ { int i; matrix rotmat; /* rotation matrix */ @@ -932,7 +932,7 @@ static void do_single_flood(FILE* edo, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, - x, + as_rvec_array(coords.data()), edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, @@ -948,7 +948,7 @@ static void do_single_flood(FILE* edo, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, - x, + as_rvec_array(coords.data()), edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, @@ -1028,14 +1028,14 @@ static void do_single_flood(FILE* edo, /* Main flooding routine, called from do_force */ -void do_flood(const t_commrec* cr, - const t_inputrec& ir, - const rvec x[], - rvec force[], - gmx_edsam* ed, - const matrix box, - int64_t step, - bool bNS) +void do_flood(const t_commrec* cr, + const t_inputrec& ir, + gmx::ArrayRef coords, + gmx::ArrayRef force, + gmx_edsam* ed, + const matrix box, + int64_t step, + bool bNS) { /* Write time to edo, when required. Output the time anyhow since we need * it in the output file for ED constraints. */ @@ -1054,7 +1054,7 @@ void do_flood(const t_commrec* cr, /* Call flooding for one matrix */ if (edi.flood.vecs.neig) { - do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS); + do_single_flood(ed->edo, coords, force, &edi, step, box, cr, bNS); } } } @@ -3142,7 +3142,13 @@ std::unique_ptr init_edsam(const gmx::MDLogger& m } -void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[], rvec v[], const matrix box, gmx_edsam* ed) +void do_edsam(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + gmx::ArrayRef coords, + gmx::ArrayRef velocities, + const matrix box, + gmx_edsam* ed) { int i, edinr, iupdate = 500; matrix rotmat; /* rotation matrix */ @@ -3191,7 +3197,7 @@ void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[] buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, - xs, + as_rvec_array(coords.data()), edi.sav.nr, edi.sav.nr_loc, edi.sav.anrs_loc, @@ -3207,7 +3213,7 @@ void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[] buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, - xs, + as_rvec_array(coords.data()), edi.sref.nr, edi.sref.nr_loc, edi.sref.anrs_loc, @@ -3313,17 +3319,17 @@ void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[] box, buf->xcoll[edi.sav.c_ind[i]], buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh); /* dx is the ED correction to the positions: */ - rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx); + rvec_sub(x_unsh, coords[edi.sav.anrs_loc[i]], dx); - if (v != nullptr) + if (!velocities.empty()) { /* dv is the ED correction to the velocity: */ svmul(dt_1, dx, dv); /* apply the velocity correction: */ - rvec_inc(v[edi.sav.anrs_loc[i]], dv); + rvec_inc(velocities[edi.sav.anrs_loc[i]], dv); } /* Finally apply the position correction due to ED: */ - copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]); + copy_rvec(x_unsh, coords[edi.sav.anrs_loc[i]]); } } } /* END of if ( bNeedDoEdsam(edi) ) */ diff --git a/src/gromacs/essentialdynamics/edsam.h b/src/gromacs/essentialdynamics/edsam.h index ab26dfaf66..957f8c8f02 100644 --- a/src/gromacs/essentialdynamics/edsam.h +++ b/src/gromacs/essentialdynamics/edsam.h @@ -52,7 +52,6 @@ #include #include "gromacs/math/vectypes.h" -#include "gromacs/utility/basedefinitions.h" /*! \brief Abstract type for essential dynamics * @@ -72,6 +71,9 @@ namespace gmx { enum class StartingBehavior; class Constraints; +template +class ArrayRef; + class EssentialDynamics { public: @@ -97,18 +99,18 @@ class MDLogger; * \param ir MD input parameter record. * \param step Number of the time step. * \param cr Data needed for MPI communication. - * \param xs The local positions on this processor. - * \param v The local velocities. + * \param coords The local positions on this processor. + * \param velocities The local velocities. * \param box The simulation box. * \param ed The essential dynamics data. */ -void do_edsam(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - rvec xs[], - rvec v[], - const matrix box, - gmx_edsam* ed); +void do_edsam(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + gmx::ArrayRef coords, + gmx::ArrayRef velocities, + const matrix box, + gmx_edsam* ed); /*! \brief Initializes the essential dynamics and flooding module. @@ -153,20 +155,20 @@ void dd_make_local_ed_indices(gmx_domdec_t* dd, gmx_edsam* ed); * * \param cr Data needed for MPI communication. * \param ir MD input parameter record. - * \param x Positions on the local processor. + * \param coords Positions on the local processor. * \param force Forcefield forces to which the flooding forces are added. * \param ed The essential dynamics data. * \param box The simulation box. * \param step Number of the time step. * \param bNS Are we in a neighbor searching step? */ -void do_flood(const t_commrec* cr, - const t_inputrec& ir, - const rvec x[], - rvec force[], - gmx_edsam* ed, - const matrix box, - int64_t step, - bool bNS); +void do_flood(const t_commrec* cr, + const t_inputrec& ir, + gmx::ArrayRef coords, + gmx::ArrayRef force, + gmx_edsam* ed, + const matrix box, + int64_t step, + bool bNS); #endif diff --git a/src/gromacs/imd/imd.cpp b/src/gromacs/imd/imd.cpp index edc6a2fed9..c75f8d7ce0 100644 --- a/src/gromacs/imd/imd.cpp +++ b/src/gromacs/imd/imd.cpp @@ -203,9 +203,9 @@ public: /*! \brief Removes shifts of molecules diffused outside of the box. */ void removeMolecularShifts(const matrix box); /*! \brief Initialize arrays used to assemble the positions from the other nodes. */ - void prepareForPositionAssembly(const t_commrec* cr, const rvec x[]); + void prepareForPositionAssembly(const t_commrec* cr, gmx::ArrayRef coords); /*! \brief Interact with any connected VMD session */ - bool run(int64_t step, bool bNS, const matrix box, const rvec x[], double t); + bool run(int64_t step, bool bNS, const matrix box, gmx::ArrayRef coords, double t); // TODO rename all the data members to have underscore suffixes @@ -1242,7 +1242,7 @@ void ImdSession::Impl::removeMolecularShifts(const matrix box) } -void ImdSession::Impl::prepareForPositionAssembly(const t_commrec* cr, const rvec x[]) +void ImdSession::Impl::prepareForPositionAssembly(const t_commrec* cr, gmx::ArrayRef coords) { snew(xa, nat); snew(xa_ind, nat); @@ -1257,7 +1257,7 @@ void ImdSession::Impl::prepareForPositionAssembly(const t_commrec* cr, const rve for (int i = 0; i < nat; i++) { int ii = ind[i]; - copy_rvec(x[ii], xa_old[i]); + copy_rvec(coords[ii], xa_old[i]); } } @@ -1297,19 +1297,19 @@ static void imd_check_integrator_parallel(const t_inputrec* ir, const t_commrec* } } -std::unique_ptr makeImdSession(const t_inputrec* ir, - const t_commrec* cr, - gmx_wallcycle* wcycle, - gmx_enerdata_t* enerd, - const gmx_multisim_t* ms, - const gmx_mtop_t& top_global, - const MDLogger& mdlog, - const rvec x[], - int nfile, - const t_filenm fnm[], - const gmx_output_env_t* oenv, - const ImdOptions& options, - const gmx::StartingBehavior startingBehavior) +std::unique_ptr makeImdSession(const t_inputrec* ir, + const t_commrec* cr, + gmx_wallcycle* wcycle, + gmx_enerdata_t* enerd, + const gmx_multisim_t* ms, + const gmx_mtop_t& top_global, + const MDLogger& mdlog, + gmx::ArrayRef coords, + int nfile, + const t_filenm fnm[], + const gmx_output_env_t* oenv, + const ImdOptions& options, + const gmx::StartingBehavior startingBehavior) { std::unique_ptr session(new ImdSession(mdlog)); auto impl = session->impl_.get(); @@ -1505,7 +1505,7 @@ std::unique_ptr makeImdSession(const t_inputrec* ir, impl->syncNodes(cr, 0); /* Initialize arrays used to assemble the positions from the other nodes */ - impl->prepareForPositionAssembly(cr, x); + impl->prepareForPositionAssembly(cr, coords); /* Initialize molecule blocks to make them whole later...*/ if (MASTER(cr)) @@ -1517,7 +1517,7 @@ std::unique_ptr makeImdSession(const t_inputrec* ir, } -bool ImdSession::Impl::run(int64_t step, bool bNS, const matrix box, const rvec x[], double t) +bool ImdSession::Impl::run(int64_t step, bool bNS, const matrix box, gmx::ArrayRef coords, double t) { /* IMD at all? */ if (!sessionPossible) @@ -1567,7 +1567,7 @@ bool ImdSession::Impl::run(int64_t step, bool bNS, const matrix box, const rvec /* Transfer the IMD positions to the master node. Every node contributes * its local positions x and stores them in the assembled xa array. */ communicate_group_positions( - cr, xa, xa_shifts, xa_eshifts, true, x, nat, nat_loc, ind_loc, xa_ind, xa_old, box); + cr, xa, xa_shifts, xa_eshifts, true, as_rvec_array(coords.data()), nat, nat_loc, ind_loc, xa_ind, xa_old, box); /* If connected and master -> remove shifts */ if ((imdstep && bConnected) && MASTER(cr)) @@ -1581,9 +1581,9 @@ bool ImdSession::Impl::run(int64_t step, bool bNS, const matrix box, const rvec return imdstep; } -bool ImdSession::run(int64_t step, bool bNS, const matrix box, const rvec x[], double t) +bool ImdSession::run(int64_t step, bool bNS, const matrix box, gmx::ArrayRef coords, double t) { - return impl_->run(step, bNS, box, x, t); + return impl_->run(step, bNS, box, coords, t); } void ImdSession::fillEnergyRecord(int64_t step, bool bHaveNewEnergies) @@ -1655,7 +1655,7 @@ void ImdSession::updateEnergyRecordAndSendPositionsAndEnergies(bool bIMDstep, in wallcycle_stop(impl_->wcycle, ewcIMD); } -void ImdSession::applyForces(rvec* f) +void ImdSession::applyForces(gmx::ArrayRef force) { if (!impl_->sessionPossible || !impl_->bForceActivated) { @@ -1677,7 +1677,7 @@ void ImdSession::applyForces(rvec* f) j = *locndx; } - rvec_inc(f[j], impl_->f[i]); + rvec_inc(force[j], impl_->f[i]); } wallcycle_stop(impl_->wcycle, ewcIMD); diff --git a/src/gromacs/imd/imd.h b/src/gromacs/imd/imd.h index 6c6517717c..1d5c4fa56d 100644 --- a/src/gromacs/imd/imd.h +++ b/src/gromacs/imd/imd.h @@ -67,7 +67,6 @@ #include #include "gromacs/math/vectypes.h" -#include "gromacs/utility/basedefinitions.h" struct gmx_domdec_t; struct gmx_enerdata_t; @@ -90,6 +89,8 @@ class InteractiveMolecularDynamics; class MDLogger; struct ImdOptions; struct MdrunOptions; +template +class ArrayRef; /*! \brief * Creates a module for interactive molecular dynamics. @@ -130,26 +131,26 @@ void write_IMDgroup_to_file(bool bIMD, * \param ms Handler for multi-simulations. * \param top_global The topology of the whole system. * \param mdlog Logger - * \param x The starting positions of the atoms. + * \param coords The starting positions of the atoms. * \param nfile Number of files. * \param fnm Struct containing file names etc. * \param oenv Output options. * \param options Options for interactive MD. * \param startingBehavior Describes whether this is a restart appending to output files */ -std::unique_ptr makeImdSession(const t_inputrec* ir, - const t_commrec* cr, - gmx_wallcycle* wcycle, - gmx_enerdata_t* enerd, - const gmx_multisim_t* ms, - const gmx_mtop_t& top_global, - const MDLogger& mdlog, - const rvec x[], - int nfile, - const t_filenm fnm[], - const gmx_output_env_t* oenv, - const ImdOptions& options, - StartingBehavior startingBehavior); +std::unique_ptr makeImdSession(const t_inputrec* ir, + const t_commrec* cr, + gmx_wallcycle* wcycle, + gmx_enerdata_t* enerd, + const gmx_multisim_t* ms, + const gmx_mtop_t& top_global, + const MDLogger& mdlog, + gmx::ArrayRef coords, + int nfile, + const t_filenm fnm[], + const gmx_output_env_t* oenv, + const ImdOptions& options, + StartingBehavior startingBehavior); class ImdSession { @@ -179,18 +180,18 @@ public: * \param step The time step. * \param bNS Is this a neighbor searching step? * \param box The simulation box. - * \param x The local atomic positions on this node. + * \param coords The local atomic positions on this node. * \param t The time. * * \returns Whether or not we have to do IMD communication at this step. */ - bool run(int64_t step, bool bNS, const matrix box, const rvec x[], double t); + bool run(int64_t step, bool bNS, const matrix box, gmx::ArrayRef coords, double t); /*! \brief Add external forces from a running interactive molecular dynamics session. * - * \param f The forces. + * \param force The forces. */ - void applyForces(rvec* f); + void applyForces(gmx::ArrayRef force); /*! \brief Copy energies and convert to float from enerdata to the IMD energy record. * @@ -220,19 +221,19 @@ private: public: // Befriend the factory function. - friend std::unique_ptr makeImdSession(const t_inputrec* ir, - const t_commrec* cr, - gmx_wallcycle* wcycle, - gmx_enerdata_t* enerd, - const gmx_multisim_t* ms, - const gmx_mtop_t& top_global, - const MDLogger& mdlog, - const rvec x[], - int nfile, - const t_filenm fnm[], - const gmx_output_env_t* oenv, - const ImdOptions& options, - StartingBehavior startingBehavior); + friend std::unique_ptr makeImdSession(const t_inputrec* ir, + const t_commrec* cr, + gmx_wallcycle* wcycle, + gmx_enerdata_t* enerd, + const gmx_multisim_t* ms, + const gmx_mtop_t& top_global, + const MDLogger& mdlog, + gmx::ArrayRef coords, + int nfile, + const t_filenm fnm[], + const gmx_output_env_t* oenv, + const ImdOptions& options, + StartingBehavior startingBehavior); }; } // namespace gmx diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index d8cddba26e..398841179c 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -778,13 +778,7 @@ bool Constraints::Impl::apply(bool bLog, if (ed && delta_step > 0) { /* apply the essential dynamics constraints here */ - do_edsam(&ir, - step, - cr, - as_rvec_array(xprime.unpaddedArrayRef().data()), - as_rvec_array(v.unpaddedArrayRef().data()), - box, - ed); + do_edsam(&ir, step, cr, xprime.unpaddedArrayRef(), v.unpaddedArrayRef(), box, ed); } } wallcycle_stop(wcycle, ewcCONSTR); diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 424e9844fc..362cbaccf4 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -693,14 +693,12 @@ static void computeSpecialForces(FILE* fplog, fplog); } } - - rvec* f = as_rvec_array(forceWithVirialMtsLevel0->force_.data()); - /* Add the forces from enforced rotation potentials (if any) */ if (inputrec.bRot) { wallcycle_start(wcycle, ewcROTadd); - enerd->term[F_COM_PULL] += add_rot_forces(enforcedRotation, f, cr, step, t); + enerd->term[F_COM_PULL] += + add_rot_forces(enforcedRotation, forceWithVirialMtsLevel0->force_, cr, step, t); wallcycle_stop(wcycle, ewcROTadd); } @@ -711,13 +709,13 @@ static void computeSpecialForces(FILE* fplog, * Thus if no other algorithm (e.g. PME) requires it, the forces * here will contribute to the virial. */ - do_flood(cr, inputrec, as_rvec_array(x.data()), f, ed, box, step, didNeighborSearch); + do_flood(cr, inputrec, x, forceWithVirialMtsLevel0->force_, ed, box, step, didNeighborSearch); } /* Add forces from interactive molecular dynamics (IMD), if any */ if (inputrec.bIMD && stepWork.computeForces) { - imdSession->applyForces(f); + imdSession->applyForces(forceWithVirialMtsLevel0->force_); } } @@ -1697,7 +1695,7 @@ void do_force(FILE* fplog, if (inputrec.bRot) { wallcycle_start(wcycle, ewcROT); - do_rotation(cr, enforcedRotation, box, as_rvec_array(x.unpaddedArrayRef().data()), t, step, stepWork.doNeighborSearch); + do_rotation(cr, enforcedRotation, box, x.unpaddedConstArrayRef(), t, step, stepWork.doNeighborSearch); wallcycle_stop(wcycle, ewcROT); } diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 5584b60232..8c58a5bbc0 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -1359,7 +1359,7 @@ void gmx::LegacySimulator::do_md() mdrunOptions.writeConfout, bSumEkinhOld); /* Check if IMD step and do IMD communication, if bIMD is TRUE. */ - bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t); + bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x, t); /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */ if (startingBehavior != StartingBehavior::NewSimulation && bFirstStep diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 7c59a56005..f1bf119794 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -1777,7 +1777,7 @@ void LegacySimulator::do_cg() } /* Send energies and positions to the IMD client if bIMD is TRUE. */ - if (MASTER(cr) && imdSession->run(step, TRUE, state_global->box, state_global->x.rvec_array(), 0)) + if (MASTER(cr) && imdSession->run(step, TRUE, state_global->box, state_global->x, 0)) { imdSession->sendPositionsAndEnergies(); } @@ -2591,7 +2591,7 @@ void LegacySimulator::do_lbfgs() } /* Send x and E to IMD client, if bIMD is TRUE. */ - if (imdSession->run(step, TRUE, state_global->box, state_global->x.rvec_array(), 0) && MASTER(cr)) + if (imdSession->run(step, TRUE, state_global->box, state_global->x, 0) && MASTER(cr)) { imdSession->sendPositionsAndEnergies(); } @@ -2940,7 +2940,7 @@ void LegacySimulator::do_steep() if (imdSession->run(count, TRUE, MASTER(cr) ? state_global->box : nullptr, - MASTER(cr) ? state_global->x.rvec_array() : nullptr, + MASTER(cr) ? state_global->x : gmx::ArrayRef(), 0) && MASTER(cr)) { diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 488951c568..f732c60ede 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -1879,7 +1879,7 @@ int Mdrunner::mdrunner() ms, mtop, mdlog, - MASTER(cr) ? globalState->x.rvec_array() : nullptr, + MASTER(cr) ? globalState->x : gmx::ArrayRef(), filenames.size(), filenames.data(), oenv, diff --git a/src/gromacs/pulling/pull_rotation.cpp b/src/gromacs/pulling/pull_rotation.cpp index ab10a7a8a8..27adb6c158 100644 --- a/src/gromacs/pulling/pull_rotation.cpp +++ b/src/gromacs/pulling/pull_rotation.cpp @@ -609,7 +609,7 @@ static void reduce_output(const t_commrec* cr, gmx_enfrot* er, real t, int64_t s /* Add the forces from enforced rotation potential to the local forces. * Should be called after the SR forces have been evaluated */ -real add_rot_forces(gmx_enfrot* er, rvec f[], const t_commrec* cr, int64_t step, real t) +real add_rot_forces(gmx_enfrot* er, gmx::ArrayRef force, const t_commrec* cr, int64_t step, real t) { real Vrot = 0.0; /* If more than one rotation group is present, Vrot assembles the local parts from all groups */ @@ -626,7 +626,7 @@ real add_rot_forces(gmx_enfrot* er, rvec f[], const t_commrec* cr, int64_t step, /* Get the right index of the local force */ int ii = localRotationGroupIndex[l]; /* Add */ - rvec_inc(f[ii], erg->f_rot_loc[l]); + rvec_inc(force[ii], erg->f_rot_loc[l]); } } @@ -2010,12 +2010,12 @@ static void flex_precalc_inner_sum(const gmx_enfrotgrp* erg) } -static real do_flex2_lowlevel(gmx_enfrotgrp* erg, - real sigma, /* The Gaussian width sigma */ - rvec x[], - gmx_bool bOutstepRot, - gmx_bool bOutstepSlab, - const matrix box) +static real do_flex2_lowlevel(gmx_enfrotgrp* erg, + real sigma, /* The Gaussian width sigma */ + gmx::ArrayRef coords, + gmx_bool bOutstepRot, + gmx_bool bOutstepSlab, + const matrix box) { int count, ii, iigrp; rvec xj; /* position in the i-sum */ @@ -2076,7 +2076,7 @@ static real do_flex2_lowlevel(gmx_enfrotgrp* erg, * Note that erg->xc_center contains the center of mass in case the flex2-t * potential was chosen. For the flex2 potential erg->xc_center must be * zero. */ - rvec_sub(x[ii], erg->xc_center, xj); + rvec_sub(coords[ii], erg->xc_center, xj); /* Shift this atom such that it is near its reference */ shift_single_coord(box, xj, erg->xc_shifts[iigrp]); @@ -2278,11 +2278,11 @@ static real do_flex2_lowlevel(gmx_enfrotgrp* erg, static real do_flex_lowlevel(gmx_enfrotgrp* erg, - real sigma, /* The Gaussian width sigma */ - rvec x[], - gmx_bool bOutstepRot, - gmx_bool bOutstepSlab, - const matrix box) + real sigma, /* The Gaussian width sigma */ + gmx::ArrayRef coords, + gmx_bool bOutstepRot, + gmx_bool bOutstepSlab, + const matrix box) { int count, iigrp; rvec xj, yj0; /* current and reference position */ @@ -2335,7 +2335,7 @@ static real do_flex_lowlevel(gmx_enfrotgrp* erg, * Note that erg->xc_center contains the center of mass in case the flex-t * potential was chosen. For the flex potential erg->xc_center must be * zero. */ - rvec_sub(x[ii], erg->xc_center, xj); + rvec_sub(coords[ii], erg->xc_center, xj); /* Shift this atom such that it is near its reference */ shift_single_coord(box, xj, erg->xc_shifts[iigrp]); @@ -2608,11 +2608,11 @@ static void get_firstlast_slab_check(gmx_enfrotgrp* erg, /* The rotation group ( static void do_flexible(gmx_bool bMaster, gmx_enfrot* enfrot, /* Other rotation data */ gmx_enfrotgrp* erg, - rvec x[], /* The local positions */ - const matrix box, - double t, /* Time in picoseconds */ - gmx_bool bOutstepRot, /* Output to main rotation output file */ - gmx_bool bOutstepSlab) /* Output per-slab data */ + gmx::ArrayRef coords, /* The local positions */ + const matrix box, + double t, /* Time in picoseconds */ + gmx_bool bOutstepRot, /* Output to main rotation output file */ + gmx_bool bOutstepSlab) /* Output per-slab data */ { int nslabs; real sigma; /* The Gaussian width sigma */ @@ -2646,12 +2646,12 @@ static void do_flexible(gmx_bool bMaster, if (erg->rotg->eType == EnforcedRotationGroupType::Flex || erg->rotg->eType == EnforcedRotationGroupType::Flext) { - erg->V = do_flex_lowlevel(erg, sigma, x, bOutstepRot, bOutstepSlab, box); + erg->V = do_flex_lowlevel(erg, sigma, coords, bOutstepRot, bOutstepSlab, box); } else if (erg->rotg->eType == EnforcedRotationGroupType::Flex2 || erg->rotg->eType == EnforcedRotationGroupType::Flex2t) { - erg->V = do_flex2_lowlevel(erg, sigma, x, bOutstepRot, bOutstepSlab, box); + erg->V = do_flex2_lowlevel(erg, sigma, coords, bOutstepRot, bOutstepSlab, box); } else { @@ -2936,11 +2936,11 @@ static void do_radial_motion(gmx_enfrotgrp* erg, /* Calculate the radial motion pivot-free potential and forces */ -static void do_radial_motion_pf(gmx_enfrotgrp* erg, - rvec x[], /* The positions */ - const matrix box, /* The simulation box */ - gmx_bool bOutstepRot, /* Output to main rotation output file */ - gmx_bool bOutstepSlab) /* Output per-slab data */ +static void do_radial_motion_pf(gmx_enfrotgrp* erg, + gmx::ArrayRef coords, /* The positions */ + const matrix box, /* The simulation box */ + gmx_bool bOutstepRot, /* Output to main rotation output file */ + gmx_bool bOutstepSlab) /* Output per-slab data */ { rvec xj; /* Current position */ rvec xj_xc; /* xj - xc */ @@ -3005,7 +3005,7 @@ static void do_radial_motion_pf(gmx_enfrotgrp* erg, wj = N_M * mj; /* Current position of this atom: x[ii][XX/YY/ZZ] */ - copy_rvec(x[ii], xj); + copy_rvec(coords[ii], xj); /* Shift this atom such that it is near its reference */ shift_single_coord(box, xj, erg->xc_shifts[iigrp]); @@ -3138,9 +3138,9 @@ static void radial_motion2_precalc_inner_sum(const gmx_enfrotgrp* erg, rvec inne /* Calculate the radial motion 2 potential and forces */ -static void do_radial_motion2(gmx_enfrotgrp* erg, - rvec x[], /* The positions */ - const matrix box, /* The simulation box */ +static void do_radial_motion2(gmx_enfrotgrp* erg, + gmx::ArrayRef coords, /* The positions */ + const matrix box, /* The simulation box */ gmx_bool bOutstepRot, /* Output to main rotation output file */ gmx_bool bOutstepSlab) /* Output per-slab data */ { @@ -3195,7 +3195,7 @@ static void do_radial_motion2(gmx_enfrotgrp* erg, mj = erg->mc[iigrp]; /* Current position of this atom: x[ii] */ - copy_rvec(x[ii], xj); + copy_rvec(coords[ii], xj); /* Shift this atom such that it is near its reference */ shift_single_coord(box, xj, erg->xc_shifts[iigrp]); @@ -3856,7 +3856,7 @@ static void rotate_local_reference(gmx_enfrotgrp* erg) /* Select the PBC representation for each local x position and store that * for later usage. We assume the right PBC image of an x is the one nearest to * its rotated reference */ -static void choose_pbc_image(rvec x[], gmx_enfrotgrp* erg, const matrix box, int npbcdim) +static void choose_pbc_image(gmx::ArrayRef coords, gmx_enfrotgrp* erg, const matrix box, int npbcdim) { const auto& localRotationGroupIndex = erg->atomSet->localIndex(); for (gmx::index i = 0; i < localRotationGroupIndex.ssize(); i++) @@ -3873,12 +3873,18 @@ static void choose_pbc_image(rvec x[], gmx_enfrotgrp* erg, const matrix box, int copy_rvec(erg->xr_loc[i], xref); rvec_inc(xref, erg->xc_ref_center); - copy_correct_pbc_image(x[ii], erg->x_loc_pbc[i], xref, box, npbcdim); + copy_correct_pbc_image(coords[ii], erg->x_loc_pbc[i], xref, box, npbcdim); } } -void do_rotation(const t_commrec* cr, gmx_enfrot* er, const matrix box, rvec x[], real t, int64_t step, gmx_bool bNS) +void do_rotation(const t_commrec* cr, + gmx_enfrot* er, + const matrix box, + gmx::ArrayRef coords, + real t, + int64_t step, + gmx_bool bNS) { gmx_bool outstep_slab, outstep_rot; gmx_bool bColl; @@ -3925,7 +3931,7 @@ void do_rotation(const t_commrec* cr, gmx_enfrot* er, const matrix box, rvec x[] erg->xc_shifts, erg->xc_eshifts, bNS, - x, + as_rvec_array(coords.data()), rotg->nat, erg->atomSet->numAtomsLocal(), erg->atomSet->localIndex().data(), @@ -3953,7 +3959,7 @@ void do_rotation(const t_commrec* cr, gmx_enfrot* er, const matrix box, rvec x[] /* Choose the nearest PBC images of the group atoms with respect * to the rotated reference positions */ - choose_pbc_image(x, erg, box, 3); + choose_pbc_image(coords, erg, box, 3); /* Get the center of the rotation group */ if ((rotg->eType == EnforcedRotationGroupType::Isopf) @@ -4016,11 +4022,11 @@ void do_rotation(const t_commrec* cr, gmx_enfrot* er, const matrix box, rvec x[] do_radial_motion(erg, outstep_rot, outstep_slab); break; case EnforcedRotationGroupType::Rmpf: - do_radial_motion_pf(erg, x, box, outstep_rot, outstep_slab); + do_radial_motion_pf(erg, coords, box, outstep_rot, outstep_slab); break; case EnforcedRotationGroupType::Rm2: case EnforcedRotationGroupType::Rm2pf: - do_radial_motion2(erg, x, box, outstep_rot, outstep_slab); + do_radial_motion2(erg, coords, box, outstep_rot, outstep_slab); break; case EnforcedRotationGroupType::Flext: case EnforcedRotationGroupType::Flex2t: @@ -4030,13 +4036,13 @@ void do_rotation(const t_commrec* cr, gmx_enfrot* er, const matrix box, rvec x[] get_center(erg->xc, erg->mc, rotg->nat, erg->xc_center); svmul(-1.0, erg->xc_center, transvec); translate_x(erg->xc, rotg->nat, transvec); - do_flexible(MASTER(cr), er, erg, x, box, t, outstep_rot, outstep_slab); + do_flexible(MASTER(cr), er, erg, coords, box, t, outstep_rot, outstep_slab); break; case EnforcedRotationGroupType::Flex: case EnforcedRotationGroupType::Flex2: /* Do NOT subtract the center of mass in the low level routines! */ clear_rvec(erg->xc_center); - do_flexible(MASTER(cr), er, erg, x, box, t, outstep_rot, outstep_slab); + do_flexible(MASTER(cr), er, erg, coords, box, t, outstep_rot, outstep_slab); break; default: gmx_fatal(FARGS, "No such rotation potential."); } diff --git a/src/gromacs/pulling/pull_rotation.h b/src/gromacs/pulling/pull_rotation.h index de5f1388a0..d367317dc7 100644 --- a/src/gromacs/pulling/pull_rotation.h +++ b/src/gromacs/pulling/pull_rotation.h @@ -54,7 +54,6 @@ #include #include "gromacs/math/vectypes.h" -#include "gromacs/utility/basedefinitions.h" struct gmx_domdec_t; struct gmx_enfrot; @@ -71,6 +70,8 @@ namespace gmx enum class StartingBehavior; class LocalAtomSetManager; struct MdrunOptions; +template +class ArrayRef; class EnforcedRotation { @@ -133,13 +134,19 @@ std::unique_ptr init_rot(FILE* fplog, * \param cr Pointer to MPI communication data. * \param er Pointer to the enforced rotation working data. * \param box Simulation box, needed to make group whole. - * \param x The positions of all the local particles. + * \param coords The positions of all the local particles. * \param t Time. * \param step The time step. * \param bNS After domain decomposition / neighbor searching several * local arrays have to be updated (masses, shifts) */ -void do_rotation(const t_commrec* cr, gmx_enfrot* er, const matrix box, rvec x[], real t, int64_t step, gmx_bool bNS); +void do_rotation(const t_commrec* cr, + gmx_enfrot* er, + const matrix box, + gmx::ArrayRef coords, + real t, + int64_t step, + bool bNS); /*! \brief Add the enforced rotation forces to the official force array. @@ -152,14 +159,14 @@ void do_rotation(const t_commrec* cr, gmx_enfrot* er, const matrix box, rvec x[] * the potential, the angle of the group(s), and torques). * * \param er Pointer to the enforced rotation working data. - * \param f The local forces to which the rotational forces have + * \param force The local forces to which the rotational forces have * to be added. * \param cr Pointer to MPI communication data. * \param step The time step, used for output. * \param t Time, used for output. * \returns The potential energy of the rotation potentials. */ -real add_rot_forces(gmx_enfrot* er, rvec f[], const t_commrec* cr, int64_t step, real t); +real add_rot_forces(gmx_enfrot* er, gmx::ArrayRef force, const t_commrec* cr, int64_t step, real t); #endif -- 2.22.0