2 * This file is part of the GROMACS molecular simulation package.
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, 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.
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.
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.
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.
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.
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.
37 #ifndef GMX_MDTYPES_STATE_H
38 #define GMX_MDTYPES_STATE_H
43 #include "gromacs/math/vectypes.h"
44 #include "gromacs/mdtypes/md_enums.h"
45 #include "gromacs/utility/basedefinitions.h"
46 #include "gromacs/utility/real.h"
48 struct energyhistory_t
;
51 * The t_state struct should contain all the (possibly) non-static
52 * information required to define the state of the system.
53 * Currently the random seeds for SD and BD are missing.
56 /* These enums are used in flags as (1<<est...).
57 * The order of these enums should not be changed,
58 * since that affects the checkpoint (.cpt) file format.
62 estBOX
, estBOX_REL
, estBOXV
, estPRES_PREV
, estNH_XI
, estTC_INT
,
63 estX
, estV
, est_SDX_NOTSUPPORTED
, estCGP
,
64 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
,
72 #define EST_DISTR(e) (!(((e) >= estLAMBDA && (e) <= estTC_INT) || ((e) >= estSVIR_PREV && (e) <= estMC_RNGI)))
74 /* The names of the state entries, defined in src/gmxlib/checkpoint.c */
75 extern const char *est_names
[estNR
];
77 /* This vector is not padded yet, padding will be added soon */
78 typedef std::vector
<gmx::RVec
> PaddedRVecVector
;
80 typedef std::shared_ptr
<energyhistory_t
> EnerHistPtr
;
82 typedef struct history_t
84 real disre_initf
; /* The scaling factor for initializing the time av. */
85 int ndisrepairs
; /* The number of distance restraints */
86 real
*disre_rm3tav
; /* The r^-3 time averaged pair distances */
87 real orire_initf
; /* The scaling factor for initializing the time av. */
88 int norire_Dtav
; /* The number of matrix element in dtav (npair*5) */
89 real
*orire_Dtav
; /* The time averaged orientation tensors */
92 /* Struct used for checkpointing only.
93 * This struct would not be required with unlimited precision.
94 * But because of limited precision, the COM motion removal implementation
95 * can cause the kinetic energy in the MD loop to differ by a few bits from
96 * the kinetic energy one would determine from state.v.
98 typedef struct ekinstate_t
106 std::vector
<double> ekinscalef_nhc
;
107 std::vector
<double> ekinscaleh_nhc
;
108 std::vector
<double> vscale_nhc
;
113 typedef struct df_history_t
115 int nlambda
; /* total number of lambda states - for history*/
117 gmx_bool bEquil
; /* Have we reached equilibration */
118 int *n_at_lam
; /* number of points observed at each lambda */
119 real
*wl_histo
; /* histogram for WL flatness determination */
120 real wl_delta
; /* current wang-landau delta */
122 real
*sum_weights
; /* weights of the states */
123 real
*sum_dg
; /* free energies of the states -- not actually used for weighting, but informational */
124 real
*sum_minvar
; /* corrections to weights for minimum variance */
125 real
*sum_variance
; /* variances of the states */
127 real
**accum_p
; /* accumulated bennett weights for n+1 */
128 real
**accum_m
; /* accumulated bennett weights for n-1 */
129 real
**accum_p2
; /* accumulated squared bennett weights for n+1 */
130 real
**accum_m2
; /* accumulated squared bennett weights for n-1 */
132 real
**Tij
; /* transition matrix */
133 real
**Tij_empirical
; /* Empirical transition matrix */
137 typedef struct edsamstate_t
139 /* If one uses essential dynamics or flooding on a group of atoms from
140 * more than one molecule, we cannot make this group whole with
141 * do_pbc_first_mtop(). We assume that the ED group has the correct PBC
142 * representation at the beginning of the simulation and keep track
143 * of the shifts to always get it into that representation.
144 * For proper restarts from a checkpoint we store the positions of the
145 * reference group at the time of checkpoint writing */
146 gmx_bool bFromCpt
; /* Did we start from a checkpoint file? */
147 int nED
; /* No. of ED/Flooding data sets, if <1 no ED */
148 int *nref
; /* No. of atoms in i'th reference structure */
149 int *nav
; /* Same for average structure */
150 rvec
**old_sref
; /* Positions of the reference atoms
151 at the last time step (with correct PBC
153 rvec
**old_sref_p
; /* Pointer to these positions */
154 rvec
**old_sav
; /* Same for the average positions */
159 typedef struct swapstateIons_t
161 int nMolReq
[eCompNR
]; /* Requested # of molecules per compartment */
162 int *nMolReq_p
[eCompNR
]; /* Pointer to this data (for .cpt writing) */
163 int inflow_net
[eCompNR
]; /* Flux determined from the # of swaps */
164 int *inflow_net_p
[eCompNR
]; /* Pointer to this data */
165 int *nMolPast
[eCompNR
]; /* Array with nAverage entries for history */
166 int *nMolPast_p
[eCompNR
]; /* Pointer points to the first entry only */
168 /* Channel flux detection, this is counting only and has no influence on whether swaps */
169 /* are performed or not: */
170 int fluxfromAtoB
[eCompNR
]; /* Flux determined from the split cylinders */
171 int *fluxfromAtoB_p
[eCompNR
]; /* Pointer to this data */
172 int nMol
; /* # of molecules, size of the following arrays */
173 unsigned char *comp_from
; /* Ion came from which compartment? */
174 unsigned char *channel_label
; /* Through which channel did this ion pass? */
177 typedef struct swapstate_t
179 int eSwapCoords
; /* Swapping along x, y, or z-direction? */
180 int nIonTypes
; /* Number of ion types, this is the size of */
181 /* the following arrays */
182 int nAverage
; /* Use average over this many swap attempt */
183 /* steps when determining the ion counts */
184 int fluxleak
; /* Ions not going through any channel (bad!) */
185 int *fluxleak_p
; /* Pointer to this data */
187 /* To also make multimeric channel proteins whole, we save the last whole configuration */
188 /* of the channels in the checkpoint file. If we have no checkpoint file, we assume */
189 /* that the starting configuration has the correct PBC representation after making the */
190 /* individual molecules whole */
191 gmx_bool bFromCpt
; /* Did we start from a checkpoint file? */
192 int nat
[eChanNR
]; /* Size of xc_old_whole, i.e. the number of */
193 /* atoms in each channel */
194 rvec
*xc_old_whole
[eChanNR
]; /* Last known whole positions of the two */
195 /* channels (important for multimeric ch.!) */
196 rvec
**xc_old_whole_p
[eChanNR
]; /* Pointer to these positions */
197 swapstateIons_t
*ionType
;
202 typedef struct t_state
207 int nhchainlength
; /* number of nose-hoover chains */
208 int flags
; /* Flags telling which entries are present */
209 int fep_state
; /* indicates which of the alchemical states we are in */
210 std::vector
<real
> lambda
; /* lambda vector */
211 matrix box
; /* box vector coordinates */
212 matrix box_rel
; /* Relitaive box vectors to preserve shape */
213 matrix boxv
; /* box velocitites for Parrinello-Rahman pcoupl */
214 matrix pres_prev
; /* Pressure of the previous step for pcoupl */
215 matrix svir_prev
; /* Shake virial for previous step for pcoupl */
216 matrix fvir_prev
; /* Force virial of the previous step for pcoupl */
217 std::vector
<double> nosehoover_xi
; /* for Nose-Hoover tcoupl (ngtc) */
218 std::vector
<double> nosehoover_vxi
; /* for N-H tcoupl (ngtc) */
219 std::vector
<double> nhpres_xi
; /* for Nose-Hoover pcoupl for barostat */
220 std::vector
<double> nhpres_vxi
; /* for Nose-Hoover pcoupl for barostat */
221 std::vector
<double> therm_integral
; /* for N-H/V-rescale tcoupl (ngtc) */
222 real veta
; /* trotter based isotropic P-coupling */
223 real vol0
; /* initial volume,required for computing NPT conserverd quantity */
224 PaddedRVecVector x
; /* the coordinates (natoms) */
225 PaddedRVecVector v
; /* the velocities (natoms) */
226 PaddedRVecVector cg_p
; /* p vector for conjugate gradient minimization */
228 history_t hist
; /* Time history for restraints */
230 ekinstate_t ekinstate
; /* The state of the kinetic energy data */
232 EnerHistPtr enerhist
; /* Energy history for statistics */
233 swapstate_t
*swapstate
; /* Position swapping */
234 df_history_t
*dfhist
; /*Free energy history for free energy analysis */
235 edsamstate_t
*edsamstate
; /* Essential dynamics / flooding history */
237 int ddp_count
; /* The DD partitioning count for this state */
238 int ddp_count_cg_gl
; /* The DD part. count for index_gl */
239 std::vector
<int> cg_gl
; /* The global cg number of the local cgs */
242 typedef struct t_extmass
244 double *Qinv
; /* inverse mass of thermostat -- computed from inputs, but a good place to store */
245 double *QPinv
; /* inverse mass of thermostat for barostat -- computed from inputs, but a good place to store */
246 double Winv
; /* Pressure mass inverse -- computed, not input, but a good place to store. Need to make a matrix later */
247 tensor Winvm
; /* inverse pressure mass tensor, computed */
261 void init_gtc_state(t_state
*state
, int ngtc
, int nnhpres
, int nhchainlength
);
263 void init_state(t_state
*state
, int natoms
, int ngtc
, int nnhpres
, int nhchainlength
, int dfhistNumLambda
);
265 void comp_state(const t_state
*st1
, const t_state
*st2
, gmx_bool bRMSD
, real ftol
, real abstol
);
267 /*! \brief Allocate an rvec pointer and copy the contents of v to it */
268 rvec
*getRvecArrayFromPaddedRVecVector(const PaddedRVecVector
*v
,