Bug Summary

File:gromacs/gmxana/gmx_bar.c
Location:line 3756, column 5
Description:Value stored to 'beta' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2010,2011,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 <ctype.h>
40#include <float.h>
41#include <math.h>
42#include <stdlib.h>
43#include <string.h>
44
45#include "typedefs.h"
46#include "gromacs/utility/smalloc.h"
47#include "gromacs/utility/futil.h"
48#include "gromacs/commandline/pargs.h"
49#include "macros.h"
50#include "gromacs/fileio/enxio.h"
51#include "physics.h"
52#include "gromacs/utility/fatalerror.h"
53#include "gromacs/fileio/xvgr.h"
54#include "viewit.h"
55#include "gmx_ana.h"
56#include "gromacs/math/utilities.h"
57#include "gromacs/utility/cstringutil.h"
58#include "names.h"
59#include "mdebin.h"
60
61
62/* Structure for the names of lambda vector components */
63typedef struct lambda_components_t
64{
65 char **names; /* Array of strings with names for the lambda vector
66 components */
67 int N; /* The number of components */
68 int Nalloc; /* The number of allocated components */
69} lambda_components_t;
70
71/* Structure for a lambda vector or a dhdl derivative direction */
72typedef struct lambda_vec_t
73{
74 double *val; /* The lambda vector component values. Only valid if
75 dhdl == -1 */
76 int dhdl; /* The coordinate index for the derivative described by this
77 structure, or -1 */
78 const lambda_components_t *lc; /* the associated lambda_components
79 structure */
80 int index; /* The state number (init-lambda-state) of this lambda
81 vector, if known. If not, it is set to -1 */
82} lambda_vec_t;
83
84/* the dhdl.xvg data from a simulation */
85typedef struct xvg_t
86{
87 const char *filename;
88 int ftp; /* file type */
89 int nset; /* number of lambdas, including dhdl */
90 int *np; /* number of data points (du or hists) per lambda */
91 int np_alloc; /* number of points (du or hists) allocated */
92 double temp; /* temperature */
93 lambda_vec_t *lambda; /* the lambdas (of first index for y). */
94 double *t; /* the times (of second index for y) */
95 double **y; /* the dU values. y[0] holds the derivative, while
96 further ones contain the energy differences between
97 the native lambda and the 'foreign' lambdas. */
98 lambda_vec_t native_lambda; /* the native lambda */
99
100 struct xvg_t *next, *prev; /*location in the global linked list of xvg_ts*/
101} xvg_t;
102
103
104typedef struct hist_t
105{
106 unsigned int *bin[2]; /* the (forward + reverse) histogram values */
107 double dx[2]; /* the histogram spacing. The reverse
108 dx is the negative of the forward dx.*/
109 gmx_int64_t x0[2]; /* the (forward + reverse) histogram start
110 point(s) as int */
111
112 int nbin[2]; /* the (forward+reverse) number of bins */
113 gmx_int64_t sum; /* the total number of counts. Must be
114 the same for forward + reverse. */
115 int nhist; /* number of hist datas (forward or reverse) */
116
117 double start_time, delta_time; /* start time, end time of histogram */
118} hist_t;
119
120
121/* an aggregate of samples for partial free energy calculation */
122typedef struct samples_t
123{
124 lambda_vec_t *native_lambda; /* pointer to native lambda vector */
125 lambda_vec_t *foreign_lambda; /* pointer to foreign lambda vector */
126 double temp; /* the temperature */
127 gmx_bool derivative; /* whether this sample is a derivative */
128
129 /* The samples come either as either delta U lists: */
130 int ndu; /* the number of delta U samples */
131 double *du; /* the delta u's */
132 double *t; /* the times associated with those samples, or: */
133 double start_time, delta_time; /*start time and delta time for linear time*/
134
135 /* or as histograms: */
136 hist_t *hist; /* a histogram */
137
138 /* allocation data: (not NULL for data 'owned' by this struct) */
139 double *du_alloc, *t_alloc; /* allocated delta u arrays */
140 size_t ndu_alloc, nt_alloc; /* pre-allocated sizes */
141 hist_t *hist_alloc; /* allocated hist */
142
143 gmx_int64_t ntot; /* total number of samples */
144 const char *filename; /* the file name this sample comes from */
145} samples_t;
146
147/* a sample range (start to end for du-style data, or boolean
148 for both du-style data and histograms */
149typedef struct sample_range_t
150{
151 int start, end; /* start and end index for du style data */
152 gmx_bool use; /* whether to use this sample */
153
154 samples_t *s; /* the samples this range belongs to */
155} sample_range_t;
156
157
158/* a collection of samples for a partial free energy calculation
159 (i.e. the collection of samples from one native lambda to one
160 foreign lambda) */
161typedef struct sample_coll_t
162{
163 lambda_vec_t *native_lambda; /* these should be the same for all samples
164 in the histogram */
165 lambda_vec_t *foreign_lambda; /* collection */
166 double temp; /* the temperature */
167
168 int nsamples; /* the number of samples */
169 samples_t **s; /* the samples themselves */
170 sample_range_t *r; /* the sample ranges */
171 int nsamples_alloc; /* number of allocated samples */
172
173 gmx_int64_t ntot; /* total number of samples in the ranges of
174 this collection */
175
176 struct sample_coll_t *next, *prev; /* next and previous in the list */
177} sample_coll_t;
178
179/* all the samples associated with a lambda point */
180typedef struct lambda_data_t
181{
182 lambda_vec_t *lambda; /* the native lambda (at start time if dynamic) */
183 double temp; /* temperature */
184
185 sample_coll_t *sc; /* the samples */
186
187 sample_coll_t sc_head; /*the pre-allocated list head for the linked list.*/
188
189 struct lambda_data_t *next, *prev; /* the next and prev in the list */
190} lambda_data_t;
191
192/* Top-level data structure of simulation data */
193typedef struct sim_data_t
194{
195 lambda_data_t *lb; /* a lambda data linked list */
196 lambda_data_t lb_head; /* The head element of the linked list */
197
198 lambda_components_t lc; /* the allowed components of the lambda
199 vectors */
200} sim_data_t;
201
202/* Top-level data structure with calculated values. */
203typedef struct {
204 sample_coll_t *a, *b; /* the simulation data */
205
206 double dg; /* the free energy difference */
207 double dg_err; /* the free energy difference */
208
209 double dg_disc_err; /* discretization error */
210 double dg_histrange_err; /* histogram range error */
211
212 double sa; /* relative entropy of b in state a */
213 double sa_err; /* error in sa */
214 double sb; /* relative entropy of a in state b */
215 double sb_err; /* error in sb */
216
217 double dg_stddev; /* expected dg stddev per sample */
218 double dg_stddev_err; /* error in dg_stddev */
219} barres_t;
220
221
222/* Initialize a lambda_components structure */
223static void lambda_components_init(lambda_components_t *lc)
224{
225 lc->N = 0;
226 lc->Nalloc = 2;
227 snew(lc->names, lc->Nalloc)(lc->names) = save_calloc("lc->names", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 227, (lc->Nalloc), sizeof(*(lc->names)));
228}
229
230/* Add a component to a lambda_components structure */
231static void lambda_components_add(lambda_components_t *lc,
232 const char *name, size_t name_length)
233{
234 while (lc->N + 1 > lc->Nalloc)
235 {
236 lc->Nalloc = (lc->Nalloc == 0) ? 2 : 2*lc->Nalloc;
237 srenew(lc->names, lc->Nalloc)(lc->names) = save_realloc("lc->names", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 237, (lc->names), (lc->Nalloc), sizeof(*(lc->names
)));
238 }
239 snew(lc->names[lc->N], name_length+1)(lc->names[lc->N]) = save_calloc("lc->names[lc->N]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 239, (name_length+1), sizeof(*(lc->names[lc->N])));
240 strncpy(lc->names[lc->N], name, name_length)__builtin_strncpy (lc->names[lc->N], name, name_length);
241 lc->N++;
242}
243
244/* check whether a component with index 'index' matches the given name, or
245 is also NULL. Returns TRUE if this is the case.
246 the string name does not need to end */
247static gmx_bool lambda_components_check(const lambda_components_t *lc,
248 int index,
249 const char *name,
250 size_t name_length)
251{
252 size_t len;
253 if (index >= lc->N)
254 {
255 return FALSE0;
256 }
257 if (name == NULL((void*)0) && lc->names[index] == NULL((void*)0))
258 {
259 return TRUE1;
260 }
261 if ((name == NULL((void*)0)) != (lc->names[index] == NULL((void*)0)))
262 {
263 return FALSE0;
264 }
265 len = strlen(lc->names[index]);
266 if (len != name_length)
267 {
268 return FALSE0;
269 }
270 if (strncmp(lc->names[index], name, name_length)(__extension__ (__builtin_constant_p (name_length) &&
((__builtin_constant_p (lc->names[index]) && strlen
(lc->names[index]) < ((size_t) (name_length))) || (__builtin_constant_p
(name) && strlen (name) < ((size_t) (name_length)
))) ? __extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(lc->names[index]) && __builtin_constant_p (name)
&& (__s1_len = strlen (lc->names[index]), __s2_len
= strlen (name), (!((size_t)(const void *)((lc->names[index
]) + 1) - (size_t)(const void *)(lc->names[index]) == 1) ||
__s1_len >= 4) && (!((size_t)(const void *)((name
) + 1) - (size_t)(const void *)(name) == 1) || __s2_len >=
4)) ? __builtin_strcmp (lc->names[index], name) : (__builtin_constant_p
(lc->names[index]) && ((size_t)(const void *)((lc
->names[index]) + 1) - (size_t)(const void *)(lc->names
[index]) == 1) && (__s1_len = strlen (lc->names[index
]), __s1_len < 4) ? (__builtin_constant_p (name) &&
((size_t)(const void *)((name) + 1) - (size_t)(const void *)
(name) == 1) ? __builtin_strcmp (lc->names[index], name) :
(__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (name); int __result = (((const unsigned
char *) (const char *) (lc->names[index]))[0] - __s2[0]);
if (__s1_len > 0 && __result == 0) { __result = (
((const unsigned char *) (const char *) (lc->names[index])
)[1] - __s2[1]); if (__s1_len > 1 && __result == 0
) { __result = (((const unsigned char *) (const char *) (lc->
names[index]))[2] - __s2[2]); if (__s1_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (lc
->names[index]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
(name) && ((size_t)(const void *)((name) + 1) - (size_t
)(const void *)(name) == 1) && (__s2_len = strlen (name
), __s2_len < 4) ? (__builtin_constant_p (lc->names[index
]) && ((size_t)(const void *)((lc->names[index]) +
1) - (size_t)(const void *)(lc->names[index]) == 1) ? __builtin_strcmp
(lc->names[index], name) : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (lc->
names[index]); int __result = (((const unsigned char *) (const
char *) (name))[0] - __s2[0]); if (__s2_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (name))[1] - __s2[1]); if (__s2_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (name))[2] - __s2[2]); if (__s2_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (name))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(lc->names[index], name)))); }) : strncmp (lc->names[index
], name, name_length))) == 0)
271 {
272 return TRUE1;
273 }
274 return FALSE0;
275}
276
277/* Find the index of a given lambda component name, or -1 if not found */
278static int lambda_components_find(const lambda_components_t *lc,
279 const char *name,
280 size_t name_length)
281{
282 int i;
283
284 for (i = 0; i < lc->N; i++)
285 {
286 if (strncmp(lc->names[i], name, name_length)(__extension__ (__builtin_constant_p (name_length) &&
((__builtin_constant_p (lc->names[i]) && strlen (
lc->names[i]) < ((size_t) (name_length))) || (__builtin_constant_p
(name) && strlen (name) < ((size_t) (name_length)
))) ? __extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(lc->names[i]) && __builtin_constant_p (name) &&
(__s1_len = strlen (lc->names[i]), __s2_len = strlen (name
), (!((size_t)(const void *)((lc->names[i]) + 1) - (size_t
)(const void *)(lc->names[i]) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)((name) + 1) - (size_t)(const void *
)(name) == 1) || __s2_len >= 4)) ? __builtin_strcmp (lc->
names[i], name) : (__builtin_constant_p (lc->names[i]) &&
((size_t)(const void *)((lc->names[i]) + 1) - (size_t)(const
void *)(lc->names[i]) == 1) && (__s1_len = strlen
(lc->names[i]), __s1_len < 4) ? (__builtin_constant_p (
name) && ((size_t)(const void *)((name) + 1) - (size_t
)(const void *)(name) == 1) ? __builtin_strcmp (lc->names[
i], name) : (__extension__ ({ const unsigned char *__s2 = (const
unsigned char *) (const char *) (name); int __result = (((const
unsigned char *) (const char *) (lc->names[i]))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (lc->names[i])
)[1] - __s2[1]); if (__s1_len > 1 && __result == 0
) { __result = (((const unsigned char *) (const char *) (lc->
names[i]))[2] - __s2[2]); if (__s1_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (lc
->names[i]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
(name) && ((size_t)(const void *)((name) + 1) - (size_t
)(const void *)(name) == 1) && (__s2_len = strlen (name
), __s2_len < 4) ? (__builtin_constant_p (lc->names[i])
&& ((size_t)(const void *)((lc->names[i]) + 1) - (
size_t)(const void *)(lc->names[i]) == 1) ? __builtin_strcmp
(lc->names[i], name) : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (lc->
names[i]); int __result = (((const unsigned char *) (const char
*) (name))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
name))[1] - __s2[1]); if (__s2_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
name))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (name
))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (lc->
names[i], name)))); }) : strncmp (lc->names[i], name, name_length
))) == 0)
287 {
288 return i;
289 }
290 }
291 return -1;
292}
293
294
295
296/* initialize a lambda vector */
297static void lambda_vec_init(lambda_vec_t *lv, const lambda_components_t *lc)
298{
299 snew(lv->val, lc->N)(lv->val) = save_calloc("lv->val", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 299, (lc->N), sizeof(*(lv->val)));
300 lv->index = -1;
301 lv->dhdl = -1;
302 lv->lc = lc;
303}
304
305static void lambda_vec_destroy(lambda_vec_t *lv)
306{
307 sfree(lv->val)save_free("lv->val", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 307, (lv->val));
308}
309
310static void lambda_vec_copy(lambda_vec_t *lv, const lambda_vec_t *orig)
311{
312 int i;
313
314 lambda_vec_init(lv, orig->lc);
315 lv->dhdl = orig->dhdl;
316 lv->index = orig->index;
317 for (i = 0; i < lv->lc->N; i++)
318 {
319 lv->val[i] = orig->val[i];
320 }
321}
322
323/* write a lambda vec to a preallocated string */
324static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
325{
326 int i;
327 size_t np;
328
329 str[0] = 0; /* reset the string */
330 if (lv->dhdl < 0)
331 {
332 if (named)
333 {
334 str += sprintf(str, "delta H to ");
335 }
336 if (lv->lc->N > 1)
337 {
338 str += sprintf(str, "(");
339 }
340 for (i = 0; i < lv->lc->N; i++)
341 {
342 str += sprintf(str, "%g", lv->val[i]);
343 if (i < lv->lc->N-1)
344 {
345 str += sprintf(str, ", ");
346 }
347 }
348 if (lv->lc->N > 1)
349 {
350 str += sprintf(str, ")");
351 }
352 }
353 else
354 {
355 /* this lambda vector describes a derivative */
356 str += sprintf(str, "dH/dl");
357 if (strlen(lv->lc->names[lv->dhdl]) > 0)
358 {
359 str += sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
360 }
361 }
362}
363
364/* write a shortened version of the lambda vec to a preallocated string */
365static void lambda_vec_print_short(const lambda_vec_t *lv, char *str)
366{
367 int i;
368 size_t np;
369
370 if (lv->index >= 0)
371 {
372 sprintf(str, "%6d", lv->index);
373 }
374 else
375 {
376 if (lv->dhdl < 0)
377 {
378 sprintf(str, "%6.3f", lv->val[0]);
379 }
380 else
381 {
382 sprintf(str, "dH/dl[%d]", lv->dhdl);
383 }
384 }
385}
386
387/* write an intermediate version of two lambda vecs to a preallocated string */
388static void lambda_vec_print_intermediate(const lambda_vec_t *a,
389 const lambda_vec_t *b, char *str)
390{
391 int i;
392 size_t np;
393
394 str[0] = 0;
395 if ( (a->index >= 0) && (b->index >= 0) )
396 {
397 sprintf(str, "%6.3f", ((double)a->index+(double)b->index)/2.);
398 }
399 else
400 {
401 if ( (a->dhdl < 0) && (b->dhdl < 0) )
402 {
403 sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.);
404 }
405 }
406}
407
408
409
410/* calculate the difference in lambda vectors: c = a-b.
411 c must be initialized already, and a and b must describe non-derivative
412 lambda points */
413static void lambda_vec_diff(const lambda_vec_t *a, const lambda_vec_t *b,
414 lambda_vec_t *c)
415{
416 int i;
417
418 if ( (a->dhdl > 0) || (b->dhdl > 0) )
419 {
420 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 420,
421 "Trying to calculate the difference between derivatives instead of lambda points");
422 }
423 if ((a->lc != b->lc) || (a->lc != c->lc) )
424 {
425 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 425,
426 "Trying to calculate the difference lambdas with differing basis set");
427 }
428 for (i = 0; i < a->lc->N; i++)
429 {
430 c->val[i] = a->val[i] - b->val[i];
431 }
432}
433
434/* calculate and return the absolute difference in lambda vectors: c = |a-b|.
435 a and b must describe non-derivative lambda points */
436static double lambda_vec_abs_diff(const lambda_vec_t *a, const lambda_vec_t *b)
437{
438 int i;
439 double ret = 0.;
440
441 if ( (a->dhdl > 0) || (b->dhdl > 0) )
442 {
443 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 443,
444 "Trying to calculate the difference between derivatives instead of lambda points");
445 }
446 if (a->lc != b->lc)
447 {
448 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 448,
449 "Trying to calculate the difference lambdas with differing basis set");
450 }
451 for (i = 0; i < a->lc->N; i++)
452 {
453 double df = a->val[i] - b->val[i];
454 ret += df*df;
455 }
456 return sqrt(ret);
457}
458
459
460/* check whether two lambda vectors are the same */
461static gmx_bool lambda_vec_same(const lambda_vec_t *a, const lambda_vec_t *b)
462{
463 int i;
464
465 if (a->lc != b->lc)
466 {
467 return FALSE0;
468 }
469 if (a->dhdl < 0)
470 {
471 for (i = 0; i < a->lc->N; i++)
472 {
473 if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS5.96046448E-08))
474 {
475 return FALSE0;
476 }
477 }
478 return TRUE1;
479 }
480 else
481 {
482 /* they're derivatives, so we check whether the indices match */
483 return (a->dhdl == b->dhdl);
484 }
485}
486
487/* Compare the sort order of two foreign lambda vectors
488
489 returns 1 if a is 'bigger' than b,
490 returns 0 if they're the same,
491 returns -1 if a is 'smaller' than b.*/
492static gmx_bool lambda_vec_cmp_foreign(const lambda_vec_t *a,
493 const lambda_vec_t *b)
494{
495 int i;
496 double norm_a = 0, norm_b = 0;
497 gmx_bool different = FALSE0;
498
499 if (a->lc != b->lc)
500 {
501 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 501, "Can't compare lambdas with differing basis sets");
502 }
503 /* if either one has an index we sort based on that */
504 if ((a->index >= 0) || (b->index >= 0))
505 {
506 if (a->index == b->index)
507 {
508 return 0;
509 }
510 return (a->index > b->index) ? 1 : -1;
511 }
512 if (a->dhdl >= 0 || b->dhdl >= 0)
513 {
514 /* lambda vectors that are derivatives always sort higher than those
515 without derivatives */
516 if ((a->dhdl >= 0) != (b->dhdl >= 0) )
517 {
518 return (a->dhdl >= 0) ? 1 : -1;
519 }
520 return a->dhdl > b->dhdl;
521 }
522
523 /* neither has an index, so we can only sort on the lambda components,
524 which is only valid if there is one component */
525 for (i = 0; i < a->lc->N; i++)
526 {
527 if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS5.96046448E-08))
528 {
529 different = TRUE1;
530 }
531 norm_a += a->val[i]*a->val[i];
532 norm_b += b->val[i]*b->val[i];
533 }
534 if (!different)
535 {
536 return 0;
537 }
538 return norm_a > norm_b;
539}
540
541/* Compare the sort order of two native lambda vectors
542
543 returns 1 if a is 'bigger' than b,
544 returns 0 if they're the same,
545 returns -1 if a is 'smaller' than b.*/
546static gmx_bool lambda_vec_cmp_native(const lambda_vec_t *a,
547 const lambda_vec_t *b)
548{
549 int i;
550
551 if (a->lc != b->lc)
552 {
553 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 553, "Can't compare lambdas with differing basis sets");
554 }
555 /* if either one has an index we sort based on that */
556 if ((a->index >= 0) || (b->index >= 0))
557 {
558 if (a->index == b->index)
559 {
560 return 0;
561 }
562 return (a->index > b->index) ? 1 : -1;
563 }
564 /* neither has an index, so we can only sort on the lambda components,
565 which is only valid if there is one component */
566 if (a->lc->N > 1)
567 {
568 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 568,
569 "Can't compare lambdas with no index and > 1 component");
570 }
571 if (a->dhdl >= 0 || b->dhdl >= 0)
572 {
573 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 573,
574 "Can't compare native lambdas that are derivatives");
575 }
576 if (gmx_within_tol(a->val[0], b->val[0], 10*GMX_REAL_EPS5.96046448E-08))
577 {
578 return 0;
579 }
580 return a->val[0] > b->val[0] ? 1 : -1;
581}
582
583
584
585
586static void hist_init(hist_t *h, int nhist, int *nbin)
587{
588 int i;
589 if (nhist > 2)
590 {
591 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 591, "histogram with more than two sets of data!");
592 }
593 for (i = 0; i < nhist; i++)
594 {
595 snew(h->bin[i], nbin[i])(h->bin[i]) = save_calloc("h->bin[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 595, (nbin[i]), sizeof(*(h->bin[i])));
596 h->x0[i] = 0;
597 h->nbin[i] = nbin[i];
598 h->start_time = h->delta_time = 0;
599 h->dx[i] = 0;
600 }
601 h->sum = 0;
602 h->nhist = nhist;
603}
604
605static void hist_destroy(hist_t *h)
606{
607 sfree(h->bin)save_free("h->bin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 607, (h->bin));
608}
609
610
611static void xvg_init(xvg_t *ba)
612{
613 ba->filename = NULL((void*)0);
614 ba->nset = 0;
615 ba->np_alloc = 0;
616 ba->np = NULL((void*)0);
617 ba->y = NULL((void*)0);
618}
619
620static void samples_init(samples_t *s, lambda_vec_t *native_lambda,
621 lambda_vec_t *foreign_lambda, double temp,
622 gmx_bool derivative, const char *filename)
623{
624 s->native_lambda = native_lambda;
625 s->foreign_lambda = foreign_lambda;
626 s->temp = temp;
627 s->derivative = derivative;
628
629 s->ndu = 0;
630 s->du = NULL((void*)0);
631 s->t = NULL((void*)0);
632 s->start_time = s->delta_time = 0;
633 s->hist = NULL((void*)0);
634 s->du_alloc = NULL((void*)0);
635 s->t_alloc = NULL((void*)0);
636 s->hist_alloc = NULL((void*)0);
637 s->ndu_alloc = 0;
638 s->nt_alloc = 0;
639
640 s->ntot = 0;
641 s->filename = filename;
642}
643
644static void sample_range_init(sample_range_t *r, samples_t *s)
645{
646 r->start = 0;
647 r->end = s->ndu;
648 r->use = TRUE1;
649 r->s = NULL((void*)0);
650}
651
652static void sample_coll_init(sample_coll_t *sc, lambda_vec_t *native_lambda,
653 lambda_vec_t *foreign_lambda, double temp)
654{
655 sc->native_lambda = native_lambda;
656 sc->foreign_lambda = foreign_lambda;
657 sc->temp = temp;
658
659 sc->nsamples = 0;
660 sc->s = NULL((void*)0);
661 sc->r = NULL((void*)0);
662 sc->nsamples_alloc = 0;
663
664 sc->ntot = 0;
665 sc->next = sc->prev = NULL((void*)0);
666}
667
668static void sample_coll_destroy(sample_coll_t *sc)
669{
670 /* don't free the samples themselves */
671 sfree(sc->r)save_free("sc->r", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 671, (sc->r));
672 sfree(sc->s)save_free("sc->s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 672, (sc->s));
673}
674
675
676static void lambda_data_init(lambda_data_t *l, lambda_vec_t *native_lambda,
677 double temp)
678{
679 l->lambda = native_lambda;
680 l->temp = temp;
681
682 l->next = NULL((void*)0);
683 l->prev = NULL((void*)0);
684
685 l->sc = &(l->sc_head);
686
687 sample_coll_init(l->sc, native_lambda, NULL((void*)0), 0.);
688 l->sc->next = l->sc;
689 l->sc->prev = l->sc;
690}
691
692static void barres_init(barres_t *br)
693{
694 br->dg = 0;
695 br->dg_err = 0;
696 br->sa = 0;
697 br->sa_err = 0;
698 br->sb = 0;
699 br->sb_err = 0;
700 br->dg_stddev = 0;
701 br->dg_stddev_err = 0;
702
703 br->a = NULL((void*)0);
704 br->b = NULL((void*)0);
705}
706
707
708/* calculate the total number of samples in a sample collection */
709static void sample_coll_calc_ntot(sample_coll_t *sc)
710{
711 int i;
712
713 sc->ntot = 0;
714 for (i = 0; i < sc->nsamples; i++)
715 {
716 if (sc->r[i].use)
717 {
718 if (sc->s[i]->hist)
719 {
720 sc->ntot += sc->s[i]->ntot;
721 }
722 else
723 {
724 sc->ntot += sc->r[i].end - sc->r[i].start;
725 }
726 }
727 }
728}
729
730
731/* find the barsamples_t associated with a lambda that corresponds to
732 a specific foreign lambda */
733static sample_coll_t *lambda_data_find_sample_coll(lambda_data_t *l,
734 lambda_vec_t *foreign_lambda)
735{
736 sample_coll_t *sc = l->sc->next;
737
738 while (sc != l->sc)
739 {
740 if (lambda_vec_same(sc->foreign_lambda, foreign_lambda))
741 {
742 return sc;
743 }
744 sc = sc->next;
745 }
746
747 return NULL((void*)0);
748}
749
750/* insert li into an ordered list of lambda_colls */
751static void lambda_data_insert_sample_coll(lambda_data_t *l, sample_coll_t *sc)
752{
753 sample_coll_t *scn = l->sc->next;
754 while ( (scn != l->sc) )
755 {
756 if (lambda_vec_cmp_foreign(scn->foreign_lambda, sc->foreign_lambda) > 0)
757 {
758 break;
759 }
760 scn = scn->next;
761 }
762 /* now insert it before the found scn */
763 sc->next = scn;
764 sc->prev = scn->prev;
765 scn->prev->next = sc;
766 scn->prev = sc;
767}
768
769/* insert li into an ordered list of lambdas */
770static void lambda_data_insert_lambda(lambda_data_t *head, lambda_data_t *li)
771{
772 lambda_data_t *lc = head->next;
773 while (lc != head)
774 {
775 if (lambda_vec_cmp_native(lc->lambda, li->lambda) > 0)
776 {
777 break;
778 }
779 lc = lc->next;
780 }
781 /* now insert ourselves before the found lc */
782 li->next = lc;
783 li->prev = lc->prev;
784 lc->prev->next = li;
785 lc->prev = li;
786}
787
788/* insert a sample and a sample_range into a sample_coll. The
789 samples are stored as a pointer, the range is copied. */
790static void sample_coll_insert_sample(sample_coll_t *sc, samples_t *s,
791 sample_range_t *r)
792{
793 /* first check if it belongs here */
794 if (sc->temp != s->temp)
795 {
796 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 796, "Temperatures in files %s and %s are not the same!",
797 s->filename, sc->next->s[0]->filename);
798 }
799 if (!lambda_vec_same(sc->native_lambda, s->native_lambda))
800 {
801 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 801, "Native lambda in files %s and %s are not the same (and they should be)!",
802 s->filename, sc->next->s[0]->filename);
803 }
804 if (!lambda_vec_same(sc->foreign_lambda, s->foreign_lambda))
805 {
806 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 806, "Foreign lambda in files %s and %s are not the same (and they should be)!",
807 s->filename, sc->next->s[0]->filename);
808 }
809
810 /* check if there's room */
811 if ( (sc->nsamples + 1) > sc->nsamples_alloc)
812 {
813 sc->nsamples_alloc = max(2*sc->nsamples_alloc, 2)(((2*sc->nsamples_alloc) > (2)) ? (2*sc->nsamples_alloc
) : (2) );
814 srenew(sc->s, sc->nsamples_alloc)(sc->s) = save_realloc("sc->s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 814, (sc->s), (sc->nsamples_alloc), sizeof(*(sc->s
)));
815 srenew(sc->r, sc->nsamples_alloc)(sc->r) = save_realloc("sc->r", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 815, (sc->r), (sc->nsamples_alloc), sizeof(*(sc->r
)));
816 }
817 sc->s[sc->nsamples] = s;
818 sc->r[sc->nsamples] = *r;
819 sc->nsamples++;
820
821 sample_coll_calc_ntot(sc);
822}
823
824/* insert a sample into a lambda_list, creating the right sample_coll if
825 neccesary */
826static void lambda_data_list_insert_sample(lambda_data_t *head, samples_t *s)
827{
828 gmx_bool found = FALSE0;
829 sample_coll_t *sc;
830 sample_range_t r;
831
832 lambda_data_t *l = head->next;
833
834 /* first search for the right lambda_data_t */
835 while (l != head)
836 {
837 if (lambda_vec_same(l->lambda, s->native_lambda) )
838 {
839 found = TRUE1;
840 break;
841 }
842 l = l->next;
843 }
844
845 if (!found)
846 {
847 snew(l, 1)(l) = save_calloc("l", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 847, (1), sizeof(*(l))); /* allocate a new one */
848 lambda_data_init(l, s->native_lambda, s->temp); /* initialize it */
849 lambda_data_insert_lambda(head, l); /* add it to the list */
850 }
851
852 /* now look for a sample collection */
853 sc = lambda_data_find_sample_coll(l, s->foreign_lambda);
854 if (!sc)
855 {
856 snew(sc, 1)(sc) = save_calloc("sc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 856, (1), sizeof(*(sc))); /* allocate a new one */
857 sample_coll_init(sc, s->native_lambda, s->foreign_lambda, s->temp);
858 lambda_data_insert_sample_coll(l, sc);
859 }
860
861 /* now insert the samples into the sample coll */
862 sample_range_init(&r, s);
863 sample_coll_insert_sample(sc, s, &r);
864}
865
866
867/* make a histogram out of a sample collection */
868static void sample_coll_make_hist(sample_coll_t *sc, int **bin,
869 int *nbin_alloc, int *nbin,
870 double *dx, double *xmin, int nbin_default)
871{
872 int i, j, k;
873 gmx_bool dx_set = FALSE0;
874 gmx_bool xmin_set = FALSE0;
875
876 gmx_bool xmax_set = FALSE0;
877 gmx_bool xmax_set_hard = FALSE0; /* whether the xmax is bounded by the
878 limits of a histogram */
879 double xmax = -1;
880
881 /* first determine dx and xmin; try the histograms */
882 for (i = 0; i < sc->nsamples; i++)
883 {
884 if (sc->s[i]->hist)
885 {
886 hist_t *hist = sc->s[i]->hist;
887 for (k = 0; k < hist->nhist; k++)
888 {
889 double hdx = hist->dx[k];
890 double xmax_now = (hist->x0[k]+hist->nbin[k])*hdx;
891
892 /* we use the biggest dx*/
893 if ( (!dx_set) || hist->dx[0] > *dx)
894 {
895 dx_set = TRUE1;
896 *dx = hist->dx[0];
897 }
898 if ( (!xmin_set) || (hist->x0[k]*hdx) < *xmin)
899 {
900 xmin_set = TRUE1;
901 *xmin = (hist->x0[k]*hdx);
902 }
903
904 if ( (!xmax_set) || (xmax_now > xmax && !xmax_set_hard) )
905 {
906 xmax_set = TRUE1;
907 xmax = xmax_now;
908 if (hist->bin[k][hist->nbin[k]-1] != 0)
909 {
910 xmax_set_hard = TRUE1;
911 }
912 }
913 if (hist->bin[k][hist->nbin[k]-1] != 0 && (xmax_now < xmax) )
914 {
915 xmax_set_hard = TRUE1;
916 xmax = xmax_now;
917 }
918 }
919 }
920 }
921 /* and the delta us */
922 for (i = 0; i < sc->nsamples; i++)
923 {
924 if (sc->s[i]->ndu > 0)
925 {
926 /* determine min and max */
927 int starti = sc->r[i].start;
928 int endi = sc->r[i].end;
929 double du_xmin = sc->s[i]->du[starti];
930 double du_xmax = sc->s[i]->du[starti];
931 for (j = starti+1; j < endi; j++)
932 {
933 if (sc->s[i]->du[j] < du_xmin)
934 {
935 du_xmin = sc->s[i]->du[j];
936 }
937 if (sc->s[i]->du[j] > du_xmax)
938 {
939 du_xmax = sc->s[i]->du[j];
940 }
941 }
942
943 /* and now change the limits */
944 if ( (!xmin_set) || (du_xmin < *xmin) )
945 {
946 xmin_set = TRUE1;
947 *xmin = du_xmin;
948 }
949 if ( (!xmax_set) || ((du_xmax > xmax) && !xmax_set_hard) )
950 {
951 xmax_set = TRUE1;
952 xmax = du_xmax;
953 }
954 }
955 }
956
957 if (!xmax_set || !xmin_set)
958 {
959 *nbin = 0;
960 return;
961 }
962
963
964 if (!dx_set)
965 {
966 *nbin = nbin_default;
967 *dx = (xmax-(*xmin))/((*nbin)-2); /* -2 because we want the last bin to
968 be 0, and we count from 0 */
969 }
970 else
971 {
972 *nbin = (xmax-(*xmin))/(*dx);
973 }
974
975 if (*nbin > *nbin_alloc)
976 {
977 *nbin_alloc = *nbin;
978 srenew(*bin, *nbin_alloc)(*bin) = save_realloc("*bin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 978, (*bin), (*nbin_alloc), sizeof(*(*bin)));
979 }
980
981 /* reset the histogram */
982 for (i = 0; i < (*nbin); i++)
983 {
984 (*bin)[i] = 0;
985 }
986
987 /* now add the actual data */
988 for (i = 0; i < sc->nsamples; i++)
989 {
990 if (sc->s[i]->hist)
991 {
992 hist_t *hist = sc->s[i]->hist;
993 for (k = 0; k < hist->nhist; k++)
994 {
995 double hdx = hist->dx[k];
996 double xmin_hist = hist->x0[k]*hdx;
997 for (j = 0; j < hist->nbin[k]; j++)
998 {
999 /* calculate the bin corresponding to the middle of the
1000 original bin */
1001 double x = hdx*(j+0.5) + xmin_hist;
1002 int binnr = (int)((x-(*xmin))/(*dx));
1003
1004 if (binnr >= *nbin || binnr < 0)
1005 {
1006 binnr = (*nbin)-1;
1007 }
1008
1009 (*bin)[binnr] += hist->bin[k][j];
1010 }
1011 }
1012 }
1013 else
1014 {
1015 int starti = sc->r[i].start;
1016 int endi = sc->r[i].end;
1017 for (j = starti; j < endi; j++)
1018 {
1019 int binnr = (int)((sc->s[i]->du[j]-(*xmin))/(*dx));
1020 if (binnr >= *nbin || binnr < 0)
1021 {
1022 binnr = (*nbin)-1;
1023 }
1024
1025 (*bin)[binnr]++;
1026 }
1027 }
1028 }
1029}
1030
1031/* write a collection of histograms to a file */
1032void sim_data_histogram(sim_data_t *sd, const char *filename,
1033 int nbin_default, const output_env_t oenv)
1034{
1035 char label_x[STRLEN4096];
1036 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1037 const char *title = "N(\\DeltaH)";
1038 const char *label_y = "Samples";
1039 FILE *fp;
1040 lambda_data_t *bl;
1041 int nsets = 0;
1042 char **setnames = NULL((void*)0);
1043 gmx_bool first_set = FALSE0;
1044 /* histogram data: */
1045 int *hist = NULL((void*)0);
1046 int nbin = 0;
1047 int nbin_alloc = 0;
1048 double dx = 0;
1049 double min = 0;
1050 int i;
1051 lambda_data_t *bl_head = sd->lb;
1052
1053 printf("\nWriting histogram to %s\n", filename);
1054 sprintf(label_x, "\\DeltaH (%s)", unit_energy"kJ/mol");
1055
1056 fp = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1057
1058 /* first get all the set names */
1059 bl = bl_head->next;
1060 /* iterate over all lambdas */
1061 while (bl != bl_head)
1062 {
1063 sample_coll_t *sc = bl->sc->next;
1064
1065 /* iterate over all samples */
1066 while (sc != bl->sc)
1067 {
1068 char buf[STRLEN4096], buf2[STRLEN4096];
1069
1070 nsets++;
1071 srenew(setnames, nsets)(setnames) = save_realloc("setnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1071, (setnames), (nsets), sizeof(*(setnames)));
1072 snew(setnames[nsets-1], STRLEN)(setnames[nsets-1]) = save_calloc("setnames[nsets-1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1072, (4096), sizeof(*(setnames[nsets-1])));
1073 if (sc->foreign_lambda->dhdl < 0)
1074 {
1075 lambda_vec_print(sc->native_lambda, buf, FALSE0);
1076 lambda_vec_print(sc->foreign_lambda, buf2, FALSE0);
1077 sprintf(setnames[nsets-1], "N(%s(%s=%s) | %s=%s)",
1078 deltag, lambda, buf2, lambda, buf);
1079 }
1080 else
1081 {
1082 lambda_vec_print(sc->native_lambda, buf, FALSE0);
1083 sprintf(setnames[nsets-1], "N(%s | %s=%s)",
1084 dhdl, lambda, buf);
1085 }
1086 sc = sc->next;
1087 }
1088
1089 bl = bl->next;
1090 }
1091 xvgr_legend(fp, nsets, (const char**)setnames, oenv);
1092
1093
1094 /* now make the histograms */
1095 bl = bl_head->next;
1096 /* iterate over all lambdas */
1097 while (bl != bl_head)
1098 {
1099 sample_coll_t *sc = bl->sc->next;
1100
1101 /* iterate over all samples */
1102 while (sc != bl->sc)
1103 {
1104 if (!first_set)
1105 {
1106 xvgr_new_dataset(fp, 0, 0, NULL((void*)0), oenv);
1107 }
1108
1109 sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &min,
1110 nbin_default);
1111
1112 for (i = 0; i < nbin; i++)
1113 {
1114 double xmin = i*dx + min;
1115 double xmax = (i+1)*dx + min;
1116
1117 fprintf(fp, "%g %d\n%g %d\n", xmin, hist[i], xmax, hist[i]);
1118 }
1119
1120 first_set = FALSE0;
1121 sc = sc->next;
1122 }
1123
1124 bl = bl->next;
1125 }
1126
1127 if (hist)
1128 {
1129 sfree(hist)save_free("hist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1129, (hist));
1130 }
1131
1132 xvgrclose(fp);
1133}
1134
1135/* create a collection (array) of barres_t object given a ordered linked list
1136 of barlamda_t sample collections */
1137static barres_t *barres_list_create(sim_data_t *sd, int *nres,
1138 gmx_bool use_dhdl)
1139{
1140 lambda_data_t *bl;
1141 int nlambda = 0;
1142 barres_t *res;
1143 int i;
1144 gmx_bool dhdl = FALSE0;
1145 gmx_bool first = TRUE1;
1146 lambda_data_t *bl_head = sd->lb;
1147
1148 /* first count the lambdas */
1149 bl = bl_head->next;
1150 while (bl != bl_head)
1151 {
1152 nlambda++;
1153 bl = bl->next;
1154 }
1155 snew(res, nlambda-1)(res) = save_calloc("res", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1155, (nlambda-1), sizeof(*(res)));
1156
1157 /* next put the right samples in the res */
1158 *nres = 0;
1159 bl = bl_head->next->next; /* we start with the second one. */
1160 while (bl != bl_head)
1161 {
1162 sample_coll_t *sc, *scprev;
1163 barres_t *br = &(res[*nres]);
1164 /* there is always a previous one. we search for that as a foreign
1165 lambda: */
1166 scprev = lambda_data_find_sample_coll(bl->prev, bl->lambda);
1167 sc = lambda_data_find_sample_coll(bl, bl->prev->lambda);
1168
1169 barres_init(br);
1170
1171 if (use_dhdl)
1172 {
1173 /* we use dhdl */
1174
1175 scprev = lambda_data_find_sample_coll(bl->prev, bl->prev->lambda);
1176 sc = lambda_data_find_sample_coll(bl, bl->lambda);
1177
1178 if (first)
1179 {
1180 printf("\nWARNING: Using the derivative data (dH/dlambda) to extrapolate delta H values.\nThis will only work if the Hamiltonian is linear in lambda.\n");
1181 dhdl = TRUE1;
1182 }
1183 if (!dhdl)
1184 {
1185 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1185, "Some dhdl files contain only one value (dH/dl), while others \ncontain multiple values (dH/dl and/or Delta H), will not proceed \nbecause of possible inconsistencies.\n");
1186 }
1187 }
1188 else if (!scprev && !sc)
1189 {
1190 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1190, "There is no path from lambda=%f -> %f that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n", bl->prev->lambda, bl->lambda);
1191 }
1192
1193 /* normal delta H */
1194 if (!scprev)
1195 {
1196 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1196, "Could not find a set for foreign lambda = %f\nin the files for lambda = %f", bl->lambda, bl->prev->lambda);
1197 }
1198 if (!sc)
1199 {
1200 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1200, "Could not find a set for foreign lambda = %f\nin the files for lambda = %f", bl->prev->lambda, bl->lambda);
1201 }
1202 br->a = scprev;
1203 br->b = sc;
1204
1205 first = FALSE0;
1206 (*nres)++;
1207 bl = bl->next;
1208 }
1209 return res;
1210}
1211
1212/* estimate the maximum discretization error */
1213static double barres_list_max_disc_err(barres_t *res, int nres)
1214{
1215 int i, j;
1216 double disc_err = 0.;
1217 double delta_lambda;
1218
1219 for (i = 0; i < nres; i++)
1220 {
1221 barres_t *br = &(res[i]);
1222
1223 delta_lambda = lambda_vec_abs_diff(br->b->native_lambda,
1224 br->a->native_lambda);
1225
1226 for (j = 0; j < br->a->nsamples; j++)
1227 {
1228 if (br->a->s[j]->hist)
1229 {
1230 double Wfac = 1.;
1231 if (br->a->s[j]->derivative)
1232 {
1233 Wfac = delta_lambda;
1234 }
1235
1236 disc_err = max(disc_err, Wfac*br->a->s[j]->hist->dx[0])(((disc_err) > (Wfac*br->a->s[j]->hist->dx[0])
) ? (disc_err) : (Wfac*br->a->s[j]->hist->dx[0]) );
1237 }
1238 }
1239 for (j = 0; j < br->b->nsamples; j++)
1240 {
1241 if (br->b->s[j]->hist)
1242 {
1243 double Wfac = 1.;
1244 if (br->b->s[j]->derivative)
1245 {
1246 Wfac = delta_lambda;
1247 }
1248 disc_err = max(disc_err, Wfac*br->b->s[j]->hist->dx[0])(((disc_err) > (Wfac*br->b->s[j]->hist->dx[0])
) ? (disc_err) : (Wfac*br->b->s[j]->hist->dx[0]) );
1249 }
1250 }
1251 }
1252 return disc_err;
1253}
1254
1255
1256/* impose start and end times on a sample collection, updating sample_ranges */
1257static void sample_coll_impose_times(sample_coll_t *sc, double begin_t,
1258 double end_t)
1259{
1260 int i;
1261 for (i = 0; i < sc->nsamples; i++)
1262 {
1263 samples_t *s = sc->s[i];
1264 sample_range_t *r = &(sc->r[i]);
1265 if (s->hist)
1266 {
1267 double end_time = s->hist->delta_time*s->hist->sum +
1268 s->hist->start_time;
1269 if (s->hist->start_time < begin_t || end_time > end_t)
1270 {
1271 r->use = FALSE0;
1272 }
1273 }
1274 else
1275 {
1276 if (!s->t)
1277 {
1278 double end_time;
1279 if (s->start_time < begin_t)
1280 {
1281 r->start = (int)((begin_t - s->start_time)/s->delta_time);
1282 }
1283 end_time = s->delta_time*s->ndu + s->start_time;
1284 if (end_time > end_t)
1285 {
1286 r->end = (int)((end_t - s->start_time)/s->delta_time);
1287 }
1288 }
1289 else
1290 {
1291 int j;
1292 for (j = 0; j < s->ndu; j++)
1293 {
1294 if (s->t[j] < begin_t)
1295 {
1296 r->start = j;
1297 }
1298
1299 if (s->t[j] >= end_t)
1300 {
1301 r->end = j;
1302 break;
1303 }
1304 }
1305 }
1306 if (r->start > r->end)
1307 {
1308 r->use = FALSE0;
1309 }
1310 }
1311 }
1312 sample_coll_calc_ntot(sc);
1313}
1314
1315static void sim_data_impose_times(sim_data_t *sd, double begin, double end)
1316{
1317 double first_t, last_t;
1318 double begin_t, end_t;
1319 lambda_data_t *lc;
1320 lambda_data_t *head = sd->lb;
1321 int j;
1322
1323 if (begin <= 0 && end < 0)
1324 {
1325 return;
1326 }
1327
1328 /* first determine the global start and end times */
1329 first_t = -1;
1330 last_t = -1;
1331 lc = head->next;
1332 while (lc != head)
1333 {
1334 sample_coll_t *sc = lc->sc->next;
1335 while (sc != lc->sc)
1336 {
1337 for (j = 0; j < sc->nsamples; j++)
1338 {
1339 double start_t, end_t;
1340
1341 start_t = sc->s[j]->start_time;
1342 end_t = sc->s[j]->start_time;
1343 if (sc->s[j]->hist)
1344 {
1345 end_t += sc->s[j]->delta_time*sc->s[j]->hist->sum;
1346 }
1347 else
1348 {
1349 if (sc->s[j]->t)
1350 {
1351 end_t = sc->s[j]->t[sc->s[j]->ndu-1];
1352 }
1353 else
1354 {
1355 end_t += sc->s[j]->delta_time*sc->s[j]->ndu;
1356 }
1357 }
1358
1359 if (start_t < first_t || first_t < 0)
1360 {
1361 first_t = start_t;
1362 }
1363 if (end_t > last_t)
1364 {
1365 last_t = end_t;
1366 }
1367 }
1368 sc = sc->next;
1369 }
1370 lc = lc->next;
1371 }
1372
1373 /* calculate the actual times */
1374 if (begin > 0)
1375 {
1376 begin_t = begin;
1377 }
1378 else
1379 {
1380 begin_t = first_t;
1381 }
1382
1383 if (end > 0)
1384 {
1385 end_t = end;
1386 }
1387 else
1388 {
1389 end_t = last_t;
1390 }
1391 printf("\n Samples in time interval: %.3f - %.3f\n", first_t, last_t);
1392
1393 if (begin_t > end_t)
1394 {
1395 return;
1396 }
1397 printf("Removing samples outside of: %.3f - %.3f\n", begin_t, end_t);
1398
1399 /* then impose them */
1400 lc = head->next;
1401 while (lc != head)
1402 {
1403 sample_coll_t *sc = lc->sc->next;
1404 while (sc != lc->sc)
1405 {
1406 sample_coll_impose_times(sc, begin_t, end_t);
1407 sc = sc->next;
1408 }
1409 lc = lc->next;
1410 }
1411}
1412
1413
1414/* create subsample i out of ni from an existing sample_coll */
1415static gmx_bool sample_coll_create_subsample(sample_coll_t *sc,
1416 sample_coll_t *sc_orig,
1417 int i, int ni)
1418{
1419 int j;
1420 int hist_start, hist_end;
1421
1422 gmx_int64_t ntot_start;
1423 gmx_int64_t ntot_end;
1424 gmx_int64_t ntot_so_far;
1425
1426 *sc = *sc_orig; /* just copy all fields */
1427
1428 /* allocate proprietary memory */
1429 snew(sc->s, sc_orig->nsamples)(sc->s) = save_calloc("sc->s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1429, (sc_orig->nsamples), sizeof(*(sc->s)));
1430 snew(sc->r, sc_orig->nsamples)(sc->r) = save_calloc("sc->r", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1430, (sc_orig->nsamples), sizeof(*(sc->r)));
1431
1432 /* copy the samples */
1433 for (j = 0; j < sc_orig->nsamples; j++)
1434 {
1435 sc->s[j] = sc_orig->s[j];
1436 sc->r[j] = sc_orig->r[j]; /* copy the ranges too */
1437 }
1438
1439 /* now fix start and end fields */
1440 /* the casts avoid possible overflows */
1441 ntot_start = (gmx_int64_t)(sc_orig->ntot*(double)i/(double)ni);
1442 ntot_end = (gmx_int64_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
1443 ntot_so_far = 0;
1444 for (j = 0; j < sc->nsamples; j++)
1445 {
1446 gmx_int64_t ntot_add;
1447 gmx_int64_t new_start, new_end;
1448
1449 if (sc->r[j].use)
1450 {
1451 if (sc->s[j]->hist)
1452 {
1453 ntot_add = sc->s[j]->hist->sum;
1454 }
1455 else
1456 {
1457 ntot_add = sc->r[j].end - sc->r[j].start;
1458 }
1459 }
1460 else
1461 {
1462 ntot_add = 0;
1463 }
1464
1465 if (!sc->s[j]->hist)
1466 {
1467 if (ntot_so_far < ntot_start)
1468 {
1469 /* adjust starting point */
1470 new_start = sc->r[j].start + (ntot_start - ntot_so_far);
1471 }
1472 else
1473 {
1474 new_start = sc->r[j].start;
1475 }
1476 /* adjust end point */
1477 new_end = sc->r[j].start + (ntot_end - ntot_so_far);
1478 if (new_end > sc->r[j].end)
1479 {
1480 new_end = sc->r[j].end;
1481 }
1482
1483 /* check if we're in range at all */
1484 if ( (new_end < new_start) || (new_start > sc->r[j].end) )
1485 {
1486 new_start = 0;
1487 new_end = 0;
1488 }
1489 /* and write the new range */
1490 sc->r[j].start = (int)new_start;
1491 sc->r[j].end = (int)new_end;
1492 }
1493 else
1494 {
1495 if (sc->r[j].use)
1496 {
1497 double overlap;
1498 double ntot_start_norm, ntot_end_norm;
1499 /* calculate the amount of overlap of the
1500 desired range (ntot_start -- ntot_end) onto
1501 the histogram range (ntot_so_far -- ntot_so_far+ntot_add)*/
1502
1503 /* first calculate normalized bounds
1504 (where 0 is the start of the hist range, and 1 the end) */
1505 ntot_start_norm = (ntot_start-ntot_so_far)/(double)ntot_add;
1506 ntot_end_norm = (ntot_end-ntot_so_far)/(double)ntot_add;
1507
1508 /* now fix the boundaries */
1509 ntot_start_norm = min(1, max(0., ntot_start_norm))(((1) < ((((0.) > (ntot_start_norm)) ? (0.) : (ntot_start_norm
) ))) ? (1) : ((((0.) > (ntot_start_norm)) ? (0.) : (ntot_start_norm
) )) );
1510 ntot_end_norm = max(0, min(1., ntot_end_norm))(((0) > ((((1.) < (ntot_end_norm)) ? (1.) : (ntot_end_norm
) ))) ? (0) : ((((1.) < (ntot_end_norm)) ? (1.) : (ntot_end_norm
) )) );
1511
1512 /* and calculate the overlap */
1513 overlap = ntot_end_norm - ntot_start_norm;
1514
1515 if (overlap > 0.95) /* we allow for 5% slack */
1516 {
1517 sc->r[j].use = TRUE1;
1518 }
1519 else if (overlap < 0.05)
1520 {
1521 sc->r[j].use = FALSE0;
1522 }
1523 else
1524 {
1525 return FALSE0;
1526 }
1527 }
1528 }
1529 ntot_so_far += ntot_add;
1530 }
1531 sample_coll_calc_ntot(sc);
1532
1533 return TRUE1;
1534}
1535
1536/* calculate minimum and maximum work values in sample collection */
1537static void sample_coll_min_max(sample_coll_t *sc, double Wfac,
1538 double *Wmin, double *Wmax)
1539{
1540 int i, j;
1541
1542 *Wmin = FLT_MAX3.40282347e+38F;
1543 *Wmax = -FLT_MAX3.40282347e+38F;
1544
1545 for (i = 0; i < sc->nsamples; i++)
1546 {
1547 samples_t *s = sc->s[i];
1548 sample_range_t *r = &(sc->r[i]);
1549 if (r->use)
1550 {
1551 if (!s->hist)
1552 {
1553 for (j = r->start; j < r->end; j++)
1554 {
1555 *Wmin = min(*Wmin, s->du[j]*Wfac)(((*Wmin) < (s->du[j]*Wfac)) ? (*Wmin) : (s->du[j]*Wfac
) );
1556 *Wmax = max(*Wmax, s->du[j]*Wfac)(((*Wmax) > (s->du[j]*Wfac)) ? (*Wmax) : (s->du[j]*Wfac
) );
1557 }
1558 }
1559 else
1560 {
1561 int hd = 0; /* determine the histogram direction: */
1562 double dx;
1563 if ( (s->hist->nhist > 1) && (Wfac < 0) )
1564 {
1565 hd = 1;
1566 }
1567 dx = s->hist->dx[hd];
1568
1569 for (j = s->hist->nbin[hd]-1; j >= 0; j--)
1570 {
1571 *Wmin = min(*Wmin, Wfac*(s->hist->x0[hd])*dx)(((*Wmin) < (Wfac*(s->hist->x0[hd])*dx)) ? (*Wmin) :
(Wfac*(s->hist->x0[hd])*dx) );
1572 *Wmax = max(*Wmax, Wfac*(s->hist->x0[hd])*dx)(((*Wmax) > (Wfac*(s->hist->x0[hd])*dx)) ? (*Wmax) :
(Wfac*(s->hist->x0[hd])*dx) );
1573 /* look for the highest value bin with values */
1574 if (s->hist->bin[hd][j] > 0)
1575 {
1576 *Wmin = min(*Wmin, Wfac*(j+s->hist->x0[hd]+1)*dx)(((*Wmin) < (Wfac*(j+s->hist->x0[hd]+1)*dx)) ? (*Wmin
) : (Wfac*(j+s->hist->x0[hd]+1)*dx) );
1577 *Wmax = max(*Wmax, Wfac*(j+s->hist->x0[hd]+1)*dx)(((*Wmax) > (Wfac*(j+s->hist->x0[hd]+1)*dx)) ? (*Wmax
) : (Wfac*(j+s->hist->x0[hd]+1)*dx) );
1578 break;
1579 }
1580 }
1581 }
1582 }
1583 }
1584}
1585
1586/* Initialize a sim_data structure */
1587static void sim_data_init(sim_data_t *sd)
1588{
1589 /* make linked list */
1590 sd->lb = &(sd->lb_head);
1591 sd->lb->next = sd->lb;
1592 sd->lb->prev = sd->lb;
1593
1594 lambda_components_init(&(sd->lc));
1595}
1596
1597
1598static double calc_bar_sum(int n, const double *W, double Wfac, double sbMmDG)
1599{
1600 int i;
1601 double sum;
1602
1603 sum = 0;
1604
1605 for (i = 0; i < n; i++)
1606 {
1607 sum += 1./(1. + exp(Wfac*W[i] + sbMmDG));
1608 }
1609
1610 return sum;
1611}
1612
1613/* calculate the BAR average given a histogram
1614
1615 if type== 0, calculate the best estimate for the average,
1616 if type==-1, calculate the minimum possible value given the histogram
1617 if type== 1, calculate the maximum possible value given the histogram */
1618static double calc_bar_sum_hist(const hist_t *hist, double Wfac, double sbMmDG,
1619 int type)
1620{
1621 double sum = 0.;
1622 int i;
1623 int max;
1624 /* normalization factor multiplied with bin width and
1625 number of samples (we normalize through M): */
1626 double normdx = 1.;
1627 int hd = 0; /* determine the histogram direction: */
1628 double dx;
1629
1630 if ( (hist->nhist > 1) && (Wfac < 0) )
1631 {
1632 hd = 1;
1633 }
1634 dx = hist->dx[hd];
1635 max = hist->nbin[hd]-1;
1636 if (type == 1)
1637 {
1638 max = hist->nbin[hd]; /* we also add whatever was out of range */
1639 }
1640
1641 for (i = 0; i < max; i++)
1642 {
1643 double x = Wfac*((i+hist->x0[hd])+0.5)*dx; /* bin middle */
1644 double pxdx = hist->bin[0][i]*normdx; /* p(x)dx */
1645
1646 sum += pxdx/(1. + exp(x + sbMmDG));
1647 }
1648
1649 return sum;
1650}
1651
1652static double calc_bar_lowlevel(sample_coll_t *ca, sample_coll_t *cb,
1653 double temp, double tol, int type)
1654{
1655 double kT, beta, M;
1656 double DG;
1657 int i, j;
1658 double Wfac1, Wfac2, Wmin, Wmax;
1659 double DG0, DG1, DG2, dDG1;
1660 double sum1, sum2;
1661 double n1, n2; /* numbers of samples as doubles */
1662
1663 kT = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp;
1664 beta = 1/kT;
1665
1666 /* count the numbers of samples */
1667 n1 = ca->ntot;
1668 n2 = cb->ntot;
1669
1670 M = log(n1/n2);
1671
1672 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1673 if (ca->foreign_lambda->dhdl < 0)
1674 {
1675 /* this is the case when the delta U were calculated directly
1676 (i.e. we're not scaling dhdl) */
1677 Wfac1 = beta;
1678 Wfac2 = beta;
1679 }
1680 else
1681 {
1682 /* we're using dhdl, so delta_lambda needs to be a
1683 multiplication factor. */
1684 /*double delta_lambda=cb->native_lambda-ca->native_lambda;*/
1685 double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1686 ca->native_lambda);
1687 if (cb->native_lambda->lc->N > 1)
1688 {
1689 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 1689,
1690 "Can't (yet) do multi-component dhdl interpolation");
1691 }
1692
1693 Wfac1 = beta*delta_lambda;
1694 Wfac2 = -beta*delta_lambda;
1695 }
1696
1697 if (beta < 1)
1698 {
1699 /* We print the output both in kT and kJ/mol.
1700 * Here we determine DG in kT, so when beta < 1
1701 * the precision has to be increased.
1702 */
1703 tol *= beta;
1704 }
1705
1706 /* Calculate minimum and maximum work to give an initial estimate of
1707 * delta G as their average.
1708 */
1709 {
1710 double Wmin1, Wmin2, Wmax1, Wmax2;
1711 sample_coll_min_max(ca, Wfac1, &Wmin1, &Wmax1);
1712 sample_coll_min_max(cb, Wfac2, &Wmin2, &Wmax2);
1713
1714 Wmin = min(Wmin1, Wmin2)(((Wmin1) < (Wmin2)) ? (Wmin1) : (Wmin2) );
1715 Wmax = max(Wmax1, Wmax2)(((Wmax1) > (Wmax2)) ? (Wmax1) : (Wmax2) );
1716 }
1717
1718 DG0 = Wmin;
1719 DG2 = Wmax;
1720
1721 if (debug)
1722 {
1723 fprintf(debug, "DG %9.5f %9.5f\n", DG0, DG2);
1724 }
1725 /* We approximate by bisection: given our initial estimates
1726 we keep checking whether the halfway point is greater or
1727 smaller than what we get out of the BAR averages.
1728
1729 For the comparison we can use twice the tolerance. */
1730 while (DG2 - DG0 > 2*tol)
1731 {
1732 DG1 = 0.5*(DG0 + DG2);
1733
1734 /* calculate the BAR averages */
1735 dDG1 = 0.;
1736
1737 for (i = 0; i < ca->nsamples; i++)
1738 {
1739 samples_t *s = ca->s[i];
1740 sample_range_t *r = &(ca->r[i]);
1741 if (r->use)
1742 {
1743 if (s->hist)
1744 {
1745 dDG1 += calc_bar_sum_hist(s->hist, Wfac1, (M-DG1), type);
1746 }
1747 else
1748 {
1749 dDG1 += calc_bar_sum(r->end - r->start, s->du + r->start,
1750 Wfac1, (M-DG1));
1751 }
1752 }
1753 }
1754 for (i = 0; i < cb->nsamples; i++)
1755 {
1756 samples_t *s = cb->s[i];
1757 sample_range_t *r = &(cb->r[i]);
1758 if (r->use)
1759 {
1760 if (s->hist)
1761 {
1762 dDG1 -= calc_bar_sum_hist(s->hist, Wfac2, -(M-DG1), type);
1763 }
1764 else
1765 {
1766 dDG1 -= calc_bar_sum(r->end - r->start, s->du + r->start,
1767 Wfac2, -(M-DG1));
1768 }
1769 }
1770 }
1771
1772 if (dDG1 < 0)
1773 {
1774 DG0 = DG1;
1775 }
1776 else
1777 {
1778 DG2 = DG1;
1779 }
1780 if (debug)
1781 {
1782 fprintf(debug, "DG %9.5f %9.5f\n", DG0, DG2);
1783 }
1784 }
1785
1786 return 0.5*(DG0 + DG2);
1787}
1788
1789static void calc_rel_entropy(sample_coll_t *ca, sample_coll_t *cb,
1790 double temp, double dg, double *sa, double *sb)
1791{
1792 int i, j;
1793 double W_ab = 0.;
1794 double W_ba = 0.;
1795 double kT, beta;
1796 double Wfac1, Wfac2;
1797 double n1, n2;
1798
1799 kT = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp;
1800 beta = 1/kT;
1801
1802 /* count the numbers of samples */
1803 n1 = ca->ntot;
1804 n2 = cb->ntot;
1805
1806 /* to ensure the work values are the same as during the delta_G */
1807 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1808 if (ca->foreign_lambda->dhdl < 0)
1809 {
1810 /* this is the case when the delta U were calculated directly
1811 (i.e. we're not scaling dhdl) */
1812 Wfac1 = beta;
1813 Wfac2 = beta;
1814 }
1815 else
1816 {
1817 /* we're using dhdl, so delta_lambda needs to be a
1818 multiplication factor. */
1819 double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1820 ca->native_lambda);
1821 Wfac1 = beta*delta_lambda;
1822 Wfac2 = -beta*delta_lambda;
1823 }
1824
1825 /* first calculate the average work in both directions */
1826 for (i = 0; i < ca->nsamples; i++)
1827 {
1828 samples_t *s = ca->s[i];
1829 sample_range_t *r = &(ca->r[i]);
1830 if (r->use)
1831 {
1832 if (!s->hist)
1833 {
1834 for (j = r->start; j < r->end; j++)
1835 {
1836 W_ab += Wfac1*s->du[j];
1837 }
1838 }
1839 else
1840 {
1841 /* normalization factor multiplied with bin width and
1842 number of samples (we normalize through M): */
1843 double normdx = 1.;
1844 double dx;
1845 int hd = 0; /* histogram direction */
1846 if ( (s->hist->nhist > 1) && (Wfac1 < 0) )
1847 {
1848 hd = 1;
1849 }
1850 dx = s->hist->dx[hd];
1851
1852 for (j = 0; j < s->hist->nbin[0]; j++)
1853 {
1854 double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1855 double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
1856 W_ab += pxdx*x;
1857 }
1858 }
1859 }
1860 }
1861 W_ab /= n1;
1862
1863 for (i = 0; i < cb->nsamples; i++)
1864 {
1865 samples_t *s = cb->s[i];
1866 sample_range_t *r = &(cb->r[i]);
1867 if (r->use)
1868 {
1869 if (!s->hist)
1870 {
1871 for (j = r->start; j < r->end; j++)
1872 {
1873 W_ba += Wfac1*s->du[j];
1874 }
1875 }
1876 else
1877 {
1878 /* normalization factor multiplied with bin width and
1879 number of samples (we normalize through M): */
1880 double normdx = 1.;
1881 double dx;
1882 int hd = 0; /* histogram direction */
1883 if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
1884 {
1885 hd = 1;
1886 }
1887 dx = s->hist->dx[hd];
1888
1889 for (j = 0; j < s->hist->nbin[0]; j++)
1890 {
1891 double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1892 double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
1893 W_ba += pxdx*x;
1894 }
1895 }
1896 }
1897 }
1898 W_ba /= n2;
1899
1900 /* then calculate the relative entropies */
1901 *sa = (W_ab - dg);
1902 *sb = (W_ba + dg);
1903}
1904
1905static void calc_dg_stddev(sample_coll_t *ca, sample_coll_t *cb,
1906 double temp, double dg, double *stddev)
1907{
1908 int i, j;
1909 double M;
1910 double sigmafact = 0.;
1911 double kT, beta;
1912 double Wfac1, Wfac2;
1913 double n1, n2;
1914
1915 kT = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp;
1916 beta = 1/kT;
1917
1918 /* count the numbers of samples */
1919 n1 = ca->ntot;
1920 n2 = cb->ntot;
1921
1922 /* to ensure the work values are the same as during the delta_G */
1923 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1924 if (ca->foreign_lambda->dhdl < 0)
1925 {
1926 /* this is the case when the delta U were calculated directly
1927 (i.e. we're not scaling dhdl) */
1928 Wfac1 = beta;
1929 Wfac2 = beta;
1930 }
1931 else
1932 {
1933 /* we're using dhdl, so delta_lambda needs to be a
1934 multiplication factor. */
1935 double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1936 ca->native_lambda);
1937 Wfac1 = beta*delta_lambda;
1938 Wfac2 = -beta*delta_lambda;
1939 }
1940
1941 M = log(n1/n2);
1942
1943
1944 /* calculate average in both directions */
1945 for (i = 0; i < ca->nsamples; i++)
1946 {
1947 samples_t *s = ca->s[i];
1948 sample_range_t *r = &(ca->r[i]);
1949 if (r->use)
1950 {
1951 if (!s->hist)
1952 {
1953 for (j = r->start; j < r->end; j++)
1954 {
1955 sigmafact += 1./(2. + 2.*cosh((M + Wfac1*s->du[j] - dg)));
1956 }
1957 }
1958 else
1959 {
1960 /* normalization factor multiplied with bin width and
1961 number of samples (we normalize through M): */
1962 double normdx = 1.;
1963 double dx;
1964 int hd = 0; /* histogram direction */
1965 if ( (s->hist->nhist > 1) && (Wfac1 < 0) )
1966 {
1967 hd = 1;
1968 }
1969 dx = s->hist->dx[hd];
1970
1971 for (j = 0; j < s->hist->nbin[0]; j++)
1972 {
1973 double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1974 double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
1975
1976 sigmafact += pxdx/(2. + 2.*cosh((M + x - dg)));
1977 }
1978 }
1979 }
1980 }
1981 for (i = 0; i < cb->nsamples; i++)
1982 {
1983 samples_t *s = cb->s[i];
1984 sample_range_t *r = &(cb->r[i]);
1985 if (r->use)
1986 {
1987 if (!s->hist)
1988 {
1989 for (j = r->start; j < r->end; j++)
1990 {
1991 sigmafact += 1./(2. + 2.*cosh((M - Wfac2*s->du[j] - dg)));
1992 }
1993 }
1994 else
1995 {
1996 /* normalization factor multiplied with bin width and
1997 number of samples (we normalize through M): */
1998 double normdx = 1.;
1999 double dx;
2000 int hd = 0; /* histogram direction */
2001 if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
2002 {
2003 hd = 1;
2004 }
2005 dx = s->hist->dx[hd];
2006
2007 for (j = 0; j < s->hist->nbin[0]; j++)
2008 {
2009 double x = Wfac2*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
2010 double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
2011
2012 sigmafact += pxdx/(2. + 2.*cosh((M - x - dg)));
2013 }
2014 }
2015 }
2016 }
2017
2018 sigmafact /= (n1 + n2);
2019
2020
2021 /* Eq. 10 from
2022 Shirts, Bair, Hooker & Pande, Phys. Rev. Lett 91, 140601 (2003): */
2023 *stddev = sqrt(((1./sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
2024}
2025
2026
2027
2028static void calc_bar(barres_t *br, double tol,
2029 int npee_min, int npee_max, gmx_bool *bEE,
2030 double *partsum)
2031{
2032 int npee, p;
2033 double dg_sig2, sa_sig2, sb_sig2, stddev_sig2; /* intermediate variance values
2034 for calculated quantities */
2035 int nsample1, nsample2;
2036 double temp = br->a->temp;
2037 int i, j;
2038 double dg_min, dg_max;
2039 gmx_bool have_hist = FALSE0;
2040
2041 br->dg = calc_bar_lowlevel(br->a, br->b, temp, tol, 0);
2042
2043 br->dg_disc_err = 0.;
2044 br->dg_histrange_err = 0.;
2045
2046 /* check if there are histograms */
2047 for (i = 0; i < br->a->nsamples; i++)
2048 {
2049 if (br->a->r[i].use && br->a->s[i]->hist)
2050 {
2051 have_hist = TRUE1;
2052 break;
2053 }
2054 }
2055 if (!have_hist)
2056 {
2057 for (i = 0; i < br->b->nsamples; i++)
2058 {
2059 if (br->b->r[i].use && br->b->s[i]->hist)
2060 {
2061 have_hist = TRUE1;
2062 break;
2063 }
2064 }
2065 }
2066
2067 /* calculate histogram-specific errors */
2068 if (have_hist)
2069 {
2070 dg_min = calc_bar_lowlevel(br->a, br->b, temp, tol, -1);
2071 dg_max = calc_bar_lowlevel(br->a, br->b, temp, tol, 1);
2072
2073 if (fabs(dg_max - dg_min) > GMX_REAL_EPS5.96046448E-08*10)
2074 {
2075 /* the histogram range error is the biggest of the differences
2076 between the best estimate and the extremes */
2077 br->dg_histrange_err = fabs(dg_max - dg_min);
2078 }
2079 br->dg_disc_err = 0.;
2080 for (i = 0; i < br->a->nsamples; i++)
2081 {
2082 if (br->a->s[i]->hist)
2083 {
2084 br->dg_disc_err = max(br->dg_disc_err, br->a->s[i]->hist->dx[0])(((br->dg_disc_err) > (br->a->s[i]->hist->dx
[0])) ? (br->dg_disc_err) : (br->a->s[i]->hist->
dx[0]) );
2085 }
2086 }
2087 for (i = 0; i < br->b->nsamples; i++)
2088 {
2089 if (br->b->s[i]->hist)
2090 {
2091 br->dg_disc_err = max(br->dg_disc_err, br->b->s[i]->hist->dx[0])(((br->dg_disc_err) > (br->b->s[i]->hist->dx
[0])) ? (br->dg_disc_err) : (br->b->s[i]->hist->
dx[0]) );
2092 }
2093 }
2094 }
2095 calc_rel_entropy(br->a, br->b, temp, br->dg, &(br->sa), &(br->sb));
2096
2097 calc_dg_stddev(br->a, br->b, temp, br->dg, &(br->dg_stddev) );
2098
2099 dg_sig2 = 0;
2100 sa_sig2 = 0;
2101 sb_sig2 = 0;
2102 stddev_sig2 = 0;
2103
2104 *bEE = TRUE1;
2105 {
2106 sample_coll_t ca, cb;
2107
2108 /* initialize the samples */
2109 sample_coll_init(&ca, br->a->native_lambda, br->a->foreign_lambda,
2110 br->a->temp);
2111 sample_coll_init(&cb, br->b->native_lambda, br->b->foreign_lambda,
2112 br->b->temp);
2113
2114 for (npee = npee_min; npee <= npee_max; npee++)
2115 {
2116 double dgs = 0;
2117 double dgs2 = 0;
2118 double dsa = 0;
2119 double dsb = 0;
2120 double dsa2 = 0;
2121 double dsb2 = 0;
2122 double dstddev = 0;
2123 double dstddev2 = 0;
2124
2125
2126 for (p = 0; p < npee; p++)
2127 {
2128 double dgp;
2129 double stddevc;
2130 double sac, sbc;
2131 gmx_bool cac, cbc;
2132
2133 cac = sample_coll_create_subsample(&ca, br->a, p, npee);
2134 cbc = sample_coll_create_subsample(&cb, br->b, p, npee);
2135
2136 if (!cac || !cbc)
2137 {
2138 printf("WARNING: histogram number incompatible with block number for averaging: can't do error estimate\n");
2139 *bEE = FALSE0;
2140 if (cac)
2141 {
2142 sample_coll_destroy(&ca);
2143 }
2144 if (cbc)
2145 {
2146 sample_coll_destroy(&cb);
2147 }
2148 return;
2149 }
2150
2151 dgp = calc_bar_lowlevel(&ca, &cb, temp, tol, 0);
2152 dgs += dgp;
2153 dgs2 += dgp*dgp;
2154
2155 partsum[npee*(npee_max+1)+p] += dgp;
2156
2157 calc_rel_entropy(&ca, &cb, temp, dgp, &sac, &sbc);
2158 dsa += sac;
2159 dsa2 += sac*sac;
2160 dsb += sbc;
2161 dsb2 += sbc*sbc;
2162 calc_dg_stddev(&ca, &cb, temp, dgp, &stddevc );
2163
2164 dstddev += stddevc;
2165 dstddev2 += stddevc*stddevc;
2166
2167 sample_coll_destroy(&ca);
2168 sample_coll_destroy(&cb);
2169 }
2170 dgs /= npee;
2171 dgs2 /= npee;
2172 dg_sig2 += (dgs2-dgs*dgs)/(npee-1);
2173
2174 dsa /= npee;
2175 dsa2 /= npee;
2176 dsb /= npee;
2177 dsb2 /= npee;
2178 sa_sig2 += (dsa2-dsa*dsa)/(npee-1);
2179 sb_sig2 += (dsb2-dsb*dsb)/(npee-1);
2180
2181 dstddev /= npee;
2182 dstddev2 /= npee;
2183 stddev_sig2 += (dstddev2-dstddev*dstddev)/(npee-1);
2184 }
2185 br->dg_err = sqrt(dg_sig2/(npee_max - npee_min + 1));
2186 br->sa_err = sqrt(sa_sig2/(npee_max - npee_min + 1));
2187 br->sb_err = sqrt(sb_sig2/(npee_max - npee_min + 1));
2188 br->dg_stddev_err = sqrt(stddev_sig2/(npee_max - npee_min + 1));
2189 }
2190}
2191
2192
2193static double bar_err(int nbmin, int nbmax, const double *partsum)
2194{
2195 int nb, b;
2196 double svar, s, s2, dg;
2197
2198 svar = 0;
2199 for (nb = nbmin; nb <= nbmax; nb++)
2200 {
2201 s = 0;
2202 s2 = 0;
2203 for (b = 0; b < nb; b++)
2204 {
2205 dg = partsum[nb*(nbmax+1)+b];
2206 s += dg;
2207 s2 += dg*dg;
2208 }
2209 s /= nb;
2210 s2 /= nb;
2211 svar += (s2 - s*s)/(nb - 1);
2212 }
2213
2214 return sqrt(svar/(nbmax + 1 - nbmin));
2215}
2216
2217
2218/* Seek the end of an identifier (consecutive non-spaces), followed by
2219 an optional number of spaces or '='-signs. Returns a pointer to the
2220 first non-space value found after that. Returns NULL if the string
2221 ends before that.
2222 */
2223static const char *find_value(const char *str)
2224{
2225 gmx_bool name_end_found = FALSE0;
2226
2227 /* if the string is a NULL pointer, return a NULL pointer. */
2228 if (str == NULL((void*)0))
2229 {
2230 return NULL((void*)0);
2231 }
2232 while (*str != '\0')
2233 {
2234 /* first find the end of the name */
2235 if (!name_end_found)
2236 {
2237 if (isspace(*str)((*__ctype_b_loc ())[(int) ((*str))] & (unsigned short int
) _ISspace) || (*str == '=') )
2238 {
2239 name_end_found = TRUE1;
2240 }
2241 }
2242 else
2243 {
2244 if (!( isspace(*str)((*__ctype_b_loc ())[(int) ((*str))] & (unsigned short int
) _ISspace) || (*str == '=') ))
2245 {
2246 return str;
2247 }
2248 }
2249 str++;
2250 }
2251 return NULL((void*)0);
2252}
2253
2254
2255
2256/* read a vector-notation description of a lambda vector */
2257static gmx_bool read_lambda_compvec(const char *str,
2258 lambda_vec_t *lv,
2259 const lambda_components_t *lc_in,
2260 lambda_components_t *lc_out,
2261 const char **end,
2262 const char *fn)
2263{
2264 gmx_bool initialize_lc = FALSE0; /* whether to initialize the lambda
2265 components, or to check them */
2266 gmx_bool start_reached = FALSE0; /* whether the start of component names
2267 has been reached */
2268 gmx_bool vector = FALSE0; /* whether there are multiple components */
2269 int n = 0; /* current component number */
2270 const char *val_start = NULL((void*)0); /* start of the component name, or NULL
2271 if not in a value */
2272 char *strtod_end;
2273 gmx_bool OK = TRUE1;
2274
2275 if (end)
2276 {
2277 *end = str;
2278 }
2279
2280
2281 if (lc_out && lc_out->N == 0)
2282 {
2283 initialize_lc = TRUE1;
2284 }
2285
2286 if (lc_in == NULL((void*)0))
2287 {
2288 lc_in = lc_out;
2289 }
2290
2291 while (1)
2292 {
2293 if (!start_reached)
2294 {
2295 if (isalnum(*str)((*__ctype_b_loc ())[(int) ((*str))] & (unsigned short int
) _ISalnum))
2296 {
2297 vector = FALSE0;
2298 start_reached = TRUE1;
2299 val_start = str;
2300 }
2301 else if (*str == '(')
2302 {
2303 vector = TRUE1;
2304 start_reached = TRUE1;
2305 }
2306 else if (!isspace(*str)((*__ctype_b_loc ())[(int) ((*str))] & (unsigned short int
) _ISspace))
2307 {
2308 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2308, "Error in lambda components in %s", fn);
2309 }
2310 }
2311 else
2312 {
2313 if (val_start)
2314 {
2315 if (isspace(*str)((*__ctype_b_loc ())[(int) ((*str))] & (unsigned short int
) _ISspace) || *str == ')' || *str == ',' || *str == '\0')
2316 {
2317 /* end of value */
2318 if (lv == NULL((void*)0))
2319 {
2320 if (initialize_lc)
2321 {
2322 lambda_components_add(lc_out, val_start,
2323 (str-val_start));
2324 }
2325 else
2326 {
2327 if (!lambda_components_check(lc_out, n, val_start,
2328 (str-val_start)))
2329 {
2330 return FALSE0;
2331 }
2332 }
2333 }
2334 else
2335 {
2336 /* add a vector component to lv */
2337 lv->val[n] = strtod(val_start, &strtod_end);
2338 if (val_start == strtod_end)
2339 {
2340 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2340,
2341 "Error reading lambda vector in %s", fn);
2342 }
2343 }
2344 /* reset for the next identifier */
2345 val_start = NULL((void*)0);
2346 n++;
2347 if (!vector)
2348 {
2349 return OK;
2350 }
2351 }
2352 }
2353 else if (isalnum(*str)((*__ctype_b_loc ())[(int) ((*str))] & (unsigned short int
) _ISalnum))
2354 {
2355 val_start = str;
2356 }
2357 if (*str == ')')
2358 {
2359 str++;
2360 if (end)
2361 {
2362 *end = str;
2363 }
2364 if (!vector)
2365 {
2366 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2366, "Error in lambda components in %s", fn);
2367 }
2368 else
2369 {
2370 if (n == lc_in->N)
2371 {
2372 return OK;
2373 }
2374 else if (lv == NULL((void*)0))
2375 {
2376 return FALSE0;
2377 }
2378 else
2379 {
2380 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2380, "Incomplete lambda vector data in %s",
2381 fn);
2382 return FALSE0;
2383 }
2384
2385 }
2386 }
2387 }
2388 if (*str == '\0')
2389 {
2390 break;
2391 }
2392 str++;
2393 if (end)
2394 {
2395 *end = str;
2396 }
2397 }
2398 if (vector)
2399 {
2400 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2400, "Incomplete lambda components data in %s", fn);
2401 return FALSE0;
2402 }
2403 return OK;
2404}
2405
2406/* read and check the component names from a string */
2407static gmx_bool read_lambda_components(const char *str,
2408 lambda_components_t *lc,
2409 const char **end,
2410 const char *fn)
2411{
2412 return read_lambda_compvec(str, NULL((void*)0), NULL((void*)0), lc, end, fn);
2413}
2414
2415/* read an initialized lambda vector from a string */
2416static gmx_bool read_lambda_vector(const char *str,
2417 lambda_vec_t *lv,
2418 const char **end,
2419 const char *fn)
2420{
2421 return read_lambda_compvec(str, lv, lv->lc, NULL((void*)0), end, fn);
2422}
2423
2424
2425
2426/* deduce lambda value from legend.
2427 fn = the file name
2428 legend = the legend string
2429 ba = the xvg data
2430 lam = the initialized lambda vector
2431 returns whether to use the data in this set.
2432 */
2433static gmx_bool legend2lambda(const char *fn,
2434 const char *legend,
2435 lambda_vec_t *lam)
2436{
2437 double lambda = 0;
2438 const char *ptr = NULL((void*)0), *ptr2 = NULL((void*)0);
2439 gmx_bool ok = FALSE0;
2440 gmx_bool bdhdl = FALSE0;
2441 const char *tostr = " to ";
2442
2443 if (legend == NULL((void*)0))
2444 {
2445 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2445, "There is no legend in file '%s', can not deduce lambda", fn);
2446 }
2447
2448 /* look for the last 'to': */
2449 ptr2 = legend;
2450 do
2451 {
2452 ptr2 = strstr(ptr2, tostr);
2453 if (ptr2 != NULL((void*)0))
2454 {
2455 ptr = ptr2;
2456 ptr2++;
2457 }
2458 }
2459 while (ptr2 != NULL((void*)0) && *ptr2 != '\0');
2460
2461 if (ptr)
2462 {
2463 ptr += strlen(tostr)-1; /* and advance past that 'to' */
2464 }
2465 else
2466 {
2467 /* look for the = sign */
2468 ptr = strrchr(legend, '=');
2469 if (!ptr)
2470 {
2471 /* otherwise look for the last space */
2472 ptr = strrchr(legend, ' ');
2473 }
2474 }
2475
2476 if (strstr(legend, "dH"))
2477 {
2478 ok = TRUE1;
2479 bdhdl = TRUE1;
2480 }
2481 else if (strchr(legend, 'D')(__extension__ (__builtin_constant_p ('D') && !__builtin_constant_p
(legend) && ('D') == '\0' ? (char *) __rawmemchr (legend
, 'D') : __builtin_strchr (legend, 'D'))) != NULL((void*)0) && strchr(legend, 'H')(__extension__ (__builtin_constant_p ('H') && !__builtin_constant_p
(legend) && ('H') == '\0' ? (char *) __rawmemchr (legend
, 'H') : __builtin_strchr (legend, 'H'))) != NULL((void*)0))
2482 {
2483 ok = TRUE1;
2484 bdhdl = FALSE0;
2485 }
2486 else /*if (strstr(legend, "pV"))*/
2487 {
2488 return FALSE0;
2489 }
2490 if (!ptr)
2491 {
2492 ok = FALSE0;
2493 }
2494
2495 if (!ok)
2496 {
2497 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2497, "There is no proper lambda legend in file '%s', can not deduce lambda", fn);
2498 }
2499 if (!bdhdl)
2500 {
2501 ptr = find_value(ptr);
2502 if (!ptr || !read_lambda_vector(ptr, lam, NULL((void*)0), fn))
2503 {
2504 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2504, "lambda vector '%s' %s faulty", legend, fn);
2505 }
2506 }
2507 else
2508 {
2509 int dhdl_index;
2510 const char *end;
2511 char buf[STRLEN4096];
2512
2513 ptr = strrchr(legend, '=');
2514 end = ptr;
2515 if (ptr)
2516 {
2517 /* there must be a component name */
2518 ptr--;
2519 if (ptr < legend)
2520 {
2521 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2521, "dhdl legend '%s' %s faulty", legend, fn);
2522 }
2523 /* now backtrack to the start of the identifier */
2524 while (isspace(*ptr)((*__ctype_b_loc ())[(int) ((*ptr))] & (unsigned short int
) _ISspace))
2525 {
2526 end = ptr;
2527 ptr--;
2528 if (ptr < legend)
2529 {
2530 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2530, "dhdl legend '%s' %s faulty", legend, fn);
2531 }
2532 }
2533 while (!isspace(*ptr)((*__ctype_b_loc ())[(int) ((*ptr))] & (unsigned short int
) _ISspace))
2534 {
2535 ptr--;
2536 if (ptr < legend)
2537 {
2538 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2538, "dhdl legend '%s' %s faulty", legend, fn);
2539 }
2540 }
2541 ptr++;
2542 strncpy(buf, ptr, (end-ptr))__builtin_strncpy (buf, ptr, (end-ptr));
2543 buf[(end-ptr)] = '\0';
2544 dhdl_index = lambda_components_find(lam->lc, ptr, (end-ptr));
2545 if (dhdl_index < 0)
2546 {
2547 char buf[STRLEN4096];
2548 strncpy(buf, ptr, (end-ptr))__builtin_strncpy (buf, ptr, (end-ptr));
2549 buf[(end-ptr)] = '\0';
2550 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2550,
2551 "Did not find lambda component for '%s' in %s",
2552 buf, fn);
2553 }
2554 }
2555 else
2556 {
2557 if (lam->lc->N > 1)
2558 {
2559 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2559,
2560 "dhdl without component name with >1 lambda component in %s",
2561 fn);
2562 }
2563 dhdl_index = 0;
2564 }
2565 lam->dhdl = dhdl_index;
2566 }
2567 return TRUE1;
2568}
2569
2570static gmx_bool subtitle2lambda(const char *subtitle, xvg_t *ba, const char *fn,
2571 lambda_components_t *lc)
2572{
2573 gmx_bool bFound;
2574 const char *ptr;
2575 char *end;
2576 double native_lambda;
2577
2578 bFound = FALSE0;
2579
2580 /* first check for a state string */
2581 ptr = strstr(subtitle, "state");
2582 if (ptr)
2583 {
2584 int index = -1;
2585 const char *val_end;
2586
2587 /* the new 4.6 style lambda vectors */
2588 ptr = find_value(ptr);
2589 if (ptr)
2590 {
2591 index = strtol(ptr, &end, 10);
2592 if (ptr == end)
2593 {
2594 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2594, "Incomplete state data in %s", fn);
2595 return FALSE0;
2596 }
2597 ptr = end;
2598 }
2599 else
2600 {
2601 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2601, "Incomplete state data in %s", fn);
2602 return FALSE0;
2603 }
2604 /* now find the lambda vector component names */
2605 while (*ptr != '(' && !isalnum(*ptr)((*__ctype_b_loc ())[(int) ((*ptr))] & (unsigned short int
) _ISalnum))
2606 {
2607 ptr++;
2608 if (*ptr == '\0')
2609 {
2610 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2610,
2611 "Incomplete lambda vector component data in %s", fn);
2612 return FALSE0;
2613 }
2614 }
2615 val_end = ptr;
2616 if (!read_lambda_components(ptr, lc, &val_end, fn))
2617 {
2618 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2618,
2619 "lambda vector components in %s don't match those previously read",
2620 fn);
2621 }
2622 ptr = find_value(val_end);
2623 if (!ptr)
2624 {
2625 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2625, "Incomplete state data in %s", fn);
2626 return FALSE0;
2627 }
2628 lambda_vec_init(&(ba->native_lambda), lc);
2629 if (!read_lambda_vector(ptr, &(ba->native_lambda), NULL((void*)0), fn))
2630 {
2631 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2631, "lambda vector in %s faulty", fn);
2632 }
2633 ba->native_lambda.index = index;
2634 bFound = TRUE1;
2635 }
2636 else
2637 {
2638 /* compatibility mode: check for lambda in other ways. */
2639 /* plain text lambda string */
2640 ptr = strstr(subtitle, "lambda");
2641 if (ptr == NULL((void*)0))
2642 {
2643 /* xmgrace formatted lambda string */
2644 ptr = strstr(subtitle, "\\xl\\f{}");
2645 }
2646 if (ptr == NULL((void*)0))
2647 {
2648 /* xmgr formatted lambda string */
2649 ptr = strstr(subtitle, "\\8l\\4");
2650 }
2651 if (ptr != NULL((void*)0))
2652 {
2653 ptr = strstr(ptr, "=");
2654 }
2655 if (ptr != NULL((void*)0))
2656 {
2657 bFound = (sscanf(ptr+1, "%lf", &(native_lambda)) == 1);
2658 /* add the lambda component name as an empty string */
2659 if (lc->N > 0)
2660 {
2661 if (!lambda_components_check(lc, 0, "", 0))
2662 {
2663 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2663,
2664 "lambda vector components in %s don't match those previously read",
2665 fn);
2666 }
2667 }
2668 else
2669 {
2670 lambda_components_add(lc, "", 0);
2671 }
2672 lambda_vec_init(&(ba->native_lambda), lc);
2673 ba->native_lambda.val[0] = native_lambda;
2674 }
2675 }
2676
2677 return bFound;
2678}
2679
2680static void filename2lambda(const char *fn)
2681{
2682 double lambda;
2683 const char *ptr, *digitptr;
2684 char *endptr;
2685 int dirsep;
2686 ptr = fn;
2687 /* go to the end of the path string and search backward to find the last
2688 directory in the path which has to contain the value of lambda
2689 */
2690 while (ptr[1] != '\0')
2691 {
2692 ptr++;
2693 }
2694 /* searching backward to find the second directory separator */
2695 dirsep = 0;
2696 digitptr = NULL((void*)0);
2697 while (ptr >= fn)
2698 {
2699 if (ptr[0] != DIR_SEPARATOR'/' && ptr[1] == DIR_SEPARATOR'/')
2700 {
2701 if (dirsep == 1)
2702 {
2703 break;
2704 }
2705 dirsep++;
2706 }
2707 /* save the last position of a digit between the last two
2708 separators = in the last dirname */
2709 if (dirsep > 0 && isdigit(*ptr)((*__ctype_b_loc ())[(int) ((*ptr))] & (unsigned short int
) _ISdigit))
2710 {
2711 digitptr = ptr;
2712 }
2713 ptr--;
2714 }
2715 if (!digitptr)
2716 {
2717 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2717, "While trying to read the lambda value from the file path:"
2718 " last directory in the path '%s' does not contain a number", fn);
2719 }
2720 if (digitptr[-1] == '-')
2721 {
2722 digitptr--;
2723 }
2724 lambda = strtod(digitptr, &endptr);
2725 if (endptr == digitptr)
2726 {
2727 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2727, "Malformed number in file path '%s'", fn);
2728 }
2729}
2730
2731static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
2732 lambda_components_t *lc)
2733{
2734 int i;
2735 char *subtitle, **legend, *ptr;
2736 int np;
2737 gmx_bool native_lambda_read = FALSE0;
2738 char buf[STRLEN4096];
2739 lambda_vec_t lv;
2740
2741 xvg_init(ba);
2742
2743 ba->filename = fn;
2744
2745 np = read_xvg_legend(fn, &ba->y, &ba->nset, &subtitle, &legend);
2746 if (!ba->y)
2747 {
2748 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2748, "File %s contains no usable data.", fn);
2749 }
2750 /* Reorder the data */
2751 ba->t = ba->y[0];
2752 for (i = 1; i < ba->nset; i++)
2753 {
2754 ba->y[i-1] = ba->y[i];
2755 }
2756 ba->nset--;
2757
2758 snew(ba->np, ba->nset)(ba->np) = save_calloc("ba->np", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2758, (ba->nset), sizeof(*(ba->np)));
2759 for (i = 0; i < ba->nset; i++)
2760 {
2761 ba->np[i] = np;
2762 }
2763
2764 ba->temp = -1;
2765 if (subtitle != NULL((void*)0))
2766 {
2767 /* try to extract temperature */
2768 ptr = strstr(subtitle, "T =");
2769 if (ptr != NULL((void*)0))
2770 {
2771 ptr += 3;
2772 if (sscanf(ptr, "%lf", &ba->temp) == 1)
2773 {
2774 if (ba->temp <= 0)
2775 {
2776 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2776, "Found temperature of %f in file '%s'",
2777 ba->temp, fn);
2778 }
2779 }
2780 }
2781 }
2782 if (ba->temp < 0)
2783 {
2784 if (*temp <= 0)
2785 {
2786 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2786, "Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]gmx bar[tt]", fn);
2787 }
2788 ba->temp = *temp;
2789 }
2790
2791 /* Try to deduce lambda from the subtitle */
2792 if (subtitle)
2793 {
2794 if (subtitle2lambda(subtitle, ba, fn, lc))
2795 {
2796 native_lambda_read = TRUE1;
2797 }
2798 }
2799 snew(ba->lambda, ba->nset)(ba->lambda) = save_calloc("ba->lambda", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2799, (ba->nset), sizeof(*(ba->lambda)));
2800 if (legend == NULL((void*)0))
2801 {
2802 /* Check if we have a single set, no legend, nset=1 means t and dH/dl */
2803 if (ba->nset == 1)
2804 {
2805 if (!native_lambda_read)
2806 {
2807 /* Deduce lambda from the file name */
2808 filename2lambda(fn);
2809 native_lambda_read = TRUE1;
2810 }
2811 ba->lambda[0] = ba->native_lambda;
2812 }
2813 else
2814 {
2815 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2815, "File %s contains multiple sets but no legends, can not determine the lambda values", fn);
2816 }
2817 }
2818 else
2819 {
2820 for (i = 0; i < ba->nset; )
2821 {
2822 gmx_bool use = FALSE0;
2823 /* Read lambda from the legend */
2824 lambda_vec_init( &(ba->lambda[i]), lc );
2825 lambda_vec_copy( &(ba->lambda[i]), &(ba->native_lambda));
2826 use = legend2lambda(fn, legend[i], &(ba->lambda[i]));
2827 if (use)
2828 {
2829 lambda_vec_print(&(ba->lambda[i]), buf, FALSE0);
2830 i++;
2831 }
2832 else
2833 {
2834 int j;
2835 printf("%s: Ignoring set '%s'.\n", fn, legend[i]);
2836 for (j = i+1; j < ba->nset; j++)
2837 {
2838 ba->y[j-1] = ba->y[j];
2839 legend[j-1] = legend[j];
2840 }
2841 ba->nset--;
2842 }
2843 }
2844 }
2845
2846 if (!native_lambda_read)
2847 {
2848 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2848, "File %s contains multiple sets but no indication of the native lambda", fn);
2849 }
2850
2851 if (legend != NULL((void*)0))
2852 {
2853 for (i = 0; i < ba->nset-1; i++)
2854 {
2855 sfree(legend[i])save_free("legend[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2855, (legend[i]));
2856 }
2857 sfree(legend)save_free("legend", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2857, (legend));
2858 }
2859}
2860
2861static void read_bar_xvg(char *fn, real *temp, sim_data_t *sd)
2862{
2863 xvg_t *barsim;
2864 samples_t *s;
2865 int i;
2866 double *lambda;
2867
2868 snew(barsim, 1)(barsim) = save_calloc("barsim", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2868, (1), sizeof(*(barsim)));
2869
2870 read_bar_xvg_lowlevel(fn, temp, barsim, &(sd->lc));
2871
2872 if (barsim->nset < 1)
2873 {
2874 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2874, "File '%s' contains fewer than two columns", fn);
2875 }
2876
2877 if (!gmx_within_tol(*temp, barsim->temp, GMX_FLOAT_EPS5.96046448E-08) && (*temp > 0) )
2878 {
2879 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2879, "Temperature in file %s different from earlier files or setting\n", fn);
2880 }
2881 *temp = barsim->temp;
2882
2883 /* now create a series of samples_t */
2884 snew(s, barsim->nset)(s) = save_calloc("s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2884, (barsim->nset), sizeof(*(s)));
2885 for (i = 0; i < barsim->nset; i++)
2886 {
2887 samples_init(s+i, &(barsim->native_lambda), &(barsim->lambda[i]),
2888 barsim->temp, lambda_vec_same(&(barsim->native_lambda),
2889 &(barsim->lambda[i])),
2890 fn);
2891 s[i].du = barsim->y[i];
2892 s[i].ndu = barsim->np[i];
2893 s[i].t = barsim->t;
2894
2895 lambda_data_list_insert_sample(sd->lb, s+i);
2896 }
2897 {
2898 char buf[STRLEN4096];
2899
2900 lambda_vec_print(s[0].native_lambda, buf, FALSE0);
2901 printf("%s: %.1f - %.1f; lambda = %s\n dH/dl & foreign lambdas:\n",
2902 fn, s[0].t[0], s[0].t[s[0].ndu-1], buf);
2903 for (i = 0; i < barsim->nset; i++)
2904 {
2905 lambda_vec_print(s[i].foreign_lambda, buf, TRUE1);
2906 printf(" %s (%d pts)\n", buf, s[i].ndu);
2907 }
2908 }
2909 printf("\n\n");
2910}
2911
2912static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
2913 double start_time, double delta_time,
2914 lambda_vec_t *native_lambda, double temp,
2915 double *last_t, const char *filename)
2916{
2917 int i, j;
2918 gmx_bool allocated;
2919 double old_foreign_lambda;
2920 lambda_vec_t *foreign_lambda;
2921 int type;
2922 samples_t *s; /* convenience pointer */
2923 int startj;
2924
2925 /* check the block types etc. */
2926 if ( (blk->nsub < 3) ||
2927 (blk->sub[0].type != xdr_datatype_int) ||
2928 (blk->sub[1].type != xdr_datatype_double) ||
2929 (
2930 (blk->sub[2].type != xdr_datatype_float) &&
2931 (blk->sub[2].type != xdr_datatype_double)
2932 ) ||
2933 (blk->sub[0].nr < 1) ||
2934 (blk->sub[1].nr < 1) )
2935 {
2936 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2936,
2937 "Unexpected/corrupted block data in file %s around time %f.",
2938 filename, start_time);
2939 }
2940
2941 snew(foreign_lambda, 1)(foreign_lambda) = save_calloc("foreign_lambda", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2941, (1), sizeof(*(foreign_lambda)));
2942 lambda_vec_init(foreign_lambda, native_lambda->lc);
2943 lambda_vec_copy(foreign_lambda, native_lambda);
2944 type = blk->sub[0].ival[0];
2945 if (type == dhbtDH)
2946 {
2947 for (i = 0; i < native_lambda->lc->N; i++)
2948 {
2949 foreign_lambda->val[i] = blk->sub[1].dval[i];
2950 }
2951 }
2952 else
2953 {
2954 if (blk->sub[0].nr > 1)
2955 {
2956 foreign_lambda->dhdl = blk->sub[0].ival[1];
2957 }
2958 else
2959 {
2960 foreign_lambda->dhdl = 0;
2961 }
2962 }
2963
2964 if (!*smp)
2965 {
2966 /* initialize the samples structure if it's empty. */
2967 snew(*smp, 1)(*smp) = save_calloc("*smp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2967, (1), sizeof(*(*smp)));
2968 samples_init(*smp, native_lambda, foreign_lambda, temp,
2969 type == dhbtDHDL, filename);
2970 (*smp)->start_time = start_time;
2971 (*smp)->delta_time = delta_time;
2972 }
2973
2974 /* set convenience pointer */
2975 s = *smp;
2976
2977 /* now double check */
2978 if (!lambda_vec_same(s->foreign_lambda, foreign_lambda) )
2979 {
2980 char buf[STRLEN4096], buf2[STRLEN4096];
2981 lambda_vec_print(foreign_lambda, buf, FALSE0);
2982 lambda_vec_print(s->foreign_lambda, buf2, FALSE0);
2983 fprintf(stderrstderr, "Got foreign lambda=%s, expected: %s\n", buf, buf2);
2984 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2984, "Corrupted data in file %s around t=%f.",
2985 filename, start_time);
2986 }
2987
2988 /* make room for the data */
2989 if (s->ndu_alloc < (size_t)(s->ndu + blk->sub[2].nr) )
2990 {
2991 s->ndu_alloc += (s->ndu_alloc < (size_t)blk->sub[2].nr) ?
2992 blk->sub[2].nr*2 : s->ndu_alloc;
2993 srenew(s->du_alloc, s->ndu_alloc)(s->du_alloc) = save_realloc("s->du_alloc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 2993, (s->du_alloc), (s->ndu_alloc), sizeof(*(s->du_alloc
)));
2994 s->du = s->du_alloc;
2995 }
2996 startj = s->ndu;
2997 s->ndu += blk->sub[2].nr;
2998 s->ntot += blk->sub[2].nr;
2999 *ndu = blk->sub[2].nr;
3000
3001 /* and copy the data*/
3002 for (j = 0; j < blk->sub[2].nr; j++)
3003 {
3004 if (blk->sub[2].type == xdr_datatype_float)
3005 {
3006 s->du[startj+j] = blk->sub[2].fval[j];
3007 }
3008 else
3009 {
3010 s->du[startj+j] = blk->sub[2].dval[j];
3011 }
3012 }
3013 if (start_time + blk->sub[2].nr*delta_time > *last_t)
3014 {
3015 *last_t = start_time + blk->sub[2].nr*delta_time;
3016 }
3017}
3018
3019static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
3020 double start_time, double delta_time,
3021 lambda_vec_t *native_lambda, double temp,
3022 double *last_t, const char *filename)
3023{
3024 int i, j;
3025 samples_t *s;
3026 int nhist;
3027 double old_foreign_lambda;
3028 lambda_vec_t *foreign_lambda;
3029 int type;
3030 int nbins[2];
3031
3032 /* check the block types etc. */
3033 if ( (blk->nsub < 2) ||
3034 (blk->sub[0].type != xdr_datatype_double) ||
3035 (blk->sub[1].type != xdr_datatype_int64) ||
3036 (blk->sub[0].nr < 2) ||
3037 (blk->sub[1].nr < 2) )
3038 {
3039 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3039,
3040 "Unexpected/corrupted block data in file %s around time %f",
3041 filename, start_time);
3042 }
3043
3044 nhist = blk->nsub-2;
3045 if (nhist == 0)
3046 {
3047 return NULL((void*)0);
3048 }
3049 if (nhist > 2)
3050 {
3051 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3051,
3052 "Unexpected/corrupted block data in file %s around time %f",
3053 filename, start_time);
3054 }
3055
3056 snew(s, 1)(s) = save_calloc("s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3056, (1), sizeof(*(s)));
3057 *nsamples = 1;
3058
3059 snew(foreign_lambda, 1)(foreign_lambda) = save_calloc("foreign_lambda", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3059, (1), sizeof(*(foreign_lambda)));
3060 lambda_vec_init(foreign_lambda, native_lambda->lc);
3061 lambda_vec_copy(foreign_lambda, native_lambda);
3062 type = (int)(blk->sub[1].lval[1]);
3063 if (type == dhbtDH)
3064 {
3065 double old_foreign_lambda;
3066
3067 old_foreign_lambda = blk->sub[0].dval[0];
3068 if (old_foreign_lambda >= 0)
3069 {
3070 foreign_lambda->val[0] = old_foreign_lambda;
3071 if (foreign_lambda->lc->N > 1)
3072 {
3073 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3073,
3074 "Single-component lambda in multi-component file %s",
3075 filename);
3076 }
3077 }
3078 else
3079 {
3080 for (i = 0; i < native_lambda->lc->N; i++)
3081 {
3082 foreign_lambda->val[i] = blk->sub[0].dval[i+2];
3083 }
3084 }
3085 }
3086 else
3087 {
3088 if (foreign_lambda->lc->N > 1)
3089 {
3090 if (blk->sub[1].nr < 3 + nhist)
3091 {
3092 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3092,
3093 "Missing derivative coord in multi-component file %s",
3094 filename);
3095 }
3096 foreign_lambda->dhdl = blk->sub[1].lval[2 + nhist];
3097 }
3098 else
3099 {
3100 foreign_lambda->dhdl = 0;
3101 }
3102 }
3103
3104 samples_init(s, native_lambda, foreign_lambda, temp, type == dhbtDHDL,
3105 filename);
3106 snew(s->hist, 1)(s->hist) = save_calloc("s->hist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3106, (1), sizeof(*(s->hist)));
3107
3108 for (i = 0; i < nhist; i++)
3109 {
3110 nbins[i] = blk->sub[i+2].nr;
3111 }
3112
3113 hist_init(s->hist, nhist, nbins);
3114
3115 for (i = 0; i < nhist; i++)
3116 {
3117 s->hist->x0[i] = blk->sub[1].lval[2+i];
3118 s->hist->dx[i] = blk->sub[0].dval[1];
3119 if (i == 1)
3120 {
3121 s->hist->dx[i] = -s->hist->dx[i];
3122 }
3123 }
3124
3125 s->hist->start_time = start_time;
3126 s->hist->delta_time = delta_time;
3127 s->start_time = start_time;
3128 s->delta_time = delta_time;
3129
3130 for (i = 0; i < nhist; i++)
3131 {
3132 int nbin;
3133 gmx_int64_t sum = 0;
3134
3135 for (j = 0; j < s->hist->nbin[i]; j++)
3136 {
3137 int binv = (int)(blk->sub[i+2].ival[j]);
3138
3139 s->hist->bin[i][j] = binv;
3140 sum += binv;
3141
3142 }
3143 if (i == 0)
3144 {
3145 s->ntot = sum;
3146 s->hist->sum = sum;
3147 }
3148 else
3149 {
3150 if (s->ntot != sum)
3151 {
3152 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3152, "Histogram counts don't match in %s",
3153 filename);
3154 }
3155 }
3156 }
3157
3158 if (start_time + s->hist->sum*delta_time > *last_t)
3159 {
3160 *last_t = start_time + s->hist->sum*delta_time;
3161 }
3162 return s;
3163}
3164
3165
3166static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
3167{
3168 int i, j;
3169 ener_file_t fp;
3170 t_enxframe *fr;
3171 int nre;
3172 gmx_enxnm_t *enm = NULL((void*)0);
3173 double first_t = -1;
3174 double last_t = -1;
3175 samples_t **samples_rawdh = NULL((void*)0); /* contains samples for raw delta_h */
3176 int *nhists = NULL((void*)0); /* array to keep count & print at end */
3177 int *npts = NULL((void*)0); /* array to keep count & print at end */
3178 lambda_vec_t **lambdas = NULL((void*)0); /* array to keep count & print at end */
3179 lambda_vec_t *native_lambda;
3180 double end_time; /* the end time of the last batch of samples */
3181 int nsamples = 0;
3182 lambda_vec_t start_lambda;
3183
3184 fp = open_enx(fn, "r");
3185 do_enxnms(fp, &nre, &enm);
3186 snew(fr, 1)(fr) = save_calloc("fr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3186, (1), sizeof(*(fr)));
3187
3188 snew(native_lambda, 1)(native_lambda) = save_calloc("native_lambda", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3188, (1), sizeof(*(native_lambda)));
3189 start_lambda.lc = NULL((void*)0);
3190
3191 while (do_enx(fp, fr))
3192 {
3193 /* count the data blocks */
3194 int nblocks_raw = 0;
3195 int nblocks_hist = 0;
3196 int nlam = 0;
3197 int k;
3198 /* DHCOLL block information: */
3199 double start_time = 0, delta_time = 0, old_start_lambda = 0, delta_lambda = 0;
3200 double rtemp = 0;
3201
3202 /* count the blocks and handle collection information: */
3203 for (i = 0; i < fr->nblock; i++)
3204 {
3205 if (fr->block[i].id == enxDHHIST)
3206 {
3207 nblocks_hist++;
3208 }
3209 if (fr->block[i].id == enxDH)
3210 {
3211 nblocks_raw++;
3212 }
3213 if (fr->block[i].id == enxDHCOLL)
3214 {
3215 nlam++;
3216 if ( (fr->block[i].nsub < 1) ||
3217 (fr->block[i].sub[0].type != xdr_datatype_double) ||
3218 (fr->block[i].sub[0].nr < 5))
3219 {
3220 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3220, "Unexpected block data in file %s", fn);
3221 }
3222
3223 /* read the data from the DHCOLL block */
3224 rtemp = fr->block[i].sub[0].dval[0];
3225 start_time = fr->block[i].sub[0].dval[1];
3226 delta_time = fr->block[i].sub[0].dval[2];
3227 old_start_lambda = fr->block[i].sub[0].dval[3];
3228 delta_lambda = fr->block[i].sub[0].dval[4];
3229
3230 if (delta_lambda > 0)
3231 {
3232 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3232, "Lambda values not constant in %s: can't apply BAR method", fn);
3233 }
3234 if ( ( *temp != rtemp) && (*temp > 0) )
3235 {
3236 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3236, "Temperature in file %s different from earlier files or setting\n", fn);
3237 }
3238 *temp = rtemp;
3239
3240 if (old_start_lambda >= 0)
3241 {
3242 if (sd->lc.N > 0)
3243 {
3244 if (!lambda_components_check(&(sd->lc), 0, "", 0))
3245 {
3246 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3246,
3247 "lambda vector components in %s don't match those previously read",
3248 fn);
3249 }
3250 }
3251 else
3252 {
3253 lambda_components_add(&(sd->lc), "", 0);
3254 }
3255 if (!start_lambda.lc)
3256 {
3257 lambda_vec_init(&start_lambda, &(sd->lc));
3258 }
3259 start_lambda.val[0] = old_start_lambda;
3260 }
3261 else
3262 {
3263 /* read lambda vector */
3264 int n_lambda_vec;
3265 gmx_bool check = (sd->lc.N > 0);
3266 if (fr->block[i].nsub < 2)
3267 {
3268 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3268,
3269 "No lambda vector, but start_lambda=%f\n",
3270 old_start_lambda);
3271 }
3272 n_lambda_vec = fr->block[i].sub[1].ival[1];
3273 for (j = 0; j < n_lambda_vec; j++)
3274 {
3275 const char *name =
3276 efpt_singular_names[fr->block[i].sub[1].ival[1+j]];
3277 if (check)
3278 {
3279 /* check the components */
3280 lambda_components_check(&(sd->lc), j, name,
3281 strlen(name));
3282 }
3283 else
3284 {
3285 lambda_components_add(&(sd->lc), name,
3286 strlen(name));
3287 }
3288 }
3289 lambda_vec_init(&start_lambda, &(sd->lc));
3290 start_lambda.index = fr->block[i].sub[1].ival[0];
3291 for (j = 0; j < n_lambda_vec; j++)
3292 {
3293 start_lambda.val[j] = fr->block[i].sub[0].dval[5+j];
3294 }
3295 }
3296 if (first_t < 0)
3297 {
3298 first_t = start_time;
3299 }
3300 }
3301 }
3302
3303 if (nlam != 1)
3304 {
3305 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3305, "Did not find delta H information in file %s", fn);
3306 }
3307 if (nblocks_raw > 0 && nblocks_hist > 0)
3308 {
3309 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3309, "Can't handle both raw delta U data and histograms in the same file %s", fn);
3310 }
3311
3312 if (nsamples > 0)
3313 {
3314 /* check the native lambda */
3315 if (!lambda_vec_same(&start_lambda, native_lambda) )
3316 {
3317 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3317, "Native lambda not constant in file %s: started at %f, and becomes %f at time %f",
3318 fn, native_lambda, start_lambda, start_time);
3319 }
3320 /* check the number of samples against the previous number */
3321 if ( ((nblocks_raw+nblocks_hist) != nsamples) || (nlam != 1 ) )
3322 {
3323 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3323, "Unexpected block count in %s: was %d, now %d\n",
3324 fn, nsamples+1, nblocks_raw+nblocks_hist+nlam);
3325 }
3326 /* check whether last iterations's end time matches with
3327 the currrent start time */
3328 if ( (fabs(last_t - start_time) > 2*delta_time) && last_t >= 0)
3329 {
3330 /* it didn't. We need to store our samples and reallocate */
3331 for (i = 0; i < nsamples; i++)
3332 {
3333 if (samples_rawdh[i])
3334 {
3335 /* insert it into the existing list */
3336 lambda_data_list_insert_sample(sd->lb,
3337 samples_rawdh[i]);
3338 /* and make sure we'll allocate a new one this time
3339 around */
3340 samples_rawdh[i] = NULL((void*)0);
3341 }
3342 }
3343 }
3344 }
3345 else
3346 {
3347 /* this is the first round; allocate the associated data
3348 structures */
3349 /*native_lambda=start_lambda;*/
3350 lambda_vec_init(native_lambda, &(sd->lc));
3351 lambda_vec_copy(native_lambda, &start_lambda);
3352 nsamples = nblocks_raw+nblocks_hist;
3353 snew(nhists, nsamples)(nhists) = save_calloc("nhists", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3353, (nsamples), sizeof(*(nhists)));
3354 snew(npts, nsamples)(npts) = save_calloc("npts", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3354, (nsamples), sizeof(*(npts)));
3355 snew(lambdas, nsamples)(lambdas) = save_calloc("lambdas", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3355, (nsamples), sizeof(*(lambdas)));
3356 snew(samples_rawdh, nsamples)(samples_rawdh) = save_calloc("samples_rawdh", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3356, (nsamples), sizeof(*(samples_rawdh)));
3357 for (i = 0; i < nsamples; i++)
3358 {
3359 nhists[i] = 0;
3360 npts[i] = 0;
3361 lambdas[i] = NULL((void*)0);
3362 samples_rawdh[i] = NULL((void*)0); /* init to NULL so we know which
3363 ones contain values */
3364 }
3365 }
3366
3367 /* and read them */
3368 k = 0; /* counter for the lambdas, etc. arrays */
3369 for (i = 0; i < fr->nblock; i++)
3370 {
3371 if (fr->block[i].id == enxDH)
3372 {
3373 int type = (fr->block[i].sub[0].ival[0]);
3374 if (type == dhbtDH || type == dhbtDHDL)
3375 {
3376 int ndu;
3377 read_edr_rawdh_block(&(samples_rawdh[k]),
3378 &ndu,
3379 &(fr->block[i]),
3380 start_time, delta_time,
3381 native_lambda, rtemp,
3382 &last_t, fn);
3383 npts[k] += ndu;
3384 if (samples_rawdh[k])
3385 {
3386 lambdas[k] = samples_rawdh[k]->foreign_lambda;
3387 }
3388 k++;
3389 }
3390 }
3391 else if (fr->block[i].id == enxDHHIST)
3392 {
3393 int type = (int)(fr->block[i].sub[1].lval[1]);
3394 if (type == dhbtDH || type == dhbtDHDL)
3395 {
3396 int j;
3397 int nb = 0;
3398 samples_t *s; /* this is where the data will go */
3399 s = read_edr_hist_block(&nb, &(fr->block[i]),
3400 start_time, delta_time,
3401 native_lambda, rtemp,
3402 &last_t, fn);
3403 nhists[k] += nb;
3404 if (nb > 0)
3405 {
3406 lambdas[k] = s->foreign_lambda;
3407 }
3408 k++;
3409 /* and insert the new sample immediately */
3410 for (j = 0; j < nb; j++)
3411 {
3412 lambda_data_list_insert_sample(sd->lb, s+j);
3413 }
3414 }
3415 }
3416 }
3417 }
3418 /* Now store all our extant sample collections */
3419 for (i = 0; i < nsamples; i++)
3420 {
3421 if (samples_rawdh[i])
3422 {
3423 /* insert it into the existing list */
3424 lambda_data_list_insert_sample(sd->lb, samples_rawdh[i]);
3425 }
3426 }
3427
3428
3429 {
3430 char buf[STRLEN4096];
3431 printf("\n");
3432 lambda_vec_print(native_lambda, buf, FALSE0);
3433 printf("%s: %.1f - %.1f; lambda = %s\n foreign lambdas:\n",
3434 fn, first_t, last_t, buf);
3435 for (i = 0; i < nsamples; i++)
3436 {
3437 if (lambdas[i])
3438 {
3439 lambda_vec_print(lambdas[i], buf, TRUE1);
3440 if (nhists[i] > 0)
3441 {
3442 printf(" %s (%d hists)\n", buf, nhists[i]);
3443 }
3444 else
3445 {
3446 printf(" %s (%d pts)\n", buf, npts[i]);
3447 }
3448 }
3449 }
3450 }
3451 printf("\n\n");
3452 sfree(npts)save_free("npts", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3452, (npts));
3453 sfree(nhists)save_free("nhists", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3453, (nhists));
3454 sfree(lambdas)save_free("lambdas", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3454, (lambdas));
3455}
3456
3457
3458int gmx_bar(int argc, char *argv[])
3459{
3460 static const char *desc[] = {
3461 "[THISMODULE] calculates free energy difference estimates through ",
3462 "Bennett's acceptance ratio method (BAR). It also automatically",
3463 "adds series of individual free energies obtained with BAR into",
3464 "a combined free energy estimate.[PAR]",
3465
3466 "Every individual BAR free energy difference relies on two ",
3467 "simulations at different states: say state A and state B, as",
3468 "controlled by a parameter, [GRK]lambda[grk] (see the [TT].mdp[tt] parameter",
3469 "[TT]init_lambda[tt]). The BAR method calculates a ratio of weighted",
3470 "average of the Hamiltonian difference of state B given state A and",
3471 "vice versa.",
3472 "The energy differences to the other state must be calculated",
3473 "explicitly during the simulation. This can be done with",
3474 "the [TT].mdp[tt] option [TT]foreign_lambda[tt].[PAR]",
3475
3476 "Input option [TT]-f[tt] expects multiple [TT]dhdl.xvg[tt] files. ",
3477 "Two types of input files are supported:[BR]",
3478 "[TT]*[tt] Files with more than one [IT]y[it]-value. ",
3479 "The files should have columns ",
3480 "with dH/d[GRK]lambda[grk] and [GRK]Delta[grk][GRK]lambda[grk]. ",
3481 "The [GRK]lambda[grk] values are inferred ",
3482 "from the legends: [GRK]lambda[grk] of the simulation from the legend of ",
3483 "dH/d[GRK]lambda[grk] and the foreign [GRK]lambda[grk] values from the ",
3484 "legends of Delta H",
3485 "[BR]",
3486 "[TT]*[tt] Files with only one [IT]y[it]-value. Using the",
3487 "[TT]-extp[tt] option for these files, it is assumed",
3488 "that the [IT]y[it]-value is dH/d[GRK]lambda[grk] and that the ",
3489 "Hamiltonian depends linearly on [GRK]lambda[grk]. ",
3490 "The [GRK]lambda[grk] value of the simulation is inferred from the ",
3491 "subtitle (if present), otherwise from a number in the subdirectory ",
3492 "in the file name.[PAR]",
3493
3494 "The [GRK]lambda[grk] of the simulation is parsed from ",
3495 "[TT]dhdl.xvg[tt] file's legend containing the string 'dH', the ",
3496 "foreign [GRK]lambda[grk] values from the legend containing the ",
3497 "capitalized letters 'D' and 'H'. The temperature is parsed from ",
3498 "the legend line containing 'T ='.[PAR]",
3499
3500 "The input option [TT]-g[tt] expects multiple [TT].edr[tt] files. ",
3501 "These can contain either lists of energy differences (see the ",
3502 "[TT].mdp[tt] option [TT]separate_dhdl_file[tt]), or a series of ",
3503 "histograms (see the [TT].mdp[tt] options [TT]dh_hist_size[tt] and ",
3504 "[TT]dh_hist_spacing[tt]).", "The temperature and [GRK]lambda[grk] ",
3505 "values are automatically deduced from the [TT]ener.edr[tt] file.[PAR]",
3506
3507 "In addition to the [TT].mdp[tt] option [TT]foreign_lambda[tt], ",
3508 "the energy difference can also be extrapolated from the ",
3509 "dH/d[GRK]lambda[grk] values. This is done with the[TT]-extp[tt]",
3510 "option, which assumes that the system's Hamiltonian depends linearly",
3511 "on [GRK]lambda[grk], which is not normally the case.[PAR]",
3512
3513 "The free energy estimates are determined using BAR with bisection, ",
3514 "with the precision of the output set with [TT]-prec[tt]. ",
3515 "An error estimate taking into account time correlations ",
3516 "is made by splitting the data into blocks and determining ",
3517 "the free energy differences over those blocks and assuming ",
3518 "the blocks are independent. ",
3519 "The final error estimate is determined from the average variance ",
3520 "over 5 blocks. A range of block numbers for error estimation can ",
3521 "be provided with the options [TT]-nbmin[tt] and [TT]-nbmax[tt].[PAR]",
3522
3523 "[THISMODULE] tries to aggregate samples with the same 'native' and ",
3524 "'foreign' [GRK]lambda[grk] values, but always assumes independent ",
3525 "samples. [BB]Note[bb] that when aggregating energy ",
3526 "differences/derivatives with different sampling intervals, this is ",
3527 "almost certainly not correct. Usually subsequent energies are ",
3528 "correlated and different time intervals mean different degrees ",
3529 "of correlation between samples.[PAR]",
3530
3531 "The results are split in two parts: the last part contains the final ",
3532 "results in kJ/mol, together with the error estimate for each part ",
3533 "and the total. The first part contains detailed free energy ",
3534 "difference estimates and phase space overlap measures in units of ",
3535 "kT (together with their computed error estimate). The printed ",
3536 "values are:[BR]",
3537 "[TT]*[tt] lam_A: the [GRK]lambda[grk] values for point A.[BR]",
3538 "[TT]*[tt] lam_B: the [GRK]lambda[grk] values for point B.[BR]",
3539 "[TT]*[tt] DG: the free energy estimate.[BR]",
3540 "[TT]*[tt] s_A: an estimate of the relative entropy of B in A.[BR]",
3541 "[TT]*[tt] s_B: an estimate of the relative entropy of A in B.[BR]",
3542 "[TT]*[tt] stdev: an estimate expected per-sample standard deviation.[PAR]",
3543
3544 "The relative entropy of both states in each other's ensemble can be ",
3545 "interpreted as a measure of phase space overlap: ",
3546 "the relative entropy s_A of the work samples of lambda_B in the ",
3547 "ensemble of lambda_A (and vice versa for s_B), is a ",
3548 "measure of the 'distance' between Boltzmann distributions of ",
3549 "the two states, that goes to zero for identical distributions. See ",
3550 "Wu & Kofke, J. Chem. Phys. 123 084109 (2005) for more information.",
3551 "[PAR]",
3552 "The estimate of the expected per-sample standard deviation, as given ",
3553 "in Bennett's original BAR paper: Bennett, J. Comp. Phys. 22, p 245 (1976).",
3554 "Eq. 10 therein gives an estimate of the quality of sampling (not directly",
3555 "of the actual statistical error, because it assumes independent samples).[PAR]",
3556
3557 "To get a visual estimate of the phase space overlap, use the ",
3558 "[TT]-oh[tt] option to write series of histograms, together with the ",
3559 "[TT]-nbin[tt] option.[PAR]"
3560 };
3561 static real begin = 0, end = -1, temp = -1;
3562 int nd = 2, nbmin = 5, nbmax = 5;
3563 int nbin = 100;
3564 gmx_bool use_dhdl = FALSE0;
3565 gmx_bool calc_s, calc_v;
3566 t_pargs pa[] = {
3567 { "-b", FALSE0, etREAL, {&begin}, "Begin time for BAR" },
3568 { "-e", FALSE0, etREAL, {&end}, "End time for BAR" },
3569 { "-temp", FALSE0, etREAL, {&temp}, "Temperature (K)" },
3570 { "-prec", FALSE0, etINT, {&nd}, "The number of digits after the decimal point" },
3571 { "-nbmin", FALSE0, etINT, {&nbmin}, "Minimum number of blocks for error estimation" },
3572 { "-nbmax", FALSE0, etINT, {&nbmax}, "Maximum number of blocks for error estimation" },
3573 { "-nbin", FALSE0, etINT, {&nbin}, "Number of bins for histogram output"},
3574 { "-extp", FALSE0, etBOOL, {&use_dhdl}, "Whether to linearly extrapolate dH/dl values to use as energies"}
3575 };
3576
3577 t_filenm fnm[] = {
3578 { efXVG, "-f", "dhdl", ffOPTRDMULT((1<<1 | 1<<5) | 1<<3) },
3579 { efEDR, "-g", "ener", ffOPTRDMULT((1<<1 | 1<<5) | 1<<3) },
3580 { efXVG, "-o", "bar", ffOPTWR(1<<2| 1<<3) },
3581 { efXVG, "-oi", "barint", ffOPTWR(1<<2| 1<<3) },
3582 { efXVG, "-oh", "histogram", ffOPTWR(1<<2| 1<<3) }
3583 };
3584#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
3585
3586 int f, i, j;
3587 int nf = 0; /* file counter */
3588 int nbs;
3589 int nfile_tot; /* total number of input files */
3590 int nxvgfile = 0;
3591 int nedrfile = 0;
3592 char **fxvgnms;
3593 char **fedrnms;
3594 sim_data_t sim_data; /* the simulation data */
3595 barres_t *results; /* the results */
3596 int nresults; /* number of results in results array */
3597
3598 double *partsum;
3599 double prec, dg_tot, dg, sig, dg_tot_max, dg_tot_min;
3600 FILE *fpb, *fpi;
3601 char dgformat[20], xvg2format[STRLEN4096], xvg3format[STRLEN4096];
3602 char buf[STRLEN4096], buf2[STRLEN4096];
3603 char ktformat[STRLEN4096], sktformat[STRLEN4096];
3604 char kteformat[STRLEN4096], skteformat[STRLEN4096];
3605 output_env_t oenv;
3606 double kT, beta;
3607 gmx_bool result_OK = TRUE1, bEE = TRUE1;
3608
3609 gmx_bool disc_err = FALSE0;
3610 double sum_disc_err = 0.; /* discretization error */
3611 gmx_bool histrange_err = FALSE0;
3612 double sum_histrange_err = 0.; /* histogram range error */
3613 double stat_err = 0.; /* statistical error */
3614
3615 if (!parse_common_args(&argc, argv,
3616 PCA_CAN_VIEW(1<<5),
3617 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv))
3618 {
3619 return 0;
3620 }
3621
3622 if (opt2bSet("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
3623 {
3624 nxvgfile = opt2fns(&fxvgnms, "-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
3625 }
3626 if (opt2bSet("-g", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
3627 {
3628 nedrfile = opt2fns(&fedrnms, "-g", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
3629 }
3630
3631 sim_data_init(&sim_data);
3632#if 0
3633 /* make linked list */
3634 lb = &lambda_head;
3635 lambda_data_init(lb, 0, 0);
3636 lb->next = lb;
3637 lb->prev = lb;
3638#endif
3639
3640
3641 nfile_tot = nxvgfile + nedrfile;
3642
3643 if (nfile_tot == 0)
3644 {
3645 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3645, "No input files!");
3646 }
3647
3648 if (nd < 0)
3649 {
3650 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3650, "Can not have negative number of digits");
3651 }
3652 prec = pow(10, -nd);
3653
3654 snew(partsum, (nbmax+1)*(nbmax+1))(partsum) = save_calloc("partsum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_bar.c"
, 3654, ((nbmax+1)*(nbmax+1)), sizeof(*(partsum)));
3655 nf = 0;
3656
3657 /* read in all files. First xvg files */
3658 for (f = 0; f < nxvgfile; f++)
3659 {
3660 read_bar_xvg(fxvgnms[f], &temp, &sim_data);
3661 nf++;
3662 }
3663 /* then .edr files */
3664 for (f = 0; f < nedrfile; f++)
3665 {
3666 read_barsim_edr(fedrnms[f], &temp, &sim_data);;
3667 nf++;
3668 }
3669
3670 /* fix the times to allow for equilibration */
3671 sim_data_impose_times(&sim_data, begin, end);
3672
3673 if (opt2bSet("-oh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
3674 {
3675 sim_data_histogram(&sim_data, opt2fn("-oh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nbin, oenv);
3676 }
3677
3678 /* assemble the output structures from the lambdas */
3679 results = barres_list_create(&sim_data, &nresults, use_dhdl);
3680
3681 sum_disc_err = barres_list_max_disc_err(results, nresults);
3682
3683 if (nresults == 0)
3684 {
3685 printf("\nNo results to calculate.\n");
3686 return 0;
3687 }
3688
3689 if (sum_disc_err > prec)
3690 {
3691 prec = sum_disc_err;
3692 nd = ceil(-log10(prec));
3693 printf("WARNING: setting the precision to %g because that is the minimum\n reasonable number, given the expected discretization error.\n", prec);
3694 }
3695
3696
3697 /*sprintf(lamformat,"%%6.3f");*/
3698 sprintf( dgformat, "%%%d.%df", 3+nd, nd);
3699 /* the format strings of the results in kT */
3700 sprintf( ktformat, "%%%d.%df", 5+nd, nd);
3701 sprintf( sktformat, "%%%ds", 6+nd);
3702 /* the format strings of the errors in kT */
3703 sprintf( kteformat, "%%%d.%df", 3+nd, nd);
3704 sprintf( skteformat, "%%%ds", 4+nd);
3705 sprintf(xvg2format, "%s %s\n", "%s", dgformat);
3706 sprintf(xvg3format, "%s %s %s\n", "%s", dgformat, dgformat);
3707
3708
3709
3710 fpb = NULL((void*)0);
3711 if (opt2bSet("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
3712 {
3713 sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
3714 fpb = xvgropen_type(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Free energy differences",
3715 "\\lambda", buf, exvggtXYDY, oenv);
3716 }
3717
3718 fpi = NULL((void*)0);
3719 if (opt2bSet("-oi", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
3720 {
3721 sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
3722 fpi = xvgropen(opt2fn("-oi", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Free energy integral",
3723 "\\lambda", buf, oenv);
3724 }
3725
3726
3727
3728 if (nbmin > nbmax)
3729 {
3730 nbmin = nbmax;
3731 }
3732
3733 /* first calculate results */
3734 bEE = TRUE1;
3735 disc_err = FALSE0;
3736 for (f = 0; f < nresults; f++)
3737 {
3738 /* Determine the free energy difference with a factor of 10
3739 * more accuracy than requested for printing.
3740 */
3741 calc_bar(&(results[f]), 0.1*prec, nbmin, nbmax,
3742 &bEE, partsum);
3743
3744 if (results[f].dg_disc_err > prec/10.)
3745 {
3746 disc_err = TRUE1;
3747 }
3748 if (results[f].dg_histrange_err > prec/10.)
3749 {
3750 histrange_err = TRUE1;
3751 }
3752 }
3753
3754 /* print results in kT */
3755 kT = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp;
3756 beta = 1/kT;
Value stored to 'beta' is never read
3757
3758 printf("\nTemperature: %g K\n", temp);
3759
3760 printf("\nDetailed results in kT (see help for explanation):\n\n");
3761 printf("%6s ", " lam_A");
3762 printf("%6s ", " lam_B");
3763 printf(sktformat, "DG ");
3764 if (bEE)
3765 {
3766 printf(skteformat, "+/- ");
3767 }
3768 if (disc_err)
3769 {
3770 printf(skteformat, "disc ");
3771 }
3772 if (histrange_err)
3773 {
3774 printf(skteformat, "range ");
3775 }
3776 printf(sktformat, "s_A ");
3777 if (bEE)
3778 {
3779 printf(skteformat, "+/- " );
3780 }
3781 printf(sktformat, "s_B ");
3782 if (bEE)
3783 {
3784 printf(skteformat, "+/- " );
3785 }
3786 printf(sktformat, "stdev ");
3787 if (bEE)
3788 {
3789 printf(skteformat, "+/- ");
3790 }
3791 printf("\n");
3792 for (f = 0; f < nresults; f++)
3793 {
3794 lambda_vec_print_short(results[f].a->native_lambda, buf);
3795 printf("%s ", buf);
3796 lambda_vec_print_short(results[f].b->native_lambda, buf);
3797 printf("%s ", buf);
3798 printf(ktformat, results[f].dg);
3799 printf(" ");
3800 if (bEE)
3801 {
3802 printf(kteformat, results[f].dg_err);
3803 printf(" ");
3804 }
3805 if (disc_err)
3806 {
3807 printf(kteformat, results[f].dg_disc_err);
3808 printf(" ");
3809 }
3810 if (histrange_err)
3811 {
3812 printf(kteformat, results[f].dg_histrange_err);
3813 printf(" ");
3814 }
3815 printf(ktformat, results[f].sa);
3816 printf(" ");
3817 if (bEE)
3818 {
3819 printf(kteformat, results[f].sa_err);
3820 printf(" ");
3821 }
3822 printf(ktformat, results[f].sb);
3823 printf(" ");
3824 if (bEE)
3825 {
3826 printf(kteformat, results[f].sb_err);
3827 printf(" ");
3828 }
3829 printf(ktformat, results[f].dg_stddev);
3830 printf(" ");
3831 if (bEE)
3832 {
3833 printf(kteformat, results[f].dg_stddev_err);
3834 }
3835 printf("\n");
3836
3837 /* Check for negative relative entropy with a 95% certainty. */
3838 if (results[f].sa < -2*results[f].sa_err ||
3839 results[f].sb < -2*results[f].sb_err)
3840 {
3841 result_OK = FALSE0;
3842 }
3843 }
3844
3845 if (!result_OK)
3846 {
3847 printf("\nWARNING: Some of these results violate the Second Law of "
3848 "Thermodynamics: \n"
3849 " This is can be the result of severe undersampling, or "
3850 "(more likely)\n"
3851 " there is something wrong with the simulations.\n");
3852 }
3853
3854
3855 /* final results in kJ/mol */
3856 printf("\n\nFinal results in kJ/mol:\n\n");
3857 dg_tot = 0;
3858 for (f = 0; f < nresults; f++)
3859 {
3860
3861 if (fpi != NULL((void*)0))
3862 {
3863 lambda_vec_print_short(results[f].a->native_lambda, buf);
3864 fprintf(fpi, xvg2format, buf, dg_tot);
3865 }
3866
3867
3868 if (fpb != NULL((void*)0))
3869 {
3870 lambda_vec_print_intermediate(results[f].a->native_lambda,
3871 results[f].b->native_lambda,
3872 buf);
3873
3874 fprintf(fpb, xvg3format, buf, results[f].dg, results[f].dg_err);
3875 }
3876
3877 printf("point ");
3878 lambda_vec_print_short(results[f].a->native_lambda, buf);
3879 lambda_vec_print_short(results[f].b->native_lambda, buf2);
3880 printf("%s - %s", buf, buf2);
3881 printf(", DG ");
3882
3883 printf(dgformat, results[f].dg*kT);
3884 if (bEE)
3885 {
3886 printf(" +/- ");
3887 printf(dgformat, results[f].dg_err*kT);
3888 }
3889 if (histrange_err)
3890 {
3891 printf(" (max. range err. = ");
3892 printf(dgformat, results[f].dg_histrange_err*kT);
3893 printf(")");
3894 sum_histrange_err += results[f].dg_histrange_err*kT;
3895 }
3896
3897 printf("\n");
3898 dg_tot += results[f].dg;
3899 }
3900 printf("\n");
3901 printf("total ");
3902 lambda_vec_print_short(results[0].a->native_lambda, buf);
3903 lambda_vec_print_short(results[nresults-1].b->native_lambda, buf2);
3904 printf("%s - %s", buf, buf2);
3905 printf(", DG ");
3906
3907 printf(dgformat, dg_tot*kT);
3908 if (bEE)
3909 {
3910 stat_err = bar_err(nbmin, nbmax, partsum)*kT;
3911 printf(" +/- ");
3912 printf(dgformat, max(max(stat_err, sum_disc_err), sum_histrange_err)((((((stat_err) > (sum_disc_err)) ? (stat_err) : (sum_disc_err
) )) > (sum_histrange_err)) ? ((((stat_err) > (sum_disc_err
)) ? (stat_err) : (sum_disc_err) )) : (sum_histrange_err) ));
3913 }
3914 printf("\n");
3915 if (disc_err)
3916 {
3917 printf("\nmaximum discretization error = ");
3918 printf(dgformat, sum_disc_err);
3919 if (bEE && stat_err < sum_disc_err)
3920 {
3921 printf("WARNING: discretization error (%g) is larger than statistical error.\n Decrease histogram spacing for more accurate results\n", stat_err);
3922 }
3923 }
3924 if (histrange_err)
3925 {
3926 printf("\nmaximum histogram range error = ");
3927 printf(dgformat, sum_histrange_err);
3928 if (bEE && stat_err < sum_histrange_err)
3929 {
3930 printf("WARNING: histogram range error (%g) is larger than statistical error.\n Increase histogram range for more accurate results\n", stat_err);
3931 }
3932
3933 }
3934 printf("\n");
3935
3936
3937 if (fpi != NULL((void*)0))
3938 {
3939 lambda_vec_print_short(results[nresults-1].b->native_lambda, buf);
3940 fprintf(fpi, xvg2format, buf, dg_tot);
3941 gmx_ffclose(fpi);
3942 }
3943
3944 do_view(oenv, opt2fn_null("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-xydy");
3945 do_view(oenv, opt2fn_null("-oi", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-xydy");
3946
3947 return 0;
3948}