2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/correlationfunctions/autocorr.h"
45 #include "gromacs/correlationfunctions/crosscorr.h"
46 #include "gromacs/correlationfunctions/expfit.h"
47 #include "gromacs/correlationfunctions/integrate.h"
48 #include "gromacs/fileio/matio.h"
49 #include "gromacs/fileio/tpxio.h"
50 #include "gromacs/fileio/trxio.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/gmxana/geminate.h"
53 #include "gromacs/gmxana/gmx_ana.h"
54 #include "gromacs/legacyheaders/copyrite.h"
55 #include "gromacs/legacyheaders/macros.h"
56 #include "gromacs/legacyheaders/txtdump.h"
57 #include "gromacs/legacyheaders/viewit.h"
58 #include "gromacs/math/units.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/topology/index.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/gmxomp.h"
66 #include "gromacs/utility/smalloc.h"
67 #include "gromacs/utility/snprintf.h"
71 /*#define HAVE_NN_LOOPS*/
73 typedef short int t_E;
76 typedef int t_hx[max_hx];
77 #define NRHXTYPES max_hx
78 const char *hxtypenames[NRHXTYPES] =
79 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
83 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == 0)
85 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == (threadNr))
88 /* -----------------------------------------*/
94 hbNo, hbDist, hbHB, hbNR, hbR2
97 noDA, ACC, DON, DA, INGROUP
100 NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR
103 static const char *grpnames[grNR] = {"0", "1", "I" };
105 static gmx_bool bDebug = FALSE;
110 #define HB_YESINS HB_YES|HB_INS
114 #define ISHB(h) (((h) & 2) == 2)
115 #define ISDIST(h) (((h) & 1) == 1)
116 #define ISDIST2(h) (((h) & 4) == 4)
117 #define ISACC(h) (((h) & 1) == 1)
118 #define ISDON(h) (((h) & 2) == 2)
119 #define ISINGRP(h) (((h) & 4) == 4)
132 typedef int t_icell[grNR];
133 typedef atom_id h_id[MAXHYDRO];
136 int history[MAXHYDRO];
137 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
138 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
140 /* Bitmask array which tells whether a hbond is present
141 * at a given time. Either of these may be NULL
143 int n0; /* First frame a HB was found */
144 int nframes, maxframes; /* Amount of frames in this hbond */
147 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
148 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
149 * acceptor distance is less than the user-specified distance (typically
156 atom_id *acc; /* Atom numbers of the acceptors */
157 int *grp; /* Group index */
158 int *aptr; /* Map atom number to acceptor index */
163 int *don; /* Atom numbers of the donors */
164 int *grp; /* Group index */
165 int *dptr; /* Map atom number to donor index */
166 int *nhydro; /* Number of hydrogens for each donor */
167 h_id *hydro; /* The atom numbers of the hydrogens */
168 h_id *nhbonds; /* The number of HBs per H at current */
171 /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/
175 int len; /* The length of frame and p. */
176 int *frame; /* The frames at which transitio*/
181 /* Periodicity history. Used for the reversible geminate recombination. */
182 t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap:
183 * pHist[d][a]. We can safely assume that the same
184 * periodic shift holds for all hydrogens of a da-pair.
186 * Nowadays it only stores TRANSITIONS, and not the shift at every frame.
187 * That saves a LOT of memory, an hopefully kills a mysterious bug where
188 * pHist gets contaminated. */
190 PSTYPE nper; /* The length of p2i */
191 ivec *p2i; /* Maps integer to periodic shift for a pair.*/
192 matrix P; /* Projection matrix to find the box shifts. */
193 int gemtype; /* enumerated type */
198 int *Etot; /* Total energy for each frame */
199 t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */
203 gmx_bool bHBmap, bDAnr, bGem;
205 /* The following arrays are nframes long */
206 int nframes, max_frames, maxhydro;
212 /* These structures are initialized from the topology at start up */
215 /* This holds a matrix with all possible hydrogen bonds */
221 /* For parallelization reasons this will have to be a pointer.
222 * Otherwise discrepancies may arise between the periodicity data
223 * seen by different threads. */
227 static void clearPshift(t_pShift *pShift)
232 sfree(pShift->frame);
237 static void calcBoxProjection(matrix B, matrix P)
239 const int vp[] = {XX, YY, ZZ};
244 for (i = 0; i < 3; i++)
247 for (j = 0; j < 3; j++)
250 U[m][n] = i == j ? 1 : 0;
254 for (i = 0; i < 3; i++)
257 mvmul(M, U[m], P[m]);
262 static void calcBoxDistance(matrix P, rvec d, ivec ibd)
264 /* returns integer distance in box coordinates.
265 * P is the projection matrix from cartesian coordinates
266 * obtained with calcBoxProjection(). */
270 /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
271 for (i = 0; i < 3; i++)
273 bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5);
275 ibd[XX] = (int)bd[XX];
276 ibd[YY] = (int)bd[YY];
277 ibd[ZZ] = (int)bd[ZZ];
280 /* Changed argument 'bMerge' into 'oneHB' below,
281 * since -contact should cause maxhydro to be 1,
283 * - Erik Marklund May 29, 2006
286 static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap)
288 /* Try to merge hbonds on the fly. That means that if the
289 * acceptor and donor are mergable, then:
290 * 1) store the hb-info so that acceptor id > donor id,
291 * 2) add the periodic shift in pairs, so that [-x,-y,-z] is
292 * stored in per.p2i[] whenever acceptor id < donor id.
293 * Note that [0,0,0] should already be the first element of per.p2i
294 * by the time this function is called. */
296 /* daSwap is TRUE if the donor and acceptor were swapped.
297 * If so, then the negative vector should be used. */
300 if (per->p2i == NULL || per->nper == 0)
302 gmx_fatal(FARGS, "'per' not initialized properly.");
304 for (i = 0; i < per->nper; i++)
306 if (r[XX] == per->p2i[i][XX] &&
307 r[YY] == per->p2i[i][YY] &&
308 r[ZZ] == per->p2i[i][ZZ])
313 /* Not found apparently. Add it to the list! */
314 /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
320 fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
325 srenew(per->p2i, per->nper+2);
327 copy_ivec(r, per->p2i[per->nper]);
330 /* Add the mirror too. It's rather likely that it'll be needed. */
331 per->p2i[per->nper][XX] = -r[XX];
332 per->p2i[per->nper][YY] = -r[YY];
333 per->p2i[per->nper][ZZ] = -r[ZZ];
336 return per->nper - 1 - (daSwap ? 0 : 1);
339 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode)
344 hb->wordlen = 8*sizeof(unsigned int);
354 hb->maxhydro = MAXHYDRO;
357 hb->per->gemtype = bGem ? gemmode : 0;
362 static void mk_hbmap(t_hbdata *hb)
366 snew(hb->hbmap, hb->d.nrd);
367 for (i = 0; (i < hb->d.nrd); i++)
369 snew(hb->hbmap[i], hb->a.nra);
370 if (hb->hbmap[i] == NULL)
372 gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
374 for (j = 0; (j > hb->a.nra); j++)
376 hb->hbmap[i][j] = NULL;
381 /* Consider redoing pHist so that is only stores transitions between
382 * periodicities and not the periodicity for all frames. This eats heaps of memory. */
383 static void mk_per(t_hbdata *hb)
388 snew(hb->per->pHist, hb->d.nrd);
389 for (i = 0; i < hb->d.nrd; i++)
391 snew(hb->per->pHist[i], hb->a.nra);
392 if (hb->per->pHist[i] == NULL)
394 gmx_fatal(FARGS, "Could not allocate enough memory for per->pHist");
396 for (j = 0; j < hb->a.nra; j++)
398 clearPshift(&(hb->per->pHist[i][j]));
401 /* add the [0,0,0] shift to element 0 of p2i. */
402 snew(hb->per->p2i, 1);
403 clear_ivec(hb->per->p2i[0]);
409 static void mk_hbEmap (t_hbdata *hb, int n0)
414 snew(hb->hbE.E, hb->d.nrd);
415 for (i = 0; i < hb->d.nrd; i++)
417 snew(hb->hbE.E[i], hb->a.nra);
418 for (j = 0; j < hb->a.nra; j++)
420 snew(hb->hbE.E[i][j], MAXHYDRO);
421 for (k = 0; k < MAXHYDRO; k++)
423 hb->hbE.E[i][j][k] = NULL;
430 static void free_hbEmap (t_hbdata *hb)
433 for (i = 0; i < hb->d.nrd; i++)
435 for (j = 0; j < hb->a.nra; j++)
437 for (k = 0; k < MAXHYDRO; k++)
439 sfree(hb->hbE.E[i][j][k]);
441 sfree(hb->hbE.E[i][j]);
449 static void addFramesNN(t_hbdata *hb, int frame)
452 #define DELTAFRAMES_HBE 10
454 int d, a, h, nframes;
456 if (frame >= hb->hbE.nframes)
458 nframes = hb->hbE.nframes + DELTAFRAMES_HBE;
459 srenew(hb->hbE.Etot, nframes);
461 for (d = 0; d < hb->d.nrd; d++)
463 for (a = 0; a < hb->a.nra; a++)
465 for (h = 0; h < hb->d.nhydro[d]; h++)
467 srenew(hb->hbE.E[d][a][h], nframes);
472 hb->hbE.nframes += DELTAFRAMES_HBE;
476 static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst,
477 matrix box, rvec hbox, t_donors *donors)
482 * alpha - angle between dipoles
483 * x[] - atomic positions
484 * EEst - the type of energy estimate (see enum in hbplugin.h)
485 * box - the box vectors \
486 * hbox - half box lengths _These two are only needed for the pbc correction
491 rvec dipole[2], xmol[3], xmean[2];
497 /* Self-interaction */
504 /* This is a simple binary existence function that sets E=1 whenever
505 * the distance between the oxygens is equal too or less than 0.35 nm.
507 rvec_sub(x[d], x[a], dist);
508 pbc_correct_gem(dist, box, hbox);
509 if (norm(dist) <= 0.35)
520 /* Negative potential energy of a dipole.
521 * E = -cos(alpha) * 1/r^3 */
523 copy_rvec(x[d], xmol[0]); /* donor */
524 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
525 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
527 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
528 rvec_inc(xmean[0], xmol[1]);
529 rvec_inc(xmean[0], xmol[2]);
530 for (i = 0; i < 3; i++)
532 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
535 /* Assumes that all acceptors are also donors. */
536 copy_rvec(x[a], xmol[0]); /* acceptor */
537 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
538 copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */
541 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
542 rvec_inc(xmean[1], xmol[1]);
543 rvec_inc(xmean[1], xmol[2]);
544 for (i = 0; i < 3; i++)
546 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
549 rvec_sub(xmean[0], xmean[1], dist);
550 pbc_correct_gem(dist, box, hbox);
553 realE = pow(r, -3.0);
554 E = (t_E)(SCALEFACTOR_E * realE);
558 /* Negative potential energy of a (unpolarizable) dipole.
559 * E = -cos(alpha) * 1/r^3 */
560 clear_rvec(dipole[1]);
561 clear_rvec(dipole[0]);
563 copy_rvec(x[d], xmol[0]); /* donor */
564 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
565 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
567 rvec_inc(dipole[0], xmol[1]);
568 rvec_inc(dipole[0], xmol[2]);
569 for (i = 0; i < 3; i++)
573 rvec_dec(dipole[0], xmol[0]);
575 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
576 rvec_inc(xmean[0], xmol[1]);
577 rvec_inc(xmean[0], xmol[2]);
578 for (i = 0; i < 3; i++)
580 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
583 /* Assumes that all acceptors are also donors. */
584 copy_rvec(x[a], xmol[0]); /* acceptor */
585 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
586 copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */
589 rvec_inc(dipole[1], xmol[1]);
590 rvec_inc(dipole[1], xmol[2]);
591 for (i = 0; i < 3; i++)
595 rvec_dec(dipole[1], xmol[0]);
597 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
598 rvec_inc(xmean[1], xmol[1]);
599 rvec_inc(xmean[1], xmol[2]);
600 for (i = 0; i < 3; i++)
602 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
605 rvec_sub(xmean[0], xmean[1], dist);
606 pbc_correct_gem(dist, box, hbox);
609 double cosalpha = cos_angle(dipole[0], dipole[1]);
610 realE = cosalpha * pow(r, -3.0);
611 E = (t_E)(SCALEFACTOR_E * realE);
615 printf("Can't do that type of energy estimate: %i\n.", EEst);
622 static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame)
624 /* hb - hbond data structure
628 E - estimate of the energy
629 frame - the current frame.
632 /* Store the estimated energy */
638 hb->hbE.E[d][a][h][frame] = E;
642 hb->hbE.Etot[frame] += E;
645 #endif /* HAVE_NN_LOOPS */
648 /* Finds -v[] in the periodicity index */
649 static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
653 for (i = 0; i < nper; i++)
655 if (v[i][XX] == -(v[p][XX]) &&
656 v[i][YY] == -(v[p][YY]) &&
657 v[i][ZZ] == -(v[p][ZZ]))
662 printf("Couldn't find mirror of [%i, %i, %i], index \n",
670 static void add_frames(t_hbdata *hb, int nframes)
674 if (nframes >= hb->max_frames)
676 hb->max_frames += 4096;
677 srenew(hb->time, hb->max_frames);
678 srenew(hb->nhb, hb->max_frames);
679 srenew(hb->ndist, hb->max_frames);
680 srenew(hb->n_bound, hb->max_frames);
681 srenew(hb->nhx, hb->max_frames);
684 srenew(hb->danr, hb->max_frames);
687 hb->nframes = nframes;
690 #define OFFSET(frame) (frame / 32)
691 #define MASK(frame) (1 << (frame % 32))
693 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
697 hbexist[OFFSET(frame)] |= MASK(frame);
701 hbexist[OFFSET(frame)] &= ~MASK(frame);
705 static gmx_bool is_hb(unsigned int hbexist[], int frame)
707 return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
710 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
712 unsigned int *ghptr = NULL;
716 ghptr = hb->hbmap[id][ia]->h[ih];
718 else if (ihb == hbDist)
720 ghptr = hb->hbmap[id][ia]->g[ih];
724 gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
727 _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
730 static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
734 snew(pHist->frame, 1);
737 pHist->frame[0] = frame;
742 if (pHist->p[pHist->len-1] != p)
745 srenew(pHist->frame, pHist->len);
746 srenew(pHist->p, pHist->len);
747 pHist->frame[pHist->len-1] = frame;
748 pHist->p[pHist->len-1] = p;
749 } /* Otherwise, there is no transition. */
753 static PSTYPE getPshift(t_pShift pHist, int frame)
758 || (pHist.len > 0 && pHist.frame[0] > frame))
763 for (i = 0; i < pHist.len; i++)
776 /* It seems that frame is after the last periodic transition. Return the last periodicity. */
777 return pHist.p[pHist.len-1];
780 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPE p)
783 t_hbond *hb = hbd->hbmap[id][ia];
784 int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]);
785 int wlen = hbd->wordlen;
787 gmx_bool bGem = hbd->bGem;
792 hb->maxframes = delta;
793 for (i = 0; (i < maxhydro); i++)
795 snew(hb->h[i], hb->maxframes/wlen);
796 snew(hb->g[i], hb->maxframes/wlen);
801 hb->nframes = frame-hb->n0;
802 /* We need a while loop here because hbonds may be returning
805 while (hb->nframes >= hb->maxframes)
807 n = hb->maxframes + delta;
808 for (i = 0; (i < maxhydro); i++)
810 srenew(hb->h[i], n/wlen);
811 srenew(hb->g[i], n/wlen);
812 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
824 set_hb(hbd, id, h, ia, frame, ihb);
827 if (p >= hbd->per->nper)
829 gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
833 addPshift(&(hbd->per->pHist[id][ia]), p, frame);
841 static void inc_nhbonds(t_donors *ddd, int d, int h)
844 int dptr = ddd->dptr[d];
846 for (j = 0; (j < ddd->nhydro[dptr]); j++)
848 if (ddd->hydro[dptr][j] == h)
850 ddd->nhbonds[dptr][j]++;
854 if (j == ddd->nhydro[dptr])
856 gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
860 static int _acceptor_index(t_acceptors *a, int grp, atom_id i,
861 const char *file, int line)
865 if (a->grp[ai] != grp)
869 fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
870 ai, a->grp[ai], grp, file, line);
879 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
881 static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line)
890 if (d->grp[di] != grp)
894 fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
895 di, d->grp[di], grp, file, line);
904 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
906 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
908 /* g_hbond doesn't allow overlapping groups */
914 donor_index(&hb->d, grpd, a) != NOTSET
915 && acceptor_index(&hb->a, grpa, d) != NOTSET;
919 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
920 int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPE p)
923 gmx_bool daSwap = FALSE;
925 if ((id = hb->d.dptr[d]) == NOTSET)
927 gmx_fatal(FARGS, "No donor atom %d", d+1);
929 else if (grpd != hb->d.grp[id])
931 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
932 grpd, hb->d.grp[id], d+1);
934 if ((ia = hb->a.aptr[a]) == NOTSET)
936 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
938 else if (grpa != hb->a.grp[ia])
940 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
941 grpa, hb->a.grp[ia], a+1);
947 if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
948 /* Then swap identity so that the id of d is lower then that of a.
950 * This should really be redundant by now, as is_hbond() now ought to return
951 * hbNo in the cases where this conditional is TRUE. */
958 /* Now repeat donor/acc check. */
959 if ((id = hb->d.dptr[d]) == NOTSET)
961 gmx_fatal(FARGS, "No donor atom %d", d+1);
963 else if (grpd != hb->d.grp[id])
965 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
966 grpd, hb->d.grp[id], d+1);
968 if ((ia = hb->a.aptr[a]) == NOTSET)
970 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
972 else if (grpa != hb->a.grp[ia])
974 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
975 grpa, hb->a.grp[ia], a+1);
982 /* Loop over hydrogens to find which hydrogen is in this particular HB */
983 if ((ihb == hbHB) && !bMerge && !bContact)
985 for (k = 0; (k < hb->d.nhydro[id]); k++)
987 if (hb->d.hydro[id][k] == h)
992 if (k == hb->d.nhydro[id])
994 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
1006 #pragma omp critical
1008 if (hb->hbmap[id][ia] == NULL)
1010 snew(hb->hbmap[id][ia], 1);
1011 snew(hb->hbmap[id][ia]->h, hb->maxhydro);
1012 snew(hb->hbmap[id][ia]->g, hb->maxhydro);
1014 add_ff(hb, id, k, ia, frame, ihb, p);
1018 /* Strange construction with frame >=0 is a relic from old code
1019 * for selected hbond analysis. It may be necessary again if that
1020 * is made to work again.
1024 hh = hb->hbmap[id][ia]->history[k];
1030 hb->hbmap[id][ia]->history[k] = hh | 2;
1041 hb->hbmap[id][ia]->history[k] = hh | 1;
1065 if (bMerge && daSwap)
1067 h = hb->d.hydro[id][0];
1069 /* Increment number if HBonds per H */
1070 if (ihb == hbHB && !bContact)
1072 inc_nhbonds(&(hb->d), d, h);
1076 static char *mkatomname(t_atoms *atoms, int i)
1078 static char buf[32];
1081 rnr = atoms->atom[i].resind;
1082 sprintf(buf, "%4s%d%-4s",
1083 *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
1088 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms)
1090 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
1093 for (i = 0; i < isize; i++)
1095 if (index[i] >= natoms)
1097 gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
1099 datable[index[i]] |= INGROUP;
1103 static void clear_datable_grp(unsigned char *datable, int size)
1105 /* Clears group information from the table */
1107 const char mask = !(char)INGROUP;
1110 for (i = 0; i < size; i++)
1117 static void add_acc(t_acceptors *a, int ia, int grp)
1119 if (a->nra >= a->max_nra)
1122 srenew(a->acc, a->max_nra);
1123 srenew(a->grp, a->max_nra);
1125 a->grp[a->nra] = grp;
1126 a->acc[a->nra++] = ia;
1129 static void search_acceptors(t_topology *top, int isize,
1130 atom_id *index, t_acceptors *a, int grp,
1132 gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
1138 for (i = 0; (i < isize); i++)
1142 (((*top->atoms.atomname[n])[0] == 'O') ||
1143 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
1144 ISINGRP(datable[n]))
1146 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
1151 snew(a->aptr, top->atoms.nr);
1152 for (i = 0; (i < top->atoms.nr); i++)
1154 a->aptr[i] = NOTSET;
1156 for (i = 0; (i < a->nra); i++)
1158 a->aptr[a->acc[i]] = i;
1162 static void add_h2d(int id, int ih, t_donors *ddd)
1166 for (i = 0; (i < ddd->nhydro[id]); i++)
1168 if (ddd->hydro[id][i] == ih)
1170 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
1175 if (i == ddd->nhydro[id])
1177 if (ddd->nhydro[id] >= MAXHYDRO)
1179 gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
1180 ddd->don[id], MAXHYDRO);
1182 ddd->hydro[id][i] = ih;
1187 static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable)
1191 if (!datable || ISDON(datable[id]))
1193 if (ddd->dptr[id] == NOTSET) /* New donor */
1205 if (ddd->nrd >= ddd->max_nrd)
1207 ddd->max_nrd += 128;
1208 srenew(ddd->don, ddd->max_nrd);
1209 srenew(ddd->nhydro, ddd->max_nrd);
1210 srenew(ddd->hydro, ddd->max_nrd);
1211 srenew(ddd->nhbonds, ddd->max_nrd);
1212 srenew(ddd->grp, ddd->max_nrd);
1214 ddd->don[ddd->nrd] = id;
1215 ddd->nhydro[ddd->nrd] = 0;
1216 ddd->grp[ddd->nrd] = grp;
1223 add_h2d(i, ih, ddd);
1228 printf("Warning: Atom %d is not in the d/a-table!\n", id);
1232 static void search_donors(t_topology *top, int isize, atom_id *index,
1233 t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
1234 unsigned char *datable)
1237 t_functype func_type;
1238 t_ilist *interaction;
1239 atom_id nr1, nr2, nr3;
1244 snew(ddd->dptr, top->atoms.nr);
1245 for (i = 0; (i < top->atoms.nr); i++)
1247 ddd->dptr[i] = NOTSET;
1255 for (i = 0; (i < isize); i++)
1257 datable[index[i]] |= DON;
1258 add_dh(ddd, index[i], -1, grp, datable);
1264 for (func_type = 0; (func_type < F_NRE); func_type++)
1266 interaction = &(top->idef.il[func_type]);
1267 if (func_type == F_POSRES || func_type == F_FBPOSRES)
1269 /* The ilist looks strange for posre. Bug in grompp?
1270 * We don't need posre interactions for hbonds anyway.*/
1273 for (i = 0; i < interaction->nr;
1274 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1277 if (func_type != top->idef.functype[interaction->iatoms[i]])
1279 fprintf(stderr, "Error in func_type %s",
1280 interaction_function[func_type].longname);
1284 /* check out this functype */
1285 if (func_type == F_SETTLE)
1287 nr1 = interaction->iatoms[i+1];
1288 nr2 = interaction->iatoms[i+2];
1289 nr3 = interaction->iatoms[i+3];
1291 if (ISINGRP(datable[nr1]))
1293 if (ISINGRP(datable[nr2]))
1295 datable[nr1] |= DON;
1296 add_dh(ddd, nr1, nr1+1, grp, datable);
1298 if (ISINGRP(datable[nr3]))
1300 datable[nr1] |= DON;
1301 add_dh(ddd, nr1, nr1+2, grp, datable);
1305 else if (IS_CHEMBOND(func_type))
1307 for (j = 0; j < 2; j++)
1309 nr1 = interaction->iatoms[i+1+j];
1310 nr2 = interaction->iatoms[i+2-j];
1311 if ((*top->atoms.atomname[nr1][0] == 'H') &&
1312 ((*top->atoms.atomname[nr2][0] == 'O') ||
1313 (*top->atoms.atomname[nr2][0] == 'N')) &&
1314 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1316 datable[nr2] |= DON;
1317 add_dh(ddd, nr2, nr1, grp, datable);
1324 for (func_type = 0; func_type < F_NRE; func_type++)
1326 interaction = &top->idef.il[func_type];
1327 for (i = 0; i < interaction->nr;
1328 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1331 if (func_type != top->idef.functype[interaction->iatoms[i]])
1333 gmx_incons("function type in search_donors");
1336 if (interaction_function[func_type].flags & IF_VSITE)
1338 nr1 = interaction->iatoms[i+1];
1339 if (*top->atoms.atomname[nr1][0] == 'H')
1343 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1354 if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1355 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1356 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1358 datable[nr2] |= DON;
1359 add_dh(ddd, nr2, nr1, grp, datable);
1369 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
1376 for (i = 0; i < DIM; i++)
1378 ngrid[i] = (box[i][i]/(1.2*rcut));
1382 if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
1384 for (i = 0; i < DIM; i++)
1391 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1392 ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
1394 snew(grid, ngrid[ZZ]);
1395 for (z = 0; z < ngrid[ZZ]; z++)
1397 snew((grid)[z], ngrid[YY]);
1398 for (y = 0; y < ngrid[YY]; y++)
1400 snew((grid)[z][y], ngrid[XX]);
1406 static void reset_nhbonds(t_donors *ddd)
1410 for (i = 0; (i < ddd->nrd); i++)
1412 for (j = 0; (j < MAXHH); j++)
1414 ddd->nhbonds[i][j] = 0;
1419 void pbc_correct_gem(rvec dx, matrix box, rvec hbox);
1421 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
1422 gmx_bool bBox, matrix box, rvec hbox,
1423 real rcut, real rshell,
1424 ivec ngrid, t_gridcell ***grid)
1426 int i, m, gr, xi, yi, zi, nr;
1429 rvec invdelta, dshell, xtemp = {0, 0, 0};
1431 gmx_bool bDoRshell, bInShell, bAcc;
1436 bDoRshell = (rshell > 0);
1437 rshell2 = sqr(rshell);
1440 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
1442 for (m = 0; m < DIM; m++)
1444 hbox[m] = box[m][m]*0.5;
1447 invdelta[m] = ngrid[m]/box[m][m];
1448 if (1/invdelta[m] < rcut)
1450 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
1451 "%s can not handle this situation, sorry.\n",
1464 /* resetting atom counts */
1465 for (gr = 0; (gr < grNR); gr++)
1467 for (zi = 0; zi < ngrid[ZZ]; zi++)
1469 for (yi = 0; yi < ngrid[YY]; yi++)
1471 for (xi = 0; xi < ngrid[XX]; xi++)
1473 grid[zi][yi][xi].d[gr].nr = 0;
1474 grid[zi][yi][xi].a[gr].nr = 0;
1480 /* put atoms in grid cells */
1481 for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
1494 for (i = 0; (i < nr); i++)
1496 /* check if we are inside the shell */
1497 /* if bDoRshell=FALSE then bInShell=TRUE always */
1502 rvec_sub(x[ad[i]], xshell, dshell);
1505 if (FALSE && !hb->bGem)
1507 for (m = DIM-1; m >= 0 && bInShell; m--)
1509 if (dshell[m] < -hbox[m])
1511 rvec_inc(dshell, box[m]);
1513 else if (dshell[m] >= hbox[m])
1515 dshell[m] -= 2*hbox[m];
1517 /* if we're outside the cube, we're outside the sphere also! */
1518 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1526 gmx_bool bDone = FALSE;
1530 for (m = DIM-1; m >= 0 && bInShell; m--)
1532 if (dshell[m] < -hbox[m])
1535 rvec_inc(dshell, box[m]);
1537 if (dshell[m] >= hbox[m])
1540 dshell[m] -= 2*hbox[m];
1544 for (m = DIM-1; m >= 0 && bInShell; m--)
1546 /* if we're outside the cube, we're outside the sphere also! */
1547 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1554 /* if we're inside the cube, check if we're inside the sphere */
1557 bInShell = norm2(dshell) < rshell2;
1567 copy_rvec(x[ad[i]], xtemp);
1569 pbc_correct_gem(x[ad[i]], box, hbox);
1571 for (m = DIM-1; m >= 0; m--)
1573 if (TRUE || !hb->bGem)
1575 /* put atom in the box */
1576 while (x[ad[i]][m] < 0)
1578 rvec_inc(x[ad[i]], box[m]);
1580 while (x[ad[i]][m] >= box[m][m])
1582 rvec_dec(x[ad[i]], box[m]);
1585 /* determine grid index of atom */
1586 grididx[m] = x[ad[i]][m]*invdelta[m];
1587 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1591 copy_rvec(xtemp, x[ad[i]]); /* copy back */
1596 range_check(gx, 0, ngrid[XX]);
1597 range_check(gy, 0, ngrid[YY]);
1598 range_check(gz, 0, ngrid[ZZ]);
1602 /* add atom to grid cell */
1605 newgrid = &(grid[gz][gy][gx].a[gr]);
1609 newgrid = &(grid[gz][gy][gx].d[gr]);
1611 if (newgrid->nr >= newgrid->maxnr)
1613 newgrid->maxnr += 10;
1614 DBB(newgrid->maxnr);
1615 srenew(newgrid->atoms, newgrid->maxnr);
1618 newgrid->atoms[newgrid->nr] = ad[i];
1626 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1630 for (gr = 0; (gr < grNR); gr++)
1633 for (zi = 0; zi < ngrid[ZZ]; zi++)
1635 for (yi = 0; yi < ngrid[YY]; yi++)
1637 for (xi = 0; xi < ngrid[XX]; xi++)
1639 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1647 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1648 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1649 * With a triclinic box all loops are 3 long, except when a cell is
1650 * located next to one of the box edges which is not parallel to the
1651 * x/y-plane, in that case all cells in a line or layer are searched.
1652 * This could be implemented slightly more efficient, but the code
1653 * would get much more complicated.
1655 static gmx_inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1657 return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
1659 static gmx_inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1661 return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1663 static gmx_inline int grid_mod(int j, int n)
1668 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1670 int gr, x, y, z, sum[grNR];
1672 fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1673 for (gr = 0; gr < grNR; gr++)
1676 fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1677 for (z = 0; z < ngrid[ZZ]; z += 2)
1679 fprintf(fp, "Z=%d,%d\n", z, z+1);
1680 for (y = 0; y < ngrid[YY]; y++)
1682 for (x = 0; x < ngrid[XX]; x++)
1684 fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1685 sum[gr] += grid[z][y][x].d[gr].nr;
1686 fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1687 sum[gr] += grid[z][y][x].a[gr].nr;
1691 if ( (z+1) < ngrid[ZZ])
1693 for (x = 0; x < ngrid[XX]; x++)
1695 fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1696 sum[gr] += grid[z+1][y][x].d[gr].nr;
1697 fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1698 sum[gr] += grid[z+1][y][x].a[gr].nr;
1705 fprintf(fp, "TOTALS:");
1706 for (gr = 0; gr < grNR; gr++)
1708 fprintf(fp, " %d=%d", gr, sum[gr]);
1713 /* New GMX record! 5 * in a row. Congratulations!
1714 * Sorry, only four left.
1716 static void free_grid(ivec ngrid, t_gridcell ****grid)
1719 t_gridcell ***g = *grid;
1721 for (z = 0; z < ngrid[ZZ]; z++)
1723 for (y = 0; y < ngrid[YY]; y++)
1733 void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1736 gmx_bool bDone = FALSE;
1740 for (m = DIM-1; m >= 0; m--)
1742 if (dx[m] < -hbox[m])
1745 rvec_inc(dx, box[m]);
1747 if (dx[m] >= hbox[m])
1750 rvec_dec(dx, box[m]);
1756 /* Added argument r2cut, changed contact and implemented
1757 * use of second cut-off.
1758 * - Erik Marklund, June 29, 2006
1760 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1761 real rcut, real r2cut, real ccut,
1762 rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1763 real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1764 gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1766 int h, hh, id, ja, ihb;
1767 rvec r_da, r_ha, r_dh, r = {0, 0, 0};
1769 real rc2, r2c2, rda2, rha2, ca;
1770 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1771 gmx_bool daSwap = FALSE;
1778 if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1779 ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
1787 rvec_sub(x[d], x[a], r_da);
1788 /* Insert projection code here */
1790 if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1792 /* Then this hbond/contact will be found again, or it has already been found. */
1797 if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1798 { /* return hbNo; */
1799 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1803 copy_rvec(r_da, r); /* Save this for later */
1804 pbc_correct_gem(r_da, box, hbox);
1808 pbc_correct_gem(r_da, box, hbox);
1811 rda2 = iprod(r_da, r_da);
1815 if (daSwap && grpa == grpd)
1823 calcBoxDistance(hb->per->P, r, ri);
1824 *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */
1828 else if (rda2 < r2c2)
1839 if (bDA && (rda2 > rc2))
1844 for (h = 0; (h < hb->d.nhydro[id]); h++)
1846 hh = hb->d.hydro[id][h];
1850 rvec_sub(x[hh], x[a], r_ha);
1853 pbc_correct_gem(r_ha, box, hbox);
1855 rha2 = iprod(r_ha, r_ha);
1860 calcBoxDistance(hb->per->P, r, ri);
1861 *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */
1864 if (bDA || (!bDA && (rha2 <= rc2)))
1866 rvec_sub(x[d], x[hh], r_dh);
1869 pbc_correct_gem(r_dh, box, hbox);
1876 ca = cos_angle(r_dh, r_da);
1877 /* if angle is smaller, cos is larger */
1881 *d_ha = sqrt(bDA ? rda2 : rha2);
1887 if (bDA || (!bDA && HAinrange))
1897 /* Fixed previously undiscovered bug in the merge
1898 code, where the last frame of each hbond disappears.
1899 - Erik Marklund, June 1, 2006 */
1900 /* Added the following arguments:
1901 * ptmp[] - temporary periodicity hisory
1902 * a1 - identity of first acceptor/donor
1903 * a2 - identity of second acceptor/donor
1904 * - Erik Marklund, FEB 20 2010 */
1906 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1907 * Will do some more testing before removing the function entirely.
1908 * - Erik Marklund, MAY 10 2010 */
1909 static void do_merge(t_hbdata *hb, int ntmp,
1910 unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[],
1911 t_hbond *hb0, t_hbond *hb1, int a1, int a2)
1913 /* Here we need to make sure we're treating periodicity in
1914 * the right way for the geminate recombination kinetics. */
1916 int m, mm, n00, n01, nn0, nnframes;
1920 /* Decide where to start from when merging */
1923 nn0 = min(n00, n01);
1924 nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1925 /* Initiate tmp arrays */
1926 for (m = 0; (m < ntmp); m++)
1932 /* Fill tmp arrays with values due to first HB */
1933 /* Once again '<' had to be replaced with '<='
1934 to catch the last frame in which the hbond
1936 - Erik Marklund, June 1, 2006 */
1937 for (m = 0; (m <= hb0->nframes); m++)
1940 htmp[mm] = is_hb(hb0->h[0], m);
1943 pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1944 if (pm > hb->per->nper)
1946 gmx_fatal(FARGS, "Illegal shift!");
1950 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1954 /* If we're doing geminate recompbination we usually don't need the distances.
1955 * Let's save some memory and time. */
1956 if (TRUE || !hb->bGem || hb->per->gemtype == gemAD)
1958 for (m = 0; (m <= hb0->nframes); m++)
1961 gtmp[mm] = is_hb(hb0->g[0], m);
1965 for (m = 0; (m <= hb1->nframes); m++)
1968 htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1969 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1970 if (hb->bGem /* && ptmp[mm] != 0 */)
1973 /* If this hbond has been seen before with donor and acceptor swapped,
1974 * then we need to find the mirrored (*-1) periodicity vector to truely
1975 * merge the hbond history. */
1976 pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper);
1977 /* Store index of mirror */
1978 if (pm > hb->per->nper)
1980 gmx_fatal(FARGS, "Illegal shift!");
1985 /* Reallocate target array */
1986 if (nnframes > hb0->maxframes)
1988 srenew(hb0->h[0], 4+nnframes/hb->wordlen);
1989 srenew(hb0->g[0], 4+nnframes/hb->wordlen);
1991 if (NULL != hb->per->pHist)
1993 clearPshift(&(hb->per->pHist[a1][a2]));
1996 /* Copy temp array to target array */
1997 for (m = 0; (m <= nnframes); m++)
1999 _set_hb(hb0->h[0], m, htmp[m]);
2000 _set_hb(hb0->g[0], m, gtmp[m]);
2003 addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
2007 /* Set scalar variables */
2009 hb0->maxframes = nnframes;
2012 /* Added argument bContact for nicer output.
2013 * Erik Marklund, June 29, 2006
2015 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
2017 int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp;
2018 unsigned int *htmp, *gtmp;
2023 indnew = hb->nrdist;
2025 /* Check whether donors are also acceptors */
2026 printf("Merging hbonds with Acceptor and Donor swapped\n");
2028 ntmp = 2*hb->max_frames;
2032 for (i = 0; (i < hb->d.nrd); i++)
2034 fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
2036 ii = hb->a.aptr[id];
2037 for (j = 0; (j < hb->a.nra); j++)
2040 jj = hb->d.dptr[ia];
2041 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
2042 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j]))))
2044 hb0 = hb->hbmap[i][j];
2045 hb1 = hb->hbmap[jj][ii];
2046 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
2048 do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j);
2049 if (ISHB(hb1->history[0]))
2053 else if (ISDIST(hb1->history[0]))
2060 gmx_incons("No contact history");
2064 gmx_incons("Neither hydrogen bond nor distance");
2070 clearPshift(&(hb->per->pHist[jj][ii]));
2074 hb1->history[0] = hbNo;
2079 fprintf(stderr, "\n");
2080 printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
2081 printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
2083 hb->nrdist = indnew;
2089 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
2091 int i, j, k, n_bound[MAXHH], nbtot;
2095 /* Set array to 0 */
2096 for (k = 0; (k < MAXHH); k++)
2100 /* Loop over possible donors */
2101 for (i = 0; (i < hb->d.nrd); i++)
2103 for (j = 0; (j < hb->d.nhydro[i]); j++)
2105 n_bound[hb->d.nhbonds[i][j]]++;
2108 fprintf(fp, "%12.5e", t);
2110 for (k = 0; (k < MAXHH); k++)
2112 fprintf(fp, " %8d", n_bound[k]);
2113 nbtot += n_bound[k]*k;
2115 fprintf(fp, " %8d\n", nbtot);
2118 /* Added argument bContact in do_hblife(...). Also
2119 * added support for -contact in function body.
2120 * - Erik Marklund, May 31, 2006 */
2121 /* Changed the contact code slightly.
2122 * - Erik Marklund, June 29, 2006
2124 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
2125 const output_env_t oenv)
2128 const char *leg[] = { "p(t)", "t p(t)" };
2130 int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
2131 int nframes = hb->nframes;
2134 double sum, integral;
2137 snew(h, hb->maxhydro);
2138 snew(histo, nframes+1);
2139 /* Total number of hbonds analyzed here */
2140 for (i = 0; (i < hb->d.nrd); i++)
2142 for (k = 0; (k < hb->a.nra); k++)
2144 hbh = hb->hbmap[i][k];
2162 for (m = 0; (m < hb->maxhydro); m++)
2166 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
2170 for (nh = 0; (nh < nhydro); nh++)
2175 /* Changed '<' into '<=' below, just like I
2176 did in the hbm-output-loop in the main code.
2177 - Erik Marklund, May 31, 2006
2179 for (j = 0; (j <= hbh->nframes); j++)
2181 ihb = is_hb(h[nh], j);
2182 if (debug && (ndump < 10))
2184 fprintf(debug, "%5d %5d\n", j, ihb);
2204 fprintf(stderr, "\n");
2207 fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
2211 fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
2215 xvgr_legend(fp, asize(leg), leg, oenv);
2217 while ((j0 > 0) && (histo[j0] == 0))
2222 for (i = 0; (i <= j0); i++)
2226 dt = hb->time[1]-hb->time[0];
2229 for (i = 1; (i <= j0); i++)
2231 t = hb->time[i] - hb->time[0] - 0.5*dt;
2232 x1 = t*histo[i]/sum;
2233 fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1);
2238 printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
2239 printf("Note that the lifetime obtained in this manner is close to useless\n");
2240 printf("Use the -ac option instead and check the Forward lifetime\n");
2241 please_cite(stdout, "Spoel2006b");
2246 /* Changed argument bMerge into oneHB to handle contacts properly.
2247 * - Erik Marklund, June 29, 2006
2249 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
2252 int i, j, k, m, nd, ihb, idist;
2253 int nframes = hb->nframes;
2261 fp = gmx_ffopen("debug-ac.xvg", "w");
2262 for (j = 0; (j < nframes); j++)
2264 fprintf(fp, "%10.3f", hb->time[j]);
2265 for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
2267 for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
2271 hbh = hb->hbmap[i][k];
2276 ihb = is_hb(hbh->h[0], j);
2277 idist = is_hb(hbh->g[0], j);
2283 for (m = 0; (m < hb->maxhydro) && !ihb; m++)
2285 ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m], j));
2286 idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
2288 /* This is not correct! */
2289 /* What isn't correct? -Erik M */
2294 fprintf(fp, " %1d-%1d", ihb, idist);
2304 static real calc_dg(real tau, real temp)
2315 return kbt*log(kbt*tau/PLANCK);
2320 int n0, n1, nparams, ndelta;
2322 real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
2325 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
2326 real kt[], real sigma_ct[], real sigma_nt[],
2327 real sigma_kt[], real *k, real *kp,
2328 real *sigma_k, real *sigma_kp,
2334 real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
2339 for (i = 0; (i < n); i++)
2341 if (t[i] >= fit_start)
2354 tl.sigma_ct = sigma_ct;
2355 tl.sigma_nt = sigma_nt;
2356 tl.sigma_kt = sigma_kt;
2360 chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
2362 *kp = tl.kkk[1] = *kp;
2364 for (j = 0; (j < NK); j++)
2366 /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
2369 kk2 += sqr(tl.kkk[0]);
2370 kp2 += sqr(tl.kkk[1]);
2373 *sigma_k = sqrt(kk2/NK - sqr(kkk/NK));
2374 *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2379 static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start,
2380 const output_env_t oenv)
2388 for (i = 0; (i < n); i++)
2404 do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(),
2405 effnEXP2, fitparm, 0);
2408 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
2409 real sigma_ct[], real sigma_nt[], real sigma_kt[],
2410 real fit_start, real temp, real smooth_tail_start,
2411 const output_env_t oenv)
2414 real k = 1, kp = 1, kow = 1;
2415 real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
2416 double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
2417 gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2419 if (smooth_tail_start >= 0)
2421 smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv);
2422 smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv);
2423 smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv);
2425 for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2431 for (i = i0; (i < n); i++)
2440 printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
2441 tmp = (sn2*sc2-sqr(scn));
2442 if ((tmp > 0) && (sn2 > 0))
2444 k = (sn2*sck-scn*snk)/tmp;
2445 kp = (k*scn-snk)/sn2;
2449 for (i = i0; (i < n); i++)
2451 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2453 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
2455 &sigma_k, &sigma_kp, fit_start);
2456 Q = 0; /* quality_of_fit(chi2, 2);*/
2457 ddg = BOLTZ*temp*sigma_k/k;
2458 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2460 printf("The Rate and Delta G are followed by an error estimate\n");
2461 printf("----------------------------------------------------------\n"
2462 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
2463 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2464 k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
2465 ddg = BOLTZ*temp*sigma_kp/kp;
2466 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2467 kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
2472 for (i = i0; (i < n); i++)
2474 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2476 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2478 printf("--------------------------------------------------\n"
2479 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
2480 printf("Forward %10.3f %8.3f %10.3f %10g\n",
2481 k, 1/k, calc_dg(1/k, temp), chi2);
2482 printf("Backward %10.3f %8.3f %10.3f\n",
2483 kp, 1/kp, calc_dg(1/kp, temp));
2489 printf("One-way %10.3f %s%8.3f %10.3f\n",
2490 kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp));
2494 printf(" - Numerical problems computing HB thermodynamics:\n"
2495 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
2496 sc2, sn2, sk2, sck, snk, scn);
2498 /* Determine integral of the correlation function */
2499 tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
2500 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
2501 bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
2503 for (i = 0; (i < n-2); i++)
2505 if ((ct[i] > e_1) && (ct[i+1] <= e_1))
2512 /* Determine tau_relax from linear interpolation */
2513 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2514 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx,
2515 tau_rlx, bError ? " " : "",
2516 calc_dg(tau_rlx, temp));
2521 printf("Correlation functions too short to compute thermodynamics\n");
2525 void compute_derivative(int nn, real x[], real y[], real dydx[])
2529 /* Compute k(t) = dc(t)/dt */
2530 for (j = 1; (j < nn-1); j++)
2532 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2534 /* Extrapolate endpoints */
2535 dydx[0] = 2*dydx[1] - dydx[2];
2536 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2539 static void parallel_print(int *data, int nThreads)
2541 /* This prints the donors on which each tread is currently working. */
2544 fprintf(stderr, "\r");
2545 for (i = 0; i < nThreads; i++)
2547 fprintf(stderr, "%-7i", data[i]);
2551 static void normalizeACF(real *ct, real *gt, int nhb, int len)
2553 real ct_fac, gt_fac = 0;
2556 /* Xu and Berne use the same normalization constant */
2561 gt_fac = 1.0/(real)nhb;
2564 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
2565 for (i = 0; i < len; i++)
2575 /* Added argument bContact in do_hbac(...). Also
2576 * added support for -contact in the actual code.
2577 * - Erik Marklund, May 31, 2006 */
2578 /* Changed contact code and added argument R2
2579 * - Erik Marklund, June 29, 2006
2581 static void do_hbac(const char *fn, t_hbdata *hb,
2582 int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
2583 real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv,
2584 const char *gemType, int nThreads,
2585 const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2588 int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
2589 const char *legNN[] = {
2593 static char **legGem;
2595 const char *legLuzar[] = {
2596 "Ac\\sfin sys\\v{}\\z{}(t)",
2598 "Cc\\scontact,hb\\v{}\\z{}(t)",
2599 "-dAc\\sfs\\v{}\\z{}/dt"
2601 gmx_bool bNorm = FALSE, bOMP = FALSE;
2604 real *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
2605 real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
2606 const real tol = 1e-3;
2607 int nframes = hb->nframes, nf;
2608 unsigned int **h = NULL, **g = NULL;
2609 int nh, nhbonds, nhydro, ngh;
2611 PSTYPE p, *pfound = NULL, np;
2613 int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
2614 real **rHbExGem = NULL;
2618 double *ctdouble, *timedouble, *fittedct;
2619 double fittolerance = 0.1;
2620 int *dondata = NULL, thisThread;
2623 AC_NONE, AC_NN, AC_GEM, AC_LUZAR
2632 printf("Doing autocorrelation ");
2634 /* Decide what kind of ACF calculations to do. */
2635 if (NN > NN_NONE && NN < NN_NR)
2637 #ifdef HAVE_NN_LOOPS
2639 printf("using the energy estimate.\n");
2642 printf("Can't do the NN-loop. Yet.\n");
2648 printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2650 nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0);
2651 snew(legGem, nSets);
2652 for (i = 0; i < nSets; i++)
2654 snew(legGem[i], 128);
2656 sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2659 sprintf(legGem[1], "Ac'(t)");
2663 sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2670 printf("according to the theory of Luzar and Chandler.\n");
2674 /* build hbexist matrix in reals for autocorr */
2675 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2677 while (n2 < nframes)
2684 if (acType != AC_NN || bOMP)
2686 snew(h, hb->maxhydro);
2687 snew(g, hb->maxhydro);
2690 /* Dump hbonds for debugging */
2691 dump_ac(hb, bMerge || bContact, nDump);
2693 /* Total number of hbonds analyzed here */
2698 if (acType != AC_LUZAR && bOMP)
2700 nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2702 gmx_omp_set_num_threads(nThreads);
2703 snew(dondata, nThreads);
2704 for (i = 0; i < nThreads; i++)
2708 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2709 "Expect close to linear scaling over this donor-loop.\n", nThreads);
2711 fprintf(stderr, "Donors: [thread no]\n");
2714 for (i = 0; i < nThreads; i++)
2716 snprintf(tmpstr, 7, "[%i]", i);
2717 fprintf(stderr, "%-7s", tmpstr);
2720 fprintf(stderr, "\n");
2724 /* Build the ACF according to acType */
2729 #ifdef HAVE_NN_LOOPS
2730 /* Here we're using the estimated energy for the hydrogen bonds. */
2733 #pragma omp parallel \
2734 private(i, j, k, nh, E, rhbex, thisThread) \
2738 thisThread = gmx_omp_get_thread_num();
2742 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2745 #pragma omp for schedule (dynamic)
2746 for (i = 0; i < hb->d.nrd; i++) /* loop over donors */
2750 #pragma omp critical
2752 dondata[thisThread] = i;
2753 parallel_print(dondata, nThreads);
2758 fprintf(stderr, "\r %i", i);
2761 for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */
2763 for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2765 E = hb->hbE.E[i][j][nh];
2768 for (k = 0; k < nframes; k++)
2770 if (E[k] != NONSENSE_E)
2772 rhbex[k] = (real)E[k];
2776 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0],
2777 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
2778 #pragma omp critical
2780 for (k = 0; (k < nn); k++)
2797 normalizeACF(ct, NULL, 0, nn);
2799 snew(timedouble, nn);
2800 for (j = 0; j < nn; j++)
2802 timedouble[j] = (double)(hb->time[j]);
2803 ctdouble[j] = (double)(ct[j]);
2806 /* Remove ballistic term */
2807 /* Ballistic component removal and fitting to the reversible geminate recombination model
2808 * will be taken out for the time being. First of all, one can remove the ballistic
2809 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2810 * problems with the robustness of the fitting to the model. More work is needed.
2811 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2812 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2813 * a BSD-licence that can do the job.
2815 * - Erik Marklund, June 18 2010.
2817 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2818 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2820 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2821 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2823 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)");
2824 xvgr_legend(fp, asize(legNN), legNN);
2826 for (j = 0; (j < nn); j++)
2828 fprintf(fp, "%10g %10g %10g\n",
2829 hb->time[j]-hb->time[0],
2837 #endif /* HAVE_NN_LOOPS */
2838 break; /* case AC_NN */
2842 memset(ct, 0, 2*n2*sizeof(real));
2844 fprintf(stderr, "Donor:\n");
2847 #define __ACDATA p_ct
2850 #pragma omp parallel \
2851 private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2852 pfound, poff, rHbExGem, p, ihb, mMax, \
2855 { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
2858 thisThread = gmx_omp_get_thread_num();
2859 snew(h, hb->maxhydro);
2860 snew(g, hb->maxhydro);
2867 memset(p_ct, 0, 2*n2*sizeof(real));
2869 /* I'm using a chunk size of 1, since I expect \
2870 * the overhead to be really small compared \
2871 * to the actual calculations \ */
2872 #pragma omp for schedule(dynamic,1) nowait
2873 for (i = 0; i < hb->d.nrd; i++)
2878 #pragma omp critical
2880 dondata[thisThread] = i;
2881 parallel_print(dondata, nThreads);
2886 fprintf(stderr, "\r %i", i);
2888 for (k = 0; k < hb->a.nra; k++)
2890 for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++)
2892 hbh = hb->hbmap[i][k];
2895 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
2896 * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
2897 * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
2898 pHist = &(hb->per->pHist[i][k]);
2899 if (ISHB(hbh->history[nh]) && pHist->len != 0)
2904 g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL;
2908 /* count the number of periodic shifts encountered and store
2909 * them in separate arrays. */
2911 for (j = 0; j < pHist->len; j++)
2914 for (m = 0; m <= np; m++)
2916 if (m == np) /* p not recognized in list. Add it and set up new array. */
2919 if (np > hb->per->nper)
2921 gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
2923 if (m >= mMax) /* Extend the arrays.
2924 * Doing it like this, using mMax to keep track of the sizes,
2925 * eleviates the need for freeing and re-allocating the arrays
2926 * when taking on the next donor-acceptor pair */
2929 srenew(pfound, np); /* The list of found periodic shifts. */
2930 srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */
2931 snew(rHbExGem[m], 2*n2);
2936 if (rHbExGem != NULL && rHbExGem[m] != NULL)
2938 /* This must be done, as this array was most likey
2939 * used to store stuff in some previous iteration. */
2940 memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
2944 fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
2956 } /* m: Loop over found shifts */
2957 } /* j: Loop over shifts */
2959 /* Now unpack and disentangle the existence funtions. */
2960 for (j = 0; j < nf; j++)
2967 * pfound: list of periodic shifts found for this pair.
2968 * poff: list of frame offsets; that is, the first
2969 * frame a hbond has a particular periodic shift. */
2970 p = getPshift(*pHist, j+n0);
2973 for (m = 0; m < np; m++)
2981 gmx_fatal(FARGS, "Shift not found, but must be there.");
2985 ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j));
2990 poff[m] = j; /* Here's where the first hbond with shift p is,
2991 * relative to the start of h[0].*/
2995 gmx_fatal(FARGS, "j<poff[m]");
2997 rHbExGem[m][j-poff[m]] += 1;
3002 /* Now, build ac. */
3003 for (m = 0; m < np; m++)
3005 if (rHbExGem[m][0] > 0 && n0+poff[m] < nn /* && m==0 */)
3007 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0],
3008 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3009 for (j = 0; (j < nn); j++)
3011 __ACDATA[j] += rHbExGem[m][j];
3014 } /* Building of ac. */
3017 } /* hydrogen loop */
3018 } /* acceptor loop */
3021 for (m = 0; m <= mMax; m++)
3034 #pragma omp critical
3036 for (i = 0; i < nn; i++)
3044 } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
3050 normalizeACF(ct, NULL, 0, nn);
3052 fprintf(stderr, "\n\nACF successfully calculated.\n");
3054 /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
3057 snew(timedouble, nn);
3060 for (j = 0; j < nn; j++)
3062 timedouble[j] = (double)(hb->time[j]);
3063 ctdouble[j] = (double)(ct[j]);
3066 /* Remove ballistic term */
3067 /* Ballistic component removal and fitting to the reversible geminate recombination model
3068 * will be taken out for the time being. First of all, one can remove the ballistic
3069 * component with g_analyze afterwards. Secondly, and more importantly, there are still
3070 * problems with the robustness of the fitting to the model. More work is needed.
3071 * A third reason is that we're currently using gsl for this and wish to reduce dependence
3072 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
3073 * a BSD-licence that can do the job.
3075 * - Erik Marklund, June 18 2010.
3077 /* if (bBallistic) { */
3078 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
3079 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
3081 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
3082 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
3085 /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
3090 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3094 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3096 xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv);
3098 for (j = 0; (j < nn); j++)
3100 fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]);
3103 fprintf(fp, " %10g", ctdouble[j]);
3107 fprintf(fp, " %10g", fittedct[j]);
3118 break; /* case AC_GEM */
3131 for (i = 0; (i < hb->d.nrd); i++)
3133 for (k = 0; (k < hb->a.nra); k++)
3136 hbh = hb->hbmap[i][k];
3140 if (bMerge || bContact)
3142 if (ISHB(hbh->history[0]))
3151 for (m = 0; (m < hb->maxhydro); m++)
3153 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
3155 g[nhydro] = hbh->g[m];
3156 h[nhydro] = hbh->h[m];
3163 for (nh = 0; (nh < nhydro); nh++)
3165 int nrint = bContact ? hb->nrdist : hb->nrhb;
3166 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
3168 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
3171 for (j = 0; (j < nframes); j++)
3173 /* Changed '<' into '<=' below, just like I did in
3174 the hbm-output-loop in the gmx_hbond() block.
3175 - Erik Marklund, May 31, 2006 */
3178 ihb = is_hb(h[nh], j);
3179 idist = is_hb(g[nh], j);
3186 /* For contacts: if a second cut-off is provided, use it,
3187 * otherwise use g(t) = 1-h(t) */
3188 if (!R2 && bContact)
3194 gt[j] = idist*(1-ihb);
3201 /* The autocorrelation function is normalized after summation only */
3202 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
3203 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
3205 /* Cross correlation analysis for thermodynamics */
3206 for (j = nframes; (j < n2); j++)
3212 cross_corr(n2, ht, gt, dght);
3214 for (j = 0; (j < nn); j++)
3223 fprintf(stderr, "\n");
3226 normalizeACF(ct, ght, nhb, nn);
3228 /* Determine tail value for statistics */
3231 for (j = nn/2; (j < nn); j++)
3234 tail2 += ct[j]*ct[j];
3236 tail /= (nn - nn/2);
3237 tail2 /= (nn - nn/2);
3238 dtail = sqrt(tail2-tail*tail);
3240 /* Check whether the ACF is long enough */
3243 printf("\nWARNING: Correlation function is probably not long enough\n"
3244 "because the standard deviation in the tail of C(t) > %g\n"
3245 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
3248 for (j = 0; (j < nn); j++)
3251 ct[j] = (cct[j]-tail)/(1-tail);
3253 /* Compute negative derivative k(t) = -dc(t)/dt */
3254 compute_derivative(nn, hb->time, ct, kt);
3255 for (j = 0; (j < nn); j++)
3263 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3267 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3269 xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
3272 for (j = 0; (j < nn); j++)
3274 fprintf(fp, "%10g %10g %10g %10g %10g\n",
3275 hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
3279 analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
3280 fit_start, temp, smooth_tail_start, oenv);
3282 do_view(oenv, fn, NULL);
3294 break; /* case AC_LUZAR */
3297 gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
3298 } /* switch (acType) */
3301 static void init_hbframe(t_hbdata *hb, int nframes, real t)
3305 hb->time[nframes] = t;
3306 hb->nhb[nframes] = 0;
3307 hb->ndist[nframes] = 0;
3308 for (i = 0; (i < max_hx); i++)
3310 hb->nhx[nframes][i] = 0;
3312 /* Loop invalidated */
3313 if (hb->bHBmap && 0)
3315 for (i = 0; (i < hb->d.nrd); i++)
3317 for (j = 0; (j < hb->a.nra); j++)
3319 for (m = 0; (m < hb->maxhydro); m++)
3321 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
3323 set_hb(hb, i, m, j, nframes, HB_NO);
3329 /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
3332 static FILE *open_donor_properties_file(const char *fn,
3334 const output_env_t oenv)
3337 const char *leg[] = { "Nbound", "Nfree" };
3344 fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3345 xvgr_legend(fp, asize(leg), leg, oenv);
3350 static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
3352 int i, j, k, nbound, nb, nhtot;
3360 for (i = 0; (i < hb->d.nrd); i++)
3362 for (k = 0; (k < hb->d.nhydro[i]); k++)
3366 for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
3368 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
3369 is_hb(hb->hbmap[i][j]->h[k], nframes))
3377 fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound);
3380 static void dump_hbmap(t_hbdata *hb,
3381 int nfile, t_filenm fnm[], gmx_bool bTwo,
3382 gmx_bool bContact, int isize[], int *index[], char *grpnames[],
3386 int ddd, hhh, aaa, i, j, k, m, grp;
3387 char ds[32], hs[32], as[32];
3390 fp = opt2FILE("-hbn", nfile, fnm, "w");
3391 if (opt2bSet("-g", nfile, fnm))
3393 fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
3394 fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
3400 for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
3402 fprintf(fp, "[ %s ]", grpnames[grp]);
3403 for (i = 0; i < isize[grp]; i++)
3405 fprintf(fp, (i%15) ? " " : "\n");
3406 fprintf(fp, " %4d", index[grp][i]+1);
3410 Added -contact support below.
3411 - Erik Marklund, May 29, 2006
3415 fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
3416 for (i = 0; (i < hb->d.nrd); i++)
3418 if (hb->d.grp[i] == grp)
3420 for (j = 0; (j < hb->d.nhydro[i]); j++)
3422 fprintf(fp, " %4d %4d", hb->d.don[i]+1,
3423 hb->d.hydro[i][j]+1);
3429 fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
3430 for (i = 0; (i < hb->a.nra); i++)
3432 if (hb->a.grp[i] == grp)
3434 fprintf(fp, (i%15 && !first) ? " " : "\n");
3435 fprintf(fp, " %4d", hb->a.acc[i]+1);
3444 fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
3445 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
3449 fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
3452 for (i = 0; (i < hb->d.nrd); i++)
3455 for (k = 0; (k < hb->a.nra); k++)
3458 for (m = 0; (m < hb->d.nhydro[i]); m++)
3460 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
3462 sprintf(ds, "%s", mkatomname(atoms, ddd));
3463 sprintf(as, "%s", mkatomname(atoms, aaa));
3466 fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
3469 fprintf(fplog, "%12s %12s\n", ds, as);
3474 hhh = hb->d.hydro[i][m];
3475 sprintf(hs, "%s", mkatomname(atoms, hhh));
3476 fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
3479 fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
3493 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3494 * It mimics add_frames() and init_frame() to some extent. */
3495 static void sync_hbdata(t_hbdata *p_hb, int nframes)
3498 if (nframes >= p_hb->max_frames)
3500 p_hb->max_frames += 4096;
3501 srenew(p_hb->nhb, p_hb->max_frames);
3502 srenew(p_hb->ndist, p_hb->max_frames);
3503 srenew(p_hb->n_bound, p_hb->max_frames);
3504 srenew(p_hb->nhx, p_hb->max_frames);
3507 srenew(p_hb->danr, p_hb->max_frames);
3509 memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3510 memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3511 p_hb->nhb[nframes] = 0;
3512 p_hb->ndist[nframes] = 0;
3515 p_hb->nframes = nframes;
3518 /* p_hb->nhx[nframes][i] */
3520 memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
3522 /* hb->per will remain constant througout the frame loop,
3523 * even though the data its members point to will change,
3524 * hence no need for re-syncing. */
3527 int gmx_hbond(int argc, char *argv[])
3529 const char *desc[] = {
3530 "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3531 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
3532 "(zero is extended) and the distance Donor - Acceptor",
3533 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
3534 "OH and NH groups are regarded as donors, O is an acceptor always,",
3535 "N is an acceptor by default, but this can be switched using",
3536 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3537 "to the first preceding non-hydrogen atom.[PAR]",
3539 "You need to specify two groups for analysis, which must be either",
3540 "identical or non-overlapping. All hydrogen bonds between the two",
3541 "groups are analyzed.[PAR]",
3543 "If you set [TT]-shell[tt], you will be asked for an additional index group",
3544 "which should contain exactly one atom. In this case, only hydrogen",
3545 "bonds between atoms within the shell distance from the one atom are",
3548 "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3549 "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3550 "If contact kinetics are analyzed by using the -contact option, then",
3551 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3552 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3553 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3554 "See mentioned literature for more details and definitions."
3557 /* "It is also possible to analyse specific hydrogen bonds with",
3558 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3559 "Donor Hydrogen Acceptor, in the following way:[PAR]",
3566 "Note that the triplets need not be on separate lines.",
3567 "Each atom triplet specifies a hydrogen bond to be analyzed,",
3568 "note also that no check is made for the types of atoms.[PAR]",
3571 "[BB]Output:[bb][BR]",
3572 "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
3573 "[TT]-ac[tt]: average over all autocorrelations of the existence",
3574 "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3575 "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3576 "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
3577 "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
3578 "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3579 "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3580 "with helices in proteins.[BR]",
3581 "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
3582 "for selected groups, all hydrogen bonded atoms from all groups and",
3583 "all solvent atoms involved in insertion.[BR]",
3584 "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
3585 "frames, this also contains information on solvent insertion",
3586 "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3588 "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3589 "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3590 "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3591 "compare results to Raman Spectroscopy.",
3593 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3594 "require an amount of memory proportional to the total numbers of donors",
3595 "times the total number of acceptors in the selected group(s)."
3598 static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
3599 static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1;
3600 static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
3601 static int nDump = 0, nFitPoints = 100;
3602 static int nThreads = 0, nBalExp = 4;
3604 static gmx_bool bContact = FALSE, bBallistic = FALSE, bGemFit = FALSE;
3605 static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
3606 static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3610 { "-a", FALSE, etREAL, {&acut},
3611 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
3612 { "-r", FALSE, etREAL, {&rcut},
3613 "Cutoff radius (nm, X - Acceptor, see next option)" },
3614 { "-da", FALSE, etBOOL, {&bDA},
3615 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3616 { "-r2", FALSE, etREAL, {&r2cut},
3617 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3618 { "-abin", FALSE, etREAL, {&abin},
3619 "Binwidth angle distribution (degrees)" },
3620 { "-rbin", FALSE, etREAL, {&rbin},
3621 "Binwidth distance distribution (nm)" },
3622 { "-nitacc", FALSE, etBOOL, {&bNitAcc},
3623 "Regard nitrogen atoms as acceptors" },
3624 { "-contact", FALSE, etBOOL, {&bContact},
3625 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3626 { "-shell", FALSE, etREAL, {&rshell},
3627 "when > 0, only calculate hydrogen bonds within # nm shell around "
3629 { "-fitstart", FALSE, etREAL, {&fit_start},
3630 "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]" },
3631 { "-fitend", FALSE, etREAL, {&fit_end},
3632 "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])" },
3633 { "-temp", FALSE, etREAL, {&temp},
3634 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3635 { "-smooth", FALSE, etREAL, {&smooth_tail_start},
3636 "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3637 { "-dump", FALSE, etINT, {&nDump},
3638 "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
3639 { "-max_hb", FALSE, etREAL, {&maxnhb},
3640 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3641 { "-merge", FALSE, etBOOL, {&bMerge},
3642 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3643 { "-geminate", FALSE, etENUM, {gemType},
3644 "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3645 { "-diff", FALSE, etREAL, {&D},
3646 "HIDDENDffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
3648 { "-nthreads", FALSE, etINT, {&nThreads},
3649 "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)"},
3652 const char *bugs[] = {
3653 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3656 { efTRX, "-f", NULL, ffREAD },
3657 { efTPR, NULL, NULL, ffREAD },
3658 { efNDX, NULL, NULL, ffOPTRD },
3659 /* { efNDX, "-sel", "select", ffOPTRD },*/
3660 { efXVG, "-num", "hbnum", ffWRITE },
3661 { efLOG, "-g", "hbond", ffOPTWR },
3662 { efXVG, "-ac", "hbac", ffOPTWR },
3663 { efXVG, "-dist", "hbdist", ffOPTWR },
3664 { efXVG, "-ang", "hbang", ffOPTWR },
3665 { efXVG, "-hx", "hbhelix", ffOPTWR },
3666 { efNDX, "-hbn", "hbond", ffOPTWR },
3667 { efXPM, "-hbm", "hbmap", ffOPTWR },
3668 { efXVG, "-don", "donor", ffOPTWR },
3669 { efXVG, "-dan", "danum", ffOPTWR },
3670 { efXVG, "-life", "hblife", ffOPTWR },
3671 { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
3674 #define NFILE asize(fnm)
3676 char hbmap [HB_NR] = { ' ', 'o', '-', '*' };
3677 const char *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
3678 t_rgb hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
3680 t_trxstatus *status;
3685 int npargs, natoms, nframes = 0, shatom;
3691 real t, ccut, dist = 0.0, ang = 0.0;
3692 double max_nhb, aver_nhb, aver_dist;
3693 int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
3695 int xj, yj, zj, aj, xjj, yjj, zjj;
3696 int xk, yk, zk, ak, xkk, ykk, zkk;
3697 gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
3698 int *adist, *rdist, *aptr, *rprt;
3699 int grp, nabin, nrbin, bin, resdist, ihb;
3701 t_hbdata *hb, *hbptr;
3702 FILE *fp, *fpins = NULL, *fpnhb = NULL, *donor_properties = NULL;
3704 t_ncell *icell, *jcell, *kcell;
3706 unsigned char *datable;
3711 int ii, jj, hh, actual_nThreads;
3713 gmx_bool bGem, bNN, bParallel;
3714 t_gemParams *params = NULL;
3715 gmx_bool bEdge_yjj, bEdge_xjj, bOMP;
3717 t_hbdata **p_hb = NULL; /* one per thread, then merge after the frame loop */
3718 int **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
3727 ppa = add_acf_pargs(&npargs, pa);
3729 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
3730 ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
3735 /* NN-loop? If so, what estimator to use ?*/
3737 /* Outcommented for now DvdS 2010-07-13
3738 while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3741 gmx_fatal(FARGS, "Invalid NN-loop type.");
3744 for (i = 2; bNN == FALSE && i < NN_NR; i++)
3746 bNN = bNN || NN == i;
3749 if (NN > NN_NONE && bMerge)
3754 /* geminate recombination? If so, which flavor? */
3756 while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0)
3760 if (gemmode == gemNR)
3762 gmx_fatal(FARGS, "Invalid recombination type.");
3766 for (i = 2; bGem == FALSE && i < gemNR; i++)
3768 bGem = bGem || gemmode == i;
3773 printf("Geminate recombination: %s\n", gemType[gemmode]);
3776 if (gemmode != gemDD)
3778 printf("Turning off -contact option...\n");
3784 if (gemmode == gemDD)
3786 printf("Turning on -contact option...\n");
3792 if (gemmode == gemAA)
3794 printf("Turning off -merge option...\n");
3800 if (gemmode != gemAA)
3802 printf("Turning on -merge option...\n");
3810 ccut = cos(acut*DEG2RAD);
3816 gmx_fatal(FARGS, "Can not analyze selected contacts.");
3820 gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
3824 /* Initiate main data structure! */
3825 bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
3826 opt2bSet("-life", NFILE, fnm) ||
3827 opt2bSet("-hbn", NFILE, fnm) ||
3828 opt2bSet("-hbm", NFILE, fnm) ||
3831 if (opt2bSet("-nhbdist", NFILE, fnm))
3833 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3834 fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
3835 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
3836 xvgr_legend(fpnhb, asize(leg), leg, oenv);
3839 hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
3842 read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
3844 snew(grpnames, grNR);
3847 /* Make Donor-Acceptor table */
3848 snew(datable, top.atoms.nr);
3849 gen_datable(index[0], isize[0], datable, top.atoms.nr);
3853 /* analyze selected hydrogen bonds */
3854 printf("Select group with selected atoms:\n");
3855 get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
3856 1, &nsel, index, grpnames);
3859 gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
3860 "and therefore cannot contain triplets of "
3861 "Donor-Hydrogen-Acceptor", grpnames[0]);
3865 for (i = 0; (i < nsel); i += 3)
3867 int dd = index[0][i];
3868 int aa = index[0][i+2];
3869 /* int */ hh = index[0][i+1];
3870 add_dh (&hb->d, dd, hh, i, datable);
3871 add_acc(&hb->a, aa, i);
3872 /* Should this be here ? */
3873 snew(hb->d.dptr, top.atoms.nr);
3874 snew(hb->a.aptr, top.atoms.nr);
3875 add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri);
3877 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
3878 isize[0], grpnames[0]);
3882 /* analyze all hydrogen bonds: get group(s) */
3883 printf("Specify 2 groups to analyze:\n");
3884 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
3885 2, isize, index, grpnames);
3887 /* check if we have two identical or two non-overlapping groups */
3888 bTwo = isize[0] != isize[1];
3889 for (i = 0; (i < isize[0]) && !bTwo; i++)
3891 bTwo = index[0][i] != index[1][i];
3895 printf("Checking for overlap in atoms between %s and %s\n",
3896 grpnames[0], grpnames[1]);
3897 for (i = 0; i < isize[1]; i++)
3899 if (ISINGRP(datable[index[1][i]]))
3901 gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
3902 grpnames[0], grpnames[1]);
3906 printf("Checking for overlap in atoms between %s and %s\n",
3907 grpnames[0],grpnames[1]);
3908 for (i=0; i<isize[0]; i++)
3909 for (j=0; j<isize[1]; j++)
3910 if (index[0][i] == index[1][j])
3911 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3912 grpnames[0],grpnames[1]);
3917 printf("Calculating %s "
3918 "between %s (%d atoms) and %s (%d atoms)\n",
3919 bContact ? "contacts" : "hydrogen bonds",
3920 grpnames[0], isize[0], grpnames[1], isize[1]);
3924 fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
3925 bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
3930 /* search donors and acceptors in groups */
3931 snew(datable, top.atoms.nr);
3932 for (i = 0; (i < grNR); i++)
3934 if ( ((i == gr0) && !bSelected ) ||
3935 ((i == gr1) && bTwo ))
3937 gen_datable(index[i], isize[i], datable, top.atoms.nr);
3940 search_acceptors(&top, isize[i], index[i], &hb->a, i,
3941 bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
3942 search_donors (&top, isize[i], index[i], &hb->d, i,
3943 TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
3947 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
3948 search_donors (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
3952 clear_datable_grp(datable, top.atoms.nr);
3957 printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
3959 snew(donors[gr0D], dons[gr0D].nrd);*/
3961 donor_properties = open_donor_properties_file(opt2fn_null("-don", NFILE, fnm), hb, oenv);
3965 printf("Making hbmap structure...");
3966 /* Generate hbond data structure */
3971 #ifdef HAVE_NN_LOOPS
3980 printf("Making per structure...");
3981 /* Generate hbond data structure */
3988 if (hb->d.nrd + hb->a.nra == 0)
3990 printf("No Donors or Acceptors found\n");
3997 printf("No Donors found\n");
4002 printf("No Acceptors found\n");
4008 gmx_fatal(FARGS, "Nothing to be done");
4017 /* get index group with atom for shell */
4020 printf("Select atom for shell (1 atom):\n");
4021 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4022 1, &shisz, &shidx, &shgrpnm);
4025 printf("group contains %d atoms, should be 1 (one)\n", shisz);
4030 printf("Will calculate hydrogen bonds within a shell "
4031 "of %g nm around atom %i\n", rshell, shatom+1);
4034 /* Analyze trajectory */
4035 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
4036 if (natoms > top.atoms.nr)
4038 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
4039 top.atoms.nr, natoms);
4042 bBox = ir.ePBC != epbcNONE;
4043 grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
4046 snew(adist, nabin+1);
4047 snew(rdist, nrbin+1);
4051 gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
4057 #define __ADIST adist
4058 #define __RDIST rdist
4060 #else /* GMX_OPENMP ================================================== \
4061 * Set up the OpenMP stuff, |
4062 * like the number of threads and such |
4063 * Also start the parallel loop. |
4065 #define __ADIST p_adist[threadNr]
4066 #define __RDIST p_rdist[threadNr]
4067 #define __HBDATA p_hb[threadNr]
4071 bParallel = !bSelected;
4075 actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
4077 gmx_omp_set_num_threads(actual_nThreads);
4078 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
4083 actual_nThreads = 1;
4086 snew(p_hb, actual_nThreads);
4087 snew(p_adist, actual_nThreads);
4088 snew(p_rdist, actual_nThreads);
4089 for (i = 0; i < actual_nThreads; i++)
4092 snew(p_adist[i], nabin+1);
4093 snew(p_rdist[i], nrbin+1);
4095 p_hb[i]->max_frames = 0;
4096 p_hb[i]->nhb = NULL;
4097 p_hb[i]->ndist = NULL;
4098 p_hb[i]->n_bound = NULL;
4099 p_hb[i]->time = NULL;
4100 p_hb[i]->nhx = NULL;
4102 p_hb[i]->bHBmap = hb->bHBmap;
4103 p_hb[i]->bDAnr = hb->bDAnr;
4104 p_hb[i]->bGem = hb->bGem;
4105 p_hb[i]->wordlen = hb->wordlen;
4106 p_hb[i]->nframes = hb->nframes;
4107 p_hb[i]->maxhydro = hb->maxhydro;
4108 p_hb[i]->danr = hb->danr;
4111 p_hb[i]->hbmap = hb->hbmap;
4112 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
4113 p_hb[i]->per = hb->per;
4115 #ifdef HAVE_NN_LOOPS
4116 p_hb[i]->hbE = hb->hbE;
4120 p_hb[i]->nrdist = 0;
4124 /* Make a thread pool here,
4125 * instead of forking anew at every frame. */
4127 #pragma omp parallel \
4129 private(j, h, ii, jj, hh, E, \
4130 xi, yi, zi, xj, yj, zj, threadNr, \
4131 dist, ang, peri, icell, jcell, \
4132 grp, ogrp, ai, aj, xjj, yjj, zjj, \
4133 xk, yk, zk, ihb, id, resdist, \
4134 xkk, ykk, zkk, kcell, ak, k, bTric, \
4135 bEdge_xjj, bEdge_yjj) \
4137 { /* Start of parallel region */
4138 threadNr = gmx_omp_get_thread_num();
4143 bTric = bBox && TRICLINIC(box);
4147 sync_hbdata(p_hb[threadNr], nframes);
4151 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
4152 rshell, ngrid, grid);
4153 reset_nhbonds(&(hb->d));
4155 if (debug && bDebug)
4157 dump_grid(debug, ngrid, grid);
4160 add_frames(hb, nframes);
4161 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
4165 count_da_grid(ngrid, grid, hb->danr[nframes]);
4171 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
4176 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
4177 /* Loop over all atom pairs and estimate interaction energy */
4181 addFramesNN(hb, nframes);
4185 #pragma omp for schedule(dynamic)
4186 for (i = 0; i < hb->d.nrd; i++)
4188 for (j = 0; j < hb->a.nra; j++)
4191 h < (bContact ? 1 : hb->d.nhydro[i]);
4194 if (i == hb->d.nrd || j == hb->a.nra)
4196 gmx_fatal(FARGS, "out of bounds");
4199 /* Get the real atom ids */
4202 hh = hb->d.hydro[i][h];
4204 /* Estimate the energy from the geometry */
4205 E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
4206 /* Store the energy */
4207 storeHbEnergy(hb, i, j, h, E, nframes);
4211 #endif /* HAVE_NN_LOOPS */
4220 /* Do not parallelize this just yet. */
4222 for (ii = 0; (ii < nsel); ii++)
4224 int dd = index[0][i];
4225 int aa = index[0][i+2];
4226 /* int */ hh = index[0][i+1];
4227 ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
4228 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4232 /* add to index if not already there */
4234 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri);
4238 } /* if (bSelected) */
4246 calcBoxProjection(box, hb->per->P);
4249 /* loop over all gridcells (xi,yi,zi) */
4250 /* Removed confusing macro, DvdS 27/12/98 */
4253 /* The outer grid loop will have to do for now. */
4254 #pragma omp for schedule(dynamic)
4255 for (xi = 0; xi < ngrid[XX]; xi++)
4257 for (yi = 0; (yi < ngrid[YY]); yi++)
4259 for (zi = 0; (zi < ngrid[ZZ]); zi++)
4262 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
4263 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
4265 icell = &(grid[zi][yi][xi].d[grp]);
4276 /* loop over all hydrogen atoms from group (grp)
4277 * in this gridcell (icell)
4279 for (ai = 0; (ai < icell->nr); ai++)
4281 i = icell->atoms[ai];
4283 /* loop over all adjacent gridcells (xj,yj,zj) */
4284 for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
4285 zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
4288 zj = grid_mod(zjj, ngrid[ZZ]);
4289 bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
4290 for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
4291 yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
4294 yj = grid_mod(yjj, ngrid[YY]);
4296 (yj == 0) || (yj == ngrid[YY] - 1) ||
4297 (zj == 0) || (zj == ngrid[ZZ] - 1);
4298 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
4299 xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
4302 xj = grid_mod(xjj, ngrid[XX]);
4303 jcell = &(grid[zj][yj][xj].a[ogrp]);
4304 /* loop over acceptor atoms from other group (ogrp)
4305 * in this adjacent gridcell (jcell)
4307 for (aj = 0; (aj < jcell->nr); aj++)
4309 j = jcell->atoms[aj];
4311 /* check if this once was a h-bond */
4313 ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
4314 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4318 /* add to index if not already there */
4320 add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri);
4322 /* make angle and distance distributions */
4323 if (ihb == hbHB && !bContact)
4327 gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
4330 __ADIST[(int)( ang/abin)]++;
4331 __RDIST[(int)(dist/rbin)]++;
4335 if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
4337 gmx_fatal(FARGS, "Invalid donor %d", i);
4339 if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
4341 gmx_fatal(FARGS, "Invalid acceptor %d", j);
4343 resdist = abs(top.atoms.atom[i].resind-
4344 top.atoms.atom[j].resind);
4345 if (resdist >= max_hx)
4349 __HBDATA->nhx[nframes][resdist]++;
4360 } /* for xi,yi,zi */
4363 } /* if (bSelected) {...} else */
4366 /* Better wait for all threads to finnish using x[] before updating it. */
4369 #pragma omp critical
4371 /* Sum up histograms and counts from p_hb[] into hb */
4374 hb->nhb[k] += p_hb[threadNr]->nhb[k];
4375 hb->ndist[k] += p_hb[threadNr]->ndist[k];
4376 for (j = 0; j < max_hx; j++)
4378 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
4383 /* Here are a handful of single constructs
4384 * to share the workload a bit. The most
4385 * important one is of course the last one,
4386 * where there's a potential bottleneck in form
4391 analyse_donor_properties(donor_properties, hb, k, t);
4398 do_nhb_dist(fpnhb, hb, t);
4401 } /* if (bNN) {...} else + */
4405 trrStatus = (read_next_x(oenv, status, &t, x, box));
4415 #pragma omp critical
4417 hb->nrhb += p_hb[threadNr]->nrhb;
4418 hb->nrdist += p_hb[threadNr]->nrdist;
4420 /* Free parallel datastructures */
4421 sfree(p_hb[threadNr]->nhb);
4422 sfree(p_hb[threadNr]->ndist);
4423 sfree(p_hb[threadNr]->nhx);
4426 for (i = 0; i < nabin; i++)
4428 for (j = 0; j < actual_nThreads; j++)
4431 adist[i] += p_adist[j][i];
4435 for (i = 0; i <= nrbin; i++)
4437 for (j = 0; j < actual_nThreads; j++)
4439 rdist[i] += p_rdist[j][i];
4443 sfree(p_adist[threadNr]);
4444 sfree(p_rdist[threadNr]);
4446 } /* End of parallel region */
4453 if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
4455 gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
4458 free_grid(ngrid, &grid);
4462 if (donor_properties)
4464 xvgrclose(donor_properties);
4472 /* Compute maximum possible number of different hbonds */
4479 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
4481 /* Added support for -contact below.
4482 * - Erik Marklund, May 29-31, 2006 */
4483 /* Changed contact code.
4484 * - Erik Marklund, June 29, 2006 */
4489 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
4493 printf("Found %d different %s in trajectory\n"
4494 "Found %d different atom-pairs within %s distance\n",
4495 hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
4496 hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
4498 /*Control the pHist.*/
4502 merge_hb(hb, bTwo, bContact);
4505 if (opt2bSet("-hbn", NFILE, fnm))
4507 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
4510 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
4511 * to make the -hbn and -hmb output match eachother.
4512 * - Erik Marklund, May 30, 2006 */
4515 /* Print out number of hbonds and distances */
4518 fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
4519 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
4521 snew(leg[0], STRLEN);
4522 snew(leg[1], STRLEN);
4523 sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
4524 sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
4525 xvgr_legend(fp, 2, (const char**)leg, oenv);
4529 for (i = 0; (i < nframes); i++)
4531 fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
4532 aver_nhb += hb->nhb[i];
4533 aver_dist += hb->ndist[i];
4536 aver_nhb /= nframes;
4537 aver_dist /= nframes;
4538 /* Print HB distance distribution */
4539 if (opt2bSet("-dist", NFILE, fnm))
4544 for (i = 0; i < nrbin; i++)
4549 fp = xvgropen(opt2fn("-dist", NFILE, fnm),
4550 "Hydrogen Bond Distribution",
4552 "Donor - Acceptor Distance (nm)" :
4553 "Hydrogen - Acceptor Distance (nm)", "", oenv);
4554 for (i = 0; i < nrbin; i++)
4556 fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
4561 /* Print HB angle distribution */
4562 if (opt2bSet("-ang", NFILE, fnm))
4567 for (i = 0; i < nabin; i++)
4572 fp = xvgropen(opt2fn("-ang", NFILE, fnm),
4573 "Hydrogen Bond Distribution",
4574 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
4575 for (i = 0; i < nabin; i++)
4577 fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
4582 /* Print HB in alpha-helix */
4583 if (opt2bSet("-hx", NFILE, fnm))
4585 fp = xvgropen(opt2fn("-hx", NFILE, fnm),
4586 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4587 xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
4588 for (i = 0; i < nframes; i++)
4590 fprintf(fp, "%10g", hb->time[i]);
4591 for (j = 0; j < max_hx; j++)
4593 fprintf(fp, " %6d", hb->nhx[i][j]);
4601 printf("Average number of %s per timeframe %.3f out of %g possible\n",
4602 bContact ? "contacts" : "hbonds",
4603 bContact ? aver_dist : aver_nhb, max_nhb);
4606 /* Do Autocorrelation etc. */
4610 Added support for -contact in ac and hbm calculations below.
4611 - Erik Marklund, May 29, 2006
4615 if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
4617 please_cite(stdout, "Spoel2006b");
4619 if (opt2bSet("-ac", NFILE, fnm))
4621 char *gemstring = NULL;
4625 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4626 gemBallistic, nBalExp);
4629 gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4632 gemstring = gmx_strdup(gemType[hb->per->gemtype]);
4633 do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
4634 bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
4635 gemstring, nThreads, NN, bBallistic, bGemFit);
4637 if (opt2bSet("-life", NFILE, fnm))
4639 do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
4641 if (opt2bSet("-hbm", NFILE, fnm))
4644 int id, ia, hh, x, y;
4645 mat.flags = mat.y0 = 0;
4647 if ((nframes > 0) && (hb->nrhb > 0))
4652 snew(mat.matrix, mat.nx);
4653 for (x = 0; (x < mat.nx); x++)
4655 snew(mat.matrix[x], mat.ny);
4658 for (id = 0; (id < hb->d.nrd); id++)
4660 for (ia = 0; (ia < hb->a.nra); ia++)
4662 for (hh = 0; (hh < hb->maxhydro); hh++)
4664 if (hb->hbmap[id][ia])
4666 if (ISHB(hb->hbmap[id][ia]->history[hh]))
4668 /* Changed '<' into '<=' in the for-statement below.
4669 * It fixed the previously undiscovered bug that caused
4670 * the last occurance of an hbond/contact to not be
4671 * set in mat.matrix. Have a look at any old -hbm-output
4672 * and you will notice that the last column is allways empty.
4673 * - Erik Marklund May 30, 2006
4675 for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
4677 int nn0 = hb->hbmap[id][ia]->n0;
4678 range_check(y, 0, mat.ny);
4679 mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
4687 mat.axis_x = hb->time;
4688 snew(mat.axis_y, mat.ny);
4689 for (j = 0; j < mat.ny; j++)
4693 sprintf(mat.title, bContact ? "Contact Existence Map" :
4694 "Hydrogen Bond Existence Map");
4695 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
4696 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
4697 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4698 mat.bDiscrete = TRUE;
4700 snew(mat.map, mat.nmap);
4701 for (i = 0; i < mat.nmap; i++)
4703 mat.map[i].code.c1 = hbmap[i];
4704 mat.map[i].desc = hbdesc[i];
4705 mat.map[i].rgb = hbrgb[i];
4707 fp = opt2FILE("-hbm", NFILE, fnm, "w");
4708 write_xpm_m(fp, mat);
4710 for (x = 0; x < mat.nx; x++)
4712 sfree(mat.matrix[x]);
4720 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
4727 fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4728 fprintf(stderr, "Freeing pHist for all donors...\n");
4729 for (i = 0; i < hb->d.nrd; i++)
4731 fprintf(stderr, "\r%i", i);
4732 if (hb->per->pHist[i] != NULL)
4734 for (j = 0; j < hb->a.nra; j++)
4736 clearPshift(&(hb->per->pHist[i][j]));
4738 sfree(hb->per->pHist[i]);
4741 sfree(hb->per->pHist);
4742 sfree(hb->per->p2i);
4744 fprintf(stderr, "...done.\n");
4747 #ifdef HAVE_NN_LOOPS
4760 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4762 fp = xvgropen(opt2fn("-dan", NFILE, fnm),
4763 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4764 nleg = (bTwo ? 2 : 1)*2;
4765 snew(legnames, nleg);
4767 for (j = 0; j < grNR; j++)
4769 if (USE_THIS_GROUP(j) )
4771 sprintf(buf, "Donors %s", grpnames[j]);
4772 legnames[i++] = gmx_strdup(buf);
4773 sprintf(buf, "Acceptors %s", grpnames[j]);
4774 legnames[i++] = gmx_strdup(buf);
4779 gmx_incons("number of legend entries");
4781 xvgr_legend(fp, nleg, (const char**)legnames, oenv);
4782 for (i = 0; i < nframes; i++)
4784 fprintf(fp, "%10g", hb->time[i]);
4785 for (j = 0; (j < grNR); j++)
4787 if (USE_THIS_GROUP(j) )
4789 fprintf(fp, " %6d", hb->danr[i][j]);