File: | gromacs/gmxana/gmx_sham.c |
Location: | line 1159, column 5 |
Description: | Function call argument is an uninitialized value |
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 "gromacs/commandline/pargs.h" | |||
46 | #include "typedefs.h" | |||
47 | #include "gromacs/utility/smalloc.h" | |||
48 | #include "macros.h" | |||
49 | #include "gromacs/utility/fatalerror.h" | |||
50 | #include "gromacs/math/vec.h" | |||
51 | #include "gromacs/utility/futil.h" | |||
52 | #include "readinp.h" | |||
53 | #include "txtdump.h" | |||
54 | #include "gstat.h" | |||
55 | #include "gromacs/fileio/xvgr.h" | |||
56 | #include "physics.h" | |||
57 | #include "gromacs/fileio/pdbio.h" | |||
58 | #include "gromacs/fileio/matio.h" | |||
59 | #include "gmx_ana.h" | |||
60 | ||||
61 | ||||
62 | static int index2(int *ibox, int x, int y) | |||
63 | { | |||
64 | return (ibox[1]*x+y); | |||
65 | } | |||
66 | ||||
67 | static int index3(int *ibox, int x, int y, int z) | |||
68 | { | |||
69 | return (ibox[2]*(ibox[1]*x+y)+z); | |||
70 | } | |||
71 | ||||
72 | static gmx_int64_t indexn(int ndim, const int *ibox, const int *nxyz) | |||
73 | { | |||
74 | gmx_int64_t d, dd; | |||
75 | int k, kk; | |||
76 | ||||
77 | /* Compute index in 1-D array */ | |||
78 | d = 0; | |||
79 | for (k = 0; (k < ndim); k++) | |||
80 | { | |||
81 | dd = nxyz[k]; | |||
82 | for (kk = k+1; (kk < ndim); kk++) | |||
83 | { | |||
84 | dd = dd*ibox[kk]; | |||
85 | } | |||
86 | d += dd; | |||
87 | } | |||
88 | return d; | |||
89 | } | |||
90 | ||||
91 | typedef struct { | |||
92 | int Nx; /* x grid points in unit cell */ | |||
93 | int Ny; /* y grid points in unit cell */ | |||
94 | int Nz; /* z grid points in unit cell */ | |||
95 | int dmin[3]; /* starting point x,y,z*/ | |||
96 | int dmax[3]; /* ending point x,y,z */ | |||
97 | real cell[6]; /* usual cell parameters */ | |||
98 | real * ed; /* data */ | |||
99 | } XplorMap; | |||
100 | ||||
101 | static void lo_write_xplor(XplorMap * map, const char * file) | |||
102 | { | |||
103 | FILE * fp; | |||
104 | int z, i, j, n; | |||
105 | ||||
106 | fp = gmx_ffopen(file, "w"); | |||
107 | /* The REMARKS part is the worst part of the XPLOR format | |||
108 | * and may cause problems with some programs | |||
109 | */ | |||
110 | fprintf(fp, "\n 2 !NTITLE\n"); | |||
111 | fprintf(fp, " REMARKS Energy Landscape from GROMACS\n"); | |||
112 | fprintf(fp, " REMARKS DATE: 2004-12-21 \n"); | |||
113 | fprintf(fp, " %7d %7d %7d %7d %7d %7d %7d %7d %7d\n", | |||
114 | map->Nx, map->dmin[0], map->dmax[0], | |||
115 | map->Ny, map->dmin[1], map->dmax[1], | |||
116 | map->Nz, map->dmin[2], map->dmax[2]); | |||
117 | fprintf(fp, "%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n", | |||
118 | map->cell[0], map->cell[1], map->cell[2], | |||
119 | map->cell[3], map->cell[4], map->cell[5]); | |||
120 | fprintf(fp, "ZYX\n"); | |||
121 | ||||
122 | z = map->dmin[2]; | |||
123 | for (n = 0; n < map->Nz; n++, z++) | |||
124 | { | |||
125 | fprintf(fp, "%8d\n", z); | |||
126 | for (i = 0; i < map->Nx*map->Ny; i += 6) | |||
127 | { | |||
128 | for (j = 0; j < 6; j++) | |||
129 | { | |||
130 | if (i+j < map->Nx*map->Ny) | |||
131 | { | |||
132 | fprintf(fp, "%12.5E", map->ed[n*map->Nx*map->Ny+i+j]); | |||
133 | } | |||
134 | } | |||
135 | fprintf(fp, "\n"); | |||
136 | } | |||
137 | } | |||
138 | fprintf(fp, " -9999\n"); | |||
139 | gmx_ffclose(fp); | |||
140 | } | |||
141 | ||||
142 | static void write_xplor(const char *file, real *data, int *ibox, real dmin[], real dmax[]) | |||
143 | { | |||
144 | XplorMap *xm; | |||
145 | int i, j, k, n; | |||
146 | ||||
147 | snew(xm, 1)(xm) = save_calloc("xm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 147, (1), sizeof(*(xm))); | |||
148 | xm->Nx = ibox[XX0]; | |||
149 | xm->Ny = ibox[YY1]; | |||
150 | xm->Nz = ibox[ZZ2]; | |||
151 | snew(xm->ed, xm->Nx*xm->Ny*xm->Nz)(xm->ed) = save_calloc("xm->ed", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 151, (xm->Nx*xm->Ny*xm->Nz), sizeof(*(xm->ed))); | |||
152 | n = 0; | |||
153 | for (k = 0; (k < xm->Nz); k++) | |||
154 | { | |||
155 | for (j = 0; (j < xm->Ny); j++) | |||
156 | { | |||
157 | for (i = 0; (i < xm->Nx); i++) | |||
158 | { | |||
159 | xm->ed[n++] = data[index3(ibox, i, j, k)]; | |||
160 | } | |||
161 | } | |||
162 | } | |||
163 | xm->cell[0] = dmax[XX0]-dmin[XX0]; | |||
164 | xm->cell[1] = dmax[YY1]-dmin[YY1]; | |||
165 | xm->cell[2] = dmax[ZZ2]-dmin[ZZ2]; | |||
166 | xm->cell[3] = xm->cell[4] = xm->cell[5] = 90; | |||
167 | ||||
168 | clear_ivec(xm->dmin); | |||
169 | xm->dmax[XX0] = ibox[XX0]-1; | |||
170 | xm->dmax[YY1] = ibox[YY1]-1; | |||
171 | xm->dmax[ZZ2] = ibox[ZZ2]-1; | |||
172 | ||||
173 | lo_write_xplor(xm, file); | |||
174 | ||||
175 | sfree(xm->ed)save_free("xm->ed", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 175, (xm->ed)); | |||
176 | sfree(xm)save_free("xm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 176, (xm)); | |||
177 | } | |||
178 | ||||
179 | static void normalize_p_e(int len, double *P, int *nbin, real *E, real pmin) | |||
180 | { | |||
181 | int i; | |||
182 | double Ptot = 0; | |||
183 | ||||
184 | for (i = 0; (i < len); i++) | |||
185 | { | |||
186 | Ptot += P[i]; | |||
187 | if (nbin[i] > 0) | |||
188 | { | |||
189 | E[i] = E[i]/nbin[i]; | |||
190 | } | |||
191 | } | |||
192 | printf("Ptot = %g\n", Ptot); | |||
193 | for (i = 0; (i < len); i++) | |||
194 | { | |||
195 | P[i] = P[i]/Ptot; | |||
196 | /* Have to check for pmin after normalizing to prevent "stretching" | |||
197 | * the energies. | |||
198 | */ | |||
199 | if (P[i] < pmin) | |||
200 | { | |||
201 | P[i] = 0; | |||
202 | } | |||
203 | } | |||
204 | } | |||
205 | ||||
206 | typedef struct { | |||
207 | gmx_int64_t index; | |||
208 | real ener; | |||
209 | } t_minimum; | |||
210 | ||||
211 | static int comp_minima(const void *a, const void *b) | |||
212 | { | |||
213 | t_minimum *ma = (t_minimum *) a; | |||
214 | t_minimum *mb = (t_minimum *) b; | |||
215 | ||||
216 | if (ma->ener < mb->ener) | |||
217 | { | |||
218 | return -1; | |||
219 | } | |||
220 | else if (ma->ener > mb->ener) | |||
221 | { | |||
222 | return 1; | |||
223 | } | |||
224 | else | |||
225 | { | |||
226 | return 0; | |||
227 | } | |||
228 | } | |||
229 | ||||
230 | static gmx_inlineinline | |||
231 | void print_minimum(FILE *fp, int num, const t_minimum *min) | |||
232 | { | |||
233 | fprintf(fp, | |||
234 | "Minimum %d at index " "%"GMX_PRId64"l" "d" " energy %10.3f\n", | |||
235 | num, min->index, min->ener); | |||
236 | } | |||
237 | ||||
238 | static gmx_inlineinline | |||
239 | void add_minimum(FILE *fp, int num, const t_minimum *min, t_minimum *mm) | |||
240 | { | |||
241 | print_minimum(fp, num, min); | |||
242 | mm[num].index = min->index; | |||
243 | mm[num].ener = min->ener; | |||
244 | } | |||
245 | ||||
246 | static gmx_inlineinline | |||
247 | gmx_bool is_local_minimum_from_below(const t_minimum *this_min, | |||
248 | int dimension_index, | |||
249 | int dimension_min, | |||
250 | int neighbour_index, | |||
251 | real *W) | |||
252 | { | |||
253 | return ((dimension_index == dimension_min) || | |||
254 | ((dimension_index > dimension_min) && | |||
255 | (this_min->ener < W[neighbour_index]))); | |||
256 | /* Note over/underflow within W cannot occur. */ | |||
257 | } | |||
258 | ||||
259 | static gmx_inlineinline | |||
260 | gmx_bool is_local_minimum_from_above(const t_minimum *this_min, | |||
261 | int dimension_index, | |||
262 | int dimension_max, | |||
263 | int neighbour_index, | |||
264 | real *W) | |||
265 | { | |||
266 | return ((dimension_index == dimension_max) || | |||
267 | ((dimension_index < dimension_max) && | |||
268 | (this_min->ener < W[neighbour_index]))); | |||
269 | /* Note over/underflow within W cannot occur. */ | |||
270 | } | |||
271 | ||||
272 | static void pick_minima(const char *logfile, int *ibox, int ndim, int len, real W[]) | |||
273 | { | |||
274 | FILE *fp; | |||
275 | int i, j, k, nmin; | |||
276 | t_minimum *mm, this_min; | |||
277 | int *this_point; | |||
278 | int loopmax, loopcounter; | |||
279 | ||||
280 | snew(mm, len)(mm) = save_calloc("mm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 280, (len), sizeof(*(mm))); | |||
281 | nmin = 0; | |||
282 | fp = gmx_ffopen(logfile, "w"); | |||
283 | /* Loop over each element in the array of dimenion ndim seeking | |||
284 | * minima with respect to every dimension. Specialized loops for | |||
285 | * speed with ndim == 2 and ndim == 3. */ | |||
286 | switch (ndim) | |||
287 | { | |||
288 | case 0: | |||
289 | /* This is probably impossible to reach anyway. */ | |||
290 | break; | |||
291 | case 2: | |||
292 | for (i = 0; (i < ibox[0]); i++) | |||
293 | { | |||
294 | for (j = 0; (j < ibox[1]); j++) | |||
295 | { | |||
296 | /* Get the index of this point in the flat array */ | |||
297 | this_min.index = index2(ibox, i, j); | |||
298 | this_min.ener = W[this_min.index]; | |||
299 | if (is_local_minimum_from_below(&this_min, i, 0, index2(ibox, i-1, j ), W) && | |||
300 | is_local_minimum_from_above(&this_min, i, ibox[0]-1, index2(ibox, i+1, j ), W) && | |||
301 | is_local_minimum_from_below(&this_min, j, 0, index2(ibox, i, j-1), W) && | |||
302 | is_local_minimum_from_above(&this_min, j, ibox[1]-1, index2(ibox, i, j+1), W)) | |||
303 | { | |||
304 | add_minimum(fp, nmin, &this_min, mm); | |||
305 | nmin++; | |||
306 | } | |||
307 | } | |||
308 | } | |||
309 | break; | |||
310 | case 3: | |||
311 | for (i = 0; (i < ibox[0]); i++) | |||
312 | { | |||
313 | for (j = 0; (j < ibox[1]); j++) | |||
314 | { | |||
315 | for (k = 0; (k < ibox[2]); k++) | |||
316 | { | |||
317 | /* Get the index of this point in the flat array */ | |||
318 | this_min.index = index3(ibox, i, j, k); | |||
319 | this_min.ener = W[this_min.index]; | |||
320 | if (is_local_minimum_from_below(&this_min, i, 0, index3(ibox, i-1, j, k ), W) && | |||
321 | is_local_minimum_from_above(&this_min, i, ibox[0]-1, index3(ibox, i+1, j, k ), W) && | |||
322 | is_local_minimum_from_below(&this_min, j, 0, index3(ibox, i, j-1, k ), W) && | |||
323 | is_local_minimum_from_above(&this_min, j, ibox[1]-1, index3(ibox, i, j+1, k ), W) && | |||
324 | is_local_minimum_from_below(&this_min, k, 0, index3(ibox, i, j, k-1), W) && | |||
325 | is_local_minimum_from_above(&this_min, k, ibox[2]-1, index3(ibox, i, j, k+1), W)) | |||
326 | { | |||
327 | add_minimum(fp, nmin, &this_min, mm); | |||
328 | nmin++; | |||
329 | } | |||
330 | } | |||
331 | } | |||
332 | } | |||
333 | break; | |||
334 | default: | |||
335 | /* Note this treats ndim == 1 and ndim > 3 */ | |||
336 | ||||
337 | /* Set up an ndim-dimensional vector to loop over the points | |||
338 | * on the grid. (0,0,0, ... 0) is an acceptable place to | |||
339 | * start. */ | |||
340 | snew(this_point, ndim)(this_point) = save_calloc("this_point", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 340, (ndim), sizeof(*(this_point))); | |||
341 | ||||
342 | /* Determine the number of points of the ndim-dimensional | |||
343 | * grid. */ | |||
344 | loopmax = ibox[0]; | |||
345 | for (i = 1; i < ndim; i++) | |||
346 | { | |||
347 | loopmax *= ibox[i]; | |||
348 | } | |||
349 | ||||
350 | loopcounter = 0; | |||
351 | while (loopmax > loopcounter) | |||
352 | { | |||
353 | gmx_bool bMin = TRUE1; | |||
354 | ||||
355 | /* Get the index of this_point in the flat array */ | |||
356 | this_min.index = indexn(ndim, ibox, this_point); | |||
357 | this_min.ener = W[this_min.index]; | |||
358 | ||||
359 | /* Is this_point a minimum from above and below in each | |||
360 | * dimension? */ | |||
361 | for (i = 0; bMin && (i < ndim); i++) | |||
362 | { | |||
363 | /* Save the index of this_point within the curent | |||
364 | * dimension so we can change that index in the | |||
365 | * this_point array for use with indexn(). */ | |||
366 | int index = this_point[i]; | |||
367 | this_point[i]--; | |||
368 | bMin = bMin && | |||
369 | is_local_minimum_from_below(&this_min, index, 0, indexn(ndim, ibox, this_point), W); | |||
370 | this_point[i] += 2; | |||
371 | bMin = bMin && | |||
372 | is_local_minimum_from_above(&this_min, index, ibox[i]-1, indexn(ndim, ibox, this_point), W); | |||
373 | this_point[i]--; | |||
374 | } | |||
375 | if (bMin) | |||
376 | { | |||
377 | add_minimum(fp, nmin, &this_min, mm); | |||
378 | nmin++; | |||
379 | } | |||
380 | ||||
381 | /* update global loop counter */ | |||
382 | loopcounter++; | |||
383 | ||||
384 | /* Avoid underflow of this_point[i] */ | |||
385 | if (loopmax > loopcounter) | |||
386 | { | |||
387 | /* update this_point non-recursively */ | |||
388 | i = ndim-1; | |||
389 | this_point[i]++; | |||
390 | while (ibox[i] == this_point[i]) | |||
391 | { | |||
392 | this_point[i] = 0; | |||
393 | i--; | |||
394 | /* this_point[i] cannot underflow because | |||
395 | * loopmax > loopcounter. */ | |||
396 | this_point[i]++; | |||
397 | } | |||
398 | } | |||
399 | } | |||
400 | ||||
401 | sfree(this_point)save_free("this_point", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 401, (this_point)); | |||
402 | break; | |||
403 | } | |||
404 | qsort(mm, nmin, sizeof(mm[0]), comp_minima); | |||
405 | fprintf(fp, "Minima sorted after energy\n"); | |||
406 | for (i = 0; (i < nmin); i++) | |||
407 | { | |||
408 | print_minimum(fp, i, &mm[i]); | |||
409 | } | |||
410 | gmx_ffclose(fp); | |||
411 | sfree(mm)save_free("mm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 411, (mm)); | |||
412 | } | |||
413 | ||||
414 | static void do_sham(const char *fn, const char *ndx, | |||
415 | const char *xpmP, const char *xpm, const char *xpm2, | |||
416 | const char *xpm3, const char *xpm4, const char *pdb, | |||
417 | const char *logf, | |||
418 | int n, int neig, real **eig, | |||
419 | gmx_bool bGE, int nenerT, real **enerT, | |||
420 | int nmap, real *mapindex, real **map, | |||
421 | real Tref, | |||
422 | real pmax, real gmax, | |||
423 | real *emin, real *emax, int nlevels, real pmin, | |||
424 | const char *mname, int *idim, int *ibox, | |||
425 | gmx_bool bXmin, real *xmin, gmx_bool bXmax, real *xmax) | |||
426 | { | |||
427 | FILE *fp; | |||
428 | real *min_eig, *max_eig; | |||
429 | real *axis_x, *axis_y, *axis_z, *axis = NULL((void*)0); | |||
430 | double *P; | |||
431 | real **PP, *W, *E, **WW, **EE, *S, **SS, *M, **MM, *bE; | |||
432 | rvec xxx; | |||
433 | char *buf; | |||
434 | double *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf, Mmin, Mmax, Minf; | |||
435 | real *delta; | |||
436 | int i, j, k, imin, len, index, d, *nbin, *bindex, bi; | |||
437 | int *nxyz, maxbox; | |||
438 | t_blocka *b; | |||
439 | gmx_bool bOutside; | |||
440 | unsigned int flags; | |||
441 | t_rgb rlo = { 0, 0, 0 }; | |||
442 | t_rgb rhi = { 1, 1, 1 }; | |||
443 | ||||
444 | /* Determine extremes for the eigenvectors */ | |||
445 | snew(min_eig, neig)(min_eig) = save_calloc("min_eig", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 445, (neig), sizeof(*(min_eig))); | |||
446 | snew(max_eig, neig)(max_eig) = save_calloc("max_eig", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 446, (neig), sizeof(*(max_eig))); | |||
447 | snew(nxyz, neig)(nxyz) = save_calloc("nxyz", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 447, (neig), sizeof(*(nxyz))); | |||
448 | snew(bfac, neig)(bfac) = save_calloc("bfac", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 448, (neig), sizeof(*(bfac))); | |||
449 | snew(delta, neig)(delta) = save_calloc("delta", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 449, (neig), sizeof(*(delta))); | |||
450 | ||||
451 | for (i = 0; (i < neig); i++) | |||
452 | { | |||
453 | /* Check for input constraints */ | |||
454 | min_eig[i] = max_eig[i] = eig[i][0]; | |||
455 | for (j = 0; (j < n); j++) | |||
456 | { | |||
457 | min_eig[i] = min(min_eig[i], eig[i][j])(((min_eig[i]) < (eig[i][j])) ? (min_eig[i]) : (eig[i][j]) ); | |||
458 | max_eig[i] = max(max_eig[i], eig[i][j])(((max_eig[i]) > (eig[i][j])) ? (max_eig[i]) : (eig[i][j]) ); | |||
459 | delta[i] = (max_eig[i]-min_eig[i])/(2.0*ibox[i]); | |||
460 | } | |||
461 | /* Add some extra space, half a bin on each side, unless the | |||
462 | * user has set the limits. | |||
463 | */ | |||
464 | if (bXmax) | |||
465 | { | |||
466 | if (max_eig[i] > xmax[i]) | |||
467 | { | |||
468 | gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f", i, xmax[i], max_eig[i]); | |||
469 | } | |||
470 | max_eig[i] = xmax[i]; | |||
471 | } | |||
472 | else | |||
473 | { | |||
474 | max_eig[i] += delta[i]; | |||
475 | } | |||
476 | ||||
477 | if (bXmin) | |||
478 | { | |||
479 | if (min_eig[i] < xmin[i]) | |||
480 | { | |||
481 | gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f", i, xmin[i], min_eig[i]); | |||
482 | } | |||
483 | min_eig[i] = xmin[i]; | |||
484 | } | |||
485 | else | |||
486 | { | |||
487 | min_eig[i] -= delta[i]; | |||
488 | } | |||
489 | bfac[i] = ibox[i]/(max_eig[i]-min_eig[i]); | |||
490 | } | |||
491 | /* Do the binning */ | |||
492 | bref = 1/(BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*Tref); | |||
493 | snew(bE, n)(bE) = save_calloc("bE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 493, (n), sizeof(*(bE))); | |||
494 | if (bGE || nenerT == 2) | |||
495 | { | |||
496 | Emin = 1e8; | |||
497 | for (j = 0; (j < n); j++) | |||
498 | { | |||
499 | if (bGE) | |||
500 | { | |||
501 | bE[j] = bref*enerT[0][j]; | |||
502 | } | |||
503 | else | |||
504 | { | |||
505 | bE[j] = (bref - 1/(BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*enerT[1][j]))*enerT[0][j]; | |||
506 | } | |||
507 | Emin = min(Emin, bE[j])(((Emin) < (bE[j])) ? (Emin) : (bE[j]) ); | |||
508 | } | |||
509 | } | |||
510 | else | |||
511 | { | |||
512 | Emin = 0; | |||
513 | } | |||
514 | len = 1; | |||
515 | for (i = 0; (i < neig); i++) | |||
516 | { | |||
517 | len = len*ibox[i]; | |||
518 | } | |||
519 | printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n", | |||
520 | len, neig, Emin); | |||
521 | snew(P, len)(P) = save_calloc("P", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 521, (len), sizeof(*(P))); | |||
522 | snew(W, len)(W) = save_calloc("W", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 522, (len), sizeof(*(W))); | |||
523 | snew(E, len)(E) = save_calloc("E", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 523, (len), sizeof(*(E))); | |||
524 | snew(S, len)(S) = save_calloc("S", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 524, (len), sizeof(*(S))); | |||
525 | snew(M, len)(M) = save_calloc("M", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 525, (len), sizeof(*(M))); | |||
526 | snew(nbin, len)(nbin) = save_calloc("nbin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 526, (len), sizeof(*(nbin))); | |||
527 | snew(bindex, n)(bindex) = save_calloc("bindex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 527, (n), sizeof(*(bindex))); | |||
528 | ||||
529 | ||||
530 | /* Loop over projections */ | |||
531 | for (j = 0; (j < n); j++) | |||
532 | { | |||
533 | /* Loop over dimensions */ | |||
534 | bOutside = FALSE0; | |||
535 | for (i = 0; (i < neig); i++) | |||
536 | { | |||
537 | nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]); | |||
538 | if (nxyz[i] < 0 || nxyz[i] >= ibox[i]) | |||
539 | { | |||
540 | bOutside = TRUE1; | |||
541 | } | |||
542 | } | |||
543 | if (!bOutside) | |||
544 | { | |||
545 | index = indexn(neig, ibox, nxyz); | |||
546 | range_check(index, 0, len)_range_check(index, 0, len, ((void*)0),"index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 546); | |||
547 | /* Compute the exponential factor */ | |||
548 | if (enerT) | |||
549 | { | |||
550 | efac = exp(-bE[j]+Emin); | |||
551 | } | |||
552 | else | |||
553 | { | |||
554 | efac = 1; | |||
555 | } | |||
556 | /* Apply the bin volume correction for a multi-dimensional distance */ | |||
557 | for (i = 0; i < neig; i++) | |||
558 | { | |||
559 | if (idim[i] == 2) | |||
560 | { | |||
561 | efac /= eig[i][j]; | |||
562 | } | |||
563 | else if (idim[i] == 3) | |||
564 | { | |||
565 | efac /= sqr(eig[i][j]); | |||
566 | } | |||
567 | else if (idim[i] == -1) | |||
568 | { | |||
569 | efac /= sin(DEG2RAD(3.14159265358979323846/180.0)*eig[i][j]); | |||
570 | } | |||
571 | } | |||
572 | /* Update the probability */ | |||
573 | P[index] += efac; | |||
574 | /* Update the energy */ | |||
575 | if (enerT) | |||
576 | { | |||
577 | E[index] += enerT[0][j]; | |||
578 | } | |||
579 | /* Statistics: which "structure" in which bin */ | |||
580 | nbin[index]++; | |||
581 | bindex[j] = index; | |||
582 | } | |||
583 | } | |||
584 | /* Normalize probability */ | |||
585 | normalize_p_e(len, P, nbin, E, pmin); | |||
586 | Pmax = 0; | |||
587 | /* Compute boundaries for the Free energy */ | |||
588 | Wmin = 1e8; | |||
589 | imin = -1; | |||
590 | Wmax = -1e8; | |||
591 | /* Recompute Emin: it may have changed due to averaging */ | |||
592 | Emin = 1e8; | |||
593 | Emax = -1e8; | |||
594 | for (i = 0; (i < len); i++) | |||
595 | { | |||
596 | if (P[i] != 0) | |||
597 | { | |||
598 | Pmax = max(P[i], Pmax)(((P[i]) > (Pmax)) ? (P[i]) : (Pmax) ); | |||
599 | W[i] = -BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*Tref*log(P[i]); | |||
600 | if (W[i] < Wmin) | |||
601 | { | |||
602 | Wmin = W[i]; | |||
603 | imin = i; | |||
604 | } | |||
605 | Emin = min(E[i], Emin)(((E[i]) < (Emin)) ? (E[i]) : (Emin) ); | |||
606 | Emax = max(E[i], Emax)(((E[i]) > (Emax)) ? (E[i]) : (Emax) ); | |||
607 | Wmax = max(W[i], Wmax)(((W[i]) > (Wmax)) ? (W[i]) : (Wmax) ); | |||
608 | } | |||
609 | } | |||
610 | if (pmax > 0) | |||
611 | { | |||
612 | Pmax = pmax; | |||
613 | } | |||
614 | if (gmax > 0) | |||
615 | { | |||
616 | Wmax = gmax; | |||
617 | } | |||
618 | else | |||
619 | { | |||
620 | Wmax -= Wmin; | |||
621 | } | |||
622 | Winf = Wmax+1; | |||
623 | Einf = Emax+1; | |||
624 | Smin = Emin-Wmax; | |||
625 | Smax = Emax-Smin; | |||
626 | Sinf = Smax+1; | |||
627 | /* Write out the free energy as a function of bin index */ | |||
628 | fp = gmx_ffopen(fn, "w"); | |||
629 | for (i = 0; (i < len); i++) | |||
630 | { | |||
631 | if (P[i] != 0) | |||
632 | { | |||
633 | W[i] -= Wmin; | |||
634 | S[i] = E[i]-W[i]-Smin; | |||
635 | fprintf(fp, "%5d %10.5e %10.5e %10.5e\n", i, W[i], E[i], S[i]); | |||
636 | } | |||
637 | else | |||
638 | { | |||
639 | W[i] = Winf; | |||
640 | E[i] = Einf; | |||
641 | S[i] = Sinf; | |||
642 | } | |||
643 | } | |||
644 | gmx_ffclose(fp); | |||
645 | /* Organize the structures in the bins */ | |||
646 | snew(b, 1)(b) = save_calloc("b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 646, (1), sizeof(*(b))); | |||
647 | snew(b->index, len+1)(b->index) = save_calloc("b->index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 647, (len+1), sizeof(*(b->index))); | |||
648 | snew(b->a, n)(b->a) = save_calloc("b->a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 648, (n), sizeof(*(b->a))); | |||
649 | b->index[0] = 0; | |||
650 | for (i = 0; (i < len); i++) | |||
651 | { | |||
652 | b->index[i+1] = b->index[i]+nbin[i]; | |||
653 | nbin[i] = 0; | |||
654 | } | |||
655 | for (i = 0; (i < n); i++) | |||
656 | { | |||
657 | bi = bindex[i]; | |||
658 | b->a[b->index[bi]+nbin[bi]] = i; | |||
659 | nbin[bi]++; | |||
660 | } | |||
661 | /* Consistency check */ | |||
662 | /* This no longer applies when we allow the plot to be smaller | |||
663 | than the sampled space. | |||
664 | for(i=0; (i<len); i++) { | |||
665 | if (nbin[i] != (b->index[i+1] - b->index[i])) | |||
666 | gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i], | |||
667 | b->index[i+1] - b->index[i]); | |||
668 | } | |||
669 | */ | |||
670 | /* Write the index file */ | |||
671 | fp = gmx_ffopen(ndx, "w"); | |||
672 | for (i = 0; (i < len); i++) | |||
673 | { | |||
674 | if (nbin[i] > 0) | |||
675 | { | |||
676 | fprintf(fp, "[ %d ]\n", i); | |||
677 | for (j = b->index[i]; (j < b->index[i+1]); j++) | |||
678 | { | |||
679 | fprintf(fp, "%d\n", b->a[j]+1); | |||
680 | } | |||
681 | } | |||
682 | } | |||
683 | gmx_ffclose(fp); | |||
684 | snew(axis_x, ibox[0]+1)(axis_x) = save_calloc("axis_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 684, (ibox[0]+1), sizeof(*(axis_x))); | |||
685 | snew(axis_y, ibox[1]+1)(axis_y) = save_calloc("axis_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 685, (ibox[1]+1), sizeof(*(axis_y))); | |||
686 | snew(axis_z, ibox[2]+1)(axis_z) = save_calloc("axis_z", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 686, (ibox[2]+1), sizeof(*(axis_z))); | |||
687 | maxbox = max(ibox[0], max(ibox[1], ibox[2]))(((ibox[0]) > ((((ibox[1]) > (ibox[2])) ? (ibox[1]) : ( ibox[2]) ))) ? (ibox[0]) : ((((ibox[1]) > (ibox[2])) ? (ibox [1]) : (ibox[2]) )) ); | |||
688 | snew(PP, maxbox*maxbox)(PP) = save_calloc("PP", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 688, (maxbox*maxbox), sizeof(*(PP))); | |||
689 | snew(WW, maxbox*maxbox)(WW) = save_calloc("WW", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 689, (maxbox*maxbox), sizeof(*(WW))); | |||
690 | snew(EE, maxbox*maxbox)(EE) = save_calloc("EE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 690, (maxbox*maxbox), sizeof(*(EE))); | |||
691 | snew(SS, maxbox*maxbox)(SS) = save_calloc("SS", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 691, (maxbox*maxbox), sizeof(*(SS))); | |||
692 | for (i = 0; (i < min(neig, 3)(((neig) < (3)) ? (neig) : (3) )); i++) | |||
693 | { | |||
694 | switch (i) | |||
695 | { | |||
696 | case 0: axis = axis_x; break; | |||
697 | case 1: axis = axis_y; break; | |||
698 | case 2: axis = axis_z; break; | |||
699 | default: break; | |||
700 | } | |||
701 | for (j = 0; j <= ibox[i]; j++) | |||
702 | { | |||
703 | axis[j] = min_eig[i] + j/bfac[i]; | |||
704 | } | |||
705 | } | |||
706 | if (map) | |||
707 | { | |||
708 | snew(M, len)(M) = save_calloc("M", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 708, (len), sizeof(*(M))); | |||
709 | snew(MM, maxbox*maxbox)(MM) = save_calloc("MM", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 709, (maxbox*maxbox), sizeof(*(MM))); | |||
710 | for (i = 0; (i < ibox[0]); i++) | |||
711 | { | |||
712 | MM[i] = &(M[i*ibox[1]]); | |||
713 | } | |||
714 | Mmin = 1e8; | |||
715 | Mmax = -1e8; | |||
716 | for (i = 0; (i < nmap); i++) | |||
717 | { | |||
718 | Mmin = min(Mmin, map[0][i])(((Mmin) < (map[0][i])) ? (Mmin) : (map[0][i]) ); | |||
719 | Mmax = max(Mmax, map[0][i])(((Mmax) > (map[0][i])) ? (Mmax) : (map[0][i]) ); | |||
720 | } | |||
721 | Minf = Mmax*1.05; | |||
722 | for (i = 0; (i < len); i++) | |||
723 | { | |||
724 | M[i] = Minf; | |||
725 | } | |||
726 | for (i = 0; (i < nmap); i++) | |||
727 | { | |||
728 | index = gmx_nint(mapindex[i]); | |||
729 | if (index >= len) | |||
730 | { | |||
731 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 731, "Number of bins in file from -mdata option does not correspond to current analysis"); | |||
732 | } | |||
733 | ||||
734 | if (P[index] != 0) | |||
735 | { | |||
736 | M[index] = map[0][i]; | |||
737 | } | |||
738 | } | |||
739 | } | |||
740 | else | |||
741 | { | |||
742 | MM = NULL((void*)0); | |||
743 | Minf = NOTSET-12345; | |||
744 | } | |||
745 | pick_minima(logf, ibox, neig, len, W); | |||
746 | if (gmax <= 0) | |||
747 | { | |||
748 | gmax = Winf; | |||
749 | } | |||
750 | flags = MAT_SPATIAL_X(1<<0) | MAT_SPATIAL_Y(1<<1); | |||
751 | if (neig == 2) | |||
752 | { | |||
753 | /* Dump to XPM file */ | |||
754 | snew(PP, ibox[0])(PP) = save_calloc("PP", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 754, (ibox[0]), sizeof(*(PP))); | |||
755 | for (i = 0; (i < ibox[0]); i++) | |||
756 | { | |||
757 | snew(PP[i], ibox[1])(PP[i]) = save_calloc("PP[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 757, (ibox[1]), sizeof(*(PP[i]))); | |||
758 | for (j = 0; j < ibox[1]; j++) | |||
759 | { | |||
760 | PP[i][j] = P[i*ibox[1]+j]; | |||
761 | } | |||
762 | WW[i] = &(W[i*ibox[1]]); | |||
763 | EE[i] = &(E[i*ibox[1]]); | |||
764 | SS[i] = &(S[i*ibox[1]]); | |||
765 | } | |||
766 | fp = gmx_ffopen(xpmP, "w"); | |||
767 | write_xpm(fp, flags, "Probability Distribution", "", "PC1", "PC2", | |||
768 | ibox[0], ibox[1], axis_x, axis_y, PP, 0, Pmax, rlo, rhi, &nlevels); | |||
769 | gmx_ffclose(fp); | |||
770 | fp = gmx_ffopen(xpm, "w"); | |||
771 | write_xpm(fp, flags, "Gibbs Energy Landscape", "G (kJ/mol)", "PC1", "PC2", | |||
772 | ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels); | |||
773 | gmx_ffclose(fp); | |||
774 | fp = gmx_ffopen(xpm2, "w"); | |||
775 | write_xpm(fp, flags, "Enthalpy Landscape", "H (kJ/mol)", "PC1", "PC2", | |||
776 | ibox[0], ibox[1], axis_x, axis_y, EE, | |||
777 | emin ? *emin : Emin, emax ? *emax : Einf, rlo, rhi, &nlevels); | |||
778 | gmx_ffclose(fp); | |||
779 | fp = gmx_ffopen(xpm3, "w"); | |||
780 | write_xpm(fp, flags, "Entropy Landscape", "TDS (kJ/mol)", "PC1", "PC2", | |||
781 | ibox[0], ibox[1], axis_x, axis_y, SS, 0, Sinf, rlo, rhi, &nlevels); | |||
782 | gmx_ffclose(fp); | |||
783 | if (map) | |||
784 | { | |||
785 | fp = gmx_ffopen(xpm4, "w"); | |||
786 | write_xpm(fp, flags, "Custom Landscape", mname, "PC1", "PC2", | |||
787 | ibox[0], ibox[1], axis_x, axis_y, MM, 0, Minf, rlo, rhi, &nlevels); | |||
788 | gmx_ffclose(fp); | |||
789 | } | |||
790 | } | |||
791 | else if (neig == 3) | |||
792 | { | |||
793 | /* Dump to PDB file */ | |||
794 | fp = gmx_ffopen(pdb, "w"); | |||
795 | for (i = 0; (i < ibox[0]); i++) | |||
796 | { | |||
797 | xxx[XX0] = 3*(i+0.5-ibox[0]/2); | |||
798 | for (j = 0; (j < ibox[1]); j++) | |||
799 | { | |||
800 | xxx[YY1] = 3*(j+0.5-ibox[1]/2); | |||
801 | for (k = 0; (k < ibox[2]); k++) | |||
802 | { | |||
803 | xxx[ZZ2] = 3*(k+0.5-ibox[2]/2); | |||
804 | index = index3(ibox, i, j, k); | |||
805 | if (P[index] > 0) | |||
806 | { | |||
807 | fprintf(fp, "%-6s%5u %-4.4s%3.3s %4d %8.3f%8.3f%8.3f%6.2f%6.2f\n", | |||
808 | "ATOM", (index+1) %10000, "H", "H", (index+1)%10000, | |||
809 | xxx[XX0], xxx[YY1], xxx[ZZ2], 1.0, W[index]); | |||
810 | } | |||
811 | } | |||
812 | } | |||
813 | } | |||
814 | gmx_ffclose(fp); | |||
815 | write_xplor("out.xplor", W, ibox, min_eig, max_eig); | |||
816 | if (map) | |||
817 | { | |||
818 | write_xplor("user.xplor", M, ibox, min_eig, max_eig); | |||
819 | } | |||
820 | nxyz[XX0] = imin/(ibox[1]*ibox[2]); | |||
821 | nxyz[YY1] = (imin-nxyz[XX0]*ibox[1]*ibox[2])/ibox[2]; | |||
822 | nxyz[ZZ2] = imin % ibox[2]; | |||
823 | for (i = 0; (i < ibox[0]); i++) | |||
824 | { | |||
825 | snew(WW[i], maxbox)(WW[i]) = save_calloc("WW[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 825, (maxbox), sizeof(*(WW[i]))); | |||
826 | for (j = 0; (j < ibox[1]); j++) | |||
827 | { | |||
828 | WW[i][j] = W[index3(ibox, i, j, nxyz[ZZ2])]; | |||
829 | } | |||
830 | } | |||
831 | snew(buf, strlen(xpm)+4)(buf) = save_calloc("buf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 831, (strlen(xpm)+4), sizeof(*(buf))); | |||
832 | sprintf(buf, "%s", xpm); | |||
833 | sprintf(&buf[strlen(xpm)-4], "12.xpm"); | |||
834 | fp = gmx_ffopen(buf, "w"); | |||
835 | write_xpm(fp, flags, "Gibbs Energy Landscape", "W (kJ/mol)", "PC1", "PC2", | |||
836 | ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels); | |||
837 | gmx_ffclose(fp); | |||
838 | for (i = 0; (i < ibox[0]); i++) | |||
839 | { | |||
840 | for (j = 0; (j < ibox[2]); j++) | |||
841 | { | |||
842 | WW[i][j] = W[index3(ibox, i, nxyz[YY1], j)]; | |||
843 | } | |||
844 | } | |||
845 | sprintf(&buf[strlen(xpm)-4], "13.xpm"); | |||
846 | fp = gmx_ffopen(buf, "w"); | |||
847 | write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC1", "PC3", | |||
848 | ibox[0], ibox[2], axis_x, axis_z, WW, 0, gmax, rlo, rhi, &nlevels); | |||
849 | gmx_ffclose(fp); | |||
850 | for (i = 0; (i < ibox[1]); i++) | |||
851 | { | |||
852 | for (j = 0; (j < ibox[2]); j++) | |||
853 | { | |||
854 | WW[i][j] = W[index3(ibox, nxyz[XX0], i, j)]; | |||
855 | } | |||
856 | } | |||
857 | sprintf(&buf[strlen(xpm)-4], "23.xpm"); | |||
858 | fp = gmx_ffopen(buf, "w"); | |||
859 | write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC2", "PC3", | |||
860 | ibox[1], ibox[2], axis_y, axis_z, WW, 0, gmax, rlo, rhi, &nlevels); | |||
861 | gmx_ffclose(fp); | |||
862 | sfree(buf)save_free("buf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 862, (buf)); | |||
863 | } | |||
864 | if (map) | |||
865 | { | |||
866 | sfree(MM)save_free("MM", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 866, (MM)); | |||
867 | sfree(M)save_free("M", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 867, (M)); | |||
868 | } | |||
869 | } | |||
870 | ||||
871 | static void ehisto(const char *fh, int n, real **enerT, const output_env_t oenv) | |||
872 | { | |||
873 | FILE *fp; | |||
874 | int i, j, k, nbin, blength; | |||
875 | int *bindex; | |||
876 | real *T, bmin, bmax, bwidth; | |||
877 | int **histo; | |||
878 | ||||
879 | bmin = 1e8; | |||
880 | bmax = -1e8; | |||
881 | snew(bindex, n)(bindex) = save_calloc("bindex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 881, (n), sizeof(*(bindex))); | |||
882 | snew(T, n)(T) = save_calloc("T", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 882, (n), sizeof(*(T))); | |||
883 | nbin = 0; | |||
884 | for (j = 1; (j < n); j++) | |||
885 | { | |||
886 | for (k = 0; (k < nbin); k++) | |||
887 | { | |||
888 | if (T[k] == enerT[1][j]) | |||
889 | { | |||
890 | bindex[j] = k; | |||
891 | break; | |||
892 | } | |||
893 | } | |||
894 | if (k == nbin) | |||
895 | { | |||
896 | bindex[j] = nbin; | |||
897 | T[nbin] = enerT[1][j]; | |||
898 | nbin++; | |||
899 | } | |||
900 | bmin = min(enerT[0][j], bmin)(((enerT[0][j]) < (bmin)) ? (enerT[0][j]) : (bmin) ); | |||
901 | bmax = max(enerT[0][j], bmax)(((enerT[0][j]) > (bmax)) ? (enerT[0][j]) : (bmax) ); | |||
902 | } | |||
903 | bwidth = 1.0; | |||
904 | blength = (bmax - bmin)/bwidth + 2; | |||
905 | snew(histo, nbin)(histo) = save_calloc("histo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 905, (nbin), sizeof(*(histo))); | |||
906 | for (i = 0; (i < nbin); i++) | |||
907 | { | |||
908 | snew(histo[i], blength)(histo[i]) = save_calloc("histo[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 908, (blength), sizeof(*(histo[i]))); | |||
909 | } | |||
910 | for (j = 0; (j < n); j++) | |||
911 | { | |||
912 | k = (enerT[0][j]-bmin)/bwidth; | |||
913 | histo[bindex[j]][k]++; | |||
914 | } | |||
915 | fp = xvgropen(fh, "Energy distribution", "E (kJ/mol)", "", oenv); | |||
916 | for (j = 0; (j < blength); j++) | |||
917 | { | |||
918 | fprintf(fp, "%8.3f", bmin+j*bwidth); | |||
919 | for (k = 0; (k < nbin); k++) | |||
920 | { | |||
921 | fprintf(fp, " %6d", histo[k][j]); | |||
922 | } | |||
923 | fprintf(fp, "\n"); | |||
924 | } | |||
925 | gmx_ffclose(fp); | |||
926 | } | |||
927 | ||||
928 | int gmx_sham(int argc, char *argv[]) | |||
929 | { | |||
930 | const char *desc[] = { | |||
931 | "[THISMODULE] makes multi-dimensional free-energy, enthalpy and entropy plots.", | |||
932 | "[THISMODULE] reads one or more [TT].xvg[tt] files and analyzes data sets.", | |||
933 | "The basic purpose of [THISMODULE] is to plot Gibbs free energy landscapes", | |||
934 | "(option [TT]-ls[tt])", | |||
935 | "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),", | |||
936 | "but it can also", | |||
937 | "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])", | |||
938 | "plots. The histograms can be made for any quantities the user supplies.", | |||
939 | "A line in the input file may start with a time", | |||
940 | "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.", | |||
941 | "Multiple sets can also be", | |||
942 | "read when they are separated by & (option [TT]-n[tt]),", | |||
943 | "in this case only one [IT]y[it]-value is read from each line.", | |||
944 | "All lines starting with # and @ are skipped.", | |||
945 | "[PAR]", | |||
946 | "Option [TT]-ge[tt] can be used to supply a file with free energies", | |||
947 | "when the ensemble is not a Boltzmann ensemble, but needs to be biased", | |||
948 | "by this free energy. One free energy value is required for each", | |||
949 | "(multi-dimensional) data point in the [TT]-f[tt] input.", | |||
950 | "[PAR]", | |||
951 | "Option [TT]-ene[tt] can be used to supply a file with energies.", | |||
952 | "These energies are used as a weighting function in the single", | |||
953 | "histogram analysis method by Kumar et al. When temperatures", | |||
954 | "are supplied (as a second column in the file), an experimental", | |||
955 | "weighting scheme is applied. In addition the vales", | |||
956 | "are used for making enthalpy and entropy plots.", | |||
957 | "[PAR]", | |||
958 | "With option [TT]-dim[tt], dimensions can be gives for distances.", | |||
959 | "When a distance is 2- or 3-dimensional, the circumference or surface", | |||
960 | "sampled by two particles increases with increasing distance.", | |||
961 | "Depending on what one would like to show, one can choose to correct", | |||
962 | "the histogram and free-energy for this volume effect.", | |||
963 | "The probability is normalized by r and r^2 for dimensions of 2 and 3, ", | |||
964 | "respectively.", | |||
965 | "A value of -1 is used to indicate an angle in degrees between two", | |||
966 | "vectors: a sin(angle) normalization will be applied.", | |||
967 | "[BB]Note[bb] that for angles between vectors the inner-product or cosine", | |||
968 | "is the natural quantity to use, as it will produce bins of the same", | |||
969 | "volume." | |||
970 | }; | |||
971 | static real tb = -1, te = -1, frac = 0.5, filtlen = 0; | |||
972 | static gmx_bool bHaveT = TRUE1, bDer = FALSE0, bSubAv = TRUE1, bAverCorr = FALSE0, bXYdy = FALSE0; | |||
973 | static gmx_bool bEESEF = FALSE0, bEENLC = FALSE0, bEeFitAc = FALSE0, bPower = FALSE0; | |||
974 | static gmx_bool bShamEner = TRUE1, bSham = TRUE1; | |||
975 | static real Tref = 298.15, pmin = 0, ttol = 0, pmax = 0, gmax = 0, emin = 0, emax = 0; | |||
976 | static rvec nrdim = {1, 1, 1}; | |||
977 | static rvec nrbox = {32, 32, 32}; | |||
978 | static rvec xmin = {0, 0, 0}, xmax = {1, 1, 1}; | |||
979 | static int nsets_in = 1, nb_min = 4, resol = 10, nlevels = 25; | |||
980 | static const char *mname = ""; | |||
981 | t_pargs pa[] = { | |||
982 | { "-time", FALSE0, etBOOL, {&bHaveT}, | |||
983 | "Expect a time in the input" }, | |||
984 | { "-b", FALSE0, etREAL, {&tb}, | |||
985 | "First time to read from set" }, | |||
986 | { "-e", FALSE0, etREAL, {&te}, | |||
987 | "Last time to read from set" }, | |||
988 | { "-ttol", FALSE0, etREAL, {&ttol}, | |||
989 | "Tolerance on time in appropriate units (usually ps)" }, | |||
990 | { "-n", FALSE0, etINT, {&nsets_in}, | |||
991 | "Read this number of sets separated by lines containing only an ampersand" }, | |||
992 | { "-d", FALSE0, etBOOL, {&bDer}, | |||
993 | "Use the derivative" }, | |||
994 | { "-sham", FALSE0, etBOOL, {&bSham}, | |||
995 | "Turn off energy weighting even if energies are given" }, | |||
996 | { "-tsham", FALSE0, etREAL, {&Tref}, | |||
997 | "Temperature for single histogram analysis" }, | |||
998 | { "-pmin", FALSE0, etREAL, {&pmin}, | |||
999 | "Minimum probability. Anything lower than this will be set to zero" }, | |||
1000 | { "-dim", FALSE0, etRVEC, {nrdim}, | |||
1001 | "Dimensions for distances, used for volume correction (max 3 values, dimensions > 3 will get the same value as the last)" }, | |||
1002 | { "-ngrid", FALSE0, etRVEC, {nrbox}, | |||
1003 | "Number of bins for energy landscapes (max 3 values, dimensions > 3 will get the same value as the last)" }, | |||
1004 | { "-xmin", FALSE0, etRVEC, {xmin}, | |||
1005 | "Minimum for the axes in energy landscape (see above for > 3 dimensions)" }, | |||
1006 | { "-xmax", FALSE0, etRVEC, {xmax}, | |||
1007 | "Maximum for the axes in energy landscape (see above for > 3 dimensions)" }, | |||
1008 | { "-pmax", FALSE0, etREAL, {&pmax}, | |||
1009 | "Maximum probability in output, default is calculate" }, | |||
1010 | { "-gmax", FALSE0, etREAL, {&gmax}, | |||
1011 | "Maximum free energy in output, default is calculate" }, | |||
1012 | { "-emin", FALSE0, etREAL, {&emin}, | |||
1013 | "Minimum enthalpy in output, default is calculate" }, | |||
1014 | { "-emax", FALSE0, etREAL, {&emax}, | |||
1015 | "Maximum enthalpy in output, default is calculate" }, | |||
1016 | { "-nlevels", FALSE0, etINT, {&nlevels}, | |||
1017 | "Number of levels for energy landscape" }, | |||
1018 | { "-mname", FALSE0, etSTR, {&mname}, | |||
1019 | "Legend label for the custom landscape" }, | |||
1020 | }; | |||
1021 | #define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))) | |||
1022 | ||||
1023 | FILE *out; | |||
1024 | int n, e_n, d_n, nlast, s, nset, e_nset, d_nset, i, j = 0, *idim, *ibox; | |||
| ||||
1025 | real **val, **et_val, **dt_val, *t, *e_t, e_dt, d_dt, *d_t, dt, tot, error; | |||
1026 | real *rmin, *rmax; | |||
1027 | double *av, *sig, cum1, cum2, cum3, cum4, db; | |||
1028 | const char *fn_ge, *fn_ene; | |||
1029 | output_env_t oenv; | |||
1030 | gmx_int64_t num_grid_points; | |||
1031 | ||||
1032 | t_filenm fnm[] = { | |||
1033 | { efXVG, "-f", "graph", ffREAD1<<1 }, | |||
1034 | { efXVG, "-ge", "gibbs", ffOPTRD(1<<1 | 1<<3) }, | |||
1035 | { efXVG, "-ene", "esham", ffOPTRD(1<<1 | 1<<3) }, | |||
1036 | { efXVG, "-dist", "ener", ffOPTWR(1<<2| 1<<3) }, | |||
1037 | { efXVG, "-histo", "edist", ffOPTWR(1<<2| 1<<3) }, | |||
1038 | { efNDX, "-bin", "bindex", ffOPTWR(1<<2| 1<<3) }, | |||
1039 | { efXPM, "-lp", "prob", ffOPTWR(1<<2| 1<<3) }, | |||
1040 | { efXPM, "-ls", "gibbs", ffOPTWR(1<<2| 1<<3) }, | |||
1041 | { efXPM, "-lsh", "enthalpy", ffOPTWR(1<<2| 1<<3) }, | |||
1042 | { efXPM, "-lss", "entropy", ffOPTWR(1<<2| 1<<3) }, | |||
1043 | { efXPM, "-map", "map", ffOPTWR(1<<2| 1<<3) }, | |||
1044 | { efPDB, "-ls3", "gibbs3", ffOPTWR(1<<2| 1<<3) }, | |||
1045 | { efXVG, "-mdata", "mapdata", ffOPTRD(1<<1 | 1<<3) }, | |||
1046 | { efLOG, "-g", "shamlog", ffOPTWR(1<<2| 1<<3) } | |||
1047 | }; | |||
1048 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) | |||
1049 | ||||
1050 | int npargs; | |||
1051 | ||||
1052 | npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))); | |||
1053 | if (!parse_common_args(&argc, argv, PCA_CAN_VIEW(1<<5) | PCA_BE_NICE(1<<13), | |||
1054 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) | |||
1055 | { | |||
1056 | return 0; | |||
1057 | } | |||
1058 | ||||
1059 | val = read_xvg_time(opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bHaveT, | |||
1060 | opt2parg_bSet("-b", npargs, pa), tb-ttol, | |||
1061 | opt2parg_bSet("-e", npargs, pa), te+ttol, | |||
1062 | nsets_in, &nset, &n, &dt, &t); | |||
1063 | printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt); | |||
1064 | ||||
1065 | fn_ge = opt2fn_null("-ge", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
1066 | fn_ene = opt2fn_null("-ene", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
1067 | ||||
1068 | if (fn_ge && fn_ene) | |||
1069 | { | |||
1070 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1070, "Can not do free energy and energy corrections at the same time"); | |||
1071 | } | |||
1072 | ||||
1073 | if (fn_ge || fn_ene) | |||
1074 | { | |||
1075 | et_val = read_xvg_time(fn_ge ? fn_ge : fn_ene, bHaveT, | |||
1076 | opt2parg_bSet("-b", npargs, pa), tb-ttol, | |||
1077 | opt2parg_bSet("-e", npargs, pa), te+ttol, | |||
1078 | 1, &e_nset, &e_n, &e_dt, &e_t); | |||
1079 | if (fn_ge) | |||
1080 | { | |||
1081 | if (e_nset != 1) | |||
1082 | { | |||
1083 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1083, "Can only handle one free energy component in %s", | |||
1084 | fn_ge); | |||
1085 | } | |||
1086 | } | |||
1087 | else | |||
1088 | { | |||
1089 | if (e_nset != 1 && e_nset != 2) | |||
1090 | { | |||
1091 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1091, "Can only handle one energy component or one energy and one T in %s", | |||
1092 | fn_ene); | |||
1093 | } | |||
1094 | } | |||
1095 | if (e_n != n) | |||
1096 | { | |||
1097 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1097, "Number of energies (%d) does not match number of entries (%d) in %s", e_n, n, opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); | |||
1098 | } | |||
1099 | } | |||
1100 | else | |||
1101 | { | |||
1102 | et_val = NULL((void*)0); | |||
1103 | } | |||
1104 | ||||
1105 | if (opt2fn_null("-mdata", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) != NULL((void*)0)) | |||
1106 | { | |||
1107 | dt_val = read_xvg_time(opt2fn("-mdata", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bHaveT, | |||
1108 | FALSE0, tb, FALSE0, te, | |||
1109 | nsets_in, &d_nset, &d_n, &d_dt, &d_t); | |||
1110 | if (d_nset != 1) | |||
1111 | { | |||
1112 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1112, "Can only handle one mapping data column in %s", | |||
1113 | opt2fn("-mdata", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); | |||
1114 | } | |||
1115 | } | |||
1116 | else | |||
1117 | { | |||
1118 | dt_val = NULL((void*)0); | |||
1119 | } | |||
1120 | ||||
1121 | if (fn_ene && et_val) | |||
1122 | { | |||
1123 | ehisto(opt2fn("-histo", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), e_n, et_val, oenv); | |||
1124 | } | |||
1125 | ||||
1126 | snew(idim, max(3, nset))(idim) = save_calloc("idim", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1126, ((((3) > (nset)) ? (3) : (nset) )), sizeof(*(idim) )); | |||
1127 | snew(ibox, max(3, nset))(ibox) = save_calloc("ibox", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1127, ((((3) > (nset)) ? (3) : (nset) )), sizeof(*(ibox) )); | |||
1128 | snew(rmin, max(3, nset))(rmin) = save_calloc("rmin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1128, ((((3) > (nset)) ? (3) : (nset) )), sizeof(*(rmin) )); | |||
1129 | snew(rmax, max(3, nset))(rmax) = save_calloc("rmax", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1129, ((((3) > (nset)) ? (3) : (nset) )), sizeof(*(rmax) )); | |||
1130 | for (i = 0; (i < min(3, nset)(((3) < (nset)) ? (3) : (nset) )); i++) | |||
1131 | { | |||
1132 | idim[i] = nrdim[i]; | |||
1133 | ibox[i] = nrbox[i]; | |||
1134 | rmin[i] = xmin[i]; | |||
1135 | rmax[i] = xmax[i]; | |||
1136 | } | |||
1137 | for (; (i < nset); i++) | |||
1138 | { | |||
1139 | idim[i] = nrdim[2]; | |||
1140 | ibox[i] = nrbox[2]; | |||
1141 | rmin[i] = xmin[2]; | |||
1142 | rmax[i] = xmax[2]; | |||
1143 | } | |||
1144 | ||||
1145 | /* Check that the grid size is manageable. */ | |||
1146 | num_grid_points = ibox[0]; | |||
1147 | for (i = 1; i < nset; i++) | |||
1148 | { | |||
1149 | gmx_int64_t result; | |||
1150 | if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result)) | |||
1151 | { | |||
1152 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sham.c" , 1152, | |||
1153 | "The number of dimensions and grid points is too large for this tool.\n"); | |||
1154 | } | |||
1155 | num_grid_points = result; | |||
1156 | } | |||
1157 | /* The number of grid points fits in a gmx_int64_t. */ | |||
1158 | ||||
1159 | do_sham(opt2fn("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-bin", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
| ||||
1160 | opt2fn("-lp", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1161 | opt2fn("-ls", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-lsh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1162 | opt2fn("-lss", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-map", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1163 | opt2fn("-ls3", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-g", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
1164 | n, nset, val, fn_ge != NULL((void*)0), e_nset, et_val, d_n, d_t, dt_val, Tref, | |||
1165 | pmax, gmax, | |||
1166 | opt2parg_bSet("-emin", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa) ? &emin : NULL((void*)0), | |||
1167 | opt2parg_bSet("-emax", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa) ? &emax : NULL((void*)0), | |||
1168 | nlevels, pmin, | |||
1169 | mname, idim, ibox, | |||
1170 | opt2parg_bSet("-xmin", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa), rmin, | |||
1171 | opt2parg_bSet("-xmax", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa), rmax); | |||
1172 | ||||
1173 | return 0; | |||
1174 | } |