*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
+ * 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.
#ifndef GMX_MDLIB_UPDATE_H
#define GMX_MDLIB_UPDATE_H
+#include <memory>
+
#include "gromacs/math/paddedvector.h"
#include "gromacs/math/vectypes.h"
+#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/utility/arrayref.h"
-#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/real.h"
class ekinstate_t;
-struct gmx_ekindata_t;
+class gmx_ekindata_t;
struct gmx_enerdata_t;
-struct t_extmass;
+enum class PbcType;
struct t_fcdata;
struct t_graph;
struct t_grpopts;
struct t_inputrec;
-struct t_mdatoms;
struct t_nrnb;
class t_state;
-
-/* Abstract type for update */
-struct gmx_update_t;
+enum class ParticleType;
namespace gmx
{
class BoxDeformation;
class Constraints;
-}
-
-/* Initialize the stochastic dynamics struct */
-gmx_update_t *init_update(const t_inputrec *ir,
- gmx::BoxDeformation *deform);
-
-/* Update pre-computed constants that depend on the reference
- * temperature for coupling.
- *
- * This could change e.g. in simulated annealing. */
-void update_temperature_constants(gmx_update_t *upd, const t_inputrec *ir);
-
-/* Update the size of per-atom arrays (e.g. after DD re-partitioning,
- which might increase the number of home atoms). */
-void update_realloc(gmx_update_t *upd, int natoms);
-
-/* Store the box at step step
- * as a reference state for simulations with box deformation.
- */
-void set_deform_reference_box(gmx_update_t *upd,
- int64_t step, matrix box);
-
-void update_tcouple(int64_t step,
- const t_inputrec *inputrec,
- t_state *state,
- gmx_ekindata_t *ekind,
- const t_extmass *MassQ,
- const t_mdatoms *md
- );
-
-/* Update Parrinello-Rahman, to be called before the coordinate update */
-void update_pcouple_before_coordinates(FILE *fplog,
- int64_t step,
- const t_inputrec *inputrec,
- t_state *state,
- matrix parrinellorahmanMu,
- matrix M,
- gmx_bool bInitStep);
-
-/* Update the box, to be called after the coordinate update.
- * For Berendsen P-coupling, also calculates the scaling factor
- * and scales the coordinates.
- * When the deform option is used, scales coordinates and box here.
- */
-void update_pcouple_after_coordinates(FILE *fplog,
- int64_t step,
- const t_inputrec *inputrec,
- const t_mdatoms *md,
- const matrix pressure,
- const matrix forceVirial,
- const matrix constraintVirial,
- const matrix parrinellorahmanMu,
- t_state *state,
- t_nrnb *nrnb,
- gmx_update_t *upd);
-
-void update_coords(int64_t step,
- const t_inputrec *inputrec, /* input record and box stuff */
- const t_mdatoms *md,
- t_state *state,
- gmx::ArrayRefWithPadding<gmx::RVec> f, /* forces on home particles */
- const t_fcdata *fcd,
- const gmx_ekindata_t *ekind,
- const matrix M,
- gmx_update_t *upd,
- int bUpdatePart,
- const t_commrec *cr, /* these shouldn't be here -- need to think about it */
- const gmx::Constraints *constr);
-
-/* Return TRUE if OK, FALSE in case of Shake Error */
-
-extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step, const t_commrec *cr,
- const t_mdatoms *md,
- gmx::ArrayRef<gmx::RVec> v,
- const gmx_update_t *upd,
- const gmx::Constraints *constr);
-
-void constrain_velocities(int64_t step,
- real *dvdlambda, /* the contribution to be added to the bonded interactions */
- t_state *state,
- tensor vir_part,
- gmx::Constraints *constr,
- gmx_bool bCalcVir,
- bool do_log,
- bool do_ene);
-
-void constrain_coordinates(int64_t step,
- real *dvdlambda, /* the contribution to be added to the bonded interactions */
- t_state *state,
- tensor vir_part,
- gmx_update_t *upd,
- gmx::Constraints *constr,
- gmx_bool bCalcVir,
- bool do_log,
- bool do_ene);
-
-void update_sd_second_half(int64_t step,
- real *dvdlambda, /* the contribution to be added to the bonded interactions */
- const t_inputrec *inputrec, /* input record and box stuff */
- const t_mdatoms *md,
- t_state *state,
- const t_commrec *cr,
- t_nrnb *nrnb,
- gmx_wallcycle_t wcycle,
- gmx_update_t *upd,
- gmx::Constraints *constr,
- bool do_log,
- bool do_ene);
-void finish_update(const t_inputrec *inputrec,
- const t_mdatoms *md,
- t_state *state,
- const t_graph *graph,
- t_nrnb *nrnb,
- gmx_wallcycle_t wcycle,
- gmx_update_t *upd,
- const gmx::Constraints *constr);
-/* Return TRUE if OK, FALSE in case of Shake Error */
+/*! \libinternal
+ * \brief Contains data for update phase */
+class Update
+{
+public:
+ /*! \brief Constructor
+ *
+ * \param[in] inputRecord Input record, used to construct SD object.
+ * \param[in] boxDeformation Periodic box deformation object.
+ */
+ Update(const t_inputrec& inputRecord, BoxDeformation* boxDeformation);
+ //! Destructor
+ ~Update();
+ /*! \brief Get the pointer to updated coordinates
+ *
+ * Update saves the updated coordinates into separate buffer, so that constraints will have
+ * access to both updated and not update coordinates. For that, update owns a separate buffer.
+ * See finish_update(...) for details.
+ *
+ * \returns The pointer to the intermediate coordinates buffer.
+ */
+ PaddedVector<gmx::RVec>* xp();
+ /*!\brief Getter to local copy of box deformation class.
+ *
+ * \returns handle to box deformation class
+ */
+ BoxDeformation* deform() const;
+ /*! \brief Sets data that changes only at domain decomposition time.
+ *
+ * \param[in] numAtoms Updated number of atoms.
+ * \param[in] cFREEZE Group index for freezing
+ * \param[in] cTC Group index for center of mass motion removal
+ * \param[in] cAcceleration Group index for constant acceleration groups
+ */
+ void updateAfterPartition(int numAtoms,
+ gmx::ArrayRef<const unsigned short> cFREEZE,
+ gmx::ArrayRef<const unsigned short> cTC,
+ gmx::ArrayRef<const unsigned short> cAcceleration);
+
+ /*! \brief Perform numerical integration step.
+ *
+ * Selects the appropriate integrator, based on the input record and performs a numerical integration step.
+ *
+ * \param[in] inputRecord Input record.
+ * \param[in] step Current timestep.
+ * \param[in] homenr The number of atoms on this processor.
+ * \param[in] havePartiallyFrozenAtoms Whether atoms are frozen along 1 or 2 (not 3) dimensions?
+ * \param[in] ptype The list of particle types.
+ * \param[in] invMass Inverse atomic mass per atom, 0 for vsites and shells.
+ * \param[in] invMassPerDim Inverse atomic mass per atom and dimension, 0 for vsites, shells and frozen dimensions
+ * \param[in] state System state object.
+ * \param[in] f Buffer with atomic forces for home particles.
+ * \param[in] fcdata Force calculation data to update distance and orientation restraints.
+ * \param[in] ekind Kinetic energy data (for temperature coupling, energy groups, etc.).
+ * \param[in] M Parrinello-Rahman velocity scaling matrix.
+ * \param[in] updatePart What should be updated, coordinates or velocities. This enum only used in VV integrator.
+ * \param[in] cr Comunication record (Old comment: these shouldn't be here -- need to think about it).
+ * \param[in] haveConstraints If the system has constraints.
+ */
+ void update_coords(const t_inputrec& inputRecord,
+ int64_t step,
+ int homenr,
+ bool havePartiallyFrozenAtoms,
+ gmx::ArrayRef<const ParticleType> ptype,
+ gmx::ArrayRef<const real> invMass,
+ gmx::ArrayRef<const rvec> invMassPerDim,
+ t_state* state,
+ const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
+ t_fcdata* fcdata,
+ const gmx_ekindata_t* ekind,
+ const matrix M,
+ int updatePart,
+ const t_commrec* cr,
+ bool haveConstraints);
+
+ /*! \brief Finalize the coordinate update.
+ *
+ * Copy the updated coordinates to the main coordinates buffer for the atoms that are not frozen.
+ *
+ * \param[in] inputRecord Input record.
+ * \param[in] havePartiallyFrozenAtoms Whether atoms are frozen along 1 or 2 (not 3) dimensions?
+ * \param[in] homenr The number of atoms on this processor.
+ * \param[in] state System state object.
+ * \param[in] wcycle Wall-clock cycle counter.
+ * \param[in] haveConstraints If the system has constraints.
+ */
+ void finish_update(const t_inputrec& inputRecord,
+ bool havePartiallyFrozenAtoms,
+ int homenr,
+ t_state* state,
+ gmx_wallcycle* wcycle,
+ bool haveConstraints);
+
+ /*! \brief Secong part of the SD integrator.
+ *
+ * The first part of integration is performed in the update_coords(...) method.
+ *
+ * \param[in] inputRecord Input record.
+ * \param[in] step Current timestep.
+ * \param[in] dvdlambda Free energy derivative. Contribution to be added to
+ * the bonded interactions.
+ * \param[in] homenr The number of atoms on this processor.
+ * \param[in] ptype The list of particle types.
+ * \param[in] invMass Inverse atomic mass per atom, 0 for vsites and shells.
+ * \param[in] state System state object.
+ * \param[in] cr Comunication record.
+ * \param[in] nrnb Cycle counters.
+ * \param[in] wcycle Wall-clock cycle counter.
+ * \param[in] constr Constraints object. The constraints are applied
+ * on coordinates after update.
+ * \param[in] do_log If this is logging step.
+ * \param[in] do_ene If this is an energy evaluation step.
+ */
+ void update_sd_second_half(const t_inputrec& inputRecord,
+ int64_t step,
+ real* dvdlambda,
+ int homenr,
+ gmx::ArrayRef<const ParticleType> ptype,
+ gmx::ArrayRef<const real> invMass,
+ t_state* state,
+ const t_commrec* cr,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle,
+ gmx::Constraints* constr,
+ bool do_log,
+ bool do_ene);
+
+ /*! \brief Performs a leap-frog update without updating \p state so the constrain virial
+ * can be computed.
+ */
+ void update_for_constraint_virial(const t_inputrec& inputRecord,
+ int homenr,
+ bool havePartiallyFrozenAtoms,
+ gmx::ArrayRef<const real> invmass,
+ gmx::ArrayRef<const rvec> invMassPerDim,
+ const t_state& state,
+ const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
+ const gmx_ekindata_t& ekind);
+
+ /*! \brief Update pre-computed constants that depend on the reference temperature for coupling.
+ *
+ * This could change e.g. in simulated annealing.
+ *
+ * \param[in] inputRecord Input record.
+ */
+ void update_temperature_constants(const t_inputrec& inputRecord);
+
+ /*!\brief Getter for the list of the randomize groups.
+ *
+ * Needed for Andersen temperature control.
+ *
+ * \returns Reference to the groups from the SD data object.
+ */
+ const std::vector<bool>& getAndersenRandomizeGroup() const;
+ /*!\brief Getter for the list of the Boltzmann factors.
+ *
+ * Needed for Andersen temperature control.
+ *
+ * \returns Reference to the Boltzmann factors from the SD data object.
+ */
+ const std::vector<real>& getBoltzmanFactor() const;
+
+private:
+ //! Implementation type.
+ class Impl;
+ //! Implementation object.
+ std::unique_ptr<Impl> impl_;
+};
+
+}; // namespace gmx
-void calc_ke_part(const t_state *state, const t_grpopts *opts, const t_mdatoms *md,
- gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel);
/*
* Compute the partial kinetic energy for home particles;
* will be accumulated in the calling routine.
*/
-void
-init_ekinstate(ekinstate_t *ekinstate, const t_inputrec *ir);
+void init_ekinstate(ekinstate_t* ekinstate, const t_inputrec* ir);
-void
-update_ekinstate(ekinstate_t *ekinstate, const gmx_ekindata_t *ekind);
+void update_ekinstate(ekinstate_t* ekinstate, const gmx_ekindata_t* ekind);
/*! \brief Restores data from \p ekinstate to \p ekind, then broadcasts it
to the rest of the simulation */
-void
-restore_ekinstate_from_state(const t_commrec *cr,
- gmx_ekindata_t *ekind, const ekinstate_t *ekinstate);
-
-void berendsen_tcoupl(const t_inputrec *ir, const gmx_ekindata_t *ekind, real dt,
- std::vector<double> &therm_integral); //NOLINT(google-runtime-references)
-
-void andersen_tcoupl(const t_inputrec *ir, int64_t step,
- const t_commrec *cr, const t_mdatoms *md,
- gmx::ArrayRef<gmx::RVec> v,
- real rate, const gmx_bool *randomize, const real *boltzfac);
-
-void nosehoover_tcoupl(const t_grpopts *opts, const gmx_ekindata_t *ekind, real dt,
- double xi[], double vxi[], const t_extmass *MassQ);
-
-void trotter_update(const t_inputrec *ir, int64_t step, gmx_ekindata_t *ekind,
- const gmx_enerdata_t *enerd, t_state *state, const tensor vir, const t_mdatoms *md,
- const t_extmass *MassQ, const int * const *trotter_seqlist, int trotter_seqno);
-
-int **init_npt_vars(const t_inputrec *ir, t_state *state, t_extmass *Mass, gmx_bool bTrotter);
-
-real NPT_energy(const t_inputrec *ir, const t_state *state, const t_extmass *MassQ);
-/* computes all the pressure/tempertature control energy terms to get a conserved energy */
-
-// TODO: This doesn't seem to be used or implemented anywhere
-void NBaroT_trotter(t_grpopts *opts, real dt,
- double xi[], double vxi[], real *veta, t_extmass *MassQ);
-
-void vrescale_tcoupl(const t_inputrec *ir, int64_t step,
- gmx_ekindata_t *ekind, real dt,
- double therm_integral[]);
-/* Compute temperature scaling. For V-rescale it is done in update. */
+void restore_ekinstate_from_state(const t_commrec* cr, gmx_ekindata_t* ekind, const ekinstate_t* ekinstate);
-void rescale_velocities(const gmx_ekindata_t *ekind, const t_mdatoms *mdatoms,
- int start, int end, rvec v[]);
-/* Rescale the velocities with the scaling factor in ekind */
-
-// TODO: This is the only function in update.h altering the inputrec
-void update_annealing_target_temp(t_inputrec *ir, real t, gmx_update_t *upd);
-/* Set reference temp for simulated annealing at time t*/
-
-real calc_temp(real ekin, real nrdf);
-/* Calculate the temperature */
-
-real calc_pres(int ePBC, int nwall, const matrix box, const tensor ekin, const tensor vir,
- tensor pres);
-/* Calculate the pressure tensor, returns the scalar pressure.
- * The unit of pressure is bar.
+/*! \brief Computes the atom range for a thread to operate on, ensuring SIMD aligned ranges
+ *
+ * \param[in] numThreads The number of threads to divide atoms over
+ * \param[in] threadIndex The thread to get the range for
+ * \param[in] numAtoms The total number of atoms (on this rank)
+ * \param[out] startAtom The start of the atom range
+ * \param[out] endAtom The end of the atom range, note that this is in general not a multiple of the SIMD width
*/
-
-void parrinellorahman_pcoupl(FILE *fplog, int64_t step,
- const t_inputrec *ir, real dt, const tensor pres,
- const tensor box, tensor box_rel, tensor boxv,
- tensor M, matrix mu,
- gmx_bool bFirstStep);
-
-void berendsen_pcoupl(FILE *fplog, int64_t step,
- const t_inputrec *ir, real dt,
- const tensor pres, const matrix box,
- const matrix force_vir, const matrix constraint_vir,
- matrix mu, double *baros_integral);
-
-void berendsen_pscale(const t_inputrec *ir, const matrix mu,
- matrix box, matrix box_rel,
- int start, int nr_atoms,
- rvec x[], const unsigned short cFREEZE[],
- t_nrnb *nrnb);
-
-// TODO: This doesn't seem to be used or implemented anywhere
-void correct_ekin(FILE *log, int start, int end, rvec v[],
- rvec vcm, real mass[], real tmass, tensor ekin);
-/* Correct ekin for vcm */
+void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* startAtom, int* endAtom);
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff