1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2008, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Groningen Machine for Chemical Simulation
41 /*#define HAVE_NN_LOOPS*/
42 /* Set environment variable CFLAGS = "-fopenmp" when running
43 * configure and define DOUSEOPENMP to make use of parallelized
44 * calculation of autocorrelation function.
45 * It also adds a new option -nthreads which sets the number of threads.
47 /*#define DOUSEOPENMP*/
62 #include "gmx_fatal.h"
75 typedef short int t_E;
78 typedef int t_hx[max_hx];
79 #define NRHXTYPES max_hx
80 const char *hxtypenames[NRHXTYPES]=
81 {"n-n","n-n+1","n-n+2","n-n+3","n-n+4","n-n+5","n-n>6"};
85 #define MASTER_THREAD_ONLY(threadNr) ((threadNr)==0)
87 #define MASTER_THREAD_ONLY(threadNr) ((threadNr)==(threadNr))
90 /* -----------------------------------------*/
92 enum { gr0, gr1, grI, grNR };
93 enum { hbNo, hbDist, hbHB, hbNR, hbR2};
94 enum { noDA, ACC, DON, DA, INGROUP};
95 enum {NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR};
97 static const char *grpnames[grNR] = {"0","1","I" };
99 static gmx_bool bDebug = FALSE;
104 #define HB_YESINS HB_YES|HB_INS
108 #define ISHB(h) (((h) & 2) == 2)
109 #define ISDIST(h) (((h) & 1) == 1)
110 #define ISDIST2(h) (((h) & 4) == 4)
111 #define ISACC(h) (((h) & 1) == 1)
112 #define ISDON(h) (((h) & 2) == 2)
113 #define ISINGRP(h) (((h) & 4) == 4)
126 typedef int t_icell[grNR];
127 typedef atom_id h_id[MAXHYDRO];
130 int history[MAXHYDRO];
131 /* Has this hbond existed ever? If so as hbDist or hbHB or both.
132 * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
134 /* Bitmask array which tells whether a hbond is present
135 * at a given time. Either of these may be NULL
137 int n0; /* First frame a HB was found */
138 int nframes,maxframes; /* Amount of frames in this hbond */
141 /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
142 * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
143 * acceptor distance is less than the user-specified distance (typically
150 atom_id *acc; /* Atom numbers of the acceptors */
151 int *grp; /* Group index */
152 int *aptr; /* Map atom number to acceptor index */
157 int *don; /* Atom numbers of the donors */
158 int *grp; /* Group index */
159 int *dptr; /* Map atom number to donor index */
160 int *nhydro; /* Number of hydrogens for each donor */
161 h_id *hydro; /* The atom numbers of the hydrogens */
162 h_id *nhbonds; /* The number of HBs per H at current */
165 /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/
169 int len; /* The length of frame and p. */
170 int *frame; /* The frames at which transitio*/
175 /* Periodicity history. Used for the reversible geminate recombination. */
176 t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap:
177 * pHist[d][a]. We can safely assume that the same
178 * periodic shift holds for all hydrogens of a da-pair.
180 * Nowadays it only stores TRANSITIONS, and not the shift at every frame.
181 * That saves a LOT of memory, an hopefully kills a mysterious bug where
182 * pHist gets contaminated. */
184 PSTYPE nper; /* The length of p2i */
185 ivec *p2i; /* Maps integer to periodic shift for a pair.*/
186 matrix P; /* Projection matrix to find the box shifts. */
187 int gemtype; /* enumerated type */
192 int *Etot; /* Total energy for each frame */
193 t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */
197 gmx_bool bHBmap,bDAnr,bGem;
199 /* The following arrays are nframes long */
200 int nframes,max_frames,maxhydro;
206 /* These structures are initialized from the topology at start up */
209 /* This holds a matrix with all possible hydrogen bonds */
215 /* For parallelization reasons this will have to be a pointer.
216 * Otherwise discrepancies may arise between the periodicity data
217 * seen by different threads. */
221 static void clearPshift(t_pShift *pShift)
223 if (pShift->len > 0) {
225 sfree(pShift->frame);
230 static void calcBoxProjection(matrix B, matrix P)
232 const int vp[] = {XX,YY,ZZ};
237 for (i=0; i<3; i++){ m = vp[i];
238 for (j=0; j<3; j++){ n = vp[j];
239 U[m][n] = i==j ? 1:0;
243 for (i=0; i<3; i++){ m = vp[i];
244 mvmul(M, U[m], P[m]);
249 static void calcBoxDistance(matrix P, rvec d, ivec ibd){
250 /* returns integer distance in box coordinates.
251 * P is the projection matrix from cartesian coordinates
252 * obtained with calcBoxProjection(). */
256 /* extend it by 0.5 in all directions since (int) rounds toward 0.*/
258 bd[i]=bd[i] + (bd[i]<0 ? -0.5 : 0.5);
259 ibd[XX] = (int)bd[XX];
260 ibd[YY] = (int)bd[YY];
261 ibd[ZZ] = (int)bd[ZZ];
264 /* Changed argument 'bMerge' into 'oneHB' below,
265 * since -contact should cause maxhydro to be 1,
267 * - Erik Marklund May 29, 2006
270 static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap) {
271 /* Try to merge hbonds on the fly. That means that if the
272 * acceptor and donor are mergable, then:
273 * 1) store the hb-info so that acceptor id > donor id,
274 * 2) add the periodic shift in pairs, so that [-x,-y,-z] is
275 * stored in per.p2i[] whenever acceptor id < donor id.
276 * Note that [0,0,0] should already be the first element of per.p2i
277 * by the time this function is called. */
279 /* daSwap is TRUE if the donor and acceptor were swapped.
280 * If so, then the negative vector should be used. */
283 if (per->p2i == NULL || per->nper == 0)
284 gmx_fatal(FARGS, "'per' not initialized properly.");
285 for (i=0; i<per->nper; i++) {
286 if (r[XX] == per->p2i[i][XX] &&
287 r[YY] == per->p2i[i][YY] &&
288 r[ZZ] == per->p2i[i][ZZ])
291 /* Not found apparently. Add it to the list! */
292 /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
294 /* Unfortunately this needs to be critical it seems. */
300 fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
304 srenew(per->p2i, per->nper+2);
305 copy_ivec(r, per->p2i[per->nper]);
308 /* Add the mirror too. It's rather likely that it'll be needed. */
309 per->p2i[per->nper][XX] = -r[XX];
310 per->p2i[per->nper][YY] = -r[YY];
311 per->p2i[per->nper][ZZ] = -r[ZZ];
314 return per->nper - 1 - (daSwap ? 0:1);
317 static t_hbdata *mk_hbdata(gmx_bool bHBmap,gmx_bool bDAnr,gmx_bool oneHB, gmx_bool bGem, int gemmode)
322 hb->wordlen = 8*sizeof(unsigned int);
329 hb->maxhydro = MAXHYDRO;
331 hb->per->gemtype = bGem ? gemmode : 0;
336 static void mk_hbmap(t_hbdata *hb,gmx_bool bTwo)
340 snew(hb->hbmap,hb->d.nrd);
341 for(i=0; (i<hb->d.nrd); i++) {
342 snew(hb->hbmap[i],hb->a.nra);
343 if (hb->hbmap[i] == NULL)
344 gmx_fatal(FARGS,"Could not allocate enough memory for hbmap");
345 for (j=0; (j>hb->a.nra); j++)
346 hb->hbmap[i][j] = NULL;
350 /* Consider redoing pHist so that is only stores transitions between
351 * periodicities and not the periodicity for all frames. This eats heaps of memory. */
352 static void mk_per(t_hbdata *hb)
356 snew(hb->per->pHist, hb->d.nrd);
357 for (i=0; i<hb->d.nrd; i++) {
358 snew(hb->per->pHist[i], hb->a.nra);
359 if (hb->per->pHist[i]==NULL)
360 gmx_fatal(FARGS,"Could not allocate enough memory for per->pHist");
361 for (j=0; j<hb->a.nra; j++) {
362 clearPshift(&(hb->per->pHist[i][j]));
365 /* add the [0,0,0] shift to element 0 of p2i. */
366 snew(hb->per->p2i, 1);
367 clear_ivec(hb->per->p2i[0]);
373 static void mk_hbEmap (t_hbdata *hb, int n0)
378 snew(hb->hbE.E, hb->d.nrd);
379 for (i=0; i<hb->d.nrd; i++)
381 snew(hb->hbE.E[i], hb->a.nra);
382 for (j=0; j<hb->a.nra; j++)
384 snew(hb->hbE.E[i][j], MAXHYDRO);
385 for (k=0; k<MAXHYDRO; k++)
386 hb->hbE.E[i][j][k] = NULL;
392 static void free_hbEmap (t_hbdata *hb)
395 for (i=0; i<hb->d.nrd; i++)
397 for (j=0; j<hb->a.nra; j++)
399 for (k=0; k<MAXHYDRO; k++)
400 sfree(hb->hbE.E[i][j][k]);
401 sfree(hb->hbE.E[i][j]);
409 static void addFramesNN(t_hbdata *hb, int frame)
412 #define DELTAFRAMES_HBE 10
416 if (frame >= hb->hbE.nframes) {
417 nframes = hb->hbE.nframes + DELTAFRAMES_HBE;
418 srenew(hb->hbE.Etot, nframes);
420 for (d=0; d<hb->d.nrd; d++)
421 for (a=0; a<hb->a.nra; a++)
422 for (h=0; h<hb->d.nhydro[d]; h++)
423 srenew(hb->hbE.E[d][a][h], nframes);
425 hb->hbE.nframes += DELTAFRAMES_HBE;
429 static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst,
430 matrix box, rvec hbox, t_donors *donors){
434 * alpha - angle between dipoles
435 * x[] - atomic positions
436 * EEst - the type of energy estimate (see enum in hbplugin.h)
437 * box - the box vectors \
438 * hbox - half box lengths _These two are only needed for the pbc correction
443 rvec dipole[2], xmol[3], xmean[2];
448 /* Self-interaction */
454 /* This is a simple binary existence function that sets E=1 whenever
455 * the distance between the oxygens is equal too or less than 0.35 nm.
457 rvec_sub(x[d], x[a], dist);
458 pbc_correct_gem(dist, box, hbox);
459 if (norm(dist) <= 0.35)
466 /* Negative potential energy of a dipole.
467 * E = -cos(alpha) * 1/r^3 */
469 copy_rvec(x[d], xmol[0]); /* donor */
470 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
471 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
473 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
474 rvec_inc(xmean[0], xmol[1]);
475 rvec_inc(xmean[0], xmol[2]);
477 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
479 /* Assumes that all acceptors are also donors. */
480 copy_rvec(x[a], xmol[0]); /* acceptor */
481 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
482 copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */
485 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
486 rvec_inc(xmean[1], xmol[1]);
487 rvec_inc(xmean[1], xmol[2]);
489 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
491 rvec_sub(xmean[0], xmean[1], dist);
492 pbc_correct_gem(dist, box, hbox);
495 realE = pow(r, -3.0);
496 E = (t_E)(SCALEFACTOR_E * realE);
500 /* Negative potential energy of a (unpolarizable) dipole.
501 * E = -cos(alpha) * 1/r^3 */
502 clear_rvec(dipole[1]);
503 clear_rvec(dipole[0]);
505 copy_rvec(x[d], xmol[0]); /* donor */
506 copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */
507 copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */
509 rvec_inc(dipole[0], xmol[1]);
510 rvec_inc(dipole[0], xmol[2]);
513 rvec_dec(dipole[0], xmol[0]);
515 svmul(15.9994*(1/1.008), xmol[0], xmean[0]);
516 rvec_inc(xmean[0], xmol[1]);
517 rvec_inc(xmean[0], xmol[2]);
519 xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008;
521 /* Assumes that all acceptors are also donors. */
522 copy_rvec(x[a], xmol[0]); /* acceptor */
523 copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */
524 copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */
527 rvec_inc(dipole[1], xmol[1]);
528 rvec_inc(dipole[1], xmol[2]);
531 rvec_dec(dipole[1], xmol[0]);
533 svmul(15.9994*(1/1.008), xmol[0], xmean[1]);
534 rvec_inc(xmean[1], xmol[1]);
535 rvec_inc(xmean[1], xmol[2]);
537 xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008;
539 rvec_sub(xmean[0], xmean[1], dist);
540 pbc_correct_gem(dist, box, hbox);
543 double cosalpha = cos_angle(dipole[0],dipole[1]);
544 realE = cosalpha * pow(r, -3.0);
545 E = (t_E)(SCALEFACTOR_E * realE);
549 printf("Can't do that type of energy estimate: %i\n.", EEst);
556 static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame){
557 /* hb - hbond data structure
561 E - estimate of the energy
562 frame - the current frame.
565 /* Store the estimated energy */
569 hb->hbE.E[d][a][h][frame] = E;
574 hb->hbE.Etot[frame] += E;
577 #endif /* HAVE_NN_LOOPS */
580 /* Finds -v[] in the periodicity index */
581 static int findMirror(PSTYPE p, ivec v[], PSTYPE nper)
585 for (i=0; i<nper; i++){
586 if (v[i][XX] == -(v[p][XX]) &&
587 v[i][YY] == -(v[p][YY]) &&
588 v[i][ZZ] == -(v[p][ZZ]))
591 printf("Couldn't find mirror of [%i, %i, %i], index \n",
599 static void add_frames(t_hbdata *hb,int nframes)
603 if (nframes >= hb->max_frames) {
604 hb->max_frames += 4096;
605 srenew(hb->time,hb->max_frames);
606 srenew(hb->nhb,hb->max_frames);
607 srenew(hb->ndist,hb->max_frames);
608 srenew(hb->n_bound,hb->max_frames);
609 srenew(hb->nhx,hb->max_frames);
611 srenew(hb->danr,hb->max_frames);
616 #define OFFSET(frame) (frame / 32)
617 #define MASK(frame) (1 << (frame % 32))
619 static void _set_hb(unsigned int hbexist[],unsigned int frame,gmx_bool bValue)
622 hbexist[OFFSET(frame)] |= MASK(frame);
624 hbexist[OFFSET(frame)] &= ~MASK(frame);
627 static gmx_bool is_hb(unsigned int hbexist[],int frame)
629 return ((hbexist[OFFSET(frame)] & MASK(frame)) != 0) ? 1 : 0;
632 static void set_hb(t_hbdata *hb,int id,int ih, int ia,int frame,int ihb)
634 unsigned int *ghptr=NULL;
637 ghptr = hb->hbmap[id][ia]->h[ih];
638 else if (ihb == hbDist)
639 ghptr = hb->hbmap[id][ia]->g[ih];
641 gmx_fatal(FARGS,"Incomprehensible iValue %d in set_hb",ihb);
643 _set_hb(ghptr,frame-hb->hbmap[id][ia]->n0,TRUE);
646 static void addPshift(t_pShift *pHist, PSTYPE p, int frame)
648 if (pHist->len == 0) {
649 snew(pHist->frame, 1);
652 pHist->frame[0] = frame;
656 if (pHist->p[pHist->len-1] != p) {
658 srenew(pHist->frame, pHist->len);
659 srenew(pHist->p, pHist->len);
660 pHist->frame[pHist->len-1] = frame;
661 pHist->p[pHist->len-1] = p;
662 } /* Otherwise, there is no transition. */
666 static PSTYPE getPshift(t_pShift pHist, int frame)
671 || (pHist.len > 0 && pHist.frame[0]>frame))
674 for (i=0; i<pHist.len; i++)
683 /* It seems that frame is after the last periodic transition. Return the last periodicity. */
684 return pHist.p[pHist.len-1];
687 static void add_ff(t_hbdata *hbd,int id,int h,int ia,int frame,int ihb, PSTYPE p)
690 t_hbond *hb = hbd->hbmap[id][ia];
691 int maxhydro = min(hbd->maxhydro,hbd->d.nhydro[id]);
692 int wlen = hbd->wordlen;
694 gmx_bool bGem = hbd->bGem;
698 hb->maxframes = delta;
699 for(i=0; (i<maxhydro); i++) {
700 snew(hb->h[i],hb->maxframes/wlen);
701 snew(hb->g[i],hb->maxframes/wlen);
704 hb->nframes = frame-hb->n0;
705 /* We need a while loop here because hbonds may be returning
708 while (hb->nframes >= hb->maxframes) {
709 n = hb->maxframes + delta;
710 for(i=0; (i<maxhydro); i++) {
711 srenew(hb->h[i],n/wlen);
712 srenew(hb->g[i],n/wlen);
713 for(j=hb->maxframes/wlen; (j<n/wlen); j++) {
723 set_hb(hbd,id,h,ia,frame,ihb);
725 if (p>=hbd->per->nper)
726 gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper);
728 addPshift(&(hbd->per->pHist[id][ia]), p, frame);
734 static void inc_nhbonds(t_donors *ddd,int d, int h)
737 int dptr = ddd->dptr[d];
739 for(j=0; (j<ddd->nhydro[dptr]); j++)
740 if (ddd->hydro[dptr][j] == h) {
741 ddd->nhbonds[dptr][j]++;
744 if (j == ddd->nhydro[dptr])
745 gmx_fatal(FARGS,"No such hydrogen %d on donor %d\n",h+1,d+1);
748 static int _acceptor_index(t_acceptors *a,int grp,atom_id i,
749 const char *file,int line)
753 if (a->grp[ai] != grp) {
755 fprintf(debug,"Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
756 ai,a->grp[ai],grp,file,line);
762 #define acceptor_index(a,grp,i) _acceptor_index(a,grp,i,__FILE__,__LINE__)
764 static int _donor_index(t_donors *d,int grp,atom_id i,const char *file,int line)
771 if (d->grp[di] != grp) {
773 fprintf(debug,"Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
774 di,d->grp[di],grp,file,line);
780 #define donor_index(d,grp,i) _donor_index(d,grp,i,__FILE__,__LINE__)
782 static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
784 /* g_hbond doesn't allow overlapping groups */
788 donor_index(&hb->d,grpd,a) != NOTSET
789 && acceptor_index(&hb->a,grpa,d) != NOTSET;
793 static void add_hbond(t_hbdata *hb,int d,int a,int h,int grpd,int grpa,
794 int frame,gmx_bool bMerge,int ihb,gmx_bool bContact, PSTYPE p)
797 gmx_bool daSwap = FALSE;
799 if ((id = hb->d.dptr[d]) == NOTSET)
800 gmx_fatal(FARGS,"No donor atom %d",d+1);
801 else if (grpd != hb->d.grp[id])
802 gmx_fatal(FARGS,"Inconsistent donor groups, %d iso %d, atom %d",
803 grpd,hb->d.grp[id],d+1);
804 if ((ia = hb->a.aptr[a]) == NOTSET)
805 gmx_fatal(FARGS,"No acceptor atom %d",a+1);
806 else if (grpa != hb->a.grp[ia])
807 gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
808 grpa,hb->a.grp[ia],a+1);
811 if ((daSwap = isInterchangable(hb, d, a, grpd, grpa) || bContact) && d>a)
812 /* Then swap identity so that the id of d is lower then that of a.
814 * This should really be redundant by now, as is_hbond() now ought to return
815 * hbNo in the cases where this conditional is TRUE. */
821 /* Now repeat donor/acc check. */
822 if ((id = hb->d.dptr[d]) == NOTSET)
823 gmx_fatal(FARGS,"No donor atom %d",d+1);
824 else if (grpd != hb->d.grp[id])
825 gmx_fatal(FARGS,"Inconsistent donor groups, %d iso %d, atom %d",
826 grpd,hb->d.grp[id],d+1);
827 if ((ia = hb->a.aptr[a]) == NOTSET)
828 gmx_fatal(FARGS,"No acceptor atom %d",a+1);
829 else if (grpa != hb->a.grp[ia])
830 gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
831 grpa,hb->a.grp[ia],a+1);
835 /* Loop over hydrogens to find which hydrogen is in this particular HB */
836 if ((ihb == hbHB) && !bMerge && !bContact) {
837 for(k=0; (k<hb->d.nhydro[id]); k++)
838 if (hb->d.hydro[id][k] == h)
840 if (k == hb->d.nhydro[id])
841 gmx_fatal(FARGS,"Donor %d does not have hydrogen %d (a = %d)",
848 if (hb->hbmap[id][ia] == NULL) {
849 snew(hb->hbmap[id][ia],1);
850 snew(hb->hbmap[id][ia]->h,hb->maxhydro);
851 snew(hb->hbmap[id][ia]->g,hb->maxhydro);
853 add_ff(hb,id,k,ia,frame,ihb,p);
856 /* Strange construction with frame >=0 is a relic from old code
857 * for selected hbond analysis. It may be necessary again if that
858 * is made to work again.
861 hh = hb->hbmap[id][ia]->history[k];
865 hb->hbmap[id][ia]->history[k] = hh | 2;
874 hb->hbmap[id][ia]->history[k] = hh | 1;
891 if (bMerge && daSwap)
892 h = hb->d.hydro[id][0];
893 /* Increment number if HBonds per H */
894 if (ihb == hbHB && !bContact)
895 inc_nhbonds(&(hb->d),d,h);
898 /* Now a redundant function. It might find use at some point though. */
899 static gmx_bool in_list(atom_id selection,int isize,atom_id *index)
905 for(i=0; (i<isize) && !bFound; i++)
906 if(selection == index[i])
912 static char *mkatomname(t_atoms *atoms,int i)
917 rnr = atoms->atom[i].resind;
918 sprintf(buf,"%4s%d%-4s",
919 *atoms->resinfo[rnr].name,atoms->resinfo[rnr].nr,*atoms->atomname[i]);
924 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms){
925 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
928 for (i=0; i<isize; i++){
929 if (index[i] >= natoms)
930 gmx_fatal(FARGS,"Atom has index %d larger than number of atoms %d.",index[i],natoms);
931 datable[index[i]] |= INGROUP;
935 static void clear_datable_grp(unsigned char *datable, int size){
936 /* Clears group information from the table */
938 const char mask = !(char)INGROUP;
944 static void add_acc(t_acceptors *a,int ia,int grp)
946 if (a->nra >= a->max_nra) {
948 srenew(a->acc,a->max_nra);
949 srenew(a->grp,a->max_nra);
951 a->grp[a->nra] = grp;
952 a->acc[a->nra++] = ia;
955 static void search_acceptors(t_topology *top,int isize,
956 atom_id *index,t_acceptors *a,int grp,
958 gmx_bool bContact,gmx_bool bDoIt, unsigned char *datable)
963 for (i=0; (i<isize); i++) {
966 (((*top->atoms.atomname[n])[0] == 'O') ||
967 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
968 ISINGRP(datable[n])) {
969 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
974 snew(a->aptr,top->atoms.nr);
975 for(i=0; (i<top->atoms.nr); i++)
977 for(i=0; (i<a->nra); i++)
978 a->aptr[a->acc[i]] = i;
981 static void add_h2d(int id,int ih,t_donors *ddd)
985 for(i=0; (i<ddd->nhydro[id]); i++)
986 if (ddd->hydro[id][i] == ih) {
987 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
991 if (i == ddd->nhydro[id]) {
992 if (ddd->nhydro[id] >= MAXHYDRO)
993 gmx_fatal(FARGS,"Donor %d has more than %d hydrogens!",
994 ddd->don[id],MAXHYDRO);
995 ddd->hydro[id][i] = ih;
1000 static void add_dh(t_donors *ddd,int id,int ih,int grp, unsigned char *datable)
1004 if (ISDON(datable[id]) || !datable) {
1005 if (ddd->dptr[id] == NOTSET) { /* New donor */
1011 if (i == ddd->nrd) {
1012 if (ddd->nrd >= ddd->max_nrd) {
1013 ddd->max_nrd += 128;
1014 srenew(ddd->don,ddd->max_nrd);
1015 srenew(ddd->nhydro,ddd->max_nrd);
1016 srenew(ddd->hydro,ddd->max_nrd);
1017 srenew(ddd->nhbonds,ddd->max_nrd);
1018 srenew(ddd->grp,ddd->max_nrd);
1020 ddd->don[ddd->nrd] = id;
1021 ddd->nhydro[ddd->nrd] = 0;
1022 ddd->grp[ddd->nrd] = grp;
1029 printf("Warning: Atom %d is not in the d/a-table!\n", id);
1032 static void search_donors(t_topology *top, int isize, atom_id *index,
1033 t_donors *ddd,int grp,gmx_bool bContact,gmx_bool bDoIt,
1034 unsigned char *datable)
1037 t_functype func_type;
1038 t_ilist *interaction;
1043 snew(ddd->dptr,top->atoms.nr);
1044 for(i=0; (i<top->atoms.nr); i++)
1045 ddd->dptr[i] = NOTSET;
1050 for(i=0; (i<isize); i++) {
1051 datable[index[i]] |= DON;
1052 add_dh(ddd,index[i],-1,grp,datable);
1056 for(func_type=0; (func_type < F_NRE); func_type++) {
1057 interaction=&(top->idef.il[func_type]);
1058 if (func_type == F_POSRES)
1059 { /* The ilist looks strange for posre. Bug in grompp?
1060 * We don't need posre interactions for hbonds anyway.*/
1063 for(i=0; i < interaction->nr;
1064 i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
1066 if (func_type != top->idef.functype[interaction->iatoms[i]]) {
1067 fprintf(stderr,"Error in func_type %s",
1068 interaction_function[func_type].longname);
1072 /* check out this functype */
1073 if (func_type == F_SETTLE) {
1074 nr1=interaction->iatoms[i+1];
1076 if (ISINGRP(datable[nr1])) {
1077 if (ISINGRP(datable[nr1+1])) {
1078 datable[nr1] |= DON;
1079 add_dh(ddd,nr1,nr1+1,grp,datable);
1081 if (ISINGRP(datable[nr1+2])) {
1082 datable[nr1] |= DON;
1083 add_dh(ddd,nr1,nr1+2,grp,datable);
1087 else if (IS_CHEMBOND(func_type)) {
1088 for (j=0; j<2; j++) {
1089 nr1=interaction->iatoms[i+1+j];
1090 nr2=interaction->iatoms[i+2-j];
1091 if ((*top->atoms.atomname[nr1][0] == 'H') &&
1092 ((*top->atoms.atomname[nr2][0] == 'O') ||
1093 (*top->atoms.atomname[nr2][0] == 'N')) &&
1094 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
1095 datable[nr2] |= DON;
1096 add_dh(ddd,nr2,nr1,grp,datable);
1103 for(func_type=0; func_type < F_NRE; func_type++) {
1104 interaction=&top->idef.il[func_type];
1105 for(i=0; i < interaction->nr;
1106 i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
1108 if (func_type != top->idef.functype[interaction->iatoms[i]])
1109 gmx_incons("function type in search_donors");
1111 if ( interaction_function[func_type].flags & IF_VSITE ) {
1112 nr1=interaction->iatoms[i+1];
1113 if ( *top->atoms.atomname[nr1][0] == 'H') {
1116 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1121 if ( !stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1122 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1123 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
1124 datable[nr2] |= DON;
1125 add_dh(ddd,nr2,nr1,grp,datable);
1135 static t_gridcell ***init_grid(gmx_bool bBox,rvec box[],real rcut,ivec ngrid)
1141 for(i=0; i<DIM; i++)
1142 ngrid[i]=(box[i][i]/(1.2*rcut));
1144 if ( !bBox || (ngrid[XX]<3) || (ngrid[YY]<3) || (ngrid[ZZ]<3) )
1145 for(i=0; i<DIM; i++)
1148 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1149 ngrid[XX],ngrid[YY],ngrid[ZZ],rcut);
1150 snew(grid,ngrid[ZZ]);
1151 for (z=0; z<ngrid[ZZ]; z++) {
1152 snew((grid)[z],ngrid[YY]);
1153 for (y=0; y<ngrid[YY]; y++)
1154 snew((grid)[z][y],ngrid[XX]);
1159 static void control_pHist(t_hbdata *hb, int nframes)
1163 for (i=0;i<hb->d.nrd;i++)
1164 for (j=0;j<hb->a.nra;j++)
1165 if (hb->per->pHist[i][j].len != 0)
1166 for (k=hb->hbmap[i][j][0].n0; k<nframes; k++) {
1167 p = getPshift(hb->per->pHist[i][j], k);
1168 if (p>hb->per->nper)
1169 fprintf(stderr, "Weird stuff in pHist[%i][%i].p at frame %i: p=%i\n",
1174 static void reset_nhbonds(t_donors *ddd)
1178 for(i=0; (i<ddd->nrd); i++)
1179 for(j=0; (j<MAXHH); j++)
1180 ddd->nhbonds[i][j] = 0;
1183 void pbc_correct_gem(rvec dx,matrix box,rvec hbox);
1185 static void build_grid(t_hbdata *hb,rvec x[], rvec xshell,
1186 gmx_bool bBox, matrix box, rvec hbox,
1187 real rcut, real rshell,
1188 ivec ngrid, t_gridcell ***grid)
1190 int i,m,gr,xi,yi,zi,nr;
1193 rvec invdelta,dshell,xtemp={0,0,0};
1195 gmx_bool bDoRshell,bInShell,bAcc;
1200 bDoRshell = (rshell > 0);
1201 rshell2 = sqr(rshell);
1204 #define DBB(x) if (debug && bDebug) fprintf(debug,"build_grid, line %d, %s = %d\n",__LINE__,#x,x)
1206 for(m=0; m<DIM; m++) {
1207 hbox[m]=box[m][m]*0.5;
1209 invdelta[m]=ngrid[m]/box[m][m];
1210 if (1/invdelta[m] < rcut)
1211 gmx_fatal(FARGS,"Your computational box has shrunk too much.\n"
1212 "%s can not handle this situation, sorry.\n",
1221 /* resetting atom counts */
1222 for(gr=0; (gr<grNR); gr++) {
1223 for (zi=0; zi<ngrid[ZZ]; zi++)
1224 for (yi=0; yi<ngrid[YY]; yi++)
1225 for (xi=0; xi<ngrid[XX]; xi++) {
1226 grid[zi][yi][xi].d[gr].nr=0;
1227 grid[zi][yi][xi].a[gr].nr=0;
1231 /* put atoms in grid cells */
1232 for(bAcc=FALSE; (bAcc<=TRUE); bAcc++) {
1242 for(i=0; (i<nr); i++) {
1243 /* check if we are inside the shell */
1244 /* if bDoRshell=FALSE then bInShell=TRUE always */
1248 rvec_sub(x[ad[i]],xshell,dshell);
1250 if (FALSE && !hb->bGem) {
1251 for(m=DIM-1; m>=0 && bInShell; m--) {
1252 if ( dshell[m] < -hbox[m] )
1253 rvec_inc(dshell,box[m]);
1254 else if ( dshell[m] >= hbox[m] )
1255 dshell[m] -= 2*hbox[m];
1256 /* if we're outside the cube, we're outside the sphere also! */
1257 if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
1261 gmx_bool bDone = FALSE;
1265 for(m=DIM-1; m>=0 && bInShell; m--) {
1266 if ( dshell[m] < -hbox[m] ) {
1268 rvec_inc(dshell,box[m]);
1270 if ( dshell[m] >= hbox[m] ) {
1272 dshell[m] -= 2*hbox[m];
1276 for(m=DIM-1; m>=0 && bInShell; m--) {
1277 /* if we're outside the cube, we're outside the sphere also! */
1278 if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
1283 /* if we're inside the cube, check if we're inside the sphere */
1285 bInShell = norm2(dshell) < rshell2;
1291 copy_rvec(x[ad[i]], xtemp);
1292 pbc_correct_gem(x[ad[i]], box, hbox);
1294 for(m=DIM-1; m>=0; m--) {
1295 if (TRUE || !hb->bGem){
1296 /* put atom in the box */
1297 while( x[ad[i]][m] < 0 )
1298 rvec_inc(x[ad[i]],box[m]);
1299 while( x[ad[i]][m] >= box[m][m] )
1300 rvec_dec(x[ad[i]],box[m]);
1302 /* determine grid index of atom */
1303 grididx[m]=x[ad[i]][m]*invdelta[m];
1304 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1307 copy_rvec(xtemp, x[ad[i]]); /* copy back */
1311 range_check(gx,0,ngrid[XX]);
1312 range_check(gy,0,ngrid[YY]);
1313 range_check(gz,0,ngrid[ZZ]);
1317 /* add atom to grid cell */
1319 newgrid=&(grid[gz][gy][gx].a[gr]);
1321 newgrid=&(grid[gz][gy][gx].d[gr]);
1322 if (newgrid->nr >= newgrid->maxnr) {
1324 DBB(newgrid->maxnr);
1325 srenew(newgrid->atoms, newgrid->maxnr);
1328 newgrid->atoms[newgrid->nr]=ad[i];
1336 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1340 for(gr=0; (gr<grNR); gr++) {
1342 for (zi=0; zi<ngrid[ZZ]; zi++)
1343 for (yi=0; yi<ngrid[YY]; yi++)
1344 for (xi=0; xi<ngrid[XX]; xi++) {
1345 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1351 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1352 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1353 * With a triclinic box all loops are 3 long, except when a cell is
1354 * located next to one of the box edges which is not parallel to the
1355 * x/y-plane, in that case all cells in a line or layer are searched.
1356 * This could be implemented slightly more efficient, but the code
1357 * would get much more complicated.
1359 #define B(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? 0 : (x-1))
1360 #define E(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? n-1 : (x+1))
1361 #define GRIDMOD(j,n) (j+n)%(n)
1362 #define LOOPGRIDINNER(x,y,z,xx,yy,zz,xo,yo,zo,n,bTric) \
1363 for(zz=B(n[ZZ],zo,bTric,FALSE); zz<=E(n[ZZ],zo,bTric,FALSE); zz++) { \
1364 z=GRIDMOD(zz,n[ZZ]); \
1365 for(yy=B(n[YY],yo,bTric,z==0||z==n[ZZ]-1); \
1366 yy<=E(n[YY],yo,bTric,z==0||z==n[ZZ]-1); yy++) { \
1367 y=GRIDMOD(yy,n[YY]); \
1368 for(xx=B(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); \
1369 xx<=E(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); xx++) { \
1370 x=GRIDMOD(xx,n[XX]);
1371 #define ENDLOOPGRIDINNER \
1377 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1379 int gr,x,y,z,sum[grNR];
1381 fprintf(fp,"grid %dx%dx%d\n",ngrid[XX],ngrid[YY],ngrid[ZZ]);
1382 for (gr=0; gr<grNR; gr++) {
1384 fprintf(fp,"GROUP %d (%s)\n",gr,grpnames[gr]);
1385 for (z=0; z<ngrid[ZZ]; z+=2) {
1386 fprintf(fp,"Z=%d,%d\n",z,z+1);
1387 for (y=0; y<ngrid[YY]; y++) {
1388 for (x=0; x<ngrid[XX]; x++) {
1389 fprintf(fp,"%3d",grid[x][y][z].d[gr].nr);
1390 sum[gr]+=grid[z][y][x].d[gr].nr;
1391 fprintf(fp,"%3d",grid[x][y][z].a[gr].nr);
1392 sum[gr]+=grid[z][y][x].a[gr].nr;
1396 if ( (z+1) < ngrid[ZZ] )
1397 for (x=0; x<ngrid[XX]; x++) {
1398 fprintf(fp,"%3d",grid[z+1][y][x].d[gr].nr);
1399 sum[gr]+=grid[z+1][y][x].d[gr].nr;
1400 fprintf(fp,"%3d",grid[z+1][y][x].a[gr].nr);
1401 sum[gr]+=grid[z+1][y][x].a[gr].nr;
1407 fprintf(fp,"TOTALS:");
1408 for (gr=0; gr<grNR; gr++)
1409 fprintf(fp," %d=%d",gr,sum[gr]);
1413 /* New GMX record! 5 * in a row. Congratulations!
1414 * Sorry, only four left.
1416 static void free_grid(ivec ngrid, t_gridcell ****grid)
1419 t_gridcell ***g = *grid;
1421 for (z=0; z<ngrid[ZZ]; z++) {
1422 for (y=0; y<ngrid[YY]; y++) {
1431 static void pbc_correct(rvec dx,matrix box,rvec hbox)
1434 for(m=DIM-1; m>=0; m--) {
1435 if ( dx[m] < -hbox[m] )
1436 rvec_inc(dx,box[m]);
1437 else if ( dx[m] >= hbox[m] )
1438 rvec_dec(dx,box[m]);
1442 void pbc_correct_gem(rvec dx,matrix box,rvec hbox)
1445 gmx_bool bDone = FALSE;
1448 for(m=DIM-1; m>=0; m--) {
1449 if ( dx[m] < -hbox[m] ) {
1451 rvec_inc(dx,box[m]);
1453 if ( dx[m] >= hbox[m] ) {
1455 rvec_dec(dx,box[m]);
1461 /* Added argument r2cut, changed contact and implemented
1462 * use of second cut-off.
1463 * - Erik Marklund, June 29, 2006
1465 static int is_hbond(t_hbdata *hb,int grpd,int grpa,int d,int a,
1466 real rcut, real r2cut, real ccut,
1467 rvec x[], gmx_bool bBox, matrix box,rvec hbox,
1468 real *d_ha, real *ang,gmx_bool bDA,int *hhh,
1469 gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1472 rvec r_da,r_ha,r_dh, r={0, 0, 0};
1474 real rc2,r2c2,rda2,rha2,ca;
1475 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1476 gmx_bool daSwap = FALSE;
1481 if (((id = donor_index(&hb->d,grpd,d)) == NOTSET) ||
1482 ((ja = acceptor_index(&hb->a,grpa,a)) == NOTSET))
1488 rvec_sub(x[d],x[a],r_da);
1489 /* Insert projection code here */
1491 /* if (d>a && ((isInterchangable(hb, d, a, grpd, grpa) && bMerge) || bContact)) */
1492 /* /\* Then this hbond will be found again, or it has already been found. *\/ */
1496 if (d>a && bMerge && (bContact || isInterchangable(hb, d, a, grpd, grpa))) { /* acceptor is also a donor and vice versa? */
1498 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1501 copy_rvec(r_da, r); /* Save this for later */
1502 pbc_correct_gem(r_da,box,hbox);
1504 pbc_correct_gem(r_da,box,hbox);
1507 rda2 = iprod(r_da,r_da);
1514 calcBoxDistance(hb->per->P, r, ri);
1515 *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */
1519 else if (rda2 < r2c2)
1526 if (bDA && (rda2 > rc2))
1529 for(h=0; (h < hb->d.nhydro[id]); h++) {
1530 hh = hb->d.hydro[id][h];
1533 rvec_sub(x[hh],x[a],r_ha);
1535 pbc_correct_gem(r_ha,box,hbox);
1537 rha2 = iprod(r_ha,r_ha);
1541 calcBoxDistance(hb->per->P, r, ri);
1542 *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */
1545 if (bDA || (!bDA && (rha2 <= rc2))) {
1546 rvec_sub(x[d],x[hh],r_dh);
1549 pbc_correct_gem(r_dh,box,hbox);
1551 pbc_correct_gem(r_dh,box,hbox);
1556 ca = cos_angle(r_dh,r_da);
1557 /* if angle is smaller, cos is larger */
1560 *d_ha = sqrt(bDA?rda2:rha2);
1566 if (bDA || (!bDA && HAinrange))
1572 /* Fixed previously undiscovered bug in the merge
1573 code, where the last frame of each hbond disappears.
1574 - Erik Marklund, June 1, 2006 */
1575 /* Added the following arguments:
1576 * ptmp[] - temporary periodicity hisory
1577 * a1 - identity of first acceptor/donor
1578 * a2 - identity of second acceptor/donor
1579 * - Erik Marklund, FEB 20 2010 */
1581 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1582 * Will do some more testing before removing the function entirely.
1583 * - Erik Marklund, MAY 10 2010 */
1584 static void do_merge(t_hbdata *hb,int ntmp,
1585 unsigned int htmp[],unsigned int gtmp[],PSTYPE ptmp[],
1586 t_hbond *hb0,t_hbond *hb1, int a1, int a2)
1588 /* Here we need to make sure we're treating periodicity in
1589 * the right way for the geminate recombination kinetics. */
1591 int m,mm,n00,n01,nn0,nnframes;
1595 /* Decide where to start from when merging */
1599 nnframes = max(n00 + hb0->nframes,n01 + hb1->nframes) - nn0;
1600 /* Initiate tmp arrays */
1601 for(m=0; (m<ntmp); m++) {
1606 /* Fill tmp arrays with values due to first HB */
1607 /* Once again '<' had to be replaced with '<='
1608 to catch the last frame in which the hbond
1610 - Erik Marklund, June 1, 2006 */
1611 for(m=0; (m<=hb0->nframes); m++) {
1613 htmp[mm] = is_hb(hb0->h[0],m);
1615 pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1616 if (pm > hb->per->nper)
1617 gmx_fatal(FARGS, "Illegal shift!");
1619 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1622 /* If we're doing geminate recompbination we usually don't need the distances.
1623 * Let's save some memory and time. */
1624 if (TRUE || !hb->bGem || hb->per->gemtype == gemAD){
1625 for(m=0; (m<=hb0->nframes); m++) {
1627 gtmp[mm] = is_hb(hb0->g[0],m);
1631 for(m=0; (m<=hb1->nframes); m++) {
1633 htmp[mm] = htmp[mm] || is_hb(hb1->h[0],m);
1634 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0],m);
1635 if (hb->bGem /* && ptmp[mm] != 0 */) {
1637 /* If this hbond has been seen before with donor and acceptor swapped,
1638 * then we need to find the mirrored (*-1) periodicity vector to truely
1639 * merge the hbond history. */
1640 pm = findMirror(getPshift(hb->per->pHist[a2][a1],m+hb1->n0), hb->per->p2i, hb->per->nper);
1641 /* Store index of mirror */
1642 if (pm > hb->per->nper)
1643 gmx_fatal(FARGS, "Illegal shift!");
1647 /* Reallocate target array */
1648 if (nnframes > hb0->maxframes) {
1649 srenew(hb0->h[0],4+nnframes/hb->wordlen);
1650 srenew(hb0->g[0],4+nnframes/hb->wordlen);
1652 clearPshift(&(hb->per->pHist[a1][a2]));
1654 /* Copy temp array to target array */
1655 for(m=0; (m<=nnframes); m++) {
1656 _set_hb(hb0->h[0],m,htmp[m]);
1657 _set_hb(hb0->g[0],m,gtmp[m]);
1659 addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
1662 /* Set scalar variables */
1664 hb0->maxframes = nnframes;
1667 /* Added argument bContact for nicer output.
1668 * Erik Marklund, June 29, 2006
1670 static void merge_hb(t_hbdata *hb,gmx_bool bTwo, gmx_bool bContact){
1671 int i,inrnew,indnew,j,ii,jj,m,id,ia,grp,ogrp,ntmp;
1672 unsigned int *htmp,*gtmp;
1677 indnew = hb->nrdist;
1679 /* Check whether donors are also acceptors */
1680 printf("Merging hbonds with Acceptor and Donor swapped\n");
1682 ntmp = 2*hb->max_frames;
1686 for(i=0; (i<hb->d.nrd); i++) {
1687 fprintf(stderr,"\r%d/%d",i+1,hb->d.nrd);
1689 ii = hb->a.aptr[id];
1690 for(j=0; (j<hb->a.nra); j++) {
1692 jj = hb->d.dptr[ia];
1693 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
1694 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j])))) {
1695 hb0 = hb->hbmap[i][j];
1696 hb1 = hb->hbmap[jj][ii];
1697 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0])) {
1698 do_merge(hb,ntmp,htmp,gtmp,ptmp,hb0,hb1,i,j);
1699 if (ISHB(hb1->history[0]))
1701 else if (ISDIST(hb1->history[0]))
1705 gmx_incons("No contact history");
1707 gmx_incons("Neither hydrogen bond nor distance");
1711 clearPshift(&(hb->per->pHist[jj][ii]));
1715 hb1->history[0] = hbNo;
1720 fprintf(stderr,"\n");
1721 printf("- Reduced number of hbonds from %d to %d\n",hb->nrhb,inrnew);
1722 printf("- Reduced number of distances from %d to %d\n",hb->nrdist,indnew);
1724 hb->nrdist = indnew;
1730 static void do_nhb_dist(FILE *fp,t_hbdata *hb,real t)
1732 int i,j,k,n_bound[MAXHH],nbtot;
1736 /* Set array to 0 */
1737 for(k=0; (k<MAXHH); k++)
1739 /* Loop over possible donors */
1740 for(i=0; (i<hb->d.nrd); i++) {
1741 for(j=0; (j<hb->d.nhydro[i]); j++)
1742 n_bound[hb->d.nhbonds[i][j]]++;
1744 fprintf(fp,"%12.5e",t);
1746 for(k=0; (k<MAXHH); k++) {
1747 fprintf(fp," %8d",n_bound[k]);
1748 nbtot += n_bound[k]*k;
1750 fprintf(fp," %8d\n",nbtot);
1753 /* Added argument bContact in do_hblife(...). Also
1754 * added support for -contact in function body.
1755 * - Erik Marklund, May 31, 2006 */
1756 /* Changed the contact code slightly.
1757 * - Erik Marklund, June 29, 2006
1759 static void do_hblife(const char *fn,t_hbdata *hb,gmx_bool bMerge,gmx_bool bContact,
1760 const output_env_t oenv)
1763 const char *leg[] = { "p(t)", "t p(t)" };
1765 int i,j,j0,k,m,nh,ihb,ohb,nhydro,ndump=0;
1766 int nframes = hb->nframes;
1769 double sum,integral;
1772 snew(h,hb->maxhydro);
1773 snew(histo,nframes+1);
1774 /* Total number of hbonds analyzed here */
1775 for(i=0; (i<hb->d.nrd); i++) {
1776 for(k=0; (k<hb->a.nra); k++) {
1777 hbh = hb->hbmap[i][k];
1789 for(m=0; (m<hb->maxhydro); m++)
1791 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
1794 for(nh=0; (nh<nhydro); nh++) {
1798 /* Changed '<' into '<=' below, just like I
1799 did in the hbm-output-loop in the main code.
1800 - Erik Marklund, May 31, 2006
1802 for(j=0; (j<=hbh->nframes); j++) {
1803 ihb = is_hb(h[nh],j);
1804 if (debug && (ndump < 10))
1805 fprintf(debug,"%5d %5d\n",j,ihb);
1821 fprintf(stderr,"\n");
1823 fp = xvgropen(fn,"Uninterrupted contact lifetime","Time (ps)","()",oenv);
1825 fp = xvgropen(fn,"Uninterrupted hydrogen bond lifetime","Time (ps)","()",
1828 xvgr_legend(fp,asize(leg),leg,oenv);
1830 while ((j0 > 0) && (histo[j0] == 0))
1833 for(i=0; (i<=j0); i++)
1835 dt = hb->time[1]-hb->time[0];
1838 for(i=1; (i<=j0); i++) {
1839 t = hb->time[i] - hb->time[0] - 0.5*dt;
1840 x1 = t*histo[i]/sum;
1841 fprintf(fp,"%8.3f %10.3e %10.3e\n",t,histo[i]/sum,x1);
1846 printf("%s lifetime = %.2f ps\n", bContact?"Contact":"HB", integral);
1847 printf("Note that the lifetime obtained in this manner is close to useless\n");
1848 printf("Use the -ac option instead and check the Forward lifetime\n");
1849 please_cite(stdout,"Spoel2006b");
1854 /* Changed argument bMerge into oneHB to handle contacts properly.
1855 * - Erik Marklund, June 29, 2006
1857 static void dump_ac(t_hbdata *hb,gmx_bool oneHB,int nDump)
1860 int i,j,k,m,nd,ihb,idist;
1861 int nframes = hb->nframes;
1867 fp = ffopen("debug-ac.xvg","w");
1868 for(j=0; (j<nframes); j++) {
1869 fprintf(fp,"%10.3f",hb->time[j]);
1870 for(i=nd=0; (i<hb->d.nrd) && (nd < nDump); i++) {
1871 for(k=0; (k<hb->a.nra) && (nd < nDump); k++) {
1874 hbh = hb->hbmap[i][k];
1877 ihb = is_hb(hbh->h[0],j);
1878 idist = is_hb(hbh->g[0],j);
1883 for(m=0; (m<hb->maxhydro) && !ihb ; m++) {
1884 ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m],j));
1885 idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m],j));
1887 /* This is not correct! */
1888 /* What isn't correct? -Erik M */
1892 fprintf(fp," %1d-%1d",ihb,idist);
1902 static real calc_dg(real tau,real temp)
1910 return kbt*log(kbt*tau/PLANCK);
1914 int n0,n1,nparams,ndelta;
1916 real *t,*ct,*nt,*kt,*sigma_ct,*sigma_nt,*sigma_kt;
1920 #include <gsl/gsl_multimin.h>
1921 #include <gsl/gsl_sf.h>
1922 #include <gsl/gsl_version.h>
1924 static double my_f(const gsl_vector *v,void *params)
1926 t_luzar *tl = (t_luzar *)params;
1928 double tol=1e-16,chi2=0;
1932 for(i=0; (i<tl->nparams); i++) {
1933 tl->kkk[i] = gsl_vector_get(v, i);
1938 for(i=tl->n0; (i<tl->n1); i+=tl->ndelta) {
1939 di=sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
1942 chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
1945 fprintf(stderr,"WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
1946 "di = %g k = %g kp = %g\n",tl->sigma_ct[i],
1947 tl->sigma_nt[i],tl->sigma_kt[i],di,k,kp);
1950 chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
1955 static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
1963 const gsl_multimin_fminimizer_type *T;
1964 gsl_multimin_fminimizer *s;
1967 gsl_multimin_function my_func;
1970 my_func.n = tl->nparams;
1971 my_func.params = (void *) tl;
1973 /* Starting point */
1974 x = gsl_vector_alloc (my_func.n);
1975 for(i=0; (i<my_func.n); i++)
1976 gsl_vector_set (x, i, tl->kkk[i]);
1978 /* Step size, different for each of the parameters */
1979 dx = gsl_vector_alloc (my_func.n);
1980 for(i=0; (i<my_func.n); i++)
1981 gsl_vector_set (dx, i, 0.01*tl->kkk[i]);
1983 T = gsl_multimin_fminimizer_nmsimplex;
1984 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
1986 gsl_multimin_fminimizer_set (s, &my_func, x, dx);
1987 gsl_vector_free (x);
1988 gsl_vector_free (dx);
1991 fprintf(fp,"%5s %12s %12s %12s %12s\n","Iter","k","kp","NM Size","Chi2");
1995 status = gsl_multimin_fminimizer_iterate (s);
1998 gmx_fatal(FARGS,"Something went wrong in the iteration in minimizer %s",
1999 gsl_multimin_fminimizer_name(s));
2001 d2 = gsl_multimin_fminimizer_minimum(s);
2002 size = gsl_multimin_fminimizer_size(s);
2003 status = gsl_multimin_test_size(size,tol);
2005 if (status == GSL_SUCCESS)
2007 fprintf(fp,"Minimum found using %s at:\n",
2008 gsl_multimin_fminimizer_name(s));
2011 fprintf(fp,"%5d", iter);
2012 for(i=0; (i<my_func.n); i++)
2013 fprintf(fp," %12.4e",gsl_vector_get (s->x,i));
2014 fprintf (fp," %12.4e %12.4e\n",size,d2);
2017 while ((status == GSL_CONTINUE) && (iter < maxiter));
2019 gsl_multimin_fminimizer_free (s);
2024 static real quality_of_fit(real chi2,int N)
2026 return gsl_sf_gamma_inc_Q((N-2)/2.0,chi2/2.0);
2030 static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
2033 fprintf(stderr,"This program needs the GNU scientific library to work.\n");
2037 static real quality_of_fit(real chi2,int N)
2039 fprintf(stderr,"This program needs the GNU scientific library to work.\n");
2046 static real compute_weighted_rates(int n,real t[],real ct[],real nt[],
2047 real kt[],real sigma_ct[],real sigma_nt[],
2048 real sigma_kt[],real *k,real *kp,
2049 real *sigma_k,real *sigma_kp,
2055 real kkk=0,kkp=0,kk2=0,kp2=0,chi2;
2060 for(i=0; (i<n); i++) {
2061 if (t[i] >= fit_start)
2072 tl.sigma_ct = sigma_ct;
2073 tl.sigma_nt = sigma_nt;
2074 tl.sigma_kt = sigma_kt;
2078 chi2 = optimize_luzar_parameters(debug,&tl,1000,1e-3);
2080 *kp = tl.kkk[1] = *kp;
2082 for(j=0; (j<NK); j++) {
2083 (void) optimize_luzar_parameters(debug,&tl,1000,1e-3);
2086 kk2 += sqr(tl.kkk[0]);
2087 kp2 += sqr(tl.kkk[1]);
2090 *sigma_k = sqrt(kk2/NK - sqr(kkk/NK));
2091 *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2096 static void smooth_tail(int n,real t[],real c[],real sigma_c[],real start,
2097 const output_env_t oenv)
2100 real e_1,fitparm[4];
2104 for(i=0; (i<n); i++)
2112 do_lmfit(n,c,sigma_c,0,t,start,t[n-1],oenv,bDebugMode(),effnEXP2,fitparm,0);
2115 void analyse_corr(int n,real t[],real ct[],real nt[],real kt[],
2116 real sigma_ct[],real sigma_nt[],real sigma_kt[],
2117 real fit_start,real temp,real smooth_tail_start,
2118 const output_env_t oenv)
2121 real k=1,kp=1,kow=1;
2122 real Q=0,chi22,chi2,dg,dgp,tau_hb,dtau,tau_rlx,e_1,dt,sigma_k,sigma_kp,ddg;
2123 double tmp,sn2=0,sc2=0,sk2=0,scn=0,sck=0,snk=0;
2124 gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2126 if (smooth_tail_start >= 0) {
2127 smooth_tail(n,t,ct,sigma_ct,smooth_tail_start,oenv);
2128 smooth_tail(n,t,nt,sigma_nt,smooth_tail_start,oenv);
2129 smooth_tail(n,t,kt,sigma_kt,smooth_tail_start,oenv);
2131 for(i0=0; (i0<n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2134 for(i=i0; (i<n); i++) {
2142 printf("Hydrogen bond thermodynamics at T = %g K\n",temp);
2143 tmp = (sn2*sc2-sqr(scn));
2144 if ((tmp > 0) && (sn2 > 0)) {
2145 k = (sn2*sck-scn*snk)/tmp;
2146 kp = (k*scn-snk)/sn2;
2149 for(i=i0; (i<n); i++) {
2150 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2152 chi22 = compute_weighted_rates(n,t,ct,nt,kt,sigma_ct,sigma_nt,
2154 &sigma_k,&sigma_kp,fit_start);
2155 Q = quality_of_fit(chi2,2);
2156 ddg = BOLTZ*temp*sigma_k/k;
2157 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2159 printf("The Rate and Delta G are followed by an error estimate\n");
2160 printf("----------------------------------------------------------\n"
2161 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
2162 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2163 k,sigma_k,1/k,calc_dg(1/k,temp),ddg);
2164 ddg = BOLTZ*temp*sigma_kp/kp;
2165 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2166 kp,sigma_kp,1/kp,calc_dg(1/kp,temp),ddg);
2170 for(i=i0; (i<n); i++) {
2171 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2173 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2175 printf("--------------------------------------------------\n"
2176 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
2177 printf("Forward %10.3f %8.3f %10.3f %10g\n",
2178 k,1/k,calc_dg(1/k,temp),chi2);
2179 printf("Backward %10.3f %8.3f %10.3f\n",
2180 kp,1/kp,calc_dg(1/kp,temp));
2185 printf("One-way %10.3f %s%8.3f %10.3f\n",
2186 kow,bError ? " " : "",1/kow,calc_dg(1/kow,temp));
2189 printf(" - Numerical problems computing HB thermodynamics:\n"
2190 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
2191 sc2,sn2,sk2,sck,snk,scn);
2192 /* Determine integral of the correlation function */
2193 tau_hb = evaluate_integral(n,t,ct,NULL,(t[n-1]-t[0])/2,&dtau);
2194 printf("Integral %10.3f %s%8.3f %10.3f\n",1/tau_hb,
2195 bError ? " " : "",tau_hb,calc_dg(tau_hb,temp));
2197 for(i=0; (i<n-2); i++) {
2198 if ((ct[i] > e_1) && (ct[i+1] <= e_1)) {
2203 /* Determine tau_relax from linear interpolation */
2204 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2205 printf("Relaxation %10.3f %8.3f %s%10.3f\n",1/tau_rlx,
2206 tau_rlx,bError ? " " : "",
2207 calc_dg(tau_rlx,temp));
2211 printf("Correlation functions too short to compute thermodynamics\n");
2214 void compute_derivative(int nn,real x[],real y[],real dydx[])
2218 /* Compute k(t) = dc(t)/dt */
2219 for(j=1; (j<nn-1); j++)
2220 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2221 /* Extrapolate endpoints */
2222 dydx[0] = 2*dydx[1] - dydx[2];
2223 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2226 static void parallel_print(int *data, int nThreads)
2228 /* This prints the donors on which each tread is currently working. */
2231 fprintf(stderr, "\r");
2232 for (i=0; i<nThreads; i++)
2233 fprintf(stderr, "%-7i",data[i]);
2236 static void normalizeACF(real *ct, real *gt, int len)
2238 real ct_fac, gt_fac;
2241 /* Xu and Berne use the same normalization constant */
2244 gt_fac = (gt!=NULL && gt[0]!=0) ? 1.0/gt[0] : 0;
2245 printf("Normalization for c(t) = %g for gh(t) = %g\n",ct_fac,gt_fac);
2246 for (i=0; i<len; i++)
2254 /* Added argument bContact in do_hbac(...). Also
2255 * added support for -contact in the actual code.
2256 * - Erik Marklund, May 31, 2006 */
2257 /* Changed contact code and added argument R2
2258 * - Erik Marklund, June 29, 2006
2260 static void do_hbac(const char *fn,t_hbdata *hb,
2261 int nDump,gmx_bool bMerge,gmx_bool bContact, real fit_start,
2262 real temp,gmx_bool R2,real smooth_tail_start, const output_env_t oenv,
2263 t_gemParams *params, const char *gemType, int nThreads,
2264 const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2267 int i,j,k,m,n,o,nd,ihb,idist,n2,nn,iter,nSets;
2268 const char *legNN[] = { "Ac(t)",
2270 static char **legGem;
2272 const char *legLuzar[] = { "Ac\\sfin sys\\v{}\\z{}(t)",
2274 "Cc\\scontact,hb\\v{}\\z{}(t)",
2275 "-dAc\\sfs\\v{}\\z{}/dt" };
2276 gmx_bool bNorm=FALSE;
2279 real *rhbex=NULL,*ht,*gt,*ght,*dght,*kt;
2280 real *ct,*p_ct,tail,tail2,dtail,ct_fac,ght_fac,*cct;
2281 const real tol = 1e-3;
2282 int nframes = hb->nframes,nf;
2283 unsigned int **h,**g;
2284 int nh,nhbonds,nhydro,ngh;
2286 PSTYPE p, *pfound = NULL, np;
2288 int *ptimes=NULL, *poff=NULL, anhb, n0, mMax=INT_MIN;
2289 real **rHbExGem = NULL;
2293 double *ctdouble, *timedouble, *fittedct;
2294 double fittolerance=0.1;
2296 enum {AC_NONE, AC_NN, AC_GEM, AC_LUZAR};
2300 int *dondata=NULL, thisThread;
2304 printf("Doing autocorrelation ");
2306 /* Decide what kind of ACF calculations to do. */
2307 if (NN > NN_NONE && NN < NN_NR) {
2308 #ifdef HAVE_NN_LOOPS
2310 printf("using the energy estimate.\n");
2313 printf("Can't do the NN-loop. Yet.\n");
2315 } else if (hb->bGem) {
2317 printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2319 nSets = 1 + (bBallistic ? 1:0) + (bGemFit ? 1:0);
2320 snew(legGem, nSets);
2321 for (i=0;i<nSets;i++)
2322 snew(legGem[i], 128);
2323 sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2325 sprintf(legGem[1], "Ac'(t)");
2327 sprintf(legGem[(bBallistic ? 3:2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2331 printf("according to the theory of Luzar and Chandler.\n");
2335 /* build hbexist matrix in reals for autocorr */
2336 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2338 while (n2 < nframes)
2343 if (acType != AC_NN ||
2350 snew(h,hb->maxhydro);
2351 snew(g,hb->maxhydro);
2354 /* Dump hbonds for debugging */
2355 dump_ac(hb,bMerge||bContact,nDump);
2357 /* Total number of hbonds analyzed here */
2362 /* ------------------------------------------------
2363 * I got tired of waiting for the acf calculations
2364 * and parallelized it with openMP
2365 * set environment variable CFLAGS = "-fopenmp" when running
2366 * configure and define DOUSEOPENMP to make use of it.
2369 #ifdef HAVE_OPENMP /* ================================================= \
2370 * Set up the OpenMP stuff, |
2371 * like the number of threads and such |
2373 if (acType != AC_LUZAR)
2375 #if (_OPENMP >= 200805) /* =====================\ */
2376 nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit());
2378 nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
2379 #endif /* _OPENMP >= 200805 ====================/ */
2381 omp_set_num_threads(nThreads);
2382 snew(dondata, nThreads);
2383 for (i=0; i<nThreads; i++)
2385 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2386 "Expect close to linear scaling over this donor-loop.\n", nThreads);
2388 fprintf(stderr, "Donors: [thread no]\n");
2391 for (i=0;i<nThreads;i++) {
2392 snprintf(tmpstr, 7, "[%i]", i);
2393 fprintf(stderr, "%-7s", tmpstr);
2396 fprintf(stderr, "\n"); /* | */
2398 #endif /* HAVE_OPENMP ===================================================/ */
2401 /* Build the ACF according to acType */
2406 #ifdef HAVE_NN_LOOPS
2407 /* Here we're using the estimated energy for the hydrogen bonds. */
2409 #ifdef HAVE_OPENMP /* ==================================\ */
2410 #pragma omp parallel \
2411 private(i, j, k, nh, E, rhbex, thisThread), \
2415 thisThread = omp_get_thread_num();
2417 #endif /* ==============================================/ */
2420 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2422 #ifdef HAVE_OPENMP /* ################################################## \
2427 #pragma omp for schedule (dynamic)
2429 for (i=0; i<hb->d.nrd; i++) /* loop over donors */
2431 #ifdef HAVE_OPENMP /* ====== Write some output ======\ */
2432 #pragma omp critical
2434 dondata[thisThread] = i;
2435 parallel_print(dondata, nThreads);
2438 fprintf(stderr, "\r %i", i);
2439 #endif /* ===========================================/ */
2441 for (j=0; j<hb->a.nra; j++) /* loop over acceptors */
2443 for (nh=0; nh<hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2445 E = hb->hbE.E[i][j][nh];
2448 for (k=0; k<nframes; k++)
2450 if (E[k] != NONSENSE_E)
2451 rhbex[k] = (real)E[k];
2454 low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rhbex),hb->time[1]-hb->time[0],
2455 eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2457 #pragma omp critical
2460 for(k=0; (k<nn); k++)
2471 } /* End of parallel block # */
2472 /* ##############################################################/ */
2475 normalizeACF(ct, NULL, nn);
2477 snew(timedouble, nn);
2478 for (j=0; j<nn; j++)
2480 timedouble[j]=(double)(hb->time[j]);
2481 ctdouble[j]=(double)(ct[j]);
2484 /* Remove ballistic term */
2485 /* Ballistic component removal and fitting to the reversible geminate recombination model
2486 * will be taken out for the time being. First of all, one can remove the ballistic
2487 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2488 * problems with the robustness of the fitting to the model. More work is needed.
2489 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2490 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2491 * a BSD-licence that can do the job.
2493 * - Erik Marklund, June 18 2010.
2495 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2496 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2498 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2499 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2501 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation","Time (ps)","C(t)");
2502 xvgr_legend(fp,asize(legNN),legNN);
2504 for(j=0; (j<nn); j++)
2505 fprintf(fp,"%10g %10g %10g\n",
2506 hb->time[j]-hb->time[0],
2513 #endif /* HAVE_NN_LOOPS */
2514 break; /* case AC_NN */
2518 memset(ct,0,2*n2*sizeof(real));
2520 fprintf(stderr, "Donor:\n");
2523 #define __ACDATA p_ct
2526 #ifdef HAVE_OPENMP /* =========================================\
2528 #pragma omp parallel default(none) \
2529 private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2530 pfound, poff, rHbExGem, p, ihb, mMax, \
2532 shared(hb, dondata, ct, nn, nThreads, n2, stderr, bNorm, \
2533 nframes, bMerge, bContact)
2534 { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
2537 thisThread = omp_get_thread_num();
2538 snew(h,hb->maxhydro);
2539 snew(g,hb->maxhydro);
2546 memset(p_ct,0,2*n2*sizeof(real));
2548 /* I'm using a chunk size of 1, since I expect \
2549 * the overhead to be really small compared \
2550 * to the actual calculations \ */
2551 #pragma omp for schedule(dynamic,1) nowait /* \ */
2552 #endif /* HAVE_OPENMP =========================================/ */
2554 for (i=0; i<hb->d.nrd; i++) {
2556 #pragma omp critical
2558 dondata[thisThread] = i;
2559 parallel_print(dondata, nThreads);
2562 fprintf(stderr, "\r %i", i);
2565 for (k=0; k<hb->a.nra; k++) {
2566 for (nh=0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) {
2567 hbh = hb->hbmap[i][k];
2569 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
2570 * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
2571 * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
2572 pHist = &(hb->per->pHist[i][k]);
2573 if (ISHB(hbh->history[nh]) && pHist->len != 0) {
2575 /* No need for a critical section */
2576 /* #ifdef HAVE_OPENMP */
2577 /* #pragma omp critical */
2581 g[nh] = hb->per->gemtype==gemAD ? hbh->g[nh] : NULL;
2585 /* count the number of periodic shifts encountered and store
2586 * them in separate arrays. */
2588 for (j=0; j<pHist->len; j++)
2591 for (m=0; m<=np; m++) {
2592 if (m == np) { /* p not recognized in list. Add it and set up new array. */
2594 if (np>hb->per->nper)
2595 gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
2596 if (m>=mMax) { /* Extend the arrays.
2597 * Doing it like this, using mMax to keep track of the sizes,
2598 * eleviates the need for freeing and re-allocating the arrays
2599 * when taking on the next donor-acceptor pair */
2601 srenew(pfound,np); /* The list of found periodic shifts. */
2602 srenew(rHbExGem,np); /* The hb existence functions (-aver_hb). */
2603 snew(rHbExGem[m],2*n2);
2607 /* This shouldn't have to be critical, right? */
2608 /* #ifdef HAVE_OPENMP */
2609 /* #pragma omp critical */
2612 if (rHbExGem != NULL && rHbExGem[m] != NULL) {
2613 /* This must be done, as this array was most likey
2614 * used to store stuff in some previous iteration. */
2615 memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
2618 fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
2627 } /* m: Loop over found shifts */
2628 } /* j: Loop over shifts */
2630 /* Now unpack and disentangle the existence funtions. */
2631 for (j=0; j<nf; j++) {
2637 * pfound: list of periodic shifts found for this pair.
2638 * poff: list of frame offsets; that is, the first
2639 * frame a hbond has a particular periodic shift. */
2640 p = getPshift(*pHist, j+n0);
2643 for (m=0; m<np; m++)
2648 gmx_fatal(FARGS,"Shift not found, but must be there.");
2651 ihb = is_hb(h[nh],j) || ((hb->per->gemtype!=gemAD || j==0) ? FALSE : is_hb(g[nh],j));
2655 poff[m] = j; /* Here's where the first hbond with shift p is,
2656 * relative to the start of h[0].*/
2658 gmx_fatal(FARGS, "j<poff[m]");
2659 rHbExGem[m][j-poff[m]] += 1;
2664 /* Now, build ac. */
2665 for (m=0; m<np; m++) {
2666 if (rHbExGem[m][0]>0 && n0+poff[m]<nn/* && m==0 */) {
2667 low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rHbExGem[m]),hb->time[1]-hb->time[0],
2668 eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2669 for(j=0; (j<nn); j++)
2670 __ACDATA[j] += rHbExGem[m][j];
2672 } /* Building of ac. */
2675 } /* hydrogen loop */
2676 } /* acceptor loop */
2679 for (m=0; m<=mMax; m++) {
2688 #ifdef HAVE_OPENMP /* =======================================\ */
2689 #pragma omp critical
2691 for (i=0; i<nn; i++)
2696 } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
2698 #endif /* HAVE_OPENMP =======================================/ */
2700 normalizeACF(ct, NULL, nn);
2702 fprintf(stderr, "\n\nACF successfully calculated.\n");
2704 /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
2707 snew(timedouble, nn);
2711 timedouble[j]=(double)(hb->time[j]);
2712 ctdouble[j]=(double)(ct[j]);
2715 /* Remove ballistic term */
2716 /* Ballistic component removal and fitting to the reversible geminate recombination model
2717 * will be taken out for the time being. First of all, one can remove the ballistic
2718 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2719 * problems with the robustness of the fitting to the model. More work is needed.
2720 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2721 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2722 * a BSD-licence that can do the job.
2724 * - Erik Marklund, June 18 2010.
2726 /* if (bBallistic) { */
2727 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2728 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2730 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2731 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2734 /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
2738 fp = xvgropen(fn, "Contact Autocorrelation","Time (ps)","C(t)",oenv);
2740 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation","Time (ps)","C(t)",oenv);
2741 xvgr_legend(fp,asize(legGem),(const char**)legGem,oenv);
2743 for(j=0; (j<nn); j++)
2745 fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0],ct[j]);
2747 fprintf(fp," %10g", ctdouble[j]);
2749 fprintf(fp," %10g", fittedct[j]);
2759 break; /* case AC_GEM */
2772 for(i=0; (i<hb->d.nrd); i++) {
2773 for(k=0; (k<hb->a.nra); k++) {
2775 hbh = hb->hbmap[i][k];
2778 if (bMerge || bContact) {
2779 if (ISHB(hbh->history[0])) {
2786 for(m=0; (m<hb->maxhydro); m++)
2787 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) {
2788 g[nhydro] = hbh->g[m];
2789 h[nhydro] = hbh->h[m];
2795 for(nh=0; (nh<nhydro); nh++) {
2796 int nrint = bContact ? hb->nrdist : hb->nrhb;
2797 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
2798 fprintf(stderr,"\rACF %d/%d",nhbonds+1,nrint);
2800 for(j=0; (j<nframes); j++) {
2801 /* Changed '<' into '<=' below, just like I did in
2802 the hbm-output-loop in the gmx_hbond() block.
2803 - Erik Marklund, May 31, 2006 */
2805 ihb = is_hb(h[nh],j);
2806 idist = is_hb(g[nh],j);
2812 /* For contacts: if a second cut-off is provided, use it,
2813 * otherwise use g(t) = 1-h(t) */
2814 if (!R2 && bContact)
2817 gt[j] = idist*(1-ihb);
2823 /* The autocorrelation function is normalized after summation only */
2824 low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&rhbex,hb->time[1]-hb->time[0],
2825 eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2827 /* Cross correlation analysis for thermodynamics */
2828 for(j=nframes; (j<n2); j++) {
2833 cross_corr(n2,ht,gt,dght);
2835 for(j=0; (j<nn); j++) {
2843 fprintf(stderr,"\n");
2846 normalizeACF(ct, gt, nn);
2848 /* Determine tail value for statistics */
2851 for(j=nn/2; (j<nn); j++) {
2853 tail2 += ct[j]*ct[j];
2855 tail /= (nn - nn/2);
2856 tail2 /= (nn - nn/2);
2857 dtail = sqrt(tail2-tail*tail);
2859 /* Check whether the ACF is long enough */
2861 printf("\nWARNING: Correlation function is probably not long enough\n"
2862 "because the standard deviation in the tail of C(t) > %g\n"
2863 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2866 for(j=0; (j<nn); j++) {
2868 ct[j] = (cct[j]-tail)/(1-tail);
2870 /* Compute negative derivative k(t) = -dc(t)/dt */
2871 compute_derivative(nn,hb->time,ct,kt);
2872 for(j=0; (j<nn); j++)
2877 fp = xvgropen(fn, "Contact Autocorrelation","Time (ps)","C(t)",oenv);
2879 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation","Time (ps)","C(t)",oenv);
2880 xvgr_legend(fp,asize(legLuzar),legLuzar, oenv);
2883 for(j=0; (j<nn); j++)
2884 fprintf(fp,"%10g %10g %10g %10g %10g\n",
2885 hb->time[j]-hb->time[0],ct[j],cct[j],ght[j],kt[j]);
2888 analyse_corr(nn,hb->time,ct,ght,kt,NULL,NULL,NULL,
2889 fit_start,temp,smooth_tail_start,oenv);
2891 do_view(oenv,fn,NULL);
2903 break; /* case AC_LUZAR */
2906 gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
2907 } /* switch (acType) */
2910 static void init_hbframe(t_hbdata *hb,int nframes,real t)
2914 hb->time[nframes] = t;
2915 hb->nhb[nframes] = 0;
2916 hb->ndist[nframes] = 0;
2917 for (i=0; (i<max_hx); i++)
2918 hb->nhx[nframes][i]=0;
2919 /* Loop invalidated */
2920 if (hb->bHBmap && 0)
2921 for (i=0; (i<hb->d.nrd); i++)
2922 for (j=0; (j<hb->a.nra); j++)
2923 for (m=0; (m<hb->maxhydro); m++)
2924 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
2925 set_hb(hb,i,m,j,nframes,HB_NO);
2926 /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
2929 static void analyse_donor_props(const char *fn,t_hbdata *hb,int nframes,real t,
2930 const output_env_t oenv)
2932 static FILE *fp = NULL;
2933 const char *leg[] = { "Nbound", "Nfree" };
2934 int i,j,k,nbound,nb,nhtot;
2939 fp = xvgropen(fn,"Donor properties","Time (ps)","Number",oenv);
2940 xvgr_legend(fp,asize(leg),leg,oenv);
2944 for(i=0; (i<hb->d.nrd); i++) {
2945 for(k=0; (k<hb->d.nhydro[i]); k++) {
2948 for(j=0; (j<hb->a.nra) && (nb == 0); j++) {
2949 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
2950 is_hb(hb->hbmap[i][j]->h[k],nframes))
2956 fprintf(fp,"%10.3e %6d %6d\n",t,nbound,nhtot-nbound);
2959 static void dump_hbmap(t_hbdata *hb,
2960 int nfile,t_filenm fnm[],gmx_bool bTwo,
2961 gmx_bool bContact, int isize[],int *index[],char *grpnames[],
2965 int ddd,hhh,aaa,i,j,k,m,grp;
2966 char ds[32],hs[32],as[32];
2969 fp = opt2FILE("-hbn",nfile,fnm,"w");
2970 if (opt2bSet("-g",nfile,fnm)) {
2971 fplog = ffopen(opt2fn("-g",nfile,fnm),"w");
2973 fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
2975 fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
2979 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp++) {
2980 fprintf(fp,"[ %s ]",grpnames[grp]);
2981 for (i=0; i<isize[grp]; i++) {
2982 fprintf(fp,(i%15)?" ":"\n");
2983 fprintf(fp," %4u",index[grp][i]+1);
2987 Added -contact support below.
2988 - Erik Marklund, May 29, 2006
2991 fprintf(fp,"[ donors_hydrogens_%s ]\n",grpnames[grp]);
2992 for (i=0; (i<hb->d.nrd); i++) {
2993 if (hb->d.grp[i] == grp) {
2994 for(j=0; (j<hb->d.nhydro[i]); j++)
2995 fprintf(fp," %4u %4u",hb->d.don[i]+1,
2996 hb->d.hydro[i][j]+1);
3001 fprintf(fp,"[ acceptors_%s ]",grpnames[grp]);
3002 for (i=0; (i<hb->a.nra); i++) {
3003 if (hb->a.grp[i] == grp) {
3004 fprintf(fp,(i%15 && !first)?" ":"\n");
3005 fprintf(fp," %4u",hb->a.acc[i]+1);
3013 fprintf(fp,bContact ? "[ contacts_%s-%s ]\n" :
3014 "[ hbonds_%s-%s ]\n",grpnames[0],grpnames[1]);
3016 fprintf(fp,bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n",grpnames[0]);
3018 for(i=0; (i<hb->d.nrd); i++) {
3020 for(k=0; (k<hb->a.nra); k++) {
3022 for(m=0; (m<hb->d.nhydro[i]); m++) {
3023 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m])) {
3024 sprintf(ds,"%s",mkatomname(atoms,ddd));
3025 sprintf(as,"%s",mkatomname(atoms,aaa));
3027 fprintf(fp," %6u %6u\n",ddd+1,aaa+1);
3029 fprintf(fplog,"%12s %12s\n",ds,as);
3031 hhh = hb->d.hydro[i][m];
3032 sprintf(hs,"%s",mkatomname(atoms,hhh));
3033 fprintf(fp," %6u %6u %6u\n",ddd+1,hhh+1,aaa+1);
3035 fprintf(fplog,"%12s %12s %12s\n",ds,hs,as);
3047 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3048 * It mimics add_frames() and init_frame() to some extent. */
3049 static void sync_hbdata(t_hbdata *hb, t_hbdata *p_hb,
3050 int nframes, real t)
3053 if (nframes >= p_hb->max_frames)
3055 p_hb->max_frames += 4096;
3056 srenew(p_hb->nhb, p_hb->max_frames);
3057 srenew(p_hb->ndist, p_hb->max_frames);
3058 srenew(p_hb->n_bound, p_hb->max_frames);
3059 srenew(p_hb->nhx, p_hb->max_frames);
3061 srenew(p_hb->danr, p_hb->max_frames);
3062 memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3063 memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3064 p_hb->nhb[nframes] = 0;
3065 p_hb->ndist[nframes] = 0;
3068 p_hb->nframes = nframes;
3071 /* p_hb->nhx[nframes][i] */
3073 memset(&(p_hb->nhx[nframes]), 0 ,sizeof(int)*max_hx); /* zero the helix count for this frame */
3075 /* hb->per will remain constant througout the frame loop,
3076 * even though the data its members point to will change,
3077 * hence no need for re-syncing. */
3081 int gmx_hbond(int argc,char *argv[])
3083 const char *desc[] = {
3084 "g_hbond computes and analyzes hydrogen bonds. Hydrogen bonds are",
3085 "determined based on cutoffs for the angle Acceptor - Donor - Hydrogen",
3086 "(zero is extended) and the distance Hydrogen - Acceptor.",
3087 "OH and NH groups are regarded as donors, O is an acceptor always,",
3088 "N is an acceptor by default, but this can be switched using",
3089 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3090 "to the first preceding non-hydrogen atom.[PAR]",
3092 "You need to specify two groups for analysis, which must be either",
3093 "identical or non-overlapping. All hydrogen bonds between the two",
3094 "groups are analyzed.[PAR]",
3096 "If you set -shell, you will be asked for an additional index group",
3097 "which should contain exactly one atom. In this case, only hydrogen",
3098 "bonds between atoms within the shell distance from the one atom are",
3101 /* "It is also possible to analyse specific hydrogen bonds with",
3102 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3103 "Donor Hydrogen Acceptor, in the following way:[PAR]",
3111 "Note that the triplets need not be on separate lines.",
3112 "Each atom triplet specifies a hydrogen bond to be analyzed,",
3113 "note also that no check is made for the types of atoms.[PAR]",
3115 "[BB]Output:[bb][BR]",
3116 "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
3117 "[TT]-ac[tt]: average over all autocorrelations of the existence",
3118 "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3119 "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3120 "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
3121 "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
3122 "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3123 "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3124 "with helices in proteins.[BR]",
3125 "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
3126 "for selected groups, all hydrogen bonded atoms from all groups and",
3127 "all solvent atoms involved in insertion.[BR]",
3128 "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
3129 "frames, this also contains information on solvent insertion",
3130 "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3132 "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3133 "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3134 "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3135 "compare results to Raman Spectroscopy.",
3137 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3138 "require an amount of memory proportional to the total numbers of donors",
3139 "times the total number of acceptors in the selected group(s)."
3142 static real acut=30, abin=1, rcut=0.35, r2cut=0, rbin=0.005, rshell=-1;
3143 static real maxnhb=0,fit_start=1,fit_end=60,temp=298.15,smooth_tail_start=-1, D=-1;
3144 static gmx_bool bNitAcc=TRUE,bDA=TRUE,bMerge=TRUE;
3145 static int nDump=0, nFitPoints=100;
3146 static int nThreads = 0, nBalExp=4;
3148 static gmx_bool bContact=FALSE, bBallistic=FALSE, bBallisticDt=FALSE, bGemFit=FALSE;
3149 static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
3150 static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3154 { "-a", FALSE, etREAL, {&acut},
3155 "Cutoff angle (degrees, Acceptor - Donor - Hydrogen)" },
3156 { "-r", FALSE, etREAL, {&rcut},
3157 "Cutoff radius (nm, X - Acceptor, see next option)" },
3158 { "-da", FALSE, etBOOL, {&bDA},
3159 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3160 { "-r2", FALSE, etREAL, {&r2cut},
3161 "Second cutoff radius. Mainly useful with -contact and -ac"},
3162 { "-abin", FALSE, etREAL, {&abin},
3163 "Binwidth angle distribution (degrees)" },
3164 { "-rbin", FALSE, etREAL, {&rbin},
3165 "Binwidth distance distribution (nm)" },
3166 { "-nitacc",FALSE, etBOOL, {&bNitAcc},
3167 "Regard nitrogen atoms as acceptors" },
3168 { "-contact",FALSE,etBOOL, {&bContact},
3169 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3170 { "-shell", FALSE, etREAL, {&rshell},
3171 "when > 0, only calculate hydrogen bonds within # nm shell around "
3173 { "-fitstart", FALSE, etREAL, {&fit_start},
3174 "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 -gemfit we suggest -fitstart 0" },
3175 { "-fitstart", FALSE, etREAL, {&fit_start},
3176 "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 -gemfit)" },
3177 { "-temp", FALSE, etREAL, {&temp},
3178 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3179 { "-smooth",FALSE, etREAL, {&smooth_tail_start},
3180 "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/tau)" },
3181 { "-dump", FALSE, etINT, {&nDump},
3182 "Dump the first N hydrogen bond ACFs in a single xvg file for debugging" },
3183 { "-max_hb",FALSE, etREAL, {&maxnhb},
3184 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3185 { "-merge", FALSE, etBOOL, {&bMerge},
3186 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3187 { "-geminate", FALSE, etENUM, {gemType},
3188 "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3189 { "-diff", FALSE, etREAL, {&D},
3190 "Dffusion coefficient to use in the rev. gem. recomb. kinetic model. If non-positive, then it will be fitted to the ACF along with ka and kd."},
3192 { "-nthreads", FALSE, etINT, {&nThreads},
3193 "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
3196 const char *bugs[] = {
3197 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3200 { efTRX, "-f", NULL, ffREAD },
3201 { efTPX, NULL, NULL, ffREAD },
3202 { efNDX, NULL, NULL, ffOPTRD },
3203 /* { efNDX, "-sel", "select", ffOPTRD },*/
3204 { efXVG, "-num", "hbnum", ffWRITE },
3205 { efLOG, "-g", "hbond", ffOPTWR },
3206 { efXVG, "-ac", "hbac", ffOPTWR },
3207 { efXVG, "-dist","hbdist", ffOPTWR },
3208 { efXVG, "-ang", "hbang", ffOPTWR },
3209 { efXVG, "-hx", "hbhelix",ffOPTWR },
3210 { efNDX, "-hbn", "hbond", ffOPTWR },
3211 { efXPM, "-hbm", "hbmap", ffOPTWR },
3212 { efXVG, "-don", "donor", ffOPTWR },
3213 { efXVG, "-dan", "danum", ffOPTWR },
3214 { efXVG, "-life","hblife", ffOPTWR },
3215 { efXVG, "-nhbdist", "nhbdist",ffOPTWR }
3218 #define NFILE asize(fnm)
3220 char hbmap [HB_NR]={ ' ', 'o', '-', '*' };
3221 const char *hbdesc[HB_NR]={ "None", "Present", "Inserted", "Present & Inserted" };
3222 t_rgb hbrgb [HB_NR]={ {1,1,1},{1,0,0}, {0,0,1}, {1,0,1} };
3224 t_trxstatus *status;
3229 int npargs,natoms,nframes=0,shatom;
3235 real t,ccut,dist=0.0,ang=0.0;
3236 double max_nhb,aver_nhb,aver_dist;
3237 int h=0,i,j,k=0,l,start,end,id,ja,ogrp,nsel;
3239 int xj,yj,zj,aj,xjj,yjj,zjj;
3240 int xk,yk,zk,ak,xkk,ykk,zkk;
3241 gmx_bool bSelected,bHBmap,bStop,bTwo,was,bBox,bTric;
3243 int grp,nabin,nrbin,bin,resdist,ihb;
3246 FILE *fp,*fpins=NULL,*fpnhb=NULL;
3248 t_ncell *icell,*jcell,*kcell;
3250 unsigned char *datable;
3255 int ii, jj, hh, actual_nThreads;
3257 gmx_bool bGem, bNN, bParallel;
3258 t_gemParams *params=NULL;
3260 CopyRight(stdout,argv[0]);
3263 ppa = add_acf_pargs(&npargs,pa);
3265 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,NFILE,fnm,npargs,
3266 ppa,asize(desc),desc,asize(bugs),bugs,&oenv);
3268 /* NN-loop? If so, what estimator to use ?*/
3270 /* Outcommented for now DvdS 2010-07-13
3271 while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3274 gmx_fatal(FARGS, "Invalid NN-loop type.");
3277 for (i=2; bNN==FALSE && i<NN_NR; i++)
3278 bNN = bNN || NN == i;
3280 if (NN > NN_NONE && bMerge)
3283 /* geminate recombination? If so, which flavor? */
3285 while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode])!=0)
3287 if (gemmode == gemNR)
3288 gmx_fatal(FARGS, "Invalid recombination type.");
3291 for (i=2; bGem==FALSE && i<gemNR; i++)
3292 bGem = bGem || gemmode == i;
3295 printf("Geminate recombination: %s\n" ,gemType[gemmode]);
3297 printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
3300 if (gemmode != gemDD) {
3301 printf("Turning off -contact option...\n");
3305 if (gemmode == gemDD) {
3306 printf("Turning on -contact option...\n");
3311 if (gemmode == gemAA) {
3312 printf("Turning off -merge option...\n");
3316 if (gemmode != gemAA) {
3317 printf("Turning on -merge option...\n");
3325 ccut = cos(acut*DEG2RAD);
3329 gmx_fatal(FARGS,"Can not analyze selected contacts.");
3331 gmx_fatal(FARGS,"Can not analyze contact between H and A: turn off -noda");
3336 /* Don't pollute stdout with information about external libraries.
3338 * printf("NO GSL! Can't find and take away ballistic term in ACF without GSL\n.");
3342 /* Initiate main data structure! */
3343 bHBmap = (opt2bSet("-ac",NFILE,fnm) ||
3344 opt2bSet("-life",NFILE,fnm) ||
3345 opt2bSet("-hbn",NFILE,fnm) ||
3346 opt2bSet("-hbm",NFILE,fnm) ||
3350 /* Same thing here. There is no reason whatsoever to write the specific version of
3351 * OpenMP used for compilation to stdout for normal usage.
3353 * printf("Compiled with OpenMP (%i)\n", _OPENMP);
3357 /* if (bContact && bGem) */
3358 /* gmx_fatal(FARGS, "Can't do reversible geminate recombination with -contact yet."); */
3360 if (opt2bSet("-nhbdist",NFILE,fnm)) {
3361 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3362 fpnhb = xvgropen(opt2fn("-nhbdist",NFILE,fnm),
3363 "Number of donor-H with N HBs","Time (ps)","N",oenv);
3364 xvgr_legend(fpnhb,asize(leg),leg,oenv);
3367 hb = mk_hbdata(bHBmap,opt2bSet("-dan",NFILE,fnm),bMerge || bContact, bGem, gemmode);
3370 read_tpx_top(ftp2fn(efTPX,NFILE,fnm),&ir,box,&natoms,NULL,NULL,NULL,&top);
3372 snew(grpnames,grNR);
3375 /* Make Donor-Acceptor table */
3376 snew(datable, top.atoms.nr);
3377 gen_datable(index[0],isize[0],datable,top.atoms.nr);
3380 /* analyze selected hydrogen bonds */
3381 printf("Select group with selected atoms:\n");
3382 get_index(&(top.atoms),opt2fn("-sel",NFILE,fnm),
3383 1,&nsel,index,grpnames);
3385 gmx_fatal(FARGS,"Number of atoms in group '%s' not a multiple of 3\n"
3386 "and therefore cannot contain triplets of "
3387 "Donor-Hydrogen-Acceptor",grpnames[0]);
3390 for(i=0; (i<nsel); i+=3) {
3391 int dd = index[0][i];
3392 int aa = index[0][i+2];
3393 /* int */ hh = index[0][i+1];
3394 add_dh (&hb->d,dd,hh,i,datable);
3395 add_acc(&hb->a,aa,i);
3396 /* Should this be here ? */
3397 snew(hb->d.dptr,top.atoms.nr);
3398 snew(hb->a.aptr,top.atoms.nr);
3399 add_hbond(hb,dd,aa,hh,gr0,gr0,0,bMerge,0,bContact,peri);
3401 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
3402 isize[0],grpnames[0]);
3405 /* analyze all hydrogen bonds: get group(s) */
3406 printf("Specify 2 groups to analyze:\n");
3407 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
3408 2,isize,index,grpnames);
3410 /* check if we have two identical or two non-overlapping groups */
3411 bTwo = isize[0] != isize[1];
3412 for (i=0; (i<isize[0]) && !bTwo; i++)
3413 bTwo = index[0][i] != index[1][i];
3415 printf("Checking for overlap in atoms between %s and %s\n",
3416 grpnames[0],grpnames[1]);
3417 for (i=0; i<isize[1];i++)
3418 if (ISINGRP(datable[index[1][i]]))
3419 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3420 grpnames[0],grpnames[1]);
3422 printf("Checking for overlap in atoms between %s and %s\n",
3423 grpnames[0],grpnames[1]);
3424 for (i=0; i<isize[0]; i++)
3425 for (j=0; j<isize[1]; j++)
3426 if (index[0][i] == index[1][j])
3427 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3428 grpnames[0],grpnames[1]);
3432 printf("Calculating %s "
3433 "between %s (%d atoms) and %s (%d atoms)\n",
3434 bContact ? "contacts" : "hydrogen bonds",
3435 grpnames[0],isize[0],grpnames[1],isize[1]);
3437 fprintf(stderr,"Calculating %s in %s (%d atoms)\n",
3438 bContact?"contacts":"hydrogen bonds",grpnames[0],isize[0]);
3442 /* search donors and acceptors in groups */
3443 snew(datable, top.atoms.nr);
3444 for (i=0; (i<grNR); i++)
3445 if ( ((i==gr0) && !bSelected ) ||
3446 ((i==gr1) && bTwo )) {
3447 gen_datable(index[i],isize[i],datable,top.atoms.nr);
3449 search_acceptors(&top,isize[i],index[i],&hb->a,i,
3450 bNitAcc,TRUE,(bTwo && (i==gr0)) || !bTwo, datable);
3451 search_donors (&top,isize[i],index[i],&hb->d,i,
3452 TRUE,(bTwo && (i==gr1)) || !bTwo, datable);
3455 search_acceptors(&top,isize[i],index[i],&hb->a,i,bNitAcc,FALSE,TRUE, datable);
3456 search_donors (&top,isize[i],index[i],&hb->d,i,FALSE,TRUE, datable);
3459 clear_datable_grp(datable,top.atoms.nr);
3462 printf("Found %d donors and %d acceptors\n",hb->d.nrd,hb->a.nra);
3464 snew(donors[gr0D], dons[gr0D].nrd);*/
3467 printf("Making hbmap structure...");
3468 /* Generate hbond data structure */
3473 #ifdef HAVE_NN_LOOPS
3479 printf("Making per structure...");
3480 /* Generate hbond data structure */
3487 if (hb->d.nrd + hb->a.nra == 0) {
3488 printf("No Donors or Acceptors found\n");
3492 if (hb->d.nrd == 0) {
3493 printf("No Donors found\n");
3496 if (hb->a.nra == 0) {
3497 printf("No Acceptors found\n");
3502 gmx_fatal(FARGS,"Nothing to be done");
3509 /* get index group with atom for shell */
3511 printf("Select atom for shell (1 atom):\n");
3512 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
3513 1,&shisz,&shidx,&shgrpnm);
3515 printf("group contains %d atoms, should be 1 (one)\n",shisz);
3516 } while(shisz != 1);
3518 printf("Will calculate hydrogen bonds within a shell "
3519 "of %g nm around atom %i\n",rshell,shatom+1);
3522 /* Analyze trajectory */
3523 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
3524 if ( natoms > top.atoms.nr )
3525 gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
3526 top.atoms.nr,natoms);
3528 bBox = ir.ePBC!=epbcNONE;
3529 grid = init_grid(bBox, box, (rcut>r2cut)?rcut:r2cut, ngrid);
3532 snew(adist,nabin+1);
3533 snew(rdist,nrbin+1);
3536 gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
3541 #define __ADIST adist
3542 #define __RDIST rdist
3544 #else /* HAVE_OPENMP ================================================== \
3545 * Set up the OpenMP stuff, |
3546 * like the number of threads and such |
3547 * Also start the parallel loop. |
3549 #define __ADIST p_adist[threadNr]
3550 #define __RDIST p_rdist[threadNr]
3551 #define __HBDATA p_hb[threadNr]
3553 bParallel = !bSelected;
3557 #if (_OPENMP > 200805)
3558 actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit());
3560 actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
3562 omp_set_num_threads(actual_nThreads);
3563 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
3568 actual_nThreads = 1;
3571 t_hbdata **p_hb; /* one per thread, then merge after the frame loop */
3572 int **p_adist, **p_rdist; /* a histogram for each thread. */
3573 snew(p_hb, actual_nThreads);
3574 snew(p_adist, actual_nThreads);
3575 snew(p_rdist, actual_nThreads);
3576 for (i=0; i<actual_nThreads; i++)
3579 snew(p_adist[i], nabin+1);
3580 snew(p_rdist[i], nrbin+1);
3582 p_hb[i]->max_frames = 0;
3583 p_hb[i]->nhb = NULL;
3584 p_hb[i]->ndist = NULL;
3585 p_hb[i]->n_bound = NULL;
3586 p_hb[i]->time = NULL;
3587 p_hb[i]->nhx = NULL;
3589 p_hb[i]->bHBmap = hb->bHBmap;
3590 p_hb[i]->bDAnr = hb->bDAnr;
3591 p_hb[i]->bGem = hb->bGem;
3592 p_hb[i]->wordlen = hb->wordlen;
3593 p_hb[i]->nframes = hb->nframes;
3594 p_hb[i]->maxhydro = hb->maxhydro;
3595 p_hb[i]->danr = hb->danr;
3598 p_hb[i]->hbmap = hb->hbmap;
3599 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
3600 p_hb[i]->per = hb->per;
3602 #ifdef HAVE_NN_LOOPS
3603 p_hb[i]->hbE = hb->hbE;
3607 p_hb[i]->nrdist = 0;
3610 /* Make a thread pool here,
3611 * instead of forking anew at every frame. */
3613 #pragma omp parallel \
3614 private(i, j, h, ii, jj, hh, E, \
3615 xi, yi, zi, xj, yj, zj, threadNr, \
3616 dist, ang, peri, icell, jcell, \
3617 grp, ogrp, ai, aj, xjj, yjj, zjj, \
3618 xk, yk, zk, ihb, id, resdist, \
3619 xkk, ykk, zkk, kcell, ak, k, bTric) \
3621 shared(hb, p_hb, p_adist, p_rdist, actual_nThreads, \
3622 x, bBox, box, hbox, rcut, r2cut, rshell, \
3623 shatom, ngrid, grid, nframes, t, \
3624 bParallel, bNN, index, bMerge, bContact, \
3625 bTwo, bDA,ccut, abin, rbin, top, \
3626 bSelected, bDebug, stderr, nsel, \
3627 bGem, oenv, fnm, fpnhb, trrStatus, natoms, \
3628 status, nabin, nrbin, adist, rdist, debug)
3629 { /* Start of parallel region */
3630 threadNr = omp_get_thread_num();
3631 #endif /* HAVE_OPENMP ================================================= */
3634 bTric = bBox && TRICLINIC(box);
3637 sync_hbdata(hb, p_hb[threadNr], nframes, t);
3641 build_grid(hb,x,x[shatom], bBox,box,hbox, (rcut>r2cut)?rcut:r2cut,
3642 rshell, ngrid,grid);
3643 reset_nhbonds(&(hb->d));
3645 if (debug && bDebug)
3646 dump_grid(debug, ngrid, grid);
3648 add_frames(hb,nframes);
3649 init_hbframe(hb,nframes,t);
3652 count_da_grid(ngrid, grid, hb->danr[nframes]);
3656 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
3660 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
3661 /* Loop over all atom pairs and estimate interaction energy */
3662 #ifdef HAVE_OPENMP /* ------- */
3664 #endif /* HAVE_OPENMP ------- */
3666 addFramesNN(hb, nframes);
3668 #ifdef HAVE_OPENMP /* ---------------- */
3670 #pragma omp for schedule(dynamic)
3671 #endif /* HAVE_OPENMP ---------------- */
3672 for (i=0; i<hb->d.nrd; i++)
3674 for(j=0;j<hb->a.nra; j++)
3677 h < (bContact ? 1 : hb->d.nhydro[i]);
3680 if (i==hb->d.nrd || j==hb->a.nra)
3681 gmx_fatal(FARGS, "out of bounds");
3683 /* Get the real atom ids */
3686 hh = hb->d.hydro[i][h];
3688 /* Estimate the energy from the geometry */
3689 E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
3690 /* Store the energy */
3691 storeHbEnergy(hb, i, j, h, E, nframes);
3695 #endif /* HAVE_NN_LOOPS */
3705 /* Do not parallelize this just yet. */
3707 for(ii=0; (ii<nsel); ii++) {
3708 int dd = index[0][i];
3709 int aa = index[0][i+2];
3710 /* int */ hh = index[0][i+1];
3711 ihb = is_hbond(hb,ii,ii,dd,aa,rcut,r2cut,ccut,x,bBox,box,
3712 hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
3715 /* add to index if not already there */
3717 add_hbond(hb,dd,aa,hh,ii,ii,nframes,bMerge,ihb,bContact,peri);
3721 } /* if (bSelected) */
3729 calcBoxProjection(box, hb->per->P);
3731 /* loop over all gridcells (xi,yi,zi) */
3732 /* Removed confusing macro, DvdS 27/12/98 */
3735 /* The outer grid loop will have to do for now. */
3736 #pragma omp for schedule(dynamic)
3738 for(xi=0; (xi<ngrid[XX]); xi++)
3739 for(yi=0; (yi<ngrid[YY]); yi++)
3740 for(zi=0; (zi<ngrid[ZZ]); zi++) {
3742 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
3743 for (grp=gr0; (grp <= (bTwo?gr1:gr0)); grp++) {
3744 icell=&(grid[zi][yi][xi].d[grp]);
3751 /* loop over all hydrogen atoms from group (grp)
3752 * in this gridcell (icell)
3754 for (ai=0; (ai<icell->nr); ai++) {
3755 i = icell->atoms[ai];
3757 /* loop over all adjacent gridcells (xj,yj,zj) */
3758 /* This is a macro!!! */
3759 LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
3760 jcell=&(grid[zj][yj][xj].a[ogrp]);
3761 /* loop over acceptor atoms from other group (ogrp)
3762 * in this adjacent gridcell (jcell)
3764 for (aj=0; (aj<jcell->nr); aj++) {
3765 j = jcell->atoms[aj];
3767 /* check if this once was a h-bond */
3769 ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
3770 hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
3773 /* add to index if not already there */
3775 add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
3777 /* make angle and distance distributions */
3778 if (ihb == hbHB && !bContact) {
3780 gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
3782 __ADIST[(int)( ang/abin)]++;
3783 __RDIST[(int)(dist/rbin)]++;
3786 if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
3787 gmx_fatal(FARGS,"Invalid donor %d",i);
3788 if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
3789 gmx_fatal(FARGS,"Invalid acceptor %d",j);
3790 resdist=abs(top.atoms.atom[i].resind-
3791 top.atoms.atom[j].resind);
3792 if (resdist >= max_hx)
3794 __HBDATA->nhx[nframes][resdist]++;
3804 } /* for xi,yi,zi */
3805 } /* if (bSelected) {...} else */
3807 #ifdef HAVE_OPENMP /* ---------------------------- */
3808 /* Better wait for all threads to finnish using x[] before updating it. */
3810 #pragma omp barrier /* */
3811 #pragma omp critical /* */
3813 /* Sum up histograms and counts from p_hb[] into hb */
3815 hb->nhb[k] += p_hb[threadNr]->nhb[k];
3816 hb->ndist[k] += p_hb[threadNr]->ndist[k];
3817 for (j=0; j<max_hx; j++) /**/
3818 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
3822 /* Here are a handful of single constructs
3823 * to share the workload a bit. The most
3824 * important one is of course the last one,
3825 * where there's a potential bottleneck in form
3827 #pragma omp single /* ++++++++++++++++, */
3828 #endif /* HAVE_OPENMP ----------------+------------*/
3830 if (hb != NULL) /* */
3832 analyse_donor_props(opt2fn_null("-don",NFILE,fnm),hb,k,t,oenv);
3835 #ifdef HAVE_OPENMP /* + */
3836 #pragma omp single /* +++ +++ */
3840 do_nhb_dist(fpnhb,hb,t);
3842 } /* if (bNN) {...} else + */
3843 #ifdef HAVE_OPENMP /* + */
3844 #pragma omp single /* +++ +++ */
3847 trrStatus = (read_next_x(oenv,status,&t,natoms,x,box));
3850 #ifdef HAVE_OPENMP /* ++++++++++++++++� */
3853 } while (trrStatus);
3856 #pragma omp critical
3858 hb->nrhb += p_hb[threadNr]->nrhb;
3859 hb->nrdist += p_hb[threadNr]->nrdist;
3861 /* Free parallel datastructures */
3862 sfree(p_hb[threadNr]->nhb);
3863 sfree(p_hb[threadNr]->ndist);
3864 sfree(p_hb[threadNr]->nhx);
3867 for (i=0; i<nabin; i++)
3868 for (j=0; j<actual_nThreads; j++)
3870 adist[i] += p_adist[j][i];
3872 for (i=0; i<=nrbin; i++)
3873 for (j=0; j<actual_nThreads; j++)
3874 rdist[i] += p_rdist[j][i];
3876 sfree(p_adist[threadNr]);
3877 sfree(p_rdist[threadNr]);
3878 } /* End of parallel region */
3883 if(nframes <2 && (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm)))
3885 gmx_fatal(FARGS,"Cannot calculate autocorrelation of life times with less than two frames");
3888 free_grid(ngrid,&grid);
3894 /* Compute maximum possible number of different hbonds */
3898 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
3900 /* Added support for -contact below.
3901 * - Erik Marklund, May 29-31, 2006 */
3902 /* Changed contact code.
3903 * - Erik Marklund, June 29, 2006 */
3904 if (bHBmap && !bNN) {
3906 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
3908 printf("Found %d different %s in trajectory\n"
3909 "Found %d different atom-pairs within %s distance\n",
3910 hb->nrhb, bContact?"contacts":"hydrogen bonds",
3911 hb->nrdist,(r2cut>0)?"second cut-off":"hydrogen bonding");
3913 /*Control the pHist.*/
3916 merge_hb(hb,bTwo,bContact);
3918 if (opt2bSet("-hbn",NFILE,fnm))
3919 dump_hbmap(hb,NFILE,fnm,bTwo,bContact,isize,index,grpnames,&top.atoms);
3921 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3922 * to make the -hbn and -hmb output match eachother.
3923 * - Erik Marklund, May 30, 2006 */
3926 /* Print out number of hbonds and distances */
3929 fp = xvgropen(opt2fn("-num",NFILE,fnm),bContact ? "Contacts" :
3930 "Hydrogen Bonds","Time","Number",oenv);
3932 snew(leg[0],STRLEN);
3933 snew(leg[1],STRLEN);
3934 sprintf(leg[0],"%s",bContact?"Contacts":"Hydrogen bonds");
3935 sprintf(leg[1],"Pairs within %g nm",(r2cut>0)?r2cut:rcut);
3936 xvgr_legend(fp,2,(const char**)leg,oenv);
3940 for(i=0; (i<nframes); i++) {
3941 fprintf(fp,"%10g %10d %10d\n",hb->time[i],hb->nhb[i],hb->ndist[i]);
3942 aver_nhb += hb->nhb[i];
3943 aver_dist += hb->ndist[i];
3946 aver_nhb /= nframes;
3947 aver_dist /= nframes;
3948 /* Print HB distance distribution */
3949 if (opt2bSet("-dist",NFILE,fnm)) {
3953 for(i=0; i<nrbin; i++)
3956 fp = xvgropen(opt2fn("-dist",NFILE,fnm),
3957 "Hydrogen Bond Distribution",
3958 "Hydrogen - Acceptor Distance (nm)","",oenv);
3959 for(i=0; i<nrbin; i++)
3960 fprintf(fp,"%10g %10g\n",(i+0.5)*rbin,rdist[i]/(rbin*(real)sum));
3964 /* Print HB angle distribution */
3965 if (opt2bSet("-ang",NFILE,fnm)) {
3969 for(i=0; i<nabin; i++)
3972 fp = xvgropen(opt2fn("-ang",NFILE,fnm),
3973 "Hydrogen Bond Distribution",
3974 "Donor - Hydrogen - Acceptor Angle (\\SO\\N)","",oenv);
3975 for(i=0; i<nabin; i++)
3976 fprintf(fp,"%10g %10g\n",(i+0.5)*abin,adist[i]/(abin*(real)sum));
3980 /* Print HB in alpha-helix */
3981 if (opt2bSet("-hx",NFILE,fnm)) {
3982 fp = xvgropen(opt2fn("-hx",NFILE,fnm),
3983 "Hydrogen Bonds","Time(ps)","Count",oenv);
3984 xvgr_legend(fp,NRHXTYPES,hxtypenames,oenv);
3985 for(i=0; i<nframes; i++) {
3986 fprintf(fp,"%10g",hb->time[i]);
3987 for (j=0; j<max_hx; j++)
3988 fprintf(fp," %6d",hb->nhx[i][j]);
3994 printf("Average number of %s per timeframe %.3f out of %g possible\n",
3995 bContact ? "contacts" : "hbonds",
3996 bContact ? aver_dist : aver_nhb, max_nhb);
3998 /* Do Autocorrelation etc. */
4001 Added support for -contact in ac and hbm calculations below.
4002 - Erik Marklund, May 29, 2006
4006 if (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm))
4007 please_cite(stdout,"Spoel2006b");
4008 if (opt2bSet("-ac",NFILE,fnm)) {
4009 char *gemstring=NULL;
4012 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
4013 gemBallistic, nBalExp, bBallisticDt);
4015 gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
4017 gemstring = strdup(gemType[hb->per->gemtype]);
4018 do_hbac(opt2fn("-ac",NFILE,fnm),hb,nDump,
4019 bMerge,bContact,fit_start,temp,r2cut>0,smooth_tail_start,oenv,
4020 params, gemstring, nThreads, NN, bBallistic, bGemFit);
4022 if (opt2bSet("-life",NFILE,fnm))
4023 do_hblife(opt2fn("-life",NFILE,fnm),hb,bMerge,bContact,oenv);
4024 if (opt2bSet("-hbm",NFILE,fnm)) {
4029 mat.ny=(bContact ? hb->nrdist : hb->nrhb);
4031 snew(mat.matrix,mat.nx);
4032 for(x=0; (x<mat.nx); x++)
4033 snew(mat.matrix[x],mat.ny);
4035 for(id=0; (id<hb->d.nrd); id++)
4036 for(ia=0; (ia<hb->a.nra); ia++) {
4037 for(hh=0; (hh<hb->maxhydro); hh++) {
4038 if (hb->hbmap[id][ia]) {
4039 if (ISHB(hb->hbmap[id][ia]->history[hh])) {
4040 /* Changed '<' into '<=' in the for-statement below.
4041 * It fixed the previously undiscovered bug that caused
4042 * the last occurance of an hbond/contact to not be
4043 * set in mat.matrix. Have a look at any old -hbm-output
4044 * and you will notice that the last column is allways empty.
4045 * - Erik Marklund May 30, 2006
4047 for(x=0; (x<=hb->hbmap[id][ia]->nframes); x++) {
4048 int nn0 = hb->hbmap[id][ia]->n0;
4049 range_check(y,0,mat.ny);
4050 mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh],x);
4057 mat.axis_x=hb->time;
4058 snew(mat.axis_y,mat.ny);
4059 for(j=0; j<mat.ny; j++)
4061 sprintf(mat.title,bContact ? "Contact Existence Map":
4062 "Hydrogen Bond Existence Map");
4063 sprintf(mat.legend,bContact ? "Contacts" : "Hydrogen Bonds");
4064 sprintf(mat.label_x,"Time (ps)");
4065 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4068 snew(mat.map,mat.nmap);
4069 for(i=0; i<mat.nmap; i++) {
4070 mat.map[i].code.c1=hbmap[i];
4071 mat.map[i].desc=hbdesc[i];
4072 mat.map[i].rgb=hbrgb[i];
4074 fp = opt2FILE("-hbm",NFILE,fnm,"w");
4075 write_xpm_m(fp, mat);
4077 for(x=0; x<mat.nx; x++)
4078 sfree(mat.matrix[x]);
4086 fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4087 fprintf(stderr, "Freeing pHist for all donors...\n");
4088 for (i=0; i<hb->d.nrd; i++) {
4089 fprintf(stderr, "\r%i",i);
4090 if (hb->per->pHist[i] != NULL) {
4091 for (j=0; j<hb->a.nra; j++)
4092 clearPshift(&(hb->per->pHist[i][j]));
4093 sfree(hb->per->pHist[i]);
4096 sfree(hb->per->pHist);
4097 sfree(hb->per->p2i);
4099 fprintf(stderr, "...done.\n");
4102 #ifdef HAVE_NN_LOOPS
4112 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4114 fp = xvgropen(opt2fn("-dan",NFILE,fnm),
4115 "Donors and Acceptors","Time(ps)","Count",oenv);
4116 nleg = (bTwo?2:1)*2;
4117 snew(legnames,nleg);
4119 for(j=0; j<grNR; j++)
4120 if ( USE_THIS_GROUP(j) ) {
4121 sprintf(buf,"Donors %s",grpnames[j]);
4122 legnames[i++] = strdup(buf);
4123 sprintf(buf,"Acceptors %s",grpnames[j]);
4124 legnames[i++] = strdup(buf);
4127 gmx_incons("number of legend entries");
4128 xvgr_legend(fp,nleg,(const char**)legnames,oenv);
4129 for(i=0; i<nframes; i++) {
4130 fprintf(fp,"%10g",hb->time[i]);
4131 for (j=0; (j<grNR); j++)
4132 if ( USE_THIS_GROUP(j) )
4133 fprintf(fp," %6d",hb->danr[i][j]);