File: | gromacs/mdlib/expanded.c |
Location: | line 1079, column 17 |
Description: | Value stored to 'dm' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by |
5 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
6 | * and including many others, as listed in the AUTHORS file in the |
7 | * top-level source directory and at http://www.gromacs.org. |
8 | * |
9 | * GROMACS is free software; you can redistribute it and/or |
10 | * modify it under the terms of the GNU Lesser General Public License |
11 | * as published by the Free Software Foundation; either version 2.1 |
12 | * of the License, or (at your option) any later version. |
13 | * |
14 | * GROMACS is distributed in the hope that it will be useful, |
15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 | * Lesser General Public License for more details. |
18 | * |
19 | * You should have received a copy of the GNU Lesser General Public |
20 | * License along with GROMACS; if not, see |
21 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
22 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
23 | * |
24 | * If you want to redistribute modifications to GROMACS, please |
25 | * consider that scientific software is very special. Version |
26 | * control is crucial - bugs must be traceable. We will be happy to |
27 | * consider code for inclusion in the official distribution, but |
28 | * derived work must not be called official GROMACS. Details are found |
29 | * in the README & COPYING files - if they are missing, get the |
30 | * official version at http://www.gromacs.org. |
31 | * |
32 | * To help us fund GROMACS development, we humbly ask that you cite |
33 | * the research papers on the package. Check out http://www.gromacs.org. |
34 | */ |
35 | #ifdef HAVE_CONFIG_H1 |
36 | #include <config.h> |
37 | #endif |
38 | |
39 | #include <stdio.h> |
40 | #include <math.h> |
41 | #include "typedefs.h" |
42 | #include "gromacs/utility/smalloc.h" |
43 | #include "names.h" |
44 | #include "gromacs/fileio/confio.h" |
45 | #include "txtdump.h" |
46 | #include "pbc.h" |
47 | #include "chargegroup.h" |
48 | #include "gromacs/math/vec.h" |
49 | #include "nrnb.h" |
50 | #include "mshift.h" |
51 | #include "mdrun.h" |
52 | #include "update.h" |
53 | #include "physics.h" |
54 | #include "mdatoms.h" |
55 | #include "force.h" |
56 | #include "bondf.h" |
57 | #include "pme.h" |
58 | #include "disre.h" |
59 | #include "orires.h" |
60 | #include "network.h" |
61 | #include "calcmu.h" |
62 | #include "constr.h" |
63 | #include "gromacs/random/random.h" |
64 | #include "domdec.h" |
65 | #include "macros.h" |
66 | |
67 | #include "gromacs/fileio/confio.h" |
68 | #include "gromacs/fileio/gmxfio.h" |
69 | #include "gromacs/fileio/trnio.h" |
70 | #include "gromacs/fileio/xtcio.h" |
71 | #include "gromacs/timing/wallcycle.h" |
72 | #include "gromacs/utility/fatalerror.h" |
73 | #include "gromacs/utility/gmxmpi.h" |
74 | |
75 | static void init_df_history_weights(df_history_t *dfhist, t_expanded *expand, int nlim) |
76 | { |
77 | int i; |
78 | dfhist->wl_delta = expand->init_wl_delta; |
79 | for (i = 0; i < nlim; i++) |
80 | { |
81 | dfhist->sum_weights[i] = expand->init_lambda_weights[i]; |
82 | dfhist->sum_dg[i] = expand->init_lambda_weights[i]; |
83 | } |
84 | } |
85 | |
86 | /* Eventually should contain all the functions needed to initialize expanded ensemble |
87 | before the md loop starts */ |
88 | extern void init_expanded_ensemble(gmx_bool bStateFromCP, t_inputrec *ir, df_history_t *dfhist) |
89 | { |
90 | if (!bStateFromCP) |
91 | { |
92 | init_df_history_weights(dfhist, ir->expandedvals, ir->fepvals->n_lambda); |
93 | } |
94 | } |
95 | |
96 | static void GenerateGibbsProbabilities(real *ene, double *p_k, double *pks, int minfep, int maxfep) |
97 | { |
98 | |
99 | int i; |
100 | real maxene; |
101 | |
102 | *pks = 0.0; |
103 | maxene = ene[minfep]; |
104 | /* find the maximum value */ |
105 | for (i = minfep; i <= maxfep; i++) |
106 | { |
107 | if (ene[i] > maxene) |
108 | { |
109 | maxene = ene[i]; |
110 | } |
111 | } |
112 | /* find the denominator */ |
113 | for (i = minfep; i <= maxfep; i++) |
114 | { |
115 | *pks += exp(ene[i]-maxene); |
116 | } |
117 | /*numerators*/ |
118 | for (i = minfep; i <= maxfep; i++) |
119 | { |
120 | p_k[i] = exp(ene[i]-maxene) / *pks; |
121 | } |
122 | } |
123 | |
124 | static void GenerateWeightedGibbsProbabilities(real *ene, double *p_k, double *pks, int nlim, real *nvals, real delta) |
125 | { |
126 | |
127 | int i; |
128 | real maxene; |
129 | real *nene; |
130 | *pks = 0.0; |
131 | |
132 | snew(nene, nlim)(nene) = save_calloc("nene", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 132, (nlim), sizeof(*(nene))); |
133 | for (i = 0; i < nlim; i++) |
134 | { |
135 | if (nvals[i] == 0) |
136 | { |
137 | /* add the delta, since we need to make sure it's greater than zero, and |
138 | we need a non-arbitrary number? */ |
139 | nene[i] = ene[i] + log(nvals[i]+delta); |
140 | } |
141 | else |
142 | { |
143 | nene[i] = ene[i] + log(nvals[i]); |
144 | } |
145 | } |
146 | |
147 | /* find the maximum value */ |
148 | maxene = nene[0]; |
149 | for (i = 0; i < nlim; i++) |
150 | { |
151 | if (nene[i] > maxene) |
152 | { |
153 | maxene = nene[i]; |
154 | } |
155 | } |
156 | |
157 | /* subtract off the maximum, avoiding overflow */ |
158 | for (i = 0; i < nlim; i++) |
159 | { |
160 | nene[i] -= maxene; |
161 | } |
162 | |
163 | /* find the denominator */ |
164 | for (i = 0; i < nlim; i++) |
165 | { |
166 | *pks += exp(nene[i]); |
167 | } |
168 | |
169 | /*numerators*/ |
170 | for (i = 0; i < nlim; i++) |
171 | { |
172 | p_k[i] = exp(nene[i]) / *pks; |
173 | } |
174 | sfree(nene)save_free("nene", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 174, (nene)); |
175 | } |
176 | |
177 | real do_logsum(int N, real *a_n) |
178 | { |
179 | |
180 | /* RETURN VALUE */ |
181 | /* log(\sum_{i=0}^(N-1) exp[a_n]) */ |
182 | real maxarg; |
183 | real sum; |
184 | int i; |
185 | real logsum; |
186 | /* compute maximum argument to exp(.) */ |
187 | |
188 | maxarg = a_n[0]; |
189 | for (i = 1; i < N; i++) |
190 | { |
191 | maxarg = max(maxarg, a_n[i])(((maxarg) > (a_n[i])) ? (maxarg) : (a_n[i]) ); |
192 | } |
193 | |
194 | /* compute sum of exp(a_n - maxarg) */ |
195 | sum = 0.0; |
196 | for (i = 0; i < N; i++) |
197 | { |
198 | sum = sum + exp(a_n[i] - maxarg); |
199 | } |
200 | |
201 | /* compute log sum */ |
202 | logsum = log(sum) + maxarg; |
203 | return logsum; |
204 | } |
205 | |
206 | int FindMinimum(real *min_metric, int N) |
207 | { |
208 | |
209 | real min_val; |
210 | int min_nval, nval; |
211 | |
212 | min_nval = 0; |
213 | min_val = min_metric[0]; |
214 | |
215 | for (nval = 0; nval < N; nval++) |
216 | { |
217 | if (min_metric[nval] < min_val) |
218 | { |
219 | min_val = min_metric[nval]; |
220 | min_nval = nval; |
221 | } |
222 | } |
223 | return min_nval; |
224 | } |
225 | |
226 | static gmx_bool CheckHistogramRatios(int nhisto, real *histo, real ratio) |
227 | { |
228 | |
229 | int i; |
230 | real nmean; |
231 | gmx_bool bIfFlat; |
232 | |
233 | nmean = 0; |
234 | for (i = 0; i < nhisto; i++) |
235 | { |
236 | nmean += histo[i]; |
237 | } |
238 | |
239 | if (nmean == 0) |
240 | { |
241 | /* no samples! is bad!*/ |
242 | bIfFlat = FALSE0; |
243 | return bIfFlat; |
244 | } |
245 | nmean /= (real)nhisto; |
246 | |
247 | bIfFlat = TRUE1; |
248 | for (i = 0; i < nhisto; i++) |
249 | { |
250 | /* make sure that all points are in the ratio < x < 1/ratio range */ |
251 | if (!((histo[i]/nmean < 1.0/ratio) && (histo[i]/nmean > ratio))) |
252 | { |
253 | bIfFlat = FALSE0; |
254 | break; |
255 | } |
256 | } |
257 | return bIfFlat; |
258 | } |
259 | |
260 | static gmx_bool CheckIfDoneEquilibrating(int nlim, t_expanded *expand, df_history_t *dfhist, gmx_int64_t step) |
261 | { |
262 | |
263 | int i, totalsamples; |
264 | gmx_bool bDoneEquilibrating = TRUE1; |
265 | gmx_bool bIfFlat; |
266 | |
267 | /* assume we have equilibrated the weights, then check to see if any of the conditions are not met */ |
268 | |
269 | /* calculate the total number of samples */ |
270 | switch (expand->elmceq) |
271 | { |
272 | case elmceqNO: |
273 | /* We have not equilibrated, and won't, ever. */ |
274 | return FALSE0; |
275 | case elmceqYES: |
276 | /* we have equilibrated -- we're done */ |
277 | return TRUE1; |
278 | case elmceqSTEPS: |
279 | /* first, check if we are equilibrating by steps, if we're still under */ |
280 | if (step < expand->equil_steps) |
281 | { |
282 | bDoneEquilibrating = FALSE0; |
283 | } |
284 | break; |
285 | case elmceqSAMPLES: |
286 | totalsamples = 0; |
287 | for (i = 0; i < nlim; i++) |
288 | { |
289 | totalsamples += dfhist->n_at_lam[i]; |
290 | } |
291 | if (totalsamples < expand->equil_samples) |
292 | { |
293 | bDoneEquilibrating = FALSE0; |
294 | } |
295 | break; |
296 | case elmceqNUMATLAM: |
297 | for (i = 0; i < nlim; i++) |
298 | { |
299 | if (dfhist->n_at_lam[i] < expand->equil_n_at_lam) /* we are still doing the initial sweep, so we're definitely not |
300 | done equilibrating*/ |
301 | { |
302 | bDoneEquilibrating = FALSE0; |
303 | break; |
304 | } |
305 | } |
306 | break; |
307 | case elmceqWLDELTA: |
308 | if (EWL(expand->elamstats)((expand->elamstats) == elamstatsWL || (expand->elamstats ) == elamstatsWWL)) /* This check is in readir as well, but |
309 | just to be sure */ |
310 | { |
311 | if (dfhist->wl_delta > expand->equil_wl_delta) |
312 | { |
313 | bDoneEquilibrating = FALSE0; |
314 | } |
315 | } |
316 | break; |
317 | case elmceqRATIO: |
318 | /* we can use the flatness as a judge of good weights, as long as |
319 | we're not doing minvar, or Wang-Landau. |
320 | But turn off for now until we figure out exactly how we do this. |
321 | */ |
322 | |
323 | if (!(EWL(expand->elamstats)((expand->elamstats) == elamstatsWL || (expand->elamstats ) == elamstatsWWL) || expand->elamstats == elamstatsMINVAR)) |
324 | { |
325 | /* we want to use flatness -avoiding- the forced-through samples. Plus, we need to convert to |
326 | floats for this histogram function. */ |
327 | |
328 | real *modhisto; |
329 | snew(modhisto, nlim)(modhisto) = save_calloc("modhisto", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 329, (nlim), sizeof(*(modhisto))); |
330 | for (i = 0; i < nlim; i++) |
331 | { |
332 | modhisto[i] = 1.0*(dfhist->n_at_lam[i]-expand->lmc_forced_nstart); |
333 | } |
334 | bIfFlat = CheckHistogramRatios(nlim, modhisto, expand->equil_ratio); |
335 | sfree(modhisto)save_free("modhisto", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 335, (modhisto)); |
336 | if (!bIfFlat) |
337 | { |
338 | bDoneEquilibrating = FALSE0; |
339 | } |
340 | } |
341 | default: |
342 | bDoneEquilibrating = TRUE1; |
343 | } |
344 | /* one last case to go though, if we are doing slow growth to get initial values, we haven't finished equilibrating */ |
345 | |
346 | if (expand->lmc_forced_nstart > 0) |
347 | { |
348 | for (i = 0; i < nlim; i++) |
349 | { |
350 | if (dfhist->n_at_lam[i] < expand->lmc_forced_nstart) /* we are still doing the initial sweep, so we're definitely not |
351 | done equilibrating*/ |
352 | { |
353 | bDoneEquilibrating = FALSE0; |
354 | break; |
355 | } |
356 | } |
357 | } |
358 | return bDoneEquilibrating; |
359 | } |
360 | |
361 | static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist, |
362 | int fep_state, real *scaled_lamee, real *weighted_lamee, gmx_int64_t step) |
363 | { |
364 | real maxdiff = 0.000000001; |
365 | gmx_bool bSufficientSamples; |
366 | int i, k, n, nz, indexi, indexk, min_n, max_n, totali; |
367 | int n0, np1, nm1, nval, min_nvalm, min_nvalp, maxc; |
368 | real omega_m1_0, omega_p1_m1, omega_m1_p1, omega_p1_0, clam_osum; |
369 | real de, de_function, dr, denom, maxdr; |
370 | real min_val, cnval, zero_sum_weights; |
371 | real *omegam_array, *weightsm_array, *omegap_array, *weightsp_array, *varm_array, *varp_array, *dwp_array, *dwm_array; |
372 | real clam_varm, clam_varp, clam_weightsm, clam_weightsp, clam_minvar; |
373 | real *lam_weights, *lam_minvar_corr, *lam_variance, *lam_dg; |
374 | double *p_k; |
375 | double pks = 0; |
376 | real *numweighted_lamee, *logfrac; |
377 | int *nonzero; |
378 | real chi_m1_0, chi_p1_0, chi_m2_0, chi_p2_0, chi_p1_m1, chi_p2_m1, chi_m1_p1, chi_m2_p1; |
379 | |
380 | /* if we have equilibrated the weights, exit now */ |
381 | if (dfhist->bEquil) |
382 | { |
383 | return FALSE0; |
384 | } |
385 | |
386 | if (CheckIfDoneEquilibrating(nlim, expand, dfhist, step)) |
387 | { |
388 | dfhist->bEquil = TRUE1; |
389 | /* zero out the visited states so we know how many equilibrated states we have |
390 | from here on out.*/ |
391 | for (i = 0; i < nlim; i++) |
392 | { |
393 | dfhist->n_at_lam[i] = 0; |
394 | } |
395 | return TRUE1; |
396 | } |
397 | |
398 | /* If we reached this far, we have not equilibrated yet, keep on |
399 | going resetting the weights */ |
400 | |
401 | if (EWL(expand->elamstats)((expand->elamstats) == elamstatsWL || (expand->elamstats ) == elamstatsWWL)) |
402 | { |
403 | if (expand->elamstats == elamstatsWL) /* Standard Wang-Landau */ |
404 | { |
405 | dfhist->sum_weights[fep_state] -= dfhist->wl_delta; |
406 | dfhist->wl_histo[fep_state] += 1.0; |
407 | } |
408 | else if (expand->elamstats == elamstatsWWL) /* Weighted Wang-Landau */ |
409 | { |
410 | snew(p_k, nlim)(p_k) = save_calloc("p_k", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 410, (nlim), sizeof(*(p_k))); |
411 | |
412 | /* first increment count */ |
413 | GenerateGibbsProbabilities(weighted_lamee, p_k, &pks, 0, nlim-1); |
414 | for (i = 0; i < nlim; i++) |
415 | { |
416 | dfhist->wl_histo[i] += (real)p_k[i]; |
417 | } |
418 | |
419 | /* then increment weights (uses count) */ |
420 | pks = 0.0; |
421 | GenerateWeightedGibbsProbabilities(weighted_lamee, p_k, &pks, nlim, dfhist->wl_histo, dfhist->wl_delta); |
422 | |
423 | for (i = 0; i < nlim; i++) |
424 | { |
425 | dfhist->sum_weights[i] -= dfhist->wl_delta*(real)p_k[i]; |
426 | } |
427 | /* Alternate definition, using logarithms. Shouldn't make very much difference! */ |
428 | /* |
429 | real di; |
430 | for (i=0;i<nlim;i++) |
431 | { |
432 | di = (real)1.0 + dfhist->wl_delta*(real)p_k[i]; |
433 | dfhist->sum_weights[i] -= log(di); |
434 | } |
435 | */ |
436 | sfree(p_k)save_free("p_k", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 436, (p_k)); |
437 | } |
438 | |
439 | zero_sum_weights = dfhist->sum_weights[0]; |
440 | for (i = 0; i < nlim; i++) |
441 | { |
442 | dfhist->sum_weights[i] -= zero_sum_weights; |
443 | } |
444 | } |
445 | |
446 | if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMETROPOLIS || expand->elamstats == elamstatsMINVAR) |
447 | { |
448 | |
449 | de_function = 0; /* to get rid of warnings, but this value will not be used because of the logic */ |
450 | maxc = 2*expand->c_range+1; |
451 | |
452 | snew(lam_dg, nlim)(lam_dg) = save_calloc("lam_dg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 452, (nlim), sizeof(*(lam_dg))); |
453 | snew(lam_variance, nlim)(lam_variance) = save_calloc("lam_variance", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 453, (nlim), sizeof(*(lam_variance))); |
454 | |
455 | snew(omegap_array, maxc)(omegap_array) = save_calloc("omegap_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 455, (maxc), sizeof(*(omegap_array))); |
456 | snew(weightsp_array, maxc)(weightsp_array) = save_calloc("weightsp_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 456, (maxc), sizeof(*(weightsp_array))); |
457 | snew(varp_array, maxc)(varp_array) = save_calloc("varp_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 457, (maxc), sizeof(*(varp_array))); |
458 | snew(dwp_array, maxc)(dwp_array) = save_calloc("dwp_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 458, (maxc), sizeof(*(dwp_array))); |
459 | |
460 | snew(omegam_array, maxc)(omegam_array) = save_calloc("omegam_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 460, (maxc), sizeof(*(omegam_array))); |
461 | snew(weightsm_array, maxc)(weightsm_array) = save_calloc("weightsm_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 461, (maxc), sizeof(*(weightsm_array))); |
462 | snew(varm_array, maxc)(varm_array) = save_calloc("varm_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 462, (maxc), sizeof(*(varm_array))); |
463 | snew(dwm_array, maxc)(dwm_array) = save_calloc("dwm_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 463, (maxc), sizeof(*(dwm_array))); |
464 | |
465 | /* unpack the current lambdas -- we will only update 2 of these */ |
466 | |
467 | for (i = 0; i < nlim-1; i++) |
468 | { /* only through the second to last */ |
469 | lam_dg[i] = dfhist->sum_dg[i+1] - dfhist->sum_dg[i]; |
470 | lam_variance[i] = pow(dfhist->sum_variance[i+1], 2) - pow(dfhist->sum_variance[i], 2); |
471 | } |
472 | |
473 | /* accumulate running averages */ |
474 | for (nval = 0; nval < maxc; nval++) |
475 | { |
476 | /* constants for later use */ |
477 | cnval = (real)(nval-expand->c_range); |
478 | /* actually, should be able to rewrite it w/o exponential, for better numerical stability */ |
479 | if (fep_state > 0) |
480 | { |
481 | de = exp(cnval - (scaled_lamee[fep_state]-scaled_lamee[fep_state-1])); |
482 | if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR) |
483 | { |
484 | de_function = 1.0/(1.0+de); |
485 | } |
486 | else if (expand->elamstats == elamstatsMETROPOLIS) |
487 | { |
488 | if (de < 1.0) |
489 | { |
490 | de_function = 1.0; |
491 | } |
492 | else |
493 | { |
494 | de_function = 1.0/de; |
495 | } |
496 | } |
497 | dfhist->accum_m[fep_state][nval] += de_function; |
498 | dfhist->accum_m2[fep_state][nval] += de_function*de_function; |
499 | } |
500 | |
501 | if (fep_state < nlim-1) |
502 | { |
503 | de = exp(-cnval + (scaled_lamee[fep_state+1]-scaled_lamee[fep_state])); |
504 | if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR) |
505 | { |
506 | de_function = 1.0/(1.0+de); |
507 | } |
508 | else if (expand->elamstats == elamstatsMETROPOLIS) |
509 | { |
510 | if (de < 1.0) |
511 | { |
512 | de_function = 1.0; |
513 | } |
514 | else |
515 | { |
516 | de_function = 1.0/de; |
517 | } |
518 | } |
519 | dfhist->accum_p[fep_state][nval] += de_function; |
520 | dfhist->accum_p2[fep_state][nval] += de_function*de_function; |
521 | } |
522 | |
523 | /* Metropolis transition and Barker transition (unoptimized Bennett) acceptance weight determination */ |
524 | |
525 | n0 = dfhist->n_at_lam[fep_state]; |
526 | if (fep_state > 0) |
527 | { |
528 | nm1 = dfhist->n_at_lam[fep_state-1]; |
529 | } |
530 | else |
531 | { |
532 | nm1 = 0; |
533 | } |
534 | if (fep_state < nlim-1) |
535 | { |
536 | np1 = dfhist->n_at_lam[fep_state+1]; |
537 | } |
538 | else |
539 | { |
540 | np1 = 0; |
541 | } |
542 | |
543 | /* logic SHOULD keep these all set correctly whatever the logic, but apparently it can't figure it out. */ |
544 | chi_m1_0 = chi_p1_0 = chi_m2_0 = chi_p2_0 = chi_p1_m1 = chi_p2_m1 = chi_m1_p1 = chi_m2_p1 = 0; |
545 | |
546 | if (n0 > 0) |
547 | { |
548 | chi_m1_0 = dfhist->accum_m[fep_state][nval]/n0; |
549 | chi_p1_0 = dfhist->accum_p[fep_state][nval]/n0; |
550 | chi_m2_0 = dfhist->accum_m2[fep_state][nval]/n0; |
551 | chi_p2_0 = dfhist->accum_p2[fep_state][nval]/n0; |
552 | } |
553 | |
554 | if ((fep_state > 0 ) && (nm1 > 0)) |
555 | { |
556 | chi_p1_m1 = dfhist->accum_p[fep_state-1][nval]/nm1; |
557 | chi_p2_m1 = dfhist->accum_p2[fep_state-1][nval]/nm1; |
558 | } |
559 | |
560 | if ((fep_state < nlim-1) && (np1 > 0)) |
561 | { |
562 | chi_m1_p1 = dfhist->accum_m[fep_state+1][nval]/np1; |
563 | chi_m2_p1 = dfhist->accum_m2[fep_state+1][nval]/np1; |
564 | } |
565 | |
566 | omega_m1_0 = 0; |
567 | omega_p1_0 = 0; |
568 | clam_weightsm = 0; |
569 | clam_weightsp = 0; |
570 | clam_varm = 0; |
571 | clam_varp = 0; |
572 | |
573 | if (fep_state > 0) |
574 | { |
575 | if (n0 > 0) |
576 | { |
577 | omega_m1_0 = chi_m2_0/(chi_m1_0*chi_m1_0) - 1.0; |
578 | } |
579 | if (nm1 > 0) |
580 | { |
581 | omega_p1_m1 = chi_p2_m1/(chi_p1_m1*chi_p1_m1) - 1.0; |
582 | } |
583 | if ((n0 > 0) && (nm1 > 0)) |
584 | { |
585 | clam_weightsm = (log(chi_m1_0) - log(chi_p1_m1)) + cnval; |
586 | clam_varm = (1.0/n0)*(omega_m1_0) + (1.0/nm1)*(omega_p1_m1); |
587 | } |
588 | } |
589 | |
590 | if (fep_state < nlim-1) |
591 | { |
592 | if (n0 > 0) |
593 | { |
594 | omega_p1_0 = chi_p2_0/(chi_p1_0*chi_p1_0) - 1.0; |
595 | } |
596 | if (np1 > 0) |
597 | { |
598 | omega_m1_p1 = chi_m2_p1/(chi_m1_p1*chi_m1_p1) - 1.0; |
599 | } |
600 | if ((n0 > 0) && (np1 > 0)) |
601 | { |
602 | clam_weightsp = (log(chi_m1_p1) - log(chi_p1_0)) + cnval; |
603 | clam_varp = (1.0/np1)*(omega_m1_p1) + (1.0/n0)*(omega_p1_0); |
604 | } |
605 | } |
606 | |
607 | if (n0 > 0) |
608 | { |
609 | omegam_array[nval] = omega_m1_0; |
610 | } |
611 | else |
612 | { |
613 | omegam_array[nval] = 0; |
614 | } |
615 | weightsm_array[nval] = clam_weightsm; |
616 | varm_array[nval] = clam_varm; |
617 | if (nm1 > 0) |
618 | { |
619 | dwm_array[nval] = fabs( (cnval + log((1.0*n0)/nm1)) - lam_dg[fep_state-1] ); |
620 | } |
621 | else |
622 | { |
623 | dwm_array[nval] = fabs( cnval - lam_dg[fep_state-1] ); |
624 | } |
625 | |
626 | if (n0 > 0) |
627 | { |
628 | omegap_array[nval] = omega_p1_0; |
629 | } |
630 | else |
631 | { |
632 | omegap_array[nval] = 0; |
633 | } |
634 | weightsp_array[nval] = clam_weightsp; |
635 | varp_array[nval] = clam_varp; |
636 | if ((np1 > 0) && (n0 > 0)) |
637 | { |
638 | dwp_array[nval] = fabs( (cnval + log((1.0*np1)/n0)) - lam_dg[fep_state] ); |
639 | } |
640 | else |
641 | { |
642 | dwp_array[nval] = fabs( cnval - lam_dg[fep_state] ); |
643 | } |
644 | |
645 | } |
646 | |
647 | /* find the C's closest to the old weights value */ |
648 | |
649 | min_nvalm = FindMinimum(dwm_array, maxc); |
650 | omega_m1_0 = omegam_array[min_nvalm]; |
651 | clam_weightsm = weightsm_array[min_nvalm]; |
652 | clam_varm = varm_array[min_nvalm]; |
653 | |
654 | min_nvalp = FindMinimum(dwp_array, maxc); |
655 | omega_p1_0 = omegap_array[min_nvalp]; |
656 | clam_weightsp = weightsp_array[min_nvalp]; |
657 | clam_varp = varp_array[min_nvalp]; |
658 | |
659 | clam_osum = omega_m1_0 + omega_p1_0; |
660 | clam_minvar = 0; |
661 | if (clam_osum > 0) |
662 | { |
663 | clam_minvar = 0.5*log(clam_osum); |
664 | } |
665 | |
666 | if (fep_state > 0) |
667 | { |
668 | lam_dg[fep_state-1] = clam_weightsm; |
669 | lam_variance[fep_state-1] = clam_varm; |
670 | } |
671 | |
672 | if (fep_state < nlim-1) |
673 | { |
674 | lam_dg[fep_state] = clam_weightsp; |
675 | lam_variance[fep_state] = clam_varp; |
676 | } |
677 | |
678 | if (expand->elamstats == elamstatsMINVAR) |
679 | { |
680 | bSufficientSamples = TRUE1; |
681 | /* make sure they are all past a threshold */ |
682 | for (i = 0; i < nlim; i++) |
683 | { |
684 | if (dfhist->n_at_lam[i] < expand->minvarmin) |
685 | { |
686 | bSufficientSamples = FALSE0; |
687 | } |
688 | } |
689 | if (bSufficientSamples) |
690 | { |
691 | dfhist->sum_minvar[fep_state] = clam_minvar; |
692 | if (fep_state == 0) |
693 | { |
694 | for (i = 0; i < nlim; i++) |
695 | { |
696 | dfhist->sum_minvar[i] += (expand->minvar_const-clam_minvar); |
697 | } |
698 | expand->minvar_const = clam_minvar; |
699 | dfhist->sum_minvar[fep_state] = 0.0; |
700 | } |
701 | else |
702 | { |
703 | dfhist->sum_minvar[fep_state] -= expand->minvar_const; |
704 | } |
705 | } |
706 | } |
707 | |
708 | /* we need to rezero minvar now, since it could change at fep_state = 0 */ |
709 | dfhist->sum_dg[0] = 0.0; |
710 | dfhist->sum_variance[0] = 0.0; |
711 | dfhist->sum_weights[0] = dfhist->sum_dg[0] + dfhist->sum_minvar[0]; /* should be zero */ |
712 | |
713 | for (i = 1; i < nlim; i++) |
714 | { |
715 | dfhist->sum_dg[i] = lam_dg[i-1] + dfhist->sum_dg[i-1]; |
716 | dfhist->sum_variance[i] = sqrt(lam_variance[i-1] + pow(dfhist->sum_variance[i-1], 2)); |
717 | dfhist->sum_weights[i] = dfhist->sum_dg[i] + dfhist->sum_minvar[i]; |
718 | } |
719 | |
720 | sfree(lam_dg)save_free("lam_dg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 720, (lam_dg)); |
721 | sfree(lam_variance)save_free("lam_variance", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 721, (lam_variance)); |
722 | |
723 | sfree(omegam_array)save_free("omegam_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 723, (omegam_array)); |
724 | sfree(weightsm_array)save_free("weightsm_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 724, (weightsm_array)); |
725 | sfree(varm_array)save_free("varm_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 725, (varm_array)); |
726 | sfree(dwm_array)save_free("dwm_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 726, (dwm_array)); |
727 | |
728 | sfree(omegap_array)save_free("omegap_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 728, (omegap_array)); |
729 | sfree(weightsp_array)save_free("weightsp_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 729, (weightsp_array)); |
730 | sfree(varp_array)save_free("varp_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 730, (varp_array)); |
731 | sfree(dwp_array)save_free("dwp_array", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 731, (dwp_array)); |
732 | } |
733 | return FALSE0; |
734 | } |
735 | |
736 | static int ChooseNewLambda(int nlim, t_expanded *expand, df_history_t *dfhist, int fep_state, real *weighted_lamee, double *p_k, |
737 | gmx_int64_t seed, gmx_int64_t step) |
738 | { |
739 | /* Choose new lambda value, and update transition matrix */ |
740 | |
741 | int i, ifep, jfep, minfep, maxfep, lamnew, lamtrial, starting_fep_state; |
742 | real r1, r2, de_old, de_new, de, trialprob, tprob = 0; |
743 | real **Tij; |
744 | double *propose, *accept, *remainder; |
745 | double pks; |
746 | real sum, pnorm; |
747 | gmx_bool bRestricted; |
748 | |
749 | starting_fep_state = fep_state; |
750 | lamnew = fep_state; /* so that there is a default setting -- stays the same */ |
751 | |
752 | if (!EWL(expand->elamstats)((expand->elamstats) == elamstatsWL || (expand->elamstats ) == elamstatsWWL)) /* ignore equilibrating the weights if using WL */ |
753 | { |
754 | if ((expand->lmc_forced_nstart > 0) && (dfhist->n_at_lam[nlim-1] <= expand->lmc_forced_nstart)) |
755 | { |
756 | /* Use a marching method to run through the lambdas and get preliminary free energy data, |
757 | before starting 'free' sampling. We start free sampling when we have enough at each lambda */ |
758 | |
759 | /* if we have enough at this lambda, move on to the next one */ |
760 | |
761 | if (dfhist->n_at_lam[fep_state] == expand->lmc_forced_nstart) |
762 | { |
763 | lamnew = fep_state+1; |
764 | if (lamnew == nlim) /* whoops, stepped too far! */ |
765 | { |
766 | lamnew -= 1; |
767 | } |
768 | } |
769 | else |
770 | { |
771 | lamnew = fep_state; |
772 | } |
773 | return lamnew; |
774 | } |
775 | } |
776 | |
777 | snew(propose, nlim)(propose) = save_calloc("propose", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 777, (nlim), sizeof(*(propose))); |
778 | snew(accept, nlim)(accept) = save_calloc("accept", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 778, (nlim), sizeof(*(accept))); |
779 | snew(remainder, nlim)(remainder) = save_calloc("remainder", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 779, (nlim), sizeof(*(remainder))); |
780 | |
781 | for (i = 0; i < expand->lmc_repeats; i++) |
782 | { |
783 | double rnd[2]; |
784 | |
785 | gmx_rng_cycle_2uniform(step, i, seed, RND_SEED_EXPANDED6, rnd); |
786 | |
787 | for (ifep = 0; ifep < nlim; ifep++) |
788 | { |
789 | propose[ifep] = 0; |
790 | accept[ifep] = 0; |
791 | } |
792 | |
793 | if ((expand->elmcmove == elmcmoveGIBBS) || (expand->elmcmove == elmcmoveMETGIBBS)) |
794 | { |
795 | bRestricted = TRUE1; |
796 | /* use the Gibbs sampler, with restricted range */ |
797 | if (expand->gibbsdeltalam < 0) |
798 | { |
799 | minfep = 0; |
800 | maxfep = nlim-1; |
801 | bRestricted = FALSE0; |
802 | } |
803 | else |
804 | { |
805 | minfep = fep_state - expand->gibbsdeltalam; |
806 | maxfep = fep_state + expand->gibbsdeltalam; |
807 | if (minfep < 0) |
808 | { |
809 | minfep = 0; |
810 | } |
811 | if (maxfep > nlim-1) |
812 | { |
813 | maxfep = nlim-1; |
814 | } |
815 | } |
816 | |
817 | GenerateGibbsProbabilities(weighted_lamee, p_k, &pks, minfep, maxfep); |
818 | |
819 | if (expand->elmcmove == elmcmoveGIBBS) |
820 | { |
821 | for (ifep = minfep; ifep <= maxfep; ifep++) |
822 | { |
823 | propose[ifep] = p_k[ifep]; |
824 | accept[ifep] = 1.0; |
825 | } |
826 | /* Gibbs sampling */ |
827 | r1 = rnd[0]; |
828 | for (lamnew = minfep; lamnew <= maxfep; lamnew++) |
829 | { |
830 | if (r1 <= p_k[lamnew]) |
831 | { |
832 | break; |
833 | } |
834 | r1 -= p_k[lamnew]; |
835 | } |
836 | } |
837 | else if (expand->elmcmove == elmcmoveMETGIBBS) |
838 | { |
839 | |
840 | /* Metropolized Gibbs sampling */ |
841 | for (ifep = minfep; ifep <= maxfep; ifep++) |
842 | { |
843 | remainder[ifep] = 1 - p_k[ifep]; |
844 | } |
845 | |
846 | /* find the proposal probabilities */ |
847 | |
848 | if (remainder[fep_state] == 0) |
849 | { |
850 | /* only the current state has any probability */ |
851 | /* we have to stay at the current state */ |
852 | lamnew = fep_state; |
853 | } |
854 | else |
855 | { |
856 | for (ifep = minfep; ifep <= maxfep; ifep++) |
857 | { |
858 | if (ifep != fep_state) |
859 | { |
860 | propose[ifep] = p_k[ifep]/remainder[fep_state]; |
861 | } |
862 | else |
863 | { |
864 | propose[ifep] = 0; |
865 | } |
866 | } |
867 | |
868 | r1 = rnd[0]; |
869 | for (lamtrial = minfep; lamtrial <= maxfep; lamtrial++) |
870 | { |
871 | pnorm = p_k[lamtrial]/remainder[fep_state]; |
872 | if (lamtrial != fep_state) |
873 | { |
874 | if (r1 <= pnorm) |
875 | { |
876 | break; |
877 | } |
878 | r1 -= pnorm; |
879 | } |
880 | } |
881 | |
882 | /* we have now selected lamtrial according to p(lamtrial)/1-p(fep_state) */ |
883 | tprob = 1.0; |
884 | /* trial probability is min{1,\frac{1 - p(old)}{1-p(new)} MRS 1/8/2008 */ |
885 | trialprob = (remainder[fep_state])/(remainder[lamtrial]); |
886 | if (trialprob < tprob) |
887 | { |
888 | tprob = trialprob; |
889 | } |
890 | r2 = rnd[1]; |
891 | if (r2 < tprob) |
892 | { |
893 | lamnew = lamtrial; |
894 | } |
895 | else |
896 | { |
897 | lamnew = fep_state; |
898 | } |
899 | } |
900 | |
901 | /* now figure out the acceptance probability for each */ |
902 | for (ifep = minfep; ifep <= maxfep; ifep++) |
903 | { |
904 | tprob = 1.0; |
905 | if (remainder[ifep] != 0) |
906 | { |
907 | trialprob = (remainder[fep_state])/(remainder[ifep]); |
908 | } |
909 | else |
910 | { |
911 | trialprob = 1.0; /* this state is the only choice! */ |
912 | } |
913 | if (trialprob < tprob) |
914 | { |
915 | tprob = trialprob; |
916 | } |
917 | /* probability for fep_state=0, but that's fine, it's never proposed! */ |
918 | accept[ifep] = tprob; |
919 | } |
920 | } |
921 | |
922 | if (lamnew > maxfep) |
923 | { |
924 | /* it's possible some rounding is failing */ |
925 | if (gmx_within_tol(remainder[fep_state], 0, 50*GMX_DOUBLE_EPS1.11022302E-16)) |
926 | { |
927 | /* numerical rounding error -- no state other than the original has weight */ |
928 | lamnew = fep_state; |
929 | } |
930 | else |
931 | { |
932 | /* probably not a numerical issue */ |
933 | int loc = 0; |
934 | int nerror = 200+(maxfep-minfep+1)*60; |
935 | char *errorstr; |
936 | snew(errorstr, nerror)(errorstr) = save_calloc("errorstr", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 936, (nerror), sizeof(*(errorstr))); |
937 | /* if its greater than maxfep, then something went wrong -- probably underflow in the calculation |
938 | of sum weights. Generated detailed info for failure */ |
939 | loc += sprintf(errorstr, "Something wrong in choosing new lambda state with a Gibbs move -- probably underflow in weight determination.\nDenominator is: %3d%17.10e\n i dE numerator weights\n", 0, pks); |
940 | for (ifep = minfep; ifep <= maxfep; ifep++) |
941 | { |
942 | loc += sprintf(&errorstr[loc], "%3d %17.10e%17.10e%17.10e\n", ifep, weighted_lamee[ifep], p_k[ifep], dfhist->sum_weights[ifep]); |
943 | } |
944 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 944, errorstr); |
945 | } |
946 | } |
947 | } |
948 | else if ((expand->elmcmove == elmcmoveMETROPOLIS) || (expand->elmcmove == elmcmoveBARKER)) |
949 | { |
950 | /* use the metropolis sampler with trial +/- 1 */ |
951 | r1 = rnd[0]; |
952 | if (r1 < 0.5) |
953 | { |
954 | if (fep_state == 0) |
955 | { |
956 | lamtrial = fep_state; |
957 | } |
958 | else |
959 | { |
960 | lamtrial = fep_state-1; |
961 | } |
962 | } |
963 | else |
964 | { |
965 | if (fep_state == nlim-1) |
966 | { |
967 | lamtrial = fep_state; |
968 | } |
969 | else |
970 | { |
971 | lamtrial = fep_state+1; |
972 | } |
973 | } |
974 | |
975 | de = weighted_lamee[lamtrial] - weighted_lamee[fep_state]; |
976 | if (expand->elmcmove == elmcmoveMETROPOLIS) |
977 | { |
978 | tprob = 1.0; |
979 | trialprob = exp(de); |
980 | if (trialprob < tprob) |
981 | { |
982 | tprob = trialprob; |
983 | } |
984 | propose[fep_state] = 0; |
985 | propose[lamtrial] = 1.0; /* note that this overwrites the above line if fep_state = ntrial, which only occurs at the ends */ |
986 | accept[fep_state] = 1.0; /* doesn't actually matter, never proposed unless fep_state = ntrial, in which case it's 1.0 anyway */ |
987 | accept[lamtrial] = tprob; |
988 | |
989 | } |
990 | else if (expand->elmcmove == elmcmoveBARKER) |
991 | { |
992 | tprob = 1.0/(1.0+exp(-de)); |
993 | |
994 | propose[fep_state] = (1-tprob); |
995 | propose[lamtrial] += tprob; /* we add, to account for the fact that at the end, they might be the same point */ |
996 | accept[fep_state] = 1.0; |
997 | accept[lamtrial] = 1.0; |
998 | } |
999 | |
1000 | r2 = rnd[1]; |
1001 | if (r2 < tprob) |
1002 | { |
1003 | lamnew = lamtrial; |
1004 | } |
1005 | else |
1006 | { |
1007 | lamnew = fep_state; |
1008 | } |
1009 | } |
1010 | |
1011 | for (ifep = 0; ifep < nlim; ifep++) |
1012 | { |
1013 | dfhist->Tij[fep_state][ifep] += propose[ifep]*accept[ifep]; |
1014 | dfhist->Tij[fep_state][fep_state] += propose[ifep]*(1.0-accept[ifep]); |
1015 | } |
1016 | fep_state = lamnew; |
1017 | } |
1018 | |
1019 | dfhist->Tij_empirical[starting_fep_state][lamnew] += 1.0; |
1020 | |
1021 | sfree(propose)save_free("propose", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1021, (propose)); |
1022 | sfree(accept)save_free("accept", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1022, (accept)); |
1023 | sfree(remainder)save_free("remainder", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1023, (remainder)); |
1024 | |
1025 | return lamnew; |
1026 | } |
1027 | |
1028 | /* print out the weights to the log, along with current state */ |
1029 | extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *expand, t_simtemp *simtemp, df_history_t *dfhist, |
1030 | int fep_state, int frequency, gmx_int64_t step) |
1031 | { |
1032 | int nlim, i, ifep, jfep; |
1033 | real dw, dg, dv, dm, Tprint; |
1034 | real *temps; |
1035 | const char *print_names[efptNR] = {" FEPL", "MassL", "CoulL", " VdwL", "BondL", "RestT", "Temp.(K)"}; |
1036 | gmx_bool bSimTemp = FALSE0; |
1037 | |
1038 | nlim = fep->n_lambda; |
1039 | if (simtemp != NULL((void*)0)) |
1040 | { |
1041 | bSimTemp = TRUE1; |
1042 | } |
1043 | |
1044 | if (mod(step, frequency)_mod((step), (frequency), "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1044) == 0) |
1045 | { |
1046 | fprintf(outfile, " MC-lambda information\n"); |
1047 | if (EWL(expand->elamstats)((expand->elamstats) == elamstatsWL || (expand->elamstats ) == elamstatsWWL) && (!(dfhist->bEquil))) |
1048 | { |
1049 | fprintf(outfile, " Wang-Landau incrementor is: %11.5g\n", dfhist->wl_delta); |
1050 | } |
1051 | fprintf(outfile, " N"); |
1052 | for (i = 0; i < efptNR; i++) |
1053 | { |
1054 | if (fep->separate_dvdl[i]) |
1055 | { |
1056 | fprintf(outfile, "%7s", print_names[i]); |
1057 | } |
1058 | else if ((i == efptTEMPERATURE) && bSimTemp) |
1059 | { |
1060 | fprintf(outfile, "%10s", print_names[i]); /* more space for temperature formats */ |
1061 | } |
1062 | } |
1063 | fprintf(outfile, " Count "); |
1064 | if (expand->elamstats == elamstatsMINVAR) |
1065 | { |
1066 | fprintf(outfile, "W(in kT) G(in kT) dG(in kT) dV(in kT)\n"); |
1067 | } |
1068 | else |
1069 | { |
1070 | fprintf(outfile, "G(in kT) dG(in kT)\n"); |
1071 | } |
1072 | for (ifep = 0; ifep < nlim; ifep++) |
1073 | { |
1074 | if (ifep == nlim-1) |
1075 | { |
1076 | dw = 0.0; |
1077 | dg = 0.0; |
1078 | dv = 0.0; |
1079 | dm = 0.0; |
Value stored to 'dm' is never read | |
1080 | } |
1081 | else |
1082 | { |
1083 | dw = dfhist->sum_weights[ifep+1] - dfhist->sum_weights[ifep]; |
1084 | dg = dfhist->sum_dg[ifep+1] - dfhist->sum_dg[ifep]; |
1085 | dv = sqrt(pow(dfhist->sum_variance[ifep+1], 2) - pow(dfhist->sum_variance[ifep], 2)); |
1086 | dm = dfhist->sum_minvar[ifep+1] - dfhist->sum_minvar[ifep]; |
1087 | |
1088 | } |
1089 | fprintf(outfile, "%3d", (ifep+1)); |
1090 | for (i = 0; i < efptNR; i++) |
1091 | { |
1092 | if (fep->separate_dvdl[i]) |
1093 | { |
1094 | fprintf(outfile, "%7.3f", fep->all_lambda[i][ifep]); |
1095 | } |
1096 | else if (i == efptTEMPERATURE && bSimTemp) |
1097 | { |
1098 | fprintf(outfile, "%9.3f", simtemp->temperatures[ifep]); |
1099 | } |
1100 | } |
1101 | if (EWL(expand->elamstats)((expand->elamstats) == elamstatsWL || (expand->elamstats ) == elamstatsWWL) && (!(dfhist->bEquil))) /* if performing WL and still haven't equilibrated */ |
1102 | { |
1103 | if (expand->elamstats == elamstatsWL) |
1104 | { |
1105 | fprintf(outfile, " %8d", (int)dfhist->wl_histo[ifep]); |
1106 | } |
1107 | else |
1108 | { |
1109 | fprintf(outfile, " %8.3f", dfhist->wl_histo[ifep]); |
1110 | } |
1111 | } |
1112 | else /* we have equilibrated weights */ |
1113 | { |
1114 | fprintf(outfile, " %8d", dfhist->n_at_lam[ifep]); |
1115 | } |
1116 | if (expand->elamstats == elamstatsMINVAR) |
1117 | { |
1118 | fprintf(outfile, " %10.5f %10.5f %10.5f %10.5f", dfhist->sum_weights[ifep], dfhist->sum_dg[ifep], dg, dv); |
1119 | } |
1120 | else |
1121 | { |
1122 | fprintf(outfile, " %10.5f %10.5f", dfhist->sum_weights[ifep], dw); |
1123 | } |
1124 | if (ifep == fep_state) |
1125 | { |
1126 | fprintf(outfile, " <<\n"); |
1127 | } |
1128 | else |
1129 | { |
1130 | fprintf(outfile, " \n"); |
1131 | } |
1132 | } |
1133 | fprintf(outfile, "\n"); |
1134 | |
1135 | if ((mod(step, expand->nstTij)_mod((step), (expand->nstTij), "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1135) == 0) && (expand->nstTij > 0) && (step > 0)) |
1136 | { |
1137 | fprintf(outfile, " Transition Matrix\n"); |
1138 | for (ifep = 0; ifep < nlim; ifep++) |
1139 | { |
1140 | fprintf(outfile, "%12d", (ifep+1)); |
1141 | } |
1142 | fprintf(outfile, "\n"); |
1143 | for (ifep = 0; ifep < nlim; ifep++) |
1144 | { |
1145 | for (jfep = 0; jfep < nlim; jfep++) |
1146 | { |
1147 | if (dfhist->n_at_lam[ifep] > 0) |
1148 | { |
1149 | if (expand->bSymmetrizedTMatrix) |
1150 | { |
1151 | Tprint = (dfhist->Tij[ifep][jfep]+dfhist->Tij[jfep][ifep])/(dfhist->n_at_lam[ifep]+dfhist->n_at_lam[jfep]); |
1152 | } |
1153 | else |
1154 | { |
1155 | Tprint = (dfhist->Tij[ifep][jfep])/(dfhist->n_at_lam[ifep]); |
1156 | } |
1157 | } |
1158 | else |
1159 | { |
1160 | Tprint = 0.0; |
1161 | } |
1162 | fprintf(outfile, "%12.8f", Tprint); |
1163 | } |
1164 | fprintf(outfile, "%3d\n", (ifep+1)); |
1165 | } |
1166 | |
1167 | fprintf(outfile, " Empirical Transition Matrix\n"); |
1168 | for (ifep = 0; ifep < nlim; ifep++) |
1169 | { |
1170 | fprintf(outfile, "%12d", (ifep+1)); |
1171 | } |
1172 | fprintf(outfile, "\n"); |
1173 | for (ifep = 0; ifep < nlim; ifep++) |
1174 | { |
1175 | for (jfep = 0; jfep < nlim; jfep++) |
1176 | { |
1177 | if (dfhist->n_at_lam[ifep] > 0) |
1178 | { |
1179 | if (expand->bSymmetrizedTMatrix) |
1180 | { |
1181 | Tprint = (dfhist->Tij_empirical[ifep][jfep]+dfhist->Tij_empirical[jfep][ifep])/(dfhist->n_at_lam[ifep]+dfhist->n_at_lam[jfep]); |
1182 | } |
1183 | else |
1184 | { |
1185 | Tprint = dfhist->Tij_empirical[ifep][jfep]/(dfhist->n_at_lam[ifep]); |
1186 | } |
1187 | } |
1188 | else |
1189 | { |
1190 | Tprint = 0.0; |
1191 | } |
1192 | fprintf(outfile, "%12.8f", Tprint); |
1193 | } |
1194 | fprintf(outfile, "%3d\n", (ifep+1)); |
1195 | } |
1196 | } |
1197 | } |
1198 | } |
1199 | |
1200 | extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *enerd, |
1201 | t_state *state, t_extmass *MassQ, int fep_state, df_history_t *dfhist, |
1202 | gmx_int64_t step, |
1203 | rvec *v, t_mdatoms *mdatoms) |
1204 | /* Note that the state variable is only needed for simulated tempering, not |
1205 | Hamiltonian expanded ensemble. May be able to remove it after integrator refactoring. */ |
1206 | { |
1207 | real *pfep_lamee, *scaled_lamee, *weighted_lamee; |
1208 | double *p_k; |
1209 | int i, nlim, lamnew, totalsamples; |
1210 | real oneovert, maxscaled = 0, maxweighted = 0; |
1211 | t_expanded *expand; |
1212 | t_simtemp *simtemp; |
1213 | double *temperature_lambdas; |
1214 | gmx_bool bIfReset, bSwitchtoOneOverT, bDoneEquilibrating = FALSE0; |
1215 | |
1216 | expand = ir->expandedvals; |
1217 | simtemp = ir->simtempvals; |
1218 | nlim = ir->fepvals->n_lambda; |
1219 | |
1220 | snew(scaled_lamee, nlim)(scaled_lamee) = save_calloc("scaled_lamee", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1220, (nlim), sizeof(*(scaled_lamee))); |
1221 | snew(weighted_lamee, nlim)(weighted_lamee) = save_calloc("weighted_lamee", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1221, (nlim), sizeof(*(weighted_lamee))); |
1222 | snew(pfep_lamee, nlim)(pfep_lamee) = save_calloc("pfep_lamee", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1222, (nlim), sizeof(*(pfep_lamee))); |
1223 | snew(p_k, nlim)(p_k) = save_calloc("p_k", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1223, (nlim), sizeof(*(p_k))); |
1224 | |
1225 | /* update the count at the current lambda*/ |
1226 | dfhist->n_at_lam[fep_state]++; |
1227 | |
1228 | /* need to calculate the PV term somewhere, but not needed here? Not until there's a lambda state that's |
1229 | pressure controlled.*/ |
1230 | /* |
1231 | pVTerm = 0; |
1232 | where does this PV term go? |
1233 | for (i=0;i<nlim;i++) |
1234 | { |
1235 | fep_lamee[i] += pVTerm; |
1236 | } |
1237 | */ |
1238 | |
1239 | /* determine the minimum value to avoid overflow. Probably a better way to do this */ |
1240 | /* we don't need to include the pressure term, since the volume is the same between the two. |
1241 | is there some term we are neglecting, however? */ |
1242 | |
1243 | if (ir->efep != efepNO) |
1244 | { |
1245 | for (i = 0; i < nlim; i++) |
1246 | { |
1247 | if (ir->bSimTemp) |
1248 | { |
1249 | /* Note -- this assumes no mass changes, since kinetic energy is not added . . . */ |
1250 | scaled_lamee[i] = (enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0])/(simtemp->temperatures[i]*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))) |
1251 | + enerd->term[F_EPOT]*(1.0/(simtemp->temperatures[i])- 1.0/(simtemp->temperatures[fep_state]))/BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3)); |
1252 | } |
1253 | else |
1254 | { |
1255 | scaled_lamee[i] = (enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0])/(expand->mc_temp*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))); |
1256 | /* mc_temp is currently set to the system reft unless otherwise defined */ |
1257 | } |
1258 | |
1259 | /* save these energies for printing, so they don't get overwritten by the next step */ |
1260 | /* they aren't overwritten in the non-free energy case, but we always print with these |
1261 | for simplicity */ |
1262 | } |
1263 | } |
1264 | else |
1265 | { |
1266 | if (ir->bSimTemp) |
1267 | { |
1268 | for (i = 0; i < nlim; i++) |
1269 | { |
1270 | scaled_lamee[i] = enerd->term[F_EPOT]*(1.0/simtemp->temperatures[i] - 1.0/simtemp->temperatures[fep_state])/BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3)); |
1271 | } |
1272 | } |
1273 | } |
1274 | |
1275 | for (i = 0; i < nlim; i++) |
1276 | { |
1277 | pfep_lamee[i] = scaled_lamee[i]; |
1278 | |
1279 | weighted_lamee[i] = dfhist->sum_weights[i] - scaled_lamee[i]; |
1280 | if (i == 0) |
1281 | { |
1282 | maxscaled = scaled_lamee[i]; |
1283 | maxweighted = weighted_lamee[i]; |
1284 | } |
1285 | else |
1286 | { |
1287 | if (scaled_lamee[i] > maxscaled) |
1288 | { |
1289 | maxscaled = scaled_lamee[i]; |
1290 | } |
1291 | if (weighted_lamee[i] > maxweighted) |
1292 | { |
1293 | maxweighted = weighted_lamee[i]; |
1294 | } |
1295 | } |
1296 | } |
1297 | |
1298 | for (i = 0; i < nlim; i++) |
1299 | { |
1300 | scaled_lamee[i] -= maxscaled; |
1301 | weighted_lamee[i] -= maxweighted; |
1302 | } |
1303 | |
1304 | /* update weights - we decide whether or not to actually do this inside */ |
1305 | |
1306 | bDoneEquilibrating = UpdateWeights(nlim, expand, dfhist, fep_state, scaled_lamee, weighted_lamee, step); |
1307 | if (bDoneEquilibrating) |
1308 | { |
1309 | if (log) |
1310 | { |
1311 | fprintf(log, "\nStep %d: Weights have equilibrated, using criteria: %s\n", (int)step, elmceq_names[expand->elmceq]); |
1312 | } |
1313 | } |
1314 | |
1315 | lamnew = ChooseNewLambda(nlim, expand, dfhist, fep_state, weighted_lamee, p_k, |
1316 | ir->expandedvals->lmc_seed, step); |
1317 | /* if using simulated tempering, we need to adjust the temperatures */ |
1318 | if (ir->bSimTemp && (lamnew != fep_state)) /* only need to change the temperatures if we change the state */ |
1319 | { |
1320 | int i, j, n, d; |
1321 | real *buf_ngtc; |
1322 | real told; |
1323 | int nstart, nend, gt; |
1324 | |
1325 | snew(buf_ngtc, ir->opts.ngtc)(buf_ngtc) = save_calloc("buf_ngtc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1325, (ir->opts.ngtc), sizeof(*(buf_ngtc))); |
1326 | |
1327 | for (i = 0; i < ir->opts.ngtc; i++) |
1328 | { |
1329 | if (ir->opts.ref_t[i] > 0) |
1330 | { |
1331 | told = ir->opts.ref_t[i]; |
1332 | ir->opts.ref_t[i] = simtemp->temperatures[lamnew]; |
1333 | buf_ngtc[i] = sqrt(ir->opts.ref_t[i]/told); /* using the buffer as temperature scaling */ |
1334 | } |
1335 | } |
1336 | |
1337 | /* we don't need to manipulate the ekind information, as it isn't due to be reset until the next step anyway */ |
1338 | |
1339 | nstart = 0; |
1340 | nend = mdatoms->homenr; |
1341 | for (n = nstart; n < nend; n++) |
1342 | { |
1343 | gt = 0; |
1344 | if (mdatoms->cTC) |
1345 | { |
1346 | gt = mdatoms->cTC[n]; |
1347 | } |
1348 | for (d = 0; d < DIM3; d++) |
1349 | { |
1350 | v[n][d] *= buf_ngtc[gt]; |
1351 | } |
1352 | } |
1353 | |
1354 | if (IR_NPT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && ((ir)->etc == etcNOSEHOOVER ))) || IR_NPH_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && (!(((ir)->etc == etcNOSEHOOVER ))))) || IR_NVT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((!((ir)->epc == epcMTTK)) && ((ir)->etc == etcNOSEHOOVER )))) |
1355 | { |
1356 | /* we need to recalculate the masses if the temperature has changed */ |
1357 | init_npt_masses(ir, state, MassQ, FALSE0); |
1358 | for (i = 0; i < state->nnhpres; i++) |
1359 | { |
1360 | for (j = 0; j < ir->opts.nhchainlength; j++) |
1361 | { |
1362 | state->nhpres_vxi[i+j] *= buf_ngtc[i]; |
1363 | } |
1364 | } |
1365 | for (i = 0; i < ir->opts.ngtc; i++) |
1366 | { |
1367 | for (j = 0; j < ir->opts.nhchainlength; j++) |
1368 | { |
1369 | state->nosehoover_vxi[i+j] *= buf_ngtc[i]; |
1370 | } |
1371 | } |
1372 | } |
1373 | sfree(buf_ngtc)save_free("buf_ngtc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1373, (buf_ngtc)); |
1374 | } |
1375 | |
1376 | /* now check on the Wang-Landau updating critera */ |
1377 | |
1378 | if (EWL(expand->elamstats)((expand->elamstats) == elamstatsWL || (expand->elamstats ) == elamstatsWWL)) |
1379 | { |
1380 | bSwitchtoOneOverT = FALSE0; |
1381 | if (expand->bWLoneovert) |
1382 | { |
1383 | totalsamples = 0; |
1384 | for (i = 0; i < nlim; i++) |
1385 | { |
1386 | totalsamples += dfhist->n_at_lam[i]; |
1387 | } |
1388 | oneovert = (1.0*nlim)/totalsamples; |
1389 | /* oneovert has decreasd by a bit since last time, so we actually make sure its within one of this number */ |
1390 | /* switch to 1/t incrementing when wl_delta has decreased at least once, and wl_delta is now less than 1/t */ |
1391 | if ((dfhist->wl_delta <= ((totalsamples)/(totalsamples-1.00001))*oneovert) && |
1392 | (dfhist->wl_delta < expand->init_wl_delta)) |
1393 | { |
1394 | bSwitchtoOneOverT = TRUE1; |
1395 | } |
1396 | } |
1397 | if (bSwitchtoOneOverT) |
1398 | { |
1399 | dfhist->wl_delta = oneovert; /* now we reduce by this each time, instead of only at flatness */ |
1400 | } |
1401 | else |
1402 | { |
1403 | bIfReset = CheckHistogramRatios(nlim, dfhist->wl_histo, expand->wl_ratio); |
1404 | if (bIfReset) |
1405 | { |
1406 | for (i = 0; i < nlim; i++) |
1407 | { |
1408 | dfhist->wl_histo[i] = 0; |
1409 | } |
1410 | dfhist->wl_delta *= expand->wl_scale; |
1411 | if (log) |
1412 | { |
1413 | fprintf(log, "\nStep %d: weights are now:", (int)step); |
1414 | for (i = 0; i < nlim; i++) |
1415 | { |
1416 | fprintf(log, " %.5f", dfhist->sum_weights[i]); |
1417 | } |
1418 | fprintf(log, "\n"); |
1419 | } |
1420 | } |
1421 | } |
1422 | } |
1423 | sfree(pfep_lamee)save_free("pfep_lamee", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1423, (pfep_lamee)); |
1424 | sfree(scaled_lamee)save_free("scaled_lamee", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1424, (scaled_lamee)); |
1425 | sfree(weighted_lamee)save_free("weighted_lamee", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1425, (weighted_lamee)); |
1426 | sfree(p_k)save_free("p_k", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/expanded.c" , 1426, (p_k)); |
1427 | |
1428 | return lamnew; |
1429 | } |