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, 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 <vector>
  59
  60 #include "gromacs/math/paddedvector.h"
  61 #include "gromacs/math/vectypes.h"
  62 #include "gromacs/mdtypes/md_enums.h"
  63 #include "gromacs/utility/basedefinitions.h"
  64 #include "gromacs/utility/real.h"
  65
  66 /*
  67  * The t_state struct should contain all the (possibly) non-static
  68  * information required to define the state of the system.
  69  * Currently the random seeds for SD and BD are missing.
  70  */
  71
  72 /* \brief Enum for all entries in \p t_state
  73  *
  74  * These enums are used in flags as (1<<est...).
  75  * The order of these enums should not be changed,
  76  * since that affects the checkpoint (.cpt) file format.
  77  */
  78 enum {
  79     estLAMBDA,
  80     estBOX, estBOX_REL, estBOXV, estPRES_PREV, estNH_XI,  estTHERM_INT,
  81     estX,   estV,       estSDX_NOTSUPPORTED,  estCGP,
  82     estLD_RNG_NOTSUPPORTED, estLD_RNGI_NOTSUPPORTED,
  83     estDISRE_INITF, estDISRE_RM3TAV,
  84     estORIRE_INITF, estORIRE_DTAV,
  85     estSVIR_PREV, estNH_VXI, estVETA, estVOL0, estNHPRES_XI, estNHPRES_VXI, estFVIR_PREV,
  86     estFEPSTATE, estMC_RNG_NOTSUPPORTED, estMC_RNGI_NOTSUPPORTED,
  87     estBAROS_INT,
  88     estNR
  89 };
  90
  91 //! \brief The names of the state entries, defined in src/gmxlib/checkpoint.c
  92 extern const char *est_names[estNR];
  93
  94 /*! \libinternal \brief History information for NMR distance and orientation restraints
  95  *
  96  * Often this is only used for reporting observables, and thus should not
  97  * actually be part of the microstate. But with time-dependent restraining
  98  * they are actually part of the (non-Markovian) microstate.
  99  * \todo Rename this with a more descriptive name.
 100  */
 101 class history_t
 102 {
 103     public:
 104         //! Constructor
 105         history_t();
 106
 107         real  disre_initf;  //!< The scaling factor for initializing the time av.
 108         int   ndisrepairs;  //!< The number of distance restraints
 109         real *disre_rm3tav; //!< The r^-3 time averaged pair distances
 110         real  orire_initf;  //!< The scaling factor for initializing the time av.
 111         int   norire_Dtav;  //!< The number of matrix element in dtav (npair*5)
 112         real *orire_Dtav;   //!< The time averaged orientation tensors
 113 };
 114
 115 /*! \libinternal \brief Struct used for checkpointing only
 116  *
 117  * This struct would not be required with unlimited precision.
 118  * But because of limited precision, the COM motion removal implementation
 119  * can cause the kinetic energy in the MD loop to differ by a few bits from
 120  * the kinetic energy one would determine from state.v.
 121  */
 122 class ekinstate_t
 123 {
 124     public:
 125         ekinstate_t();
 126
 127         gmx_bool             bUpToDate;      //!< Test if all data is up to date
 128         int                  ekin_n;         //!< The number of tensors
 129         tensor              *ekinh;          //!< Half step Ekin, size \p ekin_n
 130         tensor              *ekinf;          //!< Full step Ekin, size \p ekin_n
 131         tensor              *ekinh_old;      //!< Half step Ekin of the previous step, size \p ekin_n
 132         tensor               ekin_total;     //!< Total kinetic energy
 133         std::vector<double>  ekinscalef_nhc; //!< Nose-Hoover Ekin scaling factors for full step Ekin
 134         std::vector<double>  ekinscaleh_nhc; //!< Nose-Hoover Ekin scaling factors for half step Ekin
 135         std::vector<double>  vscale_nhc;     //!< Nose-Hoover velocity scaling factors
 136         real                 dekindl;        //!< dEkin/dlambda, with free-energy
 137         real                 mvcos;          //!< Cosine(z) component of the momentum, for viscosity calculations
 138 };
 139
 140 /*! \brief Free-energy sampling history struct
 141  *
 142  * \todo Split out into microstate and observables history.
 143  */
 144 typedef struct df_history_t
 145 {
 146     int      nlambda;        //!< total number of lambda states - for history
 147
 148     gmx_bool bEquil;         //!< Have we reached equilibration
 149     int     *n_at_lam;       //!< number of points observed at each lambda
 150     real    *wl_histo;       //!< histogram for WL flatness determination
 151     real     wl_delta;       //!< current wang-landau delta
 152
 153     real    *sum_weights;    //!< weights of the states
 154     real    *sum_dg;         //!< free energies of the states -- not actually used for weighting, but informational
 155     real    *sum_minvar;     //!< corrections to weights for minimum variance
 156     real    *sum_variance;   //!< variances of the states
 157
 158     real   **accum_p;        //!< accumulated bennett weights for n+1
 159     real   **accum_m;        //!< accumulated bennett weights for n-1
 160     real   **accum_p2;       //!< accumulated squared bennett weights for n+1
 161     real   **accum_m2;       //!< accumulated squared bennett weights for n-1
 162
 163     real   **Tij;            //!< transition matrix
 164     real   **Tij_empirical;  //!< Empirical transition matrix
 165
 166 } df_history_t;
 167
 168
 169 /*! \brief The microstate of the system
 170  *
 171  * The global state will contain complete data for all used entries.
 172  * The local state with domain decomposition will have partial entries
 173  * for which \p stateEntryIsAtomProperty() is true. Some entries that
 174  * are used in the global state might not be present in the local state.
 175  * \todo Move pure observables history to ObservablesHistory.
 176  */
 177 class t_state
 178 {
 179     public:
 180         //! Constructor
 181         t_state();
 182
 183         // All things public
 184         int                      natoms;          //!< Number of atoms, local + non-local; this is the size of \p x, \p v and \p cg_p, when used
 185         int                      ngtc;            //!< The number of temperature coupling groups
 186         int                      nnhpres;         //!< The NH-chain length for the MTTK barostat
 187         int                      nhchainlength;   //!< The NH-chain length for temperature coupling
 188         int                      flags;           //!< Set of bit-flags telling which entries are present, see enum at the top of the file
 189         int                      fep_state;       //!< indicates which of the alchemical states we are in
 190         std::array<real, efptNR> lambda;          //!< Free-energy lambda vector
 191         matrix                   box;             //!< Matrix of box vectors
 192         matrix                   box_rel;         //!< Relative box vectors to preserve box shape
 193         matrix                   boxv;            //!< Box velocities for Parrinello-Rahman P-coupling
 194         matrix                   pres_prev;       //!< Pressure of the previous step for pcoupl
 195         matrix                   svir_prev;       //!< Shake virial for previous step for pcoupl
 196         matrix                   fvir_prev;       //!< Force virial of the previous step for pcoupl
 197         std::vector<double>      nosehoover_xi;   //!< Nose-Hoover coordinates (ngtc)
 198         std::vector<double>      nosehoover_vxi;  //!< Nose-Hoover velocities (ngtc)
 199         std::vector<double>      nhpres_xi;       //!< Pressure Nose-Hoover coordinates
 200         std::vector<double>      nhpres_vxi;      //!< Pressure Nose-Hoover velocities
 201         std::vector<double>      therm_integral;  //!< Work exterted N-H/V-rescale T-coupling (ngtc)
 202         double                   baros_integral;  //!< For Berendsen P-coupling conserved quantity
 203         real                     veta;            //!< Trotter based isotropic P-coupling
 204         real                     vol0;            //!< Initial volume,required for computing MTTK conserved quantity
 205         PaddedRVecVector         x;               //!< The coordinates (natoms)
 206         PaddedRVecVector         v;               //!< The velocities (natoms)
 207         PaddedRVecVector         cg_p;            //!< p vector for conjugate gradient minimization
 208
 209         ekinstate_t              ekinstate;       //!< The state of the kinetic energy
 210
 211         /* History for special algorithms, should be moved to a history struct */
 212         history_t                hist;            //!< Time history for restraints
 213         df_history_t            *dfhist;          //!< Free-energy history for free energy analysis
 214
 215         int                      ddp_count;       //!< The DD partitioning count for this state
 216         int                      ddp_count_cg_gl; //!< The DD partitioning count for index_gl
 217         std::vector<int>         cg_gl;           //!< The global cg number of the local cgs
 218 };
 219
 220 #ifndef DOXYGEN
 221 /* We don't document the structs below, as they don't belong here.
 222  * TODO: Move the next two structs out of state.h.
 223  */
 224
 225 typedef struct t_extmass
 226 {
 227     double *Qinv;  /* inverse mass of thermostat -- computed from inputs, but a good place to store */
 228     double *QPinv; /* inverse mass of thermostat for barostat -- computed from inputs, but a good place to store */
 229     double  Winv;  /* Pressure mass inverse -- computed, not input, but a good place to store. Need to make a matrix later */
 230     tensor  Winvm; /* inverse pressure mass tensor, computed       */
 231 } t_extmass;
 232
 233
 234 typedef struct
 235 {
 236     real    veta;
 237     double  rscale;
 238     double  vscale;
 239     double  rvscale;
 240     double  alpha;
 241     double *vscale_nhc;
 242 } t_vetavars;
 243
 244 #endif // DOXYGEN
 245
 246 //! Resizes the T- and P-coupling state variables
 247 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength);
 248
 249 //! Change the number of atoms represented by this state, allocating memory as needed.
 250 void state_change_natoms(t_state *state, int natoms);
 251
 252 //! Allocates memory for free-energy history
 253 void init_dfhist_state(t_state *state, int dfhistNumLambda);
 254
 255 /*! \brief Compares two states, write the differences to stdout */
 256 void comp_state(const t_state *st1, const t_state *st2, gmx_bool bRMSD, real ftol, real abstol);
 257
 258 /*! \brief Allocates an rvec pointer and copy the contents of v to it */
 259 rvec *getRvecArrayFromPaddedRVecVector(const PaddedRVecVector *v,
 260                                        unsigned int            n);
 261
 262 #endif