File: | gromacs/gmxana/gmx_energy.c |
Location: | line 1771, column 16 |
Description: | Value stored to 'nnames' during its initialization is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2004, The GROMACS development team. |
6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
8 | * and including many others, as listed in the AUTHORS file in the |
9 | * top-level source directory and at http://www.gromacs.org. |
10 | * |
11 | * GROMACS is free software; you can redistribute it and/or |
12 | * modify it under the terms of the GNU Lesser General Public License |
13 | * as published by the Free Software Foundation; either version 2.1 |
14 | * of the License, or (at your option) any later version. |
15 | * |
16 | * GROMACS is distributed in the hope that it will be useful, |
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
19 | * Lesser General Public License for more details. |
20 | * |
21 | * You should have received a copy of the GNU Lesser General Public |
22 | * License along with GROMACS; if not, see |
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
25 | * |
26 | * If you want to redistribute modifications to GROMACS, please |
27 | * consider that scientific software is very special. Version |
28 | * control is crucial - bugs must be traceable. We will be happy to |
29 | * consider code for inclusion in the official distribution, but |
30 | * derived work must not be called official GROMACS. Details are found |
31 | * in the README & COPYING files - if they are missing, get the |
32 | * official version at http://www.gromacs.org. |
33 | * |
34 | * To help us fund GROMACS development, we humbly ask that you cite |
35 | * the research papers on the package. Check out http://www.gromacs.org. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <math.h> |
42 | #include <stdlib.h> |
43 | #include <string.h> |
44 | |
45 | #include "typedefs.h" |
46 | #include "gromacs/utility/fatalerror.h" |
47 | #include "gromacs/math/vec.h" |
48 | #include "gromacs/utility/cstringutil.h" |
49 | #include "gromacs/utility/smalloc.h" |
50 | #include "gromacs/fileio/enxio.h" |
51 | #include "gromacs/commandline/pargs.h" |
52 | #include "names.h" |
53 | #include "copyrite.h" |
54 | #include "macros.h" |
55 | #include "gromacs/fileio/xvgr.h" |
56 | #include "gstat.h" |
57 | #include "physics.h" |
58 | #include "gromacs/fileio/tpxio.h" |
59 | #include "gromacs/fileio/trxio.h" |
60 | #include "viewit.h" |
61 | #include "mtop_util.h" |
62 | #include "gmx_ana.h" |
63 | #include "mdebin.h" |
64 | |
65 | static real minthird = -1.0/3.0, minsixth = -1.0/6.0; |
66 | |
67 | typedef struct { |
68 | real sum; |
69 | real sum2; |
70 | } exactsum_t; |
71 | |
72 | typedef struct { |
73 | real *ener; |
74 | exactsum_t *es; |
75 | gmx_bool bExactStat; |
76 | double av; |
77 | double rmsd; |
78 | double ee; |
79 | double slope; |
80 | } enerdat_t; |
81 | |
82 | typedef struct { |
83 | gmx_int64_t nsteps; |
84 | gmx_int64_t npoints; |
85 | int nframes; |
86 | int *step; |
87 | int *steps; |
88 | int *points; |
89 | enerdat_t *s; |
90 | } enerdata_t; |
91 | |
92 | static double mypow(double x, double y) |
93 | { |
94 | if (x > 0) |
95 | { |
96 | return pow(x, y); |
97 | } |
98 | else |
99 | { |
100 | return 0.0; |
101 | } |
102 | } |
103 | |
104 | static int *select_it(int nre, char *nm[], int *nset) |
105 | { |
106 | gmx_bool *bE; |
107 | int n, k, j, i; |
108 | int *set; |
109 | gmx_bool bVerbose = TRUE1; |
110 | |
111 | if ((getenv("VERBOSE")) != NULL((void*)0)) |
112 | { |
113 | bVerbose = FALSE0; |
114 | } |
115 | |
116 | fprintf(stderrstderr, "Select the terms you want from the following list\n"); |
117 | fprintf(stderrstderr, "End your selection with 0\n"); |
118 | |
119 | if (bVerbose) |
120 | { |
121 | for (k = 0; (k < nre); ) |
122 | { |
123 | for (j = 0; (j < 4) && (k < nre); j++, k++) |
124 | { |
125 | fprintf(stderrstderr, " %3d=%14s", k+1, nm[k]); |
126 | } |
127 | fprintf(stderrstderr, "\n"); |
128 | } |
129 | } |
130 | |
131 | snew(bE, nre)(bE) = save_calloc("bE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 131, (nre), sizeof(*(bE))); |
132 | do |
133 | { |
134 | if (1 != scanf("%d", &n)) |
135 | { |
136 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 136, "Error reading user input"); |
137 | } |
138 | if ((n > 0) && (n <= nre)) |
139 | { |
140 | bE[n-1] = TRUE1; |
141 | } |
142 | } |
143 | while (n != 0); |
144 | |
145 | snew(set, nre)(set) = save_calloc("set", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 145, (nre), sizeof(*(set))); |
146 | for (i = (*nset) = 0; (i < nre); i++) |
147 | { |
148 | if (bE[i]) |
149 | { |
150 | set[(*nset)++] = i; |
151 | } |
152 | } |
153 | |
154 | sfree(bE)save_free("bE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 154, (bE)); |
155 | |
156 | return set; |
157 | } |
158 | |
159 | static void chomp(char *buf) |
160 | { |
161 | int len = strlen(buf); |
162 | |
163 | while ((len > 0) && (buf[len-1] == '\n')) |
164 | { |
165 | buf[len-1] = '\0'; |
166 | len--; |
167 | } |
168 | } |
169 | |
170 | static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) |
171 | { |
172 | gmx_bool *bE; |
173 | int n, k, kk, j, i, nmatch, nind, nss; |
174 | int *set; |
175 | gmx_bool bEOF, bVerbose = TRUE1, bLong = FALSE0; |
176 | char *ptr, buf[STRLEN4096]; |
177 | const char *fm4 = "%3d %-14s"; |
178 | const char *fm2 = "%3d %-34s"; |
179 | char **newnm = NULL((void*)0); |
180 | |
181 | if ((getenv("VERBOSE")) != NULL((void*)0)) |
182 | { |
183 | bVerbose = FALSE0; |
184 | } |
185 | |
186 | fprintf(stderrstderr, "\n"); |
187 | fprintf(stderrstderr, "Select the terms you want from the following list by\n"); |
188 | fprintf(stderrstderr, "selecting either (part of) the name or the number or a combination.\n"); |
189 | fprintf(stderrstderr, "End your selection with an empty line or a zero.\n"); |
190 | fprintf(stderrstderr, "-------------------------------------------------------------------\n"); |
191 | |
192 | snew(newnm, nre)(newnm) = save_calloc("newnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 192, (nre), sizeof(*(newnm))); |
193 | j = 0; |
194 | for (k = 0; k < nre; k++) |
195 | { |
196 | newnm[k] = strdup(nm[k].name)(__extension__ (__builtin_constant_p (nm[k].name) && ( (size_t)(const void *)((nm[k].name) + 1) - (size_t)(const void *)(nm[k].name) == 1) ? (((const char *) (nm[k].name))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (nm[k].name) + 1; char *__retval = (char *) malloc ( __len); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval, nm[k].name, __len); __retval; })) : __strdup (nm[ k].name))); |
197 | /* Insert dashes in all the names */ |
198 | while ((ptr = strchr(newnm[k], ' ')(__extension__ (__builtin_constant_p (' ') && !__builtin_constant_p (newnm[k]) && (' ') == '\0' ? (char *) __rawmemchr ( newnm[k], ' ') : __builtin_strchr (newnm[k], ' ')))) != NULL((void*)0)) |
199 | { |
200 | *ptr = '-'; |
201 | } |
202 | if (bVerbose) |
203 | { |
204 | if (j == 0) |
205 | { |
206 | if (k > 0) |
207 | { |
208 | fprintf(stderrstderr, "\n"); |
209 | } |
210 | bLong = FALSE0; |
211 | for (kk = k; kk < k+4; kk++) |
212 | { |
213 | if (kk < nre && strlen(nm[kk].name) > 14) |
214 | { |
215 | bLong = TRUE1; |
216 | } |
217 | } |
218 | } |
219 | else |
220 | { |
221 | fprintf(stderrstderr, " "); |
222 | } |
223 | if (!bLong) |
224 | { |
225 | fprintf(stderrstderr, fm4, k+1, newnm[k]); |
226 | j++; |
227 | if (j == 4) |
228 | { |
229 | j = 0; |
230 | } |
231 | } |
232 | else |
233 | { |
234 | fprintf(stderrstderr, fm2, k+1, newnm[k]); |
235 | j++; |
236 | if (j == 2) |
237 | { |
238 | j = 0; |
239 | } |
240 | } |
241 | } |
242 | } |
243 | if (bVerbose) |
244 | { |
245 | fprintf(stderrstderr, "\n\n"); |
246 | } |
247 | |
248 | snew(bE, nre)(bE) = save_calloc("bE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 248, (nre), sizeof(*(bE))); |
249 | |
250 | bEOF = FALSE0; |
251 | while (!bEOF && (fgets2(buf, STRLEN4096-1, stdinstdin))) |
252 | { |
253 | /* Remove newlines */ |
254 | chomp(buf); |
255 | |
256 | /* Remove spaces */ |
257 | trim(buf); |
258 | |
259 | /* Empty line means end of input */ |
260 | bEOF = (strlen(buf) == 0); |
261 | if (!bEOF) |
262 | { |
263 | ptr = buf; |
264 | do |
265 | { |
266 | if (!bEOF) |
267 | { |
268 | /* First try to read an integer */ |
269 | nss = sscanf(ptr, "%d", &nind); |
270 | if (nss == 1) |
271 | { |
272 | /* Zero means end of input */ |
273 | if (nind == 0) |
274 | { |
275 | bEOF = TRUE1; |
276 | } |
277 | else if ((1 <= nind) && (nind <= nre)) |
278 | { |
279 | bE[nind-1] = TRUE1; |
280 | } |
281 | else |
282 | { |
283 | fprintf(stderrstderr, "number %d is out of range\n", nind); |
284 | } |
285 | } |
286 | else |
287 | { |
288 | /* Now try to read a string */ |
289 | i = strlen(ptr); |
290 | nmatch = 0; |
291 | for (nind = 0; nind < nre; nind++) |
292 | { |
293 | if (gmx_strcasecmp(newnm[nind], ptr) == 0) |
294 | { |
295 | bE[nind] = TRUE1; |
296 | nmatch++; |
297 | } |
298 | } |
299 | if (nmatch == 0) |
300 | { |
301 | i = strlen(ptr); |
302 | nmatch = 0; |
303 | for (nind = 0; nind < nre; nind++) |
304 | { |
305 | if (gmx_strncasecmp(newnm[nind], ptr, i) == 0) |
306 | { |
307 | bE[nind] = TRUE1; |
308 | nmatch++; |
309 | } |
310 | } |
311 | if (nmatch == 0) |
312 | { |
313 | fprintf(stderrstderr, "String '%s' does not match anything\n", ptr); |
314 | } |
315 | } |
316 | } |
317 | } |
318 | /* Look for the first space, and remove spaces from there */ |
319 | if ((ptr = strchr(ptr, ' ')(__extension__ (__builtin_constant_p (' ') && !__builtin_constant_p (ptr) && (' ') == '\0' ? (char *) __rawmemchr (ptr, ' ' ) : __builtin_strchr (ptr, ' ')))) != NULL((void*)0)) |
320 | { |
321 | trim(ptr); |
322 | } |
323 | } |
324 | while (!bEOF && (ptr && (strlen(ptr) > 0))); |
325 | } |
326 | } |
327 | |
328 | snew(set, nre)(set) = save_calloc("set", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 328, (nre), sizeof(*(set))); |
329 | for (i = (*nset) = 0; (i < nre); i++) |
330 | { |
331 | if (bE[i]) |
332 | { |
333 | set[(*nset)++] = i; |
334 | } |
335 | } |
336 | |
337 | sfree(bE)save_free("bE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 337, (bE)); |
338 | |
339 | if (*nset == 0) |
340 | { |
341 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 341, "No energy terms selected"); |
342 | } |
343 | |
344 | for (i = 0; (i < nre); i++) |
345 | { |
346 | sfree(newnm[i])save_free("newnm[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 346, (newnm[i])); |
347 | } |
348 | sfree(newnm)save_free("newnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 348, (newnm)); |
349 | |
350 | return set; |
351 | } |
352 | |
353 | static void get_dhdl_parms(const char *topnm, t_inputrec *ir) |
354 | { |
355 | gmx_mtop_t mtop; |
356 | int natoms; |
357 | t_iatom *iatom; |
358 | matrix box; |
359 | |
360 | /* all we need is the ir to be able to write the label */ |
361 | read_tpx(topnm, ir, box, &natoms, NULL((void*)0), NULL((void*)0), NULL((void*)0), &mtop); |
362 | } |
363 | |
364 | static void get_orires_parms(const char *topnm, |
365 | int *nor, int *nex, int **label, real **obs) |
366 | { |
367 | gmx_mtop_t mtop; |
368 | gmx_localtop_t *top; |
369 | t_inputrec ir; |
370 | t_iparams *ip; |
371 | int natoms, i; |
372 | t_iatom *iatom; |
373 | int nb; |
374 | matrix box; |
375 | |
376 | read_tpx(topnm, &ir, box, &natoms, NULL((void*)0), NULL((void*)0), NULL((void*)0), &mtop); |
377 | top = gmx_mtop_generate_local_top(&mtop, &ir); |
378 | |
379 | ip = top->idef.iparams; |
380 | iatom = top->idef.il[F_ORIRES].iatoms; |
381 | |
382 | /* Count how many distance restraint there are... */ |
383 | nb = top->idef.il[F_ORIRES].nr; |
384 | if (nb == 0) |
385 | { |
386 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 386, "No orientation restraints in topology!\n"); |
387 | } |
388 | |
389 | *nor = nb/3; |
390 | *nex = 0; |
391 | snew(*label, *nor)(*label) = save_calloc("*label", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 391, (*nor), sizeof(*(*label))); |
392 | snew(*obs, *nor)(*obs) = save_calloc("*obs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 392, (*nor), sizeof(*(*obs))); |
393 | for (i = 0; i < nb; i += 3) |
394 | { |
395 | (*label)[i/3] = ip[iatom[i]].orires.label; |
396 | (*obs)[i/3] = ip[iatom[i]].orires.obs; |
397 | if (ip[iatom[i]].orires.ex >= *nex) |
398 | { |
399 | *nex = ip[iatom[i]].orires.ex+1; |
400 | } |
401 | } |
402 | fprintf(stderrstderr, "Found %d orientation restraints with %d experiments", |
403 | *nor, *nex); |
404 | } |
405 | |
406 | static int get_bounds(const char *topnm, |
407 | real **bounds, int **index, int **dr_pair, int *npairs, |
408 | gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir) |
409 | { |
410 | gmx_localtop_t *top; |
411 | t_functype *functype; |
412 | t_iparams *ip; |
413 | int natoms, i, j, k, type, ftype, natom; |
414 | t_ilist *disres; |
415 | t_iatom *iatom; |
416 | real *b; |
417 | int *ind, *pair; |
418 | int nb, label1; |
419 | matrix box; |
420 | |
421 | read_tpx(topnm, ir, box, &natoms, NULL((void*)0), NULL((void*)0), NULL((void*)0), mtop); |
422 | snew(*ltop, 1)(*ltop) = save_calloc("*ltop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 422, (1), sizeof(*(*ltop))); |
423 | top = gmx_mtop_generate_local_top(mtop, ir); |
424 | *ltop = top; |
425 | |
426 | functype = top->idef.functype; |
427 | ip = top->idef.iparams; |
428 | |
429 | /* Count how many distance restraint there are... */ |
430 | nb = top->idef.il[F_DISRES].nr; |
431 | if (nb == 0) |
432 | { |
433 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 433, "No distance restraints in topology!\n"); |
434 | } |
435 | |
436 | /* Allocate memory */ |
437 | snew(b, nb)(b) = save_calloc("b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 437, (nb), sizeof(*(b))); |
438 | snew(ind, nb)(ind) = save_calloc("ind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 438, (nb), sizeof(*(ind))); |
439 | snew(pair, nb+1)(pair) = save_calloc("pair", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 439, (nb+1), sizeof(*(pair))); |
440 | |
441 | /* Fill the bound array */ |
442 | nb = 0; |
443 | for (i = 0; (i < top->idef.ntypes); i++) |
444 | { |
445 | ftype = functype[i]; |
446 | if (ftype == F_DISRES) |
447 | { |
448 | |
449 | label1 = ip[i].disres.label; |
450 | b[nb] = ip[i].disres.up1; |
451 | ind[nb] = label1; |
452 | nb++; |
453 | } |
454 | } |
455 | *bounds = b; |
456 | |
457 | /* Fill the index array */ |
458 | label1 = -1; |
459 | disres = &(top->idef.il[F_DISRES]); |
460 | iatom = disres->iatoms; |
461 | for (i = j = k = 0; (i < disres->nr); ) |
462 | { |
463 | type = iatom[i]; |
464 | ftype = top->idef.functype[type]; |
465 | natom = interaction_function[ftype].nratoms+1; |
466 | if (label1 != top->idef.iparams[type].disres.label) |
467 | { |
468 | pair[j] = k; |
469 | label1 = top->idef.iparams[type].disres.label; |
470 | j++; |
471 | } |
472 | k++; |
473 | i += natom; |
474 | } |
475 | pair[j] = k; |
476 | *npairs = k; |
477 | if (j != nb) |
478 | { |
479 | gmx_incons("get_bounds for distance restraints")_gmx_error("incons", "get_bounds for distance restraints", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 479); |
480 | } |
481 | |
482 | *index = ind; |
483 | *dr_pair = pair; |
484 | |
485 | return nb; |
486 | } |
487 | |
488 | static void calc_violations(real rt[], real rav3[], int nb, int index[], |
489 | real bounds[], real *viol, double *st, double *sa) |
490 | { |
491 | const real sixth = 1.0/6.0; |
492 | int i, j; |
493 | double rsum, rav, sumaver, sumt; |
494 | |
495 | sumaver = 0; |
496 | sumt = 0; |
497 | for (i = 0; (i < nb); i++) |
498 | { |
499 | rsum = 0.0; |
500 | rav = 0.0; |
501 | for (j = index[i]; (j < index[i+1]); j++) |
502 | { |
503 | if (viol) |
504 | { |
505 | viol[j] += mypow(rt[j], -3.0); |
506 | } |
507 | rav += sqr(rav3[j]); |
508 | rsum += mypow(rt[j], -6); |
509 | } |
510 | rsum = max(0.0, mypow(rsum, -sixth)-bounds[i])(((0.0) > (mypow(rsum, -sixth)-bounds[i])) ? (0.0) : (mypow (rsum, -sixth)-bounds[i]) ); |
511 | rav = max(0.0, mypow(rav, -sixth)-bounds[i])(((0.0) > (mypow(rav, -sixth)-bounds[i])) ? (0.0) : (mypow (rav, -sixth)-bounds[i]) ); |
512 | |
513 | sumt += rsum; |
514 | sumaver += rav; |
515 | } |
516 | *st = sumt; |
517 | *sa = sumaver; |
518 | } |
519 | |
520 | static void analyse_disre(const char *voutfn, int nframes, |
521 | real violaver[], real bounds[], int index[], |
522 | int pair[], int nbounds, |
523 | const output_env_t oenv) |
524 | { |
525 | FILE *vout; |
526 | double sum, sumt, sumaver; |
527 | int i, j; |
528 | |
529 | /* Subtract bounds from distances, to calculate violations */ |
530 | calc_violations(violaver, violaver, |
531 | nbounds, pair, bounds, NULL((void*)0), &sumt, &sumaver); |
532 | |
533 | #ifdef DEBUG |
534 | fprintf(stdoutstdout, "\nSum of violations averaged over simulation: %g nm\n", |
535 | sumaver); |
536 | fprintf(stdoutstdout, "Largest violation averaged over simulation: %g nm\n\n", |
537 | sumt); |
538 | #endif |
539 | vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm", |
540 | oenv); |
541 | sum = 0.0; |
542 | sumt = 0.0; |
543 | for (i = 0; (i < nbounds); i++) |
544 | { |
545 | /* Do ensemble averaging */ |
546 | sumaver = 0; |
547 | for (j = pair[i]; (j < pair[i+1]); j++) |
548 | { |
549 | sumaver += sqr(violaver[j]/nframes); |
550 | } |
551 | sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i])(((0.0) > (mypow(sumaver, minsixth)-bounds[i])) ? (0.0) : ( mypow(sumaver, minsixth)-bounds[i]) ); |
552 | |
553 | sumt += sumaver; |
554 | sum = max(sum, sumaver)(((sum) > (sumaver)) ? (sum) : (sumaver) ); |
555 | fprintf(vout, "%10d %10.5e\n", index[i], sumaver); |
556 | } |
557 | #ifdef DEBUG |
558 | for (j = 0; (j < dr.ndr); j++) |
559 | { |
560 | fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird)); |
561 | } |
562 | #endif |
563 | gmx_ffclose(vout); |
564 | |
565 | fprintf(stdoutstdout, "\nSum of violations averaged over simulation: %g nm\n", |
566 | sumt); |
567 | fprintf(stdoutstdout, "Largest violation averaged over simulation: %g nm\n\n", sum); |
568 | |
569 | do_view(oenv, voutfn, "-graphtype bar"); |
570 | } |
571 | |
572 | static void einstein_visco(const char *fn, const char *fni, int nsets, |
573 | int nframes, real **sum, |
574 | real V, real T, int nsteps, double time[], |
575 | const output_env_t oenv) |
576 | { |
577 | FILE *fp0, *fp1; |
578 | double av[4], avold[4]; |
579 | double fac, dt, di; |
580 | int i, j, m, nf4; |
581 | |
582 | if (nframes < 1) |
583 | { |
584 | return; |
585 | } |
586 | |
587 | dt = (time[1]-time[0]); |
588 | nf4 = nframes/4+1; |
589 | |
590 | for (i = 0; i <= nsets; i++) |
591 | { |
592 | avold[i] = 0; |
593 | } |
594 | fp0 = xvgropen(fni, "Shear viscosity integral", |
595 | "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv); |
596 | fp1 = xvgropen(fn, "Shear viscosity using Einstein relation", |
597 | "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv); |
598 | for (i = 1; i < nf4; i++) |
599 | { |
600 | fac = dt*nframes/nsteps; |
601 | for (m = 0; m <= nsets; m++) |
602 | { |
603 | av[m] = 0; |
604 | } |
605 | for (j = 0; j < nframes-i; j++) |
606 | { |
607 | for (m = 0; m < nsets; m++) |
608 | { |
609 | di = sqr(fac*(sum[m][j+i]-sum[m][j])); |
610 | |
611 | av[m] += di; |
612 | av[nsets] += di/nsets; |
613 | } |
614 | } |
615 | /* Convert to SI for the viscosity */ |
616 | fac = (V*NANO(1e-9)*NANO(1e-9)*NANO(1e-9)*PICO(1e-12)*1e10)/(2*BOLTZMANN(1.380658e-23)*T)/(nframes-i); |
617 | fprintf(fp0, "%10g", time[i]-time[0]); |
618 | for (m = 0; (m <= nsets); m++) |
619 | { |
620 | av[m] = fac*av[m]; |
621 | fprintf(fp0, " %10g", av[m]); |
622 | } |
623 | fprintf(fp0, "\n"); |
624 | fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]); |
625 | for (m = 0; (m <= nsets); m++) |
626 | { |
627 | fprintf(fp1, " %10g", (av[m]-avold[m])/dt); |
628 | avold[m] = av[m]; |
629 | } |
630 | fprintf(fp1, "\n"); |
631 | } |
632 | gmx_ffclose(fp0); |
633 | gmx_ffclose(fp1); |
634 | } |
635 | |
636 | typedef struct { |
637 | gmx_int64_t np; |
638 | double sum; |
639 | double sav; |
640 | double sav2; |
641 | } ee_sum_t; |
642 | |
643 | typedef struct { |
644 | int b; |
645 | ee_sum_t sum; |
646 | gmx_int64_t nst; |
647 | gmx_int64_t nst_min; |
648 | } ener_ee_t; |
649 | |
650 | static void clear_ee_sum(ee_sum_t *ees) |
651 | { |
652 | ees->sav = 0; |
653 | ees->sav2 = 0; |
654 | ees->np = 0; |
655 | ees->sum = 0; |
656 | } |
657 | |
658 | static void add_ee_sum(ee_sum_t *ees, double sum, int np) |
659 | { |
660 | ees->np += np; |
661 | ees->sum += sum; |
662 | } |
663 | |
664 | static void add_ee_av(ee_sum_t *ees) |
665 | { |
666 | double av; |
667 | |
668 | av = ees->sum/ees->np; |
669 | ees->sav += av; |
670 | ees->sav2 += av*av; |
671 | ees->np = 0; |
672 | ees->sum = 0; |
673 | } |
674 | |
675 | static double calc_ee2(int nb, ee_sum_t *ees) |
676 | { |
677 | return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1); |
678 | } |
679 | |
680 | static void set_ee_av(ener_ee_t *eee) |
681 | { |
682 | if (debug) |
683 | { |
684 | char buf[STEPSTRSIZE22]; |
685 | fprintf(debug, "Storing average for err.est.: %s steps\n", |
686 | gmx_step_str(eee->nst, buf)); |
687 | } |
688 | add_ee_av(&eee->sum); |
689 | eee->b++; |
690 | if (eee->b == 1 || eee->nst < eee->nst_min) |
691 | { |
692 | eee->nst_min = eee->nst; |
693 | } |
694 | eee->nst = 0; |
695 | } |
696 | |
697 | static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax) |
698 | { |
699 | int nb, i, f, nee; |
700 | double sum, sum2, sump, see2; |
701 | gmx_int64_t steps, np, p, bound_nb; |
702 | enerdat_t *ed; |
703 | exactsum_t *es; |
704 | gmx_bool bAllZero; |
705 | double x, sx, sy, sxx, sxy; |
706 | ener_ee_t *eee; |
707 | |
708 | /* Check if we have exact statistics over all points */ |
709 | for (i = 0; i < nset; i++) |
710 | { |
711 | ed = &edat->s[i]; |
712 | ed->bExactStat = FALSE0; |
713 | if (edat->npoints > 0) |
714 | { |
715 | /* All energy file sum entries 0 signals no exact sums. |
716 | * But if all energy values are 0, we still have exact sums. |
717 | */ |
718 | bAllZero = TRUE1; |
719 | for (f = 0; f < edat->nframes && !ed->bExactStat; f++) |
720 | { |
721 | if (ed->ener[i] != 0) |
722 | { |
723 | bAllZero = FALSE0; |
724 | } |
725 | ed->bExactStat = (ed->es[f].sum != 0); |
726 | } |
727 | if (bAllZero) |
728 | { |
729 | ed->bExactStat = TRUE1; |
730 | } |
731 | } |
732 | } |
733 | |
734 | snew(eee, nbmax+1)(eee) = save_calloc("eee", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 734, (nbmax+1), sizeof(*(eee))); |
735 | for (i = 0; i < nset; i++) |
736 | { |
737 | ed = &edat->s[i]; |
738 | |
739 | sum = 0; |
740 | sum2 = 0; |
741 | np = 0; |
742 | sx = 0; |
743 | sy = 0; |
744 | sxx = 0; |
745 | sxy = 0; |
746 | for (nb = nbmin; nb <= nbmax; nb++) |
747 | { |
748 | eee[nb].b = 0; |
749 | clear_ee_sum(&eee[nb].sum); |
750 | eee[nb].nst = 0; |
751 | eee[nb].nst_min = 0; |
752 | } |
753 | for (f = 0; f < edat->nframes; f++) |
754 | { |
755 | es = &ed->es[f]; |
756 | |
757 | if (ed->bExactStat) |
758 | { |
759 | /* Add the sum and the sum of variances to the totals. */ |
760 | p = edat->points[f]; |
761 | sump = es->sum; |
762 | sum2 += es->sum2; |
763 | if (np > 0) |
764 | { |
765 | sum2 += dsqr(sum/np - (sum + es->sum)/(np + p)) |
766 | *np*(np + p)/p; |
767 | } |
768 | } |
769 | else |
770 | { |
771 | /* Add a single value to the sum and sum of squares. */ |
772 | p = 1; |
773 | sump = ed->ener[f]; |
774 | sum2 += dsqr(sump); |
775 | } |
776 | |
777 | /* sum has to be increased after sum2 */ |
778 | np += p; |
779 | sum += sump; |
780 | |
781 | /* For the linear regression use variance 1/p. |
782 | * Note that sump is the sum, not the average, so we don't need p*. |
783 | */ |
784 | x = edat->step[f] - 0.5*(edat->steps[f] - 1); |
785 | sx += p*x; |
786 | sy += sump; |
787 | sxx += p*x*x; |
788 | sxy += x*sump; |
789 | |
790 | for (nb = nbmin; nb <= nbmax; nb++) |
791 | { |
792 | /* Check if the current end step is closer to the desired |
793 | * block boundary than the next end step. |
794 | */ |
795 | bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1); |
796 | if (eee[nb].nst > 0 && |
797 | bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb) |
798 | { |
799 | set_ee_av(&eee[nb]); |
800 | } |
801 | if (f == 0) |
802 | { |
803 | eee[nb].nst = 1; |
804 | } |
805 | else |
806 | { |
807 | eee[nb].nst += edat->step[f] - edat->step[f-1]; |
808 | } |
809 | if (ed->bExactStat) |
810 | { |
811 | add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]); |
812 | } |
813 | else |
814 | { |
815 | add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1); |
816 | } |
817 | bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1); |
818 | if (edat->step[f]*nb >= bound_nb) |
819 | { |
820 | set_ee_av(&eee[nb]); |
821 | } |
822 | } |
823 | } |
824 | |
825 | edat->s[i].av = sum/np; |
826 | if (ed->bExactStat) |
827 | { |
828 | edat->s[i].rmsd = sqrt(sum2/np); |
829 | } |
830 | else |
831 | { |
832 | edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av)); |
833 | } |
834 | |
835 | if (edat->nframes > 1) |
836 | { |
837 | edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx); |
838 | } |
839 | else |
840 | { |
841 | edat->s[i].slope = 0; |
842 | } |
843 | |
844 | nee = 0; |
845 | see2 = 0; |
846 | for (nb = nbmin; nb <= nbmax; nb++) |
847 | { |
848 | /* Check if we actually got nb blocks and if the smallest |
849 | * block is not shorter than 80% of the average. |
850 | */ |
851 | if (debug) |
852 | { |
853 | char buf1[STEPSTRSIZE22], buf2[STEPSTRSIZE22]; |
854 | fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n", |
855 | nb, eee[nb].b, |
856 | gmx_step_str(eee[nb].nst_min, buf1), |
857 | gmx_step_str(edat->nsteps, buf2)); |
858 | } |
859 | if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps) |
860 | { |
861 | see2 += calc_ee2(nb, &eee[nb].sum); |
862 | nee++; |
863 | } |
864 | } |
865 | if (nee > 0) |
866 | { |
867 | edat->s[i].ee = sqrt(see2/nee); |
868 | } |
869 | else |
870 | { |
871 | edat->s[i].ee = -1; |
872 | } |
873 | } |
874 | sfree(eee)save_free("eee", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 874, (eee)); |
875 | } |
876 | |
877 | static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax) |
878 | { |
879 | enerdata_t *esum; |
880 | enerdat_t *s; |
881 | int f, i; |
882 | double sum; |
883 | |
884 | snew(esum, 1)(esum) = save_calloc("esum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 884, (1), sizeof(*(esum))); |
885 | *esum = *edat; |
886 | snew(esum->s, 1)(esum->s) = save_calloc("esum->s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 886, (1), sizeof(*(esum->s))); |
887 | s = &esum->s[0]; |
888 | snew(s->ener, esum->nframes)(s->ener) = save_calloc("s->ener", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 888, (esum->nframes), sizeof(*(s->ener))); |
889 | snew(s->es, esum->nframes)(s->es) = save_calloc("s->es", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 889, (esum->nframes), sizeof(*(s->es))); |
890 | |
891 | s->bExactStat = TRUE1; |
892 | s->slope = 0; |
893 | for (i = 0; i < nset; i++) |
894 | { |
895 | if (!edat->s[i].bExactStat) |
896 | { |
897 | s->bExactStat = FALSE0; |
898 | } |
899 | s->slope += edat->s[i].slope; |
900 | } |
901 | |
902 | for (f = 0; f < edat->nframes; f++) |
903 | { |
904 | sum = 0; |
905 | for (i = 0; i < nset; i++) |
906 | { |
907 | sum += edat->s[i].ener[f]; |
908 | } |
909 | s->ener[f] = sum; |
910 | sum = 0; |
911 | for (i = 0; i < nset; i++) |
912 | { |
913 | sum += edat->s[i].es[f].sum; |
914 | } |
915 | s->es[f].sum = sum; |
916 | s->es[f].sum2 = 0; |
917 | } |
918 | |
919 | calc_averages(1, esum, nbmin, nbmax); |
920 | |
921 | return esum; |
922 | } |
923 | |
924 | static char *ee_pr(double ee, char *buf) |
925 | { |
926 | char tmp[100]; |
927 | double rnd; |
928 | |
929 | if (ee < 0) |
930 | { |
931 | sprintf(buf, "%s", "--"); |
932 | } |
933 | else |
934 | { |
935 | /* Round to two decimals by printing. */ |
936 | sprintf(tmp, "%.1e", ee); |
937 | sscanf(tmp, "%lf", &rnd); |
938 | sprintf(buf, "%g", rnd); |
939 | } |
940 | |
941 | return buf; |
942 | } |
943 | |
944 | static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat) |
945 | { |
946 | /* Remove the drift by performing a fit to y = ax+b. |
947 | Uses 5 iterations. */ |
948 | int i, j, k; |
949 | double delta, d, sd, sd2; |
950 | |
951 | edat->npoints = edat->nframes; |
952 | edat->nsteps = edat->nframes; |
953 | |
954 | for (k = 0; (k < 5); k++) |
955 | { |
956 | for (i = 0; (i < nset); i++) |
957 | { |
958 | delta = edat->s[i].slope*dt; |
959 | |
960 | if (NULL((void*)0) != debug) |
961 | { |
962 | fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope); |
963 | } |
964 | |
965 | for (j = 0; (j < edat->nframes); j++) |
966 | { |
967 | edat->s[i].ener[j] -= j*delta; |
968 | edat->s[i].es[j].sum = 0; |
969 | edat->s[i].es[j].sum2 = 0; |
970 | } |
971 | } |
972 | calc_averages(nset, edat, nbmin, nbmax); |
973 | } |
974 | } |
975 | |
976 | static void calc_fluctuation_props(FILE *fp, |
977 | gmx_bool bDriftCorr, real dt, |
978 | int nset, int nmol, |
979 | char **leg, enerdata_t *edat, |
980 | int nbmin, int nbmax) |
981 | { |
982 | int i, j; |
983 | double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet; |
984 | double NANO3; |
985 | enum { |
986 | eVol, eEnth, eTemp, eEtot, eNR |
987 | }; |
988 | const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" }; |
989 | int ii[eNR]; |
990 | |
991 | NANO3 = NANO(1e-9)*NANO(1e-9)*NANO(1e-9); |
992 | if (!bDriftCorr) |
993 | { |
994 | fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n" |
995 | "for spurious drift in the graphs. Note that this is not\n" |
996 | "a substitute for proper equilibration and sampling!\n"); |
997 | } |
998 | else |
999 | { |
1000 | remove_drift(nset, nbmin, nbmax, dt, edat); |
1001 | } |
1002 | for (i = 0; (i < eNR); i++) |
1003 | { |
1004 | for (ii[i] = 0; (ii[i] < nset && |
1005 | (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++) |
1006 | { |
1007 | ; |
1008 | } |
1009 | /* if (ii[i] < nset) |
1010 | fprintf(fp,"Found %s data.\n",my_ener[i]); |
1011 | */ } |
1012 | /* Compute it all! */ |
1013 | alpha = kappa = cp = dcp = cv = NOTSET-12345; |
1014 | |
1015 | /* Temperature */ |
1016 | tt = NOTSET-12345; |
1017 | if (ii[eTemp] < nset) |
1018 | { |
1019 | tt = edat->s[ii[eTemp]].av; |
1020 | } |
1021 | /* Volume */ |
1022 | vv = varv = NOTSET-12345; |
1023 | if ((ii[eVol] < nset) && (ii[eTemp] < nset)) |
1024 | { |
1025 | vv = edat->s[ii[eVol]].av*NANO3; |
1026 | varv = dsqr(edat->s[ii[eVol]].rmsd*NANO3); |
1027 | kappa = (varv/vv)/(BOLTZMANN(1.380658e-23)*tt); |
1028 | } |
1029 | /* Enthalpy */ |
1030 | hh = varh = NOTSET-12345; |
1031 | if ((ii[eEnth] < nset) && (ii[eTemp] < nset)) |
1032 | { |
1033 | hh = KILO(1e3)*edat->s[ii[eEnth]].av/AVOGADRO(6.0221367e23); |
1034 | varh = dsqr(KILO(1e3)*edat->s[ii[eEnth]].rmsd/AVOGADRO(6.0221367e23)); |
1035 | cp = AVOGADRO(6.0221367e23)*((varh/nmol)/(BOLTZMANN(1.380658e-23)*tt*tt)); |
1036 | } |
1037 | /* Total energy */ |
1038 | et = varet = NOTSET-12345; |
1039 | if ((ii[eEtot] < nset) && (hh == NOTSET-12345) && (tt != NOTSET-12345)) |
1040 | { |
1041 | /* Only compute cv in constant volume runs, which we can test |
1042 | by checking whether the enthalpy was computed. |
1043 | */ |
1044 | et = edat->s[ii[eEtot]].av; |
1045 | varet = sqr(edat->s[ii[eEtot]].rmsd); |
1046 | cv = KILO(1e3)*((varet/nmol)/(BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*tt*tt)); |
1047 | } |
1048 | /* Alpha, dcp */ |
1049 | if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset)) |
1050 | { |
1051 | double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver; |
1052 | vh_sum = v_sum = h_sum = 0; |
1053 | for (j = 0; (j < edat->nframes); j++) |
1054 | { |
1055 | v = edat->s[ii[eVol]].ener[j]*NANO3; |
1056 | h = KILO(1e3)*edat->s[ii[eEnth]].ener[j]/AVOGADRO(6.0221367e23); |
1057 | v_sum += v; |
1058 | h_sum += h; |
1059 | vh_sum += (v*h); |
1060 | } |
1061 | vh_aver = vh_sum / edat->nframes; |
1062 | v_aver = v_sum / edat->nframes; |
1063 | h_aver = h_sum / edat->nframes; |
1064 | alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN(1.380658e-23)*tt*tt); |
1065 | dcp = (v_aver*AVOGADRO(6.0221367e23)/nmol)*tt*sqr(alpha)/(kappa); |
1066 | } |
1067 | |
1068 | if (tt != NOTSET-12345) |
1069 | { |
1070 | if (nmol < 2) |
1071 | { |
1072 | fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n", |
1073 | nmol); |
1074 | } |
1075 | fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt); |
1076 | fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n"); |
1077 | fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n"); |
1078 | fprintf(fp, "please use the g_dos program.\n\n"); |
1079 | fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n" |
1080 | "a block-averaging error analysis (not implemented in g_energy yet)\n"); |
1081 | |
1082 | if (debug != NULL((void*)0)) |
1083 | { |
1084 | if (varv != NOTSET-12345) |
1085 | { |
1086 | fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO(6.0221367e23)/nmol); |
1087 | } |
1088 | } |
1089 | if (vv != NOTSET-12345) |
1090 | { |
1091 | fprintf(fp, "Volume = %10g m^3/mol\n", |
1092 | vv*AVOGADRO(6.0221367e23)/nmol); |
1093 | } |
1094 | if (varh != NOTSET-12345) |
1095 | { |
1096 | fprintf(fp, "Enthalpy = %10g kJ/mol\n", |
1097 | hh*AVOGADRO(6.0221367e23)/(KILO(1e3)*nmol)); |
1098 | } |
1099 | if (alpha != NOTSET-12345) |
1100 | { |
1101 | fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n", |
1102 | alpha); |
1103 | } |
1104 | if (kappa != NOTSET-12345) |
1105 | { |
1106 | fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n", |
1107 | kappa); |
1108 | fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n", |
1109 | 1.0/kappa); |
1110 | } |
1111 | if (cp != NOTSET-12345) |
1112 | { |
1113 | fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n", |
1114 | cp); |
1115 | } |
1116 | if (cv != NOTSET-12345) |
1117 | { |
1118 | fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n", |
1119 | cv); |
1120 | } |
1121 | if (dcp != NOTSET-12345) |
1122 | { |
1123 | fprintf(fp, "Cp-Cv = %10g J/mol K\n", |
1124 | dcp); |
1125 | } |
1126 | please_cite(fp, "Allen1987a"); |
1127 | } |
1128 | else |
1129 | { |
1130 | fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n"); |
1131 | } |
1132 | } |
1133 | |
1134 | static void analyse_ener(gmx_bool bCorr, const char *corrfn, |
1135 | gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct, |
1136 | gmx_bool bVisco, const char *visfn, int nmol, |
1137 | gmx_int64_t start_step, double start_t, |
1138 | gmx_int64_t step, double t, |
1139 | double time[], real reftemp, |
1140 | enerdata_t *edat, |
1141 | int nset, int set[], gmx_bool *bIsEner, |
1142 | char **leg, gmx_enxnm_t *enm, |
1143 | real Vaver, real ezero, |
1144 | int nbmin, int nbmax, |
1145 | const output_env_t oenv) |
1146 | { |
1147 | FILE *fp; |
1148 | /* Check out the printed manual for equations! */ |
1149 | double Dt, aver, stddev, errest, delta_t, totaldrift; |
1150 | enerdata_t *esum = NULL((void*)0); |
1151 | real xxx, integral, intBulk, Temp = 0, Pres = 0; |
1152 | real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest; |
1153 | double beta = 0, expE, expEtot, *fee = NULL((void*)0); |
1154 | gmx_int64_t nsteps; |
1155 | int nexact, nnotexact; |
1156 | double x1m, x1mk; |
1157 | int i, j, k, nout; |
1158 | real chi2; |
1159 | char buf[256], eebuf[100]; |
1160 | |
1161 | nsteps = step - start_step + 1; |
1162 | if (nsteps < 1) |
1163 | { |
1164 | fprintf(stdoutstdout, "Not enough steps (%s) for statistics\n", |
1165 | gmx_step_str(nsteps, buf)); |
1166 | } |
1167 | else |
1168 | { |
1169 | /* Calculate the time difference */ |
1170 | delta_t = t - start_t; |
1171 | |
1172 | fprintf(stdoutstdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n", |
1173 | gmx_step_str(nsteps, buf), start_t, t, nset); |
1174 | |
1175 | calc_averages(nset, edat, nbmin, nbmax); |
1176 | |
1177 | if (bSum) |
1178 | { |
1179 | esum = calc_sum(nset, edat, nbmin, nbmax); |
1180 | } |
1181 | |
1182 | if (edat->npoints == 0) |
1183 | { |
1184 | nexact = 0; |
1185 | nnotexact = nset; |
1186 | } |
1187 | else |
1188 | { |
1189 | nexact = 0; |
1190 | nnotexact = 0; |
1191 | for (i = 0; (i < nset); i++) |
1192 | { |
1193 | if (edat->s[i].bExactStat) |
1194 | { |
1195 | nexact++; |
1196 | } |
1197 | else |
1198 | { |
1199 | nnotexact++; |
1200 | } |
1201 | } |
1202 | } |
1203 | |
1204 | if (nnotexact == 0) |
1205 | { |
1206 | fprintf(stdoutstdout, "All statistics are over %s points\n", |
1207 | gmx_step_str(edat->npoints, buf)); |
1208 | } |
1209 | else if (nexact == 0 || edat->npoints == edat->nframes) |
1210 | { |
1211 | fprintf(stdoutstdout, "All statistics are over %d points (frames)\n", |
1212 | edat->nframes); |
1213 | } |
1214 | else |
1215 | { |
1216 | fprintf(stdoutstdout, "The term%s", nnotexact == 1 ? "" : "s"); |
1217 | for (i = 0; (i < nset); i++) |
1218 | { |
1219 | if (!edat->s[i].bExactStat) |
1220 | { |
1221 | fprintf(stdoutstdout, " '%s'", leg[i]); |
1222 | } |
1223 | } |
1224 | fprintf(stdoutstdout, " %s has statistics over %d points (frames)\n", |
1225 | nnotexact == 1 ? "is" : "are", edat->nframes); |
1226 | fprintf(stdoutstdout, "All other statistics are over %s points\n", |
1227 | gmx_step_str(edat->npoints, buf)); |
1228 | } |
1229 | fprintf(stdoutstdout, "\n"); |
1230 | |
1231 | fprintf(stdoutstdout, "%-24s %10s %10s %10s %10s", |
1232 | "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift"); |
1233 | if (bFee) |
1234 | { |
1235 | fprintf(stdoutstdout, " %10s\n", "-kT ln<e^(E/kT)>"); |
1236 | } |
1237 | else |
1238 | { |
1239 | fprintf(stdoutstdout, "\n"); |
1240 | } |
1241 | fprintf(stdoutstdout, "-------------------------------------------------------------------------------\n"); |
1242 | |
1243 | /* Initiate locals, only used with -sum */ |
1244 | expEtot = 0; |
1245 | if (bFee) |
1246 | { |
1247 | beta = 1.0/(BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*reftemp); |
1248 | snew(fee, nset)(fee) = save_calloc("fee", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1248, (nset), sizeof(*(fee))); |
1249 | } |
1250 | for (i = 0; (i < nset); i++) |
1251 | { |
1252 | aver = edat->s[i].av; |
1253 | stddev = edat->s[i].rmsd; |
1254 | errest = edat->s[i].ee; |
1255 | |
1256 | if (bFee) |
1257 | { |
1258 | expE = 0; |
1259 | for (j = 0; (j < edat->nframes); j++) |
1260 | { |
1261 | expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol); |
1262 | } |
1263 | if (bSum) |
1264 | { |
1265 | expEtot += expE/edat->nframes; |
1266 | } |
1267 | |
1268 | fee[i] = log(expE/edat->nframes)/beta + aver/nmol; |
1269 | } |
1270 | if (strstr(leg[i], "empera") != NULL((void*)0)) |
1271 | { |
1272 | Temp = aver; |
1273 | } |
1274 | else if (strstr(leg[i], "olum") != NULL((void*)0)) |
1275 | { |
1276 | Vaver = aver; |
1277 | } |
1278 | else if (strstr(leg[i], "essure") != NULL((void*)0)) |
1279 | { |
1280 | Pres = aver; |
1281 | } |
1282 | if (bIsEner[i]) |
1283 | { |
1284 | pr_aver = aver/nmol-ezero; |
1285 | pr_stddev = stddev/nmol; |
1286 | pr_errest = errest/nmol; |
1287 | } |
1288 | else |
1289 | { |
1290 | pr_aver = aver; |
1291 | pr_stddev = stddev; |
1292 | pr_errest = errest; |
1293 | } |
1294 | |
1295 | /* Multiply the slope in steps with the number of steps taken */ |
1296 | totaldrift = (edat->nsteps - 1)*edat->s[i].slope; |
1297 | if (bIsEner[i]) |
1298 | { |
1299 | totaldrift /= nmol; |
1300 | } |
1301 | |
1302 | fprintf(stdoutstdout, "%-24s %10g %10s %10g %10g", |
1303 | leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift); |
1304 | if (bFee) |
1305 | { |
1306 | fprintf(stdoutstdout, " %10g", fee[i]); |
1307 | } |
1308 | |
1309 | fprintf(stdoutstdout, " (%s)\n", enm[set[i]].unit); |
1310 | |
1311 | if (bFluct) |
1312 | { |
1313 | for (j = 0; (j < edat->nframes); j++) |
1314 | { |
1315 | edat->s[i].ener[j] -= aver; |
1316 | } |
1317 | } |
1318 | } |
1319 | if (bSum) |
1320 | { |
1321 | totaldrift = (edat->nsteps - 1)*esum->s[0].slope; |
1322 | fprintf(stdoutstdout, "%-24s %10g %10s %10s %10g (%s)", |
1323 | "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf), |
1324 | "--", totaldrift/nmol, enm[set[0]].unit); |
1325 | /* pr_aver,pr_stddev,a,totaldrift */ |
1326 | if (bFee) |
1327 | { |
1328 | fprintf(stdoutstdout, " %10g %10g\n", |
1329 | log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta); |
1330 | } |
1331 | else |
1332 | { |
1333 | fprintf(stdoutstdout, "\n"); |
1334 | } |
1335 | } |
1336 | |
1337 | /* Do correlation function */ |
1338 | if (edat->nframes > 1) |
1339 | { |
1340 | Dt = delta_t/(edat->nframes - 1); |
1341 | } |
1342 | else |
1343 | { |
1344 | Dt = 0; |
1345 | } |
1346 | if (bVisco) |
1347 | { |
1348 | const char* leg[] = { "Shear", "Bulk" }; |
1349 | real factor; |
1350 | real **eneset; |
1351 | real **enesum; |
1352 | |
1353 | /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */ |
1354 | |
1355 | /* Symmetrise tensor! (and store in first three elements) |
1356 | * And subtract average pressure! |
1357 | */ |
1358 | snew(eneset, 12)(eneset) = save_calloc("eneset", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1358, (12), sizeof(*(eneset))); |
1359 | for (i = 0; i < 12; i++) |
1360 | { |
1361 | snew(eneset[i], edat->nframes)(eneset[i]) = save_calloc("eneset[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1361, (edat->nframes), sizeof(*(eneset[i]))); |
1362 | } |
1363 | snew(enesum, 3)(enesum) = save_calloc("enesum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1363, (3), sizeof(*(enesum))); |
1364 | for (i = 0; i < 3; i++) |
1365 | { |
1366 | snew(enesum[i], edat->nframes)(enesum[i]) = save_calloc("enesum[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1366, (edat->nframes), sizeof(*(enesum[i]))); |
1367 | } |
1368 | for (i = 0; (i < edat->nframes); i++) |
1369 | { |
1370 | eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]); |
1371 | eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]); |
1372 | eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]); |
1373 | for (j = 3; j <= 11; j++) |
1374 | { |
1375 | eneset[j][i] = edat->s[j].ener[i]; |
1376 | } |
1377 | eneset[11][i] -= Pres; |
1378 | enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum); |
1379 | enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum); |
1380 | enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum); |
1381 | } |
1382 | |
1383 | einstein_visco("evisco.xvg", "eviscoi.xvg", |
1384 | 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv); |
1385 | |
1386 | /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/ |
1387 | /* Do it for shear viscosity */ |
1388 | strcpy(buf, "Shear Viscosity"); |
1389 | low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3, |
1390 | (edat->nframes+1)/2, eneset, Dt, |
1391 | eacNormal(1<<0), 1, TRUE1, FALSE0, FALSE0, 0.0, 0.0, 0); |
1392 | |
1393 | /* Now for bulk viscosity */ |
1394 | strcpy(buf, "Bulk Viscosity"); |
1395 | low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1, |
1396 | (edat->nframes+1)/2, &(eneset[11]), Dt, |
1397 | eacNormal(1<<0), 1, TRUE1, FALSE0, FALSE0, 0.0, 0.0, 0); |
1398 | |
1399 | factor = (Vaver*1e-26/(BOLTZMANN(1.380658e-23)*Temp))*Dt; |
1400 | fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv); |
1401 | xvgr_legend(fp, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); |
1402 | |
1403 | /* Use trapezium rule for integration */ |
1404 | integral = 0; |
1405 | intBulk = 0; |
1406 | nout = get_acfnout(); |
1407 | if ((nout < 2) || (nout >= edat->nframes/2)) |
1408 | { |
1409 | nout = edat->nframes/2; |
1410 | } |
1411 | for (i = 1; (i < nout); i++) |
1412 | { |
1413 | integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor; |
1414 | intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor; |
1415 | fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk); |
1416 | } |
1417 | gmx_ffclose(fp); |
1418 | } |
1419 | else if (bCorr) |
1420 | { |
1421 | if (bFluct) |
1422 | { |
1423 | strcpy(buf, "Autocorrelation of Energy Fluctuations"); |
1424 | } |
1425 | else |
1426 | { |
1427 | strcpy(buf, "Energy Autocorrelation"); |
1428 | } |
1429 | #if 0 |
1430 | do_autocorr(corrfn, oenv, buf, edat->nframes, |
1431 | bSum ? 1 : nset, |
1432 | bSum ? &edat->s[nset-1].ener : eneset, |
1433 | (delta_t/edat->nframes), eacNormal(1<<0), FALSE0); |
1434 | #endif |
1435 | } |
1436 | } |
1437 | } |
1438 | |
1439 | static void print_time(FILE *fp, double t) |
1440 | { |
1441 | fprintf(fp, "%12.6f", t); |
1442 | } |
1443 | |
1444 | static void print1(FILE *fp, gmx_bool bDp, real e) |
1445 | { |
1446 | if (bDp) |
1447 | { |
1448 | fprintf(fp, " %16.12f", e); |
1449 | } |
1450 | else |
1451 | { |
1452 | fprintf(fp, " %10.6f", e); |
1453 | } |
1454 | } |
1455 | |
1456 | static void fec(const char *ene2fn, const char *runavgfn, |
1457 | real reftemp, int nset, int set[], char *leg[], |
1458 | enerdata_t *edat, double time[], |
1459 | const output_env_t oenv) |
1460 | { |
1461 | const char * ravgleg[] = { |
1462 | "\\8D\\4E = E\\sB\\N-E\\sA\\N", |
1463 | "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" |
1464 | }; |
1465 | FILE *fp; |
1466 | ener_file_t enx; |
1467 | int nre, timecheck, step, nenergy, nenergy2, maxenergy; |
1468 | int i, j; |
1469 | gmx_bool bCont; |
1470 | real aver, beta; |
1471 | real **eneset2; |
1472 | double dE, sum; |
1473 | gmx_enxnm_t *enm = NULL((void*)0); |
1474 | t_enxframe *fr; |
1475 | char buf[22]; |
1476 | |
1477 | /* read second energy file */ |
1478 | snew(fr, 1)(fr) = save_calloc("fr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1478, (1), sizeof(*(fr))); |
1479 | enm = NULL((void*)0); |
1480 | enx = open_enx(ene2fn, "r"); |
1481 | do_enxnms(enx, &(fr->nre), &enm); |
1482 | |
1483 | snew(eneset2, nset+1)(eneset2) = save_calloc("eneset2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1483, (nset+1), sizeof(*(eneset2))); |
1484 | nenergy2 = 0; |
1485 | maxenergy = 0; |
1486 | timecheck = 0; |
1487 | do |
1488 | { |
1489 | /* This loop searches for the first frame (when -b option is given), |
1490 | * or when this has been found it reads just one energy frame |
1491 | */ |
1492 | do |
1493 | { |
1494 | bCont = do_enx(enx, fr); |
1495 | |
1496 | if (bCont) |
1497 | { |
1498 | timecheck = check_times(fr->t); |
1499 | } |
1500 | |
1501 | } |
1502 | while (bCont && (timecheck < 0)); |
1503 | |
1504 | /* Store energies for analysis afterwards... */ |
1505 | if ((timecheck == 0) && bCont) |
1506 | { |
1507 | if (fr->nre > 0) |
1508 | { |
1509 | if (nenergy2 >= maxenergy) |
1510 | { |
1511 | maxenergy += 1000; |
1512 | for (i = 0; i <= nset; i++) |
1513 | { |
1514 | srenew(eneset2[i], maxenergy)(eneset2[i]) = save_realloc("eneset2[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1514, (eneset2[i]), (maxenergy), sizeof(*(eneset2[i]))); |
1515 | } |
1516 | } |
1517 | if (fr->t != time[nenergy2]) |
1518 | { |
1519 | fprintf(stderrstderr, "\nWARNING time mismatch %g!=%g at frame %s\n", |
1520 | fr->t, time[nenergy2], gmx_step_str(fr->step, buf)); |
1521 | } |
1522 | for (i = 0; i < nset; i++) |
1523 | { |
1524 | eneset2[i][nenergy2] = fr->ener[set[i]].e; |
1525 | } |
1526 | nenergy2++; |
1527 | } |
1528 | } |
1529 | } |
1530 | while (bCont && (timecheck == 0)); |
1531 | |
1532 | /* check */ |
1533 | if (edat->nframes != nenergy2) |
1534 | { |
1535 | fprintf(stderrstderr, "\nWARNING file length mismatch %d!=%d\n", |
1536 | edat->nframes, nenergy2); |
1537 | } |
1538 | nenergy = min(edat->nframes, nenergy2)(((edat->nframes) < (nenergy2)) ? (edat->nframes) : ( nenergy2) ); |
1539 | |
1540 | /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */ |
1541 | fp = NULL((void*)0); |
1542 | if (runavgfn) |
1543 | { |
1544 | fp = xvgropen(runavgfn, "Running average free energy difference", |
1545 | "Time (" unit_time"ps" ")", "\\8D\\4E (" unit_energy"kJ/mol" ")", oenv); |
1546 | xvgr_legend(fp, asize(ravgleg)((int)(sizeof(ravgleg)/sizeof((ravgleg)[0]))), ravgleg, oenv); |
1547 | } |
1548 | fprintf(stdoutstdout, "\n%-24s %10s\n", |
1549 | "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A"); |
1550 | sum = 0; |
1551 | beta = 1.0/(BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*reftemp); |
1552 | for (i = 0; i < nset; i++) |
1553 | { |
1554 | if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0) |
1555 | { |
1556 | fprintf(stderrstderr, "\nWARNING energy set name mismatch %s!=%s\n", |
1557 | leg[i], enm[set[i]].name); |
1558 | } |
1559 | for (j = 0; j < nenergy; j++) |
1560 | { |
1561 | dE = eneset2[i][j] - edat->s[i].ener[j]; |
1562 | sum += exp(-dE*beta); |
1563 | if (fp) |
1564 | { |
1565 | fprintf(fp, "%10g %10g %10g\n", |
1566 | time[j], dE, -BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*reftemp*log(sum/(j+1)) ); |
1567 | } |
1568 | } |
1569 | aver = -BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*reftemp*log(sum/nenergy); |
1570 | fprintf(stdoutstdout, "%-24s %10g\n", leg[i], aver); |
1571 | } |
1572 | if (fp) |
1573 | { |
1574 | gmx_ffclose(fp); |
1575 | } |
1576 | sfree(fr)save_free("fr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1576, (fr)); |
1577 | } |
1578 | |
1579 | |
1580 | static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, |
1581 | const char *filename, gmx_bool bDp, |
1582 | int *blocks, int *hists, int *samples, int *nlambdas, |
1583 | const output_env_t oenv) |
1584 | { |
1585 | const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda"; |
1586 | char title[STRLEN4096], label_x[STRLEN4096], label_y[STRLEN4096], legend[STRLEN4096]; |
1587 | char buf[STRLEN4096]; |
1588 | gmx_bool first = FALSE0; |
1589 | int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0; |
1590 | int i, j, k; |
1591 | /* coll data */ |
1592 | double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0; |
1593 | static int setnr = 0; |
1594 | double *native_lambda_vec = NULL((void*)0); |
1595 | const char **lambda_components = NULL((void*)0); |
1596 | int n_lambda_vec = 0; |
1597 | gmx_bool changing_lambda = FALSE0; |
1598 | int lambda_fep_state; |
1599 | |
1600 | /* now count the blocks & handle the global dh data */ |
1601 | for (i = 0; i < fr->nblock; i++) |
1602 | { |
1603 | if (fr->block[i].id == enxDHHIST) |
1604 | { |
1605 | nblock_hist++; |
1606 | } |
1607 | else if (fr->block[i].id == enxDH) |
1608 | { |
1609 | nblock_dh++; |
1610 | } |
1611 | else if (fr->block[i].id == enxDHCOLL) |
1612 | { |
1613 | nblock_dhcoll++; |
1614 | if ( (fr->block[i].nsub < 1) || |
1615 | (fr->block[i].sub[0].type != xdr_datatype_double) || |
1616 | (fr->block[i].sub[0].nr < 5)) |
1617 | { |
1618 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1618, "Unexpected block data"); |
1619 | } |
1620 | |
1621 | /* read the data from the DHCOLL block */ |
1622 | temp = fr->block[i].sub[0].dval[0]; |
1623 | start_time = fr->block[i].sub[0].dval[1]; |
1624 | delta_time = fr->block[i].sub[0].dval[2]; |
1625 | start_lambda = fr->block[i].sub[0].dval[3]; |
1626 | delta_lambda = fr->block[i].sub[0].dval[4]; |
1627 | changing_lambda = (delta_lambda != 0); |
1628 | if (fr->block[i].nsub > 1) |
1629 | { |
1630 | lambda_fep_state = fr->block[i].sub[1].ival[0]; |
1631 | if (n_lambda_vec == 0) |
1632 | { |
1633 | n_lambda_vec = fr->block[i].sub[1].ival[1]; |
1634 | } |
1635 | else |
1636 | { |
1637 | if (n_lambda_vec != fr->block[i].sub[1].ival[1]) |
1638 | { |
1639 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1639, |
1640 | "Unexpected change of basis set in lambda"); |
1641 | } |
1642 | } |
1643 | if (lambda_components == NULL((void*)0)) |
1644 | { |
1645 | snew(lambda_components, n_lambda_vec)(lambda_components) = save_calloc("lambda_components", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1645, (n_lambda_vec), sizeof(*(lambda_components))); |
1646 | } |
1647 | if (native_lambda_vec == NULL((void*)0)) |
1648 | { |
1649 | snew(native_lambda_vec, n_lambda_vec)(native_lambda_vec) = save_calloc("native_lambda_vec", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1649, (n_lambda_vec), sizeof(*(native_lambda_vec))); |
1650 | } |
1651 | for (j = 0; j < n_lambda_vec; j++) |
1652 | { |
1653 | native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j]; |
1654 | lambda_components[j] = |
1655 | efpt_singular_names[fr->block[i].sub[1].ival[2+j]]; |
1656 | } |
1657 | } |
1658 | } |
1659 | } |
1660 | |
1661 | if (nblock_hist == 0 && nblock_dh == 0) |
1662 | { |
1663 | /* don't do anything */ |
1664 | return; |
1665 | } |
1666 | if (nblock_hist > 0 && nblock_dh > 0) |
1667 | { |
1668 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1668, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do."); |
1669 | } |
1670 | if (!*fp_dhdl) |
1671 | { |
1672 | if (nblock_dh > 0) |
1673 | { |
1674 | /* we have standard, non-histogram data -- |
1675 | call open_dhdl to open the file */ |
1676 | /* TODO this is an ugly hack that needs to be fixed: this will only |
1677 | work if the order of data is always the same and if we're |
1678 | only using the g_energy compiled with the mdrun that produced |
1679 | the ener.edr. */ |
1680 | *fp_dhdl = open_dhdl(filename, ir, oenv); |
1681 | } |
1682 | else |
1683 | { |
1684 | sprintf(title, "N(%s)", deltag); |
1685 | sprintf(label_x, "%s (%s)", deltag, unit_energy"kJ/mol"); |
1686 | sprintf(label_y, "Samples"); |
1687 | *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv); |
1688 | sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda); |
1689 | xvgr_subtitle(*fp_dhdl, buf, oenv); |
1690 | } |
1691 | } |
1692 | |
1693 | (*hists) += nblock_hist; |
1694 | (*blocks) += nblock_dh; |
1695 | (*nlambdas) = nblock_hist+nblock_dh; |
1696 | |
1697 | /* write the data */ |
1698 | if (nblock_hist > 0) |
1699 | { |
1700 | gmx_int64_t sum = 0; |
1701 | /* histograms */ |
1702 | for (i = 0; i < fr->nblock; i++) |
1703 | { |
1704 | t_enxblock *blk = &(fr->block[i]); |
1705 | if (blk->id == enxDHHIST) |
1706 | { |
1707 | double foreign_lambda, dx; |
1708 | gmx_int64_t x0; |
1709 | int nhist, derivative; |
1710 | |
1711 | /* check the block types etc. */ |
1712 | if ( (blk->nsub < 2) || |
1713 | (blk->sub[0].type != xdr_datatype_double) || |
1714 | (blk->sub[1].type != xdr_datatype_int64) || |
1715 | (blk->sub[0].nr < 2) || |
1716 | (blk->sub[1].nr < 2) ) |
1717 | { |
1718 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1718, "Unexpected block data in file"); |
1719 | } |
1720 | foreign_lambda = blk->sub[0].dval[0]; |
1721 | dx = blk->sub[0].dval[1]; |
1722 | nhist = blk->sub[1].lval[0]; |
1723 | derivative = blk->sub[1].lval[1]; |
1724 | for (j = 0; j < nhist; j++) |
1725 | { |
1726 | const char *lg[1]; |
1727 | x0 = blk->sub[1].lval[2+j]; |
1728 | |
1729 | if (!derivative) |
1730 | { |
1731 | sprintf(legend, "N(%s(%s=%g) | %s=%g)", |
1732 | deltag, lambda, foreign_lambda, |
1733 | lambda, start_lambda); |
1734 | } |
1735 | else |
1736 | { |
1737 | sprintf(legend, "N(%s | %s=%g)", |
1738 | dhdl, lambda, start_lambda); |
1739 | } |
1740 | |
1741 | lg[0] = legend; |
1742 | xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv); |
1743 | setnr++; |
1744 | for (k = 0; k < blk->sub[j+2].nr; k++) |
1745 | { |
1746 | int hist; |
1747 | double xmin, xmax; |
1748 | |
1749 | hist = blk->sub[j+2].ival[k]; |
1750 | xmin = (x0+k)*dx; |
1751 | xmax = (x0+k+1)*dx; |
1752 | fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist, |
1753 | xmax, hist); |
1754 | sum += hist; |
1755 | } |
1756 | /* multiple histogram data blocks in one histogram |
1757 | mean that the second one is the reverse of the first one: |
1758 | for dhdl derivatives, it's important to know both the |
1759 | maximum and minimum values */ |
1760 | dx = -dx; |
1761 | } |
1762 | } |
1763 | } |
1764 | (*samples) += (int)(sum/nblock_hist); |
1765 | } |
1766 | else |
1767 | { |
1768 | /* raw dh */ |
1769 | int len = 0; |
1770 | char **setnames = NULL((void*)0); |
1771 | int nnames = nblock_dh; |
Value stored to 'nnames' during its initialization is never read | |
1772 | |
1773 | for (i = 0; i < fr->nblock; i++) |
1774 | { |
1775 | t_enxblock *blk = &(fr->block[i]); |
1776 | if (blk->id == enxDH) |
1777 | { |
1778 | if (len == 0) |
1779 | { |
1780 | len = blk->sub[2].nr; |
1781 | } |
1782 | else |
1783 | { |
1784 | if (len != blk->sub[2].nr) |
1785 | { |
1786 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 1786, "Length inconsistency in dhdl data"); |
1787 | } |
1788 | } |
1789 | } |
1790 | } |
1791 | (*samples) += len; |
1792 | |
1793 | for (i = 0; i < len; i++) |
1794 | { |
1795 | double time = start_time + delta_time*i; |
1796 | |
1797 | fprintf(*fp_dhdl, "%.4f ", time); |
1798 | |
1799 | for (j = 0; j < fr->nblock; j++) |
1800 | { |
1801 | t_enxblock *blk = &(fr->block[j]); |
1802 | if (blk->id == enxDH) |
1803 | { |
1804 | double value; |
1805 | if (blk->sub[2].type == xdr_datatype_float) |
1806 | { |
1807 | value = blk->sub[2].fval[i]; |
1808 | } |
1809 | else |
1810 | { |
1811 | value = blk->sub[2].dval[i]; |
1812 | } |
1813 | /* we need to decide which data type it is based on the count*/ |
1814 | |
1815 | if (j == 1 && ir->bExpanded) |
1816 | { |
1817 | fprintf(*fp_dhdl, "%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */ |
1818 | } |
1819 | else |
1820 | { |
1821 | if (bDp) |
1822 | { |
1823 | fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */ |
1824 | } |
1825 | else |
1826 | { |
1827 | fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */ |
1828 | } |
1829 | } |
1830 | } |
1831 | } |
1832 | fprintf(*fp_dhdl, "\n"); |
1833 | } |
1834 | } |
1835 | } |
1836 | |
1837 | |
1838 | int gmx_energy(int argc, char *argv[]) |
1839 | { |
1840 | const char *desc[] = { |
1841 | "[THISMODULE] extracts energy components or distance restraint", |
1842 | "data from an energy file. The user is prompted to interactively", |
1843 | "select the desired energy terms.[PAR]", |
1844 | |
1845 | "Average, RMSD, and drift are calculated with full precision from the", |
1846 | "simulation (see printed manual). Drift is calculated by performing", |
1847 | "a least-squares fit of the data to a straight line. The reported total drift", |
1848 | "is the difference of the fit at the first and last point.", |
1849 | "An error estimate of the average is given based on a block averages", |
1850 | "over 5 blocks using the full-precision averages. The error estimate", |
1851 | "can be performed over multiple block lengths with the options", |
1852 | "[TT]-nbmin[tt] and [TT]-nbmax[tt].", |
1853 | "[BB]Note[bb] that in most cases the energy files contains averages over all", |
1854 | "MD steps, or over many more points than the number of frames in", |
1855 | "energy file. This makes the [THISMODULE] statistics output more accurate", |
1856 | "than the [TT].xvg[tt] output. When exact averages are not present in the energy", |
1857 | "file, the statistics mentioned above are simply over the single, per-frame", |
1858 | "energy values.[PAR]", |
1859 | |
1860 | "The term fluctuation gives the RMSD around the least-squares fit.[PAR]", |
1861 | |
1862 | "Some fluctuation-dependent properties can be calculated provided", |
1863 | "the correct energy terms are selected, and that the command line option", |
1864 | "[TT]-fluct_props[tt] is given. The following properties", |
1865 | "will be computed:[BR]", |
1866 | "Property Energy terms needed[BR]", |
1867 | "---------------------------------------------------[BR]", |
1868 | "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]", |
1869 | "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]", |
1870 | "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]", |
1871 | "Isothermal compressibility: Vol, Temp [BR]", |
1872 | "Adiabatic bulk modulus: Vol, Temp [BR]", |
1873 | "---------------------------------------------------[BR]", |
1874 | "You always need to set the number of molecules [TT]-nmol[tt].", |
1875 | "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections", |
1876 | "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]" |
1877 | "When the [TT]-viol[tt] option is set, the time averaged", |
1878 | "violations are plotted and the running time-averaged and", |
1879 | "instantaneous sum of violations are recalculated. Additionally", |
1880 | "running time-averaged and instantaneous distances between", |
1881 | "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]", |
1882 | |
1883 | "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and", |
1884 | "[TT]-odt[tt] are used for analyzing orientation restraint data.", |
1885 | "The first two options plot the orientation, the last three the", |
1886 | "deviations of the orientations from the experimental values.", |
1887 | "The options that end on an 'a' plot the average over time", |
1888 | "as a function of restraint. The options that end on a 't'", |
1889 | "prompt the user for restraint label numbers and plot the data", |
1890 | "as a function of time. Option [TT]-odr[tt] plots the RMS", |
1891 | "deviation as a function of restraint.", |
1892 | "When the run used time or ensemble averaged orientation restraints,", |
1893 | "option [TT]-orinst[tt] can be used to analyse the instantaneous,", |
1894 | "not ensemble-averaged orientations and deviations instead of", |
1895 | "the time and ensemble averages.[PAR]", |
1896 | |
1897 | "Option [TT]-oten[tt] plots the eigenvalues of the molecular order", |
1898 | "tensor for each orientation restraint experiment. With option", |
1899 | "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]", |
1900 | |
1901 | "Option [TT]-odh[tt] extracts and plots the free energy data", |
1902 | "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)", |
1903 | "from the [TT]ener.edr[tt] file.[PAR]", |
1904 | |
1905 | "With [TT]-fee[tt] an estimate is calculated for the free-energy", |
1906 | "difference with an ideal gas state: [BR]", |
1907 | " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]", |
1908 | " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]", |
1909 | "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and", |
1910 | "the average is over the ensemble (or time in a trajectory).", |
1911 | "Note that this is in principle", |
1912 | "only correct when averaging over the whole (Boltzmann) ensemble", |
1913 | "and using the potential energy. This also allows for an entropy", |
1914 | "estimate using:[BR]", |
1915 | " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]", |
1916 | " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T", |
1917 | "[PAR]", |
1918 | |
1919 | "When a second energy file is specified ([TT]-f2[tt]), a free energy", |
1920 | "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,", |
1921 | "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy", |
1922 | "files, and the average is over the ensemble A. The running average", |
1923 | "of the free energy difference is printed to a file specified by [TT]-ravg[tt].", |
1924 | "[BB]Note[bb] that the energies must both be calculated from the same trajectory." |
1925 | |
1926 | }; |
1927 | static gmx_bool bSum = FALSE0, bFee = FALSE0, bPrAll = FALSE0, bFluct = FALSE0, bDriftCorr = FALSE0; |
1928 | static gmx_bool bDp = FALSE0, bMutot = FALSE0, bOrinst = FALSE0, bOvec = FALSE0, bFluctProps = FALSE0; |
1929 | static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5; |
1930 | static real reftemp = 300.0, ezero = 0; |
1931 | t_pargs pa[] = { |
1932 | { "-fee", FALSE0, etBOOL, {&bFee}, |
1933 | "Do a free energy estimate" }, |
1934 | { "-fetemp", FALSE0, etREAL, {&reftemp}, |
1935 | "Reference temperature for free energy calculation" }, |
1936 | { "-zero", FALSE0, etREAL, {&ezero}, |
1937 | "Subtract a zero-point energy" }, |
1938 | { "-sum", FALSE0, etBOOL, {&bSum}, |
1939 | "Sum the energy terms selected rather than display them all" }, |
1940 | { "-dp", FALSE0, etBOOL, {&bDp}, |
1941 | "Print energies in high precision" }, |
1942 | { "-nbmin", FALSE0, etINT, {&nbmin}, |
1943 | "Minimum number of blocks for error estimate" }, |
1944 | { "-nbmax", FALSE0, etINT, {&nbmax}, |
1945 | "Maximum number of blocks for error estimate" }, |
1946 | { "-mutot", FALSE0, etBOOL, {&bMutot}, |
1947 | "Compute the total dipole moment from the components" }, |
1948 | { "-skip", FALSE0, etINT, {&skip}, |
1949 | "Skip number of frames between data points" }, |
1950 | { "-aver", FALSE0, etBOOL, {&bPrAll}, |
1951 | "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" }, |
1952 | { "-nmol", FALSE0, etINT, {&nmol}, |
1953 | "Number of molecules in your sample: the energies are divided by this number" }, |
1954 | { "-fluct_props", FALSE0, etBOOL, {&bFluctProps}, |
1955 | "Compute properties based on energy fluctuations, like heat capacity" }, |
1956 | { "-driftcorr", FALSE0, etBOOL, {&bDriftCorr}, |
1957 | "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."}, |
1958 | { "-fluc", FALSE0, etBOOL, {&bFluct}, |
1959 | "Calculate autocorrelation of energy fluctuations rather than energy itself" }, |
1960 | { "-orinst", FALSE0, etBOOL, {&bOrinst}, |
1961 | "Analyse instantaneous orientation data" }, |
1962 | { "-ovec", FALSE0, etBOOL, {&bOvec}, |
1963 | "Also plot the eigenvectors with [TT]-oten[tt]" } |
1964 | }; |
1965 | const char * drleg[] = { |
1966 | "Running average", |
1967 | "Instantaneous" |
1968 | }; |
1969 | static const char *setnm[] = { |
1970 | "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY", |
1971 | "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature", |
1972 | "Volume", "Pressure" |
1973 | }; |
1974 | |
1975 | FILE *out = NULL((void*)0), *fp_pairs = NULL((void*)0), *fort = NULL((void*)0), *fodt = NULL((void*)0), *foten = NULL((void*)0); |
1976 | FILE *fp_dhdl = NULL((void*)0); |
1977 | FILE **drout; |
1978 | ener_file_t fp; |
1979 | int timecheck = 0; |
1980 | gmx_mtop_t mtop; |
1981 | gmx_localtop_t *top = NULL((void*)0); |
1982 | t_inputrec ir; |
1983 | t_energy **ee; |
1984 | enerdata_t edat; |
1985 | gmx_enxnm_t *enm = NULL((void*)0); |
1986 | t_enxframe *frame, *fr = NULL((void*)0); |
1987 | int cur = 0; |
1988 | #define NEXT(1-cur) (1-cur) |
1989 | int nre, teller, teller_disre, nfr; |
1990 | gmx_int64_t start_step; |
1991 | int nor = 0, nex = 0, norfr = 0, enx_i = 0; |
1992 | real start_t; |
1993 | real *bounds = NULL((void*)0), *violaver = NULL((void*)0), *oobs = NULL((void*)0), *orient = NULL((void*)0), *odrms = NULL((void*)0); |
1994 | int *index = NULL((void*)0), *pair = NULL((void*)0), norsel = 0, *orsel = NULL((void*)0), *or_label = NULL((void*)0); |
1995 | int nbounds = 0, npairs; |
1996 | gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL; |
1997 | gmx_bool bFoundStart, bCont, bEDR, bVisco; |
1998 | double sum, sumaver, sumt, ener, dbl; |
1999 | double *time = NULL((void*)0); |
2000 | real Vaver; |
2001 | int *set = NULL((void*)0), i, j, k, nset, sss; |
2002 | gmx_bool *bIsEner = NULL((void*)0); |
2003 | char **pairleg, **odtleg, **otenleg; |
2004 | char **leg = NULL((void*)0); |
2005 | char **nms; |
2006 | char *anm_j, *anm_k, *resnm_j, *resnm_k; |
2007 | int resnr_j, resnr_k; |
2008 | const char *orinst_sub = "@ subtitle \"instantaneous\"\n"; |
2009 | char buf[256]; |
2010 | output_env_t oenv; |
2011 | t_enxblock *blk = NULL((void*)0); |
2012 | t_enxblock *blk_disre = NULL((void*)0); |
2013 | int ndisre = 0; |
2014 | int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0; |
2015 | |
2016 | t_filenm fnm[] = { |
2017 | { efEDR, "-f", NULL((void*)0), ffREAD1<<1 }, |
2018 | { efEDR, "-f2", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
2019 | { efTPX, "-s", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
2020 | { efXVG, "-o", "energy", ffWRITE1<<2 }, |
2021 | { efXVG, "-viol", "violaver", ffOPTWR(1<<2| 1<<3) }, |
2022 | { efXVG, "-pairs", "pairs", ffOPTWR(1<<2| 1<<3) }, |
2023 | { efXVG, "-ora", "orienta", ffOPTWR(1<<2| 1<<3) }, |
2024 | { efXVG, "-ort", "orientt", ffOPTWR(1<<2| 1<<3) }, |
2025 | { efXVG, "-oda", "orideva", ffOPTWR(1<<2| 1<<3) }, |
2026 | { efXVG, "-odr", "oridevr", ffOPTWR(1<<2| 1<<3) }, |
2027 | { efXVG, "-odt", "oridevt", ffOPTWR(1<<2| 1<<3) }, |
2028 | { efXVG, "-oten", "oriten", ffOPTWR(1<<2| 1<<3) }, |
2029 | { efXVG, "-corr", "enecorr", ffOPTWR(1<<2| 1<<3) }, |
2030 | { efXVG, "-vis", "visco", ffOPTWR(1<<2| 1<<3) }, |
2031 | { efXVG, "-ravg", "runavgdf", ffOPTWR(1<<2| 1<<3) }, |
2032 | { efXVG, "-odh", "dhdl", ffOPTWR(1<<2| 1<<3) } |
2033 | }; |
2034 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
2035 | int npargs; |
2036 | t_pargs *ppa; |
2037 | |
2038 | npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))); |
2039 | ppa = add_acf_pargs(&npargs, pa); |
2040 | if (!parse_common_args(&argc, argv, |
2041 | PCA_CAN_VIEW(1<<5) | PCA_CAN_BEGIN(1<<6) | PCA_CAN_END(1<<7) | PCA_BE_NICE(1<<13), |
2042 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, ppa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) |
2043 | { |
2044 | return 0; |
2045 | } |
2046 | |
2047 | bDRAll = opt2bSet("-pairs", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2048 | bDisRe = opt2bSet("-viol", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || bDRAll; |
2049 | bORA = opt2bSet("-ora", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2050 | bORT = opt2bSet("-ort", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2051 | bODA = opt2bSet("-oda", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2052 | bODR = opt2bSet("-odr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2053 | bODT = opt2bSet("-odt", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2054 | bORIRE = bORA || bORT || bODA || bODR || bODT; |
2055 | bOTEN = opt2bSet("-oten", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2056 | bDHDL = opt2bSet("-odh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2057 | |
2058 | nset = 0; |
2059 | |
2060 | snew(frame, 2)(frame) = save_calloc("frame", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2060, (2), sizeof(*(frame))); |
2061 | fp = open_enx(ftp2fn(efEDR, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "r"); |
2062 | do_enxnms(fp, &nre, &enm); |
2063 | |
2064 | Vaver = -1; |
2065 | |
2066 | bVisco = opt2bSet("-vis", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
2067 | |
2068 | if ((!bDisRe) && (!bDHDL)) |
2069 | { |
2070 | if (bVisco) |
2071 | { |
2072 | nset = asize(setnm)((int)(sizeof(setnm)/sizeof((setnm)[0]))); |
2073 | snew(set, nset)(set) = save_calloc("set", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2073, (nset), sizeof(*(set))); |
2074 | /* This is nasty code... To extract Pres tensor, Volume and Temperature */ |
2075 | for (j = 0; j < nset; j++) |
2076 | { |
2077 | for (i = 0; i < nre; i++) |
2078 | { |
2079 | if (strstr(enm[i].name, setnm[j])) |
2080 | { |
2081 | set[j] = i; |
2082 | break; |
2083 | } |
2084 | } |
2085 | if (i == nre) |
2086 | { |
2087 | if (gmx_strcasecmp(setnm[j], "Volume") == 0) |
2088 | { |
2089 | printf("Enter the box volume (" unit_volume"nm" "^3" "): "); |
2090 | if (1 != scanf("%lf", &dbl)) |
2091 | { |
2092 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2092, "Error reading user input"); |
2093 | } |
2094 | Vaver = dbl; |
2095 | } |
2096 | else |
2097 | { |
2098 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2098, "Could not find term %s for viscosity calculation", |
2099 | setnm[j]); |
2100 | } |
2101 | } |
2102 | } |
2103 | } |
2104 | else |
2105 | { |
2106 | set = select_by_name(nre, enm, &nset); |
2107 | } |
2108 | /* Print all the different units once */ |
2109 | sprintf(buf, "(%s)", enm[set[0]].unit); |
2110 | for (i = 1; i < nset; i++) |
2111 | { |
2112 | for (j = 0; j < i; j++) |
2113 | { |
2114 | if (strcmp(enm[set[i]].unit, enm[set[j]].unit)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (enm[set[i]].unit) && __builtin_constant_p (enm[set[ j]].unit) && (__s1_len = strlen (enm[set[i]].unit), __s2_len = strlen (enm[set[j]].unit), (!((size_t)(const void *)((enm[ set[i]].unit) + 1) - (size_t)(const void *)(enm[set[i]].unit) == 1) || __s1_len >= 4) && (!((size_t)(const void *)((enm[set[j]].unit) + 1) - (size_t)(const void *)(enm[set[ j]].unit) == 1) || __s2_len >= 4)) ? __builtin_strcmp (enm [set[i]].unit, enm[set[j]].unit) : (__builtin_constant_p (enm [set[i]].unit) && ((size_t)(const void *)((enm[set[i] ].unit) + 1) - (size_t)(const void *)(enm[set[i]].unit) == 1) && (__s1_len = strlen (enm[set[i]].unit), __s1_len < 4) ? (__builtin_constant_p (enm[set[j]].unit) && ((size_t )(const void *)((enm[set[j]].unit) + 1) - (size_t)(const void *)(enm[set[j]].unit) == 1) ? __builtin_strcmp (enm[set[i]].unit , enm[set[j]].unit) : (__extension__ ({ const unsigned char * __s2 = (const unsigned char *) (const char *) (enm[set[j]].unit ); int __result = (((const unsigned char *) (const char *) (enm [set[i]].unit))[0] - __s2[0]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (enm[set[i]].unit))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (enm[set[i]].unit))[2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (enm[set[i]].unit))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (enm[set[j]].unit) && ((size_t)(const void *)((enm[set[j]].unit) + 1) - (size_t)(const void *)(enm[set[j]].unit) == 1) && (__s2_len = strlen (enm[set[j]].unit), __s2_len < 4) ? (__builtin_constant_p (enm[set[i]].unit) && ((size_t)(const void *)((enm[set [i]].unit) + 1) - (size_t)(const void *)(enm[set[i]].unit) == 1) ? __builtin_strcmp (enm[set[i]].unit, enm[set[j]].unit) : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (enm[set[i]].unit); int __result = (( (const unsigned char *) (const char *) (enm[set[j]].unit))[0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (enm[set[j]].unit ))[1] - __s2[1]); if (__s2_len > 1 && __result == 0 ) { __result = (((const unsigned char *) (const char *) (enm[ set[j]].unit))[2] - __s2[2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (enm [set[j]].unit))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (enm[set[i]].unit, enm[set[j]].unit)))); }) == 0) |
2115 | { |
2116 | break; |
2117 | } |
2118 | } |
2119 | if (j == i) |
2120 | { |
2121 | strcat(buf, ", ("); |
2122 | strcat(buf, enm[set[i]].unit); |
2123 | strcat(buf, ")"); |
2124 | } |
2125 | } |
2126 | out = xvgropen(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Gromacs Energies", "Time (ps)", buf, |
2127 | oenv); |
2128 | |
2129 | snew(leg, nset+1)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2129, (nset+1), sizeof(*(leg))); |
2130 | for (i = 0; (i < nset); i++) |
2131 | { |
2132 | leg[i] = enm[set[i]].name; |
2133 | } |
2134 | if (bSum) |
2135 | { |
2136 | leg[nset] = strdup("Sum")(__extension__ (__builtin_constant_p ("Sum") && ((size_t )(const void *)(("Sum") + 1) - (size_t)(const void *)("Sum") == 1) ? (((const char *) ("Sum"))[0] == '\0' ? (char *) calloc ( (size_t) 1, (size_t) 1) : ({ size_t __len = strlen ("Sum") + 1 ; char *__retval = (char *) malloc (__len); if (__retval != ( (void*)0)) __retval = (char *) memcpy (__retval, "Sum", __len ); __retval; })) : __strdup ("Sum"))); |
2137 | xvgr_legend(out, nset+1, (const char**)leg, oenv); |
2138 | } |
2139 | else |
2140 | { |
2141 | xvgr_legend(out, nset, (const char**)leg, oenv); |
2142 | } |
2143 | |
2144 | snew(bIsEner, nset)(bIsEner) = save_calloc("bIsEner", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2144, (nset), sizeof(*(bIsEner))); |
2145 | for (i = 0; (i < nset); i++) |
2146 | { |
2147 | bIsEner[i] = FALSE0; |
2148 | for (j = 0; (j <= F_ETOT); j++) |
2149 | { |
2150 | bIsEner[i] = bIsEner[i] || |
2151 | (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0); |
2152 | } |
2153 | } |
2154 | |
2155 | if (bPrAll && nset > 1) |
2156 | { |
2157 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2157, "Printing averages can only be done when a single set is selected"); |
2158 | } |
2159 | |
2160 | time = NULL((void*)0); |
2161 | |
2162 | if (bORIRE || bOTEN) |
2163 | { |
2164 | get_orires_parms(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &nor, &nex, &or_label, &oobs); |
2165 | } |
2166 | |
2167 | if (bORIRE) |
2168 | { |
2169 | if (bOrinst) |
2170 | { |
2171 | enx_i = enxORI; |
2172 | } |
2173 | else |
2174 | { |
2175 | enx_i = enxOR; |
2176 | } |
2177 | |
2178 | if (bORA || bODA) |
2179 | { |
2180 | snew(orient, nor)(orient) = save_calloc("orient", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2180, (nor), sizeof(*(orient))); |
2181 | } |
2182 | if (bODR) |
2183 | { |
2184 | snew(odrms, nor)(odrms) = save_calloc("odrms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2184, (nor), sizeof(*(odrms))); |
2185 | } |
2186 | if (bORT || bODT) |
2187 | { |
2188 | fprintf(stderrstderr, "Select the orientation restraint labels you want (-1 is all)\n"); |
2189 | fprintf(stderrstderr, "End your selection with 0\n"); |
2190 | j = -1; |
2191 | orsel = NULL((void*)0); |
2192 | do |
2193 | { |
2194 | j++; |
2195 | srenew(orsel, j+1)(orsel) = save_realloc("orsel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2195, (orsel), (j+1), sizeof(*(orsel))); |
2196 | if (1 != scanf("%d", &(orsel[j]))) |
2197 | { |
2198 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2198, "Error reading user input"); |
2199 | } |
2200 | } |
2201 | while (orsel[j] > 0); |
2202 | if (orsel[0] == -1) |
2203 | { |
2204 | fprintf(stderrstderr, "Selecting all %d orientation restraints\n", nor); |
2205 | norsel = nor; |
2206 | srenew(orsel, nor)(orsel) = save_realloc("orsel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2206, (orsel), (nor), sizeof(*(orsel))); |
2207 | for (i = 0; i < nor; i++) |
2208 | { |
2209 | orsel[i] = i; |
2210 | } |
2211 | } |
2212 | else |
2213 | { |
2214 | /* Build the selection */ |
2215 | norsel = 0; |
2216 | for (i = 0; i < j; i++) |
2217 | { |
2218 | for (k = 0; k < nor; k++) |
2219 | { |
2220 | if (or_label[k] == orsel[i]) |
2221 | { |
2222 | orsel[norsel] = k; |
2223 | norsel++; |
2224 | break; |
2225 | } |
2226 | } |
2227 | if (k == nor) |
2228 | { |
2229 | fprintf(stderrstderr, "Orientation restraint label %d not found\n", |
2230 | orsel[i]); |
2231 | } |
2232 | } |
2233 | } |
2234 | snew(odtleg, norsel)(odtleg) = save_calloc("odtleg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2234, (norsel), sizeof(*(odtleg))); |
2235 | for (i = 0; i < norsel; i++) |
2236 | { |
2237 | snew(odtleg[i], 256)(odtleg[i]) = save_calloc("odtleg[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2237, (256), sizeof(*(odtleg[i]))); |
2238 | sprintf(odtleg[i], "%d", or_label[orsel[i]]); |
2239 | } |
2240 | if (bORT) |
2241 | { |
2242 | fort = xvgropen(opt2fn("-ort", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Calculated orientations", |
2243 | "Time (ps)", "", oenv); |
2244 | if (bOrinst) |
2245 | { |
2246 | fprintf(fort, "%s", orinst_sub); |
2247 | } |
2248 | xvgr_legend(fort, norsel, (const char**)odtleg, oenv); |
2249 | } |
2250 | if (bODT) |
2251 | { |
2252 | fodt = xvgropen(opt2fn("-odt", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2253 | "Orientation restraint deviation", |
2254 | "Time (ps)", "", oenv); |
2255 | if (bOrinst) |
2256 | { |
2257 | fprintf(fodt, "%s", orinst_sub); |
2258 | } |
2259 | xvgr_legend(fodt, norsel, (const char**)odtleg, oenv); |
2260 | } |
2261 | } |
2262 | } |
2263 | if (bOTEN) |
2264 | { |
2265 | foten = xvgropen(opt2fn("-oten", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2266 | "Order tensor", "Time (ps)", "", oenv); |
2267 | snew(otenleg, bOvec ? nex*12 : nex*3)(otenleg) = save_calloc("otenleg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2267, (bOvec ? nex*12 : nex*3), sizeof(*(otenleg))); |
2268 | for (i = 0; i < nex; i++) |
2269 | { |
2270 | for (j = 0; j < 3; j++) |
2271 | { |
2272 | sprintf(buf, "eig%d", j+1); |
2273 | otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); |
2274 | } |
2275 | if (bOvec) |
2276 | { |
2277 | for (j = 0; j < 9; j++) |
2278 | { |
2279 | sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z")); |
2280 | otenleg[12*i+3+j] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); |
2281 | } |
2282 | } |
2283 | } |
2284 | xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv); |
2285 | } |
2286 | } |
2287 | else if (bDisRe) |
2288 | { |
2289 | nbounds = get_bounds(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &bounds, &index, &pair, &npairs, |
2290 | &mtop, &top, &ir); |
2291 | snew(violaver, npairs)(violaver) = save_calloc("violaver", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2291, (npairs), sizeof(*(violaver))); |
2292 | out = xvgropen(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Sum of Violations", |
2293 | "Time (ps)", "nm", oenv); |
2294 | xvgr_legend(out, 2, drleg, oenv); |
2295 | if (bDRAll) |
2296 | { |
2297 | fp_pairs = xvgropen(opt2fn("-pairs", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Pair Distances", |
2298 | "Time (ps)", "Distance (nm)", oenv); |
2299 | if (output_env_get_print_xvgr_codes(oenv)) |
2300 | { |
2301 | fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n", |
2302 | ir.dr_tau); |
2303 | } |
2304 | } |
2305 | } |
2306 | else if (bDHDL) |
2307 | { |
2308 | get_dhdl_parms(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &ir); |
2309 | } |
2310 | |
2311 | /* Initiate energies and set them to zero */ |
2312 | edat.nsteps = 0; |
2313 | edat.npoints = 0; |
2314 | edat.nframes = 0; |
2315 | edat.step = NULL((void*)0); |
2316 | edat.steps = NULL((void*)0); |
2317 | edat.points = NULL((void*)0); |
2318 | snew(edat.s, nset)(edat.s) = save_calloc("edat.s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2318, (nset), sizeof(*(edat.s))); |
2319 | |
2320 | /* Initiate counters */ |
2321 | teller = 0; |
2322 | teller_disre = 0; |
2323 | bFoundStart = FALSE0; |
2324 | start_step = 0; |
2325 | start_t = 0; |
2326 | do |
2327 | { |
2328 | /* This loop searches for the first frame (when -b option is given), |
2329 | * or when this has been found it reads just one energy frame |
2330 | */ |
2331 | do |
2332 | { |
2333 | bCont = do_enx(fp, &(frame[NEXT(1-cur)])); |
2334 | if (bCont) |
2335 | { |
2336 | timecheck = check_times(frame[NEXT(1-cur)].t); |
2337 | } |
2338 | } |
2339 | while (bCont && (timecheck < 0)); |
2340 | |
2341 | if ((timecheck == 0) && bCont) |
2342 | { |
2343 | /* We read a valid frame, so we can use it */ |
2344 | fr = &(frame[NEXT(1-cur)]); |
2345 | |
2346 | if (fr->nre > 0) |
2347 | { |
2348 | /* The frame contains energies, so update cur */ |
2349 | cur = NEXT(1-cur); |
2350 | |
2351 | if (edat.nframes % 1000 == 0) |
2352 | { |
2353 | srenew(edat.step, edat.nframes+1000)(edat.step) = save_realloc("edat.step", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2353, (edat.step), (edat.nframes+1000), sizeof(*(edat.step) )); |
2354 | memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0])); |
2355 | srenew(edat.steps, edat.nframes+1000)(edat.steps) = save_realloc("edat.steps", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2355, (edat.steps), (edat.nframes+1000), sizeof(*(edat.steps ))); |
2356 | memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0])); |
2357 | srenew(edat.points, edat.nframes+1000)(edat.points) = save_realloc("edat.points", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2357, (edat.points), (edat.nframes+1000), sizeof(*(edat.points ))); |
2358 | memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0])); |
2359 | |
2360 | for (i = 0; i < nset; i++) |
2361 | { |
2362 | srenew(edat.s[i].ener, edat.nframes+1000)(edat.s[i].ener) = save_realloc("edat.s[i].ener", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2362, (edat.s[i].ener), (edat.nframes+1000), sizeof(*(edat. s[i].ener))); |
2363 | memset(&(edat.s[i].ener[edat.nframes]), 0, |
2364 | 1000*sizeof(edat.s[i].ener[0])); |
2365 | srenew(edat.s[i].es, edat.nframes+1000)(edat.s[i].es) = save_realloc("edat.s[i].es", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2365, (edat.s[i].es), (edat.nframes+1000), sizeof(*(edat.s[ i].es))); |
2366 | memset(&(edat.s[i].es[edat.nframes]), 0, |
2367 | 1000*sizeof(edat.s[i].es[0])); |
2368 | } |
2369 | } |
2370 | |
2371 | nfr = edat.nframes; |
2372 | edat.step[nfr] = fr->step; |
2373 | |
2374 | if (!bFoundStart) |
2375 | { |
2376 | bFoundStart = TRUE1; |
2377 | /* Initiate the previous step data */ |
2378 | start_step = fr->step; |
2379 | start_t = fr->t; |
2380 | /* Initiate the energy sums */ |
2381 | edat.steps[nfr] = 1; |
2382 | edat.points[nfr] = 1; |
2383 | for (i = 0; i < nset; i++) |
2384 | { |
2385 | sss = set[i]; |
2386 | edat.s[i].es[nfr].sum = fr->ener[sss].e; |
2387 | edat.s[i].es[nfr].sum2 = 0; |
2388 | } |
2389 | edat.nsteps = 1; |
2390 | edat.npoints = 1; |
2391 | } |
2392 | else |
2393 | { |
2394 | edat.steps[nfr] = fr->nsteps; |
2395 | { |
2396 | if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps) |
2397 | { |
2398 | if (fr->nsum <= 1) |
2399 | { |
2400 | edat.points[nfr] = 1; |
2401 | for (i = 0; i < nset; i++) |
2402 | { |
2403 | sss = set[i]; |
2404 | edat.s[i].es[nfr].sum = fr->ener[sss].e; |
2405 | edat.s[i].es[nfr].sum2 = 0; |
2406 | } |
2407 | edat.npoints += 1; |
2408 | } |
2409 | else |
2410 | { |
2411 | edat.points[nfr] = fr->nsum; |
2412 | for (i = 0; i < nset; i++) |
2413 | { |
2414 | sss = set[i]; |
2415 | edat.s[i].es[nfr].sum = fr->ener[sss].esum; |
2416 | edat.s[i].es[nfr].sum2 = fr->ener[sss].eav; |
2417 | } |
2418 | edat.npoints += fr->nsum; |
2419 | } |
2420 | } |
2421 | else |
2422 | { |
2423 | /* The interval does not match fr->nsteps: |
2424 | * can not do exact averages. |
2425 | */ |
2426 | edat.npoints = 0; |
2427 | } |
2428 | edat.nsteps = fr->step - start_step + 1; |
2429 | } |
2430 | } |
2431 | for (i = 0; i < nset; i++) |
2432 | { |
2433 | edat.s[i].ener[nfr] = fr->ener[set[i]].e; |
2434 | } |
2435 | } |
2436 | /* |
2437 | * Define distance restraint legends. Can only be done after |
2438 | * the first frame has been read... (Then we know how many there are) |
2439 | */ |
2440 | blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL((void*)0)); |
2441 | if (bDisRe && bDRAll && !leg && blk_disre) |
2442 | { |
2443 | t_iatom *fa; |
2444 | t_iparams *ip; |
2445 | |
2446 | fa = top->idef.il[F_DISRES].iatoms; |
2447 | ip = top->idef.iparams; |
2448 | if (blk_disre->nsub != 2 || |
2449 | (blk_disre->sub[0].nr != blk_disre->sub[1].nr) ) |
2450 | { |
2451 | gmx_incons("Number of disre sub-blocks not equal to 2")_gmx_error("incons", "Number of disre sub-blocks not equal to 2" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2451); |
2452 | } |
2453 | |
2454 | ndisre = blk_disre->sub[0].nr; |
2455 | if (ndisre != top->idef.il[F_DISRES].nr/3) |
2456 | { |
2457 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2457, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n", |
2458 | ndisre, top->idef.il[F_DISRES].nr/3); |
2459 | } |
2460 | snew(pairleg, ndisre)(pairleg) = save_calloc("pairleg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2460, (ndisre), sizeof(*(pairleg))); |
2461 | for (i = 0; i < ndisre; i++) |
2462 | { |
2463 | snew(pairleg[i], 30)(pairleg[i]) = save_calloc("pairleg[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2463, (30), sizeof(*(pairleg[i]))); |
2464 | j = fa[3*i+1]; |
2465 | k = fa[3*i+2]; |
2466 | gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j); |
2467 | gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k); |
2468 | sprintf(pairleg[i], "%d %s %d %s (%d)", |
2469 | resnr_j, anm_j, resnr_k, anm_k, |
2470 | ip[fa[3*i]].disres.label); |
2471 | } |
2472 | set = select_it(ndisre, pairleg, &nset); |
2473 | snew(leg, 2*nset)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2473, (2*nset), sizeof(*(leg))); |
2474 | for (i = 0; (i < nset); i++) |
2475 | { |
2476 | snew(leg[2*i], 32)(leg[2*i]) = save_calloc("leg[2*i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2476, (32), sizeof(*(leg[2*i]))); |
2477 | sprintf(leg[2*i], "a %s", pairleg[set[i]]); |
2478 | snew(leg[2*i+1], 32)(leg[2*i+1]) = save_calloc("leg[2*i+1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2478, (32), sizeof(*(leg[2*i+1]))); |
2479 | sprintf(leg[2*i+1], "i %s", pairleg[set[i]]); |
2480 | } |
2481 | xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv); |
2482 | } |
2483 | |
2484 | /* |
2485 | * Store energies for analysis afterwards... |
2486 | */ |
2487 | if (!bDisRe && !bDHDL && (fr->nre > 0)) |
2488 | { |
2489 | if (edat.nframes % 1000 == 0) |
2490 | { |
2491 | srenew(time, edat.nframes+1000)(time) = save_realloc("time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2491, (time), (edat.nframes+1000), sizeof(*(time))); |
2492 | } |
2493 | time[edat.nframes] = fr->t; |
2494 | edat.nframes++; |
2495 | } |
2496 | /* |
2497 | * Printing time, only when we do not want to skip frames |
2498 | */ |
2499 | if (!skip || teller % skip == 0) |
2500 | { |
2501 | if (bDisRe) |
2502 | { |
2503 | /******************************************* |
2504 | * D I S T A N C E R E S T R A I N T S |
2505 | *******************************************/ |
2506 | if (ndisre > 0) |
2507 | { |
2508 | #ifndef GMX_DOUBLE |
2509 | float *disre_rt = blk_disre->sub[0].fval; |
2510 | float *disre_rm3tav = blk_disre->sub[1].fval; |
2511 | #else |
2512 | double *disre_rt = blk_disre->sub[0].dval; |
2513 | double *disre_rm3tav = blk_disre->sub[1].dval; |
2514 | #endif |
2515 | |
2516 | print_time(out, fr->t); |
2517 | if (violaver == NULL((void*)0)) |
2518 | { |
2519 | snew(violaver, ndisre)(violaver) = save_calloc("violaver", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2519, (ndisre), sizeof(*(violaver))); |
2520 | } |
2521 | |
2522 | /* Subtract bounds from distances, to calculate violations */ |
2523 | calc_violations(disre_rt, disre_rm3tav, |
2524 | nbounds, pair, bounds, violaver, &sumt, &sumaver); |
2525 | |
2526 | fprintf(out, " %8.4f %8.4f\n", sumaver, sumt); |
2527 | if (bDRAll) |
2528 | { |
2529 | print_time(fp_pairs, fr->t); |
2530 | for (i = 0; (i < nset); i++) |
2531 | { |
2532 | sss = set[i]; |
2533 | fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird)); |
2534 | fprintf(fp_pairs, " %8.4f", disre_rt[sss]); |
2535 | } |
2536 | fprintf(fp_pairs, "\n"); |
2537 | } |
2538 | teller_disre++; |
2539 | } |
2540 | } |
2541 | else if (bDHDL) |
2542 | { |
2543 | do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv); |
2544 | } |
2545 | |
2546 | /******************************************* |
2547 | * E N E R G I E S |
2548 | *******************************************/ |
2549 | else |
2550 | { |
2551 | if (fr->nre > 0) |
2552 | { |
2553 | if (bPrAll) |
2554 | { |
2555 | /* We skip frames with single points (usually only the first frame), |
2556 | * since they would result in an average plot with outliers. |
2557 | */ |
2558 | if (fr->nsum > 1) |
2559 | { |
2560 | print_time(out, fr->t); |
2561 | print1(out, bDp, fr->ener[set[0]].e); |
2562 | print1(out, bDp, fr->ener[set[0]].esum/fr->nsum); |
2563 | print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum)); |
2564 | fprintf(out, "\n"); |
2565 | } |
2566 | } |
2567 | else |
2568 | { |
2569 | print_time(out, fr->t); |
2570 | if (bSum) |
2571 | { |
2572 | sum = 0; |
2573 | for (i = 0; i < nset; i++) |
2574 | { |
2575 | sum += fr->ener[set[i]].e; |
2576 | } |
2577 | print1(out, bDp, sum/nmol-ezero); |
2578 | } |
2579 | else |
2580 | { |
2581 | for (i = 0; (i < nset); i++) |
2582 | { |
2583 | if (bIsEner[i]) |
2584 | { |
2585 | print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero); |
2586 | } |
2587 | else |
2588 | { |
2589 | print1(out, bDp, fr->ener[set[i]].e); |
2590 | } |
2591 | } |
2592 | } |
2593 | fprintf(out, "\n"); |
2594 | } |
2595 | } |
2596 | blk = find_block_id_enxframe(fr, enx_i, NULL((void*)0)); |
2597 | if (bORIRE && blk) |
2598 | { |
2599 | #ifndef GMX_DOUBLE |
2600 | xdr_datatype dt = xdr_datatype_float; |
2601 | #else |
2602 | xdr_datatype dt = xdr_datatype_double; |
2603 | #endif |
2604 | real *vals; |
2605 | |
2606 | if ( (blk->nsub != 1) || (blk->sub[0].type != dt) ) |
2607 | { |
2608 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2608, "Orientational restraints read in incorrectly"); |
2609 | } |
2610 | #ifndef GMX_DOUBLE |
2611 | vals = blk->sub[0].fval; |
2612 | #else |
2613 | vals = blk->sub[0].dval; |
2614 | #endif |
2615 | |
2616 | if (blk->sub[0].nr != (size_t)nor) |
2617 | { |
2618 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2618, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr); |
2619 | } |
2620 | if (bORA || bODA) |
2621 | { |
2622 | for (i = 0; i < nor; i++) |
2623 | { |
2624 | orient[i] += vals[i]; |
2625 | } |
2626 | } |
2627 | if (bODR) |
2628 | { |
2629 | for (i = 0; i < nor; i++) |
2630 | { |
2631 | odrms[i] += sqr(vals[i]-oobs[i]); |
2632 | } |
2633 | } |
2634 | if (bORT) |
2635 | { |
2636 | fprintf(fort, " %10f", fr->t); |
2637 | for (i = 0; i < norsel; i++) |
2638 | { |
2639 | fprintf(fort, " %g", vals[orsel[i]]); |
2640 | } |
2641 | fprintf(fort, "\n"); |
2642 | } |
2643 | if (bODT) |
2644 | { |
2645 | fprintf(fodt, " %10f", fr->t); |
2646 | for (i = 0; i < norsel; i++) |
2647 | { |
2648 | fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]); |
2649 | } |
2650 | fprintf(fodt, "\n"); |
2651 | } |
2652 | norfr++; |
2653 | } |
2654 | blk = find_block_id_enxframe(fr, enxORT, NULL((void*)0)); |
2655 | if (bOTEN && blk) |
2656 | { |
2657 | #ifndef GMX_DOUBLE |
2658 | xdr_datatype dt = xdr_datatype_float; |
2659 | #else |
2660 | xdr_datatype dt = xdr_datatype_double; |
2661 | #endif |
2662 | real *vals; |
2663 | |
2664 | if ( (blk->nsub != 1) || (blk->sub[0].type != dt) ) |
2665 | { |
2666 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2666, "Orientational restraints read in incorrectly"); |
2667 | } |
2668 | #ifndef GMX_DOUBLE |
2669 | vals = blk->sub[0].fval; |
2670 | #else |
2671 | vals = blk->sub[0].dval; |
2672 | #endif |
2673 | |
2674 | if (blk->sub[0].nr != (size_t)(nex*12)) |
2675 | { |
2676 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2676, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)", |
2677 | blk->sub[0].nr/12, nex); |
2678 | } |
2679 | fprintf(foten, " %10f", fr->t); |
2680 | for (i = 0; i < nex; i++) |
2681 | { |
2682 | for (j = 0; j < (bOvec ? 12 : 3); j++) |
2683 | { |
2684 | fprintf(foten, " %g", vals[i*12+j]); |
2685 | } |
2686 | } |
2687 | fprintf(foten, "\n"); |
2688 | } |
2689 | } |
2690 | } |
2691 | teller++; |
2692 | } |
2693 | } |
2694 | while (bCont && (timecheck == 0)); |
2695 | |
2696 | fprintf(stderrstderr, "\n"); |
2697 | close_enx(fp); |
2698 | if (out) |
2699 | { |
2700 | gmx_ffclose(out); |
2701 | } |
2702 | |
2703 | if (bDRAll) |
2704 | { |
2705 | gmx_ffclose(fp_pairs); |
2706 | } |
2707 | |
2708 | if (bORT) |
2709 | { |
2710 | gmx_ffclose(fort); |
2711 | } |
2712 | if (bODT) |
2713 | { |
2714 | gmx_ffclose(fodt); |
2715 | } |
2716 | if (bORA) |
2717 | { |
2718 | out = xvgropen(opt2fn("-ora", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2719 | "Average calculated orientations", |
2720 | "Restraint label", "", oenv); |
2721 | if (bOrinst) |
2722 | { |
2723 | fprintf(out, "%s", orinst_sub); |
2724 | } |
2725 | for (i = 0; i < nor; i++) |
2726 | { |
2727 | fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr); |
2728 | } |
2729 | gmx_ffclose(out); |
2730 | } |
2731 | if (bODA) |
2732 | { |
2733 | out = xvgropen(opt2fn("-oda", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2734 | "Average restraint deviation", |
2735 | "Restraint label", "", oenv); |
2736 | if (bOrinst) |
2737 | { |
2738 | fprintf(out, "%s", orinst_sub); |
2739 | } |
2740 | for (i = 0; i < nor; i++) |
2741 | { |
2742 | fprintf(out, "%5d %g\n", or_label[i], orient[i]/norfr-oobs[i]); |
2743 | } |
2744 | gmx_ffclose(out); |
2745 | } |
2746 | if (bODR) |
2747 | { |
2748 | out = xvgropen(opt2fn("-odr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2749 | "RMS orientation restraint deviations", |
2750 | "Restraint label", "", oenv); |
2751 | if (bOrinst) |
2752 | { |
2753 | fprintf(out, "%s", orinst_sub); |
2754 | } |
2755 | for (i = 0; i < nor; i++) |
2756 | { |
2757 | fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr)); |
2758 | } |
2759 | gmx_ffclose(out); |
2760 | } |
2761 | if (bOTEN) |
2762 | { |
2763 | gmx_ffclose(foten); |
2764 | } |
2765 | |
2766 | if (bDisRe) |
2767 | { |
2768 | analyse_disre(opt2fn("-viol", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2769 | teller_disre, violaver, bounds, index, pair, nbounds, oenv); |
2770 | } |
2771 | else if (bDHDL) |
2772 | { |
2773 | if (fp_dhdl) |
2774 | { |
2775 | gmx_ffclose(fp_dhdl); |
2776 | printf("\n\nWrote %d lambda values with %d samples as ", |
2777 | dh_lambdas, dh_samples); |
2778 | if (dh_hists > 0) |
2779 | { |
2780 | printf("%d dH histograms ", dh_hists); |
2781 | } |
2782 | if (dh_blocks > 0) |
2783 | { |
2784 | printf("%d dH data blocks ", dh_blocks); |
2785 | } |
2786 | printf("to %s\n", opt2fn("-odh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
2787 | |
2788 | } |
2789 | else |
2790 | { |
2791 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_energy.c" , 2791, "No dH data in %s\n", opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
2792 | } |
2793 | |
2794 | } |
2795 | else |
2796 | { |
2797 | double dt = (frame[cur].t-start_t)/(edat.nframes-1); |
2798 | analyse_ener(opt2bSet("-corr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-corr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2799 | bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa), |
2800 | bVisco, opt2fn("-vis", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2801 | nmol, |
2802 | start_step, start_t, frame[cur].step, frame[cur].t, |
2803 | time, reftemp, &edat, |
2804 | nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax, |
2805 | oenv); |
2806 | if (bFluctProps) |
2807 | { |
2808 | calc_fluctuation_props(stdoutstdout, bDriftCorr, dt, nset, nmol, leg, &edat, |
2809 | nbmin, nbmax); |
2810 | } |
2811 | } |
2812 | if (opt2bSet("-f2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
2813 | { |
2814 | fec(opt2fn("-f2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-ravg", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
2815 | reftemp, nset, set, leg, &edat, time, oenv); |
2816 | } |
2817 | |
2818 | { |
2819 | const char *nxy = "-nxy"; |
2820 | |
2821 | do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2822 | do_view(oenv, opt2fn_null("-ravg", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2823 | do_view(oenv, opt2fn_null("-ora", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2824 | do_view(oenv, opt2fn_null("-ort", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2825 | do_view(oenv, opt2fn_null("-oda", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2826 | do_view(oenv, opt2fn_null("-odr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2827 | do_view(oenv, opt2fn_null("-odt", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2828 | do_view(oenv, opt2fn_null("-oten", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2829 | do_view(oenv, opt2fn_null("-odh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), nxy); |
2830 | } |
2831 | |
2832 | return 0; |
2833 | } |