File: | gromacs/gmxana/gmx_energy.c |
Location: | line 2509, column 52 |
Description: | Access to field 'sub' results in a dereference of a null pointer (loaded from variable 'blk_disre') |
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; | |||
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 | } |