File: | gromacs/gmxana/gmx_hbond.c |
Location: | line 1998, column 25 |
Description: | Array access (via field 'pHist') results in a null pointer dereference |
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, | |||
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 | } |