File: | gromacs/gmxana/gmx_bar.c |
Location: | line 2370, column 30 |
Description: | Access to field 'N' results in a dereference of a null pointer (loaded from variable 'lc_in') |
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 */ | |||
63 | typedef 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 */ | |||
72 | typedef 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 */ | |||
85 | typedef 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 | ||||
104 | typedef 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 */ | |||
122 | typedef 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 */ | |||
149 | typedef 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) */ | |||
161 | typedef 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 */ | |||
180 | typedef 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 */ | |||
193 | typedef 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. */ | |||
203 | typedef 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 */ | |||
223 | static 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 */ | |||
231 | static 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 */ | |||
247 | static 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 */ | |||
278 | static 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 */ | |||
297 | static 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 | ||||
305 | static 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 | ||||
310 | static 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 */ | |||
324 | static 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 */ | |||
365 | static 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 */ | |||
388 | static 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 */ | |||
413 | static 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 */ | |||
436 | static 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 */ | |||
461 | static 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.*/ | |||
492 | static 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.*/ | |||
546 | static 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 | ||||
586 | static 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 | ||||
605 | static 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 | ||||
611 | static 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 | ||||
620 | static 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 | ||||
644 | static 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 | ||||
652 | static 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 | ||||
668 | static 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 | ||||
676 | static 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 | ||||
692 | static 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 */ | |||
709 | static 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 */ | |||
733 | static 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 */ | |||
751 | static 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 */ | |||
770 | static 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. */ | |||
790 | static 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 */ | |||
826 | static 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 */ | |||
868 | static 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 */ | |||
1032 | void 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 */ | |||
1137 | static 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 */ | |||
1213 | static 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 */ | |||
1257 | static 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 | ||||
1315 | static 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 */ | |||
1415 | static 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 */ | |||
1537 | static 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 */ | |||
1587 | static 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 | ||||
1598 | static 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 */ | |||
1618 | static 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 | ||||
1652 | static 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 | ||||
1789 | static 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 | ||||
1905 | static 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 | ||||
2028 | static 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 | ||||
2193 | static 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 | */ | |||
2223 | static 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 */ | |||
2257 | static 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 */ | |||
2407 | static 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 */ | |||
2416 | static 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 | */ | |||
2433 | static 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 | ||||
2570 | static 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 | ||||
2680 | static 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 | ||||
2731 | static 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 | ||||
2861 | static 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 | ||||
2912 | static 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 | ||||
3019 | static 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 | ||||
3166 | static 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 | ||||
3458 | int 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; | |||
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 | } |