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, 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.
42 /*#define HAVE_NN_LOOPS*/
44 #include "gromacs/commandline/pargs.h"
50 #include "gmx_fatal.h"
52 #include "gromacs/utility/smalloc.h"
56 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/fileio/futil.h"
63 #include "gromacs/fileio/matio.h"
64 #include "gromacs/fileio/tpxio.h"
65 #include "gromacs/fileio/trxio.h"
66 #include "gromacs/utility/gmxomp.h"
68 typedef short int t_E;
71 typedef int t_hx[max_hx];
72 #define NRHXTYPES max_hx
73 const char *hxtypenames[NRHXTYPES] =
74 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
78 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == 0)
80 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == (threadNr))
83 /* -----------------------------------------*/
89 hbNo, hbDist, hbHB, hbNR, hbR2
92 noDA, ACC, DON, DA, INGROUP
95 NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR
98 static const char *grpnames[grNR] = {"0", "1", "I" };
100 static gmx_bool bDebug = FALSE;
105 #define HB_YESINS HB_YES|HB_INS
109 #define ISHB(h) (((h) & 2) == 2)
110 #define ISDIST(h) (((h) & 1) == 1)
111 #define ISDIST2(h) (((h) & 4) == 4)
112 #define ISACC(h) (((h) & 1) == 1)
113 #define ISDON(h) (((h) & 2) == 2)
114 #define ISINGRP(h) (((h) & 4) == 4)
127 typedef int t_icell[grNR];
128 typedef atom_id h_id[MAXHYDRO];
131 int history[MAXHYDRO];
132 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
133 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
135 /* Bitmask array which tells whether a hbond is present
136 * at a given time. Either of these may be NULL
138 int n0; /* First frame a HB was found */
139 int nframes, maxframes; /* Amount of frames in this hbond */
142 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
143 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
144 * acceptor distance is less than the user-specified distance (typically
151 atom_id *acc; /* Atom numbers of the acceptors */
152 int *grp; /* Group index */
153 int *aptr; /* Map atom number to acceptor index */
158 int *don; /* Atom numbers of the donors */
159 int *grp; /* Group index */
160 int *dptr; /* Map atom number to donor index */
161 int *nhydro; /* Number of hydrogens for each donor */
162 h_id *hydro; /* The atom numbers of the hydrogens */
163 h_id *nhbonds; /* The number of HBs per H at current */
166 /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/
170 int len; /* The length of frame and p. */
171 int *frame; /* The frames at which transitio*/
176 /* Periodicity history. Used for the reversible geminate recombination. */
177 t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap:
178 * pHist[d][a]. We can safely assume that the same
179 * periodic shift holds for all hydrogens of a da-pair.
181 * Nowadays it only stores TRANSITIONS, and not the shift at every frame.
182 * That saves a LOT of memory, an hopefully kills a mysterious bug where
183 * pHist gets contaminated. */
185 PSTYPE nper; /* The length of p2i */
186 ivec *p2i; /* Maps integer to periodic shift for a pair.*/
187 matrix P; /* Projection matrix to find the box shifts. */
188 int gemtype; /* enumerated type */
193 int *Etot; /* Total energy for each frame */
194 t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */
198 gmx_bool bHBmap, bDAnr, bGem;
200 /* The following arrays are nframes long */
201 int nframes, max_frames, maxhydro;
207 /* These structures are initialized from the topology at start up */
210 /* This holds a matrix with all possible hydrogen bonds */
216 /* For parallelization reasons this will have to be a pointer.
217 * Otherwise discrepancies may arise between the periodicity data
218 * seen by different threads. */
222 static void clearPshift(t_pShift *pShift)
227 sfree(pShift->frame);
232 static void calcBoxProjection(matrix B, matrix P)
234 const int vp[] = {XX, YY, ZZ};
239 for (i = 0; i < 3; i++)
242 for (j = 0; j < 3; j++)
245 U[m][n] = i == j ? 1 : 0;
249 for (i = 0; i < 3; i++)
252 mvmul(M, U[m], P[m]);
257 static void calcBoxDistance(matrix P, rvec d, ivec ibd)
259 /* returns integer distance in box coordinates.
260 * P is the projection matrix from cartesian coordinates
261 * obtained with calcBoxProjection(). */
265 /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
266 for (i = 0; i < 3; i++)
268 bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5);
270 ibd[XX] = (int)bd[XX];
271 ibd[YY] = (int)bd[YY];
272 ibd[ZZ] = (int)bd[ZZ];
275 /* Changed argument 'bMerge' into 'oneHB' below,
276 * since -contact should cause maxhydro to be 1,
278 * - Erik Marklund May 29, 2006
281 static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap)
283 /* Try to merge hbonds on the fly. That means that if the
284 * acceptor and donor are mergable, then:
285 * 1) store the hb-info so that acceptor id > donor id,
286 * 2) add the periodic shift in pairs, so that [-x,-y,-z] is
287 * stored in per.p2i[] whenever acceptor id < donor id.
288 * Note that [0,0,0] should already be the first element of per.p2i
289 * by the time this function is called. */
291 /* daSwap is TRUE if the donor and acceptor were swapped.
292 * If so, then the negative vector should be used. */
295 if (per->p2i == NULL || per->nper == 0)
297 gmx_fatal(FARGS, "'per' not initialized properly.");
299 for (i = 0; i < per->nper; i++)
301 if (r[XX] == per->p2i[i][XX] &&
302 r[YY] == per->p2i[i][YY] &&
303 r[ZZ] == per->p2i[i][ZZ])
308 /* Not found apparently. Add it to the list! */
309 /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
315 fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
320 srenew(per->p2i, per->nper+2);
322 copy_ivec(r, per->p2i[per->nper]);
325 /* Add the mirror too. It's rather likely that it'll be needed. */
326 per->p2i[per->nper][XX] = -r[XX];
327 per->p2i[per->nper][YY] = -r[YY];
328 per->p2i[per->nper][ZZ] = -r[ZZ];
331 return per->nper - 1 - (daSwap ? 0 : 1);
334 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode)
339 hb->wordlen = 8*sizeof(unsigned int);
349 hb->maxhydro = MAXHYDRO;
352 hb->per->gemtype = bGem ? gemmode : 0;
357 static void mk_hbmap(t_hbdata *hb)
361 snew(hb->hbmap, hb->d.nrd);
362 for (i = 0; (i < hb->d.nrd); i++)
364 snew(hb->hbmap[i], hb->a.nra);
365 if (hb->hbmap[i] == NULL)
367 gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
369 for (j = 0; (j > hb->a.nra); j++)
371 hb->hbmap[i][j] = NULL;
376 /* Consider redoing pHist so that is only stores transitions between
377 * periodicities and not the periodicity for all frames. This eats heaps of memory. */
378 static void mk_per(t_hbdata *hb)
383 snew(hb->per->pHist, hb->d.nrd);
384 for (i = 0; i < hb->d.nrd; i++)
386 snew(hb->per->pHist[i], hb->a.nra);
387 if (hb->per->pHist[i] == NULL)
389 gmx_fatal(FARGS, "Could not allocate enough memory for per->pHist");
391 for (j = 0; j < hb->a.nra; j++)
393 clearPshift(&(hb->per->pHist[i][j]));
396 /* add the [0,0,0] shift to element 0 of p2i. */
397 snew(hb->per->p2i, 1);
398 clear_ivec(hb->per->p2i[0]);
404 static void mk_hbEmap (t_hbdata *hb, int n0)
409 snew(hb->hbE.E, hb->d.nrd);
410 for (i = 0; i < hb->d.nrd; i++)
412 snew(hb->hbE.E[i], hb->a.nra);
413 for (j = 0; j < hb->a.nra; j++)
415 snew(hb->hbE.E[i][j], MAXHYDRO);
416 for (k = 0; k < MAXHYDRO; k++)
418 hb->hbE.E[i][j][k] = NULL;
425 static void free_hbEmap (t_hbdata *hb)
428 for (i = 0; i < hb->d.nrd; i++)
430 for (j = 0; j < hb->a.nra; j++)
432 for (k = 0; k < MAXHYDRO; k++)
434 sfree(hb->hbE.E[i][j][k]);
436 sfree(hb->hbE.E[i][j]);
444 static void addFramesNN(t_hbdata *hb, int frame)
447 #define DELTAFRAMES_HBE 10
449 int d, a, h, nframes;
451 if (frame >= hb->hbE.nframes)
453 nframes = hb->hbE.nframes + DELTAFRAMES_HBE;
454 srenew(hb->hbE.Etot, nframes);
456 for (d = 0; d < hb->d.nrd; d++)
458 for (a = 0; a < hb->a.nra; a++)
460 for (h = 0; h < hb->d.nhydro[d]; h++)
462 srenew(hb->hbE.E[d][a][h], nframes);
467 hb->hbE.nframes += DELTAFRAMES_HBE;
471 static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst,
472 matrix box, rvec hbox, t_donors *donors)
477 * alpha - angle between dipoles
478 * x[] - atomic positions
479 * EEst - the type of energy estimate (see enum in hbplugin.h)
480 * box - the box vectors \
481 * hbox - half box lengths _These two are only needed for the pbc correction
486 rvec dipole[2], xmol[3], xmean[2];
492 /* Self-interaction */
499 /* This is a simple binary existence function that sets E=1 whenever
500 * the distance between the oxygens is equal too or less than 0.35 nm.
502 rvec_sub(x[d], x[a], dist);
503 pbc_correct_gem(dist, box, hbox);
504 if (norm(dist) <= 0.35)
515 /* Negative potential energy of a dipole.
516 * E = -cos(alpha) * 1/r^3 */
518 copy_rvec(x[d], xmol[0]); /* donor */
519 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
520 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
522 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
523 rvec_inc(xmean[0], xmol[1]);
524 rvec_inc(xmean[0], xmol[2]);
525 for (i = 0; i < 3; i++)
527 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
530 /* Assumes that all acceptors are also donors. */
531 copy_rvec(x[a], xmol[0]); /* acceptor */
532 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
533 copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */
536 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
537 rvec_inc(xmean[1], xmol[1]);
538 rvec_inc(xmean[1], xmol[2]);
539 for (i = 0; i < 3; i++)
541 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
544 rvec_sub(xmean[0], xmean[1], dist);
545 pbc_correct_gem(dist, box, hbox);
548 realE = pow(r, -3.0);
549 E = (t_E)(SCALEFACTOR_E * realE);
553 /* Negative potential energy of a (unpolarizable) dipole.
554 * E = -cos(alpha) * 1/r^3 */
555 clear_rvec(dipole[1]);
556 clear_rvec(dipole[0]);
558 copy_rvec(x[d], xmol[0]); /* donor */
559 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
560 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
562 rvec_inc(dipole[0], xmol[1]);
563 rvec_inc(dipole[0], xmol[2]);
564 for (i = 0; i < 3; i++)
568 rvec_dec(dipole[0], xmol[0]);
570 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
571 rvec_inc(xmean[0], xmol[1]);
572 rvec_inc(xmean[0], xmol[2]);
573 for (i = 0; i < 3; i++)
575 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
578 /* Assumes that all acceptors are also donors. */
579 copy_rvec(x[a], xmol[0]); /* acceptor */
580 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
581 copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */
584 rvec_inc(dipole[1], xmol[1]);
585 rvec_inc(dipole[1], xmol[2]);
586 for (i = 0; i < 3; i++)
590 rvec_dec(dipole[1], xmol[0]);
592 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
593 rvec_inc(xmean[1], xmol[1]);
594 rvec_inc(xmean[1], xmol[2]);
595 for (i = 0; i < 3; i++)
597 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
600 rvec_sub(xmean[0], xmean[1], dist);
601 pbc_correct_gem(dist, box, hbox);
604 double cosalpha = cos_angle(dipole[0], dipole[1]);
605 realE = cosalpha * pow(r, -3.0);
606 E = (t_E)(SCALEFACTOR_E * realE);
610 printf("Can't do that type of energy estimate: %i\n.", EEst);
617 static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame)
619 /* hb - hbond data structure
623 E - estimate of the energy
624 frame - the current frame.
627 /* Store the estimated energy */
633 hb->hbE.E[d][a][h][frame] = E;
637 hb->hbE.Etot[frame] += E;
640 #endif /* HAVE_NN_LOOPS */
643 /* Finds -v[] in the periodicity index */
644 static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
648 for (i = 0; i < nper; i++)
650 if (v[i][XX] == -(v[p][XX]) &&
651 v[i][YY] == -(v[p][YY]) &&
652 v[i][ZZ] == -(v[p][ZZ]))
657 printf("Couldn't find mirror of [%i, %i, %i], index \n",
665 static void add_frames(t_hbdata *hb, int nframes)
669 if (nframes >= hb->max_frames)
671 hb->max_frames += 4096;
672 srenew(hb->time, hb->max_frames);
673 srenew(hb->nhb, hb->max_frames);
674 srenew(hb->ndist, hb->max_frames);
675 srenew(hb->n_bound, hb->max_frames);
676 srenew(hb->nhx, hb->max_frames);
679 srenew(hb->danr, hb->max_frames);
682 hb->nframes = nframes;
685 #define OFFSET(frame) (frame / 32)
686 #define MASK(frame) (1 << (frame % 32))
688 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
692 hbexist[OFFSET(frame)] |= MASK(frame);
696 hbexist[OFFSET(frame)] &= ~MASK(frame);
700 static gmx_bool is_hb(unsigned int hbexist[], int frame)
702 return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
705 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
707 unsigned int *ghptr = NULL;
711 ghptr = hb->hbmap[id][ia]->h[ih];
713 else if (ihb == hbDist)
715 ghptr = hb->hbmap[id][ia]->g[ih];
719 gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
722 _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
725 static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
729 snew(pHist->frame, 1);
732 pHist->frame[0] = frame;
737 if (pHist->p[pHist->len-1] != p)
740 srenew(pHist->frame, pHist->len);
741 srenew(pHist->p, pHist->len);
742 pHist->frame[pHist->len-1] = frame;
743 pHist->p[pHist->len-1] = p;
744 } /* Otherwise, there is no transition. */
748 static PSTYPE getPshift(t_pShift pHist, int frame)
753 || (pHist.len > 0 && pHist.frame[0] > frame))
758 for (i = 0; i < pHist.len; i++)
771 /* It seems that frame is after the last periodic transition. Return the last periodicity. */
772 return pHist.p[pHist.len-1];
775 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPE p)
778 t_hbond *hb = hbd->hbmap[id][ia];
779 int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]);
780 int wlen = hbd->wordlen;
782 gmx_bool bGem = hbd->bGem;
787 hb->maxframes = delta;
788 for (i = 0; (i < maxhydro); i++)
790 snew(hb->h[i], hb->maxframes/wlen);
791 snew(hb->g[i], hb->maxframes/wlen);
796 hb->nframes = frame-hb->n0;
797 /* We need a while loop here because hbonds may be returning
800 while (hb->nframes >= hb->maxframes)
802 n = hb->maxframes + delta;
803 for (i = 0; (i < maxhydro); i++)
805 srenew(hb->h[i], n/wlen);
806 srenew(hb->g[i], n/wlen);
807 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
819 set_hb(hbd, id, h, ia, frame, ihb);
822 if (p >= hbd->per->nper)
824 gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
828 addPshift(&(hbd->per->pHist[id][ia]), p, frame);
836 static void inc_nhbonds(t_donors *ddd, int d, int h)
839 int dptr = ddd->dptr[d];
841 for (j = 0; (j < ddd->nhydro[dptr]); j++)
843 if (ddd->hydro[dptr][j] == h)
845 ddd->nhbonds[dptr][j]++;
849 if (j == ddd->nhydro[dptr])
851 gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
855 static int _acceptor_index(t_acceptors *a, int grp, atom_id i,
856 const char *file, int line)
860 if (a->grp[ai] != grp)
864 fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
865 ai, a->grp[ai], grp, file, line);
874 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
876 static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line)
885 if (d->grp[di] != grp)
889 fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
890 di, d->grp[di], grp, file, line);
899 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
901 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
903 /* g_hbond doesn't allow overlapping groups */
909 donor_index(&hb->d, grpd, a) != NOTSET
910 && acceptor_index(&hb->a, grpa, d) != NOTSET;
914 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
915 int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPE p)
918 gmx_bool daSwap = FALSE;
920 if ((id = hb->d.dptr[d]) == NOTSET)
922 gmx_fatal(FARGS, "No donor atom %d", d+1);
924 else if (grpd != hb->d.grp[id])
926 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
927 grpd, hb->d.grp[id], d+1);
929 if ((ia = hb->a.aptr[a]) == NOTSET)
931 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
933 else if (grpa != hb->a.grp[ia])
935 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
936 grpa, hb->a.grp[ia], a+1);
942 if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
943 /* Then swap identity so that the id of d is lower then that of a.
945 * This should really be redundant by now, as is_hbond() now ought to return
946 * hbNo in the cases where this conditional is TRUE. */
953 /* Now repeat donor/acc check. */
954 if ((id = hb->d.dptr[d]) == NOTSET)
956 gmx_fatal(FARGS, "No donor atom %d", d+1);
958 else if (grpd != hb->d.grp[id])
960 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
961 grpd, hb->d.grp[id], d+1);
963 if ((ia = hb->a.aptr[a]) == NOTSET)
965 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
967 else if (grpa != hb->a.grp[ia])
969 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
970 grpa, hb->a.grp[ia], a+1);
977 /* Loop over hydrogens to find which hydrogen is in this particular HB */
978 if ((ihb == hbHB) && !bMerge && !bContact)
980 for (k = 0; (k < hb->d.nhydro[id]); k++)
982 if (hb->d.hydro[id][k] == h)
987 if (k == hb->d.nhydro[id])
989 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
1001 #pragma omp critical
1003 if (hb->hbmap[id][ia] == NULL)
1005 snew(hb->hbmap[id][ia], 1);
1006 snew(hb->hbmap[id][ia]->h, hb->maxhydro);
1007 snew(hb->hbmap[id][ia]->g, hb->maxhydro);
1009 add_ff(hb, id, k, ia, frame, ihb, p);
1013 /* Strange construction with frame >=0 is a relic from old code
1014 * for selected hbond analysis. It may be necessary again if that
1015 * is made to work again.
1019 hh = hb->hbmap[id][ia]->history[k];
1025 hb->hbmap[id][ia]->history[k] = hh | 2;
1036 hb->hbmap[id][ia]->history[k] = hh | 1;
1060 if (bMerge && daSwap)
1062 h = hb->d.hydro[id][0];
1064 /* Increment number if HBonds per H */
1065 if (ihb == hbHB && !bContact)
1067 inc_nhbonds(&(hb->d), d, h);
1071 static char *mkatomname(t_atoms *atoms, int i)
1073 static char buf[32];
1076 rnr = atoms->atom[i].resind;
1077 sprintf(buf, "%4s%d%-4s",
1078 *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
1083 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms)
1085 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
1088 for (i = 0; i < isize; i++)
1090 if (index[i] >= natoms)
1092 gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
1094 datable[index[i]] |= INGROUP;
1098 static void clear_datable_grp(unsigned char *datable, int size)
1100 /* Clears group information from the table */
1102 const char mask = !(char)INGROUP;
1105 for (i = 0; i < size; i++)
1112 static void add_acc(t_acceptors *a, int ia, int grp)
1114 if (a->nra >= a->max_nra)
1117 srenew(a->acc, a->max_nra);
1118 srenew(a->grp, a->max_nra);
1120 a->grp[a->nra] = grp;
1121 a->acc[a->nra++] = ia;
1124 static void search_acceptors(t_topology *top, int isize,
1125 atom_id *index, t_acceptors *a, int grp,
1127 gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
1133 for (i = 0; (i < isize); i++)
1137 (((*top->atoms.atomname[n])[0] == 'O') ||
1138 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
1139 ISINGRP(datable[n]))
1141 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
1146 snew(a->aptr, top->atoms.nr);
1147 for (i = 0; (i < top->atoms.nr); i++)
1149 a->aptr[i] = NOTSET;
1151 for (i = 0; (i < a->nra); i++)
1153 a->aptr[a->acc[i]] = i;
1157 static void add_h2d(int id, int ih, t_donors *ddd)
1161 for (i = 0; (i < ddd->nhydro[id]); i++)
1163 if (ddd->hydro[id][i] == ih)
1165 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
1170 if (i == ddd->nhydro[id])
1172 if (ddd->nhydro[id] >= MAXHYDRO)
1174 gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
1175 ddd->don[id], MAXHYDRO);
1177 ddd->hydro[id][i] = ih;
1182 static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable)
1186 if (ISDON(datable[id]) || !datable)
1188 if (ddd->dptr[id] == NOTSET) /* New donor */
1200 if (ddd->nrd >= ddd->max_nrd)
1202 ddd->max_nrd += 128;
1203 srenew(ddd->don, ddd->max_nrd);
1204 srenew(ddd->nhydro, ddd->max_nrd);
1205 srenew(ddd->hydro, ddd->max_nrd);
1206 srenew(ddd->nhbonds, ddd->max_nrd);
1207 srenew(ddd->grp, ddd->max_nrd);
1209 ddd->don[ddd->nrd] = id;
1210 ddd->nhydro[ddd->nrd] = 0;
1211 ddd->grp[ddd->nrd] = grp;
1218 add_h2d(i, ih, ddd);
1223 printf("Warning: Atom %d is not in the d/a-table!\n", id);
1227 static void search_donors(t_topology *top, int isize, atom_id *index,
1228 t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
1229 unsigned char *datable)
1232 t_functype func_type;
1233 t_ilist *interaction;
1234 atom_id nr1, nr2, nr3;
1239 snew(ddd->dptr, top->atoms.nr);
1240 for (i = 0; (i < top->atoms.nr); i++)
1242 ddd->dptr[i] = NOTSET;
1250 for (i = 0; (i < isize); i++)
1252 datable[index[i]] |= DON;
1253 add_dh(ddd, index[i], -1, grp, datable);
1259 for (func_type = 0; (func_type < F_NRE); func_type++)
1261 interaction = &(top->idef.il[func_type]);
1262 if (func_type == F_POSRES || func_type == F_FBPOSRES)
1264 /* The ilist looks strange for posre. Bug in grompp?
1265 * We don't need posre interactions for hbonds anyway.*/
1268 for (i = 0; i < interaction->nr;
1269 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1272 if (func_type != top->idef.functype[interaction->iatoms[i]])
1274 fprintf(stderr, "Error in func_type %s",
1275 interaction_function[func_type].longname);
1279 /* check out this functype */
1280 if (func_type == F_SETTLE)
1282 nr1 = interaction->iatoms[i+1];
1283 nr2 = interaction->iatoms[i+2];
1284 nr3 = interaction->iatoms[i+3];
1286 if (ISINGRP(datable[nr1]))
1288 if (ISINGRP(datable[nr2]))
1290 datable[nr1] |= DON;
1291 add_dh(ddd, nr1, nr1+1, grp, datable);
1293 if (ISINGRP(datable[nr3]))
1295 datable[nr1] |= DON;
1296 add_dh(ddd, nr1, nr1+2, grp, datable);
1300 else if (IS_CHEMBOND(func_type))
1302 for (j = 0; j < 2; j++)
1304 nr1 = interaction->iatoms[i+1+j];
1305 nr2 = interaction->iatoms[i+2-j];
1306 if ((*top->atoms.atomname[nr1][0] == 'H') &&
1307 ((*top->atoms.atomname[nr2][0] == 'O') ||
1308 (*top->atoms.atomname[nr2][0] == 'N')) &&
1309 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1311 datable[nr2] |= DON;
1312 add_dh(ddd, nr2, nr1, grp, datable);
1319 for (func_type = 0; func_type < F_NRE; func_type++)
1321 interaction = &top->idef.il[func_type];
1322 for (i = 0; i < interaction->nr;
1323 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1326 if (func_type != top->idef.functype[interaction->iatoms[i]])
1328 gmx_incons("function type in search_donors");
1331 if (interaction_function[func_type].flags & IF_VSITE)
1333 nr1 = interaction->iatoms[i+1];
1334 if (*top->atoms.atomname[nr1][0] == 'H')
1338 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1349 if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1350 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1351 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1353 datable[nr2] |= DON;
1354 add_dh(ddd, nr2, nr1, grp, datable);
1364 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
1371 for (i = 0; i < DIM; i++)
1373 ngrid[i] = (box[i][i]/(1.2*rcut));
1377 if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
1379 for (i = 0; i < DIM; i++)
1386 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1387 ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
1389 snew(grid, ngrid[ZZ]);
1390 for (z = 0; z < ngrid[ZZ]; z++)
1392 snew((grid)[z], ngrid[YY]);
1393 for (y = 0; y < ngrid[YY]; y++)
1395 snew((grid)[z][y], ngrid[XX]);
1401 static void reset_nhbonds(t_donors *ddd)
1405 for (i = 0; (i < ddd->nrd); i++)
1407 for (j = 0; (j < MAXHH); j++)
1409 ddd->nhbonds[i][j] = 0;
1414 void pbc_correct_gem(rvec dx, matrix box, rvec hbox);
1416 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
1417 gmx_bool bBox, matrix box, rvec hbox,
1418 real rcut, real rshell,
1419 ivec ngrid, t_gridcell ***grid)
1421 int i, m, gr, xi, yi, zi, nr;
1424 rvec invdelta, dshell, xtemp = {0, 0, 0};
1426 gmx_bool bDoRshell, bInShell, bAcc;
1431 bDoRshell = (rshell > 0);
1432 rshell2 = sqr(rshell);
1435 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
1437 for (m = 0; m < DIM; m++)
1439 hbox[m] = box[m][m]*0.5;
1442 invdelta[m] = ngrid[m]/box[m][m];
1443 if (1/invdelta[m] < rcut)
1445 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
1446 "%s can not handle this situation, sorry.\n",
1459 /* resetting atom counts */
1460 for (gr = 0; (gr < grNR); gr++)
1462 for (zi = 0; zi < ngrid[ZZ]; zi++)
1464 for (yi = 0; yi < ngrid[YY]; yi++)
1466 for (xi = 0; xi < ngrid[XX]; xi++)
1468 grid[zi][yi][xi].d[gr].nr = 0;
1469 grid[zi][yi][xi].a[gr].nr = 0;
1475 /* put atoms in grid cells */
1476 for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
1489 for (i = 0; (i < nr); i++)
1491 /* check if we are inside the shell */
1492 /* if bDoRshell=FALSE then bInShell=TRUE always */
1497 rvec_sub(x[ad[i]], xshell, dshell);
1500 if (FALSE && !hb->bGem)
1502 for (m = DIM-1; m >= 0 && bInShell; m--)
1504 if (dshell[m] < -hbox[m])
1506 rvec_inc(dshell, box[m]);
1508 else if (dshell[m] >= hbox[m])
1510 dshell[m] -= 2*hbox[m];
1512 /* if we're outside the cube, we're outside the sphere also! */
1513 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1521 gmx_bool bDone = FALSE;
1525 for (m = DIM-1; m >= 0 && bInShell; m--)
1527 if (dshell[m] < -hbox[m])
1530 rvec_inc(dshell, box[m]);
1532 if (dshell[m] >= hbox[m])
1535 dshell[m] -= 2*hbox[m];
1539 for (m = DIM-1; m >= 0 && bInShell; m--)
1541 /* if we're outside the cube, we're outside the sphere also! */
1542 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1549 /* if we're inside the cube, check if we're inside the sphere */
1552 bInShell = norm2(dshell) < rshell2;
1562 copy_rvec(x[ad[i]], xtemp);
1564 pbc_correct_gem(x[ad[i]], box, hbox);
1566 for (m = DIM-1; m >= 0; m--)
1568 if (TRUE || !hb->bGem)
1570 /* put atom in the box */
1571 while (x[ad[i]][m] < 0)
1573 rvec_inc(x[ad[i]], box[m]);
1575 while (x[ad[i]][m] >= box[m][m])
1577 rvec_dec(x[ad[i]], box[m]);
1580 /* determine grid index of atom */
1581 grididx[m] = x[ad[i]][m]*invdelta[m];
1582 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1586 copy_rvec(xtemp, x[ad[i]]); /* copy back */
1591 range_check(gx, 0, ngrid[XX]);
1592 range_check(gy, 0, ngrid[YY]);
1593 range_check(gz, 0, ngrid[ZZ]);
1597 /* add atom to grid cell */
1600 newgrid = &(grid[gz][gy][gx].a[gr]);
1604 newgrid = &(grid[gz][gy][gx].d[gr]);
1606 if (newgrid->nr >= newgrid->maxnr)
1608 newgrid->maxnr += 10;
1609 DBB(newgrid->maxnr);
1610 srenew(newgrid->atoms, newgrid->maxnr);
1613 newgrid->atoms[newgrid->nr] = ad[i];
1621 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1625 for (gr = 0; (gr < grNR); gr++)
1628 for (zi = 0; zi < ngrid[ZZ]; zi++)
1630 for (yi = 0; yi < ngrid[YY]; yi++)
1632 for (xi = 0; xi < ngrid[XX]; xi++)
1634 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1642 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1643 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1644 * With a triclinic box all loops are 3 long, except when a cell is
1645 * located next to one of the box edges which is not parallel to the
1646 * x/y-plane, in that case all cells in a line or layer are searched.
1647 * This could be implemented slightly more efficient, but the code
1648 * would get much more complicated.
1650 static gmx_inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1652 return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
1654 static gmx_inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1656 return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1658 static gmx_inline int grid_mod(int j, int n)
1663 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1665 int gr, x, y, z, sum[grNR];
1667 fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1668 for (gr = 0; gr < grNR; gr++)
1671 fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1672 for (z = 0; z < ngrid[ZZ]; z += 2)
1674 fprintf(fp, "Z=%d,%d\n", z, z+1);
1675 for (y = 0; y < ngrid[YY]; y++)
1677 for (x = 0; x < ngrid[XX]; x++)
1679 fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1680 sum[gr] += grid[z][y][x].d[gr].nr;
1681 fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1682 sum[gr] += grid[z][y][x].a[gr].nr;
1686 if ( (z+1) < ngrid[ZZ])
1688 for (x = 0; x < ngrid[XX]; x++)
1690 fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1691 sum[gr] += grid[z+1][y][x].d[gr].nr;
1692 fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1693 sum[gr] += grid[z+1][y][x].a[gr].nr;
1700 fprintf(fp, "TOTALS:");
1701 for (gr = 0; gr < grNR; gr++)
1703 fprintf(fp, " %d=%d", gr, sum[gr]);
1708 /* New GMX record! 5 * in a row. Congratulations!
1709 * Sorry, only four left.
1711 static void free_grid(ivec ngrid, t_gridcell ****grid)
1714 t_gridcell ***g = *grid;
1716 for (z = 0; z < ngrid[ZZ]; z++)
1718 for (y = 0; y < ngrid[YY]; y++)
1728 void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1731 gmx_bool bDone = FALSE;
1735 for (m = DIM-1; m >= 0; m--)
1737 if (dx[m] < -hbox[m])
1740 rvec_inc(dx, box[m]);
1742 if (dx[m] >= hbox[m])
1745 rvec_dec(dx, box[m]);
1751 /* Added argument r2cut, changed contact and implemented
1752 * use of second cut-off.
1753 * - Erik Marklund, June 29, 2006
1755 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1756 real rcut, real r2cut, real ccut,
1757 rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1758 real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1759 gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1761 int h, hh, id, ja, ihb;
1762 rvec r_da, r_ha, r_dh, r = {0, 0, 0};
1764 real rc2, r2c2, rda2, rha2, ca;
1765 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1766 gmx_bool daSwap = FALSE;
1773 if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1774 ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
1782 rvec_sub(x[d], x[a], r_da);
1783 /* Insert projection code here */
1785 if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1787 /* Then this hbond/contact will be found again, or it has already been found. */
1792 if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1793 { /* return hbNo; */
1794 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1798 copy_rvec(r_da, r); /* Save this for later */
1799 pbc_correct_gem(r_da, box, hbox);
1803 pbc_correct_gem(r_da, box, hbox);
1806 rda2 = iprod(r_da, r_da);
1810 if (daSwap && grpa == grpd)
1818 calcBoxDistance(hb->per->P, r, ri);
1819 *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */
1823 else if (rda2 < r2c2)
1834 if (bDA && (rda2 > rc2))
1839 for (h = 0; (h < hb->d.nhydro[id]); h++)
1841 hh = hb->d.hydro[id][h];
1845 rvec_sub(x[hh], x[a], r_ha);
1848 pbc_correct_gem(r_ha, box, hbox);
1850 rha2 = iprod(r_ha, r_ha);
1855 calcBoxDistance(hb->per->P, r, ri);
1856 *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */
1859 if (bDA || (!bDA && (rha2 <= rc2)))
1861 rvec_sub(x[d], x[hh], r_dh);
1864 pbc_correct_gem(r_dh, box, hbox);
1871 ca = cos_angle(r_dh, r_da);
1872 /* if angle is smaller, cos is larger */
1876 *d_ha = sqrt(bDA ? rda2 : rha2);
1882 if (bDA || (!bDA && HAinrange))
1892 /* Fixed previously undiscovered bug in the merge
1893 code, where the last frame of each hbond disappears.
1894 - Erik Marklund, June 1, 2006 */
1895 /* Added the following arguments:
1896 * ptmp[] - temporary periodicity hisory
1897 * a1 - identity of first acceptor/donor
1898 * a2 - identity of second acceptor/donor
1899 * - Erik Marklund, FEB 20 2010 */
1901 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1902 * Will do some more testing before removing the function entirely.
1903 * - Erik Marklund, MAY 10 2010 */
1904 static void do_merge(t_hbdata *hb, int ntmp,
1905 unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[],
1906 t_hbond *hb0, t_hbond *hb1, int a1, int a2)
1908 /* Here we need to make sure we're treating periodicity in
1909 * the right way for the geminate recombination kinetics. */
1911 int m, mm, n00, n01, nn0, nnframes;
1915 /* Decide where to start from when merging */
1918 nn0 = min(n00, n01);
1919 nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1920 /* Initiate tmp arrays */
1921 for (m = 0; (m < ntmp); m++)
1927 /* Fill tmp arrays with values due to first HB */
1928 /* Once again '<' had to be replaced with '<='
1929 to catch the last frame in which the hbond
1931 - Erik Marklund, June 1, 2006 */
1932 for (m = 0; (m <= hb0->nframes); m++)
1935 htmp[mm] = is_hb(hb0->h[0], m);
1938 pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1939 if (pm > hb->per->nper)
1941 gmx_fatal(FARGS, "Illegal shift!");
1945 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1949 /* If we're doing geminate recompbination we usually don't need the distances.
1950 * Let's save some memory and time. */
1951 if (TRUE || !hb->bGem || hb->per->gemtype == gemAD)
1953 for (m = 0; (m <= hb0->nframes); m++)
1956 gtmp[mm] = is_hb(hb0->g[0], m);
1960 for (m = 0; (m <= hb1->nframes); m++)
1963 htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1964 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1965 if (hb->bGem /* && ptmp[mm] != 0 */)
1968 /* If this hbond has been seen before with donor and acceptor swapped,
1969 * then we need to find the mirrored (*-1) periodicity vector to truely
1970 * merge the hbond history. */
1971 pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper);
1972 /* Store index of mirror */
1973 if (pm > hb->per->nper)
1975 gmx_fatal(FARGS, "Illegal shift!");
1980 /* Reallocate target array */
1981 if (nnframes > hb0->maxframes)
1983 srenew(hb0->h[0], 4+nnframes/hb->wordlen);
1984 srenew(hb0->g[0], 4+nnframes/hb->wordlen);
1986 if (NULL != hb->per->pHist)
1988 clearPshift(&(hb->per->pHist[a1][a2]));
1991 /* Copy temp array to target array */
1992 for (m = 0; (m <= nnframes); m++)
1994 _set_hb(hb0->h[0], m, htmp[m]);
1995 _set_hb(hb0->g[0], m, gtmp[m]);
1998 addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
2002 /* Set scalar variables */
2004 hb0->maxframes = nnframes;
2007 /* Added argument bContact for nicer output.
2008 * Erik Marklund, June 29, 2006
2010 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
2012 int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp;
2013 unsigned int *htmp, *gtmp;
2018 indnew = hb->nrdist;
2020 /* Check whether donors are also acceptors */
2021 printf("Merging hbonds with Acceptor and Donor swapped\n");
2023 ntmp = 2*hb->max_frames;
2027 for (i = 0; (i < hb->d.nrd); i++)
2029 fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
2031 ii = hb->a.aptr[id];
2032 for (j = 0; (j < hb->a.nra); j++)
2035 jj = hb->d.dptr[ia];
2036 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
2037 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j]))))
2039 hb0 = hb->hbmap[i][j];
2040 hb1 = hb->hbmap[jj][ii];
2041 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
2043 do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j);
2044 if (ISHB(hb1->history[0]))
2048 else if (ISDIST(hb1->history[0]))
2055 gmx_incons("No contact history");
2059 gmx_incons("Neither hydrogen bond nor distance");
2065 clearPshift(&(hb->per->pHist[jj][ii]));
2069 hb1->history[0] = hbNo;
2074 fprintf(stderr, "\n");
2075 printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
2076 printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
2078 hb->nrdist = indnew;
2084 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
2086 int i, j, k, n_bound[MAXHH], nbtot;
2090 /* Set array to 0 */
2091 for (k = 0; (k < MAXHH); k++)
2095 /* Loop over possible donors */
2096 for (i = 0; (i < hb->d.nrd); i++)
2098 for (j = 0; (j < hb->d.nhydro[i]); j++)
2100 n_bound[hb->d.nhbonds[i][j]]++;
2103 fprintf(fp, "%12.5e", t);
2105 for (k = 0; (k < MAXHH); k++)
2107 fprintf(fp, " %8d", n_bound[k]);
2108 nbtot += n_bound[k]*k;
2110 fprintf(fp, " %8d\n", nbtot);
2113 /* Added argument bContact in do_hblife(...). Also
2114 * added support for -contact in function body.
2115 * - Erik Marklund, May 31, 2006 */
2116 /* Changed the contact code slightly.
2117 * - Erik Marklund, June 29, 2006
2119 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
2120 const output_env_t oenv)
2123 const char *leg[] = { "p(t)", "t p(t)" };
2125 int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
2126 int nframes = hb->nframes;
2129 double sum, integral;
2132 snew(h, hb->maxhydro);
2133 snew(histo, nframes+1);
2134 /* Total number of hbonds analyzed here */
2135 for (i = 0; (i < hb->d.nrd); i++)
2137 for (k = 0; (k < hb->a.nra); k++)
2139 hbh = hb->hbmap[i][k];
2157 for (m = 0; (m < hb->maxhydro); m++)
2161 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
2165 for (nh = 0; (nh < nhydro); nh++)
2170 /* Changed '<' into '<=' below, just like I
2171 did in the hbm-output-loop in the main code.
2172 - Erik Marklund, May 31, 2006
2174 for (j = 0; (j <= hbh->nframes); j++)
2176 ihb = is_hb(h[nh], j);
2177 if (debug && (ndump < 10))
2179 fprintf(debug, "%5d %5d\n", j, ihb);
2199 fprintf(stderr, "\n");
2202 fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
2206 fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
2210 xvgr_legend(fp, asize(leg), leg, oenv);
2212 while ((j0 > 0) && (histo[j0] == 0))
2217 for (i = 0; (i <= j0); i++)
2221 dt = hb->time[1]-hb->time[0];
2224 for (i = 1; (i <= j0); i++)
2226 t = hb->time[i] - hb->time[0] - 0.5*dt;
2227 x1 = t*histo[i]/sum;
2228 fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1);
2233 printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
2234 printf("Note that the lifetime obtained in this manner is close to useless\n");
2235 printf("Use the -ac option instead and check the Forward lifetime\n");
2236 please_cite(stdout, "Spoel2006b");
2241 /* Changed argument bMerge into oneHB to handle contacts properly.
2242 * - Erik Marklund, June 29, 2006
2244 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
2247 int i, j, k, m, nd, ihb, idist;
2248 int nframes = hb->nframes;
2256 fp = gmx_ffopen("debug-ac.xvg", "w");
2257 for (j = 0; (j < nframes); j++)
2259 fprintf(fp, "%10.3f", hb->time[j]);
2260 for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
2262 for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
2266 hbh = hb->hbmap[i][k];
2271 ihb = is_hb(hbh->h[0], j);
2272 idist = is_hb(hbh->g[0], j);
2278 for (m = 0; (m < hb->maxhydro) && !ihb; m++)
2280 ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m], j));
2281 idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
2283 /* This is not correct! */
2284 /* What isn't correct? -Erik M */
2289 fprintf(fp, " %1d-%1d", ihb, idist);
2299 static real calc_dg(real tau, real temp)
2310 return kbt*log(kbt*tau/PLANCK);
2315 int n0, n1, nparams, ndelta;
2317 real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
2321 #include <gsl/gsl_multimin.h>
2322 #include <gsl/gsl_sf.h>
2323 #include <gsl/gsl_version.h>
2325 static double my_f(const gsl_vector *v, void *params)
2327 t_luzar *tl = (t_luzar *)params;
2329 double tol = 1e-16, chi2 = 0;
2333 for (i = 0; (i < tl->nparams); i++)
2335 tl->kkk[i] = gsl_vector_get(v, i);
2340 for (i = tl->n0; (i < tl->n1); i += tl->ndelta)
2342 di = sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
2346 chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
2351 fprintf(stderr, "WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
2352 "di = %g k = %g kp = %g\n", tl->sigma_ct[i],
2353 tl->sigma_nt[i], tl->sigma_kt[i], di, k, kp);
2357 chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
2362 static real optimize_luzar_parameters(FILE *fp, t_luzar *tl, int maxiter,
2370 const gsl_multimin_fminimizer_type *T;
2371 gsl_multimin_fminimizer *s;
2374 gsl_multimin_function my_func;
2377 my_func.n = tl->nparams;
2378 my_func.params = (void *) tl;
2380 /* Starting point */
2381 x = gsl_vector_alloc (my_func.n);
2382 for (i = 0; (i < my_func.n); i++)
2384 gsl_vector_set (x, i, tl->kkk[i]);
2387 /* Step size, different for each of the parameters */
2388 dx = gsl_vector_alloc (my_func.n);
2389 for (i = 0; (i < my_func.n); i++)
2391 gsl_vector_set (dx, i, 0.01*tl->kkk[i]);
2394 T = gsl_multimin_fminimizer_nmsimplex;
2395 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
2397 gsl_multimin_fminimizer_set (s, &my_func, x, dx);
2398 gsl_vector_free (x);
2399 gsl_vector_free (dx);
2403 fprintf(fp, "%5s %12s %12s %12s %12s\n", "Iter", "k", "kp", "NM Size", "Chi2");
2409 status = gsl_multimin_fminimizer_iterate (s);
2413 gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s",
2414 gsl_multimin_fminimizer_name(s));
2417 d2 = gsl_multimin_fminimizer_minimum(s);
2418 size = gsl_multimin_fminimizer_size(s);
2419 status = gsl_multimin_test_size(size, tol);
2421 if (status == GSL_SUCCESS)
2425 fprintf(fp, "Minimum found using %s at:\n",
2426 gsl_multimin_fminimizer_name(s));
2432 fprintf(fp, "%5d", iter);
2433 for (i = 0; (i < my_func.n); i++)
2435 fprintf(fp, " %12.4e", gsl_vector_get (s->x, i));
2437 fprintf (fp, " %12.4e %12.4e\n", size, d2);
2440 while ((status == GSL_CONTINUE) && (iter < maxiter));
2442 gsl_multimin_fminimizer_free (s);
2447 static real quality_of_fit(real chi2, int N)
2449 return gsl_sf_gamma_inc_Q((N-2)/2.0, chi2/2.0);
2453 static real optimize_luzar_parameters(FILE gmx_unused *fp,
2454 t_luzar gmx_unused *tl,
2455 int gmx_unused maxiter,
2456 real gmx_unused tol)
2458 fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2462 static real quality_of_fit(real gmx_unused chi2, int gmx_unused N)
2464 fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2471 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
2472 real kt[], real sigma_ct[], real sigma_nt[],
2473 real sigma_kt[], real *k, real *kp,
2474 real *sigma_k, real *sigma_kp,
2480 real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
2485 for (i = 0; (i < n); i++)
2487 if (t[i] >= fit_start)
2500 tl.sigma_ct = sigma_ct;
2501 tl.sigma_nt = sigma_nt;
2502 tl.sigma_kt = sigma_kt;
2506 chi2 = optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2508 *kp = tl.kkk[1] = *kp;
2510 for (j = 0; (j < NK); j++)
2512 (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2515 kk2 += sqr(tl.kkk[0]);
2516 kp2 += sqr(tl.kkk[1]);
2519 *sigma_k = sqrt(kk2/NK - sqr(kkk/NK));
2520 *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2525 static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start,
2526 const output_env_t oenv)
2529 real e_1, fitparm[4];
2533 for (i = 0; (i < n); i++)
2549 do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(), effnEXP2, fitparm, 0);
2552 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
2553 real sigma_ct[], real sigma_nt[], real sigma_kt[],
2554 real fit_start, real temp, real smooth_tail_start,
2555 const output_env_t oenv)
2558 real k = 1, kp = 1, kow = 1;
2559 real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
2560 double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
2561 gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2563 if (smooth_tail_start >= 0)
2565 smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv);
2566 smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv);
2567 smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv);
2569 for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2575 for (i = i0; (i < n); i++)
2584 printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
2585 tmp = (sn2*sc2-sqr(scn));
2586 if ((tmp > 0) && (sn2 > 0))
2588 k = (sn2*sck-scn*snk)/tmp;
2589 kp = (k*scn-snk)/sn2;
2593 for (i = i0; (i < n); i++)
2595 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2597 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
2599 &sigma_k, &sigma_kp, fit_start);
2600 Q = quality_of_fit(chi2, 2);
2601 ddg = BOLTZ*temp*sigma_k/k;
2602 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2604 printf("The Rate and Delta G are followed by an error estimate\n");
2605 printf("----------------------------------------------------------\n"
2606 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
2607 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2608 k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
2609 ddg = BOLTZ*temp*sigma_kp/kp;
2610 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2611 kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
2616 for (i = i0; (i < n); i++)
2618 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2620 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2622 printf("--------------------------------------------------\n"
2623 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
2624 printf("Forward %10.3f %8.3f %10.3f %10g\n",
2625 k, 1/k, calc_dg(1/k, temp), chi2);
2626 printf("Backward %10.3f %8.3f %10.3f\n",
2627 kp, 1/kp, calc_dg(1/kp, temp));
2633 printf("One-way %10.3f %s%8.3f %10.3f\n",
2634 kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp));
2638 printf(" - Numerical problems computing HB thermodynamics:\n"
2639 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
2640 sc2, sn2, sk2, sck, snk, scn);
2642 /* Determine integral of the correlation function */
2643 tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
2644 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
2645 bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
2647 for (i = 0; (i < n-2); i++)
2649 if ((ct[i] > e_1) && (ct[i+1] <= e_1))
2656 /* Determine tau_relax from linear interpolation */
2657 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2658 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx,
2659 tau_rlx, bError ? " " : "",
2660 calc_dg(tau_rlx, temp));
2665 printf("Correlation functions too short to compute thermodynamics\n");
2669 void compute_derivative(int nn, real x[], real y[], real dydx[])
2673 /* Compute k(t) = dc(t)/dt */
2674 for (j = 1; (j < nn-1); j++)
2676 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2678 /* Extrapolate endpoints */
2679 dydx[0] = 2*dydx[1] - dydx[2];
2680 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2683 static void parallel_print(int *data, int nThreads)
2685 /* This prints the donors on which each tread is currently working. */
2688 fprintf(stderr, "\r");
2689 for (i = 0; i < nThreads; i++)
2691 fprintf(stderr, "%-7i", data[i]);
2695 static void normalizeACF(real *ct, real *gt, int nhb, int len)
2697 real ct_fac, gt_fac;
2700 /* Xu and Berne use the same normalization constant */
2703 gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
2705 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
2706 for (i = 0; i < len; i++)
2716 /* Added argument bContact in do_hbac(...). Also
2717 * added support for -contact in the actual code.
2718 * - Erik Marklund, May 31, 2006 */
2719 /* Changed contact code and added argument R2
2720 * - Erik Marklund, June 29, 2006
2722 static void do_hbac(const char *fn, t_hbdata *hb,
2723 int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
2724 real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv,
2725 const char *gemType, int nThreads,
2726 const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2729 int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
2730 const char *legNN[] = {
2734 static char **legGem;
2736 const char *legLuzar[] = {
2737 "Ac\\sfin sys\\v{}\\z{}(t)",
2739 "Cc\\scontact,hb\\v{}\\z{}(t)",
2740 "-dAc\\sfs\\v{}\\z{}/dt"
2742 gmx_bool bNorm = FALSE, bOMP = FALSE;
2745 real *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
2746 real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
2747 const real tol = 1e-3;
2748 int nframes = hb->nframes, nf;
2749 unsigned int **h = NULL, **g = NULL;
2750 int nh, nhbonds, nhydro, ngh;
2752 PSTYPE p, *pfound = NULL, np;
2754 int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
2755 real **rHbExGem = NULL;
2759 double *ctdouble, *timedouble, *fittedct;
2760 double fittolerance = 0.1;
2761 int *dondata = NULL, thisThread;
2764 AC_NONE, AC_NN, AC_GEM, AC_LUZAR
2773 printf("Doing autocorrelation ");
2775 /* Decide what kind of ACF calculations to do. */
2776 if (NN > NN_NONE && NN < NN_NR)
2778 #ifdef HAVE_NN_LOOPS
2780 printf("using the energy estimate.\n");
2783 printf("Can't do the NN-loop. Yet.\n");
2789 printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2791 nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0);
2792 snew(legGem, nSets);
2793 for (i = 0; i < nSets; i++)
2795 snew(legGem[i], 128);
2797 sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2800 sprintf(legGem[1], "Ac'(t)");
2804 sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2811 printf("according to the theory of Luzar and Chandler.\n");
2815 /* build hbexist matrix in reals for autocorr */
2816 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2818 while (n2 < nframes)
2825 if (acType != AC_NN || bOMP)
2827 snew(h, hb->maxhydro);
2828 snew(g, hb->maxhydro);
2831 /* Dump hbonds for debugging */
2832 dump_ac(hb, bMerge || bContact, nDump);
2834 /* Total number of hbonds analyzed here */
2839 if (acType != AC_LUZAR && bOMP)
2841 nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2843 gmx_omp_set_num_threads(nThreads);
2844 snew(dondata, nThreads);
2845 for (i = 0; i < nThreads; i++)
2849 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2850 "Expect close to linear scaling over this donor-loop.\n", nThreads);
2852 fprintf(stderr, "Donors: [thread no]\n");
2855 for (i = 0; i < nThreads; i++)
2857 snprintf(tmpstr, 7, "[%i]", i);
2858 fprintf(stderr, "%-7s", tmpstr);
2861 fprintf(stderr, "\n");
2865 /* Build the ACF according to acType */
2870 #ifdef HAVE_NN_LOOPS
2871 /* Here we're using the estimated energy for the hydrogen bonds. */
2874 #pragma omp parallel \
2875 private(i, j, k, nh, E, rhbex, thisThread) \
2879 thisThread = gmx_omp_get_thread_num();
2883 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2886 #pragma omp for schedule (dynamic)
2887 for (i = 0; i < hb->d.nrd; i++) /* loop over donors */
2891 #pragma omp critical
2893 dondata[thisThread] = i;
2894 parallel_print(dondata, nThreads);
2899 fprintf(stderr, "\r %i", i);
2902 for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */
2904 for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2906 E = hb->hbE.E[i][j][nh];
2909 for (k = 0; k < nframes; k++)
2911 if (E[k] != NONSENSE_E)
2913 rhbex[k] = (real)E[k];
2917 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0],
2918 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
2919 #pragma omp critical
2921 for (k = 0; (k < nn); k++)
2938 normalizeACF(ct, NULL, 0, nn);
2940 snew(timedouble, nn);
2941 for (j = 0; j < nn; j++)
2943 timedouble[j] = (double)(hb->time[j]);
2944 ctdouble[j] = (double)(ct[j]);
2947 /* Remove ballistic term */
2948 /* Ballistic component removal and fitting to the reversible geminate recombination model
2949 * will be taken out for the time being. First of all, one can remove the ballistic
2950 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2951 * problems with the robustness of the fitting to the model. More work is needed.
2952 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2953 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2954 * a BSD-licence that can do the job.
2956 * - Erik Marklund, June 18 2010.
2958 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2959 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2961 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2962 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2964 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)");
2965 xvgr_legend(fp, asize(legNN), legNN);
2967 for (j = 0; (j < nn); j++)
2969 fprintf(fp, "%10g %10g %10g\n",
2970 hb->time[j]-hb->time[0],
2978 #endif /* HAVE_NN_LOOPS */
2979 break; /* case AC_NN */
2983 memset(ct, 0, 2*n2*sizeof(real));
2985 fprintf(stderr, "Donor:\n");
2988 #define __ACDATA p_ct
2991 #pragma omp parallel \
2992 private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2993 pfound, poff, rHbExGem, p, ihb, mMax, \
2996 { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
2999 thisThread = gmx_omp_get_thread_num();
3000 snew(h, hb->maxhydro);
3001 snew(g, hb->maxhydro);
3008 memset(p_ct, 0, 2*n2*sizeof(real));
3010 /* I'm using a chunk size of 1, since I expect \
3011 * the overhead to be really small compared \
3012 * to the actual calculations \ */
3013 #pragma omp for schedule(dynamic,1) nowait
3014 for (i = 0; i < hb->d.nrd; i++)
3019 #pragma omp critical
3021 dondata[thisThread] = i;
3022 parallel_print(dondata, nThreads);
3027 fprintf(stderr, "\r %i", i);
3029 for (k = 0; k < hb->a.nra; k++)
3031 for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++)
3033 hbh = hb->hbmap[i][k];
3036 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
3037 * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
3038 * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
3039 pHist = &(hb->per->pHist[i][k]);
3040 if (ISHB(hbh->history[nh]) && pHist->len != 0)
3045 g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL;
3049 /* count the number of periodic shifts encountered and store
3050 * them in separate arrays. */
3052 for (j = 0; j < pHist->len; j++)
3055 for (m = 0; m <= np; m++)
3057 if (m == np) /* p not recognized in list. Add it and set up new array. */
3060 if (np > hb->per->nper)
3062 gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
3064 if (m >= mMax) /* Extend the arrays.
3065 * Doing it like this, using mMax to keep track of the sizes,
3066 * eleviates the need for freeing and re-allocating the arrays
3067 * when taking on the next donor-acceptor pair */
3070 srenew(pfound, np); /* The list of found periodic shifts. */
3071 srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */
3072 snew(rHbExGem[m], 2*n2);
3077 if (rHbExGem != NULL && rHbExGem[m] != NULL)
3079 /* This must be done, as this array was most likey
3080 * used to store stuff in some previous iteration. */
3081 memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
3085 fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
3097 } /* m: Loop over found shifts */
3098 } /* j: Loop over shifts */
3100 /* Now unpack and disentangle the existence funtions. */
3101 for (j = 0; j < nf; j++)
3108 * pfound: list of periodic shifts found for this pair.
3109 * poff: list of frame offsets; that is, the first
3110 * frame a hbond has a particular periodic shift. */
3111 p = getPshift(*pHist, j+n0);
3114 for (m = 0; m < np; m++)
3122 gmx_fatal(FARGS, "Shift not found, but must be there.");
3126 ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j));
3131 poff[m] = j; /* Here's where the first hbond with shift p is,
3132 * relative to the start of h[0].*/
3136 gmx_fatal(FARGS, "j<poff[m]");
3138 rHbExGem[m][j-poff[m]] += 1;
3143 /* Now, build ac. */
3144 for (m = 0; m < np; m++)
3146 if (rHbExGem[m][0] > 0 && n0+poff[m] < nn /* && m==0 */)
3148 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0],
3149 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3150 for (j = 0; (j < nn); j++)
3152 __ACDATA[j] += rHbExGem[m][j];
3155 } /* Building of ac. */
3158 } /* hydrogen loop */
3159 } /* acceptor loop */
3162 for (m = 0; m <= mMax; m++)
3175 #pragma omp critical
3177 for (i = 0; i < nn; i++)
3185 } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
3191 normalizeACF(ct, NULL, 0, nn);
3193 fprintf(stderr, "\n\nACF successfully calculated.\n");
3195 /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
3198 snew(timedouble, nn);
3201 for (j = 0; j < nn; j++)
3203 timedouble[j] = (double)(hb->time[j]);
3204 ctdouble[j] = (double)(ct[j]);
3207 /* Remove ballistic term */
3208 /* Ballistic component removal and fitting to the reversible geminate recombination model
3209 * will be taken out for the time being. First of all, one can remove the ballistic
3210 * component with g_analyze afterwards. Secondly, and more importantly, there are still
3211 * problems with the robustness of the fitting to the model. More work is needed.
3212 * A third reason is that we're currently using gsl for this and wish to reduce dependence
3213 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
3214 * a BSD-licence that can do the job.
3216 * - Erik Marklund, June 18 2010.
3218 /* if (bBallistic) { */
3219 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
3220 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
3222 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
3223 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
3226 /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
3231 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3235 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3237 xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv);
3239 for (j = 0; (j < nn); j++)
3241 fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]);
3244 fprintf(fp, " %10g", ctdouble[j]);
3248 fprintf(fp, " %10g", fittedct[j]);
3259 break; /* case AC_GEM */
3272 for (i = 0; (i < hb->d.nrd); i++)
3274 for (k = 0; (k < hb->a.nra); k++)
3277 hbh = hb->hbmap[i][k];
3281 if (bMerge || bContact)
3283 if (ISHB(hbh->history[0]))
3292 for (m = 0; (m < hb->maxhydro); m++)
3294 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
3296 g[nhydro] = hbh->g[m];
3297 h[nhydro] = hbh->h[m];
3304 for (nh = 0; (nh < nhydro); nh++)
3306 int nrint = bContact ? hb->nrdist : hb->nrhb;
3307 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
3309 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
3312 for (j = 0; (j < nframes); j++)
3314 /* Changed '<' into '<=' below, just like I did in
3315 the hbm-output-loop in the gmx_hbond() block.
3316 - Erik Marklund, May 31, 2006 */
3319 ihb = is_hb(h[nh], j);
3320 idist = is_hb(g[nh], j);
3327 /* For contacts: if a second cut-off is provided, use it,
3328 * otherwise use g(t) = 1-h(t) */
3329 if (!R2 && bContact)
3335 gt[j] = idist*(1-ihb);
3342 /* The autocorrelation function is normalized after summation only */
3343 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
3344 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3346 /* Cross correlation analysis for thermodynamics */
3347 for (j = nframes; (j < n2); j++)
3353 cross_corr(n2, ht, gt, dght);
3355 for (j = 0; (j < nn); j++)
3364 fprintf(stderr, "\n");
3367 normalizeACF(ct, ght, nhb, nn);
3369 /* Determine tail value for statistics */
3372 for (j = nn/2; (j < nn); j++)
3375 tail2 += ct[j]*ct[j];
3377 tail /= (nn - nn/2);
3378 tail2 /= (nn - nn/2);
3379 dtail = sqrt(tail2-tail*tail);
3381 /* Check whether the ACF is long enough */
3384 printf("\nWARNING: Correlation function is probably not long enough\n"
3385 "because the standard deviation in the tail of C(t) > %g\n"
3386 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
3389 for (j = 0; (j < nn); j++)
3392 ct[j] = (cct[j]-tail)/(1-tail);
3394 /* Compute negative derivative k(t) = -dc(t)/dt */
3395 compute_derivative(nn, hb->time, ct, kt);
3396 for (j = 0; (j < nn); j++)
3404 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3408 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3410 xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
3413 for (j = 0; (j < nn); j++)
3415 fprintf(fp, "%10g %10g %10g %10g %10g\n",
3416 hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
3420 analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
3421 fit_start, temp, smooth_tail_start, oenv);
3423 do_view(oenv, fn, NULL);
3435 break; /* case AC_LUZAR */
3438 gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
3439 } /* switch (acType) */
3442 static void init_hbframe(t_hbdata *hb, int nframes, real t)
3446 hb->time[nframes] = t;
3447 hb->nhb[nframes] = 0;
3448 hb->ndist[nframes] = 0;
3449 for (i = 0; (i < max_hx); i++)
3451 hb->nhx[nframes][i] = 0;
3453 /* Loop invalidated */
3454 if (hb->bHBmap && 0)
3456 for (i = 0; (i < hb->d.nrd); i++)
3458 for (j = 0; (j < hb->a.nra); j++)
3460 for (m = 0; (m < hb->maxhydro); m++)
3462 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
3464 set_hb(hb, i, m, j, nframes, HB_NO);
3470 /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
3473 static void analyse_donor_props(const char *fn, t_hbdata *hb, int nframes, real t,
3474 const output_env_t oenv)
3476 static FILE *fp = NULL;
3477 const char *leg[] = { "Nbound", "Nfree" };
3478 int i, j, k, nbound, nb, nhtot;
3486 fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3487 xvgr_legend(fp, asize(leg), leg, oenv);
3491 for (i = 0; (i < hb->d.nrd); i++)
3493 for (k = 0; (k < hb->d.nhydro[i]); k++)
3497 for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
3499 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
3500 is_hb(hb->hbmap[i][j]->h[k], nframes))
3508 fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound);
3511 static void dump_hbmap(t_hbdata *hb,
3512 int nfile, t_filenm fnm[], gmx_bool bTwo,
3513 gmx_bool bContact, int isize[], int *index[], char *grpnames[],
3517 int ddd, hhh, aaa, i, j, k, m, grp;
3518 char ds[32], hs[32], as[32];
3521 fp = opt2FILE("-hbn", nfile, fnm, "w");
3522 if (opt2bSet("-g", nfile, fnm))
3524 fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
3525 fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
3531 for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
3533 fprintf(fp, "[ %s ]", grpnames[grp]);
3534 for (i = 0; i < isize[grp]; i++)
3536 fprintf(fp, (i%15) ? " " : "\n");
3537 fprintf(fp, " %4u", index[grp][i]+1);
3541 Added -contact support below.
3542 - Erik Marklund, May 29, 2006
3546 fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
3547 for (i = 0; (i < hb->d.nrd); i++)
3549 if (hb->d.grp[i] == grp)
3551 for (j = 0; (j < hb->d.nhydro[i]); j++)
3553 fprintf(fp, " %4u %4u", hb->d.don[i]+1,
3554 hb->d.hydro[i][j]+1);
3560 fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
3561 for (i = 0; (i < hb->a.nra); i++)
3563 if (hb->a.grp[i] == grp)
3565 fprintf(fp, (i%15 && !first) ? " " : "\n");
3566 fprintf(fp, " %4u", hb->a.acc[i]+1);
3575 fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
3576 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
3580 fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
3583 for (i = 0; (i < hb->d.nrd); i++)
3586 for (k = 0; (k < hb->a.nra); k++)
3589 for (m = 0; (m < hb->d.nhydro[i]); m++)
3591 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
3593 sprintf(ds, "%s", mkatomname(atoms, ddd));
3594 sprintf(as, "%s", mkatomname(atoms, aaa));
3597 fprintf(fp, " %6u %6u\n", ddd+1, aaa+1);
3600 fprintf(fplog, "%12s %12s\n", ds, as);
3605 hhh = hb->d.hydro[i][m];
3606 sprintf(hs, "%s", mkatomname(atoms, hhh));
3607 fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1);
3610 fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
3624 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3625 * It mimics add_frames() and init_frame() to some extent. */
3626 static void sync_hbdata(t_hbdata *p_hb, int nframes)
3629 if (nframes >= p_hb->max_frames)
3631 p_hb->max_frames += 4096;
3632 srenew(p_hb->nhb, p_hb->max_frames);
3633 srenew(p_hb->ndist, p_hb->max_frames);
3634 srenew(p_hb->n_bound, p_hb->max_frames);
3635 srenew(p_hb->nhx, p_hb->max_frames);
3638 srenew(p_hb->danr, p_hb->max_frames);
3640 memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3641 memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3642 p_hb->nhb[nframes] = 0;
3643 p_hb->ndist[nframes] = 0;
3646 p_hb->nframes = nframes;
3649 /* p_hb->nhx[nframes][i] */
3651 memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
3653 /* hb->per will remain constant througout the frame loop,
3654 * even though the data its members point to will change,
3655 * hence no need for re-syncing. */
3658 int gmx_hbond(int argc, char *argv[])
3660 const char *desc[] = {
3661 "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3662 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
3663 "(zero is extended) and the distance Donor - Acceptor",
3664 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
3665 "OH and NH groups are regarded as donors, O is an acceptor always,",
3666 "N is an acceptor by default, but this can be switched using",
3667 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3668 "to the first preceding non-hydrogen atom.[PAR]",
3670 "You need to specify two groups for analysis, which must be either",
3671 "identical or non-overlapping. All hydrogen bonds between the two",
3672 "groups are analyzed.[PAR]",
3674 "If you set [TT]-shell[tt], you will be asked for an additional index group",
3675 "which should contain exactly one atom. In this case, only hydrogen",
3676 "bonds between atoms within the shell distance from the one atom are",
3679 "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3680 "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3681 "If contact kinetics are analyzed by using the -contact option, then",
3682 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3683 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3684 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3685 "See mentioned literature for more details and definitions."
3688 /* "It is also possible to analyse specific hydrogen bonds with",
3689 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3690 "Donor Hydrogen Acceptor, in the following way:[PAR]",
3698 "Note that the triplets need not be on separate lines.",
3699 "Each atom triplet specifies a hydrogen bond to be analyzed,",
3700 "note also that no check is made for the types of atoms.[PAR]",
3702 "[BB]Output:[bb][BR]",
3703 "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
3704 "[TT]-ac[tt]: average over all autocorrelations of the existence",
3705 "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3706 "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3707 "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
3708 "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
3709 "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3710 "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3711 "with helices in proteins.[BR]",
3712 "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
3713 "for selected groups, all hydrogen bonded atoms from all groups and",
3714 "all solvent atoms involved in insertion.[BR]",
3715 "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
3716 "frames, this also contains information on solvent insertion",
3717 "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3719 "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3720 "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3721 "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3722 "compare results to Raman Spectroscopy.",
3724 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3725 "require an amount of memory proportional to the total numbers of donors",
3726 "times the total number of acceptors in the selected group(s)."
3729 static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
3730 static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1;
3731 static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
3732 static int nDump = 0, nFitPoints = 100;
3733 static int nThreads = 0, nBalExp = 4;
3735 static gmx_bool bContact = FALSE, bBallistic = FALSE, bGemFit = FALSE;
3736 static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
3737 static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3741 { "-a", FALSE, etREAL, {&acut},
3742 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
3743 { "-r", FALSE, etREAL, {&rcut},
3744 "Cutoff radius (nm, X - Acceptor, see next option)" },
3745 { "-da", FALSE, etBOOL, {&bDA},
3746 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3747 { "-r2", FALSE, etREAL, {&r2cut},
3748 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3749 { "-abin", FALSE, etREAL, {&abin},
3750 "Binwidth angle distribution (degrees)" },
3751 { "-rbin", FALSE, etREAL, {&rbin},
3752 "Binwidth distance distribution (nm)" },
3753 { "-nitacc", FALSE, etBOOL, {&bNitAcc},
3754 "Regard nitrogen atoms as acceptors" },
3755 { "-contact", FALSE, etBOOL, {&bContact},
3756 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3757 { "-shell", FALSE, etREAL, {&rshell},
3758 "when > 0, only calculate hydrogen bonds within # nm shell around "
3760 { "-fitstart", FALSE, etREAL, {&fit_start},
3761 "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]" },
3762 { "-fitend", FALSE, etREAL, {&fit_end},
3763 "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])" },
3764 { "-temp", FALSE, etREAL, {&temp},
3765 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3766 { "-smooth", FALSE, etREAL, {&smooth_tail_start},
3767 "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3768 { "-dump", FALSE, etINT, {&nDump},
3769 "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
3770 { "-max_hb", FALSE, etREAL, {&maxnhb},
3771 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3772 { "-merge", FALSE, etBOOL, {&bMerge},
3773 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3774 { "-geminate", FALSE, etENUM, {gemType},
3775 "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3776 { "-diff", FALSE, etREAL, {&D},
3777 "Dffusion 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."},
3779 { "-nthreads", FALSE, etINT, {&nThreads},
3780 "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
3783 const char *bugs[] = {
3784 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3787 { efTRX, "-f", NULL, ffREAD },
3788 { efTPX, NULL, NULL, ffREAD },
3789 { efNDX, NULL, NULL, ffOPTRD },
3790 /* { efNDX, "-sel", "select", ffOPTRD },*/
3791 { efXVG, "-num", "hbnum", ffWRITE },
3792 { efLOG, "-g", "hbond", ffOPTWR },
3793 { efXVG, "-ac", "hbac", ffOPTWR },
3794 { efXVG, "-dist", "hbdist", ffOPTWR },
3795 { efXVG, "-ang", "hbang", ffOPTWR },
3796 { efXVG, "-hx", "hbhelix", ffOPTWR },
3797 { efNDX, "-hbn", "hbond", ffOPTWR },
3798 { efXPM, "-hbm", "hbmap", ffOPTWR },
3799 { efXVG, "-don", "donor", ffOPTWR },
3800 { efXVG, "-dan", "danum", ffOPTWR },
3801 { efXVG, "-life", "hblife", ffOPTWR },
3802 { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
3805 #define NFILE asize(fnm)
3807 char hbmap [HB_NR] = { ' ', 'o', '-', '*' };
3808 const char *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
3809 t_rgb hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
3811 t_trxstatus *status;
3816 int npargs, natoms, nframes = 0, shatom;
3822 real t, ccut, dist = 0.0, ang = 0.0;
3823 double max_nhb, aver_nhb, aver_dist;
3824 int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
3826 int xj, yj, zj, aj, xjj, yjj, zjj;
3827 int xk, yk, zk, ak, xkk, ykk, zkk;
3828 gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
3829 int *adist, *rdist, *aptr, *rprt;
3830 int grp, nabin, nrbin, bin, resdist, ihb;
3832 t_hbdata *hb, *hbptr;
3833 FILE *fp, *fpins = NULL, *fpnhb = NULL;
3835 t_ncell *icell, *jcell, *kcell;
3837 unsigned char *datable;
3842 int ii, jj, hh, actual_nThreads;
3844 gmx_bool bGem, bNN, bParallel;
3845 t_gemParams *params = NULL;
3846 gmx_bool bEdge_yjj, bEdge_xjj, bOMP;
3848 t_hbdata **p_hb = NULL; /* one per thread, then merge after the frame loop */
3849 int **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
3858 ppa = add_acf_pargs(&npargs, pa);
3860 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE, NFILE, fnm, npargs,
3861 ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
3866 /* NN-loop? If so, what estimator to use ?*/
3868 /* Outcommented for now DvdS 2010-07-13
3869 while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3872 gmx_fatal(FARGS, "Invalid NN-loop type.");
3875 for (i = 2; bNN == FALSE && i < NN_NR; i++)
3877 bNN = bNN || NN == i;
3880 if (NN > NN_NONE && bMerge)
3885 /* geminate recombination? If so, which flavor? */
3887 while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0)
3891 if (gemmode == gemNR)
3893 gmx_fatal(FARGS, "Invalid recombination type.");
3897 for (i = 2; bGem == FALSE && i < gemNR; i++)
3899 bGem = bGem || gemmode == i;
3904 printf("Geminate recombination: %s\n", gemType[gemmode]);
3906 printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
3910 if (gemmode != gemDD)
3912 printf("Turning off -contact option...\n");
3918 if (gemmode == gemDD)
3920 printf("Turning on -contact option...\n");
3926 if (gemmode == gemAA)
3928 printf("Turning off -merge option...\n");
3934 if (gemmode != gemAA)
3936 printf("Turning on -merge option...\n");
3944 ccut = cos(acut*DEG2RAD);
3950 gmx_fatal(FARGS, "Can not analyze selected contacts.");
3954 gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
3958 /* Initiate main data structure! */
3959 bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
3960 opt2bSet("-life", NFILE, fnm) ||
3961 opt2bSet("-hbn", NFILE, fnm) ||
3962 opt2bSet("-hbm", NFILE, fnm) ||
3965 if (opt2bSet("-nhbdist", NFILE, fnm))
3967 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3968 fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
3969 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
3970 xvgr_legend(fpnhb, asize(leg), leg, oenv);
3973 hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
3976 read_tpx_top(ftp2fn(efTPX, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
3978 snew(grpnames, grNR);
3981 /* Make Donor-Acceptor table */
3982 snew(datable, top.atoms.nr);
3983 gen_datable(index[0], isize[0], datable, top.atoms.nr);
3987 /* analyze selected hydrogen bonds */
3988 printf("Select group with selected atoms:\n");
3989 get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
3990 1, &nsel, index, grpnames);
3993 gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
3994 "and therefore cannot contain triplets of "
3995 "Donor-Hydrogen-Acceptor", grpnames[0]);
3999 for (i = 0; (i < nsel); i += 3)
4001 int dd = index[0][i];
4002 int aa = index[0][i+2];
4003 /* int */ hh = index[0][i+1];
4004 add_dh (&hb->d, dd, hh, i, datable);
4005 add_acc(&hb->a, aa, i);
4006 /* Should this be here ? */
4007 snew(hb->d.dptr, top.atoms.nr);
4008 snew(hb->a.aptr, top.atoms.nr);
4009 add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri);
4011 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
4012 isize[0], grpnames[0]);
4016 /* analyze all hydrogen bonds: get group(s) */
4017 printf("Specify 2 groups to analyze:\n");
4018 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4019 2, isize, index, grpnames);
4021 /* check if we have two identical or two non-overlapping groups */
4022 bTwo = isize[0] != isize[1];
4023 for (i = 0; (i < isize[0]) && !bTwo; i++)
4025 bTwo = index[0][i] != index[1][i];
4029 printf("Checking for overlap in atoms between %s and %s\n",
4030 grpnames[0], grpnames[1]);
4031 for (i = 0; i < isize[1]; i++)
4033 if (ISINGRP(datable[index[1][i]]))
4035 gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
4036 grpnames[0], grpnames[1]);
4040 printf("Checking for overlap in atoms between %s and %s\n",
4041 grpnames[0],grpnames[1]);
4042 for (i=0; i<isize[0]; i++)
4043 for (j=0; j<isize[1]; j++)
4044 if (index[0][i] == index[1][j])
4045 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
4046 grpnames[0],grpnames[1]);
4051 printf("Calculating %s "
4052 "between %s (%d atoms) and %s (%d atoms)\n",
4053 bContact ? "contacts" : "hydrogen bonds",
4054 grpnames[0], isize[0], grpnames[1], isize[1]);
4058 fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
4059 bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
4064 /* search donors and acceptors in groups */
4065 snew(datable, top.atoms.nr);
4066 for (i = 0; (i < grNR); i++)
4068 if ( ((i == gr0) && !bSelected ) ||
4069 ((i == gr1) && bTwo ))
4071 gen_datable(index[i], isize[i], datable, top.atoms.nr);
4074 search_acceptors(&top, isize[i], index[i], &hb->a, i,
4075 bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
4076 search_donors (&top, isize[i], index[i], &hb->d, i,
4077 TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
4081 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
4082 search_donors (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
4086 clear_datable_grp(datable, top.atoms.nr);
4091 printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
4093 snew(donors[gr0D], dons[gr0D].nrd);*/
4097 printf("Making hbmap structure...");
4098 /* Generate hbond data structure */
4103 #ifdef HAVE_NN_LOOPS
4112 printf("Making per structure...");
4113 /* Generate hbond data structure */
4120 if (hb->d.nrd + hb->a.nra == 0)
4122 printf("No Donors or Acceptors found\n");
4129 printf("No Donors found\n");
4134 printf("No Acceptors found\n");
4140 gmx_fatal(FARGS, "Nothing to be done");
4149 /* get index group with atom for shell */
4152 printf("Select atom for shell (1 atom):\n");
4153 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4154 1, &shisz, &shidx, &shgrpnm);
4157 printf("group contains %d atoms, should be 1 (one)\n", shisz);
4162 printf("Will calculate hydrogen bonds within a shell "
4163 "of %g nm around atom %i\n", rshell, shatom+1);
4166 /* Analyze trajectory */
4167 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
4168 if (natoms > top.atoms.nr)
4170 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
4171 top.atoms.nr, natoms);
4174 bBox = ir.ePBC != epbcNONE;
4175 grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
4178 snew(adist, nabin+1);
4179 snew(rdist, nrbin+1);
4183 gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
4189 #define __ADIST adist
4190 #define __RDIST rdist
4192 #else /* GMX_OPENMP ================================================== \
4193 * Set up the OpenMP stuff, |
4194 * like the number of threads and such |
4195 * Also start the parallel loop. |
4197 #define __ADIST p_adist[threadNr]
4198 #define __RDIST p_rdist[threadNr]
4199 #define __HBDATA p_hb[threadNr]
4203 bParallel = !bSelected;
4207 actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
4209 gmx_omp_set_num_threads(actual_nThreads);
4210 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
4215 actual_nThreads = 1;
4218 snew(p_hb, actual_nThreads);
4219 snew(p_adist, actual_nThreads);
4220 snew(p_rdist, actual_nThreads);
4221 for (i = 0; i < actual_nThreads; i++)
4224 snew(p_adist[i], nabin+1);
4225 snew(p_rdist[i], nrbin+1);
4227 p_hb[i]->max_frames = 0;
4228 p_hb[i]->nhb = NULL;
4229 p_hb[i]->ndist = NULL;
4230 p_hb[i]->n_bound = NULL;
4231 p_hb[i]->time = NULL;
4232 p_hb[i]->nhx = NULL;
4234 p_hb[i]->bHBmap = hb->bHBmap;
4235 p_hb[i]->bDAnr = hb->bDAnr;
4236 p_hb[i]->bGem = hb->bGem;
4237 p_hb[i]->wordlen = hb->wordlen;
4238 p_hb[i]->nframes = hb->nframes;
4239 p_hb[i]->maxhydro = hb->maxhydro;
4240 p_hb[i]->danr = hb->danr;
4243 p_hb[i]->hbmap = hb->hbmap;
4244 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
4245 p_hb[i]->per = hb->per;
4247 #ifdef HAVE_NN_LOOPS
4248 p_hb[i]->hbE = hb->hbE;
4252 p_hb[i]->nrdist = 0;
4256 /* Make a thread pool here,
4257 * instead of forking anew at every frame. */
4259 #pragma omp parallel \
4261 private(j, h, ii, jj, hh, E, \
4262 xi, yi, zi, xj, yj, zj, threadNr, \
4263 dist, ang, peri, icell, jcell, \
4264 grp, ogrp, ai, aj, xjj, yjj, zjj, \
4265 xk, yk, zk, ihb, id, resdist, \
4266 xkk, ykk, zkk, kcell, ak, k, bTric, \
4267 bEdge_xjj, bEdge_yjj) \
4269 { /* Start of parallel region */
4270 threadNr = gmx_omp_get_thread_num();
4275 bTric = bBox && TRICLINIC(box);
4279 sync_hbdata(p_hb[threadNr], nframes);
4283 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
4284 rshell, ngrid, grid);
4285 reset_nhbonds(&(hb->d));
4287 if (debug && bDebug)
4289 dump_grid(debug, ngrid, grid);
4292 add_frames(hb, nframes);
4293 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
4297 count_da_grid(ngrid, grid, hb->danr[nframes]);
4303 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
4308 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
4309 /* Loop over all atom pairs and estimate interaction energy */
4313 addFramesNN(hb, nframes);
4317 #pragma omp for schedule(dynamic)
4318 for (i = 0; i < hb->d.nrd; i++)
4320 for (j = 0; j < hb->a.nra; j++)
4323 h < (bContact ? 1 : hb->d.nhydro[i]);
4326 if (i == hb->d.nrd || j == hb->a.nra)
4328 gmx_fatal(FARGS, "out of bounds");
4331 /* Get the real atom ids */
4334 hh = hb->d.hydro[i][h];
4336 /* Estimate the energy from the geometry */
4337 E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
4338 /* Store the energy */
4339 storeHbEnergy(hb, i, j, h, E, nframes);
4343 #endif /* HAVE_NN_LOOPS */
4352 /* Do not parallelize this just yet. */
4354 for (ii = 0; (ii < nsel); ii++)
4356 int dd = index[0][i];
4357 int aa = index[0][i+2];
4358 /* int */ hh = index[0][i+1];
4359 ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
4360 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4364 /* add to index if not already there */
4366 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri);
4370 } /* if (bSelected) */
4378 calcBoxProjection(box, hb->per->P);
4381 /* loop over all gridcells (xi,yi,zi) */
4382 /* Removed confusing macro, DvdS 27/12/98 */
4385 /* The outer grid loop will have to do for now. */
4386 #pragma omp for schedule(dynamic)
4387 for (xi = 0; xi < ngrid[XX]; xi++)
4389 for (yi = 0; (yi < ngrid[YY]); yi++)
4391 for (zi = 0; (zi < ngrid[ZZ]); zi++)
4394 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
4395 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
4397 icell = &(grid[zi][yi][xi].d[grp]);
4408 /* loop over all hydrogen atoms from group (grp)
4409 * in this gridcell (icell)
4411 for (ai = 0; (ai < icell->nr); ai++)
4413 i = icell->atoms[ai];
4415 /* loop over all adjacent gridcells (xj,yj,zj) */
4416 for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
4417 zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
4420 zj = grid_mod(zjj, ngrid[ZZ]);
4421 bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
4422 for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
4423 yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
4426 yj = grid_mod(yjj, ngrid[YY]);
4428 (yj == 0) || (yj == ngrid[YY] - 1) ||
4429 (zj == 0) || (zj == ngrid[ZZ] - 1);
4430 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
4431 xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
4434 xj = grid_mod(xjj, ngrid[XX]);
4435 jcell = &(grid[zj][yj][xj].a[ogrp]);
4436 /* loop over acceptor atoms from other group (ogrp)
4437 * in this adjacent gridcell (jcell)
4439 for (aj = 0; (aj < jcell->nr); aj++)
4441 j = jcell->atoms[aj];
4443 /* check if this once was a h-bond */
4445 ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
4446 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4450 /* add to index if not already there */
4452 add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri);
4454 /* make angle and distance distributions */
4455 if (ihb == hbHB && !bContact)
4459 gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
4462 __ADIST[(int)( ang/abin)]++;
4463 __RDIST[(int)(dist/rbin)]++;
4467 if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
4469 gmx_fatal(FARGS, "Invalid donor %d", i);
4471 if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
4473 gmx_fatal(FARGS, "Invalid acceptor %d", j);
4475 resdist = abs(top.atoms.atom[i].resind-
4476 top.atoms.atom[j].resind);
4477 if (resdist >= max_hx)
4481 __HBDATA->nhx[nframes][resdist]++;
4492 } /* for xi,yi,zi */
4495 } /* if (bSelected) {...} else */
4498 /* Better wait for all threads to finnish using x[] before updating it. */
4501 #pragma omp critical
4503 /* Sum up histograms and counts from p_hb[] into hb */
4506 hb->nhb[k] += p_hb[threadNr]->nhb[k];
4507 hb->ndist[k] += p_hb[threadNr]->ndist[k];
4508 for (j = 0; j < max_hx; j++)
4510 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
4515 /* Here are a handful of single constructs
4516 * to share the workload a bit. The most
4517 * important one is of course the last one,
4518 * where there's a potential bottleneck in form
4525 analyse_donor_props(opt2fn_null("-don", NFILE, fnm), hb, k, t, oenv);
4533 do_nhb_dist(fpnhb, hb, t);
4536 } /* if (bNN) {...} else + */
4540 trrStatus = (read_next_x(oenv, status, &t, x, box));
4550 #pragma omp critical
4552 hb->nrhb += p_hb[threadNr]->nrhb;
4553 hb->nrdist += p_hb[threadNr]->nrdist;
4555 /* Free parallel datastructures */
4556 sfree(p_hb[threadNr]->nhb);
4557 sfree(p_hb[threadNr]->ndist);
4558 sfree(p_hb[threadNr]->nhx);
4561 for (i = 0; i < nabin; i++)
4563 for (j = 0; j < actual_nThreads; j++)
4566 adist[i] += p_adist[j][i];
4570 for (i = 0; i <= nrbin; i++)
4572 for (j = 0; j < actual_nThreads; j++)
4574 rdist[i] += p_rdist[j][i];
4578 sfree(p_adist[threadNr]);
4579 sfree(p_rdist[threadNr]);
4581 } /* End of parallel region */
4588 if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
4590 gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
4593 free_grid(ngrid, &grid);
4601 /* Compute maximum possible number of different hbonds */
4608 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
4610 /* Added support for -contact below.
4611 * - Erik Marklund, May 29-31, 2006 */
4612 /* Changed contact code.
4613 * - Erik Marklund, June 29, 2006 */
4618 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
4622 printf("Found %d different %s in trajectory\n"
4623 "Found %d different atom-pairs within %s distance\n",
4624 hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
4625 hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
4627 /*Control the pHist.*/
4631 merge_hb(hb, bTwo, bContact);
4634 if (opt2bSet("-hbn", NFILE, fnm))
4636 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
4639 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
4640 * to make the -hbn and -hmb output match eachother.
4641 * - Erik Marklund, May 30, 2006 */
4644 /* Print out number of hbonds and distances */
4647 fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
4648 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
4650 snew(leg[0], STRLEN);
4651 snew(leg[1], STRLEN);
4652 sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
4653 sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
4654 xvgr_legend(fp, 2, (const char**)leg, oenv);
4658 for (i = 0; (i < nframes); i++)
4660 fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
4661 aver_nhb += hb->nhb[i];
4662 aver_dist += hb->ndist[i];
4665 aver_nhb /= nframes;
4666 aver_dist /= nframes;
4667 /* Print HB distance distribution */
4668 if (opt2bSet("-dist", NFILE, fnm))
4673 for (i = 0; i < nrbin; i++)
4678 fp = xvgropen(opt2fn("-dist", NFILE, fnm),
4679 "Hydrogen Bond Distribution",
4681 "Donor - Acceptor Distance (nm)" :
4682 "Hydrogen - Acceptor Distance (nm)", "", oenv);
4683 for (i = 0; i < nrbin; i++)
4685 fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
4690 /* Print HB angle distribution */
4691 if (opt2bSet("-ang", NFILE, fnm))
4696 for (i = 0; i < nabin; i++)
4701 fp = xvgropen(opt2fn("-ang", NFILE, fnm),
4702 "Hydrogen Bond Distribution",
4703 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
4704 for (i = 0; i < nabin; i++)
4706 fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
4711 /* Print HB in alpha-helix */
4712 if (opt2bSet("-hx", NFILE, fnm))
4714 fp = xvgropen(opt2fn("-hx", NFILE, fnm),
4715 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4716 xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
4717 for (i = 0; i < nframes; i++)
4719 fprintf(fp, "%10g", hb->time[i]);
4720 for (j = 0; j < max_hx; j++)
4722 fprintf(fp, " %6d", hb->nhx[i][j]);
4730 printf("Average number of %s per timeframe %.3f out of %g possible\n",
4731 bContact ? "contacts" : "hbonds",
4732 bContact ? aver_dist : aver_nhb, max_nhb);
4735 /* Do Autocorrelation etc. */
4739 Added support for -contact in ac and hbm calculations below.
4740 - Erik Marklund, May 29, 2006
4744 if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
4746 please_cite(stdout, "Spoel2006b");
4748 if (opt2bSet("-ac", NFILE, fnm))
4750 char *gemstring = NULL;
4754 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4755 gemBallistic, nBalExp);
4758 gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4761 gemstring = strdup(gemType[hb->per->gemtype]);
4762 do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
4763 bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
4764 gemstring, nThreads, NN, bBallistic, bGemFit);
4766 if (opt2bSet("-life", NFILE, fnm))
4768 do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
4770 if (opt2bSet("-hbm", NFILE, fnm))
4773 int id, ia, hh, x, y;
4775 if ((nframes > 0) && (hb->nrhb > 0))
4780 snew(mat.matrix, mat.nx);
4781 for (x = 0; (x < mat.nx); x++)
4783 snew(mat.matrix[x], mat.ny);
4786 for (id = 0; (id < hb->d.nrd); id++)
4788 for (ia = 0; (ia < hb->a.nra); ia++)
4790 for (hh = 0; (hh < hb->maxhydro); hh++)
4792 if (hb->hbmap[id][ia])
4794 if (ISHB(hb->hbmap[id][ia]->history[hh]))
4796 /* Changed '<' into '<=' in the for-statement below.
4797 * It fixed the previously undiscovered bug that caused
4798 * the last occurance of an hbond/contact to not be
4799 * set in mat.matrix. Have a look at any old -hbm-output
4800 * and you will notice that the last column is allways empty.
4801 * - Erik Marklund May 30, 2006
4803 for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
4805 int nn0 = hb->hbmap[id][ia]->n0;
4806 range_check(y, 0, mat.ny);
4807 mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
4815 mat.axis_x = hb->time;
4816 snew(mat.axis_y, mat.ny);
4817 for (j = 0; j < mat.ny; j++)
4821 sprintf(mat.title, bContact ? "Contact Existence Map" :
4822 "Hydrogen Bond Existence Map");
4823 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
4824 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
4825 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4826 mat.bDiscrete = TRUE;
4828 snew(mat.map, mat.nmap);
4829 for (i = 0; i < mat.nmap; i++)
4831 mat.map[i].code.c1 = hbmap[i];
4832 mat.map[i].desc = hbdesc[i];
4833 mat.map[i].rgb = hbrgb[i];
4835 fp = opt2FILE("-hbm", NFILE, fnm, "w");
4836 write_xpm_m(fp, mat);
4838 for (x = 0; x < mat.nx; x++)
4840 sfree(mat.matrix[x]);
4848 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
4855 fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4856 fprintf(stderr, "Freeing pHist for all donors...\n");
4857 for (i = 0; i < hb->d.nrd; i++)
4859 fprintf(stderr, "\r%i", i);
4860 if (hb->per->pHist[i] != NULL)
4862 for (j = 0; j < hb->a.nra; j++)
4864 clearPshift(&(hb->per->pHist[i][j]));
4866 sfree(hb->per->pHist[i]);
4869 sfree(hb->per->pHist);
4870 sfree(hb->per->p2i);
4872 fprintf(stderr, "...done.\n");
4875 #ifdef HAVE_NN_LOOPS
4888 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4890 fp = xvgropen(opt2fn("-dan", NFILE, fnm),
4891 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4892 nleg = (bTwo ? 2 : 1)*2;
4893 snew(legnames, nleg);
4895 for (j = 0; j < grNR; j++)
4897 if (USE_THIS_GROUP(j) )
4899 sprintf(buf, "Donors %s", grpnames[j]);
4900 legnames[i++] = strdup(buf);
4901 sprintf(buf, "Acceptors %s", grpnames[j]);
4902 legnames[i++] = strdup(buf);
4907 gmx_incons("number of legend entries");
4909 xvgr_legend(fp, nleg, (const char**)legnames, oenv);
4910 for (i = 0; i < nframes; i++)
4912 fprintf(fp, "%10g", hb->time[i]);
4913 for (j = 0; (j < grNR); j++)
4915 if (USE_THIS_GROUP(j) )
4917 fprintf(fp, " %6d", hb->danr[i][j]);