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     int   ndisrepairs;  //!< The number of distance restraints
 156     real* disre_rm3tav; //!< The r^-3 time averaged pair distances
 157     real  orire_initf;  //!< The scaling factor for initializing the time av.
 158     int   norire_Dtav;  //!< The number of matrix element in dtav (npair*5)
 159     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 typedef 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 } df_history_t;
 228
 229
 230 /*! \brief The microstate of the system
 231  *
 232  * The global state will contain complete data for all used entries.
 233  * The local state with domain decomposition will have partial entries
 234  * for which \p stateEntryIsAtomProperty() is true. Some entries that
 235  * are used in the global state might not be present in the local state.
 236  * \todo Move pure observables history to ObservablesHistory.
 237  */
 238 class t_state
 239 {
 240 public:
 241     t_state();
 242
 243     // All things public
 244     int natoms; //!< Number of atoms, local + non-local; this is the size of \p x, \p v and \p cg_p, when used
 245     int ngtc;          //!< The number of temperature coupling groups
 246     int nnhpres;       //!< The number of NH-chains for the MTTK barostat (always 1 or 0)
 247     int nhchainlength; //!< The NH-chain length for temperature coupling and MTTK barostat
 248     int flags; //!< Set of bit-flags telling which entries are present, see enum at the top of the file
 249     int fep_state; //!< indicates which of the alchemical states we are in
 250     gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, real> lambda; //!< Free-energy lambda vector
 251     matrix                                                          box; //!< Matrix of box vectors
 252     matrix              box_rel;        //!< Relative box vectors to preserve box shape
 253     matrix              boxv;           //!< Box velocities for Parrinello-Rahman P-coupling
 254     matrix              pres_prev;      //!< Pressure of the previous step for pcoupl
 255     matrix              svir_prev;      //!< Shake virial for previous step for pcoupl
 256     matrix              fvir_prev;      //!< Force virial of the previous step for pcoupl
 257     std::vector<double> nosehoover_xi;  //!< Nose-Hoover coordinates (ngtc)
 258     std::vector<double> nosehoover_vxi; //!< Nose-Hoover velocities (ngtc)
 259     std::vector<double> nhpres_xi;      //!< Pressure Nose-Hoover coordinates
 260     std::vector<double> nhpres_vxi;     //!< Pressure Nose-Hoover velocities
 261     std::vector<double> therm_integral; //!< Work exterted N-H/V-rescale T-coupling (ngtc)
 262     double              baros_integral; //!< For Berendsen P-coupling conserved quantity
 263     real                veta;           //!< Trotter based isotropic P-coupling
 264     real                vol0; //!< Initial volume,required for computing MTTK conserved quantity
 265     PaddedHostVector<gmx::RVec> x;    //!< The coordinates (natoms)
 266     PaddedHostVector<gmx::RVec> v;    //!< The velocities (natoms)
 267     PaddedHostVector<gmx::RVec> cg_p; //!< p vector for conjugate gradient minimization
 268
 269     ekinstate_t ekinstate; //!< The state of the kinetic energy
 270
 271     /* History for special algorithms, should be moved to a history struct */
 272     history_t                        hist;       //!< Time history for restraints
 273     df_history_t*                    dfhist;     //!< Free-energy history for free energy analysis
 274     std::shared_ptr<gmx::AwhHistory> awhHistory; //!< Accelerated weight histogram history
 275
 276     int              ddp_count;       //!< The DD partitioning count for this state
 277     int              ddp_count_cg_gl; //!< The DD partitioning count for index_gl
 278     std::vector<int> cg_gl;           //!< The global cg number of the local cgs
 279
 280     std::vector<double> pull_com_prev_step; //!< The COM of the previous step of each pull group
 281 };
 282
 283 #ifndef DOXYGEN
 284 /* We don't document the structs below, as they don't belong here.
 285  * TODO: Move the next two structs out of state.h.
 286  */
 287
 288 struct t_extmass
 289 {
 290     std::vector<double> Qinv; /* inverse mass of thermostat -- computed from inputs, but a good place to store */
 291     std::vector<double> QPinv; /* inverse mass of thermostat for barostat -- computed from inputs, but a good place to store */
 292     double              Winv; /* Pressure mass inverse -- computed, not input, but a good place to store. Need to make a matrix later */
 293     tensor              Winvm; /* inverse pressure mass tensor, computed       */
 294 };
 295
 296
 297 typedef struct
 298 {
 299     real    veta;
 300     double  rscale;
 301     double  vscale;
 302     double  rvscale;
 303     double  alpha;
 304     double* vscale_nhc;
 305 } t_vetavars;
 306
 307 #endif // DOXYGEN
 308
 309 //! Resizes the T- and P-coupling state variables
 310 void init_gtc_state(t_state* state, int ngtc, int nnhpres, int nhchainlength);
 311
 312 //! Change the number of atoms represented by this state, allocating memory as needed.
 313 void state_change_natoms(t_state* state, int natoms);
 314
 315 //! Allocates memory for free-energy history
 316 void init_dfhist_state(t_state* state, int dfhistNumLambda);
 317
 318 /*! \brief Compares two states, write the differences to stdout */
 319 void comp_state(const t_state* st1, const t_state* st2, bool bRMSD, real ftol, real abstol);
 320
 321 /*! \brief Allocates an rvec pointer and copy the contents of v to it */
 322 rvec* makeRvecArray(gmx::ArrayRef<const gmx::RVec> v, gmx::index n);
 323
 324 /*! \brief Determine the relative box components
 325  *
 326  * Set box_rel e.g. used in mdrun state, used to preserve the box shape
 327  * \param[in]    ir      Input record
 328  * \param[inout] state   State
 329  */
 330 void set_box_rel(const t_inputrec* ir, t_state* state);
 331
 332 /*! \brief Make sure the relative box shape remains the same
 333  *
 334  * This function ensures that the relative box dimensions are
 335  * preserved, which otherwise might diffuse away due to rounding
 336  * errors in pressure coupling or the deform option.
 337  *
 338  * \param[in]    ir      Input record
 339  * \param[in]    box_rel Relative box dimensions
 340  * \param[inout] box     The corrected actual box dimensions
 341  */
 342 void preserve_box_shape(const t_inputrec* ir, matrix box_rel, matrix box);
 343
 344 /*! \brief Returns an arrayRef to the positions in \p state when \p state!=null
 345  *
 346  * When \p state=nullptr, returns an empty arrayRef.
 347  *
 348  * \note The size returned is the number of atoms, without padding.
 349  *
 350  * \param[in] state  The state, can be nullptr
 351  */
 352 static inline gmx::ArrayRef<const gmx::RVec> positionsFromStatePointer(const t_state* state)
 353 {
 354     if (state)
 355     {
 356         return gmx::makeConstArrayRef(state->x).subArray(0, state->natoms);
 357     }
 358     else
 359     {
 360         return {};
 361     }
 362 };
 363
 364 /*! \brief Prints the current lambda state to the log file.
 365  *
 366  * \param[in] fplog  The log file. If fplog == nullptr there will be no output.
 367  * \param[in] lambda The array of lambda values.
 368  * \param[in] isInitialOutput Whether this output is the initial lambda state or not.
 369  */
 370 void printLambdaStateToLog(FILE* fplog, gmx::ArrayRef<const real> lambda, bool isInitialOutput);
 371
 372
 373 /*! \brief Fills fep_state and lambda if needed
 374  *
 375  * If FEP or simulated tempering is in use,  fills fep_state
 376  * and lambda on master rank.
 377  *
 378  * Reports the initial lambda state to the log file. */
 379 void initialize_lambdas(FILE* fplog, const t_inputrec& ir, bool isMaster, int* fep_state, gmx::ArrayRef<real> lambda);
 380
 381 #endif