state.h
1 
/*


2 
* This file is part of the GROMACS molecular simulation package.

3 
*

4 
* Copyright (c) 19912000, University of Groningen, The Netherlands.

5 
* Copyright (c) 20012004, 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 
* toplevel 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 021101301 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) nonstatic

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 backwardcompatibility. */

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 wanglandau 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 n1 */

139 
real **accum_p2; /* accumulated squared bennett weights for n+1 */

140 
real **accum_m2; /* accumulated squared bennett weights for n1 */

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 zdirection? */ 
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 nosehoover 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 ParrinelloRahman 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 NoseHoover tcoupl (ngtc) */ 
236 
double *nosehoover_vxi; /* for NH tcoupl (ngtc) */ 
237 
double *nhpres_xi; /* for NoseHoover pcoupl for barostat */ 
238 
double *nhpres_vxi; /* for NoseHoover pcoupl for barostat */ 
239 
double *therm_integral; /* for NH/Vrescale tcoupl (ngtc) */ 
240 
real veta; /* trotter based isotropic Pcoupling */

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_ */ 