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