src/gromacs/mdlib/update.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifndef GMX_MDLIB_UPDATE_H
  39 #define GMX_MDLIB_UPDATE_H
  40
  41 #include "gromacs/math/paddedvector.h"
  42 #include "gromacs/math/vectypes.h"
  43 #include "gromacs/mdtypes/md_enums.h"
  44 #include "gromacs/timing/wallcycle.h"
  45 #include "gromacs/utility/arrayref.h"
  46 #include "gromacs/utility/basedefinitions.h"
  47 #include "gromacs/utility/classhelpers.h"
  48 #include "gromacs/utility/real.h"
  49
  50 class ekinstate_t;
  51 struct gmx_ekindata_t;
  52 struct gmx_enerdata_t;
  53 enum class PbcType;
  54 struct t_extmass;
  55 struct t_fcdata;
  56 struct t_graph;
  57 struct t_grpopts;
  58 struct t_inputrec;
  59 struct t_mdatoms;
  60 struct t_nrnb;
  61 class t_state;
  62
  63 namespace gmx
  64 {
  65 class BoxDeformation;
  66 class Constraints;
  67
  68
  69 /*! \libinternal
  70  * \brief Contains data for update phase */
  71 class Update
  72 {
  73 public:
  74     /*! \brief Constructor
  75      *
  76      * \param[in] inputRecord     Input record, used to construct SD object.
  77      * \param[in] boxDeformation  Periodic box deformation object.
  78      */
  79     Update(const t_inputrec& inputRecord, BoxDeformation* boxDeformation);
  80     //! Destructor
  81     ~Update();
  82     /*! \brief Get the pointer to updated coordinates
  83      *
  84      * Update saves the updated coordinates into separate buffer, so that constraints will have
  85      * access to both updated and not update coordinates. For that, update owns a separate buffer.
  86      * See finish_update(...) for details.
  87      *
  88      * \returns The pointer to the intermediate coordinates buffer.
  89      */
  90     PaddedVector<gmx::RVec>* xp();
  91     /*!\brief Getter to local copy of box deformation class.
  92      *
  93      * \returns handle to box deformation class
  94      */
  95     BoxDeformation* deform() const;
  96     /*! \brief Resizes buffer that stores intermediate coordinates.
  97      *
  98      * \param[in] numAtoms  Updated number of atoms.
  99      */
 100     void setNumAtoms(int numAtoms);
 101
 102     /*! \brief Perform numerical integration step.
 103      *
 104      * Selects the appropriate integrator, based on the input record and performs a numerical integration step.
 105      *
 106      * \param[in]  inputRecord      Input record.
 107      * \param[in]  step             Current timestep.
 108      * \param[in]  md               MD atoms data.
 109      * \param[in]  state            System state object.
 110      * \param[in]  f                Buffer with atomic forces for home particles.
 111      * \param[in]  fcd              Force calculation data to update distance and orientation restraints.
 112      * \param[in]  ekind            Kinetic energy data (for temperature coupling, energy groups, etc.).
 113      * \param[in]  M                Parrinello-Rahman velocity scaling matrix.
 114      * \param[in]  updatePart       What should be updated, coordinates or velocities. This enum only used in VV integrator.
 115      * \param[in]  cr               Comunication record  (Old comment: these shouldn't be here -- need to think about it).
 116      * \param[in]  haveConstraints  If the system has constraints.
 117      */
 118     void update_coords(const t_inputrec&                                inputRecord,
 119                        int64_t                                          step,
 120                        const t_mdatoms*                                 md,
 121                        t_state*                                         state,
 122                        const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 123                        const t_fcdata*                                  fcd,
 124                        const gmx_ekindata_t*                            ekind,
 125                        const matrix                                     M,
 126                        int                                              updatePart,
 127                        const t_commrec*                                 cr,
 128                        bool                                             haveConstraints);
 129
 130     /*! \brief Finalize the coordinate update.
 131      *
 132      * Copy the updated coordinates to the main coordinates buffer for the atoms that are not frozen.
 133      *
 134      * \param[in]  inputRecord      Input record.
 135      * \param[in]  md               MD atoms data.
 136      * \param[in]  state            System state object.
 137      * \param[in]  wcycle           Wall-clock cycle counter.
 138      * \param[in]  haveConstraints  If the system has constraints.
 139      */
 140     void finish_update(const t_inputrec& inputRecord,
 141                        const t_mdatoms*  md,
 142                        t_state*          state,
 143                        gmx_wallcycle_t   wcycle,
 144                        bool              haveConstraints);
 145
 146     /*! \brief Secong part of the SD integrator.
 147      *
 148      * The first part of integration is performed in the update_coords(...) method.
 149      *
 150      * \param[in]  inputRecord  Input record.
 151      * \param[in]  step         Current timestep.
 152      * \param[in]  dvdlambda    Free energy derivative. Contribution to be added to the bonded
 153      * interactions. \param[in]  md           MD atoms data. \param[in]  state        System state
 154      * object. \param[in]  cr           Comunication record. \param[in]  nrnb         Cycle
 155      * counters. \param[in]  wcycle       Wall-clock cycle counter. \param[in]  constr Constraints
 156      * object. The constraints are applied on coordinates after update. \param[in]  do_log       If
 157      * this is logging step. \param[in]  do_ene       If this is an energy evaluation step.
 158      */
 159     void update_sd_second_half(const t_inputrec& inputRecord,
 160                                int64_t           step,
 161                                real*             dvdlambda,
 162                                const t_mdatoms*  md,
 163                                t_state*          state,
 164                                const t_commrec*  cr,
 165                                t_nrnb*           nrnb,
 166                                gmx_wallcycle_t   wcycle,
 167                                gmx::Constraints* constr,
 168                                bool              do_log,
 169                                bool              do_ene);
 170     /*! \brief Update pre-computed constants that depend on the reference temperature for coupling.
 171      *
 172      * This could change e.g. in simulated annealing.
 173      *
 174      * \param[in]  inputRecord  Input record.
 175      */
 176     void update_temperature_constants(const t_inputrec& inputRecord);
 177
 178     /*!\brief Getter for the list of the randomize groups.
 179      *
 180      *  Needed for Andersen temperature control.
 181      *
 182      * \returns Reference to the groups from the SD data object.
 183      */
 184     const std::vector<bool>& getAndersenRandomizeGroup() const;
 185     /*!\brief Getter for the list of the Boltzmann factors.
 186      *
 187      *  Needed for Andersen temperature control.
 188      *
 189      * \returns Reference to the Boltzmann factors from the SD data object.
 190      */
 191     const std::vector<real>& getBoltzmanFactor() const;
 192
 193 private:
 194     //! Implementation type.
 195     class Impl;
 196     //! Implementation object.
 197     PrivateImplPointer<Impl> impl_;
 198 };
 199
 200 }; // namespace gmx
 201
 202
 203 /* Update the size of per-atom arrays (e.g. after DD re-partitioning,
 204    which might increase the number of home atoms). */
 205
 206 void update_tcouple(int64_t           step,
 207                     const t_inputrec* inputrec,
 208                     t_state*          state,
 209                     gmx_ekindata_t*   ekind,
 210                     const t_extmass*  MassQ,
 211                     const t_mdatoms*  md);
 212
 213 /* Update Parrinello-Rahman, to be called before the coordinate update */
 214 void update_pcouple_before_coordinates(FILE*             fplog,
 215                                        int64_t           step,
 216                                        const t_inputrec* inputrec,
 217                                        t_state*          state,
 218                                        matrix            parrinellorahmanMu,
 219                                        matrix            M,
 220                                        gmx_bool          bInitStep);
 221
 222 /* Update the box, to be called after the coordinate update.
 223  * For Berendsen P-coupling, also calculates the scaling factor
 224  * and scales the coordinates.
 225  * When the deform option is used, scales coordinates and box here.
 226  */
 227 void update_pcouple_after_coordinates(FILE*                fplog,
 228                                       int64_t              step,
 229                                       const t_inputrec*    inputrec,
 230                                       const t_mdatoms*     md,
 231                                       const matrix         pressure,
 232                                       const matrix         forceVirial,
 233                                       const matrix         constraintVirial,
 234                                       matrix               pressureCouplingMu,
 235                                       t_state*             state,
 236                                       t_nrnb*              nrnb,
 237                                       gmx::BoxDeformation* boxDeformation,
 238                                       bool                 scaleCoordinates);
 239
 240 /* Return TRUE if OK, FALSE in case of Shake Error */
 241
 242 extern gmx_bool update_randomize_velocities(const t_inputrec*        ir,
 243                                             int64_t                  step,
 244                                             const t_commrec*         cr,
 245                                             const t_mdatoms*         md,
 246                                             gmx::ArrayRef<gmx::RVec> v,
 247                                             const gmx::Update*       upd,
 248                                             const gmx::Constraints*  constr);
 249
 250 /*
 251  * Compute the partial kinetic energy for home particles;
 252  * will be accumulated in the calling routine.
 253  * The tensor is
 254  *
 255  * Ekin = SUM(i) 0.5 m[i] v[i] (x) v[i]
 256  *
 257  *     use v[i] = v[i] - u[i] when calculating temperature
 258  *
 259  * u must be accumulated already.
 260  *
 261  * Now also computes the contribution of the kinetic energy to the
 262  * free energy
 263  *
 264  */
 265
 266
 267 void init_ekinstate(ekinstate_t* ekinstate, const t_inputrec* ir);
 268
 269 void update_ekinstate(ekinstate_t* ekinstate, const gmx_ekindata_t* ekind);
 270
 271 /*! \brief Restores data from \p ekinstate to \p ekind, then broadcasts it
 272    to the rest of the simulation */
 273 void restore_ekinstate_from_state(const t_commrec* cr, gmx_ekindata_t* ekind, const ekinstate_t* ekinstate);
 274
 275 void berendsen_tcoupl(const t_inputrec*    ir,
 276                       gmx_ekindata_t*      ekind,
 277                       real                 dt,
 278                       std::vector<double>& therm_integral); //NOLINT(google-runtime-references)
 279
 280 void andersen_tcoupl(const t_inputrec*         ir,
 281                      int64_t                   step,
 282                      const t_commrec*          cr,
 283                      const t_mdatoms*          md,
 284                      gmx::ArrayRef<gmx::RVec>  v,
 285                      real                      rate,
 286                      const std::vector<bool>&  randomize,
 287                      gmx::ArrayRef<const real> boltzfac);
 288
 289 void nosehoover_tcoupl(const t_grpopts*      opts,
 290                        const gmx_ekindata_t* ekind,
 291                        real                  dt,
 292                        double                xi[],
 293                        double                vxi[],
 294                        const t_extmass*      MassQ);
 295
 296 void trotter_update(const t_inputrec*               ir,
 297                     int64_t                         step,
 298                     gmx_ekindata_t*                 ekind,
 299                     const gmx_enerdata_t*           enerd,
 300                     t_state*                        state,
 301                     const tensor                    vir,
 302                     const t_mdatoms*                md,
 303                     const t_extmass*                MassQ,
 304                     gmx::ArrayRef<std::vector<int>> trotter_seqlist,
 305                     int                             trotter_seqno);
 306
 307 std::array<std::vector<int>, ettTSEQMAX>
 308 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* Mass, gmx_bool bTrotter);
 309
 310 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ);
 311 /* computes all the pressure/tempertature control energy terms to get a conserved energy */
 312
 313 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[]);
 314 /* Compute temperature scaling. For V-rescale it is done in update. */
 315
 316 void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[]);
 317 /* Rescale the velocities with the scaling factor in ekind */
 318
 319 //! Check whether we do simulated annealing.
 320 bool doSimulatedAnnealing(const t_inputrec* ir);
 321
 322 //! Initialize simulated annealing.
 323 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd);
 324
 325 // TODO: This is the only function in update.h altering the inputrec
 326 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd);
 327 /* Set reference temp for simulated annealing at time t*/
 328
 329 real calc_temp(real ekin, real nrdf);
 330 /* Calculate the temperature */
 331
 332 real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres);
 333 /* Calculate the pressure tensor, returns the scalar pressure.
 334  * The unit of pressure is bar.
 335  */
 336
 337 void parrinellorahman_pcoupl(FILE*             fplog,
 338                              int64_t           step,
 339                              const t_inputrec* ir,
 340                              real              dt,
 341                              const tensor      pres,
 342                              const tensor      box,
 343                              tensor            box_rel,
 344                              tensor            boxv,
 345                              tensor            M,
 346                              matrix            mu,
 347                              gmx_bool          bFirstStep);
 348
 349 void berendsen_pcoupl(FILE*             fplog,
 350                       int64_t           step,
 351                       const t_inputrec* ir,
 352                       real              dt,
 353                       const tensor      pres,
 354                       const matrix      box,
 355                       const matrix      force_vir,
 356                       const matrix      constraint_vir,
 357                       matrix            mu,
 358                       double*           baros_integral);
 359
 360 void berendsen_pscale(const t_inputrec*    ir,
 361                       const matrix         mu,
 362                       matrix               box,
 363                       matrix               box_rel,
 364                       int                  start,
 365                       int                  nr_atoms,
 366                       rvec                 x[],
 367                       const unsigned short cFREEZE[],
 368                       t_nrnb*              nrnb,
 369                       bool                 scaleCoordinates);
 370
 371 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir);
 372
 373 /*! \brief Computes the atom range for a thread to operate on, ensuring SIMD aligned ranges
 374  *
 375  * \param[in]  numThreads   The number of threads to divide atoms over
 376  * \param[in]  threadIndex  The thread to get the range for
 377  * \param[in]  numAtoms     The total number of atoms (on this rank)
 378  * \param[out] startAtom    The start of the atom range
 379  * \param[out] endAtom      The end of the atom range, note that this is in general not a multiple of the SIMD width
 380  */
 381 void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* startAtom, int* endAtom);
 382
 383 /*! \brief Generate a new kinetic energy for the v-rescale thermostat
 384  *
 385  * Generates a new value for the kinetic energy, according to
 386  * Bussi et al JCP (2007), Eq. (A7)
 387  *
 388  * This is used by update_tcoupl(), and by the VRescaleThermostat of the modular
 389  * simulator.
 390  * TODO: Move this to the VRescaleThermostat once the modular simulator becomes
 391  *       the default code path.
 392  *
 393  * @param kk     present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
 394  * @param sigma  target average value of the kinetic energy (ndeg k_b T/2)  (in the same units as kk)
 395  * @param ndeg   number of degrees of freedom of the atoms to be thermalized
 396  * @param taut   relaxation time of the thermostat, in units of 'how often this routine is called'
 397  * @param step   the time step this routine is called on
 398  * @param seed   the random number generator seed
 399  * @return  the new kinetic energy
 400  */
 401 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed);
 402
 403 #endif