#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"
#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
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" };
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;
/* 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;
hb->wordlen = 8*sizeof(unsigned int);
hb->bHBmap = bHBmap;
hb->bDAnr = bDAnr;
- hb->bGem = bGem;
if (oneHB)
{
hb->maxhydro = 1;
{
hb->maxhydro = MAXHYDRO;
}
- snew(hb->per, 1);
- hb->per->gemtype = bGem ? gemmode : 0;
-
return 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;
_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])
{
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);
- }
-
- }
}
}
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;
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);
}
}
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;
}
}
}
{
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--)
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];
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};
{ /* 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);
}
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)
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);
}
}
-/* 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;
{
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 '<='
{
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++)
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)
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 */
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;
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);
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--;
}
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;
hb->nrdist = indnew;
sfree(gtmp);
sfree(htmp);
- sfree(ptmp);
}
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)
{
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);
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;
}
}
-/* 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;
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;
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 */
}
- /* 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<poff[m]");
- }
- rHbExGem[m][j-poff[m]] += 1;
- }
- }
- }
-
- /* Now, build ac. */
- for (m = 0; m < np; m++)
- {
- if (rHbExGem[m][0] > 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)
{
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,
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]);
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};
"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)"},
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 */
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);
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))
{
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);
/* 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]);
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)
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
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;
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;
#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, \
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
{
{
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)
{
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);
}
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)
{
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))
{
{
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;
}
}
- 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;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff