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,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 /*! \file
  39  *
  40  * \brief
  41  * This file contains the definition of the microstate of the simulated system
  42  *
  43  * History of observables that needs to be checkpointed should be stored
  44  * in ObservablesHistory.
  45  * The state of the mdrun machinery that needs to be checkpointed is also
  46  * stored elsewhere.
  47  *
  48  * \author Berk Hess
  49  *
  50  * \inlibraryapi
  51  * \ingroup module_mdtypes
  52  */
  53
  54 #ifndef GMX_MDTYPES_STATE_H
  55 #define GMX_MDTYPES_STATE_H
  56
  57 #include <array>
  58 #include <memory>
  59 #include <vector>
  60
  61 #include "gromacs/gpu_utils/hostallocator.h"
  62 #include "gromacs/math/paddedvector.h"
  63 #include "gromacs/math/vectypes.h"
  64 #include "gromacs/mdtypes/md_enums.h"
  65 #include "gromacs/utility/arrayref.h"
  66 #include "gromacs/utility/basedefinitions.h"
  67 #include "gromacs/utility/real.h"
  68
  69 struct t_inputrec;
  70
  71 namespace gmx
  72 {
  73 struct AwhHistory;
  74 }
  75
  76 //! Convenience alias for until all is moved in the gmx namespace
  77 template <class T>
  78 using PaddedHostVector = gmx::PaddedHostVector<T>;
  79
  80 /*
  81  * The t_state struct should contain all the (possibly) non-static
  82  * information required to define the state of the system.
  83  * Currently the random seeds for SD and BD are missing.
  84  */
  85
  86 /* \brief Enum for all entries in \p t_state
  87  *
  88  * These enums are used in flags as (1<<est...).
  89  * The order of these enums should not be changed,
  90  * since that affects the checkpoint (.cpt) file format.
  91  */
  92 enum {
  93     estLAMBDA,
  94     estBOX, estBOX_REL, estBOXV, estPRES_PREV, estNH_XI,  estTHERM_INT,
  95     estX,   estV,       estSDX_NOTSUPPORTED,  estCGP,
  96     estLD_RNG_NOTSUPPORTED, estLD_RNGI_NOTSUPPORTED,
  97     estDISRE_INITF, estDISRE_RM3TAV,
  98     estORIRE_INITF, estORIRE_DTAV,
  99     estSVIR_PREV, estNH_VXI, estVETA, estVOL0, estNHPRES_XI, estNHPRES_VXI, estFVIR_PREV,
 100     estFEPSTATE, estMC_RNG_NOTSUPPORTED, estMC_RNGI_NOTSUPPORTED,
 101     estBAROS_INT, estPULLCOMPREVSTEP,
 102     estNR
 103 };
 104
 105 //! \brief The names of the state entries, defined in src/gmxlib/checkpoint.c
 106 extern const char *est_names[estNR];
 107
 108 /*! \libinternal \brief History information for NMR distance and orientation restraints
 109  *
 110  * Often this is only used for reporting observables, and thus should not
 111  * actually be part of the microstate. But with time-dependent restraining
 112  * they are actually part of the (non-Markovian) microstate.
 113  * \todo Rename this with a more descriptive name.
 114  */
 115 class history_t
 116 {
 117     public:
 118         history_t();
 119
 120         real  disre_initf;  //!< The scaling factor for initializing the time av.
 121         int   ndisrepairs;  //!< The number of distance restraints
 122         real *disre_rm3tav; //!< The r^-3 time averaged pair distances
 123         real  orire_initf;  //!< The scaling factor for initializing the time av.
 124         int   norire_Dtav;  //!< The number of matrix element in dtav (npair*5)
 125         real *orire_Dtav;   //!< The time averaged orientation tensors
 126 };
 127
 128 /*! \libinternal \brief Struct used for checkpointing only
 129  *
 130  * This struct would not be required with unlimited precision.
 131  * But because of limited precision, the COM motion removal implementation
 132  * can cause the kinetic energy in the MD loop to differ by a few bits from
 133  * the kinetic energy one would determine from state.v.
 134  */
 135 class ekinstate_t
 136 {
 137     public:
 138         ekinstate_t();
 139
 140         gmx_bool             bUpToDate;      //!< Test if all data is up to date
 141         int                  ekin_n;         //!< The number of tensors
 142         tensor              *ekinh;          //!< Half step Ekin, size \p ekin_n
 143         tensor              *ekinf;          //!< Full step Ekin, size \p ekin_n
 144         tensor              *ekinh_old;      //!< Half step Ekin of the previous step, size \p ekin_n
 145         tensor               ekin_total;     //!< Total kinetic energy
 146         std::vector<double>  ekinscalef_nhc; //!< Nose-Hoover Ekin scaling factors for full step Ekin
 147         std::vector<double>  ekinscaleh_nhc; //!< Nose-Hoover Ekin scaling factors for half step Ekin
 148         std::vector<double>  vscale_nhc;     //!< Nose-Hoover velocity scaling factors
 149         real                 dekindl;        //!< dEkin/dlambda, with free-energy
 150         real                 mvcos;          //!< Cosine(z) component of the momentum, for viscosity calculations
 151         /*! \brief Whether KE terms have been read from the checkpoint.
 152          *
 153          * Only used for managing whether the call to compute_globals
 154          * before we enter the MD loop should compute these quantities
 155          * fresh, or not. */
 156         bool hasReadEkinState;
 157 };
 158
 159 /*! \brief Free-energy sampling history struct
 160  *
 161  * \todo Split out into microstate and observables history.
 162  */
 163 typedef struct df_history_t
 164 {
 165     int      nlambda;        //!< total number of lambda states - for history
 166
 167     gmx_bool bEquil;         //!< Have we reached equilibration
 168     int     *n_at_lam;       //!< number of points observed at each lambda
 169     real    *wl_histo;       //!< histogram for WL flatness determination
 170     real     wl_delta;       //!< current wang-landau delta
 171
 172     real    *sum_weights;    //!< weights of the states
 173     real    *sum_dg;         //!< free energies of the states -- not actually used for weighting, but informational
 174     real    *sum_minvar;     //!< corrections to weights for minimum variance
 175     real    *sum_variance;   //!< variances of the states
 176
 177     real   **accum_p;        //!< accumulated bennett weights for n+1
 178     real   **accum_m;        //!< accumulated bennett weights for n-1
 179     real   **accum_p2;       //!< accumulated squared bennett weights for n+1
 180     real   **accum_m2;       //!< accumulated squared bennett weights for n-1
 181
 182     real   **Tij;            //!< transition matrix
 183     real   **Tij_empirical;  //!< Empirical transition matrix
 184
 185 } df_history_t;
 186
 187
 188 /*! \brief The microstate of the system
 189  *
 190  * The global state will contain complete data for all used entries.
 191  * The local state with domain decomposition will have partial entries
 192  * for which \p stateEntryIsAtomProperty() is true. Some entries that
 193  * are used in the global state might not be present in the local state.
 194  * \todo Move pure observables history to ObservablesHistory.
 195  */
 196 class t_state
 197 {
 198     public:
 199         t_state();
 200
 201         // All things public
 202         int                             natoms;         //!< Number of atoms, local + non-local; this is the size of \p x, \p v and \p cg_p, when used
 203         int                             ngtc;           //!< The number of temperature coupling groups
 204         int                             nnhpres;        //!< The NH-chain length for the MTTK barostat
 205         int                             nhchainlength;  //!< The NH-chain length for temperature coupling
 206         int                             flags;          //!< Set of bit-flags telling which entries are present, see enum at the top of the file
 207         int                             fep_state;      //!< indicates which of the alchemical states we are in
 208         std::array<real, efptNR>        lambda;         //!< Free-energy lambda vector
 209         matrix                          box;            //!< Matrix of box vectors
 210         matrix                          box_rel;        //!< Relative box vectors to preserve box shape
 211         matrix                          boxv;           //!< Box velocities for Parrinello-Rahman P-coupling
 212         matrix                          pres_prev;      //!< Pressure of the previous step for pcoupl
 213         matrix                          svir_prev;      //!< Shake virial for previous step for pcoupl
 214         matrix                          fvir_prev;      //!< Force virial of the previous step for pcoupl
 215         std::vector<double>             nosehoover_xi;  //!< Nose-Hoover coordinates (ngtc)
 216         std::vector<double>             nosehoover_vxi; //!< Nose-Hoover velocities (ngtc)
 217         std::vector<double>             nhpres_xi;      //!< Pressure Nose-Hoover coordinates
 218         std::vector<double>             nhpres_vxi;     //!< Pressure Nose-Hoover velocities
 219         std::vector<double>             therm_integral; //!< Work exterted N-H/V-rescale T-coupling (ngtc)
 220         double                          baros_integral; //!< For Berendsen P-coupling conserved quantity
 221         real                            veta;           //!< Trotter based isotropic P-coupling
 222         real                            vol0;           //!< Initial volume,required for computing MTTK conserved quantity
 223         PaddedHostVector<gmx::RVec>     x;              //!< The coordinates (natoms)
 224         PaddedHostVector<gmx::RVec>     v;              //!< The velocities (natoms)
 225         PaddedHostVector<gmx::RVec>     cg_p;           //!< p vector for conjugate gradient minimization
 226
 227         ekinstate_t                     ekinstate;      //!< The state of the kinetic energy
 228
 229         /* History for special algorithms, should be moved to a history struct */
 230         history_t                         hist;               //!< Time history for restraints
 231         df_history_t                     *dfhist;             //!< Free-energy history for free energy analysis
 232         std::shared_ptr<gmx::AwhHistory>  awhHistory;         //!< Accelerated weight histogram history
 233
 234         int                               ddp_count;          //!< The DD partitioning count for this state
 235         int                               ddp_count_cg_gl;    //!< The DD partitioning count for index_gl
 236         std::vector<int>                  cg_gl;              //!< The global cg number of the local cgs
 237
 238         std::vector<double>               pull_com_prev_step; //!< The COM of the previous step of each pull group
 239 };
 240
 241 #ifndef DOXYGEN
 242 /* We don't document the structs below, as they don't belong here.
 243  * TODO: Move the next two structs out of state.h.
 244  */
 245
 246 struct t_extmass
 247 {
 248     std::vector<double> Qinv;  /* inverse mass of thermostat -- computed from inputs, but a good place to store */
 249     std::vector<double> QPinv; /* inverse mass of thermostat for barostat -- computed from inputs, but a good place to store */
 250     double              Winv;  /* Pressure mass inverse -- computed, not input, but a good place to store. Need to make a matrix later */
 251     tensor              Winvm; /* inverse pressure mass tensor, computed       */
 252 };
 253
 254
 255 typedef struct
 256 {
 257     real    veta;
 258     double  rscale;
 259     double  vscale;
 260     double  rvscale;
 261     double  alpha;
 262     double *vscale_nhc;
 263 } t_vetavars;
 264
 265 #endif // DOXYGEN
 266
 267 //! Resizes the T- and P-coupling state variables
 268 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength);
 269
 270 //! Change the number of atoms represented by this state, allocating memory as needed.
 271 void state_change_natoms(t_state *state, int natoms);
 272
 273 //! Allocates memory for free-energy history
 274 void init_dfhist_state(t_state *state, int dfhistNumLambda);
 275
 276 /*! \brief Compares two states, write the differences to stdout */
 277 void comp_state(const t_state *st1, const t_state *st2, gmx_bool bRMSD, real ftol, real abstol);
 278
 279 /*! \brief Allocates an rvec pointer and copy the contents of v to it */
 280 rvec *makeRvecArray(gmx::ArrayRef<const gmx::RVec> v,
 281                     gmx::index                     n);
 282
 283 /*! \brief Determine the relative box components
 284  *
 285  * Set box_rel e.g. used in mdrun state, used to preserve the box shape
 286  * \param[in]    ir      Input record
 287  * \param[inout] state   State
 288  */
 289 void set_box_rel(const t_inputrec *ir, t_state *state);
 290
 291 /*! \brief Make sure the relative box shape remains the same
 292  *
 293  * This function ensures that the relative box dimensions are
 294  * preserved, which otherwise might diffuse away due to rounding
 295  * errors in pressure coupling or the deform option.
 296  *
 297  * \param[in]    ir      Input record
 298  * \param[in]    box_rel Relative box dimensions
 299  * \param[inout] box     The corrected actual box dimensions
 300  */
 301 void preserve_box_shape(const t_inputrec *ir, matrix box_rel, matrix box);
 302
 303 /*! \brief Returns an arrayRef to the positions in \p state when \p state!=null
 304  *
 305  * When \p state=nullptr, returns an empty arrayRef.
 306  *
 307  * \note The size returned is the number of atoms, without padding.
 308  *
 309  * \param[in] state  The state, can be nullptr
 310  */
 311 static inline gmx::ArrayRef<const gmx::RVec>
 312 positionsFromStatePointer(const t_state *state)
 313 {
 314     if (state)
 315     {
 316         return gmx::makeConstArrayRef(state->x).subArray(0, state->natoms);
 317     }
 318     else
 319     {
 320         return {};
 321     }
 322 };
 323
 324 /*! \brief Fills fep_state, lambda, and lam0 if needed
 325  *
 326  * If FEP or simulated tempering is in use:
 327  *
 328  *    fills non-null lam0 with the initial lambda values, and
 329  *    on master rank fills fep_state and lambda.
 330  *
 331  * Reports the initial lambda state to the log file. */
 332 void initialize_lambdas(FILE               *fplog,
 333                         const t_inputrec   &ir,
 334                         bool                isMaster,
 335                         int                *fep_state,
 336                         gmx::ArrayRef<real> lambda,
 337                         double             *lam0);
 338
 339 #endif