File: | gromacs/gmxana/gmx_hbond.c |
Location: | line 2446, column 17 |
Description: | Value stored to 'chi22' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2008, 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 | #include <math.h> |
41 | |
42 | /*#define HAVE_NN_LOOPS*/ |
43 | |
44 | #include "gromacs/commandline/pargs.h" |
45 | #include "copyrite.h" |
46 | #include "txtdump.h" |
47 | #include "physics.h" |
48 | #include "macros.h" |
49 | #include "gromacs/utility/fatalerror.h" |
50 | #include "index.h" |
51 | #include "gromacs/utility/smalloc.h" |
52 | #include "gromacs/math/vec.h" |
53 | #include "gromacs/fileio/xvgr.h" |
54 | #include "viewit.h" |
55 | #include "gstat.h" |
56 | #include "gromacs/utility/cstringutil.h" |
57 | #include "pbc.h" |
58 | #include "correl.h" |
59 | #include "gmx_ana.h" |
60 | #include "geminate.h" |
61 | |
62 | #include "gromacs/utility/futil.h" |
63 | #include "gromacs/fileio/matio.h" |
64 | #include "gromacs/fileio/tpxio.h" |
65 | #include "gromacs/fileio/trxio.h" |
66 | #include "gromacs/utility/gmxomp.h" |
67 | |
68 | typedef short int t_E; |
69 | typedef int t_EEst; |
70 | #define max_hx7 7 |
71 | typedef int t_hx[max_hx7]; |
72 | #define NRHXTYPES7 max_hx7 |
73 | const char *hxtypenames[NRHXTYPES7] = |
74 | {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"}; |
75 | #define MAXHH4 4 |
76 | |
77 | #ifdef GMX_OPENMP |
78 | #define MASTER_THREAD_ONLY(threadNr)((threadNr) == (threadNr)) ((threadNr) == 0) |
79 | #else |
80 | #define MASTER_THREAD_ONLY(threadNr)((threadNr) == (threadNr)) ((threadNr) == (threadNr)) |
81 | #endif |
82 | |
83 | /* -----------------------------------------*/ |
84 | |
85 | enum { |
86 | gr0, gr1, grI, grNR |
87 | }; |
88 | enum { |
89 | hbNo, hbDist, hbHB, hbNR, hbR2 |
90 | }; |
91 | enum { |
92 | noDA, ACC, DON, DA, INGROUP |
93 | }; |
94 | enum { |
95 | NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR |
96 | }; |
97 | |
98 | static const char *grpnames[grNR] = {"0", "1", "I" }; |
99 | |
100 | static gmx_bool bDebug = FALSE0; |
101 | |
102 | #define HB_NO0 0 |
103 | #define HB_YES1<<0 1<<0 |
104 | #define HB_INS1<<1 1<<1 |
105 | #define HB_YESINS1<<0|1<<1 HB_YES1<<0|HB_INS1<<1 |
106 | #define HB_NR(1<<2) (1<<2) |
107 | #define MAXHYDRO4 4 |
108 | |
109 | #define ISHB(h)(((h) & 2) == 2) (((h) & 2) == 2) |
110 | #define ISDIST(h)(((h) & 1) == 1) (((h) & 1) == 1) |
111 | #define ISDIST2(h)(((h) & 4) == 4) (((h) & 4) == 4) |
112 | #define ISACC(h)(((h) & 1) == 1) (((h) & 1) == 1) |
113 | #define ISDON(h)(((h) & 2) == 2) (((h) & 2) == 2) |
114 | #define ISINGRP(h)(((h) & 4) == 4) (((h) & 4) == 4) |
115 | |
116 | typedef struct { |
117 | int nr; |
118 | int maxnr; |
119 | atom_id *atoms; |
120 | } t_ncell; |
121 | |
122 | typedef struct { |
123 | t_ncell d[grNR]; |
124 | t_ncell a[grNR]; |
125 | } t_gridcell; |
126 | |
127 | typedef int t_icell[grNR]; |
128 | typedef atom_id h_id[MAXHYDRO4]; |
129 | |
130 | typedef struct { |
131 | int history[MAXHYDRO4]; |
132 | /* Has this hbond existed ever? If so as hbDist or hbHB or both. |
133 | * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB) |
134 | */ |
135 | /* Bitmask array which tells whether a hbond is present |
136 | * at a given time. Either of these may be NULL |
137 | */ |
138 | int n0; /* First frame a HB was found */ |
139 | int nframes, maxframes; /* Amount of frames in this hbond */ |
140 | unsigned int **h; |
141 | unsigned int **g; |
142 | /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the |
143 | * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor- |
144 | * acceptor distance is less than the user-specified distance (typically |
145 | * 0.35 nm). |
146 | */ |
147 | } t_hbond; |
148 | |
149 | typedef struct { |
150 | int nra, max_nra; |
151 | atom_id *acc; /* Atom numbers of the acceptors */ |
152 | int *grp; /* Group index */ |
153 | int *aptr; /* Map atom number to acceptor index */ |
154 | } t_acceptors; |
155 | |
156 | typedef struct { |
157 | int nrd, max_nrd; |
158 | int *don; /* Atom numbers of the donors */ |
159 | int *grp; /* Group index */ |
160 | int *dptr; /* Map atom number to donor index */ |
161 | int *nhydro; /* Number of hydrogens for each donor */ |
162 | h_id *hydro; /* The atom numbers of the hydrogens */ |
163 | h_id *nhbonds; /* The number of HBs per H at current */ |
164 | } t_donors; |
165 | |
166 | /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/ |
167 | #define PSTYPEint int |
168 | |
169 | typedef struct { |
170 | int len; /* The length of frame and p. */ |
171 | int *frame; /* The frames at which transitio*/ |
172 | PSTYPEint *p; |
173 | } t_pShift; |
174 | |
175 | typedef struct { |
176 | /* Periodicity history. Used for the reversible geminate recombination. */ |
177 | t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap: |
178 | * pHist[d][a]. We can safely assume that the same |
179 | * periodic shift holds for all hydrogens of a da-pair. |
180 | * |
181 | * Nowadays it only stores TRANSITIONS, and not the shift at every frame. |
182 | * That saves a LOT of memory, an hopefully kills a mysterious bug where |
183 | * pHist gets contaminated. */ |
184 | |
185 | PSTYPEint nper; /* The length of p2i */ |
186 | ivec *p2i; /* Maps integer to periodic shift for a pair.*/ |
187 | matrix P; /* Projection matrix to find the box shifts. */ |
188 | int gemtype; /* enumerated type */ |
189 | } t_gemPeriod; |
190 | |
191 | typedef struct { |
192 | int nframes; |
193 | int *Etot; /* Total energy for each frame */ |
194 | t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */ |
195 | } t_hbEmap; |
196 | |
197 | typedef struct { |
198 | gmx_bool bHBmap, bDAnr, bGem; |
199 | int wordlen; |
200 | /* The following arrays are nframes long */ |
201 | int nframes, max_frames, maxhydro; |
202 | int *nhb, *ndist; |
203 | h_id *n_bound; |
204 | real *time; |
205 | t_icell *danr; |
206 | t_hx *nhx; |
207 | /* These structures are initialized from the topology at start up */ |
208 | t_donors d; |
209 | t_acceptors a; |
210 | /* This holds a matrix with all possible hydrogen bonds */ |
211 | int nrhb, nrdist; |
212 | t_hbond ***hbmap; |
213 | #ifdef HAVE_NN_LOOPS |
214 | t_hbEmap hbE; |
215 | #endif |
216 | /* For parallelization reasons this will have to be a pointer. |
217 | * Otherwise discrepancies may arise between the periodicity data |
218 | * seen by different threads. */ |
219 | t_gemPeriod *per; |
220 | } t_hbdata; |
221 | |
222 | static void clearPshift(t_pShift *pShift) |
223 | { |
224 | if (pShift->len > 0) |
225 | { |
226 | sfree(pShift->p)save_free("pShift->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 226, (pShift->p)); |
227 | sfree(pShift->frame)save_free("pShift->frame", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 227, (pShift->frame)); |
228 | pShift->len = 0; |
229 | } |
230 | } |
231 | |
232 | static void calcBoxProjection(matrix B, matrix P) |
233 | { |
234 | const int vp[] = {XX0, YY1, ZZ2}; |
235 | int i, j; |
236 | int m, n; |
237 | matrix M, N, U; |
238 | |
239 | for (i = 0; i < 3; i++) |
240 | { |
241 | m = vp[i]; |
242 | for (j = 0; j < 3; j++) |
243 | { |
244 | n = vp[j]; |
245 | U[m][n] = i == j ? 1 : 0; |
246 | } |
247 | } |
248 | m_inv(B, M); |
249 | for (i = 0; i < 3; i++) |
250 | { |
251 | m = vp[i]; |
252 | mvmul(M, U[m], P[m]); |
253 | } |
254 | transpose(P, N); |
255 | } |
256 | |
257 | static void calcBoxDistance(matrix P, rvec d, ivec ibd) |
258 | { |
259 | /* returns integer distance in box coordinates. |
260 | * P is the projection matrix from cartesian coordinates |
261 | * obtained with calcBoxProjection(). */ |
262 | int i; |
263 | rvec bd; |
264 | mvmul(P, d, bd); |
265 | /* extend it by 0.5 in all directions since (int) rounds toward 0.*/ |
266 | for (i = 0; i < 3; i++) |
267 | { |
268 | bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5); |
269 | } |
270 | ibd[XX0] = (int)bd[XX0]; |
271 | ibd[YY1] = (int)bd[YY1]; |
272 | ibd[ZZ2] = (int)bd[ZZ2]; |
273 | } |
274 | |
275 | /* Changed argument 'bMerge' into 'oneHB' below, |
276 | * since -contact should cause maxhydro to be 1, |
277 | * not just -merge. |
278 | * - Erik Marklund May 29, 2006 |
279 | */ |
280 | |
281 | static PSTYPEint periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap) |
282 | { |
283 | /* Try to merge hbonds on the fly. That means that if the |
284 | * acceptor and donor are mergable, then: |
285 | * 1) store the hb-info so that acceptor id > donor id, |
286 | * 2) add the periodic shift in pairs, so that [-x,-y,-z] is |
287 | * stored in per.p2i[] whenever acceptor id < donor id. |
288 | * Note that [0,0,0] should already be the first element of per.p2i |
289 | * by the time this function is called. */ |
290 | |
291 | /* daSwap is TRUE if the donor and acceptor were swapped. |
292 | * If so, then the negative vector should be used. */ |
293 | PSTYPEint i; |
294 | |
295 | if (per->p2i == NULL((void*)0) || per->nper == 0) |
296 | { |
297 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 297, "'per' not initialized properly."); |
298 | } |
299 | for (i = 0; i < per->nper; i++) |
300 | { |
301 | if (r[XX0] == per->p2i[i][XX0] && |
302 | r[YY1] == per->p2i[i][YY1] && |
303 | r[ZZ2] == per->p2i[i][ZZ2]) |
304 | { |
305 | return i; |
306 | } |
307 | } |
308 | /* Not found apparently. Add it to the list! */ |
309 | /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */ |
310 | |
311 | #pragma omp critical |
312 | { |
313 | if (!per->p2i) |
314 | { |
315 | fprintf(stderrstderr, "p2i not initialized. This shouldn't happen!\n"); |
316 | snew(per->p2i, 1)(per->p2i) = save_calloc("per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 316, (1), sizeof(*(per->p2i))); |
317 | } |
318 | else |
319 | { |
320 | srenew(per->p2i, per->nper+2)(per->p2i) = save_realloc("per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 320, (per->p2i), (per->nper+2), sizeof(*(per->p2i) )); |
321 | } |
322 | copy_ivec(r, per->p2i[per->nper]); |
323 | (per->nper)++; |
324 | |
325 | /* Add the mirror too. It's rather likely that it'll be needed. */ |
326 | per->p2i[per->nper][XX0] = -r[XX0]; |
327 | per->p2i[per->nper][YY1] = -r[YY1]; |
328 | per->p2i[per->nper][ZZ2] = -r[ZZ2]; |
329 | (per->nper)++; |
330 | } /* omp critical */ |
331 | return per->nper - 1 - (daSwap ? 0 : 1); |
332 | } |
333 | |
334 | static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode) |
335 | { |
336 | t_hbdata *hb; |
337 | |
338 | snew(hb, 1)(hb) = save_calloc("hb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 338, (1), sizeof(*(hb))); |
339 | hb->wordlen = 8*sizeof(unsigned int); |
340 | hb->bHBmap = bHBmap; |
341 | hb->bDAnr = bDAnr; |
342 | hb->bGem = bGem; |
343 | if (oneHB) |
344 | { |
345 | hb->maxhydro = 1; |
346 | } |
347 | else |
348 | { |
349 | hb->maxhydro = MAXHYDRO4; |
350 | } |
351 | snew(hb->per, 1)(hb->per) = save_calloc("hb->per", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 351, (1), sizeof(*(hb->per))); |
352 | hb->per->gemtype = bGem ? gemmode : 0; |
353 | |
354 | return hb; |
355 | } |
356 | |
357 | static void mk_hbmap(t_hbdata *hb) |
358 | { |
359 | int i, j; |
360 | |
361 | snew(hb->hbmap, hb->d.nrd)(hb->hbmap) = save_calloc("hb->hbmap", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 361, (hb->d.nrd), sizeof(*(hb->hbmap))); |
362 | for (i = 0; (i < hb->d.nrd); i++) |
363 | { |
364 | snew(hb->hbmap[i], hb->a.nra)(hb->hbmap[i]) = save_calloc("hb->hbmap[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 364, (hb->a.nra), sizeof(*(hb->hbmap[i]))); |
365 | if (hb->hbmap[i] == NULL((void*)0)) |
366 | { |
367 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 367, "Could not allocate enough memory for hbmap"); |
368 | } |
369 | for (j = 0; (j > hb->a.nra); j++) |
370 | { |
371 | hb->hbmap[i][j] = NULL((void*)0); |
372 | } |
373 | } |
374 | } |
375 | |
376 | /* Consider redoing pHist so that is only stores transitions between |
377 | * periodicities and not the periodicity for all frames. This eats heaps of memory. */ |
378 | static void mk_per(t_hbdata *hb) |
379 | { |
380 | int i, j; |
381 | if (hb->bGem) |
382 | { |
383 | snew(hb->per->pHist, hb->d.nrd)(hb->per->pHist) = save_calloc("hb->per->pHist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 383, (hb->d.nrd), sizeof(*(hb->per->pHist))); |
384 | for (i = 0; i < hb->d.nrd; i++) |
385 | { |
386 | snew(hb->per->pHist[i], hb->a.nra)(hb->per->pHist[i]) = save_calloc("hb->per->pHist[i]" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 386, (hb->a.nra), sizeof(*(hb->per->pHist[i]))); |
387 | if (hb->per->pHist[i] == NULL((void*)0)) |
388 | { |
389 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 389, "Could not allocate enough memory for per->pHist"); |
390 | } |
391 | for (j = 0; j < hb->a.nra; j++) |
392 | { |
393 | clearPshift(&(hb->per->pHist[i][j])); |
394 | } |
395 | } |
396 | /* add the [0,0,0] shift to element 0 of p2i. */ |
397 | snew(hb->per->p2i, 1)(hb->per->p2i) = save_calloc("hb->per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 397, (1), sizeof(*(hb->per->p2i))); |
398 | clear_ivec(hb->per->p2i[0]); |
399 | hb->per->nper = 1; |
400 | } |
401 | } |
402 | |
403 | #ifdef HAVE_NN_LOOPS |
404 | static void mk_hbEmap (t_hbdata *hb, int n0) |
405 | { |
406 | int i, j, k; |
407 | hb->hbE.E = NULL((void*)0); |
408 | hb->hbE.nframes = 0; |
409 | snew(hb->hbE.E, hb->d.nrd)(hb->hbE.E) = save_calloc("hb->hbE.E", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 409, (hb->d.nrd), sizeof(*(hb->hbE.E))); |
410 | for (i = 0; i < hb->d.nrd; i++) |
411 | { |
412 | snew(hb->hbE.E[i], hb->a.nra)(hb->hbE.E[i]) = save_calloc("hb->hbE.E[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 412, (hb->a.nra), sizeof(*(hb->hbE.E[i]))); |
413 | for (j = 0; j < hb->a.nra; j++) |
414 | { |
415 | snew(hb->hbE.E[i][j], MAXHYDRO)(hb->hbE.E[i][j]) = save_calloc("hb->hbE.E[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 415, (4), sizeof(*(hb->hbE.E[i][j]))); |
416 | for (k = 0; k < MAXHYDRO4; k++) |
417 | { |
418 | hb->hbE.E[i][j][k] = NULL((void*)0); |
419 | } |
420 | } |
421 | } |
422 | hb->hbE.Etot = NULL((void*)0); |
423 | } |
424 | |
425 | static void free_hbEmap (t_hbdata *hb) |
426 | { |
427 | int i, j, k; |
428 | for (i = 0; i < hb->d.nrd; i++) |
429 | { |
430 | for (j = 0; j < hb->a.nra; j++) |
431 | { |
432 | for (k = 0; k < MAXHYDRO4; k++) |
433 | { |
434 | sfree(hb->hbE.E[i][j][k])save_free("hb->hbE.E[i][j][k]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 434, (hb->hbE.E[i][j][k])); |
435 | } |
436 | sfree(hb->hbE.E[i][j])save_free("hb->hbE.E[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 436, (hb->hbE.E[i][j])); |
437 | } |
438 | sfree(hb->hbE.E[i])save_free("hb->hbE.E[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 438, (hb->hbE.E[i])); |
439 | } |
440 | sfree(hb->hbE.E)save_free("hb->hbE.E", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 440, (hb->hbE.E)); |
441 | sfree(hb->hbE.Etot)save_free("hb->hbE.Etot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 441, (hb->hbE.Etot)); |
442 | } |
443 | |
444 | static void addFramesNN(t_hbdata *hb, int frame) |
445 | { |
446 | |
447 | #define DELTAFRAMES_HBE 10 |
448 | |
449 | int d, a, h, nframes; |
450 | |
451 | if (frame >= hb->hbE.nframes) |
452 | { |
453 | nframes = hb->hbE.nframes + DELTAFRAMES_HBE; |
454 | srenew(hb->hbE.Etot, nframes)(hb->hbE.Etot) = save_realloc("hb->hbE.Etot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 454, (hb->hbE.Etot), (nframes), sizeof(*(hb->hbE.Etot ))); |
455 | |
456 | for (d = 0; d < hb->d.nrd; d++) |
457 | { |
458 | for (a = 0; a < hb->a.nra; a++) |
459 | { |
460 | for (h = 0; h < hb->d.nhydro[d]; h++) |
461 | { |
462 | srenew(hb->hbE.E[d][a][h], nframes)(hb->hbE.E[d][a][h]) = save_realloc("hb->hbE.E[d][a][h]" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 462, (hb->hbE.E[d][a][h]), (nframes), sizeof(*(hb->hbE .E[d][a][h]))); |
463 | } |
464 | } |
465 | } |
466 | |
467 | hb->hbE.nframes += DELTAFRAMES_HBE; |
468 | } |
469 | } |
470 | |
471 | static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst, |
472 | matrix box, rvec hbox, t_donors *donors) |
473 | { |
474 | /* d - donor atom |
475 | * a - acceptor atom |
476 | * h - hydrogen |
477 | * alpha - angle between dipoles |
478 | * x[] - atomic positions |
479 | * EEst - the type of energy estimate (see enum in hbplugin.h) |
480 | * box - the box vectors \ |
481 | * hbox - half box lengths _These two are only needed for the pbc correction |
482 | */ |
483 | |
484 | t_E E; |
485 | rvec dist; |
486 | rvec dipole[2], xmol[3], xmean[2]; |
487 | int i; |
488 | real r, realE; |
489 | |
490 | if (d == a) |
491 | { |
492 | /* Self-interaction */ |
493 | return NONSENSE_E; |
494 | } |
495 | |
496 | switch (EEst) |
497 | { |
498 | case NN_BINARY: |
499 | /* This is a simple binary existence function that sets E=1 whenever |
500 | * the distance between the oxygens is equal too or less than 0.35 nm. |
501 | */ |
502 | rvec_sub(x[d], x[a], dist); |
503 | pbc_correct_gem(dist, box, hbox); |
504 | if (norm(dist) <= 0.35) |
505 | { |
506 | E = 1; |
507 | } |
508 | else |
509 | { |
510 | E = 0; |
511 | } |
512 | break; |
513 | |
514 | case NN_1_over_r3: |
515 | /* Negative potential energy of a dipole. |
516 | * E = -cos(alpha) * 1/r^3 */ |
517 | |
518 | copy_rvec(x[d], xmol[0]); /* donor */ |
519 | copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ |
520 | copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ |
521 | |
522 | svmul(15.9994*(1/1.008), xmol[0], xmean[0]); |
523 | rvec_inc(xmean[0], xmol[1]); |
524 | rvec_inc(xmean[0], xmol[2]); |
525 | for (i = 0; i < 3; i++) |
526 | { |
527 | xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
528 | } |
529 | |
530 | /* Assumes that all acceptors are also donors. */ |
531 | copy_rvec(x[a], xmol[0]); /* acceptor */ |
532 | copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ |
533 | copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */ |
534 | |
535 | |
536 | svmul(15.9994*(1/1.008), xmol[0], xmean[1]); |
537 | rvec_inc(xmean[1], xmol[1]); |
538 | rvec_inc(xmean[1], xmol[2]); |
539 | for (i = 0; i < 3; i++) |
540 | { |
541 | xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
542 | } |
543 | |
544 | rvec_sub(xmean[0], xmean[1], dist); |
545 | pbc_correct_gem(dist, box, hbox); |
546 | r = norm(dist); |
547 | |
548 | realE = pow(r, -3.0); |
549 | E = (t_E)(SCALEFACTOR_E * realE); |
550 | break; |
551 | |
552 | case NN_dipole: |
553 | /* Negative potential energy of a (unpolarizable) dipole. |
554 | * E = -cos(alpha) * 1/r^3 */ |
555 | clear_rvec(dipole[1]); |
556 | clear_rvec(dipole[0]); |
557 | |
558 | copy_rvec(x[d], xmol[0]); /* donor */ |
559 | copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ |
560 | copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ |
561 | |
562 | rvec_inc(dipole[0], xmol[1]); |
563 | rvec_inc(dipole[0], xmol[2]); |
564 | for (i = 0; i < 3; i++) |
565 | { |
566 | dipole[0][i] *= 0.5; |
567 | } |
568 | rvec_dec(dipole[0], xmol[0]); |
569 | |
570 | svmul(15.9994*(1/1.008), xmol[0], xmean[0]); |
571 | rvec_inc(xmean[0], xmol[1]); |
572 | rvec_inc(xmean[0], xmol[2]); |
573 | for (i = 0; i < 3; i++) |
574 | { |
575 | xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
576 | } |
577 | |
578 | /* Assumes that all acceptors are also donors. */ |
579 | copy_rvec(x[a], xmol[0]); /* acceptor */ |
580 | copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ |
581 | copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */ |
582 | |
583 | |
584 | rvec_inc(dipole[1], xmol[1]); |
585 | rvec_inc(dipole[1], xmol[2]); |
586 | for (i = 0; i < 3; i++) |
587 | { |
588 | dipole[1][i] *= 0.5; |
589 | } |
590 | rvec_dec(dipole[1], xmol[0]); |
591 | |
592 | svmul(15.9994*(1/1.008), xmol[0], xmean[1]); |
593 | rvec_inc(xmean[1], xmol[1]); |
594 | rvec_inc(xmean[1], xmol[2]); |
595 | for (i = 0; i < 3; i++) |
596 | { |
597 | xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
598 | } |
599 | |
600 | rvec_sub(xmean[0], xmean[1], dist); |
601 | pbc_correct_gem(dist, box, hbox); |
602 | r = norm(dist); |
603 | |
604 | double cosalpha = cos_angle(dipole[0], dipole[1]); |
605 | realE = cosalpha * pow(r, -3.0); |
606 | E = (t_E)(SCALEFACTOR_E * realE); |
607 | break; |
608 | |
609 | default: |
610 | printf("Can't do that type of energy estimate: %i\n.", EEst); |
611 | E = NONSENSE_E; |
612 | } |
613 | |
614 | return E; |
615 | } |
616 | |
617 | static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame) |
618 | { |
619 | /* hb - hbond data structure |
620 | d - donor |
621 | a - acceptor |
622 | h - hydrogen |
623 | E - estimate of the energy |
624 | frame - the current frame. |
625 | */ |
626 | |
627 | /* Store the estimated energy */ |
628 | if (E == NONSENSE_E) |
629 | { |
630 | E = 0; |
631 | } |
632 | |
633 | hb->hbE.E[d][a][h][frame] = E; |
634 | |
635 | #pragma omp critical |
636 | { |
637 | hb->hbE.Etot[frame] += E; |
638 | } |
639 | } |
640 | #endif /* HAVE_NN_LOOPS */ |
641 | |
642 | |
643 | /* Finds -v[] in the periodicity index */ |
644 | static int findMirror(PSTYPEint p, ivec v[], PSTYPEint nper) |
645 | { |
646 | PSTYPEint i; |
647 | ivec u; |
648 | for (i = 0; i < nper; i++) |
649 | { |
650 | if (v[i][XX0] == -(v[p][XX0]) && |
651 | v[i][YY1] == -(v[p][YY1]) && |
652 | v[i][ZZ2] == -(v[p][ZZ2])) |
653 | { |
654 | return (int)i; |
655 | } |
656 | } |
657 | printf("Couldn't find mirror of [%i, %i, %i], index \n", |
658 | v[p][XX0], |
659 | v[p][YY1], |
660 | v[p][ZZ2]); |
661 | return -1; |
662 | } |
663 | |
664 | |
665 | static void add_frames(t_hbdata *hb, int nframes) |
666 | { |
667 | int i, j, k, l; |
668 | |
669 | if (nframes >= hb->max_frames) |
670 | { |
671 | hb->max_frames += 4096; |
672 | srenew(hb->time, hb->max_frames)(hb->time) = save_realloc("hb->time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 672, (hb->time), (hb->max_frames), sizeof(*(hb->time ))); |
673 | srenew(hb->nhb, hb->max_frames)(hb->nhb) = save_realloc("hb->nhb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 673, (hb->nhb), (hb->max_frames), sizeof(*(hb->nhb ))); |
674 | srenew(hb->ndist, hb->max_frames)(hb->ndist) = save_realloc("hb->ndist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 674, (hb->ndist), (hb->max_frames), sizeof(*(hb->ndist ))); |
675 | srenew(hb->n_bound, hb->max_frames)(hb->n_bound) = save_realloc("hb->n_bound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 675, (hb->n_bound), (hb->max_frames), sizeof(*(hb-> n_bound))); |
676 | srenew(hb->nhx, hb->max_frames)(hb->nhx) = save_realloc("hb->nhx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 676, (hb->nhx), (hb->max_frames), sizeof(*(hb->nhx ))); |
677 | if (hb->bDAnr) |
678 | { |
679 | srenew(hb->danr, hb->max_frames)(hb->danr) = save_realloc("hb->danr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 679, (hb->danr), (hb->max_frames), sizeof(*(hb->danr ))); |
680 | } |
681 | } |
682 | hb->nframes = nframes; |
683 | } |
684 | |
685 | #define OFFSET(frame)(frame / 32) (frame / 32) |
686 | #define MASK(frame)(1 << (frame % 32)) (1 << (frame % 32)) |
687 | |
688 | static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue) |
689 | { |
690 | if (bValue) |
691 | { |
692 | hbexist[OFFSET(frame)(frame / 32)] |= MASK(frame)(1 << (frame % 32)); |
693 | } |
694 | else |
695 | { |
696 | hbexist[OFFSET(frame)(frame / 32)] &= ~MASK(frame)(1 << (frame % 32)); |
697 | } |
698 | } |
699 | |
700 | static gmx_bool is_hb(unsigned int hbexist[], int frame) |
701 | { |
702 | return ((hbexist[OFFSET(frame)(frame / 32)] & MASK(frame)(1 << (frame % 32))) != 0) ? 1 : 0; |
703 | } |
704 | |
705 | static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb) |
706 | { |
707 | unsigned int *ghptr = NULL((void*)0); |
708 | |
709 | if (ihb == hbHB) |
710 | { |
711 | ghptr = hb->hbmap[id][ia]->h[ih]; |
712 | } |
713 | else if (ihb == hbDist) |
714 | { |
715 | ghptr = hb->hbmap[id][ia]->g[ih]; |
716 | } |
717 | else |
718 | { |
719 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 719, "Incomprehensible iValue %d in set_hb", ihb); |
720 | } |
721 | |
722 | _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE1); |
723 | } |
724 | |
725 | static void addPshift(t_pShift *pHist, PSTYPEint p, int frame) |
726 | { |
727 | if (pHist->len == 0) |
728 | { |
729 | snew(pHist->frame, 1)(pHist->frame) = save_calloc("pHist->frame", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 729, (1), sizeof(*(pHist->frame))); |
730 | snew(pHist->p, 1)(pHist->p) = save_calloc("pHist->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 730, (1), sizeof(*(pHist->p))); |
731 | pHist->len = 1; |
732 | pHist->frame[0] = frame; |
733 | pHist->p[0] = p; |
734 | return; |
735 | } |
736 | else |
737 | if (pHist->p[pHist->len-1] != p) |
738 | { |
739 | pHist->len++; |
740 | srenew(pHist->frame, pHist->len)(pHist->frame) = save_realloc("pHist->frame", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 740, (pHist->frame), (pHist->len), sizeof(*(pHist-> frame))); |
741 | srenew(pHist->p, pHist->len)(pHist->p) = save_realloc("pHist->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 741, (pHist->p), (pHist->len), sizeof(*(pHist->p)) ); |
742 | pHist->frame[pHist->len-1] = frame; |
743 | pHist->p[pHist->len-1] = p; |
744 | } /* Otherwise, there is no transition. */ |
745 | return; |
746 | } |
747 | |
748 | static PSTYPEint getPshift(t_pShift pHist, int frame) |
749 | { |
750 | int f, i; |
751 | |
752 | if (pHist.len == 0 |
753 | || (pHist.len > 0 && pHist.frame[0] > frame)) |
754 | { |
755 | return -1; |
756 | } |
757 | |
758 | for (i = 0; i < pHist.len; i++) |
759 | { |
760 | f = pHist.frame[i]; |
761 | if (f == frame) |
762 | { |
763 | return pHist.p[i]; |
764 | } |
765 | if (f > frame) |
766 | { |
767 | return pHist.p[i-1]; |
768 | } |
769 | } |
770 | |
771 | /* It seems that frame is after the last periodic transition. Return the last periodicity. */ |
772 | return pHist.p[pHist.len-1]; |
773 | } |
774 | |
775 | static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPEint p) |
776 | { |
777 | int i, j, n; |
778 | t_hbond *hb = hbd->hbmap[id][ia]; |
779 | int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id])(((hbd->maxhydro) < (hbd->d.nhydro[id])) ? (hbd-> maxhydro) : (hbd->d.nhydro[id]) ); |
780 | int wlen = hbd->wordlen; |
781 | int delta = 32*wlen; |
782 | gmx_bool bGem = hbd->bGem; |
783 | |
784 | if (!hb->h[0]) |
785 | { |
786 | hb->n0 = frame; |
787 | hb->maxframes = delta; |
788 | for (i = 0; (i < maxhydro); i++) |
789 | { |
790 | snew(hb->h[i], hb->maxframes/wlen)(hb->h[i]) = save_calloc("hb->h[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 790, (hb->maxframes/wlen), sizeof(*(hb->h[i]))); |
791 | snew(hb->g[i], hb->maxframes/wlen)(hb->g[i]) = save_calloc("hb->g[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 791, (hb->maxframes/wlen), sizeof(*(hb->g[i]))); |
792 | } |
793 | } |
794 | else |
795 | { |
796 | hb->nframes = frame-hb->n0; |
797 | /* We need a while loop here because hbonds may be returning |
798 | * after a long time. |
799 | */ |
800 | while (hb->nframes >= hb->maxframes) |
801 | { |
802 | n = hb->maxframes + delta; |
803 | for (i = 0; (i < maxhydro); i++) |
804 | { |
805 | srenew(hb->h[i], n/wlen)(hb->h[i]) = save_realloc("hb->h[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 805, (hb->h[i]), (n/wlen), sizeof(*(hb->h[i]))); |
806 | srenew(hb->g[i], n/wlen)(hb->g[i]) = save_realloc("hb->g[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 806, (hb->g[i]), (n/wlen), sizeof(*(hb->g[i]))); |
807 | for (j = hb->maxframes/wlen; (j < n/wlen); j++) |
808 | { |
809 | hb->h[i][j] = 0; |
810 | hb->g[i][j] = 0; |
811 | } |
812 | } |
813 | |
814 | hb->maxframes = n; |
815 | } |
816 | } |
817 | if (frame >= 0) |
818 | { |
819 | set_hb(hbd, id, h, ia, frame, ihb); |
820 | if (bGem) |
821 | { |
822 | if (p >= hbd->per->nper) |
823 | { |
824 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 824, "invalid shift: p=%u, nper=%u", p, hbd->per->nper); |
825 | } |
826 | else |
827 | { |
828 | addPshift(&(hbd->per->pHist[id][ia]), p, frame); |
829 | } |
830 | |
831 | } |
832 | } |
833 | |
834 | } |
835 | |
836 | static void inc_nhbonds(t_donors *ddd, int d, int h) |
837 | { |
838 | int j; |
839 | int dptr = ddd->dptr[d]; |
840 | |
841 | for (j = 0; (j < ddd->nhydro[dptr]); j++) |
842 | { |
843 | if (ddd->hydro[dptr][j] == h) |
844 | { |
845 | ddd->nhbonds[dptr][j]++; |
846 | break; |
847 | } |
848 | } |
849 | if (j == ddd->nhydro[dptr]) |
850 | { |
851 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 851, "No such hydrogen %d on donor %d\n", h+1, d+1); |
852 | } |
853 | } |
854 | |
855 | static int _acceptor_index(t_acceptors *a, int grp, atom_id i, |
856 | const char *file, int line) |
857 | { |
858 | int ai = a->aptr[i]; |
859 | |
860 | if (a->grp[ai] != grp) |
861 | { |
862 | if (debug && bDebug) |
863 | { |
864 | fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n", |
865 | ai, a->grp[ai], grp, file, line); |
866 | } |
867 | return NOTSET-12345; |
868 | } |
869 | else |
870 | { |
871 | return ai; |
872 | } |
873 | } |
874 | #define acceptor_index(a, grp, i)_acceptor_index(a, grp, i, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 874) _acceptor_index(a, grp, i, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c", __LINE__874) |
875 | |
876 | static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line) |
877 | { |
878 | int di = d->dptr[i]; |
879 | |
880 | if (di == NOTSET-12345) |
881 | { |
882 | return NOTSET-12345; |
883 | } |
884 | |
885 | if (d->grp[di] != grp) |
886 | { |
887 | if (debug && bDebug) |
888 | { |
889 | fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n", |
890 | di, d->grp[di], grp, file, line); |
891 | } |
892 | return NOTSET-12345; |
893 | } |
894 | else |
895 | { |
896 | return di; |
897 | } |
898 | } |
899 | #define donor_index(d, grp, i)_donor_index(d, grp, i, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 899) _donor_index(d, grp, i, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c", __LINE__899) |
900 | |
901 | static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd) |
902 | { |
903 | /* g_hbond doesn't allow overlapping groups */ |
904 | if (grpa != grpd) |
905 | { |
906 | return FALSE0; |
907 | } |
908 | return |
909 | donor_index(&hb->d, grpd, a)_donor_index(&hb->d, grpd, a, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 909) != NOTSET-12345 |
910 | && acceptor_index(&hb->a, grpa, d)_acceptor_index(&hb->a, grpa, d, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 910) != NOTSET-12345; |
911 | } |
912 | |
913 | |
914 | static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa, |
915 | int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPEint p) |
916 | { |
917 | int k, id, ia, hh; |
918 | gmx_bool daSwap = FALSE0; |
919 | |
920 | if ((id = hb->d.dptr[d]) == NOTSET-12345) |
921 | { |
922 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 922, "No donor atom %d", d+1); |
923 | } |
924 | else if (grpd != hb->d.grp[id]) |
925 | { |
926 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 926, "Inconsistent donor groups, %d iso %d, atom %d", |
927 | grpd, hb->d.grp[id], d+1); |
928 | } |
929 | if ((ia = hb->a.aptr[a]) == NOTSET-12345) |
930 | { |
931 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 931, "No acceptor atom %d", a+1); |
932 | } |
933 | else if (grpa != hb->a.grp[ia]) |
934 | { |
935 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 935, "Inconsistent acceptor groups, %d iso %d, atom %d", |
936 | grpa, hb->a.grp[ia], a+1); |
937 | } |
938 | |
939 | if (bMerge) |
940 | { |
941 | |
942 | if (isInterchangable(hb, d, a, grpd, grpa) && d > a) |
943 | /* Then swap identity so that the id of d is lower then that of a. |
944 | * |
945 | * This should really be redundant by now, as is_hbond() now ought to return |
946 | * hbNo in the cases where this conditional is TRUE. */ |
947 | { |
948 | daSwap = TRUE1; |
949 | k = d; |
950 | d = a; |
951 | a = k; |
952 | |
953 | /* Now repeat donor/acc check. */ |
954 | if ((id = hb->d.dptr[d]) == NOTSET-12345) |
955 | { |
956 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 956, "No donor atom %d", d+1); |
957 | } |
958 | else if (grpd != hb->d.grp[id]) |
959 | { |
960 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 960, "Inconsistent donor groups, %d iso %d, atom %d", |
961 | grpd, hb->d.grp[id], d+1); |
962 | } |
963 | if ((ia = hb->a.aptr[a]) == NOTSET-12345) |
964 | { |
965 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 965, "No acceptor atom %d", a+1); |
966 | } |
967 | else if (grpa != hb->a.grp[ia]) |
968 | { |
969 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 969, "Inconsistent acceptor groups, %d iso %d, atom %d", |
970 | grpa, hb->a.grp[ia], a+1); |
971 | } |
972 | } |
973 | } |
974 | |
975 | if (hb->hbmap) |
976 | { |
977 | /* Loop over hydrogens to find which hydrogen is in this particular HB */ |
978 | if ((ihb == hbHB) && !bMerge && !bContact) |
979 | { |
980 | for (k = 0; (k < hb->d.nhydro[id]); k++) |
981 | { |
982 | if (hb->d.hydro[id][k] == h) |
983 | { |
984 | break; |
985 | } |
986 | } |
987 | if (k == hb->d.nhydro[id]) |
988 | { |
989 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 989, "Donor %d does not have hydrogen %d (a = %d)", |
990 | d+1, h+1, a+1); |
991 | } |
992 | } |
993 | else |
994 | { |
995 | k = 0; |
996 | } |
997 | |
998 | if (hb->bHBmap) |
999 | { |
1000 | |
1001 | #pragma omp critical |
1002 | { |
1003 | if (hb->hbmap[id][ia] == NULL((void*)0)) |
1004 | { |
1005 | snew(hb->hbmap[id][ia], 1)(hb->hbmap[id][ia]) = save_calloc("hb->hbmap[id][ia]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1005, (1), sizeof(*(hb->hbmap[id][ia]))); |
1006 | snew(hb->hbmap[id][ia]->h, hb->maxhydro)(hb->hbmap[id][ia]->h) = save_calloc("hb->hbmap[id][ia]->h" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1006, (hb->maxhydro), sizeof(*(hb->hbmap[id][ia]-> h))); |
1007 | snew(hb->hbmap[id][ia]->g, hb->maxhydro)(hb->hbmap[id][ia]->g) = save_calloc("hb->hbmap[id][ia]->g" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1007, (hb->maxhydro), sizeof(*(hb->hbmap[id][ia]-> g))); |
1008 | } |
1009 | add_ff(hb, id, k, ia, frame, ihb, p); |
1010 | } |
1011 | } |
1012 | |
1013 | /* Strange construction with frame >=0 is a relic from old code |
1014 | * for selected hbond analysis. It may be necessary again if that |
1015 | * is made to work again. |
1016 | */ |
1017 | if (frame >= 0) |
1018 | { |
1019 | hh = hb->hbmap[id][ia]->history[k]; |
1020 | if (ihb == hbHB) |
1021 | { |
1022 | hb->nhb[frame]++; |
1023 | if (!(ISHB(hh)(((hh) & 2) == 2))) |
1024 | { |
1025 | hb->hbmap[id][ia]->history[k] = hh | 2; |
1026 | hb->nrhb++; |
1027 | } |
1028 | } |
1029 | else |
1030 | { |
1031 | if (ihb == hbDist) |
1032 | { |
1033 | hb->ndist[frame]++; |
1034 | if (!(ISDIST(hh)(((hh) & 1) == 1))) |
1035 | { |
1036 | hb->hbmap[id][ia]->history[k] = hh | 1; |
1037 | hb->nrdist++; |
1038 | } |
1039 | } |
1040 | } |
1041 | } |
1042 | } |
1043 | else |
1044 | { |
1045 | if (frame >= 0) |
1046 | { |
1047 | if (ihb == hbHB) |
1048 | { |
1049 | hb->nhb[frame]++; |
1050 | } |
1051 | else |
1052 | { |
1053 | if (ihb == hbDist) |
1054 | { |
1055 | hb->ndist[frame]++; |
1056 | } |
1057 | } |
1058 | } |
1059 | } |
1060 | if (bMerge && daSwap) |
1061 | { |
1062 | h = hb->d.hydro[id][0]; |
1063 | } |
1064 | /* Increment number if HBonds per H */ |
1065 | if (ihb == hbHB && !bContact) |
1066 | { |
1067 | inc_nhbonds(&(hb->d), d, h); |
1068 | } |
1069 | } |
1070 | |
1071 | static char *mkatomname(t_atoms *atoms, int i) |
1072 | { |
1073 | static char buf[32]; |
1074 | int rnr; |
1075 | |
1076 | rnr = atoms->atom[i].resind; |
1077 | sprintf(buf, "%4s%d%-4s", |
1078 | *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]); |
1079 | |
1080 | return buf; |
1081 | } |
1082 | |
1083 | static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms) |
1084 | { |
1085 | /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */ |
1086 | int i; |
1087 | |
1088 | for (i = 0; i < isize; i++) |
1089 | { |
1090 | if (index[i] >= natoms) |
1091 | { |
1092 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1092, "Atom has index %d larger than number of atoms %d.", index[i], natoms); |
1093 | } |
1094 | datable[index[i]] |= INGROUP; |
1095 | } |
1096 | } |
1097 | |
1098 | static void clear_datable_grp(unsigned char *datable, int size) |
1099 | { |
1100 | /* Clears group information from the table */ |
1101 | int i; |
1102 | const char mask = !(char)INGROUP; |
1103 | if (size > 0) |
1104 | { |
1105 | for (i = 0; i < size; i++) |
1106 | { |
1107 | datable[i] &= mask; |
1108 | } |
1109 | } |
1110 | } |
1111 | |
1112 | static void add_acc(t_acceptors *a, int ia, int grp) |
1113 | { |
1114 | if (a->nra >= a->max_nra) |
1115 | { |
1116 | a->max_nra += 16; |
1117 | srenew(a->acc, a->max_nra)(a->acc) = save_realloc("a->acc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1117, (a->acc), (a->max_nra), sizeof(*(a->acc))); |
1118 | srenew(a->grp, a->max_nra)(a->grp) = save_realloc("a->grp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1118, (a->grp), (a->max_nra), sizeof(*(a->grp))); |
1119 | } |
1120 | a->grp[a->nra] = grp; |
1121 | a->acc[a->nra++] = ia; |
1122 | } |
1123 | |
1124 | static void search_acceptors(t_topology *top, int isize, |
1125 | atom_id *index, t_acceptors *a, int grp, |
1126 | gmx_bool bNitAcc, |
1127 | gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable) |
1128 | { |
1129 | int i, n; |
1130 | |
1131 | if (bDoIt) |
1132 | { |
1133 | for (i = 0; (i < isize); i++) |
1134 | { |
1135 | n = index[i]; |
1136 | if ((bContact || |
1137 | (((*top->atoms.atomname[n])[0] == 'O') || |
1138 | (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) && |
1139 | ISINGRP(datable[n])(((datable[n]) & 4) == 4)) |
1140 | { |
1141 | datable[n] |= ACC; /* set the atom's acceptor flag in datable. */ |
1142 | add_acc(a, n, grp); |
1143 | } |
1144 | } |
1145 | } |
1146 | snew(a->aptr, top->atoms.nr)(a->aptr) = save_calloc("a->aptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1146, (top->atoms.nr), sizeof(*(a->aptr))); |
1147 | for (i = 0; (i < top->atoms.nr); i++) |
1148 | { |
1149 | a->aptr[i] = NOTSET-12345; |
1150 | } |
1151 | for (i = 0; (i < a->nra); i++) |
1152 | { |
1153 | a->aptr[a->acc[i]] = i; |
1154 | } |
1155 | } |
1156 | |
1157 | static void add_h2d(int id, int ih, t_donors *ddd) |
1158 | { |
1159 | int i; |
1160 | |
1161 | for (i = 0; (i < ddd->nhydro[id]); i++) |
1162 | { |
1163 | if (ddd->hydro[id][i] == ih) |
1164 | { |
1165 | printf("Hm. This isn't the first time I found this donor (%d,%d)\n", |
1166 | ddd->don[id], ih); |
1167 | break; |
1168 | } |
1169 | } |
1170 | if (i == ddd->nhydro[id]) |
1171 | { |
1172 | if (ddd->nhydro[id] >= MAXHYDRO4) |
1173 | { |
1174 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1174, "Donor %d has more than %d hydrogens!", |
1175 | ddd->don[id], MAXHYDRO4); |
1176 | } |
1177 | ddd->hydro[id][i] = ih; |
1178 | ddd->nhydro[id]++; |
1179 | } |
1180 | } |
1181 | |
1182 | static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable) |
1183 | { |
1184 | int i; |
1185 | |
1186 | if (ISDON(datable[id])(((datable[id]) & 2) == 2) || !datable) |
1187 | { |
1188 | if (ddd->dptr[id] == NOTSET-12345) /* New donor */ |
1189 | { |
1190 | i = ddd->nrd; |
1191 | ddd->dptr[id] = i; |
1192 | } |
1193 | else |
1194 | { |
1195 | i = ddd->dptr[id]; |
1196 | } |
1197 | |
1198 | if (i == ddd->nrd) |
1199 | { |
1200 | if (ddd->nrd >= ddd->max_nrd) |
1201 | { |
1202 | ddd->max_nrd += 128; |
1203 | srenew(ddd->don, ddd->max_nrd)(ddd->don) = save_realloc("ddd->don", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1203, (ddd->don), (ddd->max_nrd), sizeof(*(ddd->don ))); |
1204 | srenew(ddd->nhydro, ddd->max_nrd)(ddd->nhydro) = save_realloc("ddd->nhydro", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1204, (ddd->nhydro), (ddd->max_nrd), sizeof(*(ddd-> nhydro))); |
1205 | srenew(ddd->hydro, ddd->max_nrd)(ddd->hydro) = save_realloc("ddd->hydro", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1205, (ddd->hydro), (ddd->max_nrd), sizeof(*(ddd-> hydro))); |
1206 | srenew(ddd->nhbonds, ddd->max_nrd)(ddd->nhbonds) = save_realloc("ddd->nhbonds", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1206, (ddd->nhbonds), (ddd->max_nrd), sizeof(*(ddd-> nhbonds))); |
1207 | srenew(ddd->grp, ddd->max_nrd)(ddd->grp) = save_realloc("ddd->grp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1207, (ddd->grp), (ddd->max_nrd), sizeof(*(ddd->grp ))); |
1208 | } |
1209 | ddd->don[ddd->nrd] = id; |
1210 | ddd->nhydro[ddd->nrd] = 0; |
1211 | ddd->grp[ddd->nrd] = grp; |
1212 | ddd->nrd++; |
1213 | } |
1214 | else |
1215 | { |
1216 | ddd->don[i] = id; |
1217 | } |
1218 | add_h2d(i, ih, ddd); |
1219 | } |
1220 | else |
1221 | if (datable) |
1222 | { |
1223 | printf("Warning: Atom %d is not in the d/a-table!\n", id); |
1224 | } |
1225 | } |
1226 | |
1227 | static void search_donors(t_topology *top, int isize, atom_id *index, |
1228 | t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt, |
1229 | unsigned char *datable) |
1230 | { |
1231 | int i, j, nra, n; |
1232 | t_functype func_type; |
1233 | t_ilist *interaction; |
1234 | atom_id nr1, nr2, nr3; |
1235 | gmx_bool stop; |
1236 | |
1237 | if (!ddd->dptr) |
1238 | { |
1239 | snew(ddd->dptr, top->atoms.nr)(ddd->dptr) = save_calloc("ddd->dptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1239, (top->atoms.nr), sizeof(*(ddd->dptr))); |
1240 | for (i = 0; (i < top->atoms.nr); i++) |
1241 | { |
1242 | ddd->dptr[i] = NOTSET-12345; |
1243 | } |
1244 | } |
1245 | |
1246 | if (bContact) |
1247 | { |
1248 | if (bDoIt) |
1249 | { |
1250 | for (i = 0; (i < isize); i++) |
1251 | { |
1252 | datable[index[i]] |= DON; |
1253 | add_dh(ddd, index[i], -1, grp, datable); |
1254 | } |
1255 | } |
1256 | } |
1257 | else |
1258 | { |
1259 | for (func_type = 0; (func_type < F_NRE); func_type++) |
1260 | { |
1261 | interaction = &(top->idef.il[func_type]); |
1262 | if (func_type == F_POSRES || func_type == F_FBPOSRES) |
1263 | { |
1264 | /* The ilist looks strange for posre. Bug in grompp? |
1265 | * We don't need posre interactions for hbonds anyway.*/ |
1266 | continue; |
1267 | } |
1268 | for (i = 0; i < interaction->nr; |
1269 | i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) |
1270 | { |
1271 | /* next function */ |
1272 | if (func_type != top->idef.functype[interaction->iatoms[i]]) |
1273 | { |
1274 | fprintf(stderrstderr, "Error in func_type %s", |
1275 | interaction_function[func_type].longname); |
1276 | continue; |
1277 | } |
1278 | |
1279 | /* check out this functype */ |
1280 | if (func_type == F_SETTLE) |
1281 | { |
1282 | nr1 = interaction->iatoms[i+1]; |
1283 | nr2 = interaction->iatoms[i+2]; |
1284 | nr3 = interaction->iatoms[i+3]; |
1285 | |
1286 | if (ISINGRP(datable[nr1])(((datable[nr1]) & 4) == 4)) |
1287 | { |
1288 | if (ISINGRP(datable[nr2])(((datable[nr2]) & 4) == 4)) |
1289 | { |
1290 | datable[nr1] |= DON; |
1291 | add_dh(ddd, nr1, nr1+1, grp, datable); |
1292 | } |
1293 | if (ISINGRP(datable[nr3])(((datable[nr3]) & 4) == 4)) |
1294 | { |
1295 | datable[nr1] |= DON; |
1296 | add_dh(ddd, nr1, nr1+2, grp, datable); |
1297 | } |
1298 | } |
1299 | } |
1300 | else if (IS_CHEMBOND(func_type)(interaction_function[(func_type)].nratoms == 2 && (interaction_function [(func_type)].flags & 1<<3))) |
1301 | { |
1302 | for (j = 0; j < 2; j++) |
1303 | { |
1304 | nr1 = interaction->iatoms[i+1+j]; |
1305 | nr2 = interaction->iatoms[i+2-j]; |
1306 | if ((*top->atoms.atomname[nr1][0] == 'H') && |
1307 | ((*top->atoms.atomname[nr2][0] == 'O') || |
1308 | (*top->atoms.atomname[nr2][0] == 'N')) && |
1309 | ISINGRP(datable[nr1])(((datable[nr1]) & 4) == 4) && ISINGRP(datable[nr2])(((datable[nr2]) & 4) == 4)) |
1310 | { |
1311 | datable[nr2] |= DON; |
1312 | add_dh(ddd, nr2, nr1, grp, datable); |
1313 | } |
1314 | } |
1315 | } |
1316 | } |
1317 | } |
1318 | #ifdef SAFEVSITES |
1319 | for (func_type = 0; func_type < F_NRE; func_type++) |
1320 | { |
1321 | interaction = &top->idef.il[func_type]; |
1322 | for (i = 0; i < interaction->nr; |
1323 | i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) |
1324 | { |
1325 | /* next function */ |
1326 | if (func_type != top->idef.functype[interaction->iatoms[i]]) |
1327 | { |
1328 | gmx_incons("function type in search_donors")_gmx_error("incons", "function type in search_donors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1328); |
1329 | } |
1330 | |
1331 | if (interaction_function[func_type].flags & IF_VSITE1<<1) |
1332 | { |
1333 | nr1 = interaction->iatoms[i+1]; |
1334 | if (*top->atoms.atomname[nr1][0] == 'H') |
1335 | { |
1336 | nr2 = nr1-1; |
1337 | stop = FALSE0; |
1338 | while (!stop && ( *top->atoms.atomname[nr2][0] == 'H')) |
1339 | { |
1340 | if (nr2) |
1341 | { |
1342 | nr2--; |
1343 | } |
1344 | else |
1345 | { |
1346 | stop = TRUE1; |
1347 | } |
1348 | } |
1349 | if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') || |
1350 | ( *top->atoms.atomname[nr2][0] == 'N') ) && |
1351 | ISINGRP(datable[nr1])(((datable[nr1]) & 4) == 4) && ISINGRP(datable[nr2])(((datable[nr2]) & 4) == 4)) |
1352 | { |
1353 | datable[nr2] |= DON; |
1354 | add_dh(ddd, nr2, nr1, grp, datable); |
1355 | } |
1356 | } |
1357 | } |
1358 | } |
1359 | } |
1360 | #endif |
1361 | } |
1362 | } |
1363 | |
1364 | static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid) |
1365 | { |
1366 | t_gridcell ***grid; |
1367 | int i, y, z; |
1368 | |
1369 | if (bBox) |
1370 | { |
1371 | for (i = 0; i < DIM3; i++) |
1372 | { |
1373 | ngrid[i] = (box[i][i]/(1.2*rcut)); |
1374 | } |
1375 | } |
1376 | |
1377 | if (!bBox || (ngrid[XX0] < 3) || (ngrid[YY1] < 3) || (ngrid[ZZ2] < 3) ) |
1378 | { |
1379 | for (i = 0; i < DIM3; i++) |
1380 | { |
1381 | ngrid[i] = 1; |
1382 | } |
1383 | } |
1384 | else |
1385 | { |
1386 | printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n", |
1387 | ngrid[XX0], ngrid[YY1], ngrid[ZZ2], rcut); |
1388 | } |
1389 | snew(grid, ngrid[ZZ])(grid) = save_calloc("grid", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1389, (ngrid[2]), sizeof(*(grid))); |
1390 | for (z = 0; z < ngrid[ZZ2]; z++) |
1391 | { |
1392 | snew((grid)[z], ngrid[YY])((grid)[z]) = save_calloc("(grid)[z]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1392, (ngrid[1]), sizeof(*((grid)[z]))); |
1393 | for (y = 0; y < ngrid[YY1]; y++) |
1394 | { |
1395 | snew((grid)[z][y], ngrid[XX])((grid)[z][y]) = save_calloc("(grid)[z][y]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1395, (ngrid[0]), sizeof(*((grid)[z][y]))); |
1396 | } |
1397 | } |
1398 | return grid; |
1399 | } |
1400 | |
1401 | static void reset_nhbonds(t_donors *ddd) |
1402 | { |
1403 | int i, j; |
1404 | |
1405 | for (i = 0; (i < ddd->nrd); i++) |
1406 | { |
1407 | for (j = 0; (j < MAXHH4); j++) |
1408 | { |
1409 | ddd->nhbonds[i][j] = 0; |
1410 | } |
1411 | } |
1412 | } |
1413 | |
1414 | void pbc_correct_gem(rvec dx, matrix box, rvec hbox); |
1415 | |
1416 | static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, |
1417 | gmx_bool bBox, matrix box, rvec hbox, |
1418 | real rcut, real rshell, |
1419 | ivec ngrid, t_gridcell ***grid) |
1420 | { |
1421 | int i, m, gr, xi, yi, zi, nr; |
1422 | atom_id *ad; |
1423 | ivec grididx; |
1424 | rvec invdelta, dshell, xtemp = {0, 0, 0}; |
1425 | t_ncell *newgrid; |
1426 | gmx_bool bDoRshell, bInShell, bAcc; |
1427 | real rshell2 = 0; |
1428 | int gx, gy, gz; |
1429 | int dum = -1; |
1430 | |
1431 | bDoRshell = (rshell > 0); |
1432 | rshell2 = sqr(rshell); |
1433 | bInShell = TRUE1; |
1434 | |
1435 | #define DBB(x)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1435,"x", x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__1435,#x, x) |
1436 | DBB(dum)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1436,"dum", dum); |
1437 | for (m = 0; m < DIM3; m++) |
1438 | { |
1439 | hbox[m] = box[m][m]*0.5; |
1440 | if (bBox) |
1441 | { |
1442 | invdelta[m] = ngrid[m]/box[m][m]; |
1443 | if (1/invdelta[m] < rcut) |
1444 | { |
1445 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1445, "Your computational box has shrunk too much.\n" |
1446 | "%s can not handle this situation, sorry.\n", |
1447 | ShortProgram()); |
1448 | } |
1449 | } |
1450 | else |
1451 | { |
1452 | invdelta[m] = 0; |
1453 | } |
1454 | } |
1455 | grididx[XX0] = 0; |
1456 | grididx[YY1] = 0; |
1457 | grididx[ZZ2] = 0; |
1458 | DBB(dum)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1458,"dum", dum); |
1459 | /* resetting atom counts */ |
1460 | for (gr = 0; (gr < grNR); gr++) |
1461 | { |
1462 | for (zi = 0; zi < ngrid[ZZ2]; zi++) |
1463 | { |
1464 | for (yi = 0; yi < ngrid[YY1]; yi++) |
1465 | { |
1466 | for (xi = 0; xi < ngrid[XX0]; xi++) |
1467 | { |
1468 | grid[zi][yi][xi].d[gr].nr = 0; |
1469 | grid[zi][yi][xi].a[gr].nr = 0; |
1470 | } |
1471 | } |
1472 | } |
1473 | DBB(dum)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1473,"dum", dum); |
1474 | |
1475 | /* put atoms in grid cells */ |
1476 | for (bAcc = FALSE0; (bAcc <= TRUE1); bAcc++) |
1477 | { |
1478 | if (bAcc) |
1479 | { |
1480 | nr = hb->a.nra; |
1481 | ad = hb->a.acc; |
1482 | } |
1483 | else |
1484 | { |
1485 | nr = hb->d.nrd; |
1486 | ad = hb->d.don; |
1487 | } |
1488 | DBB(bAcc)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1488,"bAcc", bAcc); |
1489 | for (i = 0; (i < nr); i++) |
1490 | { |
1491 | /* check if we are inside the shell */ |
1492 | /* if bDoRshell=FALSE then bInShell=TRUE always */ |
1493 | DBB(i)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1493,"i", i); |
1494 | if (bDoRshell) |
1495 | { |
1496 | bInShell = TRUE1; |
1497 | rvec_sub(x[ad[i]], xshell, dshell); |
1498 | if (bBox) |
1499 | { |
1500 | if (FALSE0 && !hb->bGem) |
1501 | { |
1502 | for (m = DIM3-1; m >= 0 && bInShell; m--) |
1503 | { |
1504 | if (dshell[m] < -hbox[m]) |
1505 | { |
1506 | rvec_inc(dshell, box[m]); |
1507 | } |
1508 | else if (dshell[m] >= hbox[m]) |
1509 | { |
1510 | dshell[m] -= 2*hbox[m]; |
1511 | } |
1512 | /* if we're outside the cube, we're outside the sphere also! */ |
1513 | if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) |
1514 | { |
1515 | bInShell = FALSE0; |
1516 | } |
1517 | } |
1518 | } |
1519 | else |
1520 | { |
1521 | gmx_bool bDone = FALSE0; |
1522 | while (!bDone) |
1523 | { |
1524 | bDone = TRUE1; |
1525 | for (m = DIM3-1; m >= 0 && bInShell; m--) |
1526 | { |
1527 | if (dshell[m] < -hbox[m]) |
1528 | { |
1529 | bDone = FALSE0; |
1530 | rvec_inc(dshell, box[m]); |
1531 | } |
1532 | if (dshell[m] >= hbox[m]) |
1533 | { |
1534 | bDone = FALSE0; |
1535 | dshell[m] -= 2*hbox[m]; |
1536 | } |
1537 | } |
1538 | } |
1539 | for (m = DIM3-1; m >= 0 && bInShell; m--) |
1540 | { |
1541 | /* if we're outside the cube, we're outside the sphere also! */ |
1542 | if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) |
1543 | { |
1544 | bInShell = FALSE0; |
1545 | } |
1546 | } |
1547 | } |
1548 | } |
1549 | /* if we're inside the cube, check if we're inside the sphere */ |
1550 | if (bInShell) |
1551 | { |
1552 | bInShell = norm2(dshell) < rshell2; |
1553 | } |
1554 | } |
1555 | DBB(i)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1555,"i", i); |
1556 | if (bInShell) |
1557 | { |
1558 | if (bBox) |
1559 | { |
1560 | if (hb->bGem) |
1561 | { |
1562 | copy_rvec(x[ad[i]], xtemp); |
1563 | } |
1564 | pbc_correct_gem(x[ad[i]], box, hbox); |
1565 | } |
1566 | for (m = DIM3-1; m >= 0; m--) |
1567 | { |
1568 | if (TRUE1 || !hb->bGem) |
1569 | { |
1570 | /* put atom in the box */ |
1571 | while (x[ad[i]][m] < 0) |
1572 | { |
1573 | rvec_inc(x[ad[i]], box[m]); |
1574 | } |
1575 | while (x[ad[i]][m] >= box[m][m]) |
1576 | { |
1577 | rvec_dec(x[ad[i]], box[m]); |
1578 | } |
1579 | } |
1580 | /* determine grid index of atom */ |
1581 | grididx[m] = x[ad[i]][m]*invdelta[m]; |
1582 | grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m]; |
1583 | } |
1584 | if (hb->bGem) |
1585 | { |
1586 | copy_rvec(xtemp, x[ad[i]]); /* copy back */ |
1587 | } |
1588 | gx = grididx[XX0]; |
1589 | gy = grididx[YY1]; |
1590 | gz = grididx[ZZ2]; |
1591 | range_check(gx, 0, ngrid[XX])_range_check(gx, 0, ngrid[0], ((void*)0),"gx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1591); |
1592 | range_check(gy, 0, ngrid[YY])_range_check(gy, 0, ngrid[1], ((void*)0),"gy", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1592); |
1593 | range_check(gz, 0, ngrid[ZZ])_range_check(gz, 0, ngrid[2], ((void*)0),"gz", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1593); |
1594 | DBB(gx)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1594,"gx", gx); |
1595 | DBB(gy)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1595,"gy", gy); |
1596 | DBB(gz)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1596,"gz", gz); |
1597 | /* add atom to grid cell */ |
1598 | if (bAcc) |
1599 | { |
1600 | newgrid = &(grid[gz][gy][gx].a[gr]); |
1601 | } |
1602 | else |
1603 | { |
1604 | newgrid = &(grid[gz][gy][gx].d[gr]); |
1605 | } |
1606 | if (newgrid->nr >= newgrid->maxnr) |
1607 | { |
1608 | newgrid->maxnr += 10; |
1609 | DBB(newgrid->maxnr)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1609,"newgrid->maxnr", newgrid->maxnr); |
1610 | srenew(newgrid->atoms, newgrid->maxnr)(newgrid->atoms) = save_realloc("newgrid->atoms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1610, (newgrid->atoms), (newgrid->maxnr), sizeof(*(newgrid ->atoms))); |
1611 | } |
1612 | DBB(newgrid->nr)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1612,"newgrid->nr", newgrid->nr); |
1613 | newgrid->atoms[newgrid->nr] = ad[i]; |
1614 | newgrid->nr++; |
1615 | } |
1616 | } |
1617 | } |
1618 | } |
1619 | } |
1620 | |
1621 | static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr) |
1622 | { |
1623 | int gr, xi, yi, zi; |
1624 | |
1625 | for (gr = 0; (gr < grNR); gr++) |
1626 | { |
1627 | danr[gr] = 0; |
1628 | for (zi = 0; zi < ngrid[ZZ2]; zi++) |
1629 | { |
1630 | for (yi = 0; yi < ngrid[YY1]; yi++) |
1631 | { |
1632 | for (xi = 0; xi < ngrid[XX0]; xi++) |
1633 | { |
1634 | danr[gr] += grid[zi][yi][xi].d[gr].nr; |
1635 | } |
1636 | } |
1637 | } |
1638 | } |
1639 | } |
1640 | |
1641 | /* The grid loop. |
1642 | * Without a box, the grid is 1x1x1, so all loops are 1 long. |
1643 | * With a rectangular box (bTric==FALSE) all loops are 3 long. |
1644 | * With a triclinic box all loops are 3 long, except when a cell is |
1645 | * located next to one of the box edges which is not parallel to the |
1646 | * x/y-plane, in that case all cells in a line or layer are searched. |
1647 | * This could be implemented slightly more efficient, but the code |
1648 | * would get much more complicated. |
1649 | */ |
1650 | static gmx_inlineinline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge) |
1651 | { |
1652 | return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1)); |
1653 | } |
1654 | static gmx_inlineinline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge) |
1655 | { |
1656 | return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1)); |
1657 | } |
1658 | static gmx_inlineinline int grid_mod(int j, int n) |
1659 | { |
1660 | return (j+n) % (n); |
1661 | } |
1662 | |
1663 | static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid) |
1664 | { |
1665 | int gr, x, y, z, sum[grNR]; |
1666 | |
1667 | fprintf(fp, "grid %dx%dx%d\n", ngrid[XX0], ngrid[YY1], ngrid[ZZ2]); |
1668 | for (gr = 0; gr < grNR; gr++) |
1669 | { |
1670 | sum[gr] = 0; |
1671 | fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]); |
1672 | for (z = 0; z < ngrid[ZZ2]; z += 2) |
1673 | { |
1674 | fprintf(fp, "Z=%d,%d\n", z, z+1); |
1675 | for (y = 0; y < ngrid[YY1]; y++) |
1676 | { |
1677 | for (x = 0; x < ngrid[XX0]; x++) |
1678 | { |
1679 | fprintf(fp, "%3d", grid[x][y][z].d[gr].nr); |
1680 | sum[gr] += grid[z][y][x].d[gr].nr; |
1681 | fprintf(fp, "%3d", grid[x][y][z].a[gr].nr); |
1682 | sum[gr] += grid[z][y][x].a[gr].nr; |
1683 | |
1684 | } |
1685 | fprintf(fp, " | "); |
1686 | if ( (z+1) < ngrid[ZZ2]) |
1687 | { |
1688 | for (x = 0; x < ngrid[XX0]; x++) |
1689 | { |
1690 | fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr); |
1691 | sum[gr] += grid[z+1][y][x].d[gr].nr; |
1692 | fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr); |
1693 | sum[gr] += grid[z+1][y][x].a[gr].nr; |
1694 | } |
1695 | } |
1696 | fprintf(fp, "\n"); |
1697 | } |
1698 | } |
1699 | } |
1700 | fprintf(fp, "TOTALS:"); |
1701 | for (gr = 0; gr < grNR; gr++) |
1702 | { |
1703 | fprintf(fp, " %d=%d", gr, sum[gr]); |
1704 | } |
1705 | fprintf(fp, "\n"); |
1706 | } |
1707 | |
1708 | /* New GMX record! 5 * in a row. Congratulations! |
1709 | * Sorry, only four left. |
1710 | */ |
1711 | static void free_grid(ivec ngrid, t_gridcell ****grid) |
1712 | { |
1713 | int y, z; |
1714 | t_gridcell ***g = *grid; |
1715 | |
1716 | for (z = 0; z < ngrid[ZZ2]; z++) |
1717 | { |
1718 | for (y = 0; y < ngrid[YY1]; y++) |
1719 | { |
1720 | sfree(g[z][y])save_free("g[z][y]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1720, (g[z][y])); |
1721 | } |
1722 | sfree(g[z])save_free("g[z]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1722, (g[z])); |
1723 | } |
1724 | sfree(g)save_free("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1724, (g)); |
1725 | g = NULL((void*)0); |
1726 | } |
1727 | |
1728 | void pbc_correct_gem(rvec dx, matrix box, rvec hbox) |
1729 | { |
1730 | int m; |
1731 | gmx_bool bDone = FALSE0; |
1732 | while (!bDone) |
1733 | { |
1734 | bDone = TRUE1; |
1735 | for (m = DIM3-1; m >= 0; m--) |
1736 | { |
1737 | if (dx[m] < -hbox[m]) |
1738 | { |
1739 | bDone = FALSE0; |
1740 | rvec_inc(dx, box[m]); |
1741 | } |
1742 | if (dx[m] >= hbox[m]) |
1743 | { |
1744 | bDone = FALSE0; |
1745 | rvec_dec(dx, box[m]); |
1746 | } |
1747 | } |
1748 | } |
1749 | } |
1750 | |
1751 | /* Added argument r2cut, changed contact and implemented |
1752 | * use of second cut-off. |
1753 | * - Erik Marklund, June 29, 2006 |
1754 | */ |
1755 | static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, |
1756 | real rcut, real r2cut, real ccut, |
1757 | rvec x[], gmx_bool bBox, matrix box, rvec hbox, |
1758 | real *d_ha, real *ang, gmx_bool bDA, int *hhh, |
1759 | gmx_bool bContact, gmx_bool bMerge, PSTYPEint *p) |
1760 | { |
1761 | int h, hh, id, ja, ihb; |
1762 | rvec r_da, r_ha, r_dh, r = {0, 0, 0}; |
1763 | ivec ri; |
1764 | real rc2, r2c2, rda2, rha2, ca; |
1765 | gmx_bool HAinrange = FALSE0; /* If !bDA. Needed for returning hbDist in a correct way. */ |
1766 | gmx_bool daSwap = FALSE0; |
1767 | |
1768 | if (d == a) |
1769 | { |
1770 | return hbNo; |
1771 | } |
1772 | |
1773 | if (((id = donor_index(&hb->d, grpd, d)_donor_index(&hb->d, grpd, d, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1773)) == NOTSET-12345) || |
1774 | ((ja = acceptor_index(&hb->a, grpa, a)_acceptor_index(&hb->a, grpa, a, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1774)) == NOTSET-12345)) |
1775 | { |
1776 | return hbNo; |
1777 | } |
1778 | |
1779 | rc2 = rcut*rcut; |
1780 | r2c2 = r2cut*r2cut; |
1781 | |
1782 | rvec_sub(x[d], x[a], r_da); |
1783 | /* Insert projection code here */ |
1784 | |
1785 | if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa)) |
1786 | { |
1787 | /* Then this hbond/contact will be found again, or it has already been found. */ |
1788 | /*return hbNo;*/ |
1789 | } |
1790 | if (bBox) |
1791 | { |
1792 | if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */ |
1793 | { /* return hbNo; */ |
1794 | daSwap = TRUE1; /* If so, then their history should be filed with donor and acceptor swapped. */ |
1795 | } |
1796 | if (hb->bGem) |
1797 | { |
1798 | copy_rvec(r_da, r); /* Save this for later */ |
1799 | pbc_correct_gem(r_da, box, hbox); |
1800 | } |
1801 | else |
1802 | { |
1803 | pbc_correct_gem(r_da, box, hbox); |
1804 | } |
1805 | } |
1806 | rda2 = iprod(r_da, r_da); |
1807 | |
1808 | if (bContact) |
1809 | { |
1810 | if (daSwap && grpa == grpd) |
1811 | { |
1812 | return hbNo; |
1813 | } |
1814 | if (rda2 <= rc2) |
1815 | { |
1816 | if (hb->bGem) |
1817 | { |
1818 | calcBoxDistance(hb->per->P, r, ri); |
1819 | *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */ |
1820 | } |
1821 | return hbHB; |
1822 | } |
1823 | else if (rda2 < r2c2) |
1824 | { |
1825 | return hbDist; |
1826 | } |
1827 | else |
1828 | { |
1829 | return hbNo; |
1830 | } |
1831 | } |
1832 | *hhh = NOTSET-12345; |
1833 | |
1834 | if (bDA && (rda2 > rc2)) |
1835 | { |
1836 | return hbNo; |
1837 | } |
1838 | |
1839 | for (h = 0; (h < hb->d.nhydro[id]); h++) |
1840 | { |
1841 | hh = hb->d.hydro[id][h]; |
1842 | rha2 = rc2+1; |
1843 | if (!bDA) |
1844 | { |
1845 | rvec_sub(x[hh], x[a], r_ha); |
1846 | if (bBox) |
1847 | { |
1848 | pbc_correct_gem(r_ha, box, hbox); |
1849 | } |
1850 | rha2 = iprod(r_ha, r_ha); |
1851 | } |
1852 | |
1853 | if (hb->bGem) |
1854 | { |
1855 | calcBoxDistance(hb->per->P, r, ri); |
1856 | *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */ |
1857 | } |
1858 | |
1859 | if (bDA || (!bDA && (rha2 <= rc2))) |
1860 | { |
1861 | rvec_sub(x[d], x[hh], r_dh); |
1862 | if (bBox) |
1863 | { |
1864 | pbc_correct_gem(r_dh, box, hbox); |
1865 | } |
1866 | |
1867 | if (!bDA) |
1868 | { |
1869 | HAinrange = TRUE1; |
1870 | } |
1871 | ca = cos_angle(r_dh, r_da); |
1872 | /* if angle is smaller, cos is larger */ |
1873 | if (ca >= ccut) |
1874 | { |
1875 | *hhh = hh; |
1876 | *d_ha = sqrt(bDA ? rda2 : rha2); |
1877 | *ang = acos(ca); |
1878 | return hbHB; |
1879 | } |
1880 | } |
1881 | } |
1882 | if (bDA || (!bDA && HAinrange)) |
1883 | { |
1884 | return hbDist; |
1885 | } |
1886 | else |
1887 | { |
1888 | return hbNo; |
1889 | } |
1890 | } |
1891 | |
1892 | /* Fixed previously undiscovered bug in the merge |
1893 | code, where the last frame of each hbond disappears. |
1894 | - Erik Marklund, June 1, 2006 */ |
1895 | /* Added the following arguments: |
1896 | * ptmp[] - temporary periodicity hisory |
1897 | * a1 - identity of first acceptor/donor |
1898 | * a2 - identity of second acceptor/donor |
1899 | * - Erik Marklund, FEB 20 2010 */ |
1900 | |
1901 | /* Merging is now done on the fly, so do_merge is most likely obsolete now. |
1902 | * Will do some more testing before removing the function entirely. |
1903 | * - Erik Marklund, MAY 10 2010 */ |
1904 | static void do_merge(t_hbdata *hb, int ntmp, |
1905 | unsigned int htmp[], unsigned int gtmp[], PSTYPEint ptmp[], |
1906 | t_hbond *hb0, t_hbond *hb1, int a1, int a2) |
1907 | { |
1908 | /* Here we need to make sure we're treating periodicity in |
1909 | * the right way for the geminate recombination kinetics. */ |
1910 | |
1911 | int m, mm, n00, n01, nn0, nnframes; |
1912 | PSTYPEint pm; |
1913 | t_pShift *pShift; |
1914 | |
1915 | /* Decide where to start from when merging */ |
1916 | n00 = hb0->n0; |
1917 | n01 = hb1->n0; |
1918 | nn0 = min(n00, n01)(((n00) < (n01)) ? (n00) : (n01) ); |
1919 | nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes)(((n00 + hb0->nframes) > (n01 + hb1->nframes)) ? (n00 + hb0->nframes) : (n01 + hb1->nframes) ) - nn0; |
1920 | /* Initiate tmp arrays */ |
1921 | for (m = 0; (m < ntmp); m++) |
1922 | { |
1923 | htmp[m] = 0; |
1924 | gtmp[m] = 0; |
1925 | ptmp[m] = 0; |
1926 | } |
1927 | /* Fill tmp arrays with values due to first HB */ |
1928 | /* Once again '<' had to be replaced with '<=' |
1929 | to catch the last frame in which the hbond |
1930 | appears. |
1931 | - Erik Marklund, June 1, 2006 */ |
1932 | for (m = 0; (m <= hb0->nframes); m++) |
1933 | { |
1934 | mm = m+n00-nn0; |
1935 | htmp[mm] = is_hb(hb0->h[0], m); |
1936 | if (hb->bGem) |
1937 | { |
1938 | pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0); |
1939 | if (pm > hb->per->nper) |
1940 | { |
1941 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1941, "Illegal shift!"); |
1942 | } |
1943 | else |
1944 | { |
1945 | ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/ |
1946 | } |
1947 | } |
1948 | } |
1949 | /* If we're doing geminate recompbination we usually don't need the distances. |
1950 | * Let's save some memory and time. */ |
1951 | if (TRUE1 || !hb->bGem || hb->per->gemtype == gemAD) |
1952 | { |
1953 | for (m = 0; (m <= hb0->nframes); m++) |
1954 | { |
1955 | mm = m+n00-nn0; |
1956 | gtmp[mm] = is_hb(hb0->g[0], m); |
1957 | } |
1958 | } |
1959 | /* Next HB */ |
1960 | for (m = 0; (m <= hb1->nframes); m++) |
1961 | { |
1962 | mm = m+n01-nn0; |
1963 | htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m); |
1964 | gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m); |
1965 | if (hb->bGem /* && ptmp[mm] != 0 */) |
1966 | { |
1967 | |
1968 | /* If this hbond has been seen before with donor and acceptor swapped, |
1969 | * then we need to find the mirrored (*-1) periodicity vector to truely |
1970 | * merge the hbond history. */ |
1971 | pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper); |
1972 | /* Store index of mirror */ |
1973 | if (pm > hb->per->nper) |
1974 | { |
1975 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1975, "Illegal shift!"); |
1976 | } |
1977 | ptmp[mm] = pm; |
1978 | } |
1979 | } |
1980 | /* Reallocate target array */ |
1981 | if (nnframes > hb0->maxframes) |
1982 | { |
1983 | srenew(hb0->h[0], 4+nnframes/hb->wordlen)(hb0->h[0]) = save_realloc("hb0->h[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1983, (hb0->h[0]), (4+nnframes/hb->wordlen), sizeof(* (hb0->h[0]))); |
1984 | srenew(hb0->g[0], 4+nnframes/hb->wordlen)(hb0->g[0]) = save_realloc("hb0->g[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1984, (hb0->g[0]), (4+nnframes/hb->wordlen), sizeof(* (hb0->g[0]))); |
1985 | } |
1986 | if (NULL((void*)0) != hb->per->pHist) |
1987 | { |
1988 | clearPshift(&(hb->per->pHist[a1][a2])); |
1989 | } |
1990 | |
1991 | /* Copy temp array to target array */ |
1992 | for (m = 0; (m <= nnframes); m++) |
1993 | { |
1994 | _set_hb(hb0->h[0], m, htmp[m]); |
1995 | _set_hb(hb0->g[0], m, gtmp[m]); |
1996 | if (hb->bGem) |
1997 | { |
1998 | addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0); |
1999 | } |
2000 | } |
2001 | |
2002 | /* Set scalar variables */ |
2003 | hb0->n0 = nn0; |
2004 | hb0->maxframes = nnframes; |
2005 | } |
2006 | |
2007 | /* Added argument bContact for nicer output. |
2008 | * Erik Marklund, June 29, 2006 |
2009 | */ |
2010 | static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) |
2011 | { |
2012 | int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp; |
2013 | unsigned int *htmp, *gtmp; |
2014 | PSTYPEint *ptmp; |
2015 | t_hbond *hb0, *hb1; |
2016 | |
2017 | inrnew = hb->nrhb; |
2018 | indnew = hb->nrdist; |
2019 | |
2020 | /* Check whether donors are also acceptors */ |
2021 | printf("Merging hbonds with Acceptor and Donor swapped\n"); |
2022 | |
2023 | ntmp = 2*hb->max_frames; |
2024 | snew(gtmp, ntmp)(gtmp) = save_calloc("gtmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2024, (ntmp), sizeof(*(gtmp))); |
2025 | snew(htmp, ntmp)(htmp) = save_calloc("htmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2025, (ntmp), sizeof(*(htmp))); |
2026 | snew(ptmp, ntmp)(ptmp) = save_calloc("ptmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2026, (ntmp), sizeof(*(ptmp))); |
2027 | for (i = 0; (i < hb->d.nrd); i++) |
2028 | { |
2029 | fprintf(stderrstderr, "\r%d/%d", i+1, hb->d.nrd); |
2030 | id = hb->d.don[i]; |
2031 | ii = hb->a.aptr[id]; |
2032 | for (j = 0; (j < hb->a.nra); j++) |
2033 | { |
2034 | ia = hb->a.acc[j]; |
2035 | jj = hb->d.dptr[ia]; |
2036 | if ((id != ia) && (ii != NOTSET-12345) && (jj != NOTSET-12345) && |
2037 | (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j])))) |
2038 | { |
2039 | hb0 = hb->hbmap[i][j]; |
2040 | hb1 = hb->hbmap[jj][ii]; |
2041 | if (hb0 && hb1 && ISHB(hb0->history[0])(((hb0->history[0]) & 2) == 2) && ISHB(hb1->history[0])(((hb1->history[0]) & 2) == 2)) |
2042 | { |
2043 | do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j); |
2044 | if (ISHB(hb1->history[0])(((hb1->history[0]) & 2) == 2)) |
2045 | { |
2046 | inrnew--; |
2047 | } |
2048 | else if (ISDIST(hb1->history[0])(((hb1->history[0]) & 1) == 1)) |
2049 | { |
2050 | indnew--; |
2051 | } |
2052 | else |
2053 | if (bContact) |
2054 | { |
2055 | gmx_incons("No contact history")_gmx_error("incons", "No contact history", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2055); |
2056 | } |
2057 | else |
2058 | { |
2059 | gmx_incons("Neither hydrogen bond nor distance")_gmx_error("incons", "Neither hydrogen bond nor distance", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2059); |
2060 | } |
2061 | sfree(hb1->h[0])save_free("hb1->h[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2061, (hb1->h[0])); |
2062 | sfree(hb1->g[0])save_free("hb1->g[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2062, (hb1->g[0])); |
2063 | if (hb->bGem) |
2064 | { |
2065 | clearPshift(&(hb->per->pHist[jj][ii])); |
2066 | } |
2067 | hb1->h[0] = NULL((void*)0); |
2068 | hb1->g[0] = NULL((void*)0); |
2069 | hb1->history[0] = hbNo; |
2070 | } |
2071 | } |
2072 | } |
2073 | } |
2074 | fprintf(stderrstderr, "\n"); |
2075 | printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew); |
2076 | printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew); |
2077 | hb->nrhb = inrnew; |
2078 | hb->nrdist = indnew; |
2079 | sfree(gtmp)save_free("gtmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2079, (gtmp)); |
2080 | sfree(htmp)save_free("htmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2080, (htmp)); |
2081 | sfree(ptmp)save_free("ptmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2081, (ptmp)); |
2082 | } |
2083 | |
2084 | static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t) |
2085 | { |
2086 | int i, j, k, n_bound[MAXHH4], nbtot; |
2087 | h_id nhb; |
2088 | |
2089 | |
2090 | /* Set array to 0 */ |
2091 | for (k = 0; (k < MAXHH4); k++) |
2092 | { |
2093 | n_bound[k] = 0; |
2094 | } |
2095 | /* Loop over possible donors */ |
2096 | for (i = 0; (i < hb->d.nrd); i++) |
2097 | { |
2098 | for (j = 0; (j < hb->d.nhydro[i]); j++) |
2099 | { |
2100 | n_bound[hb->d.nhbonds[i][j]]++; |
2101 | } |
2102 | } |
2103 | fprintf(fp, "%12.5e", t); |
2104 | nbtot = 0; |
2105 | for (k = 0; (k < MAXHH4); k++) |
2106 | { |
2107 | fprintf(fp, " %8d", n_bound[k]); |
2108 | nbtot += n_bound[k]*k; |
2109 | } |
2110 | fprintf(fp, " %8d\n", nbtot); |
2111 | } |
2112 | |
2113 | /* Added argument bContact in do_hblife(...). Also |
2114 | * added support for -contact in function body. |
2115 | * - Erik Marklund, May 31, 2006 */ |
2116 | /* Changed the contact code slightly. |
2117 | * - Erik Marklund, June 29, 2006 |
2118 | */ |
2119 | static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact, |
2120 | const output_env_t oenv) |
2121 | { |
2122 | FILE *fp; |
2123 | const char *leg[] = { "p(t)", "t p(t)" }; |
2124 | int *histo; |
2125 | int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0; |
2126 | int nframes = hb->nframes; |
2127 | unsigned int **h; |
2128 | real t, x1, dt; |
2129 | double sum, integral; |
2130 | t_hbond *hbh; |
2131 | |
2132 | snew(h, hb->maxhydro)(h) = save_calloc("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2132, (hb->maxhydro), sizeof(*(h))); |
2133 | snew(histo, nframes+1)(histo) = save_calloc("histo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2133, (nframes+1), sizeof(*(histo))); |
2134 | /* Total number of hbonds analyzed here */ |
2135 | for (i = 0; (i < hb->d.nrd); i++) |
2136 | { |
2137 | for (k = 0; (k < hb->a.nra); k++) |
2138 | { |
2139 | hbh = hb->hbmap[i][k]; |
2140 | if (hbh) |
2141 | { |
2142 | if (bMerge) |
2143 | { |
2144 | if (hbh->h[0]) |
2145 | { |
2146 | h[0] = hbh->h[0]; |
2147 | nhydro = 1; |
2148 | } |
2149 | else |
2150 | { |
2151 | nhydro = 0; |
2152 | } |
2153 | } |
2154 | else |
2155 | { |
2156 | nhydro = 0; |
2157 | for (m = 0; (m < hb->maxhydro); m++) |
2158 | { |
2159 | if (hbh->h[m]) |
2160 | { |
2161 | h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m]; |
2162 | } |
2163 | } |
2164 | } |
2165 | for (nh = 0; (nh < nhydro); nh++) |
2166 | { |
2167 | ohb = 0; |
2168 | j0 = 0; |
2169 | |
2170 | /* Changed '<' into '<=' below, just like I |
2171 | did in the hbm-output-loop in the main code. |
2172 | - Erik Marklund, May 31, 2006 |
2173 | */ |
2174 | for (j = 0; (j <= hbh->nframes); j++) |
2175 | { |
2176 | ihb = is_hb(h[nh], j); |
2177 | if (debug && (ndump < 10)) |
2178 | { |
2179 | fprintf(debug, "%5d %5d\n", j, ihb); |
2180 | } |
2181 | if (ihb != ohb) |
2182 | { |
2183 | if (ihb) |
2184 | { |
2185 | j0 = j; |
2186 | } |
2187 | else |
2188 | { |
2189 | histo[j-j0]++; |
2190 | } |
2191 | ohb = ihb; |
2192 | } |
2193 | } |
2194 | ndump++; |
2195 | } |
2196 | } |
2197 | } |
2198 | } |
2199 | fprintf(stderrstderr, "\n"); |
2200 | if (bContact) |
2201 | { |
2202 | fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv); |
2203 | } |
2204 | else |
2205 | { |
2206 | fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()", |
2207 | oenv); |
2208 | } |
2209 | |
2210 | xvgr_legend(fp, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); |
2211 | j0 = nframes-1; |
2212 | while ((j0 > 0) && (histo[j0] == 0)) |
2213 | { |
2214 | j0--; |
2215 | } |
2216 | sum = 0; |
2217 | for (i = 0; (i <= j0); i++) |
2218 | { |
2219 | sum += histo[i]; |
2220 | } |
2221 | dt = hb->time[1]-hb->time[0]; |
2222 | sum = dt*sum; |
2223 | integral = 0; |
2224 | for (i = 1; (i <= j0); i++) |
2225 | { |
2226 | t = hb->time[i] - hb->time[0] - 0.5*dt; |
2227 | x1 = t*histo[i]/sum; |
2228 | fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1); |
2229 | integral += x1; |
2230 | } |
2231 | integral *= dt; |
2232 | gmx_ffclose(fp); |
2233 | printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral); |
2234 | printf("Note that the lifetime obtained in this manner is close to useless\n"); |
2235 | printf("Use the -ac option instead and check the Forward lifetime\n"); |
2236 | please_cite(stdoutstdout, "Spoel2006b"); |
2237 | sfree(h)save_free("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2237, (h)); |
2238 | sfree(histo)save_free("histo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2238, (histo)); |
2239 | } |
2240 | |
2241 | /* Changed argument bMerge into oneHB to handle contacts properly. |
2242 | * - Erik Marklund, June 29, 2006 |
2243 | */ |
2244 | static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump) |
2245 | { |
2246 | FILE *fp; |
2247 | int i, j, k, m, nd, ihb, idist; |
2248 | int nframes = hb->nframes; |
2249 | gmx_bool bPrint; |
2250 | t_hbond *hbh; |
2251 | |
2252 | if (nDump <= 0) |
2253 | { |
2254 | return; |
2255 | } |
2256 | fp = gmx_ffopen("debug-ac.xvg", "w"); |
2257 | for (j = 0; (j < nframes); j++) |
2258 | { |
2259 | fprintf(fp, "%10.3f", hb->time[j]); |
2260 | for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++) |
2261 | { |
2262 | for (k = 0; (k < hb->a.nra) && (nd < nDump); k++) |
2263 | { |
2264 | bPrint = FALSE0; |
2265 | ihb = idist = 0; |
2266 | hbh = hb->hbmap[i][k]; |
2267 | if (oneHB) |
2268 | { |
2269 | if (hbh->h[0]) |
2270 | { |
2271 | ihb = is_hb(hbh->h[0], j); |
2272 | idist = is_hb(hbh->g[0], j); |
2273 | bPrint = TRUE1; |
2274 | } |
2275 | } |
2276 | else |
2277 | { |
2278 | for (m = 0; (m < hb->maxhydro) && !ihb; m++) |
2279 | { |
2280 | ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m], j)); |
2281 | idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j)); |
2282 | } |
2283 | /* This is not correct! */ |
2284 | /* What isn't correct? -Erik M */ |
2285 | bPrint = TRUE1; |
2286 | } |
2287 | if (bPrint) |
2288 | { |
2289 | fprintf(fp, " %1d-%1d", ihb, idist); |
2290 | nd++; |
2291 | } |
2292 | } |
2293 | } |
2294 | fprintf(fp, "\n"); |
2295 | } |
2296 | gmx_ffclose(fp); |
2297 | } |
2298 | |
2299 | static real calc_dg(real tau, real temp) |
2300 | { |
2301 | real kbt; |
2302 | |
2303 | kbt = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp; |
2304 | if (tau <= 0) |
2305 | { |
2306 | return -666; |
2307 | } |
2308 | else |
2309 | { |
2310 | return kbt*log(kbt*tau/PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))); |
2311 | } |
2312 | } |
2313 | |
2314 | typedef struct { |
2315 | int n0, n1, nparams, ndelta; |
2316 | real kkk[2]; |
2317 | real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt; |
2318 | } t_luzar; |
2319 | |
2320 | static real compute_weighted_rates(int n, real t[], real ct[], real nt[], |
2321 | real kt[], real sigma_ct[], real sigma_nt[], |
2322 | real sigma_kt[], real *k, real *kp, |
2323 | real *sigma_k, real *sigma_kp, |
2324 | real fit_start) |
2325 | { |
2326 | #define NK10 10 |
2327 | int i, j; |
2328 | t_luzar tl; |
2329 | real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2; |
2330 | |
2331 | *sigma_k = 0; |
2332 | *sigma_kp = 0; |
2333 | |
2334 | for (i = 0; (i < n); i++) |
2335 | { |
2336 | if (t[i] >= fit_start) |
2337 | { |
2338 | break; |
2339 | } |
2340 | } |
2341 | tl.n0 = i; |
2342 | tl.n1 = n; |
2343 | tl.nparams = 2; |
2344 | tl.ndelta = 1; |
2345 | tl.t = t; |
2346 | tl.ct = ct; |
2347 | tl.nt = nt; |
2348 | tl.kt = kt; |
2349 | tl.sigma_ct = sigma_ct; |
2350 | tl.sigma_nt = sigma_nt; |
2351 | tl.sigma_kt = sigma_kt; |
2352 | tl.kkk[0] = *k; |
2353 | tl.kkk[1] = *kp; |
2354 | |
2355 | chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */ |
2356 | *k = tl.kkk[0]; |
2357 | *kp = tl.kkk[1] = *kp; |
2358 | tl.ndelta = NK10; |
2359 | for (j = 0; (j < NK10); j++) |
2360 | { |
2361 | /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */ |
2362 | kkk += tl.kkk[0]; |
2363 | kkp += tl.kkk[1]; |
2364 | kk2 += sqr(tl.kkk[0]); |
2365 | kp2 += sqr(tl.kkk[1]); |
2366 | tl.n0++; |
2367 | } |
2368 | *sigma_k = sqrt(kk2/NK10 - sqr(kkk/NK10)); |
2369 | *sigma_kp = sqrt(kp2/NK10 - sqr(kkp/NK10)); |
2370 | |
2371 | return chi2; |
2372 | } |
2373 | |
2374 | static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start, |
2375 | const output_env_t oenv) |
2376 | { |
2377 | FILE *fp; |
2378 | real e_1, fitparm[4]; |
2379 | int i; |
2380 | |
2381 | e_1 = exp(-1); |
2382 | for (i = 0; (i < n); i++) |
2383 | { |
2384 | if (c[i] < e_1) |
2385 | { |
2386 | break; |
2387 | } |
2388 | } |
2389 | if (i < n) |
2390 | { |
2391 | fitparm[0] = t[i]; |
2392 | } |
2393 | else |
2394 | { |
2395 | fitparm[0] = 10; |
2396 | } |
2397 | fitparm[1] = 0.95; |
2398 | do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(), effnEXP2, fitparm, 0); |
2399 | } |
2400 | |
2401 | void analyse_corr(int n, real t[], real ct[], real nt[], real kt[], |
2402 | real sigma_ct[], real sigma_nt[], real sigma_kt[], |
2403 | real fit_start, real temp, real smooth_tail_start, |
2404 | const output_env_t oenv) |
2405 | { |
2406 | int i0, i; |
2407 | real k = 1, kp = 1, kow = 1; |
2408 | real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg; |
2409 | double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0; |
2410 | gmx_bool bError = (sigma_ct != NULL((void*)0)) && (sigma_nt != NULL((void*)0)) && (sigma_kt != NULL((void*)0)); |
2411 | |
2412 | if (smooth_tail_start >= 0) |
2413 | { |
2414 | smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv); |
2415 | smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv); |
2416 | smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv); |
2417 | } |
2418 | for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++) |
2419 | { |
2420 | ; |
2421 | } |
2422 | if (i0 < n-2) |
2423 | { |
2424 | for (i = i0; (i < n); i++) |
2425 | { |
2426 | sc2 += sqr(ct[i]); |
2427 | sn2 += sqr(nt[i]); |
2428 | sk2 += sqr(kt[i]); |
2429 | sck += ct[i]*kt[i]; |
2430 | snk += nt[i]*kt[i]; |
2431 | scn += ct[i]*nt[i]; |
2432 | } |
2433 | printf("Hydrogen bond thermodynamics at T = %g K\n", temp); |
2434 | tmp = (sn2*sc2-sqr(scn)); |
2435 | if ((tmp > 0) && (sn2 > 0)) |
2436 | { |
2437 | k = (sn2*sck-scn*snk)/tmp; |
2438 | kp = (k*scn-snk)/sn2; |
2439 | if (bError) |
2440 | { |
2441 | chi2 = 0; |
2442 | for (i = i0; (i < n); i++) |
2443 | { |
2444 | chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]); |
2445 | } |
2446 | chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt, |
Value stored to 'chi22' is never read | |
2447 | sigma_kt, &k, &kp, |
2448 | &sigma_k, &sigma_kp, fit_start); |
2449 | Q = 0; /* quality_of_fit(chi2, 2);*/ |
2450 | ddg = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp*sigma_k/k; |
2451 | printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n", |
2452 | chi2, Q); |
2453 | printf("The Rate and Delta G are followed by an error estimate\n"); |
2454 | printf("----------------------------------------------------------\n" |
2455 | "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n"); |
2456 | printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n", |
2457 | k, sigma_k, 1/k, calc_dg(1/k, temp), ddg); |
2458 | ddg = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp*sigma_kp/kp; |
2459 | printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n", |
2460 | kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg); |
2461 | } |
2462 | else |
2463 | { |
2464 | chi2 = 0; |
2465 | for (i = i0; (i < n); i++) |
2466 | { |
2467 | chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]); |
2468 | } |
2469 | printf("Fitting parameters chi^2 = %10g\nQ = %10g\n", |
2470 | chi2, Q); |
2471 | printf("--------------------------------------------------\n" |
2472 | "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n"); |
2473 | printf("Forward %10.3f %8.3f %10.3f %10g\n", |
2474 | k, 1/k, calc_dg(1/k, temp), chi2); |
2475 | printf("Backward %10.3f %8.3f %10.3f\n", |
2476 | kp, 1/kp, calc_dg(1/kp, temp)); |
2477 | } |
2478 | } |
2479 | if (sc2 > 0) |
2480 | { |
2481 | kow = 2*sck/sc2; |
2482 | printf("One-way %10.3f %s%8.3f %10.3f\n", |
2483 | kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp)); |
2484 | } |
2485 | else |
2486 | { |
2487 | printf(" - Numerical problems computing HB thermodynamics:\n" |
2488 | "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n", |
2489 | sc2, sn2, sk2, sck, snk, scn); |
2490 | } |
2491 | /* Determine integral of the correlation function */ |
2492 | tau_hb = evaluate_integral(n, t, ct, NULL((void*)0), (t[n-1]-t[0])/2, &dtau); |
2493 | printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb, |
2494 | bError ? " " : "", tau_hb, calc_dg(tau_hb, temp)); |
2495 | e_1 = exp(-1); |
2496 | for (i = 0; (i < n-2); i++) |
2497 | { |
2498 | if ((ct[i] > e_1) && (ct[i+1] <= e_1)) |
2499 | { |
2500 | break; |
2501 | } |
2502 | } |
2503 | if (i < n-2) |
2504 | { |
2505 | /* Determine tau_relax from linear interpolation */ |
2506 | tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]); |
2507 | printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx, |
2508 | tau_rlx, bError ? " " : "", |
2509 | calc_dg(tau_rlx, temp)); |
2510 | } |
2511 | } |
2512 | else |
2513 | { |
2514 | printf("Correlation functions too short to compute thermodynamics\n"); |
2515 | } |
2516 | } |
2517 | |
2518 | void compute_derivative(int nn, real x[], real y[], real dydx[]) |
2519 | { |
2520 | int j; |
2521 | |
2522 | /* Compute k(t) = dc(t)/dt */ |
2523 | for (j = 1; (j < nn-1); j++) |
2524 | { |
2525 | dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]); |
2526 | } |
2527 | /* Extrapolate endpoints */ |
2528 | dydx[0] = 2*dydx[1] - dydx[2]; |
2529 | dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3]; |
2530 | } |
2531 | |
2532 | static void parallel_print(int *data, int nThreads) |
2533 | { |
2534 | /* This prints the donors on which each tread is currently working. */ |
2535 | int i; |
2536 | |
2537 | fprintf(stderrstderr, "\r"); |
2538 | for (i = 0; i < nThreads; i++) |
2539 | { |
2540 | fprintf(stderrstderr, "%-7i", data[i]); |
2541 | } |
2542 | } |
2543 | |
2544 | static void normalizeACF(real *ct, real *gt, int nhb, int len) |
2545 | { |
2546 | real ct_fac, gt_fac; |
2547 | int i; |
2548 | |
2549 | /* Xu and Berne use the same normalization constant */ |
2550 | |
2551 | ct_fac = 1.0/ct[0]; |
2552 | gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb; |
2553 | |
2554 | printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac); |
2555 | for (i = 0; i < len; i++) |
2556 | { |
2557 | ct[i] *= ct_fac; |
2558 | if (gt != NULL((void*)0)) |
2559 | { |
2560 | gt[i] *= gt_fac; |
2561 | } |
2562 | } |
2563 | } |
2564 | |
2565 | /* Added argument bContact in do_hbac(...). Also |
2566 | * added support for -contact in the actual code. |
2567 | * - Erik Marklund, May 31, 2006 */ |
2568 | /* Changed contact code and added argument R2 |
2569 | * - Erik Marklund, June 29, 2006 |
2570 | */ |
2571 | static void do_hbac(const char *fn, t_hbdata *hb, |
2572 | int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start, |
2573 | real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv, |
2574 | const char *gemType, int nThreads, |
2575 | const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit) |
2576 | { |
2577 | FILE *fp; |
2578 | int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets; |
2579 | const char *legNN[] = { |
2580 | "Ac(t)", |
2581 | "Ac'(t)" |
2582 | }; |
2583 | static char **legGem; |
2584 | |
2585 | const char *legLuzar[] = { |
2586 | "Ac\\sfin sys\\v{}\\z{}(t)", |
2587 | "Ac(t)", |
2588 | "Cc\\scontact,hb\\v{}\\z{}(t)", |
2589 | "-dAc\\sfs\\v{}\\z{}/dt" |
2590 | }; |
2591 | gmx_bool bNorm = FALSE0, bOMP = FALSE0; |
2592 | double nhb = 0; |
2593 | int nhbi = 0; |
2594 | real *rhbex = NULL((void*)0), *ht, *gt, *ght, *dght, *kt; |
2595 | real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct; |
2596 | const real tol = 1e-3; |
2597 | int nframes = hb->nframes, nf; |
2598 | unsigned int **h = NULL((void*)0), **g = NULL((void*)0); |
2599 | int nh, nhbonds, nhydro, ngh; |
2600 | t_hbond *hbh; |
2601 | PSTYPEint p, *pfound = NULL((void*)0), np; |
2602 | t_pShift *pHist; |
2603 | int *ptimes = NULL((void*)0), *poff = NULL((void*)0), anhb, n0, mMax = INT_MIN(-2147483647 -1); |
2604 | real **rHbExGem = NULL((void*)0); |
2605 | gmx_bool c; |
2606 | int acType; |
2607 | t_E *E; |
2608 | double *ctdouble, *timedouble, *fittedct; |
2609 | double fittolerance = 0.1; |
2610 | int *dondata = NULL((void*)0), thisThread; |
2611 | |
2612 | enum { |
2613 | AC_NONE, AC_NN, AC_GEM, AC_LUZAR |
2614 | }; |
2615 | |
2616 | #ifdef GMX_OPENMP |
2617 | bOMP = TRUE1; |
2618 | #else |
2619 | bOMP = FALSE0; |
2620 | #endif |
2621 | |
2622 | printf("Doing autocorrelation "); |
2623 | |
2624 | /* Decide what kind of ACF calculations to do. */ |
2625 | if (NN > NN_NONE && NN < NN_NR) |
2626 | { |
2627 | #ifdef HAVE_NN_LOOPS |
2628 | acType = AC_NN; |
2629 | printf("using the energy estimate.\n"); |
2630 | #else |
2631 | acType = AC_NONE; |
2632 | printf("Can't do the NN-loop. Yet.\n"); |
2633 | #endif |
2634 | } |
2635 | else if (hb->bGem) |
2636 | { |
2637 | acType = AC_GEM; |
2638 | printf("according to the reversible geminate recombination model by Omer Markowitch.\n"); |
2639 | |
2640 | nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0); |
2641 | snew(legGem, nSets)(legGem) = save_calloc("legGem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2641, (nSets), sizeof(*(legGem))); |
2642 | for (i = 0; i < nSets; i++) |
2643 | { |
2644 | snew(legGem[i], 128)(legGem[i]) = save_calloc("legGem[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2644, (128), sizeof(*(legGem[i]))); |
2645 | } |
2646 | sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType); |
2647 | if (bBallistic) |
2648 | { |
2649 | sprintf(legGem[1], "Ac'(t)"); |
2650 | } |
2651 | if (bGemFit) |
2652 | { |
2653 | sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType); |
2654 | } |
2655 | |
2656 | } |
2657 | else |
2658 | { |
2659 | acType = AC_LUZAR; |
2660 | printf("according to the theory of Luzar and Chandler.\n"); |
2661 | } |
2662 | fflush(stdoutstdout); |
2663 | |
2664 | /* build hbexist matrix in reals for autocorr */ |
2665 | /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */ |
2666 | n2 = 1; |
2667 | while (n2 < nframes) |
2668 | { |
2669 | n2 *= 2; |
2670 | } |
2671 | |
2672 | nn = nframes/2; |
2673 | |
2674 | if (acType != AC_NN || bOMP) |
2675 | { |
2676 | snew(h, hb->maxhydro)(h) = save_calloc("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2676, (hb->maxhydro), sizeof(*(h))); |
2677 | snew(g, hb->maxhydro)(g) = save_calloc("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2677, (hb->maxhydro), sizeof(*(g))); |
2678 | } |
2679 | |
2680 | /* Dump hbonds for debugging */ |
2681 | dump_ac(hb, bMerge || bContact, nDump); |
2682 | |
2683 | /* Total number of hbonds analyzed here */ |
2684 | nhbonds = 0; |
2685 | ngh = 0; |
2686 | anhb = 0; |
2687 | |
2688 | if (acType != AC_LUZAR && bOMP) |
2689 | { |
2690 | nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads())((((nThreads <= 0) ? 2147483647 : nThreads) < (gmx_omp_get_max_threads ())) ? ((nThreads <= 0) ? 2147483647 : nThreads) : (gmx_omp_get_max_threads ()) ); |
2691 | |
2692 | gmx_omp_set_num_threads(nThreads); |
2693 | snew(dondata, nThreads)(dondata) = save_calloc("dondata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2693, (nThreads), sizeof(*(dondata))); |
2694 | for (i = 0; i < nThreads; i++) |
2695 | { |
2696 | dondata[i] = -1; |
2697 | } |
2698 | printf("ACF calculations parallelized with OpenMP using %i threads.\n" |
2699 | "Expect close to linear scaling over this donor-loop.\n", nThreads); |
2700 | fflush(stdoutstdout); |
2701 | fprintf(stderrstderr, "Donors: [thread no]\n"); |
2702 | { |
2703 | char tmpstr[7]; |
2704 | for (i = 0; i < nThreads; i++) |
2705 | { |
2706 | snprintf(tmpstr, 7, "[%i]", i); |
2707 | fprintf(stderrstderr, "%-7s", tmpstr); |
2708 | } |
2709 | } |
2710 | fprintf(stderrstderr, "\n"); |
2711 | } |
2712 | |
2713 | |
2714 | /* Build the ACF according to acType */ |
2715 | switch (acType) |
2716 | { |
2717 | |
2718 | case AC_NN: |
2719 | #ifdef HAVE_NN_LOOPS |
2720 | /* Here we're using the estimated energy for the hydrogen bonds. */ |
2721 | snew(ct, nn)(ct) = save_calloc("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2721, (nn), sizeof(*(ct))); |
2722 | |
2723 | #pragma omp parallel \ |
2724 | private(i, j, k, nh, E, rhbex, thisThread) \ |
2725 | default(shared) |
2726 | { |
2727 | #pragma omp barrier |
2728 | thisThread = gmx_omp_get_thread_num(); |
2729 | rhbex = NULL((void*)0); |
2730 | |
2731 | snew(rhbex, n2)(rhbex) = save_calloc("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2731, (n2), sizeof(*(rhbex))); |
2732 | memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */ |
2733 | |
2734 | #pragma omp barrier |
2735 | #pragma omp for schedule (dynamic) |
2736 | for (i = 0; i < hb->d.nrd; i++) /* loop over donors */ |
2737 | { |
2738 | if (bOMP) |
2739 | { |
2740 | #pragma omp critical |
2741 | { |
2742 | dondata[thisThread] = i; |
2743 | parallel_print(dondata, nThreads); |
2744 | } |
2745 | } |
2746 | else |
2747 | { |
2748 | fprintf(stderrstderr, "\r %i", i); |
2749 | } |
2750 | |
2751 | for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */ |
2752 | { |
2753 | for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */ |
2754 | { |
2755 | E = hb->hbE.E[i][j][nh]; |
2756 | if (E != NULL((void*)0)) |
2757 | { |
2758 | for (k = 0; k < nframes; k++) |
2759 | { |
2760 | if (E[k] != NONSENSE_E) |
2761 | { |
2762 | rhbex[k] = (real)E[k]; |
2763 | } |
2764 | } |
2765 | |
2766 | low_do_autocorr(NULL((void*)0), oenv, NULL((void*)0), nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0], |
2767 | eacNormal(1<<0), 1, FALSE0, bNorm, FALSE0, 0, -1, 0, 1); |
2768 | #pragma omp critical |
2769 | { |
2770 | for (k = 0; (k < nn); k++) |
2771 | { |
2772 | ct[k] += rhbex[k]; |
2773 | } |
2774 | } |
2775 | } |
2776 | } /* k loop */ |
2777 | } /* j loop */ |
2778 | } /* i loop */ |
2779 | sfree(rhbex)save_free("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2779, (rhbex)); |
2780 | #pragma omp barrier |
2781 | } |
2782 | |
2783 | if (bOMP) |
2784 | { |
2785 | sfree(dondata)save_free("dondata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2785, (dondata)); |
2786 | } |
2787 | normalizeACF(ct, NULL((void*)0), 0, nn); |
2788 | snew(ctdouble, nn)(ctdouble) = save_calloc("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2788, (nn), sizeof(*(ctdouble))); |
2789 | snew(timedouble, nn)(timedouble) = save_calloc("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2789, (nn), sizeof(*(timedouble))); |
2790 | for (j = 0; j < nn; j++) |
2791 | { |
2792 | timedouble[j] = (double)(hb->time[j]); |
2793 | ctdouble[j] = (double)(ct[j]); |
2794 | } |
2795 | |
2796 | /* Remove ballistic term */ |
2797 | /* Ballistic component removal and fitting to the reversible geminate recombination model |
2798 | * will be taken out for the time being. First of all, one can remove the ballistic |
2799 | * component with g_analyze afterwards. Secondly, and more importantly, there are still |
2800 | * problems with the robustness of the fitting to the model. More work is needed. |
2801 | * A third reason is that we're currently using gsl for this and wish to reduce dependence |
2802 | * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with |
2803 | * a BSD-licence that can do the job. |
2804 | * |
2805 | * - Erik Marklund, June 18 2010. |
2806 | */ |
2807 | /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ |
2808 | /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ |
2809 | /* else */ |
2810 | /* printf("\nNumber of data points is less than the number of parameters to fit\n." */ |
2811 | /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ |
2812 | |
2813 | fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)"); |
2814 | xvgr_legend(fp, asize(legNN)((int)(sizeof(legNN)/sizeof((legNN)[0]))), legNN); |
2815 | |
2816 | for (j = 0; (j < nn); j++) |
2817 | { |
2818 | fprintf(fp, "%10g %10g %10g\n", |
2819 | hb->time[j]-hb->time[0], |
2820 | ct[j], |
2821 | ctdouble[j]); |
2822 | } |
2823 | xvgrclose(fp); |
2824 | sfree(ct)save_free("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2824, (ct)); |
2825 | sfree(ctdouble)save_free("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2825, (ctdouble)); |
2826 | sfree(timedouble)save_free("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2826, (timedouble)); |
2827 | #endif /* HAVE_NN_LOOPS */ |
2828 | break; /* case AC_NN */ |
2829 | |
2830 | case AC_GEM: |
2831 | snew(ct, 2*n2)(ct) = save_calloc("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2831, (2*n2), sizeof(*(ct))); |
2832 | memset(ct, 0, 2*n2*sizeof(real)); |
2833 | #ifndef GMX_OPENMP |
2834 | fprintf(stderrstderr, "Donor:\n"); |
2835 | #define __ACDATAct ct |
2836 | #else |
2837 | #define __ACDATAct p_ct |
2838 | #endif |
2839 | |
2840 | #pragma omp parallel \ |
2841 | private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \ |
2842 | pfound, poff, rHbExGem, p, ihb, mMax, \ |
2843 | thisThread, p_ct) \ |
2844 | default(shared) |
2845 | { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */ |
2846 | h = NULL((void*)0); |
2847 | g = NULL((void*)0); |
2848 | thisThread = gmx_omp_get_thread_num(); |
2849 | snew(h, hb->maxhydro)(h) = save_calloc("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2849, (hb->maxhydro), sizeof(*(h))); |
2850 | snew(g, hb->maxhydro)(g) = save_calloc("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2850, (hb->maxhydro), sizeof(*(g))); |
2851 | mMax = INT_MIN(-2147483647 -1); |
2852 | rHbExGem = NULL((void*)0); |
2853 | poff = NULL((void*)0); |
2854 | pfound = NULL((void*)0); |
2855 | p_ct = NULL((void*)0); |
2856 | snew(p_ct, 2*n2)(p_ct) = save_calloc("p_ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2856, (2*n2), sizeof(*(p_ct))); |
2857 | memset(p_ct, 0, 2*n2*sizeof(real)); |
2858 | |
2859 | /* I'm using a chunk size of 1, since I expect \ |
2860 | * the overhead to be really small compared \ |
2861 | * to the actual calculations \ */ |
2862 | #pragma omp for schedule(dynamic,1) nowait |
2863 | for (i = 0; i < hb->d.nrd; i++) |
2864 | { |
2865 | |
2866 | if (bOMP) |
2867 | { |
2868 | #pragma omp critical |
2869 | { |
2870 | dondata[thisThread] = i; |
2871 | parallel_print(dondata, nThreads); |
2872 | } |
2873 | } |
2874 | else |
2875 | { |
2876 | fprintf(stderrstderr, "\r %i", i); |
2877 | } |
2878 | for (k = 0; k < hb->a.nra; k++) |
2879 | { |
2880 | for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) |
2881 | { |
2882 | hbh = hb->hbmap[i][k]; |
2883 | if (hbh) |
2884 | { |
2885 | /* Note that if hb->per->gemtype==gemDD, then distances will be stored in |
2886 | * hb->hbmap[d][a].h array anyway, because the contact flag will be set. |
2887 | * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */ |
2888 | pHist = &(hb->per->pHist[i][k]); |
2889 | if (ISHB(hbh->history[nh])(((hbh->history[nh]) & 2) == 2) && pHist->len != 0) |
2890 | { |
2891 | |
2892 | { |
2893 | h[nh] = hbh->h[nh]; |
2894 | g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL((void*)0); |
2895 | } |
2896 | n0 = hbh->n0; |
2897 | nf = hbh->nframes; |
2898 | /* count the number of periodic shifts encountered and store |
2899 | * them in separate arrays. */ |
2900 | np = 0; |
2901 | for (j = 0; j < pHist->len; j++) |
2902 | { |
2903 | p = pHist->p[j]; |
2904 | for (m = 0; m <= np; m++) |
2905 | { |
2906 | if (m == np) /* p not recognized in list. Add it and set up new array. */ |
2907 | { |
2908 | np++; |
2909 | if (np > hb->per->nper) |
2910 | { |
2911 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2911, "Too many pshifts. Something's utterly wrong here."); |
2912 | } |
2913 | if (m >= mMax) /* Extend the arrays. |
2914 | * Doing it like this, using mMax to keep track of the sizes, |
2915 | * eleviates the need for freeing and re-allocating the arrays |
2916 | * when taking on the next donor-acceptor pair */ |
2917 | { |
2918 | mMax = m; |
2919 | srenew(pfound, np)(pfound) = save_realloc("pfound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2919, (pfound), (np), sizeof(*(pfound))); /* The list of found periodic shifts. */ |
2920 | srenew(rHbExGem, np)(rHbExGem) = save_realloc("rHbExGem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2920, (rHbExGem), (np), sizeof(*(rHbExGem))); /* The hb existence functions (-aver_hb). */ |
2921 | snew(rHbExGem[m], 2*n2)(rHbExGem[m]) = save_calloc("rHbExGem[m]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2921, (2*n2), sizeof(*(rHbExGem[m]))); |
2922 | srenew(poff, np)(poff) = save_realloc("poff", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2922, (poff), (np), sizeof(*(poff))); |
2923 | } |
2924 | |
2925 | { |
2926 | if (rHbExGem != NULL((void*)0) && rHbExGem[m] != NULL((void*)0)) |
2927 | { |
2928 | /* This must be done, as this array was most likey |
2929 | * used to store stuff in some previous iteration. */ |
2930 | memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2)); |
2931 | } |
2932 | else |
2933 | { |
2934 | fprintf(stderrstderr, "rHbExGem not initialized! m = %i\n", m); |
2935 | } |
2936 | } |
2937 | pfound[m] = p; |
2938 | poff[m] = -1; |
2939 | |
2940 | break; |
2941 | } /* m==np */ |
2942 | if (p == pfound[m]) |
2943 | { |
2944 | break; |
2945 | } |
2946 | } /* m: Loop over found shifts */ |
2947 | } /* j: Loop over shifts */ |
2948 | |
2949 | /* Now unpack and disentangle the existence funtions. */ |
2950 | for (j = 0; j < nf; j++) |
2951 | { |
2952 | /* i: donor, |
2953 | * k: acceptor |
2954 | * nh: hydrogen |
2955 | * j: time |
2956 | * p: periodic shift |
2957 | * pfound: list of periodic shifts found for this pair. |
2958 | * poff: list of frame offsets; that is, the first |
2959 | * frame a hbond has a particular periodic shift. */ |
2960 | p = getPshift(*pHist, j+n0); |
2961 | if (p != -1) |
2962 | { |
2963 | for (m = 0; m < np; m++) |
2964 | { |
2965 | if (pfound[m] == p) |
2966 | { |
2967 | break; |
2968 | } |
2969 | if (m == (np-1)) |
2970 | { |
2971 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2971, "Shift not found, but must be there."); |
2972 | } |
2973 | } |
2974 | |
2975 | ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE0 : is_hb(g[nh], j)); |
2976 | if (ihb) |
2977 | { |
2978 | if (poff[m] == -1) |
2979 | { |
2980 | poff[m] = j; /* Here's where the first hbond with shift p is, |
2981 | * relative to the start of h[0].*/ |
2982 | } |
2983 | if (j < poff[m]) |
2984 | { |
2985 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2985, "j<poff[m]"); |
2986 | } |
2987 | rHbExGem[m][j-poff[m]] += 1; |
2988 | } |
2989 | } |
2990 | } |
2991 | |
2992 | /* Now, build ac. */ |
2993 | for (m = 0; m < np; m++) |
2994 | { |
2995 | if (rHbExGem[m][0] > 0 && n0+poff[m] < nn /* && m==0 */) |
2996 | { |
2997 | low_do_autocorr(NULL((void*)0), oenv, NULL((void*)0), nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0], |
2998 | eacNormal(1<<0), 1, FALSE0, bNorm, FALSE0, 0, -1, 0); |
2999 | for (j = 0; (j < nn); j++) |
3000 | { |
3001 | __ACDATAct[j] += rHbExGem[m][j]; |
3002 | } |
3003 | } |
3004 | } /* Building of ac. */ |
3005 | } /* if (ISHB(...*/ |
3006 | } /* if (hbh) */ |
3007 | } /* hydrogen loop */ |
3008 | } /* acceptor loop */ |
3009 | } /* donor loop */ |
3010 | |
3011 | for (m = 0; m <= mMax; m++) |
3012 | { |
3013 | sfree(rHbExGem[m])save_free("rHbExGem[m]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3013, (rHbExGem[m])); |
3014 | } |
3015 | sfree(pfound)save_free("pfound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3015, (pfound)); |
3016 | sfree(poff)save_free("poff", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3016, (poff)); |
3017 | sfree(rHbExGem)save_free("rHbExGem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3017, (rHbExGem)); |
3018 | |
3019 | sfree(h)save_free("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3019, (h)); |
3020 | sfree(g)save_free("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3020, (g)); |
3021 | |
3022 | if (bOMP) |
3023 | { |
3024 | #pragma omp critical |
3025 | { |
3026 | for (i = 0; i < nn; i++) |
3027 | { |
3028 | ct[i] += p_ct[i]; |
3029 | } |
3030 | } |
3031 | sfree(p_ct)save_free("p_ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3031, (p_ct)); |
3032 | } |
3033 | |
3034 | } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */ |
3035 | if (bOMP) |
3036 | { |
3037 | sfree(dondata)save_free("dondata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3037, (dondata)); |
3038 | } |
3039 | |
3040 | normalizeACF(ct, NULL((void*)0), 0, nn); |
3041 | |
3042 | fprintf(stderrstderr, "\n\nACF successfully calculated.\n"); |
3043 | |
3044 | /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */ |
3045 | |
3046 | snew(ctdouble, nn)(ctdouble) = save_calloc("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3046, (nn), sizeof(*(ctdouble))); |
3047 | snew(timedouble, nn)(timedouble) = save_calloc("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3047, (nn), sizeof(*(timedouble))); |
3048 | snew(fittedct, nn)(fittedct) = save_calloc("fittedct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3048, (nn), sizeof(*(fittedct))); |
3049 | |
3050 | for (j = 0; j < nn; j++) |
3051 | { |
3052 | timedouble[j] = (double)(hb->time[j]); |
3053 | ctdouble[j] = (double)(ct[j]); |
3054 | } |
3055 | |
3056 | /* Remove ballistic term */ |
3057 | /* Ballistic component removal and fitting to the reversible geminate recombination model |
3058 | * will be taken out for the time being. First of all, one can remove the ballistic |
3059 | * component with g_analyze afterwards. Secondly, and more importantly, there are still |
3060 | * problems with the robustness of the fitting to the model. More work is needed. |
3061 | * A third reason is that we're currently using gsl for this and wish to reduce dependence |
3062 | * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with |
3063 | * a BSD-licence that can do the job. |
3064 | * |
3065 | * - Erik Marklund, June 18 2010. |
3066 | */ |
3067 | /* if (bBallistic) { */ |
3068 | /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ |
3069 | /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ |
3070 | /* else */ |
3071 | /* printf("\nNumber of data points is less than the number of parameters to fit\n." */ |
3072 | /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ |
3073 | /* } */ |
3074 | /* if (bGemFit) */ |
3075 | /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */ |
3076 | |
3077 | |
3078 | if (bContact) |
3079 | { |
3080 | fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
3081 | } |
3082 | else |
3083 | { |
3084 | fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
3085 | } |
3086 | xvgr_legend(fp, asize(legGem)((int)(sizeof(legGem)/sizeof((legGem)[0]))), (const char**)legGem, oenv); |
3087 | |
3088 | for (j = 0; (j < nn); j++) |
3089 | { |
3090 | fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]); |
3091 | if (bBallistic) |
3092 | { |
3093 | fprintf(fp, " %10g", ctdouble[j]); |
3094 | } |
3095 | if (bGemFit) |
3096 | { |
3097 | fprintf(fp, " %10g", fittedct[j]); |
3098 | } |
3099 | fprintf(fp, "\n"); |
3100 | } |
3101 | xvgrclose(fp); |
3102 | |
3103 | sfree(ctdouble)save_free("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3103, (ctdouble)); |
3104 | sfree(timedouble)save_free("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3104, (timedouble)); |
3105 | sfree(fittedct)save_free("fittedct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3105, (fittedct)); |
3106 | sfree(ct)save_free("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3106, (ct)); |
3107 | |
3108 | break; /* case AC_GEM */ |
3109 | |
3110 | case AC_LUZAR: |
3111 | snew(rhbex, 2*n2)(rhbex) = save_calloc("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3111, (2*n2), sizeof(*(rhbex))); |
3112 | snew(ct, 2*n2)(ct) = save_calloc("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3112, (2*n2), sizeof(*(ct))); |
3113 | snew(gt, 2*n2)(gt) = save_calloc("gt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3113, (2*n2), sizeof(*(gt))); |
3114 | snew(ht, 2*n2)(ht) = save_calloc("ht", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3114, (2*n2), sizeof(*(ht))); |
3115 | snew(ght, 2*n2)(ght) = save_calloc("ght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3115, (2*n2), sizeof(*(ght))); |
3116 | snew(dght, 2*n2)(dght) = save_calloc("dght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3116, (2*n2), sizeof(*(dght))); |
3117 | |
3118 | snew(kt, nn)(kt) = save_calloc("kt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3118, (nn), sizeof(*(kt))); |
3119 | snew(cct, nn)(cct) = save_calloc("cct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3119, (nn), sizeof(*(cct))); |
3120 | |
3121 | for (i = 0; (i < hb->d.nrd); i++) |
3122 | { |
3123 | for (k = 0; (k < hb->a.nra); k++) |
3124 | { |
3125 | nhydro = 0; |
3126 | hbh = hb->hbmap[i][k]; |
3127 | |
3128 | if (hbh) |
3129 | { |
3130 | if (bMerge || bContact) |
3131 | { |
3132 | if (ISHB(hbh->history[0])(((hbh->history[0]) & 2) == 2)) |
3133 | { |
3134 | h[0] = hbh->h[0]; |
3135 | g[0] = hbh->g[0]; |
3136 | nhydro = 1; |
3137 | } |
3138 | } |
3139 | else |
3140 | { |
3141 | for (m = 0; (m < hb->maxhydro); m++) |
3142 | { |
3143 | if (bContact ? ISDIST(hbh->history[m])(((hbh->history[m]) & 1) == 1) : ISHB(hbh->history[m])(((hbh->history[m]) & 2) == 2)) |
3144 | { |
3145 | g[nhydro] = hbh->g[m]; |
3146 | h[nhydro] = hbh->h[m]; |
3147 | nhydro++; |
3148 | } |
3149 | } |
3150 | } |
3151 | |
3152 | nf = hbh->nframes; |
3153 | for (nh = 0; (nh < nhydro); nh++) |
3154 | { |
3155 | int nrint = bContact ? hb->nrdist : hb->nrhb; |
3156 | if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint)) |
3157 | { |
3158 | fprintf(stderrstderr, "\rACF %d/%d", nhbonds+1, nrint); |
3159 | } |
3160 | nhbonds++; |
3161 | for (j = 0; (j < nframes); j++) |
3162 | { |
3163 | /* Changed '<' into '<=' below, just like I did in |
3164 | the hbm-output-loop in the gmx_hbond() block. |
3165 | - Erik Marklund, May 31, 2006 */ |
3166 | if (j <= nf) |
3167 | { |
3168 | ihb = is_hb(h[nh], j); |
3169 | idist = is_hb(g[nh], j); |
3170 | } |
3171 | else |
3172 | { |
3173 | ihb = idist = 0; |
3174 | } |
3175 | rhbex[j] = ihb; |
3176 | /* For contacts: if a second cut-off is provided, use it, |
3177 | * otherwise use g(t) = 1-h(t) */ |
3178 | if (!R2 && bContact) |
3179 | { |
3180 | gt[j] = 1-ihb; |
3181 | } |
3182 | else |
3183 | { |
3184 | gt[j] = idist*(1-ihb); |
3185 | } |
3186 | ht[j] = rhbex[j]; |
3187 | nhb += ihb; |
3188 | } |
3189 | |
3190 | |
3191 | /* The autocorrelation function is normalized after summation only */ |
3192 | low_do_autocorr(NULL((void*)0), oenv, NULL((void*)0), nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0], |
3193 | eacNormal(1<<0), 1, FALSE0, bNorm, FALSE0, 0, -1, 0); |
3194 | |
3195 | /* Cross correlation analysis for thermodynamics */ |
3196 | for (j = nframes; (j < n2); j++) |
3197 | { |
3198 | ht[j] = 0; |
3199 | gt[j] = 0; |
3200 | } |
3201 | |
3202 | cross_corr(n2, ht, gt, dght); |
3203 | |
3204 | for (j = 0; (j < nn); j++) |
3205 | { |
3206 | ct[j] += rhbex[j]; |
3207 | ght[j] += dght[j]; |
3208 | } |
3209 | } |
3210 | } |
3211 | } |
3212 | } |
3213 | fprintf(stderrstderr, "\n"); |
3214 | sfree(h)save_free("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3214, (h)); |
3215 | sfree(g)save_free("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3215, (g)); |
3216 | normalizeACF(ct, ght, nhb, nn); |
3217 | |
3218 | /* Determine tail value for statistics */ |
3219 | tail = 0; |
3220 | tail2 = 0; |
3221 | for (j = nn/2; (j < nn); j++) |
3222 | { |
3223 | tail += ct[j]; |
3224 | tail2 += ct[j]*ct[j]; |
3225 | } |
3226 | tail /= (nn - nn/2); |
3227 | tail2 /= (nn - nn/2); |
3228 | dtail = sqrt(tail2-tail*tail); |
3229 | |
3230 | /* Check whether the ACF is long enough */ |
3231 | if (dtail > tol) |
3232 | { |
3233 | printf("\nWARNING: Correlation function is probably not long enough\n" |
3234 | "because the standard deviation in the tail of C(t) > %g\n" |
3235 | "Tail value (average C(t) over second half of acf): %g +/- %g\n", |
3236 | tol, tail, dtail); |
3237 | } |
3238 | for (j = 0; (j < nn); j++) |
3239 | { |
3240 | cct[j] = ct[j]; |
3241 | ct[j] = (cct[j]-tail)/(1-tail); |
3242 | } |
3243 | /* Compute negative derivative k(t) = -dc(t)/dt */ |
3244 | compute_derivative(nn, hb->time, ct, kt); |
3245 | for (j = 0; (j < nn); j++) |
3246 | { |
3247 | kt[j] = -kt[j]; |
3248 | } |
3249 | |
3250 | |
3251 | if (bContact) |
3252 | { |
3253 | fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
3254 | } |
3255 | else |
3256 | { |
3257 | fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
3258 | } |
3259 | xvgr_legend(fp, asize(legLuzar)((int)(sizeof(legLuzar)/sizeof((legLuzar)[0]))), legLuzar, oenv); |
3260 | |
3261 | |
3262 | for (j = 0; (j < nn); j++) |
3263 | { |
3264 | fprintf(fp, "%10g %10g %10g %10g %10g\n", |
3265 | hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]); |
3266 | } |
3267 | gmx_ffclose(fp); |
3268 | |
3269 | analyse_corr(nn, hb->time, ct, ght, kt, NULL((void*)0), NULL((void*)0), NULL((void*)0), |
3270 | fit_start, temp, smooth_tail_start, oenv); |
3271 | |
3272 | do_view(oenv, fn, NULL((void*)0)); |
3273 | sfree(rhbex)save_free("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3273, (rhbex)); |
3274 | sfree(ct)save_free("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3274, (ct)); |
3275 | sfree(gt)save_free("gt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3275, (gt)); |
3276 | sfree(ht)save_free("ht", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3276, (ht)); |
3277 | sfree(ght)save_free("ght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3277, (ght)); |
3278 | sfree(dght)save_free("dght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3278, (dght)); |
3279 | sfree(cct)save_free("cct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3279, (cct)); |
3280 | sfree(kt)save_free("kt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3280, (kt)); |
3281 | /* sfree(h); */ |
3282 | /* sfree(g); */ |
3283 | |
3284 | break; /* case AC_LUZAR */ |
3285 | |
3286 | default: |
3287 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3287, "Unrecognized type of ACF-calulation. acType = %i.", acType); |
3288 | } /* switch (acType) */ |
3289 | } |
3290 | |
3291 | static void init_hbframe(t_hbdata *hb, int nframes, real t) |
3292 | { |
3293 | int i, j, m; |
3294 | |
3295 | hb->time[nframes] = t; |
3296 | hb->nhb[nframes] = 0; |
3297 | hb->ndist[nframes] = 0; |
3298 | for (i = 0; (i < max_hx7); i++) |
3299 | { |
3300 | hb->nhx[nframes][i] = 0; |
3301 | } |
3302 | /* Loop invalidated */ |
3303 | if (hb->bHBmap && 0) |
3304 | { |
3305 | for (i = 0; (i < hb->d.nrd); i++) |
3306 | { |
3307 | for (j = 0; (j < hb->a.nra); j++) |
3308 | { |
3309 | for (m = 0; (m < hb->maxhydro); m++) |
3310 | { |
3311 | if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m]) |
3312 | { |
3313 | set_hb(hb, i, m, j, nframes, HB_NO0); |
3314 | } |
3315 | } |
3316 | } |
3317 | } |
3318 | } |
3319 | /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/ |
3320 | } |
3321 | |
3322 | static void analyse_donor_props(const char *fn, t_hbdata *hb, int nframes, real t, |
3323 | const output_env_t oenv) |
3324 | { |
3325 | static FILE *fp = NULL((void*)0); |
3326 | const char *leg[] = { "Nbound", "Nfree" }; |
3327 | int i, j, k, nbound, nb, nhtot; |
3328 | |
3329 | if (!fn) |
3330 | { |
3331 | return; |
3332 | } |
3333 | if (!fp) |
3334 | { |
3335 | fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv); |
3336 | xvgr_legend(fp, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); |
3337 | } |
3338 | nbound = 0; |
3339 | nhtot = 0; |
3340 | for (i = 0; (i < hb->d.nrd); i++) |
3341 | { |
3342 | for (k = 0; (k < hb->d.nhydro[i]); k++) |
3343 | { |
3344 | nb = 0; |
3345 | nhtot++; |
3346 | for (j = 0; (j < hb->a.nra) && (nb == 0); j++) |
3347 | { |
3348 | if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] && |
3349 | is_hb(hb->hbmap[i][j]->h[k], nframes)) |
3350 | { |
3351 | nb = 1; |
3352 | } |
3353 | } |
3354 | nbound += nb; |
3355 | } |
3356 | } |
3357 | fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound); |
3358 | } |
3359 | |
3360 | static void dump_hbmap(t_hbdata *hb, |
3361 | int nfile, t_filenm fnm[], gmx_bool bTwo, |
3362 | gmx_bool bContact, int isize[], int *index[], char *grpnames[], |
3363 | t_atoms *atoms) |
3364 | { |
3365 | FILE *fp, *fplog; |
3366 | int ddd, hhh, aaa, i, j, k, m, grp; |
3367 | char ds[32], hs[32], as[32]; |
3368 | gmx_bool first; |
3369 | |
3370 | fp = opt2FILE("-hbn", nfile, fnm, "w")gmx_ffopen(opt2fn("-hbn", nfile, fnm), "w"); |
3371 | if (opt2bSet("-g", nfile, fnm)) |
3372 | { |
3373 | fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w"); |
3374 | fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor"); |
3375 | } |
3376 | else |
3377 | { |
3378 | fplog = NULL((void*)0); |
3379 | } |
3380 | for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++) |
3381 | { |
3382 | fprintf(fp, "[ %s ]", grpnames[grp]); |
3383 | for (i = 0; i < isize[grp]; i++) |
3384 | { |
3385 | fprintf(fp, (i%15) ? " " : "\n"); |
3386 | fprintf(fp, " %4u", index[grp][i]+1); |
3387 | } |
3388 | fprintf(fp, "\n"); |
3389 | /* |
3390 | Added -contact support below. |
3391 | - Erik Marklund, May 29, 2006 |
3392 | */ |
3393 | if (!bContact) |
3394 | { |
3395 | fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]); |
3396 | for (i = 0; (i < hb->d.nrd); i++) |
3397 | { |
3398 | if (hb->d.grp[i] == grp) |
3399 | { |
3400 | for (j = 0; (j < hb->d.nhydro[i]); j++) |
3401 | { |
3402 | fprintf(fp, " %4u %4u", hb->d.don[i]+1, |
3403 | hb->d.hydro[i][j]+1); |
3404 | } |
3405 | fprintf(fp, "\n"); |
3406 | } |
3407 | } |
3408 | first = TRUE1; |
3409 | fprintf(fp, "[ acceptors_%s ]", grpnames[grp]); |
3410 | for (i = 0; (i < hb->a.nra); i++) |
3411 | { |
3412 | if (hb->a.grp[i] == grp) |
3413 | { |
3414 | fprintf(fp, (i%15 && !first) ? " " : "\n"); |
3415 | fprintf(fp, " %4u", hb->a.acc[i]+1); |
3416 | first = FALSE0; |
3417 | } |
3418 | } |
3419 | fprintf(fp, "\n"); |
3420 | } |
3421 | } |
3422 | if (bTwo) |
3423 | { |
3424 | fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" : |
3425 | "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]); |
3426 | } |
3427 | else |
3428 | { |
3429 | fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]); |
3430 | } |
3431 | |
3432 | for (i = 0; (i < hb->d.nrd); i++) |
3433 | { |
3434 | ddd = hb->d.don[i]; |
3435 | for (k = 0; (k < hb->a.nra); k++) |
3436 | { |
3437 | aaa = hb->a.acc[k]; |
3438 | for (m = 0; (m < hb->d.nhydro[i]); m++) |
3439 | { |
3440 | if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m])(((hb->hbmap[i][k]->history[m]) & 2) == 2)) |
3441 | { |
3442 | sprintf(ds, "%s", mkatomname(atoms, ddd)); |
3443 | sprintf(as, "%s", mkatomname(atoms, aaa)); |
3444 | if (bContact) |
3445 | { |
3446 | fprintf(fp, " %6u %6u\n", ddd+1, aaa+1); |
3447 | if (fplog) |
3448 | { |
3449 | fprintf(fplog, "%12s %12s\n", ds, as); |
3450 | } |
3451 | } |
3452 | else |
3453 | { |
3454 | hhh = hb->d.hydro[i][m]; |
3455 | sprintf(hs, "%s", mkatomname(atoms, hhh)); |
3456 | fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1); |
3457 | if (fplog) |
3458 | { |
3459 | fprintf(fplog, "%12s %12s %12s\n", ds, hs, as); |
3460 | } |
3461 | } |
3462 | } |
3463 | } |
3464 | } |
3465 | } |
3466 | gmx_ffclose(fp); |
3467 | if (fplog) |
3468 | { |
3469 | gmx_ffclose(fplog); |
3470 | } |
3471 | } |
3472 | |
3473 | /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template. |
3474 | * It mimics add_frames() and init_frame() to some extent. */ |
3475 | static void sync_hbdata(t_hbdata *p_hb, int nframes) |
3476 | { |
3477 | int i; |
3478 | if (nframes >= p_hb->max_frames) |
3479 | { |
3480 | p_hb->max_frames += 4096; |
3481 | srenew(p_hb->nhb, p_hb->max_frames)(p_hb->nhb) = save_realloc("p_hb->nhb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3481, (p_hb->nhb), (p_hb->max_frames), sizeof(*(p_hb-> nhb))); |
3482 | srenew(p_hb->ndist, p_hb->max_frames)(p_hb->ndist) = save_realloc("p_hb->ndist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3482, (p_hb->ndist), (p_hb->max_frames), sizeof(*(p_hb ->ndist))); |
3483 | srenew(p_hb->n_bound, p_hb->max_frames)(p_hb->n_bound) = save_realloc("p_hb->n_bound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3483, (p_hb->n_bound), (p_hb->max_frames), sizeof(*(p_hb ->n_bound))); |
3484 | srenew(p_hb->nhx, p_hb->max_frames)(p_hb->nhx) = save_realloc("p_hb->nhx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3484, (p_hb->nhx), (p_hb->max_frames), sizeof(*(p_hb-> nhx))); |
3485 | if (p_hb->bDAnr) |
3486 | { |
3487 | srenew(p_hb->danr, p_hb->max_frames)(p_hb->danr) = save_realloc("p_hb->danr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3487, (p_hb->danr), (p_hb->max_frames), sizeof(*(p_hb ->danr))); |
3488 | } |
3489 | memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); |
3490 | memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); |
3491 | p_hb->nhb[nframes] = 0; |
3492 | p_hb->ndist[nframes] = 0; |
3493 | |
3494 | } |
3495 | p_hb->nframes = nframes; |
3496 | /* for (i=0;) */ |
3497 | /* { */ |
3498 | /* p_hb->nhx[nframes][i] */ |
3499 | /* } */ |
3500 | memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx7); /* zero the helix count for this frame */ |
3501 | |
3502 | /* hb->per will remain constant througout the frame loop, |
3503 | * even though the data its members point to will change, |
3504 | * hence no need for re-syncing. */ |
3505 | } |
3506 | |
3507 | int gmx_hbond(int argc, char *argv[]) |
3508 | { |
3509 | const char *desc[] = { |
3510 | "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are", |
3511 | "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor", |
3512 | "(zero is extended) and the distance Donor - Acceptor", |
3513 | "(or Hydrogen - Acceptor using [TT]-noda[tt]).", |
3514 | "OH and NH groups are regarded as donors, O is an acceptor always,", |
3515 | "N is an acceptor by default, but this can be switched using", |
3516 | "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected", |
3517 | "to the first preceding non-hydrogen atom.[PAR]", |
3518 | |
3519 | "You need to specify two groups for analysis, which must be either", |
3520 | "identical or non-overlapping. All hydrogen bonds between the two", |
3521 | "groups are analyzed.[PAR]", |
3522 | |
3523 | "If you set [TT]-shell[tt], you will be asked for an additional index group", |
3524 | "which should contain exactly one atom. In this case, only hydrogen", |
3525 | "bonds between atoms within the shell distance from the one atom are", |
3526 | "considered.[PAR]", |
3527 | |
3528 | "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler", |
3529 | "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).", |
3530 | "If contact kinetics are analyzed by using the -contact option, then", |
3531 | "n(t) can be defined as either all pairs that are not within contact distance r at time t", |
3532 | "(corresponding to leaving the -r2 option at the default value 0) or all pairs that", |
3533 | "are within distance r2 (corresponding to setting a second cut-off value with option -r2).", |
3534 | "See mentioned literature for more details and definitions." |
3535 | "[PAR]", |
3536 | |
3537 | /* "It is also possible to analyse specific hydrogen bonds with", |
3538 | "[TT]-sel[tt]. This index file must contain a group of atom triplets", |
3539 | "Donor Hydrogen Acceptor, in the following way:[PAR]", |
3540 | */ |
3541 | "[TT]", |
3542 | "[ selected ][BR]", |
3543 | " 20 21 24[BR]", |
3544 | " 25 26 29[BR]", |
3545 | " 1 3 6[BR]", |
3546 | "[tt][BR]", |
3547 | "Note that the triplets need not be on separate lines.", |
3548 | "Each atom triplet specifies a hydrogen bond to be analyzed,", |
3549 | "note also that no check is made for the types of atoms.[PAR]", |
3550 | |
3551 | "[BB]Output:[bb][BR]", |
3552 | "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]", |
3553 | "[TT]-ac[tt]: average over all autocorrelations of the existence", |
3554 | "functions (either 0 or 1) of all hydrogen bonds.[BR]", |
3555 | "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]", |
3556 | "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]", |
3557 | "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time", |
3558 | "where n and n+i stand for residue numbers and i ranges from 0 to 6.", |
3559 | "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated", |
3560 | "with helices in proteins.[BR]", |
3561 | "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors", |
3562 | "for selected groups, all hydrogen bonded atoms from all groups and", |
3563 | "all solvent atoms involved in insertion.[BR]", |
3564 | "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all", |
3565 | "frames, this also contains information on solvent insertion", |
3566 | "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]", |
3567 | "index file.[BR]", |
3568 | "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for", |
3569 | "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]", |
3570 | "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to", |
3571 | "compare results to Raman Spectroscopy.", |
3572 | "[PAR]", |
3573 | "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]", |
3574 | "require an amount of memory proportional to the total numbers of donors", |
3575 | "times the total number of acceptors in the selected group(s)." |
3576 | }; |
3577 | |
3578 | static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1; |
3579 | static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1; |
3580 | static gmx_bool bNitAcc = TRUE1, bDA = TRUE1, bMerge = TRUE1; |
3581 | static int nDump = 0, nFitPoints = 100; |
3582 | static int nThreads = 0, nBalExp = 4; |
3583 | |
3584 | static gmx_bool bContact = FALSE0, bBallistic = FALSE0, bGemFit = FALSE0; |
3585 | static real logAfterTime = 10, gemBallistic = 0.2; /* ps */ |
3586 | static const char *NNtype[] = {NULL((void*)0), "none", "binary", "oneOverR3", "dipole", NULL((void*)0)}; |
3587 | |
3588 | /* options */ |
3589 | t_pargs pa [] = { |
3590 | { "-a", FALSE0, etREAL, {&acut}, |
3591 | "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" }, |
3592 | { "-r", FALSE0, etREAL, {&rcut}, |
3593 | "Cutoff radius (nm, X - Acceptor, see next option)" }, |
3594 | { "-da", FALSE0, etBOOL, {&bDA}, |
3595 | "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" }, |
3596 | { "-r2", FALSE0, etREAL, {&r2cut}, |
3597 | "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"}, |
3598 | { "-abin", FALSE0, etREAL, {&abin}, |
3599 | "Binwidth angle distribution (degrees)" }, |
3600 | { "-rbin", FALSE0, etREAL, {&rbin}, |
3601 | "Binwidth distance distribution (nm)" }, |
3602 | { "-nitacc", FALSE0, etBOOL, {&bNitAcc}, |
3603 | "Regard nitrogen atoms as acceptors" }, |
3604 | { "-contact", FALSE0, etBOOL, {&bContact}, |
3605 | "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" }, |
3606 | { "-shell", FALSE0, etREAL, {&rshell}, |
3607 | "when > 0, only calculate hydrogen bonds within # nm shell around " |
3608 | "one particle" }, |
3609 | { "-fitstart", FALSE0, etREAL, {&fit_start}, |
3610 | "Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation. With [TT]-gemfit[tt] we suggest [TT]-fitstart 0[tt]" }, |
3611 | { "-fitend", FALSE0, etREAL, {&fit_end}, |
3612 | "Time (ps) to which to stop fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation (only with [TT]-gemfit[tt])" }, |
3613 | { "-temp", FALSE0, etREAL, {&temp}, |
3614 | "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" }, |
3615 | { "-smooth", FALSE0, etREAL, {&smooth_tail_start}, |
3616 | "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" }, |
3617 | { "-dump", FALSE0, etINT, {&nDump}, |
3618 | "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" }, |
3619 | { "-max_hb", FALSE0, etREAL, {&maxnhb}, |
3620 | "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" }, |
3621 | { "-merge", FALSE0, etBOOL, {&bMerge}, |
3622 | "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." }, |
3623 | { "-geminate", FALSE0, etENUM, {gemType}, |
3624 | "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."}, |
3625 | { "-diff", FALSE0, etREAL, {&D}, |
3626 | "HIDDENDffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."}, |
3627 | #ifdef GMX_OPENMP |
3628 | { "-nthreads", FALSE0, etINT, {&nThreads}, |
3629 | "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"}, |
3630 | #endif |
3631 | }; |
3632 | const char *bugs[] = { |
3633 | "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being." |
3634 | }; |
3635 | t_filenm fnm[] = { |
3636 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, |
3637 | { efTPX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
3638 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
3639 | /* { efNDX, "-sel", "select", ffOPTRD },*/ |
3640 | { efXVG, "-num", "hbnum", ffWRITE1<<2 }, |
3641 | { efLOG, "-g", "hbond", ffOPTWR(1<<2| 1<<3) }, |
3642 | { efXVG, "-ac", "hbac", ffOPTWR(1<<2| 1<<3) }, |
3643 | { efXVG, "-dist", "hbdist", ffOPTWR(1<<2| 1<<3) }, |
3644 | { efXVG, "-ang", "hbang", ffOPTWR(1<<2| 1<<3) }, |
3645 | { efXVG, "-hx", "hbhelix", ffOPTWR(1<<2| 1<<3) }, |
3646 | { efNDX, "-hbn", "hbond", ffOPTWR(1<<2| 1<<3) }, |
3647 | { efXPM, "-hbm", "hbmap", ffOPTWR(1<<2| 1<<3) }, |
3648 | { efXVG, "-don", "donor", ffOPTWR(1<<2| 1<<3) }, |
3649 | { efXVG, "-dan", "danum", ffOPTWR(1<<2| 1<<3) }, |
3650 | { efXVG, "-life", "hblife", ffOPTWR(1<<2| 1<<3) }, |
3651 | { efXVG, "-nhbdist", "nhbdist", ffOPTWR(1<<2| 1<<3) } |
3652 | |
3653 | }; |
3654 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
3655 | |
3656 | char hbmap [HB_NR(1<<2)] = { ' ', 'o', '-', '*' }; |
3657 | const char *hbdesc[HB_NR(1<<2)] = { "None", "Present", "Inserted", "Present & Inserted" }; |
3658 | t_rgb hbrgb [HB_NR(1<<2)] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} }; |
3659 | |
3660 | t_trxstatus *status; |
3661 | int trrStatus = 1; |
3662 | t_topology top; |
3663 | t_inputrec ir; |
3664 | t_pargs *ppa; |
3665 | int npargs, natoms, nframes = 0, shatom; |
3666 | int *isize; |
3667 | char **grpnames; |
3668 | atom_id **index; |
3669 | rvec *x, hbox; |
3670 | matrix box; |
3671 | real t, ccut, dist = 0.0, ang = 0.0; |
3672 | double max_nhb, aver_nhb, aver_dist; |
3673 | int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel; |
3674 | int xi, yi, zi, ai; |
3675 | int xj, yj, zj, aj, xjj, yjj, zjj; |
3676 | int xk, yk, zk, ak, xkk, ykk, zkk; |
3677 | gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric; |
3678 | int *adist, *rdist, *aptr, *rprt; |
3679 | int grp, nabin, nrbin, bin, resdist, ihb; |
3680 | char **leg; |
3681 | t_hbdata *hb, *hbptr; |
3682 | FILE *fp, *fpins = NULL((void*)0), *fpnhb = NULL((void*)0); |
3683 | t_gridcell ***grid; |
3684 | t_ncell *icell, *jcell, *kcell; |
3685 | ivec ngrid; |
3686 | unsigned char *datable; |
3687 | output_env_t oenv; |
3688 | int gemmode, NN; |
3689 | PSTYPEint peri = 0; |
3690 | t_E E; |
3691 | int ii, jj, hh, actual_nThreads; |
3692 | int threadNr = 0; |
3693 | gmx_bool bGem, bNN, bParallel; |
3694 | t_gemParams *params = NULL((void*)0); |
3695 | gmx_bool bEdge_yjj, bEdge_xjj, bOMP; |
3696 | |
3697 | t_hbdata **p_hb = NULL((void*)0); /* one per thread, then merge after the frame loop */ |
3698 | int **p_adist = NULL((void*)0), **p_rdist = NULL((void*)0); /* a histogram for each thread. */ |
3699 | |
3700 | #ifdef GMX_OPENMP |
3701 | bOMP = TRUE1; |
3702 | #else |
3703 | bOMP = FALSE0; |
3704 | #endif |
3705 | |
3706 | npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))); |
3707 | ppa = add_acf_pargs(&npargs, pa); |
3708 | |
3709 | if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13), NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, |
3710 | ppa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, asize(bugs)((int)(sizeof(bugs)/sizeof((bugs)[0]))), bugs, &oenv)) |
3711 | { |
3712 | return 0; |
3713 | } |
3714 | |
3715 | /* NN-loop? If so, what estimator to use ?*/ |
3716 | NN = 1; |
3717 | /* Outcommented for now DvdS 2010-07-13 |
3718 | while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0) |
3719 | NN++; |
3720 | if (NN == NN_NR) |
3721 | gmx_fatal(FARGS, "Invalid NN-loop type."); |
3722 | */ |
3723 | bNN = FALSE0; |
3724 | for (i = 2; bNN == FALSE0 && i < NN_NR; i++) |
3725 | { |
3726 | bNN = bNN || NN == i; |
3727 | } |
3728 | |
3729 | if (NN > NN_NONE && bMerge) |
3730 | { |
3731 | bMerge = FALSE0; |
3732 | } |
3733 | |
3734 | /* geminate recombination? If so, which flavor? */ |
3735 | gemmode = 1; |
3736 | while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0) |
3737 | { |
3738 | gemmode++; |
3739 | } |
3740 | if (gemmode == gemNR) |
3741 | { |
3742 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3742, "Invalid recombination type."); |
3743 | } |
3744 | |
3745 | bGem = FALSE0; |
3746 | for (i = 2; bGem == FALSE0 && i < gemNR; i++) |
3747 | { |
3748 | bGem = bGem || gemmode == i; |
3749 | } |
3750 | |
3751 | if (bGem) |
3752 | { |
3753 | printf("Geminate recombination: %s\n", gemType[gemmode]); |
3754 | if (bContact) |
3755 | { |
3756 | if (gemmode != gemDD) |
3757 | { |
3758 | printf("Turning off -contact option...\n"); |
3759 | bContact = FALSE0; |
3760 | } |
3761 | } |
3762 | else |
3763 | { |
3764 | if (gemmode == gemDD) |
3765 | { |
3766 | printf("Turning on -contact option...\n"); |
3767 | bContact = TRUE1; |
3768 | } |
3769 | } |
3770 | if (bMerge) |
3771 | { |
3772 | if (gemmode == gemAA) |
3773 | { |
3774 | printf("Turning off -merge option...\n"); |
3775 | bMerge = FALSE0; |
3776 | } |
3777 | } |
3778 | else |
3779 | { |
3780 | if (gemmode != gemAA) |
3781 | { |
3782 | printf("Turning on -merge option...\n"); |
3783 | bMerge = TRUE1; |
3784 | } |
3785 | } |
3786 | } |
3787 | |
3788 | /* process input */ |
3789 | bSelected = FALSE0; |
3790 | ccut = cos(acut*DEG2RAD(3.14159265358979323846/180.0)); |
3791 | |
3792 | if (bContact) |
3793 | { |
3794 | if (bSelected) |
3795 | { |
3796 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3796, "Can not analyze selected contacts."); |
3797 | } |
3798 | if (!bDA) |
3799 | { |
3800 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3800, "Can not analyze contact between H and A: turn off -noda"); |
3801 | } |
3802 | } |
3803 | |
3804 | /* Initiate main data structure! */ |
3805 | bHBmap = (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
3806 | opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
3807 | opt2bSet("-hbn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
3808 | opt2bSet("-hbm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
3809 | bGem); |
3810 | |
3811 | if (opt2bSet("-nhbdist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
3812 | { |
3813 | const char *leg[MAXHH4+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" }; |
3814 | fpnhb = xvgropen(opt2fn("-nhbdist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
3815 | "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv); |
3816 | xvgr_legend(fpnhb, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); |
3817 | } |
3818 | |
3819 | hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bMerge || bContact, bGem, gemmode); |
3820 | |
3821 | /* get topology */ |
3822 | read_tpx_top(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &ir, box, &natoms, NULL((void*)0), NULL((void*)0), NULL((void*)0), &top); |
3823 | |
3824 | snew(grpnames, grNR)(grpnames) = save_calloc("grpnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3824, (grNR), sizeof(*(grpnames))); |
3825 | snew(index, grNR)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3825, (grNR), sizeof(*(index))); |
3826 | snew(isize, grNR)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3826, (grNR), sizeof(*(isize))); |
3827 | /* Make Donor-Acceptor table */ |
3828 | snew(datable, top.atoms.nr)(datable) = save_calloc("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3828, (top.atoms.nr), sizeof(*(datable))); |
3829 | gen_datable(index[0], isize[0], datable, top.atoms.nr); |
3830 | |
3831 | if (bSelected) |
3832 | { |
3833 | /* analyze selected hydrogen bonds */ |
3834 | printf("Select group with selected atoms:\n"); |
3835 | get_index(&(top.atoms), opt2fn("-sel", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
3836 | 1, &nsel, index, grpnames); |
3837 | if (nsel % 3) |
3838 | { |
3839 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3839, "Number of atoms in group '%s' not a multiple of 3\n" |
3840 | "and therefore cannot contain triplets of " |
3841 | "Donor-Hydrogen-Acceptor", grpnames[0]); |
3842 | } |
3843 | bTwo = FALSE0; |
3844 | |
3845 | for (i = 0; (i < nsel); i += 3) |
3846 | { |
3847 | int dd = index[0][i]; |
3848 | int aa = index[0][i+2]; |
3849 | /* int */ hh = index[0][i+1]; |
3850 | add_dh (&hb->d, dd, hh, i, datable); |
3851 | add_acc(&hb->a, aa, i); |
3852 | /* Should this be here ? */ |
3853 | snew(hb->d.dptr, top.atoms.nr)(hb->d.dptr) = save_calloc("hb->d.dptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3853, (top.atoms.nr), sizeof(*(hb->d.dptr))); |
3854 | snew(hb->a.aptr, top.atoms.nr)(hb->a.aptr) = save_calloc("hb->a.aptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3854, (top.atoms.nr), sizeof(*(hb->a.aptr))); |
3855 | add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri); |
3856 | } |
3857 | printf("Analyzing %d selected hydrogen bonds from '%s'\n", |
3858 | isize[0], grpnames[0]); |
3859 | } |
3860 | else |
3861 | { |
3862 | /* analyze all hydrogen bonds: get group(s) */ |
3863 | printf("Specify 2 groups to analyze:\n"); |
3864 | get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
3865 | 2, isize, index, grpnames); |
3866 | |
3867 | /* check if we have two identical or two non-overlapping groups */ |
3868 | bTwo = isize[0] != isize[1]; |
3869 | for (i = 0; (i < isize[0]) && !bTwo; i++) |
3870 | { |
3871 | bTwo = index[0][i] != index[1][i]; |
3872 | } |
3873 | if (bTwo) |
3874 | { |
3875 | printf("Checking for overlap in atoms between %s and %s\n", |
3876 | grpnames[0], grpnames[1]); |
3877 | for (i = 0; i < isize[1]; i++) |
3878 | { |
3879 | if (ISINGRP(datable[index[1][i]])(((datable[index[1][i]]) & 4) == 4)) |
3880 | { |
3881 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3881, "Partial overlap between groups '%s' and '%s'", |
3882 | grpnames[0], grpnames[1]); |
3883 | } |
3884 | } |
3885 | /* |
3886 | printf("Checking for overlap in atoms between %s and %s\n", |
3887 | grpnames[0],grpnames[1]); |
3888 | for (i=0; i<isize[0]; i++) |
3889 | for (j=0; j<isize[1]; j++) |
3890 | if (index[0][i] == index[1][j]) |
3891 | gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'", |
3892 | grpnames[0],grpnames[1]); |
3893 | */ |
3894 | } |
3895 | if (bTwo) |
3896 | { |
3897 | printf("Calculating %s " |
3898 | "between %s (%d atoms) and %s (%d atoms)\n", |
3899 | bContact ? "contacts" : "hydrogen bonds", |
3900 | grpnames[0], isize[0], grpnames[1], isize[1]); |
3901 | } |
3902 | else |
3903 | { |
3904 | fprintf(stderrstderr, "Calculating %s in %s (%d atoms)\n", |
3905 | bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]); |
3906 | } |
3907 | } |
3908 | sfree(datable)save_free("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3908, (datable)); |
3909 | |
3910 | /* search donors and acceptors in groups */ |
3911 | snew(datable, top.atoms.nr)(datable) = save_calloc("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3911, (top.atoms.nr), sizeof(*(datable))); |
3912 | for (i = 0; (i < grNR); i++) |
3913 | { |
3914 | if ( ((i == gr0) && !bSelected ) || |
3915 | ((i == gr1) && bTwo )) |
3916 | { |
3917 | gen_datable(index[i], isize[i], datable, top.atoms.nr); |
3918 | if (bContact) |
3919 | { |
3920 | search_acceptors(&top, isize[i], index[i], &hb->a, i, |
3921 | bNitAcc, TRUE1, (bTwo && (i == gr0)) || !bTwo, datable); |
3922 | search_donors (&top, isize[i], index[i], &hb->d, i, |
3923 | TRUE1, (bTwo && (i == gr1)) || !bTwo, datable); |
3924 | } |
3925 | else |
3926 | { |
3927 | search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE0, TRUE1, datable); |
3928 | search_donors (&top, isize[i], index[i], &hb->d, i, FALSE0, TRUE1, datable); |
3929 | } |
3930 | if (bTwo) |
3931 | { |
3932 | clear_datable_grp(datable, top.atoms.nr); |
3933 | } |
3934 | } |
3935 | } |
3936 | sfree(datable)save_free("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3936, (datable)); |
3937 | printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra); |
3938 | /*if (bSelected) |
3939 | snew(donors[gr0D], dons[gr0D].nrd);*/ |
3940 | |
3941 | if (bHBmap) |
3942 | { |
3943 | printf("Making hbmap structure..."); |
3944 | /* Generate hbond data structure */ |
3945 | mk_hbmap(hb); |
3946 | printf("done.\n"); |
3947 | } |
3948 | |
3949 | #ifdef HAVE_NN_LOOPS |
3950 | if (bNN) |
3951 | { |
3952 | mk_hbEmap(hb, 0); |
3953 | } |
3954 | #endif |
3955 | |
3956 | if (bGem) |
3957 | { |
3958 | printf("Making per structure..."); |
3959 | /* Generate hbond data structure */ |
3960 | mk_per(hb); |
3961 | printf("done.\n"); |
3962 | } |
3963 | |
3964 | /* check input */ |
3965 | bStop = FALSE0; |
3966 | if (hb->d.nrd + hb->a.nra == 0) |
3967 | { |
3968 | printf("No Donors or Acceptors found\n"); |
3969 | bStop = TRUE1; |
3970 | } |
3971 | if (!bStop) |
3972 | { |
3973 | if (hb->d.nrd == 0) |
3974 | { |
3975 | printf("No Donors found\n"); |
3976 | bStop = TRUE1; |
3977 | } |
3978 | if (hb->a.nra == 0) |
3979 | { |
3980 | printf("No Acceptors found\n"); |
3981 | bStop = TRUE1; |
3982 | } |
3983 | } |
3984 | if (bStop) |
3985 | { |
3986 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3986, "Nothing to be done"); |
3987 | } |
3988 | |
3989 | shatom = 0; |
3990 | if (rshell > 0) |
3991 | { |
3992 | int shisz; |
3993 | atom_id *shidx; |
3994 | char *shgrpnm; |
3995 | /* get index group with atom for shell */ |
3996 | do |
3997 | { |
3998 | printf("Select atom for shell (1 atom):\n"); |
3999 | get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
4000 | 1, &shisz, &shidx, &shgrpnm); |
4001 | if (shisz != 1) |
4002 | { |
4003 | printf("group contains %d atoms, should be 1 (one)\n", shisz); |
4004 | } |
4005 | } |
4006 | while (shisz != 1); |
4007 | shatom = shidx[0]; |
4008 | printf("Will calculate hydrogen bonds within a shell " |
4009 | "of %g nm around atom %i\n", rshell, shatom+1); |
4010 | } |
4011 | |
4012 | /* Analyze trajectory */ |
4013 | natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box); |
4014 | if (natoms > top.atoms.nr) |
4015 | { |
4016 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4016, "Topology (%d atoms) does not match trajectory (%d atoms)", |
4017 | top.atoms.nr, natoms); |
4018 | } |
4019 | |
4020 | bBox = ir.ePBC != epbcNONE; |
4021 | grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid); |
4022 | nabin = acut/abin; |
4023 | nrbin = rcut/rbin; |
4024 | snew(adist, nabin+1)(adist) = save_calloc("adist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4024, (nabin+1), sizeof(*(adist))); |
4025 | snew(rdist, nrbin+1)(rdist) = save_calloc("rdist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4025, (nrbin+1), sizeof(*(rdist))); |
4026 | |
4027 | if (bGem && !bBox) |
4028 | { |
4029 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4029, "Can't do geminate recombination without periodic box."); |
4030 | } |
4031 | |
4032 | bParallel = FALSE0; |
4033 | |
4034 | #ifndef GMX_OPENMP |
4035 | #define __ADISTadist adist |
4036 | #define __RDISTrdist rdist |
4037 | #define __HBDATAhb hb |
4038 | #else /* GMX_OPENMP ================================================== \ |
4039 | * Set up the OpenMP stuff, | |
4040 | * like the number of threads and such | |
4041 | * Also start the parallel loop. | |
4042 | */ |
4043 | #define __ADISTadist p_adist[threadNr] |
4044 | #define __RDISTrdist p_rdist[threadNr] |
4045 | #define __HBDATAhb p_hb[threadNr] |
4046 | #endif |
4047 | if (bOMP) |
4048 | { |
4049 | bParallel = !bSelected; |
4050 | |
4051 | if (bParallel) |
4052 | { |
4053 | actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads())((((nThreads <= 0) ? 2147483647 : nThreads) < (gmx_omp_get_max_threads ())) ? ((nThreads <= 0) ? 2147483647 : nThreads) : (gmx_omp_get_max_threads ()) ); |
4054 | |
4055 | gmx_omp_set_num_threads(actual_nThreads); |
4056 | printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads); |
4057 | fflush(stdoutstdout); |
4058 | } |
4059 | else |
4060 | { |
4061 | actual_nThreads = 1; |
4062 | } |
4063 | |
4064 | snew(p_hb, actual_nThreads)(p_hb) = save_calloc("p_hb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4064, (actual_nThreads), sizeof(*(p_hb))); |
4065 | snew(p_adist, actual_nThreads)(p_adist) = save_calloc("p_adist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4065, (actual_nThreads), sizeof(*(p_adist))); |
4066 | snew(p_rdist, actual_nThreads)(p_rdist) = save_calloc("p_rdist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4066, (actual_nThreads), sizeof(*(p_rdist))); |
4067 | for (i = 0; i < actual_nThreads; i++) |
4068 | { |
4069 | snew(p_hb[i], 1)(p_hb[i]) = save_calloc("p_hb[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4069, (1), sizeof(*(p_hb[i]))); |
4070 | snew(p_adist[i], nabin+1)(p_adist[i]) = save_calloc("p_adist[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4070, (nabin+1), sizeof(*(p_adist[i]))); |
4071 | snew(p_rdist[i], nrbin+1)(p_rdist[i]) = save_calloc("p_rdist[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4071, (nrbin+1), sizeof(*(p_rdist[i]))); |
4072 | |
4073 | p_hb[i]->max_frames = 0; |
4074 | p_hb[i]->nhb = NULL((void*)0); |
4075 | p_hb[i]->ndist = NULL((void*)0); |
4076 | p_hb[i]->n_bound = NULL((void*)0); |
4077 | p_hb[i]->time = NULL((void*)0); |
4078 | p_hb[i]->nhx = NULL((void*)0); |
4079 | |
4080 | p_hb[i]->bHBmap = hb->bHBmap; |
4081 | p_hb[i]->bDAnr = hb->bDAnr; |
4082 | p_hb[i]->bGem = hb->bGem; |
4083 | p_hb[i]->wordlen = hb->wordlen; |
4084 | p_hb[i]->nframes = hb->nframes; |
4085 | p_hb[i]->maxhydro = hb->maxhydro; |
4086 | p_hb[i]->danr = hb->danr; |
4087 | p_hb[i]->d = hb->d; |
4088 | p_hb[i]->a = hb->a; |
4089 | p_hb[i]->hbmap = hb->hbmap; |
4090 | p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */ |
4091 | p_hb[i]->per = hb->per; |
4092 | |
4093 | #ifdef HAVE_NN_LOOPS |
4094 | p_hb[i]->hbE = hb->hbE; |
4095 | #endif |
4096 | |
4097 | p_hb[i]->nrhb = 0; |
4098 | p_hb[i]->nrdist = 0; |
4099 | } |
4100 | } |
4101 | |
4102 | /* Make a thread pool here, |
4103 | * instead of forking anew at every frame. */ |
4104 | |
4105 | #pragma omp parallel \ |
4106 | firstprivate(i) \ |
4107 | private(j, h, ii, jj, hh, E, \ |
4108 | xi, yi, zi, xj, yj, zj, threadNr, \ |
4109 | dist, ang, peri, icell, jcell, \ |
4110 | grp, ogrp, ai, aj, xjj, yjj, zjj, \ |
4111 | xk, yk, zk, ihb, id, resdist, \ |
4112 | xkk, ykk, zkk, kcell, ak, k, bTric, \ |
4113 | bEdge_xjj, bEdge_yjj) \ |
4114 | default(shared) |
4115 | { /* Start of parallel region */ |
4116 | threadNr = gmx_omp_get_thread_num(); |
4117 | |
4118 | do |
4119 | { |
4120 | |
4121 | bTric = bBox && TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0); |
4122 | |
4123 | if (bOMP) |
4124 | { |
4125 | sync_hbdata(p_hb[threadNr], nframes); |
4126 | } |
4127 | #pragma omp single |
4128 | { |
4129 | build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut, |
4130 | rshell, ngrid, grid); |
4131 | reset_nhbonds(&(hb->d)); |
4132 | |
4133 | if (debug && bDebug) |
4134 | { |
4135 | dump_grid(debug, ngrid, grid); |
4136 | } |
4137 | |
4138 | add_frames(hb, nframes); |
4139 | init_hbframe(hb, nframes, output_env_conv_time(oenv, t)); |
4140 | |
4141 | if (hb->bDAnr) |
4142 | { |
4143 | count_da_grid(ngrid, grid, hb->danr[nframes]); |
4144 | } |
4145 | } /* omp single */ |
4146 | |
4147 | if (bOMP) |
4148 | { |
4149 | p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */ |
4150 | } |
4151 | |
4152 | if (bNN) |
4153 | { |
4154 | #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */ |
4155 | /* Loop over all atom pairs and estimate interaction energy */ |
4156 | |
4157 | #pragma omp single |
4158 | { |
4159 | addFramesNN(hb, nframes); |
4160 | } |
4161 | |
4162 | #pragma omp barrier |
4163 | #pragma omp for schedule(dynamic) |
4164 | for (i = 0; i < hb->d.nrd; i++) |
4165 | { |
4166 | for (j = 0; j < hb->a.nra; j++) |
4167 | { |
4168 | for (h = 0; |
4169 | h < (bContact ? 1 : hb->d.nhydro[i]); |
4170 | h++) |
4171 | { |
4172 | if (i == hb->d.nrd || j == hb->a.nra) |
4173 | { |
4174 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4174, "out of bounds"); |
4175 | } |
4176 | |
4177 | /* Get the real atom ids */ |
4178 | ii = hb->d.don[i]; |
4179 | jj = hb->a.acc[j]; |
4180 | hh = hb->d.hydro[i][h]; |
4181 | |
4182 | /* Estimate the energy from the geometry */ |
4183 | E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d)); |
4184 | /* Store the energy */ |
4185 | storeHbEnergy(hb, i, j, h, E, nframes); |
4186 | } |
4187 | } |
4188 | } |
4189 | #endif /* HAVE_NN_LOOPS */ |
4190 | } /* if (bNN)*/ |
4191 | else |
4192 | { |
4193 | if (bSelected) |
4194 | { |
4195 | |
4196 | #pragma omp single |
4197 | { |
4198 | /* Do not parallelize this just yet. */ |
4199 | /* int ii; */ |
4200 | for (ii = 0; (ii < nsel); ii++) |
4201 | { |
4202 | int dd = index[0][i]; |
4203 | int aa = index[0][i+2]; |
4204 | /* int */ hh = index[0][i+1]; |
4205 | ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box, |
4206 | hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); |
4207 | |
4208 | if (ihb) |
4209 | { |
4210 | /* add to index if not already there */ |
4211 | /* Add a hbond */ |
4212 | add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri); |
4213 | } |
4214 | } |
4215 | } /* omp single */ |
4216 | } /* if (bSelected) */ |
4217 | else |
4218 | { |
4219 | |
4220 | #pragma omp single |
4221 | { |
4222 | if (bGem) |
4223 | { |
4224 | calcBoxProjection(box, hb->per->P); |
4225 | } |
4226 | |
4227 | /* loop over all gridcells (xi,yi,zi) */ |
4228 | /* Removed confusing macro, DvdS 27/12/98 */ |
4229 | |
4230 | } |
4231 | /* The outer grid loop will have to do for now. */ |
4232 | #pragma omp for schedule(dynamic) |
4233 | for (xi = 0; xi < ngrid[XX0]; xi++) |
4234 | { |
4235 | for (yi = 0; (yi < ngrid[YY1]); yi++) |
4236 | { |
4237 | for (zi = 0; (zi < ngrid[ZZ2]); zi++) |
4238 | { |
4239 | |
4240 | /* loop over donor groups gr0 (always) and gr1 (if necessary) */ |
4241 | for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++) |
4242 | { |
4243 | icell = &(grid[zi][yi][xi].d[grp]); |
4244 | |
4245 | if (bTwo) |
4246 | { |
4247 | ogrp = 1-grp; |
4248 | } |
4249 | else |
4250 | { |
4251 | ogrp = grp; |
4252 | } |
4253 | |
4254 | /* loop over all hydrogen atoms from group (grp) |
4255 | * in this gridcell (icell) |
4256 | */ |
4257 | for (ai = 0; (ai < icell->nr); ai++) |
4258 | { |
4259 | i = icell->atoms[ai]; |
4260 | |
4261 | /* loop over all adjacent gridcells (xj,yj,zj) */ |
4262 | for (zjj = grid_loop_begin(ngrid[ZZ2], zi, bTric, FALSE0); |
4263 | zjj <= grid_loop_end(ngrid[ZZ2], zi, bTric, FALSE0); |
4264 | zjj++) |
4265 | { |
4266 | zj = grid_mod(zjj, ngrid[ZZ2]); |
4267 | bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ2] - 1); |
4268 | for (yjj = grid_loop_begin(ngrid[YY1], yi, bTric, bEdge_yjj); |
4269 | yjj <= grid_loop_end(ngrid[YY1], yi, bTric, bEdge_yjj); |
4270 | yjj++) |
4271 | { |
4272 | yj = grid_mod(yjj, ngrid[YY1]); |
4273 | bEdge_xjj = |
4274 | (yj == 0) || (yj == ngrid[YY1] - 1) || |
4275 | (zj == 0) || (zj == ngrid[ZZ2] - 1); |
4276 | for (xjj = grid_loop_begin(ngrid[XX0], xi, bTric, bEdge_xjj); |
4277 | xjj <= grid_loop_end(ngrid[XX0], xi, bTric, bEdge_xjj); |
4278 | xjj++) |
4279 | { |
4280 | xj = grid_mod(xjj, ngrid[XX0]); |
4281 | jcell = &(grid[zj][yj][xj].a[ogrp]); |
4282 | /* loop over acceptor atoms from other group (ogrp) |
4283 | * in this adjacent gridcell (jcell) |
4284 | */ |
4285 | for (aj = 0; (aj < jcell->nr); aj++) |
4286 | { |
4287 | j = jcell->atoms[aj]; |
4288 | |
4289 | /* check if this once was a h-bond */ |
4290 | peri = -1; |
4291 | ihb = is_hbond(__HBDATAhb, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box, |
4292 | hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); |
4293 | |
4294 | if (ihb) |
4295 | { |
4296 | /* add to index if not already there */ |
4297 | /* Add a hbond */ |
4298 | add_hbond(__HBDATAhb, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri); |
4299 | |
4300 | /* make angle and distance distributions */ |
4301 | if (ihb == hbHB && !bContact) |
4302 | { |
4303 | if (dist > rcut) |
4304 | { |
4305 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4305, "distance is higher than what is allowed for an hbond: %f", dist); |
4306 | } |
4307 | ang *= RAD2DEG(180.0/3.14159265358979323846); |
4308 | __ADISTadist[(int)( ang/abin)]++; |
4309 | __RDISTrdist[(int)(dist/rbin)]++; |
4310 | if (!bTwo) |
4311 | { |
4312 | int id, ia; |
4313 | if ((id = donor_index(&hb->d, grp, i)_donor_index(&hb->d, grp, i, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4313)) == NOTSET-12345) |
4314 | { |
4315 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4315, "Invalid donor %d", i); |
4316 | } |
4317 | if ((ia = acceptor_index(&hb->a, ogrp, j)_acceptor_index(&hb->a, ogrp, j, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4317)) == NOTSET-12345) |
4318 | { |
4319 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4319, "Invalid acceptor %d", j); |
4320 | } |
4321 | resdist = abs(top.atoms.atom[i].resind- |
4322 | top.atoms.atom[j].resind); |
4323 | if (resdist >= max_hx7) |
4324 | { |
4325 | resdist = max_hx7-1; |
4326 | } |
4327 | __HBDATAhb->nhx[nframes][resdist]++; |
4328 | } |
4329 | } |
4330 | |
4331 | } |
4332 | } /* for aj */ |
4333 | } /* for xjj */ |
4334 | } /* for yjj */ |
4335 | } /* for zjj */ |
4336 | } /* for ai */ |
4337 | } /* for grp */ |
4338 | } /* for xi,yi,zi */ |
4339 | } |
4340 | } |
4341 | } /* if (bSelected) {...} else */ |
4342 | |
4343 | |
4344 | /* Better wait for all threads to finnish using x[] before updating it. */ |
4345 | k = nframes; |
4346 | #pragma omp barrier |
4347 | #pragma omp critical |
4348 | { |
4349 | /* Sum up histograms and counts from p_hb[] into hb */ |
4350 | if (bOMP) |
4351 | { |
4352 | hb->nhb[k] += p_hb[threadNr]->nhb[k]; |
4353 | hb->ndist[k] += p_hb[threadNr]->ndist[k]; |
4354 | for (j = 0; j < max_hx7; j++) |
4355 | { |
4356 | hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j]; |
4357 | } |
4358 | } |
4359 | } |
4360 | |
4361 | /* Here are a handful of single constructs |
4362 | * to share the workload a bit. The most |
4363 | * important one is of course the last one, |
4364 | * where there's a potential bottleneck in form |
4365 | * of slow I/O. */ |
4366 | #pragma omp barrier |
4367 | #pragma omp single |
4368 | { |
4369 | if (hb != NULL((void*)0)) |
4370 | { |
4371 | analyse_donor_props(opt2fn_null("-don", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), hb, k, t, oenv); |
4372 | } |
4373 | } |
4374 | |
4375 | #pragma omp single |
4376 | { |
4377 | if (fpnhb) |
4378 | { |
4379 | do_nhb_dist(fpnhb, hb, t); |
4380 | } |
4381 | } |
4382 | } /* if (bNN) {...} else + */ |
4383 | |
4384 | #pragma omp single |
4385 | { |
4386 | trrStatus = (read_next_x(oenv, status, &t, x, box)); |
4387 | nframes++; |
4388 | } |
4389 | |
4390 | #pragma omp barrier |
4391 | } |
4392 | while (trrStatus); |
4393 | |
4394 | if (bOMP) |
4395 | { |
4396 | #pragma omp critical |
4397 | { |
4398 | hb->nrhb += p_hb[threadNr]->nrhb; |
4399 | hb->nrdist += p_hb[threadNr]->nrdist; |
4400 | } |
4401 | /* Free parallel datastructures */ |
4402 | sfree(p_hb[threadNr]->nhb)save_free("p_hb[threadNr]->nhb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4402, (p_hb[threadNr]->nhb)); |
4403 | sfree(p_hb[threadNr]->ndist)save_free("p_hb[threadNr]->ndist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4403, (p_hb[threadNr]->ndist)); |
4404 | sfree(p_hb[threadNr]->nhx)save_free("p_hb[threadNr]->nhx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4404, (p_hb[threadNr]->nhx)); |
4405 | |
4406 | #pragma omp for |
4407 | for (i = 0; i < nabin; i++) |
4408 | { |
4409 | for (j = 0; j < actual_nThreads; j++) |
4410 | { |
4411 | |
4412 | adist[i] += p_adist[j][i]; |
4413 | } |
4414 | } |
4415 | #pragma omp for |
4416 | for (i = 0; i <= nrbin; i++) |
4417 | { |
4418 | for (j = 0; j < actual_nThreads; j++) |
4419 | { |
4420 | rdist[i] += p_rdist[j][i]; |
4421 | } |
4422 | } |
4423 | |
4424 | sfree(p_adist[threadNr])save_free("p_adist[threadNr]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4424, (p_adist[threadNr])); |
4425 | sfree(p_rdist[threadNr])save_free("p_rdist[threadNr]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4425, (p_rdist[threadNr])); |
4426 | } |
4427 | } /* End of parallel region */ |
4428 | if (bOMP) |
4429 | { |
4430 | sfree(p_adist)save_free("p_adist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4430, (p_adist)); |
4431 | sfree(p_rdist)save_free("p_rdist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4431, (p_rdist)); |
4432 | } |
4433 | |
4434 | if (nframes < 2 && (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))) |
4435 | { |
4436 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4436, "Cannot calculate autocorrelation of life times with less than two frames"); |
4437 | } |
4438 | |
4439 | free_grid(ngrid, &grid); |
4440 | |
4441 | close_trj(status); |
4442 | if (fpnhb) |
4443 | { |
4444 | gmx_ffclose(fpnhb); |
4445 | } |
4446 | |
4447 | /* Compute maximum possible number of different hbonds */ |
4448 | if (maxnhb > 0) |
4449 | { |
4450 | max_nhb = maxnhb; |
4451 | } |
4452 | else |
4453 | { |
4454 | max_nhb = 0.5*(hb->d.nrd*hb->a.nra); |
4455 | } |
4456 | /* Added support for -contact below. |
4457 | * - Erik Marklund, May 29-31, 2006 */ |
4458 | /* Changed contact code. |
4459 | * - Erik Marklund, June 29, 2006 */ |
4460 | if (bHBmap && !bNN) |
4461 | { |
4462 | if (hb->nrhb == 0) |
4463 | { |
4464 | printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds"); |
4465 | } |
4466 | else |
4467 | { |
4468 | printf("Found %d different %s in trajectory\n" |
4469 | "Found %d different atom-pairs within %s distance\n", |
4470 | hb->nrhb, bContact ? "contacts" : "hydrogen bonds", |
4471 | hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding"); |
4472 | |
4473 | /*Control the pHist.*/ |
4474 | |
4475 | if (bMerge) |
4476 | { |
4477 | merge_hb(hb, bTwo, bContact); |
4478 | } |
4479 | |
4480 | if (opt2bSet("-hbn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4481 | { |
4482 | dump_hbmap(hb, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, bTwo, bContact, isize, index, grpnames, &top.atoms); |
4483 | } |
4484 | |
4485 | /* Moved the call to merge_hb() to a line BEFORE dump_hbmap |
4486 | * to make the -hbn and -hmb output match eachother. |
4487 | * - Erik Marklund, May 30, 2006 */ |
4488 | } |
4489 | } |
4490 | /* Print out number of hbonds and distances */ |
4491 | aver_nhb = 0; |
4492 | aver_dist = 0; |
4493 | fp = xvgropen(opt2fn("-num", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bContact ? "Contacts" : |
4494 | "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv); |
4495 | snew(leg, 2)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4495, (2), sizeof(*(leg))); |
4496 | snew(leg[0], STRLEN)(leg[0]) = save_calloc("leg[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4496, (4096), sizeof(*(leg[0]))); |
4497 | snew(leg[1], STRLEN)(leg[1]) = save_calloc("leg[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4497, (4096), sizeof(*(leg[1]))); |
4498 | sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds"); |
4499 | sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut); |
4500 | xvgr_legend(fp, 2, (const char**)leg, oenv); |
4501 | sfree(leg[1])save_free("leg[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4501, (leg[1])); |
4502 | sfree(leg[0])save_free("leg[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4502, (leg[0])); |
4503 | sfree(leg)save_free("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4503, (leg)); |
4504 | for (i = 0; (i < nframes); i++) |
4505 | { |
4506 | fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]); |
4507 | aver_nhb += hb->nhb[i]; |
4508 | aver_dist += hb->ndist[i]; |
4509 | } |
4510 | gmx_ffclose(fp); |
4511 | aver_nhb /= nframes; |
4512 | aver_dist /= nframes; |
4513 | /* Print HB distance distribution */ |
4514 | if (opt2bSet("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4515 | { |
4516 | long sum; |
4517 | |
4518 | sum = 0; |
4519 | for (i = 0; i < nrbin; i++) |
4520 | { |
4521 | sum += rdist[i]; |
4522 | } |
4523 | |
4524 | fp = xvgropen(opt2fn("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
4525 | "Hydrogen Bond Distribution", |
4526 | bDA ? |
4527 | "Donor - Acceptor Distance (nm)" : |
4528 | "Hydrogen - Acceptor Distance (nm)", "", oenv); |
4529 | for (i = 0; i < nrbin; i++) |
4530 | { |
4531 | fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum)); |
4532 | } |
4533 | gmx_ffclose(fp); |
4534 | } |
4535 | |
4536 | /* Print HB angle distribution */ |
4537 | if (opt2bSet("-ang", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4538 | { |
4539 | long sum; |
4540 | |
4541 | sum = 0; |
4542 | for (i = 0; i < nabin; i++) |
4543 | { |
4544 | sum += adist[i]; |
4545 | } |
4546 | |
4547 | fp = xvgropen(opt2fn("-ang", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
4548 | "Hydrogen Bond Distribution", |
4549 | "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv); |
4550 | for (i = 0; i < nabin; i++) |
4551 | { |
4552 | fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum)); |
4553 | } |
4554 | gmx_ffclose(fp); |
4555 | } |
4556 | |
4557 | /* Print HB in alpha-helix */ |
4558 | if (opt2bSet("-hx", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4559 | { |
4560 | fp = xvgropen(opt2fn("-hx", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
4561 | "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv); |
4562 | xvgr_legend(fp, NRHXTYPES7, hxtypenames, oenv); |
4563 | for (i = 0; i < nframes; i++) |
4564 | { |
4565 | fprintf(fp, "%10g", hb->time[i]); |
4566 | for (j = 0; j < max_hx7; j++) |
4567 | { |
4568 | fprintf(fp, " %6d", hb->nhx[i][j]); |
4569 | } |
4570 | fprintf(fp, "\n"); |
4571 | } |
4572 | gmx_ffclose(fp); |
4573 | } |
4574 | if (!bNN) |
4575 | { |
4576 | printf("Average number of %s per timeframe %.3f out of %g possible\n", |
4577 | bContact ? "contacts" : "hbonds", |
4578 | bContact ? aver_dist : aver_nhb, max_nhb); |
4579 | } |
4580 | |
4581 | /* Do Autocorrelation etc. */ |
4582 | if (hb->bHBmap) |
4583 | { |
4584 | /* |
4585 | Added support for -contact in ac and hbm calculations below. |
4586 | - Erik Marklund, May 29, 2006 |
4587 | */ |
4588 | ivec itmp; |
4589 | rvec rtmp; |
4590 | if (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4591 | { |
4592 | please_cite(stdoutstdout, "Spoel2006b"); |
4593 | } |
4594 | if (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4595 | { |
4596 | char *gemstring = NULL((void*)0); |
4597 | |
4598 | if (bGem || bNN) |
4599 | { |
4600 | params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end, |
4601 | gemBallistic, nBalExp); |
4602 | if (params == NULL((void*)0)) |
4603 | { |
4604 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4604, "Could not initiate t_gemParams params."); |
4605 | } |
4606 | } |
4607 | gemstring = strdup(gemType[hb->per->gemtype])(__extension__ (__builtin_constant_p (gemType[hb->per-> gemtype]) && ((size_t)(const void *)((gemType[hb-> per->gemtype]) + 1) - (size_t)(const void *)(gemType[hb-> per->gemtype]) == 1) ? (((const char *) (gemType[hb->per ->gemtype]))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t ) 1) : ({ size_t __len = strlen (gemType[hb->per->gemtype ]) + 1; char *__retval = (char *) malloc (__len); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval, gemType [hb->per->gemtype], __len); __retval; })) : __strdup (gemType [hb->per->gemtype]))); |
4608 | do_hbac(opt2fn("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), hb, nDump, |
4609 | bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv, |
4610 | gemstring, nThreads, NN, bBallistic, bGemFit); |
4611 | } |
4612 | if (opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4613 | { |
4614 | do_hblife(opt2fn("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), hb, bMerge, bContact, oenv); |
4615 | } |
4616 | if (opt2bSet("-hbm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
4617 | { |
4618 | t_matrix mat; |
4619 | int id, ia, hh, x, y; |
4620 | |
4621 | if ((nframes > 0) && (hb->nrhb > 0)) |
4622 | { |
4623 | mat.nx = nframes; |
4624 | mat.ny = hb->nrhb; |
4625 | |
4626 | snew(mat.matrix, mat.nx)(mat.matrix) = save_calloc("mat.matrix", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4626, (mat.nx), sizeof(*(mat.matrix))); |
4627 | for (x = 0; (x < mat.nx); x++) |
4628 | { |
4629 | snew(mat.matrix[x], mat.ny)(mat.matrix[x]) = save_calloc("mat.matrix[x]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4629, (mat.ny), sizeof(*(mat.matrix[x]))); |
4630 | } |
4631 | y = 0; |
4632 | for (id = 0; (id < hb->d.nrd); id++) |
4633 | { |
4634 | for (ia = 0; (ia < hb->a.nra); ia++) |
4635 | { |
4636 | for (hh = 0; (hh < hb->maxhydro); hh++) |
4637 | { |
4638 | if (hb->hbmap[id][ia]) |
4639 | { |
4640 | if (ISHB(hb->hbmap[id][ia]->history[hh])(((hb->hbmap[id][ia]->history[hh]) & 2) == 2)) |
4641 | { |
4642 | /* Changed '<' into '<=' in the for-statement below. |
4643 | * It fixed the previously undiscovered bug that caused |
4644 | * the last occurance of an hbond/contact to not be |
4645 | * set in mat.matrix. Have a look at any old -hbm-output |
4646 | * and you will notice that the last column is allways empty. |
4647 | * - Erik Marklund May 30, 2006 |
4648 | */ |
4649 | for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++) |
4650 | { |
4651 | int nn0 = hb->hbmap[id][ia]->n0; |
4652 | range_check(y, 0, mat.ny)_range_check(y, 0, mat.ny, ((void*)0),"y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4652); |
4653 | mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x); |
4654 | } |
4655 | y++; |
4656 | } |
4657 | } |
4658 | } |
4659 | } |
4660 | } |
4661 | mat.axis_x = hb->time; |
4662 | snew(mat.axis_y, mat.ny)(mat.axis_y) = save_calloc("mat.axis_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4662, (mat.ny), sizeof(*(mat.axis_y))); |
4663 | for (j = 0; j < mat.ny; j++) |
4664 | { |
4665 | mat.axis_y[j] = j; |
4666 | } |
4667 | sprintf(mat.title, bContact ? "Contact Existence Map" : |
4668 | "Hydrogen Bond Existence Map"); |
4669 | sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds"); |
4670 | sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv)); |
4671 | sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index"); |
4672 | mat.bDiscrete = TRUE1; |
4673 | mat.nmap = 2; |
4674 | snew(mat.map, mat.nmap)(mat.map) = save_calloc("mat.map", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4674, (mat.nmap), sizeof(*(mat.map))); |
4675 | for (i = 0; i < mat.nmap; i++) |
4676 | { |
4677 | mat.map[i].code.c1 = hbmap[i]; |
4678 | mat.map[i].desc = hbdesc[i]; |
4679 | mat.map[i].rgb = hbrgb[i]; |
4680 | } |
4681 | fp = opt2FILE("-hbm", NFILE, fnm, "w")gmx_ffopen(opt2fn("-hbm", ((int)(sizeof(fnm)/sizeof((fnm)[0]) )), fnm), "w"); |
4682 | write_xpm_m(fp, mat); |
4683 | gmx_ffclose(fp); |
4684 | for (x = 0; x < mat.nx; x++) |
4685 | { |
4686 | sfree(mat.matrix[x])save_free("mat.matrix[x]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4686, (mat.matrix[x])); |
4687 | } |
4688 | sfree(mat.axis_y)save_free("mat.axis_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4688, (mat.axis_y)); |
4689 | sfree(mat.matrix)save_free("mat.matrix", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4689, (mat.matrix)); |
4690 | sfree(mat.map)save_free("mat.map", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4690, (mat.map)); |
4691 | } |
4692 | else |
4693 | { |
4694 | fprintf(stderrstderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n"); |
4695 | } |
4696 | } |
4697 | } |
4698 | |
4699 | if (bGem) |
4700 | { |
4701 | fprintf(stderrstderr, "There were %i periodic shifts\n", hb->per->nper); |
4702 | fprintf(stderrstderr, "Freeing pHist for all donors...\n"); |
4703 | for (i = 0; i < hb->d.nrd; i++) |
4704 | { |
4705 | fprintf(stderrstderr, "\r%i", i); |
4706 | if (hb->per->pHist[i] != NULL((void*)0)) |
4707 | { |
4708 | for (j = 0; j < hb->a.nra; j++) |
4709 | { |
4710 | clearPshift(&(hb->per->pHist[i][j])); |
4711 | } |
4712 | sfree(hb->per->pHist[i])save_free("hb->per->pHist[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4712, (hb->per->pHist[i])); |
4713 | } |
4714 | } |
4715 | sfree(hb->per->pHist)save_free("hb->per->pHist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4715, (hb->per->pHist)); |
4716 | sfree(hb->per->p2i)save_free("hb->per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4716, (hb->per->p2i)); |
4717 | sfree(hb->per)save_free("hb->per", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4717, (hb->per)); |
4718 | fprintf(stderrstderr, "...done.\n"); |
4719 | } |
4720 | |
4721 | #ifdef HAVE_NN_LOOPS |
4722 | if (bNN) |
4723 | { |
4724 | free_hbEmap(hb); |
4725 | } |
4726 | #endif |
4727 | |
4728 | if (hb->bDAnr) |
4729 | { |
4730 | int i, j, nleg; |
4731 | char **legnames; |
4732 | char buf[STRLEN4096]; |
4733 | |
4734 | #define USE_THIS_GROUP(j)( (j == gr0) || (bTwo && (j == gr1)) ) ( (j == gr0) || (bTwo && (j == gr1)) ) |
4735 | |
4736 | fp = xvgropen(opt2fn("-dan", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
4737 | "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv); |
4738 | nleg = (bTwo ? 2 : 1)*2; |
4739 | snew(legnames, nleg)(legnames) = save_calloc("legnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4739, (nleg), sizeof(*(legnames))); |
4740 | i = 0; |
4741 | for (j = 0; j < grNR; j++) |
4742 | { |
4743 | if (USE_THIS_GROUP(j)( (j == gr0) || (bTwo && (j == gr1)) ) ) |
4744 | { |
4745 | sprintf(buf, "Donors %s", grpnames[j]); |
4746 | legnames[i++] = 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))); |
4747 | sprintf(buf, "Acceptors %s", grpnames[j]); |
4748 | legnames[i++] = 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))); |
4749 | } |
4750 | } |
4751 | if (i != nleg) |
4752 | { |
4753 | gmx_incons("number of legend entries")_gmx_error("incons", "number of legend entries", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4753); |
4754 | } |
4755 | xvgr_legend(fp, nleg, (const char**)legnames, oenv); |
4756 | for (i = 0; i < nframes; i++) |
4757 | { |
4758 | fprintf(fp, "%10g", hb->time[i]); |
4759 | for (j = 0; (j < grNR); j++) |
4760 | { |
4761 | if (USE_THIS_GROUP(j)( (j == gr0) || (bTwo && (j == gr1)) ) ) |
4762 | { |
4763 | fprintf(fp, " %6d", hb->danr[i][j]); |
4764 | } |
4765 | } |
4766 | fprintf(fp, "\n"); |
4767 | } |
4768 | gmx_ffclose(fp); |
4769 | } |
4770 | |
4771 | return 0; |
4772 | } |