2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016, 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.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/commandline/viewit.h"
49 #include "gromacs/correlationfunctions/autocorr.h"
50 #include "gromacs/correlationfunctions/crosscorr.h"
51 #include "gromacs/correlationfunctions/expfit.h"
52 #include "gromacs/correlationfunctions/integrate.h"
53 #include "gromacs/fileio/matio.h"
54 #include "gromacs/fileio/tpxio.h"
55 #include "gromacs/fileio/trxio.h"
56 #include "gromacs/fileio/xvgr.h"
57 #include "gromacs/gmxana/gmx_ana.h"
58 #include "gromacs/gmxana/gstat.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/mdrunutility/mdmodules.h"
63 #include "gromacs/mdtypes/inputrec.h"
64 #include "gromacs/pbcutil/pbc.h"
65 #include "gromacs/topology/ifunc.h"
66 #include "gromacs/topology/index.h"
67 #include "gromacs/topology/topology.h"
68 #include "gromacs/utility/arraysize.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/exceptions.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/futil.h"
73 #include "gromacs/utility/gmxomp.h"
74 #include "gromacs/utility/pleasecite.h"
75 #include "gromacs/utility/programcontext.h"
76 #include "gromacs/utility/smalloc.h"
77 #include "gromacs/utility/snprintf.h"
80 typedef int t_hx[max_hx];
81 #define NRHXTYPES max_hx
82 const char *hxtypenames[NRHXTYPES] =
83 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
86 static const int NOTSET = -49297;
88 /* -----------------------------------------*/
94 hbNo, hbDist, hbHB, hbNR, hbR2
96 static const unsigned char c_acceptorMask = (1 << 0);
97 static const unsigned char c_donorMask = (1 << 1);
98 static const unsigned char c_inGroupMask = (1 << 2);
101 static const char *grpnames[grNR] = {"0", "1", "I" };
103 static gmx_bool bDebug = FALSE;
108 #define HB_YESINS HB_YES|HB_INS
112 #define ISHB(h) ((h) & 2)
113 #define ISDIST(h) ((h) & 1)
114 #define ISDIST2(h) ((h) & 4)
115 #define ISACC(h) ((h) & c_acceptorMask)
116 #define ISDON(h) ((h) & c_donorMask)
117 #define ISINGRP(h) ((h) & c_inGroupMask)
130 typedef int t_icell[grNR];
131 typedef int h_id[MAXHYDRO];
134 int history[MAXHYDRO];
135 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
136 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
138 /* Bitmask array which tells whether a hbond is present
139 * at a given time. Either of these may be NULL
141 int n0; /* First frame a HB was found */
142 int nframes, maxframes; /* Amount of frames in this hbond */
145 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
146 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
147 * acceptor distance is less than the user-specified distance (typically
154 int *acc; /* Atom numbers of the acceptors */
155 int *grp; /* Group index */
156 int *aptr; /* Map atom number to acceptor index */
161 int *don; /* Atom numbers of the donors */
162 int *grp; /* Group index */
163 int *dptr; /* Map atom number to donor index */
164 int *nhydro; /* Number of hydrogens for each donor */
165 h_id *hydro; /* The atom numbers of the hydrogens */
166 h_id *nhbonds; /* The number of HBs per H at current */
170 gmx_bool bHBmap, bDAnr;
172 /* The following arrays are nframes long */
173 int nframes, max_frames, maxhydro;
179 /* These structures are initialized from the topology at start up */
182 /* This holds a matrix with all possible hydrogen bonds */
187 /* Changed argument 'bMerge' into 'oneHB' below,
188 * since -contact should cause maxhydro to be 1,
190 * - Erik Marklund May 29, 2006
193 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB)
198 hb->wordlen = 8*sizeof(unsigned int);
207 hb->maxhydro = MAXHYDRO;
212 static void mk_hbmap(t_hbdata *hb)
216 snew(hb->hbmap, hb->d.nrd);
217 for (i = 0; (i < hb->d.nrd); i++)
219 snew(hb->hbmap[i], hb->a.nra);
220 if (hb->hbmap[i] == NULL)
222 gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
224 for (j = 0; (j > hb->a.nra); j++)
226 hb->hbmap[i][j] = NULL;
231 static void add_frames(t_hbdata *hb, int nframes)
233 if (nframes >= hb->max_frames)
235 hb->max_frames += 4096;
236 srenew(hb->time, hb->max_frames);
237 srenew(hb->nhb, hb->max_frames);
238 srenew(hb->ndist, hb->max_frames);
239 srenew(hb->n_bound, hb->max_frames);
240 srenew(hb->nhx, hb->max_frames);
243 srenew(hb->danr, hb->max_frames);
246 hb->nframes = nframes;
249 #define OFFSET(frame) (frame / 32)
250 #define MASK(frame) (1 << (frame % 32))
252 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
256 hbexist[OFFSET(frame)] |= MASK(frame);
260 hbexist[OFFSET(frame)] &= ~MASK(frame);
264 static gmx_bool is_hb(unsigned int hbexist[], int frame)
266 return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
269 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
271 unsigned int *ghptr = NULL;
275 ghptr = hb->hbmap[id][ia]->h[ih];
277 else if (ihb == hbDist)
279 ghptr = hb->hbmap[id][ia]->g[ih];
283 gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
286 _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
289 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb)
292 t_hbond *hb = hbd->hbmap[id][ia];
293 int maxhydro = std::min(hbd->maxhydro, hbd->d.nhydro[id]);
294 int wlen = hbd->wordlen;
300 hb->maxframes = delta;
301 for (i = 0; (i < maxhydro); i++)
303 snew(hb->h[i], hb->maxframes/wlen);
304 snew(hb->g[i], hb->maxframes/wlen);
309 hb->nframes = frame-hb->n0;
310 /* We need a while loop here because hbonds may be returning
313 while (hb->nframes >= hb->maxframes)
315 n = hb->maxframes + delta;
316 for (i = 0; (i < maxhydro); i++)
318 srenew(hb->h[i], n/wlen);
319 srenew(hb->g[i], n/wlen);
320 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
332 set_hb(hbd, id, h, ia, frame, ihb);
337 static void inc_nhbonds(t_donors *ddd, int d, int h)
340 int dptr = ddd->dptr[d];
342 for (j = 0; (j < ddd->nhydro[dptr]); j++)
344 if (ddd->hydro[dptr][j] == h)
346 ddd->nhbonds[dptr][j]++;
350 if (j == ddd->nhydro[dptr])
352 gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
356 static int _acceptor_index(t_acceptors *a, int grp, int i,
357 const char *file, int line)
361 if (a->grp[ai] != grp)
365 fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
366 ai, a->grp[ai], grp, file, line);
375 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
377 static int _donor_index(t_donors *d, int grp, int i, const char *file, int line)
386 if (d->grp[di] != grp)
390 fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
391 di, d->grp[di], grp, file, line);
400 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
402 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
404 /* g_hbond doesn't allow overlapping groups */
410 donor_index(&hb->d, grpd, a) != NOTSET
411 && acceptor_index(&hb->a, grpa, d) != NOTSET;
415 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
416 int frame, gmx_bool bMerge, int ihb, gmx_bool bContact)
419 gmx_bool daSwap = FALSE;
421 if ((id = hb->d.dptr[d]) == NOTSET)
423 gmx_fatal(FARGS, "No donor atom %d", d+1);
425 else if (grpd != hb->d.grp[id])
427 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
428 grpd, hb->d.grp[id], d+1);
430 if ((ia = hb->a.aptr[a]) == NOTSET)
432 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
434 else if (grpa != hb->a.grp[ia])
436 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
437 grpa, hb->a.grp[ia], a+1);
443 if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
444 /* Then swap identity so that the id of d is lower then that of a.
446 * This should really be redundant by now, as is_hbond() now ought to return
447 * hbNo in the cases where this conditional is TRUE. */
454 /* Now repeat donor/acc check. */
455 if ((id = hb->d.dptr[d]) == NOTSET)
457 gmx_fatal(FARGS, "No donor atom %d", d+1);
459 else if (grpd != hb->d.grp[id])
461 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
462 grpd, hb->d.grp[id], d+1);
464 if ((ia = hb->a.aptr[a]) == NOTSET)
466 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
468 else if (grpa != hb->a.grp[ia])
470 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
471 grpa, hb->a.grp[ia], a+1);
478 /* Loop over hydrogens to find which hydrogen is in this particular HB */
479 if ((ihb == hbHB) && !bMerge && !bContact)
481 for (k = 0; (k < hb->d.nhydro[id]); k++)
483 if (hb->d.hydro[id][k] == h)
488 if (k == hb->d.nhydro[id])
490 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
506 if (hb->hbmap[id][ia] == NULL)
508 snew(hb->hbmap[id][ia], 1);
509 snew(hb->hbmap[id][ia]->h, hb->maxhydro);
510 snew(hb->hbmap[id][ia]->g, hb->maxhydro);
512 add_ff(hb, id, k, ia, frame, ihb);
514 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
518 /* Strange construction with frame >=0 is a relic from old code
519 * for selected hbond analysis. It may be necessary again if that
520 * is made to work again.
524 hh = hb->hbmap[id][ia]->history[k];
530 hb->hbmap[id][ia]->history[k] = hh | 2;
541 hb->hbmap[id][ia]->history[k] = hh | 1;
565 if (bMerge && daSwap)
567 h = hb->d.hydro[id][0];
569 /* Increment number if HBonds per H */
570 if (ihb == hbHB && !bContact)
572 inc_nhbonds(&(hb->d), d, h);
576 static char *mkatomname(const t_atoms *atoms, int i)
581 rnr = atoms->atom[i].resind;
582 sprintf(buf, "%4s%d%-4s",
583 *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
588 static void gen_datable(int *index, int isize, unsigned char *datable, int natoms)
590 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
593 for (i = 0; i < isize; i++)
595 if (index[i] >= natoms)
597 gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
599 datable[index[i]] |= c_inGroupMask;
603 static void clear_datable_grp(unsigned char *datable, int size)
605 /* Clears group information from the table */
609 for (i = 0; i < size; i++)
611 datable[i] &= ~c_inGroupMask;
616 static void add_acc(t_acceptors *a, int ia, int grp)
618 if (a->nra >= a->max_nra)
621 srenew(a->acc, a->max_nra);
622 srenew(a->grp, a->max_nra);
624 a->grp[a->nra] = grp;
625 a->acc[a->nra++] = ia;
628 static void search_acceptors(const t_topology *top, int isize,
629 const int *index, t_acceptors *a, int grp,
631 gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
637 for (i = 0; (i < isize); i++)
641 (((*top->atoms.atomname[n])[0] == 'O') ||
642 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
645 datable[n] |= c_acceptorMask;
650 snew(a->aptr, top->atoms.nr);
651 for (i = 0; (i < top->atoms.nr); i++)
655 for (i = 0; (i < a->nra); i++)
657 a->aptr[a->acc[i]] = i;
661 static void add_h2d(int id, int ih, t_donors *ddd)
665 for (i = 0; (i < ddd->nhydro[id]); i++)
667 if (ddd->hydro[id][i] == ih)
669 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
674 if (i == ddd->nhydro[id])
676 if (ddd->nhydro[id] >= MAXHYDRO)
678 gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
679 ddd->don[id], MAXHYDRO);
681 ddd->hydro[id][i] = ih;
686 static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable)
690 if (!datable || ISDON(datable[id]))
692 if (ddd->dptr[id] == NOTSET) /* New donor */
704 if (ddd->nrd >= ddd->max_nrd)
707 srenew(ddd->don, ddd->max_nrd);
708 srenew(ddd->nhydro, ddd->max_nrd);
709 srenew(ddd->hydro, ddd->max_nrd);
710 srenew(ddd->nhbonds, ddd->max_nrd);
711 srenew(ddd->grp, ddd->max_nrd);
713 ddd->don[ddd->nrd] = id;
714 ddd->nhydro[ddd->nrd] = 0;
715 ddd->grp[ddd->nrd] = grp;
727 printf("Warning: Atom %d is not in the d/a-table!\n", id);
731 static void search_donors(const t_topology *top, int isize, const int *index,
732 t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
733 unsigned char *datable)
736 t_functype func_type;
741 snew(ddd->dptr, top->atoms.nr);
742 for (i = 0; (i < top->atoms.nr); i++)
744 ddd->dptr[i] = NOTSET;
752 for (i = 0; (i < isize); i++)
754 datable[index[i]] |= c_donorMask;
755 add_dh(ddd, index[i], -1, grp, datable);
761 for (func_type = 0; (func_type < F_NRE); func_type++)
763 const t_ilist *interaction = &(top->idef.il[func_type]);
764 if (func_type == F_POSRES || func_type == F_FBPOSRES)
766 /* The ilist looks strange for posre. Bug in grompp?
767 * We don't need posre interactions for hbonds anyway.*/
770 for (i = 0; i < interaction->nr;
771 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
774 if (func_type != top->idef.functype[interaction->iatoms[i]])
776 fprintf(stderr, "Error in func_type %s",
777 interaction_function[func_type].longname);
781 /* check out this functype */
782 if (func_type == F_SETTLE)
784 nr1 = interaction->iatoms[i+1];
785 nr2 = interaction->iatoms[i+2];
786 nr3 = interaction->iatoms[i+3];
788 if (ISINGRP(datable[nr1]))
790 if (ISINGRP(datable[nr2]))
792 datable[nr1] |= c_donorMask;
793 add_dh(ddd, nr1, nr1+1, grp, datable);
795 if (ISINGRP(datable[nr3]))
797 datable[nr1] |= c_donorMask;
798 add_dh(ddd, nr1, nr1+2, grp, datable);
802 else if (IS_CHEMBOND(func_type))
804 for (j = 0; j < 2; j++)
806 nr1 = interaction->iatoms[i+1+j];
807 nr2 = interaction->iatoms[i+2-j];
808 if ((*top->atoms.atomname[nr1][0] == 'H') &&
809 ((*top->atoms.atomname[nr2][0] == 'O') ||
810 (*top->atoms.atomname[nr2][0] == 'N')) &&
811 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
813 datable[nr2] |= c_donorMask;
814 add_dh(ddd, nr2, nr1, grp, datable);
821 for (func_type = 0; func_type < F_NRE; func_type++)
823 interaction = &top->idef.il[func_type];
824 for (i = 0; i < interaction->nr;
825 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
828 if (func_type != top->idef.functype[interaction->iatoms[i]])
830 gmx_incons("function type in search_donors");
833 if (interaction_function[func_type].flags & IF_VSITE)
835 nr1 = interaction->iatoms[i+1];
836 if (*top->atoms.atomname[nr1][0] == 'H')
840 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
851 if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
852 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
853 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
855 datable[nr2] |= c_donorMask;
856 add_dh(ddd, nr2, nr1, grp, datable);
866 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
873 for (i = 0; i < DIM; i++)
875 ngrid[i] = static_cast<int>(box[i][i]/(1.2*rcut));
879 if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
881 for (i = 0; i < DIM; i++)
888 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
889 ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
891 snew(grid, ngrid[ZZ]);
892 for (z = 0; z < ngrid[ZZ]; z++)
894 snew((grid)[z], ngrid[YY]);
895 for (y = 0; y < ngrid[YY]; y++)
897 snew((grid)[z][y], ngrid[XX]);
903 static void reset_nhbonds(t_donors *ddd)
907 for (i = 0; (i < ddd->nrd); i++)
909 for (j = 0; (j < MAXHH); j++)
911 ddd->nhbonds[i][j] = 0;
916 static void pbc_correct_gem(rvec dx, matrix box, rvec hbox);
917 static void pbc_in_gridbox(rvec dx, matrix box);
919 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
920 gmx_bool bBox, matrix box, rvec hbox,
921 real rcut, real rshell,
922 ivec ngrid, t_gridcell ***grid)
924 int i, m, gr, xi, yi, zi, nr;
927 rvec invdelta, dshell;
929 gmx_bool bDoRshell, bInShell, bAcc;
934 bDoRshell = (rshell > 0);
935 rshell2 = gmx::square(rshell);
938 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
940 for (m = 0; m < DIM; m++)
942 hbox[m] = box[m][m]*0.5;
945 invdelta[m] = ngrid[m]/box[m][m];
946 if (1/invdelta[m] < rcut)
948 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
949 "%s can not handle this situation, sorry.\n",
950 gmx::getProgramContext().displayName());
962 /* resetting atom counts */
963 for (gr = 0; (gr < grNR); gr++)
965 for (zi = 0; zi < ngrid[ZZ]; zi++)
967 for (yi = 0; yi < ngrid[YY]; yi++)
969 for (xi = 0; xi < ngrid[XX]; xi++)
971 grid[zi][yi][xi].d[gr].nr = 0;
972 grid[zi][yi][xi].a[gr].nr = 0;
978 /* put atoms in grid cells */
979 for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
992 for (i = 0; (i < nr); i++)
994 /* check if we are inside the shell */
995 /* if bDoRshell=FALSE then bInShell=TRUE always */
1000 rvec_sub(x[ad[i]], xshell, dshell);
1003 gmx_bool bDone = FALSE;
1007 for (m = DIM-1; m >= 0 && bInShell; m--)
1009 if (dshell[m] < -hbox[m])
1012 rvec_inc(dshell, box[m]);
1014 if (dshell[m] >= hbox[m])
1017 dshell[m] -= 2*hbox[m];
1021 for (m = DIM-1; m >= 0 && bInShell; m--)
1023 /* if we're outside the cube, we're outside the sphere also! */
1024 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1030 /* if we're inside the cube, check if we're inside the sphere */
1033 bInShell = norm2(dshell) < rshell2;
1041 pbc_in_gridbox(x[ad[i]], box);
1043 for (m = DIM-1; m >= 0; m--)
1044 { /* determine grid index of atom */
1045 grididx[m] = static_cast<int>(x[ad[i]][m]*invdelta[m]);
1046 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1053 range_check(gx, 0, ngrid[XX]);
1054 range_check(gy, 0, ngrid[YY]);
1055 range_check(gz, 0, ngrid[ZZ]);
1059 /* add atom to grid cell */
1062 newgrid = &(grid[gz][gy][gx].a[gr]);
1066 newgrid = &(grid[gz][gy][gx].d[gr]);
1068 if (newgrid->nr >= newgrid->maxnr)
1070 newgrid->maxnr += 10;
1071 DBB(newgrid->maxnr);
1072 srenew(newgrid->atoms, newgrid->maxnr);
1075 newgrid->atoms[newgrid->nr] = ad[i];
1083 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1087 for (gr = 0; (gr < grNR); gr++)
1090 for (zi = 0; zi < ngrid[ZZ]; zi++)
1092 for (yi = 0; yi < ngrid[YY]; yi++)
1094 for (xi = 0; xi < ngrid[XX]; xi++)
1096 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1104 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1105 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1106 * With a triclinic box all loops are 3 long, except when a cell is
1107 * located next to one of the box edges which is not parallel to the
1108 * x/y-plane, in that case all cells in a line or layer are searched.
1109 * This could be implemented slightly more efficient, but the code
1110 * would get much more complicated.
1112 static gmx_inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1114 return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
1116 static gmx_inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1118 return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1120 static gmx_inline int grid_mod(int j, int n)
1125 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1127 int gr, x, y, z, sum[grNR];
1129 fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1130 for (gr = 0; gr < grNR; gr++)
1133 fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1134 for (z = 0; z < ngrid[ZZ]; z += 2)
1136 fprintf(fp, "Z=%d,%d\n", z, z+1);
1137 for (y = 0; y < ngrid[YY]; y++)
1139 for (x = 0; x < ngrid[XX]; x++)
1141 fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1142 sum[gr] += grid[z][y][x].d[gr].nr;
1143 fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1144 sum[gr] += grid[z][y][x].a[gr].nr;
1148 if ( (z+1) < ngrid[ZZ])
1150 for (x = 0; x < ngrid[XX]; x++)
1152 fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1153 sum[gr] += grid[z+1][y][x].d[gr].nr;
1154 fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1155 sum[gr] += grid[z+1][y][x].a[gr].nr;
1162 fprintf(fp, "TOTALS:");
1163 for (gr = 0; gr < grNR; gr++)
1165 fprintf(fp, " %d=%d", gr, sum[gr]);
1170 /* New GMX record! 5 * in a row. Congratulations!
1171 * Sorry, only four left.
1173 static void free_grid(ivec ngrid, t_gridcell ****grid)
1176 t_gridcell ***g = *grid;
1178 for (z = 0; z < ngrid[ZZ]; z++)
1180 for (y = 0; y < ngrid[YY]; y++)
1190 static void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1193 gmx_bool bDone = FALSE;
1197 for (m = DIM-1; m >= 0; m--)
1199 if (dx[m] < -hbox[m])
1202 rvec_inc(dx, box[m]);
1204 if (dx[m] >= hbox[m])
1207 rvec_dec(dx, box[m]);
1213 static void pbc_in_gridbox(rvec dx, matrix box)
1216 gmx_bool bDone = FALSE;
1220 for (m = DIM-1; m >= 0; m--)
1225 rvec_inc(dx, box[m]);
1227 if (dx[m] >= box[m][m])
1230 rvec_dec(dx, box[m]);
1236 /* Added argument r2cut, changed contact and implemented
1237 * use of second cut-off.
1238 * - Erik Marklund, June 29, 2006
1240 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1241 real rcut, real r2cut, real ccut,
1242 rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1243 real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1244 gmx_bool bContact, gmx_bool bMerge)
1247 rvec r_da, r_ha, r_dh;
1248 real rc2, r2c2, rda2, rha2, ca;
1249 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1250 gmx_bool daSwap = FALSE;
1257 if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1258 (acceptor_index(&hb->a, grpa, a) == NOTSET))
1266 rvec_sub(x[d], x[a], r_da);
1267 /* Insert projection code here */
1269 if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1271 /* Then this hbond/contact will be found again, or it has already been found. */
1276 if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1277 { /* return hbNo; */
1278 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1280 pbc_correct_gem(r_da, box, hbox);
1282 rda2 = iprod(r_da, r_da);
1286 if (daSwap && grpa == grpd)
1294 else if (rda2 < r2c2)
1305 if (bDA && (rda2 > rc2))
1310 for (h = 0; (h < hb->d.nhydro[id]); h++)
1312 hh = hb->d.hydro[id][h];
1316 rvec_sub(x[hh], x[a], r_ha);
1319 pbc_correct_gem(r_ha, box, hbox);
1321 rha2 = iprod(r_ha, r_ha);
1324 if (bDA || (rha2 <= rc2))
1326 rvec_sub(x[d], x[hh], r_dh);
1329 pbc_correct_gem(r_dh, box, hbox);
1336 ca = cos_angle(r_dh, r_da);
1337 /* if angle is smaller, cos is larger */
1341 *d_ha = std::sqrt(bDA ? rda2 : rha2);
1342 *ang = std::acos(ca);
1347 if (bDA || HAinrange)
1357 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1358 * Will do some more testing before removing the function entirely.
1359 * - Erik Marklund, MAY 10 2010 */
1360 static void do_merge(t_hbdata *hb, int ntmp,
1361 unsigned int htmp[], unsigned int gtmp[],
1362 t_hbond *hb0, t_hbond *hb1)
1364 /* Here we need to make sure we're treating periodicity in
1365 * the right way for the geminate recombination kinetics. */
1367 int m, mm, n00, n01, nn0, nnframes;
1369 /* Decide where to start from when merging */
1372 nn0 = std::min(n00, n01);
1373 nnframes = std::max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1374 /* Initiate tmp arrays */
1375 for (m = 0; (m < ntmp); m++)
1380 /* Fill tmp arrays with values due to first HB */
1381 /* Once again '<' had to be replaced with '<='
1382 to catch the last frame in which the hbond
1384 - Erik Marklund, June 1, 2006 */
1385 for (m = 0; (m <= hb0->nframes); m++)
1388 htmp[mm] = is_hb(hb0->h[0], m);
1390 for (m = 0; (m <= hb0->nframes); m++)
1393 gtmp[mm] = is_hb(hb0->g[0], m);
1396 for (m = 0; (m <= hb1->nframes); m++)
1399 htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1400 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1402 /* Reallocate target array */
1403 if (nnframes > hb0->maxframes)
1405 srenew(hb0->h[0], 4+nnframes/hb->wordlen);
1406 srenew(hb0->g[0], 4+nnframes/hb->wordlen);
1409 /* Copy temp array to target array */
1410 for (m = 0; (m <= nnframes); m++)
1412 _set_hb(hb0->h[0], m, htmp[m]);
1413 _set_hb(hb0->g[0], m, gtmp[m]);
1416 /* Set scalar variables */
1418 hb0->maxframes = nnframes;
1421 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
1423 int i, inrnew, indnew, j, ii, jj, id, ia, ntmp;
1424 unsigned int *htmp, *gtmp;
1428 indnew = hb->nrdist;
1430 /* Check whether donors are also acceptors */
1431 printf("Merging hbonds with Acceptor and Donor swapped\n");
1433 ntmp = 2*hb->max_frames;
1436 for (i = 0; (i < hb->d.nrd); i++)
1438 fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
1441 ii = hb->a.aptr[id];
1442 for (j = 0; (j < hb->a.nra); j++)
1445 jj = hb->d.dptr[ia];
1446 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
1447 (!bTwo || (hb->d.grp[i] != hb->a.grp[j])))
1449 hb0 = hb->hbmap[i][j];
1450 hb1 = hb->hbmap[jj][ii];
1451 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
1453 do_merge(hb, ntmp, htmp, gtmp, hb0, hb1);
1454 if (ISHB(hb1->history[0]))
1458 else if (ISDIST(hb1->history[0]))
1465 gmx_incons("No contact history");
1469 gmx_incons("Neither hydrogen bond nor distance");
1475 hb1->history[0] = hbNo;
1480 fprintf(stderr, "\n");
1481 printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
1482 printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
1484 hb->nrdist = indnew;
1489 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
1491 int i, j, k, n_bound[MAXHH], nbtot;
1493 /* Set array to 0 */
1494 for (k = 0; (k < MAXHH); k++)
1498 /* Loop over possible donors */
1499 for (i = 0; (i < hb->d.nrd); i++)
1501 for (j = 0; (j < hb->d.nhydro[i]); j++)
1503 n_bound[hb->d.nhbonds[i][j]]++;
1506 fprintf(fp, "%12.5e", t);
1508 for (k = 0; (k < MAXHH); k++)
1510 fprintf(fp, " %8d", n_bound[k]);
1511 nbtot += n_bound[k]*k;
1513 fprintf(fp, " %8d\n", nbtot);
1516 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
1517 const gmx_output_env_t *oenv)
1520 const char *leg[] = { "p(t)", "t p(t)" };
1522 int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
1523 int nframes = hb->nframes;
1526 double sum, integral;
1529 snew(h, hb->maxhydro);
1530 snew(histo, nframes+1);
1531 /* Total number of hbonds analyzed here */
1532 for (i = 0; (i < hb->d.nrd); i++)
1534 for (k = 0; (k < hb->a.nra); k++)
1536 hbh = hb->hbmap[i][k];
1554 for (m = 0; (m < hb->maxhydro); m++)
1558 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
1562 for (nh = 0; (nh < nhydro); nh++)
1567 for (j = 0; (j <= hbh->nframes); j++)
1569 ihb = is_hb(h[nh], j);
1570 if (debug && (ndump < 10))
1572 fprintf(debug, "%5d %5d\n", j, ihb);
1592 fprintf(stderr, "\n");
1595 fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
1599 fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
1603 xvgr_legend(fp, asize(leg), leg, oenv);
1605 while ((j0 > 0) && (histo[j0] == 0))
1610 for (i = 0; (i <= j0); i++)
1614 dt = hb->time[1]-hb->time[0];
1617 for (i = 1; (i <= j0); i++)
1619 t = hb->time[i] - hb->time[0] - 0.5*dt;
1620 x1 = t*histo[i]/sum;
1621 fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1);
1626 printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
1627 printf("Note that the lifetime obtained in this manner is close to useless\n");
1628 printf("Use the -ac option instead and check the Forward lifetime\n");
1629 please_cite(stdout, "Spoel2006b");
1634 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
1637 int i, j, k, m, nd, ihb, idist;
1638 int nframes = hb->nframes;
1646 fp = gmx_ffopen("debug-ac.xvg", "w");
1647 for (j = 0; (j < nframes); j++)
1649 fprintf(fp, "%10.3f", hb->time[j]);
1650 for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
1652 for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
1656 hbh = hb->hbmap[i][k];
1661 ihb = is_hb(hbh->h[0], j);
1662 idist = is_hb(hbh->g[0], j);
1668 for (m = 0; (m < hb->maxhydro) && !ihb; m++)
1670 ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m], j));
1671 idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
1673 /* This is not correct! */
1674 /* What isn't correct? -Erik M */
1679 fprintf(fp, " %1d-%1d", ihb, idist);
1689 static real calc_dg(real tau, real temp)
1700 return kbt*std::log(kbt*tau/PLANCK);
1705 int n0, n1, nparams, ndelta;
1707 real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
1710 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
1711 real kt[], real sigma_ct[], real sigma_nt[],
1712 real sigma_kt[], real *k, real *kp,
1713 real *sigma_k, real *sigma_kp,
1719 real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
1724 for (i = 0; (i < n); i++)
1726 if (t[i] >= fit_start)
1739 tl.sigma_ct = sigma_ct;
1740 tl.sigma_nt = sigma_nt;
1741 tl.sigma_kt = sigma_kt;
1745 chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
1747 *kp = tl.kkk[1] = *kp;
1749 for (j = 0; (j < NK); j++)
1753 kk2 += gmx::square(tl.kkk[0]);
1754 kp2 += gmx::square(tl.kkk[1]);
1757 *sigma_k = std::sqrt(kk2/NK - gmx::square(kkk/NK));
1758 *sigma_kp = std::sqrt(kp2/NK - gmx::square(kkp/NK));
1763 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
1764 real sigma_ct[], real sigma_nt[], real sigma_kt[],
1765 real fit_start, real temp)
1768 real k = 1, kp = 1, kow = 1;
1769 real Q = 0, chi2, tau_hb, dtau, tau_rlx, e_1, sigma_k, sigma_kp, ddg;
1770 double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
1771 gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
1773 for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
1779 for (i = i0; (i < n); i++)
1781 sc2 += gmx::square(ct[i]);
1782 sn2 += gmx::square(nt[i]);
1783 sk2 += gmx::square(kt[i]);
1788 printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
1789 tmp = (sn2*sc2-gmx::square(scn));
1790 if ((tmp > 0) && (sn2 > 0))
1792 k = (sn2*sck-scn*snk)/tmp;
1793 kp = (k*scn-snk)/sn2;
1797 for (i = i0; (i < n); i++)
1799 chi2 += gmx::square(k*ct[i]-kp*nt[i]-kt[i]);
1801 compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
1803 &sigma_k, &sigma_kp, fit_start);
1804 Q = 0; /* quality_of_fit(chi2, 2);*/
1805 ddg = BOLTZ*temp*sigma_k/k;
1806 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
1808 printf("The Rate and Delta G are followed by an error estimate\n");
1809 printf("----------------------------------------------------------\n"
1810 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
1811 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1812 k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
1813 ddg = BOLTZ*temp*sigma_kp/kp;
1814 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1815 kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
1820 for (i = i0; (i < n); i++)
1822 chi2 += gmx::square(k*ct[i]-kp*nt[i]-kt[i]);
1824 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
1826 printf("--------------------------------------------------\n"
1827 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
1828 printf("Forward %10.3f %8.3f %10.3f %10g\n",
1829 k, 1/k, calc_dg(1/k, temp), chi2);
1830 printf("Backward %10.3f %8.3f %10.3f\n",
1831 kp, 1/kp, calc_dg(1/kp, temp));
1837 printf("One-way %10.3f %s%8.3f %10.3f\n",
1838 kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp));
1842 printf(" - Numerical problems computing HB thermodynamics:\n"
1843 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
1844 sc2, sn2, sk2, sck, snk, scn);
1846 /* Determine integral of the correlation function */
1847 tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
1848 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
1849 bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
1850 e_1 = std::exp(-1.0);
1851 for (i = 0; (i < n-2); i++)
1853 if ((ct[i] > e_1) && (ct[i+1] <= e_1))
1860 /* Determine tau_relax from linear interpolation */
1861 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
1862 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx,
1863 tau_rlx, bError ? " " : "",
1864 calc_dg(tau_rlx, temp));
1869 printf("Correlation functions too short to compute thermodynamics\n");
1873 void compute_derivative(int nn, real x[], real y[], real dydx[])
1877 /* Compute k(t) = dc(t)/dt */
1878 for (j = 1; (j < nn-1); j++)
1880 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
1882 /* Extrapolate endpoints */
1883 dydx[0] = 2*dydx[1] - dydx[2];
1884 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
1887 static void normalizeACF(real *ct, real *gt, int nhb, int len)
1889 real ct_fac, gt_fac = 0;
1892 /* Xu and Berne use the same normalization constant */
1900 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
1901 for (i = 0; i < len; i++)
1911 static void do_hbac(const char *fn, t_hbdata *hb,
1912 int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
1913 real temp, gmx_bool R2, const gmx_output_env_t *oenv,
1917 int i, j, k, m, ihb, idist, n2, nn;
1919 const char *legLuzar[] = {
1920 "Ac\\sfin sys\\v{}\\z{}(t)",
1922 "Cc\\scontact,hb\\v{}\\z{}(t)",
1923 "-dAc\\sfs\\v{}\\z{}/dt"
1925 gmx_bool bNorm = FALSE;
1927 real *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
1928 real *ct, tail, tail2, dtail, *cct;
1929 const real tol = 1e-3;
1930 int nframes = hb->nframes;
1931 unsigned int **h = NULL, **g = NULL;
1932 int nh, nhbonds, nhydro;
1935 int *dondata = NULL;
1938 AC_NONE, AC_NN, AC_GEM, AC_LUZAR
1941 const bool bOMP = GMX_OPENMP;
1943 printf("Doing autocorrelation ");
1946 printf("according to the theory of Luzar and Chandler.\n");
1949 /* build hbexist matrix in reals for autocorr */
1950 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
1952 while (n2 < nframes)
1959 if (acType != AC_NN || bOMP)
1961 snew(h, hb->maxhydro);
1962 snew(g, hb->maxhydro);
1965 /* Dump hbonds for debugging */
1966 dump_ac(hb, bMerge || bContact, nDump);
1968 /* Total number of hbonds analyzed here */
1971 if (acType != AC_LUZAR && bOMP)
1973 nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
1975 gmx_omp_set_num_threads(nThreads);
1976 snew(dondata, nThreads);
1977 for (i = 0; i < nThreads; i++)
1981 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
1982 "Expect close to linear scaling over this donor-loop.\n", nThreads);
1984 fprintf(stderr, "Donors: [thread no]\n");
1987 for (i = 0; i < nThreads; i++)
1989 snprintf(tmpstr, 7, "[%i]", i);
1990 fprintf(stderr, "%-7s", tmpstr);
1993 fprintf(stderr, "\n");
2008 for (i = 0; (i < hb->d.nrd); i++)
2010 for (k = 0; (k < hb->a.nra); k++)
2013 hbh = hb->hbmap[i][k];
2017 if (bMerge || bContact)
2019 if (ISHB(hbh->history[0]))
2028 for (m = 0; (m < hb->maxhydro); m++)
2030 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
2032 g[nhydro] = hbh->g[m];
2033 h[nhydro] = hbh->h[m];
2039 int nf = hbh->nframes;
2040 for (nh = 0; (nh < nhydro); nh++)
2042 int nrint = bContact ? hb->nrdist : hb->nrhb;
2043 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
2045 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
2049 for (j = 0; (j < nframes); j++)
2053 ihb = is_hb(h[nh], j);
2054 idist = is_hb(g[nh], j);
2061 /* For contacts: if a second cut-off is provided, use it,
2062 * otherwise use g(t) = 1-h(t) */
2063 if (!R2 && bContact)
2069 gt[j] = idist*(1-ihb);
2075 /* The autocorrelation function is normalized after summation only */
2076 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
2077 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
2079 /* Cross correlation analysis for thermodynamics */
2080 for (j = nframes; (j < n2); j++)
2086 cross_corr(n2, ht, gt, dght);
2088 for (j = 0; (j < nn); j++)
2097 fprintf(stderr, "\n");
2100 normalizeACF(ct, ght, static_cast<int>(nhb), nn);
2102 /* Determine tail value for statistics */
2105 for (j = nn/2; (j < nn); j++)
2108 tail2 += ct[j]*ct[j];
2110 tail /= (nn - nn/2);
2111 tail2 /= (nn - nn/2);
2112 dtail = std::sqrt(tail2-tail*tail);
2114 /* Check whether the ACF is long enough */
2117 printf("\nWARNING: Correlation function is probably not long enough\n"
2118 "because the standard deviation in the tail of C(t) > %g\n"
2119 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2122 for (j = 0; (j < nn); j++)
2125 ct[j] = (cct[j]-tail)/(1-tail);
2127 /* Compute negative derivative k(t) = -dc(t)/dt */
2128 compute_derivative(nn, hb->time, ct, kt);
2129 for (j = 0; (j < nn); j++)
2137 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
2141 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
2143 xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
2146 for (j = 0; (j < nn); j++)
2148 fprintf(fp, "%10g %10g %10g %10g %10g\n",
2149 hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
2153 analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
2156 do_view(oenv, fn, NULL);
2167 static void init_hbframe(t_hbdata *hb, int nframes, real t)
2171 hb->time[nframes] = t;
2172 hb->nhb[nframes] = 0;
2173 hb->ndist[nframes] = 0;
2174 for (i = 0; (i < max_hx); i++)
2176 hb->nhx[nframes][i] = 0;
2180 static FILE *open_donor_properties_file(const char *fn,
2182 const gmx_output_env_t *oenv)
2185 const char *leg[] = { "Nbound", "Nfree" };
2192 fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
2193 xvgr_legend(fp, asize(leg), leg, oenv);
2198 static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
2200 int i, j, k, nbound, nb, nhtot;
2208 for (i = 0; (i < hb->d.nrd); i++)
2210 for (k = 0; (k < hb->d.nhydro[i]); k++)
2214 for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
2216 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
2217 is_hb(hb->hbmap[i][j]->h[k], nframes))
2225 fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound);
2228 static void dump_hbmap(t_hbdata *hb,
2229 int nfile, t_filenm fnm[], gmx_bool bTwo,
2230 gmx_bool bContact, int isize[], int *index[], char *grpnames[],
2231 const t_atoms *atoms)
2234 int ddd, hhh, aaa, i, j, k, m, grp;
2235 char ds[32], hs[32], as[32];
2238 fp = opt2FILE("-hbn", nfile, fnm, "w");
2239 if (opt2bSet("-g", nfile, fnm))
2241 fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
2242 fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
2248 for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
2250 fprintf(fp, "[ %s ]", grpnames[grp]);
2251 for (i = 0; i < isize[grp]; i++)
2253 fprintf(fp, (i%15) ? " " : "\n");
2254 fprintf(fp, " %4d", index[grp][i]+1);
2260 fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
2261 for (i = 0; (i < hb->d.nrd); i++)
2263 if (hb->d.grp[i] == grp)
2265 for (j = 0; (j < hb->d.nhydro[i]); j++)
2267 fprintf(fp, " %4d %4d", hb->d.don[i]+1,
2268 hb->d.hydro[i][j]+1);
2274 fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
2275 for (i = 0; (i < hb->a.nra); i++)
2277 if (hb->a.grp[i] == grp)
2279 fprintf(fp, (i%15 && !first) ? " " : "\n");
2280 fprintf(fp, " %4d", hb->a.acc[i]+1);
2289 fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
2290 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
2294 fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
2297 for (i = 0; (i < hb->d.nrd); i++)
2300 for (k = 0; (k < hb->a.nra); k++)
2303 for (m = 0; (m < hb->d.nhydro[i]); m++)
2305 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
2307 sprintf(ds, "%s", mkatomname(atoms, ddd));
2308 sprintf(as, "%s", mkatomname(atoms, aaa));
2311 fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
2314 fprintf(fplog, "%12s %12s\n", ds, as);
2319 hhh = hb->d.hydro[i][m];
2320 sprintf(hs, "%s", mkatomname(atoms, hhh));
2321 fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
2324 fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
2338 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
2339 * It mimics add_frames() and init_frame() to some extent. */
2340 static void sync_hbdata(t_hbdata *p_hb, int nframes)
2342 if (nframes >= p_hb->max_frames)
2344 p_hb->max_frames += 4096;
2345 srenew(p_hb->nhb, p_hb->max_frames);
2346 srenew(p_hb->ndist, p_hb->max_frames);
2347 srenew(p_hb->n_bound, p_hb->max_frames);
2348 srenew(p_hb->nhx, p_hb->max_frames);
2351 srenew(p_hb->danr, p_hb->max_frames);
2353 std::memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
2354 std::memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
2355 p_hb->nhb[nframes] = 0;
2356 p_hb->ndist[nframes] = 0;
2359 p_hb->nframes = nframes;
2361 std::memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
2364 int gmx_hbond(int argc, char *argv[])
2366 const char *desc[] = {
2367 "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
2368 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
2369 "(zero is extended) and the distance Donor - Acceptor",
2370 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
2371 "OH and NH groups are regarded as donors, O is an acceptor always,",
2372 "N is an acceptor by default, but this can be switched using",
2373 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
2374 "to the first preceding non-hydrogen atom.[PAR]",
2376 "You need to specify two groups for analysis, which must be either",
2377 "identical or non-overlapping. All hydrogen bonds between the two",
2378 "groups are analyzed.[PAR]",
2380 "If you set [TT]-shell[tt], you will be asked for an additional index group",
2381 "which should contain exactly one atom. In this case, only hydrogen",
2382 "bonds between atoms within the shell distance from the one atom are",
2385 "With option -ac, rate constants for hydrogen bonding can be derived with the",
2386 "model of Luzar and Chandler (Nature 379:55, 1996; J. Chem. Phys. 113:23, 2000).",
2387 "If contact kinetics are analyzed by using the -contact option, then",
2388 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
2389 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
2390 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
2391 "See mentioned literature for more details and definitions."
2394 /* "It is also possible to analyse specific hydrogen bonds with",
2395 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
2396 "Donor Hydrogen Acceptor, in the following way::",
2403 "Note that the triplets need not be on separate lines.",
2404 "Each atom triplet specifies a hydrogen bond to be analyzed,",
2405 "note also that no check is made for the types of atoms.[PAR]",
2410 " * [TT]-num[tt]: number of hydrogen bonds as a function of time.",
2411 " * [TT]-ac[tt]: average over all autocorrelations of the existence",
2412 " functions (either 0 or 1) of all hydrogen bonds.",
2413 " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
2414 " * [TT]-ang[tt]: angle distribution of all hydrogen bonds.",
2415 " * [TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
2416 " where n and n+i stand for residue numbers and i ranges from 0 to 6.",
2417 " This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
2418 " with helices in proteins.",
2419 " * [TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
2420 " for selected groups, all hydrogen bonded atoms from all groups and",
2421 " all solvent atoms involved in insertion.",
2422 " * [TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
2423 " frames, this also contains information on solvent insertion",
2424 " into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
2426 " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
2427 " each timeframe. This is especially useful when using [TT]-shell[tt].",
2428 " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
2429 " compare results to Raman Spectroscopy.",
2431 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
2432 "require an amount of memory proportional to the total numbers of donors",
2433 "times the total number of acceptors in the selected group(s)."
2436 static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
2437 static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15;
2438 static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
2439 static int nDump = 0;
2440 static int nThreads = 0;
2442 static gmx_bool bContact = FALSE;
2446 { "-a", FALSE, etREAL, {&acut},
2447 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
2448 { "-r", FALSE, etREAL, {&rcut},
2449 "Cutoff radius (nm, X - Acceptor, see next option)" },
2450 { "-da", FALSE, etBOOL, {&bDA},
2451 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
2452 { "-r2", FALSE, etREAL, {&r2cut},
2453 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
2454 { "-abin", FALSE, etREAL, {&abin},
2455 "Binwidth angle distribution (degrees)" },
2456 { "-rbin", FALSE, etREAL, {&rbin},
2457 "Binwidth distance distribution (nm)" },
2458 { "-nitacc", FALSE, etBOOL, {&bNitAcc},
2459 "Regard nitrogen atoms as acceptors" },
2460 { "-contact", FALSE, etBOOL, {&bContact},
2461 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
2462 { "-shell", FALSE, etREAL, {&rshell},
2463 "when > 0, only calculate hydrogen bonds within # nm shell around "
2465 { "-fitstart", FALSE, etREAL, {&fit_start},
2466 "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]" },
2467 { "-fitend", FALSE, etREAL, {&fit_end},
2468 "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])" },
2469 { "-temp", FALSE, etREAL, {&temp},
2470 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
2471 { "-dump", FALSE, etINT, {&nDump},
2472 "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
2473 { "-max_hb", FALSE, etREAL, {&maxnhb},
2474 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
2475 { "-merge", FALSE, etBOOL, {&bMerge},
2476 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
2478 { "-nthreads", FALSE, etINT, {&nThreads},
2479 "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
2482 const char *bugs[] = {
2483 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
2486 { efTRX, "-f", NULL, ffREAD },
2487 { efTPR, NULL, NULL, ffREAD },
2488 { efNDX, NULL, NULL, ffOPTRD },
2489 /* { efNDX, "-sel", "select", ffOPTRD },*/
2490 { efXVG, "-num", "hbnum", ffWRITE },
2491 { efLOG, "-g", "hbond", ffOPTWR },
2492 { efXVG, "-ac", "hbac", ffOPTWR },
2493 { efXVG, "-dist", "hbdist", ffOPTWR },
2494 { efXVG, "-ang", "hbang", ffOPTWR },
2495 { efXVG, "-hx", "hbhelix", ffOPTWR },
2496 { efNDX, "-hbn", "hbond", ffOPTWR },
2497 { efXPM, "-hbm", "hbmap", ffOPTWR },
2498 { efXVG, "-don", "donor", ffOPTWR },
2499 { efXVG, "-dan", "danum", ffOPTWR },
2500 { efXVG, "-life", "hblife", ffOPTWR },
2501 { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
2504 #define NFILE asize(fnm)
2506 char hbmap [HB_NR] = { ' ', 'o', '-', '*' };
2507 const char *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
2508 t_rgb hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
2510 t_trxstatus *status;
2514 int npargs, natoms, nframes = 0, shatom;
2520 real t, ccut, dist = 0.0, ang = 0.0;
2521 double max_nhb, aver_nhb, aver_dist;
2522 int h = 0, i = 0, j, k = 0, ogrp, nsel;
2524 int xj, yj, zj, aj, xjj, yjj, zjj;
2525 gmx_bool bSelected, bHBmap, bStop, bTwo, bBox, bTric;
2527 int grp, nabin, nrbin, resdist, ihb;
2530 FILE *fp, *fpnhb = NULL, *donor_properties = NULL;
2532 t_ncell *icell, *jcell;
2534 unsigned char *datable;
2535 gmx_output_env_t *oenv;
2536 int ii, hh, actual_nThreads;
2539 gmx_bool bEdge_yjj, bEdge_xjj;
2541 t_hbdata **p_hb = NULL; /* one per thread, then merge after the frame loop */
2542 int **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
2544 const bool bOMP = GMX_OPENMP;
2547 ppa = add_acf_pargs(&npargs, pa);
2549 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
2550 ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
2558 ccut = std::cos(acut*DEG2RAD);
2564 gmx_fatal(FARGS, "Can not analyze selected contacts.");
2568 gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
2572 /* Initiate main data structure! */
2573 bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
2574 opt2bSet("-life", NFILE, fnm) ||
2575 opt2bSet("-hbn", NFILE, fnm) ||
2576 opt2bSet("-hbm", NFILE, fnm));
2578 if (opt2bSet("-nhbdist", NFILE, fnm))
2580 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
2581 fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
2582 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
2583 xvgr_legend(fpnhb, asize(leg), leg, oenv);
2586 hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact);
2589 gmx::MDModules mdModules;
2590 t_inputrec *ir = mdModules.inputrec();
2591 read_tpx_top(ftp2fn(efTPR, NFILE, fnm), ir, box, &natoms, NULL, NULL, &top);
2593 snew(grpnames, grNR);
2596 /* Make Donor-Acceptor table */
2597 snew(datable, top.atoms.nr);
2601 /* analyze selected hydrogen bonds */
2602 printf("Select group with selected atoms:\n");
2603 get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
2604 1, &nsel, index, grpnames);
2607 gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
2608 "and therefore cannot contain triplets of "
2609 "Donor-Hydrogen-Acceptor", grpnames[0]);
2613 for (i = 0; (i < nsel); i += 3)
2615 int dd = index[0][i];
2616 int aa = index[0][i+2];
2617 /* int */ hh = index[0][i+1];
2618 add_dh (&hb->d, dd, hh, i, datable);
2619 add_acc(&hb->a, aa, i);
2620 /* Should this be here ? */
2621 snew(hb->d.dptr, top.atoms.nr);
2622 snew(hb->a.aptr, top.atoms.nr);
2623 add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact);
2625 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
2626 isize[0], grpnames[0]);
2630 /* analyze all hydrogen bonds: get group(s) */
2631 printf("Specify 2 groups to analyze:\n");
2632 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
2633 2, isize, index, grpnames);
2635 /* check if we have two identical or two non-overlapping groups */
2636 bTwo = isize[0] != isize[1];
2637 for (i = 0; (i < isize[0]) && !bTwo; i++)
2639 bTwo = index[0][i] != index[1][i];
2643 printf("Checking for overlap in atoms between %s and %s\n",
2644 grpnames[0], grpnames[1]);
2646 gen_datable(index[0], isize[0], datable, top.atoms.nr);
2648 for (i = 0; i < isize[1]; i++)
2650 if (ISINGRP(datable[index[1][i]]))
2652 gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
2653 grpnames[0], grpnames[1]);
2659 printf("Calculating %s "
2660 "between %s (%d atoms) and %s (%d atoms)\n",
2661 bContact ? "contacts" : "hydrogen bonds",
2662 grpnames[0], isize[0], grpnames[1], isize[1]);
2666 fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
2667 bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
2672 /* search donors and acceptors in groups */
2673 snew(datable, top.atoms.nr);
2674 for (i = 0; (i < grNR); i++)
2676 if ( ((i == gr0) && !bSelected ) ||
2677 ((i == gr1) && bTwo ))
2679 gen_datable(index[i], isize[i], datable, top.atoms.nr);
2682 search_acceptors(&top, isize[i], index[i], &hb->a, i,
2683 bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
2684 search_donors (&top, isize[i], index[i], &hb->d, i,
2685 TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
2689 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
2690 search_donors (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
2694 clear_datable_grp(datable, top.atoms.nr);
2699 printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
2701 snew(donors[gr0D], dons[gr0D].nrd);*/
2703 donor_properties = open_donor_properties_file(opt2fn_null("-don", NFILE, fnm), hb, oenv);
2707 printf("Making hbmap structure...");
2708 /* Generate hbond data structure */
2715 if (hb->d.nrd + hb->a.nra == 0)
2717 printf("No Donors or Acceptors found\n");
2724 printf("No Donors found\n");
2729 printf("No Acceptors found\n");
2735 gmx_fatal(FARGS, "Nothing to be done");
2744 /* get index group with atom for shell */
2747 printf("Select atom for shell (1 atom):\n");
2748 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
2749 1, &shisz, &shidx, &shgrpnm);
2752 printf("group contains %d atoms, should be 1 (one)\n", shisz);
2757 printf("Will calculate hydrogen bonds within a shell "
2758 "of %g nm around atom %i\n", rshell, shatom+1);
2761 /* Analyze trajectory */
2762 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
2763 if (natoms > top.atoms.nr)
2765 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
2766 top.atoms.nr, natoms);
2769 bBox = (ir->ePBC != epbcNONE);
2770 grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
2771 nabin = static_cast<int>(acut/abin);
2772 nrbin = static_cast<int>(rcut/rbin);
2773 snew(adist, nabin+1);
2774 snew(rdist, nrbin+1);
2777 #define __ADIST adist
2778 #define __RDIST rdist
2780 #else /* GMX_OPENMP ================================================== \
2781 * Set up the OpenMP stuff, |
2782 * like the number of threads and such |
2783 * Also start the parallel loop. |
2785 #define __ADIST p_adist[threadNr]
2786 #define __RDIST p_rdist[threadNr]
2787 #define __HBDATA p_hb[threadNr]
2791 bParallel = !bSelected;
2795 actual_nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2797 gmx_omp_set_num_threads(actual_nThreads);
2798 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
2803 actual_nThreads = 1;
2806 snew(p_hb, actual_nThreads);
2807 snew(p_adist, actual_nThreads);
2808 snew(p_rdist, actual_nThreads);
2809 for (i = 0; i < actual_nThreads; i++)
2812 snew(p_adist[i], nabin+1);
2813 snew(p_rdist[i], nrbin+1);
2815 p_hb[i]->max_frames = 0;
2816 p_hb[i]->nhb = NULL;
2817 p_hb[i]->ndist = NULL;
2818 p_hb[i]->n_bound = NULL;
2819 p_hb[i]->time = NULL;
2820 p_hb[i]->nhx = NULL;
2822 p_hb[i]->bHBmap = hb->bHBmap;
2823 p_hb[i]->bDAnr = hb->bDAnr;
2824 p_hb[i]->wordlen = hb->wordlen;
2825 p_hb[i]->nframes = hb->nframes;
2826 p_hb[i]->maxhydro = hb->maxhydro;
2827 p_hb[i]->danr = hb->danr;
2830 p_hb[i]->hbmap = hb->hbmap;
2831 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
2834 p_hb[i]->nrdist = 0;
2838 /* Make a thread pool here,
2839 * instead of forking anew at every frame. */
2841 #pragma omp parallel \
2843 private(j, h, ii, hh, \
2844 xi, yi, zi, xj, yj, zj, threadNr, \
2845 dist, ang, icell, jcell, \
2846 grp, ogrp, ai, aj, xjj, yjj, zjj, \
2849 bEdge_xjj, bEdge_yjj) \
2851 { /* Start of parallel region */
2852 #if !defined __clang_analyzer__ // clang complains about unused value.
2853 threadNr = gmx_omp_get_thread_num();
2859 bTric = bBox && TRICLINIC(box);
2865 sync_hbdata(p_hb[threadNr], nframes);
2867 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2873 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
2874 rshell, ngrid, grid);
2875 reset_nhbonds(&(hb->d));
2877 if (debug && bDebug)
2879 dump_grid(debug, ngrid, grid);
2882 add_frames(hb, nframes);
2883 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
2887 count_da_grid(ngrid, grid, hb->danr[nframes]);
2890 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2895 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
2905 /* Do not parallelize this just yet. */
2907 for (ii = 0; (ii < nsel); ii++)
2909 int dd = index[0][i];
2910 int aa = index[0][i+2];
2911 /* int */ hh = index[0][i+1];
2912 ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
2913 hbox, &dist, &ang, bDA, &h, bContact, bMerge);
2917 /* add to index if not already there */
2919 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact);
2923 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2925 } /* if (bSelected) */
2928 /* The outer grid loop will have to do for now. */
2929 #pragma omp for schedule(dynamic)
2930 for (xi = 0; xi < ngrid[XX]; xi++)
2934 for (yi = 0; (yi < ngrid[YY]); yi++)
2936 for (zi = 0; (zi < ngrid[ZZ]); zi++)
2939 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
2940 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
2942 icell = &(grid[zi][yi][xi].d[grp]);
2953 /* loop over all hydrogen atoms from group (grp)
2954 * in this gridcell (icell)
2956 for (ai = 0; (ai < icell->nr); ai++)
2958 i = icell->atoms[ai];
2960 /* loop over all adjacent gridcells (xj,yj,zj) */
2961 for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
2962 zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
2965 zj = grid_mod(zjj, ngrid[ZZ]);
2966 bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
2967 for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
2968 yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
2971 yj = grid_mod(yjj, ngrid[YY]);
2973 (yj == 0) || (yj == ngrid[YY] - 1) ||
2974 (zj == 0) || (zj == ngrid[ZZ] - 1);
2975 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
2976 xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
2979 xj = grid_mod(xjj, ngrid[XX]);
2980 jcell = &(grid[zj][yj][xj].a[ogrp]);
2981 /* loop over acceptor atoms from other group (ogrp)
2982 * in this adjacent gridcell (jcell)
2984 for (aj = 0; (aj < jcell->nr); aj++)
2986 j = jcell->atoms[aj];
2988 /* check if this once was a h-bond */
2989 ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
2990 hbox, &dist, &ang, bDA, &h, bContact, bMerge);
2994 /* add to index if not already there */
2996 add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact);
2998 /* make angle and distance distributions */
2999 if (ihb == hbHB && !bContact)
3003 gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
3006 __ADIST[static_cast<int>( ang/abin)]++;
3007 __RDIST[static_cast<int>(dist/rbin)]++;
3010 if (donor_index(&hb->d, grp, i) == NOTSET)
3012 gmx_fatal(FARGS, "Invalid donor %d", i);
3014 if (acceptor_index(&hb->a, ogrp, j) == NOTSET)
3016 gmx_fatal(FARGS, "Invalid acceptor %d", j);
3018 resdist = std::abs(top.atoms.atom[i].resind-top.atoms.atom[j].resind);
3019 if (resdist >= max_hx)
3023 __HBDATA->nhx[nframes][resdist]++;
3034 } /* for xi,yi,zi */
3037 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3039 } /* if (bSelected) {...} else */
3042 /* Better wait for all threads to finnish using x[] before updating it. */
3045 #pragma omp critical
3049 /* Sum up histograms and counts from p_hb[] into hb */
3052 hb->nhb[k] += p_hb[threadNr]->nhb[k];
3053 hb->ndist[k] += p_hb[threadNr]->ndist[k];
3054 for (j = 0; j < max_hx; j++)
3056 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
3060 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3063 /* Here are a handful of single constructs
3064 * to share the workload a bit. The most
3065 * important one is of course the last one,
3066 * where there's a potential bottleneck in form
3073 analyse_donor_properties(donor_properties, hb, k, t);
3075 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3084 do_nhb_dist(fpnhb, hb, t);
3087 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3094 trrStatus = (read_next_x(oenv, status, &t, x, box));
3097 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3106 #pragma omp critical
3108 hb->nrhb += p_hb[threadNr]->nrhb;
3109 hb->nrdist += p_hb[threadNr]->nrdist;
3112 /* Free parallel datastructures */
3113 sfree(p_hb[threadNr]->nhb);
3114 sfree(p_hb[threadNr]->ndist);
3115 sfree(p_hb[threadNr]->nhx);
3118 for (i = 0; i < nabin; i++)
3122 for (j = 0; j < actual_nThreads; j++)
3125 adist[i] += p_adist[j][i];
3128 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3131 for (i = 0; i <= nrbin; i++)
3135 for (j = 0; j < actual_nThreads; j++)
3137 rdist[i] += p_rdist[j][i];
3140 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3143 sfree(p_adist[threadNr]);
3144 sfree(p_rdist[threadNr]);
3146 } /* End of parallel region */
3153 if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
3155 gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
3158 free_grid(ngrid, &grid);
3162 if (donor_properties)
3164 xvgrclose(donor_properties);
3172 /* Compute maximum possible number of different hbonds */
3179 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
3186 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
3190 printf("Found %d different %s in trajectory\n"
3191 "Found %d different atom-pairs within %s distance\n",
3192 hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
3193 hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
3197 merge_hb(hb, bTwo, bContact);
3200 if (opt2bSet("-hbn", NFILE, fnm))
3202 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
3205 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3206 * to make the -hbn and -hmb output match eachother.
3207 * - Erik Marklund, May 30, 2006 */
3210 /* Print out number of hbonds and distances */
3213 fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
3214 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3216 snew(leg[0], STRLEN);
3217 snew(leg[1], STRLEN);
3218 sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
3219 sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
3220 xvgr_legend(fp, 2, (const char**)leg, oenv);
3224 for (i = 0; (i < nframes); i++)
3226 fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
3227 aver_nhb += hb->nhb[i];
3228 aver_dist += hb->ndist[i];
3231 aver_nhb /= nframes;
3232 aver_dist /= nframes;
3233 /* Print HB distance distribution */
3234 if (opt2bSet("-dist", NFILE, fnm))
3239 for (i = 0; i < nrbin; i++)
3244 fp = xvgropen(opt2fn("-dist", NFILE, fnm),
3245 "Hydrogen Bond Distribution",
3247 "Donor - Acceptor Distance (nm)" :
3248 "Hydrogen - Acceptor Distance (nm)", "", oenv);
3249 for (i = 0; i < nrbin; i++)
3251 fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*sum));
3256 /* Print HB angle distribution */
3257 if (opt2bSet("-ang", NFILE, fnm))
3262 for (i = 0; i < nabin; i++)
3267 fp = xvgropen(opt2fn("-ang", NFILE, fnm),
3268 "Hydrogen Bond Distribution",
3269 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
3270 for (i = 0; i < nabin; i++)
3272 fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*sum));
3277 /* Print HB in alpha-helix */
3278 if (opt2bSet("-hx", NFILE, fnm))
3280 fp = xvgropen(opt2fn("-hx", NFILE, fnm),
3281 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
3282 xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
3283 for (i = 0; i < nframes; i++)
3285 fprintf(fp, "%10g", hb->time[i]);
3286 for (j = 0; j < max_hx; j++)
3288 fprintf(fp, " %6d", hb->nhx[i][j]);
3295 printf("Average number of %s per timeframe %.3f out of %g possible\n",
3296 bContact ? "contacts" : "hbonds",
3297 bContact ? aver_dist : aver_nhb, max_nhb);
3299 /* Do Autocorrelation etc. */
3303 Added support for -contact in ac and hbm calculations below.
3304 - Erik Marklund, May 29, 2006
3306 if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
3308 please_cite(stdout, "Spoel2006b");
3310 if (opt2bSet("-ac", NFILE, fnm))
3312 do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
3313 bMerge, bContact, fit_start, temp, r2cut > 0, oenv,
3316 if (opt2bSet("-life", NFILE, fnm))
3318 do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
3320 if (opt2bSet("-hbm", NFILE, fnm))
3323 int id, ia, hh, x, y;
3324 mat.flags = mat.y0 = 0;
3326 if ((nframes > 0) && (hb->nrhb > 0))
3331 snew(mat.matrix, mat.nx);
3332 for (x = 0; (x < mat.nx); x++)
3334 snew(mat.matrix[x], mat.ny);
3337 for (id = 0; (id < hb->d.nrd); id++)
3339 for (ia = 0; (ia < hb->a.nra); ia++)
3341 for (hh = 0; (hh < hb->maxhydro); hh++)
3343 if (hb->hbmap[id][ia])
3345 if (ISHB(hb->hbmap[id][ia]->history[hh]))
3347 for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
3349 int nn0 = hb->hbmap[id][ia]->n0;
3350 range_check(y, 0, mat.ny);
3351 mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
3359 mat.axis_x = hb->time;
3360 snew(mat.axis_y, mat.ny);
3361 for (j = 0; j < mat.ny; j++)
3365 sprintf(mat.title, bContact ? "Contact Existence Map" :
3366 "Hydrogen Bond Existence Map");
3367 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
3368 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
3369 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
3370 mat.bDiscrete = TRUE;
3372 snew(mat.map, mat.nmap);
3373 for (i = 0; i < mat.nmap; i++)
3375 mat.map[i].code.c1 = hbmap[i];
3376 mat.map[i].desc = hbdesc[i];
3377 mat.map[i].rgb = hbrgb[i];
3379 fp = opt2FILE("-hbm", NFILE, fnm, "w");
3380 write_xpm_m(fp, mat);
3382 for (x = 0; x < mat.nx; x++)
3384 sfree(mat.matrix[x]);
3392 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
3403 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
3405 fp = xvgropen(opt2fn("-dan", NFILE, fnm),
3406 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
3407 nleg = (bTwo ? 2 : 1)*2;
3408 snew(legnames, nleg);
3410 for (j = 0; j < grNR; j++)
3412 if (USE_THIS_GROUP(j) )
3414 sprintf(buf, "Donors %s", grpnames[j]);
3415 legnames[i++] = gmx_strdup(buf);
3416 sprintf(buf, "Acceptors %s", grpnames[j]);
3417 legnames[i++] = gmx_strdup(buf);
3422 gmx_incons("number of legend entries");
3424 xvgr_legend(fp, nleg, (const char**)legnames, oenv);
3425 for (i = 0; i < nframes; i++)
3427 fprintf(fp, "%10g", hb->time[i]);
3428 for (j = 0; (j < grNR); j++)
3430 if (USE_THIS_GROUP(j) )
3432 fprintf(fp, " %6d", hb->danr[i][j]);