2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/commandline/viewit.h"
49 #include "gromacs/correlationfunctions/autocorr.h"
50 #include "gromacs/correlationfunctions/crosscorr.h"
51 #include "gromacs/correlationfunctions/expfit.h"
52 #include "gromacs/correlationfunctions/integrate.h"
53 #include "gromacs/fileio/matio.h"
54 #include "gromacs/fileio/tpxio.h"
55 #include "gromacs/fileio/trxio.h"
56 #include "gromacs/fileio/xvgr.h"
57 #include "gromacs/gmxana/gmx_ana.h"
58 #include "gromacs/gmxana/gstat.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/topology/ifunc.h"
65 #include "gromacs/topology/index.h"
66 #include "gromacs/topology/topology.h"
67 #include "gromacs/utility/arraysize.h"
68 #include "gromacs/utility/cstringutil.h"
69 #include "gromacs/utility/exceptions.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/futil.h"
72 #include "gromacs/utility/gmxomp.h"
73 #include "gromacs/utility/pleasecite.h"
74 #include "gromacs/utility/programcontext.h"
75 #include "gromacs/utility/smalloc.h"
76 #include "gromacs/utility/snprintf.h"
79 typedef int t_hx
[max_hx
];
80 #define NRHXTYPES max_hx
81 const char *hxtypenames
[NRHXTYPES
] =
82 {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
85 static const int NOTSET
= -49297;
87 /* -----------------------------------------*/
93 hbNo
, hbDist
, hbHB
, hbNR
, hbR2
96 noDA
, ACC
, DON
, DA
, INGROUP
99 static const char *grpnames
[grNR
] = {"0", "1", "I" };
101 static gmx_bool bDebug
= FALSE
;
106 #define HB_YESINS HB_YES|HB_INS
110 #define ISHB(h) (((h) & 2) == 2)
111 #define ISDIST(h) (((h) & 1) == 1)
112 #define ISDIST2(h) (((h) & 4) == 4)
113 #define ISACC(h) (((h) & 1) == 1)
114 #define ISDON(h) (((h) & 2) == 2)
115 #define ISINGRP(h) (((h) & 4) == 4)
128 typedef int t_icell
[grNR
];
129 typedef int h_id
[MAXHYDRO
];
132 int history
[MAXHYDRO
];
133 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
134 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
136 /* Bitmask array which tells whether a hbond is present
137 * at a given time. Either of these may be NULL
139 int n0
; /* First frame a HB was found */
140 int nframes
, maxframes
; /* Amount of frames in this hbond */
143 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
144 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
145 * acceptor distance is less than the user-specified distance (typically
152 int *acc
; /* Atom numbers of the acceptors */
153 int *grp
; /* Group index */
154 int *aptr
; /* Map atom number to acceptor index */
159 int *don
; /* Atom numbers of the donors */
160 int *grp
; /* Group index */
161 int *dptr
; /* Map atom number to donor index */
162 int *nhydro
; /* Number of hydrogens for each donor */
163 h_id
*hydro
; /* The atom numbers of the hydrogens */
164 h_id
*nhbonds
; /* The number of HBs per H at current */
168 gmx_bool bHBmap
, bDAnr
;
170 /* The following arrays are nframes long */
171 int nframes
, max_frames
, maxhydro
;
177 /* These structures are initialized from the topology at start up */
180 /* This holds a matrix with all possible hydrogen bonds */
185 /* Changed argument 'bMerge' into 'oneHB' below,
186 * since -contact should cause maxhydro to be 1,
188 * - Erik Marklund May 29, 2006
191 static t_hbdata
*mk_hbdata(gmx_bool bHBmap
, gmx_bool bDAnr
, gmx_bool oneHB
)
196 hb
->wordlen
= 8*sizeof(unsigned int);
205 hb
->maxhydro
= MAXHYDRO
;
210 static void mk_hbmap(t_hbdata
*hb
)
214 snew(hb
->hbmap
, hb
->d
.nrd
);
215 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
217 snew(hb
->hbmap
[i
], hb
->a
.nra
);
218 if (hb
->hbmap
[i
] == NULL
)
220 gmx_fatal(FARGS
, "Could not allocate enough memory for hbmap");
222 for (j
= 0; (j
> hb
->a
.nra
); j
++)
224 hb
->hbmap
[i
][j
] = NULL
;
229 static void add_frames(t_hbdata
*hb
, int nframes
)
231 if (nframes
>= hb
->max_frames
)
233 hb
->max_frames
+= 4096;
234 srenew(hb
->time
, hb
->max_frames
);
235 srenew(hb
->nhb
, hb
->max_frames
);
236 srenew(hb
->ndist
, hb
->max_frames
);
237 srenew(hb
->n_bound
, hb
->max_frames
);
238 srenew(hb
->nhx
, hb
->max_frames
);
241 srenew(hb
->danr
, hb
->max_frames
);
244 hb
->nframes
= nframes
;
247 #define OFFSET(frame) (frame / 32)
248 #define MASK(frame) (1 << (frame % 32))
250 static void _set_hb(unsigned int hbexist
[], unsigned int frame
, gmx_bool bValue
)
254 hbexist
[OFFSET(frame
)] |= MASK(frame
);
258 hbexist
[OFFSET(frame
)] &= ~MASK(frame
);
262 static gmx_bool
is_hb(unsigned int hbexist
[], int frame
)
264 return ((hbexist
[OFFSET(frame
)] & MASK(frame
)) != 0) ? 1 : 0;
267 static void set_hb(t_hbdata
*hb
, int id
, int ih
, int ia
, int frame
, int ihb
)
269 unsigned int *ghptr
= NULL
;
273 ghptr
= hb
->hbmap
[id
][ia
]->h
[ih
];
275 else if (ihb
== hbDist
)
277 ghptr
= hb
->hbmap
[id
][ia
]->g
[ih
];
281 gmx_fatal(FARGS
, "Incomprehensible iValue %d in set_hb", ihb
);
284 _set_hb(ghptr
, frame
-hb
->hbmap
[id
][ia
]->n0
, TRUE
);
287 static void add_ff(t_hbdata
*hbd
, int id
, int h
, int ia
, int frame
, int ihb
)
290 t_hbond
*hb
= hbd
->hbmap
[id
][ia
];
291 int maxhydro
= std::min(hbd
->maxhydro
, hbd
->d
.nhydro
[id
]);
292 int wlen
= hbd
->wordlen
;
298 hb
->maxframes
= delta
;
299 for (i
= 0; (i
< maxhydro
); i
++)
301 snew(hb
->h
[i
], hb
->maxframes
/wlen
);
302 snew(hb
->g
[i
], hb
->maxframes
/wlen
);
307 hb
->nframes
= frame
-hb
->n0
;
308 /* We need a while loop here because hbonds may be returning
311 while (hb
->nframes
>= hb
->maxframes
)
313 n
= hb
->maxframes
+ delta
;
314 for (i
= 0; (i
< maxhydro
); i
++)
316 srenew(hb
->h
[i
], n
/wlen
);
317 srenew(hb
->g
[i
], n
/wlen
);
318 for (j
= hb
->maxframes
/wlen
; (j
< n
/wlen
); j
++)
330 set_hb(hbd
, id
, h
, ia
, frame
, ihb
);
335 static void inc_nhbonds(t_donors
*ddd
, int d
, int h
)
338 int dptr
= ddd
->dptr
[d
];
340 for (j
= 0; (j
< ddd
->nhydro
[dptr
]); j
++)
342 if (ddd
->hydro
[dptr
][j
] == h
)
344 ddd
->nhbonds
[dptr
][j
]++;
348 if (j
== ddd
->nhydro
[dptr
])
350 gmx_fatal(FARGS
, "No such hydrogen %d on donor %d\n", h
+1, d
+1);
354 static int _acceptor_index(t_acceptors
*a
, int grp
, int i
,
355 const char *file
, int line
)
359 if (a
->grp
[ai
] != grp
)
363 fprintf(debug
, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
364 ai
, a
->grp
[ai
], grp
, file
, line
);
373 #define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
375 static int _donor_index(t_donors
*d
, int grp
, int i
, const char *file
, int line
)
384 if (d
->grp
[di
] != grp
)
388 fprintf(debug
, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
389 di
, d
->grp
[di
], grp
, file
, line
);
398 #define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
400 static gmx_bool
isInterchangable(t_hbdata
*hb
, int d
, int a
, int grpa
, int grpd
)
402 /* g_hbond doesn't allow overlapping groups */
408 donor_index(&hb
->d
, grpd
, a
) != NOTSET
409 && acceptor_index(&hb
->a
, grpa
, d
) != NOTSET
;
413 static void add_hbond(t_hbdata
*hb
, int d
, int a
, int h
, int grpd
, int grpa
,
414 int frame
, gmx_bool bMerge
, int ihb
, gmx_bool bContact
)
417 gmx_bool daSwap
= FALSE
;
419 if ((id
= hb
->d
.dptr
[d
]) == NOTSET
)
421 gmx_fatal(FARGS
, "No donor atom %d", d
+1);
423 else if (grpd
!= hb
->d
.grp
[id
])
425 gmx_fatal(FARGS
, "Inconsistent donor groups, %d iso %d, atom %d",
426 grpd
, hb
->d
.grp
[id
], d
+1);
428 if ((ia
= hb
->a
.aptr
[a
]) == NOTSET
)
430 gmx_fatal(FARGS
, "No acceptor atom %d", a
+1);
432 else if (grpa
!= hb
->a
.grp
[ia
])
434 gmx_fatal(FARGS
, "Inconsistent acceptor groups, %d iso %d, atom %d",
435 grpa
, hb
->a
.grp
[ia
], a
+1);
441 if (isInterchangable(hb
, d
, a
, grpd
, grpa
) && d
> a
)
442 /* Then swap identity so that the id of d is lower then that of a.
444 * This should really be redundant by now, as is_hbond() now ought to return
445 * hbNo in the cases where this conditional is TRUE. */
452 /* Now repeat donor/acc check. */
453 if ((id
= hb
->d
.dptr
[d
]) == NOTSET
)
455 gmx_fatal(FARGS
, "No donor atom %d", d
+1);
457 else if (grpd
!= hb
->d
.grp
[id
])
459 gmx_fatal(FARGS
, "Inconsistent donor groups, %d iso %d, atom %d",
460 grpd
, hb
->d
.grp
[id
], d
+1);
462 if ((ia
= hb
->a
.aptr
[a
]) == NOTSET
)
464 gmx_fatal(FARGS
, "No acceptor atom %d", a
+1);
466 else if (grpa
!= hb
->a
.grp
[ia
])
468 gmx_fatal(FARGS
, "Inconsistent acceptor groups, %d iso %d, atom %d",
469 grpa
, hb
->a
.grp
[ia
], a
+1);
476 /* Loop over hydrogens to find which hydrogen is in this particular HB */
477 if ((ihb
== hbHB
) && !bMerge
&& !bContact
)
479 for (k
= 0; (k
< hb
->d
.nhydro
[id
]); k
++)
481 if (hb
->d
.hydro
[id
][k
] == h
)
486 if (k
== hb
->d
.nhydro
[id
])
488 gmx_fatal(FARGS
, "Donor %d does not have hydrogen %d (a = %d)",
504 if (hb
->hbmap
[id
][ia
] == NULL
)
506 snew(hb
->hbmap
[id
][ia
], 1);
507 snew(hb
->hbmap
[id
][ia
]->h
, hb
->maxhydro
);
508 snew(hb
->hbmap
[id
][ia
]->g
, hb
->maxhydro
);
510 add_ff(hb
, id
, k
, ia
, frame
, ihb
);
512 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
516 /* Strange construction with frame >=0 is a relic from old code
517 * for selected hbond analysis. It may be necessary again if that
518 * is made to work again.
522 hh
= hb
->hbmap
[id
][ia
]->history
[k
];
528 hb
->hbmap
[id
][ia
]->history
[k
] = hh
| 2;
539 hb
->hbmap
[id
][ia
]->history
[k
] = hh
| 1;
563 if (bMerge
&& daSwap
)
565 h
= hb
->d
.hydro
[id
][0];
567 /* Increment number if HBonds per H */
568 if (ihb
== hbHB
&& !bContact
)
570 inc_nhbonds(&(hb
->d
), d
, h
);
574 static char *mkatomname(const t_atoms
*atoms
, int i
)
579 rnr
= atoms
->atom
[i
].resind
;
580 sprintf(buf
, "%4s%d%-4s",
581 *atoms
->resinfo
[rnr
].name
, atoms
->resinfo
[rnr
].nr
, *atoms
->atomname
[i
]);
586 static void gen_datable(int *index
, int isize
, unsigned char *datable
, int natoms
)
588 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
591 for (i
= 0; i
< isize
; i
++)
593 if (index
[i
] >= natoms
)
595 gmx_fatal(FARGS
, "Atom has index %d larger than number of atoms %d.", index
[i
], natoms
);
597 datable
[index
[i
]] |= INGROUP
;
601 static void clear_datable_grp(unsigned char *datable
, int size
)
603 /* Clears group information from the table */
605 const char mask
= !(char)INGROUP
;
608 for (i
= 0; i
< size
; i
++)
615 static void add_acc(t_acceptors
*a
, int ia
, int grp
)
617 if (a
->nra
>= a
->max_nra
)
620 srenew(a
->acc
, a
->max_nra
);
621 srenew(a
->grp
, a
->max_nra
);
623 a
->grp
[a
->nra
] = grp
;
624 a
->acc
[a
->nra
++] = ia
;
627 static void search_acceptors(const t_topology
*top
, int isize
,
628 const int *index
, t_acceptors
*a
, int grp
,
630 gmx_bool bContact
, gmx_bool bDoIt
, unsigned char *datable
)
636 for (i
= 0; (i
< isize
); i
++)
640 (((*top
->atoms
.atomname
[n
])[0] == 'O') ||
641 (bNitAcc
&& ((*top
->atoms
.atomname
[n
])[0] == 'N')))) &&
644 datable
[n
] |= ACC
; /* set the atom's acceptor flag in datable. */
649 snew(a
->aptr
, top
->atoms
.nr
);
650 for (i
= 0; (i
< top
->atoms
.nr
); i
++)
654 for (i
= 0; (i
< a
->nra
); i
++)
656 a
->aptr
[a
->acc
[i
]] = i
;
660 static void add_h2d(int id
, int ih
, t_donors
*ddd
)
664 for (i
= 0; (i
< ddd
->nhydro
[id
]); i
++)
666 if (ddd
->hydro
[id
][i
] == ih
)
668 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
673 if (i
== ddd
->nhydro
[id
])
675 if (ddd
->nhydro
[id
] >= MAXHYDRO
)
677 gmx_fatal(FARGS
, "Donor %d has more than %d hydrogens!",
678 ddd
->don
[id
], MAXHYDRO
);
680 ddd
->hydro
[id
][i
] = ih
;
685 static void add_dh(t_donors
*ddd
, int id
, int ih
, int grp
, unsigned char *datable
)
689 if (!datable
|| ISDON(datable
[id
]))
691 if (ddd
->dptr
[id
] == NOTSET
) /* New donor */
703 if (ddd
->nrd
>= ddd
->max_nrd
)
706 srenew(ddd
->don
, ddd
->max_nrd
);
707 srenew(ddd
->nhydro
, ddd
->max_nrd
);
708 srenew(ddd
->hydro
, ddd
->max_nrd
);
709 srenew(ddd
->nhbonds
, ddd
->max_nrd
);
710 srenew(ddd
->grp
, ddd
->max_nrd
);
712 ddd
->don
[ddd
->nrd
] = id
;
713 ddd
->nhydro
[ddd
->nrd
] = 0;
714 ddd
->grp
[ddd
->nrd
] = grp
;
726 printf("Warning: Atom %d is not in the d/a-table!\n", id
);
730 static void search_donors(const t_topology
*top
, int isize
, const int *index
,
731 t_donors
*ddd
, int grp
, gmx_bool bContact
, gmx_bool bDoIt
,
732 unsigned char *datable
)
735 t_functype func_type
;
740 snew(ddd
->dptr
, top
->atoms
.nr
);
741 for (i
= 0; (i
< top
->atoms
.nr
); i
++)
743 ddd
->dptr
[i
] = NOTSET
;
751 for (i
= 0; (i
< isize
); i
++)
753 datable
[index
[i
]] |= DON
;
754 add_dh(ddd
, index
[i
], -1, grp
, datable
);
760 for (func_type
= 0; (func_type
< F_NRE
); func_type
++)
762 const t_ilist
*interaction
= &(top
->idef
.il
[func_type
]);
763 if (func_type
== F_POSRES
|| func_type
== F_FBPOSRES
)
765 /* The ilist looks strange for posre. Bug in grompp?
766 * We don't need posre interactions for hbonds anyway.*/
769 for (i
= 0; i
< interaction
->nr
;
770 i
+= interaction_function
[top
->idef
.functype
[interaction
->iatoms
[i
]]].nratoms
+1)
773 if (func_type
!= top
->idef
.functype
[interaction
->iatoms
[i
]])
775 fprintf(stderr
, "Error in func_type %s",
776 interaction_function
[func_type
].longname
);
780 /* check out this functype */
781 if (func_type
== F_SETTLE
)
783 nr1
= interaction
->iatoms
[i
+1];
784 nr2
= interaction
->iatoms
[i
+2];
785 nr3
= interaction
->iatoms
[i
+3];
787 if (ISINGRP(datable
[nr1
]))
789 if (ISINGRP(datable
[nr2
]))
792 add_dh(ddd
, nr1
, nr1
+1, grp
, datable
);
794 if (ISINGRP(datable
[nr3
]))
797 add_dh(ddd
, nr1
, nr1
+2, grp
, datable
);
801 else if (IS_CHEMBOND(func_type
))
803 for (j
= 0; j
< 2; j
++)
805 nr1
= interaction
->iatoms
[i
+1+j
];
806 nr2
= interaction
->iatoms
[i
+2-j
];
807 if ((*top
->atoms
.atomname
[nr1
][0] == 'H') &&
808 ((*top
->atoms
.atomname
[nr2
][0] == 'O') ||
809 (*top
->atoms
.atomname
[nr2
][0] == 'N')) &&
810 ISINGRP(datable
[nr1
]) && ISINGRP(datable
[nr2
]))
813 add_dh(ddd
, nr2
, nr1
, grp
, datable
);
820 for (func_type
= 0; func_type
< F_NRE
; func_type
++)
822 interaction
= &top
->idef
.il
[func_type
];
823 for (i
= 0; i
< interaction
->nr
;
824 i
+= interaction_function
[top
->idef
.functype
[interaction
->iatoms
[i
]]].nratoms
+1)
827 if (func_type
!= top
->idef
.functype
[interaction
->iatoms
[i
]])
829 gmx_incons("function type in search_donors");
832 if (interaction_function
[func_type
].flags
& IF_VSITE
)
834 nr1
= interaction
->iatoms
[i
+1];
835 if (*top
->atoms
.atomname
[nr1
][0] == 'H')
839 while (!stop
&& ( *top
->atoms
.atomname
[nr2
][0] == 'H'))
850 if (!stop
&& ( ( *top
->atoms
.atomname
[nr2
][0] == 'O') ||
851 ( *top
->atoms
.atomname
[nr2
][0] == 'N') ) &&
852 ISINGRP(datable
[nr1
]) && ISINGRP(datable
[nr2
]))
855 add_dh(ddd
, nr2
, nr1
, grp
, datable
);
865 static t_gridcell
***init_grid(gmx_bool bBox
, rvec box
[], real rcut
, ivec ngrid
)
872 for (i
= 0; i
< DIM
; i
++)
874 ngrid
[i
] = static_cast<int>(box
[i
][i
]/(1.2*rcut
));
878 if (!bBox
|| (ngrid
[XX
] < 3) || (ngrid
[YY
] < 3) || (ngrid
[ZZ
] < 3) )
880 for (i
= 0; i
< DIM
; i
++)
887 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
888 ngrid
[XX
], ngrid
[YY
], ngrid
[ZZ
], rcut
);
890 snew(grid
, ngrid
[ZZ
]);
891 for (z
= 0; z
< ngrid
[ZZ
]; z
++)
893 snew((grid
)[z
], ngrid
[YY
]);
894 for (y
= 0; y
< ngrid
[YY
]; y
++)
896 snew((grid
)[z
][y
], ngrid
[XX
]);
902 static void reset_nhbonds(t_donors
*ddd
)
906 for (i
= 0; (i
< ddd
->nrd
); i
++)
908 for (j
= 0; (j
< MAXHH
); j
++)
910 ddd
->nhbonds
[i
][j
] = 0;
915 void pbc_correct_gem(rvec dx
, matrix box
, rvec hbox
);
916 void pbc_in_gridbox(rvec dx
, matrix box
);
918 static void build_grid(t_hbdata
*hb
, rvec x
[], rvec xshell
,
919 gmx_bool bBox
, matrix box
, rvec hbox
,
920 real rcut
, real rshell
,
921 ivec ngrid
, t_gridcell
***grid
)
923 int i
, m
, gr
, xi
, yi
, zi
, nr
;
926 rvec invdelta
, dshell
;
928 gmx_bool bDoRshell
, bInShell
, bAcc
;
933 bDoRshell
= (rshell
> 0);
934 rshell2
= gmx::square(rshell
);
937 #define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
939 for (m
= 0; m
< DIM
; m
++)
941 hbox
[m
] = box
[m
][m
]*0.5;
944 invdelta
[m
] = ngrid
[m
]/box
[m
][m
];
945 if (1/invdelta
[m
] < rcut
)
947 gmx_fatal(FARGS
, "Your computational box has shrunk too much.\n"
948 "%s can not handle this situation, sorry.\n",
949 gmx::getProgramContext().displayName());
961 /* resetting atom counts */
962 for (gr
= 0; (gr
< grNR
); gr
++)
964 for (zi
= 0; zi
< ngrid
[ZZ
]; zi
++)
966 for (yi
= 0; yi
< ngrid
[YY
]; yi
++)
968 for (xi
= 0; xi
< ngrid
[XX
]; xi
++)
970 grid
[zi
][yi
][xi
].d
[gr
].nr
= 0;
971 grid
[zi
][yi
][xi
].a
[gr
].nr
= 0;
977 /* put atoms in grid cells */
978 for (bAcc
= FALSE
; (bAcc
<= TRUE
); bAcc
++)
991 for (i
= 0; (i
< nr
); i
++)
993 /* check if we are inside the shell */
994 /* if bDoRshell=FALSE then bInShell=TRUE always */
999 rvec_sub(x
[ad
[i
]], xshell
, dshell
);
1002 gmx_bool bDone
= FALSE
;
1006 for (m
= DIM
-1; m
>= 0 && bInShell
; m
--)
1008 if (dshell
[m
] < -hbox
[m
])
1011 rvec_inc(dshell
, box
[m
]);
1013 if (dshell
[m
] >= hbox
[m
])
1016 dshell
[m
] -= 2*hbox
[m
];
1020 for (m
= DIM
-1; m
>= 0 && bInShell
; m
--)
1022 /* if we're outside the cube, we're outside the sphere also! */
1023 if ( (dshell
[m
] > rshell
) || (-dshell
[m
] > rshell
) )
1029 /* if we're inside the cube, check if we're inside the sphere */
1032 bInShell
= norm2(dshell
) < rshell2
;
1040 pbc_in_gridbox(x
[ad
[i
]], box
);
1042 for (m
= DIM
-1; m
>= 0; m
--)
1043 { /* determine grid index of atom */
1044 grididx
[m
] = static_cast<int>(x
[ad
[i
]][m
]*invdelta
[m
]);
1045 grididx
[m
] = (grididx
[m
]+ngrid
[m
]) % ngrid
[m
];
1052 range_check(gx
, 0, ngrid
[XX
]);
1053 range_check(gy
, 0, ngrid
[YY
]);
1054 range_check(gz
, 0, ngrid
[ZZ
]);
1058 /* add atom to grid cell */
1061 newgrid
= &(grid
[gz
][gy
][gx
].a
[gr
]);
1065 newgrid
= &(grid
[gz
][gy
][gx
].d
[gr
]);
1067 if (newgrid
->nr
>= newgrid
->maxnr
)
1069 newgrid
->maxnr
+= 10;
1070 DBB(newgrid
->maxnr
);
1071 srenew(newgrid
->atoms
, newgrid
->maxnr
);
1074 newgrid
->atoms
[newgrid
->nr
] = ad
[i
];
1082 static void count_da_grid(ivec ngrid
, t_gridcell
***grid
, t_icell danr
)
1086 for (gr
= 0; (gr
< grNR
); gr
++)
1089 for (zi
= 0; zi
< ngrid
[ZZ
]; zi
++)
1091 for (yi
= 0; yi
< ngrid
[YY
]; yi
++)
1093 for (xi
= 0; xi
< ngrid
[XX
]; xi
++)
1095 danr
[gr
] += grid
[zi
][yi
][xi
].d
[gr
].nr
;
1103 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1104 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1105 * With a triclinic box all loops are 3 long, except when a cell is
1106 * located next to one of the box edges which is not parallel to the
1107 * x/y-plane, in that case all cells in a line or layer are searched.
1108 * This could be implemented slightly more efficient, but the code
1109 * would get much more complicated.
1111 static gmx_inline gmx_bool
grid_loop_begin(int n
, int x
, gmx_bool bTric
, gmx_bool bEdge
)
1113 return ((n
== 1) ? x
: bTric
&& bEdge
? 0 : (x
-1));
1115 static gmx_inline gmx_bool
grid_loop_end(int n
, int x
, gmx_bool bTric
, gmx_bool bEdge
)
1117 return ((n
== 1) ? x
: bTric
&& bEdge
? (n
-1) : (x
+1));
1119 static gmx_inline
int grid_mod(int j
, int n
)
1124 static void dump_grid(FILE *fp
, ivec ngrid
, t_gridcell
***grid
)
1126 int gr
, x
, y
, z
, sum
[grNR
];
1128 fprintf(fp
, "grid %dx%dx%d\n", ngrid
[XX
], ngrid
[YY
], ngrid
[ZZ
]);
1129 for (gr
= 0; gr
< grNR
; gr
++)
1132 fprintf(fp
, "GROUP %d (%s)\n", gr
, grpnames
[gr
]);
1133 for (z
= 0; z
< ngrid
[ZZ
]; z
+= 2)
1135 fprintf(fp
, "Z=%d,%d\n", z
, z
+1);
1136 for (y
= 0; y
< ngrid
[YY
]; y
++)
1138 for (x
= 0; x
< ngrid
[XX
]; x
++)
1140 fprintf(fp
, "%3d", grid
[x
][y
][z
].d
[gr
].nr
);
1141 sum
[gr
] += grid
[z
][y
][x
].d
[gr
].nr
;
1142 fprintf(fp
, "%3d", grid
[x
][y
][z
].a
[gr
].nr
);
1143 sum
[gr
] += grid
[z
][y
][x
].a
[gr
].nr
;
1147 if ( (z
+1) < ngrid
[ZZ
])
1149 for (x
= 0; x
< ngrid
[XX
]; x
++)
1151 fprintf(fp
, "%3d", grid
[z
+1][y
][x
].d
[gr
].nr
);
1152 sum
[gr
] += grid
[z
+1][y
][x
].d
[gr
].nr
;
1153 fprintf(fp
, "%3d", grid
[z
+1][y
][x
].a
[gr
].nr
);
1154 sum
[gr
] += grid
[z
+1][y
][x
].a
[gr
].nr
;
1161 fprintf(fp
, "TOTALS:");
1162 for (gr
= 0; gr
< grNR
; gr
++)
1164 fprintf(fp
, " %d=%d", gr
, sum
[gr
]);
1169 /* New GMX record! 5 * in a row. Congratulations!
1170 * Sorry, only four left.
1172 static void free_grid(ivec ngrid
, t_gridcell
****grid
)
1175 t_gridcell
***g
= *grid
;
1177 for (z
= 0; z
< ngrid
[ZZ
]; z
++)
1179 for (y
= 0; y
< ngrid
[YY
]; y
++)
1189 void pbc_correct_gem(rvec dx
, matrix box
, rvec hbox
)
1192 gmx_bool bDone
= FALSE
;
1196 for (m
= DIM
-1; m
>= 0; m
--)
1198 if (dx
[m
] < -hbox
[m
])
1201 rvec_inc(dx
, box
[m
]);
1203 if (dx
[m
] >= hbox
[m
])
1206 rvec_dec(dx
, box
[m
]);
1212 void pbc_in_gridbox(rvec dx
, matrix box
)
1215 gmx_bool bDone
= FALSE
;
1219 for (m
= DIM
-1; m
>= 0; m
--)
1224 rvec_inc(dx
, box
[m
]);
1226 if (dx
[m
] >= box
[m
][m
])
1229 rvec_dec(dx
, box
[m
]);
1235 /* Added argument r2cut, changed contact and implemented
1236 * use of second cut-off.
1237 * - Erik Marklund, June 29, 2006
1239 static int is_hbond(t_hbdata
*hb
, int grpd
, int grpa
, int d
, int a
,
1240 real rcut
, real r2cut
, real ccut
,
1241 rvec x
[], gmx_bool bBox
, matrix box
, rvec hbox
,
1242 real
*d_ha
, real
*ang
, gmx_bool bDA
, int *hhh
,
1243 gmx_bool bContact
, gmx_bool bMerge
)
1246 rvec r_da
, r_ha
, r_dh
;
1247 real rc2
, r2c2
, rda2
, rha2
, ca
;
1248 gmx_bool HAinrange
= FALSE
; /* If !bDA. Needed for returning hbDist in a correct way. */
1249 gmx_bool daSwap
= FALSE
;
1256 if (((id
= donor_index(&hb
->d
, grpd
, d
)) == NOTSET
) ||
1257 (acceptor_index(&hb
->a
, grpa
, a
) == NOTSET
))
1265 rvec_sub(x
[d
], x
[a
], r_da
);
1266 /* Insert projection code here */
1268 if (bMerge
&& d
> a
&& isInterchangable(hb
, d
, a
, grpd
, grpa
))
1270 /* Then this hbond/contact will be found again, or it has already been found. */
1275 if (d
> a
&& bMerge
&& isInterchangable(hb
, d
, a
, grpd
, grpa
)) /* acceptor is also a donor and vice versa? */
1276 { /* return hbNo; */
1277 daSwap
= TRUE
; /* If so, then their history should be filed with donor and acceptor swapped. */
1279 pbc_correct_gem(r_da
, box
, hbox
);
1281 rda2
= iprod(r_da
, r_da
);
1285 if (daSwap
&& grpa
== grpd
)
1293 else if (rda2
< r2c2
)
1304 if (bDA
&& (rda2
> rc2
))
1309 for (h
= 0; (h
< hb
->d
.nhydro
[id
]); h
++)
1311 hh
= hb
->d
.hydro
[id
][h
];
1315 rvec_sub(x
[hh
], x
[a
], r_ha
);
1318 pbc_correct_gem(r_ha
, box
, hbox
);
1320 rha2
= iprod(r_ha
, r_ha
);
1323 if (bDA
|| (rha2
<= rc2
))
1325 rvec_sub(x
[d
], x
[hh
], r_dh
);
1328 pbc_correct_gem(r_dh
, box
, hbox
);
1335 ca
= cos_angle(r_dh
, r_da
);
1336 /* if angle is smaller, cos is larger */
1340 *d_ha
= std::sqrt(bDA
? rda2
: rha2
);
1341 *ang
= std::acos(ca
);
1346 if (bDA
|| HAinrange
)
1356 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1357 * Will do some more testing before removing the function entirely.
1358 * - Erik Marklund, MAY 10 2010 */
1359 static void do_merge(t_hbdata
*hb
, int ntmp
,
1360 unsigned int htmp
[], unsigned int gtmp
[],
1361 t_hbond
*hb0
, t_hbond
*hb1
)
1363 /* Here we need to make sure we're treating periodicity in
1364 * the right way for the geminate recombination kinetics. */
1366 int m
, mm
, n00
, n01
, nn0
, nnframes
;
1368 /* Decide where to start from when merging */
1371 nn0
= std::min(n00
, n01
);
1372 nnframes
= std::max(n00
+ hb0
->nframes
, n01
+ hb1
->nframes
) - nn0
;
1373 /* Initiate tmp arrays */
1374 for (m
= 0; (m
< ntmp
); m
++)
1379 /* Fill tmp arrays with values due to first HB */
1380 /* Once again '<' had to be replaced with '<='
1381 to catch the last frame in which the hbond
1383 - Erik Marklund, June 1, 2006 */
1384 for (m
= 0; (m
<= hb0
->nframes
); m
++)
1387 htmp
[mm
] = is_hb(hb0
->h
[0], m
);
1389 for (m
= 0; (m
<= hb0
->nframes
); m
++)
1392 gtmp
[mm
] = is_hb(hb0
->g
[0], m
);
1395 for (m
= 0; (m
<= hb1
->nframes
); m
++)
1398 htmp
[mm
] = htmp
[mm
] || is_hb(hb1
->h
[0], m
);
1399 gtmp
[mm
] = gtmp
[mm
] || is_hb(hb1
->g
[0], m
);
1401 /* Reallocate target array */
1402 if (nnframes
> hb0
->maxframes
)
1404 srenew(hb0
->h
[0], 4+nnframes
/hb
->wordlen
);
1405 srenew(hb0
->g
[0], 4+nnframes
/hb
->wordlen
);
1408 /* Copy temp array to target array */
1409 for (m
= 0; (m
<= nnframes
); m
++)
1411 _set_hb(hb0
->h
[0], m
, htmp
[m
]);
1412 _set_hb(hb0
->g
[0], m
, gtmp
[m
]);
1415 /* Set scalar variables */
1417 hb0
->maxframes
= nnframes
;
1420 static void merge_hb(t_hbdata
*hb
, gmx_bool bTwo
, gmx_bool bContact
)
1422 int i
, inrnew
, indnew
, j
, ii
, jj
, id
, ia
, ntmp
;
1423 unsigned int *htmp
, *gtmp
;
1427 indnew
= hb
->nrdist
;
1429 /* Check whether donors are also acceptors */
1430 printf("Merging hbonds with Acceptor and Donor swapped\n");
1432 ntmp
= 2*hb
->max_frames
;
1435 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
1437 fprintf(stderr
, "\r%d/%d", i
+1, hb
->d
.nrd
);
1439 ii
= hb
->a
.aptr
[id
];
1440 for (j
= 0; (j
< hb
->a
.nra
); j
++)
1443 jj
= hb
->d
.dptr
[ia
];
1444 if ((id
!= ia
) && (ii
!= NOTSET
) && (jj
!= NOTSET
) &&
1445 (!bTwo
|| (hb
->d
.grp
[i
] != hb
->a
.grp
[j
])))
1447 hb0
= hb
->hbmap
[i
][j
];
1448 hb1
= hb
->hbmap
[jj
][ii
];
1449 if (hb0
&& hb1
&& ISHB(hb0
->history
[0]) && ISHB(hb1
->history
[0]))
1451 do_merge(hb
, ntmp
, htmp
, gtmp
, hb0
, hb1
);
1452 if (ISHB(hb1
->history
[0]))
1456 else if (ISDIST(hb1
->history
[0]))
1463 gmx_incons("No contact history");
1467 gmx_incons("Neither hydrogen bond nor distance");
1473 hb1
->history
[0] = hbNo
;
1478 fprintf(stderr
, "\n");
1479 printf("- Reduced number of hbonds from %d to %d\n", hb
->nrhb
, inrnew
);
1480 printf("- Reduced number of distances from %d to %d\n", hb
->nrdist
, indnew
);
1482 hb
->nrdist
= indnew
;
1487 static void do_nhb_dist(FILE *fp
, t_hbdata
*hb
, real t
)
1489 int i
, j
, k
, n_bound
[MAXHH
], nbtot
;
1491 /* Set array to 0 */
1492 for (k
= 0; (k
< MAXHH
); k
++)
1496 /* Loop over possible donors */
1497 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
1499 for (j
= 0; (j
< hb
->d
.nhydro
[i
]); j
++)
1501 n_bound
[hb
->d
.nhbonds
[i
][j
]]++;
1504 fprintf(fp
, "%12.5e", t
);
1506 for (k
= 0; (k
< MAXHH
); k
++)
1508 fprintf(fp
, " %8d", n_bound
[k
]);
1509 nbtot
+= n_bound
[k
]*k
;
1511 fprintf(fp
, " %8d\n", nbtot
);
1514 static void do_hblife(const char *fn
, t_hbdata
*hb
, gmx_bool bMerge
, gmx_bool bContact
,
1515 const gmx_output_env_t
*oenv
)
1518 const char *leg
[] = { "p(t)", "t p(t)" };
1520 int i
, j
, j0
, k
, m
, nh
, ihb
, ohb
, nhydro
, ndump
= 0;
1521 int nframes
= hb
->nframes
;
1524 double sum
, integral
;
1527 snew(h
, hb
->maxhydro
);
1528 snew(histo
, nframes
+1);
1529 /* Total number of hbonds analyzed here */
1530 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
1532 for (k
= 0; (k
< hb
->a
.nra
); k
++)
1534 hbh
= hb
->hbmap
[i
][k
];
1552 for (m
= 0; (m
< hb
->maxhydro
); m
++)
1556 h
[nhydro
++] = bContact
? hbh
->g
[m
] : hbh
->h
[m
];
1560 for (nh
= 0; (nh
< nhydro
); nh
++)
1565 for (j
= 0; (j
<= hbh
->nframes
); j
++)
1567 ihb
= is_hb(h
[nh
], j
);
1568 if (debug
&& (ndump
< 10))
1570 fprintf(debug
, "%5d %5d\n", j
, ihb
);
1590 fprintf(stderr
, "\n");
1593 fp
= xvgropen(fn
, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv
), "()", oenv
);
1597 fp
= xvgropen(fn
, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv
), "()",
1601 xvgr_legend(fp
, asize(leg
), leg
, oenv
);
1603 while ((j0
> 0) && (histo
[j0
] == 0))
1608 for (i
= 0; (i
<= j0
); i
++)
1612 dt
= hb
->time
[1]-hb
->time
[0];
1615 for (i
= 1; (i
<= j0
); i
++)
1617 t
= hb
->time
[i
] - hb
->time
[0] - 0.5*dt
;
1618 x1
= t
*histo
[i
]/sum
;
1619 fprintf(fp
, "%8.3f %10.3e %10.3e\n", t
, histo
[i
]/sum
, x1
);
1624 printf("%s lifetime = %.2f ps\n", bContact
? "Contact" : "HB", integral
);
1625 printf("Note that the lifetime obtained in this manner is close to useless\n");
1626 printf("Use the -ac option instead and check the Forward lifetime\n");
1627 please_cite(stdout
, "Spoel2006b");
1632 static void dump_ac(t_hbdata
*hb
, gmx_bool oneHB
, int nDump
)
1635 int i
, j
, k
, m
, nd
, ihb
, idist
;
1636 int nframes
= hb
->nframes
;
1644 fp
= gmx_ffopen("debug-ac.xvg", "w");
1645 for (j
= 0; (j
< nframes
); j
++)
1647 fprintf(fp
, "%10.3f", hb
->time
[j
]);
1648 for (i
= nd
= 0; (i
< hb
->d
.nrd
) && (nd
< nDump
); i
++)
1650 for (k
= 0; (k
< hb
->a
.nra
) && (nd
< nDump
); k
++)
1654 hbh
= hb
->hbmap
[i
][k
];
1659 ihb
= is_hb(hbh
->h
[0], j
);
1660 idist
= is_hb(hbh
->g
[0], j
);
1666 for (m
= 0; (m
< hb
->maxhydro
) && !ihb
; m
++)
1668 ihb
= ihb
|| ((hbh
->h
[m
]) && is_hb(hbh
->h
[m
], j
));
1669 idist
= idist
|| ((hbh
->g
[m
]) && is_hb(hbh
->g
[m
], j
));
1671 /* This is not correct! */
1672 /* What isn't correct? -Erik M */
1677 fprintf(fp
, " %1d-%1d", ihb
, idist
);
1687 static real
calc_dg(real tau
, real temp
)
1698 return kbt
*std::log(kbt
*tau
/PLANCK
);
1703 int n0
, n1
, nparams
, ndelta
;
1705 real
*t
, *ct
, *nt
, *kt
, *sigma_ct
, *sigma_nt
, *sigma_kt
;
1708 static real
compute_weighted_rates(int n
, real t
[], real ct
[], real nt
[],
1709 real kt
[], real sigma_ct
[], real sigma_nt
[],
1710 real sigma_kt
[], real
*k
, real
*kp
,
1711 real
*sigma_k
, real
*sigma_kp
,
1717 real kkk
= 0, kkp
= 0, kk2
= 0, kp2
= 0, chi2
;
1722 for (i
= 0; (i
< n
); i
++)
1724 if (t
[i
] >= fit_start
)
1737 tl
.sigma_ct
= sigma_ct
;
1738 tl
.sigma_nt
= sigma_nt
;
1739 tl
.sigma_kt
= sigma_kt
;
1743 chi2
= 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
1745 *kp
= tl
.kkk
[1] = *kp
;
1747 for (j
= 0; (j
< NK
); j
++)
1751 kk2
+= gmx::square(tl
.kkk
[0]);
1752 kp2
+= gmx::square(tl
.kkk
[1]);
1755 *sigma_k
= std::sqrt(kk2
/NK
- gmx::square(kkk
/NK
));
1756 *sigma_kp
= std::sqrt(kp2
/NK
- gmx::square(kkp
/NK
));
1761 void analyse_corr(int n
, real t
[], real ct
[], real nt
[], real kt
[],
1762 real sigma_ct
[], real sigma_nt
[], real sigma_kt
[],
1763 real fit_start
, real temp
)
1766 real k
= 1, kp
= 1, kow
= 1;
1767 real Q
= 0, chi2
, tau_hb
, dtau
, tau_rlx
, e_1
, sigma_k
, sigma_kp
, ddg
;
1768 double tmp
, sn2
= 0, sc2
= 0, sk2
= 0, scn
= 0, sck
= 0, snk
= 0;
1769 gmx_bool bError
= (sigma_ct
!= NULL
) && (sigma_nt
!= NULL
) && (sigma_kt
!= NULL
);
1771 for (i0
= 0; (i0
< n
-2) && ((t
[i0
]-t
[0]) < fit_start
); i0
++)
1777 for (i
= i0
; (i
< n
); i
++)
1779 sc2
+= gmx::square(ct
[i
]);
1780 sn2
+= gmx::square(nt
[i
]);
1781 sk2
+= gmx::square(kt
[i
]);
1786 printf("Hydrogen bond thermodynamics at T = %g K\n", temp
);
1787 tmp
= (sn2
*sc2
-gmx::square(scn
));
1788 if ((tmp
> 0) && (sn2
> 0))
1790 k
= (sn2
*sck
-scn
*snk
)/tmp
;
1791 kp
= (k
*scn
-snk
)/sn2
;
1795 for (i
= i0
; (i
< n
); i
++)
1797 chi2
+= gmx::square(k
*ct
[i
]-kp
*nt
[i
]-kt
[i
]);
1799 compute_weighted_rates(n
, t
, ct
, nt
, kt
, sigma_ct
, sigma_nt
,
1801 &sigma_k
, &sigma_kp
, fit_start
);
1802 Q
= 0; /* quality_of_fit(chi2, 2);*/
1803 ddg
= BOLTZ
*temp
*sigma_k
/k
;
1804 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
1806 printf("The Rate and Delta G are followed by an error estimate\n");
1807 printf("----------------------------------------------------------\n"
1808 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
1809 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1810 k
, sigma_k
, 1/k
, calc_dg(1/k
, temp
), ddg
);
1811 ddg
= BOLTZ
*temp
*sigma_kp
/kp
;
1812 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
1813 kp
, sigma_kp
, 1/kp
, calc_dg(1/kp
, temp
), ddg
);
1818 for (i
= i0
; (i
< n
); i
++)
1820 chi2
+= gmx::square(k
*ct
[i
]-kp
*nt
[i
]-kt
[i
]);
1822 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
1824 printf("--------------------------------------------------\n"
1825 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
1826 printf("Forward %10.3f %8.3f %10.3f %10g\n",
1827 k
, 1/k
, calc_dg(1/k
, temp
), chi2
);
1828 printf("Backward %10.3f %8.3f %10.3f\n",
1829 kp
, 1/kp
, calc_dg(1/kp
, temp
));
1835 printf("One-way %10.3f %s%8.3f %10.3f\n",
1836 kow
, bError
? " " : "", 1/kow
, calc_dg(1/kow
, temp
));
1840 printf(" - Numerical problems computing HB thermodynamics:\n"
1841 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
1842 sc2
, sn2
, sk2
, sck
, snk
, scn
);
1844 /* Determine integral of the correlation function */
1845 tau_hb
= evaluate_integral(n
, t
, ct
, NULL
, (t
[n
-1]-t
[0])/2, &dtau
);
1846 printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb
,
1847 bError
? " " : "", tau_hb
, calc_dg(tau_hb
, temp
));
1848 e_1
= std::exp(-1.0);
1849 for (i
= 0; (i
< n
-2); i
++)
1851 if ((ct
[i
] > e_1
) && (ct
[i
+1] <= e_1
))
1858 /* Determine tau_relax from linear interpolation */
1859 tau_rlx
= t
[i
]-t
[0] + (e_1
-ct
[i
])*(t
[i
+1]-t
[i
])/(ct
[i
+1]-ct
[i
]);
1860 printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx
,
1861 tau_rlx
, bError
? " " : "",
1862 calc_dg(tau_rlx
, temp
));
1867 printf("Correlation functions too short to compute thermodynamics\n");
1871 void compute_derivative(int nn
, real x
[], real y
[], real dydx
[])
1875 /* Compute k(t) = dc(t)/dt */
1876 for (j
= 1; (j
< nn
-1); j
++)
1878 dydx
[j
] = (y
[j
+1]-y
[j
-1])/(x
[j
+1]-x
[j
-1]);
1880 /* Extrapolate endpoints */
1881 dydx
[0] = 2*dydx
[1] - dydx
[2];
1882 dydx
[nn
-1] = 2*dydx
[nn
-2] - dydx
[nn
-3];
1885 static void normalizeACF(real
*ct
, real
*gt
, int nhb
, int len
)
1887 real ct_fac
, gt_fac
= 0;
1890 /* Xu and Berne use the same normalization constant */
1898 printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac
, gt_fac
);
1899 for (i
= 0; i
< len
; i
++)
1909 static void do_hbac(const char *fn
, t_hbdata
*hb
,
1910 int nDump
, gmx_bool bMerge
, gmx_bool bContact
, real fit_start
,
1911 real temp
, gmx_bool R2
, const gmx_output_env_t
*oenv
,
1915 int i
, j
, k
, m
, ihb
, idist
, n2
, nn
;
1917 const char *legLuzar
[] = {
1918 "Ac\\sfin sys\\v{}\\z{}(t)",
1920 "Cc\\scontact,hb\\v{}\\z{}(t)",
1921 "-dAc\\sfs\\v{}\\z{}/dt"
1923 gmx_bool bNorm
= FALSE
;
1925 real
*rhbex
= NULL
, *ht
, *gt
, *ght
, *dght
, *kt
;
1926 real
*ct
, tail
, tail2
, dtail
, *cct
;
1927 const real tol
= 1e-3;
1928 int nframes
= hb
->nframes
;
1929 unsigned int **h
= NULL
, **g
= NULL
;
1930 int nh
, nhbonds
, nhydro
;
1933 int *dondata
= NULL
;
1936 AC_NONE
, AC_NN
, AC_GEM
, AC_LUZAR
1939 const bool bOMP
= GMX_OPENMP
;
1941 printf("Doing autocorrelation ");
1944 printf("according to the theory of Luzar and Chandler.\n");
1947 /* build hbexist matrix in reals for autocorr */
1948 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
1950 while (n2
< nframes
)
1957 if (acType
!= AC_NN
|| bOMP
)
1959 snew(h
, hb
->maxhydro
);
1960 snew(g
, hb
->maxhydro
);
1963 /* Dump hbonds for debugging */
1964 dump_ac(hb
, bMerge
|| bContact
, nDump
);
1966 /* Total number of hbonds analyzed here */
1969 if (acType
!= AC_LUZAR
&& bOMP
)
1971 nThreads
= std::min((nThreads
<= 0) ? INT_MAX
: nThreads
, gmx_omp_get_max_threads());
1973 gmx_omp_set_num_threads(nThreads
);
1974 snew(dondata
, nThreads
);
1975 for (i
= 0; i
< nThreads
; i
++)
1979 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
1980 "Expect close to linear scaling over this donor-loop.\n", nThreads
);
1982 fprintf(stderr
, "Donors: [thread no]\n");
1985 for (i
= 0; i
< nThreads
; i
++)
1987 snprintf(tmpstr
, 7, "[%i]", i
);
1988 fprintf(stderr
, "%-7s", tmpstr
);
1991 fprintf(stderr
, "\n");
2006 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2008 for (k
= 0; (k
< hb
->a
.nra
); k
++)
2011 hbh
= hb
->hbmap
[i
][k
];
2015 if (bMerge
|| bContact
)
2017 if (ISHB(hbh
->history
[0]))
2026 for (m
= 0; (m
< hb
->maxhydro
); m
++)
2028 if (bContact
? ISDIST(hbh
->history
[m
]) : ISHB(hbh
->history
[m
]))
2030 g
[nhydro
] = hbh
->g
[m
];
2031 h
[nhydro
] = hbh
->h
[m
];
2037 int nf
= hbh
->nframes
;
2038 for (nh
= 0; (nh
< nhydro
); nh
++)
2040 int nrint
= bContact
? hb
->nrdist
: hb
->nrhb
;
2041 if ((((nhbonds
+1) % 10) == 0) || (nhbonds
+1 == nrint
))
2043 fprintf(stderr
, "\rACF %d/%d", nhbonds
+1, nrint
);
2046 for (j
= 0; (j
< nframes
); j
++)
2050 ihb
= is_hb(h
[nh
], j
);
2051 idist
= is_hb(g
[nh
], j
);
2058 /* For contacts: if a second cut-off is provided, use it,
2059 * otherwise use g(t) = 1-h(t) */
2060 if (!R2
&& bContact
)
2066 gt
[j
] = idist
*(1-ihb
);
2072 /* The autocorrelation function is normalized after summation only */
2073 low_do_autocorr(NULL
, oenv
, NULL
, nframes
, 1, -1, &rhbex
, hb
->time
[1]-hb
->time
[0],
2074 eacNormal
, 1, FALSE
, bNorm
, FALSE
, 0, -1, 0);
2076 /* Cross correlation analysis for thermodynamics */
2077 for (j
= nframes
; (j
< n2
); j
++)
2083 cross_corr(n2
, ht
, gt
, dght
);
2085 for (j
= 0; (j
< nn
); j
++)
2094 fprintf(stderr
, "\n");
2097 normalizeACF(ct
, ght
, static_cast<int>(nhb
), nn
);
2099 /* Determine tail value for statistics */
2102 for (j
= nn
/2; (j
< nn
); j
++)
2105 tail2
+= ct
[j
]*ct
[j
];
2107 tail
/= (nn
- nn
/2);
2108 tail2
/= (nn
- nn
/2);
2109 dtail
= std::sqrt(tail2
-tail
*tail
);
2111 /* Check whether the ACF is long enough */
2114 printf("\nWARNING: Correlation function is probably not long enough\n"
2115 "because the standard deviation in the tail of C(t) > %g\n"
2116 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2119 for (j
= 0; (j
< nn
); j
++)
2122 ct
[j
] = (cct
[j
]-tail
)/(1-tail
);
2124 /* Compute negative derivative k(t) = -dc(t)/dt */
2125 compute_derivative(nn
, hb
->time
, ct
, kt
);
2126 for (j
= 0; (j
< nn
); j
++)
2134 fp
= xvgropen(fn
, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv
), "C(t)", oenv
);
2138 fp
= xvgropen(fn
, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv
), "C(t)", oenv
);
2140 xvgr_legend(fp
, asize(legLuzar
), legLuzar
, oenv
);
2143 for (j
= 0; (j
< nn
); j
++)
2145 fprintf(fp
, "%10g %10g %10g %10g %10g\n",
2146 hb
->time
[j
]-hb
->time
[0], ct
[j
], cct
[j
], ght
[j
], kt
[j
]);
2150 analyse_corr(nn
, hb
->time
, ct
, ght
, kt
, NULL
, NULL
, NULL
,
2153 do_view(oenv
, fn
, NULL
);
2164 static void init_hbframe(t_hbdata
*hb
, int nframes
, real t
)
2168 hb
->time
[nframes
] = t
;
2169 hb
->nhb
[nframes
] = 0;
2170 hb
->ndist
[nframes
] = 0;
2171 for (i
= 0; (i
< max_hx
); i
++)
2173 hb
->nhx
[nframes
][i
] = 0;
2177 static FILE *open_donor_properties_file(const char *fn
,
2179 const gmx_output_env_t
*oenv
)
2182 const char *leg
[] = { "Nbound", "Nfree" };
2189 fp
= xvgropen(fn
, "Donor properties", output_env_get_xvgr_tlabel(oenv
), "Number", oenv
);
2190 xvgr_legend(fp
, asize(leg
), leg
, oenv
);
2195 static void analyse_donor_properties(FILE *fp
, t_hbdata
*hb
, int nframes
, real t
)
2197 int i
, j
, k
, nbound
, nb
, nhtot
;
2205 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2207 for (k
= 0; (k
< hb
->d
.nhydro
[i
]); k
++)
2211 for (j
= 0; (j
< hb
->a
.nra
) && (nb
== 0); j
++)
2213 if (hb
->hbmap
[i
][j
] && hb
->hbmap
[i
][j
]->h
[k
] &&
2214 is_hb(hb
->hbmap
[i
][j
]->h
[k
], nframes
))
2222 fprintf(fp
, "%10.3e %6d %6d\n", t
, nbound
, nhtot
-nbound
);
2225 static void dump_hbmap(t_hbdata
*hb
,
2226 int nfile
, t_filenm fnm
[], gmx_bool bTwo
,
2227 gmx_bool bContact
, int isize
[], int *index
[], char *grpnames
[],
2228 const t_atoms
*atoms
)
2231 int ddd
, hhh
, aaa
, i
, j
, k
, m
, grp
;
2232 char ds
[32], hs
[32], as
[32];
2235 fp
= opt2FILE("-hbn", nfile
, fnm
, "w");
2236 if (opt2bSet("-g", nfile
, fnm
))
2238 fplog
= gmx_ffopen(opt2fn("-g", nfile
, fnm
), "w");
2239 fprintf(fplog
, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
2245 for (grp
= gr0
; grp
<= (bTwo
? gr1
: gr0
); grp
++)
2247 fprintf(fp
, "[ %s ]", grpnames
[grp
]);
2248 for (i
= 0; i
< isize
[grp
]; i
++)
2250 fprintf(fp
, (i
%15) ? " " : "\n");
2251 fprintf(fp
, " %4d", index
[grp
][i
]+1);
2257 fprintf(fp
, "[ donors_hydrogens_%s ]\n", grpnames
[grp
]);
2258 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2260 if (hb
->d
.grp
[i
] == grp
)
2262 for (j
= 0; (j
< hb
->d
.nhydro
[i
]); j
++)
2264 fprintf(fp
, " %4d %4d", hb
->d
.don
[i
]+1,
2265 hb
->d
.hydro
[i
][j
]+1);
2271 fprintf(fp
, "[ acceptors_%s ]", grpnames
[grp
]);
2272 for (i
= 0; (i
< hb
->a
.nra
); i
++)
2274 if (hb
->a
.grp
[i
] == grp
)
2276 fprintf(fp
, (i
%15 && !first
) ? " " : "\n");
2277 fprintf(fp
, " %4d", hb
->a
.acc
[i
]+1);
2286 fprintf(fp
, bContact
? "[ contacts_%s-%s ]\n" :
2287 "[ hbonds_%s-%s ]\n", grpnames
[0], grpnames
[1]);
2291 fprintf(fp
, bContact
? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames
[0]);
2294 for (i
= 0; (i
< hb
->d
.nrd
); i
++)
2297 for (k
= 0; (k
< hb
->a
.nra
); k
++)
2300 for (m
= 0; (m
< hb
->d
.nhydro
[i
]); m
++)
2302 if (hb
->hbmap
[i
][k
] && ISHB(hb
->hbmap
[i
][k
]->history
[m
]))
2304 sprintf(ds
, "%s", mkatomname(atoms
, ddd
));
2305 sprintf(as
, "%s", mkatomname(atoms
, aaa
));
2308 fprintf(fp
, " %6d %6d\n", ddd
+1, aaa
+1);
2311 fprintf(fplog
, "%12s %12s\n", ds
, as
);
2316 hhh
= hb
->d
.hydro
[i
][m
];
2317 sprintf(hs
, "%s", mkatomname(atoms
, hhh
));
2318 fprintf(fp
, " %6d %6d %6d\n", ddd
+1, hhh
+1, aaa
+1);
2321 fprintf(fplog
, "%12s %12s %12s\n", ds
, hs
, as
);
2335 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
2336 * It mimics add_frames() and init_frame() to some extent. */
2337 static void sync_hbdata(t_hbdata
*p_hb
, int nframes
)
2339 if (nframes
>= p_hb
->max_frames
)
2341 p_hb
->max_frames
+= 4096;
2342 srenew(p_hb
->nhb
, p_hb
->max_frames
);
2343 srenew(p_hb
->ndist
, p_hb
->max_frames
);
2344 srenew(p_hb
->n_bound
, p_hb
->max_frames
);
2345 srenew(p_hb
->nhx
, p_hb
->max_frames
);
2348 srenew(p_hb
->danr
, p_hb
->max_frames
);
2350 std::memset(&(p_hb
->nhb
[nframes
]), 0, sizeof(int) * (p_hb
->max_frames
-nframes
));
2351 std::memset(&(p_hb
->ndist
[nframes
]), 0, sizeof(int) * (p_hb
->max_frames
-nframes
));
2352 p_hb
->nhb
[nframes
] = 0;
2353 p_hb
->ndist
[nframes
] = 0;
2356 p_hb
->nframes
= nframes
;
2358 std::memset(&(p_hb
->nhx
[nframes
]), 0, sizeof(int)*max_hx
); /* zero the helix count for this frame */
2361 int gmx_hbond(int argc
, char *argv
[])
2363 const char *desc
[] = {
2364 "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
2365 "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
2366 "(zero is extended) and the distance Donor - Acceptor",
2367 "(or Hydrogen - Acceptor using [TT]-noda[tt]).",
2368 "OH and NH groups are regarded as donors, O is an acceptor always,",
2369 "N is an acceptor by default, but this can be switched using",
2370 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
2371 "to the first preceding non-hydrogen atom.[PAR]",
2373 "You need to specify two groups for analysis, which must be either",
2374 "identical or non-overlapping. All hydrogen bonds between the two",
2375 "groups are analyzed.[PAR]",
2377 "If you set [TT]-shell[tt], you will be asked for an additional index group",
2378 "which should contain exactly one atom. In this case, only hydrogen",
2379 "bonds between atoms within the shell distance from the one atom are",
2382 "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
2383 "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
2384 "If contact kinetics are analyzed by using the -contact option, then",
2385 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
2386 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
2387 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
2388 "See mentioned literature for more details and definitions."
2391 /* "It is also possible to analyse specific hydrogen bonds with",
2392 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
2393 "Donor Hydrogen Acceptor, in the following way::",
2400 "Note that the triplets need not be on separate lines.",
2401 "Each atom triplet specifies a hydrogen bond to be analyzed,",
2402 "note also that no check is made for the types of atoms.[PAR]",
2407 " * [TT]-num[tt]: number of hydrogen bonds as a function of time.",
2408 " * [TT]-ac[tt]: average over all autocorrelations of the existence",
2409 " functions (either 0 or 1) of all hydrogen bonds.",
2410 " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
2411 " * [TT]-ang[tt]: angle distribution of all hydrogen bonds.",
2412 " * [TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
2413 " where n and n+i stand for residue numbers and i ranges from 0 to 6.",
2414 " This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
2415 " with helices in proteins.",
2416 " * [TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
2417 " for selected groups, all hydrogen bonded atoms from all groups and",
2418 " all solvent atoms involved in insertion.",
2419 " * [TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
2420 " frames, this also contains information on solvent insertion",
2421 " into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
2423 " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
2424 " each timeframe. This is especially useful when using [TT]-shell[tt].",
2425 " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
2426 " compare results to Raman Spectroscopy.",
2428 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
2429 "require an amount of memory proportional to the total numbers of donors",
2430 "times the total number of acceptors in the selected group(s)."
2433 static real acut
= 30, abin
= 1, rcut
= 0.35, r2cut
= 0, rbin
= 0.005, rshell
= -1;
2434 static real maxnhb
= 0, fit_start
= 1, fit_end
= 60, temp
= 298.15;
2435 static gmx_bool bNitAcc
= TRUE
, bDA
= TRUE
, bMerge
= TRUE
;
2436 static int nDump
= 0;
2437 static int nThreads
= 0;
2439 static gmx_bool bContact
= FALSE
;
2443 { "-a", FALSE
, etREAL
, {&acut
},
2444 "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
2445 { "-r", FALSE
, etREAL
, {&rcut
},
2446 "Cutoff radius (nm, X - Acceptor, see next option)" },
2447 { "-da", FALSE
, etBOOL
, {&bDA
},
2448 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
2449 { "-r2", FALSE
, etREAL
, {&r2cut
},
2450 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
2451 { "-abin", FALSE
, etREAL
, {&abin
},
2452 "Binwidth angle distribution (degrees)" },
2453 { "-rbin", FALSE
, etREAL
, {&rbin
},
2454 "Binwidth distance distribution (nm)" },
2455 { "-nitacc", FALSE
, etBOOL
, {&bNitAcc
},
2456 "Regard nitrogen atoms as acceptors" },
2457 { "-contact", FALSE
, etBOOL
, {&bContact
},
2458 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
2459 { "-shell", FALSE
, etREAL
, {&rshell
},
2460 "when > 0, only calculate hydrogen bonds within # nm shell around "
2462 { "-fitstart", FALSE
, etREAL
, {&fit_start
},
2463 "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]" },
2464 { "-fitend", FALSE
, etREAL
, {&fit_end
},
2465 "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])" },
2466 { "-temp", FALSE
, etREAL
, {&temp
},
2467 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
2468 { "-dump", FALSE
, etINT
, {&nDump
},
2469 "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
2470 { "-max_hb", FALSE
, etREAL
, {&maxnhb
},
2471 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
2472 { "-merge", FALSE
, etBOOL
, {&bMerge
},
2473 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
2475 { "-nthreads", FALSE
, etINT
, {&nThreads
},
2476 "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)"},
2479 const char *bugs
[] = {
2480 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
2483 { efTRX
, "-f", NULL
, ffREAD
},
2484 { efTPR
, NULL
, NULL
, ffREAD
},
2485 { efNDX
, NULL
, NULL
, ffOPTRD
},
2486 /* { efNDX, "-sel", "select", ffOPTRD },*/
2487 { efXVG
, "-num", "hbnum", ffWRITE
},
2488 { efLOG
, "-g", "hbond", ffOPTWR
},
2489 { efXVG
, "-ac", "hbac", ffOPTWR
},
2490 { efXVG
, "-dist", "hbdist", ffOPTWR
},
2491 { efXVG
, "-ang", "hbang", ffOPTWR
},
2492 { efXVG
, "-hx", "hbhelix", ffOPTWR
},
2493 { efNDX
, "-hbn", "hbond", ffOPTWR
},
2494 { efXPM
, "-hbm", "hbmap", ffOPTWR
},
2495 { efXVG
, "-don", "donor", ffOPTWR
},
2496 { efXVG
, "-dan", "danum", ffOPTWR
},
2497 { efXVG
, "-life", "hblife", ffOPTWR
},
2498 { efXVG
, "-nhbdist", "nhbdist", ffOPTWR
}
2501 #define NFILE asize(fnm)
2503 char hbmap
[HB_NR
] = { ' ', 'o', '-', '*' };
2504 const char *hbdesc
[HB_NR
] = { "None", "Present", "Inserted", "Present & Inserted" };
2505 t_rgb hbrgb
[HB_NR
] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
2507 t_trxstatus
*status
;
2512 int npargs
, natoms
, nframes
= 0, shatom
;
2518 real t
, ccut
, dist
= 0.0, ang
= 0.0;
2519 double max_nhb
, aver_nhb
, aver_dist
;
2520 int h
= 0, i
= 0, j
, k
= 0, ogrp
, nsel
;
2522 int xj
, yj
, zj
, aj
, xjj
, yjj
, zjj
;
2523 gmx_bool bSelected
, bHBmap
, bStop
, bTwo
, bBox
, bTric
;
2525 int grp
, nabin
, nrbin
, resdist
, ihb
;
2528 FILE *fp
, *fpnhb
= NULL
, *donor_properties
= NULL
;
2530 t_ncell
*icell
, *jcell
;
2532 unsigned char *datable
;
2533 gmx_output_env_t
*oenv
;
2534 int ii
, hh
, actual_nThreads
;
2537 gmx_bool bEdge_yjj
, bEdge_xjj
;
2539 t_hbdata
**p_hb
= NULL
; /* one per thread, then merge after the frame loop */
2540 int **p_adist
= NULL
, **p_rdist
= NULL
; /* a histogram for each thread. */
2542 const bool bOMP
= GMX_OPENMP
;
2545 ppa
= add_acf_pargs(&npargs
, pa
);
2547 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_TIME_UNIT
, NFILE
, fnm
, npargs
,
2548 ppa
, asize(desc
), desc
, asize(bugs
), bugs
, &oenv
))
2555 ccut
= std::cos(acut
*DEG2RAD
);
2561 gmx_fatal(FARGS
, "Can not analyze selected contacts.");
2565 gmx_fatal(FARGS
, "Can not analyze contact between H and A: turn off -noda");
2569 /* Initiate main data structure! */
2570 bHBmap
= (opt2bSet("-ac", NFILE
, fnm
) ||
2571 opt2bSet("-life", NFILE
, fnm
) ||
2572 opt2bSet("-hbn", NFILE
, fnm
) ||
2573 opt2bSet("-hbm", NFILE
, fnm
));
2575 if (opt2bSet("-nhbdist", NFILE
, fnm
))
2577 const char *leg
[MAXHH
+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
2578 fpnhb
= xvgropen(opt2fn("-nhbdist", NFILE
, fnm
),
2579 "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv
), "N", oenv
);
2580 xvgr_legend(fpnhb
, asize(leg
), leg
, oenv
);
2583 hb
= mk_hbdata(bHBmap
, opt2bSet("-dan", NFILE
, fnm
), bMerge
|| bContact
);
2586 read_tpx_top(ftp2fn(efTPR
, NFILE
, fnm
), &ir
, box
, &natoms
, NULL
, NULL
, &top
);
2588 snew(grpnames
, grNR
);
2591 /* Make Donor-Acceptor table */
2592 snew(datable
, top
.atoms
.nr
);
2593 gen_datable(index
[0], isize
[0], datable
, top
.atoms
.nr
);
2597 /* analyze selected hydrogen bonds */
2598 printf("Select group with selected atoms:\n");
2599 get_index(&(top
.atoms
), opt2fn("-sel", NFILE
, fnm
),
2600 1, &nsel
, index
, grpnames
);
2603 gmx_fatal(FARGS
, "Number of atoms in group '%s' not a multiple of 3\n"
2604 "and therefore cannot contain triplets of "
2605 "Donor-Hydrogen-Acceptor", grpnames
[0]);
2609 for (i
= 0; (i
< nsel
); i
+= 3)
2611 int dd
= index
[0][i
];
2612 int aa
= index
[0][i
+2];
2613 /* int */ hh
= index
[0][i
+1];
2614 add_dh (&hb
->d
, dd
, hh
, i
, datable
);
2615 add_acc(&hb
->a
, aa
, i
);
2616 /* Should this be here ? */
2617 snew(hb
->d
.dptr
, top
.atoms
.nr
);
2618 snew(hb
->a
.aptr
, top
.atoms
.nr
);
2619 add_hbond(hb
, dd
, aa
, hh
, gr0
, gr0
, 0, bMerge
, 0, bContact
);
2621 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
2622 isize
[0], grpnames
[0]);
2626 /* analyze all hydrogen bonds: get group(s) */
2627 printf("Specify 2 groups to analyze:\n");
2628 get_index(&(top
.atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
),
2629 2, isize
, index
, grpnames
);
2631 /* check if we have two identical or two non-overlapping groups */
2632 bTwo
= isize
[0] != isize
[1];
2633 for (i
= 0; (i
< isize
[0]) && !bTwo
; i
++)
2635 bTwo
= index
[0][i
] != index
[1][i
];
2639 printf("Checking for overlap in atoms between %s and %s\n",
2640 grpnames
[0], grpnames
[1]);
2641 for (i
= 0; i
< isize
[1]; i
++)
2643 if (ISINGRP(datable
[index
[1][i
]]))
2645 gmx_fatal(FARGS
, "Partial overlap between groups '%s' and '%s'",
2646 grpnames
[0], grpnames
[1]);
2650 printf("Checking for overlap in atoms between %s and %s\n",
2651 grpnames[0],grpnames[1]);
2652 for (i=0; i<isize[0]; i++)
2653 for (j=0; j<isize[1]; j++)
2654 if (index[0][i] == index[1][j])
2655 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
2656 grpnames[0],grpnames[1]);
2661 printf("Calculating %s "
2662 "between %s (%d atoms) and %s (%d atoms)\n",
2663 bContact
? "contacts" : "hydrogen bonds",
2664 grpnames
[0], isize
[0], grpnames
[1], isize
[1]);
2668 fprintf(stderr
, "Calculating %s in %s (%d atoms)\n",
2669 bContact
? "contacts" : "hydrogen bonds", grpnames
[0], isize
[0]);
2674 /* search donors and acceptors in groups */
2675 snew(datable
, top
.atoms
.nr
);
2676 for (i
= 0; (i
< grNR
); i
++)
2678 if ( ((i
== gr0
) && !bSelected
) ||
2679 ((i
== gr1
) && bTwo
))
2681 gen_datable(index
[i
], isize
[i
], datable
, top
.atoms
.nr
);
2684 search_acceptors(&top
, isize
[i
], index
[i
], &hb
->a
, i
,
2685 bNitAcc
, TRUE
, (bTwo
&& (i
== gr0
)) || !bTwo
, datable
);
2686 search_donors (&top
, isize
[i
], index
[i
], &hb
->d
, i
,
2687 TRUE
, (bTwo
&& (i
== gr1
)) || !bTwo
, datable
);
2691 search_acceptors(&top
, isize
[i
], index
[i
], &hb
->a
, i
, bNitAcc
, FALSE
, TRUE
, datable
);
2692 search_donors (&top
, isize
[i
], index
[i
], &hb
->d
, i
, FALSE
, TRUE
, datable
);
2696 clear_datable_grp(datable
, top
.atoms
.nr
);
2701 printf("Found %d donors and %d acceptors\n", hb
->d
.nrd
, hb
->a
.nra
);
2703 snew(donors[gr0D], dons[gr0D].nrd);*/
2705 donor_properties
= open_donor_properties_file(opt2fn_null("-don", NFILE
, fnm
), hb
, oenv
);
2709 printf("Making hbmap structure...");
2710 /* Generate hbond data structure */
2717 if (hb
->d
.nrd
+ hb
->a
.nra
== 0)
2719 printf("No Donors or Acceptors found\n");
2726 printf("No Donors found\n");
2731 printf("No Acceptors found\n");
2737 gmx_fatal(FARGS
, "Nothing to be done");
2746 /* get index group with atom for shell */
2749 printf("Select atom for shell (1 atom):\n");
2750 get_index(&(top
.atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
),
2751 1, &shisz
, &shidx
, &shgrpnm
);
2754 printf("group contains %d atoms, should be 1 (one)\n", shisz
);
2759 printf("Will calculate hydrogen bonds within a shell "
2760 "of %g nm around atom %i\n", rshell
, shatom
+1);
2763 /* Analyze trajectory */
2764 natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
2765 if (natoms
> top
.atoms
.nr
)
2767 gmx_fatal(FARGS
, "Topology (%d atoms) does not match trajectory (%d atoms)",
2768 top
.atoms
.nr
, natoms
);
2771 bBox
= ir
.ePBC
!= epbcNONE
;
2772 grid
= init_grid(bBox
, box
, (rcut
> r2cut
) ? rcut
: r2cut
, ngrid
);
2773 nabin
= static_cast<int>(acut
/abin
);
2774 nrbin
= static_cast<int>(rcut
/rbin
);
2775 snew(adist
, nabin
+1);
2776 snew(rdist
, nrbin
+1);
2779 #define __ADIST adist
2780 #define __RDIST rdist
2782 #else /* GMX_OPENMP ================================================== \
2783 * Set up the OpenMP stuff, |
2784 * like the number of threads and such |
2785 * Also start the parallel loop. |
2787 #define __ADIST p_adist[threadNr]
2788 #define __RDIST p_rdist[threadNr]
2789 #define __HBDATA p_hb[threadNr]
2793 bParallel
= !bSelected
;
2797 actual_nThreads
= std::min((nThreads
<= 0) ? INT_MAX
: nThreads
, gmx_omp_get_max_threads());
2799 gmx_omp_set_num_threads(actual_nThreads
);
2800 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads
);
2805 actual_nThreads
= 1;
2808 snew(p_hb
, actual_nThreads
);
2809 snew(p_adist
, actual_nThreads
);
2810 snew(p_rdist
, actual_nThreads
);
2811 for (i
= 0; i
< actual_nThreads
; i
++)
2814 snew(p_adist
[i
], nabin
+1);
2815 snew(p_rdist
[i
], nrbin
+1);
2817 p_hb
[i
]->max_frames
= 0;
2818 p_hb
[i
]->nhb
= NULL
;
2819 p_hb
[i
]->ndist
= NULL
;
2820 p_hb
[i
]->n_bound
= NULL
;
2821 p_hb
[i
]->time
= NULL
;
2822 p_hb
[i
]->nhx
= NULL
;
2824 p_hb
[i
]->bHBmap
= hb
->bHBmap
;
2825 p_hb
[i
]->bDAnr
= hb
->bDAnr
;
2826 p_hb
[i
]->wordlen
= hb
->wordlen
;
2827 p_hb
[i
]->nframes
= hb
->nframes
;
2828 p_hb
[i
]->maxhydro
= hb
->maxhydro
;
2829 p_hb
[i
]->danr
= hb
->danr
;
2832 p_hb
[i
]->hbmap
= hb
->hbmap
;
2833 p_hb
[i
]->time
= hb
->time
; /* This may need re-syncing at every frame. */
2836 p_hb
[i
]->nrdist
= 0;
2840 /* Make a thread pool here,
2841 * instead of forking anew at every frame. */
2843 #pragma omp parallel \
2845 private(j, h, ii, hh, \
2846 xi, yi, zi, xj, yj, zj, threadNr, \
2847 dist, ang, icell, jcell, \
2848 grp, ogrp, ai, aj, xjj, yjj, zjj, \
2851 bEdge_xjj, bEdge_yjj) \
2853 { /* Start of parallel region */
2854 #if !defined __clang_analyzer__ // clang complains about unused value.
2855 threadNr
= gmx_omp_get_thread_num();
2861 bTric
= bBox
&& TRICLINIC(box
);
2867 sync_hbdata(p_hb
[threadNr
], nframes
);
2869 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2875 build_grid(hb
, x
, x
[shatom
], bBox
, box
, hbox
, (rcut
> r2cut
) ? rcut
: r2cut
,
2876 rshell
, ngrid
, grid
);
2877 reset_nhbonds(&(hb
->d
));
2879 if (debug
&& bDebug
)
2881 dump_grid(debug
, ngrid
, grid
);
2884 add_frames(hb
, nframes
);
2885 init_hbframe(hb
, nframes
, output_env_conv_time(oenv
, t
));
2889 count_da_grid(ngrid
, grid
, hb
->danr
[nframes
]);
2892 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2897 p_hb
[threadNr
]->time
= hb
->time
; /* This pointer may have changed. */
2907 /* Do not parallelize this just yet. */
2909 for (ii
= 0; (ii
< nsel
); ii
++)
2911 int dd
= index
[0][i
];
2912 int aa
= index
[0][i
+2];
2913 /* int */ hh
= index
[0][i
+1];
2914 ihb
= is_hbond(hb
, ii
, ii
, dd
, aa
, rcut
, r2cut
, ccut
, x
, bBox
, box
,
2915 hbox
, &dist
, &ang
, bDA
, &h
, bContact
, bMerge
);
2919 /* add to index if not already there */
2921 add_hbond(hb
, dd
, aa
, hh
, ii
, ii
, nframes
, bMerge
, ihb
, bContact
);
2925 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2927 } /* if (bSelected) */
2930 /* The outer grid loop will have to do for now. */
2931 #pragma omp for schedule(dynamic)
2932 for (xi
= 0; xi
< ngrid
[XX
]; xi
++)
2936 for (yi
= 0; (yi
< ngrid
[YY
]); yi
++)
2938 for (zi
= 0; (zi
< ngrid
[ZZ
]); zi
++)
2941 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
2942 for (grp
= gr0
; (grp
<= (bTwo
? gr1
: gr0
)); grp
++)
2944 icell
= &(grid
[zi
][yi
][xi
].d
[grp
]);
2955 /* loop over all hydrogen atoms from group (grp)
2956 * in this gridcell (icell)
2958 for (ai
= 0; (ai
< icell
->nr
); ai
++)
2960 i
= icell
->atoms
[ai
];
2962 /* loop over all adjacent gridcells (xj,yj,zj) */
2963 for (zjj
= grid_loop_begin(ngrid
[ZZ
], zi
, bTric
, FALSE
);
2964 zjj
<= grid_loop_end(ngrid
[ZZ
], zi
, bTric
, FALSE
);
2967 zj
= grid_mod(zjj
, ngrid
[ZZ
]);
2968 bEdge_yjj
= (zj
== 0) || (zj
== ngrid
[ZZ
] - 1);
2969 for (yjj
= grid_loop_begin(ngrid
[YY
], yi
, bTric
, bEdge_yjj
);
2970 yjj
<= grid_loop_end(ngrid
[YY
], yi
, bTric
, bEdge_yjj
);
2973 yj
= grid_mod(yjj
, ngrid
[YY
]);
2975 (yj
== 0) || (yj
== ngrid
[YY
] - 1) ||
2976 (zj
== 0) || (zj
== ngrid
[ZZ
] - 1);
2977 for (xjj
= grid_loop_begin(ngrid
[XX
], xi
, bTric
, bEdge_xjj
);
2978 xjj
<= grid_loop_end(ngrid
[XX
], xi
, bTric
, bEdge_xjj
);
2981 xj
= grid_mod(xjj
, ngrid
[XX
]);
2982 jcell
= &(grid
[zj
][yj
][xj
].a
[ogrp
]);
2983 /* loop over acceptor atoms from other group (ogrp)
2984 * in this adjacent gridcell (jcell)
2986 for (aj
= 0; (aj
< jcell
->nr
); aj
++)
2988 j
= jcell
->atoms
[aj
];
2990 /* check if this once was a h-bond */
2991 ihb
= is_hbond(__HBDATA
, grp
, ogrp
, i
, j
, rcut
, r2cut
, ccut
, x
, bBox
, box
,
2992 hbox
, &dist
, &ang
, bDA
, &h
, bContact
, bMerge
);
2996 /* add to index if not already there */
2998 add_hbond(__HBDATA
, i
, j
, h
, grp
, ogrp
, nframes
, bMerge
, ihb
, bContact
);
3000 /* make angle and distance distributions */
3001 if (ihb
== hbHB
&& !bContact
)
3005 gmx_fatal(FARGS
, "distance is higher than what is allowed for an hbond: %f", dist
);
3008 __ADIST
[static_cast<int>( ang
/abin
)]++;
3009 __RDIST
[static_cast<int>(dist
/rbin
)]++;
3012 if (donor_index(&hb
->d
, grp
, i
) == NOTSET
)
3014 gmx_fatal(FARGS
, "Invalid donor %d", i
);
3016 if (acceptor_index(&hb
->a
, ogrp
, j
) == NOTSET
)
3018 gmx_fatal(FARGS
, "Invalid acceptor %d", j
);
3020 resdist
= std::abs(top
.atoms
.atom
[i
].resind
-top
.atoms
.atom
[j
].resind
);
3021 if (resdist
>= max_hx
)
3025 __HBDATA
->nhx
[nframes
][resdist
]++;
3036 } /* for xi,yi,zi */
3039 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3041 } /* if (bSelected) {...} else */
3044 /* Better wait for all threads to finnish using x[] before updating it. */
3047 #pragma omp critical
3051 /* Sum up histograms and counts from p_hb[] into hb */
3054 hb
->nhb
[k
] += p_hb
[threadNr
]->nhb
[k
];
3055 hb
->ndist
[k
] += p_hb
[threadNr
]->ndist
[k
];
3056 for (j
= 0; j
< max_hx
; j
++)
3058 hb
->nhx
[k
][j
] += p_hb
[threadNr
]->nhx
[k
][j
];
3062 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3065 /* Here are a handful of single constructs
3066 * to share the workload a bit. The most
3067 * important one is of course the last one,
3068 * where there's a potential bottleneck in form
3075 analyse_donor_properties(donor_properties
, hb
, k
, t
);
3077 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3086 do_nhb_dist(fpnhb
, hb
, t
);
3089 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3096 trrStatus
= (read_next_x(oenv
, status
, &t
, x
, box
));
3099 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3108 #pragma omp critical
3110 hb
->nrhb
+= p_hb
[threadNr
]->nrhb
;
3111 hb
->nrdist
+= p_hb
[threadNr
]->nrdist
;
3114 /* Free parallel datastructures */
3115 sfree(p_hb
[threadNr
]->nhb
);
3116 sfree(p_hb
[threadNr
]->ndist
);
3117 sfree(p_hb
[threadNr
]->nhx
);
3120 for (i
= 0; i
< nabin
; i
++)
3124 for (j
= 0; j
< actual_nThreads
; j
++)
3127 adist
[i
] += p_adist
[j
][i
];
3130 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3133 for (i
= 0; i
<= nrbin
; i
++)
3137 for (j
= 0; j
< actual_nThreads
; j
++)
3139 rdist
[i
] += p_rdist
[j
][i
];
3142 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
3145 sfree(p_adist
[threadNr
]);
3146 sfree(p_rdist
[threadNr
]);
3148 } /* End of parallel region */
3155 if (nframes
< 2 && (opt2bSet("-ac", NFILE
, fnm
) || opt2bSet("-life", NFILE
, fnm
)))
3157 gmx_fatal(FARGS
, "Cannot calculate autocorrelation of life times with less than two frames");
3160 free_grid(ngrid
, &grid
);
3164 if (donor_properties
)
3166 xvgrclose(donor_properties
);
3174 /* Compute maximum possible number of different hbonds */
3181 max_nhb
= 0.5*(hb
->d
.nrd
*hb
->a
.nra
);
3188 printf("No %s found!!\n", bContact
? "contacts" : "hydrogen bonds");
3192 printf("Found %d different %s in trajectory\n"
3193 "Found %d different atom-pairs within %s distance\n",
3194 hb
->nrhb
, bContact
? "contacts" : "hydrogen bonds",
3195 hb
->nrdist
, (r2cut
> 0) ? "second cut-off" : "hydrogen bonding");
3199 merge_hb(hb
, bTwo
, bContact
);
3202 if (opt2bSet("-hbn", NFILE
, fnm
))
3204 dump_hbmap(hb
, NFILE
, fnm
, bTwo
, bContact
, isize
, index
, grpnames
, &top
.atoms
);
3207 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3208 * to make the -hbn and -hmb output match eachother.
3209 * - Erik Marklund, May 30, 2006 */
3212 /* Print out number of hbonds and distances */
3215 fp
= xvgropen(opt2fn("-num", NFILE
, fnm
), bContact
? "Contacts" :
3216 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv
), "Number", oenv
);
3218 snew(leg
[0], STRLEN
);
3219 snew(leg
[1], STRLEN
);
3220 sprintf(leg
[0], "%s", bContact
? "Contacts" : "Hydrogen bonds");
3221 sprintf(leg
[1], "Pairs within %g nm", (r2cut
> 0) ? r2cut
: rcut
);
3222 xvgr_legend(fp
, 2, (const char**)leg
, oenv
);
3226 for (i
= 0; (i
< nframes
); i
++)
3228 fprintf(fp
, "%10g %10d %10d\n", hb
->time
[i
], hb
->nhb
[i
], hb
->ndist
[i
]);
3229 aver_nhb
+= hb
->nhb
[i
];
3230 aver_dist
+= hb
->ndist
[i
];
3233 aver_nhb
/= nframes
;
3234 aver_dist
/= nframes
;
3235 /* Print HB distance distribution */
3236 if (opt2bSet("-dist", NFILE
, fnm
))
3241 for (i
= 0; i
< nrbin
; i
++)
3246 fp
= xvgropen(opt2fn("-dist", NFILE
, fnm
),
3247 "Hydrogen Bond Distribution",
3249 "Donor - Acceptor Distance (nm)" :
3250 "Hydrogen - Acceptor Distance (nm)", "", oenv
);
3251 for (i
= 0; i
< nrbin
; i
++)
3253 fprintf(fp
, "%10g %10g\n", (i
+0.5)*rbin
, rdist
[i
]/(rbin
*sum
));
3258 /* Print HB angle distribution */
3259 if (opt2bSet("-ang", NFILE
, fnm
))
3264 for (i
= 0; i
< nabin
; i
++)
3269 fp
= xvgropen(opt2fn("-ang", NFILE
, fnm
),
3270 "Hydrogen Bond Distribution",
3271 "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv
);
3272 for (i
= 0; i
< nabin
; i
++)
3274 fprintf(fp
, "%10g %10g\n", (i
+0.5)*abin
, adist
[i
]/(abin
*sum
));
3279 /* Print HB in alpha-helix */
3280 if (opt2bSet("-hx", NFILE
, fnm
))
3282 fp
= xvgropen(opt2fn("-hx", NFILE
, fnm
),
3283 "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv
), "Count", oenv
);
3284 xvgr_legend(fp
, NRHXTYPES
, hxtypenames
, oenv
);
3285 for (i
= 0; i
< nframes
; i
++)
3287 fprintf(fp
, "%10g", hb
->time
[i
]);
3288 for (j
= 0; j
< max_hx
; j
++)
3290 fprintf(fp
, " %6d", hb
->nhx
[i
][j
]);
3297 printf("Average number of %s per timeframe %.3f out of %g possible\n",
3298 bContact
? "contacts" : "hbonds",
3299 bContact
? aver_dist
: aver_nhb
, max_nhb
);
3301 /* Do Autocorrelation etc. */
3305 Added support for -contact in ac and hbm calculations below.
3306 - Erik Marklund, May 29, 2006
3308 if (opt2bSet("-ac", NFILE
, fnm
) || opt2bSet("-life", NFILE
, fnm
))
3310 please_cite(stdout
, "Spoel2006b");
3312 if (opt2bSet("-ac", NFILE
, fnm
))
3314 do_hbac(opt2fn("-ac", NFILE
, fnm
), hb
, nDump
,
3315 bMerge
, bContact
, fit_start
, temp
, r2cut
> 0, oenv
,
3318 if (opt2bSet("-life", NFILE
, fnm
))
3320 do_hblife(opt2fn("-life", NFILE
, fnm
), hb
, bMerge
, bContact
, oenv
);
3322 if (opt2bSet("-hbm", NFILE
, fnm
))
3325 int id
, ia
, hh
, x
, y
;
3326 mat
.flags
= mat
.y0
= 0;
3328 if ((nframes
> 0) && (hb
->nrhb
> 0))
3333 snew(mat
.matrix
, mat
.nx
);
3334 for (x
= 0; (x
< mat
.nx
); x
++)
3336 snew(mat
.matrix
[x
], mat
.ny
);
3339 for (id
= 0; (id
< hb
->d
.nrd
); id
++)
3341 for (ia
= 0; (ia
< hb
->a
.nra
); ia
++)
3343 for (hh
= 0; (hh
< hb
->maxhydro
); hh
++)
3345 if (hb
->hbmap
[id
][ia
])
3347 if (ISHB(hb
->hbmap
[id
][ia
]->history
[hh
]))
3349 for (x
= 0; (x
<= hb
->hbmap
[id
][ia
]->nframes
); x
++)
3351 int nn0
= hb
->hbmap
[id
][ia
]->n0
;
3352 range_check(y
, 0, mat
.ny
);
3353 mat
.matrix
[x
+nn0
][y
] = is_hb(hb
->hbmap
[id
][ia
]->h
[hh
], x
);
3361 mat
.axis_x
= hb
->time
;
3362 snew(mat
.axis_y
, mat
.ny
);
3363 for (j
= 0; j
< mat
.ny
; j
++)
3367 sprintf(mat
.title
, bContact
? "Contact Existence Map" :
3368 "Hydrogen Bond Existence Map");
3369 sprintf(mat
.legend
, bContact
? "Contacts" : "Hydrogen Bonds");
3370 sprintf(mat
.label_x
, "%s", output_env_get_xvgr_tlabel(oenv
));
3371 sprintf(mat
.label_y
, bContact
? "Contact Index" : "Hydrogen Bond Index");
3372 mat
.bDiscrete
= TRUE
;
3374 snew(mat
.map
, mat
.nmap
);
3375 for (i
= 0; i
< mat
.nmap
; i
++)
3377 mat
.map
[i
].code
.c1
= hbmap
[i
];
3378 mat
.map
[i
].desc
= hbdesc
[i
];
3379 mat
.map
[i
].rgb
= hbrgb
[i
];
3381 fp
= opt2FILE("-hbm", NFILE
, fnm
, "w");
3382 write_xpm_m(fp
, mat
);
3384 for (x
= 0; x
< mat
.nx
; x
++)
3386 sfree(mat
.matrix
[x
]);
3394 fprintf(stderr
, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
3405 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
3407 fp
= xvgropen(opt2fn("-dan", NFILE
, fnm
),
3408 "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv
), "Count", oenv
);
3409 nleg
= (bTwo
? 2 : 1)*2;
3410 snew(legnames
, nleg
);
3412 for (j
= 0; j
< grNR
; j
++)
3414 if (USE_THIS_GROUP(j
) )
3416 sprintf(buf
, "Donors %s", grpnames
[j
]);
3417 legnames
[i
++] = gmx_strdup(buf
);
3418 sprintf(buf
, "Acceptors %s", grpnames
[j
]);
3419 legnames
[i
++] = gmx_strdup(buf
);
3424 gmx_incons("number of legend entries");
3426 xvgr_legend(fp
, nleg
, (const char**)legnames
, oenv
);
3427 for (i
= 0; i
< nframes
; i
++)
3429 fprintf(fp
, "%10g", hb
->time
[i
]);
3430 for (j
= 0; (j
< grNR
); j
++)
3432 if (USE_THIS_GROUP(j
) )
3434 fprintf(fp
, " %6d", hb
->danr
[i
][j
]);