2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/commandline/viewit.h"
49 #include "gromacs/correlationfunctions/autocorr.h"
50 #include "gromacs/correlationfunctions/crosscorr.h"
51 #include "gromacs/correlationfunctions/expfit.h"
52 #include "gromacs/correlationfunctions/integrate.h"
53 #include "gromacs/fileio/copyrite.h"
54 #include "gromacs/fileio/matio.h"
55 #include "gromacs/fileio/tpxio.h"
56 #include "gromacs/fileio/trxio.h"
57 #include "gromacs/fileio/xvgr.h"
58 #include "gromacs/gmxana/gmx_ana.h"
59 #include "gromacs/gmxana/gstat.h"
60 #include "gromacs/gmxlib/ifunc.h"
61 #include "gromacs/math/units.h"
62 #include "gromacs/math/vec.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/topology/index.h"
65 #include "gromacs/topology/topology.h"
66 #include "gromacs/utility/arraysize.h"
67 #include "gromacs/utility/cstringutil.h"
68 #include "gromacs/utility/exceptions.h"
69 #include "gromacs/utility/fatalerror.h"
70 #include "gromacs/utility/futil.h"
71 #include "gromacs/utility/gmxomp.h"
72 #include "gromacs/utility/programcontext.h"
73 #include "gromacs/utility/smalloc.h"
74 #include "gromacs/utility/snprintf.h"
77 typedef int t_hx
[max_hx
];
78 #define NRHXTYPES max_hx
79 const char *hxtypenames
[NRHXTYPES
] =
80 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
83 static const int NOTSET
= -49297;
86 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == 0)
88 #define MASTER_THREAD_ONLY(threadNr) ((threadNr) == (threadNr))
91 /* -----------------------------------------*/
97 hbNo
, hbDist
, hbHB
, hbNR
, hbR2
100 noDA
, ACC
, DON
, DA
, INGROUP
103 static const char *grpnames
[grNR
] = {"0", "1", "I" };
105 static gmx_bool bDebug
= FALSE
;
110 #define HB_YESINS HB_YES|HB_INS
114 #define ISHB(h) (((h) & 2) == 2)
115 #define ISDIST(h) (((h) & 1) == 1)
116 #define ISDIST2(h) (((h) & 4) == 4)
117 #define ISACC(h) (((h) & 1) == 1)
118 #define ISDON(h) (((h) & 2) == 2)
119 #define ISINGRP(h) (((h) & 4) == 4)
132 typedef int t_icell
[grNR
];
133 typedef int h_id
[MAXHYDRO
];
136 int history
[MAXHYDRO
];
137 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
138 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
140 /* Bitmask array which tells whether a hbond is present
141 * at a given time. Either of these may be NULL
143 int n0
; /* First frame a HB was found */
144 int nframes
, maxframes
; /* Amount of frames in this hbond */
147 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
148 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
149 * acceptor distance is less than the user-specified distance (typically
156 int *acc
; /* Atom numbers of the acceptors */
157 int *grp
; /* Group index */
158 int *aptr
; /* Map atom number to acceptor index */
163 int *don
; /* Atom numbers of the donors */
164 int *grp
; /* Group index */
165 int *dptr
; /* Map atom number to donor index */
166 int *nhydro
; /* Number of hydrogens for each donor */
167 h_id
*hydro
; /* The atom numbers of the hydrogens */
168 h_id
*nhbonds
; /* The number of HBs per H at current */
172 gmx_bool bHBmap
, bDAnr
;
174 /* The following arrays are nframes long */
175 int nframes
, max_frames
, maxhydro
;
181 /* These structures are initialized from the topology at start up */
184 /* This holds a matrix with all possible hydrogen bonds */
189 /* Changed argument 'bMerge' into 'oneHB' below,
190 * since -contact should cause maxhydro to be 1,
192 * - Erik Marklund May 29, 2006
195 static t_hbdata
*mk_hbdata(gmx_bool bHBmap
, gmx_bool bDAnr
, gmx_bool oneHB
)
200 hb
->wordlen
= 8*sizeof(unsigned int);
209 hb
->maxhydro
= MAXHYDRO
;
214 static void mk_hbmap(t_hbdata
*hb
)
218 snew(hb
->hbmap
, hb
->d
.nrd
);
219 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
221 snew(hb
->hbmap
[i
], hb
->a
.nra
);
222 if (hb
->hbmap
[i
] == NULL
)
224 gmx_fatal(FARGS
, "Could not allocate enough memory for hbmap");
226 for (j
= 0; (j
> hb
->a
.nra
); j
++)
228 hb
->hbmap
[i
][j
] = NULL
;
233 static void add_frames(t_hbdata
*hb
, int nframes
)
235 if (nframes
>= hb
->max_frames
)
237 hb
->max_frames
+= 4096;
238 srenew(hb
->time
, hb
->max_frames
);
239 srenew(hb
->nhb
, hb
->max_frames
);
240 srenew(hb
->ndist
, hb
->max_frames
);
241 srenew(hb
->n_bound
, hb
->max_frames
);
242 srenew(hb
->nhx
, hb
->max_frames
);
245 srenew(hb
->danr
, hb
->max_frames
);
248 hb
->nframes
= nframes
;
251 #define OFFSET(frame) (frame / 32)
252 #define MASK(frame) (1 << (frame % 32))
254 static void _set_hb(unsigned int hbexist
[], unsigned int frame
, gmx_bool bValue
)
258 hbexist
[OFFSET(frame
)] |= MASK(frame
);
262 hbexist
[OFFSET(frame
)] &= ~MASK(frame
);
266 static gmx_bool
is_hb(unsigned int hbexist
[], int frame
)
268 return ((hbexist
[OFFSET(frame
)] & MASK(frame
)) != 0) ? 1 : 0;
271 static void set_hb(t_hbdata
*hb
, int id
, int ih
, int ia
, int frame
, int ihb
)
273 unsigned int *ghptr
= NULL
;
277 ghptr
= hb
->hbmap
[id
][ia
]->h
[ih
];
279 else if (ihb
== hbDist
)
281 ghptr
= hb
->hbmap
[id
][ia
]->g
[ih
];
285 gmx_fatal(FARGS
, "Incomprehensible iValue %d in set_hb", ihb
);
288 _set_hb(ghptr
, frame
-hb
->hbmap
[id
][ia
]->n0
, TRUE
);
291 static void add_ff(t_hbdata
*hbd
, int id
, int h
, int ia
, int frame
, int ihb
)
294 t_hbond
*hb
= hbd
->hbmap
[id
][ia
];
295 int maxhydro
= std::min(hbd
->maxhydro
, hbd
->d
.nhydro
[id
]);
296 int wlen
= hbd
->wordlen
;
302 hb
->maxframes
= delta
;
303 for (i
= 0; (i
< maxhydro
); i
++)
305 snew(hb
->h
[i
], hb
->maxframes
/wlen
);
306 snew(hb
->g
[i
], hb
->maxframes
/wlen
);
311 hb
->nframes
= frame
-hb
->n0
;
312 /* We need a while loop here because hbonds may be returning
315 while (hb
->nframes
>= hb
->maxframes
)
317 n
= hb
->maxframes
+ delta
;
318 for (i
= 0; (i
< maxhydro
); i
++)
320 srenew(hb
->h
[i
], n
/wlen
);
321 srenew(hb
->g
[i
], n
/wlen
);
322 for (j
= hb
->maxframes
/wlen
; (j
< n
/wlen
); j
++)
334 set_hb(hbd
, id
, h
, ia
, frame
, ihb
);
339 static void inc_nhbonds(t_donors
*ddd
, int d
, int h
)
342 int dptr
= ddd
->dptr
[d
];
344 for (j
= 0; (j
< ddd
->nhydro
[dptr
]); j
++)
346 if (ddd
->hydro
[dptr
][j
] == h
)
348 ddd
->nhbonds
[dptr
][j
]++;
352 if (j
== ddd
->nhydro
[dptr
])
354 gmx_fatal(FARGS
, "No such hydrogen %d on donor %d\n", h
+1, d
+1);
358 static int _acceptor_index(t_acceptors
*a
, int grp
, int i
,
359 const char *file
, int line
)
363 if (a
->grp
[ai
] != grp
)
367 fprintf(debug
, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
368 ai
, a
->grp
[ai
], grp
, file
, line
);
377 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
379 static int _donor_index(t_donors
*d
, int grp
, int i
, const char *file
, int line
)
388 if (d
->grp
[di
] != grp
)
392 fprintf(debug
, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
393 di
, d
->grp
[di
], grp
, file
, line
);
402 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
404 static gmx_bool
isInterchangable(t_hbdata
*hb
, int d
, int a
, int grpa
, int grpd
)
406 /* g_hbond doesn't allow overlapping groups */
412 donor_index(&hb
->d
, grpd
, a
) != NOTSET
413 && acceptor_index(&hb
->a
, grpa
, d
) != NOTSET
;
417 static void add_hbond(t_hbdata
*hb
, int d
, int a
, int h
, int grpd
, int grpa
,
418 int frame
, gmx_bool bMerge
, int ihb
, gmx_bool bContact
)
421 gmx_bool daSwap
= FALSE
;
423 if ((id
= hb
->d
.dptr
[d
]) == NOTSET
)
425 gmx_fatal(FARGS
, "No donor atom %d", d
+1);
427 else if (grpd
!= hb
->d
.grp
[id
])
429 gmx_fatal(FARGS
, "Inconsistent donor groups, %d iso %d, atom %d",
430 grpd
, hb
->d
.grp
[id
], d
+1);
432 if ((ia
= hb
->a
.aptr
[a
]) == NOTSET
)
434 gmx_fatal(FARGS
, "No acceptor atom %d", a
+1);
436 else if (grpa
!= hb
->a
.grp
[ia
])
438 gmx_fatal(FARGS
, "Inconsistent acceptor groups, %d iso %d, atom %d",
439 grpa
, hb
->a
.grp
[ia
], a
+1);
445 if (isInterchangable(hb
, d
, a
, grpd
, grpa
) && d
> a
)
446 /* Then swap identity so that the id of d is lower then that of a.
448 * This should really be redundant by now, as is_hbond() now ought to return
449 * hbNo in the cases where this conditional is TRUE. */
456 /* Now repeat donor/acc check. */
457 if ((id
= hb
->d
.dptr
[d
]) == NOTSET
)
459 gmx_fatal(FARGS
, "No donor atom %d", d
+1);
461 else if (grpd
!= hb
->d
.grp
[id
])
463 gmx_fatal(FARGS
, "Inconsistent donor groups, %d iso %d, atom %d",
464 grpd
, hb
->d
.grp
[id
], d
+1);
466 if ((ia
= hb
->a
.aptr
[a
]) == NOTSET
)
468 gmx_fatal(FARGS
, "No acceptor atom %d", a
+1);
470 else if (grpa
!= hb
->a
.grp
[ia
])
472 gmx_fatal(FARGS
, "Inconsistent acceptor groups, %d iso %d, atom %d",
473 grpa
, hb
->a
.grp
[ia
], a
+1);
480 /* Loop over hydrogens to find which hydrogen is in this particular HB */
481 if ((ihb
== hbHB
) && !bMerge
&& !bContact
)
483 for (k
= 0; (k
< hb
->d
.nhydro
[id
]); k
++)
485 if (hb
->d
.hydro
[id
][k
] == h
)
490 if (k
== hb
->d
.nhydro
[id
])
492 gmx_fatal(FARGS
, "Donor %d does not have hydrogen %d (a = %d)",
508 if (hb
->hbmap
[id
][ia
] == NULL
)
510 snew(hb
->hbmap
[id
][ia
], 1);
511 snew(hb
->hbmap
[id
][ia
]->h
, hb
->maxhydro
);
512 snew(hb
->hbmap
[id
][ia
]->g
, hb
->maxhydro
);
514 add_ff(hb
, id
, k
, ia
, frame
, ihb
);
516 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
520 /* Strange construction with frame >=0 is a relic from old code
521 * for selected hbond analysis. It may be necessary again if that
522 * is made to work again.
526 hh
= hb
->hbmap
[id
][ia
]->history
[k
];
532 hb
->hbmap
[id
][ia
]->history
[k
] = hh
| 2;
543 hb
->hbmap
[id
][ia
]->history
[k
] = hh
| 1;
567 if (bMerge
&& daSwap
)
569 h
= hb
->d
.hydro
[id
][0];
571 /* Increment number if HBonds per H */
572 if (ihb
== hbHB
&& !bContact
)
574 inc_nhbonds(&(hb
->d
), d
, h
);
578 static char *mkatomname(t_atoms
*atoms
, int i
)
583 rnr
= atoms
->atom
[i
].resind
;
584 sprintf(buf
, "%4s%d%-4s",
585 *atoms
->resinfo
[rnr
].name
, atoms
->resinfo
[rnr
].nr
, *atoms
->atomname
[i
]);
590 static void gen_datable(int *index
, int isize
, unsigned char *datable
, int natoms
)
592 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
595 for (i
= 0; i
< isize
; i
++)
597 if (index
[i
] >= natoms
)
599 gmx_fatal(FARGS
, "Atom has index %d larger than number of atoms %d.", index
[i
], natoms
);
601 datable
[index
[i
]] |= INGROUP
;
605 static void clear_datable_grp(unsigned char *datable
, int size
)
607 /* Clears group information from the table */
609 const char mask
= !(char)INGROUP
;
612 for (i
= 0; i
< size
; i
++)
619 static void add_acc(t_acceptors
*a
, int ia
, int grp
)
621 if (a
->nra
>= a
->max_nra
)
624 srenew(a
->acc
, a
->max_nra
);
625 srenew(a
->grp
, a
->max_nra
);
627 a
->grp
[a
->nra
] = grp
;
628 a
->acc
[a
->nra
++] = ia
;
631 static void search_acceptors(t_topology
*top
, int isize
,
632 int *index
, t_acceptors
*a
, int grp
,
634 gmx_bool bContact
, gmx_bool bDoIt
, unsigned char *datable
)
640 for (i
= 0; (i
< isize
); i
++)
644 (((*top
->atoms
.atomname
[n
])[0] == 'O') ||
645 (bNitAcc
&& ((*top
->atoms
.atomname
[n
])[0] == 'N')))) &&
648 datable
[n
] |= ACC
; /* set the atom's acceptor flag in datable. */
653 snew(a
->aptr
, top
->atoms
.nr
);
654 for (i
= 0; (i
< top
->atoms
.nr
); i
++)
658 for (i
= 0; (i
< a
->nra
); i
++)
660 a
->aptr
[a
->acc
[i
]] = i
;
664 static void add_h2d(int id
, int ih
, t_donors
*ddd
)
668 for (i
= 0; (i
< ddd
->nhydro
[id
]); i
++)
670 if (ddd
->hydro
[id
][i
] == ih
)
672 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
677 if (i
== ddd
->nhydro
[id
])
679 if (ddd
->nhydro
[id
] >= MAXHYDRO
)
681 gmx_fatal(FARGS
, "Donor %d has more than %d hydrogens!",
682 ddd
->don
[id
], MAXHYDRO
);
684 ddd
->hydro
[id
][i
] = ih
;
689 static void add_dh(t_donors
*ddd
, int id
, int ih
, int grp
, unsigned char *datable
)
693 if (!datable
|| ISDON(datable
[id
]))
695 if (ddd
->dptr
[id
] == NOTSET
) /* New donor */
707 if (ddd
->nrd
>= ddd
->max_nrd
)
710 srenew(ddd
->don
, ddd
->max_nrd
);
711 srenew(ddd
->nhydro
, ddd
->max_nrd
);
712 srenew(ddd
->hydro
, ddd
->max_nrd
);
713 srenew(ddd
->nhbonds
, ddd
->max_nrd
);
714 srenew(ddd
->grp
, ddd
->max_nrd
);
716 ddd
->don
[ddd
->nrd
] = id
;
717 ddd
->nhydro
[ddd
->nrd
] = 0;
718 ddd
->grp
[ddd
->nrd
] = grp
;
730 printf("Warning: Atom %d is not in the d/a-table!\n", id
);
734 static void search_donors(t_topology
*top
, int isize
, int *index
,
735 t_donors
*ddd
, int grp
, gmx_bool bContact
, gmx_bool bDoIt
,
736 unsigned char *datable
)
739 t_functype func_type
;
740 t_ilist
*interaction
;
745 snew(ddd
->dptr
, top
->atoms
.nr
);
746 for (i
= 0; (i
< top
->atoms
.nr
); i
++)
748 ddd
->dptr
[i
] = NOTSET
;
756 for (i
= 0; (i
< isize
); i
++)
758 datable
[index
[i
]] |= DON
;
759 add_dh(ddd
, index
[i
], -1, grp
, datable
);
765 for (func_type
= 0; (func_type
< F_NRE
); func_type
++)
767 interaction
= &(top
->idef
.il
[func_type
]);
768 if (func_type
== F_POSRES
|| func_type
== F_FBPOSRES
)
770 /* The ilist looks strange for posre. Bug in grompp?
771 * We don't need posre interactions for hbonds anyway.*/
774 for (i
= 0; i
< interaction
->nr
;
775 i
+= interaction_function
[top
->idef
.functype
[interaction
->iatoms
[i
]]].nratoms
+1)
778 if (func_type
!= top
->idef
.functype
[interaction
->iatoms
[i
]])
780 fprintf(stderr
, "Error in func_type %s",
781 interaction_function
[func_type
].longname
);
785 /* check out this functype */
786 if (func_type
== F_SETTLE
)
788 nr1
= interaction
->iatoms
[i
+1];
789 nr2
= interaction
->iatoms
[i
+2];
790 nr3
= interaction
->iatoms
[i
+3];
792 if (ISINGRP(datable
[nr1
]))
794 if (ISINGRP(datable
[nr2
]))
797 add_dh(ddd
, nr1
, nr1
+1, grp
, datable
);
799 if (ISINGRP(datable
[nr3
]))
802 add_dh(ddd
, nr1
, nr1
+2, grp
, datable
);
806 else if (IS_CHEMBOND(func_type
))
808 for (j
= 0; j
< 2; j
++)
810 nr1
= interaction
->iatoms
[i
+1+j
];
811 nr2
= interaction
->iatoms
[i
+2-j
];
812 if ((*top
->atoms
.atomname
[nr1
][0] == 'H') &&
813 ((*top
->atoms
.atomname
[nr2
][0] == 'O') ||
814 (*top
->atoms
.atomname
[nr2
][0] == 'N')) &&
815 ISINGRP(datable
[nr1
]) && ISINGRP(datable
[nr2
]))
818 add_dh(ddd
, nr2
, nr1
, grp
, datable
);
825 for (func_type
= 0; func_type
< F_NRE
; func_type
++)
827 interaction
= &top
->idef
.il
[func_type
];
828 for (i
= 0; i
< interaction
->nr
;
829 i
+= interaction_function
[top
->idef
.functype
[interaction
->iatoms
[i
]]].nratoms
+1)
832 if (func_type
!= top
->idef
.functype
[interaction
->iatoms
[i
]])
834 gmx_incons("function type in search_donors");
837 if (interaction_function
[func_type
].flags
& IF_VSITE
)
839 nr1
= interaction
->iatoms
[i
+1];
840 if (*top
->atoms
.atomname
[nr1
][0] == 'H')
844 while (!stop
&& ( *top
->atoms
.atomname
[nr2
][0] == 'H'))
855 if (!stop
&& ( ( *top
->atoms
.atomname
[nr2
][0] == 'O') ||
856 ( *top
->atoms
.atomname
[nr2
][0] == 'N') ) &&
857 ISINGRP(datable
[nr1
]) && ISINGRP(datable
[nr2
]))
860 add_dh(ddd
, nr2
, nr1
, grp
, datable
);
870 static t_gridcell
***init_grid(gmx_bool bBox
, rvec box
[], real rcut
, ivec ngrid
)
877 for (i
= 0; i
< DIM
; i
++)
879 ngrid
[i
] = static_cast<int>(box
[i
][i
]/(1.2*rcut
));
883 if (!bBox
|| (ngrid
[XX
] < 3) || (ngrid
[YY
] < 3) || (ngrid
[ZZ
] < 3) )
885 for (i
= 0; i
< DIM
; i
++)
892 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
893 ngrid
[XX
], ngrid
[YY
], ngrid
[ZZ
], rcut
);
895 snew(grid
, ngrid
[ZZ
]);
896 for (z
= 0; z
< ngrid
[ZZ
]; z
++)
898 snew((grid
)[z
], ngrid
[YY
]);
899 for (y
= 0; y
< ngrid
[YY
]; y
++)
901 snew((grid
)[z
][y
], ngrid
[XX
]);
907 static void reset_nhbonds(t_donors
*ddd
)
911 for (i
= 0; (i
< ddd
->nrd
); i
++)
913 for (j
= 0; (j
< MAXHH
); j
++)
915 ddd
->nhbonds
[i
][j
] = 0;
920 void pbc_correct_gem(rvec dx
, matrix box
, rvec hbox
);
921 void pbc_in_gridbox(rvec dx
, matrix box
);
923 static void build_grid(t_hbdata
*hb
, rvec x
[], rvec xshell
,
924 gmx_bool bBox
, matrix box
, rvec hbox
,
925 real rcut
, real rshell
,
926 ivec ngrid
, t_gridcell
***grid
)
928 int i
, m
, gr
, xi
, yi
, zi
, nr
;
931 rvec invdelta
, dshell
;
933 gmx_bool bDoRshell
, bInShell
, bAcc
;
938 bDoRshell
= (rshell
> 0);
939 rshell2
= sqr(rshell
);
942 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
944 for (m
= 0; m
< DIM
; m
++)
946 hbox
[m
] = box
[m
][m
]*0.5;
949 invdelta
[m
] = ngrid
[m
]/box
[m
][m
];
950 if (1/invdelta
[m
] < rcut
)
952 gmx_fatal(FARGS
, "Your computational box has shrunk too much.\n"
953 "%s can not handle this situation, sorry.\n",
954 gmx::getProgramContext().displayName());
966 /* resetting atom counts */
967 for (gr
= 0; (gr
< grNR
); gr
++)
969 for (zi
= 0; zi
< ngrid
[ZZ
]; zi
++)
971 for (yi
= 0; yi
< ngrid
[YY
]; yi
++)
973 for (xi
= 0; xi
< ngrid
[XX
]; xi
++)
975 grid
[zi
][yi
][xi
].d
[gr
].nr
= 0;
976 grid
[zi
][yi
][xi
].a
[gr
].nr
= 0;
982 /* put atoms in grid cells */
983 for (bAcc
= FALSE
; (bAcc
<= TRUE
); bAcc
++)
996 for (i
= 0; (i
< nr
); i
++)
998 /* check if we are inside the shell */
999 /* if bDoRshell=FALSE then bInShell=TRUE always */
1004 rvec_sub(x
[ad
[i
]], xshell
, dshell
);
1007 gmx_bool bDone
= FALSE
;
1011 for (m
= DIM
-1; m
>= 0 && bInShell
; m
--)
1013 if (dshell
[m
] < -hbox
[m
])
1016 rvec_inc(dshell
, box
[m
]);
1018 if (dshell
[m
] >= hbox
[m
])
1021 dshell
[m
] -= 2*hbox
[m
];
1025 for (m
= DIM
-1; m
>= 0 && bInShell
; m
--)
1027 /* if we're outside the cube, we're outside the sphere also! */
1028 if ( (dshell
[m
] > rshell
) || (-dshell
[m
] > rshell
) )
1034 /* if we're inside the cube, check if we're inside the sphere */
1037 bInShell
= norm2(dshell
) < rshell2
;
1045 pbc_in_gridbox(x
[ad
[i
]], box
);
1047 for (m
= DIM
-1; m
>= 0; m
--)
1048 { /* determine grid index of atom */
1049 grididx
[m
] = static_cast<int>(x
[ad
[i
]][m
]*invdelta
[m
]);
1050 grididx
[m
] = (grididx
[m
]+ngrid
[m
]) % ngrid
[m
];
1057 range_check(gx
, 0, ngrid
[XX
]);
1058 range_check(gy
, 0, ngrid
[YY
]);
1059 range_check(gz
, 0, ngrid
[ZZ
]);
1063 /* add atom to grid cell */
1066 newgrid
= &(grid
[gz
][gy
][gx
].a
[gr
]);
1070 newgrid
= &(grid
[gz
][gy
][gx
].d
[gr
]);
1072 if (newgrid
->nr
>= newgrid
->maxnr
)
1074 newgrid
->maxnr
+= 10;
1075 DBB(newgrid
->maxnr
);
1076 srenew(newgrid
->atoms
, newgrid
->maxnr
);
1079 newgrid
->atoms
[newgrid
->nr
] = ad
[i
];
1087 static void count_da_grid(ivec ngrid
, t_gridcell
***grid
, t_icell danr
)
1091 for (gr
= 0; (gr
< grNR
); gr
++)
1094 for (zi
= 0; zi
< ngrid
[ZZ
]; zi
++)
1096 for (yi
= 0; yi
< ngrid
[YY
]; yi
++)
1098 for (xi
= 0; xi
< ngrid
[XX
]; xi
++)
1100 danr
[gr
] += grid
[zi
][yi
][xi
].d
[gr
].nr
;
1108 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1109 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1110 * With a triclinic box all loops are 3 long, except when a cell is
1111 * located next to one of the box edges which is not parallel to the
1112 * x/y-plane, in that case all cells in a line or layer are searched.
1113 * This could be implemented slightly more efficient, but the code
1114 * would get much more complicated.
1116 static gmx_inline gmx_bool
grid_loop_begin(int n
, int x
, gmx_bool bTric
, gmx_bool bEdge
)
1118 return ((n
== 1) ? x
: bTric
&& bEdge
? 0 : (x
-1));
1120 static gmx_inline gmx_bool
grid_loop_end(int n
, int x
, gmx_bool bTric
, gmx_bool bEdge
)
1122 return ((n
== 1) ? x
: bTric
&& bEdge
? (n
-1) : (x
+1));
1124 static gmx_inline
int grid_mod(int j
, int n
)
1129 static void dump_grid(FILE *fp
, ivec ngrid
, t_gridcell
***grid
)
1131 int gr
, x
, y
, z
, sum
[grNR
];
1133 fprintf(fp
, "grid %dx%dx%d\n", ngrid
[XX
], ngrid
[YY
], ngrid
[ZZ
]);
1134 for (gr
= 0; gr
< grNR
; gr
++)
1137 fprintf(fp
, "GROUP %d (%s)\n", gr
, grpnames
[gr
]);
1138 for (z
= 0; z
< ngrid
[ZZ
]; z
+= 2)
1140 fprintf(fp
, "Z=%d,%d\n", z
, z
+1);
1141 for (y
= 0; y
< ngrid
[YY
]; y
++)
1143 for (x
= 0; x
< ngrid
[XX
]; x
++)
1145 fprintf(fp
, "%3d", grid
[x
][y
][z
].d
[gr
].nr
);
1146 sum
[gr
] += grid
[z
][y
][x
].d
[gr
].nr
;
1147 fprintf(fp
, "%3d", grid
[x
][y
][z
].a
[gr
].nr
);
1148 sum
[gr
] += grid
[z
][y
][x
].a
[gr
].nr
;
1152 if ( (z
+1) < ngrid
[ZZ
])
1154 for (x
= 0; x
< ngrid
[XX
]; x
++)
1156 fprintf(fp
, "%3d", grid
[z
+1][y
][x
].d
[gr
].nr
);
1157 sum
[gr
] += grid
[z
+1][y
][x
].d
[gr
].nr
;
1158 fprintf(fp
, "%3d", grid
[z
+1][y
][x
].a
[gr
].nr
);
1159 sum
[gr
] += grid
[z
+1][y
][x
].a
[gr
].nr
;
1166 fprintf(fp
, "TOTALS:");
1167 for (gr
= 0; gr
< grNR
; gr
++)
1169 fprintf(fp
, " %d=%d", gr
, sum
[gr
]);
1174 /* New GMX record! 5 * in a row. Congratulations!
1175 * Sorry, only four left.
1177 static void free_grid(ivec ngrid
, t_gridcell
****grid
)
1180 t_gridcell
***g
= *grid
;
1182 for (z
= 0; z
< ngrid
[ZZ
]; z
++)
1184 for (y
= 0; y
< ngrid
[YY
]; y
++)
1194 void pbc_correct_gem(rvec dx
, matrix box
, rvec hbox
)
1197 gmx_bool bDone
= FALSE
;
1201 for (m
= DIM
-1; m
>= 0; m
--)
1203 if (dx
[m
] < -hbox
[m
])
1206 rvec_inc(dx
, box
[m
]);
1208 if (dx
[m
] >= hbox
[m
])
1211 rvec_dec(dx
, box
[m
]);
1217 void pbc_in_gridbox(rvec dx
, matrix box
)
1220 gmx_bool bDone
= FALSE
;
1224 for (m
= DIM
-1; m
>= 0; m
--)
1229 rvec_inc(dx
, box
[m
]);
1231 if (dx
[m
] >= box
[m
][m
])
1234 rvec_dec(dx
, box
[m
]);
1240 /* Added argument r2cut, changed contact and implemented
1241 * use of second cut-off.
1242 * - Erik Marklund, June 29, 2006
1244 static int is_hbond(t_hbdata
*hb
, int grpd
, int grpa
, int d
, int a
,
1245 real rcut
, real r2cut
, real ccut
,
1246 rvec x
[], gmx_bool bBox
, matrix box
, rvec hbox
,
1247 real
*d_ha
, real
*ang
, gmx_bool bDA
, int *hhh
,
1248 gmx_bool bContact
, gmx_bool bMerge
)
1251 rvec r_da
, r_ha
, r_dh
;
1252 real rc2
, r2c2
, rda2
, rha2
, ca
;
1253 gmx_bool HAinrange
= FALSE
; /* If !bDA. Needed for returning hbDist in a correct way. */
1254 gmx_bool daSwap
= FALSE
;
1261 if (((id
= donor_index(&hb
->d
, grpd
, d
)) == NOTSET
) ||
1262 (acceptor_index(&hb
->a
, grpa
, a
) == NOTSET
))
1270 rvec_sub(x
[d
], x
[a
], r_da
);
1271 /* Insert projection code here */
1273 if (bMerge
&& d
> a
&& isInterchangable(hb
, d
, a
, grpd
, grpa
))
1275 /* Then this hbond/contact will be found again, or it has already been found. */
1280 if (d
> a
&& bMerge
&& isInterchangable(hb
, d
, a
, grpd
, grpa
)) /* acceptor is also a donor and vice versa? */
1281 { /* return hbNo; */
1282 daSwap
= TRUE
; /* If so, then their history should be filed with donor and acceptor swapped. */
1284 pbc_correct_gem(r_da
, box
, hbox
);
1286 rda2
= iprod(r_da
, r_da
);
1290 if (daSwap
&& grpa
== grpd
)
1298 else if (rda2
< r2c2
)
1309 if (bDA
&& (rda2
> rc2
))
1314 for (h
= 0; (h
< hb
->d
.nhydro
[id
]); h
++)
1316 hh
= hb
->d
.hydro
[id
][h
];
1320 rvec_sub(x
[hh
], x
[a
], r_ha
);
1323 pbc_correct_gem(r_ha
, box
, hbox
);
1325 rha2
= iprod(r_ha
, r_ha
);
1328 if (bDA
|| (!bDA
&& (rha2
<= rc2
)))
1330 rvec_sub(x
[d
], x
[hh
], r_dh
);
1333 pbc_correct_gem(r_dh
, box
, hbox
);
1340 ca
= cos_angle(r_dh
, r_da
);
1341 /* if angle is smaller, cos is larger */
1345 *d_ha
= std::sqrt(bDA
? rda2
: rha2
);
1346 *ang
= std::acos(ca
);
1351 if (bDA
|| (!bDA
&& HAinrange
))
1361 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1362 * Will do some more testing before removing the function entirely.
1363 * - Erik Marklund, MAY 10 2010 */
1364 static void do_merge(t_hbdata
*hb
, int ntmp
,
1365 unsigned int htmp
[], unsigned int gtmp
[],
1366 t_hbond
*hb0
, t_hbond
*hb1
)
1368 /* Here we need to make sure we're treating periodicity in
1369 * the right way for the geminate recombination kinetics. */
1371 int m
, mm
, n00
, n01
, nn0
, nnframes
;
1373 /* Decide where to start from when merging */
1376 nn0
= std::min(n00
, n01
);
1377 nnframes
= std::max(n00
+ hb0
->nframes
, n01
+ hb1
->nframes
) - nn0
;
1378 /* Initiate tmp arrays */
1379 for (m
= 0; (m
< ntmp
); m
++)
1384 /* Fill tmp arrays with values due to first HB */
1385 /* Once again '<' had to be replaced with '<='
1386 to catch the last frame in which the hbond
1388 - Erik Marklund, June 1, 2006 */
1389 for (m
= 0; (m
<= hb0
->nframes
); m
++)
1392 htmp
[mm
] = is_hb(hb0
->h
[0], m
);
1394 for (m
= 0; (m
<= hb0
->nframes
); m
++)
1397 gtmp
[mm
] = is_hb(hb0
->g
[0], m
);
1400 for (m
= 0; (m
<= hb1
->nframes
); m
++)
1403 htmp
[mm
] = htmp
[mm
] || is_hb(hb1
->h
[0], m
);
1404 gtmp
[mm
] = gtmp
[mm
] || is_hb(hb1
->g
[0], m
);
1406 /* Reallocate target array */
1407 if (nnframes
> hb0
->maxframes
)
1409 srenew(hb0
->h
[0], 4+nnframes
/hb
->wordlen
);
1410 srenew(hb0
->g
[0], 4+nnframes
/hb
->wordlen
);
1413 /* Copy temp array to target array */
1414 for (m
= 0; (m
<= nnframes
); m
++)
1416 _set_hb(hb0
->h
[0], m
, htmp
[m
]);
1417 _set_hb(hb0
->g
[0], m
, gtmp
[m
]);
1420 /* Set scalar variables */
1422 hb0
->maxframes
= nnframes
;
1425 static void merge_hb(t_hbdata
*hb
, gmx_bool bTwo
, gmx_bool bContact
)
1427 int i
, inrnew
, indnew
, j
, ii
, jj
, id
, ia
, ntmp
;
1428 unsigned int *htmp
, *gtmp
;
1432 indnew
= hb
->nrdist
;
1434 /* Check whether donors are also acceptors */
1435 printf("Merging hbonds with Acceptor and Donor swapped\n");
1437 ntmp
= 2*hb
->max_frames
;
1440 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
1442 fprintf(stderr
, "\r%d/%d", i
+1, hb
->d
.nrd
);
1444 ii
= hb
->a
.aptr
[id
];
1445 for (j
= 0; (j
< hb
->a
.nra
); j
++)
1448 jj
= hb
->d
.dptr
[ia
];
1449 if ((id
!= ia
) && (ii
!= NOTSET
) && (jj
!= NOTSET
) &&
1450 (!bTwo
|| (bTwo
&& (hb
->d
.grp
[i
] != hb
->a
.grp
[j
]))))
1452 hb0
= hb
->hbmap
[i
][j
];
1453 hb1
= hb
->hbmap
[jj
][ii
];
1454 if (hb0
&& hb1
&& ISHB(hb0
->history
[0]) && ISHB(hb1
->history
[0]))
1456 do_merge(hb
, ntmp
, htmp
, gtmp
, hb0
, hb1
);
1457 if (ISHB(hb1
->history
[0]))
1461 else if (ISDIST(hb1
->history
[0]))
1468 gmx_incons("No contact history");
1472 gmx_incons("Neither hydrogen bond nor distance");
1478 hb1
->history
[0] = hbNo
;
1483 fprintf(stderr
, "\n");
1484 printf("- Reduced number of hbonds from %d to %d\n", hb
->nrhb
, inrnew
);
1485 printf("- Reduced number of distances from %d to %d\n", hb
->nrdist
, indnew
);
1487 hb
->nrdist
= indnew
;
1492 static void do_nhb_dist(FILE *fp
, t_hbdata
*hb
, real t
)
1494 int i
, j
, k
, n_bound
[MAXHH
], nbtot
;
1496 /* Set array to 0 */
1497 for (k
= 0; (k
< MAXHH
); k
++)
1501 /* Loop over possible donors */
1502 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
1504 for (j
= 0; (j
< hb
->d
.nhydro
[i
]); j
++)
1506 n_bound
[hb
->d
.nhbonds
[i
][j
]]++;
1509 fprintf(fp
, "%12.5e", t
);
1511 for (k
= 0; (k
< MAXHH
); k
++)
1513 fprintf(fp
, " %8d", n_bound
[k
]);
1514 nbtot
+= n_bound
[k
]*k
;
1516 fprintf(fp
, " %8d\n", nbtot
);
1519 static void do_hblife(const char *fn
, t_hbdata
*hb
, gmx_bool bMerge
, gmx_bool bContact
,
1520 const gmx_output_env_t
*oenv
)
1523 const char *leg
[] = { "p(t)", "t p(t)" };
1525 int i
, j
, j0
, k
, m
, nh
, ihb
, ohb
, nhydro
, ndump
= 0;
1526 int nframes
= hb
->nframes
;
1529 double sum
, integral
;
1532 snew(h
, hb
->maxhydro
);
1533 snew(histo
, nframes
+1);
1534 /* Total number of hbonds analyzed here */
1535 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
1537 for (k
= 0; (k
< hb
->a
.nra
); k
++)
1539 hbh
= hb
->hbmap
[i
][k
];
1557 for (m
= 0; (m
< hb
->maxhydro
); m
++)
1561 h
[nhydro
++] = bContact
? hbh
->g
[m
] : hbh
->h
[m
];
1565 for (nh
= 0; (nh
< nhydro
); nh
++)
1570 for (j
= 0; (j
<= hbh
->nframes
); j
++)
1572 ihb
= is_hb(h
[nh
], j
);
1573 if (debug
&& (ndump
< 10))
1575 fprintf(debug
, "%5d %5d\n", j
, ihb
);
1595 fprintf(stderr
, "\n");
1598 fp
= xvgropen(fn
, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv
), "()", oenv
);
1602 fp
= xvgropen(fn
, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv
), "()",
1606 xvgr_legend(fp
, asize(leg
), leg
, oenv
);
1608 while ((j0
> 0) && (histo
[j0
] == 0))
1613 for (i
= 0; (i
<= j0
); i
++)
1617 dt
= hb
->time
[1]-hb
->time
[0];
1620 for (i
= 1; (i
<= j0
); i
++)
1622 t
= hb
->time
[i
] - hb
->time
[0] - 0.5*dt
;
1623 x1
= t
*histo
[i
]/sum
;
1624 fprintf(fp
, "%8.3f %10.3e %10.3e\n", t
, histo
[i
]/sum
, x1
);
1629 printf("%s lifetime = %.2f ps\n", bContact
? "Contact" : "HB", integral
);
1630 printf("Note that the lifetime obtained in this manner is close to useless\n");
1631 printf("Use the -ac option instead and check the Forward lifetime\n");
1632 please_cite(stdout
, "Spoel2006b");
1637 static void dump_ac(t_hbdata
*hb
, gmx_bool oneHB
, int nDump
)
1640 int i
, j
, k
, m
, nd
, ihb
, idist
;
1641 int nframes
= hb
->nframes
;
1649 fp
= gmx_ffopen("debug-ac.xvg", "w");
1650 for (j
= 0; (j
< nframes
); j
++)
1652 fprintf(fp
, "%10.3f", hb
->time
[j
]);
1653 for (i
= nd
= 0; (i
< hb
->d
.nrd
) && (nd
< nDump
); i
++)
1655 for (k
= 0; (k
< hb
->a
.nra
) && (nd
< nDump
); k
++)
1659 hbh
= hb
->hbmap
[i
][k
];
1664 ihb
= is_hb(hbh
->h
[0], j
);
1665 idist
= is_hb(hbh
->g
[0], j
);
1671 for (m
= 0; (m
< hb
->maxhydro
) && !ihb
; m
++)
1673 ihb
= ihb
|| ((hbh
->h
[m
]) && is_hb(hbh
->h
[m
], j
));
1674 idist
= idist
|| ((hbh
->g
[m
]) && is_hb(hbh
->g
[m
], j
));
1676 /* This is not correct! */
1677 /* What isn't correct? -Erik M */
1682 fprintf(fp
, " %1d-%1d", ihb
, idist
);
1692 static real
calc_dg(real tau
, real temp
)
1703 return kbt
*std::log(kbt
*tau
/PLANCK
);
1708 int n0
, n1
, nparams
, ndelta
;
1710 real
*t
, *ct
, *nt
, *kt
, *sigma_ct
, *sigma_nt
, *sigma_kt
;
1713 static real
compute_weighted_rates(int n
, real t
[], real ct
[], real nt
[],
1714 real kt
[], real sigma_ct
[], real sigma_nt
[],
1715 real sigma_kt
[], real
*k
, real
*kp
,
1716 real
*sigma_k
, real
*sigma_kp
,
1722 real kkk
= 0, kkp
= 0, kk2
= 0, kp2
= 0, chi2
;
1727 for (i
= 0; (i
< n
); i
++)
1729 if (t
[i
] >= fit_start
)
1742 tl
.sigma_ct
= sigma_ct
;
1743 tl
.sigma_nt
= sigma_nt
;
1744 tl
.sigma_kt
= sigma_kt
;
1748 chi2
= 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
1750 *kp
= tl
.kkk
[1] = *kp
;
1752 for (j
= 0; (j
< NK
); j
++)
1756 kk2
+= sqr(tl
.kkk
[0]);
1757 kp2
+= sqr(tl
.kkk
[1]);
1760 *sigma_k
= std::sqrt(kk2
/NK
- sqr(kkk
/NK
));
1761 *sigma_kp
= std::sqrt(kp2
/NK
- sqr(kkp
/NK
));
1766 void analyse_corr(int n
, real t
[], real ct
[], real nt
[], real kt
[],
1767 real sigma_ct
[], real sigma_nt
[], real sigma_kt
[],
1768 real fit_start
, real temp
)
1771 real k
= 1, kp
= 1, kow
= 1;
1772 real Q
= 0, chi2
, tau_hb
, dtau
, tau_rlx
, e_1
, sigma_k
, sigma_kp
, ddg
;
1773 double tmp
, sn2
= 0, sc2
= 0, sk2
= 0, scn
= 0, sck
= 0, snk
= 0;
1774 gmx_bool bError
= (sigma_ct
!= NULL
) && (sigma_nt
!= NULL
) && (sigma_kt
!= NULL
);
1776 for (i0
= 0; (i0
< n
-2) && ((t
[i0
]-t
[0]) < fit_start
); i0
++)
1782 for (i
= i0
; (i
< n
); i
++)
1791 printf("Hydrogen bond thermodynamics at T = %g K\n", temp
);
1792 tmp
= (sn2
*sc2
-sqr(scn
));
1793 if ((tmp
> 0) && (sn2
> 0))
1795 k
= (sn2
*sck
-scn
*snk
)/tmp
;
1796 kp
= (k
*scn
-snk
)/sn2
;
1800 for (i
= i0
; (i
< n
); i
++)
1802 chi2
+= sqr(k
*ct
[i
]-kp
*nt
[i
]-kt
[i
]);
1804 compute_weighted_rates(n
, t
, ct
, nt
, kt
, sigma_ct
, sigma_nt
,
1806 &sigma_k
, &sigma_kp
, fit_start
);
1807 Q
= 0; /* quality_of_fit(chi2, 2);*/
1808 ddg
= BOLTZ
*temp
*sigma_k
/k
;
1809 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
1811 printf("The Rate and Delta G are followed by an error estimate\n");
1812 printf("----------------------------------------------------------\n"
1813 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
1814 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1815 k
, sigma_k
, 1/k
, calc_dg(1/k
, temp
), ddg
);
1816 ddg
= BOLTZ
*temp
*sigma_kp
/kp
;
1817 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1818 kp
, sigma_kp
, 1/kp
, calc_dg(1/kp
, temp
), ddg
);
1823 for (i
= i0
; (i
< n
); i
++)
1825 chi2
+= sqr(k
*ct
[i
]-kp
*nt
[i
]-kt
[i
]);
1827 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
1829 printf("--------------------------------------------------\n"
1830 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
1831 printf("Forward %10.3f %8.3f %10.3f %10g\n",
1832 k
, 1/k
, calc_dg(1/k
, temp
), chi2
);
1833 printf("Backward %10.3f %8.3f %10.3f\n",
1834 kp
, 1/kp
, calc_dg(1/kp
, temp
));
1840 printf("One-way %10.3f %s%8.3f %10.3f\n",
1841 kow
, bError
? " " : "", 1/kow
, calc_dg(1/kow
, temp
));
1845 printf(" - Numerical problems computing HB thermodynamics:\n"
1846 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
1847 sc2
, sn2
, sk2
, sck
, snk
, scn
);
1849 /* Determine integral of the correlation function */
1850 tau_hb
= evaluate_integral(n
, t
, ct
, NULL
, (t
[n
-1]-t
[0])/2, &dtau
);
1851 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb
,
1852 bError
? " " : "", tau_hb
, calc_dg(tau_hb
, temp
));
1853 e_1
= std::exp(-1.0);
1854 for (i
= 0; (i
< n
-2); i
++)
1856 if ((ct
[i
] > e_1
) && (ct
[i
+1] <= e_1
))
1863 /* Determine tau_relax from linear interpolation */
1864 tau_rlx
= t
[i
]-t
[0] + (e_1
-ct
[i
])*(t
[i
+1]-t
[i
])/(ct
[i
+1]-ct
[i
]);
1865 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx
,
1866 tau_rlx
, bError
? " " : "",
1867 calc_dg(tau_rlx
, temp
));
1872 printf("Correlation functions too short to compute thermodynamics\n");
1876 void compute_derivative(int nn
, real x
[], real y
[], real dydx
[])
1880 /* Compute k(t) = dc(t)/dt */
1881 for (j
= 1; (j
< nn
-1); j
++)
1883 dydx
[j
] = (y
[j
+1]-y
[j
-1])/(x
[j
+1]-x
[j
-1]);
1885 /* Extrapolate endpoints */
1886 dydx
[0] = 2*dydx
[1] - dydx
[2];
1887 dydx
[nn
-1] = 2*dydx
[nn
-2] - dydx
[nn
-3];
1890 static void normalizeACF(real
*ct
, real
*gt
, int nhb
, int len
)
1892 real ct_fac
, gt_fac
= 0;
1895 /* Xu and Berne use the same normalization constant */
1903 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac
, gt_fac
);
1904 for (i
= 0; i
< len
; i
++)
1914 static void do_hbac(const char *fn
, t_hbdata
*hb
,
1915 int nDump
, gmx_bool bMerge
, gmx_bool bContact
, real fit_start
,
1916 real temp
, gmx_bool R2
, const gmx_output_env_t
*oenv
,
1920 int i
, j
, k
, m
, ihb
, idist
, n2
, nn
;
1922 const char *legLuzar
[] = {
1923 "Ac\\sfin sys\\v{}\\z{}(t)",
1925 "Cc\\scontact,hb\\v{}\\z{}(t)",
1926 "-dAc\\sfs\\v{}\\z{}/dt"
1928 gmx_bool bNorm
= FALSE
, bOMP
= FALSE
;
1930 real
*rhbex
= NULL
, *ht
, *gt
, *ght
, *dght
, *kt
;
1931 real
*ct
, tail
, tail2
, dtail
, *cct
;
1932 const real tol
= 1e-3;
1933 int nframes
= hb
->nframes
;
1934 unsigned int **h
= NULL
, **g
= NULL
;
1935 int nh
, nhbonds
, nhydro
;
1938 int *dondata
= NULL
;
1941 AC_NONE
, AC_NN
, AC_GEM
, AC_LUZAR
1950 printf("Doing autocorrelation ");
1953 printf("according to the theory of Luzar and Chandler.\n");
1956 /* build hbexist matrix in reals for autocorr */
1957 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
1959 while (n2
< nframes
)
1966 if (acType
!= AC_NN
|| bOMP
)
1968 snew(h
, hb
->maxhydro
);
1969 snew(g
, hb
->maxhydro
);
1972 /* Dump hbonds for debugging */
1973 dump_ac(hb
, bMerge
|| bContact
, nDump
);
1975 /* Total number of hbonds analyzed here */
1978 if (acType
!= AC_LUZAR
&& bOMP
)
1980 nThreads
= std::min((nThreads
<= 0) ? INT_MAX
: nThreads
, gmx_omp_get_max_threads());
1982 gmx_omp_set_num_threads(nThreads
);
1983 snew(dondata
, nThreads
);
1984 for (i
= 0; i
< nThreads
; i
++)
1988 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
1989 "Expect close to linear scaling over this donor-loop.\n", nThreads
);
1991 fprintf(stderr
, "Donors: [thread no]\n");
1994 for (i
= 0; i
< nThreads
; i
++)
1996 snprintf(tmpstr
, 7, "[%i]", i
);
1997 fprintf(stderr
, "%-7s", tmpstr
);
2000 fprintf(stderr
, "\n");
2015 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2017 for (k
= 0; (k
< hb
->a
.nra
); k
++)
2020 hbh
= hb
->hbmap
[i
][k
];
2024 if (bMerge
|| bContact
)
2026 if (ISHB(hbh
->history
[0]))
2035 for (m
= 0; (m
< hb
->maxhydro
); m
++)
2037 if (bContact
? ISDIST(hbh
->history
[m
]) : ISHB(hbh
->history
[m
]))
2039 g
[nhydro
] = hbh
->g
[m
];
2040 h
[nhydro
] = hbh
->h
[m
];
2046 int nf
= hbh
->nframes
;
2047 for (nh
= 0; (nh
< nhydro
); nh
++)
2049 int nrint
= bContact
? hb
->nrdist
: hb
->nrhb
;
2050 if ((((nhbonds
+1) % 10) == 0) || (nhbonds
+1 == nrint
))
2052 fprintf(stderr
, "\rACF %d/%d", nhbonds
+1, nrint
);
2055 for (j
= 0; (j
< nframes
); j
++)
2059 ihb
= is_hb(h
[nh
], j
);
2060 idist
= is_hb(g
[nh
], j
);
2067 /* For contacts: if a second cut-off is provided, use it,
2068 * otherwise use g(t) = 1-h(t) */
2069 if (!R2
&& bContact
)
2075 gt
[j
] = idist
*(1-ihb
);
2081 /* The autocorrelation function is normalized after summation only */
2082 low_do_autocorr(NULL
, oenv
, NULL
, nframes
, 1, -1, &rhbex
, hb
->time
[1]-hb
->time
[0],
2083 eacNormal
, 1, FALSE
, bNorm
, FALSE
, 0, -1, 0);
2085 /* Cross correlation analysis for thermodynamics */
2086 for (j
= nframes
; (j
< n2
); j
++)
2092 cross_corr(n2
, ht
, gt
, dght
);
2094 for (j
= 0; (j
< nn
); j
++)
2103 fprintf(stderr
, "\n");
2106 normalizeACF(ct
, ght
, static_cast<int>(nhb
), nn
);
2108 /* Determine tail value for statistics */
2111 for (j
= nn
/2; (j
< nn
); j
++)
2114 tail2
+= ct
[j
]*ct
[j
];
2116 tail
/= (nn
- nn
/2);
2117 tail2
/= (nn
- nn
/2);
2118 dtail
= std::sqrt(tail2
-tail
*tail
);
2120 /* Check whether the ACF is long enough */
2123 printf("\nWARNING: Correlation function is probably not long enough\n"
2124 "because the standard deviation in the tail of C(t) > %g\n"
2125 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2128 for (j
= 0; (j
< nn
); j
++)
2131 ct
[j
] = (cct
[j
]-tail
)/(1-tail
);
2133 /* Compute negative derivative k(t) = -dc(t)/dt */
2134 compute_derivative(nn
, hb
->time
, ct
, kt
);
2135 for (j
= 0; (j
< nn
); j
++)
2143 fp
= xvgropen(fn
, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv
), "C(t)", oenv
);
2147 fp
= xvgropen(fn
, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv
), "C(t)", oenv
);
2149 xvgr_legend(fp
, asize(legLuzar
), legLuzar
, oenv
);
2152 for (j
= 0; (j
< nn
); j
++)
2154 fprintf(fp
, "%10g %10g %10g %10g %10g\n",
2155 hb
->time
[j
]-hb
->time
[0], ct
[j
], cct
[j
], ght
[j
], kt
[j
]);
2159 analyse_corr(nn
, hb
->time
, ct
, ght
, kt
, NULL
, NULL
, NULL
,
2162 do_view(oenv
, fn
, NULL
);
2173 static void init_hbframe(t_hbdata
*hb
, int nframes
, real t
)
2177 hb
->time
[nframes
] = t
;
2178 hb
->nhb
[nframes
] = 0;
2179 hb
->ndist
[nframes
] = 0;
2180 for (i
= 0; (i
< max_hx
); i
++)
2182 hb
->nhx
[nframes
][i
] = 0;
2186 static FILE *open_donor_properties_file(const char *fn
,
2188 const gmx_output_env_t
*oenv
)
2191 const char *leg
[] = { "Nbound", "Nfree" };
2198 fp
= xvgropen(fn
, "Donor properties", output_env_get_xvgr_tlabel(oenv
), "Number", oenv
);
2199 xvgr_legend(fp
, asize(leg
), leg
, oenv
);
2204 static void analyse_donor_properties(FILE *fp
, t_hbdata
*hb
, int nframes
, real t
)
2206 int i
, j
, k
, nbound
, nb
, nhtot
;
2214 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2216 for (k
= 0; (k
< hb
->d
.nhydro
[i
]); k
++)
2220 for (j
= 0; (j
< hb
->a
.nra
) && (nb
== 0); j
++)
2222 if (hb
->hbmap
[i
][j
] && hb
->hbmap
[i
][j
]->h
[k
] &&
2223 is_hb(hb
->hbmap
[i
][j
]->h
[k
], nframes
))
2231 fprintf(fp
, "%10.3e %6d %6d\n", t
, nbound
, nhtot
-nbound
);
2234 static void dump_hbmap(t_hbdata
*hb
,
2235 int nfile
, t_filenm fnm
[], gmx_bool bTwo
,
2236 gmx_bool bContact
, int isize
[], int *index
[], char *grpnames
[],
2240 int ddd
, hhh
, aaa
, i
, j
, k
, m
, grp
;
2241 char ds
[32], hs
[32], as
[32];
2244 fp
= opt2FILE("-hbn", nfile
, fnm
, "w");
2245 if (opt2bSet("-g", nfile
, fnm
))
2247 fplog
= gmx_ffopen(opt2fn("-g", nfile
, fnm
), "w");
2248 fprintf(fplog
, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
2254 for (grp
= gr0
; grp
<= (bTwo
? gr1
: gr0
); grp
++)
2256 fprintf(fp
, "[ %s ]", grpnames
[grp
]);
2257 for (i
= 0; i
< isize
[grp
]; i
++)
2259 fprintf(fp
, (i
%15) ? " " : "\n");
2260 fprintf(fp
, " %4d", index
[grp
][i
]+1);
2266 fprintf(fp
, "[ donors_hydrogens_%s ]\n", grpnames
[grp
]);
2267 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2269 if (hb
->d
.grp
[i
] == grp
)
2271 for (j
= 0; (j
< hb
->d
.nhydro
[i
]); j
++)
2273 fprintf(fp
, " %4d %4d", hb
->d
.don
[i
]+1,
2274 hb
->d
.hydro
[i
][j
]+1);
2280 fprintf(fp
, "[ acceptors_%s ]", grpnames
[grp
]);
2281 for (i
= 0; (i
< hb
->a
.nra
); i
++)
2283 if (hb
->a
.grp
[i
] == grp
)
2285 fprintf(fp
, (i
%15 && !first
) ? " " : "\n");
2286 fprintf(fp
, " %4d", hb
->a
.acc
[i
]+1);
2295 fprintf(fp
, bContact
? "[ contacts_%s-%s ]\n" :
2296 "[ hbonds_%s-%s ]\n", grpnames
[0], grpnames
[1]);
2300 fprintf(fp
, bContact
? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames
[0]);
2303 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2306 for (k
= 0; (k
< hb
->a
.nra
); k
++)
2309 for (m
= 0; (m
< hb
->d
.nhydro
[i
]); m
++)
2311 if (hb
->hbmap
[i
][k
] && ISHB(hb
->hbmap
[i
][k
]->history
[m
]))
2313 sprintf(ds
, "%s", mkatomname(atoms
, ddd
));
2314 sprintf(as
, "%s", mkatomname(atoms
, aaa
));
2317 fprintf(fp
, " %6d %6d\n", ddd
+1, aaa
+1);
2320 fprintf(fplog
, "%12s %12s\n", ds
, as
);
2325 hhh
= hb
->d
.hydro
[i
][m
];
2326 sprintf(hs
, "%s", mkatomname(atoms
, hhh
));
2327 fprintf(fp
, " %6d %6d %6d\n", ddd
+1, hhh
+1, aaa
+1);
2330 fprintf(fplog
, "%12s %12s %12s\n", ds
, hs
, as
);
2344 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
2345 * It mimics add_frames() and init_frame() to some extent. */
2346 static void sync_hbdata(t_hbdata
*p_hb
, int nframes
)
2348 if (nframes
>= p_hb
->max_frames
)
2350 p_hb
->max_frames
+= 4096;
2351 srenew(p_hb
->nhb
, p_hb
->max_frames
);
2352 srenew(p_hb
->ndist
, p_hb
->max_frames
);
2353 srenew(p_hb
->n_bound
, p_hb
->max_frames
);
2354 srenew(p_hb
->nhx
, p_hb
->max_frames
);
2357 srenew(p_hb
->danr
, p_hb
->max_frames
);
2359 std::memset(&(p_hb
->nhb
[nframes
]), 0, sizeof(int) * (p_hb
->max_frames
-nframes
));
2360 std::memset(&(p_hb
->ndist
[nframes
]), 0, sizeof(int) * (p_hb
->max_frames
-nframes
));
2361 p_hb
->nhb
[nframes
] = 0;
2362 p_hb
->ndist
[nframes
] = 0;
2365 p_hb
->nframes
= nframes
;
2367 std::memset(&(p_hb
->nhx
[nframes
]), 0, sizeof(int)*max_hx
); /* zero the helix count for this frame */
2370 int gmx_hbond(int argc
, char *argv
[])
2372 const char *desc
[] = {
2373 "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
2374 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
2375 "(zero is extended) and the distance Donor - Acceptor",
2376 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
2377 "OH and NH groups are regarded as donors, O is an acceptor always,",
2378 "N is an acceptor by default, but this can be switched using",
2379 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
2380 "to the first preceding non-hydrogen atom.[PAR]",
2382 "You need to specify two groups for analysis, which must be either",
2383 "identical or non-overlapping. All hydrogen bonds between the two",
2384 "groups are analyzed.[PAR]",
2386 "If you set [TT]-shell[tt], you will be asked for an additional index group",
2387 "which should contain exactly one atom. In this case, only hydrogen",
2388 "bonds between atoms within the shell distance from the one atom are",
2391 "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
2392 "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
2393 "If contact kinetics are analyzed by using the -contact option, then",
2394 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
2395 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
2396 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
2397 "See mentioned literature for more details and definitions."
2400 /* "It is also possible to analyse specific hydrogen bonds with",
2401 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
2402 "Donor Hydrogen Acceptor, in the following way::",
2409 "Note that the triplets need not be on separate lines.",
2410 "Each atom triplet specifies a hydrogen bond to be analyzed,",
2411 "note also that no check is made for the types of atoms.[PAR]",
2416 " * [TT]-num[tt]: number of hydrogen bonds as a function of time.",
2417 " * [TT]-ac[tt]: average over all autocorrelations of the existence",
2418 " functions (either 0 or 1) of all hydrogen bonds.",
2419 " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
2420 " * [TT]-ang[tt]: angle distribution of all hydrogen bonds.",
2421 " * [TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
2422 " where n and n+i stand for residue numbers and i ranges from 0 to 6.",
2423 " This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
2424 " with helices in proteins.",
2425 " * [TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
2426 " for selected groups, all hydrogen bonded atoms from all groups and",
2427 " all solvent atoms involved in insertion.",
2428 " * [TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
2429 " frames, this also contains information on solvent insertion",
2430 " into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
2432 " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
2433 " each timeframe. This is especially useful when using [TT]-shell[tt].",
2434 " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
2435 " compare results to Raman Spectroscopy.",
2437 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
2438 "require an amount of memory proportional to the total numbers of donors",
2439 "times the total number of acceptors in the selected group(s)."
2442 static real acut
= 30, abin
= 1, rcut
= 0.35, r2cut
= 0, rbin
= 0.005, rshell
= -1;
2443 static real maxnhb
= 0, fit_start
= 1, fit_end
= 60, temp
= 298.15;
2444 static gmx_bool bNitAcc
= TRUE
, bDA
= TRUE
, bMerge
= TRUE
;
2445 static int nDump
= 0;
2446 static int nThreads
= 0;
2448 static gmx_bool bContact
= FALSE
;
2452 { "-a", FALSE
, etREAL
, {&acut
},
2453 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
2454 { "-r", FALSE
, etREAL
, {&rcut
},
2455 "Cutoff radius (nm, X - Acceptor, see next option)" },
2456 { "-da", FALSE
, etBOOL
, {&bDA
},
2457 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
2458 { "-r2", FALSE
, etREAL
, {&r2cut
},
2459 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
2460 { "-abin", FALSE
, etREAL
, {&abin
},
2461 "Binwidth angle distribution (degrees)" },
2462 { "-rbin", FALSE
, etREAL
, {&rbin
},
2463 "Binwidth distance distribution (nm)" },
2464 { "-nitacc", FALSE
, etBOOL
, {&bNitAcc
},
2465 "Regard nitrogen atoms as acceptors" },
2466 { "-contact", FALSE
, etBOOL
, {&bContact
},
2467 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
2468 { "-shell", FALSE
, etREAL
, {&rshell
},
2469 "when > 0, only calculate hydrogen bonds within # nm shell around "
2471 { "-fitstart", FALSE
, etREAL
, {&fit_start
},
2472 "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]" },
2473 { "-fitend", FALSE
, etREAL
, {&fit_end
},
2474 "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])" },
2475 { "-temp", FALSE
, etREAL
, {&temp
},
2476 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
2477 { "-dump", FALSE
, etINT
, {&nDump
},
2478 "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
2479 { "-max_hb", FALSE
, etREAL
, {&maxnhb
},
2480 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
2481 { "-merge", FALSE
, etBOOL
, {&bMerge
},
2482 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
2484 { "-nthreads", FALSE
, etINT
, {&nThreads
},
2485 "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
2488 const char *bugs
[] = {
2489 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
2492 { efTRX
, "-f", NULL
, ffREAD
},
2493 { efTPR
, NULL
, NULL
, ffREAD
},
2494 { efNDX
, NULL
, NULL
, ffOPTRD
},
2495 /* { efNDX, "-sel", "select", ffOPTRD },*/
2496 { efXVG
, "-num", "hbnum", ffWRITE
},
2497 { efLOG
, "-g", "hbond", ffOPTWR
},
2498 { efXVG
, "-ac", "hbac", ffOPTWR
},
2499 { efXVG
, "-dist", "hbdist", ffOPTWR
},
2500 { efXVG
, "-ang", "hbang", ffOPTWR
},
2501 { efXVG
, "-hx", "hbhelix", ffOPTWR
},
2502 { efNDX
, "-hbn", "hbond", ffOPTWR
},
2503 { efXPM
, "-hbm", "hbmap", ffOPTWR
},
2504 { efXVG
, "-don", "donor", ffOPTWR
},
2505 { efXVG
, "-dan", "danum", ffOPTWR
},
2506 { efXVG
, "-life", "hblife", ffOPTWR
},
2507 { efXVG
, "-nhbdist", "nhbdist", ffOPTWR
}
2510 #define NFILE asize(fnm)
2512 char hbmap
[HB_NR
] = { ' ', 'o', '-', '*' };
2513 const char *hbdesc
[HB_NR
] = { "None", "Present", "Inserted", "Present & Inserted" };
2514 t_rgb hbrgb
[HB_NR
] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
2516 t_trxstatus
*status
;
2521 int npargs
, natoms
, nframes
= 0, shatom
;
2527 real t
, ccut
, dist
= 0.0, ang
= 0.0;
2528 double max_nhb
, aver_nhb
, aver_dist
;
2529 int h
= 0, i
= 0, j
, k
= 0, ogrp
, nsel
;
2531 int xj
, yj
, zj
, aj
, xjj
, yjj
, zjj
;
2532 gmx_bool bSelected
, bHBmap
, bStop
, bTwo
, bBox
, bTric
;
2534 int grp
, nabin
, nrbin
, resdist
, ihb
;
2537 FILE *fp
, *fpnhb
= NULL
, *donor_properties
= NULL
;
2539 t_ncell
*icell
, *jcell
;
2541 unsigned char *datable
;
2542 gmx_output_env_t
*oenv
;
2543 int ii
, hh
, actual_nThreads
;
2546 gmx_bool bEdge_yjj
, bEdge_xjj
, bOMP
;
2548 t_hbdata
**p_hb
= NULL
; /* one per thread, then merge after the frame loop */
2549 int **p_adist
= NULL
, **p_rdist
= NULL
; /* a histogram for each thread. */
2558 ppa
= add_acf_pargs(&npargs
, pa
);
2560 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_TIME_UNIT
, NFILE
, fnm
, npargs
,
2561 ppa
, asize(desc
), desc
, asize(bugs
), bugs
, &oenv
))
2568 ccut
= std::cos(acut
*DEG2RAD
);
2574 gmx_fatal(FARGS
, "Can not analyze selected contacts.");
2578 gmx_fatal(FARGS
, "Can not analyze contact between H and A: turn off -noda");
2582 /* Initiate main data structure! */
2583 bHBmap
= (opt2bSet("-ac", NFILE
, fnm
) ||
2584 opt2bSet("-life", NFILE
, fnm
) ||
2585 opt2bSet("-hbn", NFILE
, fnm
) ||
2586 opt2bSet("-hbm", NFILE
, fnm
));
2588 if (opt2bSet("-nhbdist", NFILE
, fnm
))
2590 const char *leg
[MAXHH
+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
2591 fpnhb
= xvgropen(opt2fn("-nhbdist", NFILE
, fnm
),
2592 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv
), "N", oenv
);
2593 xvgr_legend(fpnhb
, asize(leg
), leg
, oenv
);
2596 hb
= mk_hbdata(bHBmap
, opt2bSet("-dan", NFILE
, fnm
), bMerge
|| bContact
);
2599 read_tpx_top(ftp2fn(efTPR
, NFILE
, fnm
), &ir
, box
, &natoms
, NULL
, NULL
, &top
);
2601 snew(grpnames
, grNR
);
2604 /* Make Donor-Acceptor table */
2605 snew(datable
, top
.atoms
.nr
);
2606 gen_datable(index
[0], isize
[0], datable
, top
.atoms
.nr
);
2610 /* analyze selected hydrogen bonds */
2611 printf("Select group with selected atoms:\n");
2612 get_index(&(top
.atoms
), opt2fn("-sel", NFILE
, fnm
),
2613 1, &nsel
, index
, grpnames
);
2616 gmx_fatal(FARGS
, "Number of atoms in group '%s' not a multiple of 3\n"
2617 "and therefore cannot contain triplets of "
2618 "Donor-Hydrogen-Acceptor", grpnames
[0]);
2622 for (i
= 0; (i
< nsel
); i
+= 3)
2624 int dd
= index
[0][i
];
2625 int aa
= index
[0][i
+2];
2626 /* int */ hh
= index
[0][i
+1];
2627 add_dh (&hb
->d
, dd
, hh
, i
, datable
);
2628 add_acc(&hb
->a
, aa
, i
);
2629 /* Should this be here ? */
2630 snew(hb
->d
.dptr
, top
.atoms
.nr
);
2631 snew(hb
->a
.aptr
, top
.atoms
.nr
);
2632 add_hbond(hb
, dd
, aa
, hh
, gr0
, gr0
, 0, bMerge
, 0, bContact
);
2634 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
2635 isize
[0], grpnames
[0]);
2639 /* analyze all hydrogen bonds: get group(s) */
2640 printf("Specify 2 groups to analyze:\n");
2641 get_index(&(top
.atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
),
2642 2, isize
, index
, grpnames
);
2644 /* check if we have two identical or two non-overlapping groups */
2645 bTwo
= isize
[0] != isize
[1];
2646 for (i
= 0; (i
< isize
[0]) && !bTwo
; i
++)
2648 bTwo
= index
[0][i
] != index
[1][i
];
2652 printf("Checking for overlap in atoms between %s and %s\n",
2653 grpnames
[0], grpnames
[1]);
2654 for (i
= 0; i
< isize
[1]; i
++)
2656 if (ISINGRP(datable
[index
[1][i
]]))
2658 gmx_fatal(FARGS
, "Partial overlap between groups '%s' and '%s'",
2659 grpnames
[0], grpnames
[1]);
2663 printf("Checking for overlap in atoms between %s and %s\n",
2664 grpnames[0],grpnames[1]);
2665 for (i=0; i<isize[0]; i++)
2666 for (j=0; j<isize[1]; j++)
2667 if (index[0][i] == index[1][j])
2668 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
2669 grpnames[0],grpnames[1]);
2674 printf("Calculating %s "
2675 "between %s (%d atoms) and %s (%d atoms)\n",
2676 bContact
? "contacts" : "hydrogen bonds",
2677 grpnames
[0], isize
[0], grpnames
[1], isize
[1]);
2681 fprintf(stderr
, "Calculating %s in %s (%d atoms)\n",
2682 bContact
? "contacts" : "hydrogen bonds", grpnames
[0], isize
[0]);
2687 /* search donors and acceptors in groups */
2688 snew(datable
, top
.atoms
.nr
);
2689 for (i
= 0; (i
< grNR
); i
++)
2691 if ( ((i
== gr0
) && !bSelected
) ||
2692 ((i
== gr1
) && bTwo
))
2694 gen_datable(index
[i
], isize
[i
], datable
, top
.atoms
.nr
);
2697 search_acceptors(&top
, isize
[i
], index
[i
], &hb
->a
, i
,
2698 bNitAcc
, TRUE
, (bTwo
&& (i
== gr0
)) || !bTwo
, datable
);
2699 search_donors (&top
, isize
[i
], index
[i
], &hb
->d
, i
,
2700 TRUE
, (bTwo
&& (i
== gr1
)) || !bTwo
, datable
);
2704 search_acceptors(&top
, isize
[i
], index
[i
], &hb
->a
, i
, bNitAcc
, FALSE
, TRUE
, datable
);
2705 search_donors (&top
, isize
[i
], index
[i
], &hb
->d
, i
, FALSE
, TRUE
, datable
);
2709 clear_datable_grp(datable
, top
.atoms
.nr
);
2714 printf("Found %d donors and %d acceptors\n", hb
->d
.nrd
, hb
->a
.nra
);
2716 snew(donors[gr0D], dons[gr0D].nrd);*/
2718 donor_properties
= open_donor_properties_file(opt2fn_null("-don", NFILE
, fnm
), hb
, oenv
);
2722 printf("Making hbmap structure...");
2723 /* Generate hbond data structure */
2730 if (hb
->d
.nrd
+ hb
->a
.nra
== 0)
2732 printf("No Donors or Acceptors found\n");
2739 printf("No Donors found\n");
2744 printf("No Acceptors found\n");
2750 gmx_fatal(FARGS
, "Nothing to be done");
2759 /* get index group with atom for shell */
2762 printf("Select atom for shell (1 atom):\n");
2763 get_index(&(top
.atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
),
2764 1, &shisz
, &shidx
, &shgrpnm
);
2767 printf("group contains %d atoms, should be 1 (one)\n", shisz
);
2772 printf("Will calculate hydrogen bonds within a shell "
2773 "of %g nm around atom %i\n", rshell
, shatom
+1);
2776 /* Analyze trajectory */
2777 natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
2778 if (natoms
> top
.atoms
.nr
)
2780 gmx_fatal(FARGS
, "Topology (%d atoms) does not match trajectory (%d atoms)",
2781 top
.atoms
.nr
, natoms
);
2784 bBox
= ir
.ePBC
!= epbcNONE
;
2785 grid
= init_grid(bBox
, box
, (rcut
> r2cut
) ? rcut
: r2cut
, ngrid
);
2786 nabin
= static_cast<int>(acut
/abin
);
2787 nrbin
= static_cast<int>(rcut
/rbin
);
2788 snew(adist
, nabin
+1);
2789 snew(rdist
, nrbin
+1);
2792 #define __ADIST adist
2793 #define __RDIST rdist
2795 #else /* GMX_OPENMP ================================================== \
2796 * Set up the OpenMP stuff, |
2797 * like the number of threads and such |
2798 * Also start the parallel loop. |
2800 #define __ADIST p_adist[threadNr]
2801 #define __RDIST p_rdist[threadNr]
2802 #define __HBDATA p_hb[threadNr]
2806 bParallel
= !bSelected
;
2810 actual_nThreads
= std::min((nThreads
<= 0) ? INT_MAX
: nThreads
, gmx_omp_get_max_threads());
2812 gmx_omp_set_num_threads(actual_nThreads
);
2813 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads
);
2818 actual_nThreads
= 1;
2821 snew(p_hb
, actual_nThreads
);
2822 snew(p_adist
, actual_nThreads
);
2823 snew(p_rdist
, actual_nThreads
);
2824 for (i
= 0; i
< actual_nThreads
; i
++)
2827 snew(p_adist
[i
], nabin
+1);
2828 snew(p_rdist
[i
], nrbin
+1);
2830 p_hb
[i
]->max_frames
= 0;
2831 p_hb
[i
]->nhb
= NULL
;
2832 p_hb
[i
]->ndist
= NULL
;
2833 p_hb
[i
]->n_bound
= NULL
;
2834 p_hb
[i
]->time
= NULL
;
2835 p_hb
[i
]->nhx
= NULL
;
2837 p_hb
[i
]->bHBmap
= hb
->bHBmap
;
2838 p_hb
[i
]->bDAnr
= hb
->bDAnr
;
2839 p_hb
[i
]->wordlen
= hb
->wordlen
;
2840 p_hb
[i
]->nframes
= hb
->nframes
;
2841 p_hb
[i
]->maxhydro
= hb
->maxhydro
;
2842 p_hb
[i
]->danr
= hb
->danr
;
2845 p_hb
[i
]->hbmap
= hb
->hbmap
;
2846 p_hb
[i
]->time
= hb
->time
; /* This may need re-syncing at every frame. */
2849 p_hb
[i
]->nrdist
= 0;
2853 /* Make a thread pool here,
2854 * instead of forking anew at every frame. */
2856 #pragma omp parallel \
2858 private(j, h, ii, hh, \
2859 xi, yi, zi, xj, yj, zj, threadNr, \
2860 dist, ang, icell, jcell, \
2861 grp, ogrp, ai, aj, xjj, yjj, zjj, \
2864 bEdge_xjj, bEdge_yjj) \
2866 { /* Start of parallel region */
2867 threadNr
= gmx_omp_get_thread_num();
2872 bTric
= bBox
&& TRICLINIC(box
);
2878 sync_hbdata(p_hb
[threadNr
], nframes
);
2880 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2886 build_grid(hb
, x
, x
[shatom
], bBox
, box
, hbox
, (rcut
> r2cut
) ? rcut
: r2cut
,
2887 rshell
, ngrid
, grid
);
2888 reset_nhbonds(&(hb
->d
));
2890 if (debug
&& bDebug
)
2892 dump_grid(debug
, ngrid
, grid
);
2895 add_frames(hb
, nframes
);
2896 init_hbframe(hb
, nframes
, output_env_conv_time(oenv
, t
));
2900 count_da_grid(ngrid
, grid
, hb
->danr
[nframes
]);
2903 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2908 p_hb
[threadNr
]->time
= hb
->time
; /* This pointer may have changed. */
2918 /* Do not parallelize this just yet. */
2920 for (ii
= 0; (ii
< nsel
); ii
++)
2922 int dd
= index
[0][i
];
2923 int aa
= index
[0][i
+2];
2924 /* int */ hh
= index
[0][i
+1];
2925 ihb
= is_hbond(hb
, ii
, ii
, dd
, aa
, rcut
, r2cut
, ccut
, x
, bBox
, box
,
2926 hbox
, &dist
, &ang
, bDA
, &h
, bContact
, bMerge
);
2930 /* add to index if not already there */
2932 add_hbond(hb
, dd
, aa
, hh
, ii
, ii
, nframes
, bMerge
, ihb
, bContact
);
2936 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2938 } /* if (bSelected) */
2941 /* The outer grid loop will have to do for now. */
2942 #pragma omp for schedule(dynamic)
2943 for (xi
= 0; xi
< ngrid
[XX
]; xi
++)
2947 for (yi
= 0; (yi
< ngrid
[YY
]); yi
++)
2949 for (zi
= 0; (zi
< ngrid
[ZZ
]); zi
++)
2952 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
2953 for (grp
= gr0
; (grp
<= (bTwo
? gr1
: gr0
)); grp
++)
2955 icell
= &(grid
[zi
][yi
][xi
].d
[grp
]);
2966 /* loop over all hydrogen atoms from group (grp)
2967 * in this gridcell (icell)
2969 for (ai
= 0; (ai
< icell
->nr
); ai
++)
2971 i
= icell
->atoms
[ai
];
2973 /* loop over all adjacent gridcells (xj,yj,zj) */
2974 for (zjj
= grid_loop_begin(ngrid
[ZZ
], zi
, bTric
, FALSE
);
2975 zjj
<= grid_loop_end(ngrid
[ZZ
], zi
, bTric
, FALSE
);
2978 zj
= grid_mod(zjj
, ngrid
[ZZ
]);
2979 bEdge_yjj
= (zj
== 0) || (zj
== ngrid
[ZZ
] - 1);
2980 for (yjj
= grid_loop_begin(ngrid
[YY
], yi
, bTric
, bEdge_yjj
);
2981 yjj
<= grid_loop_end(ngrid
[YY
], yi
, bTric
, bEdge_yjj
);
2984 yj
= grid_mod(yjj
, ngrid
[YY
]);
2986 (yj
== 0) || (yj
== ngrid
[YY
] - 1) ||
2987 (zj
== 0) || (zj
== ngrid
[ZZ
] - 1);
2988 for (xjj
= grid_loop_begin(ngrid
[XX
], xi
, bTric
, bEdge_xjj
);
2989 xjj
<= grid_loop_end(ngrid
[XX
], xi
, bTric
, bEdge_xjj
);
2992 xj
= grid_mod(xjj
, ngrid
[XX
]);
2993 jcell
= &(grid
[zj
][yj
][xj
].a
[ogrp
]);
2994 /* loop over acceptor atoms from other group (ogrp)
2995 * in this adjacent gridcell (jcell)
2997 for (aj
= 0; (aj
< jcell
->nr
); aj
++)
2999 j
= jcell
->atoms
[aj
];
3001 /* check if this once was a h-bond */
3002 ihb
= is_hbond(__HBDATA
, grp
, ogrp
, i
, j
, rcut
, r2cut
, ccut
, x
, bBox
, box
,
3003 hbox
, &dist
, &ang
, bDA
, &h
, bContact
, bMerge
);
3007 /* add to index if not already there */
3009 add_hbond(__HBDATA
, i
, j
, h
, grp
, ogrp
, nframes
, bMerge
, ihb
, bContact
);
3011 /* make angle and distance distributions */
3012 if (ihb
== hbHB
&& !bContact
)
3016 gmx_fatal(FARGS
, "distance is higher than what is allowed for an hbond: %f", dist
);
3019 __ADIST
[static_cast<int>( ang
/abin
)]++;
3020 __RDIST
[static_cast<int>(dist
/rbin
)]++;
3023 if (donor_index(&hb
->d
, grp
, i
) == NOTSET
)
3025 gmx_fatal(FARGS
, "Invalid donor %d", i
);
3027 if (acceptor_index(&hb
->a
, ogrp
, j
) == NOTSET
)
3029 gmx_fatal(FARGS
, "Invalid acceptor %d", j
);
3031 resdist
= std::abs(top
.atoms
.atom
[i
].resind
-top
.atoms
.atom
[j
].resind
);
3032 if (resdist
>= max_hx
)
3036 __HBDATA
->nhx
[nframes
][resdist
]++;
3047 } /* for xi,yi,zi */
3050 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3052 } /* if (bSelected) {...} else */
3055 /* Better wait for all threads to finnish using x[] before updating it. */
3058 #pragma omp critical
3062 /* Sum up histograms and counts from p_hb[] into hb */
3065 hb
->nhb
[k
] += p_hb
[threadNr
]->nhb
[k
];
3066 hb
->ndist
[k
] += p_hb
[threadNr
]->ndist
[k
];
3067 for (j
= 0; j
< max_hx
; j
++)
3069 hb
->nhx
[k
][j
] += p_hb
[threadNr
]->nhx
[k
][j
];
3073 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3076 /* Here are a handful of single constructs
3077 * to share the workload a bit. The most
3078 * important one is of course the last one,
3079 * where there's a potential bottleneck in form
3086 analyse_donor_properties(donor_properties
, hb
, k
, t
);
3088 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3097 do_nhb_dist(fpnhb
, hb
, t
);
3100 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3107 trrStatus
= (read_next_x(oenv
, status
, &t
, x
, box
));
3110 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3119 #pragma omp critical
3121 hb
->nrhb
+= p_hb
[threadNr
]->nrhb
;
3122 hb
->nrdist
+= p_hb
[threadNr
]->nrdist
;
3125 /* Free parallel datastructures */
3126 sfree(p_hb
[threadNr
]->nhb
);
3127 sfree(p_hb
[threadNr
]->ndist
);
3128 sfree(p_hb
[threadNr
]->nhx
);
3131 for (i
= 0; i
< nabin
; i
++)
3135 for (j
= 0; j
< actual_nThreads
; j
++)
3138 adist
[i
] += p_adist
[j
][i
];
3141 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3144 for (i
= 0; i
<= nrbin
; i
++)
3148 for (j
= 0; j
< actual_nThreads
; j
++)
3150 rdist
[i
] += p_rdist
[j
][i
];
3153 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3156 sfree(p_adist
[threadNr
]);
3157 sfree(p_rdist
[threadNr
]);
3159 } /* End of parallel region */
3166 if (nframes
< 2 && (opt2bSet("-ac", NFILE
, fnm
) || opt2bSet("-life", NFILE
, fnm
)))
3168 gmx_fatal(FARGS
, "Cannot calculate autocorrelation of life times with less than two frames");
3171 free_grid(ngrid
, &grid
);
3175 if (donor_properties
)
3177 xvgrclose(donor_properties
);
3185 /* Compute maximum possible number of different hbonds */
3192 max_nhb
= 0.5*(hb
->d
.nrd
*hb
->a
.nra
);
3199 printf("No %s found!!\n", bContact
? "contacts" : "hydrogen bonds");
3203 printf("Found %d different %s in trajectory\n"
3204 "Found %d different atom-pairs within %s distance\n",
3205 hb
->nrhb
, bContact
? "contacts" : "hydrogen bonds",
3206 hb
->nrdist
, (r2cut
> 0) ? "second cut-off" : "hydrogen bonding");
3210 merge_hb(hb
, bTwo
, bContact
);
3213 if (opt2bSet("-hbn", NFILE
, fnm
))
3215 dump_hbmap(hb
, NFILE
, fnm
, bTwo
, bContact
, isize
, index
, grpnames
, &top
.atoms
);
3218 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3219 * to make the -hbn and -hmb output match eachother.
3220 * - Erik Marklund, May 30, 2006 */
3223 /* Print out number of hbonds and distances */
3226 fp
= xvgropen(opt2fn("-num", NFILE
, fnm
), bContact
? "Contacts" :
3227 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv
), "Number", oenv
);
3229 snew(leg
[0], STRLEN
);
3230 snew(leg
[1], STRLEN
);
3231 sprintf(leg
[0], "%s", bContact
? "Contacts" : "Hydrogen bonds");
3232 sprintf(leg
[1], "Pairs within %g nm", (r2cut
> 0) ? r2cut
: rcut
);
3233 xvgr_legend(fp
, 2, (const char**)leg
, oenv
);
3237 for (i
= 0; (i
< nframes
); i
++)
3239 fprintf(fp
, "%10g %10d %10d\n", hb
->time
[i
], hb
->nhb
[i
], hb
->ndist
[i
]);
3240 aver_nhb
+= hb
->nhb
[i
];
3241 aver_dist
+= hb
->ndist
[i
];
3244 aver_nhb
/= nframes
;
3245 aver_dist
/= nframes
;
3246 /* Print HB distance distribution */
3247 if (opt2bSet("-dist", NFILE
, fnm
))
3252 for (i
= 0; i
< nrbin
; i
++)
3257 fp
= xvgropen(opt2fn("-dist", NFILE
, fnm
),
3258 "Hydrogen Bond Distribution",
3260 "Donor - Acceptor Distance (nm)" :
3261 "Hydrogen - Acceptor Distance (nm)", "", oenv
);
3262 for (i
= 0; i
< nrbin
; i
++)
3264 fprintf(fp
, "%10g %10g\n", (i
+0.5)*rbin
, rdist
[i
]/(rbin
*sum
));
3269 /* Print HB angle distribution */
3270 if (opt2bSet("-ang", NFILE
, fnm
))
3275 for (i
= 0; i
< nabin
; i
++)
3280 fp
= xvgropen(opt2fn("-ang", NFILE
, fnm
),
3281 "Hydrogen Bond Distribution",
3282 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv
);
3283 for (i
= 0; i
< nabin
; i
++)
3285 fprintf(fp
, "%10g %10g\n", (i
+0.5)*abin
, adist
[i
]/(abin
*sum
));
3290 /* Print HB in alpha-helix */
3291 if (opt2bSet("-hx", NFILE
, fnm
))
3293 fp
= xvgropen(opt2fn("-hx", NFILE
, fnm
),
3294 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv
), "Count", oenv
);
3295 xvgr_legend(fp
, NRHXTYPES
, hxtypenames
, oenv
);
3296 for (i
= 0; i
< nframes
; i
++)
3298 fprintf(fp
, "%10g", hb
->time
[i
]);
3299 for (j
= 0; j
< max_hx
; j
++)
3301 fprintf(fp
, " %6d", hb
->nhx
[i
][j
]);
3308 printf("Average number of %s per timeframe %.3f out of %g possible\n",
3309 bContact
? "contacts" : "hbonds",
3310 bContact
? aver_dist
: aver_nhb
, max_nhb
);
3312 /* Do Autocorrelation etc. */
3316 Added support for -contact in ac and hbm calculations below.
3317 - Erik Marklund, May 29, 2006
3319 if (opt2bSet("-ac", NFILE
, fnm
) || opt2bSet("-life", NFILE
, fnm
))
3321 please_cite(stdout
, "Spoel2006b");
3323 if (opt2bSet("-ac", NFILE
, fnm
))
3325 do_hbac(opt2fn("-ac", NFILE
, fnm
), hb
, nDump
,
3326 bMerge
, bContact
, fit_start
, temp
, r2cut
> 0, oenv
,
3329 if (opt2bSet("-life", NFILE
, fnm
))
3331 do_hblife(opt2fn("-life", NFILE
, fnm
), hb
, bMerge
, bContact
, oenv
);
3333 if (opt2bSet("-hbm", NFILE
, fnm
))
3336 int id
, ia
, hh
, x
, y
;
3337 mat
.flags
= mat
.y0
= 0;
3339 if ((nframes
> 0) && (hb
->nrhb
> 0))
3344 snew(mat
.matrix
, mat
.nx
);
3345 for (x
= 0; (x
< mat
.nx
); x
++)
3347 snew(mat
.matrix
[x
], mat
.ny
);
3350 for (id
= 0; (id
< hb
->d
.nrd
); id
++)
3352 for (ia
= 0; (ia
< hb
->a
.nra
); ia
++)
3354 for (hh
= 0; (hh
< hb
->maxhydro
); hh
++)
3356 if (hb
->hbmap
[id
][ia
])
3358 if (ISHB(hb
->hbmap
[id
][ia
]->history
[hh
]))
3360 for (x
= 0; (x
<= hb
->hbmap
[id
][ia
]->nframes
); x
++)
3362 int nn0
= hb
->hbmap
[id
][ia
]->n0
;
3363 range_check(y
, 0, mat
.ny
);
3364 mat
.matrix
[x
+nn0
][y
] = is_hb(hb
->hbmap
[id
][ia
]->h
[hh
], x
);
3372 mat
.axis_x
= hb
->time
;
3373 snew(mat
.axis_y
, mat
.ny
);
3374 for (j
= 0; j
< mat
.ny
; j
++)
3378 sprintf(mat
.title
, bContact
? "Contact Existence Map" :
3379 "Hydrogen Bond Existence Map");
3380 sprintf(mat
.legend
, bContact
? "Contacts" : "Hydrogen Bonds");
3381 sprintf(mat
.label_x
, "%s", output_env_get_xvgr_tlabel(oenv
));
3382 sprintf(mat
.label_y
, bContact
? "Contact Index" : "Hydrogen Bond Index");
3383 mat
.bDiscrete
= TRUE
;
3385 snew(mat
.map
, mat
.nmap
);
3386 for (i
= 0; i
< mat
.nmap
; i
++)
3388 mat
.map
[i
].code
.c1
= hbmap
[i
];
3389 mat
.map
[i
].desc
= hbdesc
[i
];
3390 mat
.map
[i
].rgb
= hbrgb
[i
];
3392 fp
= opt2FILE("-hbm", NFILE
, fnm
, "w");
3393 write_xpm_m(fp
, mat
);
3395 for (x
= 0; x
< mat
.nx
; x
++)
3397 sfree(mat
.matrix
[x
]);
3405 fprintf(stderr
, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
3416 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
3418 fp
= xvgropen(opt2fn("-dan", NFILE
, fnm
),
3419 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv
), "Count", oenv
);
3420 nleg
= (bTwo
? 2 : 1)*2;
3421 snew(legnames
, nleg
);
3423 for (j
= 0; j
< grNR
; j
++)
3425 if (USE_THIS_GROUP(j
) )
3427 sprintf(buf
, "Donors %s", grpnames
[j
]);
3428 legnames
[i
++] = gmx_strdup(buf
);
3429 sprintf(buf
, "Acceptors %s", grpnames
[j
]);
3430 legnames
[i
++] = gmx_strdup(buf
);
3435 gmx_incons("number of legend entries");
3437 xvgr_legend(fp
, nleg
, (const char**)legnames
, oenv
);
3438 for (i
= 0; i
< nframes
; i
++)
3440 fprintf(fp
, "%10g", hb
->time
[i
]);
3441 for (j
= 0; (j
< grNR
); j
++)
3443 if (USE_THIS_GROUP(j
) )
3445 fprintf(fp
, " %6d", hb
->danr
[i
][j
]);