src/gromacs/mdtypes/state.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,2021, 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
  39 /*! \file
  40  *
  41  * \brief
  42  * This file contains the definition of the microstate of the simulated system
  43  *
  44  * History of observables that needs to be checkpointed should be stored
  45  * in ObservablesHistory.
  46  * The state of the mdrun machinery that needs to be checkpointed is also
  47  * stored elsewhere.
  48  *
  49  * \author Berk Hess
  50  *
  51  * \inlibraryapi
  52  * \ingroup module_mdtypes
  53  */
  54
  55 #ifndef GMX_MDTYPES_STATE_H
  56 #define GMX_MDTYPES_STATE_H
  57
  58 #include <array>
  59 #include <memory>
  60 #include <vector>
  61
  62 #include "gromacs/gpu_utils/hostallocator.h"
  63 #include "gromacs/math/paddedvector.h"
  64 #include "gromacs/math/vectypes.h"
  65 #include "gromacs/mdtypes/md_enums.h"
  66 #include "gromacs/utility/enumerationhelpers.h"
  67 #include "gromacs/utility/real.h"
  68
  69 struct t_inputrec;
  70 struct t_lambda;
  71 enum class FreeEnergyPerturbationType;
  72
  73 namespace gmx
  74 {
  75 struct AwhHistory;
  76 enum class CheckpointDataOperation;
  77 template<CheckpointDataOperation operation>
  78 class CheckpointData;
  79 } // namespace gmx
  80
  81 //! Convenience alias for until all is moved in the gmx namespace
  82 template<class T>
  83 using PaddedHostVector = gmx::PaddedHostVector<T>;
  84
  85 /*
  86  * The t_state struct should contain all the (possibly) non-static
  87  * information required to define the state of the system.
  88  * Currently the random seeds for SD and BD are missing.
  89  */
  90
  91 /*! \brief Enum for all entries in \p t_state
  92  *
  93  * These enums are used in flags as (1<<est...).
  94  * The order of these enums should not be changed,
  95  * since that affects the checkpoint (.cpt) file format.
  96  */
  97 enum class StateEntry : int
  98 {
  99     Lambda,
 100     Box,
 101     BoxRel,
 102     BoxV,
 103     PressurePrevious,
 104     Nhxi,
 105     ThermInt,
 106     X,
 107     V,
 108     SDxNotSupported,
 109     Cgp,
 110     LDRngNotSupported,
 111     LDRngINotSupported,
 112     DisreInitF,
 113     DisreRm3Tav,
 114     OrireInitF,
 115     OrireDtav,
 116     SVirPrev,
 117     Nhvxi,
 118     Veta,
 119     Vol0,
 120     Nhpresxi,
 121     Nhpresvxi,
 122     FVirPrev,
 123     FepState,
 124     MCRngNotSupported,
 125     MCRngINotSupported,
 126     BarosInt,
 127     PullComPrevStep,
 128     Count
 129 };
 130
 131 //! \brief The names of the state entries, defined in src/gromacs/fileio/checkpoint.cpp
 132 const char* enumValueToString(StateEntry enumValue);
 133 /*! \brief Convert enum to bitmask value.
 134  *
 135  * Used for setting flags in checkpoint header and verifying which flags are set.
 136  */
 137 template<typename Enum>
 138 inline int enumValueToBitMask(Enum enumValue)
 139 {
 140     static_assert(static_cast<int>(Enum::Count) <= std::numeric_limits<int>::digits);
 141     return 1 << static_cast<int>(enumValue);
 142 }
 143
 144 /*! \libinternal \brief History information for NMR distance and orientation restraints
 145  *
 146  * Often this is only used for reporting observables, and thus should not
 147  * actually be part of the microstate. But with time-dependent restraining
 148  * they are actually part of the (non-Markovian) microstate.
 149  * \todo Rename this with a more descriptive name.
 150  */
 151 class history_t
 152 {
 153 public:
 154     history_t();
 155
 156     real              disre_initf;  //!< The scaling factor for initializing the time av.
 157     std::vector<real> disre_rm3tav; //!< The r^-3 time averaged pair distances
 158     real              orire_initf;  //!< The scaling factor for initializing the time av.
 159     std::vector<real> orire_Dtav;   //!< The time averaged orientation tensors
 160 };
 161
 162 /*! \libinternal \brief Struct used for checkpointing only
 163  *
 164  * This struct would not be required with unlimited precision.
 165  * But because of limited precision, the COM motion removal implementation
 166  * can cause the kinetic energy in the MD loop to differ by a few bits from
 167  * the kinetic energy one would determine from state.v.
 168  */
 169 class ekinstate_t
 170 {
 171 public:
 172     ekinstate_t();
 173
 174     bool                bUpToDate;      //!< Test if all data is up to date
 175     int                 ekin_n;         //!< The number of tensors
 176     tensor*             ekinh;          //!< Half step Ekin, size \p ekin_n
 177     tensor*             ekinf;          //!< Full step Ekin, size \p ekin_n
 178     tensor*             ekinh_old;      //!< Half step Ekin of the previous step, size \p ekin_n
 179     tensor              ekin_total;     //!< Total kinetic energy
 180     std::vector<double> ekinscalef_nhc; //!< Nose-Hoover Ekin scaling factors for full step Ekin
 181     std::vector<double> ekinscaleh_nhc; //!< Nose-Hoover Ekin scaling factors for half step Ekin
 182     std::vector<double> vscale_nhc;     //!< Nose-Hoover velocity scaling factors
 183     real                dekindl;        //!< dEkin/dlambda, with free-energy
 184     real                mvcos; //!< Cosine(z) component of the momentum, for viscosity calculations
 185     /*! \brief Whether KE terms have been read from the checkpoint.
 186      *
 187      * Only used for managing whether the call to compute_globals
 188      * before we enter the MD loop should compute these quantities
 189      * fresh, or not. */
 190     bool hasReadEkinState;
 191
 192     /*!
 193      * \brief Allows to read and write checkpoint within modular simulator
 194      * \tparam operation  Whether we're reading or writing
 195      * \param checkpointData  The CheckpointData object
 196      */
 197     template<gmx::CheckpointDataOperation operation>
 198     void doCheckpoint(gmx::CheckpointData<operation> checkpointData);
 199 };
 200
 201 /*! \brief Free-energy sampling history struct
 202  *
 203  * \todo Split out into microstate and observables history.
 204  */
 205 struct df_history_t
 206 {
 207     int nlambda; //!< total number of lambda states - for history
 208
 209     bool  bEquil;   //!< Have we reached equilibration
 210     int*  n_at_lam; //!< number of points observed at each lambda
 211     real* wl_histo; //!< histogram for WL flatness determination
 212     real  wl_delta; //!< current wang-landau delta
 213
 214     real* sum_weights; //!< weights of the states
 215     real* sum_dg; //!< free energies of the states -- not actually used for weighting, but informational
 216     real* sum_minvar;   //!< corrections to weights for minimum variance
 217     real* sum_variance; //!< variances of the states
 218
 219     real** accum_p;  //!< accumulated bennett weights for n+1
 220     real** accum_m;  //!< accumulated bennett weights for n-1
 221     real** accum_p2; //!< accumulated squared bennett weights for n+1
 222     real** accum_m2; //!< accumulated squared bennett weights for n-1
 223
 224     real** Tij;           //!< transition matrix
 225     real** Tij_empirical; //!< Empirical transition matrix
 226
 227     /*! \brief Allows to read and write checkpoint within modular simulator
 228      *
 229      * \tparam operation  Whether we're reading or writing
 230      * \param checkpointData  The CheckpointData object
 231      * \param elamstats  How the lambda weights are calculated
 232      */
 233     template<gmx::CheckpointDataOperation operation>
 234     void doCheckpoint(gmx::CheckpointData<operation> checkpointData, LambdaWeightCalculation elamstats);
 235 };
 236
 237
 238 /*! \brief The microstate of the system
 239  *
 240  * The global state will contain complete data for all used entries.
 241  * The local state with domain decomposition will have partial entries
 242  * for which \p stateEntryIsAtomProperty() is true. Some entries that
 243  * are used in the global state might not be present in the local state.
 244  * \todo Move pure observables history to ObservablesHistory.
 245  */
 246 class t_state
 247 {
 248 public:
 249     t_state();
 250
 251     // All things public
 252     int natoms; //!< Number of atoms, local + non-local; this is the size of \p x, \p v and \p cg_p, when used
 253     int ngtc;          //!< The number of temperature coupling groups
 254     int nnhpres;       //!< The number of NH-chains for the MTTK barostat (always 1 or 0)
 255     int nhchainlength; //!< The NH-chain length for temperature coupling and MTTK barostat
 256     int flags; //!< Set of bit-flags telling which entries are present, see enum at the top of the file
 257     int fep_state; //!< indicates which of the alchemical states we are in
 258     gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> lambda; //!< Free-energy lambda vector
 259     matrix                                                          box; //!< Matrix of box vectors
 260     matrix              box_rel;        //!< Relative box vectors to preserve box shape
 261     matrix              boxv;           //!< Box velocities for Parrinello-Rahman P-coupling
 262     matrix              pres_prev;      //!< Pressure of the previous step for pcoupl
 263     matrix              svir_prev;      //!< Shake virial for previous step for pcoupl
 264     matrix              fvir_prev;      //!< Force virial of the previous step for pcoupl
 265     std::vector<double> nosehoover_xi;  //!< Nose-Hoover coordinates (ngtc)
 266     std::vector<double> nosehoover_vxi; //!< Nose-Hoover velocities (ngtc)
 267     std::vector<double> nhpres_xi;      //!< Pressure Nose-Hoover coordinates
 268     std::vector<double> nhpres_vxi;     //!< Pressure Nose-Hoover velocities
 269     std::vector<double> therm_integral; //!< Work exterted N-H/V-rescale T-coupling (ngtc)
 270     double              baros_integral; //!< For Berendsen P-coupling conserved quantity
 271     real                veta;           //!< Trotter based isotropic P-coupling
 272     real                vol0; //!< Initial volume,required for computing MTTK conserved quantity
 273     PaddedHostVector<gmx::RVec> x;    //!< The coordinates (natoms)
 274     PaddedHostVector<gmx::RVec> v;    //!< The velocities (natoms)
 275     PaddedHostVector<gmx::RVec> cg_p; //!< p vector for conjugate gradient minimization
 276
 277     ekinstate_t ekinstate; //!< The state of the kinetic energy
 278
 279     /* History for special algorithms, should be moved to a history struct */
 280     history_t                        hist;       //!< Time history for restraints
 281     df_history_t*                    dfhist;     //!< Free-energy history for free energy analysis
 282     std::shared_ptr<gmx::AwhHistory> awhHistory; //!< Accelerated weight histogram history
 283
 284     int              ddp_count;       //!< The DD partitioning count for this state
 285     int              ddp_count_cg_gl; //!< The DD partitioning count for index_gl
 286     std::vector<int> cg_gl;           //!< The global cg number of the local cgs
 287
 288     std::vector<double> pull_com_prev_step; //!< The COM of the previous step of each pull group
 289 };
 290
 291 #ifndef DOXYGEN
 292 /* We don't document the structs below, as they don't belong here.
 293  * TODO: Move the next two structs out of state.h.
 294  */
 295
 296 struct t_extmass
 297 {
 298     std::vector<double> Qinv; /* inverse mass of thermostat -- computed from inputs, but a good place to store */
 299     std::vector<double> QPinv; /* inverse mass of thermostat for barostat -- computed from inputs, but a good place to store */
 300     double              Winv; /* Pressure mass inverse -- computed, not input, but a good place to store. Need to make a matrix later */
 301     tensor              Winvm; /* inverse pressure mass tensor, computed       */
 302 };
 303
 304 #endif // DOXYGEN
 305
 306 //! Resizes the T- and P-coupling state variables
 307 void init_gtc_state(t_state* state, int ngtc, int nnhpres, int nhchainlength);
 308
 309 //! Change the number of atoms represented by this state, allocating memory as needed.
 310 void state_change_natoms(t_state* state, int natoms);
 311
 312 //! Allocates memory for free-energy history
 313 void init_dfhist_state(t_state* state, int dfhistNumLambda);
 314
 315 /*! \brief Compares two states, write the differences to stdout */
 316 void comp_state(const t_state* st1, const t_state* st2, bool bRMSD, real ftol, real abstol);
 317
 318 /*! \brief Allocates an rvec pointer and copy the contents of v to it */
 319 rvec* makeRvecArray(gmx::ArrayRef<const gmx::RVec> v, gmx::index n);
 320
 321 /*! \brief Determine the relative box components
 322  *
 323  * Set box_rel e.g. used in mdrun state, used to preserve the box shape
 324  * \param[in]    ir      Input record
 325  * \param[inout] state   State
 326  */
 327 void set_box_rel(const t_inputrec* ir, t_state* state);
 328
 329 /*! \brief Make sure the relative box shape remains the same
 330  *
 331  * This function ensures that the relative box dimensions are
 332  * preserved, which otherwise might diffuse away due to rounding
 333  * errors in pressure coupling or the deform option.
 334  *
 335  * \param[in]    ir      Input record
 336  * \param[in]    box_rel Relative box dimensions
 337  * \param[inout] box     The corrected actual box dimensions
 338  */
 339 void preserve_box_shape(const t_inputrec* ir, matrix box_rel, matrix box);
 340
 341 /*! \brief Returns an arrayRef to the positions in \p state when \p state!=null
 342  *
 343  * When \p state=nullptr, returns an empty arrayRef.
 344  *
 345  * \note The size returned is the number of atoms, without padding.
 346  *
 347  * \param[in] state  The state, can be nullptr
 348  */
 349 static inline gmx::ArrayRef<const gmx::RVec> positionsFromStatePointer(const t_state* state)
 350 {
 351     if (state)
 352     {
 353         return gmx::makeConstArrayRef(state->x).subArray(0, state->natoms);
 354     }
 355     else
 356     {
 357         return {};
 358     }
 359 };
 360
 361 /*! \brief Prints the current lambda state to the log file.
 362  *
 363  * \param[in] fplog  The log file. If fplog == nullptr there will be no output.
 364  * \param[in] lambda The array of lambda values.
 365  * \param[in] isInitialOutput Whether this output is the initial lambda state or not.
 366  */
 367 void printLambdaStateToLog(FILE* fplog, gmx::ArrayRef<const real> lambda, bool isInitialOutput);
 368
 369
 370 /*! \brief Fills fep_state and lambda if needed
 371  *
 372  * If FEP or simulated tempering is in use,  fills fep_state
 373  * and lambda on master rank.
 374  *
 375  * Reports the initial lambda state to the log file. */
 376 void initialize_lambdas(FILE*                      fplog,
 377                         FreeEnergyPerturbationType freeEnergyPerturbationType,
 378                         bool                       haveSimulatedTempering,
 379                         const t_lambda&            fep,
 380                         gmx::ArrayRef<const real>  simulatedTemperingTemps,
 381                         gmx::ArrayRef<real>        ref_t,
 382                         bool                       isMaster,
 383                         int*                       fep_state,
 384                         gmx::ArrayRef<real>        lambda);
 385
 386 #endif