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