-/* 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;
-}
-
-