2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, 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.
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.
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.
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.
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.
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.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/correlationfunctions/autocorr.h"
45 #include "gromacs/correlationfunctions/crosscorr.h"
46 #include "gromacs/correlationfunctions/expfit.h"
47 #include "gromacs/correlationfunctions/integrate.h"
48 #include "gromacs/fileio/matio.h"
49 #include "gromacs/fileio/tpxio.h"
50 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/gmxana/geminate.h"
53 #include "gromacs/gmxana/gmx_ana.h"
54 #include "gromacs/legacyheaders/copyrite.h"
55 #include "gromacs/legacyheaders/macros.h"
56 #include "gromacs/legacyheaders/txtdump.h"
57 #include "gromacs/legacyheaders/viewit.h"
58 #include "gromacs/math/units.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/topology/index.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/gmxomp.h"
66 #include "gromacs/utility/smalloc.h"
67 #include "gromacs/utility/snprintf.h"
71 /*#define HAVE_NN_LOOPS*/
73 typedef short int t_E;
76 typedef int t_hx[max_hx];
77 #define NRHXTYPES max_hx
78 const char *hxtypenames[NRHXTYPES] =
79 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
83 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == 0)
85 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == (threadNr))
88 /* -----------------------------------------*/
94 hbNo, hbDist, hbHB, hbNR, hbR2
97 noDA, ACC, DON, DA, INGROUP
100 NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR
103 static const char *grpnames[grNR] = {"0", "1", "I" };
105 static gmx_bool bDebug = FALSE;
110 #define HB_YESINS HB_YES|HB_INS
114 #define ISHB(h) (((h) & 2) == 2)
115 #define ISDIST(h) (((h) & 1) == 1)
116 #define ISDIST2(h) (((h) & 4) == 4)
117 #define ISACC(h) (((h) & 1) == 1)
118 #define ISDON(h) (((h) & 2) == 2)
119 #define ISINGRP(h) (((h) & 4) == 4)
132 typedef int t_icell[grNR];
133 typedef atom_id h_id[MAXHYDRO];
136 int history[MAXHYDRO];
137 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
138 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
140 /* Bitmask array which tells whether a hbond is present
141 * at a given time. Either of these may be NULL
143 int n0; /* First frame a HB was found */
144 int nframes, maxframes; /* Amount of frames in this hbond */
147 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
148 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
149 * acceptor distance is less than the user-specified distance (typically
156 atom_id *acc; /* Atom numbers of the acceptors */
157 int *grp; /* Group index */
158 int *aptr; /* Map atom number to acceptor index */
163 int *don; /* Atom numbers of the donors */
164 int *grp; /* Group index */
165 int *dptr; /* Map atom number to donor index */
166 int *nhydro; /* Number of hydrogens for each donor */
167 h_id *hydro; /* The atom numbers of the hydrogens */
168 h_id *nhbonds; /* The number of HBs per H at current */
171 /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/
175 int len; /* The length of frame and p. */
176 int *frame; /* The frames at which transitio*/
181 /* Periodicity history. Used for the reversible geminate recombination. */
182 t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap:
183 * pHist[d][a]. We can safely assume that the same
184 * periodic shift holds for all hydrogens of a da-pair.
186 * Nowadays it only stores TRANSITIONS, and not the shift at every frame.
187 * That saves a LOT of memory, an hopefully kills a mysterious bug where
188 * pHist gets contaminated. */
190 PSTYPE nper; /* The length of p2i */
191 ivec *p2i; /* Maps integer to periodic shift for a pair.*/
192 matrix P; /* Projection matrix to find the box shifts. */
193 int gemtype; /* enumerated type */
198 int *Etot; /* Total energy for each frame */
199 t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */
203 gmx_bool bHBmap, bDAnr, bGem;
205 /* The following arrays are nframes long */
206 int nframes, max_frames, maxhydro;
212 /* These structures are initialized from the topology at start up */
215 /* This holds a matrix with all possible hydrogen bonds */
221 /* For parallelization reasons this will have to be a pointer.
222 * Otherwise discrepancies may arise between the periodicity data
223 * seen by different threads. */
227 static void clearPshift(t_pShift *pShift)
232 sfree(pShift->frame);
237 static void calcBoxProjection(matrix B, matrix P)
239 const int vp[] = {XX, YY, ZZ};
244 for (i = 0; i < 3; i++)
247 for (j = 0; j < 3; j++)
250 U[m][n] = i == j ? 1 : 0;
254 for (i = 0; i < 3; i++)
257 mvmul(M, U[m], P[m]);
262 static void calcBoxDistance(matrix P, rvec d, ivec ibd)
264 /* returns integer distance in box coordinates.
265 * P is the projection matrix from cartesian coordinates
266 * obtained with calcBoxProjection(). */
270 /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
271 for (i = 0; i < 3; i++)
273 bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5);
275 ibd[XX] = (int)bd[XX];
276 ibd[YY] = (int)bd[YY];
277 ibd[ZZ] = (int)bd[ZZ];
280 /* Changed argument 'bMerge' into 'oneHB' below,
281 * since -contact should cause maxhydro to be 1,
283 * - Erik Marklund May 29, 2006
286 static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap)
288 /* Try to merge hbonds on the fly. That means that if the
289 * acceptor and donor are mergable, then:
290 * 1) store the hb-info so that acceptor id > donor id,
291 * 2) add the periodic shift in pairs, so that [-x,-y,-z] is
292 * stored in per.p2i[] whenever acceptor id < donor id.
293 * Note that [0,0,0] should already be the first element of per.p2i
294 * by the time this function is called. */
296 /* daSwap is TRUE if the donor and acceptor were swapped.
297 * If so, then the negative vector should be used. */
300 if (per->p2i == NULL || per->nper == 0)
302 gmx_fatal(FARGS, "'per' not initialized properly.");
304 for (i = 0; i < per->nper; i++)
306 if (r[XX] == per->p2i[i][XX] &&
307 r[YY] == per->p2i[i][YY] &&
308 r[ZZ] == per->p2i[i][ZZ])
313 /* Not found apparently. Add it to the list! */
314 /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
320 fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
325 srenew(per->p2i, per->nper+2);
327 copy_ivec(r, per->p2i[per->nper]);
330 /* Add the mirror too. It's rather likely that it'll be needed. */
331 per->p2i[per->nper][XX] = -r[XX];
332 per->p2i[per->nper][YY] = -r[YY];
333 per->p2i[per->nper][ZZ] = -r[ZZ];
336 return per->nper - 1 - (daSwap ? 0 : 1);
339 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode)
344 hb->wordlen = 8*sizeof(unsigned int);
354 hb->maxhydro = MAXHYDRO;
357 hb->per->gemtype = bGem ? gemmode : 0;
362 static void mk_hbmap(t_hbdata *hb)
366 snew(hb->hbmap, hb->d.nrd);
367 for (i = 0; (i < hb->d.nrd); i++)
369 snew(hb->hbmap[i], hb->a.nra);
370 if (hb->hbmap[i] == NULL)
372 gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
374 for (j = 0; (j > hb->a.nra); j++)
376 hb->hbmap[i][j] = NULL;
381 /* Consider redoing pHist so that is only stores transitions between
382 * periodicities and not the periodicity for all frames. This eats heaps of memory. */
383 static void mk_per(t_hbdata *hb)
388 snew(hb->per->pHist, hb->d.nrd);
389 for (i = 0; i < hb->d.nrd; i++)
391 snew(hb->per->pHist[i], hb->a.nra);
392 if (hb->per->pHist[i] == NULL)
394 gmx_fatal(FARGS, "Could not allocate enough memory for per->pHist");
396 for (j = 0; j < hb->a.nra; j++)
398 clearPshift(&(hb->per->pHist[i][j]));
401 /* add the [0,0,0] shift to element 0 of p2i. */
402 snew(hb->per->p2i, 1);
403 clear_ivec(hb->per->p2i[0]);
409 static void mk_hbEmap (t_hbdata *hb, int n0)
414 snew(hb->hbE.E, hb->d.nrd);
415 for (i = 0; i < hb->d.nrd; i++)
417 snew(hb->hbE.E[i], hb->a.nra);
418 for (j = 0; j < hb->a.nra; j++)
420 snew(hb->hbE.E[i][j], MAXHYDRO);
421 for (k = 0; k < MAXHYDRO; k++)
423 hb->hbE.E[i][j][k] = NULL;
430 static void free_hbEmap (t_hbdata *hb)
433 for (i = 0; i < hb->d.nrd; i++)
435 for (j = 0; j < hb->a.nra; j++)
437 for (k = 0; k < MAXHYDRO; k++)
439 sfree(hb->hbE.E[i][j][k]);
441 sfree(hb->hbE.E[i][j]);
449 static void addFramesNN(t_hbdata *hb, int frame)
452 #define DELTAFRAMES_HBE 10
454 int d, a, h, nframes;
456 if (frame >= hb->hbE.nframes)
458 nframes = hb->hbE.nframes + DELTAFRAMES_HBE;
459 srenew(hb->hbE.Etot, nframes);
461 for (d = 0; d < hb->d.nrd; d++)
463 for (a = 0; a < hb->a.nra; a++)
465 for (h = 0; h < hb->d.nhydro[d]; h++)
467 srenew(hb->hbE.E[d][a][h], nframes);
472 hb->hbE.nframes += DELTAFRAMES_HBE;
476 static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst,
477 matrix box, rvec hbox, t_donors *donors)
482 * alpha - angle between dipoles
483 * x[] - atomic positions
484 * EEst - the type of energy estimate (see enum in hbplugin.h)
485 * box - the box vectors \
486 * hbox - half box lengths _These two are only needed for the pbc correction
491 rvec dipole[2], xmol[3], xmean[2];
497 /* Self-interaction */
504 /* This is a simple binary existence function that sets E=1 whenever
505 * the distance between the oxygens is equal too or less than 0.35 nm.
507 rvec_sub(x[d], x[a], dist);
508 pbc_correct_gem(dist, box, hbox);
509 if (norm(dist) <= 0.35)
520 /* Negative potential energy of a dipole.
521 * E = -cos(alpha) * 1/r^3 */
523 copy_rvec(x[d], xmol[0]); /* donor */
524 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
525 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
527 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
528 rvec_inc(xmean[0], xmol[1]);
529 rvec_inc(xmean[0], xmol[2]);
530 for (i = 0; i < 3; i++)
532 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
535 /* Assumes that all acceptors are also donors. */
536 copy_rvec(x[a], xmol[0]); /* acceptor */
537 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
538 copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */
541 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
542 rvec_inc(xmean[1], xmol[1]);
543 rvec_inc(xmean[1], xmol[2]);
544 for (i = 0; i < 3; i++)
546 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
549 rvec_sub(xmean[0], xmean[1], dist);
550 pbc_correct_gem(dist, box, hbox);
553 realE = pow(r, -3.0);
554 E = (t_E)(SCALEFACTOR_E * realE);
558 /* Negative potential energy of a (unpolarizable) dipole.
559 * E = -cos(alpha) * 1/r^3 */
560 clear_rvec(dipole[1]);
561 clear_rvec(dipole[0]);
563 copy_rvec(x[d], xmol[0]); /* donor */
564 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
565 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
567 rvec_inc(dipole[0], xmol[1]);
568 rvec_inc(dipole[0], xmol[2]);
569 for (i = 0; i < 3; i++)
573 rvec_dec(dipole[0], xmol[0]);
575 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
576 rvec_inc(xmean[0], xmol[1]);
577 rvec_inc(xmean[0], xmol[2]);
578 for (i = 0; i < 3; i++)
580 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
583 /* Assumes that all acceptors are also donors. */
584 copy_rvec(x[a], xmol[0]); /* acceptor */
585 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
586 copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */
589 rvec_inc(dipole[1], xmol[1]);
590 rvec_inc(dipole[1], xmol[2]);
591 for (i = 0; i < 3; i++)
595 rvec_dec(dipole[1], xmol[0]);
597 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
598 rvec_inc(xmean[1], xmol[1]);
599 rvec_inc(xmean[1], xmol[2]);
600 for (i = 0; i < 3; i++)
602 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
605 rvec_sub(xmean[0], xmean[1], dist);
606 pbc_correct_gem(dist, box, hbox);
609 double cosalpha = cos_angle(dipole[0], dipole[1]);
610 realE = cosalpha * pow(r, -3.0);
611 E = (t_E)(SCALEFACTOR_E * realE);
615 printf("Can't do that type of energy estimate: %i\n.", EEst);
622 static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame)
624 /* hb - hbond data structure
628 E - estimate of the energy
629 frame - the current frame.
632 /* Store the estimated energy */
638 hb->hbE.E[d][a][h][frame] = E;
642 hb->hbE.Etot[frame] += E;
645 #endif /* HAVE_NN_LOOPS */
648 /* Finds -v[] in the periodicity index */
649 static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
653 for (i = 0; i < nper; i++)
655 if (v[i][XX] == -(v[p][XX]) &&
656 v[i][YY] == -(v[p][YY]) &&
657 v[i][ZZ] == -(v[p][ZZ]))
662 printf("Couldn't find mirror of [%i, %i, %i], index \n",
670 static void add_frames(t_hbdata *hb, int nframes)
674 if (nframes >= hb->max_frames)
676 hb->max_frames += 4096;
677 srenew(hb->time, hb->max_frames);
678 srenew(hb->nhb, hb->max_frames);
679 srenew(hb->ndist, hb->max_frames);
680 srenew(hb->n_bound, hb->max_frames);
681 srenew(hb->nhx, hb->max_frames);
684 srenew(hb->danr, hb->max_frames);
687 hb->nframes = nframes;
690 #define OFFSET(frame) (frame / 32)
691 #define MASK(frame) (1 << (frame % 32))
693 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
697 hbexist[OFFSET(frame)] |= MASK(frame);
701 hbexist[OFFSET(frame)] &= ~MASK(frame);
705 static gmx_bool is_hb(unsigned int hbexist[], int frame)
707 return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
710 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
712 unsigned int *ghptr = NULL;
716 ghptr = hb->hbmap[id][ia]->h[ih];
718 else if (ihb == hbDist)
720 ghptr = hb->hbmap[id][ia]->g[ih];
724 gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
727 _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
730 static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
734 snew(pHist->frame, 1);
737 pHist->frame[0] = frame;
742 if (pHist->p[pHist->len-1] != p)
745 srenew(pHist->frame, pHist->len);
746 srenew(pHist->p, pHist->len);
747 pHist->frame[pHist->len-1] = frame;
748 pHist->p[pHist->len-1] = p;
749 } /* Otherwise, there is no transition. */
753 static PSTYPE getPshift(t_pShift pHist, int frame)
758 || (pHist.len > 0 && pHist.frame[0] > frame))
763 for (i = 0; i < pHist.len; i++)
776 /* It seems that frame is after the last periodic transition. Return the last periodicity. */
777 return pHist.p[pHist.len-1];
780 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPE p)
783 t_hbond *hb = hbd->hbmap[id][ia];
784 int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]);
785 int wlen = hbd->wordlen;
787 gmx_bool bGem = hbd->bGem;
792 hb->maxframes = delta;
793 for (i = 0; (i < maxhydro); i++)
795 snew(hb->h[i], hb->maxframes/wlen);
796 snew(hb->g[i], hb->maxframes/wlen);
801 hb->nframes = frame-hb->n0;
802 /* We need a while loop here because hbonds may be returning
805 while (hb->nframes >= hb->maxframes)
807 n = hb->maxframes + delta;
808 for (i = 0; (i < maxhydro); i++)
810 srenew(hb->h[i], n/wlen);
811 srenew(hb->g[i], n/wlen);
812 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
824 set_hb(hbd, id, h, ia, frame, ihb);
827 if (p >= hbd->per->nper)
829 gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
833 addPshift(&(hbd->per->pHist[id][ia]), p, frame);
841 static void inc_nhbonds(t_donors *ddd, int d, int h)
844 int dptr = ddd->dptr[d];
846 for (j = 0; (j < ddd->nhydro[dptr]); j++)
848 if (ddd->hydro[dptr][j] == h)
850 ddd->nhbonds[dptr][j]++;
854 if (j == ddd->nhydro[dptr])
856 gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
860 static int _acceptor_index(t_acceptors *a, int grp, atom_id i,
861 const char *file, int line)
865 if (a->grp[ai] != grp)
869 fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
870 ai, a->grp[ai], grp, file, line);
879 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
881 static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line)
890 if (d->grp[di] != grp)
894 fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
895 di, d->grp[di], grp, file, line);
904 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
906 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
908 /* g_hbond doesn't allow overlapping groups */
914 donor_index(&hb->d, grpd, a) != NOTSET
915 && acceptor_index(&hb->a, grpa, d) != NOTSET;
919 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
920 int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPE p)
923 gmx_bool daSwap = FALSE;
925 if ((id = hb->d.dptr[d]) == NOTSET)
927 gmx_fatal(FARGS, "No donor atom %d", d+1);
929 else if (grpd != hb->d.grp[id])
931 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
932 grpd, hb->d.grp[id], d+1);
934 if ((ia = hb->a.aptr[a]) == NOTSET)
936 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
938 else if (grpa != hb->a.grp[ia])
940 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
941 grpa, hb->a.grp[ia], a+1);
947 if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
948 /* Then swap identity so that the id of d is lower then that of a.
950 * This should really be redundant by now, as is_hbond() now ought to return
951 * hbNo in the cases where this conditional is TRUE. */
958 /* Now repeat donor/acc check. */
959 if ((id = hb->d.dptr[d]) == NOTSET)
961 gmx_fatal(FARGS, "No donor atom %d", d+1);
963 else if (grpd != hb->d.grp[id])
965 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
966 grpd, hb->d.grp[id], d+1);
968 if ((ia = hb->a.aptr[a]) == NOTSET)
970 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
972 else if (grpa != hb->a.grp[ia])
974 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
975 grpa, hb->a.grp[ia], a+1);
982 /* Loop over hydrogens to find which hydrogen is in this particular HB */
983 if ((ihb == hbHB) && !bMerge && !bContact)
985 for (k = 0; (k < hb->d.nhydro[id]); k++)
987 if (hb->d.hydro[id][k] == h)
992 if (k == hb->d.nhydro[id])
994 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
1006 #pragma omp critical
1008 if (hb->hbmap[id][ia] == NULL)
1010 snew(hb->hbmap[id][ia], 1);
1011 snew(hb->hbmap[id][ia]->h, hb->maxhydro);
1012 snew(hb->hbmap[id][ia]->g, hb->maxhydro);
1014 add_ff(hb, id, k, ia, frame, ihb, p);
1018 /* Strange construction with frame >=0 is a relic from old code
1019 * for selected hbond analysis. It may be necessary again if that
1020 * is made to work again.
1024 hh = hb->hbmap[id][ia]->history[k];
1030 hb->hbmap[id][ia]->history[k] = hh | 2;
1041 hb->hbmap[id][ia]->history[k] = hh | 1;
1065 if (bMerge && daSwap)
1067 h = hb->d.hydro[id][0];
1069 /* Increment number if HBonds per H */
1070 if (ihb == hbHB && !bContact)
1072 inc_nhbonds(&(hb->d), d, h);
1076 static char *mkatomname(t_atoms *atoms, int i)
1078 static char buf[32];
1081 rnr = atoms->atom[i].resind;
1082 sprintf(buf, "%4s%d%-4s",
1083 *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
1088 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms)
1090 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
1093 for (i = 0; i < isize; i++)
1095 if (index[i] >= natoms)
1097 gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
1099 datable[index[i]] |= INGROUP;
1103 static void clear_datable_grp(unsigned char *datable, int size)
1105 /* Clears group information from the table */
1107 const char mask = !(char)INGROUP;
1110 for (i = 0; i < size; i++)
1117 static void add_acc(t_acceptors *a, int ia, int grp)
1119 if (a->nra >= a->max_nra)
1122 srenew(a->acc, a->max_nra);
1123 srenew(a->grp, a->max_nra);
1125 a->grp[a->nra] = grp;
1126 a->acc[a->nra++] = ia;
1129 static void search_acceptors(t_topology *top, int isize,
1130 atom_id *index, t_acceptors *a, int grp,
1132 gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
1138 for (i = 0; (i < isize); i++)
1142 (((*top->atoms.atomname[n])[0] == 'O') ||
1143 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
1144 ISINGRP(datable[n]))
1146 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
1151 snew(a->aptr, top->atoms.nr);
1152 for (i = 0; (i < top->atoms.nr); i++)
1154 a->aptr[i] = NOTSET;
1156 for (i = 0; (i < a->nra); i++)
1158 a->aptr[a->acc[i]] = i;
1162 static void add_h2d(int id, int ih, t_donors *ddd)
1166 for (i = 0; (i < ddd->nhydro[id]); i++)
1168 if (ddd->hydro[id][i] == ih)
1170 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
1175 if (i == ddd->nhydro[id])
1177 if (ddd->nhydro[id] >= MAXHYDRO)
1179 gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
1180 ddd->don[id], MAXHYDRO);
1182 ddd->hydro[id][i] = ih;
1187 static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable)
1191 if (!datable || ISDON(datable[id]))
1193 if (ddd->dptr[id] == NOTSET) /* New donor */
1205 if (ddd->nrd >= ddd->max_nrd)
1207 ddd->max_nrd += 128;
1208 srenew(ddd->don, ddd->max_nrd);
1209 srenew(ddd->nhydro, ddd->max_nrd);
1210 srenew(ddd->hydro, ddd->max_nrd);
1211 srenew(ddd->nhbonds, ddd->max_nrd);
1212 srenew(ddd->grp, ddd->max_nrd);
1214 ddd->don[ddd->nrd] = id;
1215 ddd->nhydro[ddd->nrd] = 0;
1216 ddd->grp[ddd->nrd] = grp;
1223 add_h2d(i, ih, ddd);
1228 printf("Warning: Atom %d is not in the d/a-table!\n", id);
1232 static void search_donors(t_topology *top, int isize, atom_id *index,
1233 t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
1234 unsigned char *datable)
1237 t_functype func_type;
1238 t_ilist *interaction;
1239 atom_id nr1, nr2, nr3;
1244 snew(ddd->dptr, top->atoms.nr);
1245 for (i = 0; (i < top->atoms.nr); i++)
1247 ddd->dptr[i] = NOTSET;
1255 for (i = 0; (i < isize); i++)
1257 datable[index[i]] |= DON;
1258 add_dh(ddd, index[i], -1, grp, datable);
1264 for (func_type = 0; (func_type < F_NRE); func_type++)
1266 interaction = &(top->idef.il[func_type]);
1267 if (func_type == F_POSRES || func_type == F_FBPOSRES)
1269 /* The ilist looks strange for posre. Bug in grompp?
1270 * We don't need posre interactions for hbonds anyway.*/
1273 for (i = 0; i < interaction->nr;
1274 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1277 if (func_type != top->idef.functype[interaction->iatoms[i]])
1279 fprintf(stderr, "Error in func_type %s",
1280 interaction_function[func_type].longname);
1284 /* check out this functype */
1285 if (func_type == F_SETTLE)
1287 nr1 = interaction->iatoms[i+1];
1288 nr2 = interaction->iatoms[i+2];
1289 nr3 = interaction->iatoms[i+3];
1291 if (ISINGRP(datable[nr1]))
1293 if (ISINGRP(datable[nr2]))
1295 datable[nr1] |= DON;
1296 add_dh(ddd, nr1, nr1+1, grp, datable);
1298 if (ISINGRP(datable[nr3]))
1300 datable[nr1] |= DON;
1301 add_dh(ddd, nr1, nr1+2, grp, datable);
1305 else if (IS_CHEMBOND(func_type))
1307 for (j = 0; j < 2; j++)
1309 nr1 = interaction->iatoms[i+1+j];
1310 nr2 = interaction->iatoms[i+2-j];
1311 if ((*top->atoms.atomname[nr1][0] == 'H') &&
1312 ((*top->atoms.atomname[nr2][0] == 'O') ||
1313 (*top->atoms.atomname[nr2][0] == 'N')) &&
1314 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1316 datable[nr2] |= DON;
1317 add_dh(ddd, nr2, nr1, grp, datable);
1324 for (func_type = 0; func_type < F_NRE; func_type++)
1326 interaction = &top->idef.il[func_type];
1327 for (i = 0; i < interaction->nr;
1328 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1331 if (func_type != top->idef.functype[interaction->iatoms[i]])
1333 gmx_incons("function type in search_donors");
1336 if (interaction_function[func_type].flags & IF_VSITE)
1338 nr1 = interaction->iatoms[i+1];
1339 if (*top->atoms.atomname[nr1][0] == 'H')
1343 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1354 if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1355 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1356 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1358 datable[nr2] |= DON;
1359 add_dh(ddd, nr2, nr1, grp, datable);
1369 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
1376 for (i = 0; i < DIM; i++)
1378 ngrid[i] = (box[i][i]/(1.2*rcut));
1382 if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
1384 for (i = 0; i < DIM; i++)
1391 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1392 ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
1394 snew(grid, ngrid[ZZ]);
1395 for (z = 0; z < ngrid[ZZ]; z++)
1397 snew((grid)[z], ngrid[YY]);
1398 for (y = 0; y < ngrid[YY]; y++)
1400 snew((grid)[z][y], ngrid[XX]);
1406 static void reset_nhbonds(t_donors *ddd)
1410 for (i = 0; (i < ddd->nrd); i++)
1412 for (j = 0; (j < MAXHH); j++)
1414 ddd->nhbonds[i][j] = 0;
1419 void pbc_correct_gem(rvec dx, matrix box, rvec hbox);
1420 void pbc_in_gridbox(rvec dx, matrix box);
1422 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
1423 gmx_bool bBox, matrix box, rvec hbox,
1424 real rcut, real rshell,
1425 ivec ngrid, t_gridcell ***grid)
1427 int i, m, gr, xi, yi, zi, nr;
1430 rvec invdelta, dshell, xtemp = {0, 0, 0};
1432 gmx_bool bDoRshell, bInShell, bAcc;
1437 bDoRshell = (rshell > 0);
1438 rshell2 = sqr(rshell);
1441 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
1443 for (m = 0; m < DIM; m++)
1445 hbox[m] = box[m][m]*0.5;
1448 invdelta[m] = ngrid[m]/box[m][m];
1449 if (1/invdelta[m] < rcut)
1451 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
1452 "%s can not handle this situation, sorry.\n",
1465 /* resetting atom counts */
1466 for (gr = 0; (gr < grNR); gr++)
1468 for (zi = 0; zi < ngrid[ZZ]; zi++)
1470 for (yi = 0; yi < ngrid[YY]; yi++)
1472 for (xi = 0; xi < ngrid[XX]; xi++)
1474 grid[zi][yi][xi].d[gr].nr = 0;
1475 grid[zi][yi][xi].a[gr].nr = 0;
1481 /* put atoms in grid cells */
1482 for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
1495 for (i = 0; (i < nr); i++)
1497 /* check if we are inside the shell */
1498 /* if bDoRshell=FALSE then bInShell=TRUE always */
1503 rvec_sub(x[ad[i]], xshell, dshell);
1506 if (FALSE && !hb->bGem)
1508 for (m = DIM-1; m >= 0 && bInShell; m--)
1510 if (dshell[m] < -hbox[m])
1512 rvec_inc(dshell, box[m]);
1514 else if (dshell[m] >= hbox[m])
1516 dshell[m] -= 2*hbox[m];
1518 /* if we're outside the cube, we're outside the sphere also! */
1519 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1527 gmx_bool bDone = FALSE;
1531 for (m = DIM-1; m >= 0 && bInShell; m--)
1533 if (dshell[m] < -hbox[m])
1536 rvec_inc(dshell, box[m]);
1538 if (dshell[m] >= hbox[m])
1541 dshell[m] -= 2*hbox[m];
1545 for (m = DIM-1; m >= 0 && bInShell; m--)
1547 /* if we're outside the cube, we're outside the sphere also! */
1548 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1555 /* if we're inside the cube, check if we're inside the sphere */
1558 bInShell = norm2(dshell) < rshell2;
1568 copy_rvec(x[ad[i]], xtemp);
1570 pbc_in_gridbox(x[ad[i]], box);
1572 for (m = DIM-1; m >= 0; m--)
1573 { /* determine grid index of atom */
1574 grididx[m] = x[ad[i]][m]*invdelta[m];
1575 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1579 copy_rvec(xtemp, x[ad[i]]); /* copy back */
1586 range_check(gx, 0, ngrid[XX]);
1587 range_check(gy, 0, ngrid[YY]);
1588 range_check(gz, 0, ngrid[ZZ]);
1592 /* add atom to grid cell */
1595 newgrid = &(grid[gz][gy][gx].a[gr]);
1599 newgrid = &(grid[gz][gy][gx].d[gr]);
1601 if (newgrid->nr >= newgrid->maxnr)
1603 newgrid->maxnr += 10;
1604 DBB(newgrid->maxnr);
1605 srenew(newgrid->atoms, newgrid->maxnr);
1608 newgrid->atoms[newgrid->nr] = ad[i];
1616 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1620 for (gr = 0; (gr < grNR); gr++)
1623 for (zi = 0; zi < ngrid[ZZ]; zi++)
1625 for (yi = 0; yi < ngrid[YY]; yi++)
1627 for (xi = 0; xi < ngrid[XX]; xi++)
1629 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1637 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1638 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1639 * With a triclinic box all loops are 3 long, except when a cell is
1640 * located next to one of the box edges which is not parallel to the
1641 * x/y-plane, in that case all cells in a line or layer are searched.
1642 * This could be implemented slightly more efficient, but the code
1643 * would get much more complicated.
1645 static gmx_inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1647 return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
1649 static gmx_inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1651 return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1653 static gmx_inline int grid_mod(int j, int n)
1658 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1660 int gr, x, y, z, sum[grNR];
1662 fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1663 for (gr = 0; gr < grNR; gr++)
1666 fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1667 for (z = 0; z < ngrid[ZZ]; z += 2)
1669 fprintf(fp, "Z=%d,%d\n", z, z+1);
1670 for (y = 0; y < ngrid[YY]; y++)
1672 for (x = 0; x < ngrid[XX]; x++)
1674 fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1675 sum[gr] += grid[z][y][x].d[gr].nr;
1676 fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1677 sum[gr] += grid[z][y][x].a[gr].nr;
1681 if ( (z+1) < ngrid[ZZ])
1683 for (x = 0; x < ngrid[XX]; x++)
1685 fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1686 sum[gr] += grid[z+1][y][x].d[gr].nr;
1687 fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1688 sum[gr] += grid[z+1][y][x].a[gr].nr;
1695 fprintf(fp, "TOTALS:");
1696 for (gr = 0; gr < grNR; gr++)
1698 fprintf(fp, " %d=%d", gr, sum[gr]);
1703 /* New GMX record! 5 * in a row. Congratulations!
1704 * Sorry, only four left.
1706 static void free_grid(ivec ngrid, t_gridcell ****grid)
1709 t_gridcell ***g = *grid;
1711 for (z = 0; z < ngrid[ZZ]; z++)
1713 for (y = 0; y < ngrid[YY]; y++)
1723 void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1726 gmx_bool bDone = FALSE;
1730 for (m = DIM-1; m >= 0; m--)
1732 if (dx[m] < -hbox[m])
1735 rvec_inc(dx, box[m]);
1737 if (dx[m] >= hbox[m])
1740 rvec_dec(dx, box[m]);
1746 void pbc_in_gridbox(rvec dx, matrix box)
1749 gmx_bool bDone = FALSE;
1753 for (m = DIM-1; m >= 0; m--)
1758 rvec_inc(dx, box[m]);
1760 if (dx[m] >= box[m][m])
1763 rvec_dec(dx, box[m]);
1769 /* Added argument r2cut, changed contact and implemented
1770 * use of second cut-off.
1771 * - Erik Marklund, June 29, 2006
1773 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1774 real rcut, real r2cut, real ccut,
1775 rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1776 real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1777 gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1779 int h, hh, id, ja, ihb;
1780 rvec r_da, r_ha, r_dh, r = {0, 0, 0};
1782 real rc2, r2c2, rda2, rha2, ca;
1783 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1784 gmx_bool daSwap = FALSE;
1791 if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1792 ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
1800 rvec_sub(x[d], x[a], r_da);
1801 /* Insert projection code here */
1803 if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1805 /* Then this hbond/contact will be found again, or it has already been found. */
1810 if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1811 { /* return hbNo; */
1812 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1816 copy_rvec(r_da, r); /* Save this for later */
1817 pbc_correct_gem(r_da, box, hbox);
1821 pbc_correct_gem(r_da, box, hbox);
1824 rda2 = iprod(r_da, r_da);
1828 if (daSwap && grpa == grpd)
1836 calcBoxDistance(hb->per->P, r, ri);
1837 *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */
1841 else if (rda2 < r2c2)
1852 if (bDA && (rda2 > rc2))
1857 for (h = 0; (h < hb->d.nhydro[id]); h++)
1859 hh = hb->d.hydro[id][h];
1863 rvec_sub(x[hh], x[a], r_ha);
1866 pbc_correct_gem(r_ha, box, hbox);
1868 rha2 = iprod(r_ha, r_ha);
1873 calcBoxDistance(hb->per->P, r, ri);
1874 *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */
1877 if (bDA || (!bDA && (rha2 <= rc2)))
1879 rvec_sub(x[d], x[hh], r_dh);
1882 pbc_correct_gem(r_dh, box, hbox);
1889 ca = cos_angle(r_dh, r_da);
1890 /* if angle is smaller, cos is larger */
1894 *d_ha = sqrt(bDA ? rda2 : rha2);
1900 if (bDA || (!bDA && HAinrange))
1910 /* Fixed previously undiscovered bug in the merge
1911 code, where the last frame of each hbond disappears.
1912 - Erik Marklund, June 1, 2006 */
1913 /* Added the following arguments:
1914 * ptmp[] - temporary periodicity hisory
1915 * a1 - identity of first acceptor/donor
1916 * a2 - identity of second acceptor/donor
1917 * - Erik Marklund, FEB 20 2010 */
1919 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1920 * Will do some more testing before removing the function entirely.
1921 * - Erik Marklund, MAY 10 2010 */
1922 static void do_merge(t_hbdata *hb, int ntmp,
1923 unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[],
1924 t_hbond *hb0, t_hbond *hb1, int a1, int a2)
1926 /* Here we need to make sure we're treating periodicity in
1927 * the right way for the geminate recombination kinetics. */
1929 int m, mm, n00, n01, nn0, nnframes;
1933 /* Decide where to start from when merging */
1936 nn0 = min(n00, n01);
1937 nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1938 /* Initiate tmp arrays */
1939 for (m = 0; (m < ntmp); m++)
1945 /* Fill tmp arrays with values due to first HB */
1946 /* Once again '<' had to be replaced with '<='
1947 to catch the last frame in which the hbond
1949 - Erik Marklund, June 1, 2006 */
1950 for (m = 0; (m <= hb0->nframes); m++)
1953 htmp[mm] = is_hb(hb0->h[0], m);
1956 pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1957 if (pm > hb->per->nper)
1959 gmx_fatal(FARGS, "Illegal shift!");
1963 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1967 /* If we're doing geminate recompbination we usually don't need the distances.
1968 * Let's save some memory and time. */
1969 if (TRUE || !hb->bGem || hb->per->gemtype == gemAD)
1971 for (m = 0; (m <= hb0->nframes); m++)
1974 gtmp[mm] = is_hb(hb0->g[0], m);
1978 for (m = 0; (m <= hb1->nframes); m++)
1981 htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1982 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1983 if (hb->bGem /* && ptmp[mm] != 0 */)
1986 /* If this hbond has been seen before with donor and acceptor swapped,
1987 * then we need to find the mirrored (*-1) periodicity vector to truely
1988 * merge the hbond history. */
1989 pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper);
1990 /* Store index of mirror */
1991 if (pm > hb->per->nper)
1993 gmx_fatal(FARGS, "Illegal shift!");
1998 /* Reallocate target array */
1999 if (nnframes > hb0->maxframes)
2001 srenew(hb0->h[0], 4+nnframes/hb->wordlen);
2002 srenew(hb0->g[0], 4+nnframes/hb->wordlen);
2004 if (NULL != hb->per->pHist)
2006 clearPshift(&(hb->per->pHist[a1][a2]));
2009 /* Copy temp array to target array */
2010 for (m = 0; (m <= nnframes); m++)
2012 _set_hb(hb0->h[0], m, htmp[m]);
2013 _set_hb(hb0->g[0], m, gtmp[m]);
2016 addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
2020 /* Set scalar variables */
2022 hb0->maxframes = nnframes;
2025 /* Added argument bContact for nicer output.
2026 * Erik Marklund, June 29, 2006
2028 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
2030 int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp;
2031 unsigned int *htmp, *gtmp;
2036 indnew = hb->nrdist;
2038 /* Check whether donors are also acceptors */
2039 printf("Merging hbonds with Acceptor and Donor swapped\n");
2041 ntmp = 2*hb->max_frames;
2045 for (i = 0; (i < hb->d.nrd); i++)
2047 fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
2049 ii = hb->a.aptr[id];
2050 for (j = 0; (j < hb->a.nra); j++)
2053 jj = hb->d.dptr[ia];
2054 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
2055 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j]))))
2057 hb0 = hb->hbmap[i][j];
2058 hb1 = hb->hbmap[jj][ii];
2059 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
2061 do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j);
2062 if (ISHB(hb1->history[0]))
2066 else if (ISDIST(hb1->history[0]))
2073 gmx_incons("No contact history");
2077 gmx_incons("Neither hydrogen bond nor distance");
2083 clearPshift(&(hb->per->pHist[jj][ii]));
2087 hb1->history[0] = hbNo;
2092 fprintf(stderr, "\n");
2093 printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
2094 printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
2096 hb->nrdist = indnew;
2102 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
2104 int i, j, k, n_bound[MAXHH], nbtot;
2108 /* Set array to 0 */
2109 for (k = 0; (k < MAXHH); k++)
2113 /* Loop over possible donors */
2114 for (i = 0; (i < hb->d.nrd); i++)
2116 for (j = 0; (j < hb->d.nhydro[i]); j++)
2118 n_bound[hb->d.nhbonds[i][j]]++;
2121 fprintf(fp, "%12.5e", t);
2123 for (k = 0; (k < MAXHH); k++)
2125 fprintf(fp, " %8d", n_bound[k]);
2126 nbtot += n_bound[k]*k;
2128 fprintf(fp, " %8d\n", nbtot);
2131 /* Added argument bContact in do_hblife(...). Also
2132 * added support for -contact in function body.
2133 * - Erik Marklund, May 31, 2006 */
2134 /* Changed the contact code slightly.
2135 * - Erik Marklund, June 29, 2006
2137 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
2138 const output_env_t oenv)
2141 const char *leg[] = { "p(t)", "t p(t)" };
2143 int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
2144 int nframes = hb->nframes;
2147 double sum, integral;
2150 snew(h, hb->maxhydro);
2151 snew(histo, nframes+1);
2152 /* Total number of hbonds analyzed here */
2153 for (i = 0; (i < hb->d.nrd); i++)
2155 for (k = 0; (k < hb->a.nra); k++)
2157 hbh = hb->hbmap[i][k];
2175 for (m = 0; (m < hb->maxhydro); m++)
2179 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
2183 for (nh = 0; (nh < nhydro); nh++)
2188 /* Changed '<' into '<=' below, just like I
2189 did in the hbm-output-loop in the main code.
2190 - Erik Marklund, May 31, 2006
2192 for (j = 0; (j <= hbh->nframes); j++)
2194 ihb = is_hb(h[nh], j);
2195 if (debug && (ndump < 10))
2197 fprintf(debug, "%5d %5d\n", j, ihb);
2217 fprintf(stderr, "\n");
2220 fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
2224 fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
2228 xvgr_legend(fp, asize(leg), leg, oenv);
2230 while ((j0 > 0) && (histo[j0] == 0))
2235 for (i = 0; (i <= j0); i++)
2239 dt = hb->time[1]-hb->time[0];
2242 for (i = 1; (i <= j0); i++)
2244 t = hb->time[i] - hb->time[0] - 0.5*dt;
2245 x1 = t*histo[i]/sum;
2246 fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1);
2251 printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
2252 printf("Note that the lifetime obtained in this manner is close to useless\n");
2253 printf("Use the -ac option instead and check the Forward lifetime\n");
2254 please_cite(stdout, "Spoel2006b");
2259 /* Changed argument bMerge into oneHB to handle contacts properly.
2260 * - Erik Marklund, June 29, 2006
2262 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
2265 int i, j, k, m, nd, ihb, idist;
2266 int nframes = hb->nframes;
2274 fp = gmx_ffopen("debug-ac.xvg", "w");
2275 for (j = 0; (j < nframes); j++)
2277 fprintf(fp, "%10.3f", hb->time[j]);
2278 for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
2280 for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
2284 hbh = hb->hbmap[i][k];
2289 ihb = is_hb(hbh->h[0], j);
2290 idist = is_hb(hbh->g[0], j);
2296 for (m = 0; (m < hb->maxhydro) && !ihb; m++)
2298 ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m], j));
2299 idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
2301 /* This is not correct! */
2302 /* What isn't correct? -Erik M */
2307 fprintf(fp, " %1d-%1d", ihb, idist);
2317 static real calc_dg(real tau, real temp)
2328 return kbt*log(kbt*tau/PLANCK);
2333 int n0, n1, nparams, ndelta;
2335 real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
2338 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
2339 real kt[], real sigma_ct[], real sigma_nt[],
2340 real sigma_kt[], real *k, real *kp,
2341 real *sigma_k, real *sigma_kp,
2347 real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
2352 for (i = 0; (i < n); i++)
2354 if (t[i] >= fit_start)
2367 tl.sigma_ct = sigma_ct;
2368 tl.sigma_nt = sigma_nt;
2369 tl.sigma_kt = sigma_kt;
2373 chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
2375 *kp = tl.kkk[1] = *kp;
2377 for (j = 0; (j < NK); j++)
2379 /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
2382 kk2 += sqr(tl.kkk[0]);
2383 kp2 += sqr(tl.kkk[1]);
2386 *sigma_k = sqrt(kk2/NK - sqr(kkk/NK));
2387 *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2392 static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start,
2393 const output_env_t oenv)
2401 for (i = 0; (i < n); i++)
2417 do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(),
2418 effnEXP2, fitparm, 0, NULL);
2421 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
2422 real sigma_ct[], real sigma_nt[], real sigma_kt[],
2423 real fit_start, real temp, real smooth_tail_start,
2424 const output_env_t oenv)
2427 real k = 1, kp = 1, kow = 1;
2428 real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
2429 double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
2430 gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2432 if (smooth_tail_start >= 0)
2434 smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv);
2435 smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv);
2436 smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv);
2438 for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2444 for (i = i0; (i < n); i++)
2453 printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
2454 tmp = (sn2*sc2-sqr(scn));
2455 if ((tmp > 0) && (sn2 > 0))
2457 k = (sn2*sck-scn*snk)/tmp;
2458 kp = (k*scn-snk)/sn2;
2462 for (i = i0; (i < n); i++)
2464 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2466 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
2468 &sigma_k, &sigma_kp, fit_start);
2469 Q = 0; /* quality_of_fit(chi2, 2);*/
2470 ddg = BOLTZ*temp*sigma_k/k;
2471 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2473 printf("The Rate and Delta G are followed by an error estimate\n");
2474 printf("----------------------------------------------------------\n"
2475 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
2476 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2477 k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
2478 ddg = BOLTZ*temp*sigma_kp/kp;
2479 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2480 kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
2485 for (i = i0; (i < n); i++)
2487 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2489 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2491 printf("--------------------------------------------------\n"
2492 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
2493 printf("Forward %10.3f %8.3f %10.3f %10g\n",
2494 k, 1/k, calc_dg(1/k, temp), chi2);
2495 printf("Backward %10.3f %8.3f %10.3f\n",
2496 kp, 1/kp, calc_dg(1/kp, temp));
2502 printf("One-way %10.3f %s%8.3f %10.3f\n",
2503 kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp));
2507 printf(" - Numerical problems computing HB thermodynamics:\n"
2508 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
2509 sc2, sn2, sk2, sck, snk, scn);
2511 /* Determine integral of the correlation function */
2512 tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
2513 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
2514 bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
2516 for (i = 0; (i < n-2); i++)
2518 if ((ct[i] > e_1) && (ct[i+1] <= e_1))
2525 /* Determine tau_relax from linear interpolation */
2526 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2527 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx,
2528 tau_rlx, bError ? " " : "",
2529 calc_dg(tau_rlx, temp));
2534 printf("Correlation functions too short to compute thermodynamics\n");
2538 void compute_derivative(int nn, real x[], real y[], real dydx[])
2542 /* Compute k(t) = dc(t)/dt */
2543 for (j = 1; (j < nn-1); j++)
2545 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2547 /* Extrapolate endpoints */
2548 dydx[0] = 2*dydx[1] - dydx[2];
2549 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2552 static void parallel_print(int *data, int nThreads)
2554 /* This prints the donors on which each tread is currently working. */
2557 fprintf(stderr, "\r");
2558 for (i = 0; i < nThreads; i++)
2560 fprintf(stderr, "%-7i", data[i]);
2564 static void normalizeACF(real *ct, real *gt, int nhb, int len)
2566 real ct_fac, gt_fac = 0;
2569 /* Xu and Berne use the same normalization constant */
2574 gt_fac = 1.0/(real)nhb;
2577 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
2578 for (i = 0; i < len; i++)
2588 /* Added argument bContact in do_hbac(...). Also
2589 * added support for -contact in the actual code.
2590 * - Erik Marklund, May 31, 2006 */
2591 /* Changed contact code and added argument R2
2592 * - Erik Marklund, June 29, 2006
2594 static void do_hbac(const char *fn, t_hbdata *hb,
2595 int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
2596 real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv,
2597 const char *gemType, int nThreads,
2598 const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2601 int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
2602 const char *legNN[] = {
2606 static char **legGem;
2608 const char *legLuzar[] = {
2609 "Ac\\sfin sys\\v{}\\z{}(t)",
2611 "Cc\\scontact,hb\\v{}\\z{}(t)",
2612 "-dAc\\sfs\\v{}\\z{}/dt"
2614 gmx_bool bNorm = FALSE, bOMP = FALSE;
2617 real *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
2618 real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
2619 const real tol = 1e-3;
2620 int nframes = hb->nframes, nf;
2621 unsigned int **h = NULL, **g = NULL;
2622 int nh, nhbonds, nhydro, ngh;
2624 PSTYPE p, *pfound = NULL, np;
2626 int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
2627 real **rHbExGem = NULL;
2631 double *ctdouble, *timedouble, *fittedct;
2632 double fittolerance = 0.1;
2633 int *dondata = NULL, thisThread;
2636 AC_NONE, AC_NN, AC_GEM, AC_LUZAR
2645 printf("Doing autocorrelation ");
2647 /* Decide what kind of ACF calculations to do. */
2648 if (NN > NN_NONE && NN < NN_NR)
2650 #ifdef HAVE_NN_LOOPS
2652 printf("using the energy estimate.\n");
2655 printf("Can't do the NN-loop. Yet.\n");
2661 printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2663 nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0);
2664 snew(legGem, nSets);
2665 for (i = 0; i < nSets; i++)
2667 snew(legGem[i], 128);
2669 sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2672 sprintf(legGem[1], "Ac'(t)");
2676 sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2683 printf("according to the theory of Luzar and Chandler.\n");
2687 /* build hbexist matrix in reals for autocorr */
2688 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2690 while (n2 < nframes)
2697 if (acType != AC_NN || bOMP)
2699 snew(h, hb->maxhydro);
2700 snew(g, hb->maxhydro);
2703 /* Dump hbonds for debugging */
2704 dump_ac(hb, bMerge || bContact, nDump);
2706 /* Total number of hbonds analyzed here */
2711 if (acType != AC_LUZAR && bOMP)
2713 nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2715 gmx_omp_set_num_threads(nThreads);
2716 snew(dondata, nThreads);
2717 for (i = 0; i < nThreads; i++)
2721 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2722 "Expect close to linear scaling over this donor-loop.\n", nThreads);
2724 fprintf(stderr, "Donors: [thread no]\n");
2727 for (i = 0; i < nThreads; i++)
2729 snprintf(tmpstr, 7, "[%i]", i);
2730 fprintf(stderr, "%-7s", tmpstr);
2733 fprintf(stderr, "\n");
2737 /* Build the ACF according to acType */
2742 #ifdef HAVE_NN_LOOPS
2743 /* Here we're using the estimated energy for the hydrogen bonds. */
2746 #pragma omp parallel \
2747 private(i, j, k, nh, E, rhbex, thisThread) \
2751 thisThread = gmx_omp_get_thread_num();
2755 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2758 #pragma omp for schedule (dynamic)
2759 for (i = 0; i < hb->d.nrd; i++) /* loop over donors */
2763 #pragma omp critical
2765 dondata[thisThread] = i;
2766 parallel_print(dondata, nThreads);
2771 fprintf(stderr, "\r %i", i);
2774 for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */
2776 for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2778 E = hb->hbE.E[i][j][nh];
2781 for (k = 0; k < nframes; k++)
2783 if (E[k] != NONSENSE_E)
2785 rhbex[k] = (real)E[k];
2789 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0],
2790 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
2791 #pragma omp critical
2793 for (k = 0; (k < nn); k++)
2810 normalizeACF(ct, NULL, 0, nn);
2812 snew(timedouble, nn);
2813 for (j = 0; j < nn; j++)
2815 timedouble[j] = (double)(hb->time[j]);
2816 ctdouble[j] = (double)(ct[j]);
2819 /* Remove ballistic term */
2820 /* Ballistic component removal and fitting to the reversible geminate recombination model
2821 * will be taken out for the time being. First of all, one can remove the ballistic
2822 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2823 * problems with the robustness of the fitting to the model. More work is needed.
2824 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2825 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2826 * a BSD-licence that can do the job.
2828 * - Erik Marklund, June 18 2010.
2830 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2831 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2833 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2834 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2836 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)");
2837 xvgr_legend(fp, asize(legNN), legNN);
2839 for (j = 0; (j < nn); j++)
2841 fprintf(fp, "%10g %10g %10g\n",
2842 hb->time[j]-hb->time[0],
2850 #endif /* HAVE_NN_LOOPS */
2851 break; /* case AC_NN */
2855 memset(ct, 0, 2*n2*sizeof(real));
2857 fprintf(stderr, "Donor:\n");
2860 #define __ACDATA p_ct
2863 #pragma omp parallel \
2864 private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2865 pfound, poff, rHbExGem, p, ihb, mMax, \
2868 { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
2871 thisThread = gmx_omp_get_thread_num();
2872 snew(h, hb->maxhydro);
2873 snew(g, hb->maxhydro);
2880 memset(p_ct, 0, 2*n2*sizeof(real));
2882 /* I'm using a chunk size of 1, since I expect \
2883 * the overhead to be really small compared \
2884 * to the actual calculations \ */
2885 #pragma omp for schedule(dynamic,1) nowait
2886 for (i = 0; i < hb->d.nrd; i++)
2891 #pragma omp critical
2893 dondata[thisThread] = i;
2894 parallel_print(dondata, nThreads);
2899 fprintf(stderr, "\r %i", i);
2901 for (k = 0; k < hb->a.nra; k++)
2903 for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++)
2905 hbh = hb->hbmap[i][k];
2908 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
2909 * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
2910 * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
2911 pHist = &(hb->per->pHist[i][k]);
2912 if (ISHB(hbh->history[nh]) && pHist->len != 0)
2917 g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL;
2921 /* count the number of periodic shifts encountered and store
2922 * them in separate arrays. */
2924 for (j = 0; j < pHist->len; j++)
2927 for (m = 0; m <= np; m++)
2929 if (m == np) /* p not recognized in list. Add it and set up new array. */
2932 if (np > hb->per->nper)
2934 gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
2936 if (m >= mMax) /* Extend the arrays.
2937 * Doing it like this, using mMax to keep track of the sizes,
2938 * eleviates the need for freeing and re-allocating the arrays
2939 * when taking on the next donor-acceptor pair */
2942 srenew(pfound, np); /* The list of found periodic shifts. */
2943 srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */
2944 snew(rHbExGem[m], 2*n2);
2949 if (rHbExGem != NULL && rHbExGem[m] != NULL)
2951 /* This must be done, as this array was most likey
2952 * used to store stuff in some previous iteration. */
2953 memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
2957 fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
2969 } /* m: Loop over found shifts */
2970 } /* j: Loop over shifts */
2972 /* Now unpack and disentangle the existence funtions. */
2973 for (j = 0; j < nf; j++)
2980 * pfound: list of periodic shifts found for this pair.
2981 * poff: list of frame offsets; that is, the first
2982 * frame a hbond has a particular periodic shift. */
2983 p = getPshift(*pHist, j+n0);
2986 for (m = 0; m < np; m++)
2994 gmx_fatal(FARGS, "Shift not found, but must be there.");
2998 ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j));
3003 poff[m] = j; /* Here's where the first hbond with shift p is,
3004 * relative to the start of h[0].*/
3008 gmx_fatal(FARGS, "j<poff[m]");
3010 rHbExGem[m][j-poff[m]] += 1;
3015 /* Now, build ac. */
3016 for (m = 0; m < np; m++)
3018 if (rHbExGem[m][0] > 0 && n0+poff[m] < nn /* && m==0 */)
3020 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0],
3021 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3022 for (j = 0; (j < nn); j++)
3024 __ACDATA[j] += rHbExGem[m][j];
3027 } /* Building of ac. */
3030 } /* hydrogen loop */
3031 } /* acceptor loop */
3034 for (m = 0; m <= mMax; m++)
3047 #pragma omp critical
3049 for (i = 0; i < nn; i++)
3057 } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
3063 normalizeACF(ct, NULL, 0, nn);
3065 fprintf(stderr, "\n\nACF successfully calculated.\n");
3067 /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
3070 snew(timedouble, nn);
3073 for (j = 0; j < nn; j++)
3075 timedouble[j] = (double)(hb->time[j]);
3076 ctdouble[j] = (double)(ct[j]);
3079 /* Remove ballistic term */
3080 /* Ballistic component removal and fitting to the reversible geminate recombination model
3081 * will be taken out for the time being. First of all, one can remove the ballistic
3082 * component with g_analyze afterwards. Secondly, and more importantly, there are still
3083 * problems with the robustness of the fitting to the model. More work is needed.
3084 * A third reason is that we're currently using gsl for this and wish to reduce dependence
3085 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
3086 * a BSD-licence that can do the job.
3088 * - Erik Marklund, June 18 2010.
3090 /* if (bBallistic) { */
3091 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
3092 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
3094 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
3095 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
3098 /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
3103 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3107 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3109 xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv);
3111 for (j = 0; (j < nn); j++)
3113 fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]);
3116 fprintf(fp, " %10g", ctdouble[j]);
3120 fprintf(fp, " %10g", fittedct[j]);
3131 break; /* case AC_GEM */
3144 for (i = 0; (i < hb->d.nrd); i++)
3146 for (k = 0; (k < hb->a.nra); k++)
3149 hbh = hb->hbmap[i][k];
3153 if (bMerge || bContact)
3155 if (ISHB(hbh->history[0]))
3164 for (m = 0; (m < hb->maxhydro); m++)
3166 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
3168 g[nhydro] = hbh->g[m];
3169 h[nhydro] = hbh->h[m];
3176 for (nh = 0; (nh < nhydro); nh++)
3178 int nrint = bContact ? hb->nrdist : hb->nrhb;
3179 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
3181 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
3184 for (j = 0; (j < nframes); j++)
3186 /* Changed '<' into '<=' below, just like I did in
3187 the hbm-output-loop in the gmx_hbond() block.
3188 - Erik Marklund, May 31, 2006 */
3191 ihb = is_hb(h[nh], j);
3192 idist = is_hb(g[nh], j);
3199 /* For contacts: if a second cut-off is provided, use it,
3200 * otherwise use g(t) = 1-h(t) */
3201 if (!R2 && bContact)
3207 gt[j] = idist*(1-ihb);
3214 /* The autocorrelation function is normalized after summation only */
3215 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
3216 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3218 /* Cross correlation analysis for thermodynamics */
3219 for (j = nframes; (j < n2); j++)
3225 cross_corr(n2, ht, gt, dght);
3227 for (j = 0; (j < nn); j++)
3236 fprintf(stderr, "\n");
3239 normalizeACF(ct, ght, nhb, nn);
3241 /* Determine tail value for statistics */
3244 for (j = nn/2; (j < nn); j++)
3247 tail2 += ct[j]*ct[j];
3249 tail /= (nn - nn/2);
3250 tail2 /= (nn - nn/2);
3251 dtail = sqrt(tail2-tail*tail);
3253 /* Check whether the ACF is long enough */
3256 printf("\nWARNING: Correlation function is probably not long enough\n"
3257 "because the standard deviation in the tail of C(t) > %g\n"
3258 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
3261 for (j = 0; (j < nn); j++)
3264 ct[j] = (cct[j]-tail)/(1-tail);
3266 /* Compute negative derivative k(t) = -dc(t)/dt */
3267 compute_derivative(nn, hb->time, ct, kt);
3268 for (j = 0; (j < nn); j++)
3276 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3280 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3282 xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
3285 for (j = 0; (j < nn); j++)
3287 fprintf(fp, "%10g %10g %10g %10g %10g\n",
3288 hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
3292 analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
3293 fit_start, temp, smooth_tail_start, oenv);
3295 do_view(oenv, fn, NULL);
3307 break; /* case AC_LUZAR */
3310 gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
3311 } /* switch (acType) */
3314 static void init_hbframe(t_hbdata *hb, int nframes, real t)
3318 hb->time[nframes] = t;
3319 hb->nhb[nframes] = 0;
3320 hb->ndist[nframes] = 0;
3321 for (i = 0; (i < max_hx); i++)
3323 hb->nhx[nframes][i] = 0;
3325 /* Loop invalidated */
3326 if (hb->bHBmap && 0)
3328 for (i = 0; (i < hb->d.nrd); i++)
3330 for (j = 0; (j < hb->a.nra); j++)
3332 for (m = 0; (m < hb->maxhydro); m++)
3334 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
3336 set_hb(hb, i, m, j, nframes, HB_NO);
3342 /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
3345 static FILE *open_donor_properties_file(const char *fn,
3347 const output_env_t oenv)
3350 const char *leg[] = { "Nbound", "Nfree" };
3357 fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3358 xvgr_legend(fp, asize(leg), leg, oenv);
3363 static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
3365 int i, j, k, nbound, nb, nhtot;
3373 for (i = 0; (i < hb->d.nrd); i++)
3375 for (k = 0; (k < hb->d.nhydro[i]); k++)
3379 for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
3381 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
3382 is_hb(hb->hbmap[i][j]->h[k], nframes))
3390 fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound);
3393 static void dump_hbmap(t_hbdata *hb,
3394 int nfile, t_filenm fnm[], gmx_bool bTwo,
3395 gmx_bool bContact, int isize[], int *index[], char *grpnames[],
3399 int ddd, hhh, aaa, i, j, k, m, grp;
3400 char ds[32], hs[32], as[32];
3403 fp = opt2FILE("-hbn", nfile, fnm, "w");
3404 if (opt2bSet("-g", nfile, fnm))
3406 fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
3407 fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
3413 for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
3415 fprintf(fp, "[ %s ]", grpnames[grp]);
3416 for (i = 0; i < isize[grp]; i++)
3418 fprintf(fp, (i%15) ? " " : "\n");
3419 fprintf(fp, " %4d", index[grp][i]+1);
3423 Added -contact support below.
3424 - Erik Marklund, May 29, 2006
3428 fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
3429 for (i = 0; (i < hb->d.nrd); i++)
3431 if (hb->d.grp[i] == grp)
3433 for (j = 0; (j < hb->d.nhydro[i]); j++)
3435 fprintf(fp, " %4d %4d", hb->d.don[i]+1,
3436 hb->d.hydro[i][j]+1);
3442 fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
3443 for (i = 0; (i < hb->a.nra); i++)
3445 if (hb->a.grp[i] == grp)
3447 fprintf(fp, (i%15 && !first) ? " " : "\n");
3448 fprintf(fp, " %4d", hb->a.acc[i]+1);
3457 fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
3458 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
3462 fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
3465 for (i = 0; (i < hb->d.nrd); i++)
3468 for (k = 0; (k < hb->a.nra); k++)
3471 for (m = 0; (m < hb->d.nhydro[i]); m++)
3473 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
3475 sprintf(ds, "%s", mkatomname(atoms, ddd));
3476 sprintf(as, "%s", mkatomname(atoms, aaa));
3479 fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
3482 fprintf(fplog, "%12s %12s\n", ds, as);
3487 hhh = hb->d.hydro[i][m];
3488 sprintf(hs, "%s", mkatomname(atoms, hhh));
3489 fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
3492 fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
3506 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3507 * It mimics add_frames() and init_frame() to some extent. */
3508 static void sync_hbdata(t_hbdata *p_hb, int nframes)
3511 if (nframes >= p_hb->max_frames)
3513 p_hb->max_frames += 4096;
3514 srenew(p_hb->nhb, p_hb->max_frames);
3515 srenew(p_hb->ndist, p_hb->max_frames);
3516 srenew(p_hb->n_bound, p_hb->max_frames);
3517 srenew(p_hb->nhx, p_hb->max_frames);
3520 srenew(p_hb->danr, p_hb->max_frames);
3522 memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3523 memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3524 p_hb->nhb[nframes] = 0;
3525 p_hb->ndist[nframes] = 0;
3528 p_hb->nframes = nframes;
3531 /* p_hb->nhx[nframes][i] */
3533 memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
3535 /* hb->per will remain constant througout the frame loop,
3536 * even though the data its members point to will change,
3537 * hence no need for re-syncing. */
3540 int gmx_hbond(int argc, char *argv[])
3542 const char *desc[] = {
3543 "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3544 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
3545 "(zero is extended) and the distance Donor - Acceptor",
3546 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
3547 "OH and NH groups are regarded as donors, O is an acceptor always,",
3548 "N is an acceptor by default, but this can be switched using",
3549 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3550 "to the first preceding non-hydrogen atom.[PAR]",
3552 "You need to specify two groups for analysis, which must be either",
3553 "identical or non-overlapping. All hydrogen bonds between the two",
3554 "groups are analyzed.[PAR]",
3556 "If you set [TT]-shell[tt], you will be asked for an additional index group",
3557 "which should contain exactly one atom. In this case, only hydrogen",
3558 "bonds between atoms within the shell distance from the one atom are",
3561 "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3562 "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3563 "If contact kinetics are analyzed by using the -contact option, then",
3564 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3565 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3566 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3567 "See mentioned literature for more details and definitions."
3570 /* "It is also possible to analyse specific hydrogen bonds with",
3571 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3572 "Donor Hydrogen Acceptor, in the following way::",
3579 "Note that the triplets need not be on separate lines.",
3580 "Each atom triplet specifies a hydrogen bond to be analyzed,",
3581 "note also that no check is made for the types of atoms.[PAR]",
3586 " * [TT]-num[tt]: number of hydrogen bonds as a function of time.",
3587 " * [TT]-ac[tt]: average over all autocorrelations of the existence",
3588 " functions (either 0 or 1) of all hydrogen bonds.",
3589 " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
3590 " * [TT]-ang[tt]: angle distribution of all hydrogen bonds.",
3591 " * [TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
3592 " where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3593 " This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3594 " with helices in proteins.",
3595 " * [TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
3596 " for selected groups, all hydrogen bonded atoms from all groups and",
3597 " all solvent atoms involved in insertion.",
3598 " * [TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
3599 " frames, this also contains information on solvent insertion",
3600 " into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3602 " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3603 " each timeframe. This is especially useful when using [TT]-shell[tt].",
3604 " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3605 " compare results to Raman Spectroscopy.",
3607 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3608 "require an amount of memory proportional to the total numbers of donors",
3609 "times the total number of acceptors in the selected group(s)."
3612 static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
3613 static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1;
3614 static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
3615 static int nDump = 0, nFitPoints = 100;
3616 static int nThreads = 0, nBalExp = 4;
3618 static gmx_bool bContact = FALSE, bBallistic = FALSE, bGemFit = FALSE;
3619 static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
3620 static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3624 { "-a", FALSE, etREAL, {&acut},
3625 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
3626 { "-r", FALSE, etREAL, {&rcut},
3627 "Cutoff radius (nm, X - Acceptor, see next option)" },
3628 { "-da", FALSE, etBOOL, {&bDA},
3629 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3630 { "-r2", FALSE, etREAL, {&r2cut},
3631 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3632 { "-abin", FALSE, etREAL, {&abin},
3633 "Binwidth angle distribution (degrees)" },
3634 { "-rbin", FALSE, etREAL, {&rbin},
3635 "Binwidth distance distribution (nm)" },
3636 { "-nitacc", FALSE, etBOOL, {&bNitAcc},
3637 "Regard nitrogen atoms as acceptors" },
3638 { "-contact", FALSE, etBOOL, {&bContact},
3639 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3640 { "-shell", FALSE, etREAL, {&rshell},
3641 "when > 0, only calculate hydrogen bonds within # nm shell around "
3643 { "-fitstart", FALSE, etREAL, {&fit_start},
3644 "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]" },
3645 { "-fitend", FALSE, etREAL, {&fit_end},
3646 "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])" },
3647 { "-temp", FALSE, etREAL, {&temp},
3648 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3649 { "-smooth", FALSE, etREAL, {&smooth_tail_start},
3650 "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3651 { "-dump", FALSE, etINT, {&nDump},
3652 "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
3653 { "-max_hb", FALSE, etREAL, {&maxnhb},
3654 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3655 { "-merge", FALSE, etBOOL, {&bMerge},
3656 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3657 { "-geminate", FALSE, etENUM, {gemType},
3658 "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3659 { "-diff", FALSE, etREAL, {&D},
3660 "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."},
3662 { "-nthreads", FALSE, etINT, {&nThreads},
3663 "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 cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
3666 const char *bugs[] = {
3667 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3670 { efTRX, "-f", NULL, ffREAD },
3671 { efTPR, NULL, NULL, ffREAD },
3672 { efNDX, NULL, NULL, ffOPTRD },
3673 /* { efNDX, "-sel", "select", ffOPTRD },*/
3674 { efXVG, "-num", "hbnum", ffWRITE },
3675 { efLOG, "-g", "hbond", ffOPTWR },
3676 { efXVG, "-ac", "hbac", ffOPTWR },
3677 { efXVG, "-dist", "hbdist", ffOPTWR },
3678 { efXVG, "-ang", "hbang", ffOPTWR },
3679 { efXVG, "-hx", "hbhelix", ffOPTWR },
3680 { efNDX, "-hbn", "hbond", ffOPTWR },
3681 { efXPM, "-hbm", "hbmap", ffOPTWR },
3682 { efXVG, "-don", "donor", ffOPTWR },
3683 { efXVG, "-dan", "danum", ffOPTWR },
3684 { efXVG, "-life", "hblife", ffOPTWR },
3685 { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
3688 #define NFILE asize(fnm)
3690 char hbmap [HB_NR] = { ' ', 'o', '-', '*' };
3691 const char *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
3692 t_rgb hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
3694 t_trxstatus *status;
3699 int npargs, natoms, nframes = 0, shatom;
3705 real t, ccut, dist = 0.0, ang = 0.0;
3706 double max_nhb, aver_nhb, aver_dist;
3707 int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
3709 int xj, yj, zj, aj, xjj, yjj, zjj;
3710 int xk, yk, zk, ak, xkk, ykk, zkk;
3711 gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
3712 int *adist, *rdist, *aptr, *rprt;
3713 int grp, nabin, nrbin, bin, resdist, ihb;
3715 t_hbdata *hb, *hbptr;
3716 FILE *fp, *fpins = NULL, *fpnhb = NULL, *donor_properties = NULL;
3718 t_ncell *icell, *jcell, *kcell;
3720 unsigned char *datable;
3725 int ii, jj, hh, actual_nThreads;
3727 gmx_bool bGem, bNN, bParallel;
3728 t_gemParams *params = NULL;
3729 gmx_bool bEdge_yjj, bEdge_xjj, bOMP;
3731 t_hbdata **p_hb = NULL; /* one per thread, then merge after the frame loop */
3732 int **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
3741 ppa = add_acf_pargs(&npargs, pa);
3743 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
3744 ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
3749 /* NN-loop? If so, what estimator to use ?*/
3751 /* Outcommented for now DvdS 2010-07-13
3752 while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3755 gmx_fatal(FARGS, "Invalid NN-loop type.");
3758 for (i = 2; bNN == FALSE && i < NN_NR; i++)
3760 bNN = bNN || NN == i;
3763 if (NN > NN_NONE && bMerge)
3768 /* geminate recombination? If so, which flavor? */
3770 while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0)
3774 if (gemmode == gemNR)
3776 gmx_fatal(FARGS, "Invalid recombination type.");
3780 for (i = 2; bGem == FALSE && i < gemNR; i++)
3782 bGem = bGem || gemmode == i;
3787 printf("Geminate recombination: %s\n", gemType[gemmode]);
3790 if (gemmode != gemDD)
3792 printf("Turning off -contact option...\n");
3798 if (gemmode == gemDD)
3800 printf("Turning on -contact option...\n");
3806 if (gemmode == gemAA)
3808 printf("Turning off -merge option...\n");
3814 if (gemmode != gemAA)
3816 printf("Turning on -merge option...\n");
3824 ccut = cos(acut*DEG2RAD);
3830 gmx_fatal(FARGS, "Can not analyze selected contacts.");
3834 gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
3838 /* Initiate main data structure! */
3839 bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
3840 opt2bSet("-life", NFILE, fnm) ||
3841 opt2bSet("-hbn", NFILE, fnm) ||
3842 opt2bSet("-hbm", NFILE, fnm) ||
3845 if (opt2bSet("-nhbdist", NFILE, fnm))
3847 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3848 fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
3849 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
3850 xvgr_legend(fpnhb, asize(leg), leg, oenv);
3853 hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
3856 read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
3858 snew(grpnames, grNR);
3861 /* Make Donor-Acceptor table */
3862 snew(datable, top.atoms.nr);
3863 gen_datable(index[0], isize[0], datable, top.atoms.nr);
3867 /* analyze selected hydrogen bonds */
3868 printf("Select group with selected atoms:\n");
3869 get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
3870 1, &nsel, index, grpnames);
3873 gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
3874 "and therefore cannot contain triplets of "
3875 "Donor-Hydrogen-Acceptor", grpnames[0]);
3879 for (i = 0; (i < nsel); i += 3)
3881 int dd = index[0][i];
3882 int aa = index[0][i+2];
3883 /* int */ hh = index[0][i+1];
3884 add_dh (&hb->d, dd, hh, i, datable);
3885 add_acc(&hb->a, aa, i);
3886 /* Should this be here ? */
3887 snew(hb->d.dptr, top.atoms.nr);
3888 snew(hb->a.aptr, top.atoms.nr);
3889 add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri);
3891 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
3892 isize[0], grpnames[0]);
3896 /* analyze all hydrogen bonds: get group(s) */
3897 printf("Specify 2 groups to analyze:\n");
3898 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
3899 2, isize, index, grpnames);
3901 /* check if we have two identical or two non-overlapping groups */
3902 bTwo = isize[0] != isize[1];
3903 for (i = 0; (i < isize[0]) && !bTwo; i++)
3905 bTwo = index[0][i] != index[1][i];
3909 printf("Checking for overlap in atoms between %s and %s\n",
3910 grpnames[0], grpnames[1]);
3911 for (i = 0; i < isize[1]; i++)
3913 if (ISINGRP(datable[index[1][i]]))
3915 gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
3916 grpnames[0], grpnames[1]);
3920 printf("Checking for overlap in atoms between %s and %s\n",
3921 grpnames[0],grpnames[1]);
3922 for (i=0; i<isize[0]; i++)
3923 for (j=0; j<isize[1]; j++)
3924 if (index[0][i] == index[1][j])
3925 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3926 grpnames[0],grpnames[1]);
3931 printf("Calculating %s "
3932 "between %s (%d atoms) and %s (%d atoms)\n",
3933 bContact ? "contacts" : "hydrogen bonds",
3934 grpnames[0], isize[0], grpnames[1], isize[1]);
3938 fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
3939 bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
3944 /* search donors and acceptors in groups */
3945 snew(datable, top.atoms.nr);
3946 for (i = 0; (i < grNR); i++)
3948 if ( ((i == gr0) && !bSelected ) ||
3949 ((i == gr1) && bTwo ))
3951 gen_datable(index[i], isize[i], datable, top.atoms.nr);
3954 search_acceptors(&top, isize[i], index[i], &hb->a, i,
3955 bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
3956 search_donors (&top, isize[i], index[i], &hb->d, i,
3957 TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
3961 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
3962 search_donors (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
3966 clear_datable_grp(datable, top.atoms.nr);
3971 printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
3973 snew(donors[gr0D], dons[gr0D].nrd);*/
3975 donor_properties = open_donor_properties_file(opt2fn_null("-don", NFILE, fnm), hb, oenv);
3979 printf("Making hbmap structure...");
3980 /* Generate hbond data structure */
3985 #ifdef HAVE_NN_LOOPS
3994 printf("Making per structure...");
3995 /* Generate hbond data structure */
4002 if (hb->d.nrd + hb->a.nra == 0)
4004 printf("No Donors or Acceptors found\n");
4011 printf("No Donors found\n");
4016 printf("No Acceptors found\n");
4022 gmx_fatal(FARGS, "Nothing to be done");
4031 /* get index group with atom for shell */
4034 printf("Select atom for shell (1 atom):\n");
4035 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4036 1, &shisz, &shidx, &shgrpnm);
4039 printf("group contains %d atoms, should be 1 (one)\n", shisz);
4044 printf("Will calculate hydrogen bonds within a shell "
4045 "of %g nm around atom %i\n", rshell, shatom+1);
4048 /* Analyze trajectory */
4049 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
4050 if (natoms > top.atoms.nr)
4052 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
4053 top.atoms.nr, natoms);
4056 bBox = ir.ePBC != epbcNONE;
4057 grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
4060 snew(adist, nabin+1);
4061 snew(rdist, nrbin+1);
4065 gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
4071 #define __ADIST adist
4072 #define __RDIST rdist
4074 #else /* GMX_OPENMP ================================================== \
4075 * Set up the OpenMP stuff, |
4076 * like the number of threads and such |
4077 * Also start the parallel loop. |
4079 #define __ADIST p_adist[threadNr]
4080 #define __RDIST p_rdist[threadNr]
4081 #define __HBDATA p_hb[threadNr]
4085 bParallel = !bSelected;
4089 actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
4091 gmx_omp_set_num_threads(actual_nThreads);
4092 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
4097 actual_nThreads = 1;
4100 snew(p_hb, actual_nThreads);
4101 snew(p_adist, actual_nThreads);
4102 snew(p_rdist, actual_nThreads);
4103 for (i = 0; i < actual_nThreads; i++)
4106 snew(p_adist[i], nabin+1);
4107 snew(p_rdist[i], nrbin+1);
4109 p_hb[i]->max_frames = 0;
4110 p_hb[i]->nhb = NULL;
4111 p_hb[i]->ndist = NULL;
4112 p_hb[i]->n_bound = NULL;
4113 p_hb[i]->time = NULL;
4114 p_hb[i]->nhx = NULL;
4116 p_hb[i]->bHBmap = hb->bHBmap;
4117 p_hb[i]->bDAnr = hb->bDAnr;
4118 p_hb[i]->bGem = hb->bGem;
4119 p_hb[i]->wordlen = hb->wordlen;
4120 p_hb[i]->nframes = hb->nframes;
4121 p_hb[i]->maxhydro = hb->maxhydro;
4122 p_hb[i]->danr = hb->danr;
4125 p_hb[i]->hbmap = hb->hbmap;
4126 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
4127 p_hb[i]->per = hb->per;
4129 #ifdef HAVE_NN_LOOPS
4130 p_hb[i]->hbE = hb->hbE;
4134 p_hb[i]->nrdist = 0;
4138 /* Make a thread pool here,
4139 * instead of forking anew at every frame. */
4141 #pragma omp parallel \
4143 private(j, h, ii, jj, hh, E, \
4144 xi, yi, zi, xj, yj, zj, threadNr, \
4145 dist, ang, peri, icell, jcell, \
4146 grp, ogrp, ai, aj, xjj, yjj, zjj, \
4147 xk, yk, zk, ihb, id, resdist, \
4148 xkk, ykk, zkk, kcell, ak, k, bTric, \
4149 bEdge_xjj, bEdge_yjj) \
4151 { /* Start of parallel region */
4152 threadNr = gmx_omp_get_thread_num();
4157 bTric = bBox && TRICLINIC(box);
4161 sync_hbdata(p_hb[threadNr], nframes);
4165 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
4166 rshell, ngrid, grid);
4167 reset_nhbonds(&(hb->d));
4169 if (debug && bDebug)
4171 dump_grid(debug, ngrid, grid);
4174 add_frames(hb, nframes);
4175 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
4179 count_da_grid(ngrid, grid, hb->danr[nframes]);
4185 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
4190 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
4191 /* Loop over all atom pairs and estimate interaction energy */
4195 addFramesNN(hb, nframes);
4199 #pragma omp for schedule(dynamic)
4200 for (i = 0; i < hb->d.nrd; i++)
4202 for (j = 0; j < hb->a.nra; j++)
4205 h < (bContact ? 1 : hb->d.nhydro[i]);
4208 if (i == hb->d.nrd || j == hb->a.nra)
4210 gmx_fatal(FARGS, "out of bounds");
4213 /* Get the real atom ids */
4216 hh = hb->d.hydro[i][h];
4218 /* Estimate the energy from the geometry */
4219 E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
4220 /* Store the energy */
4221 storeHbEnergy(hb, i, j, h, E, nframes);
4225 #endif /* HAVE_NN_LOOPS */
4234 /* Do not parallelize this just yet. */
4236 for (ii = 0; (ii < nsel); ii++)
4238 int dd = index[0][i];
4239 int aa = index[0][i+2];
4240 /* int */ hh = index[0][i+1];
4241 ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
4242 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4246 /* add to index if not already there */
4248 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri);
4252 } /* if (bSelected) */
4260 calcBoxProjection(box, hb->per->P);
4263 /* loop over all gridcells (xi,yi,zi) */
4264 /* Removed confusing macro, DvdS 27/12/98 */
4267 /* The outer grid loop will have to do for now. */
4268 #pragma omp for schedule(dynamic)
4269 for (xi = 0; xi < ngrid[XX]; xi++)
4271 for (yi = 0; (yi < ngrid[YY]); yi++)
4273 for (zi = 0; (zi < ngrid[ZZ]); zi++)
4276 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
4277 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
4279 icell = &(grid[zi][yi][xi].d[grp]);
4290 /* loop over all hydrogen atoms from group (grp)
4291 * in this gridcell (icell)
4293 for (ai = 0; (ai < icell->nr); ai++)
4295 i = icell->atoms[ai];
4297 /* loop over all adjacent gridcells (xj,yj,zj) */
4298 for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
4299 zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
4302 zj = grid_mod(zjj, ngrid[ZZ]);
4303 bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
4304 for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
4305 yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
4308 yj = grid_mod(yjj, ngrid[YY]);
4310 (yj == 0) || (yj == ngrid[YY] - 1) ||
4311 (zj == 0) || (zj == ngrid[ZZ] - 1);
4312 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
4313 xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
4316 xj = grid_mod(xjj, ngrid[XX]);
4317 jcell = &(grid[zj][yj][xj].a[ogrp]);
4318 /* loop over acceptor atoms from other group (ogrp)
4319 * in this adjacent gridcell (jcell)
4321 for (aj = 0; (aj < jcell->nr); aj++)
4323 j = jcell->atoms[aj];
4325 /* check if this once was a h-bond */
4327 ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
4328 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4332 /* add to index if not already there */
4334 add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri);
4336 /* make angle and distance distributions */
4337 if (ihb == hbHB && !bContact)
4341 gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
4344 __ADIST[(int)( ang/abin)]++;
4345 __RDIST[(int)(dist/rbin)]++;
4349 if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
4351 gmx_fatal(FARGS, "Invalid donor %d", i);
4353 if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
4355 gmx_fatal(FARGS, "Invalid acceptor %d", j);
4357 resdist = abs(top.atoms.atom[i].resind-
4358 top.atoms.atom[j].resind);
4359 if (resdist >= max_hx)
4363 __HBDATA->nhx[nframes][resdist]++;
4374 } /* for xi,yi,zi */
4377 } /* if (bSelected) {...} else */
4380 /* Better wait for all threads to finnish using x[] before updating it. */
4383 #pragma omp critical
4385 /* Sum up histograms and counts from p_hb[] into hb */
4388 hb->nhb[k] += p_hb[threadNr]->nhb[k];
4389 hb->ndist[k] += p_hb[threadNr]->ndist[k];
4390 for (j = 0; j < max_hx; j++)
4392 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
4397 /* Here are a handful of single constructs
4398 * to share the workload a bit. The most
4399 * important one is of course the last one,
4400 * where there's a potential bottleneck in form
4405 analyse_donor_properties(donor_properties, hb, k, t);
4412 do_nhb_dist(fpnhb, hb, t);
4415 } /* if (bNN) {...} else + */
4419 trrStatus = (read_next_x(oenv, status, &t, x, box));
4429 #pragma omp critical
4431 hb->nrhb += p_hb[threadNr]->nrhb;
4432 hb->nrdist += p_hb[threadNr]->nrdist;
4434 /* Free parallel datastructures */
4435 sfree(p_hb[threadNr]->nhb);
4436 sfree(p_hb[threadNr]->ndist);
4437 sfree(p_hb[threadNr]->nhx);
4440 for (i = 0; i < nabin; i++)
4442 for (j = 0; j < actual_nThreads; j++)
4445 adist[i] += p_adist[j][i];
4449 for (i = 0; i <= nrbin; i++)
4451 for (j = 0; j < actual_nThreads; j++)
4453 rdist[i] += p_rdist[j][i];
4457 sfree(p_adist[threadNr]);
4458 sfree(p_rdist[threadNr]);
4460 } /* End of parallel region */
4467 if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
4469 gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
4472 free_grid(ngrid, &grid);
4476 if (donor_properties)
4478 xvgrclose(donor_properties);
4486 /* Compute maximum possible number of different hbonds */
4493 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
4495 /* Added support for -contact below.
4496 * - Erik Marklund, May 29-31, 2006 */
4497 /* Changed contact code.
4498 * - Erik Marklund, June 29, 2006 */
4503 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
4507 printf("Found %d different %s in trajectory\n"
4508 "Found %d different atom-pairs within %s distance\n",
4509 hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
4510 hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
4512 /*Control the pHist.*/
4516 merge_hb(hb, bTwo, bContact);
4519 if (opt2bSet("-hbn", NFILE, fnm))
4521 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
4524 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
4525 * to make the -hbn and -hmb output match eachother.
4526 * - Erik Marklund, May 30, 2006 */
4529 /* Print out number of hbonds and distances */
4532 fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
4533 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
4535 snew(leg[0], STRLEN);
4536 snew(leg[1], STRLEN);
4537 sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
4538 sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
4539 xvgr_legend(fp, 2, (const char**)leg, oenv);
4543 for (i = 0; (i < nframes); i++)
4545 fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
4546 aver_nhb += hb->nhb[i];
4547 aver_dist += hb->ndist[i];
4550 aver_nhb /= nframes;
4551 aver_dist /= nframes;
4552 /* Print HB distance distribution */
4553 if (opt2bSet("-dist", NFILE, fnm))
4558 for (i = 0; i < nrbin; i++)
4563 fp = xvgropen(opt2fn("-dist", NFILE, fnm),
4564 "Hydrogen Bond Distribution",
4566 "Donor - Acceptor Distance (nm)" :
4567 "Hydrogen - Acceptor Distance (nm)", "", oenv);
4568 for (i = 0; i < nrbin; i++)
4570 fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
4575 /* Print HB angle distribution */
4576 if (opt2bSet("-ang", NFILE, fnm))
4581 for (i = 0; i < nabin; i++)
4586 fp = xvgropen(opt2fn("-ang", NFILE, fnm),
4587 "Hydrogen Bond Distribution",
4588 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
4589 for (i = 0; i < nabin; i++)
4591 fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
4596 /* Print HB in alpha-helix */
4597 if (opt2bSet("-hx", NFILE, fnm))
4599 fp = xvgropen(opt2fn("-hx", NFILE, fnm),
4600 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4601 xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
4602 for (i = 0; i < nframes; i++)
4604 fprintf(fp, "%10g", hb->time[i]);
4605 for (j = 0; j < max_hx; j++)
4607 fprintf(fp, " %6d", hb->nhx[i][j]);
4615 printf("Average number of %s per timeframe %.3f out of %g possible\n",
4616 bContact ? "contacts" : "hbonds",
4617 bContact ? aver_dist : aver_nhb, max_nhb);
4620 /* Do Autocorrelation etc. */
4624 Added support for -contact in ac and hbm calculations below.
4625 - Erik Marklund, May 29, 2006
4629 if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
4631 please_cite(stdout, "Spoel2006b");
4633 if (opt2bSet("-ac", NFILE, fnm))
4635 char *gemstring = NULL;
4639 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4640 gemBallistic, nBalExp);
4643 gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4646 gemstring = gmx_strdup(gemType[hb->per->gemtype]);
4647 do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
4648 bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
4649 gemstring, nThreads, NN, bBallistic, bGemFit);
4651 if (opt2bSet("-life", NFILE, fnm))
4653 do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
4655 if (opt2bSet("-hbm", NFILE, fnm))
4658 int id, ia, hh, x, y;
4659 mat.flags = mat.y0 = 0;
4661 if ((nframes > 0) && (hb->nrhb > 0))
4666 snew(mat.matrix, mat.nx);
4667 for (x = 0; (x < mat.nx); x++)
4669 snew(mat.matrix[x], mat.ny);
4672 for (id = 0; (id < hb->d.nrd); id++)
4674 for (ia = 0; (ia < hb->a.nra); ia++)
4676 for (hh = 0; (hh < hb->maxhydro); hh++)
4678 if (hb->hbmap[id][ia])
4680 if (ISHB(hb->hbmap[id][ia]->history[hh]))
4682 /* Changed '<' into '<=' in the for-statement below.
4683 * It fixed the previously undiscovered bug that caused
4684 * the last occurance of an hbond/contact to not be
4685 * set in mat.matrix. Have a look at any old -hbm-output
4686 * and you will notice that the last column is allways empty.
4687 * - Erik Marklund May 30, 2006
4689 for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
4691 int nn0 = hb->hbmap[id][ia]->n0;
4692 range_check(y, 0, mat.ny);
4693 mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
4701 mat.axis_x = hb->time;
4702 snew(mat.axis_y, mat.ny);
4703 for (j = 0; j < mat.ny; j++)
4707 sprintf(mat.title, bContact ? "Contact Existence Map" :
4708 "Hydrogen Bond Existence Map");
4709 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
4710 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
4711 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4712 mat.bDiscrete = TRUE;
4714 snew(mat.map, mat.nmap);
4715 for (i = 0; i < mat.nmap; i++)
4717 mat.map[i].code.c1 = hbmap[i];
4718 mat.map[i].desc = hbdesc[i];
4719 mat.map[i].rgb = hbrgb[i];
4721 fp = opt2FILE("-hbm", NFILE, fnm, "w");
4722 write_xpm_m(fp, mat);
4724 for (x = 0; x < mat.nx; x++)
4726 sfree(mat.matrix[x]);
4734 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
4741 fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4742 fprintf(stderr, "Freeing pHist for all donors...\n");
4743 for (i = 0; i < hb->d.nrd; i++)
4745 fprintf(stderr, "\r%i", i);
4746 if (hb->per->pHist[i] != NULL)
4748 for (j = 0; j < hb->a.nra; j++)
4750 clearPshift(&(hb->per->pHist[i][j]));
4752 sfree(hb->per->pHist[i]);
4755 sfree(hb->per->pHist);
4756 sfree(hb->per->p2i);
4758 fprintf(stderr, "...done.\n");
4761 #ifdef HAVE_NN_LOOPS
4774 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4776 fp = xvgropen(opt2fn("-dan", NFILE, fnm),
4777 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4778 nleg = (bTwo ? 2 : 1)*2;
4779 snew(legnames, nleg);
4781 for (j = 0; j < grNR; j++)
4783 if (USE_THIS_GROUP(j) )
4785 sprintf(buf, "Donors %s", grpnames[j]);
4786 legnames[i++] = gmx_strdup(buf);
4787 sprintf(buf, "Acceptors %s", grpnames[j]);
4788 legnames[i++] = gmx_strdup(buf);
4793 gmx_incons("number of legend entries");
4795 xvgr_legend(fp, nleg, (const char**)legnames, oenv);
4796 for (i = 0; i < nframes; i++)
4798 fprintf(fp, "%10g", hb->time[i]);
4799 for (j = 0; (j < grNR); j++)
4801 if (USE_THIS_GROUP(j) )
4803 fprintf(fp, " %6d", hb->danr[i][j]);