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, 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 |
#ifndef _state_h_
|
38 |
#define _state_h_
|
39 |
|
40 |
|
41 |
#include "simple.h" |
42 |
#include "../../swap/enums.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 |
/* These enums are used in flags as (1<<est...).
|
55 |
* The order of these enums should not be changed,
|
56 |
* since that affects the checkpoint (.cpt) file format.
|
57 |
*/
|
58 |
enum {
|
59 |
estLAMBDA, |
60 |
estBOX, estBOX_REL, estBOXV, estPRES_PREV, estNH_XI, estTC_INT, |
61 |
estX, estV, estSDX, estCGP, estLD_RNG, estLD_RNGI, |
62 |
estDISRE_INITF, estDISRE_RM3TAV, |
63 |
estORIRE_INITF, estORIRE_DTAV, |
64 |
estSVIR_PREV, estNH_VXI, estVETA, estVOL0, estNHPRES_XI, estNHPRES_VXI, estFVIR_PREV, |
65 |
estFEPSTATE, estMC_RNG, estMC_RNGI, |
66 |
estNR |
67 |
}; |
68 |
|
69 |
#define EST_DISTR(e) (!(((e) >= estLAMBDA && (e) <= estTC_INT) || ((e) >= estSVIR_PREV && (e) <= estMC_RNGI)))
|
70 |
|
71 |
/* The names of the state entries, defined in src/gmxlib/checkpoint.c */
|
72 |
extern const char *est_names[estNR]; |
73 |
|
74 |
typedef struct |
75 |
{ |
76 |
real disre_initf; /* The scaling factor for initializing the time av. */
|
77 |
int ndisrepairs; /* The number of distance restraints */ |
78 |
real *disre_rm3tav; /* The r^-3 time averaged pair distances */
|
79 |
real orire_initf; /* The scaling factor for initializing the time av. */
|
80 |
int norire_Dtav; /* The number of matrix element in dtav (npair*5) */ |
81 |
real *orire_Dtav; /* The time averaged orientation tensors */
|
82 |
} history_t; |
83 |
|
84 |
/* Struct used for checkpointing only.
|
85 |
* This struct would not be required with unlimited precision.
|
86 |
* But because of limited precision, the COM motion removal implementation
|
87 |
* can cause the kinetic energy in the MD loop to differ by a few bits from
|
88 |
* the kinetic energy one would determine from state.v.
|
89 |
*/
|
90 |
typedef struct |
91 |
{ |
92 |
gmx_bool bUpToDate; |
93 |
int ekin_n;
|
94 |
tensor *ekinh; |
95 |
tensor *ekinf; |
96 |
tensor *ekinh_old; |
97 |
tensor ekin_total; |
98 |
double *ekinscalef_nhc;
|
99 |
double *ekinscaleh_nhc;
|
100 |
double *vscale_nhc;
|
101 |
real dekindl; |
102 |
real mvcos; |
103 |
} ekinstate_t; |
104 |
|
105 |
/* energy history for delta_h histograms */
|
106 |
typedef struct |
107 |
{ |
108 |
int nndh; /* the number of energy difference lists */ |
109 |
int *ndh; /* the number in each energy difference list */ |
110 |
real **dh; /* the energy difference lists */
|
111 |
|
112 |
double start_time; /* the start time of these energy diff blocks */ |
113 |
double start_lambda; /* lambda at start time */ |
114 |
|
115 |
gmx_bool start_lambda_set; /* whether the lambda value is set. Here
|
116 |
For backward-compatibility. */
|
117 |
} delta_h_history_t; |
118 |
|
119 |
typedef struct |
120 |
{ |
121 |
int nlambda; /* total number of lambda states - for history*/ |
122 |
|
123 |
gmx_bool bEquil; /* Have we reached equilibration */
|
124 |
int *n_at_lam; /* number of points observed at each lambda */ |
125 |
real *wl_histo; /* histogram for WL flatness determination */
|
126 |
real wl_delta; /* current wang-landau delta */
|
127 |
real *dfavg; /* used to update the free energy estimate for AIM */
|
128 |
real *sum_dvdl; /* used to update the sum of the derivative of the potential energy term at lambda */
|
129 |
real *store_fepot; /* store the potential energy at lambda */
|
130 |
int *aim_at_lam; /* number of points observed at each lambda after chosen*/ |
131 |
int *laccept; /* track lambda acceptance */ |
132 |
real *sum_weights; /* weights of the states */
|
133 |
real *sum_dg; /* free energies of the states -- not actually used for weighting, but informational */
|
134 |
real *sum_minvar; /* corrections to weights for minimum variance */
|
135 |
real *sum_variance; /* variances of the states */
|
136 |
|
137 |
real **accum_p; /* accumulated bennett weights for n+1 */
|
138 |
real **accum_m; /* accumulated bennett weights for n-1 */
|
139 |
real **accum_p2; /* accumulated squared bennett weights for n+1 */
|
140 |
real **accum_m2; /* accumulated squared bennett weights for n-1 */
|
141 |
|
142 |
real **Tij; /* transition matrix */
|
143 |
real **Tij_empirical; /* Empirical transition matrix */
|
144 |
|
145 |
} df_history_t; |
146 |
|
147 |
typedef struct |
148 |
{ |
149 |
gmx_int64_t nsteps; /* The number of steps in the history */
|
150 |
gmx_int64_t nsum; /* The nr. of steps in the ener_ave and ener_sum */
|
151 |
double * ener_ave; /* Energy term history sum to get fluctuations */ |
152 |
double * ener_sum; /* Energy term history sum to get fluctuations */ |
153 |
int nener; /* Number of energy terms in two previous arrays */ |
154 |
gmx_int64_t nsteps_sim; /* The number of steps in ener_sum_sim */
|
155 |
gmx_int64_t nsum_sim; /* The number of frames in ener_sum_sim */
|
156 |
double * ener_sum_sim; /* Energy term history sum of the whole sim */ |
157 |
|
158 |
delta_h_history_t *dht; /* The BAR energy differences */
|
159 |
} |
160 |
energyhistory_t; |
161 |
|
162 |
typedef struct |
163 |
{ |
164 |
/* If one uses essential dynamics or flooding on a group of atoms from
|
165 |
* more than one molecule, we cannot make this group whole with
|
166 |
* do_pbc_first_mtop(). We assume that the ED group has the correct PBC
|
167 |
* representation at the beginning of the simulation and keep track
|
168 |
* of the shifts to always get it into that representation.
|
169 |
* For proper restarts from a checkpoint we store the positions of the
|
170 |
* reference group at the time of checkpoint writing */
|
171 |
gmx_bool bFromCpt; /* Did we start from a checkpoint file? */
|
172 |
int nED; /* No. of ED/Flooding data sets, if <1 no ED */ |
173 |
int *nref; /* No. of atoms in i'th reference structure */ |
174 |
int *nav; /* Same for average structure */ |
175 |
rvec **old_sref; /* Positions of the reference atoms
|
176 |
at the last time step (with correct PBC
|
177 |
representation) */
|
178 |
rvec **old_sref_p; /* Pointer to these positions */
|
179 |
rvec **old_sav; /* Same for the average positions */
|
180 |
rvec **old_sav_p; |
181 |
} |
182 |
edsamstate_t; |
183 |
|
184 |
|
185 |
typedef struct |
186 |
{ |
187 |
int eSwapCoords; /* Swapping along x, y, or z-direction? */ |
188 |
int nat_req[eCompNR][eIonNR]; /* Requested ion numbers per type an comp. */ |
189 |
int *nat_req_p[eCompNR][eIonNR]; /* Pointer to this data (for .cpt writing) */ |
190 |
int nAverage; /* Use average over this many swap attempt |
191 |
steps when determining the ion counts */
|
192 |
int inflow_netto[eCompNR][eIonNR]; /* Flux determined from the # of swaps */ |
193 |
int *inflow_netto_p[eCompNR][eIonNR]; /* Pointer to this data */ |
194 |
int *nat_past[eCompNR][eIonNR]; /* Array with nAverage entries for history */ |
195 |
int *nat_past_p[eCompNR][eIonNR]; /* Pointer points to the first entry only */ |
196 |
|
197 |
/* Channel flux detection, this is counting only and has no influence on whether swaps
|
198 |
* are performed or not: */
|
199 |
int fluxfromAtoB[eCompNR][eIonNR]; /* Flux determined from the split cylinders */ |
200 |
int *fluxfromAtoB_p[eCompNR][eIonNR]; /* Pointer to this data */ |
201 |
int *fluxleak; /* Flux not going through any channel */ |
202 |
int nions; /* Size of the following arrays */ |
203 |
unsigned char *comp_from; /* Ion came from which compartment? */ |
204 |
unsigned char *channel_label; /* Through which channel did this ion pass? */ |
205 |
|
206 |
/* To also make multimeric channel proteins whole, we save the last whole configuration of
|
207 |
* the channels in the checkpoint file. If we have no checkpoint file, we assume that the
|
208 |
* starting configuration hast the correct PBC representation after making the individual
|
209 |
* molecules whole */
|
210 |
gmx_bool bFromCpt; /* Did we started from a checkpoint file? */
|
211 |
int nat[eChanNR]; /* Size of xc_old_whole, i.e. the number of |
212 |
atoms in each channel */
|
213 |
rvec *xc_old_whole[eChanNR]; /* Last known whole positions of the two
|
214 |
channels (important for multimeric ch.!) */
|
215 |
rvec **xc_old_whole_p[eChanNR]; /* Pointer to these positions */
|
216 |
} |
217 |
swapstate_t; |
218 |
|
219 |
|
220 |
typedef struct |
221 |
{ |
222 |
int natoms;
|
223 |
int ngtc;
|
224 |
int nnhpres;
|
225 |
int nhchainlength; /* number of nose-hoover chains */ |
226 |
int flags; /* Flags telling which entries are present */ |
227 |
int fep_state; /* indicates which of the alchemical states we are in */ |
228 |
real *lambda; /* lambda vector */
|
229 |
matrix box; /* box vector coordinates */
|
230 |
matrix box_rel; /* Relitaive box vectors to preserve shape */
|
231 |
matrix boxv; /* box velocitites for Parrinello-Rahman pcoupl */
|
232 |
matrix pres_prev; /* Pressure of the previous step for pcoupl */
|
233 |
matrix svir_prev; /* Shake virial for previous step for pcoupl */
|
234 |
matrix fvir_prev; /* Force virial of the previous step for pcoupl */
|
235 |
double *nosehoover_xi; /* for Nose-Hoover tcoupl (ngtc) */ |
236 |
double *nosehoover_vxi; /* for N-H tcoupl (ngtc) */ |
237 |
double *nhpres_xi; /* for Nose-Hoover pcoupl for barostat */ |
238 |
double *nhpres_vxi; /* for Nose-Hoover pcoupl for barostat */ |
239 |
double *therm_integral; /* for N-H/V-rescale tcoupl (ngtc) */ |
240 |
real veta; /* trotter based isotropic P-coupling */
|
241 |
real vol0; /* initial volume,required for computing NPT conserverd quantity */
|
242 |
int nalloc; /* Allocation size for x, v and sd_x when !=NULL*/ |
243 |
rvec *x; /* the coordinates (natoms) */
|
244 |
rvec *v; /* the velocities (natoms) */
|
245 |
rvec *sd_X; /* random part of the x update for stoch. dyn. */
|
246 |
rvec *cg_p; /* p vector for conjugate gradient minimization */
|
247 |
|
248 |
history_t hist; /* Time history for restraints */
|
249 |
|
250 |
ekinstate_t ekinstate; /* The state of the kinetic energy data */
|
251 |
|
252 |
energyhistory_t enerhist; /* Energy history for statistics */
|
253 |
swapstate_t swapstate; /* Position swapping */
|
254 |
df_history_t dfhist; /*Free energy history for free energy analysis */
|
255 |
edsamstate_t edsamstate; /* Essential dynamics / flooding history */
|
256 |
|
257 |
int ddp_count; /* The DD partitioning count for this state */ |
258 |
int ddp_count_cg_gl; /* The DD part. count for index_gl */ |
259 |
int ncg_gl; /* The number of local charge groups */ |
260 |
int *cg_gl; /* The global cg number of the local cgs */ |
261 |
int cg_gl_nalloc; /* Allocation size of cg_gl; */ |
262 |
} t_state; |
263 |
|
264 |
typedef struct |
265 |
{ |
266 |
double *Qinv; /* inverse mass of thermostat -- computed from inputs, but a good place to store */ |
267 |
double *QPinv; /* inverse mass of thermostat for barostat -- computed from inputs, but a good place to store */ |
268 |
double Winv; /* Pressure mass inverse -- computed, not input, but a good place to store. Need to make a matrix later */ |
269 |
tensor Winvm; /* inverse pressure mass tensor, computed */
|
270 |
} t_extmass; |
271 |
|
272 |
|
273 |
typedef struct |
274 |
{ |
275 |
real veta; |
276 |
double rscale;
|
277 |
double vscale;
|
278 |
double rvscale;
|
279 |
double alpha;
|
280 |
double *vscale_nhc;
|
281 |
} t_vetavars; |
282 |
|
283 |
#ifdef __cplusplus
|
284 |
} |
285 |
#endif
|
286 |
|
287 |
|
288 |
#endif /* _state_h_ */ |