X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_hbond.c;fp=src%2Fgromacs%2Fgmxana%2Fgmx_hbond.c;h=8bde9f1ac0724d8889a63ae5ef96dc37036d1acc;hb=8a9ed358b11c7181e46129f811c6f9a9805cbdd3;hp=73b17d1a05180fd8f32226f48718799e035c281c;hpb=87823c5ebe2ecab05a07cfdfd5165e53cb07c30d;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_hbond.c b/src/gromacs/gmxana/gmx_hbond.c index 73b17d1a05..8bde9f1ac0 100644 --- a/src/gromacs/gmxana/gmx_hbond.c +++ b/src/gromacs/gmxana/gmx_hbond.c @@ -49,7 +49,6 @@ #include "gromacs/fileio/tpxio.h" #include "gromacs/fileio/trxio.h" #include "gromacs/fileio/xvgr.h" -#include "gromacs/gmxana/geminate.h" #include "gromacs/gmxana/gmx_ana.h" #include "gromacs/legacyheaders/copyrite.h" #include "gromacs/legacyheaders/macros.h" @@ -66,12 +65,6 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/snprintf.h" -#include "geminate.h" - -/*#define HAVE_NN_LOOPS*/ - -typedef short int t_E; -typedef int t_EEst; #define max_hx 7 typedef int t_hx[max_hx]; #define NRHXTYPES max_hx @@ -96,9 +89,6 @@ enum { enum { noDA, ACC, DON, DA, INGROUP }; -enum { - NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR -}; static const char *grpnames[grNR] = {"0", "1", "I" }; @@ -168,39 +158,8 @@ typedef struct { h_id *nhbonds; /* The number of HBs per H at current */ } t_donors; -/* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/ -#define PSTYPE int - -typedef struct { - int len; /* The length of frame and p. */ - int *frame; /* The frames at which transitio*/ - PSTYPE *p; -} t_pShift; - -typedef struct { - /* Periodicity history. Used for the reversible geminate recombination. */ - t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap: - * pHist[d][a]. We can safely assume that the same - * periodic shift holds for all hydrogens of a da-pair. - * - * Nowadays it only stores TRANSITIONS, and not the shift at every frame. - * That saves a LOT of memory, an hopefully kills a mysterious bug where - * pHist gets contaminated. */ - - PSTYPE nper; /* The length of p2i */ - ivec *p2i; /* Maps integer to periodic shift for a pair.*/ - matrix P; /* Projection matrix to find the box shifts. */ - int gemtype; /* enumerated type */ -} t_gemPeriod; - -typedef struct { - int nframes; - int *Etot; /* Total energy for each frame */ - t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */ -} t_hbEmap; - typedef struct { - gmx_bool bHBmap, bDAnr, bGem; + gmx_bool bHBmap, bDAnr; int wordlen; /* The following arrays are nframes long */ int nframes, max_frames, maxhydro; @@ -215,128 +174,15 @@ typedef struct { /* This holds a matrix with all possible hydrogen bonds */ int nrhb, nrdist; t_hbond ***hbmap; -#ifdef HAVE_NN_LOOPS - t_hbEmap hbE; -#endif - /* For parallelization reasons this will have to be a pointer. - * Otherwise discrepancies may arise between the periodicity data - * seen by different threads. */ - t_gemPeriod *per; } t_hbdata; -static void clearPshift(t_pShift *pShift) -{ - if (pShift->len > 0) - { - sfree(pShift->p); - sfree(pShift->frame); - pShift->len = 0; - } -} - -static void calcBoxProjection(matrix B, matrix P) -{ - const int vp[] = {XX, YY, ZZ}; - int i, j; - int m, n; - matrix M, N, U; - - for (i = 0; i < 3; i++) - { - m = vp[i]; - for (j = 0; j < 3; j++) - { - n = vp[j]; - U[m][n] = i == j ? 1 : 0; - } - } - m_inv(B, M); - for (i = 0; i < 3; i++) - { - m = vp[i]; - mvmul(M, U[m], P[m]); - } - transpose(P, N); -} - -static void calcBoxDistance(matrix P, rvec d, ivec ibd) -{ - /* returns integer distance in box coordinates. - * P is the projection matrix from cartesian coordinates - * obtained with calcBoxProjection(). */ - int i; - rvec bd; - mvmul(P, d, bd); - /* extend it by 0.5 in all directions since (int) rounds toward 0.*/ - for (i = 0; i < 3; i++) - { - bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5); - } - ibd[XX] = (int)bd[XX]; - ibd[YY] = (int)bd[YY]; - ibd[ZZ] = (int)bd[ZZ]; -} - /* Changed argument 'bMerge' into 'oneHB' below, * since -contact should cause maxhydro to be 1, * not just -merge. * - Erik Marklund May 29, 2006 */ -static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap) -{ - /* Try to merge hbonds on the fly. That means that if the - * acceptor and donor are mergable, then: - * 1) store the hb-info so that acceptor id > donor id, - * 2) add the periodic shift in pairs, so that [-x,-y,-z] is - * stored in per.p2i[] whenever acceptor id < donor id. - * Note that [0,0,0] should already be the first element of per.p2i - * by the time this function is called. */ - - /* daSwap is TRUE if the donor and acceptor were swapped. - * If so, then the negative vector should be used. */ - PSTYPE i; - - if (per->p2i == NULL || per->nper == 0) - { - gmx_fatal(FARGS, "'per' not initialized properly."); - } - for (i = 0; i < per->nper; i++) - { - if (r[XX] == per->p2i[i][XX] && - r[YY] == per->p2i[i][YY] && - r[ZZ] == per->p2i[i][ZZ]) - { - return i; - } - } - /* Not found apparently. Add it to the list! */ - /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */ - -#pragma omp critical - { - if (!per->p2i) - { - fprintf(stderr, "p2i not initialized. This shouldn't happen!\n"); - snew(per->p2i, 1); - } - else - { - srenew(per->p2i, per->nper+2); - } - copy_ivec(r, per->p2i[per->nper]); - (per->nper)++; - - /* Add the mirror too. It's rather likely that it'll be needed. */ - per->p2i[per->nper][XX] = -r[XX]; - per->p2i[per->nper][YY] = -r[YY]; - per->p2i[per->nper][ZZ] = -r[ZZ]; - (per->nper)++; - } /* omp critical */ - return per->nper - 1 - (daSwap ? 0 : 1); -} - -static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode) +static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB) { t_hbdata *hb; @@ -344,7 +190,6 @@ static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_ hb->wordlen = 8*sizeof(unsigned int); hb->bHBmap = bHBmap; hb->bDAnr = bDAnr; - hb->bGem = bGem; if (oneHB) { hb->maxhydro = 1; @@ -353,9 +198,6 @@ static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_ { hb->maxhydro = MAXHYDRO; } - snew(hb->per, 1); - hb->per->gemtype = bGem ? gemmode : 0; - return hb; } @@ -378,295 +220,6 @@ static void mk_hbmap(t_hbdata *hb) } } -/* Consider redoing pHist so that is only stores transitions between - * periodicities and not the periodicity for all frames. This eats heaps of memory. */ -static void mk_per(t_hbdata *hb) -{ - int i, j; - if (hb->bGem) - { - snew(hb->per->pHist, hb->d.nrd); - for (i = 0; i < hb->d.nrd; i++) - { - snew(hb->per->pHist[i], hb->a.nra); - if (hb->per->pHist[i] == NULL) - { - gmx_fatal(FARGS, "Could not allocate enough memory for per->pHist"); - } - for (j = 0; j < hb->a.nra; j++) - { - clearPshift(&(hb->per->pHist[i][j])); - } - } - /* add the [0,0,0] shift to element 0 of p2i. */ - snew(hb->per->p2i, 1); - clear_ivec(hb->per->p2i[0]); - hb->per->nper = 1; - } -} - -#ifdef HAVE_NN_LOOPS -static void mk_hbEmap (t_hbdata *hb, int n0) -{ - int i, j, k; - hb->hbE.E = NULL; - hb->hbE.nframes = 0; - snew(hb->hbE.E, hb->d.nrd); - for (i = 0; i < hb->d.nrd; i++) - { - snew(hb->hbE.E[i], hb->a.nra); - for (j = 0; j < hb->a.nra; j++) - { - snew(hb->hbE.E[i][j], MAXHYDRO); - for (k = 0; k < MAXHYDRO; k++) - { - hb->hbE.E[i][j][k] = NULL; - } - } - } - hb->hbE.Etot = NULL; -} - -static void free_hbEmap (t_hbdata *hb) -{ - int i, j, k; - for (i = 0; i < hb->d.nrd; i++) - { - for (j = 0; j < hb->a.nra; j++) - { - for (k = 0; k < MAXHYDRO; k++) - { - sfree(hb->hbE.E[i][j][k]); - } - sfree(hb->hbE.E[i][j]); - } - sfree(hb->hbE.E[i]); - } - sfree(hb->hbE.E); - sfree(hb->hbE.Etot); -} - -static void addFramesNN(t_hbdata *hb, int frame) -{ - -#define DELTAFRAMES_HBE 10 - - int d, a, h, nframes; - - if (frame >= hb->hbE.nframes) - { - nframes = hb->hbE.nframes + DELTAFRAMES_HBE; - srenew(hb->hbE.Etot, nframes); - - for (d = 0; d < hb->d.nrd; d++) - { - for (a = 0; a < hb->a.nra; a++) - { - for (h = 0; h < hb->d.nhydro[d]; h++) - { - srenew(hb->hbE.E[d][a][h], nframes); - } - } - } - - hb->hbE.nframes += DELTAFRAMES_HBE; - } -} - -static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst, - matrix box, rvec hbox, t_donors *donors) -{ - /* d - donor atom - * a - acceptor atom - * h - hydrogen - * alpha - angle between dipoles - * x[] - atomic positions - * EEst - the type of energy estimate (see enum in hbplugin.h) - * box - the box vectors \ - * hbox - half box lengths _These two are only needed for the pbc correction - */ - - t_E E; - rvec dist; - rvec dipole[2], xmol[3], xmean[2]; - int i; - real r, realE; - - if (d == a) - { - /* Self-interaction */ - return NONSENSE_E; - } - - switch (EEst) - { - case NN_BINARY: - /* This is a simple binary existence function that sets E=1 whenever - * the distance between the oxygens is equal too or less than 0.35 nm. - */ - rvec_sub(x[d], x[a], dist); - pbc_correct_gem(dist, box, hbox); - if (norm(dist) <= 0.35) - { - E = 1; - } - else - { - E = 0; - } - break; - - case NN_1_over_r3: - /* Negative potential energy of a dipole. - * E = -cos(alpha) * 1/r^3 */ - - copy_rvec(x[d], xmol[0]); /* donor */ - copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ - - svmul(15.9994*(1/1.008), xmol[0], xmean[0]); - rvec_inc(xmean[0], xmol[1]); - rvec_inc(xmean[0], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - /* Assumes that all acceptors are also donors. */ - copy_rvec(x[a], xmol[0]); /* acceptor */ - copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */ - - - svmul(15.9994*(1/1.008), xmol[0], xmean[1]); - rvec_inc(xmean[1], xmol[1]); - rvec_inc(xmean[1], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - rvec_sub(xmean[0], xmean[1], dist); - pbc_correct_gem(dist, box, hbox); - r = norm(dist); - - realE = pow(r, -3.0); - E = (t_E)(SCALEFACTOR_E * realE); - break; - - case NN_dipole: - /* Negative potential energy of a (unpolarizable) dipole. - * E = -cos(alpha) * 1/r^3 */ - clear_rvec(dipole[1]); - clear_rvec(dipole[0]); - - copy_rvec(x[d], xmol[0]); /* donor */ - copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ - - rvec_inc(dipole[0], xmol[1]); - rvec_inc(dipole[0], xmol[2]); - for (i = 0; i < 3; i++) - { - dipole[0][i] *= 0.5; - } - rvec_dec(dipole[0], xmol[0]); - - svmul(15.9994*(1/1.008), xmol[0], xmean[0]); - rvec_inc(xmean[0], xmol[1]); - rvec_inc(xmean[0], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - /* Assumes that all acceptors are also donors. */ - copy_rvec(x[a], xmol[0]); /* acceptor */ - copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */ - - - rvec_inc(dipole[1], xmol[1]); - rvec_inc(dipole[1], xmol[2]); - for (i = 0; i < 3; i++) - { - dipole[1][i] *= 0.5; - } - rvec_dec(dipole[1], xmol[0]); - - svmul(15.9994*(1/1.008), xmol[0], xmean[1]); - rvec_inc(xmean[1], xmol[1]); - rvec_inc(xmean[1], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - rvec_sub(xmean[0], xmean[1], dist); - pbc_correct_gem(dist, box, hbox); - r = norm(dist); - - double cosalpha = cos_angle(dipole[0], dipole[1]); - realE = cosalpha * pow(r, -3.0); - E = (t_E)(SCALEFACTOR_E * realE); - break; - - default: - printf("Can't do that type of energy estimate: %i\n.", EEst); - E = NONSENSE_E; - } - - return E; -} - -static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame) -{ - /* hb - hbond data structure - d - donor - a - acceptor - h - hydrogen - E - estimate of the energy - frame - the current frame. - */ - - /* Store the estimated energy */ - if (E == NONSENSE_E) - { - E = 0; - } - - hb->hbE.E[d][a][h][frame] = E; - -#pragma omp critical - { - hb->hbE.Etot[frame] += E; - } -} -#endif /* HAVE_NN_LOOPS */ - - -/* Finds -v[] in the periodicity index */ -static int findMirror(PSTYPE p, ivec v[], PSTYPE nper) -{ - PSTYPE i; - ivec u; - for (i = 0; i < nper; i++) - { - if (v[i][XX] == -(v[p][XX]) && - v[i][YY] == -(v[p][YY]) && - v[i][ZZ] == -(v[p][ZZ])) - { - return (int)i; - } - } - printf("Couldn't find mirror of [%i, %i, %i], index \n", - v[p][XX], - v[p][YY], - v[p][ZZ]); - return -1; -} - - static void add_frames(t_hbdata *hb, int nframes) { int i, j, k, l; @@ -727,64 +280,13 @@ static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb) _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE); } -static void addPshift(t_pShift *pHist, PSTYPE p, int frame) -{ - if (pHist->len == 0) - { - snew(pHist->frame, 1); - snew(pHist->p, 1); - pHist->len = 1; - pHist->frame[0] = frame; - pHist->p[0] = p; - return; - } - else - if (pHist->p[pHist->len-1] != p) - { - pHist->len++; - srenew(pHist->frame, pHist->len); - srenew(pHist->p, pHist->len); - pHist->frame[pHist->len-1] = frame; - pHist->p[pHist->len-1] = p; - } /* Otherwise, there is no transition. */ - return; -} - -static PSTYPE getPshift(t_pShift pHist, int frame) -{ - int f, i; - - if (pHist.len == 0 - || (pHist.len > 0 && pHist.frame[0] > frame)) - { - return -1; - } - - for (i = 0; i < pHist.len; i++) - { - f = pHist.frame[i]; - if (f == frame) - { - return pHist.p[i]; - } - if (f > frame) - { - return pHist.p[i-1]; - } - } - - /* It seems that frame is after the last periodic transition. Return the last periodicity. */ - return pHist.p[pHist.len-1]; -} - -static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPE p) +static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb) { int i, j, n; t_hbond *hb = hbd->hbmap[id][ia]; int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]); int wlen = hbd->wordlen; int delta = 32*wlen; - gmx_bool bGem = hbd->bGem; if (!hb->h[0]) { @@ -822,18 +324,6 @@ static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PST if (frame >= 0) { set_hb(hbd, id, h, ia, frame, ihb); - if (bGem) - { - if (p >= hbd->per->nper) - { - gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper); - } - else - { - addPshift(&(hbd->per->pHist[id][ia]), p, frame); - } - - } } } @@ -917,7 +407,7 @@ static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd) static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa, - int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPE p) + int frame, gmx_bool bMerge, int ihb, gmx_bool bContact) { int k, id, ia, hh; gmx_bool daSwap = FALSE; @@ -1011,7 +501,7 @@ static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa, snew(hb->hbmap[id][ia]->h, hb->maxhydro); snew(hb->hbmap[id][ia]->g, hb->maxhydro); } - add_ff(hb, id, k, ia, frame, ihb, p); + add_ff(hb, id, k, ia, frame, ihb); } } @@ -1503,52 +993,30 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, rvec_sub(x[ad[i]], xshell, dshell); if (bBox) { - if (FALSE && !hb->bGem) + gmx_bool bDone = FALSE; + while (!bDone) { + bDone = TRUE; for (m = DIM-1; m >= 0 && bInShell; m--) { if (dshell[m] < -hbox[m]) { + bDone = FALSE; rvec_inc(dshell, box[m]); } - else if (dshell[m] >= hbox[m]) + if (dshell[m] >= hbox[m]) { + bDone = FALSE; dshell[m] -= 2*hbox[m]; } - /* if we're outside the cube, we're outside the sphere also! */ - if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) - { - bInShell = FALSE; - } } } - else + for (m = DIM-1; m >= 0 && bInShell; m--) { - gmx_bool bDone = FALSE; - while (!bDone) - { - bDone = TRUE; - for (m = DIM-1; m >= 0 && bInShell; m--) - { - if (dshell[m] < -hbox[m]) - { - bDone = FALSE; - rvec_inc(dshell, box[m]); - } - if (dshell[m] >= hbox[m]) - { - bDone = FALSE; - dshell[m] -= 2*hbox[m]; - } - } - } - for (m = DIM-1; m >= 0 && bInShell; m--) + /* if we're outside the cube, we're outside the sphere also! */ + if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) { - /* if we're outside the cube, we're outside the sphere also! */ - if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) - { - bInShell = FALSE; - } + bInShell = FALSE; } } } @@ -1563,10 +1031,6 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, { if (bBox) { - if (hb->bGem) - { - copy_rvec(x[ad[i]], xtemp); - } pbc_in_gridbox(x[ad[i]], box); for (m = DIM-1; m >= 0; m--) @@ -1574,10 +1038,6 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, grididx[m] = x[ad[i]][m]*invdelta[m]; grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m]; } - if (hb->bGem) - { - copy_rvec(xtemp, x[ad[i]]); /* copy back */ - } } gx = grididx[XX]; @@ -1774,7 +1234,7 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, real rcut, real r2cut, real ccut, rvec x[], gmx_bool bBox, matrix box, rvec hbox, real *d_ha, real *ang, gmx_bool bDA, int *hhh, - gmx_bool bContact, gmx_bool bMerge, PSTYPE *p) + gmx_bool bContact, gmx_bool bMerge) { int h, hh, id, ja, ihb; rvec r_da, r_ha, r_dh, r = {0, 0, 0}; @@ -1811,15 +1271,7 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, { /* return hbNo; */ daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */ } - if (hb->bGem) - { - copy_rvec(r_da, r); /* Save this for later */ - pbc_correct_gem(r_da, box, hbox); - } - else - { - pbc_correct_gem(r_da, box, hbox); - } + pbc_correct_gem(r_da, box, hbox); } rda2 = iprod(r_da, r_da); @@ -1831,11 +1283,6 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, } if (rda2 <= rc2) { - if (hb->bGem) - { - calcBoxDistance(hb->per->P, r, ri); - *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */ - } return hbHB; } else if (rda2 < r2c2) @@ -1868,12 +1315,6 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, rha2 = iprod(r_ha, r_ha); } - if (hb->bGem) - { - calcBoxDistance(hb->per->P, r, ri); - *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */ - } - if (bDA || (!bDA && (rha2 <= rc2))) { rvec_sub(x[d], x[hh], r_dh); @@ -1907,28 +1348,17 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, } } -/* Fixed previously undiscovered bug in the merge - code, where the last frame of each hbond disappears. - - Erik Marklund, June 1, 2006 */ -/* Added the following arguments: - * ptmp[] - temporary periodicity hisory - * a1 - identity of first acceptor/donor - * a2 - identity of second acceptor/donor - * - Erik Marklund, FEB 20 2010 */ - /* Merging is now done on the fly, so do_merge is most likely obsolete now. * Will do some more testing before removing the function entirely. * - Erik Marklund, MAY 10 2010 */ static void do_merge(t_hbdata *hb, int ntmp, - unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[], - t_hbond *hb0, t_hbond *hb1, int a1, int a2) + unsigned int htmp[], unsigned int gtmp[], + t_hbond *hb0, t_hbond *hb1) { /* Here we need to make sure we're treating periodicity in * the right way for the geminate recombination kinetics. */ int m, mm, n00, n01, nn0, nnframes; - PSTYPE pm; - t_pShift *pShift; /* Decide where to start from when merging */ n00 = hb0->n0; @@ -1940,7 +1370,6 @@ static void do_merge(t_hbdata *hb, int ntmp, { htmp[m] = 0; gtmp[m] = 0; - ptmp[m] = 0; } /* Fill tmp arrays with values due to first HB */ /* Once again '<' had to be replaced with '<=' @@ -1951,28 +1380,11 @@ static void do_merge(t_hbdata *hb, int ntmp, { mm = m+n00-nn0; htmp[mm] = is_hb(hb0->h[0], m); - if (hb->bGem) - { - pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0); - if (pm > hb->per->nper) - { - gmx_fatal(FARGS, "Illegal shift!"); - } - else - { - ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/ - } - } } - /* If we're doing geminate recompbination we usually don't need the distances. - * Let's save some memory and time. */ - if (TRUE || !hb->bGem || hb->per->gemtype == gemAD) + for (m = 0; (m <= hb0->nframes); m++) { - for (m = 0; (m <= hb0->nframes); m++) - { - mm = m+n00-nn0; - gtmp[mm] = is_hb(hb0->g[0], m); - } + mm = m+n00-nn0; + gtmp[mm] = is_hb(hb0->g[0], m); } /* Next HB */ for (m = 0; (m <= hb1->nframes); m++) @@ -1980,20 +1392,6 @@ static void do_merge(t_hbdata *hb, int ntmp, mm = m+n01-nn0; htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m); gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m); - if (hb->bGem /* && ptmp[mm] != 0 */) - { - - /* If this hbond has been seen before with donor and acceptor swapped, - * then we need to find the mirrored (*-1) periodicity vector to truely - * merge the hbond history. */ - pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper); - /* Store index of mirror */ - if (pm > hb->per->nper) - { - gmx_fatal(FARGS, "Illegal shift!"); - } - ptmp[mm] = pm; - } } /* Reallocate target array */ if (nnframes > hb0->maxframes) @@ -2001,20 +1399,12 @@ static void do_merge(t_hbdata *hb, int ntmp, srenew(hb0->h[0], 4+nnframes/hb->wordlen); srenew(hb0->g[0], 4+nnframes/hb->wordlen); } - if (NULL != hb->per->pHist) - { - clearPshift(&(hb->per->pHist[a1][a2])); - } /* Copy temp array to target array */ for (m = 0; (m <= nnframes); m++) { _set_hb(hb0->h[0], m, htmp[m]); _set_hb(hb0->g[0], m, gtmp[m]); - if (hb->bGem) - { - addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0); - } } /* Set scalar variables */ @@ -2022,14 +1412,10 @@ static void do_merge(t_hbdata *hb, int ntmp, hb0->maxframes = nnframes; } -/* Added argument bContact for nicer output. - * Erik Marklund, June 29, 2006 - */ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) { int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp; unsigned int *htmp, *gtmp; - PSTYPE *ptmp; t_hbond *hb0, *hb1; inrnew = hb->nrhb; @@ -2041,7 +1427,6 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) ntmp = 2*hb->max_frames; snew(gtmp, ntmp); snew(htmp, ntmp); - snew(ptmp, ntmp); for (i = 0; (i < hb->d.nrd); i++) { fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd); @@ -2058,7 +1443,7 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) hb1 = hb->hbmap[jj][ii]; if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0])) { - do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j); + do_merge(hb, ntmp, htmp, gtmp, hb0, hb1); if (ISHB(hb1->history[0])) { inrnew--; @@ -2078,10 +1463,6 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) } sfree(hb1->h[0]); sfree(hb1->g[0]); - if (hb->bGem) - { - clearPshift(&(hb->per->pHist[jj][ii])); - } hb1->h[0] = NULL; hb1->g[0] = NULL; hb1->history[0] = hbNo; @@ -2096,7 +1477,6 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) hb->nrdist = indnew; sfree(gtmp); sfree(htmp); - sfree(ptmp); } static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t) @@ -2128,12 +1508,6 @@ static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t) fprintf(fp, " %8d\n", nbtot); } -/* Added argument bContact in do_hblife(...). Also - * added support for -contact in function body. - * - Erik Marklund, May 31, 2006 */ -/* Changed the contact code slightly. - * - Erik Marklund, June 29, 2006 - */ static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact, const output_env_t oenv) { @@ -2185,10 +1559,6 @@ static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bC ohb = 0; j0 = 0; - /* Changed '<' into '<=' below, just like I - did in the hbm-output-loop in the main code. - - Erik Marklund, May 31, 2006 - */ for (j = 0; (j <= hbh->nframes); j++) { ihb = is_hb(h[nh], j); @@ -2256,9 +1626,6 @@ static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bC sfree(histo); } -/* Changed argument bMerge into oneHB to handle contacts properly. - * - Erik Marklund, June 29, 2006 - */ static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump) { FILE *fp; @@ -2549,17 +1916,10 @@ static void normalizeACF(real *ct, real *gt, int nhb, int len) } } -/* Added argument bContact in do_hbac(...). Also - * added support for -contact in the actual code. - * - Erik Marklund, May 31, 2006 */ -/* Changed contact code and added argument R2 - * - Erik Marklund, June 29, 2006 - */ static void do_hbac(const char *fn, t_hbdata *hb, int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start, real temp, gmx_bool R2, const output_env_t oenv, - const char *gemType, int nThreads, - const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit) + int nThreads) { FILE *fp; int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets; @@ -2585,13 +1945,9 @@ static void do_hbac(const char *fn, t_hbdata *hb, unsigned int **h = NULL, **g = NULL; int nh, nhbonds, nhydro, ngh; t_hbond *hbh; - PSTYPE p, *pfound = NULL, np; - t_pShift *pHist; int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN; - real **rHbExGem = NULL; gmx_bool c; int acType; - t_E *E; double *ctdouble, *timedouble, *fittedct; double fittolerance = 0.1; int *dondata = NULL, thisThread; @@ -2608,44 +1964,8 @@ static void do_hbac(const char *fn, t_hbdata *hb, printf("Doing autocorrelation "); - /* Decide what kind of ACF calculations to do. */ - if (NN > NN_NONE && NN < NN_NR) - { -#ifdef HAVE_NN_LOOPS - acType = AC_NN; - printf("using the energy estimate.\n"); -#else - acType = AC_NONE; - printf("Can't do the NN-loop. Yet.\n"); -#endif - } - else if (hb->bGem) - { - acType = AC_GEM; - printf("according to the reversible geminate recombination model by Omer Markowitch.\n"); - - nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0); - snew(legGem, nSets); - for (i = 0; i < nSets; i++) - { - snew(legGem[i], 128); - } - sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType); - if (bBallistic) - { - sprintf(legGem[1], "Ac'(t)"); - } - if (bGemFit) - { - sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType); - } - - } - else - { - acType = AC_LUZAR; - printf("according to the theory of Luzar and Chandler.\n"); - } + acType = AC_LUZAR; + printf("according to the theory of Luzar and Chandler.\n"); fflush(stdout); /* build hbexist matrix in reals for autocorr */ @@ -2698,581 +2018,173 @@ static void do_hbac(const char *fn, t_hbdata *hb, } - /* Build the ACF according to acType */ - switch (acType) - { - - case AC_NN: -#ifdef HAVE_NN_LOOPS - /* Here we're using the estimated energy for the hydrogen bonds. */ - snew(ct, nn); - -#pragma omp parallel \ - private(i, j, k, nh, E, rhbex, thisThread) \ - default(shared) - { -#pragma omp barrier - thisThread = gmx_omp_get_thread_num(); - rhbex = NULL; - - snew(rhbex, n2); - memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */ + /* Build the ACF */ + snew(rhbex, 2*n2); + snew(ct, 2*n2); + snew(gt, 2*n2); + snew(ht, 2*n2); + snew(ght, 2*n2); + snew(dght, 2*n2); -#pragma omp barrier -#pragma omp for schedule (dynamic) - for (i = 0; i < hb->d.nrd; i++) /* loop over donors */ - { - if (bOMP) - { -#pragma omp critical - { - dondata[thisThread] = i; - parallel_print(dondata, nThreads); - } - } - else - { - fprintf(stderr, "\r %i", i); - } + snew(kt, nn); + snew(cct, nn); - for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */ - { - for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */ - { - E = hb->hbE.E[i][j][nh]; - if (E != NULL) - { - for (k = 0; k < nframes; k++) - { - if (E[k] != NONSENSE_E) - { - rhbex[k] = (real)E[k]; - } - } - - low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0], - eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1); -#pragma omp critical - { - for (k = 0; (k < nn); k++) - { - ct[k] += rhbex[k]; - } - } - } - } /* k loop */ - } /* j loop */ - } /* i loop */ - sfree(rhbex); -#pragma omp barrier - } - - if (bOMP) - { - sfree(dondata); - } - normalizeACF(ct, NULL, 0, nn); - snew(ctdouble, nn); - snew(timedouble, nn); - for (j = 0; j < nn; j++) - { - timedouble[j] = (double)(hb->time[j]); - ctdouble[j] = (double)(ct[j]); - } - - /* Remove ballistic term */ - /* Ballistic component removal and fitting to the reversible geminate recombination model - * will be taken out for the time being. First of all, one can remove the ballistic - * component with g_analyze afterwards. Secondly, and more importantly, there are still - * problems with the robustness of the fitting to the model. More work is needed. - * A third reason is that we're currently using gsl for this and wish to reduce dependence - * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with - * a BSD-licence that can do the job. - * - * - Erik Marklund, June 18 2010. - */ -/* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ -/* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ -/* else */ -/* printf("\nNumber of data points is less than the number of parameters to fit\n." */ -/* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ - - fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)"); - xvgr_legend(fp, asize(legNN), legNN); + for (i = 0; (i < hb->d.nrd); i++) + { + for (k = 0; (k < hb->a.nra); k++) + { + nhydro = 0; + hbh = hb->hbmap[i][k]; - for (j = 0; (j < nn); j++) + if (hbh) { - fprintf(fp, "%10g %10g %10g\n", - hb->time[j]-hb->time[0], - ct[j], - ctdouble[j]); - } - xvgrclose(fp); - sfree(ct); - sfree(ctdouble); - sfree(timedouble); -#endif /* HAVE_NN_LOOPS */ - break; /* case AC_NN */ - - case AC_GEM: - snew(ct, 2*n2); - memset(ct, 0, 2*n2*sizeof(real)); -#ifndef GMX_OPENMP - fprintf(stderr, "Donor:\n"); -#define __ACDATA ct -#else -#define __ACDATA p_ct -#endif - -#pragma omp parallel \ - private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \ - pfound, poff, rHbExGem, p, ihb, mMax, \ - thisThread, p_ct) \ - default(shared) - { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */ - h = NULL; - g = NULL; - thisThread = gmx_omp_get_thread_num(); - snew(h, hb->maxhydro); - snew(g, hb->maxhydro); - mMax = INT_MIN; - rHbExGem = NULL; - poff = NULL; - pfound = NULL; - p_ct = NULL; - snew(p_ct, 2*n2); - memset(p_ct, 0, 2*n2*sizeof(real)); - - /* I'm using a chunk size of 1, since I expect \ - * the overhead to be really small compared \ - * to the actual calculations \ */ -#pragma omp for schedule(dynamic,1) nowait - for (i = 0; i < hb->d.nrd; i++) + if (bMerge || bContact) { - - if (bOMP) + if (ISHB(hbh->history[0])) { -#pragma omp critical - { - dondata[thisThread] = i; - parallel_print(dondata, nThreads); - } - } - else - { - fprintf(stderr, "\r %i", i); + h[0] = hbh->h[0]; + g[0] = hbh->g[0]; + nhydro = 1; } - for (k = 0; k < hb->a.nra; k++) - { - for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) - { - hbh = hb->hbmap[i][k]; - if (hbh) - { - /* Note that if hb->per->gemtype==gemDD, then distances will be stored in - * hb->hbmap[d][a].h array anyway, because the contact flag will be set. - * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */ - pHist = &(hb->per->pHist[i][k]); - if (ISHB(hbh->history[nh]) && pHist->len != 0) - { - - { - h[nh] = hbh->h[nh]; - g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL; - } - n0 = hbh->n0; - nf = hbh->nframes; - /* count the number of periodic shifts encountered and store - * them in separate arrays. */ - np = 0; - for (j = 0; j < pHist->len; j++) - { - p = pHist->p[j]; - for (m = 0; m <= np; m++) - { - if (m == np) /* p not recognized in list. Add it and set up new array. */ - { - np++; - if (np > hb->per->nper) - { - gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here."); - } - if (m >= mMax) /* Extend the arrays. - * Doing it like this, using mMax to keep track of the sizes, - * eleviates the need for freeing and re-allocating the arrays - * when taking on the next donor-acceptor pair */ - { - mMax = m; - srenew(pfound, np); /* The list of found periodic shifts. */ - srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */ - snew(rHbExGem[m], 2*n2); - srenew(poff, np); - } - - { - if (rHbExGem != NULL && rHbExGem[m] != NULL) - { - /* This must be done, as this array was most likey - * used to store stuff in some previous iteration. */ - memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2)); - } - else - { - fprintf(stderr, "rHbExGem not initialized! m = %i\n", m); - } - } - pfound[m] = p; - poff[m] = -1; - - break; - } /* m==np */ - if (p == pfound[m]) - { - break; - } - } /* m: Loop over found shifts */ - } /* j: Loop over shifts */ - - /* Now unpack and disentangle the existence funtions. */ - for (j = 0; j < nf; j++) - { - /* i: donor, - * k: acceptor - * nh: hydrogen - * j: time - * p: periodic shift - * pfound: list of periodic shifts found for this pair. - * poff: list of frame offsets; that is, the first - * frame a hbond has a particular periodic shift. */ - p = getPshift(*pHist, j+n0); - if (p != -1) - { - for (m = 0; m < np; m++) - { - if (pfound[m] == p) - { - break; - } - if (m == (np-1)) - { - gmx_fatal(FARGS, "Shift not found, but must be there."); - } - } - - ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j)); - if (ihb) - { - if (poff[m] == -1) - { - poff[m] = j; /* Here's where the first hbond with shift p is, - * relative to the start of h[0].*/ - } - if (j < poff[m]) - { - gmx_fatal(FARGS, "j 0 && n0+poff[m] < nn /* && m==0 */) - { - low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0], - eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0); - for (j = 0; (j < nn); j++) - { - __ACDATA[j] += rHbExGem[m][j]; - } - } - } /* Building of ac. */ - } /* if (ISHB(...*/ - } /* if (hbh) */ - } /* hydrogen loop */ - } /* acceptor loop */ - } /* donor loop */ - - for (m = 0; m <= mMax; m++) - { - sfree(rHbExGem[m]); } - sfree(pfound); - sfree(poff); - sfree(rHbExGem); - - sfree(h); - sfree(g); - - if (bOMP) + else { -#pragma omp critical + for (m = 0; (m < hb->maxhydro); m++) { - for (i = 0; i < nn; i++) + if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) { - ct[i] += p_ct[i]; + g[nhydro] = hbh->g[m]; + h[nhydro] = hbh->h[m]; + nhydro++; } } - sfree(p_ct); } - } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */ - if (bOMP) - { - sfree(dondata); - } - - normalizeACF(ct, NULL, 0, nn); - - fprintf(stderr, "\n\nACF successfully calculated.\n"); - - /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */ - - snew(ctdouble, nn); - snew(timedouble, nn); - snew(fittedct, nn); - - for (j = 0; j < nn; j++) - { - timedouble[j] = (double)(hb->time[j]); - ctdouble[j] = (double)(ct[j]); - } - - /* Remove ballistic term */ - /* Ballistic component removal and fitting to the reversible geminate recombination model - * will be taken out for the time being. First of all, one can remove the ballistic - * component with g_analyze afterwards. Secondly, and more importantly, there are still - * problems with the robustness of the fitting to the model. More work is needed. - * A third reason is that we're currently using gsl for this and wish to reduce dependence - * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with - * a BSD-licence that can do the job. - * - * - Erik Marklund, June 18 2010. - */ -/* if (bBallistic) { */ -/* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ -/* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ -/* else */ -/* printf("\nNumber of data points is less than the number of parameters to fit\n." */ -/* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ -/* } */ -/* if (bGemFit) */ -/* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */ - - - if (bContact) - { - fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - else - { - fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv); - - for (j = 0; (j < nn); j++) - { - fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]); - if (bBallistic) - { - fprintf(fp, " %10g", ctdouble[j]); - } - if (bGemFit) - { - fprintf(fp, " %10g", fittedct[j]); - } - fprintf(fp, "\n"); - } - xvgrclose(fp); - - sfree(ctdouble); - sfree(timedouble); - sfree(fittedct); - sfree(ct); - - break; /* case AC_GEM */ - - case AC_LUZAR: - snew(rhbex, 2*n2); - snew(ct, 2*n2); - snew(gt, 2*n2); - snew(ht, 2*n2); - snew(ght, 2*n2); - snew(dght, 2*n2); - - snew(kt, nn); - snew(cct, nn); - - for (i = 0; (i < hb->d.nrd); i++) - { - for (k = 0; (k < hb->a.nra); k++) + nf = hbh->nframes; + for (nh = 0; (nh < nhydro); nh++) { - nhydro = 0; - hbh = hb->hbmap[i][k]; - - if (hbh) + int nrint = bContact ? hb->nrdist : hb->nrhb; + if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint)) + { + fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint); + } + nhbonds++; + for (j = 0; (j < nframes); j++) { - if (bMerge || bContact) + if (j <= nf) { - if (ISHB(hbh->history[0])) - { - h[0] = hbh->h[0]; - g[0] = hbh->g[0]; - nhydro = 1; - } + ihb = is_hb(h[nh], j); + idist = is_hb(g[nh], j); } else { - for (m = 0; (m < hb->maxhydro); m++) - { - if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) - { - g[nhydro] = hbh->g[m]; - h[nhydro] = hbh->h[m]; - nhydro++; - } - } + ihb = idist = 0; } - - nf = hbh->nframes; - for (nh = 0; (nh < nhydro); nh++) + rhbex[j] = ihb; + /* For contacts: if a second cut-off is provided, use it, + * otherwise use g(t) = 1-h(t) */ + if (!R2 && bContact) { - int nrint = bContact ? hb->nrdist : hb->nrhb; - if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint)) - { - fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint); - } - nhbonds++; - for (j = 0; (j < nframes); j++) - { - /* Changed '<' into '<=' below, just like I did in - the hbm-output-loop in the gmx_hbond() block. - - Erik Marklund, May 31, 2006 */ - if (j <= nf) - { - ihb = is_hb(h[nh], j); - idist = is_hb(g[nh], j); - } - else - { - ihb = idist = 0; - } - rhbex[j] = ihb; - /* For contacts: if a second cut-off is provided, use it, - * otherwise use g(t) = 1-h(t) */ - if (!R2 && bContact) - { - gt[j] = 1-ihb; - } - else - { - gt[j] = idist*(1-ihb); - } - ht[j] = rhbex[j]; - nhb += ihb; - } - + gt[j] = 1-ihb; + } + else + { + gt[j] = idist*(1-ihb); + } + ht[j] = rhbex[j]; + nhb += ihb; + } - /* The autocorrelation function is normalized after summation only */ - low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0], - eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0); + /* The autocorrelation function is normalized after summation only */ + low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0], + eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0); - /* Cross correlation analysis for thermodynamics */ - for (j = nframes; (j < n2); j++) - { - ht[j] = 0; - gt[j] = 0; - } + /* Cross correlation analysis for thermodynamics */ + for (j = nframes; (j < n2); j++) + { + ht[j] = 0; + gt[j] = 0; + } - cross_corr(n2, ht, gt, dght); + cross_corr(n2, ht, gt, dght); - for (j = 0; (j < nn); j++) - { - ct[j] += rhbex[j]; - ght[j] += dght[j]; - } - } + for (j = 0; (j < nn); j++) + { + ct[j] += rhbex[j]; + ght[j] += dght[j]; } } } - fprintf(stderr, "\n"); - sfree(h); - sfree(g); - normalizeACF(ct, ght, nhb, nn); - - /* Determine tail value for statistics */ - tail = 0; - tail2 = 0; - for (j = nn/2; (j < nn); j++) - { - tail += ct[j]; - tail2 += ct[j]*ct[j]; - } - tail /= (nn - nn/2); - tail2 /= (nn - nn/2); - dtail = sqrt(tail2-tail*tail); + } + } + fprintf(stderr, "\n"); + sfree(h); + sfree(g); + normalizeACF(ct, ght, nhb, nn); - /* Check whether the ACF is long enough */ - if (dtail > tol) - { - printf("\nWARNING: Correlation function is probably not long enough\n" - "because the standard deviation in the tail of C(t) > %g\n" - "Tail value (average C(t) over second half of acf): %g +/- %g\n", - tol, tail, dtail); - } - for (j = 0; (j < nn); j++) - { - cct[j] = ct[j]; - ct[j] = (cct[j]-tail)/(1-tail); - } - /* Compute negative derivative k(t) = -dc(t)/dt */ - compute_derivative(nn, hb->time, ct, kt); - for (j = 0; (j < nn); j++) - { - kt[j] = -kt[j]; - } + /* Determine tail value for statistics */ + tail = 0; + tail2 = 0; + for (j = nn/2; (j < nn); j++) + { + tail += ct[j]; + tail2 += ct[j]*ct[j]; + } + tail /= (nn - nn/2); + tail2 /= (nn - nn/2); + dtail = sqrt(tail2-tail*tail); + /* Check whether the ACF is long enough */ + if (dtail > tol) + { + printf("\nWARNING: Correlation function is probably not long enough\n" + "because the standard deviation in the tail of C(t) > %g\n" + "Tail value (average C(t) over second half of acf): %g +/- %g\n", + tol, tail, dtail); + } + for (j = 0; (j < nn); j++) + { + cct[j] = ct[j]; + ct[j] = (cct[j]-tail)/(1-tail); + } + /* Compute negative derivative k(t) = -dc(t)/dt */ + compute_derivative(nn, hb->time, ct, kt); + for (j = 0; (j < nn); j++) + { + kt[j] = -kt[j]; + } - if (bContact) - { - fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - else - { - fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - xvgr_legend(fp, asize(legLuzar), legLuzar, oenv); + if (bContact) + { + fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); + } + else + { + fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); + } + xvgr_legend(fp, asize(legLuzar), legLuzar, oenv); - for (j = 0; (j < nn); j++) - { - fprintf(fp, "%10g %10g %10g %10g %10g\n", - hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]); - } - xvgrclose(fp); - - analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL, - fit_start, temp); - - do_view(oenv, fn, NULL); - sfree(rhbex); - sfree(ct); - sfree(gt); - sfree(ht); - sfree(ght); - sfree(dght); - sfree(cct); - sfree(kt); - /* sfree(h); */ -/* sfree(g); */ - - break; /* case AC_LUZAR */ - - default: - gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType); - } /* switch (acType) */ + + for (j = 0; (j < nn); j++) + { + fprintf(fp, "%10g %10g %10g %10g %10g\n", + hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]); + } + xvgrclose(fp); + + analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL, + fit_start, temp); + + do_view(oenv, fn, NULL); + sfree(rhbex); + sfree(ct); + sfree(gt); + sfree(ht); + sfree(ght); + sfree(dght); + sfree(cct); + sfree(kt); } static void init_hbframe(t_hbdata *hb, int nframes, real t) @@ -3286,24 +2198,6 @@ static void init_hbframe(t_hbdata *hb, int nframes, real t) { hb->nhx[nframes][i] = 0; } - /* Loop invalidated */ - if (hb->bHBmap && 0) - { - for (i = 0; (i < hb->d.nrd); i++) - { - for (j = 0; (j < hb->a.nra); j++) - { - for (m = 0; (m < hb->maxhydro); m++) - { - if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m]) - { - set_hb(hb, i, m, j, nframes, HB_NO); - } - } - } - } - } - /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/ } static FILE *open_donor_properties_file(const char *fn, @@ -3383,10 +2277,7 @@ static void dump_hbmap(t_hbdata *hb, fprintf(fp, " %4d", index[grp][i]+1); } fprintf(fp, "\n"); - /* - Added -contact support below. - - Erik Marklund, May 29, 2006 - */ + if (!bContact) { fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]); @@ -3579,7 +2470,7 @@ int gmx_hbond(int argc, char *argv[]) static int nDump = 0, nFitPoints = 100; static int nThreads = 0, nBalExp = 4; - static gmx_bool bContact = FALSE, bBallistic = FALSE, bGemFit = FALSE; + static gmx_bool bContact = FALSE; static real logAfterTime = 10, gemBallistic = 0.2; /* ps */ static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL}; @@ -3616,10 +2507,6 @@ int gmx_hbond(int argc, char *argv[]) "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" }, { "-merge", FALSE, etBOOL, {&bMerge}, "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." }, - { "-geminate", FALSE, etENUM, {gemType}, - "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."}, - { "-diff", FALSE, etREAL, {&D}, - "HIDDENDffusion 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."}, #ifdef GMX_OPENMP { "-nthreads", FALSE, etINT, {&nThreads}, "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"}, @@ -3681,13 +2568,10 @@ int gmx_hbond(int argc, char *argv[]) ivec ngrid; unsigned char *datable; output_env_t oenv; - int gemmode, NN; - PSTYPE peri = 0; - t_E E; + int NN; int ii, jj, hh, actual_nThreads; int threadNr = 0; - gmx_bool bGem, bNN, bParallel; - t_gemParams *params = NULL; + gmx_bool bParallel; gmx_bool bEdge_yjj, bEdge_xjj, bOMP; t_hbdata **p_hb = NULL; /* one per thread, then merge after the frame loop */ @@ -3708,79 +2592,6 @@ int gmx_hbond(int argc, char *argv[]) return 0; } - /* NN-loop? If so, what estimator to use ?*/ - NN = 1; - /* Outcommented for now DvdS 2010-07-13 - while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0) - NN++; - if (NN == NN_NR) - gmx_fatal(FARGS, "Invalid NN-loop type."); - */ - bNN = FALSE; - for (i = 2; bNN == FALSE && i < NN_NR; i++) - { - bNN = bNN || NN == i; - } - - if (NN > NN_NONE && bMerge) - { - bMerge = FALSE; - } - - /* geminate recombination? If so, which flavor? */ - gemmode = 1; - while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0) - { - gemmode++; - } - if (gemmode == gemNR) - { - gmx_fatal(FARGS, "Invalid recombination type."); - } - - bGem = FALSE; - for (i = 2; bGem == FALSE && i < gemNR; i++) - { - bGem = bGem || gemmode == i; - } - - if (bGem) - { - printf("Geminate recombination: %s\n", gemType[gemmode]); - if (bContact) - { - if (gemmode != gemDD) - { - printf("Turning off -contact option...\n"); - bContact = FALSE; - } - } - else - { - if (gemmode == gemDD) - { - printf("Turning on -contact option...\n"); - bContact = TRUE; - } - } - if (bMerge) - { - if (gemmode == gemAA) - { - printf("Turning off -merge option...\n"); - bMerge = FALSE; - } - } - else - { - if (gemmode != gemAA) - { - printf("Turning on -merge option...\n"); - bMerge = TRUE; - } - } - } - /* process input */ bSelected = FALSE; ccut = cos(acut*DEG2RAD); @@ -3801,8 +2612,7 @@ int gmx_hbond(int argc, char *argv[]) bHBmap = (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm) || opt2bSet("-hbn", NFILE, fnm) || - opt2bSet("-hbm", NFILE, fnm) || - bGem); + opt2bSet("-hbm", NFILE, fnm)); if (opt2bSet("-nhbdist", NFILE, fnm)) { @@ -3812,7 +2622,7 @@ int gmx_hbond(int argc, char *argv[]) xvgr_legend(fpnhb, asize(leg), leg, oenv); } - hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode); + hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact); /* get topology */ read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top); @@ -3848,7 +2658,7 @@ int gmx_hbond(int argc, char *argv[]) /* Should this be here ? */ snew(hb->d.dptr, top.atoms.nr); snew(hb->a.aptr, top.atoms.nr); - add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri); + add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact); } printf("Analyzing %d selected hydrogen bonds from '%s'\n", isize[0], grpnames[0]); @@ -3944,21 +2754,6 @@ int gmx_hbond(int argc, char *argv[]) printf("done.\n"); } -#ifdef HAVE_NN_LOOPS - if (bNN) - { - mk_hbEmap(hb, 0); - } -#endif - - if (bGem) - { - printf("Making per structure..."); - /* Generate hbond data structure */ - mk_per(hb); - printf("done.\n"); - } - /* check input */ bStop = FALSE; if (hb->d.nrd + hb->a.nra == 0) @@ -4022,11 +2817,6 @@ int gmx_hbond(int argc, char *argv[]) snew(adist, nabin+1); snew(rdist, nrbin+1); - if (bGem && !bBox) - { - gmx_fatal(FARGS, "Can't do geminate recombination without periodic box."); - } - bParallel = FALSE; #ifndef GMX_OPENMP @@ -4077,7 +2867,6 @@ int gmx_hbond(int argc, char *argv[]) p_hb[i]->bHBmap = hb->bHBmap; p_hb[i]->bDAnr = hb->bDAnr; - p_hb[i]->bGem = hb->bGem; p_hb[i]->wordlen = hb->wordlen; p_hb[i]->nframes = hb->nframes; p_hb[i]->maxhydro = hb->maxhydro; @@ -4086,11 +2875,6 @@ int gmx_hbond(int argc, char *argv[]) p_hb[i]->a = hb->a; p_hb[i]->hbmap = hb->hbmap; p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */ - p_hb[i]->per = hb->per; - -#ifdef HAVE_NN_LOOPS - p_hb[i]->hbE = hb->hbE; -#endif p_hb[i]->nrhb = 0; p_hb[i]->nrdist = 0; @@ -4102,9 +2886,9 @@ int gmx_hbond(int argc, char *argv[]) #pragma omp parallel \ firstprivate(i) \ - private(j, h, ii, jj, hh, E, \ + private(j, h, ii, jj, hh, \ xi, yi, zi, xj, yj, zj, threadNr, \ - dist, ang, peri, icell, jcell, \ + dist, ang, icell, jcell, \ grp, ogrp, ai, aj, xjj, yjj, zjj, \ xk, yk, zk, ihb, id, resdist, \ xkk, ykk, zkk, kcell, ak, k, bTric, \ @@ -4147,234 +2931,179 @@ int gmx_hbond(int argc, char *argv[]) p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */ } - if (bNN) + if (bSelected) { -#ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */ - /* Loop over all atom pairs and estimate interaction energy */ #pragma omp single { - addFramesNN(hb, nframes); - } - -#pragma omp barrier -#pragma omp for schedule(dynamic) - for (i = 0; i < hb->d.nrd; i++) - { - for (j = 0; j < hb->a.nra; j++) + /* Do not parallelize this just yet. */ + /* int ii; */ + for (ii = 0; (ii < nsel); ii++) { - for (h = 0; - h < (bContact ? 1 : hb->d.nhydro[i]); - h++) - { - if (i == hb->d.nrd || j == hb->a.nra) - { - gmx_fatal(FARGS, "out of bounds"); - } - - /* Get the real atom ids */ - ii = hb->d.don[i]; - jj = hb->a.acc[j]; - hh = hb->d.hydro[i][h]; + int dd = index[0][i]; + int aa = index[0][i+2]; + /* int */ hh = index[0][i+1]; + ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box, + hbox, &dist, &ang, bDA, &h, bContact, bMerge); - /* Estimate the energy from the geometry */ - E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d)); - /* Store the energy */ - storeHbEnergy(hb, i, j, h, E, nframes); + if (ihb) + { + /* add to index if not already there */ + /* Add a hbond */ + add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact); } } - } -#endif /* HAVE_NN_LOOPS */ - } /* if (bNN)*/ + } /* omp single */ + } /* if (bSelected) */ else { - if (bSelected) - { - -#pragma omp single - { - /* Do not parallelize this just yet. */ - /* int ii; */ - for (ii = 0; (ii < nsel); ii++) - { - int dd = index[0][i]; - int aa = index[0][i+2]; - /* int */ hh = index[0][i+1]; - ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box, - hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); - - if (ihb) - { - /* add to index if not already there */ - /* Add a hbond */ - add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri); - } - } - } /* omp single */ - } /* if (bSelected) */ - else + /* The outer grid loop will have to do for now. */ +#pragma omp for schedule(dynamic) + for (xi = 0; xi < ngrid[XX]; xi++) { - -#pragma omp single + for (yi = 0; (yi < ngrid[YY]); yi++) { - if (bGem) + for (zi = 0; (zi < ngrid[ZZ]); zi++) { - calcBoxProjection(box, hb->per->P); - } - /* loop over all gridcells (xi,yi,zi) */ - /* Removed confusing macro, DvdS 27/12/98 */ - - } - /* The outer grid loop will have to do for now. */ -#pragma omp for schedule(dynamic) - for (xi = 0; xi < ngrid[XX]; xi++) - { - for (yi = 0; (yi < ngrid[YY]); yi++) - { - for (zi = 0; (zi < ngrid[ZZ]); zi++) + /* loop over donor groups gr0 (always) and gr1 (if necessary) */ + for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++) { + icell = &(grid[zi][yi][xi].d[grp]); - /* loop over donor groups gr0 (always) and gr1 (if necessary) */ - for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++) + if (bTwo) + { + ogrp = 1-grp; + } + else { - icell = &(grid[zi][yi][xi].d[grp]); + ogrp = grp; + } - if (bTwo) - { - ogrp = 1-grp; - } - else - { - ogrp = grp; - } + /* loop over all hydrogen atoms from group (grp) + * in this gridcell (icell) + */ + for (ai = 0; (ai < icell->nr); ai++) + { + i = icell->atoms[ai]; - /* loop over all hydrogen atoms from group (grp) - * in this gridcell (icell) - */ - for (ai = 0; (ai < icell->nr); ai++) + /* loop over all adjacent gridcells (xj,yj,zj) */ + for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE); + zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE); + zjj++) { - i = icell->atoms[ai]; - - /* loop over all adjacent gridcells (xj,yj,zj) */ - for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE); - zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE); - zjj++) + zj = grid_mod(zjj, ngrid[ZZ]); + bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1); + for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj); + yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj); + yjj++) { - zj = grid_mod(zjj, ngrid[ZZ]); - bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1); - for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj); - yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj); - yjj++) + yj = grid_mod(yjj, ngrid[YY]); + bEdge_xjj = + (yj == 0) || (yj == ngrid[YY] - 1) || + (zj == 0) || (zj == ngrid[ZZ] - 1); + for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj); + xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj); + xjj++) { - yj = grid_mod(yjj, ngrid[YY]); - bEdge_xjj = - (yj == 0) || (yj == ngrid[YY] - 1) || - (zj == 0) || (zj == ngrid[ZZ] - 1); - for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj); - xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj); - xjj++) + xj = grid_mod(xjj, ngrid[XX]); + jcell = &(grid[zj][yj][xj].a[ogrp]); + /* loop over acceptor atoms from other group (ogrp) + * in this adjacent gridcell (jcell) + */ + for (aj = 0; (aj < jcell->nr); aj++) { - xj = grid_mod(xjj, ngrid[XX]); - jcell = &(grid[zj][yj][xj].a[ogrp]); - /* loop over acceptor atoms from other group (ogrp) - * in this adjacent gridcell (jcell) - */ - for (aj = 0; (aj < jcell->nr); aj++) - { - j = jcell->atoms[aj]; + j = jcell->atoms[aj]; - /* check if this once was a h-bond */ - peri = -1; - ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box, - hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); + /* check if this once was a h-bond */ + ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box, + hbox, &dist, &ang, bDA, &h, bContact, bMerge); - if (ihb) - { - /* add to index if not already there */ - /* Add a hbond */ - add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri); + if (ihb) + { + /* add to index if not already there */ + /* Add a hbond */ + add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact); - /* make angle and distance distributions */ - if (ihb == hbHB && !bContact) + /* make angle and distance distributions */ + if (ihb == hbHB && !bContact) + { + if (dist > rcut) { - if (dist > rcut) + gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist); + } + ang *= RAD2DEG; + __ADIST[(int)( ang/abin)]++; + __RDIST[(int)(dist/rbin)]++; + if (!bTwo) + { + int id, ia; + if ((id = donor_index(&hb->d, grp, i)) == NOTSET) { - gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist); + gmx_fatal(FARGS, "Invalid donor %d", i); } - ang *= RAD2DEG; - __ADIST[(int)( ang/abin)]++; - __RDIST[(int)(dist/rbin)]++; - if (!bTwo) + if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET) { - int id, ia; - if ((id = donor_index(&hb->d, grp, i)) == NOTSET) - { - gmx_fatal(FARGS, "Invalid donor %d", i); - } - if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET) - { - gmx_fatal(FARGS, "Invalid acceptor %d", j); - } - resdist = abs(top.atoms.atom[i].resind- - top.atoms.atom[j].resind); - if (resdist >= max_hx) - { - resdist = max_hx-1; - } - __HBDATA->nhx[nframes][resdist]++; + gmx_fatal(FARGS, "Invalid acceptor %d", j); } + resdist = abs(top.atoms.atom[i].resind- + top.atoms.atom[j].resind); + if (resdist >= max_hx) + { + resdist = max_hx-1; + } + __HBDATA->nhx[nframes][resdist]++; } - } - } /* for aj */ - } /* for xjj */ - } /* for yjj */ - } /* for zjj */ - } /* for ai */ - } /* for grp */ - } /* for xi,yi,zi */ - } + + } + } /* for aj */ + } /* for xjj */ + } /* for yjj */ + } /* for zjj */ + } /* for ai */ + } /* for grp */ + } /* for xi,yi,zi */ } - } /* if (bSelected) {...} else */ + } + } /* if (bSelected) {...} else */ - /* Better wait for all threads to finnish using x[] before updating it. */ - k = nframes; + /* Better wait for all threads to finnish using x[] before updating it. */ + k = nframes; #pragma omp barrier #pragma omp critical + { + /* Sum up histograms and counts from p_hb[] into hb */ + if (bOMP) { - /* Sum up histograms and counts from p_hb[] into hb */ - if (bOMP) + hb->nhb[k] += p_hb[threadNr]->nhb[k]; + hb->ndist[k] += p_hb[threadNr]->ndist[k]; + for (j = 0; j < max_hx; j++) { - hb->nhb[k] += p_hb[threadNr]->nhb[k]; - hb->ndist[k] += p_hb[threadNr]->ndist[k]; - for (j = 0; j < max_hx; j++) - { - hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j]; - } + hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j]; } } + } - /* Here are a handful of single constructs - * to share the workload a bit. The most - * important one is of course the last one, - * where there's a potential bottleneck in form - * of slow I/O. */ + /* Here are a handful of single constructs + * to share the workload a bit. The most + * important one is of course the last one, + * where there's a potential bottleneck in form + * of slow I/O. */ #pragma omp barrier #pragma omp single - { - analyse_donor_properties(donor_properties, hb, k, t); - } + { + analyse_donor_properties(donor_properties, hb, k, t); + } #pragma omp single + { + if (fpnhb) { - if (fpnhb) - { - do_nhb_dist(fpnhb, hb, t); - } + do_nhb_dist(fpnhb, hb, t); } - } /* if (bNN) {...} else + */ + } #pragma omp single { @@ -4454,11 +3183,8 @@ int gmx_hbond(int argc, char *argv[]) { max_nhb = 0.5*(hb->d.nrd*hb->a.nra); } - /* Added support for -contact below. - * - Erik Marklund, May 29-31, 2006 */ - /* Changed contact code. - * - Erik Marklund, June 29, 2006 */ - if (bHBmap && !bNN) + + if (bHBmap) { if (hb->nrhb == 0) { @@ -4471,8 +3197,6 @@ int gmx_hbond(int argc, char *argv[]) hb->nrhb, bContact ? "contacts" : "hydrogen bonds", hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding"); - /*Control the pHist.*/ - if (bMerge) { merge_hb(hb, bTwo, bContact); @@ -4572,12 +3296,10 @@ int gmx_hbond(int argc, char *argv[]) } xvgrclose(fp); } - if (!bNN) - { - printf("Average number of %s per timeframe %.3f out of %g possible\n", - bContact ? "contacts" : "hbonds", - bContact ? aver_dist : aver_nhb, max_nhb); - } + + printf("Average number of %s per timeframe %.3f out of %g possible\n", + bContact ? "contacts" : "hbonds", + bContact ? aver_dist : aver_nhb, max_nhb); /* Do Autocorrelation etc. */ if (hb->bHBmap) @@ -4596,19 +3318,9 @@ int gmx_hbond(int argc, char *argv[]) { char *gemstring = NULL; - if (bGem || bNN) - { - params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end, - gemBallistic, nBalExp); - if (params == NULL) - { - gmx_fatal(FARGS, "Could not initiate t_gemParams params."); - } - } - gemstring = gmx_strdup(gemType[hb->per->gemtype]); do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump, bMerge, bContact, fit_start, temp, r2cut > 0, oenv, - gemstring, nThreads, NN, bBallistic, bGemFit); + nThreads); } if (opt2bSet("-life", NFILE, fnm)) { @@ -4641,13 +3353,6 @@ int gmx_hbond(int argc, char *argv[]) { if (ISHB(hb->hbmap[id][ia]->history[hh])) { - /* Changed '<' into '<=' in the for-statement below. - * It fixed the previously undiscovered bug that caused - * the last occurance of an hbond/contact to not be - * set in mat.matrix. Have a look at any old -hbm-output - * and you will notice that the last column is allways empty. - * - Erik Marklund May 30, 2006 - */ for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++) { int nn0 = hb->hbmap[id][ia]->n0; @@ -4698,35 +3403,6 @@ int gmx_hbond(int argc, char *argv[]) } } - if (bGem) - { - fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper); - fprintf(stderr, "Freeing pHist for all donors...\n"); - for (i = 0; i < hb->d.nrd; i++) - { - fprintf(stderr, "\r%i", i); - if (hb->per->pHist[i] != NULL) - { - for (j = 0; j < hb->a.nra; j++) - { - clearPshift(&(hb->per->pHist[i][j])); - } - sfree(hb->per->pHist[i]); - } - } - sfree(hb->per->pHist); - sfree(hb->per->p2i); - sfree(hb->per); - fprintf(stderr, "...done.\n"); - } - -#ifdef HAVE_NN_LOOPS - if (bNN) - { - free_hbEmap(hb); - } -#endif - if (hb->bDAnr) { int i, j, nleg;