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