include/types/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * 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 #ifndef _state_h_
  39 #define _state_h_
  40
  41
  42 #include "simple.h"
  43
  44 #ifdef __cplusplus
  45 extern "C" {
  46 #endif
  47
  48 /*
  49  * The t_state struct should contain all the (possibly) non-static
  50  * information required to define the state of the system.
  51  * Currently the random seeds for SD and BD are missing.
  52  */
  53
  54 /* for now, define the length of the NH chains here */
  55 #define NHCHAINLENGTH 10
  56 #define MAXLAMBDAS 1024
  57
  58 /* These enums are used in flags as (1<<est...).
  59  * The order of these enums should not be changed,
  60  * since that affects the checkpoint (.cpt) file format.
  61  */
  62   enum { estLAMBDA,
  63          estBOX, estBOX_REL, estBOXV, estPRES_PREV, estNH_XI,  estTC_INT,
  64          estX,   estV,       estSDX,  estCGP,       estLD_RNG, estLD_RNGI,
  65          estDISRE_INITF, estDISRE_RM3TAV,
  66          estORIRE_INITF, estORIRE_DTAV,
  67          estSVIR_PREV, estNH_VXI, estVETA, estVOL0, estNHPRES_XI, estNHPRES_VXI, estFVIR_PREV,
  68          estFEPSTATE, estMC_RNG, estMC_RNGI,
  69          estNR };
  70
  71 #define EST_DISTR(e) (!(((e) >= estLAMBDA && (e) <= estTC_INT) || ((e) >= estSVIR_PREV && (e) <= estMC_RNGI)))
  72
  73 /* The names of the state entries, defined in src/gmxlib/checkpoint.c */
  74 extern const char *est_names[estNR];
  75
  76 typedef struct
  77 {
  78   real disre_initf;    /* The scaling factor for initializing the time av. */
  79   int  ndisrepairs;    /* The number of distance restraints                */
  80   real *disre_rm3tav;  /* The r^-3 time averaged pair distances            */
  81   real orire_initf;    /* The scaling factor for initializing the time av. */
  82   int  norire_Dtav;    /* The number of matrix element in dtav (npair*5)   */
  83   real *orire_Dtav;    /* The time averaged orientation tensors            */
  84 } history_t;
  85
  86 /* Struct used for checkpointing only.
  87  * This struct would not be required with unlimited precision.
  88  * But because of limited precision, the COM motion removal implementation
  89  * can cause the kinetic energy in the MD loop to differ by a few bits from
  90  * the kinetic energy one would determine from state.v.
  91  */
  92 typedef struct
  93 {
  94   gmx_bool     bUpToDate;
  95   int      ekin_n;
  96   tensor  *ekinh;
  97   tensor  *ekinf;
  98   tensor  *ekinh_old;
  99   tensor   ekin_total;
 100   double  *ekinscalef_nhc;
 101   double  *ekinscaleh_nhc;
 102   double  *vscale_nhc;
 103   real     dekindl;
 104   real     mvcos;
 105 } ekinstate_t;
 106
 107 /* energy history for delta_h histograms */
 108 typedef struct
 109 {
 110     int nndh;           /* the number of energy difference lists */
 111     int  *ndh;          /* the number in each energy difference list */
 112     real **dh;          /* the energy difference lists */
 113
 114     double start_time;     /* the start time of these energy diff blocks */
 115     double start_lambda;   /* lambda at start time */
 116
 117     gmx_bool start_lambda_set; /* whether the lambda value is set. Here
 118                                   For backward-compatibility. */
 119 } delta_h_history_t;
 120
 121 typedef struct
 122 {
 123   int nlambda;               /* total number of lambda states - for history*/
 124
 125   gmx_bool bEquil;           /* reached equilibration */
 126   int  *n_at_lam;            /* number of points observed at each lambda */
 127   real *wl_histo;            /* histogram for WL flatness determination */
 128   real wl_delta;             /* current wang-landau delta */
 129
 130   real *sum_weights;         /* weights of the states */
 131   real *sum_dg;              /* free energies of the states -- not actually used for weighting, but informational */
 132   real *sum_minvar;          /* corrections to weights for minimum variance */
 133   real *sum_variance;        /* variances of the states */
 134
 135   real **accum_p;            /* accumulated bennett weights for n+1 */
 136   real **accum_m;            /* accumulated bennett weights for n-1 */
 137   real **accum_p2;           /* accumulated squared bennett weights for n+1 */
 138   real **accum_m2;           /* accumulated squared bennett weights for n-1 */
 139
 140   real **Tij;                /* transition matrix */
 141   real **Tij_empirical;      /* Empirical transition matrix */
 142 } df_history_t;
 143
 144 typedef struct
 145 {
 146   gmx_large_int_t nsteps;       /* The number of steps in the history            */
 147   gmx_large_int_t nsum;         /* The nr. of steps in the ener_ave and ener_sum */
 148   double *   ener_ave;     /* Energy term history sum to get fluctuations   */
 149   double *   ener_sum;     /* Energy term history sum to get fluctuations   */
 150   int        nener;        /* Number of energy terms in two previous arrays */
 151   gmx_large_int_t nsteps_sim;   /* The number of steps in ener_sum_sim      */
 152   gmx_large_int_t nsum_sim;     /* The number of frames in ener_sum_sim     */
 153   double *   ener_sum_sim; /* Energy term history sum of the whole sim      */
 154
 155   delta_h_history_t *dht;  /* The BAR energy differences */
 156 }
 157 energyhistory_t;
 158
 159 typedef struct
 160 {
 161     /* If one uses essential dynamics or flooding on a group of atoms from
 162      * more than one molecule, we cannot make this group whole with
 163      * do_pbc_first_mtop(). We assume that the ED group has the correct PBC
 164      * representation at the beginning of the simulation and keep track
 165      * of the shifts to always get it into that representation.
 166      * For proper restarts from a checkpoint we store the positions of the
 167      * reference group at the time of checkpoint writing */
 168     gmx_bool    bFromCpt;       /* Did we start from a checkpoint file?       */
 169     int         nED;            /* No. of ED/Flooding data sets, if <1 no ED  */
 170     int         *nref;          /* No. of atoms in i'th reference structure   */
 171     int         *nav;           /* Same for average structure                 */
 172     rvec        **old_sref;     /* Positions of the reference atoms
 173                                    at the last time step (with correct PBC
 174                                    representation)                            */
 175     rvec        **old_sref_p;   /* Pointer to these positions                 */
 176     rvec        **old_sav;      /* Same for the average positions             */
 177     rvec        **old_sav_p;
 178 }
 179 edsamstate_t;
 180
 181 typedef struct
 182 {
 183   int           natoms;
 184   int           ngtc;
 185   int           nnhpres;
 186   int           nhchainlength; /* number of nose-hoover chains               */
 187   int           nrng;
 188   int           nrngi;
 189   int           flags;  /* Flags telling which entries are present      */
 190   int           fep_state; /* indicates which of the alchemical states we are in                 */
 191   real          *lambda; /* lambda vector                               */
 192   matrix        box;    /* box vector coordinates                       */
 193   matrix        box_rel; /* Relitaive box vectors to preserve shape     */
 194   matrix        boxv;   /* box velocitites for Parrinello-Rahman pcoupl */
 195   matrix        pres_prev; /* Pressure of the previous step for pcoupl  */
 196   matrix        svir_prev; /* Shake virial for previous step for pcoupl */
 197   matrix        fvir_prev; /* Force virial of the previous step for pcoupl  */
 198   double        *nosehoover_xi;  /* for Nose-Hoover tcoupl (ngtc)       */
 199   double        *nosehoover_vxi; /* for N-H tcoupl (ngtc)               */
 200   double        *nhpres_xi;  /* for Nose-Hoover pcoupl for barostat     */
 201   double        *nhpres_vxi; /* for Nose-Hoover pcoupl for barostat     */
 202   double        *therm_integral; /* for N-H/V-rescale tcoupl (ngtc)     */
 203   real          veta; /* trotter based isotropic P-coupling             */
 204   real          vol0; /* initial volume,required for computing NPT conserverd quantity */
 205   int           nalloc; /* Allocation size for x, v and sd_x when !=NULL*/
 206   rvec          *x;     /* the coordinates (natoms)                     */
 207   rvec          *v;     /* the velocities (natoms)                      */
 208   rvec          *sd_X;  /* random part of the x update for stoch. dyn.  */
 209   rvec          *cg_p;  /* p vector for conjugate gradient minimization */
 210
 211   unsigned int  *ld_rng;  /* RNG random state                           */
 212   int           *ld_rngi; /* RNG index                                  */
 213
 214   int           nmcrng;   /* number of RNG states                       */
 215   unsigned int  *mc_rng;  /* lambda MC RNG random state                 */
 216   int           *mc_rngi; /* lambda MC RNG index                        */
 217
 218   history_t     hist;   /* Time history for restraints                  */
 219
 220   ekinstate_t   ekinstate; /* The state of the kinetic energy data      */
 221
 222   energyhistory_t  enerhist; /* Energy history for statistics           */
 223   df_history_t  dfhist; /*Free energy history for free energy analysis  */
 224   edsamstate_t  edsamstate;    /* Essential dynamics / flooding history */
 225
 226   int           ddp_count; /* The DD partitioning count for this state  */
 227   int           ddp_count_cg_gl; /* The DD part. count for index_gl     */
 228   int           ncg_gl; /* The number of local charge groups            */
 229   int           *cg_gl; /* The global cg number of the local cgs        */
 230   int           cg_gl_nalloc; /* Allocation size of cg_gl;              */
 231 } t_state;
 232
 233 typedef struct
 234 {
 235   double *Qinv;  /* inverse mass of thermostat -- computed from inputs, but a good place to store */
 236   double *QPinv; /* inverse mass of thermostat for barostat -- computed from inputs, but a good place to store */
 237   double Winv;   /* Pressure mass inverse -- computed, not input, but a good place to store. Need to make a matrix later */
 238   tensor Winvm;  /* inverse pressure mass tensor, computed       */
 239 } t_extmass;
 240
 241
 242 typedef struct
 243 {
 244   real veta;
 245   double rscale;
 246   double vscale;
 247   double rvscale;
 248   double alpha;
 249   double *vscale_nhc;
 250 } t_vetavars;
 251
 252 #ifdef __cplusplus
 253 }
 254 #endif
 255
 256
 257 #endif /* _state_h_ */