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 static char *mkatomname(t_atoms *atoms,int i)
903 rnr = atoms->atom[i].resind;
904 sprintf(buf,"%4s%d%-4s",
905 *atoms->resinfo[rnr].name,atoms->resinfo[rnr].nr,*atoms->atomname[i]);
910 static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms){
911 /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
914 for (i=0; i<isize; i++){
915 if (index[i] >= natoms)
916 gmx_fatal(FARGS,"Atom has index %d larger than number of atoms %d.",index[i],natoms);
917 datable[index[i]] |= INGROUP;
921 static void clear_datable_grp(unsigned char *datable, int size){
922 /* Clears group information from the table */
924 const char mask = !(char)INGROUP;
930 static void add_acc(t_acceptors *a,int ia,int grp)
932 if (a->nra >= a->max_nra) {
934 srenew(a->acc,a->max_nra);
935 srenew(a->grp,a->max_nra);
937 a->grp[a->nra] = grp;
938 a->acc[a->nra++] = ia;
941 static void search_acceptors(t_topology *top,int isize,
942 atom_id *index,t_acceptors *a,int grp,
944 gmx_bool bContact,gmx_bool bDoIt, unsigned char *datable)
949 for (i=0; (i<isize); i++) {
952 (((*top->atoms.atomname[n])[0] == 'O') ||
953 (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
954 ISINGRP(datable[n])) {
955 datable[n] |= ACC; /* set the atom's acceptor flag in datable. */
960 snew(a->aptr,top->atoms.nr);
961 for(i=0; (i<top->atoms.nr); i++)
963 for(i=0; (i<a->nra); i++)
964 a->aptr[a->acc[i]] = i;
967 static void add_h2d(int id,int ih,t_donors *ddd)
971 for(i=0; (i<ddd->nhydro[id]); i++)
972 if (ddd->hydro[id][i] == ih) {
973 printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
977 if (i == ddd->nhydro[id]) {
978 if (ddd->nhydro[id] >= MAXHYDRO)
979 gmx_fatal(FARGS,"Donor %d has more than %d hydrogens!",
980 ddd->don[id],MAXHYDRO);
981 ddd->hydro[id][i] = ih;
986 static void add_dh(t_donors *ddd,int id,int ih,int grp, unsigned char *datable)
990 if (ISDON(datable[id]) || !datable) {
991 if (ddd->dptr[id] == NOTSET) { /* New donor */
998 if (ddd->nrd >= ddd->max_nrd) {
1000 srenew(ddd->don,ddd->max_nrd);
1001 srenew(ddd->nhydro,ddd->max_nrd);
1002 srenew(ddd->hydro,ddd->max_nrd);
1003 srenew(ddd->nhbonds,ddd->max_nrd);
1004 srenew(ddd->grp,ddd->max_nrd);
1006 ddd->don[ddd->nrd] = id;
1007 ddd->nhydro[ddd->nrd] = 0;
1008 ddd->grp[ddd->nrd] = grp;
1015 printf("Warning: Atom %d is not in the d/a-table!\n", id);
1018 static void search_donors(t_topology *top, int isize, atom_id *index,
1019 t_donors *ddd,int grp,gmx_bool bContact,gmx_bool bDoIt,
1020 unsigned char *datable)
1023 t_functype func_type;
1024 t_ilist *interaction;
1025 atom_id nr1,nr2,nr3;
1029 snew(ddd->dptr,top->atoms.nr);
1030 for(i=0; (i<top->atoms.nr); i++)
1031 ddd->dptr[i] = NOTSET;
1036 for(i=0; (i<isize); i++) {
1037 datable[index[i]] |= DON;
1038 add_dh(ddd,index[i],-1,grp,datable);
1042 for(func_type=0; (func_type < F_NRE); func_type++) {
1043 interaction=&(top->idef.il[func_type]);
1044 if (func_type == F_POSRES)
1045 { /* The ilist looks strange for posre. Bug in grompp?
1046 * We don't need posre interactions for hbonds anyway.*/
1049 for(i=0; i < interaction->nr;
1050 i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
1052 if (func_type != top->idef.functype[interaction->iatoms[i]]) {
1053 fprintf(stderr,"Error in func_type %s",
1054 interaction_function[func_type].longname);
1058 /* check out this functype */
1059 if (func_type == F_SETTLE) {
1060 nr1 = interaction->iatoms[i+1];
1061 nr2 = interaction->iatoms[i+2];
1062 nr3 = interaction->iatoms[i+3];
1064 if (ISINGRP(datable[nr1])) {
1065 if (ISINGRP(datable[nr2])) {
1066 datable[nr1] |= DON;
1067 add_dh(ddd,nr1,nr1+1,grp,datable);
1069 if (ISINGRP(datable[nr3])) {
1070 datable[nr1] |= DON;
1071 add_dh(ddd,nr1,nr1+2,grp,datable);
1075 else if (IS_CHEMBOND(func_type)) {
1076 for (j=0; j<2; j++) {
1077 nr1=interaction->iatoms[i+1+j];
1078 nr2=interaction->iatoms[i+2-j];
1079 if ((*top->atoms.atomname[nr1][0] == 'H') &&
1080 ((*top->atoms.atomname[nr2][0] == 'O') ||
1081 (*top->atoms.atomname[nr2][0] == 'N')) &&
1082 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
1083 datable[nr2] |= DON;
1084 add_dh(ddd,nr2,nr1,grp,datable);
1091 for(func_type=0; func_type < F_NRE; func_type++) {
1092 interaction=&top->idef.il[func_type];
1093 for(i=0; i < interaction->nr;
1094 i+=interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) {
1096 if (func_type != top->idef.functype[interaction->iatoms[i]])
1097 gmx_incons("function type in search_donors");
1099 if ( interaction_function[func_type].flags & IF_VSITE ) {
1100 nr1=interaction->iatoms[i+1];
1101 if ( *top->atoms.atomname[nr1][0] == 'H') {
1104 while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
1109 if ( !stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
1110 ( *top->atoms.atomname[nr2][0] == 'N') ) &&
1111 ISINGRP(datable[nr1]) && ISINGRP(datable[nr2])) {
1112 datable[nr2] |= DON;
1113 add_dh(ddd,nr2,nr1,grp,datable);
1123 static t_gridcell ***init_grid(gmx_bool bBox,rvec box[],real rcut,ivec ngrid)
1129 for(i=0; i<DIM; i++)
1130 ngrid[i]=(box[i][i]/(1.2*rcut));
1132 if ( !bBox || (ngrid[XX]<3) || (ngrid[YY]<3) || (ngrid[ZZ]<3) )
1133 for(i=0; i<DIM; i++)
1136 printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n",
1137 ngrid[XX],ngrid[YY],ngrid[ZZ],rcut);
1138 snew(grid,ngrid[ZZ]);
1139 for (z=0; z<ngrid[ZZ]; z++) {
1140 snew((grid)[z],ngrid[YY]);
1141 for (y=0; y<ngrid[YY]; y++)
1142 snew((grid)[z][y],ngrid[XX]);
1147 static void reset_nhbonds(t_donors *ddd)
1151 for(i=0; (i<ddd->nrd); i++)
1152 for(j=0; (j<MAXHH); j++)
1153 ddd->nhbonds[i][j] = 0;
1156 void pbc_correct_gem(rvec dx,matrix box,rvec hbox);
1158 static void build_grid(t_hbdata *hb,rvec x[], rvec xshell,
1159 gmx_bool bBox, matrix box, rvec hbox,
1160 real rcut, real rshell,
1161 ivec ngrid, t_gridcell ***grid)
1163 int i,m,gr,xi,yi,zi,nr;
1166 rvec invdelta,dshell,xtemp={0,0,0};
1168 gmx_bool bDoRshell,bInShell,bAcc;
1173 bDoRshell = (rshell > 0);
1174 rshell2 = sqr(rshell);
1177 #define DBB(x) if (debug && bDebug) fprintf(debug,"build_grid, line %d, %s = %d\n",__LINE__,#x,x)
1179 for(m=0; m<DIM; m++) {
1180 hbox[m]=box[m][m]*0.5;
1182 invdelta[m]=ngrid[m]/box[m][m];
1183 if (1/invdelta[m] < rcut)
1184 gmx_fatal(FARGS,"Your computational box has shrunk too much.\n"
1185 "%s can not handle this situation, sorry.\n",
1194 /* resetting atom counts */
1195 for(gr=0; (gr<grNR); gr++) {
1196 for (zi=0; zi<ngrid[ZZ]; zi++)
1197 for (yi=0; yi<ngrid[YY]; yi++)
1198 for (xi=0; xi<ngrid[XX]; xi++) {
1199 grid[zi][yi][xi].d[gr].nr=0;
1200 grid[zi][yi][xi].a[gr].nr=0;
1204 /* put atoms in grid cells */
1205 for(bAcc=FALSE; (bAcc<=TRUE); bAcc++) {
1215 for(i=0; (i<nr); i++) {
1216 /* check if we are inside the shell */
1217 /* if bDoRshell=FALSE then bInShell=TRUE always */
1221 rvec_sub(x[ad[i]],xshell,dshell);
1223 if (FALSE && !hb->bGem) {
1224 for(m=DIM-1; m>=0 && bInShell; m--) {
1225 if ( dshell[m] < -hbox[m] )
1226 rvec_inc(dshell,box[m]);
1227 else if ( dshell[m] >= hbox[m] )
1228 dshell[m] -= 2*hbox[m];
1229 /* if we're outside the cube, we're outside the sphere also! */
1230 if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
1234 gmx_bool bDone = FALSE;
1238 for(m=DIM-1; m>=0 && bInShell; m--) {
1239 if ( dshell[m] < -hbox[m] ) {
1241 rvec_inc(dshell,box[m]);
1243 if ( dshell[m] >= hbox[m] ) {
1245 dshell[m] -= 2*hbox[m];
1249 for(m=DIM-1; m>=0 && bInShell; m--) {
1250 /* if we're outside the cube, we're outside the sphere also! */
1251 if ( (dshell[m]>rshell) || (-dshell[m]>rshell) )
1256 /* if we're inside the cube, check if we're inside the sphere */
1258 bInShell = norm2(dshell) < rshell2;
1264 copy_rvec(x[ad[i]], xtemp);
1265 pbc_correct_gem(x[ad[i]], box, hbox);
1267 for(m=DIM-1; m>=0; m--) {
1268 if (TRUE || !hb->bGem){
1269 /* put atom in the box */
1270 while( x[ad[i]][m] < 0 )
1271 rvec_inc(x[ad[i]],box[m]);
1272 while( x[ad[i]][m] >= box[m][m] )
1273 rvec_dec(x[ad[i]],box[m]);
1275 /* determine grid index of atom */
1276 grididx[m]=x[ad[i]][m]*invdelta[m];
1277 grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
1280 copy_rvec(xtemp, x[ad[i]]); /* copy back */
1284 range_check(gx,0,ngrid[XX]);
1285 range_check(gy,0,ngrid[YY]);
1286 range_check(gz,0,ngrid[ZZ]);
1290 /* add atom to grid cell */
1292 newgrid=&(grid[gz][gy][gx].a[gr]);
1294 newgrid=&(grid[gz][gy][gx].d[gr]);
1295 if (newgrid->nr >= newgrid->maxnr) {
1297 DBB(newgrid->maxnr);
1298 srenew(newgrid->atoms, newgrid->maxnr);
1301 newgrid->atoms[newgrid->nr]=ad[i];
1309 static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr)
1313 for(gr=0; (gr<grNR); gr++) {
1315 for (zi=0; zi<ngrid[ZZ]; zi++)
1316 for (yi=0; yi<ngrid[YY]; yi++)
1317 for (xi=0; xi<ngrid[XX]; xi++) {
1318 danr[gr] += grid[zi][yi][xi].d[gr].nr;
1324 * Without a box, the grid is 1x1x1, so all loops are 1 long.
1325 * With a rectangular box (bTric==FALSE) all loops are 3 long.
1326 * With a triclinic box all loops are 3 long, except when a cell is
1327 * located next to one of the box edges which is not parallel to the
1328 * x/y-plane, in that case all cells in a line or layer are searched.
1329 * This could be implemented slightly more efficient, but the code
1330 * would get much more complicated.
1332 static inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1334 return ((n==1) ? x : bTric && bEdge ? 0 : (x-1));
1336 static inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
1338 return ((n==1) ? x : bTric && bEdge ? (n-1) : (x+1));
1340 static inline int grid_mod(int j, int n)
1345 static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
1347 int gr,x,y,z,sum[grNR];
1349 fprintf(fp,"grid %dx%dx%d\n",ngrid[XX],ngrid[YY],ngrid[ZZ]);
1350 for (gr=0; gr<grNR; gr++) {
1352 fprintf(fp,"GROUP %d (%s)\n",gr,grpnames[gr]);
1353 for (z=0; z<ngrid[ZZ]; z+=2) {
1354 fprintf(fp,"Z=%d,%d\n",z,z+1);
1355 for (y=0; y<ngrid[YY]; y++) {
1356 for (x=0; x<ngrid[XX]; x++) {
1357 fprintf(fp,"%3d",grid[x][y][z].d[gr].nr);
1358 sum[gr]+=grid[z][y][x].d[gr].nr;
1359 fprintf(fp,"%3d",grid[x][y][z].a[gr].nr);
1360 sum[gr]+=grid[z][y][x].a[gr].nr;
1364 if ( (z+1) < ngrid[ZZ] )
1365 for (x=0; x<ngrid[XX]; x++) {
1366 fprintf(fp,"%3d",grid[z+1][y][x].d[gr].nr);
1367 sum[gr]+=grid[z+1][y][x].d[gr].nr;
1368 fprintf(fp,"%3d",grid[z+1][y][x].a[gr].nr);
1369 sum[gr]+=grid[z+1][y][x].a[gr].nr;
1375 fprintf(fp,"TOTALS:");
1376 for (gr=0; gr<grNR; gr++)
1377 fprintf(fp," %d=%d",gr,sum[gr]);
1381 /* New GMX record! 5 * in a row. Congratulations!
1382 * Sorry, only four left.
1384 static void free_grid(ivec ngrid, t_gridcell ****grid)
1387 t_gridcell ***g = *grid;
1389 for (z=0; z<ngrid[ZZ]; z++) {
1390 for (y=0; y<ngrid[YY]; y++) {
1399 void pbc_correct_gem(rvec dx,matrix box,rvec hbox)
1402 gmx_bool bDone = FALSE;
1405 for(m=DIM-1; m>=0; m--) {
1406 if ( dx[m] < -hbox[m] ) {
1408 rvec_inc(dx,box[m]);
1410 if ( dx[m] >= hbox[m] ) {
1412 rvec_dec(dx,box[m]);
1418 /* Added argument r2cut, changed contact and implemented
1419 * use of second cut-off.
1420 * - Erik Marklund, June 29, 2006
1422 static int is_hbond(t_hbdata *hb,int grpd,int grpa,int d,int a,
1423 real rcut, real r2cut, real ccut,
1424 rvec x[], gmx_bool bBox, matrix box,rvec hbox,
1425 real *d_ha, real *ang,gmx_bool bDA,int *hhh,
1426 gmx_bool bContact, gmx_bool bMerge, PSTYPE *p)
1429 rvec r_da,r_ha,r_dh, r={0, 0, 0};
1431 real rc2,r2c2,rda2,rha2,ca;
1432 gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
1433 gmx_bool daSwap = FALSE;
1438 if (((id = donor_index(&hb->d,grpd,d)) == NOTSET) ||
1439 ((ja = acceptor_index(&hb->a,grpa,a)) == NOTSET))
1445 rvec_sub(x[d],x[a],r_da);
1446 /* Insert projection code here */
1448 /* if (d>a && ((isInterchangable(hb, d, a, grpd, grpa) && bMerge) || bContact)) */
1449 /* /\* Then this hbond will be found again, or it has already been found. *\/ */
1453 if (d>a && bMerge && (bContact || isInterchangable(hb, d, a, grpd, grpa))) { /* acceptor is also a donor and vice versa? */
1455 daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
1458 copy_rvec(r_da, r); /* Save this for later */
1459 pbc_correct_gem(r_da,box,hbox);
1461 pbc_correct_gem(r_da,box,hbox);
1464 rda2 = iprod(r_da,r_da);
1471 calcBoxDistance(hb->per->P, r, ri);
1472 *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */
1476 else if (rda2 < r2c2)
1483 if (bDA && (rda2 > rc2))
1486 for(h=0; (h < hb->d.nhydro[id]); h++) {
1487 hh = hb->d.hydro[id][h];
1490 rvec_sub(x[hh],x[a],r_ha);
1492 pbc_correct_gem(r_ha,box,hbox);
1494 rha2 = iprod(r_ha,r_ha);
1498 calcBoxDistance(hb->per->P, r, ri);
1499 *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */
1502 if (bDA || (!bDA && (rha2 <= rc2))) {
1503 rvec_sub(x[d],x[hh],r_dh);
1505 pbc_correct_gem(r_dh,box,hbox);
1510 ca = cos_angle(r_dh,r_da);
1511 /* if angle is smaller, cos is larger */
1514 *d_ha = sqrt(bDA?rda2:rha2);
1520 if (bDA || (!bDA && HAinrange))
1526 /* Fixed previously undiscovered bug in the merge
1527 code, where the last frame of each hbond disappears.
1528 - Erik Marklund, June 1, 2006 */
1529 /* Added the following arguments:
1530 * ptmp[] - temporary periodicity hisory
1531 * a1 - identity of first acceptor/donor
1532 * a2 - identity of second acceptor/donor
1533 * - Erik Marklund, FEB 20 2010 */
1535 /* Merging is now done on the fly, so do_merge is most likely obsolete now.
1536 * Will do some more testing before removing the function entirely.
1537 * - Erik Marklund, MAY 10 2010 */
1538 static void do_merge(t_hbdata *hb,int ntmp,
1539 unsigned int htmp[],unsigned int gtmp[],PSTYPE ptmp[],
1540 t_hbond *hb0,t_hbond *hb1, int a1, int a2)
1542 /* Here we need to make sure we're treating periodicity in
1543 * the right way for the geminate recombination kinetics. */
1545 int m,mm,n00,n01,nn0,nnframes;
1549 /* Decide where to start from when merging */
1553 nnframes = max(n00 + hb0->nframes,n01 + hb1->nframes) - nn0;
1554 /* Initiate tmp arrays */
1555 for(m=0; (m<ntmp); m++) {
1560 /* Fill tmp arrays with values due to first HB */
1561 /* Once again '<' had to be replaced with '<='
1562 to catch the last frame in which the hbond
1564 - Erik Marklund, June 1, 2006 */
1565 for(m=0; (m<=hb0->nframes); m++) {
1567 htmp[mm] = is_hb(hb0->h[0],m);
1569 pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0);
1570 if (pm > hb->per->nper)
1571 gmx_fatal(FARGS, "Illegal shift!");
1573 ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/
1576 /* If we're doing geminate recompbination we usually don't need the distances.
1577 * Let's save some memory and time. */
1578 if (TRUE || !hb->bGem || hb->per->gemtype == gemAD){
1579 for(m=0; (m<=hb0->nframes); m++) {
1581 gtmp[mm] = is_hb(hb0->g[0],m);
1585 for(m=0; (m<=hb1->nframes); m++) {
1587 htmp[mm] = htmp[mm] || is_hb(hb1->h[0],m);
1588 gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0],m);
1589 if (hb->bGem /* && ptmp[mm] != 0 */) {
1591 /* If this hbond has been seen before with donor and acceptor swapped,
1592 * then we need to find the mirrored (*-1) periodicity vector to truely
1593 * merge the hbond history. */
1594 pm = findMirror(getPshift(hb->per->pHist[a2][a1],m+hb1->n0), hb->per->p2i, hb->per->nper);
1595 /* Store index of mirror */
1596 if (pm > hb->per->nper)
1597 gmx_fatal(FARGS, "Illegal shift!");
1601 /* Reallocate target array */
1602 if (nnframes > hb0->maxframes) {
1603 srenew(hb0->h[0],4+nnframes/hb->wordlen);
1604 srenew(hb0->g[0],4+nnframes/hb->wordlen);
1606 if (NULL != hb->per->pHist)
1608 clearPshift(&(hb->per->pHist[a1][a2]));
1611 /* Copy temp array to target array */
1612 for(m=0; (m<=nnframes); m++) {
1613 _set_hb(hb0->h[0],m,htmp[m]);
1614 _set_hb(hb0->g[0],m,gtmp[m]);
1616 addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0);
1619 /* Set scalar variables */
1621 hb0->maxframes = nnframes;
1624 /* Added argument bContact for nicer output.
1625 * Erik Marklund, June 29, 2006
1627 static void merge_hb(t_hbdata *hb,gmx_bool bTwo, gmx_bool bContact){
1628 int i,inrnew,indnew,j,ii,jj,m,id,ia,grp,ogrp,ntmp;
1629 unsigned int *htmp,*gtmp;
1634 indnew = hb->nrdist;
1636 /* Check whether donors are also acceptors */
1637 printf("Merging hbonds with Acceptor and Donor swapped\n");
1639 ntmp = 2*hb->max_frames;
1643 for(i=0; (i<hb->d.nrd); i++) {
1644 fprintf(stderr,"\r%d/%d",i+1,hb->d.nrd);
1646 ii = hb->a.aptr[id];
1647 for(j=0; (j<hb->a.nra); j++) {
1649 jj = hb->d.dptr[ia];
1650 if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
1651 (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j])))) {
1652 hb0 = hb->hbmap[i][j];
1653 hb1 = hb->hbmap[jj][ii];
1654 if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0])) {
1655 do_merge(hb,ntmp,htmp,gtmp,ptmp,hb0,hb1,i,j);
1656 if (ISHB(hb1->history[0]))
1658 else if (ISDIST(hb1->history[0]))
1662 gmx_incons("No contact history");
1664 gmx_incons("Neither hydrogen bond nor distance");
1668 clearPshift(&(hb->per->pHist[jj][ii]));
1672 hb1->history[0] = hbNo;
1677 fprintf(stderr,"\n");
1678 printf("- Reduced number of hbonds from %d to %d\n",hb->nrhb,inrnew);
1679 printf("- Reduced number of distances from %d to %d\n",hb->nrdist,indnew);
1681 hb->nrdist = indnew;
1687 static void do_nhb_dist(FILE *fp,t_hbdata *hb,real t)
1689 int i,j,k,n_bound[MAXHH],nbtot;
1693 /* Set array to 0 */
1694 for(k=0; (k<MAXHH); k++)
1696 /* Loop over possible donors */
1697 for(i=0; (i<hb->d.nrd); i++) {
1698 for(j=0; (j<hb->d.nhydro[i]); j++)
1699 n_bound[hb->d.nhbonds[i][j]]++;
1701 fprintf(fp,"%12.5e",t);
1703 for(k=0; (k<MAXHH); k++) {
1704 fprintf(fp," %8d",n_bound[k]);
1705 nbtot += n_bound[k]*k;
1707 fprintf(fp," %8d\n",nbtot);
1710 /* Added argument bContact in do_hblife(...). Also
1711 * added support for -contact in function body.
1712 * - Erik Marklund, May 31, 2006 */
1713 /* Changed the contact code slightly.
1714 * - Erik Marklund, June 29, 2006
1716 static void do_hblife(const char *fn,t_hbdata *hb,gmx_bool bMerge,gmx_bool bContact,
1717 const output_env_t oenv)
1720 const char *leg[] = { "p(t)", "t p(t)" };
1722 int i,j,j0,k,m,nh,ihb,ohb,nhydro,ndump=0;
1723 int nframes = hb->nframes;
1726 double sum,integral;
1729 snew(h,hb->maxhydro);
1730 snew(histo,nframes+1);
1731 /* Total number of hbonds analyzed here */
1732 for(i=0; (i<hb->d.nrd); i++) {
1733 for(k=0; (k<hb->a.nra); k++) {
1734 hbh = hb->hbmap[i][k];
1746 for(m=0; (m<hb->maxhydro); m++)
1748 h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m];
1751 for(nh=0; (nh<nhydro); nh++) {
1755 /* Changed '<' into '<=' below, just like I
1756 did in the hbm-output-loop in the main code.
1757 - Erik Marklund, May 31, 2006
1759 for(j=0; (j<=hbh->nframes); j++) {
1760 ihb = is_hb(h[nh],j);
1761 if (debug && (ndump < 10))
1762 fprintf(debug,"%5d %5d\n",j,ihb);
1778 fprintf(stderr,"\n");
1780 fp = xvgropen(fn,"Uninterrupted contact lifetime",output_env_get_xvgr_tlabel(oenv),"()",oenv);
1782 fp = xvgropen(fn,"Uninterrupted hydrogen bond lifetime",output_env_get_xvgr_tlabel(oenv),"()",
1785 xvgr_legend(fp,asize(leg),leg,oenv);
1787 while ((j0 > 0) && (histo[j0] == 0))
1790 for(i=0; (i<=j0); i++)
1792 dt = hb->time[1]-hb->time[0];
1795 for(i=1; (i<=j0); i++) {
1796 t = hb->time[i] - hb->time[0] - 0.5*dt;
1797 x1 = t*histo[i]/sum;
1798 fprintf(fp,"%8.3f %10.3e %10.3e\n",t,histo[i]/sum,x1);
1803 printf("%s lifetime = %.2f ps\n", bContact?"Contact":"HB", integral);
1804 printf("Note that the lifetime obtained in this manner is close to useless\n");
1805 printf("Use the -ac option instead and check the Forward lifetime\n");
1806 please_cite(stdout,"Spoel2006b");
1811 /* Changed argument bMerge into oneHB to handle contacts properly.
1812 * - Erik Marklund, June 29, 2006
1814 static void dump_ac(t_hbdata *hb,gmx_bool oneHB,int nDump)
1817 int i,j,k,m,nd,ihb,idist;
1818 int nframes = hb->nframes;
1824 fp = ffopen("debug-ac.xvg","w");
1825 for(j=0; (j<nframes); j++) {
1826 fprintf(fp,"%10.3f",hb->time[j]);
1827 for(i=nd=0; (i<hb->d.nrd) && (nd < nDump); i++) {
1828 for(k=0; (k<hb->a.nra) && (nd < nDump); k++) {
1831 hbh = hb->hbmap[i][k];
1834 ihb = is_hb(hbh->h[0],j);
1835 idist = is_hb(hbh->g[0],j);
1840 for(m=0; (m<hb->maxhydro) && !ihb ; m++) {
1841 ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m],j));
1842 idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m],j));
1844 /* This is not correct! */
1845 /* What isn't correct? -Erik M */
1849 fprintf(fp," %1d-%1d",ihb,idist);
1859 static real calc_dg(real tau,real temp)
1867 return kbt*log(kbt*tau/PLANCK);
1871 int n0,n1,nparams,ndelta;
1873 real *t,*ct,*nt,*kt,*sigma_ct,*sigma_nt,*sigma_kt;
1877 #include <gsl/gsl_multimin.h>
1878 #include <gsl/gsl_sf.h>
1879 #include <gsl/gsl_version.h>
1881 static double my_f(const gsl_vector *v,void *params)
1883 t_luzar *tl = (t_luzar *)params;
1885 double tol=1e-16,chi2=0;
1889 for(i=0; (i<tl->nparams); i++) {
1890 tl->kkk[i] = gsl_vector_get(v, i);
1895 for(i=tl->n0; (i<tl->n1); i+=tl->ndelta) {
1896 di=sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
1899 chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
1902 fprintf(stderr,"WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
1903 "di = %g k = %g kp = %g\n",tl->sigma_ct[i],
1904 tl->sigma_nt[i],tl->sigma_kt[i],di,k,kp);
1907 chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
1912 static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
1920 const gsl_multimin_fminimizer_type *T;
1921 gsl_multimin_fminimizer *s;
1924 gsl_multimin_function my_func;
1927 my_func.n = tl->nparams;
1928 my_func.params = (void *) tl;
1930 /* Starting point */
1931 x = gsl_vector_alloc (my_func.n);
1932 for(i=0; (i<my_func.n); i++)
1933 gsl_vector_set (x, i, tl->kkk[i]);
1935 /* Step size, different for each of the parameters */
1936 dx = gsl_vector_alloc (my_func.n);
1937 for(i=0; (i<my_func.n); i++)
1938 gsl_vector_set (dx, i, 0.01*tl->kkk[i]);
1940 T = gsl_multimin_fminimizer_nmsimplex;
1941 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
1943 gsl_multimin_fminimizer_set (s, &my_func, x, dx);
1944 gsl_vector_free (x);
1945 gsl_vector_free (dx);
1948 fprintf(fp,"%5s %12s %12s %12s %12s\n","Iter","k","kp","NM Size","Chi2");
1952 status = gsl_multimin_fminimizer_iterate (s);
1955 gmx_fatal(FARGS,"Something went wrong in the iteration in minimizer %s",
1956 gsl_multimin_fminimizer_name(s));
1958 d2 = gsl_multimin_fminimizer_minimum(s);
1959 size = gsl_multimin_fminimizer_size(s);
1960 status = gsl_multimin_test_size(size,tol);
1962 if (status == GSL_SUCCESS)
1964 fprintf(fp,"Minimum found using %s at:\n",
1965 gsl_multimin_fminimizer_name(s));
1968 fprintf(fp,"%5d", iter);
1969 for(i=0; (i<my_func.n); i++)
1970 fprintf(fp," %12.4e",gsl_vector_get (s->x,i));
1971 fprintf (fp," %12.4e %12.4e\n",size,d2);
1974 while ((status == GSL_CONTINUE) && (iter < maxiter));
1976 gsl_multimin_fminimizer_free (s);
1981 static real quality_of_fit(real chi2,int N)
1983 return gsl_sf_gamma_inc_Q((N-2)/2.0,chi2/2.0);
1987 static real optimize_luzar_parameters(FILE *fp,t_luzar *tl,int maxiter,
1990 fprintf(stderr,"This program needs the GNU scientific library to work.\n");
1994 static real quality_of_fit(real chi2,int N)
1996 fprintf(stderr,"This program needs the GNU scientific library to work.\n");
2003 static real compute_weighted_rates(int n,real t[],real ct[],real nt[],
2004 real kt[],real sigma_ct[],real sigma_nt[],
2005 real sigma_kt[],real *k,real *kp,
2006 real *sigma_k,real *sigma_kp,
2012 real kkk=0,kkp=0,kk2=0,kp2=0,chi2;
2017 for(i=0; (i<n); i++) {
2018 if (t[i] >= fit_start)
2029 tl.sigma_ct = sigma_ct;
2030 tl.sigma_nt = sigma_nt;
2031 tl.sigma_kt = sigma_kt;
2035 chi2 = optimize_luzar_parameters(debug,&tl,1000,1e-3);
2037 *kp = tl.kkk[1] = *kp;
2039 for(j=0; (j<NK); j++) {
2040 (void) optimize_luzar_parameters(debug,&tl,1000,1e-3);
2043 kk2 += sqr(tl.kkk[0]);
2044 kp2 += sqr(tl.kkk[1]);
2047 *sigma_k = sqrt(kk2/NK - sqr(kkk/NK));
2048 *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
2053 static void smooth_tail(int n,real t[],real c[],real sigma_c[],real start,
2054 const output_env_t oenv)
2057 real e_1,fitparm[4];
2061 for(i=0; (i<n); i++)
2069 do_lmfit(n,c,sigma_c,0,t,start,t[n-1],oenv,bDebugMode(),effnEXP2,fitparm,0);
2072 void analyse_corr(int n,real t[],real ct[],real nt[],real kt[],
2073 real sigma_ct[],real sigma_nt[],real sigma_kt[],
2074 real fit_start,real temp,real smooth_tail_start,
2075 const output_env_t oenv)
2078 real k=1,kp=1,kow=1;
2079 real Q=0,chi22,chi2,dg,dgp,tau_hb,dtau,tau_rlx,e_1,dt,sigma_k,sigma_kp,ddg;
2080 double tmp,sn2=0,sc2=0,sk2=0,scn=0,sck=0,snk=0;
2081 gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
2083 if (smooth_tail_start >= 0) {
2084 smooth_tail(n,t,ct,sigma_ct,smooth_tail_start,oenv);
2085 smooth_tail(n,t,nt,sigma_nt,smooth_tail_start,oenv);
2086 smooth_tail(n,t,kt,sigma_kt,smooth_tail_start,oenv);
2088 for(i0=0; (i0<n-2) && ((t[i0]-t[0]) < fit_start); i0++)
2091 for(i=i0; (i<n); i++) {
2099 printf("Hydrogen bond thermodynamics at T = %g K\n",temp);
2100 tmp = (sn2*sc2-sqr(scn));
2101 if ((tmp > 0) && (sn2 > 0)) {
2102 k = (sn2*sck-scn*snk)/tmp;
2103 kp = (k*scn-snk)/sn2;
2106 for(i=i0; (i<n); i++) {
2107 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2109 chi22 = compute_weighted_rates(n,t,ct,nt,kt,sigma_ct,sigma_nt,
2111 &sigma_k,&sigma_kp,fit_start);
2112 Q = quality_of_fit(chi2,2);
2113 ddg = BOLTZ*temp*sigma_k/k;
2114 printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
2116 printf("The Rate and Delta G are followed by an error estimate\n");
2117 printf("----------------------------------------------------------\n"
2118 "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
2119 printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2120 k,sigma_k,1/k,calc_dg(1/k,temp),ddg);
2121 ddg = BOLTZ*temp*sigma_kp/kp;
2122 printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
2123 kp,sigma_kp,1/kp,calc_dg(1/kp,temp),ddg);
2127 for(i=i0; (i<n); i++) {
2128 chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
2130 printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
2132 printf("--------------------------------------------------\n"
2133 "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
2134 printf("Forward %10.3f %8.3f %10.3f %10g\n",
2135 k,1/k,calc_dg(1/k,temp),chi2);
2136 printf("Backward %10.3f %8.3f %10.3f\n",
2137 kp,1/kp,calc_dg(1/kp,temp));
2142 printf("One-way %10.3f %s%8.3f %10.3f\n",
2143 kow,bError ? " " : "",1/kow,calc_dg(1/kow,temp));
2146 printf(" - Numerical problems computing HB thermodynamics:\n"
2147 "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n",
2148 sc2,sn2,sk2,sck,snk,scn);
2149 /* Determine integral of the correlation function */
2150 tau_hb = evaluate_integral(n,t,ct,NULL,(t[n-1]-t[0])/2,&dtau);
2151 printf("Integral %10.3f %s%8.3f %10.3f\n",1/tau_hb,
2152 bError ? " " : "",tau_hb,calc_dg(tau_hb,temp));
2154 for(i=0; (i<n-2); i++) {
2155 if ((ct[i] > e_1) && (ct[i+1] <= e_1)) {
2160 /* Determine tau_relax from linear interpolation */
2161 tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
2162 printf("Relaxation %10.3f %8.3f %s%10.3f\n",1/tau_rlx,
2163 tau_rlx,bError ? " " : "",
2164 calc_dg(tau_rlx,temp));
2168 printf("Correlation functions too short to compute thermodynamics\n");
2171 void compute_derivative(int nn,real x[],real y[],real dydx[])
2175 /* Compute k(t) = dc(t)/dt */
2176 for(j=1; (j<nn-1); j++)
2177 dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
2178 /* Extrapolate endpoints */
2179 dydx[0] = 2*dydx[1] - dydx[2];
2180 dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
2183 static void parallel_print(int *data, int nThreads)
2185 /* This prints the donors on which each tread is currently working. */
2188 fprintf(stderr, "\r");
2189 for (i=0; i<nThreads; i++)
2190 fprintf(stderr, "%-7i",data[i]);
2193 static void normalizeACF(real *ct, real *gt, int len)
2195 real ct_fac, gt_fac;
2198 /* Xu and Berne use the same normalization constant */
2201 gt_fac = (gt!=NULL && gt[0]!=0) ? 1.0/gt[0] : 0;
2202 printf("Normalization for c(t) = %g for gh(t) = %g\n",ct_fac,gt_fac);
2203 for (i=0; i<len; i++)
2211 /* Added argument bContact in do_hbac(...). Also
2212 * added support for -contact in the actual code.
2213 * - Erik Marklund, May 31, 2006 */
2214 /* Changed contact code and added argument R2
2215 * - Erik Marklund, June 29, 2006
2217 static void do_hbac(const char *fn,t_hbdata *hb,
2218 int nDump,gmx_bool bMerge,gmx_bool bContact, real fit_start,
2219 real temp,gmx_bool R2,real smooth_tail_start, const output_env_t oenv,
2220 t_gemParams *params, const char *gemType, int nThreads,
2221 const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit)
2224 int i,j,k,m,n,o,nd,ihb,idist,n2,nn,iter,nSets;
2225 const char *legNN[] = { "Ac(t)",
2227 static char **legGem;
2229 const char *legLuzar[] = { "Ac\\sfin sys\\v{}\\z{}(t)",
2231 "Cc\\scontact,hb\\v{}\\z{}(t)",
2232 "-dAc\\sfs\\v{}\\z{}/dt" };
2233 gmx_bool bNorm=FALSE;
2236 real *rhbex=NULL,*ht,*gt,*ght,*dght,*kt;
2237 real *ct,*p_ct,tail,tail2,dtail,ct_fac,ght_fac,*cct;
2238 const real tol = 1e-3;
2239 int nframes = hb->nframes,nf;
2240 unsigned int **h,**g;
2241 int nh,nhbonds,nhydro,ngh;
2243 PSTYPE p, *pfound = NULL, np;
2245 int *ptimes=NULL, *poff=NULL, anhb, n0, mMax=INT_MIN;
2246 real **rHbExGem = NULL;
2250 double *ctdouble, *timedouble, *fittedct;
2251 double fittolerance=0.1;
2253 enum {AC_NONE, AC_NN, AC_GEM, AC_LUZAR};
2257 int *dondata=NULL, thisThread;
2261 printf("Doing autocorrelation ");
2263 /* Decide what kind of ACF calculations to do. */
2264 if (NN > NN_NONE && NN < NN_NR) {
2265 #ifdef HAVE_NN_LOOPS
2267 printf("using the energy estimate.\n");
2270 printf("Can't do the NN-loop. Yet.\n");
2272 } else if (hb->bGem) {
2274 printf("according to the reversible geminate recombination model by Omer Markowitch.\n");
2276 nSets = 1 + (bBallistic ? 1:0) + (bGemFit ? 1:0);
2277 snew(legGem, nSets);
2278 for (i=0;i<nSets;i++)
2279 snew(legGem[i], 128);
2280 sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType);
2282 sprintf(legGem[1], "Ac'(t)");
2284 sprintf(legGem[(bBallistic ? 3:2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
2288 printf("according to the theory of Luzar and Chandler.\n");
2292 /* build hbexist matrix in reals for autocorr */
2293 /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */
2295 while (n2 < nframes)
2300 if (acType != AC_NN ||
2307 snew(h,hb->maxhydro);
2308 snew(g,hb->maxhydro);
2311 /* Dump hbonds for debugging */
2312 dump_ac(hb,bMerge||bContact,nDump);
2314 /* Total number of hbonds analyzed here */
2319 /* ------------------------------------------------
2320 * I got tired of waiting for the acf calculations
2321 * and parallelized it with openMP
2322 * set environment variable CFLAGS = "-fopenmp" when running
2323 * configure and define DOUSEOPENMP to make use of it.
2326 #ifdef HAVE_OPENMP /* ================================================= \
2327 * Set up the OpenMP stuff, |
2328 * like the number of threads and such |
2330 if (acType != AC_LUZAR)
2332 /* #if (_OPENMP >= 200805) /\* =====================\ *\/ */
2333 /* nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit()); */
2335 nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_num_procs());
2336 /* #endif /\* _OPENMP >= 200805 ====================/ *\/ */
2338 gmx_omp_set_num_threads(nThreads);
2339 snew(dondata, nThreads);
2340 for (i=0; i<nThreads; i++)
2342 printf("ACF calculations parallelized with OpenMP using %i threads.\n"
2343 "Expect close to linear scaling over this donor-loop.\n", nThreads);
2345 fprintf(stderr, "Donors: [thread no]\n");
2348 for (i=0;i<nThreads;i++) {
2349 snprintf(tmpstr, 7, "[%i]", i);
2350 fprintf(stderr, "%-7s", tmpstr);
2353 fprintf(stderr, "\n"); /* | */
2355 #endif /* HAVE_OPENMP ===================================================/ */
2358 /* Build the ACF according to acType */
2363 #ifdef HAVE_NN_LOOPS
2364 /* Here we're using the estimated energy for the hydrogen bonds. */
2366 #ifdef HAVE_OPENMP /* ==================================\ */
2367 #pragma omp parallel \
2368 private(i, j, k, nh, E, rhbex, thisThread), \
2372 thisThread = gmx_omp_get_thread_num();
2374 #endif /* ==============================================/ */
2377 memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
2379 #ifdef HAVE_OPENMP /* ################################################## \
2384 #pragma omp for schedule (dynamic)
2386 for (i=0; i<hb->d.nrd; i++) /* loop over donors */
2388 #ifdef HAVE_OPENMP /* ====== Write some output ======\ */
2389 #pragma omp critical
2391 dondata[thisThread] = i;
2392 parallel_print(dondata, nThreads);
2395 fprintf(stderr, "\r %i", i);
2396 #endif /* ===========================================/ */
2398 for (j=0; j<hb->a.nra; j++) /* loop over acceptors */
2400 for (nh=0; nh<hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */
2402 E = hb->hbE.E[i][j][nh];
2405 for (k=0; k<nframes; k++)
2407 if (E[k] != NONSENSE_E)
2408 rhbex[k] = (real)E[k];
2411 low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rhbex),hb->time[1]-hb->time[0],
2412 eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2414 #pragma omp critical
2417 for(k=0; (k<nn); k++)
2428 } /* End of parallel block # */
2429 /* ##############################################################/ */
2432 normalizeACF(ct, NULL, nn);
2434 snew(timedouble, nn);
2435 for (j=0; j<nn; j++)
2437 timedouble[j]=(double)(hb->time[j]);
2438 ctdouble[j]=(double)(ct[j]);
2441 /* Remove ballistic term */
2442 /* Ballistic component removal and fitting to the reversible geminate recombination model
2443 * will be taken out for the time being. First of all, one can remove the ballistic
2444 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2445 * problems with the robustness of the fitting to the model. More work is needed.
2446 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2447 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2448 * a BSD-licence that can do the job.
2450 * - Erik Marklund, June 18 2010.
2452 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2453 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2455 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2456 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2458 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)");
2459 xvgr_legend(fp,asize(legNN),legNN);
2461 for(j=0; (j<nn); j++)
2462 fprintf(fp,"%10g %10g %10g\n",
2463 hb->time[j]-hb->time[0],
2470 #endif /* HAVE_NN_LOOPS */
2471 break; /* case AC_NN */
2475 memset(ct,0,2*n2*sizeof(real));
2477 fprintf(stderr, "Donor:\n");
2480 #define __ACDATA p_ct
2483 #ifdef HAVE_OPENMP /* =========================================\
2485 #pragma omp parallel default(none) \
2486 private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
2487 pfound, poff, rHbExGem, p, ihb, mMax, \
2489 shared(hb, dondata, ct, nn, nThreads, n2, stderr, bNorm, \
2490 nframes, bMerge, bContact)
2491 { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
2494 thisThread = gmx_omp_get_thread_num();
2495 snew(h,hb->maxhydro);
2496 snew(g,hb->maxhydro);
2503 memset(p_ct,0,2*n2*sizeof(real));
2505 /* I'm using a chunk size of 1, since I expect \
2506 * the overhead to be really small compared \
2507 * to the actual calculations \ */
2508 #pragma omp for schedule(dynamic,1) nowait /* \ */
2509 #endif /* HAVE_OPENMP =========================================/ */
2511 for (i=0; i<hb->d.nrd; i++) {
2513 #pragma omp critical
2515 dondata[thisThread] = i;
2516 parallel_print(dondata, nThreads);
2519 fprintf(stderr, "\r %i", i);
2522 for (k=0; k<hb->a.nra; k++) {
2523 for (nh=0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) {
2524 hbh = hb->hbmap[i][k];
2526 /* Note that if hb->per->gemtype==gemDD, then distances will be stored in
2527 * hb->hbmap[d][a].h array anyway, because the contact flag will be set.
2528 * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */
2529 pHist = &(hb->per->pHist[i][k]);
2530 if (ISHB(hbh->history[nh]) && pHist->len != 0) {
2532 /* No need for a critical section */
2533 /* #ifdef HAVE_OPENMP */
2534 /* #pragma omp critical */
2538 g[nh] = hb->per->gemtype==gemAD ? hbh->g[nh] : NULL;
2542 /* count the number of periodic shifts encountered and store
2543 * them in separate arrays. */
2545 for (j=0; j<pHist->len; j++)
2548 for (m=0; m<=np; m++) {
2549 if (m == np) { /* p not recognized in list. Add it and set up new array. */
2551 if (np>hb->per->nper)
2552 gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here.");
2553 if (m>=mMax) { /* Extend the arrays.
2554 * Doing it like this, using mMax to keep track of the sizes,
2555 * eleviates the need for freeing and re-allocating the arrays
2556 * when taking on the next donor-acceptor pair */
2558 srenew(pfound,np); /* The list of found periodic shifts. */
2559 srenew(rHbExGem,np); /* The hb existence functions (-aver_hb). */
2560 snew(rHbExGem[m],2*n2);
2564 /* This shouldn't have to be critical, right? */
2565 /* #ifdef HAVE_OPENMP */
2566 /* #pragma omp critical */
2569 if (rHbExGem != NULL && rHbExGem[m] != NULL) {
2570 /* This must be done, as this array was most likey
2571 * used to store stuff in some previous iteration. */
2572 memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2));
2575 fprintf(stderr, "rHbExGem not initialized! m = %i\n", m);
2584 } /* m: Loop over found shifts */
2585 } /* j: Loop over shifts */
2587 /* Now unpack and disentangle the existence funtions. */
2588 for (j=0; j<nf; j++) {
2594 * pfound: list of periodic shifts found for this pair.
2595 * poff: list of frame offsets; that is, the first
2596 * frame a hbond has a particular periodic shift. */
2597 p = getPshift(*pHist, j+n0);
2600 for (m=0; m<np; m++)
2605 gmx_fatal(FARGS,"Shift not found, but must be there.");
2608 ihb = is_hb(h[nh],j) || ((hb->per->gemtype!=gemAD || j==0) ? FALSE : is_hb(g[nh],j));
2612 poff[m] = j; /* Here's where the first hbond with shift p is,
2613 * relative to the start of h[0].*/
2615 gmx_fatal(FARGS, "j<poff[m]");
2616 rHbExGem[m][j-poff[m]] += 1;
2621 /* Now, build ac. */
2622 for (m=0; m<np; m++) {
2623 if (rHbExGem[m][0]>0 && n0+poff[m]<nn/* && m==0 */) {
2624 low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rHbExGem[m]),hb->time[1]-hb->time[0],
2625 eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2626 for(j=0; (j<nn); j++)
2627 __ACDATA[j] += rHbExGem[m][j];
2629 } /* Building of ac. */
2632 } /* hydrogen loop */
2633 } /* acceptor loop */
2636 for (m=0; m<=mMax; m++) {
2645 #ifdef HAVE_OPENMP /* =======================================\ */
2646 #pragma omp critical
2648 for (i=0; i<nn; i++)
2653 } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
2655 #endif /* HAVE_OPENMP =======================================/ */
2657 normalizeACF(ct, NULL, nn);
2659 fprintf(stderr, "\n\nACF successfully calculated.\n");
2661 /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */
2664 snew(timedouble, nn);
2668 timedouble[j]=(double)(hb->time[j]);
2669 ctdouble[j]=(double)(ct[j]);
2672 /* Remove ballistic term */
2673 /* Ballistic component removal and fitting to the reversible geminate recombination model
2674 * will be taken out for the time being. First of all, one can remove the ballistic
2675 * component with g_analyze afterwards. Secondly, and more importantly, there are still
2676 * problems with the robustness of the fitting to the model. More work is needed.
2677 * A third reason is that we're currently using gsl for this and wish to reduce dependence
2678 * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with
2679 * a BSD-licence that can do the job.
2681 * - Erik Marklund, June 18 2010.
2683 /* if (bBallistic) { */
2684 /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */
2685 /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */
2687 /* printf("\nNumber of data points is less than the number of parameters to fit\n." */
2688 /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */
2691 /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */
2695 fp = xvgropen(fn, "Contact Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
2697 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
2698 xvgr_legend(fp,asize(legGem),(const char**)legGem,oenv);
2700 for(j=0; (j<nn); j++)
2702 fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0],ct[j]);
2704 fprintf(fp," %10g", ctdouble[j]);
2706 fprintf(fp," %10g", fittedct[j]);
2716 break; /* case AC_GEM */
2729 for(i=0; (i<hb->d.nrd); i++) {
2730 for(k=0; (k<hb->a.nra); k++) {
2732 hbh = hb->hbmap[i][k];
2735 if (bMerge || bContact) {
2736 if (ISHB(hbh->history[0])) {
2743 for(m=0; (m<hb->maxhydro); m++)
2744 if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) {
2745 g[nhydro] = hbh->g[m];
2746 h[nhydro] = hbh->h[m];
2752 for(nh=0; (nh<nhydro); nh++) {
2753 int nrint = bContact ? hb->nrdist : hb->nrhb;
2754 if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
2755 fprintf(stderr,"\rACF %d/%d",nhbonds+1,nrint);
2757 for(j=0; (j<nframes); j++) {
2758 /* Changed '<' into '<=' below, just like I did in
2759 the hbm-output-loop in the gmx_hbond() block.
2760 - Erik Marklund, May 31, 2006 */
2762 ihb = is_hb(h[nh],j);
2763 idist = is_hb(g[nh],j);
2769 /* For contacts: if a second cut-off is provided, use it,
2770 * otherwise use g(t) = 1-h(t) */
2771 if (!R2 && bContact)
2774 gt[j] = idist*(1-ihb);
2780 /* The autocorrelation function is normalized after summation only */
2781 low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&rhbex,hb->time[1]-hb->time[0],
2782 eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
2784 /* Cross correlation analysis for thermodynamics */
2785 for(j=nframes; (j<n2); j++) {
2790 cross_corr(n2,ht,gt,dght);
2792 for(j=0; (j<nn); j++) {
2800 fprintf(stderr,"\n");
2803 normalizeACF(ct, gt, nn);
2805 /* Determine tail value for statistics */
2808 for(j=nn/2; (j<nn); j++) {
2810 tail2 += ct[j]*ct[j];
2812 tail /= (nn - nn/2);
2813 tail2 /= (nn - nn/2);
2814 dtail = sqrt(tail2-tail*tail);
2816 /* Check whether the ACF is long enough */
2818 printf("\nWARNING: Correlation function is probably not long enough\n"
2819 "because the standard deviation in the tail of C(t) > %g\n"
2820 "Tail value (average C(t) over second half of acf): %g +/- %g\n",
2823 for(j=0; (j<nn); j++) {
2825 ct[j] = (cct[j]-tail)/(1-tail);
2827 /* Compute negative derivative k(t) = -dc(t)/dt */
2828 compute_derivative(nn,hb->time,ct,kt);
2829 for(j=0; (j<nn); j++)
2834 fp = xvgropen(fn, "Contact Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
2836 fp = xvgropen(fn, "Hydrogen Bond Autocorrelation",output_env_get_xvgr_tlabel(oenv),"C(t)",oenv);
2837 xvgr_legend(fp,asize(legLuzar),legLuzar, oenv);
2840 for(j=0; (j<nn); j++)
2841 fprintf(fp,"%10g %10g %10g %10g %10g\n",
2842 hb->time[j]-hb->time[0],ct[j],cct[j],ght[j],kt[j]);
2845 analyse_corr(nn,hb->time,ct,ght,kt,NULL,NULL,NULL,
2846 fit_start,temp,smooth_tail_start,oenv);
2848 do_view(oenv,fn,NULL);
2860 break; /* case AC_LUZAR */
2863 gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType);
2864 } /* switch (acType) */
2867 static void init_hbframe(t_hbdata *hb,int nframes,real t)
2871 hb->time[nframes] = t;
2872 hb->nhb[nframes] = 0;
2873 hb->ndist[nframes] = 0;
2874 for (i=0; (i<max_hx); i++)
2875 hb->nhx[nframes][i]=0;
2876 /* Loop invalidated */
2877 if (hb->bHBmap && 0)
2878 for (i=0; (i<hb->d.nrd); i++)
2879 for (j=0; (j<hb->a.nra); j++)
2880 for (m=0; (m<hb->maxhydro); m++)
2881 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m])
2882 set_hb(hb,i,m,j,nframes,HB_NO);
2883 /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
2886 static void analyse_donor_props(const char *fn,t_hbdata *hb,int nframes,real t,
2887 const output_env_t oenv)
2889 static FILE *fp = NULL;
2890 const char *leg[] = { "Nbound", "Nfree" };
2891 int i,j,k,nbound,nb,nhtot;
2896 fp = xvgropen(fn,"Donor properties",output_env_get_xvgr_tlabel(oenv),"Number",oenv);
2897 xvgr_legend(fp,asize(leg),leg,oenv);
2901 for(i=0; (i<hb->d.nrd); i++) {
2902 for(k=0; (k<hb->d.nhydro[i]); k++) {
2905 for(j=0; (j<hb->a.nra) && (nb == 0); j++) {
2906 if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
2907 is_hb(hb->hbmap[i][j]->h[k],nframes))
2913 fprintf(fp,"%10.3e %6d %6d\n",t,nbound,nhtot-nbound);
2916 static void dump_hbmap(t_hbdata *hb,
2917 int nfile,t_filenm fnm[],gmx_bool bTwo,
2918 gmx_bool bContact, int isize[],int *index[],char *grpnames[],
2922 int ddd,hhh,aaa,i,j,k,m,grp;
2923 char ds[32],hs[32],as[32];
2926 fp = opt2FILE("-hbn",nfile,fnm,"w");
2927 if (opt2bSet("-g",nfile,fnm)) {
2928 fplog = ffopen(opt2fn("-g",nfile,fnm),"w");
2929 fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
2933 for (grp=gr0; grp<=(bTwo?gr1:gr0); grp++) {
2934 fprintf(fp,"[ %s ]",grpnames[grp]);
2935 for (i=0; i<isize[grp]; i++) {
2936 fprintf(fp,(i%15)?" ":"\n");
2937 fprintf(fp," %4u",index[grp][i]+1);
2941 Added -contact support below.
2942 - Erik Marklund, May 29, 2006
2945 fprintf(fp,"[ donors_hydrogens_%s ]\n",grpnames[grp]);
2946 for (i=0; (i<hb->d.nrd); i++) {
2947 if (hb->d.grp[i] == grp) {
2948 for(j=0; (j<hb->d.nhydro[i]); j++)
2949 fprintf(fp," %4u %4u",hb->d.don[i]+1,
2950 hb->d.hydro[i][j]+1);
2955 fprintf(fp,"[ acceptors_%s ]",grpnames[grp]);
2956 for (i=0; (i<hb->a.nra); i++) {
2957 if (hb->a.grp[i] == grp) {
2958 fprintf(fp,(i%15 && !first)?" ":"\n");
2959 fprintf(fp," %4u",hb->a.acc[i]+1);
2967 fprintf(fp,bContact ? "[ contacts_%s-%s ]\n" :
2968 "[ hbonds_%s-%s ]\n",grpnames[0],grpnames[1]);
2970 fprintf(fp,bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n",grpnames[0]);
2972 for(i=0; (i<hb->d.nrd); i++) {
2974 for(k=0; (k<hb->a.nra); k++) {
2976 for(m=0; (m<hb->d.nhydro[i]); m++) {
2977 if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m])) {
2978 sprintf(ds,"%s",mkatomname(atoms,ddd));
2979 sprintf(as,"%s",mkatomname(atoms,aaa));
2981 fprintf(fp," %6u %6u\n",ddd+1,aaa+1);
2983 fprintf(fplog,"%12s %12s\n",ds,as);
2985 hhh = hb->d.hydro[i][m];
2986 sprintf(hs,"%s",mkatomname(atoms,hhh));
2987 fprintf(fp," %6u %6u %6u\n",ddd+1,hhh+1,aaa+1);
2989 fprintf(fplog,"%12s %12s %12s\n",ds,hs,as);
3001 /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
3002 * It mimics add_frames() and init_frame() to some extent. */
3003 static void sync_hbdata(t_hbdata *hb, t_hbdata *p_hb,
3004 int nframes, real t)
3007 if (nframes >= p_hb->max_frames)
3009 p_hb->max_frames += 4096;
3010 srenew(p_hb->nhb, p_hb->max_frames);
3011 srenew(p_hb->ndist, p_hb->max_frames);
3012 srenew(p_hb->n_bound, p_hb->max_frames);
3013 srenew(p_hb->nhx, p_hb->max_frames);
3015 srenew(p_hb->danr, p_hb->max_frames);
3016 memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3017 memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
3018 p_hb->nhb[nframes] = 0;
3019 p_hb->ndist[nframes] = 0;
3022 p_hb->nframes = nframes;
3025 /* p_hb->nhx[nframes][i] */
3027 memset(&(p_hb->nhx[nframes]), 0 ,sizeof(int)*max_hx); /* zero the helix count for this frame */
3029 /* hb->per will remain constant througout the frame loop,
3030 * even though the data its members point to will change,
3031 * hence no need for re-syncing. */
3035 int gmx_hbond(int argc,char *argv[])
3037 const char *desc[] = {
3038 "[TT]g_hbond[tt] computes and analyzes hydrogen bonds. Hydrogen bonds are",
3039 "determined based on cutoffs for the angle Acceptor - Donor - Hydrogen",
3040 "(zero is extended) and the distance Hydrogen - Acceptor.",
3041 "OH and NH groups are regarded as donors, O is an acceptor always,",
3042 "N is an acceptor by default, but this can be switched using",
3043 "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected",
3044 "to the first preceding non-hydrogen atom.[PAR]",
3046 "You need to specify two groups for analysis, which must be either",
3047 "identical or non-overlapping. All hydrogen bonds between the two",
3048 "groups are analyzed.[PAR]",
3050 "If you set [TT]-shell[tt], you will be asked for an additional index group",
3051 "which should contain exactly one atom. In this case, only hydrogen",
3052 "bonds between atoms within the shell distance from the one atom are",
3055 "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
3056 "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
3057 "If contact kinetics are analyzed by using the -contact option, then",
3058 "n(t) can be defined as either all pairs that are not within contact distance r at time t",
3059 "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
3060 "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
3061 "See mentioned literature for more details and definitions."
3064 /* "It is also possible to analyse specific hydrogen bonds with",
3065 "[TT]-sel[tt]. This index file must contain a group of atom triplets",
3066 "Donor Hydrogen Acceptor, in the following way:[PAR]",
3074 "Note that the triplets need not be on separate lines.",
3075 "Each atom triplet specifies a hydrogen bond to be analyzed,",
3076 "note also that no check is made for the types of atoms.[PAR]",
3078 "[BB]Output:[bb][BR]",
3079 "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
3080 "[TT]-ac[tt]: average over all autocorrelations of the existence",
3081 "functions (either 0 or 1) of all hydrogen bonds.[BR]",
3082 "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
3083 "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
3084 "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
3085 "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
3086 "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
3087 "with helices in proteins.[BR]",
3088 "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
3089 "for selected groups, all hydrogen bonded atoms from all groups and",
3090 "all solvent atoms involved in insertion.[BR]",
3091 "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
3092 "frames, this also contains information on solvent insertion",
3093 "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
3095 "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
3096 "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
3097 "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
3098 "compare results to Raman Spectroscopy.",
3100 "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
3101 "require an amount of memory proportional to the total numbers of donors",
3102 "times the total number of acceptors in the selected group(s)."
3105 static real acut=30, abin=1, rcut=0.35, r2cut=0, rbin=0.005, rshell=-1;
3106 static real maxnhb=0,fit_start=1,fit_end=60,temp=298.15,smooth_tail_start=-1, D=-1;
3107 static gmx_bool bNitAcc=TRUE,bDA=TRUE,bMerge=TRUE;
3108 static int nDump=0, nFitPoints=100;
3109 static int nThreads = 0, nBalExp=4;
3111 static gmx_bool bContact=FALSE, bBallistic=FALSE, bBallisticDt=FALSE, bGemFit=FALSE;
3112 static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
3113 static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
3117 { "-a", FALSE, etREAL, {&acut},
3118 "Cutoff angle (degrees, Acceptor - Donor - Hydrogen)" },
3119 { "-r", FALSE, etREAL, {&rcut},
3120 "Cutoff radius (nm, X - Acceptor, see next option)" },
3121 { "-da", FALSE, etBOOL, {&bDA},
3122 "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
3123 { "-r2", FALSE, etREAL, {&r2cut},
3124 "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
3125 { "-abin", FALSE, etREAL, {&abin},
3126 "Binwidth angle distribution (degrees)" },
3127 { "-rbin", FALSE, etREAL, {&rbin},
3128 "Binwidth distance distribution (nm)" },
3129 { "-nitacc",FALSE, etBOOL, {&bNitAcc},
3130 "Regard nitrogen atoms as acceptors" },
3131 { "-contact",FALSE,etBOOL, {&bContact},
3132 "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
3133 { "-shell", FALSE, etREAL, {&rshell},
3134 "when > 0, only calculate hydrogen bonds within # nm shell around "
3136 { "-fitstart", FALSE, etREAL, {&fit_start},
3137 "Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation. With [TT]-gemfit[tt] we suggest [TT]-fitstart 0[tt]" },
3138 { "-fitstart", FALSE, etREAL, {&fit_start},
3139 "Time (ps) to which to stop fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation (only with [TT]-gemfit[tt])" },
3140 { "-temp", FALSE, etREAL, {&temp},
3141 "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
3142 { "-smooth",FALSE, etREAL, {&smooth_tail_start},
3143 "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" },
3144 { "-dump", FALSE, etINT, {&nDump},
3145 "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
3146 { "-max_hb",FALSE, etREAL, {&maxnhb},
3147 "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
3148 { "-merge", FALSE, etBOOL, {&bMerge},
3149 "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
3150 { "-geminate", FALSE, etENUM, {gemType},
3151 "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
3152 { "-diff", FALSE, etREAL, {&D},
3153 "Dffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
3155 { "-nthreads", FALSE, etINT, {&nThreads},
3156 "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)"},
3159 const char *bugs[] = {
3160 "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
3163 { efTRX, "-f", NULL, ffREAD },
3164 { efTPX, NULL, NULL, ffREAD },
3165 { efNDX, NULL, NULL, ffOPTRD },
3166 /* { efNDX, "-sel", "select", ffOPTRD },*/
3167 { efXVG, "-num", "hbnum", ffWRITE },
3168 { efLOG, "-g", "hbond", ffOPTWR },
3169 { efXVG, "-ac", "hbac", ffOPTWR },
3170 { efXVG, "-dist","hbdist", ffOPTWR },
3171 { efXVG, "-ang", "hbang", ffOPTWR },
3172 { efXVG, "-hx", "hbhelix",ffOPTWR },
3173 { efNDX, "-hbn", "hbond", ffOPTWR },
3174 { efXPM, "-hbm", "hbmap", ffOPTWR },
3175 { efXVG, "-don", "donor", ffOPTWR },
3176 { efXVG, "-dan", "danum", ffOPTWR },
3177 { efXVG, "-life","hblife", ffOPTWR },
3178 { efXVG, "-nhbdist", "nhbdist",ffOPTWR }
3181 #define NFILE asize(fnm)
3183 char hbmap [HB_NR]={ ' ', 'o', '-', '*' };
3184 const char *hbdesc[HB_NR]={ "None", "Present", "Inserted", "Present & Inserted" };
3185 t_rgb hbrgb [HB_NR]={ {1,1,1},{1,0,0}, {0,0,1}, {1,0,1} };
3187 t_trxstatus *status;
3192 int npargs,natoms,nframes=0,shatom;
3198 real t,ccut,dist=0.0,ang=0.0;
3199 double max_nhb,aver_nhb,aver_dist;
3200 int h=0,i,j,k=0,l,start,end,id,ja,ogrp,nsel;
3202 int xj,yj,zj,aj,xjj,yjj,zjj;
3203 int xk,yk,zk,ak,xkk,ykk,zkk;
3204 gmx_bool bSelected,bHBmap,bStop,bTwo,was,bBox,bTric;
3206 int grp,nabin,nrbin,bin,resdist,ihb;
3209 FILE *fp,*fpins=NULL,*fpnhb=NULL;
3211 t_ncell *icell,*jcell,*kcell;
3213 unsigned char *datable;
3218 int ii, jj, hh, actual_nThreads;
3220 gmx_bool bGem, bNN, bParallel;
3221 t_gemParams *params=NULL;
3222 gmx_bool bEdge_yjj, bEdge_xjj;
3224 CopyRight(stdout,argv[0]);
3227 ppa = add_acf_pargs(&npargs,pa);
3229 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,NFILE,fnm,npargs,
3230 ppa,asize(desc),desc,asize(bugs),bugs,&oenv);
3232 /* NN-loop? If so, what estimator to use ?*/
3234 /* Outcommented for now DvdS 2010-07-13
3235 while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0)
3238 gmx_fatal(FARGS, "Invalid NN-loop type.");
3241 for (i=2; bNN==FALSE && i<NN_NR; i++)
3242 bNN = bNN || NN == i;
3244 if (NN > NN_NONE && bMerge)
3247 /* geminate recombination? If so, which flavor? */
3249 while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode])!=0)
3251 if (gemmode == gemNR)
3252 gmx_fatal(FARGS, "Invalid recombination type.");
3255 for (i=2; bGem==FALSE && i<gemNR; i++)
3256 bGem = bGem || gemmode == i;
3259 printf("Geminate recombination: %s\n" ,gemType[gemmode]);
3261 printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
3264 if (gemmode != gemDD) {
3265 printf("Turning off -contact option...\n");
3269 if (gemmode == gemDD) {
3270 printf("Turning on -contact option...\n");
3275 if (gemmode == gemAA) {
3276 printf("Turning off -merge option...\n");
3280 if (gemmode != gemAA) {
3281 printf("Turning on -merge option...\n");
3289 ccut = cos(acut*DEG2RAD);
3293 gmx_fatal(FARGS,"Can not analyze selected contacts.");
3295 gmx_fatal(FARGS,"Can not analyze contact between H and A: turn off -noda");
3300 /* Don't pollute stdout with information about external libraries.
3302 * printf("NO GSL! Can't find and take away ballistic term in ACF without GSL\n.");
3306 /* Initiate main data structure! */
3307 bHBmap = (opt2bSet("-ac",NFILE,fnm) ||
3308 opt2bSet("-life",NFILE,fnm) ||
3309 opt2bSet("-hbn",NFILE,fnm) ||
3310 opt2bSet("-hbm",NFILE,fnm) ||
3314 /* Same thing here. There is no reason whatsoever to write the specific version of
3315 * OpenMP used for compilation to stdout for normal usage.
3317 * printf("Compiled with OpenMP (%i)\n", _OPENMP);
3321 /* if (bContact && bGem) */
3322 /* gmx_fatal(FARGS, "Can't do reversible geminate recombination with -contact yet."); */
3324 if (opt2bSet("-nhbdist",NFILE,fnm)) {
3325 const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
3326 fpnhb = xvgropen(opt2fn("-nhbdist",NFILE,fnm),
3327 "Number of donor-H with N HBs",output_env_get_xvgr_tlabel(oenv),"N",oenv);
3328 xvgr_legend(fpnhb,asize(leg),leg,oenv);
3331 hb = mk_hbdata(bHBmap,opt2bSet("-dan",NFILE,fnm),bMerge || bContact, bGem, gemmode);
3334 read_tpx_top(ftp2fn(efTPX,NFILE,fnm),&ir,box,&natoms,NULL,NULL,NULL,&top);
3336 snew(grpnames,grNR);
3339 /* Make Donor-Acceptor table */
3340 snew(datable, top.atoms.nr);
3341 gen_datable(index[0],isize[0],datable,top.atoms.nr);
3344 /* analyze selected hydrogen bonds */
3345 printf("Select group with selected atoms:\n");
3346 get_index(&(top.atoms),opt2fn("-sel",NFILE,fnm),
3347 1,&nsel,index,grpnames);
3349 gmx_fatal(FARGS,"Number of atoms in group '%s' not a multiple of 3\n"
3350 "and therefore cannot contain triplets of "
3351 "Donor-Hydrogen-Acceptor",grpnames[0]);
3354 for(i=0; (i<nsel); i+=3) {
3355 int dd = index[0][i];
3356 int aa = index[0][i+2];
3357 /* int */ hh = index[0][i+1];
3358 add_dh (&hb->d,dd,hh,i,datable);
3359 add_acc(&hb->a,aa,i);
3360 /* Should this be here ? */
3361 snew(hb->d.dptr,top.atoms.nr);
3362 snew(hb->a.aptr,top.atoms.nr);
3363 add_hbond(hb,dd,aa,hh,gr0,gr0,0,bMerge,0,bContact,peri);
3365 printf("Analyzing %d selected hydrogen bonds from '%s'\n",
3366 isize[0],grpnames[0]);
3369 /* analyze all hydrogen bonds: get group(s) */
3370 printf("Specify 2 groups to analyze:\n");
3371 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
3372 2,isize,index,grpnames);
3374 /* check if we have two identical or two non-overlapping groups */
3375 bTwo = isize[0] != isize[1];
3376 for (i=0; (i<isize[0]) && !bTwo; i++)
3377 bTwo = index[0][i] != index[1][i];
3379 printf("Checking for overlap in atoms between %s and %s\n",
3380 grpnames[0],grpnames[1]);
3381 for (i=0; i<isize[1];i++)
3382 if (ISINGRP(datable[index[1][i]]))
3383 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3384 grpnames[0],grpnames[1]);
3386 printf("Checking for overlap in atoms between %s and %s\n",
3387 grpnames[0],grpnames[1]);
3388 for (i=0; i<isize[0]; i++)
3389 for (j=0; j<isize[1]; j++)
3390 if (index[0][i] == index[1][j])
3391 gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'",
3392 grpnames[0],grpnames[1]);
3396 printf("Calculating %s "
3397 "between %s (%d atoms) and %s (%d atoms)\n",
3398 bContact ? "contacts" : "hydrogen bonds",
3399 grpnames[0],isize[0],grpnames[1],isize[1]);
3401 fprintf(stderr,"Calculating %s in %s (%d atoms)\n",
3402 bContact?"contacts":"hydrogen bonds",grpnames[0],isize[0]);
3406 /* search donors and acceptors in groups */
3407 snew(datable, top.atoms.nr);
3408 for (i=0; (i<grNR); i++)
3409 if ( ((i==gr0) && !bSelected ) ||
3410 ((i==gr1) && bTwo )) {
3411 gen_datable(index[i],isize[i],datable,top.atoms.nr);
3413 search_acceptors(&top,isize[i],index[i],&hb->a,i,
3414 bNitAcc,TRUE,(bTwo && (i==gr0)) || !bTwo, datable);
3415 search_donors (&top,isize[i],index[i],&hb->d,i,
3416 TRUE,(bTwo && (i==gr1)) || !bTwo, datable);
3419 search_acceptors(&top,isize[i],index[i],&hb->a,i,bNitAcc,FALSE,TRUE, datable);
3420 search_donors (&top,isize[i],index[i],&hb->d,i,FALSE,TRUE, datable);
3423 clear_datable_grp(datable,top.atoms.nr);
3426 printf("Found %d donors and %d acceptors\n",hb->d.nrd,hb->a.nra);
3428 snew(donors[gr0D], dons[gr0D].nrd);*/
3431 printf("Making hbmap structure...");
3432 /* Generate hbond data structure */
3437 #ifdef HAVE_NN_LOOPS
3443 printf("Making per structure...");
3444 /* Generate hbond data structure */
3451 if (hb->d.nrd + hb->a.nra == 0) {
3452 printf("No Donors or Acceptors found\n");
3456 if (hb->d.nrd == 0) {
3457 printf("No Donors found\n");
3460 if (hb->a.nra == 0) {
3461 printf("No Acceptors found\n");
3466 gmx_fatal(FARGS,"Nothing to be done");
3473 /* get index group with atom for shell */
3475 printf("Select atom for shell (1 atom):\n");
3476 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
3477 1,&shisz,&shidx,&shgrpnm);
3479 printf("group contains %d atoms, should be 1 (one)\n",shisz);
3480 } while(shisz != 1);
3482 printf("Will calculate hydrogen bonds within a shell "
3483 "of %g nm around atom %i\n",rshell,shatom+1);
3486 /* Analyze trajectory */
3487 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
3488 if ( natoms > top.atoms.nr )
3489 gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
3490 top.atoms.nr,natoms);
3492 bBox = ir.ePBC!=epbcNONE;
3493 grid = init_grid(bBox, box, (rcut>r2cut)?rcut:r2cut, ngrid);
3496 snew(adist,nabin+1);
3497 snew(rdist,nrbin+1);
3500 gmx_fatal(FARGS, "Can't do geminate recombination without periodic box.");
3505 #define __ADIST adist
3506 #define __RDIST rdist
3508 #else /* HAVE_OPENMP ================================================== \
3509 * Set up the OpenMP stuff, |
3510 * like the number of threads and such |
3511 * Also start the parallel loop. |
3513 #define __ADIST p_adist[threadNr]
3514 #define __RDIST p_rdist[threadNr]
3515 #define __HBDATA p_hb[threadNr]
3517 bParallel = !bSelected;
3521 /* #if (_OPENMP > 200805) */
3522 /* actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit()); */
3524 actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_num_procs());
3526 gmx_omp_set_num_threads(actual_nThreads);
3527 printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
3532 actual_nThreads = 1;
3535 t_hbdata **p_hb; /* one per thread, then merge after the frame loop */
3536 int **p_adist, **p_rdist; /* a histogram for each thread. */
3537 snew(p_hb, actual_nThreads);
3538 snew(p_adist, actual_nThreads);
3539 snew(p_rdist, actual_nThreads);
3540 for (i=0; i<actual_nThreads; i++)
3543 snew(p_adist[i], nabin+1);
3544 snew(p_rdist[i], nrbin+1);
3546 p_hb[i]->max_frames = 0;
3547 p_hb[i]->nhb = NULL;
3548 p_hb[i]->ndist = NULL;
3549 p_hb[i]->n_bound = NULL;
3550 p_hb[i]->time = NULL;
3551 p_hb[i]->nhx = NULL;
3553 p_hb[i]->bHBmap = hb->bHBmap;
3554 p_hb[i]->bDAnr = hb->bDAnr;
3555 p_hb[i]->bGem = hb->bGem;
3556 p_hb[i]->wordlen = hb->wordlen;
3557 p_hb[i]->nframes = hb->nframes;
3558 p_hb[i]->maxhydro = hb->maxhydro;
3559 p_hb[i]->danr = hb->danr;
3562 p_hb[i]->hbmap = hb->hbmap;
3563 p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
3564 p_hb[i]->per = hb->per;
3566 #ifdef HAVE_NN_LOOPS
3567 p_hb[i]->hbE = hb->hbE;
3571 p_hb[i]->nrdist = 0;
3574 /* Make a thread pool here,
3575 * instead of forking anew at every frame. */
3577 #pragma omp parallel \
3578 private(i, j, h, ii, jj, hh, E, \
3579 xi, yi, zi, xj, yj, zj, threadNr, \
3580 dist, ang, peri, icell, jcell, \
3581 grp, ogrp, ai, aj, xjj, yjj, zjj, \
3582 xk, yk, zk, ihb, id, resdist, \
3583 xkk, ykk, zkk, kcell, ak, k, bTric) \
3585 shared(hb, p_hb, p_adist, p_rdist, actual_nThreads, \
3586 x, bBox, box, hbox, rcut, r2cut, rshell, \
3587 shatom, ngrid, grid, nframes, t, \
3588 bParallel, bNN, index, bMerge, bContact, \
3589 bTwo, bDA,ccut, abin, rbin, top, \
3590 bSelected, bDebug, stderr, nsel, \
3591 bGem, oenv, fnm, fpnhb, trrStatus, natoms, \
3592 status, nabin, nrbin, adist, rdist, debug)
3593 { /* Start of parallel region */
3594 threadNr = gmx_omp_get_thread_num();
3595 #endif /* HAVE_OPENMP ================================================= */
3598 bTric = bBox && TRICLINIC(box);
3601 sync_hbdata(hb, p_hb[threadNr], nframes, t);
3605 build_grid(hb,x,x[shatom], bBox,box,hbox, (rcut>r2cut)?rcut:r2cut,
3606 rshell, ngrid,grid);
3607 reset_nhbonds(&(hb->d));
3609 if (debug && bDebug)
3610 dump_grid(debug, ngrid, grid);
3612 add_frames(hb,nframes);
3613 init_hbframe(hb,nframes,output_env_conv_time(oenv,t));
3616 count_da_grid(ngrid, grid, hb->danr[nframes]);
3619 p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
3623 #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
3624 /* Loop over all atom pairs and estimate interaction energy */
3625 #ifdef HAVE_OPENMP /* ------- */
3627 #endif /* HAVE_OPENMP ------- */
3629 addFramesNN(hb, nframes);
3631 #ifdef HAVE_OPENMP /* ---------------- */
3633 #pragma omp for schedule(dynamic)
3634 #endif /* HAVE_OPENMP ---------------- */
3635 for (i=0; i<hb->d.nrd; i++)
3637 for(j=0;j<hb->a.nra; j++)
3640 h < (bContact ? 1 : hb->d.nhydro[i]);
3643 if (i==hb->d.nrd || j==hb->a.nra)
3644 gmx_fatal(FARGS, "out of bounds");
3646 /* Get the real atom ids */
3649 hh = hb->d.hydro[i][h];
3651 /* Estimate the energy from the geometry */
3652 E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d));
3653 /* Store the energy */
3654 storeHbEnergy(hb, i, j, h, E, nframes);
3658 #endif /* HAVE_NN_LOOPS */
3668 /* Do not parallelize this just yet. */
3670 for(ii=0; (ii<nsel); ii++) {
3671 int dd = index[0][i];
3672 int aa = index[0][i+2];
3673 /* int */ hh = index[0][i+1];
3674 ihb = is_hbond(hb,ii,ii,dd,aa,rcut,r2cut,ccut,x,bBox,box,
3675 hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
3678 /* add to index if not already there */
3680 add_hbond(hb,dd,aa,hh,ii,ii,nframes,bMerge,ihb,bContact,peri);
3684 } /* if (bSelected) */
3692 calcBoxProjection(box, hb->per->P);
3694 /* loop over all gridcells (xi,yi,zi) */
3695 /* Removed confusing macro, DvdS 27/12/98 */
3698 /* The outer grid loop will have to do for now. */
3699 #pragma omp for schedule(dynamic)
3701 for(xi=0; xi<ngrid[XX]; xi++)
3702 for(yi=0; (yi<ngrid[YY]); yi++)
3703 for(zi=0; (zi<ngrid[ZZ]); zi++) {
3705 /* loop over donor groups gr0 (always) and gr1 (if necessary) */
3706 for (grp=gr0; (grp <= (bTwo?gr1:gr0)); grp++) {
3707 icell=&(grid[zi][yi][xi].d[grp]);
3714 /* loop over all hydrogen atoms from group (grp)
3715 * in this gridcell (icell)
3717 for (ai=0; (ai<icell->nr); ai++) {
3718 i = icell->atoms[ai];
3720 /* loop over all adjacent gridcells (xj,yj,zj) */
3721 for(zjj = grid_loop_begin(ngrid[ZZ],zi,bTric,FALSE);
3722 zjj <= grid_loop_end(ngrid[ZZ],zi,bTric,FALSE);
3725 zj = grid_mod(zjj,ngrid[ZZ]);
3726 bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
3727 for(yjj = grid_loop_begin(ngrid[YY],yi,bTric,bEdge_yjj);
3728 yjj <= grid_loop_end(ngrid[YY],yi,bTric,bEdge_yjj);
3731 yj = grid_mod(yjj,ngrid[YY]);
3733 (yj == 0) || (yj == ngrid[YY] - 1) ||
3734 (zj == 0) || (zj == ngrid[ZZ] - 1);
3735 for(xjj = grid_loop_begin(ngrid[XX],xi,bTric,bEdge_xjj);
3736 xjj <= grid_loop_end(ngrid[XX],xi,bTric,bEdge_xjj);
3739 xj = grid_mod(xjj,ngrid[XX]);
3740 jcell=&(grid[zj][yj][xj].a[ogrp]);
3741 /* loop over acceptor atoms from other group (ogrp)
3742 * in this adjacent gridcell (jcell)
3744 for (aj=0; (aj<jcell->nr); aj++) {
3745 j = jcell->atoms[aj];
3747 /* check if this once was a h-bond */
3749 ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
3750 hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
3753 /* add to index if not already there */
3755 add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
3757 /* make angle and distance distributions */
3758 if (ihb == hbHB && !bContact) {
3760 gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
3762 __ADIST[(int)( ang/abin)]++;
3763 __RDIST[(int)(dist/rbin)]++;
3766 if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
3767 gmx_fatal(FARGS,"Invalid donor %d",i);
3768 if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
3769 gmx_fatal(FARGS,"Invalid acceptor %d",j);
3770 resdist=abs(top.atoms.atom[i].resind-
3771 top.atoms.atom[j].resind);
3772 if (resdist >= max_hx)
3774 __HBDATA->nhx[nframes][resdist]++;
3785 } /* for xi,yi,zi */
3786 } /* if (bSelected) {...} else */
3788 #ifdef HAVE_OPENMP /* ---------------------------- */
3789 /* Better wait for all threads to finnish using x[] before updating it. */
3791 #pragma omp barrier /* */
3792 #pragma omp critical /* */
3794 /* Sum up histograms and counts from p_hb[] into hb */
3796 hb->nhb[k] += p_hb[threadNr]->nhb[k];
3797 hb->ndist[k] += p_hb[threadNr]->ndist[k];
3798 for (j=0; j<max_hx; j++) /**/
3799 hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
3803 /* Here are a handful of single constructs
3804 * to share the workload a bit. The most
3805 * important one is of course the last one,
3806 * where there's a potential bottleneck in form
3808 #pragma omp single /* ++++++++++++++++, */
3809 #endif /* HAVE_OPENMP ----------------+------------*/
3811 if (hb != NULL) /* */
3813 analyse_donor_props(opt2fn_null("-don",NFILE,fnm),hb,k,t,oenv);
3816 #ifdef HAVE_OPENMP /* + */
3817 #pragma omp single /* +++ +++ */
3821 do_nhb_dist(fpnhb,hb,t);
3823 } /* if (bNN) {...} else + */
3824 #ifdef HAVE_OPENMP /* + */
3825 #pragma omp single /* +++ +++ */
3828 trrStatus = (read_next_x(oenv,status,&t,natoms,x,box));
3831 #ifdef HAVE_OPENMP /* +++++++++++++++++ */
3834 } while (trrStatus);
3837 #pragma omp critical
3839 hb->nrhb += p_hb[threadNr]->nrhb;
3840 hb->nrdist += p_hb[threadNr]->nrdist;
3842 /* Free parallel datastructures */
3843 sfree(p_hb[threadNr]->nhb);
3844 sfree(p_hb[threadNr]->ndist);
3845 sfree(p_hb[threadNr]->nhx);
3848 for (i=0; i<nabin; i++)
3849 for (j=0; j<actual_nThreads; j++)
3851 adist[i] += p_adist[j][i];
3853 for (i=0; i<=nrbin; i++)
3854 for (j=0; j<actual_nThreads; j++)
3855 rdist[i] += p_rdist[j][i];
3857 sfree(p_adist[threadNr]);
3858 sfree(p_rdist[threadNr]);
3859 } /* End of parallel region */
3864 if(nframes <2 && (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm)))
3866 gmx_fatal(FARGS,"Cannot calculate autocorrelation of life times with less than two frames");
3869 free_grid(ngrid,&grid);
3875 /* Compute maximum possible number of different hbonds */
3879 max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
3881 /* Added support for -contact below.
3882 * - Erik Marklund, May 29-31, 2006 */
3883 /* Changed contact code.
3884 * - Erik Marklund, June 29, 2006 */
3885 if (bHBmap && !bNN) {
3887 printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds");
3889 printf("Found %d different %s in trajectory\n"
3890 "Found %d different atom-pairs within %s distance\n",
3891 hb->nrhb, bContact?"contacts":"hydrogen bonds",
3892 hb->nrdist,(r2cut>0)?"second cut-off":"hydrogen bonding");
3894 /*Control the pHist.*/
3897 merge_hb(hb,bTwo,bContact);
3899 if (opt2bSet("-hbn",NFILE,fnm))
3900 dump_hbmap(hb,NFILE,fnm,bTwo,bContact,isize,index,grpnames,&top.atoms);
3902 /* Moved the call to merge_hb() to a line BEFORE dump_hbmap
3903 * to make the -hbn and -hmb output match eachother.
3904 * - Erik Marklund, May 30, 2006 */
3907 /* Print out number of hbonds and distances */
3910 fp = xvgropen(opt2fn("-num",NFILE,fnm),bContact ? "Contacts" :
3911 "Hydrogen Bonds",output_env_get_xvgr_tlabel(oenv),"Number",oenv);
3913 snew(leg[0],STRLEN);
3914 snew(leg[1],STRLEN);
3915 sprintf(leg[0],"%s",bContact?"Contacts":"Hydrogen bonds");
3916 sprintf(leg[1],"Pairs within %g nm",(r2cut>0)?r2cut:rcut);
3917 xvgr_legend(fp,2,(const char**)leg,oenv);
3921 for(i=0; (i<nframes); i++) {
3922 fprintf(fp,"%10g %10d %10d\n",hb->time[i],hb->nhb[i],hb->ndist[i]);
3923 aver_nhb += hb->nhb[i];
3924 aver_dist += hb->ndist[i];
3927 aver_nhb /= nframes;
3928 aver_dist /= nframes;
3929 /* Print HB distance distribution */
3930 if (opt2bSet("-dist",NFILE,fnm)) {
3934 for(i=0; i<nrbin; i++)
3937 fp = xvgropen(opt2fn("-dist",NFILE,fnm),
3938 "Hydrogen Bond Distribution",
3939 "Hydrogen - Acceptor Distance (nm)","",oenv);
3940 for(i=0; i<nrbin; i++)
3941 fprintf(fp,"%10g %10g\n",(i+0.5)*rbin,rdist[i]/(rbin*(real)sum));
3945 /* Print HB angle distribution */
3946 if (opt2bSet("-ang",NFILE,fnm)) {
3950 for(i=0; i<nabin; i++)
3953 fp = xvgropen(opt2fn("-ang",NFILE,fnm),
3954 "Hydrogen Bond Distribution",
3955 "Donor - Hydrogen - Acceptor Angle (\\SO\\N)","",oenv);
3956 for(i=0; i<nabin; i++)
3957 fprintf(fp,"%10g %10g\n",(i+0.5)*abin,adist[i]/(abin*(real)sum));
3961 /* Print HB in alpha-helix */
3962 if (opt2bSet("-hx",NFILE,fnm)) {
3963 fp = xvgropen(opt2fn("-hx",NFILE,fnm),
3964 "Hydrogen Bonds",output_env_get_xvgr_tlabel(oenv),"Count",oenv);
3965 xvgr_legend(fp,NRHXTYPES,hxtypenames,oenv);
3966 for(i=0; i<nframes; i++) {
3967 fprintf(fp,"%10g",hb->time[i]);
3968 for (j=0; j<max_hx; j++)
3969 fprintf(fp," %6d",hb->nhx[i][j]);
3975 printf("Average number of %s per timeframe %.3f out of %g possible\n",
3976 bContact ? "contacts" : "hbonds",
3977 bContact ? aver_dist : aver_nhb, max_nhb);
3979 /* Do Autocorrelation etc. */
3982 Added support for -contact in ac and hbm calculations below.
3983 - Erik Marklund, May 29, 2006
3987 if (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm))
3988 please_cite(stdout,"Spoel2006b");
3989 if (opt2bSet("-ac",NFILE,fnm)) {
3990 char *gemstring=NULL;
3993 params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end,
3994 gemBallistic, nBalExp, bBallisticDt);
3996 gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
3998 gemstring = strdup(gemType[hb->per->gemtype]);
3999 do_hbac(opt2fn("-ac",NFILE,fnm),hb,nDump,
4000 bMerge,bContact,fit_start,temp,r2cut>0,smooth_tail_start,oenv,
4001 params, gemstring, nThreads, NN, bBallistic, bGemFit);
4003 if (opt2bSet("-life",NFILE,fnm))
4004 do_hblife(opt2fn("-life",NFILE,fnm),hb,bMerge,bContact,oenv);
4005 if (opt2bSet("-hbm",NFILE,fnm)) {
4010 mat.ny=(bContact ? hb->nrdist : hb->nrhb);
4012 snew(mat.matrix,mat.nx);
4013 for(x=0; (x<mat.nx); x++)
4014 snew(mat.matrix[x],mat.ny);
4016 for(id=0; (id<hb->d.nrd); id++)
4017 for(ia=0; (ia<hb->a.nra); ia++) {
4018 for(hh=0; (hh<hb->maxhydro); hh++) {
4019 if (hb->hbmap[id][ia]) {
4020 if (ISHB(hb->hbmap[id][ia]->history[hh])) {
4021 /* Changed '<' into '<=' in the for-statement below.
4022 * It fixed the previously undiscovered bug that caused
4023 * the last occurance of an hbond/contact to not be
4024 * set in mat.matrix. Have a look at any old -hbm-output
4025 * and you will notice that the last column is allways empty.
4026 * - Erik Marklund May 30, 2006
4028 for(x=0; (x<=hb->hbmap[id][ia]->nframes); x++) {
4029 int nn0 = hb->hbmap[id][ia]->n0;
4030 range_check(y,0,mat.ny);
4031 mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh],x);
4038 mat.axis_x=hb->time;
4039 snew(mat.axis_y,mat.ny);
4040 for(j=0; j<mat.ny; j++)
4042 sprintf(mat.title,bContact ? "Contact Existence Map":
4043 "Hydrogen Bond Existence Map");
4044 sprintf(mat.legend,bContact ? "Contacts" : "Hydrogen Bonds");
4045 sprintf(mat.label_x,"%s",output_env_get_xvgr_tlabel(oenv));
4046 sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index");
4049 snew(mat.map,mat.nmap);
4050 for(i=0; i<mat.nmap; i++) {
4051 mat.map[i].code.c1=hbmap[i];
4052 mat.map[i].desc=hbdesc[i];
4053 mat.map[i].rgb=hbrgb[i];
4055 fp = opt2FILE("-hbm",NFILE,fnm,"w");
4056 write_xpm_m(fp, mat);
4058 for(x=0; x<mat.nx; x++)
4059 sfree(mat.matrix[x]);
4067 fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper);
4068 fprintf(stderr, "Freeing pHist for all donors...\n");
4069 for (i=0; i<hb->d.nrd; i++) {
4070 fprintf(stderr, "\r%i",i);
4071 if (hb->per->pHist[i] != NULL) {
4072 for (j=0; j<hb->a.nra; j++)
4073 clearPshift(&(hb->per->pHist[i][j]));
4074 sfree(hb->per->pHist[i]);
4077 sfree(hb->per->pHist);
4078 sfree(hb->per->p2i);
4080 fprintf(stderr, "...done.\n");
4083 #ifdef HAVE_NN_LOOPS
4093 #define USE_THIS_GROUP(j) ( (j == gr0) || (bTwo && (j == gr1)) )
4095 fp = xvgropen(opt2fn("-dan",NFILE,fnm),
4096 "Donors and Acceptors",output_env_get_xvgr_tlabel(oenv),"Count",oenv);
4097 nleg = (bTwo?2:1)*2;
4098 snew(legnames,nleg);
4100 for(j=0; j<grNR; j++)
4101 if ( USE_THIS_GROUP(j) ) {
4102 sprintf(buf,"Donors %s",grpnames[j]);
4103 legnames[i++] = strdup(buf);
4104 sprintf(buf,"Acceptors %s",grpnames[j]);
4105 legnames[i++] = strdup(buf);
4108 gmx_incons("number of legend entries");
4109 xvgr_legend(fp,nleg,(const char**)legnames,oenv);
4110 for(i=0; i<nframes; i++) {
4111 fprintf(fp,"%10g",hb->time[i]);
4112 for (j=0; (j<grNR); j++)
4113 if ( USE_THIS_GROUP(j) )
4114 fprintf(fp," %6d",hb->danr[i][j]);