1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2008, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Groningen Machine for Chemical Simulation
41 /*#define HAVE_NN_LOOPS*/
53 #include "gmx_fatal.h"
66 typedef short int t_E;
69 typedef int t_hx[max_hx];
70 #define NRHXTYPES max_hx
71 const char *hxtypenames[NRHXTYPES] =
72 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
76 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == 0)
78 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == (threadNr))
81 /* -----------------------------------------*/
87 hbNo, hbDist, hbHB, hbNR, hbR2
90 noDA, ACC, DON, DA, INGROUP
93 NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR
96 static const char *grpnames[grNR] = {"0", "1", "I" };
98 static gmx_bool bDebug = FALSE;
103 #define HB_YESINS HB_YES|HB_INS
107 #define ISHB(h) (((h) & 2) == 2)
108 #define ISDIST(h) (((h) & 1) == 1)
109 #define ISDIST2(h) (((h) & 4) == 4)
110 #define ISACC(h) (((h) & 1) == 1)
111 #define ISDON(h) (((h) & 2) == 2)
112 #define ISINGRP(h) (((h) & 4) == 4)
125 typedef int t_icell[grNR];
126 typedef atom_id h_id[MAXHYDRO];
129 int history[MAXHYDRO];
130 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
131 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
133 /* Bitmask array which tells whether a hbond is present
134 * at a given time. Either of these may be NULL
136 int n0; /* First frame a HB was found */
137 int nframes, maxframes; /* Amount of frames in this hbond */
140 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
141 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
142 * acceptor distance is less than the user-specified distance (typically
149 atom_id *acc; /* Atom numbers of the acceptors */
150 int *grp; /* Group index */
151 int *aptr; /* Map atom number to acceptor index */
156 int *don; /* Atom numbers of the donors */
157 int *grp; /* Group index */
158 int *dptr; /* Map atom number to donor index */
159 int *nhydro; /* Number of hydrogens for each donor */
160 h_id *hydro; /* The atom numbers of the hydrogens */
161 h_id *nhbonds; /* The number of HBs per H at current */
164 /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/
168 int len; /* The length of frame and p. */
169 int *frame; /* The frames at which transitio*/
174 /* Periodicity history. Used for the reversible geminate recombination. */
175 t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap:
176 * pHist[d][a]. We can safely assume that the same
177 * periodic shift holds for all hydrogens of a da-pair.
179 * Nowadays it only stores TRANSITIONS, and not the shift at every frame.
180 * That saves a LOT of memory, an hopefully kills a mysterious bug where
181 * pHist gets contaminated. */
183 PSTYPE nper; /* The length of p2i */
184 ivec *p2i; /* Maps integer to periodic shift for a pair.*/
185 matrix P; /* Projection matrix to find the box shifts. */
186 int gemtype; /* enumerated type */
191 int *Etot; /* Total energy for each frame */
192 t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */
196 gmx_bool bHBmap, bDAnr, bGem;
198 /* The following arrays are nframes long */
199 int nframes, max_frames, maxhydro;
205 /* These structures are initialized from the topology at start up */
208 /* This holds a matrix with all possible hydrogen bonds */
214 /* For parallelization reasons this will have to be a pointer.
215 * Otherwise discrepancies may arise between the periodicity data
216 * seen by different threads. */
220 static void clearPshift(t_pShift *pShift)
225 sfree(pShift->frame);
230 static void calcBoxProjection(matrix B, matrix P)
232 const int vp[] = {XX, YY, ZZ};
237 for (i = 0; i < 3; i++)
240 for (j = 0; j < 3; j++)
243 U[m][n] = i == j ? 1 : 0;
247 for (i = 0; i < 3; i++)
250 mvmul(M, U[m], P[m]);
255 static void calcBoxDistance(matrix P, rvec d, ivec ibd)
257 /* returns integer distance in box coordinates.
258 * P is the projection matrix from cartesian coordinates
259 * obtained with calcBoxProjection(). */
263 /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
264 for (i = 0; i < 3; i++)
266 bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5);
268 ibd[XX] = (int)bd[XX];
269 ibd[YY] = (int)bd[YY];
270 ibd[ZZ] = (int)bd[ZZ];
273 /* Changed argument 'bMerge' into 'oneHB' below,
274 * since -contact should cause maxhydro to be 1,
276 * - Erik Marklund May 29, 2006
279 static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap)
281 /* Try to merge hbonds on the fly. That means that if the
282 * acceptor and donor are mergable, then:
283 * 1) store the hb-info so that acceptor id > donor id,
284 * 2) add the periodic shift in pairs, so that [-x,-y,-z] is
285 * stored in per.p2i[] whenever acceptor id < donor id.
286 * Note that [0,0,0] should already be the first element of per.p2i
287 * by the time this function is called. */
289 /* daSwap is TRUE if the donor and acceptor were swapped.
290 * If so, then the negative vector should be used. */
293 if (per->p2i == NULL || per->nper == 0)
295 gmx_fatal(FARGS, "'per' not initialized properly.");
297 for (i = 0; i < per->nper; i++)
299 if (r[XX] == per->p2i[i][XX] &&
300 r[YY] == per->p2i[i][YY] &&
301 r[ZZ] == per->p2i[i][ZZ])
306 /* Not found apparently. Add it to the list! */
307 /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
313 fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
318 srenew(per->p2i, per->nper+2);
320 copy_ivec(r, per->p2i[per->nper]);
323 /* Add the mirror too. It's rather likely that it'll be needed. */
324 per->p2i[per->nper][XX] = -r[XX];
325 per->p2i[per->nper][YY] = -r[YY];
326 per->p2i[per->nper][ZZ] = -r[ZZ];
329 return per->nper - 1 - (daSwap ? 0 : 1);
332 static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode)
337 hb->wordlen = 8*sizeof(unsigned int);
347 hb->maxhydro = MAXHYDRO;
350 hb->per->gemtype = bGem ? gemmode : 0;
355 static void mk_hbmap(t_hbdata *hb, gmx_bool bTwo)
359 snew(hb->hbmap, hb->d.nrd);
360 for (i = 0; (i < hb->d.nrd); i++)
362 snew(hb->hbmap[i], hb->a.nra);
363 if (hb->hbmap[i] == NULL)
365 gmx_fatal(FARGS, "Could not allocate enough memory for hbmap");
367 for (j = 0; (j > hb->a.nra); j++)
369 hb->hbmap[i][j] = NULL;
374 /* Consider redoing pHist so that is only stores transitions between
375 * periodicities and not the periodicity for all frames. This eats heaps of memory. */
376 static void mk_per(t_hbdata *hb)
381 snew(hb->per->pHist, hb->d.nrd);
382 for (i = 0; i < hb->d.nrd; i++)
384 snew(hb->per->pHist[i], hb->a.nra);
385 if (hb->per->pHist[i] == NULL)
387 gmx_fatal(FARGS, "Could not allocate enough memory for per->pHist");
389 for (j = 0; j < hb->a.nra; j++)
391 clearPshift(&(hb->per->pHist[i][j]));
394 /* add the [0,0,0] shift to element 0 of p2i. */
395 snew(hb->per->p2i, 1);
396 clear_ivec(hb->per->p2i[0]);
402 static void mk_hbEmap (t_hbdata *hb, int n0)
407 snew(hb->hbE.E, hb->d.nrd);
408 for (i = 0; i < hb->d.nrd; i++)
410 snew(hb->hbE.E[i], hb->a.nra);
411 for (j = 0; j < hb->a.nra; j++)
413 snew(hb->hbE.E[i][j], MAXHYDRO);
414 for (k = 0; k < MAXHYDRO; k++)
416 hb->hbE.E[i][j][k] = NULL;
423 static void free_hbEmap (t_hbdata *hb)
426 for (i = 0; i < hb->d.nrd; i++)
428 for (j = 0; j < hb->a.nra; j++)
430 for (k = 0; k < MAXHYDRO; k++)
432 sfree(hb->hbE.E[i][j][k]);
434 sfree(hb->hbE.E[i][j]);
442 static void addFramesNN(t_hbdata *hb, int frame)
445 #define DELTAFRAMES_HBE 10
447 int d, a, h, nframes;
449 if (frame >= hb->hbE.nframes)
451 nframes = hb->hbE.nframes + DELTAFRAMES_HBE;
452 srenew(hb->hbE.Etot, nframes);
454 for (d = 0; d < hb->d.nrd; d++)
456 for (a = 0; a < hb->a.nra; a++)
458 for (h = 0; h < hb->d.nhydro[d]; h++)
460 srenew(hb->hbE.E[d][a][h], nframes);
465 hb->hbE.nframes += DELTAFRAMES_HBE;
469 static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst,
470 matrix box, rvec hbox, t_donors *donors)
475 * alpha - angle between dipoles
476 * x[] - atomic positions
477 * EEst - the type of energy estimate (see enum in hbplugin.h)
478 * box - the box vectors \
479 * hbox - half box lengths _These two are only needed for the pbc correction
484 rvec dipole[2], xmol[3], xmean[2];
490 /* Self-interaction */
497 /* This is a simple binary existence function that sets E=1 whenever
498 * the distance between the oxygens is equal too or less than 0.35 nm.
500 rvec_sub(x[d], x[a], dist);
501 pbc_correct_gem(dist, box, hbox);
502 if (norm(dist) <= 0.35)
513 /* Negative potential energy of a dipole.
514 * E = -cos(alpha) * 1/r^3 */
516 copy_rvec(x[d], xmol[0]); /* donor */
517 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
518 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
520 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
521 rvec_inc(xmean[0], xmol[1]);
522 rvec_inc(xmean[0], xmol[2]);
523 for (i = 0; i < 3; i++)
525 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
528 /* Assumes that all acceptors are also donors. */
529 copy_rvec(x[a], xmol[0]); /* acceptor */
530 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
531 copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */
534 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
535 rvec_inc(xmean[1], xmol[1]);
536 rvec_inc(xmean[1], xmol[2]);
537 for (i = 0; i < 3; i++)
539 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
542 rvec_sub(xmean[0], xmean[1], dist);
543 pbc_correct_gem(dist, box, hbox);
546 realE = pow(r, -3.0);
547 E = (t_E)(SCALEFACTOR_E * realE);
551 /* Negative potential energy of a (unpolarizable) dipole.
552 * E = -cos(alpha) * 1/r^3 */
553 clear_rvec(dipole[1]);
554 clear_rvec(dipole[0]);
556 copy_rvec(x[d], xmol[0]); /* donor */
557 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
558 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
560 rvec_inc(dipole[0], xmol[1]);
561 rvec_inc(dipole[0], xmol[2]);
562 for (i = 0; i < 3; i++)
566 rvec_dec(dipole[0], xmol[0]);
568 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
569 rvec_inc(xmean[0], xmol[1]);
570 rvec_inc(xmean[0], xmol[2]);
571 for (i = 0; i < 3; i++)
573 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
576 /* Assumes that all acceptors are also donors. */
577 copy_rvec(x[a], xmol[0]); /* acceptor */
578 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
579 copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */
582 rvec_inc(dipole[1], xmol[1]);
583 rvec_inc(dipole[1], xmol[2]);
584 for (i = 0; i < 3; i++)
588 rvec_dec(dipole[1], xmol[0]);
590 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
591 rvec_inc(xmean[1], xmol[1]);
592 rvec_inc(xmean[1], xmol[2]);
593 for (i = 0; i < 3; i++)
595 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
598 rvec_sub(xmean[0], xmean[1], dist);
599 pbc_correct_gem(dist, box, hbox);
602 double cosalpha = cos_angle(dipole[0], dipole[1]);
603 realE = cosalpha * pow(r, -3.0);
604 E = (t_E)(SCALEFACTOR_E * realE);
608 printf("Can't do that type of energy estimate: %i\n.", EEst);
615 static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame)
617 /* hb - hbond data structure
621 E - estimate of the energy
622 frame - the current frame.
625 /* Store the estimated energy */
631 hb->hbE.E[d][a][h][frame] = E;
635 hb->hbE.Etot[frame] += E;
638 #endif /* HAVE_NN_LOOPS */
641 /* Finds -v[] in the periodicity index */
642 static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
646 for (i = 0; i < nper; i++)
648 if (v[i][XX] == -(v[p][XX]) &&
649 v[i][YY] == -(v[p][YY]) &&
650 v[i][ZZ] == -(v[p][ZZ]))
655 printf("Couldn't find mirror of [%i, %i, %i], index \n",
663 static void add_frames(t_hbdata *hb, int nframes)
667 if (nframes >= hb->max_frames)
669 hb->max_frames += 4096;
670 srenew(hb->time, hb->max_frames);
671 srenew(hb->nhb, hb->max_frames);
672 srenew(hb->ndist, hb->max_frames);
673 srenew(hb->n_bound, hb->max_frames);
674 srenew(hb->nhx, hb->max_frames);
677 srenew(hb->danr, hb->max_frames);
680 hb->nframes = nframes;
683 #define OFFSET(frame) (frame / 32)
684 #define MASK(frame) (1 << (frame % 32))
686 static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
690 hbexist[OFFSET(frame)] |= MASK(frame);
694 hbexist[OFFSET(frame)] &= ~MASK(frame);
698 static gmx_bool is_hb(unsigned int hbexist[], int frame)
700 return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
703 static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
705 unsigned int *ghptr = NULL;
709 ghptr = hb->hbmap[id][ia]->h[ih];
711 else if (ihb == hbDist)
713 ghptr = hb->hbmap[id][ia]->g[ih];
717 gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
720 _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
723 static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
727 snew(pHist->frame, 1);
730 pHist->frame[0] = frame;
735 if (pHist->p[pHist->len-1] != p)
738 srenew(pHist->frame, pHist->len);
739 srenew(pHist->p, pHist->len);
740 pHist->frame[pHist->len-1] = frame;
741 pHist->p[pHist->len-1] = p;
742 } /* Otherwise, there is no transition. */
746 static PSTYPE getPshift(t_pShift pHist, int frame)
751 || (pHist.len > 0 && pHist.frame[0] > frame))
756 for (i = 0; i < pHist.len; i++)
769 /* It seems that frame is after the last periodic transition. Return the last periodicity. */
770 return pHist.p[pHist.len-1];
773 static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPE p)
776 t_hbond *hb = hbd->hbmap[id][ia];
777 int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]);
778 int wlen = hbd->wordlen;
780 gmx_bool bGem = hbd->bGem;
785 hb->maxframes = delta;
786 for (i = 0; (i < maxhydro); i++)
788 snew(hb->h[i], hb->maxframes/wlen);
789 snew(hb->g[i], hb->maxframes/wlen);
794 hb->nframes = frame-hb->n0;
795 /* We need a while loop here because hbonds may be returning
798 while (hb->nframes >= hb->maxframes)
800 n = hb->maxframes + delta;
801 for (i = 0; (i < maxhydro); i++)
803 srenew(hb->h[i], n/wlen);
804 srenew(hb->g[i], n/wlen);
805 for (j = hb->maxframes/wlen; (j < n/wlen); j++)
817 set_hb(hbd, id, h, ia, frame, ihb);
820 if (p >= hbd->per->nper)
822 gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
826 addPshift(&(hbd->per->pHist[id][ia]), p, frame);
834 static void inc_nhbonds(t_donors *ddd, int d, int h)
837 int dptr = ddd->dptr[d];
839 for (j = 0; (j < ddd->nhydro[dptr]); j++)
841 if (ddd->hydro[dptr][j] == h)
843 ddd->nhbonds[dptr][j]++;
847 if (j == ddd->nhydro[dptr])
849 gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
853 static int _acceptor_index(t_acceptors *a, int grp, atom_id i,
854 const char *file, int line)
858 if (a->grp[ai] != grp)
862 fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
863 ai, a->grp[ai], grp, file, line);
872 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
874 static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line)
883 if (d->grp[di] != grp)
887 fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
888 di, d->grp[di], grp, file, line);
897 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
899 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
901 /* g_hbond doesn't allow overlapping groups */
907 donor_index(&hb->d, grpd, a) != NOTSET
908 && acceptor_index(&hb->a, grpa, d) != NOTSET;
912 static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
913 int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPE p)
916 gmx_bool daSwap = FALSE;
918 if ((id = hb->d.dptr[d]) == NOTSET)
920 gmx_fatal(FARGS, "No donor atom %d", d+1);
922 else if (grpd != hb->d.grp[id])
924 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
925 grpd, hb->d.grp[id], d+1);
927 if ((ia = hb->a.aptr[a]) == NOTSET)
929 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
931 else if (grpa != hb->a.grp[ia])
933 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
934 grpa, hb->a.grp[ia], a+1);
940 if (isInterchangable(hb, d, a, grpd, grpa) && d > a)
941 /* Then swap identity so that the id of d is lower then that of a.
943 * This should really be redundant by now, as is_hbond() now ought to return
944 * hbNo in the cases where this conditional is TRUE. */
951 /* Now repeat donor/acc check. */
952 if ((id = hb->d.dptr[d]) == NOTSET)
954 gmx_fatal(FARGS, "No donor atom %d", d+1);
956 else if (grpd != hb->d.grp[id])
958 gmx_fatal(FARGS, "Inconsistent donor groups, %d iso %d, atom %d",
959 grpd, hb->d.grp[id], d+1);
961 if ((ia = hb->a.aptr[a]) == NOTSET)
963 gmx_fatal(FARGS, "No acceptor atom %d", a+1);
965 else if (grpa != hb->a.grp[ia])
967 gmx_fatal(FARGS, "Inconsistent acceptor groups, %d iso %d, atom %d",
968 grpa, hb->a.grp[ia], a+1);
975 /* Loop over hydrogens to find which hydrogen is in this particular HB */
976 if ((ihb == hbHB) && !bMerge && !bContact)
978 for (k = 0; (k < hb->d.nhydro[id]); k++)
980 if (hb->d.hydro[id][k] == h)
985 if (k == hb->d.nhydro[id])
987 gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
1001 if (hb->hbmap[id][ia] == NULL)
1003 snew(hb->hbmap[id][ia], 1);
1004 snew(hb->hbmap[id][ia]->h, hb->maxhydro);
1005 snew(hb->hbmap[id][ia]->g, hb->maxhydro);
1007 add_ff(hb, id, k, ia, frame, ihb, p);
1011 /* Strange construction with frame >=0 is a relic from old code
1012 * for selected hbond analysis. It may be necessary again if that
1013 * is made to work again.
1017 hh = hb->hbmap[id][ia]->history[k];
1023 hb->hbmap[id][ia]->history[k] = hh | 2;
1034 hb->hbmap[id][ia]->history[k] = hh | 1;
1058 if (bMerge && daSwap)
1060 h = hb->d.hydro[id][0];
1062 /* Increment number if HBonds per H */
1063 if (ihb == hbHB && !bContact)
1065 inc_nhbonds(&(hb->d), d, h);
1069 static char *mkatomname(t_atoms *atoms, int i)
1071 static char buf[32];
1074 rnr = atoms->atom[i].resind;
1075 sprintf(buf, "%4s%d%-4s",
1076 *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
1081 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms)
1083 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
1086 for (i = 0; i < isize; i++)
1088 if (index[i] >= natoms)
1090 gmx_fatal(FARGS, "Atom has index %d larger than number of atoms %d.", index[i], natoms);
1092 datable[index[i]] |= INGROUP;
1096 static void clear_datable_grp(unsigned char *datable, int size)
1098 /* Clears group information from the table */
1100 const char mask = !(char)INGROUP;
1103 for (i = 0; i < size; i++)
1110 static void add_acc(t_acceptors *a, int ia, int grp)
1112 if (a->nra >= a->max_nra)
1115 srenew(a->acc, a->max_nra);
1116 srenew(a->grp, a->max_nra);
1118 a->grp[a->nra] = grp;
1119 a->acc[a->nra++] = ia;
1122 static void search_acceptors(t_topology *top, int isize,
1123 atom_id *index, t_acceptors *a, int grp,
1125 gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
1131 for (i = 0; (i < isize); i++)
1135 (((*top->atoms.atomname[n])[0] == 'O') ||
1136 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
1137 ISINGRP(datable[n]))
1139 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
1144 snew(a->aptr, top->atoms.nr);
1145 for (i = 0; (i < top->atoms.nr); i++)
1147 a->aptr[i] = NOTSET;
1149 for (i = 0; (i < a->nra); i++)
1151 a->aptr[a->acc[i]] = i;
1155 static void add_h2d(int id, int ih, t_donors *ddd)
1159 for (i = 0; (i < ddd->nhydro[id]); i++)
1161 if (ddd->hydro[id][i] == ih)
1163 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
1168 if (i == ddd->nhydro[id])
1170 if (ddd->nhydro[id] >= MAXHYDRO)
1172 gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
1173 ddd->don[id], MAXHYDRO);
1175 ddd->hydro[id][i] = ih;
1180 static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable)
1184 if (ISDON(datable[id]) || !datable)
1186 if (ddd->dptr[id] == NOTSET) /* New donor */
1198 if (ddd->nrd >= ddd->max_nrd)
1200 ddd->max_nrd += 128;
1201 srenew(ddd->don, ddd->max_nrd);
1202 srenew(ddd->nhydro, ddd->max_nrd);
1203 srenew(ddd->hydro, ddd->max_nrd);
1204 srenew(ddd->nhbonds, ddd->max_nrd);
1205 srenew(ddd->grp, ddd->max_nrd);
1207 ddd->don[ddd->nrd] = id;
1208 ddd->nhydro[ddd->nrd] = 0;
1209 ddd->grp[ddd->nrd] = grp;
1216 add_h2d(i, ih, ddd);
1221 printf("Warning: Atom %d is not in the d/a-table!\n", id);
1225 static void search_donors(t_topology *top, int isize, atom_id *index,
1226 t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
1227 unsigned char *datable)
1230 t_functype func_type;
1231 t_ilist *interaction;
1232 atom_id nr1, nr2, nr3;
1237 snew(ddd->dptr, top->atoms.nr);
1238 for (i = 0; (i < top->atoms.nr); i++)
1240 ddd->dptr[i] = NOTSET;
1248 for (i = 0; (i < isize); i++)
1250 datable[index[i]] |= DON;
1251 add_dh(ddd, index[i], -1, grp, datable);
1257 for (func_type = 0; (func_type < F_NRE); func_type++)
1259 interaction = &(top->idef.il[func_type]);
1260 if (func_type == F_POSRES || func_type == F_FBPOSRES)
1262 /* The ilist looks strange for posre. Bug in grompp?
1263 * We don't need posre interactions for hbonds anyway.*/
1266 for (i = 0; i < interaction->nr;
1267 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1270 if (func_type != top->idef.functype[interaction->iatoms[i]])
1272 fprintf(stderr, "Error in func_type %s",
1273 interaction_function[func_type].longname);
1277 /* check out this functype */
1278 if (func_type == F_SETTLE)
1280 nr1 = interaction->iatoms[i+1];
1281 nr2 = interaction->iatoms[i+2];
1282 nr3 = interaction->iatoms[i+3];
1284 if (ISINGRP(datable[nr1]))
1286 if (ISINGRP(datable[nr2]))
1288 datable[nr1] |= DON;
1289 add_dh(ddd, nr1, nr1+1, grp, datable);
1291 if (ISINGRP(datable[nr3]))
1293 datable[nr1] |= DON;
1294 add_dh(ddd, nr1, nr1+2, grp, datable);
1298 else if (IS_CHEMBOND(func_type))
1300 for (j = 0; j < 2; j++)
1302 nr1 = interaction->iatoms[i+1+j];
1303 nr2 = interaction->iatoms[i+2-j];
1304 if ((*top->atoms.atomname[nr1][0] == 'H') &&
1305 ((*top->atoms.atomname[nr2][0] == 'O') ||
1306 (*top->atoms.atomname[nr2][0] == 'N')) &&
1307 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1309 datable[nr2] |= DON;
1310 add_dh(ddd, nr2, nr1, grp, datable);
1317 for (func_type = 0; func_type < F_NRE; func_type++)
1319 interaction = &top->idef.il[func_type];
1320 for (i = 0; i < interaction->nr;
1321 i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
1324 if (func_type != top->idef.functype[interaction->iatoms[i]])
1326 gmx_incons("function type in search_donors");
1329 if (interaction_function[func_type].flags & IF_VSITE)
1331 nr1 = interaction->iatoms[i+1];
1332 if (*top->atoms.atomname[nr1][0] == 'H')
1336 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1347 if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1348 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1349 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
1351 datable[nr2] |= DON;
1352 add_dh(ddd, nr2, nr1, grp, datable);
1362 static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
1369 for (i = 0; i < DIM; i++)
1371 ngrid[i] = (box[i][i]/(1.2*rcut));
1375 if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
1377 for (i = 0; i < DIM; i++)
1384 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1385 ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
1387 snew(grid, ngrid[ZZ]);
1388 for (z = 0; z < ngrid[ZZ]; z++)
1390 snew((grid)[z], ngrid[YY]);
1391 for (y = 0; y < ngrid[YY]; y++)
1393 snew((grid)[z][y], ngrid[XX]);
1399 static void reset_nhbonds(t_donors *ddd)
1403 for (i = 0; (i < ddd->nrd); i++)
1405 for (j = 0; (j < MAXHH); j++)
1407 ddd->nhbonds[i][j] = 0;
1412 void pbc_correct_gem(rvec dx, matrix box, rvec hbox);
1414 static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
1415 gmx_bool bBox, matrix box, rvec hbox,
1416 real rcut, real rshell,
1417 ivec ngrid, t_gridcell ***grid)
1419 int i, m, gr, xi, yi, zi, nr;
1422 rvec invdelta, dshell, xtemp = {0, 0, 0};
1424 gmx_bool bDoRshell, bInShell, bAcc;
1429 bDoRshell = (rshell > 0);
1430 rshell2 = sqr(rshell);
1433 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
1435 for (m = 0; m < DIM; m++)
1437 hbox[m] = box[m][m]*0.5;
1440 invdelta[m] = ngrid[m]/box[m][m];
1441 if (1/invdelta[m] < rcut)
1443 gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
1444 "%s can not handle this situation, sorry.\n",
1457 /* resetting atom counts */
1458 for (gr = 0; (gr < grNR); gr++)
1460 for (zi = 0; zi < ngrid[ZZ]; zi++)
1462 for (yi = 0; yi < ngrid[YY]; yi++)
1464 for (xi = 0; xi < ngrid[XX]; xi++)
1466 grid[zi][yi][xi].d[gr].nr = 0;
1467 grid[zi][yi][xi].a[gr].nr = 0;
1473 /* put atoms in grid cells */
1474 for (bAcc = FALSE; (bAcc <= TRUE); bAcc++)
1487 for (i = 0; (i < nr); i++)
1489 /* check if we are inside the shell */
1490 /* if bDoRshell=FALSE then bInShell=TRUE always */
1495 rvec_sub(x[ad[i]], xshell, dshell);
1498 if (FALSE && !hb->bGem)
1500 for (m = DIM-1; m >= 0 && bInShell; m--)
1502 if (dshell[m] < -hbox[m])
1504 rvec_inc(dshell, box[m]);
1506 else if (dshell[m] >= hbox[m])
1508 dshell[m] -= 2*hbox[m];
1510 /* if we're outside the cube, we're outside the sphere also! */
1511 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1519 gmx_bool bDone = FALSE;
1523 for (m = DIM-1; m >= 0 && bInShell; m--)
1525 if (dshell[m] < -hbox[m])
1528 rvec_inc(dshell, box[m]);
1530 if (dshell[m] >= hbox[m])
1533 dshell[m] -= 2*hbox[m];
1537 for (m = DIM-1; m >= 0 && bInShell; m--)
1539 /* if we're outside the cube, we're outside the sphere also! */
1540 if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
1547 /* if we're inside the cube, check if we're inside the sphere */
1550 bInShell = norm2(dshell) < rshell2;
1560 copy_rvec(x[ad[i]], xtemp);
1562 pbc_correct_gem(x[ad[i]], box, hbox);
1564 for (m = DIM-1; m >= 0; m--)
1566 if (TRUE || !hb->bGem)
1568 /* put atom in the box */
1569 while (x[ad[i]][m] < 0)
1571 rvec_inc(x[ad[i]], box[m]);
1573 while (x[ad[i]][m] >= box[m][m])
1575 rvec_dec(x[ad[i]], box[m]);
1578 /* determine grid index of atom */
1579 grididx[m] = x[ad[i]][m]*invdelta[m];
1580 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1584 copy_rvec(xtemp, x[ad[i]]); /* copy back */
1589 range_check(gx, 0, ngrid[XX]);
1590 range_check(gy, 0, ngrid[YY]);
1591 range_check(gz, 0, ngrid[ZZ]);
1595 /* add atom to grid cell */
1598 newgrid = &(grid[gz][gy][gx].a[gr]);
1602 newgrid = &(grid[gz][gy][gx].d[gr]);
1604 if (newgrid->nr >= newgrid->maxnr)
1606 newgrid->maxnr += 10;
1607 DBB(newgrid->maxnr);
1608 srenew(newgrid->atoms, newgrid->maxnr);
1611 newgrid->atoms[newgrid->nr] = ad[i];
1619 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1623 for (gr = 0; (gr < grNR); gr++)
1626 for (zi = 0; zi < ngrid[ZZ]; zi++)
1628 for (yi = 0; yi < ngrid[YY]; yi++)
1630 for (xi = 0; xi < ngrid[XX]; xi++)
1632 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1640 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1641 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1642 * With a triclinic box all loops are 3 long, except when a cell is
1643 * located next to one of the box edges which is not parallel to the
1644 * x/y-plane, in that case all cells in a line or layer are searched.
1645 * This could be implemented slightly more efficient, but the code
1646 * would get much more complicated.
1648 static inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1650 return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
1652 static inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1654 return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
1656 static inline int grid_mod(int j, int n)
1661 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1663 int gr, x, y, z, sum[grNR];
1665 fprintf(fp, "grid %dx%dx%d\n", ngrid[XX], ngrid[YY], ngrid[ZZ]);
1666 for (gr = 0; gr < grNR; gr++)
1669 fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
1670 for (z = 0; z < ngrid[ZZ]; z += 2)
1672 fprintf(fp, "Z=%d,%d\n", z, z+1);
1673 for (y = 0; y < ngrid[YY]; y++)
1675 for (x = 0; x < ngrid[XX]; x++)
1677 fprintf(fp, "%3d", grid[x][y][z].d[gr].nr);
1678 sum[gr] += grid[z][y][x].d[gr].nr;
1679 fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
1680 sum[gr] += grid[z][y][x].a[gr].nr;
1684 if ( (z+1) < ngrid[ZZ])
1686 for (x = 0; x < ngrid[XX]; x++)
1688 fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
1689 sum[gr] += grid[z+1][y][x].d[gr].nr;
1690 fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
1691 sum[gr] += grid[z+1][y][x].a[gr].nr;
1698 fprintf(fp, "TOTALS:");
1699 for (gr = 0; gr < grNR; gr++)
1701 fprintf(fp, " %d=%d", gr, sum[gr]);
1706 /* New GMX record! 5 * in a row. Congratulations!
1707 * Sorry, only four left.
1709 static void free_grid(ivec ngrid, t_gridcell ****grid)
1712 t_gridcell ***g = *grid;
1714 for (z = 0; z < ngrid[ZZ]; z++)
1716 for (y = 0; y < ngrid[YY]; y++)
1726 void pbc_correct_gem(rvec dx, matrix box, rvec hbox)
1729 gmx_bool bDone = FALSE;
1733 for (m = DIM-1; m >= 0; m--)
1735 if (dx[m] < -hbox[m])
1738 rvec_inc(dx, box[m]);
1740 if (dx[m] >= hbox[m])
1743 rvec_dec(dx, box[m]);
1749 /* Added argument r2cut, changed contact and implemented
1750 * use of second cut-off.
1751 * - Erik Marklund, June 29, 2006
1753 static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
1754 real rcut, real r2cut, real ccut,
1755 rvec x[], gmx_bool bBox, matrix box, rvec hbox,
1756 real *d_ha, real *ang, gmx_bool bDA, int *hhh,
1757 gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1759 int h, hh, id, ja, ihb;
1760 rvec r_da, r_ha, r_dh, r = {0, 0, 0};
1762 real rc2, r2c2, rda2, rha2, ca;
1763 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1764 gmx_bool daSwap = FALSE;
1771 if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
1772 ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
1780 rvec_sub(x[d], x[a], r_da);
1781 /* Insert projection code here */
1783 if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa))
1785 /* Then this hbond/contact will be found again, or it has already been found. */
1790 if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
1791 { /* return hbNo; */
1792 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1796 copy_rvec(r_da, r); /* Save this for later */
1797 pbc_correct_gem(r_da, box, hbox);
1801 pbc_correct_gem(r_da, box, hbox);
1804 rda2 = iprod(r_da, r_da);
1808 if (daSwap && grpa == grpd)
1816 calcBoxDistance(hb->per->P, r, ri);
1817 *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */
1821 else if (rda2 < r2c2)
1832 if (bDA && (rda2 > rc2))
1837 for (h = 0; (h < hb->d.nhydro[id]); h++)
1839 hh = hb->d.hydro[id][h];
1843 rvec_sub(x[hh], x[a], r_ha);
1846 pbc_correct_gem(r_ha, box, hbox);
1848 rha2 = iprod(r_ha, r_ha);
1853 calcBoxDistance(hb->per->P, r, ri);
1854 *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */
1857 if (bDA || (!bDA && (rha2 <= rc2)))
1859 rvec_sub(x[d], x[hh], r_dh);
1862 pbc_correct_gem(r_dh, box, hbox);
1869 ca = cos_angle(r_dh, r_da);
1870 /* if angle is smaller, cos is larger */
1874 *d_ha = sqrt(bDA ? rda2 : rha2);
1880 if (bDA || (!bDA && HAinrange))
1890 /* Fixed previously undiscovered bug in the merge
1891 code, where the last frame of each hbond disappears.
1892 - Erik Marklund, June 1, 2006 */
1893 /* Added the following arguments:
1894 * ptmp[] - temporary periodicity hisory
1895 * a1 - identity of first acceptor/donor
1896 * a2 - identity of second acceptor/donor
1897 * - Erik Marklund, FEB 20 2010 */
1899 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1900 * Will do some more testing before removing the function entirely.
1901 * - Erik Marklund, MAY 10 2010 */
1902 static void do_merge(t_hbdata *hb, int ntmp,
1903 unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[],
1904 t_hbond *hb0, t_hbond *hb1, int a1, int a2)
1906 /* Here we need to make sure we're treating periodicity in
1907 * the right way for the geminate recombination kinetics. */
1909 int m, mm, n00, n01, nn0, nnframes;
1913 /* Decide where to start from when merging */
1916 nn0 = min(n00, n01);
1917 nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
1918 /* Initiate tmp arrays */
1919 for (m = 0; (m < ntmp); m++)
1925 /* Fill tmp arrays with values due to first HB */
1926 /* Once again '<' had to be replaced with '<='
1927 to catch the last frame in which the hbond
1929 - Erik Marklund, June 1, 2006 */
1930 for (m = 0; (m <= hb0->nframes); m++)
1933 htmp[mm] = is_hb(hb0->h[0], m);
1936 pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1937 if (pm > hb->per->nper)
1939 gmx_fatal(FARGS, "Illegal shift!");
1943 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1947 /* If we're doing geminate recompbination we usually don't need the distances.
1948 * Let's save some memory and time. */
1949 if (TRUE || !hb->bGem || hb->per->gemtype == gemAD)
1951 for (m = 0; (m <= hb0->nframes); m++)
1954 gtmp[mm] = is_hb(hb0->g[0], m);
1958 for (m = 0; (m <= hb1->nframes); m++)
1961 htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
1962 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
1963 if (hb->bGem /* && ptmp[mm] != 0 */)
1966 /* If this hbond has been seen before with donor and acceptor swapped,
1967 * then we need to find the mirrored (*-1) periodicity vector to truely
1968 * merge the hbond history. */
1969 pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper);
1970 /* Store index of mirror */
1971 if (pm > hb->per->nper)
1973 gmx_fatal(FARGS, "Illegal shift!");
1978 /* Reallocate target array */
1979 if (nnframes > hb0->maxframes)
1981 srenew(hb0->h[0], 4+nnframes/hb->wordlen);
1982 srenew(hb0->g[0], 4+nnframes/hb->wordlen);
1984 if (NULL != hb->per->pHist)
1986 clearPshift(&(hb->per->pHist[a1][a2]));
1989 /* Copy temp array to target array */
1990 for (m = 0; (m <= nnframes); m++)
1992 _set_hb(hb0->h[0], m, htmp[m]);
1993 _set_hb(hb0->g[0], m, gtmp[m]);
1996 addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
2000 /* Set scalar variables */
2002 hb0->maxframes = nnframes;
2005 /* Added argument bContact for nicer output.
2006 * Erik Marklund, June 29, 2006
2008 static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
2010 int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp;
2011 unsigned int *htmp, *gtmp;
2016 indnew = hb->nrdist;
2018 /* Check whether donors are also acceptors */
2019 printf("Merging hbonds with Acceptor and Donor swapped\n");
2021 ntmp = 2*hb->max_frames;
2025 for (i = 0; (i < hb->d.nrd); i++)
2027 fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
2029 ii = hb->a.aptr[id];
2030 for (j = 0; (j < hb->a.nra); j++)
2033 jj = hb->d.dptr[ia];
2034 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
2035 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j]))))
2037 hb0 = hb->hbmap[i][j];
2038 hb1 = hb->hbmap[jj][ii];
2039 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0]))
2041 do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j);
2042 if (ISHB(hb1->history[0]))
2046 else if (ISDIST(hb1->history[0]))
2053 gmx_incons("No contact history");
2057 gmx_incons("Neither hydrogen bond nor distance");
2063 clearPshift(&(hb->per->pHist[jj][ii]));
2067 hb1->history[0] = hbNo;
2072 fprintf(stderr, "\n");
2073 printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew);
2074 printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew);
2076 hb->nrdist = indnew;
2082 static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
2084 int i, j, k, n_bound[MAXHH], nbtot;
2088 /* Set array to 0 */
2089 for (k = 0; (k < MAXHH); k++)
2093 /* Loop over possible donors */
2094 for (i = 0; (i < hb->d.nrd); i++)
2096 for (j = 0; (j < hb->d.nhydro[i]); j++)
2098 n_bound[hb->d.nhbonds[i][j]]++;
2101 fprintf(fp, "%12.5e", t);
2103 for (k = 0; (k < MAXHH); k++)
2105 fprintf(fp, " %8d", n_bound[k]);
2106 nbtot += n_bound[k]*k;
2108 fprintf(fp, " %8d\n", nbtot);
2111 /* Added argument bContact in do_hblife(...). Also
2112 * added support for -contact in function body.
2113 * - Erik Marklund, May 31, 2006 */
2114 /* Changed the contact code slightly.
2115 * - Erik Marklund, June 29, 2006
2117 static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
2118 const output_env_t oenv)
2121 const char *leg[] = { "p(t)", "t p(t)" };
2123 int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
2124 int nframes = hb->nframes;
2127 double sum, integral;
2130 snew(h, hb->maxhydro);
2131 snew(histo, nframes+1);
2132 /* Total number of hbonds analyzed here */
2133 for (i = 0; (i < hb->d.nrd); i++)
2135 for (k = 0; (k < hb->a.nra); k++)
2137 hbh = hb->hbmap[i][k];
2155 for (m = 0; (m < hb->maxhydro); m++)
2159 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
2163 for (nh = 0; (nh < nhydro); nh++)
2168 /* Changed '<' into '<=' below, just like I
2169 did in the hbm-output-loop in the main code.
2170 - Erik Marklund, May 31, 2006
2172 for (j = 0; (j <= hbh->nframes); j++)
2174 ihb = is_hb(h[nh], j);
2175 if (debug && (ndump < 10))
2177 fprintf(debug, "%5d %5d\n", j, ihb);
2197 fprintf(stderr, "\n");
2200 fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv);
2204 fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
2208 xvgr_legend(fp, asize(leg), leg, oenv);
2210 while ((j0 > 0) && (histo[j0] == 0))
2215 for (i = 0; (i <= j0); i++)
2219 dt = hb->time[1]-hb->time[0];
2222 for (i = 1; (i <= j0); i++)
2224 t = hb->time[i] - hb->time[0] - 0.5*dt;
2225 x1 = t*histo[i]/sum;
2226 fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1);
2231 printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
2232 printf("Note that the lifetime obtained in this manner is close to useless\n");
2233 printf("Use the -ac option instead and check the Forward lifetime\n");
2234 please_cite(stdout, "Spoel2006b");
2239 /* Changed argument bMerge into oneHB to handle contacts properly.
2240 * - Erik Marklund, June 29, 2006
2242 static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
2245 int i, j, k, m, nd, ihb, idist;
2246 int nframes = hb->nframes;
2254 fp = ffopen("debug-ac.xvg", "w");
2255 for (j = 0; (j < nframes); j++)
2257 fprintf(fp, "%10.3f", hb->time[j]);
2258 for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++)
2260 for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
2264 hbh = hb->hbmap[i][k];
2269 ihb = is_hb(hbh->h[0], j);
2270 idist = is_hb(hbh->g[0], j);
2276 for (m = 0; (m < hb->maxhydro) && !ihb; m++)
2278 ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m], j));
2279 idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j));
2281 /* This is not correct! */
2282 /* What isn't correct? -Erik M */
2287 fprintf(fp, " %1d-%1d", ihb, idist);
2297 static real calc_dg(real tau, real temp)
2308 return kbt*log(kbt*tau/PLANCK);
2313 int n0, n1, nparams, ndelta;
2315 real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
2319 #include <gsl/gsl_multimin.h>
2320 #include <gsl/gsl_sf.h>
2321 #include <gsl/gsl_version.h>
2323 static double my_f(const gsl_vector *v, void *params)
2325 t_luzar *tl = (t_luzar *)params;
2327 double tol = 1e-16, chi2 = 0;
2331 for (i = 0; (i < tl->nparams); i++)
2333 tl->kkk[i] = gsl_vector_get(v, i);
2338 for (i = tl->n0; (i < tl->n1); i += tl->ndelta)
2340 di = sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
2344 chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
2349 fprintf(stderr, "WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
2350 "di = %g k = %g kp = %g\n", tl->sigma_ct[i],
2351 tl->sigma_nt[i], tl->sigma_kt[i], di, k, kp);
2355 chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
2360 static real optimize_luzar_parameters(FILE *fp, t_luzar *tl, int maxiter,
2368 const gsl_multimin_fminimizer_type *T;
2369 gsl_multimin_fminimizer *s;
2372 gsl_multimin_function my_func;
2375 my_func.n = tl->nparams;
2376 my_func.params = (void *) tl;
2378 /* Starting point */
2379 x = gsl_vector_alloc (my_func.n);
2380 for (i = 0; (i < my_func.n); i++)
2382 gsl_vector_set (x, i, tl->kkk[i]);
2385 /* Step size, different for each of the parameters */
2386 dx = gsl_vector_alloc (my_func.n);
2387 for (i = 0; (i < my_func.n); i++)
2389 gsl_vector_set (dx, i, 0.01*tl->kkk[i]);
2392 T = gsl_multimin_fminimizer_nmsimplex;
2393 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
2395 gsl_multimin_fminimizer_set (s, &my_func, x, dx);
2396 gsl_vector_free (x);
2397 gsl_vector_free (dx);
2401 fprintf(fp, "%5s %12s %12s %12s %12s\n", "Iter", "k", "kp", "NM Size", "Chi2");
2407 status = gsl_multimin_fminimizer_iterate (s);
2411 gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s",
2412 gsl_multimin_fminimizer_name(s));
2415 d2 = gsl_multimin_fminimizer_minimum(s);
2416 size = gsl_multimin_fminimizer_size(s);
2417 status = gsl_multimin_test_size(size, tol);
2419 if (status == GSL_SUCCESS)
2423 fprintf(fp, "Minimum found using %s at:\n",
2424 gsl_multimin_fminimizer_name(s));
2430 fprintf(fp, "%5d", iter);
2431 for (i = 0; (i < my_func.n); i++)
2433 fprintf(fp, " %12.4e", gsl_vector_get (s->x, i));
2435 fprintf (fp, " %12.4e %12.4e\n", size, d2);
2438 while ((status == GSL_CONTINUE) && (iter < maxiter));
2440 gsl_multimin_fminimizer_free (s);
2445 static real quality_of_fit(real chi2, int N)
2447 return gsl_sf_gamma_inc_Q((N-2)/2.0, chi2/2.0);
2451 static real optimize_luzar_parameters(FILE *fp, t_luzar *tl, int maxiter,
2454 fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2458 static real quality_of_fit(real chi2, int N)
2460 fprintf(stderr, "This program needs the GNU scientific library to work.\n");
2467 static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
2468 real kt[], real sigma_ct[], real sigma_nt[],
2469 real sigma_kt[], real *k, real *kp,
2470 real *sigma_k, real *sigma_kp,
2476 real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
2481 for (i = 0; (i < n); i++)
2483 if (t[i] >= fit_start)
2496 tl.sigma_ct = sigma_ct;
2497 tl.sigma_nt = sigma_nt;
2498 tl.sigma_kt = sigma_kt;
2502 chi2 = optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2504 *kp = tl.kkk[1] = *kp;
2506 for (j = 0; (j < NK); j++)
2508 (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
2511 kk2 += sqr(tl.kkk[0]);
2512 kp2 += sqr(tl.kkk[1]);
2515 *sigma_k = sqrt(kk2/NK - sqr(kkk/NK));
2516 *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2521 static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start,
2522 const output_env_t oenv)
2525 real e_1, fitparm[4];
2529 for (i = 0; (i < n); i++)
2545 do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(), effnEXP2, fitparm, 0);
2548 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
2549 real sigma_ct[], real sigma_nt[], real sigma_kt[],
2550 real fit_start, real temp, real smooth_tail_start,
2551 const output_env_t oenv)
2554 real k = 1, kp = 1, kow = 1;
2555 real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
2556 double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
2557 gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2559 if (smooth_tail_start >= 0)
2561 smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv);
2562 smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv);
2563 smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv);
2565 for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2571 for (i = i0; (i < n); i++)
2580 printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
2581 tmp = (sn2*sc2-sqr(scn));
2582 if ((tmp > 0) && (sn2 > 0))
2584 k = (sn2*sck-scn*snk)/tmp;
2585 kp = (k*scn-snk)/sn2;
2589 for (i = i0; (i < n); i++)
2591 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2593 chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
2595 &sigma_k, &sigma_kp, fit_start);
2596 Q = quality_of_fit(chi2, 2);
2597 ddg = BOLTZ*temp*sigma_k/k;
2598 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2600 printf("The Rate and Delta G are followed by an error estimate\n");
2601 printf("----------------------------------------------------------\n"
2602 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
2603 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2604 k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
2605 ddg = BOLTZ*temp*sigma_kp/kp;
2606 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2607 kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
2612 for (i = i0; (i < n); i++)
2614 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2616 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2618 printf("--------------------------------------------------\n"
2619 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
2620 printf("Forward %10.3f %8.3f %10.3f %10g\n",
2621 k, 1/k, calc_dg(1/k, temp), chi2);
2622 printf("Backward %10.3f %8.3f %10.3f\n",
2623 kp, 1/kp, calc_dg(1/kp, temp));
2629 printf("One-way %10.3f %s%8.3f %10.3f\n",
2630 kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp));
2634 printf(" - Numerical problems computing HB thermodynamics:\n"
2635 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
2636 sc2, sn2, sk2, sck, snk, scn);
2638 /* Determine integral of the correlation function */
2639 tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
2640 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
2641 bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
2643 for (i = 0; (i < n-2); i++)
2645 if ((ct[i] > e_1) && (ct[i+1] <= e_1))
2652 /* Determine tau_relax from linear interpolation */
2653 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2654 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx,
2655 tau_rlx, bError ? " " : "",
2656 calc_dg(tau_rlx, temp));
2661 printf("Correlation functions too short to compute thermodynamics\n");
2665 void compute_derivative(int nn, real x[], real y[], real dydx[])
2669 /* Compute k(t) = dc(t)/dt */
2670 for (j = 1; (j < nn-1); j++)
2672 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2674 /* Extrapolate endpoints */
2675 dydx[0] = 2*dydx[1] - dydx[2];
2676 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2679 static void parallel_print(int *data, int nThreads)
2681 /* This prints the donors on which each tread is currently working. */
2684 fprintf(stderr, "\r");
2685 for (i = 0; i < nThreads; i++)
2687 fprintf(stderr, "%-7i", data[i]);
2691 static void normalizeACF(real *ct, real *gt, int nhb, int len)
2693 real ct_fac, gt_fac;
2696 /* Xu and Berne use the same normalization constant */
2699 gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
2701 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
2702 for (i = 0; i < len; i++)
2712 /* Added argument bContact in do_hbac(...). Also
2713 * added support for -contact in the actual code.
2714 * - Erik Marklund, May 31, 2006 */
2715 /* Changed contact code and added argument R2
2716 * - Erik Marklund, June 29, 2006
2718 static void do_hbac(const char *fn, t_hbdata *hb,
2719 int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
2720 real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv,
2721 t_gemParams *params, const char *gemType, int nThreads,
2722 const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2725 int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
2726 const char *legNN[] = {
2730 static char **legGem;
2732 const char *legLuzar[] = {
2733 "Ac\\sfin sys\\v{}\\z{}(t)",
2735 "Cc\\scontact,hb\\v{}\\z{}(t)",
2736 "-dAc\\sfs\\v{}\\z{}/dt"
2738 gmx_bool bNorm = FALSE, bOMP = FALSE;
2741 real *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
2742 real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
2743 const real tol = 1e-3;
2744 int nframes = hb->nframes, nf;
2745 unsigned int **h = NULL, **g = NULL;
2746 int nh, nhbonds, nhydro, ngh;
2748 PSTYPE p, *pfound = NULL, np;
2750 int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
2751 real **rHbExGem = NULL;
2755 double *ctdouble, *timedouble, *fittedct;
2756 double fittolerance = 0.1;
2757 int *dondata = NULL, thisThread;
2760 AC_NONE, AC_NN, AC_GEM, AC_LUZAR
2769 printf("Doing autocorrelation ");
2771 /* Decide what kind of ACF calculations to do. */
2772 if (NN > NN_NONE && NN < NN_NR)
2774 #ifdef HAVE_NN_LOOPS
2776 printf("using the energy estimate.\n");
2779 printf("Can't do the NN-loop. Yet.\n");
2785 printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2787 nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0);
2788 snew(legGem, nSets);
2789 for (i = 0; i < nSets; i++)
2791 snew(legGem[i], 128);
2793 sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2796 sprintf(legGem[1], "Ac'(t)");
2800 sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2807 printf("according to the theory of Luzar and Chandler.\n");
2811 /* build hbexist matrix in reals for autocorr */
2812 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2814 while (n2 < nframes)
2821 if (acType != AC_NN || bOMP)
2823 snew(h, hb->maxhydro);
2824 snew(g, hb->maxhydro);
2827 /* Dump hbonds for debugging */
2828 dump_ac(hb, bMerge || bContact, nDump);
2830 /* Total number of hbonds analyzed here */
2835 if (acType != AC_LUZAR && bOMP)
2837 nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
2839 gmx_omp_set_num_threads(nThreads);
2840 snew(dondata, nThreads);
2841 for (i = 0; i < nThreads; i++)
2845 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2846 "Expect close to linear scaling over this donor-loop.\n", nThreads);
2848 fprintf(stderr, "Donors: [thread no]\n");
2851 for (i = 0; i < nThreads; i++)
2853 snprintf(tmpstr, 7, "[%i]", i);
2854 fprintf(stderr, "%-7s", tmpstr);
2857 fprintf(stderr, "\n");
2861 /* Build the ACF according to acType */
2866 #ifdef HAVE_NN_LOOPS
2867 /* Here we're using the estimated energy for the hydrogen bonds. */
2870 #pragma omp parallel \
2871 private(i, j, k, nh, E, rhbex, thisThread) \
2875 thisThread = gmx_omp_get_thread_num();
2879 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2882 #pragma omp for schedule (dynamic)
2883 for (i = 0; i < hb->d.nrd; i++) /* loop over donors */
2887 #pragma omp critical
2889 dondata[thisThread] = i;
2890 parallel_print(dondata, nThreads);
2895 fprintf(stderr, "\r %i", i);
2898 for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */
2900 for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2902 E = hb->hbE.E[i][j][nh];
2905 for (k = 0; k < nframes; k++)
2907 if (E[k] != NONSENSE_E)
2909 rhbex[k] = (real)E[k];
2913 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0],
2914 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
2915 #pragma omp critical
2917 for (k = 0; (k < nn); k++)
2934 normalizeACF(ct, NULL, 0, nn);
2936 snew(timedouble, nn);
2937 for (j = 0; j < nn; j++)
2939 timedouble[j] = (double)(hb->time[j]);
2940 ctdouble[j] = (double)(ct[j]);
2943 /* Remove ballistic term */
2944 /* Ballistic component removal and fitting to the reversible geminate recombination model
2945 * will be taken out for the time being. First of all, one can remove the ballistic
2946 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2947 * problems with the robustness of the fitting to the model. More work is needed.
2948 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2949 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2950 * a BSD-licence that can do the job.
2952 * - Erik Marklund, June 18 2010.
2954 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2955 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2957 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2958 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2960 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)");
2961 xvgr_legend(fp, asize(legNN), legNN);
2963 for (j = 0; (j < nn); j++)
2965 fprintf(fp, "%10g %10g %10g\n",
2966 hb->time[j]-hb->time[0],
2974 #endif /* HAVE_NN_LOOPS */
2975 break; /* case AC_NN */
2979 memset(ct, 0, 2*n2*sizeof(real));
2981 fprintf(stderr, "Donor:\n");
2984 #define __ACDATA p_ct
2987 #pragma omp parallel \
2988 private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2989 pfound, poff, rHbExGem, p, ihb, mMax, \
2992 { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
2995 thisThread = gmx_omp_get_thread_num();
2996 snew(h, hb->maxhydro);
2997 snew(g, hb->maxhydro);
3004 memset(p_ct, 0, 2*n2*sizeof(real));
3006 /* I'm using a chunk size of 1, since I expect \
3007 * the overhead to be really small compared \
3008 * to the actual calculations \ */
3009 #pragma omp for schedule(dynamic,1) nowait
3010 for (i = 0; i < hb->d.nrd; i++)
3015 #pragma omp critical
3017 dondata[thisThread] = i;
3018 parallel_print(dondata, nThreads);
3023 fprintf(stderr, "\r %i", i);
3025 for (k = 0; k < hb->a.nra; k++)
3027 for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++)
3029 hbh = hb->hbmap[i][k];
3032 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
3033 * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
3034 * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
3035 pHist = &(hb->per->pHist[i][k]);
3036 if (ISHB(hbh->history[nh]) && pHist->len != 0)
3041 g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL;
3045 /* count the number of periodic shifts encountered and store
3046 * them in separate arrays. */
3048 for (j = 0; j < pHist->len; j++)
3051 for (m = 0; m <= np; m++)
3053 if (m == np) /* p not recognized in list. Add it and set up new array. */
3056 if (np > hb->per->nper)
3058 gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
3060 if (m >= mMax) /* Extend the arrays.
3061 * Doing it like this, using mMax to keep track of the sizes,
3062 * eleviates the need for freeing and re-allocating the arrays
3063 * when taking on the next donor-acceptor pair */
3066 srenew(pfound, np); /* The list of found periodic shifts. */
3067 srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */
3068 snew(rHbExGem[m], 2*n2);
3073 if (rHbExGem != NULL && rHbExGem[m] != NULL)
3075 /* This must be done, as this array was most likey
3076 * used to store stuff in some previous iteration. */
3077 memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
3081 fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
3093 } /* m: Loop over found shifts */
3094 } /* j: Loop over shifts */
3096 /* Now unpack and disentangle the existence funtions. */
3097 for (j = 0; j < nf; j++)
3104 * pfound: list of periodic shifts found for this pair.
3105 * poff: list of frame offsets; that is, the first
3106 * frame a hbond has a particular periodic shift. */
3107 p = getPshift(*pHist, j+n0);
3110 for (m = 0; m < np; m++)
3118 gmx_fatal(FARGS, "Shift not found, but must be there.");
3122 ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j));
3127 poff[m] = j; /* Here's where the first hbond with shift p is,
3128 * relative to the start of h[0].*/
3132 gmx_fatal(FARGS, "j<poff[m]");
3134 rHbExGem[m][j-poff[m]] += 1;
3139 /* Now, build ac. */
3140 for (m = 0; m < np; m++)
3142 if (rHbExGem[m][0] > 0 && n0+poff[m] < nn /* && m==0 */)
3144 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0],
3145 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
3146 for (j = 0; (j < nn); j++)
3148 __ACDATA[j] += rHbExGem[m][j];
3151 } /* Building of ac. */
3154 } /* hydrogen loop */
3155 } /* acceptor loop */
3158 for (m = 0; m <= mMax; m++)
3171 #pragma omp critical
3173 for (i = 0; i < nn; i++)
3181 } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
3187 normalizeACF(ct, NULL, 0, nn);
3189 fprintf(stderr, "\n\nACF successfully calculated.\n");
3191 /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
3194 snew(timedouble, nn);
3197 for (j = 0; j < nn; j++)
3199 timedouble[j] = (double)(hb->time[j]);
3200 ctdouble[j] = (double)(ct[j]);
3203 /* Remove ballistic term */
3204 /* Ballistic component removal and fitting to the reversible geminate recombination model
3205 * will be taken out for the time being. First of all, one can remove the ballistic
3206 * component with g_analyze afterwards. Secondly, and more importantly, there are still
3207 * problems with the robustness of the fitting to the model. More work is needed.
3208 * A third reason is that we're currently using gsl for this and wish to reduce dependence
3209 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
3210 * a BSD-licence that can do the job.
3212 * - Erik Marklund, June 18 2010.
3214 /* if (bBallistic) { */
3215 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
3216 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
3218 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
3219 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
3222 /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
3227 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3231 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3233 xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv);
3235 for (j = 0; (j < nn); j++)
3237 fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]);
3240 fprintf(fp, " %10g", ctdouble[j]);
3244 fprintf(fp, " %10g", fittedct[j]);
3255 break; /* case AC_GEM */
3268 for (i = 0; (i < hb->d.nrd); i++)
3270 for (k = 0; (k < hb->a.nra); k++)
3273 hbh = hb->hbmap[i][k];
3277 if (bMerge || bContact)
3279 if (ISHB(hbh->history[0]))
3288 for (m = 0; (m < hb->maxhydro); m++)
3290 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m]))
3292 g[nhydro] = hbh->g[m];
3293 h[nhydro] = hbh->h[m];
3300 for (nh = 0; (nh < nhydro); nh++)
3302 int nrint = bContact ? hb->nrdist : hb->nrhb;
3303 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
3305 fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
3308 for (j = 0; (j < nframes); j++)
3310 /* Changed '<' into '<=' below, just like I did in
3311 the hbm-output-loop in the gmx_hbond() block.
3312 - Erik Marklund, May 31, 2006 */
3315 ihb = is_hb(h[nh], j);
3316 idist = is_hb(g[nh], j);
3323 /* For contacts: if a second cut-off is provided, use it,
3324 * otherwise use g(t) = 1-h(t) */
3325 if (!R2 && bContact)
3331 gt[j] = idist*(1-ihb);
3338 /* The autocorrelation function is normalized after summation only */
3339 low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
3340 eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1);
3342 /* Cross correlation analysis for thermodynamics */
3343 for (j = nframes; (j < n2); j++)
3349 cross_corr(n2, ht, gt, dght);
3351 for (j = 0; (j < nn); j++)
3360 fprintf(stderr, "\n");
3363 normalizeACF(ct, ght, nhb, nn);
3365 /* Determine tail value for statistics */
3368 for (j = nn/2; (j < nn); j++)
3371 tail2 += ct[j]*ct[j];
3373 tail /= (nn - nn/2);
3374 tail2 /= (nn - nn/2);
3375 dtail = sqrt(tail2-tail*tail);
3377 /* Check whether the ACF is long enough */
3380 printf("\nWARNING: Correlation function is probably not long enough\n"
3381 "because the standard deviation in the tail of C(t) > %g\n"
3382 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
3385 for (j = 0; (j < nn); j++)
3388 ct[j] = (cct[j]-tail)/(1-tail);
3390 /* Compute negative derivative k(t) = -dc(t)/dt */
3391 compute_derivative(nn, hb->time, ct, kt);
3392 for (j = 0; (j < nn); j++)
3400 fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3404 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv);
3406 xvgr_legend(fp, asize(legLuzar), legLuzar, oenv);
3409 for (j = 0; (j < nn); j++)
3411 fprintf(fp, "%10g %10g %10g %10g %10g\n",
3412 hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
3416 analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
3417 fit_start, temp, smooth_tail_start, oenv);
3419 do_view(oenv, fn, NULL);
3431 break; /* case AC_LUZAR */
3434 gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
3435 } /* switch (acType) */
3438 static void init_hbframe(t_hbdata *hb, int nframes, real t)
3442 hb->time[nframes] = t;
3443 hb->nhb[nframes] = 0;
3444 hb->ndist[nframes] = 0;
3445 for (i = 0; (i < max_hx); i++)
3447 hb->nhx[nframes][i] = 0;
3449 /* Loop invalidated */
3450 if (hb->bHBmap && 0)
3452 for (i = 0; (i < hb->d.nrd); i++)
3454 for (j = 0; (j < hb->a.nra); j++)
3456 for (m = 0; (m < hb->maxhydro); m++)
3458 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
3460 set_hb(hb, i, m, j, nframes, HB_NO);
3466 /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
3469 static void analyse_donor_props(const char *fn, t_hbdata *hb, int nframes, real t,
3470 const output_env_t oenv)
3472 static FILE *fp = NULL;
3473 const char *leg[] = { "Nbound", "Nfree" };
3474 int i, j, k, nbound, nb, nhtot;
3482 fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
3483 xvgr_legend(fp, asize(leg), leg, oenv);
3487 for (i = 0; (i < hb->d.nrd); i++)
3489 for (k = 0; (k < hb->d.nhydro[i]); k++)
3493 for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
3495 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
3496 is_hb(hb->hbmap[i][j]->h[k], nframes))
3504 fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound);
3507 static void dump_hbmap(t_hbdata *hb,
3508 int nfile, t_filenm fnm[], gmx_bool bTwo,
3509 gmx_bool bContact, int isize[], int *index[], char *grpnames[],
3513 int ddd, hhh, aaa, i, j, k, m, grp;
3514 char ds[32], hs[32], as[32];
3517 fp = opt2FILE("-hbn", nfile, fnm, "w");
3518 if (opt2bSet("-g", nfile, fnm))
3520 fplog = ffopen(opt2fn("-g", nfile, fnm), "w");
3521 fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
3527 for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++)
3529 fprintf(fp, "[ %s ]", grpnames[grp]);
3530 for (i = 0; i < isize[grp]; i++)
3532 fprintf(fp, (i%15) ? " " : "\n");
3533 fprintf(fp, " %4u", index[grp][i]+1);
3537 Added -contact support below.
3538 - Erik Marklund, May 29, 2006
3542 fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]);
3543 for (i = 0; (i < hb->d.nrd); i++)
3545 if (hb->d.grp[i] == grp)
3547 for (j = 0; (j < hb->d.nhydro[i]); j++)
3549 fprintf(fp, " %4u %4u", hb->d.don[i]+1,
3550 hb->d.hydro[i][j]+1);
3556 fprintf(fp, "[ acceptors_%s ]", grpnames[grp]);
3557 for (i = 0; (i < hb->a.nra); i++)
3559 if (hb->a.grp[i] == grp)
3561 fprintf(fp, (i%15 && !first) ? " " : "\n");
3562 fprintf(fp, " %4u", hb->a.acc[i]+1);
3571 fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
3572 "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
3576 fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]);
3579 for (i = 0; (i < hb->d.nrd); i++)
3582 for (k = 0; (k < hb->a.nra); k++)
3585 for (m = 0; (m < hb->d.nhydro[i]); m++)
3587 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m]))
3589 sprintf(ds, "%s", mkatomname(atoms, ddd));
3590 sprintf(as, "%s", mkatomname(atoms, aaa));
3593 fprintf(fp, " %6u %6u\n", ddd+1, aaa+1);
3596 fprintf(fplog, "%12s %12s\n", ds, as);
3601 hhh = hb->d.hydro[i][m];
3602 sprintf(hs, "%s", mkatomname(atoms, hhh));
3603 fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1);
3606 fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
3620 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3621 * It mimics add_frames() and init_frame() to some extent. */
3622 static void sync_hbdata(t_hbdata *hb, t_hbdata *p_hb,
3623 int nframes, real t)
3626 if (nframes >= p_hb->max_frames)
3628 p_hb->max_frames += 4096;
3629 srenew(p_hb->nhb, p_hb->max_frames);
3630 srenew(p_hb->ndist, p_hb->max_frames);
3631 srenew(p_hb->n_bound, p_hb->max_frames);
3632 srenew(p_hb->nhx, p_hb->max_frames);
3635 srenew(p_hb->danr, p_hb->max_frames);
3637 memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3638 memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3639 p_hb->nhb[nframes] = 0;
3640 p_hb->ndist[nframes] = 0;
3643 p_hb->nframes = nframes;
3646 /* p_hb->nhx[nframes][i] */
3648 memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
3650 /* hb->per will remain constant througout the frame loop,
3651 * even though the data its members point to will change,
3652 * hence no need for re-syncing. */
3655 int gmx_hbond(int argc, char *argv[])
3657 const char *desc[] = {
3658 "[TT]g_hbond[tt] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3659 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
3660 "(zero is extended) and the distance Donor - Acceptor",
3661 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
3662 "OH and NH groups are regarded as donors, O is an acceptor always,",
3663 "N is an acceptor by default, but this can be switched using",
3664 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3665 "to the first preceding non-hydrogen atom.[PAR]",
3667 "You need to specify two groups for analysis, which must be either",
3668 "identical or non-overlapping. All hydrogen bonds between the two",
3669 "groups are analyzed.[PAR]",
3671 "If you set [TT]-shell[tt], you will be asked for an additional index group",
3672 "which should contain exactly one atom. In this case, only hydrogen",
3673 "bonds between atoms within the shell distance from the one atom are",
3676 "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3677 "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3678 "If contact kinetics are analyzed by using the -contact option, then",
3679 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3680 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3681 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3682 "See mentioned literature for more details and definitions."
3685 /* "It is also possible to analyse specific hydrogen bonds with",
3686 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3687 "Donor Hydrogen Acceptor, in the following way:[PAR]",
3695 "Note that the triplets need not be on separate lines.",
3696 "Each atom triplet specifies a hydrogen bond to be analyzed,",
3697 "note also that no check is made for the types of atoms.[PAR]",
3699 "[BB]Output:[bb][BR]",
3700 "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
3701 "[TT]-ac[tt]: average over all autocorrelations of the existence",
3702 "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3703 "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3704 "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
3705 "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
3706 "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3707 "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3708 "with helices in proteins.[BR]",
3709 "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
3710 "for selected groups, all hydrogen bonded atoms from all groups and",
3711 "all solvent atoms involved in insertion.[BR]",
3712 "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
3713 "frames, this also contains information on solvent insertion",
3714 "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3716 "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3717 "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3718 "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3719 "compare results to Raman Spectroscopy.",
3721 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3722 "require an amount of memory proportional to the total numbers of donors",
3723 "times the total number of acceptors in the selected group(s)."
3726 static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
3727 static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1;
3728 static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
3729 static int nDump = 0, nFitPoints = 100;
3730 static int nThreads = 0, nBalExp = 4;
3732 static gmx_bool bContact = FALSE, bBallistic = FALSE, bBallisticDt = FALSE, bGemFit = FALSE;
3733 static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
3734 static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3738 { "-a", FALSE, etREAL, {&acut},
3739 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
3740 { "-r", FALSE, etREAL, {&rcut},
3741 "Cutoff radius (nm, X - Acceptor, see next option)" },
3742 { "-da", FALSE, etBOOL, {&bDA},
3743 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3744 { "-r2", FALSE, etREAL, {&r2cut},
3745 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3746 { "-abin", FALSE, etREAL, {&abin},
3747 "Binwidth angle distribution (degrees)" },
3748 { "-rbin", FALSE, etREAL, {&rbin},
3749 "Binwidth distance distribution (nm)" },
3750 { "-nitacc", FALSE, etBOOL, {&bNitAcc},
3751 "Regard nitrogen atoms as acceptors" },
3752 { "-contact", FALSE, etBOOL, {&bContact},
3753 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3754 { "-shell", FALSE, etREAL, {&rshell},
3755 "when > 0, only calculate hydrogen bonds within # nm shell around "
3757 { "-fitstart", FALSE, etREAL, {&fit_start},
3758 "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]" },
3759 { "-fitstart", FALSE, etREAL, {&fit_start},
3760 "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])" },
3761 { "-temp", FALSE, etREAL, {&temp},
3762 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3763 { "-smooth", FALSE, etREAL, {&smooth_tail_start},
3764 "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3765 { "-dump", FALSE, etINT, {&nDump},
3766 "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
3767 { "-max_hb", FALSE, etREAL, {&maxnhb},
3768 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3769 { "-merge", FALSE, etBOOL, {&bMerge},
3770 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3771 { "-geminate", FALSE, etENUM, {gemType},
3772 "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3773 { "-diff", FALSE, etREAL, {&D},
3774 "Dffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
3776 { "-nthreads", FALSE, etINT, {&nThreads},
3777 "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
3780 const char *bugs[] = {
3781 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3784 { efTRX, "-f", NULL, ffREAD },
3785 { efTPX, NULL, NULL, ffREAD },
3786 { efNDX, NULL, NULL, ffOPTRD },
3787 /* { efNDX, "-sel", "select", ffOPTRD },*/
3788 { efXVG, "-num", "hbnum", ffWRITE },
3789 { efLOG, "-g", "hbond", ffOPTWR },
3790 { efXVG, "-ac", "hbac", ffOPTWR },
3791 { efXVG, "-dist", "hbdist", ffOPTWR },
3792 { efXVG, "-ang", "hbang", ffOPTWR },
3793 { efXVG, "-hx", "hbhelix", ffOPTWR },
3794 { efNDX, "-hbn", "hbond", ffOPTWR },
3795 { efXPM, "-hbm", "hbmap", ffOPTWR },
3796 { efXVG, "-don", "donor", ffOPTWR },
3797 { efXVG, "-dan", "danum", ffOPTWR },
3798 { efXVG, "-life", "hblife", ffOPTWR },
3799 { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
3802 #define NFILE asize(fnm)
3804 char hbmap [HB_NR] = { ' ', 'o', '-', '*' };
3805 const char *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
3806 t_rgb hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
3808 t_trxstatus *status;
3813 int npargs, natoms, nframes = 0, shatom;
3819 real t, ccut, dist = 0.0, ang = 0.0;
3820 double max_nhb, aver_nhb, aver_dist;
3821 int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
3823 int xj, yj, zj, aj, xjj, yjj, zjj;
3824 int xk, yk, zk, ak, xkk, ykk, zkk;
3825 gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
3826 int *adist, *rdist, *aptr, *rprt;
3827 int grp, nabin, nrbin, bin, resdist, ihb;
3829 t_hbdata *hb, *hbptr;
3830 FILE *fp, *fpins = NULL, *fpnhb = NULL;
3832 t_ncell *icell, *jcell, *kcell;
3834 unsigned char *datable;
3839 int ii, jj, hh, actual_nThreads;
3841 gmx_bool bGem, bNN, bParallel;
3842 t_gemParams *params = NULL;
3843 gmx_bool bEdge_yjj, bEdge_xjj, bOMP;
3845 t_hbdata **p_hb = NULL; /* one per thread, then merge after the frame loop */
3846 int **p_adist = NULL, **p_rdist = NULL; /* a histogram for each thread. */
3855 ppa = add_acf_pargs(&npargs, pa);
3857 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE, NFILE, fnm, npargs,
3858 ppa, asize(desc), desc, asize(bugs), bugs, &oenv);
3860 /* NN-loop? If so, what estimator to use ?*/
3862 /* Outcommented for now DvdS 2010-07-13
3863 while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3866 gmx_fatal(FARGS, "Invalid NN-loop type.");
3869 for (i = 2; bNN == FALSE && i < NN_NR; i++)
3871 bNN = bNN || NN == i;
3874 if (NN > NN_NONE && bMerge)
3879 /* geminate recombination? If so, which flavor? */
3881 while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0)
3885 if (gemmode == gemNR)
3887 gmx_fatal(FARGS, "Invalid recombination type.");
3891 for (i = 2; bGem == FALSE && i < gemNR; i++)
3893 bGem = bGem || gemmode == i;
3898 printf("Geminate recombination: %s\n", gemType[gemmode]);
3900 printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
3904 if (gemmode != gemDD)
3906 printf("Turning off -contact option...\n");
3912 if (gemmode == gemDD)
3914 printf("Turning on -contact option...\n");
3920 if (gemmode == gemAA)
3922 printf("Turning off -merge option...\n");
3928 if (gemmode != gemAA)
3930 printf("Turning on -merge option...\n");
3938 ccut = cos(acut*DEG2RAD);
3944 gmx_fatal(FARGS, "Can not analyze selected contacts.");
3948 gmx_fatal(FARGS, "Can not analyze contact between H and A: turn off -noda");
3952 /* Initiate main data structure! */
3953 bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
3954 opt2bSet("-life", NFILE, fnm) ||
3955 opt2bSet("-hbn", NFILE, fnm) ||
3956 opt2bSet("-hbm", NFILE, fnm) ||
3959 if (opt2bSet("-nhbdist", NFILE, fnm))
3961 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3962 fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
3963 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
3964 xvgr_legend(fpnhb, asize(leg), leg, oenv);
3967 hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
3970 read_tpx_top(ftp2fn(efTPX, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
3972 snew(grpnames, grNR);
3975 /* Make Donor-Acceptor table */
3976 snew(datable, top.atoms.nr);
3977 gen_datable(index[0], isize[0], datable, top.atoms.nr);
3981 /* analyze selected hydrogen bonds */
3982 printf("Select group with selected atoms:\n");
3983 get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
3984 1, &nsel, index, grpnames);
3987 gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
3988 "and therefore cannot contain triplets of "
3989 "Donor-Hydrogen-Acceptor", grpnames[0]);
3993 for (i = 0; (i < nsel); i += 3)
3995 int dd = index[0][i];
3996 int aa = index[0][i+2];
3997 /* int */ hh = index[0][i+1];
3998 add_dh (&hb->d, dd, hh, i, datable);
3999 add_acc(&hb->a, aa, i);
4000 /* Should this be here ? */
4001 snew(hb->d.dptr, top.atoms.nr);
4002 snew(hb->a.aptr, top.atoms.nr);
4003 add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri);
4005 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
4006 isize[0], grpnames[0]);
4010 /* analyze all hydrogen bonds: get group(s) */
4011 printf("Specify 2 groups to analyze:\n");
4012 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4013 2, isize, index, grpnames);
4015 /* check if we have two identical or two non-overlapping groups */
4016 bTwo = isize[0] != isize[1];
4017 for (i = 0; (i < isize[0]) && !bTwo; i++)
4019 bTwo = index[0][i] != index[1][i];
4023 printf("Checking for overlap in atoms between %s and %s\n",
4024 grpnames[0], grpnames[1]);
4025 for (i = 0; i < isize[1]; i++)
4027 if (ISINGRP(datable[index[1][i]]))
4029 gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
4030 grpnames[0], grpnames[1]);
4034 printf("Checking for overlap in atoms between %s and %s\n",
4035 grpnames[0],grpnames[1]);
4036 for (i=0; i<isize[0]; i++)
4037 for (j=0; j<isize[1]; j++)
4038 if (index[0][i] == index[1][j])
4039 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
4040 grpnames[0],grpnames[1]);
4045 printf("Calculating %s "
4046 "between %s (%d atoms) and %s (%d atoms)\n",
4047 bContact ? "contacts" : "hydrogen bonds",
4048 grpnames[0], isize[0], grpnames[1], isize[1]);
4052 fprintf(stderr, "Calculating %s in %s (%d atoms)\n",
4053 bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]);
4058 /* search donors and acceptors in groups */
4059 snew(datable, top.atoms.nr);
4060 for (i = 0; (i < grNR); i++)
4062 if ( ((i == gr0) && !bSelected ) ||
4063 ((i == gr1) && bTwo ))
4065 gen_datable(index[i], isize[i], datable, top.atoms.nr);
4068 search_acceptors(&top, isize[i], index[i], &hb->a, i,
4069 bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
4070 search_donors (&top, isize[i], index[i], &hb->d, i,
4071 TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
4075 search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
4076 search_donors (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
4080 clear_datable_grp(datable, top.atoms.nr);
4085 printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra);
4087 snew(donors[gr0D], dons[gr0D].nrd);*/
4091 printf("Making hbmap structure...");
4092 /* Generate hbond data structure */
4097 #ifdef HAVE_NN_LOOPS
4106 printf("Making per structure...");
4107 /* Generate hbond data structure */
4114 if (hb->d.nrd + hb->a.nra == 0)
4116 printf("No Donors or Acceptors found\n");
4123 printf("No Donors found\n");
4128 printf("No Acceptors found\n");
4134 gmx_fatal(FARGS, "Nothing to be done");
4143 /* get index group with atom for shell */
4146 printf("Select atom for shell (1 atom):\n");
4147 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
4148 1, &shisz, &shidx, &shgrpnm);
4151 printf("group contains %d atoms, should be 1 (one)\n", shisz);
4156 printf("Will calculate hydrogen bonds within a shell "
4157 "of %g nm around atom %i\n", rshell, shatom+1);
4160 /* Analyze trajectory */
4161 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
4162 if (natoms > top.atoms.nr)
4164 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
4165 top.atoms.nr, natoms);
4168 bBox = ir.ePBC != epbcNONE;
4169 grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
4172 snew(adist, nabin+1);
4173 snew(rdist, nrbin+1);
4177 gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
4183 #define __ADIST adist
4184 #define __RDIST rdist
4186 #else /* GMX_OPENMP ================================================== \
4187 * Set up the OpenMP stuff, |
4188 * like the number of threads and such |
4189 * Also start the parallel loop. |
4191 #define __ADIST p_adist[threadNr]
4192 #define __RDIST p_rdist[threadNr]
4193 #define __HBDATA p_hb[threadNr]
4197 bParallel = !bSelected;
4201 actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
4203 gmx_omp_set_num_threads(actual_nThreads);
4204 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
4209 actual_nThreads = 1;
4212 snew(p_hb, actual_nThreads);
4213 snew(p_adist, actual_nThreads);
4214 snew(p_rdist, actual_nThreads);
4215 for (i = 0; i < actual_nThreads; i++)
4218 snew(p_adist[i], nabin+1);
4219 snew(p_rdist[i], nrbin+1);
4221 p_hb[i]->max_frames = 0;
4222 p_hb[i]->nhb = NULL;
4223 p_hb[i]->ndist = NULL;
4224 p_hb[i]->n_bound = NULL;
4225 p_hb[i]->time = NULL;
4226 p_hb[i]->nhx = NULL;
4228 p_hb[i]->bHBmap = hb->bHBmap;
4229 p_hb[i]->bDAnr = hb->bDAnr;
4230 p_hb[i]->bGem = hb->bGem;
4231 p_hb[i]->wordlen = hb->wordlen;
4232 p_hb[i]->nframes = hb->nframes;
4233 p_hb[i]->maxhydro = hb->maxhydro;
4234 p_hb[i]->danr = hb->danr;
4237 p_hb[i]->hbmap = hb->hbmap;
4238 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
4239 p_hb[i]->per = hb->per;
4241 #ifdef HAVE_NN_LOOPS
4242 p_hb[i]->hbE = hb->hbE;
4246 p_hb[i]->nrdist = 0;
4250 /* Make a thread pool here,
4251 * instead of forking anew at every frame. */
4253 #pragma omp parallel \
4255 private(j, h, ii, jj, hh, E, \
4256 xi, yi, zi, xj, yj, zj, threadNr, \
4257 dist, ang, peri, icell, jcell, \
4258 grp, ogrp, ai, aj, xjj, yjj, zjj, \
4259 xk, yk, zk, ihb, id, resdist, \
4260 xkk, ykk, zkk, kcell, ak, k, bTric, \
4261 bEdge_xjj, bEdge_yjj) \
4263 { /* Start of parallel region */
4264 threadNr = gmx_omp_get_thread_num();
4269 bTric = bBox && TRICLINIC(box);
4273 sync_hbdata(hb, p_hb[threadNr], nframes, t);
4277 build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut,
4278 rshell, ngrid, grid);
4279 reset_nhbonds(&(hb->d));
4281 if (debug && bDebug)
4283 dump_grid(debug, ngrid, grid);
4286 add_frames(hb, nframes);
4287 init_hbframe(hb, nframes, output_env_conv_time(oenv, t));
4291 count_da_grid(ngrid, grid, hb->danr[nframes]);
4297 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
4302 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
4303 /* Loop over all atom pairs and estimate interaction energy */
4307 addFramesNN(hb, nframes);
4311 #pragma omp for schedule(dynamic)
4312 for (i = 0; i < hb->d.nrd; i++)
4314 for (j = 0; j < hb->a.nra; j++)
4317 h < (bContact ? 1 : hb->d.nhydro[i]);
4320 if (i == hb->d.nrd || j == hb->a.nra)
4322 gmx_fatal(FARGS, "out of bounds");
4325 /* Get the real atom ids */
4328 hh = hb->d.hydro[i][h];
4330 /* Estimate the energy from the geometry */
4331 E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
4332 /* Store the energy */
4333 storeHbEnergy(hb, i, j, h, E, nframes);
4337 #endif /* HAVE_NN_LOOPS */
4346 /* Do not parallelize this just yet. */
4348 for (ii = 0; (ii < nsel); ii++)
4350 int dd = index[0][i];
4351 int aa = index[0][i+2];
4352 /* int */ hh = index[0][i+1];
4353 ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
4354 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4358 /* add to index if not already there */
4360 add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri);
4364 } /* if (bSelected) */
4372 calcBoxProjection(box, hb->per->P);
4375 /* loop over all gridcells (xi,yi,zi) */
4376 /* Removed confusing macro, DvdS 27/12/98 */
4379 /* The outer grid loop will have to do for now. */
4380 #pragma omp for schedule(dynamic)
4381 for (xi = 0; xi < ngrid[XX]; xi++)
4383 for (yi = 0; (yi < ngrid[YY]); yi++)
4385 for (zi = 0; (zi < ngrid[ZZ]); zi++)
4388 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
4389 for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++)
4391 icell = &(grid[zi][yi][xi].d[grp]);
4402 /* loop over all hydrogen atoms from group (grp)
4403 * in this gridcell (icell)
4405 for (ai = 0; (ai < icell->nr); ai++)
4407 i = icell->atoms[ai];
4409 /* loop over all adjacent gridcells (xj,yj,zj) */
4410 for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
4411 zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
4414 zj = grid_mod(zjj, ngrid[ZZ]);
4415 bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
4416 for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj);
4417 yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj);
4420 yj = grid_mod(yjj, ngrid[YY]);
4422 (yj == 0) || (yj == ngrid[YY] - 1) ||
4423 (zj == 0) || (zj == ngrid[ZZ] - 1);
4424 for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
4425 xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
4428 xj = grid_mod(xjj, ngrid[XX]);
4429 jcell = &(grid[zj][yj][xj].a[ogrp]);
4430 /* loop over acceptor atoms from other group (ogrp)
4431 * in this adjacent gridcell (jcell)
4433 for (aj = 0; (aj < jcell->nr); aj++)
4435 j = jcell->atoms[aj];
4437 /* check if this once was a h-bond */
4439 ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
4440 hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri);
4444 /* add to index if not already there */
4446 add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri);
4448 /* make angle and distance distributions */
4449 if (ihb == hbHB && !bContact)
4453 gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
4456 __ADIST[(int)( ang/abin)]++;
4457 __RDIST[(int)(dist/rbin)]++;
4461 if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
4463 gmx_fatal(FARGS, "Invalid donor %d", i);
4465 if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
4467 gmx_fatal(FARGS, "Invalid acceptor %d", j);
4469 resdist = abs(top.atoms.atom[i].resind-
4470 top.atoms.atom[j].resind);
4471 if (resdist >= max_hx)
4475 __HBDATA->nhx[nframes][resdist]++;
4486 } /* for xi,yi,zi */
4489 } /* if (bSelected) {...} else */
4492 /* Better wait for all threads to finnish using x[] before updating it. */
4495 #pragma omp critical
4497 /* Sum up histograms and counts from p_hb[] into hb */
4500 hb->nhb[k] += p_hb[threadNr]->nhb[k];
4501 hb->ndist[k] += p_hb[threadNr]->ndist[k];
4502 for (j = 0; j < max_hx; j++)
4504 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
4509 /* Here are a handful of single constructs
4510 * to share the workload a bit. The most
4511 * important one is of course the last one,
4512 * where there's a potential bottleneck in form
4519 analyse_donor_props(opt2fn_null("-don", NFILE, fnm), hb, k, t, oenv);
4527 do_nhb_dist(fpnhb, hb, t);
4530 } /* if (bNN) {...} else + */
4534 trrStatus = (read_next_x(oenv, status, &t, natoms, x, box));
4544 #pragma omp critical
4546 hb->nrhb += p_hb[threadNr]->nrhb;
4547 hb->nrdist += p_hb[threadNr]->nrdist;
4549 /* Free parallel datastructures */
4550 sfree(p_hb[threadNr]->nhb);
4551 sfree(p_hb[threadNr]->ndist);
4552 sfree(p_hb[threadNr]->nhx);
4555 for (i = 0; i < nabin; i++)
4557 for (j = 0; j < actual_nThreads; j++)
4560 adist[i] += p_adist[j][i];
4564 for (i = 0; i <= nrbin; i++)
4566 for (j = 0; j < actual_nThreads; j++)
4568 rdist[i] += p_rdist[j][i];
4572 sfree(p_adist[threadNr]);
4573 sfree(p_rdist[threadNr]);
4575 } /* End of parallel region */
4582 if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
4584 gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
4587 free_grid(ngrid, &grid);
4595 /* Compute maximum possible number of different hbonds */
4602 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
4604 /* Added support for -contact below.
4605 * - Erik Marklund, May 29-31, 2006 */
4606 /* Changed contact code.
4607 * - Erik Marklund, June 29, 2006 */
4612 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
4616 printf("Found %d different %s in trajectory\n"
4617 "Found %d different atom-pairs within %s distance\n",
4618 hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
4619 hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
4621 /*Control the pHist.*/
4625 merge_hb(hb, bTwo, bContact);
4628 if (opt2bSet("-hbn", NFILE, fnm))
4630 dump_hbmap(hb, NFILE, fnm, bTwo, bContact, isize, index, grpnames, &top.atoms);
4633 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
4634 * to make the -hbn and -hmb output match eachother.
4635 * - Erik Marklund, May 30, 2006 */
4638 /* Print out number of hbonds and distances */
4641 fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
4642 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
4644 snew(leg[0], STRLEN);
4645 snew(leg[1], STRLEN);
4646 sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds");
4647 sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut);
4648 xvgr_legend(fp, 2, (const char**)leg, oenv);
4652 for (i = 0; (i < nframes); i++)
4654 fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
4655 aver_nhb += hb->nhb[i];
4656 aver_dist += hb->ndist[i];
4659 aver_nhb /= nframes;
4660 aver_dist /= nframes;
4661 /* Print HB distance distribution */
4662 if (opt2bSet("-dist", NFILE, fnm))
4667 for (i = 0; i < nrbin; i++)
4672 fp = xvgropen(opt2fn("-dist", NFILE, fnm),
4673 "Hydrogen Bond Distribution",
4675 "Donor - Acceptor Distance (nm)" :
4676 "Hydrogen - Acceptor Distance (nm)", "", oenv);
4677 for (i = 0; i < nrbin; i++)
4679 fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
4684 /* Print HB angle distribution */
4685 if (opt2bSet("-ang", NFILE, fnm))
4690 for (i = 0; i < nabin; i++)
4695 fp = xvgropen(opt2fn("-ang", NFILE, fnm),
4696 "Hydrogen Bond Distribution",
4697 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
4698 for (i = 0; i < nabin; i++)
4700 fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
4705 /* Print HB in alpha-helix */
4706 if (opt2bSet("-hx", NFILE, fnm))
4708 fp = xvgropen(opt2fn("-hx", NFILE, fnm),
4709 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4710 xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
4711 for (i = 0; i < nframes; i++)
4713 fprintf(fp, "%10g", hb->time[i]);
4714 for (j = 0; j < max_hx; j++)
4716 fprintf(fp, " %6d", hb->nhx[i][j]);
4724 printf("Average number of %s per timeframe %.3f out of %g possible\n",
4725 bContact ? "contacts" : "hbonds",
4726 bContact ? aver_dist : aver_nhb, max_nhb);
4729 /* Do Autocorrelation etc. */
4733 Added support for -contact in ac and hbm calculations below.
4734 - Erik Marklund, May 29, 2006
4738 if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
4740 please_cite(stdout, "Spoel2006b");
4742 if (opt2bSet("-ac", NFILE, fnm))
4744 char *gemstring = NULL;
4748 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4749 gemBallistic, nBalExp, bBallisticDt);
4752 gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4755 gemstring = strdup(gemType[hb->per->gemtype]);
4756 do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
4757 bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
4758 params, gemstring, nThreads, NN, bBallistic, bGemFit);
4760 if (opt2bSet("-life", NFILE, fnm))
4762 do_hblife(opt2fn("-life", NFILE, fnm), hb, bMerge, bContact, oenv);
4764 if (opt2bSet("-hbm", NFILE, fnm))
4767 int id, ia, hh, x, y;
4769 if ((nframes > 0) && (hb->nrhb > 0))
4774 snew(mat.matrix, mat.nx);
4775 for (x = 0; (x < mat.nx); x++)
4777 snew(mat.matrix[x], mat.ny);
4780 for (id = 0; (id < hb->d.nrd); id++)
4782 for (ia = 0; (ia < hb->a.nra); ia++)
4784 for (hh = 0; (hh < hb->maxhydro); hh++)
4786 if (hb->hbmap[id][ia])
4788 if (ISHB(hb->hbmap[id][ia]->history[hh]))
4790 /* Changed '<' into '<=' in the for-statement below.
4791 * It fixed the previously undiscovered bug that caused
4792 * the last occurance of an hbond/contact to not be
4793 * set in mat.matrix. Have a look at any old -hbm-output
4794 * and you will notice that the last column is allways empty.
4795 * - Erik Marklund May 30, 2006
4797 for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++)
4799 int nn0 = hb->hbmap[id][ia]->n0;
4800 range_check(y, 0, mat.ny);
4801 mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x);
4809 mat.axis_x = hb->time;
4810 snew(mat.axis_y, mat.ny);
4811 for (j = 0; j < mat.ny; j++)
4815 sprintf(mat.title, bContact ? "Contact Existence Map" :
4816 "Hydrogen Bond Existence Map");
4817 sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds");
4818 sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv));
4819 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4820 mat.bDiscrete = TRUE;
4822 snew(mat.map, mat.nmap);
4823 for (i = 0; i < mat.nmap; i++)
4825 mat.map[i].code.c1 = hbmap[i];
4826 mat.map[i].desc = hbdesc[i];
4827 mat.map[i].rgb = hbrgb[i];
4829 fp = opt2FILE("-hbm", NFILE, fnm, "w");
4830 write_xpm_m(fp, mat);
4832 for (x = 0; x < mat.nx; x++)
4834 sfree(mat.matrix[x]);
4842 fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
4849 fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4850 fprintf(stderr, "Freeing pHist for all donors...\n");
4851 for (i = 0; i < hb->d.nrd; i++)
4853 fprintf(stderr, "\r%i", i);
4854 if (hb->per->pHist[i] != NULL)
4856 for (j = 0; j < hb->a.nra; j++)
4858 clearPshift(&(hb->per->pHist[i][j]));
4860 sfree(hb->per->pHist[i]);
4863 sfree(hb->per->pHist);
4864 sfree(hb->per->p2i);
4866 fprintf(stderr, "...done.\n");
4869 #ifdef HAVE_NN_LOOPS
4882 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4884 fp = xvgropen(opt2fn("-dan", NFILE, fnm),
4885 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
4886 nleg = (bTwo ? 2 : 1)*2;
4887 snew(legnames, nleg);
4889 for (j = 0; j < grNR; j++)
4891 if (USE_THIS_GROUP(j) )
4893 sprintf(buf, "Donors %s", grpnames[j]);
4894 legnames[i++] = strdup(buf);
4895 sprintf(buf, "Acceptors %s", grpnames[j]);
4896 legnames[i++] = strdup(buf);
4901 gmx_incons("number of legend entries");
4903 xvgr_legend(fp, nleg, (const char**)legnames, oenv);
4904 for (i = 0; i < nframes; i++)
4906 fprintf(fp, "%10g", hb->time[i]);
4907 for (j = 0; (j < grNR); j++)
4909 if (USE_THIS_GROUP(j) )
4911 fprintf(fp, " %6d", hb->danr[i][j]);