Introduce SimulatorBuilder
[gromacs.git] / src / gromacs / gmxana / gmx_hbond.cpp
blob9fa96015cd0d27c8824df92938669262c9254f6b
1 /*
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,2017,2018,2019, 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.
37 #include "gmxpre.h"
39 #include "config.h"
41 #include <cmath>
42 #include <cstdlib>
43 #include <cstring>
45 #include <algorithm>
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/mdtypes/inputrec.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/topology/ifunc.h"
65 #include "gromacs/topology/index.h"
66 #include "gromacs/topology/topology.h"
67 #include "gromacs/utility/arraysize.h"
68 #include "gromacs/utility/cstringutil.h"
69 #include "gromacs/utility/exceptions.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/futil.h"
72 #include "gromacs/utility/gmxomp.h"
73 #include "gromacs/utility/pleasecite.h"
74 #include "gromacs/utility/programcontext.h"
75 #include "gromacs/utility/smalloc.h"
76 #include "gromacs/utility/snprintf.h"
78 #define max_hx 7
79 typedef int t_hx[max_hx];
80 #define NRHXTYPES max_hx
81 static const char *hxtypenames[NRHXTYPES] =
82 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
83 #define MAXHH 4
85 static const int NOTSET = -49297;
87 /* -----------------------------------------*/
89 enum {
90 gr0, gr1, grI, grNR
92 enum {
93 hbNo, hbDist, hbHB, hbNR, hbR2
95 static const unsigned char c_acceptorMask = (1 << 0);
96 static const unsigned char c_donorMask = (1 << 1);
97 static const unsigned char c_inGroupMask = (1 << 2);
100 static const char *grpnames[grNR] = {"0", "1", "I" };
102 static gmx_bool bDebug = FALSE;
104 #define HB_NO 0
105 #define HB_YES 1<<0
106 #define HB_INS 1<<1
107 #define HB_YESINS (HB_YES|HB_INS)
108 #define HB_NR (1<<2)
109 #define MAXHYDRO 4
111 #define ISHB(h) ((h) & 2)
112 #define ISDIST(h) ((h) & 1)
113 #define ISDIST2(h) ((h) & 4)
114 #define ISACC(h) ((h) & c_acceptorMask)
115 #define ISDON(h) ((h) & c_donorMask)
116 #define ISINGRP(h) ((h) & c_inGroupMask)
118 typedef struct {
119 int nr;
120 int maxnr;
121 int *atoms;
122 } t_ncell;
124 typedef struct {
125 t_ncell d[grNR];
126 t_ncell a[grNR];
127 } t_gridcell;
129 typedef int t_icell[grNR];
130 typedef int h_id[MAXHYDRO];
132 typedef struct {
133 int history[MAXHYDRO];
134 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
135 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
137 /* Bitmask array which tells whether a hbond is present
138 * at a given time. Either of these may be NULL
140 int n0; /* First frame a HB was found */
141 int nframes, maxframes; /* Amount of frames in this hbond */
142 unsigned int **h;
143 unsigned int **g;
144 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
145 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
146 * acceptor distance is less than the user-specified distance (typically
147 * 0.35 nm).
149 } t_hbond;
151 typedef struct {
152 int nra, max_nra;
153 int *acc; /* Atom numbers of the acceptors */
154 int *grp; /* Group index */
155 int *aptr; /* Map atom number to acceptor index */
156 } t_acceptors;
158 typedef struct {
159 int nrd, max_nrd;
160 int *don; /* Atom numbers of the donors */
161 int *grp; /* Group index */
162 int *dptr; /* Map atom number to donor index */
163 int *nhydro; /* Number of hydrogens for each donor */
164 h_id *hydro; /* The atom numbers of the hydrogens */
165 h_id *nhbonds; /* The number of HBs per H at current */
166 } t_donors;
168 typedef struct {
169 gmx_bool bHBmap, bDAnr;
170 int wordlen;
171 /* The following arrays are nframes long */
172 int nframes, max_frames, maxhydro;
173 int *nhb, *ndist;
174 h_id *n_bound;
175 real *time;
176 t_icell *danr;
177 t_hx *nhx;
178 /* These structures are initialized from the topology at start up */
179 t_donors d;
180 t_acceptors a;
181 /* This holds a matrix with all possible hydrogen bonds */
182 int nrhb, nrdist;
183 t_hbond ***hbmap;
184 } t_hbdata;
186 /* Changed argument 'bMerge' into 'oneHB' below,
187 * since -contact should cause maxhydro to be 1,
188 * not just -merge.
189 * - Erik Marklund May 29, 2006
192 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB)
194 t_hbdata *hb;
196 snew(hb, 1);
197 hb->wordlen = 8*sizeof(unsigned int);
198 hb->bHBmap = bHBmap;
199 hb->bDAnr = bDAnr;
200 if (oneHB)
202 hb->maxhydro = 1;
204 else
206 hb->maxhydro = MAXHYDRO;
208 return hb;
211 static void mk_hbmap(t_hbdata *hb)
213 int i, j;
215 snew(hb->hbmap, hb->d.nrd);
216 for (i = 0; (i < hb->d.nrd); i++)
218 snew(hb->hbmap[i], hb->a.nra);
219 if (hb->hbmap[i] == nullptr)
221 gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
223 for (j = 0; j < hb->a.nra; j++)
225 hb->hbmap[i][j] = nullptr;
230 static void add_frames(t_hbdata *hb, int nframes)
232 if (nframes >= hb->max_frames)
234 hb->max_frames += 4096;
235 srenew(hb->time, hb->max_frames);
236 srenew(hb->nhb, hb->max_frames);
237 srenew(hb->ndist, hb->max_frames);
238 srenew(hb->n_bound, hb->max_frames);
239 srenew(hb->nhx, hb->max_frames);
240 if (hb->bDAnr)
242 srenew(hb->danr, hb->max_frames);
245 hb->nframes = nframes;
248 #define OFFSET(frame) ((frame) / 32)
249 #define MASK(frame) (1 << ((frame) % 32))
251 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
253 if (bValue)
255 hbexist[OFFSET(frame)] |= MASK(frame);
257 else
259 hbexist[OFFSET(frame)] &= ~MASK(frame);
263 static gmx_bool is_hb(const unsigned int hbexist[], int frame)
265 return (hbexist[OFFSET(frame)] & MASK(frame)) != 0;
268 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
270 unsigned int *ghptr = nullptr;
272 if (ihb == hbHB)
274 ghptr = hb->hbmap[id][ia]->h[ih];
276 else if (ihb == hbDist)
278 ghptr = hb->hbmap[id][ia]->g[ih];
280 else
282 gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
285 _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
288 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb)
290 int i, j, n;
291 t_hbond *hb = hbd->hbmap[id][ia];
292 int maxhydro = std::min(hbd->maxhydro, hbd->d.nhydro[id]);
293 int wlen = hbd->wordlen;
294 int delta = 32*wlen;
296 if (!hb->h[0])
298 hb->n0 = frame;
299 hb->maxframes = delta;
300 for (i = 0; (i < maxhydro); i++)
302 snew(hb->h[i], hb->maxframes/wlen);
303 snew(hb->g[i], hb->maxframes/wlen);
306 else
308 hb->nframes = frame-hb->n0;
309 /* We need a while loop here because hbonds may be returning
310 * after a long time.
312 while (hb->nframes >= hb->maxframes)
314 n = hb->maxframes + delta;
315 for (i = 0; (i < maxhydro); i++)
317 srenew(hb->h[i], n/wlen);
318 srenew(hb->g[i], n/wlen);
319 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
321 hb->h[i][j] = 0;
322 hb->g[i][j] = 0;
326 hb->maxframes = n;
329 if (frame >= 0)
331 set_hb(hbd, id, h, ia, frame, ihb);
336 static void inc_nhbonds(t_donors *ddd, int d, int h)
338 int j;
339 int dptr = ddd->dptr[d];
341 for (j = 0; (j < ddd->nhydro[dptr]); j++)
343 if (ddd->hydro[dptr][j] == h)
345 ddd->nhbonds[dptr][j]++;
346 break;
349 if (j == ddd->nhydro[dptr])
351 gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
355 static int _acceptor_index(t_acceptors *a, int grp, int i,
356 const char *file, int line)
358 int ai = a->aptr[i];
360 if (a->grp[ai] != grp)
362 if (debug && bDebug)
364 fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
365 ai, a->grp[ai], grp, file, line);
367 return NOTSET;
369 else
371 return ai;
374 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
376 static int _donor_index(t_donors *d, int grp, int i, const char *file, int line)
378 int di = d->dptr[i];
380 if (di == NOTSET)
382 return NOTSET;
385 if (d->grp[di] != grp)
387 if (debug && bDebug)
389 fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
390 di, d->grp[di], grp, file, line);
392 return NOTSET;
394 else
396 return di;
399 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
401 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
403 /* g_hbond doesn't allow overlapping groups */
404 if (grpa != grpd)
406 return FALSE;
408 return
409 donor_index(&hb->d, grpd, a) != NOTSET
410 && acceptor_index(&hb->a, grpa, d) != NOTSET;
414 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
415 int frame, gmx_bool bMerge, int ihb, gmx_bool bContact)
417 int k, id, ia, hh;
418 gmx_bool daSwap = FALSE;
420 if ((id = hb->d.dptr[d]) == NOTSET)
422 gmx_fatal(FARGS, "No donor atom %d", d+1);
424 else if (grpd != hb->d.grp[id])
426 gmx_fatal(FARGS, "Inconsistent donor groups, %d instead of %d, atom %d",
427 grpd, hb->d.grp[id], d+1);
429 if ((ia = hb->a.aptr[a]) == NOTSET)
431 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
433 else if (grpa != hb->a.grp[ia])
435 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d instead of %d, atom %d",
436 grpa, hb->a.grp[ia], a+1);
439 if (bMerge)
442 if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
443 /* Then swap identity so that the id of d is lower then that of a.
445 * This should really be redundant by now, as is_hbond() now ought to return
446 * hbNo in the cases where this conditional is TRUE. */
448 daSwap = TRUE;
449 k = d;
450 d = a;
451 a = k;
453 /* Now repeat donor/acc check. */
454 if ((id = hb->d.dptr[d]) == NOTSET)
456 gmx_fatal(FARGS, "No donor atom %d", d+1);
458 else if (grpd != hb->d.grp[id])
460 gmx_fatal(FARGS, "Inconsistent donor groups, %d instead of %d, atom %d",
461 grpd, hb->d.grp[id], d+1);
463 if ((ia = hb->a.aptr[a]) == NOTSET)
465 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
467 else if (grpa != hb->a.grp[ia])
469 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d instead of %d, atom %d",
470 grpa, hb->a.grp[ia], a+1);
475 if (hb->hbmap)
477 /* Loop over hydrogens to find which hydrogen is in this particular HB */
478 if ((ihb == hbHB) && !bMerge && !bContact)
480 for (k = 0; (k < hb->d.nhydro[id]); k++)
482 if (hb->d.hydro[id][k] == h)
484 break;
487 if (k == hb->d.nhydro[id])
489 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
490 d+1, h+1, a+1);
493 else
495 k = 0;
498 if (hb->bHBmap)
501 #pragma omp critical
505 if (hb->hbmap[id][ia] == nullptr)
507 snew(hb->hbmap[id][ia], 1);
508 snew(hb->hbmap[id][ia]->h, hb->maxhydro);
509 snew(hb->hbmap[id][ia]->g, hb->maxhydro);
511 add_ff(hb, id, k, ia, frame, ihb);
513 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
517 /* Strange construction with frame >=0 is a relic from old code
518 * for selected hbond analysis. It may be necessary again if that
519 * is made to work again.
521 if (frame >= 0)
523 hh = hb->hbmap[id][ia]->history[k];
524 if (ihb == hbHB)
526 hb->nhb[frame]++;
527 if (!(ISHB(hh)))
529 hb->hbmap[id][ia]->history[k] = hh | 2;
530 hb->nrhb++;
533 else
535 if (ihb == hbDist)
537 hb->ndist[frame]++;
538 if (!(ISDIST(hh)))
540 hb->hbmap[id][ia]->history[k] = hh | 1;
541 hb->nrdist++;
547 else
549 if (frame >= 0)
551 if (ihb == hbHB)
553 hb->nhb[frame]++;
555 else
557 if (ihb == hbDist)
559 hb->ndist[frame]++;
564 if (bMerge && daSwap)
566 h = hb->d.hydro[id][0];
568 /* Increment number if HBonds per H */
569 if (ihb == hbHB && !bContact)
571 inc_nhbonds(&(hb->d), d, h);
575 static char *mkatomname(const t_atoms *atoms, int i)
577 static char buf[32];
578 int rnr;
580 rnr = atoms->atom[i].resind;
581 sprintf(buf, "%4s%d%-4s",
582 *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
584 return buf;
587 static void gen_datable(int *index, int isize, unsigned char *datable, int natoms)
589 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
590 int i;
592 for (i = 0; i < isize; i++)
594 if (index[i] >= natoms)
596 gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
598 datable[index[i]] |= c_inGroupMask;
602 static void clear_datable_grp(unsigned char *datable, int size)
604 /* Clears group information from the table */
605 int i;
606 if (size > 0)
608 for (i = 0; i < size; i++)
610 datable[i] &= ~c_inGroupMask;
615 static void add_acc(t_acceptors *a, int ia, int grp)
617 if (a->nra >= a->max_nra)
619 a->max_nra += 16;
620 srenew(a->acc, a->max_nra);
621 srenew(a->grp, a->max_nra);
623 a->grp[a->nra] = grp;
624 a->acc[a->nra++] = ia;
627 static void search_acceptors(const t_topology *top, int isize,
628 const int *index, t_acceptors *a, int grp,
629 gmx_bool bNitAcc,
630 gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
632 int i, n;
634 if (bDoIt)
636 for (i = 0; (i < isize); i++)
638 n = index[i];
639 if ((bContact ||
640 (((*top->atoms.atomname[n])[0] == 'O') ||
641 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
642 ISINGRP(datable[n]))
644 datable[n] |= c_acceptorMask;
645 add_acc(a, n, grp);
649 snew(a->aptr, top->atoms.nr);
650 for (i = 0; (i < top->atoms.nr); i++)
652 a->aptr[i] = NOTSET;
654 for (i = 0; (i < a->nra); i++)
656 a->aptr[a->acc[i]] = i;
660 static void add_h2d(int id, int ih, t_donors *ddd)
662 int i;
664 for (i = 0; (i < ddd->nhydro[id]); i++)
666 if (ddd->hydro[id][i] == ih)
668 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
669 ddd->don[id], ih);
670 break;
673 if (i == ddd->nhydro[id])
675 if (ddd->nhydro[id] >= MAXHYDRO)
677 gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
678 ddd->don[id], MAXHYDRO);
680 ddd->hydro[id][i] = ih;
681 ddd->nhydro[id]++;
685 static void add_dh(t_donors *ddd, int id, int ih, int grp, const unsigned char *datable)
687 int i;
689 if (!datable || ISDON(datable[id]))
691 if (ddd->dptr[id] == NOTSET) /* New donor */
693 i = ddd->nrd;
694 ddd->dptr[id] = i;
696 else
698 i = ddd->dptr[id];
701 if (i == ddd->nrd)
703 if (ddd->nrd >= ddd->max_nrd)
705 ddd->max_nrd += 128;
706 srenew(ddd->don, ddd->max_nrd);
707 srenew(ddd->nhydro, ddd->max_nrd);
708 srenew(ddd->hydro, ddd->max_nrd);
709 srenew(ddd->nhbonds, ddd->max_nrd);
710 srenew(ddd->grp, ddd->max_nrd);
712 ddd->don[ddd->nrd] = id;
713 ddd->nhydro[ddd->nrd] = 0;
714 ddd->grp[ddd->nrd] = grp;
715 ddd->nrd++;
717 else
719 ddd->don[i] = id;
721 add_h2d(i, ih, ddd);
723 else
724 if (datable)
726 printf("Warning: Atom %d is not in the d/a-table!\n", id);
730 static void search_donors(const t_topology *top, int isize, const int *index,
731 t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
732 unsigned char *datable)
734 int i, j;
735 t_functype func_type;
736 int nr1, nr2, nr3;
738 if (!ddd->dptr)
740 snew(ddd->dptr, top->atoms.nr);
741 for (i = 0; (i < top->atoms.nr); i++)
743 ddd->dptr[i] = NOTSET;
747 if (bContact)
749 if (bDoIt)
751 for (i = 0; (i < isize); i++)
753 datable[index[i]] |= c_donorMask;
754 add_dh(ddd, index[i], -1, grp, datable);
758 else
760 for (func_type = 0; (func_type < F_NRE); func_type++)
762 const t_ilist *interaction = &(top->idef.il[func_type]);
763 if (func_type == F_POSRES || func_type == F_FBPOSRES)
765 /* The ilist looks strange for posre. Bug in grompp?
766 * We don't need posre interactions for hbonds anyway.*/
767 continue;
769 for (i = 0; i < interaction->nr;
770 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
772 /* next function */
773 if (func_type != top->idef.functype[interaction->iatoms[i]])
775 fprintf(stderr, "Error in func_type %s",
776 interaction_function[func_type].longname);
777 continue;
780 /* check out this functype */
781 if (func_type == F_SETTLE)
783 nr1 = interaction->iatoms[i+1];
784 nr2 = interaction->iatoms[i+2];
785 nr3 = interaction->iatoms[i+3];
787 if (ISINGRP(datable[nr1]))
789 if (ISINGRP(datable[nr2]))
791 datable[nr1] |= c_donorMask;
792 add_dh(ddd, nr1, nr1+1, grp, datable);
794 if (ISINGRP(datable[nr3]))
796 datable[nr1] |= c_donorMask;
797 add_dh(ddd, nr1, nr1+2, grp, datable);
801 else if (IS_CHEMBOND(func_type))
803 for (j = 0; j < 2; j++)
805 nr1 = interaction->iatoms[i+1+j];
806 nr2 = interaction->iatoms[i+2-j];
807 if ((*top->atoms.atomname[nr1][0] == 'H') &&
808 ((*top->atoms.atomname[nr2][0] == 'O') ||
809 (*top->atoms.atomname[nr2][0] == 'N')) &&
810 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
812 datable[nr2] |= c_donorMask;
813 add_dh(ddd, nr2, nr1, grp, datable);
819 #ifdef SAFEVSITES
820 for (func_type = 0; func_type < F_NRE; func_type++)
822 interaction = &top->idef.il[func_type];
823 for (i = 0; i < interaction->nr;
824 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
826 /* next function */
827 if (func_type != top->idef.functype[interaction->iatoms[i]])
829 gmx_incons("function type in search_donors");
832 if (interaction_function[func_type].flags & IF_VSITE)
834 nr1 = interaction->iatoms[i+1];
835 if (*top->atoms.atomname[nr1][0] == 'H')
837 nr2 = nr1-1;
838 stop = FALSE;
839 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
841 if (nr2)
843 nr2--;
845 else
847 stop = TRUE;
850 if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
851 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
852 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
854 datable[nr2] |= c_donorMask;
855 add_dh(ddd, nr2, nr1, grp, datable);
861 #endif
865 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
867 t_gridcell ***grid;
868 int i, y, z;
870 if (bBox)
872 for (i = 0; i < DIM; i++)
874 ngrid[i] = static_cast<int>(box[i][i]/(1.2*rcut));
878 if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
880 for (i = 0; i < DIM; i++)
882 ngrid[i] = 1;
885 else
887 printf("\nWill do grid-search on %dx%dx%d grid, rcut=%3.8f\n",
888 ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
890 if (((ngrid[XX]*ngrid[YY]*ngrid[ZZ]) * sizeof(grid)) > INT_MAX)
892 gmx_fatal(FARGS, "Failed to allocate memory for %d x %d x %d grid points, which is larger than the maximum of %zu. "
893 "You are likely either using a box that is too large (box dimensions are %3.8f nm x%3.8f nm x%3.8f nm) or a cutoff (%3.8f nm) that is too small.",
894 ngrid[XX], ngrid[YY], ngrid[ZZ], INT_MAX/sizeof(grid), box[XX][XX], box[YY][YY], box[ZZ][ZZ], rcut);
896 snew(grid, ngrid[ZZ]);
897 for (z = 0; z < ngrid[ZZ]; z++)
899 snew((grid)[z], ngrid[YY]);
900 for (y = 0; y < ngrid[YY]; y++)
902 snew((grid)[z][y], ngrid[XX]);
905 return grid;
908 static void reset_nhbonds(t_donors *ddd)
910 int i, j;
912 for (i = 0; (i < ddd->nrd); i++)
914 for (j = 0; (j < MAXHH); j++)
916 ddd->nhbonds[i][j] = 0;
921 static void pbc_correct_gem(rvec dx, matrix box, const rvec hbox);
922 static void pbc_in_gridbox(rvec dx, matrix box);
924 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
925 gmx_bool bBox, matrix box, rvec hbox,
926 real rcut, real rshell,
927 ivec ngrid, t_gridcell ***grid)
929 int i, m, gr, xi, yi, zi, nr;
930 int *ad;
931 ivec grididx;
932 rvec invdelta, dshell;
933 t_ncell *newgrid;
934 gmx_bool bDoRshell, bInShell;
935 real rshell2 = 0;
936 int gx, gy, gz;
937 int dum = -1;
939 bDoRshell = (rshell > 0);
940 rshell2 = gmx::square(rshell);
941 bInShell = TRUE;
943 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
944 DBB(dum);
945 for (m = 0; m < DIM; m++)
947 hbox[m] = box[m][m]*0.5;
948 if (bBox)
950 invdelta[m] = ngrid[m]/box[m][m];
951 if (1/invdelta[m] < rcut)
953 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
954 "%s can not handle this situation, sorry.\n",
955 gmx::getProgramContext().displayName());
958 else
960 invdelta[m] = 0;
963 grididx[XX] = 0;
964 grididx[YY] = 0;
965 grididx[ZZ] = 0;
966 DBB(dum);
967 /* resetting atom counts */
968 for (gr = 0; (gr < grNR); gr++)
970 for (zi = 0; zi < ngrid[ZZ]; zi++)
972 for (yi = 0; yi < ngrid[YY]; yi++)
974 for (xi = 0; xi < ngrid[XX]; xi++)
976 grid[zi][yi][xi].d[gr].nr = 0;
977 grid[zi][yi][xi].a[gr].nr = 0;
981 DBB(dum);
983 /* put atoms in grid cells */
984 for (int acc = 0; acc < 2; acc++)
986 if (acc == 1)
988 nr = hb->a.nra;
989 ad = hb->a.acc;
991 else
993 nr = hb->d.nrd;
994 ad = hb->d.don;
996 DBB(acc);
997 for (i = 0; (i < nr); i++)
999 /* check if we are inside the shell */
1000 /* if bDoRshell=FALSE then bInShell=TRUE always */
1001 DBB(i);
1002 if (bDoRshell)
1004 bInShell = TRUE;
1005 rvec_sub(x[ad[i]], xshell, dshell);
1006 if (bBox)
1008 gmx_bool bDone = FALSE;
1009 while (!bDone)
1011 bDone = TRUE;
1012 for (m = DIM-1; m >= 0 && bInShell; m--)
1014 if (dshell[m] < -hbox[m])
1016 bDone = FALSE;
1017 rvec_inc(dshell, box[m]);
1019 if (dshell[m] >= hbox[m])
1021 bDone = FALSE;
1022 dshell[m] -= 2*hbox[m];
1026 for (m = DIM-1; m >= 0 && bInShell; m--)
1028 /* if we're outside the cube, we're outside the sphere also! */
1029 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1031 bInShell = FALSE;
1035 /* if we're inside the cube, check if we're inside the sphere */
1036 if (bInShell)
1038 bInShell = norm2(dshell) < rshell2;
1041 DBB(i);
1042 if (bInShell)
1044 if (bBox)
1046 pbc_in_gridbox(x[ad[i]], box);
1048 for (m = DIM-1; m >= 0; m--)
1049 { /* determine grid index of atom */
1050 grididx[m] = static_cast<int>(x[ad[i]][m]*invdelta[m]);
1051 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1055 gx = grididx[XX];
1056 gy = grididx[YY];
1057 gz = grididx[ZZ];
1058 range_check(gx, 0, ngrid[XX]);
1059 range_check(gy, 0, ngrid[YY]);
1060 range_check(gz, 0, ngrid[ZZ]);
1061 DBB(gx);
1062 DBB(gy);
1063 DBB(gz);
1064 /* add atom to grid cell */
1065 if (acc == 1)
1067 newgrid = &(grid[gz][gy][gx].a[gr]);
1069 else
1071 newgrid = &(grid[gz][gy][gx].d[gr]);
1073 if (newgrid->nr >= newgrid->maxnr)
1075 newgrid->maxnr += 10;
1076 DBB(newgrid->maxnr);
1077 srenew(newgrid->atoms, newgrid->maxnr);
1079 DBB(newgrid->nr);
1080 newgrid->atoms[newgrid->nr] = ad[i];
1081 newgrid->nr++;
1088 static void count_da_grid(const ivec ngrid, t_gridcell ***grid, t_icell danr)
1090 int gr, xi, yi, zi;
1092 for (gr = 0; (gr < grNR); gr++)
1094 danr[gr] = 0;
1095 for (zi = 0; zi < ngrid[ZZ]; zi++)
1097 for (yi = 0; yi < ngrid[YY]; yi++)
1099 for (xi = 0; xi < ngrid[XX]; xi++)
1101 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1108 /* The grid loop.
1109 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1110 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1111 * With a triclinic box all loops are 3 long, except when a cell is
1112 * located next to one of the box edges which is not parallel to the
1113 * x/y-plane, in that case all cells in a line or layer are searched.
1114 * This could be implemented slightly more efficient, but the code
1115 * would get much more complicated.
1117 static inline int grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1119 return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
1121 static inline int grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1123 return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1125 static inline int grid_mod(int j, int n)
1127 return (j+n) % (n);
1130 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1132 int gr, x, y, z, sum[grNR];
1134 fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1135 for (gr = 0; gr < grNR; gr++)
1137 sum[gr] = 0;
1138 fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1139 for (z = 0; z < ngrid[ZZ]; z += 2)
1141 fprintf(fp, "Z=%d,%d\n", z, z+1);
1142 for (y = 0; y < ngrid[YY]; y++)
1144 for (x = 0; x < ngrid[XX]; x++)
1146 fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1147 sum[gr] += grid[z][y][x].d[gr].nr;
1148 fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1149 sum[gr] += grid[z][y][x].a[gr].nr;
1152 fprintf(fp, " | ");
1153 if ( (z+1) < ngrid[ZZ])
1155 for (x = 0; x < ngrid[XX]; x++)
1157 fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1158 sum[gr] += grid[z+1][y][x].d[gr].nr;
1159 fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1160 sum[gr] += grid[z+1][y][x].a[gr].nr;
1163 fprintf(fp, "\n");
1167 fprintf(fp, "TOTALS:");
1168 for (gr = 0; gr < grNR; gr++)
1170 fprintf(fp, " %d=%d", gr, sum[gr]);
1172 fprintf(fp, "\n");
1175 /* New GMX record! 5 * in a row. Congratulations!
1176 * Sorry, only four left.
1178 static void free_grid(const ivec ngrid, t_gridcell ****grid)
1180 int y, z;
1181 t_gridcell ***g = *grid;
1183 for (z = 0; z < ngrid[ZZ]; z++)
1185 for (y = 0; y < ngrid[YY]; y++)
1187 sfree(g[z][y]);
1189 sfree(g[z]);
1191 sfree(g);
1192 g = nullptr;
1195 static void pbc_correct_gem(rvec dx, matrix box, const rvec hbox)
1197 int m;
1198 gmx_bool bDone = FALSE;
1199 while (!bDone)
1201 bDone = TRUE;
1202 for (m = DIM-1; m >= 0; m--)
1204 if (dx[m] < -hbox[m])
1206 bDone = FALSE;
1207 rvec_inc(dx, box[m]);
1209 if (dx[m] >= hbox[m])
1211 bDone = FALSE;
1212 rvec_dec(dx, box[m]);
1218 static void pbc_in_gridbox(rvec dx, matrix box)
1220 int m;
1221 gmx_bool bDone = FALSE;
1222 while (!bDone)
1224 bDone = TRUE;
1225 for (m = DIM-1; m >= 0; m--)
1227 if (dx[m] < 0)
1229 bDone = FALSE;
1230 rvec_inc(dx, box[m]);
1232 if (dx[m] >= box[m][m])
1234 bDone = FALSE;
1235 rvec_dec(dx, box[m]);
1241 /* Added argument r2cut, changed contact and implemented
1242 * use of second cut-off.
1243 * - Erik Marklund, June 29, 2006
1245 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1246 real rcut, real r2cut, real ccut,
1247 rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1248 real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1249 gmx_bool bContact, gmx_bool bMerge)
1251 int h, hh, id;
1252 rvec r_da, r_ha, r_dh;
1253 real rc2, r2c2, rda2, rha2, ca;
1254 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1255 gmx_bool daSwap = FALSE;
1257 if (d == a)
1259 return hbNo;
1262 if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1263 (acceptor_index(&hb->a, grpa, a) == NOTSET))
1265 return hbNo;
1268 rc2 = rcut*rcut;
1269 r2c2 = r2cut*r2cut;
1271 rvec_sub(x[d], x[a], r_da);
1272 /* Insert projection code here */
1274 if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1276 /* Then this hbond/contact will be found again, or it has already been found. */
1277 /*return hbNo;*/
1279 if (bBox)
1281 if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1282 { /* return hbNo; */
1283 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1285 pbc_correct_gem(r_da, box, hbox);
1287 rda2 = iprod(r_da, r_da);
1289 if (bContact)
1291 if (daSwap && grpa == grpd)
1293 return hbNo;
1295 if (rda2 <= rc2)
1297 return hbHB;
1299 else if (rda2 < r2c2)
1301 return hbDist;
1303 else
1305 return hbNo;
1308 *hhh = NOTSET;
1310 if (bDA && (rda2 > rc2))
1312 return hbNo;
1315 for (h = 0; (h < hb->d.nhydro[id]); h++)
1317 hh = hb->d.hydro[id][h];
1318 rha2 = rc2+1;
1319 if (!bDA)
1321 rvec_sub(x[hh], x[a], r_ha);
1322 if (bBox)
1324 pbc_correct_gem(r_ha, box, hbox);
1326 rha2 = iprod(r_ha, r_ha);
1329 if (bDA || (rha2 <= rc2))
1331 rvec_sub(x[d], x[hh], r_dh);
1332 if (bBox)
1334 pbc_correct_gem(r_dh, box, hbox);
1337 if (!bDA)
1339 HAinrange = TRUE;
1341 ca = cos_angle(r_dh, r_da);
1342 /* if angle is smaller, cos is larger */
1343 if (ca >= ccut)
1345 *hhh = hh;
1346 *d_ha = std::sqrt(bDA ? rda2 : rha2);
1347 *ang = std::acos(ca);
1348 return hbHB;
1352 if (bDA || HAinrange)
1354 return hbDist;
1356 else
1358 return hbNo;
1362 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1363 * Will do some more testing before removing the function entirely.
1364 * - Erik Marklund, MAY 10 2010 */
1365 static void do_merge(t_hbdata *hb, int ntmp,
1366 bool htmp[], bool gtmp[],
1367 t_hbond *hb0, t_hbond *hb1)
1369 /* Here we need to make sure we're treating periodicity in
1370 * the right way for the geminate recombination kinetics. */
1372 int m, mm, n00, n01, nn0, nnframes;
1374 /* Decide where to start from when merging */
1375 n00 = hb0->n0;
1376 n01 = hb1->n0;
1377 nn0 = std::min(n00, n01);
1378 nnframes = std::max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1379 /* Initiate tmp arrays */
1380 for (m = 0; (m < ntmp); m++)
1382 htmp[m] = false;
1383 gtmp[m] = false;
1385 /* Fill tmp arrays with values due to first HB */
1386 /* Once again '<' had to be replaced with '<='
1387 to catch the last frame in which the hbond
1388 appears.
1389 - Erik Marklund, June 1, 2006 */
1390 for (m = 0; (m <= hb0->nframes); m++)
1392 mm = m+n00-nn0;
1393 htmp[mm] = is_hb(hb0->h[0], m);
1395 for (m = 0; (m <= hb0->nframes); m++)
1397 mm = m+n00-nn0;
1398 gtmp[mm] = is_hb(hb0->g[0], m);
1400 /* Next HB */
1401 for (m = 0; (m <= hb1->nframes); m++)
1403 mm = m+n01-nn0;
1404 htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1405 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1407 /* Reallocate target array */
1408 if (nnframes > hb0->maxframes)
1410 srenew(hb0->h[0], 4+nnframes/hb->wordlen);
1411 srenew(hb0->g[0], 4+nnframes/hb->wordlen);
1414 /* Copy temp array to target array */
1415 for (m = 0; (m <= nnframes); m++)
1417 _set_hb(hb0->h[0], m, htmp[m]);
1418 _set_hb(hb0->g[0], m, gtmp[m]);
1421 /* Set scalar variables */
1422 hb0->n0 = nn0;
1423 hb0->maxframes = nnframes;
1426 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
1428 int i, inrnew, indnew, j, ii, jj, id, ia, ntmp;
1429 bool *htmp, *gtmp;
1430 t_hbond *hb0, *hb1;
1432 inrnew = hb->nrhb;
1433 indnew = hb->nrdist;
1435 /* Check whether donors are also acceptors */
1436 printf("Merging hbonds with Acceptor and Donor swapped\n");
1438 ntmp = 2*hb->max_frames;
1439 snew(gtmp, ntmp);
1440 snew(htmp, ntmp);
1441 for (i = 0; (i < hb->d.nrd); i++)
1443 fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
1444 fflush(stderr);
1445 id = hb->d.don[i];
1446 ii = hb->a.aptr[id];
1447 for (j = 0; (j < hb->a.nra); j++)
1449 ia = hb->a.acc[j];
1450 jj = hb->d.dptr[ia];
1451 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
1452 (!bTwo || (hb->d.grp[i] != hb->a.grp[j])))
1454 hb0 = hb->hbmap[i][j];
1455 hb1 = hb->hbmap[jj][ii];
1456 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
1458 do_merge(hb, ntmp, htmp, gtmp, hb0, hb1);
1459 if (ISHB(hb1->history[0]))
1461 inrnew--;
1463 else if (ISDIST(hb1->history[0]))
1465 indnew--;
1467 else
1468 if (bContact)
1470 gmx_incons("No contact history");
1472 else
1474 gmx_incons("Neither hydrogen bond nor distance");
1476 sfree(hb1->h[0]);
1477 sfree(hb1->g[0]);
1478 hb1->h[0] = nullptr;
1479 hb1->g[0] = nullptr;
1480 hb1->history[0] = hbNo;
1485 fprintf(stderr, "\n");
1486 printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
1487 printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
1488 hb->nrhb = inrnew;
1489 hb->nrdist = indnew;
1490 sfree(gtmp);
1491 sfree(htmp);
1494 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
1496 int i, j, k, n_bound[MAXHH], nbtot;
1498 /* Set array to 0 */
1499 for (k = 0; (k < MAXHH); k++)
1501 n_bound[k] = 0;
1503 /* Loop over possible donors */
1504 for (i = 0; (i < hb->d.nrd); i++)
1506 for (j = 0; (j < hb->d.nhydro[i]); j++)
1508 n_bound[hb->d.nhbonds[i][j]]++;
1511 fprintf(fp, "%12.5e", t);
1512 nbtot = 0;
1513 for (k = 0; (k < MAXHH); k++)
1515 fprintf(fp, " %8d", n_bound[k]);
1516 nbtot += n_bound[k]*k;
1518 fprintf(fp, " %8d\n", nbtot);
1521 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
1522 const gmx_output_env_t *oenv)
1524 FILE *fp;
1525 const char *leg[] = { "p(t)", "t p(t)" };
1526 int *histo;
1527 int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
1528 int nframes = hb->nframes;
1529 unsigned int **h;
1530 real t, x1, dt;
1531 double sum, integral;
1532 t_hbond *hbh;
1534 snew(h, hb->maxhydro);
1535 snew(histo, nframes+1);
1536 /* Total number of hbonds analyzed here */
1537 for (i = 0; (i < hb->d.nrd); i++)
1539 for (k = 0; (k < hb->a.nra); k++)
1541 hbh = hb->hbmap[i][k];
1542 if (hbh)
1544 if (bMerge)
1546 if (hbh->h[0])
1548 h[0] = hbh->h[0];
1549 nhydro = 1;
1551 else
1553 nhydro = 0;
1556 else
1558 nhydro = 0;
1559 for (m = 0; (m < hb->maxhydro); m++)
1561 if (hbh->h[m])
1563 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
1567 for (nh = 0; (nh < nhydro); nh++)
1569 ohb = 0;
1570 j0 = 0;
1572 for (j = 0; (j <= hbh->nframes); j++)
1574 ihb = static_cast<int>(is_hb(h[nh], j));
1575 if (debug && (ndump < 10))
1577 fprintf(debug, "%5d %5d\n", j, ihb);
1579 if (ihb != ohb)
1581 if (ihb)
1583 j0 = j;
1585 else
1587 histo[j-j0]++;
1589 ohb = ihb;
1592 ndump++;
1597 fprintf(stderr, "\n");
1598 if (bContact)
1600 fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
1602 else
1604 fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
1605 oenv);
1608 xvgr_legend(fp, asize(leg), leg, oenv);
1609 j0 = nframes-1;
1610 while ((j0 > 0) && (histo[j0] == 0))
1612 j0--;
1614 sum = 0;
1615 for (i = 0; (i <= j0); i++)
1617 sum += histo[i];
1619 dt = hb->time[1]-hb->time[0];
1620 sum = dt*sum;
1621 integral = 0;
1622 for (i = 1; (i <= j0); i++)
1624 t = hb->time[i] - hb->time[0] - 0.5*dt;
1625 x1 = t*histo[i]/sum;
1626 fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1);
1627 integral += x1;
1629 integral *= dt;
1630 xvgrclose(fp);
1631 printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
1632 printf("Note that the lifetime obtained in this manner is close to useless\n");
1633 printf("Use the -ac option instead and check the Forward lifetime\n");
1634 please_cite(stdout, "Spoel2006b");
1635 sfree(h);
1636 sfree(histo);
1639 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
1641 FILE *fp;
1642 int i, j, k, m, nd, ihb, idist;
1643 int nframes = hb->nframes;
1644 gmx_bool bPrint;
1645 t_hbond *hbh;
1647 if (nDump <= 0)
1649 return;
1651 fp = gmx_ffopen("debug-ac.xvg", "w");
1652 for (j = 0; (j < nframes); j++)
1654 fprintf(fp, "%10.3f", hb->time[j]);
1655 for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
1657 for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
1659 bPrint = FALSE;
1660 ihb = idist = 0;
1661 hbh = hb->hbmap[i][k];
1662 if (oneHB)
1664 if (hbh->h[0])
1666 ihb = static_cast<int>(is_hb(hbh->h[0], j));
1667 idist = static_cast<int>(is_hb(hbh->g[0], j));
1668 bPrint = TRUE;
1671 else
1673 for (m = 0; (m < hb->maxhydro) && !ihb; m++)
1675 ihb = static_cast<int>((ihb != 0) || (((hbh->h[m]) != nullptr) && is_hb(hbh->h[m], j)));
1676 idist = static_cast<int>((idist != 0) || (((hbh->g[m]) != nullptr) && is_hb(hbh->g[m], j)));
1678 /* This is not correct! */
1679 /* What isn't correct? -Erik M */
1680 bPrint = TRUE;
1682 if (bPrint)
1684 fprintf(fp, " %1d-%1d", ihb, idist);
1685 nd++;
1689 fprintf(fp, "\n");
1691 gmx_ffclose(fp);
1694 static real calc_dg(real tau, real temp)
1696 real kbt;
1698 kbt = BOLTZ*temp;
1699 if (tau <= 0)
1701 return -666;
1703 else
1705 return kbt*std::log(kbt*tau/PLANCK);
1709 typedef struct {
1710 int n0, n1, nparams, ndelta;
1711 real kkk[2];
1712 real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
1713 } t_luzar;
1715 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
1716 real kt[], real sigma_ct[], real sigma_nt[],
1717 real sigma_kt[], real *k, real *kp,
1718 real *sigma_k, real *sigma_kp,
1719 real fit_start)
1721 #define NK 10
1722 int i, j;
1723 t_luzar tl;
1724 real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
1726 *sigma_k = 0;
1727 *sigma_kp = 0;
1729 for (i = 0; (i < n); i++)
1731 if (t[i] >= fit_start)
1733 break;
1736 tl.n0 = i;
1737 tl.n1 = n;
1738 tl.nparams = 2;
1739 tl.ndelta = 1;
1740 tl.t = t;
1741 tl.ct = ct;
1742 tl.nt = nt;
1743 tl.kt = kt;
1744 tl.sigma_ct = sigma_ct;
1745 tl.sigma_nt = sigma_nt;
1746 tl.sigma_kt = sigma_kt;
1747 tl.kkk[0] = *k;
1748 tl.kkk[1] = *kp;
1750 chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
1751 *k = tl.kkk[0];
1752 *kp = tl.kkk[1] = *kp;
1753 tl.ndelta = NK;
1754 for (j = 0; (j < NK); j++)
1756 kkk += tl.kkk[0];
1757 kkp += tl.kkk[1];
1758 kk2 += gmx::square(tl.kkk[0]);
1759 kp2 += gmx::square(tl.kkk[1]);
1760 tl.n0++;
1762 *sigma_k = std::sqrt(kk2/NK - gmx::square(kkk/NK));
1763 *sigma_kp = std::sqrt(kp2/NK - gmx::square(kkp/NK));
1765 return chi2;
1768 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
1769 real sigma_ct[], real sigma_nt[], real sigma_kt[],
1770 real fit_start, real temp)
1772 int i0, i;
1773 real k = 1, kp = 1, kow = 1;
1774 real Q = 0, chi2, tau_hb, dtau, tau_rlx, e_1, sigma_k, sigma_kp, ddg;
1775 double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
1776 gmx_bool bError = (sigma_ct != nullptr) && (sigma_nt != nullptr) && (sigma_kt != nullptr);
1778 for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
1782 if (i0 < n-2)
1784 for (i = i0; (i < n); i++)
1786 sc2 += gmx::square(ct[i]);
1787 sn2 += gmx::square(nt[i]);
1788 sk2 += gmx::square(kt[i]);
1789 sck += ct[i]*kt[i];
1790 snk += nt[i]*kt[i];
1791 scn += ct[i]*nt[i];
1793 printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
1794 tmp = (sn2*sc2-gmx::square(scn));
1795 if ((tmp > 0) && (sn2 > 0))
1797 k = (sn2*sck-scn*snk)/tmp;
1798 kp = (k*scn-snk)/sn2;
1799 if (bError)
1801 chi2 = 0;
1802 for (i = i0; (i < n); i++)
1804 chi2 += gmx::square(k*ct[i]-kp*nt[i]-kt[i]);
1806 compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
1807 sigma_kt, &k, &kp,
1808 &sigma_k, &sigma_kp, fit_start);
1809 Q = 0; /* quality_of_fit(chi2, 2);*/
1810 ddg = BOLTZ*temp*sigma_k/k;
1811 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
1812 chi2, Q);
1813 printf("The Rate and Delta G are followed by an error estimate\n");
1814 printf("----------------------------------------------------------\n"
1815 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
1816 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1817 k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
1818 ddg = BOLTZ*temp*sigma_kp/kp;
1819 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1820 kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
1822 else
1824 chi2 = 0;
1825 for (i = i0; (i < n); i++)
1827 chi2 += gmx::square(k*ct[i]-kp*nt[i]-kt[i]);
1829 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
1830 chi2, Q);
1831 printf("--------------------------------------------------\n"
1832 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
1833 printf("Forward %10.3f %8.3f %10.3f %10g\n",
1834 k, 1/k, calc_dg(1/k, temp), chi2);
1835 printf("Backward %10.3f %8.3f %10.3f\n",
1836 kp, 1/kp, calc_dg(1/kp, temp));
1839 if (sc2 > 0)
1841 kow = 2*sck/sc2;
1842 printf("One-way %10.3f %s%8.3f %10.3f\n",
1843 kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp));
1845 else
1847 printf(" - Numerical problems computing HB thermodynamics:\n"
1848 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
1849 sc2, sn2, sk2, sck, snk, scn);
1851 /* Determine integral of the correlation function */
1852 tau_hb = evaluate_integral(n, t, ct, nullptr, (t[n-1]-t[0])/2, &dtau);
1853 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
1854 bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
1855 e_1 = std::exp(-1.0);
1856 for (i = 0; (i < n-2); i++)
1858 if ((ct[i] > e_1) && (ct[i+1] <= e_1))
1860 break;
1863 if (i < n-2)
1865 /* Determine tau_relax from linear interpolation */
1866 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
1867 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx,
1868 tau_rlx, bError ? " " : "",
1869 calc_dg(tau_rlx, temp));
1872 else
1874 printf("Correlation functions too short to compute thermodynamics\n");
1878 void compute_derivative(int nn, const real x[], const real y[], real dydx[])
1880 int j;
1882 /* Compute k(t) = dc(t)/dt */
1883 for (j = 1; (j < nn-1); j++)
1885 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
1887 /* Extrapolate endpoints */
1888 dydx[0] = 2*dydx[1] - dydx[2];
1889 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
1892 static void normalizeACF(real *ct, real *gt, int nhb, int len)
1894 real ct_fac, gt_fac = 0;
1895 int i;
1897 /* Xu and Berne use the same normalization constant */
1899 ct_fac = 1.0/ct[0];
1900 if (nhb != 0)
1902 gt_fac = 1.0/nhb;
1905 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
1906 for (i = 0; i < len; i++)
1908 ct[i] *= ct_fac;
1909 if (gt != nullptr)
1911 gt[i] *= gt_fac;
1916 static void do_hbac(const char *fn, t_hbdata *hb,
1917 int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
1918 real temp, gmx_bool R2, const gmx_output_env_t *oenv,
1919 int nThreads)
1921 FILE *fp;
1922 int i, j, k, m, ihb, idist, n2, nn;
1924 const char *legLuzar[] = {
1925 "Ac\\sfin sys\\v{}\\z{}(t)",
1926 "Ac(t)",
1927 "Cc\\scontact,hb\\v{}\\z{}(t)",
1928 "-dAc\\sfs\\v{}\\z{}/dt"
1930 gmx_bool bNorm = FALSE;
1931 double nhb = 0;
1932 real *rhbex = nullptr, *ht, *gt, *ght, *dght, *kt;
1933 real *ct, tail, tail2, dtail, *cct;
1934 const real tol = 1e-3;
1935 int nframes = hb->nframes;
1936 unsigned int **h = nullptr, **g = nullptr;
1937 int nh, nhbonds, nhydro;
1938 t_hbond *hbh;
1939 int acType;
1940 int *dondata = nullptr;
1942 enum {
1943 AC_NONE, AC_NN, AC_GEM, AC_LUZAR
1946 const bool bOMP = GMX_OPENMP;
1948 printf("Doing autocorrelation ");
1950 acType = AC_LUZAR;
1951 printf("according to the theory of Luzar and Chandler.\n");
1952 fflush(stdout);
1954 /* build hbexist matrix in reals for autocorr */
1955 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
1956 n2 = 1;
1957 while (n2 < nframes)
1959 n2 *= 2;
1962 nn = nframes/2;
1964 if (acType != AC_NN || bOMP)
1966 snew(h, hb->maxhydro);
1967 snew(g, hb->maxhydro);
1970 /* Dump hbonds for debugging */
1971 dump_ac(hb, bMerge || bContact, nDump);
1973 /* Total number of hbonds analyzed here */
1974 nhbonds = 0;
1976 if (acType != AC_LUZAR && bOMP)
1978 nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
1980 gmx_omp_set_num_threads(nThreads);
1981 snew(dondata, nThreads);
1982 for (i = 0; i < nThreads; i++)
1984 dondata[i] = -1;
1986 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
1987 "Expect close to linear scaling over this donor-loop.\n", nThreads);
1988 fflush(stdout);
1989 fprintf(stderr, "Donors: [thread no]\n");
1991 char tmpstr[7];
1992 for (i = 0; i < nThreads; i++)
1994 snprintf(tmpstr, 7, "[%i]", i);
1995 fprintf(stderr, "%-7s", tmpstr);
1998 fprintf(stderr, "\n");
2002 /* Build the ACF */
2003 snew(rhbex, 2*n2);
2004 snew(ct, 2*n2);
2005 snew(gt, 2*n2);
2006 snew(ht, 2*n2);
2007 snew(ght, 2*n2);
2008 snew(dght, 2*n2);
2010 snew(kt, nn);
2011 snew(cct, nn);
2013 for (i = 0; (i < hb->d.nrd); i++)
2015 for (k = 0; (k < hb->a.nra); k++)
2017 nhydro = 0;
2018 hbh = hb->hbmap[i][k];
2020 if (hbh)
2022 if (bMerge || bContact)
2024 if (ISHB(hbh->history[0]))
2026 h[0] = hbh->h[0];
2027 g[0] = hbh->g[0];
2028 nhydro = 1;
2031 else
2033 for (m = 0; (m < hb->maxhydro); m++)
2035 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
2037 g[nhydro] = hbh->g[m];
2038 h[nhydro] = hbh->h[m];
2039 nhydro++;
2044 int nf = hbh->nframes;
2045 for (nh = 0; (nh < nhydro); nh++)
2047 int nrint = bContact ? hb->nrdist : hb->nrhb;
2048 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
2050 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
2051 fflush(stderr);
2053 nhbonds++;
2054 for (j = 0; (j < nframes); j++)
2056 if (j <= nf)
2058 ihb = static_cast<int>(is_hb(h[nh], j));
2059 idist = static_cast<int>(is_hb(g[nh], j));
2061 else
2063 ihb = idist = 0;
2065 rhbex[j] = ihb;
2066 /* For contacts: if a second cut-off is provided, use it,
2067 * otherwise use g(t) = 1-h(t) */
2068 if (!R2 && bContact)
2070 gt[j] = 1-ihb;
2072 else
2074 gt[j] = idist*(1-ihb);
2076 ht[j] = rhbex[j];
2077 nhb += ihb;
2080 /* The autocorrelation function is normalized after summation only */
2081 low_do_autocorr(nullptr, oenv, nullptr, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
2082 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
2084 /* Cross correlation analysis for thermodynamics */
2085 for (j = nframes; (j < n2); j++)
2087 ht[j] = 0;
2088 gt[j] = 0;
2091 cross_corr(n2, ht, gt, dght);
2093 for (j = 0; (j < nn); j++)
2095 ct[j] += rhbex[j];
2096 ght[j] += dght[j];
2102 fprintf(stderr, "\n");
2103 sfree(h);
2104 sfree(g);
2105 normalizeACF(ct, ght, static_cast<int>(nhb), nn);
2107 /* Determine tail value for statistics */
2108 tail = 0;
2109 tail2 = 0;
2110 for (j = nn/2; (j < nn); j++)
2112 tail += ct[j];
2113 tail2 += ct[j]*ct[j];
2115 tail /= (nn - int{nn/2});
2116 tail2 /= (nn - int{nn/2});
2117 dtail = std::sqrt(tail2-tail*tail);
2119 /* Check whether the ACF is long enough */
2120 if (dtail > tol)
2122 printf("\nWARNING: Correlation function is probably not long enough\n"
2123 "because the standard deviation in the tail of C(t) > %g\n"
2124 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2125 tol, tail, dtail);
2127 for (j = 0; (j < nn); j++)
2129 cct[j] = ct[j];
2130 ct[j] = (cct[j]-tail)/(1-tail);
2132 /* Compute negative derivative k(t) = -dc(t)/dt */
2133 compute_derivative(nn, hb->time, ct, kt);
2134 for (j = 0; (j < nn); j++)
2136 kt[j] = -kt[j];
2140 if (bContact)
2142 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
2144 else
2146 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
2148 xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
2151 for (j = 0; (j < nn); j++)
2153 fprintf(fp, "%10g %10g %10g %10g %10g\n",
2154 hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
2156 xvgrclose(fp);
2158 analyse_corr(nn, hb->time, ct, ght, kt, nullptr, nullptr, nullptr,
2159 fit_start, temp);
2161 do_view(oenv, fn, nullptr);
2162 sfree(rhbex);
2163 sfree(ct);
2164 sfree(gt);
2165 sfree(ht);
2166 sfree(ght);
2167 sfree(dght);
2168 sfree(cct);
2169 sfree(kt);
2172 static void init_hbframe(t_hbdata *hb, int nframes, real t)
2174 int i;
2176 hb->time[nframes] = t;
2177 hb->nhb[nframes] = 0;
2178 hb->ndist[nframes] = 0;
2179 for (i = 0; (i < max_hx); i++)
2181 hb->nhx[nframes][i] = 0;
2185 static FILE *open_donor_properties_file(const char *fn,
2186 t_hbdata *hb,
2187 const gmx_output_env_t *oenv)
2189 FILE *fp = nullptr;
2190 const char *leg[] = { "Nbound", "Nfree" };
2192 if (!fn || !hb)
2194 return nullptr;
2197 fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
2198 xvgr_legend(fp, asize(leg), leg, oenv);
2200 return fp;
2203 static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
2205 int i, j, k, nbound, nb, nhtot;
2207 if (!fp || !hb)
2209 return;
2211 nbound = 0;
2212 nhtot = 0;
2213 for (i = 0; (i < hb->d.nrd); i++)
2215 for (k = 0; (k < hb->d.nhydro[i]); k++)
2217 nb = 0;
2218 nhtot++;
2219 for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
2221 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
2222 is_hb(hb->hbmap[i][j]->h[k], nframes))
2224 nb = 1;
2227 nbound += nb;
2230 fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound);
2233 static void dump_hbmap(t_hbdata *hb,
2234 int nfile, t_filenm fnm[], gmx_bool bTwo,
2235 gmx_bool bContact, const int isize[], int *index[], char *grpnames[],
2236 const t_atoms *atoms)
2238 FILE *fp, *fplog;
2239 int ddd, hhh, aaa, i, j, k, m, grp;
2240 char ds[32], hs[32], as[32];
2241 gmx_bool first;
2243 fp = opt2FILE("-hbn", nfile, fnm, "w");
2244 if (opt2bSet("-g", nfile, fnm))
2246 fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
2247 fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
2249 else
2251 fplog = nullptr;
2253 for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
2255 fprintf(fp, "[ %s ]", grpnames[grp]);
2256 for (i = 0; i < isize[grp]; i++)
2258 fprintf(fp, (i%15) ? " " : "\n");
2259 fprintf(fp, " %4d", index[grp][i]+1);
2261 fprintf(fp, "\n");
2263 if (!bContact)
2265 fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
2266 for (i = 0; (i < hb->d.nrd); i++)
2268 if (hb->d.grp[i] == grp)
2270 for (j = 0; (j < hb->d.nhydro[i]); j++)
2272 fprintf(fp, " %4d %4d", hb->d.don[i]+1,
2273 hb->d.hydro[i][j]+1);
2275 fprintf(fp, "\n");
2278 first = TRUE;
2279 fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
2280 for (i = 0; (i < hb->a.nra); i++)
2282 if (hb->a.grp[i] == grp)
2284 fprintf(fp, (i%15 && !first) ? " " : "\n");
2285 fprintf(fp, " %4d", hb->a.acc[i]+1);
2286 first = FALSE;
2289 fprintf(fp, "\n");
2292 if (bTwo)
2294 fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
2295 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
2297 else
2299 fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
2302 for (i = 0; (i < hb->d.nrd); i++)
2304 ddd = hb->d.don[i];
2305 for (k = 0; (k < hb->a.nra); k++)
2307 aaa = hb->a.acc[k];
2308 for (m = 0; (m < hb->d.nhydro[i]); m++)
2310 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
2312 sprintf(ds, "%s", mkatomname(atoms, ddd));
2313 sprintf(as, "%s", mkatomname(atoms, aaa));
2314 if (bContact)
2316 fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
2317 if (fplog)
2319 fprintf(fplog, "%12s %12s\n", ds, as);
2322 else
2324 hhh = hb->d.hydro[i][m];
2325 sprintf(hs, "%s", mkatomname(atoms, hhh));
2326 fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
2327 if (fplog)
2329 fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
2336 gmx_ffclose(fp);
2337 if (fplog)
2339 gmx_ffclose(fplog);
2343 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
2344 * It mimics add_frames() and init_frame() to some extent. */
2345 static void sync_hbdata(t_hbdata *p_hb, int nframes)
2347 if (nframes >= p_hb->max_frames)
2349 p_hb->max_frames += 4096;
2350 srenew(p_hb->nhb, p_hb->max_frames);
2351 srenew(p_hb->ndist, p_hb->max_frames);
2352 srenew(p_hb->n_bound, p_hb->max_frames);
2353 srenew(p_hb->nhx, p_hb->max_frames);
2354 if (p_hb->bDAnr)
2356 srenew(p_hb->danr, p_hb->max_frames);
2358 std::memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
2359 std::memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
2360 p_hb->nhb[nframes] = 0;
2361 p_hb->ndist[nframes] = 0;
2364 p_hb->nframes = nframes;
2366 std::memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
2369 int gmx_hbond(int argc, char *argv[])
2371 const char *desc[] = {
2372 "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
2373 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
2374 "(zero is extended) and the distance Donor - Acceptor",
2375 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
2376 "OH and NH groups are regarded as donors, O is an acceptor always,",
2377 "N is an acceptor by default, but this can be switched using",
2378 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
2379 "to the first preceding non-hydrogen atom.[PAR]",
2381 "You need to specify two groups for analysis, which must be either",
2382 "identical or non-overlapping. All hydrogen bonds between the two",
2383 "groups are analyzed.[PAR]",
2385 "If you set [TT]-shell[tt], you will be asked for an additional index group",
2386 "which should contain exactly one atom. In this case, only hydrogen",
2387 "bonds between atoms within the shell distance from the one atom are",
2388 "considered.[PAR]",
2390 "With option -ac, rate constants for hydrogen bonding can be derived with the",
2391 "model of Luzar and Chandler (Nature 379:55, 1996; J. Chem. Phys. 113:23, 2000).",
2392 "If contact kinetics are analyzed by using the -contact option, then",
2393 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
2394 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
2395 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
2396 "See mentioned literature for more details and definitions.",
2397 "[PAR]",
2399 /* "It is also possible to analyse specific hydrogen bonds with",
2400 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
2401 "Donor Hydrogen Acceptor, in the following way::",
2403 "[ selected ]",
2404 " 20 21 24",
2405 " 25 26 29",
2406 " 1 3 6",
2408 "Note that the triplets need not be on separate lines.",
2409 "Each atom triplet specifies a hydrogen bond to be analyzed,",
2410 "note also that no check is made for the types of atoms.[PAR]",
2413 "[BB]Output:[bb]",
2415 " * [TT]-num[tt]: number of hydrogen bonds as a function of time.",
2416 " * [TT]-ac[tt]: average over all autocorrelations of the existence",
2417 " functions (either 0 or 1) of all hydrogen bonds.",
2418 " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
2419 " * [TT]-ang[tt]: angle distribution of all hydrogen bonds.",
2420 " * [TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
2421 " where n and n+i stand for residue numbers and i ranges from 0 to 6.",
2422 " This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
2423 " with helices in proteins.",
2424 " * [TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
2425 " for selected groups, all hydrogen bonded atoms from all groups and",
2426 " all solvent atoms involved in insertion.",
2427 " * [TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
2428 " frames, this also contains information on solvent insertion",
2429 " into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
2430 " index file.",
2431 " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
2432 " each timeframe. This is especially useful when using [TT]-shell[tt].",
2433 " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
2434 " compare results to Raman Spectroscopy.",
2436 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
2437 "require an amount of memory proportional to the total numbers of donors",
2438 "times the total number of acceptors in the selected group(s)."
2441 static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
2442 static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15;
2443 static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
2444 static int nDump = 0;
2445 static int nThreads = 0;
2447 static gmx_bool bContact = FALSE;
2449 /* options */
2450 t_pargs pa [] = {
2451 { "-a", FALSE, etREAL, {&acut},
2452 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
2453 { "-r", FALSE, etREAL, {&rcut},
2454 "Cutoff radius (nm, X - Acceptor, see next option)" },
2455 { "-da", FALSE, etBOOL, {&bDA},
2456 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
2457 { "-r2", FALSE, etREAL, {&r2cut},
2458 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
2459 { "-abin", FALSE, etREAL, {&abin},
2460 "Binwidth angle distribution (degrees)" },
2461 { "-rbin", FALSE, etREAL, {&rbin},
2462 "Binwidth distance distribution (nm)" },
2463 { "-nitacc", FALSE, etBOOL, {&bNitAcc},
2464 "Regard nitrogen atoms as acceptors" },
2465 { "-contact", FALSE, etBOOL, {&bContact},
2466 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
2467 { "-shell", FALSE, etREAL, {&rshell},
2468 "when > 0, only calculate hydrogen bonds within # nm shell around "
2469 "one particle" },
2470 { "-fitstart", FALSE, etREAL, {&fit_start},
2471 "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]" },
2472 { "-fitend", FALSE, etREAL, {&fit_end},
2473 "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])" },
2474 { "-temp", FALSE, etREAL, {&temp},
2475 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
2476 { "-dump", FALSE, etINT, {&nDump},
2477 "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
2478 { "-max_hb", FALSE, etREAL, {&maxnhb},
2479 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
2480 { "-merge", FALSE, etBOOL, {&bMerge},
2481 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
2482 #if GMX_OPENMP
2483 { "-nthreads", FALSE, etINT, {&nThreads},
2484 "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)"},
2485 #endif
2487 const char *bugs[] = {
2488 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
2490 t_filenm fnm[] = {
2491 { efTRX, "-f", nullptr, ffREAD },
2492 { efTPR, nullptr, nullptr, ffREAD },
2493 { efNDX, nullptr, nullptr, ffOPTRD },
2494 /* { efNDX, "-sel", "select", ffOPTRD },*/
2495 { efXVG, "-num", "hbnum", ffWRITE },
2496 { efLOG, "-g", "hbond", ffOPTWR },
2497 { efXVG, "-ac", "hbac", ffOPTWR },
2498 { efXVG, "-dist", "hbdist", ffOPTWR },
2499 { efXVG, "-ang", "hbang", ffOPTWR },
2500 { efXVG, "-hx", "hbhelix", ffOPTWR },
2501 { efNDX, "-hbn", "hbond", ffOPTWR },
2502 { efXPM, "-hbm", "hbmap", ffOPTWR },
2503 { efXVG, "-don", "donor", ffOPTWR },
2504 { efXVG, "-dan", "danum", ffOPTWR },
2505 { efXVG, "-life", "hblife", ffOPTWR },
2506 { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
2509 #define NFILE asize(fnm)
2511 char hbmap [HB_NR] = { ' ', 'o', '-', '*' };
2512 const char *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
2513 t_rgb hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
2515 t_trxstatus *status;
2516 bool trrStatus = true;
2517 t_topology top;
2518 t_pargs *ppa;
2519 int npargs, natoms, nframes = 0, shatom;
2520 int *isize;
2521 char **grpnames;
2522 int **index;
2523 rvec *x, hbox;
2524 matrix box;
2525 real t, ccut, dist = 0.0, ang = 0.0;
2526 double max_nhb, aver_nhb, aver_dist;
2527 int h = 0, i = 0, j, k = 0, ogrp, nsel;
2528 int xi, yi, zi, ai;
2529 int xj, yj, zj, aj, xjj, yjj, zjj;
2530 gmx_bool bSelected, bHBmap, bStop, bTwo, bBox, bTric;
2531 int *adist, *rdist;
2532 int grp, nabin, nrbin, resdist, ihb;
2533 char **leg;
2534 t_hbdata *hb;
2535 FILE *fp, *fpnhb = nullptr, *donor_properties = nullptr;
2536 t_gridcell ***grid;
2537 t_ncell *icell, *jcell;
2538 ivec ngrid;
2539 unsigned char *datable;
2540 gmx_output_env_t *oenv;
2541 int ii, hh, actual_nThreads;
2542 int threadNr = 0;
2543 gmx_bool bParallel;
2544 gmx_bool bEdge_yjj, bEdge_xjj;
2546 t_hbdata **p_hb = nullptr; /* one per thread, then merge after the frame loop */
2547 int **p_adist = nullptr, **p_rdist = nullptr; /* a histogram for each thread. */
2549 const bool bOMP = GMX_OPENMP;
2551 npargs = asize(pa);
2552 ppa = add_acf_pargs(&npargs, pa);
2554 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
2555 ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
2557 sfree(ppa);
2558 return 0;
2561 /* process input */
2562 bSelected = FALSE;
2563 ccut = std::cos(acut*DEG2RAD);
2565 if (bContact)
2567 if (bSelected)
2569 gmx_fatal(FARGS, "Can not analyze selected contacts.");
2571 if (!bDA)
2573 gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
2577 /* Initiate main data structure! */
2578 bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
2579 opt2bSet("-life", NFILE, fnm) ||
2580 opt2bSet("-hbn", NFILE, fnm) ||
2581 opt2bSet("-hbm", NFILE, fnm));
2583 if (opt2bSet("-nhbdist", NFILE, fnm))
2585 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
2586 fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
2587 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
2588 xvgr_legend(fpnhb, asize(leg), leg, oenv);
2591 hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact);
2593 /* get topology */
2594 t_inputrec irInstance;
2595 t_inputrec *ir = &irInstance;
2596 read_tpx_top(ftp2fn(efTPR, NFILE, fnm), ir, box, &natoms, nullptr, nullptr, &top);
2598 snew(grpnames, grNR);
2599 snew(index, grNR);
2600 snew(isize, grNR);
2601 /* Make Donor-Acceptor table */
2602 snew(datable, top.atoms.nr);
2604 if (bSelected)
2606 /* analyze selected hydrogen bonds */
2607 printf("Select group with selected atoms:\n");
2608 get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
2609 1, &nsel, index, grpnames);
2610 if (nsel % 3)
2612 gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
2613 "and therefore cannot contain triplets of "
2614 "Donor-Hydrogen-Acceptor", grpnames[0]);
2616 bTwo = FALSE;
2618 for (i = 0; (i < nsel); i += 3)
2620 int dd = index[0][i];
2621 int aa = index[0][i+2];
2622 /* int */ hh = index[0][i+1];
2623 add_dh (&hb->d, dd, hh, i, datable);
2624 add_acc(&hb->a, aa, i);
2625 /* Should this be here ? */
2626 snew(hb->d.dptr, top.atoms.nr);
2627 snew(hb->a.aptr, top.atoms.nr);
2628 add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact);
2630 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
2631 isize[0], grpnames[0]);
2633 else
2635 /* analyze all hydrogen bonds: get group(s) */
2636 printf("Specify 2 groups to analyze:\n");
2637 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
2638 2, isize, index, grpnames);
2640 /* check if we have two identical or two non-overlapping groups */
2641 bTwo = isize[0] != isize[1];
2642 for (i = 0; (i < isize[0]) && !bTwo; i++)
2644 bTwo = index[0][i] != index[1][i];
2646 if (bTwo)
2648 printf("Checking for overlap in atoms between %s and %s\n",
2649 grpnames[0], grpnames[1]);
2651 gen_datable(index[0], isize[0], datable, top.atoms.nr);
2653 for (i = 0; i < isize[1]; i++)
2655 if (ISINGRP(datable[index[1][i]]))
2657 gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
2658 grpnames[0], grpnames[1]);
2662 if (bTwo)
2664 printf("Calculating %s "
2665 "between %s (%d atoms) and %s (%d atoms)\n",
2666 bContact ? "contacts" : "hydrogen bonds",
2667 grpnames[0], isize[0], grpnames[1], isize[1]);
2669 else
2671 fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
2672 bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
2675 sfree(datable);
2677 /* search donors and acceptors in groups */
2678 snew(datable, top.atoms.nr);
2679 for (i = 0; (i < grNR); i++)
2681 if ( ((i == gr0) && !bSelected ) ||
2682 ((i == gr1) && bTwo ))
2684 gen_datable(index[i], isize[i], datable, top.atoms.nr);
2685 if (bContact)
2687 search_acceptors(&top, isize[i], index[i], &hb->a, i,
2688 bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
2689 search_donors (&top, isize[i], index[i], &hb->d, i,
2690 TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
2692 else
2694 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
2695 search_donors (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
2697 if (bTwo)
2699 clear_datable_grp(datable, top.atoms.nr);
2703 sfree(datable);
2704 printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
2705 /*if (bSelected)
2706 snew(donors[gr0D], dons[gr0D].nrd);*/
2708 donor_properties = open_donor_properties_file(opt2fn_null("-don", NFILE, fnm), hb, oenv);
2710 if (bHBmap)
2712 printf("Making hbmap structure...");
2713 /* Generate hbond data structure */
2714 mk_hbmap(hb);
2715 printf("done.\n");
2718 /* check input */
2719 bStop = FALSE;
2720 if (hb->d.nrd + hb->a.nra == 0)
2722 printf("No Donors or Acceptors found\n");
2723 bStop = TRUE;
2725 if (!bStop)
2727 if (hb->d.nrd == 0)
2729 printf("No Donors found\n");
2730 bStop = TRUE;
2732 if (hb->a.nra == 0)
2734 printf("No Acceptors found\n");
2735 bStop = TRUE;
2738 if (bStop)
2740 gmx_fatal(FARGS, "Nothing to be done");
2743 shatom = 0;
2744 if (rshell > 0)
2746 int shisz;
2747 int *shidx;
2748 char *shgrpnm;
2749 /* get index group with atom for shell */
2752 printf("Select atom for shell (1 atom):\n");
2753 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
2754 1, &shisz, &shidx, &shgrpnm);
2755 if (shisz != 1)
2757 printf("group contains %d atoms, should be 1 (one)\n", shisz);
2760 while (shisz != 1);
2761 shatom = shidx[0];
2762 printf("Will calculate hydrogen bonds within a shell "
2763 "of %g nm around atom %i\n", rshell, shatom+1);
2766 /* Analyze trajectory */
2767 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
2768 if (natoms > top.atoms.nr)
2770 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
2771 top.atoms.nr, natoms);
2774 bBox = (ir->ePBC != epbcNONE);
2775 grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
2776 nabin = static_cast<int>(acut/abin);
2777 nrbin = static_cast<int>(rcut/rbin);
2778 snew(adist, nabin+1);
2779 snew(rdist, nrbin+1);
2781 #if !GMX_OPENMP
2782 #define __ADIST adist
2783 #define __RDIST rdist
2784 #define __HBDATA hb
2785 #else /* GMX_OPENMP ================================================== \
2786 * Set up the OpenMP stuff, |
2787 * like the number of threads and such |
2788 * Also start the parallel loop. |
2790 #define __ADIST p_adist[threadNr]
2791 #define __RDIST p_rdist[threadNr]
2792 #define __HBDATA p_hb[threadNr]
2793 #endif
2794 if (bOMP)
2796 bParallel = !bSelected;
2798 if (bParallel)
2800 actual_nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2802 gmx_omp_set_num_threads(actual_nThreads);
2803 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
2804 fflush(stdout);
2806 else
2808 actual_nThreads = 1;
2811 snew(p_hb, actual_nThreads);
2812 snew(p_adist, actual_nThreads);
2813 snew(p_rdist, actual_nThreads);
2814 for (i = 0; i < actual_nThreads; i++)
2816 snew(p_hb[i], 1);
2817 snew(p_adist[i], nabin+1);
2818 snew(p_rdist[i], nrbin+1);
2820 p_hb[i]->max_frames = 0;
2821 p_hb[i]->nhb = nullptr;
2822 p_hb[i]->ndist = nullptr;
2823 p_hb[i]->n_bound = nullptr;
2824 p_hb[i]->time = nullptr;
2825 p_hb[i]->nhx = nullptr;
2827 p_hb[i]->bHBmap = hb->bHBmap;
2828 p_hb[i]->bDAnr = hb->bDAnr;
2829 p_hb[i]->wordlen = hb->wordlen;
2830 p_hb[i]->nframes = hb->nframes;
2831 p_hb[i]->maxhydro = hb->maxhydro;
2832 p_hb[i]->danr = hb->danr;
2833 p_hb[i]->d = hb->d;
2834 p_hb[i]->a = hb->a;
2835 p_hb[i]->hbmap = hb->hbmap;
2836 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
2838 p_hb[i]->nrhb = 0;
2839 p_hb[i]->nrdist = 0;
2843 /* Make a thread pool here,
2844 * instead of forking anew at every frame. */
2846 #pragma omp parallel \
2847 firstprivate(i) \
2848 private(j, h, ii, hh, \
2849 xi, yi, zi, xj, yj, zj, threadNr, \
2850 dist, ang, icell, jcell, \
2851 grp, ogrp, ai, aj, xjj, yjj, zjj, \
2852 ihb, resdist, \
2853 k, bTric, \
2854 bEdge_xjj, bEdge_yjj) \
2855 default(shared)
2856 { /* Start of parallel region */
2857 if (bOMP)
2859 threadNr = gmx_omp_get_thread_num();
2864 bTric = bBox && TRICLINIC(box);
2866 if (bOMP)
2870 sync_hbdata(p_hb[threadNr], nframes);
2872 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2874 #pragma omp single
2878 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
2879 rshell, ngrid, grid);
2880 reset_nhbonds(&(hb->d));
2882 if (debug && bDebug)
2884 dump_grid(debug, ngrid, grid);
2887 add_frames(hb, nframes);
2888 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
2890 if (hb->bDAnr)
2892 count_da_grid(ngrid, grid, hb->danr[nframes]);
2895 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2896 } /* omp single */
2898 if (bOMP)
2900 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
2903 if (bSelected)
2906 #pragma omp single
2910 /* Do not parallelize this just yet. */
2911 /* int ii; */
2912 for (ii = 0; (ii < nsel); ii++)
2914 int dd = index[0][i];
2915 int aa = index[0][i+2];
2916 /* int */ hh = index[0][i+1];
2917 ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
2918 hbox, &dist, &ang, bDA, &h, bContact, bMerge);
2920 if (ihb)
2922 /* add to index if not already there */
2923 /* Add a hbond */
2924 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact);
2928 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2929 } /* omp single */
2930 } /* if (bSelected) */
2931 else
2933 /* The outer grid loop will have to do for now. */
2934 #pragma omp for schedule(dynamic)
2935 for (xi = 0; xi < ngrid[XX]; xi++)
2939 for (yi = 0; (yi < ngrid[YY]); yi++)
2941 for (zi = 0; (zi < ngrid[ZZ]); zi++)
2944 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
2945 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
2947 icell = &(grid[zi][yi][xi].d[grp]);
2949 if (bTwo)
2951 ogrp = 1-grp;
2953 else
2955 ogrp = grp;
2958 /* loop over all hydrogen atoms from group (grp)
2959 * in this gridcell (icell)
2961 for (ai = 0; (ai < icell->nr); ai++)
2963 i = icell->atoms[ai];
2965 /* loop over all adjacent gridcells (xj,yj,zj) */
2966 for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
2967 zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
2968 zjj++)
2970 zj = grid_mod(zjj, ngrid[ZZ]);
2971 bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
2972 for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
2973 yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
2974 yjj++)
2976 yj = grid_mod(yjj, ngrid[YY]);
2977 bEdge_xjj =
2978 (yj == 0) || (yj == ngrid[YY] - 1) ||
2979 (zj == 0) || (zj == ngrid[ZZ] - 1);
2980 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
2981 xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
2982 xjj++)
2984 xj = grid_mod(xjj, ngrid[XX]);
2985 jcell = &(grid[zj][yj][xj].a[ogrp]);
2986 /* loop over acceptor atoms from other group (ogrp)
2987 * in this adjacent gridcell (jcell)
2989 for (aj = 0; (aj < jcell->nr); aj++)
2991 j = jcell->atoms[aj];
2993 /* check if this once was a h-bond */
2994 ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
2995 hbox, &dist, &ang, bDA, &h, bContact, bMerge);
2997 if (ihb)
2999 /* add to index if not already there */
3000 /* Add a hbond */
3001 add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact);
3003 /* make angle and distance distributions */
3004 if (ihb == hbHB && !bContact)
3006 if (dist > rcut)
3008 gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
3010 ang *= RAD2DEG;
3011 __ADIST[static_cast<int>( ang/abin)]++;
3012 __RDIST[static_cast<int>(dist/rbin)]++;
3013 if (!bTwo)
3015 if (donor_index(&hb->d, grp, i) == NOTSET)
3017 gmx_fatal(FARGS, "Invalid donor %d", i);
3019 if (acceptor_index(&hb->a, ogrp, j) == NOTSET)
3021 gmx_fatal(FARGS, "Invalid acceptor %d", j);
3023 resdist = std::abs(top.atoms.atom[i].resind-top.atoms.atom[j].resind);
3024 if (resdist >= max_hx)
3026 resdist = max_hx-1;
3028 __HBDATA->nhx[nframes][resdist]++;
3033 } /* for aj */
3034 } /* for xjj */
3035 } /* for yjj */
3036 } /* for zjj */
3037 } /* for ai */
3038 } /* for grp */
3039 } /* for xi,yi,zi */
3042 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3044 } /* if (bSelected) {...} else */
3047 /* Better wait for all threads to finnish using x[] before updating it. */
3048 k = nframes;
3049 #pragma omp barrier
3050 #pragma omp critical
3054 /* Sum up histograms and counts from p_hb[] into hb */
3055 if (bOMP)
3057 hb->nhb[k] += p_hb[threadNr]->nhb[k];
3058 hb->ndist[k] += p_hb[threadNr]->ndist[k];
3059 for (j = 0; j < max_hx; j++)
3061 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
3065 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3068 /* Here are a handful of single constructs
3069 * to share the workload a bit. The most
3070 * important one is of course the last one,
3071 * where there's a potential bottleneck in form
3072 * of slow I/O. */
3073 #pragma omp barrier
3074 #pragma omp single
3078 analyse_donor_properties(donor_properties, hb, k, t);
3080 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3083 #pragma omp single
3087 if (fpnhb)
3089 do_nhb_dist(fpnhb, hb, t);
3092 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3095 #pragma omp single
3099 trrStatus = (read_next_x(oenv, status, &t, x, box));
3100 nframes++;
3102 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3105 #pragma omp barrier
3107 while (trrStatus);
3109 if (bOMP)
3111 #pragma omp critical
3113 hb->nrhb += p_hb[threadNr]->nrhb;
3114 hb->nrdist += p_hb[threadNr]->nrdist;
3117 /* Free parallel datastructures */
3118 sfree(p_hb[threadNr]->nhb);
3119 sfree(p_hb[threadNr]->ndist);
3120 sfree(p_hb[threadNr]->nhx);
3122 #pragma omp for
3123 for (i = 0; i < nabin; i++)
3127 for (j = 0; j < actual_nThreads; j++)
3130 adist[i] += p_adist[j][i];
3133 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3135 #pragma omp for
3136 for (i = 0; i <= nrbin; i++)
3140 for (j = 0; j < actual_nThreads; j++)
3142 rdist[i] += p_rdist[j][i];
3145 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
3148 sfree(p_adist[threadNr]);
3149 sfree(p_rdist[threadNr]);
3151 } /* End of parallel region */
3152 if (bOMP)
3154 sfree(p_adist);
3155 sfree(p_rdist);
3158 if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
3160 gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
3163 free_grid(ngrid, &grid);
3165 close_trx(status);
3167 if (donor_properties)
3169 xvgrclose(donor_properties);
3172 if (fpnhb)
3174 xvgrclose(fpnhb);
3177 /* Compute maximum possible number of different hbonds */
3178 if (maxnhb > 0)
3180 max_nhb = maxnhb;
3182 else
3184 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
3187 if (bHBmap)
3189 if (hb->nrhb == 0)
3191 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
3193 else
3195 printf("Found %d different %s in trajectory\n"
3196 "Found %d different atom-pairs within %s distance\n",
3197 hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
3198 hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
3200 if (bMerge)
3202 merge_hb(hb, bTwo, bContact);
3205 if (opt2bSet("-hbn", NFILE, fnm))
3207 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
3210 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3211 * to make the -hbn and -hmb output match eachother.
3212 * - Erik Marklund, May 30, 2006 */
3215 /* Print out number of hbonds and distances */
3216 aver_nhb = 0;
3217 aver_dist = 0;
3218 fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
3219 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3220 snew(leg, 2);
3221 snew(leg[0], STRLEN);
3222 snew(leg[1], STRLEN);
3223 sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
3224 sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
3225 xvgr_legend(fp, 2, leg, oenv);
3226 sfree(leg[1]);
3227 sfree(leg[0]);
3228 sfree(leg);
3229 for (i = 0; (i < nframes); i++)
3231 fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
3232 aver_nhb += hb->nhb[i];
3233 aver_dist += hb->ndist[i];
3235 xvgrclose(fp);
3236 aver_nhb /= nframes;
3237 aver_dist /= nframes;
3238 /* Print HB distance distribution */
3239 if (opt2bSet("-dist", NFILE, fnm))
3241 int sum;
3243 sum = 0;
3244 for (i = 0; i < nrbin; i++)
3246 sum += rdist[i];
3249 fp = xvgropen(opt2fn("-dist", NFILE, fnm),
3250 "Hydrogen Bond Distribution",
3251 bDA ?
3252 "Donor - Acceptor Distance (nm)" :
3253 "Hydrogen - Acceptor Distance (nm)", "", oenv);
3254 for (i = 0; i < nrbin; i++)
3256 fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*sum));
3258 xvgrclose(fp);
3261 /* Print HB angle distribution */
3262 if (opt2bSet("-ang", NFILE, fnm))
3264 long sum;
3266 sum = 0;
3267 for (i = 0; i < nabin; i++)
3269 sum += adist[i];
3272 fp = xvgropen(opt2fn("-ang", NFILE, fnm),
3273 "Hydrogen Bond Distribution",
3274 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
3275 for (i = 0; i < nabin; i++)
3277 fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*sum));
3279 xvgrclose(fp);
3282 /* Print HB in alpha-helix */
3283 if (opt2bSet("-hx", NFILE, fnm))
3285 fp = xvgropen(opt2fn("-hx", NFILE, fnm),
3286 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
3287 xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
3288 for (i = 0; i < nframes; i++)
3290 fprintf(fp, "%10g", hb->time[i]);
3291 for (j = 0; j < max_hx; j++)
3293 fprintf(fp, " %6d", hb->nhx[i][j]);
3295 fprintf(fp, "\n");
3297 xvgrclose(fp);
3300 printf("Average number of %s per timeframe %.3f out of %g possible\n",
3301 bContact ? "contacts" : "hbonds",
3302 bContact ? aver_dist : aver_nhb, max_nhb);
3304 /* Do Autocorrelation etc. */
3305 if (hb->bHBmap)
3308 Added support for -contact in ac and hbm calculations below.
3309 - Erik Marklund, May 29, 2006
3311 if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
3313 please_cite(stdout, "Spoel2006b");
3315 if (opt2bSet("-ac", NFILE, fnm))
3317 do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
3318 bMerge, bContact, fit_start, temp, r2cut > 0, oenv,
3319 nThreads);
3321 if (opt2bSet("-life", NFILE, fnm))
3323 do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
3325 if (opt2bSet("-hbm", NFILE, fnm))
3327 t_matrix mat;
3328 int id, ia, hh, x, y;
3329 mat.flags = mat.y0 = 0;
3331 if ((nframes > 0) && (hb->nrhb > 0))
3333 mat.nx = nframes;
3334 mat.ny = hb->nrhb;
3336 snew(mat.matrix, mat.nx);
3337 for (x = 0; (x < mat.nx); x++)
3339 snew(mat.matrix[x], mat.ny);
3341 y = 0;
3342 for (id = 0; (id < hb->d.nrd); id++)
3344 for (ia = 0; (ia < hb->a.nra); ia++)
3346 for (hh = 0; (hh < hb->maxhydro); hh++)
3348 if (hb->hbmap[id][ia])
3350 if (ISHB(hb->hbmap[id][ia]->history[hh]))
3352 for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
3354 int nn0 = hb->hbmap[id][ia]->n0;
3355 range_check(y, 0, mat.ny);
3356 mat.matrix[x+nn0][y] = static_cast<t_matelmt>(is_hb(hb->hbmap[id][ia]->h[hh], x));
3358 y++;
3364 mat.axis_x = hb->time;
3365 snew(mat.axis_y, mat.ny);
3366 for (j = 0; j < mat.ny; j++)
3368 mat.axis_y[j] = j;
3370 sprintf(mat.title, bContact ? "Contact Existence Map" :
3371 "Hydrogen Bond Existence Map");
3372 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
3373 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv).c_str());
3374 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
3375 mat.bDiscrete = TRUE;
3376 mat.nmap = 2;
3377 snew(mat.map, mat.nmap);
3378 for (i = 0; i < mat.nmap; i++)
3380 mat.map[i].code.c1 = hbmap[i];
3381 mat.map[i].desc = hbdesc[i];
3382 mat.map[i].rgb = hbrgb[i];
3384 fp = opt2FILE("-hbm", NFILE, fnm, "w");
3385 write_xpm_m(fp, mat);
3386 gmx_ffclose(fp);
3387 for (x = 0; x < mat.nx; x++)
3389 sfree(mat.matrix[x]);
3391 sfree(mat.axis_y);
3392 sfree(mat.matrix);
3393 sfree(mat.map);
3395 else
3397 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
3402 if (hb->bDAnr)
3404 int i, j, nleg;
3405 char **legnames;
3406 char buf[STRLEN];
3408 #define USE_THIS_GROUP(j) ( ((j) == gr0) || (bTwo && ((j) == gr1)) )
3410 fp = xvgropen(opt2fn("-dan", NFILE, fnm),
3411 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
3412 nleg = (bTwo ? 2 : 1)*2;
3413 snew(legnames, nleg);
3414 i = 0;
3415 for (j = 0; j < grNR; j++)
3417 if (USE_THIS_GROUP(j) )
3419 sprintf(buf, "Donors %s", grpnames[j]);
3420 legnames[i++] = gmx_strdup(buf);
3421 sprintf(buf, "Acceptors %s", grpnames[j]);
3422 legnames[i++] = gmx_strdup(buf);
3425 if (i != nleg)
3427 gmx_incons("number of legend entries");
3429 xvgr_legend(fp, nleg, legnames, oenv);
3430 for (i = 0; i < nframes; i++)
3432 fprintf(fp, "%10g", hb->time[i]);
3433 for (j = 0; (j < grNR); j++)
3435 if (USE_THIS_GROUP(j) )
3437 fprintf(fp, " %6d", hb->danr[i][j]);
3440 fprintf(fp, "\n");
3442 xvgrclose(fp);
3445 return 0;