File: | gromacs/gmxana/gmx_bar.c |
Location: | line 350, column 13 |
Description: | Value stored to 'str' is never read |
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, ")"); |
Value stored to 'str' is never read | |
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 | } |