*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
+#include <cmath>
+#include <cstdio>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/trxio.h"
char title[256];
int i, j, k, Dih, ntrans;
int cur_bin, new_bin;
- real ttime, tt;
+ real ttime;
real *rot_occ[NROT];
int (*calc_bin)(real, int, real);
real dt;
prev = dih[i][j];
- mind = min(mind, dih[i][j]);
- maxd = max(maxd, dih[i][j]);
+ mind = std::min(mind, dih[i][j]);
+ maxd = std::max(maxd, dih[i][j]);
if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming */
{ /* multiplicity 3. Not so general.*/
tr_f[j]++;
{
for (i = 0; (i < nlist); i++)
{
- strncpy(name, dlist[i].name, 3);
+ std::strncpy(name, dlist[i].name, 3);
name[3] = '\0';
if (((Dih < edOmega) ) ||
((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
/* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))
{
- if ( ((strstr(name, "PHE") != NULL) && (Dih == edChi2)) ||
- ((strstr(name, "TYR") != NULL) && (Dih == edChi2)) ||
- ((strstr(name, "PTR") != NULL) && (Dih == edChi2)) ||
- ((strstr(name, "TRP") != NULL) && (Dih == edChi2)) ||
- ((strstr(name, "HIS") != NULL) && (Dih == edChi2)) ||
- ((strstr(name, "GLU") != NULL) && (Dih == edChi3)) ||
- ((strstr(name, "ASP") != NULL) && (Dih == edChi2)) ||
- ((strstr(name, "GLN") != NULL) && (Dih == edChi3)) ||
- ((strstr(name, "ASN") != NULL) && (Dih == edChi2)) ||
- ((strstr(name, "ARG") != NULL) && (Dih == edChi4)) )
+ if ( ((std::strstr(name, "PHE") != NULL) && (Dih == edChi2)) ||
+ ((std::strstr(name, "TYR") != NULL) && (Dih == edChi2)) ||
+ ((std::strstr(name, "PTR") != NULL) && (Dih == edChi2)) ||
+ ((std::strstr(name, "TRP") != NULL) && (Dih == edChi2)) ||
+ ((std::strstr(name, "HIS") != NULL) && (Dih == edChi2)) ||
+ ((std::strstr(name, "GLU") != NULL) && (Dih == edChi3)) ||
+ ((std::strstr(name, "ASP") != NULL) && (Dih == edChi2)) ||
+ ((std::strstr(name, "GLN") != NULL) && (Dih == edChi3)) ||
+ ((std::strstr(name, "ASN") != NULL) && (Dih == edChi2)) ||
+ ((std::strstr(name, "ARG") != NULL) && (Dih == edChi4)) )
{
multiplicity[j] = 2;
}
index = lookup[i][0]; /* index into dih of chi1 of res i */
if (index == -1)
{
- b = 0;
bRotZero = TRUE;
bHaveChi = FALSE;
}
{
d = invth*histo[j];
ang = j*fac-start;
- c1 = cos(ang);
- c2 = c1*c1;
+ c1 = std::cos(ang);
dc = d*c1;
- ds = d*sin(ang);
+ ds = d*std::sin(ang);
tdc += dc;
tds += ds;
for (i = 0; (i < nkkk); i++)
{
- c1 = cos(ang+kkk[i].offset);
+ c1 = std::cos(ang+kkk[i].offset);
c2 = c1*c1;
Jc = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
kkk[i].Jc += histo[j]*Jc;
for (i = 0; (i < nkkk); i++)
{
kkk[i].Jc /= th;
- kkk[i].Jcsig = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
+ kkk[i].Jcsig = std::sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
}
*S2 = tdc*tdc+tds*tds;
}
{
int i, ix, t1, t2;
rvec r_ij, r_kj;
- real costh;
+ real costh = 0.0;
for (i = ix = 0; (ix < n3); i++, ix += 3)
{
minx = maxx = data[0];
for (i = 1; (i < ndata); i++)
{
- minx = min(minx, data[i]);
- maxx = max(maxx, data[i]);
+ minx = std::min(minx, data[i]);
+ maxx = std::max(maxx, data[i]);
}
fprintf(log, "Min data: %10g Max data: %10g\n", minx, maxx);
}
}
for (i = 0; (i < ndata); i++)
{
- ind = (data[i]-minx)*dx;
+ ind = static_cast<int>((data[i]-minx)*dx);
if ((ind >= 0) && (ind < npoints))
{
histo[ind]++;
{
struct t_pbc *pbc;
t_trxstatus *status;
- int i, angind, natoms, total, teller;
+ int i, angind, total, teller;
int nangles, n_alloc;
real t, fraction, pifac, aa, angle;
real *angles[2];
#define prev (1-cur)
snew(pbc, 1);
- natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box);
+ read_first_x(oenv, &status, trj_fn, &t, &x, box);
if (bAngles)
{
for (i = 0; (i < nangles); i++)
{
real dd = angles[cur][i];
- angles[cur][i] = atan2(sin(dd), cos(dd));
+ angles[cur][i] = std::atan2(sin(dd), cos(dd));
}
}
else
}
/* Update the distribution histogram */
- angind = (int) ((angle*maxangstat)/pifac + 0.5);
+ angind = static_cast<int>((angle*maxangstat)/pifac + 0.5);
if (angind == maxangstat)
{
angind = 0;
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
*/
#include "gmxpre.h"
-#include <stdio.h>
+#include "binsearch.h"
+
+#include <cstdio>
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/real.h"
int BinarySearch (real *array, int low, int high, real key, int direction)
{
- int mid, max, min;
- max = high+2;
- min = low+1;
+ int iMid, iMax, iMin;
+ iMax = high+2;
+ iMin = low+1;
/*Iterative implementation*/
if (direction >= 0)
{
- while (max-min > 1)
+ while (iMax-iMin > 1)
{
- mid = (min+max)>>1;
- if (key < array[mid-1])
+ iMid = (iMin+iMax)>>1;
+ if (key < array[iMid-1])
{
- max = mid;
+ iMax = iMid;
}
else
{
- min = mid;
+ iMin = iMid;
}
}
- return min;
+ return iMin;
}
else if (direction < 0)
{
- while (max-min > 1)
+ while (iMax-iMin > 1)
{
- mid = (min+max)>>1;
- if (key > array[mid-1])
+ iMid = (iMin+iMax)>>1;
+ if (key > array[iMid-1])
{
- max = mid;
+ iMax = iMid;
}
else
{
- min = mid;
+ iMin = iMid;
}
}
- return min-1;
+ return iMin-1;
} /*end -ifelse direction*/
return -1;
if (direction >= 0)
{
- keyindex = startindx;
for (i = startindx; i <= stopindx; i++)
{
(*count)++;
}
else if (direction < 0)
{
- keyindex = stopindx;
for (i = stopindx; i >= startindx; i--)
{
(*count)++;
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
{
#endif
-extern void rangeArray(int *ar, int size);
+void rangeArray(int *ar, int size);
-extern void insertionSort(real *ar, int *perm, int start, int end, int direction);
+void insertionSort(real *ar, int *perm, int start, int end, int direction);
-extern int BinarySearch(real *ar, int start, int end, real key, int direction);
+int BinarySearch(real *ar, int start, int end, real key, int direction);
-extern int start_binsearch(real *array, int *perm, int low, int high,
- real key, int direction);
+int start_binsearch(real *array, int *perm, int low, int high,
+ real key, int direction);
-extern int LinearSearch(double *array, int startindx, int stopindx,
- double key, int *count, int direction);
+int LinearSearch(double *array, int startindx, int stopindx,
+ double key, int *count, int direction);
#ifdef __cplusplus
}
#include "cmat.h"
+#include <algorithm>
+
#include "gromacs/fileio/matio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/legacyheaders/macros.h"
void set_mat_entry(t_mat *m, int i, int j, real val)
{
m->mat[i][j] = m->mat[j][i] = val;
- m->maxrms = max(m->maxrms, val);
+ m->maxrms = std::max(m->maxrms, val);
if (j != i)
{
- m->minrms = min(m->minrms, val);
+ m->minrms = std::min(m->minrms, val);
}
m->sumrms += val;
- m->nn = max(m->nn, max(j+1, i+1));
+ m->nn = std::max(m->nn, std::max(j+1, i+1));
}
void done_mat(t_mat **m)
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/legacyheaders/typedefs.h"
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
typedef struct {
int i, j;
real dist;
extern t_clustid *new_clustid(int n1);
+#ifdef __cplusplus
+}
+#endif
+
#endif
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2011,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2011,2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "dens_filter.h"
-#include <math.h>
+#include <cmath>
#include "gromacs/legacyheaders/typedefs.h"
#include "gromacs/math/vec.h"
#include "gromacs/utility/smalloc.h"
-#ifdef GMX_DOUBLE
-#define EXP(x) (exp(x))
-#else
-#define EXP(x) (expf(x))
-#endif
-
-/* Modulo function assuming y>0 but with arbitrary INTEGER x */
-static int MOD(int x, int y)
-{
- if (x < 0)
- {
- x = x+y;
- }
- return (mod(x, y));
-}
-
-
gmx_bool convolution(int dataSize, real *x, int kernelSize, real* kernel)
{
int i, j, k;
gmx_bool periodic_convolution(int datasize, real *x, int kernelsize,
real *kernel)
{
- int i, j, k, nj;
+ int i, j, idx;
real *filtered;
if (!x || !kernel)
{
for (j = 0; (j < kernelsize); j++)
{
- filtered[i] += kernel[j]*x[MOD((i-j), datasize)];
+ // add datasize in case i-j is <0
+ idx = i-j + datasize;
+ filtered[i] += kernel[j]*x[idx % datasize];
}
}
for (i = 0; i < datasize; i++)
for (i = -k; i <= k; i++)
{
arg = (i*i)/(2*var);
- tot += out[j++] = EXP(-arg);
+ tot += out[j++] = std::exp(-arg);
}
for (i = 0; i < j; i++)
{
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
*/
#include "gmxpre.h"
-#include <stdlib.h>
-#include <string.h>
+#include <cstdlib>
+#include <cstring>
#include "gromacs/gmxana/gstat.h"
#include "gromacs/topology/residuetypes.h"
gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bHChi,
int maxchi, int r0, gmx_residuetype_t *rt)
{
- int ires, i, j, k, ii;
+ int i, j, ii;
t_dihatms atm, prev;
int nl = 0, nc[edMax];
char *thisres;
{
nc[i] = 0;
}
- ires = -1;
i = 0;
while (i < atoms->nr)
{
- ires = atoms->atom[i].resind;
+ int ires = atoms->atom[i].resind;
/* Initiate all atom numbers to -1 */
atm.minC = atm.H = atm.N = atm.C = atm.O = atm.minCalpha = -1;
/* maybe should allow for chis to hydrogens? */
while ((i < atoms->nr) && (atoms->atom[i].resind == ires))
{
- if ((strcmp(*(atoms->atomname[i]), "H") == 0) ||
- (strcmp(*(atoms->atomname[i]), "H1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "HN") == 0) )
+ if ((std::strcmp(*(atoms->atomname[i]), "H") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "H1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "HN") == 0) )
{
atm.H = i;
}
- else if (strcmp(*(atoms->atomname[i]), "N") == 0)
+ else if (std::strcmp(*(atoms->atomname[i]), "N") == 0)
{
atm.N = i;
}
- else if (strcmp(*(atoms->atomname[i]), "C") == 0)
+ else if (std::strcmp(*(atoms->atomname[i]), "C") == 0)
{
atm.C = i;
}
- else if ((strcmp(*(atoms->atomname[i]), "O") == 0) ||
- (strcmp(*(atoms->atomname[i]), "O1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "OC1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "OT1") == 0))
+ else if ((std::strcmp(*(atoms->atomname[i]), "O") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "O1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "OC1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "OT1") == 0))
{
atm.O = i;
}
- else if (strcmp(*(atoms->atomname[i]), "CA") == 0)
+ else if (std::strcmp(*(atoms->atomname[i]), "CA") == 0)
{
atm.Cn[1] = i;
}
- else if (strcmp(*(atoms->atomname[i]), "CB") == 0)
+ else if (std::strcmp(*(atoms->atomname[i]), "CB") == 0)
{
atm.Cn[2] = i;
}
- else if ((strcmp(*(atoms->atomname[i]), "CG") == 0) ||
- (strcmp(*(atoms->atomname[i]), "CG1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "OG") == 0) ||
- (strcmp(*(atoms->atomname[i]), "OG1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "SG") == 0))
+ else if ((std::strcmp(*(atoms->atomname[i]), "CG") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "CG1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "OG") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "OG1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "SG") == 0))
{
atm.Cn[3] = i;
}
- else if ((strcmp(*(atoms->atomname[i]), "CD") == 0) ||
- (strcmp(*(atoms->atomname[i]), "CD1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "SD") == 0) ||
- (strcmp(*(atoms->atomname[i]), "OD1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "ND1") == 0))
+ else if ((std::strcmp(*(atoms->atomname[i]), "CD") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "CD1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "SD") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "OD1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "ND1") == 0))
{
atm.Cn[4] = i;
}
/* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
- else if (bHChi && ((strcmp(*(atoms->atomname[i]), "HG") == 0) ||
- (strcmp(*(atoms->atomname[i]), "HG1") == 0)) )
+ else if (bHChi && ((std::strcmp(*(atoms->atomname[i]), "HG") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "HG1") == 0)) )
{
atm.Cn[4] = i;
}
- else if ((strcmp(*(atoms->atomname[i]), "CE") == 0) ||
- (strcmp(*(atoms->atomname[i]), "CE1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "OE1") == 0) ||
- (strcmp(*(atoms->atomname[i]), "NE") == 0))
+ else if ((std::strcmp(*(atoms->atomname[i]), "CE") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "CE1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "OE1") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "NE") == 0))
{
atm.Cn[5] = i;
}
- else if ((strcmp(*(atoms->atomname[i]), "CZ") == 0) ||
- (strcmp(*(atoms->atomname[i]), "NZ") == 0))
+ else if ((std::strcmp(*(atoms->atomname[i]), "CZ") == 0) ||
+ (std::strcmp(*(atoms->atomname[i]), "NZ") == 0))
{
atm.Cn[6] = i;
}
/* HChi flag here too */
- else if (bHChi && (strcmp(*(atoms->atomname[i]), "NH1") == 0))
+ else if (bHChi && (std::strcmp(*(atoms->atomname[i]), "NH1") == 0))
{
atm.Cn[7] = i;
}
/* added by grs - special case for aromatics, whose chis above 2 are
not real and produce rubbish output - so set back to -1 */
- if (strcmp(thisres, "PHE") == 0 ||
- strcmp(thisres, "TYR") == 0 ||
- strcmp(thisres, "PTR") == 0 ||
- strcmp(thisres, "TRP") == 0 ||
- strcmp(thisres, "HIS") == 0 ||
- strcmp(thisres, "HISA") == 0 ||
- strcmp(thisres, "HISB") == 0)
+ if (std::strcmp(thisres, "PHE") == 0 ||
+ std::strcmp(thisres, "TYR") == 0 ||
+ std::strcmp(thisres, "PTR") == 0 ||
+ std::strcmp(thisres, "TRP") == 0 ||
+ std::strcmp(thisres, "HIS") == 0 ||
+ std::strcmp(thisres, "HISA") == 0 ||
+ std::strcmp(thisres, "HISB") == 0)
{
for (ii = 5; ii <= 7; ii++)
{
fprintf(fp, "Residue\t&$\\chi_%d$\t\\\\\n", Xi+1);
for (i = 0; (i < nl); i++)
{
- nn = dl[i].ntr[Xi]/dt;
+ nn = static_cast<int>(dl[i].ntr[Xi]/dt);
if (nn == 0)
{
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
static void delete_from_block(t_block *block, int inr)
{
/* Update block data structure */
- int i, i1, j1, j, k;
+ int i, i1, j;
for (i = 0; (i < block->nr); i++)
{
void delete_atom(t_topology *top, int inr)
{
- int k;
-
if ((inr < 0) || (inr >= top->atoms.nr))
{
gmx_fatal(FARGS, "Delete_atom: inr (%d) not in %d .. %d", inr, 0,
#include "gromacs/legacyheaders/typedefs.h"
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
enum {
eWXR_NO, eWXR_YES, eWXR_NOFIT
};
int eigvalnr[],
real eigval[]);
+#ifdef __cplusplus
+}
+#endif
+
#endif
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "fitahx.h"
+#include <cmath>
+
#include "gromacs/math/do_fit.h"
#include "gromacs/math/vec.h"
#include "gromacs/utility/fatalerror.h"
for (i = 0; (i < nca); i++)
{
ai = caindex[i];
- xref[ai][XX] = rad*cos(phi0);
- xref[ai][YY] = rad*sin(phi0);
+ xref[ai][XX] = rad*std::cos(phi0);
+ xref[ai][YY] = rad*std::sin(phi0);
xref[ai][ZZ] = d*i;
/* Set the mass to select that this Calpha contributes to fitting */
{
rvec_sub(x[ai], xref[ai], dx);
rms = iprod(dx, dx);
- bb[i].rmsa += sqrt(rms);
+ bb[i].rmsa += std::sqrt(rms);
bb[i].nrms++;
trms += rms;
mass[ai] = 0.0;
}
}
- return sqrt(trms/nca);
+ return std::sqrt(trms/nca);
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/gmxana/hxprops.h"
#include "gromacs/legacyheaders/typedefs.h"
-extern real fit_ahx(int nres, t_bb bb[], int natoms, int nall, atom_id allindex[],
- rvec x[], int nca, atom_id caindex[], gmx_bool bFit);
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
+real fit_ahx(int nres, t_bb bb[], int natoms, int nall, atom_id allindex[],
+ rvec x[], int nca, atom_id caindex[], gmx_bool bFit);
+
+#ifdef __cplusplus
+}
+#endif
#endif
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "geminate.h"
-
-#include <math.h>
-#include <stdlib.h>
-
-#include "gromacs/legacyheaders/typedefs.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/gmxomp.h"
-#include "gromacs/utility/smalloc.h"
-
-static void missing_code_message()
-{
- fprintf(stderr, "You have requested code to run that is deprecated.\n");
- fprintf(stderr, "Revert to an older GROMACS version or help in porting the code.\n");
-}
-
-/* The first few sections of this file contain functions that were adopted,
- * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
- * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
- * This is also the case with the function eq10v2().
- *
- * The parts menetioned in the previous paragraph were contributed under the BSD license.
- */
-
-
-/* This first part is from complex.c which I recieved from Omer Markowitch.
- * - Erik Marklund
- *
- * ------------- from complex.c ------------- */
-
-/* Complexation of a paired number (r,i) */
-static gem_complex gem_cmplx(double r, double i)
-{
- gem_complex value;
- value.r = r;
- value.i = i;
- return value;
-}
-
-/* Complexation of a real number, x */
-static gem_complex gem_c(double x)
-{
- gem_complex value;
- value.r = x;
- value.i = 0;
- return value;
-}
-
-/* Magnitude of a complex number z */
-static double gem_cx_abs(gem_complex z)
-{
- return (sqrt(z.r*z.r+z.i*z.i));
-}
-
-/* Addition of two complex numbers z1 and z2 */
-static gem_complex gem_cxadd(gem_complex z1, gem_complex z2)
-{
- gem_complex value;
- value.r = z1.r+z2.r;
- value.i = z1.i+z2.i;
- return value;
-}
-
-/* Addition of a complex number z1 and a real number r */
-static gem_complex gem_cxradd(gem_complex z, double r)
-{
- gem_complex value;
- value.r = z.r + r;
- value.i = z.i;
- return value;
-}
-
-/* Subtraction of two complex numbers z1 and z2 */
-static gem_complex gem_cxsub(gem_complex z1, gem_complex z2)
-{
- gem_complex value;
- value.r = z1.r-z2.r;
- value.i = z1.i-z2.i;
- return value;
-}
-
-/* Multiplication of two complex numbers z1 and z2 */
-static gem_complex gem_cxmul(gem_complex z1, gem_complex z2)
-{
- gem_complex value;
- value.r = z1.r*z2.r-z1.i*z2.i;
- value.i = z1.r*z2.i+z1.i*z2.r;
- return value;
-}
-
-/* Square of a complex number z */
-static gem_complex gem_cxsq(gem_complex z)
-{
- gem_complex value;
- value.r = z.r*z.r-z.i*z.i;
- value.i = z.r*z.i*2.;
- return value;
-}
-
-/* multiplication of a complex number z and a real number r */
-static gem_complex gem_cxrmul(gem_complex z, double r)
-{
- gem_complex value;
- value.r = z.r*r;
- value.i = z.i*r;
- return value;
-}
-
-/* Division of two complex numbers z1 and z2 */
-static gem_complex gem_cxdiv(gem_complex z1, gem_complex z2)
-{
- gem_complex value;
- double num;
- num = z2.r*z2.r+z2.i*z2.i;
- if (num == 0.)
- {
- fprintf(stderr, "ERROR in gem_cxdiv function\n");
- }
- value.r = (z1.r*z2.r+z1.i*z2.i)/num; value.i = (z1.i*z2.r-z1.r*z2.i)/num;
- return value;
-}
-
-/* division of a complex z number by a real number x */
-static gem_complex gem_cxrdiv(gem_complex z, double r)
-{
- gem_complex value;
- value.r = z.r/r;
- value.i = z.i/r;
- return value;
-}
-
-/* division of a real number r by a complex number x */
-static gem_complex gem_rxcdiv(double r, gem_complex z)
-{
- gem_complex value;
- double f;
- f = r/(z.r*z.r+z.i*z.i);
- value.r = f*z.r;
- value.i = -f*z.i;
- return value;
-}
-
-/* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/
-static gem_complex gem_cxdexp(gem_complex z)
-{
- gem_complex value;
- double exp_z_r;
- exp_z_r = exp(z.r);
- value.r = exp_z_r*cos(z.i);
- value.i = exp_z_r*sin(z.i);
- return value;
-}
-
-/* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z), */
-/* where -PI < Arg(z) < PI */
-static gem_complex gem_cxlog(gem_complex z)
-{
- gem_complex value;
- double mag2;
- mag2 = z.r*z.r+z.i*z.i;
- if (mag2 < 0.)
- {
- fprintf(stderr, "ERROR in gem_cxlog func\n");
- }
- value.r = log(sqrt(mag2));
- if (z.r == 0.)
- {
- value.i = PI/2.;
- if (z.i < 0.)
- {
- value.i = -value.i;
- }
- }
- else
- {
- value.i = atan2(z.i, z.r);
- }
- return value;
-}
-
-/* Square root of a complex number z = |z| exp(I*the) -- z^(1/2) */
-/* z^(1/2)=|z|^(1/2)*[cos(the/2)+I*sin(the/2)] */
-/* where 0 < the < 2*PI */
-static gem_complex gem_cxdsqrt(gem_complex z)
-{
- gem_complex value;
- double sq;
- sq = gem_cx_abs(z);
- value.r = sqrt(fabs((sq+z.r)*0.5)); /* z'.r={|z|*[1+cos(the)]/2}^(1/2) */
- value.i = sqrt(fabs((sq-z.r)*0.5)); /* z'.i={|z|*[1-cos(the)]/2}^(1/2) */
- if (z.i < 0.)
- {
- value.r = -value.r;
- }
- return value;
-}
-
-/* Complex power of a complex number z1^z2 */
-static gem_complex gem_cxdpow(gem_complex z1, gem_complex z2)
-{
- gem_complex value;
- value = gem_cxdexp(gem_cxmul(gem_cxlog(z1), z2));
- return value;
-}
-
-/* ------------ end of complex.c ------------ */
-
-/* This next part was derived from cubic.c, also received from Omer Markovitch.
- * ------------- from cubic.c ------------- */
-
-/* Solver for a cubic equation: x^3-a*x^2+b*x-c=0 */
-static void gem_solve(gem_complex *al, gem_complex *be, gem_complex *gam,
- double a, double b, double c)
-{
- double t1, t2, two_3, temp;
- gem_complex ctemp, ct3;
-
- two_3 = pow(2., 1./3.); t1 = -a*a+3.*b; t2 = 2.*a*a*a-9.*a*b+27.*c;
- temp = 4.*t1*t1*t1+t2*t2;
-
- ctemp = gem_cmplx(temp, 0.); ctemp = gem_cxadd(gem_cmplx(t2, 0.), gem_cxdsqrt(ctemp));
- ct3 = gem_cxdpow(ctemp, gem_cmplx(1./3., 0.));
-
- ctemp = gem_rxcdiv(-two_3*t1/3., ct3);
- ctemp = gem_cxadd(ctemp, gem_cxrdiv(ct3, 3.*two_3));
-
- *gam = gem_cxadd(gem_cmplx(a/3., 0.), ctemp);
-
- ctemp = gem_cxmul(gem_cxsq(*gam), gem_cxsq(gem_cxsub(*gam, gem_cmplx(a, 0.))));
- ctemp = gem_cxadd(ctemp, gem_cxmul(gem_cmplx(-4.*c, 0.), *gam));
- ctemp = gem_cxdiv(gem_cxdsqrt(ctemp), *gam);
- *al = gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
- *be = gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
-}
-
-/* ------------ end of cubic.c ------------ */
-
-/* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch.
- * ------------- from [cr]error.c ------------- */
-
-/************************************************************/
-/* Real valued error function and related functions */
-/* x, y : real variables */
-/* erf(x) : error function */
-/* erfc(x) : complementary error function */
-/* omega(x) : exp(x*x)*erfc(x) */
-/* W(x,y) : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */
-/************************************************************/
-
-/*---------------------------------------------------------------------------*/
-/* Utilzed the series approximation of erf(x) */
-/* Relative error=|err(x)|/erf(x)<EPS */
-/* Handbook of Mathematical functions, Abramowitz, p 297 */
-/* Note: When x>=6 series sum deteriorates -> Asymptotic series used instead */
-/*---------------------------------------------------------------------------*/
-/* This gives erfc(z) correctly only upto >10-15 */
-
-static double gem_erf(double x)
-{
- double n, sum, temp, exp2, xm, x2, x4, x6, x8, x10, x12;
- temp = x;
- sum = temp;
- xm = 26.;
- x2 = x*x;
- x4 = x2*x2;
- x6 = x4*x2;
- x8 = x6*x2;
- x10 = x8*x2;
- x12 = x10*x2;
- exp2 = exp(-x2);
- if (x <= xm)
- {
- for (n = 1.; n <= 2000.; n += 1.)
- {
- temp *= 2.*x2/(2.*n+1.);
- sum += temp;
- if (fabs(temp/sum) < 1.E-16)
- {
- break;
- }
- }
-
- if (n >= 2000.)
- {
- fprintf(stderr, "In Erf calc - iteration exceeds %lg\n", n);
- }
- sum *= 2./sPI*exp2;
- }
- else
- {
- /* from the asymptotic expansion of experfc(x) */
- sum = (1. - 0.5/x2 + 0.75/x4
- - 1.875/x6 + 6.5625/x8
- - 29.53125/x10 + 162.421875/x12)
- / sPI/x;
- sum *= exp2; /* now sum is erfc(x) */
- sum = -sum+1.;
- }
- return x >= 0.0 ? sum : -sum;
-}
-
-/* Result --> Alex's code for erfc and experfc behaves better than this */
-/* complementray error function. Returns 1.-erf(x) */
-static double gem_erfc(double x)
-{
- double t, z, ans;
- z = fabs(x);
- t = 1.0/(1.0+0.5*z);
-
- ans = t * exp(-z*z - 1.26551223 +
- t*(1.00002368 +
- t*(0.37409196 +
- t*(0.09678418 +
- t*(-0.18628806 +
- t*(0.27886807 +
- t*(-1.13520398 +
- t*(1.48851587 +
- t*(-0.82215223 +
- t*0.17087277)))))))));
-
- return x >= 0.0 ? ans : 2.0-ans;
-}
-
-/* omega(x)=exp(x^2)erfc(x) */
-static double gem_omega(double x)
-{
- double xm, ans, xx, x4, x6, x8, x10, x12;
- xm = 26;
- xx = x*x;
- x4 = xx*xx;
- x6 = x4*xx;
- x8 = x6*xx;
- x10 = x8*xx;
- x12 = x10*xx;
-
- if (x <= xm)
- {
- ans = exp(xx)*gem_erfc(x);
- }
- else
- {
- /* Asymptotic expansion */
- ans = (1. - 0.5/xx + 0.75/x4 - 1.875/x6 + 6.5625/x8 - 29.53125/x10 + 162.421875/x12) / sPI/x;
- }
- return ans;
-}
-
-/*---------------------------------------------------------------------------*/
-/* Utilzed the series approximation of erf(z=x+iy) */
-/* Relative error=|err(z)|/|erf(z)|<EPS */
-/* Handbook of Mathematical functions, Abramowitz, p 299 */
-/* comega(z=x+iy)=cexp(z^2)*cerfc(z) */
-/*---------------------------------------------------------------------------*/
-static gem_complex gem_comega(gem_complex z)
-{
- gem_complex value;
- double x, y;
- double sumr, sumi, n, n2, f, temp, temp1;
- double x2, y2, cos_2xy, sin_2xy, cosh_2xy, sinh_2xy, cosh_ny, sinh_ny, exp_y2;
-
- x = z.r;
- y = z.i;
- x2 = x*x;
- y2 = y*y;
- sumr = 0.;
- sumi = 0.;
- cos_2xy = cos(2.*x*y);
- sin_2xy = sin(2.*x*y);
- cosh_2xy = cosh(2.*x*y);
- sinh_2xy = sinh(2.*x*y);
- exp_y2 = exp(-y2);
-
- for (n = 1.0, temp = 0.; n <= 2000.; n += 1.0)
- {
- n2 = n*n;
- cosh_ny = cosh(n*y);
- sinh_ny = sinh(n*y);
- f = exp(-n2/4.)/(n2+4.*x2);
- /* if(f<1.E-200) break; */
- sumr += (2.*x - 2.*x*cosh_ny*cos_2xy + n*sinh_ny*sin_2xy)*f;
- sumi += (2.*x*cosh_ny*sin_2xy + n*sinh_ny*cos_2xy)*f;
- temp1 = sqrt(sumr*sumr+sumi*sumi);
- if (fabs((temp1-temp)/temp1) < 1.E-16)
- {
- break;
- }
- temp = temp1;
- }
- if (n == 2000.)
- {
- fprintf(stderr, "iteration exceeds %lg\n", n);
- }
- sumr *= 2./PI;
- sumi *= 2./PI;
-
- if (x != 0.)
- {
- f = 1./2./PI/x;
- }
- else
- {
- f = 0.;
- }
- value.r = gem_omega(x)-(f*(1.-cos_2xy)+sumr);
- value.i = -(f*sin_2xy+sumi);
- value = gem_cxmul(value, gem_cmplx(exp_y2*cos_2xy, exp_y2*sin_2xy));
- return (value);
-}
-
-/* ------------ end of [cr]error.c ------------ */
-
-/*_ REVERSIBLE GEMINATE RECOMBINATION
- *
- * Here are the functions for reversible geminate recombination. */
-
-/* Changes the unit from square cm per s to square Ångström per ps,
- * since Omers code uses the latter units while g_mds outputs the former.
- * g_hbond expects a diffusion coefficent given in square cm per s. */
-static double sqcm_per_s_to_sqA_per_ps (real D)
-{
- fprintf(stdout, "Diffusion coefficient is %f A^2/ps\n", D*1e4);
- return (double)(D*1e4);
-}
-
-
-static double eq10v2(double theoryCt[], double time[], int manytimes,
- double ka, double kd, t_gemParams *params)
-{
- /* Finding the 3 roots */
- int i;
- double kD, D, r, a, b, c, tsqrt, sumimaginary;
- gem_complex
- alpha, beta, gamma,
- c1, c2, c3, c4,
- oma, omb, omc,
- part1, part2, part3, part4;
-
- kD = params->kD;
- D = params->D;
- r = params->sigma;
- a = (1.0 + ka/kD) * sqrt(D)/r;
- b = kd;
- c = kd * sqrt(D)/r;
-
- gem_solve(&alpha, &beta, &gamma, a, b, c);
- /* Finding the 3 roots */
-
- sumimaginary = 0;
- part1 = gem_cxmul(alpha, gem_cxmul(gem_cxadd(beta, gamma), gem_cxsub(beta, gamma))); /* 1(2+3)(2-3) */
- part2 = gem_cxmul(beta, gem_cxmul(gem_cxadd(gamma, alpha), gem_cxsub(gamma, alpha))); /* 2(3+1)(3-1) */
- part3 = gem_cxmul(gamma, gem_cxmul(gem_cxadd(alpha, beta), gem_cxsub(alpha, beta))); /* 3(1+2)(1-2) */
- part4 = gem_cxmul(gem_cxsub(gamma, alpha), gem_cxmul(gem_cxsub(alpha, beta), gem_cxsub(beta, gamma))); /* (3-1)(1-2)(2-3) */
-
-#pragma omp parallel for \
- private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \
- reduction(+:sumimaginary) \
- default(shared) \
- schedule(guided)
- for (i = 0; i < manytimes; i++)
- {
- tsqrt = sqrt(time[i]);
- oma = gem_comega(gem_cxrmul(alpha, tsqrt));
- c1 = gem_cxmul(oma, gem_cxdiv(part1, part4));
- omb = gem_comega(gem_cxrmul(beta, tsqrt));
- c2 = gem_cxmul(omb, gem_cxdiv(part2, part4));
- omc = gem_comega(gem_cxrmul(gamma, tsqrt));
- c3 = gem_cxmul(omc, gem_cxdiv(part3, part4));
- c4.r = c1.r+c2.r+c3.r;
- c4.i = c1.i+c2.i+c3.i;
- theoryCt[i] = c4.r;
- sumimaginary += c4.i * c4.i;
- }
-
- return sumimaginary;
-
-} /* eq10v2 */
-
-/* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */
-static double getLogIndex(const int i, const t_gemParams *params)
-{
- return gmx_expm1(((double)(i)) * params->logQuota);
-}
-
-extern t_gemParams *init_gemParams(const double sigma, const double D,
- const real *t, const int len, const int nFitPoints,
- const real begFit, const real endFit,
- const real ballistic, const int nBalExp)
-{
- double tDelta;
- t_gemParams *p;
-
- snew(p, 1);
-
- /* A few hardcoded things here. For now anyway. */
-/* p->ka_min = 0; */
-/* p->ka_max = 100; */
-/* p->dka = 10; */
-/* p->kd_min = 0; */
-/* p->kd_max = 2; */
-/* p->dkd = 0.2; */
- p->ka = 0;
- p->kd = 0;
-/* p->lsq = -1; */
-/* p->lifetime = 0; */
- p->sigma = sigma*10; /* Omer uses Å, not nm */
-/* p->lsq_old = 99999; */
- p->D = sqcm_per_s_to_sqA_per_ps(D);
- p->kD = 4*acos(-1.0)*sigma*p->D;
-
-
- /* Parameters used by calcsquare(). Better to calculate them
- * here than in calcsquare every time it's called. */
- p->len = len;
-/* p->logAfterTime = logAfterTime; */
- tDelta = (t[len-1]-t[0]) / len;
- if (tDelta <= 0)
- {
- gmx_fatal(FARGS, "Time between frames is non-positive!");
- }
- else
- {
- p->tDelta = tDelta;
- }
-
- p->nExpFit = nBalExp;
-/* p->nLin = logAfterTime / tDelta; */
- p->nFitPoints = nFitPoints;
- p->begFit = begFit;
- p->endFit = endFit;
- p->logQuota = (double)(log(p->len))/(p->nFitPoints-1);
-/* if (p->nLin <= 0) { */
-/* fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */
-/* sfree(p); */
-/* return NULL; */
-/* } */
-/* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */
-/* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */
-/* p->logPF = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */
-/* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */
-
-/* p->logMult = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */
- p->ballistic = ballistic;
- return p;
-}
-
-/* There was a misunderstanding regarding the fitting. From our
- * recent correspondence it appears that Omer's code require
- * the ACF data on a log-scale and does not operate on the raw data.
- * This needs to be redone in gemFunc_residual() as well as in the
- * t_gemParams structure. */
-
-static real* d2r(const double *d, const int nn);
-
-extern real fitGemRecomb(double gmx_unused *ct,
- double gmx_unused *time,
- double gmx_unused **ctFit,
- const int gmx_unused nData,
- t_gemParams gmx_unused *params)
-{
-
- int nThreads, i, iter, status, maxiter;
- real size, d2, tol, *dumpdata;
- size_t p, n;
- gemFitData *GD;
- char *dumpstr, dumpname[128];
-
- missing_code_message();
- return -1;
-
-}
-
-
-/* Removes the ballistic term from the beginning of the ACF,
- * just like in Omer's paper.
- */
-extern void takeAwayBallistic(double gmx_unused *ct, double *t, int len, real tMax, int nexp, gmx_bool gmx_unused bDerivative)
-{
-
- /* Fit with 4 exponentials and one constant term,
- * subtract the fatest exponential. */
-
- int nData, i, status, iter;
- balData *BD;
- double *guess, /* Initial guess. */
- *A, /* The fitted parameters. (A1, B1, A2, B2,... C) */
- a[2],
- ddt[2];
- gmx_bool sorted;
- size_t n;
- size_t p;
-
- nData = 0;
- do
- {
- nData++;
- }
- while (t[nData] < tMax+t[0] && nData < len);
-
- p = nexp*2+1; /* Number of parameters. */
-
- missing_code_message();
- return;
-}
-
-extern void dumpN(const real *e, const int nn, char *fn)
-{
- /* For debugging only */
- int i;
- FILE *f;
- char standardName[] = "Nt.xvg";
- if (fn == NULL)
- {
- fn = standardName;
- }
-
- f = fopen(fn, "w");
- fprintf(f,
- "@ type XY\n"
- "@ xaxis label \"Frame\"\n"
- "@ yaxis label \"N\"\n"
- "@ s0 line type 3\n");
-
- for (i = 0; i < nn; i++)
- {
- fprintf(f, "%-10i %-g\n", i, e[i]);
- }
-
- fclose(f);
-}
-
-static real* d2r(const double *d, const int nn)
-{
- real *r;
- int i;
-
- snew(r, nn);
- for (i = 0; i < nn; i++)
- {
- r[i] = (real)d[i];
- }
-
- return r;
-}
-
-static void _patchBad(double *ct, int n, double dy)
-{
- /* Just do lin. interpolation for now. */
- int i;
-
- for (i = 1; i < n; i++)
- {
- ct[i] = ct[0]+i*dy;
- }
-}
-
-static void patchBadPart(double *ct, int n)
-{
- _patchBad(ct, n, (ct[n] - ct[0])/n);
-}
-
-static void patchBadTail(double *ct, int n)
-{
- _patchBad(ct+1, n-1, ct[1]-ct[0]);
-
-}
-
-extern void fixGemACF(double *ct, int len)
-{
- int i, j, b, e;
- gmx_bool bBad;
-
- /* Let's separate two things:
- * - identification of bad parts
- * - patching of bad parts.
- */
-
- b = 0; /* Start of a bad stretch */
- e = 0; /* End of a bad stretch */
- bBad = FALSE;
-
- /* An acf of binary data must be one at t=0. */
- if (fabs(ct[0]-1.0) > 1e-6)
- {
- ct[0] = 1.0;
- fprintf(stderr, "|ct[0]-1.0| = %1.6f. Setting ct[0] to 1.0.\n", fabs(ct[0]-1.0));
- }
-
- for (i = 0; i < len; i++)
- {
-
-#ifdef HAS_ISFINITE
- if (isfinite(ct[i]))
-#elif defined(HAS__ISFINITE)
- if (_isfinite(ct[i]))
-#else
- if (1)
-#endif
- {
- if (!bBad)
- {
- /* Still on a good stretch. Proceed.*/
- continue;
- }
-
- /* Patch up preceding bad stretch. */
- if (i == (len-1))
- {
- /* It's the tail */
- if (b <= 1)
- {
- gmx_fatal(FARGS, "The ACF is mostly NaN or Inf. Aborting.");
- }
- patchBadTail(&(ct[b-2]), (len-b)+1);
- }
-
- e = i;
- patchBadPart(&(ct[b-1]), (e-b)+1);
- bBad = FALSE;
- }
- else
- {
- if (!bBad)
- {
- b = i;
-
- bBad = TRUE;
- }
- }
- }
-}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifndef _GEMINATE_H
-#define _GEMINATE_H
-
-#include <stddef.h>
-
-#include "gromacs/utility/basedefinitions.h"
-#include "gromacs/utility/real.h"
-
-enum {
- gemNULL, gemNONE, gemDD, gemAD, gemAA, gemA4, gemNR
-};
-static const char *gemType[] = {NULL, "none", "dd", "ad", "aa", "a4", NULL};
-
-/* The first few sections of this file contain functions that were adopted,
- * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
- * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
- * This is also the case with the function eq10v2() in geminate.c.
- *
- * The parts menetioned in the previous paragraph were contributed under a BSD license.
- */
-
-/* This first part is derived from complex.c which I recieved from Omer Markowitch.
- * - Erik Marklund
- *
- * ------------- from complex.c ------------- */
-
-#include <math.h>
-/* definition of PI */
-#ifndef PI
-#define PI (acos(-1.0))
-#endif
-
-/* definition of the type complex */
-typedef struct
-{
- double r, i;
-} gem_complex;
-
-
-/* ------------ end of complex.c ------------ */
-
-/* This next part was derived from cerror.c and rerror.c,
- * also received from Omer Markovitch.
- * ------------- from [cr]error.c ------------- */
-
-#ifndef sPI
-#define sPI (sqrt(PI))
-#endif
-
-/* ------------ end of [cr]error.c ------------ */
-
-/* ///////////////// REVERSIBLE GEMINATE RECOMBINATION ///////////////////
- * Here follow routines and structs for reversible geminate recombination.
- */
-
-typedef struct {
- size_t n;
- double *y;
- double tDelta;
- int nexp;
-} balData;
-
-
-typedef struct {
- /* Used in the rewritten version of Omer's gem. recomb. analysis */
- double ka, kd; /* -- || -- results[] */
- double sigma; /* -- || -- sigma */
- double D; /* The diffusion coefficient */
- double kD; /* Replaces kD in analyse_corr_gem3d() */
-
- /* The following members are for calcsquare() and takeAwayBallistic() */
- double tDelta; /* Time between frames */
- /* double logAfterTime; /\* Time after which we do the lsq calculations on a logarithmic timescale. *\/ */
- int nFitPoints; /* Number of points to take from the ACF for fitting */
- double begFit; /* Fit from this time (ps) */
- double endFit; /* Fit up to this time (ps) */
-/* double logDelta; */
-/* double logPF; */
-/* To get an equal number of points in the lin and log regime,
- * we'll use logDelta to calculate where to sample the ACF.
- * if i and j are indices in the linear and log regime, then:
- * j = Cexp(A(i+nLin)),
- * where C = (nLin**2 / len) and A = log(len/nLin) / nLin.
- * This expands to j = (nLin**2 / len) * exp((i+nLin) * log(len/nLin) / nLin).
- * We'll save part of our brains and some computing time if we pre-calculate
- * 1) logDelta = log(len/nLin) / nLin
- * 2) logPF = nLin**2 / len
- * and get j = logPF * exp((i+nLin) * logDelta). */
-
- /* I will redo this for a fit done entirely in log-log.
- * j' = j+1
- * nFitPoints' = nFitPoints-1
- *
- * j' = Cexp(Ai)
- * (j'= 1 <=> i=0)
- * => C=1
- * (j'=len <=> i=nFitPoints')
- * => A=log(len)/nFitPoints'
- * => j = exp(i*log(len)/(nFitPoints-1)) -1
- **/
-/* #define GETLOGINDEX(i,params) (params)->logPF * exp(((i)+(params)->nLin) * (params)->logDelta)
- */
- double logQuota;
- int nLin; /* Number of timepoints in the linear regime */
- int len; /* Length of time and ct arrays */
- int nExpFit; /* Number of exponentials to fit */
- real ballistic; /* Time before which the ballistic term should be fitted */
- gmx_bool bDt; /* TRUE => use time derivative at time 0
- * to find fastest component.
- * FALSE => use coefficient in exponenetial
- * to find fastest component. */
-} t_gemParams;
-
-
-typedef struct {
- size_t n; /* Number of data points (lin-log) */
- double *y; /* The datapoints */
- double *ctTheory; /* The theoretical ct which will be built by gemFunc_f. */
- double *LinLog;
- double *time;
- double ka;
- double kd;
- double tDelta; /* time difference between subsequent datapoints */
- size_t nData; /* real size of the data */
- int *logtime;
- double *doubleLogTime;
- t_gemParams *params;
-} gemFitData;
-
-extern void takeAwayBallistic(double *ct, double *t,
- int len, real tMax,
- int nexp, gmx_bool bDerivative);
-
-
-extern t_gemParams *init_gemParams(const double sigma, const double D,
- const real *t, const int len, const int nFitPoints,
- const real begFit, const real endFit,
- const real ballistic, const int nBalExp);
-
-/* Fit to geminate recombination model.
- Returns root mean square error of fit. */
-extern real fitGemRecomb(double *ct, double *time, double **ctFit,
- const int nData, t_gemParams *params);
-
-extern void dumpN(const real *e, const int nn, char *fn);
-
-/* Fix NaN that might appear in the theoretical acf. */
-extern void fixGemACF(double *ct, int len);
-
-#endif
#include "gromacs/correlationfunctions/expfit.h"
#include "gromacs/correlationfunctions/integrate.h"
#include "gromacs/fileio/xvgr.h"
-#include "gromacs/gmxana/geminate.h"
#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/gmxana/gstat.h"
#include "gromacs/legacyheaders/copyrite.h"
}
}
-static void do_ballistic(const char *balFile, int nData,
- real *t, real **val, int nSet,
- real balTime, int nBalExp,
- const output_env_t oenv)
-{
- double **ctd = NULL, *td = NULL;
- t_gemParams *GP = init_gemParams(0, 0, t, nData, 0, 0, 0, balTime, nBalExp);
- static char *leg[] = {"Ac'(t)"};
- FILE *fp;
- int i, set;
-
- if (GP->ballistic/GP->tDelta >= GP->nExpFit*2+1)
- {
- snew(ctd, nSet);
- snew(td, nData);
-
- fp = xvgropen(balFile, "Hydrogen Bond Autocorrelation", "Time (ps)", "C'(t)", oenv);
- xvgr_legend(fp, asize(leg), (const char**)leg, oenv);
-
- for (set = 0; set < nSet; set++)
- {
- snew(ctd[set], nData);
- for (i = 0; i < nData; i++)
- {
- ctd[set][i] = (double)val[set][i];
- if (set == 0)
- {
- td[i] = (double)t[i];
- }
- }
-
- takeAwayBallistic(ctd[set], td, nData, GP->ballistic, GP->nExpFit, GP->bDt);
- }
-
- for (i = 0; i < nData; i++)
- {
- fprintf(fp, " %g", t[i]);
- for (set = 0; set < nSet; set++)
- {
- fprintf(fp, " %g", ctd[set][i]);
- }
- fprintf(fp, "\n");
- }
-
-
- for (set = 0; set < nSet; set++)
- {
- sfree(ctd[set]);
- }
- sfree(ctd);
- sfree(td);
- xvgrclose(fp);
- }
- else
- {
- printf("Number 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");
- }
-}
-
-static void do_geminate(const char *gemFile, int nData,
- real *t, real **val, int nSet,
- const real D, const real rcut, const real balTime,
- const int nFitPoints, const real begFit, const real endFit,
- const output_env_t oenv)
-{
- double **ctd = NULL, **ctdGem = NULL, *td = NULL;
- t_gemParams *GP = init_gemParams(rcut, D, t, nData, nFitPoints,
- begFit, endFit, balTime, 1);
- const char *leg[] = {"Ac\\sgem\\N(t)"};
- FILE *fp;
- int i, set;
-
- snew(ctd, nSet);
- snew(ctdGem, nSet);
- snew(td, nData);
-
- fp = xvgropen(gemFile, "Hydrogen Bond Autocorrelation", "Time (ps)", "C'(t)", oenv);
- xvgr_legend(fp, asize(leg), leg, oenv);
-
- for (set = 0; set < nSet; set++)
- {
- snew(ctd[set], nData);
- snew(ctdGem[set], nData);
- for (i = 0; i < nData; i++)
- {
- ctd[set][i] = (double)val[set][i];
- if (set == 0)
- {
- td[i] = (double)t[i];
- }
- }
- fitGemRecomb(ctd[set], td, &(ctd[set]), nData, GP);
- }
-
- for (i = 0; i < nData; i++)
- {
- fprintf(fp, " %g", t[i]);
- for (set = 0; set < nSet; set++)
- {
- fprintf(fp, " %g", ctdGem[set][i]);
- }
- fprintf(fp, "\n");
- }
-
- for (set = 0; set < nSet; set++)
- {
- sfree(ctd[set]);
- sfree(ctdGem[set]);
- }
- sfree(ctd);
- sfree(ctdGem);
- sfree(td);
- xvgrclose(fp);
-}
-
static void print_fitted_function(const char *fitfile,
const char *fn_fitted,
gmx_bool bXYdy,
"The complete derivation is given in",
"B. Hess, J. Chem. Phys. 116:209-217, 2002.[PAR]",
- "Option [TT]-bal[tt] finds and subtracts the ultrafast \"ballistic\"",
- "component from a hydrogen bond autocorrelation function by the fitting",
- "of a sum of exponentials, as described in e.g.",
- "O. Markovitch, J. Chem. Phys. 129:084505, 2008. The fastest term",
- "is the one with the most negative coefficient in the exponential,",
- "or with [TT]-d[tt], the one with most negative time derivative at time 0.",
- "[TT]-nbalexp[tt] sets the number of exponentials to fit.[PAR]",
-
- "Option [TT]-gem[tt] fits bimolecular rate constants ka and kb",
- "(and optionally kD) to the hydrogen bond autocorrelation function",
- "according to the reversible geminate recombination model. Removal of",
- "the ballistic component first is strongly advised. The model is presented in",
- "O. Markovitch, J. Chem. Phys. 129:084505, 2008.[PAR]",
-
"Option [TT]-filter[tt] prints the RMS high-frequency fluctuation",
"of each set and over all sets with respect to a filtered average.",
"The filter is proportional to cos([GRK]pi[grk] t/len) where t goes from -len/2",
static gmx_bool bHaveT = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE;
static gmx_bool bEESEF = FALSE, bEENLC = FALSE, bEeFitAc = FALSE, bPower = FALSE;
static gmx_bool bIntegrate = FALSE, bRegression = FALSE, bLuzar = FALSE, bLuzarError = FALSE;
- static int nsets_in = 1, d = 1, nb_min = 4, resol = 10, nBalExp = 4, nFitPoints = 100;
- static real temp = 298.15, fit_start = 1, fit_end = 60, balTime = 0.2, diffusion = 5e-5, rcut = 0.35;
+ static int nsets_in = 1, d = 1, nb_min = 4, resol = 10, nFitPoints = 100;
+ static real temp = 298.15, fit_start = 1, fit_end = 60, diffusion = 5e-5, rcut = 0.35;
/* must correspond to enum avbar* declared at beginning of file */
static const char *avbar_opt[avbarNR+1] = {
"Subtract the average before autocorrelating" },
{ "-oneacf", FALSE, etBOOL, {&bAverCorr},
"Calculate one ACF over all sets" },
- { "-nbalexp", FALSE, etINT, {&nBalExp},
- "HIDDENNumber of exponentials to fit to the ultrafast component" },
- { "-baltime", FALSE, etREAL, {&balTime},
- "HIDDENTime up to which the ballistic component will be fitted" },
-/* { "-gemnp", FALSE, etINT, {&nFitPoints}, */
-/* "HIDDENNumber of data points taken from the ACF to use for fitting to rev. gem. recomb. model."}, */
-/* { "-rcut", FALSE, etREAL, {&rcut}, */
-/* "Cut-off for hydrogen bonds in geminate algorithms" }, */
-/* { "-gemtype", FALSE, etENUM, {gemType}, */
-/* "What type of gminate recombination to use"}, */
-/* { "-D", FALSE, etREAL, {&diffusion}, */
-/* "The self diffusion coefficient which is used for the reversible geminate recombination model."} */
};
#define NPA asize(pa)
int n, nlast, s, nset, i, j = 0;
real **val, *t, dt, tot, error;
double *av, *sig, cum1, cum2, cum3, cum4, db;
- const char *acfile, *msdfile, *ccfile, *distfile, *avfile, *eefile, *balfile, *gemfile, *fitfile;
+ const char *acfile, *msdfile, *ccfile, *distfile, *avfile, *eefile, *fitfile;
output_env_t oenv;
t_filenm fnm[] = {
{ efXVG, "-dist", "distr", ffOPTWR },
{ efXVG, "-av", "average", ffOPTWR },
{ efXVG, "-ee", "errest", ffOPTWR },
- { efXVG, "-bal", "ballisitc", ffOPTWR },
{ efXVG, "-fitted", "fitted", ffOPTWR },
-/* { efXVG, "-gem", "geminate", ffOPTWR }, */
{ efLOG, "-g", "fitlog", ffOPTWR }
};
#define NFILE asize(fnm)
distfile = opt2fn_null("-dist", NFILE, fnm);
avfile = opt2fn_null("-av", NFILE, fnm);
eefile = opt2fn_null("-ee", NFILE, fnm);
- balfile = opt2fn_null("-bal", NFILE, fnm);
-/* gemfile = opt2fn_null("-gem",NFILE,fnm); */
- /* When doing autocorrelation we don't want a fitlog for fitting
- * the function itself (not the acf) when the user did not ask for it.
- */
if (opt2parg_bSet("-fitfn", npargs, ppa) && acfile == NULL)
{
fitfile = opt2fn("-g", NFILE, fnm);
estimate_error(eefile, nb_min, resol, n, nset, av, sig, val, dt,
bEeFitAc, bEESEF, bEENLC, oenv);
}
- if (balfile)
- {
- do_ballistic(balfile, n, t, val, nset, balTime, nBalExp, oenv);
- }
-/* if (gemfile) */
-/* do_geminate(gemfile,n,t,val,nset,diffusion,rcut,balTime, */
-/* nFitPoints, fit_start, fit_end, oenv); */
if (bPower)
{
power_fit(n, nset, val, t);
#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;
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "hxprops.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/legacyheaders/macros.h"
#include "gromacs/legacyheaders/typedefs.h"
#define NPPW asize(ppw)
int i, j;
- const real Xi = 1.0;
real ell, pp2, phi, psi;
ell = 0;
fprintf(fp, "\n");
}
- return sqrt(dlt/nca);
+ return std::sqrt(dlt/nca);
}
static real rot(rvec x1, rvec x2)
real phi1, dphi, cp, sp;
real xx, yy;
- phi1 = atan2(x1[YY], x1[XX]);
- cp = cos(phi1);
- sp = sin(phi1);
+ phi1 = std::atan2(x1[YY], x1[XX]);
+ cp = std::cos(phi1);
+ sp = std::sin(phi1);
xx = cp*x2[XX]+sp*x2[YY];
yy = -sp*x2[XX]+cp*x2[YY];
- dphi = RAD2DEG*atan2(yy, xx);
+ dphi = RAD2DEG*std::atan2(yy, xx);
return dphi;
}
r0 = r1 = atom[(*index)[0]].resind;
for (i = 1; (i < gnx); i++)
{
- r0 = min(r0, atom[(*index)[i]].resind);
- r1 = max(r1, atom[(*index)[i]].resind);
+ r0 = std::min(r0, atom[(*index)[i]].resind);
+ r1 = std::max(r1, atom[(*index)[i]].resind);
}
rnr = r1-r0+1;
fprintf(stderr, "There are %d residues\n", rnr);
{
ai = (*index)[i];
ri = atom[ai].resind-r0;
- if (strcmp(*(resinfo[ri].name), "PRO") == 0)
+ if (std::strcmp(*(resinfo[ri].name), "PRO") == 0)
{
- if (strcmp(*(atomname[ai]), "CD") == 0)
+ if (std::strcmp(*(atomname[ai]), "CD") == 0)
{
bb[ri].H = ai;
}
}
for (k = 0; (k < NBB); k++)
{
- if (strcmp(bb_nm[k], *(atomname[ai])) == 0)
+ if (std::strcmp(bb_nm[k], *(atomname[ai])) == 0)
{
j++;
break;
bb[i].Cprev = bb[i-1].C;
bb[i].Nnext = bb[i+1].N;
}
- rnr = max(0, i1-i0+1);
+ rnr = std::max(0, i1-i0+1);
fprintf(stderr, "There are %d complete backbone residues (from %d to %d)\n",
rnr, bb[i0].resno, bb[i1].resno);
if (rnr == 0)
{
if (bb[i].bHelix)
{
- rms = sqrt(bb[i].pprms2);
+ rms = std::sqrt(bb[i].pprms2);
rmst += rms;
rms2 += bb[i].pprms2;
fprintf(fp, "%10g ", rms);
}
}
fprintf(fp, "\n");
- rms = sqrt(rms2/n-sqr(rmst/n));
+ rms = std::sqrt(rms2/n-sqr(rmst/n));
return rms;
}
bb[i].pprms2 = sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX);
bb[i].jcaha +=
- 1.4*sin((bb[i].psi+138.0)*DEG2RAD) -
- 4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
- 2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0));
+ 1.4*std::sin((bb[i].psi+138.0)*DEG2RAD) -
+ 4.1*std::cos(2.0*DEG2RAD*(bb[i].psi+138.0)) +
+ 2.0*std::cos(2.0*DEG2RAD*(bb[i].phi+30.0));
}
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/legacyheaders/typedefs.h"
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+
#define PHI_AHX (-55.0)
#define PSI_AHX (-45.0)
/* Canonical values of the helix phi/psi angles */
extern void pr_bb(FILE *fp, int nres, t_bb bb[]);
+#ifdef __cplusplus
+}
+#endif
+
#endif
#include "nrama.h"
-#include <math.h>
-#include <stdlib.h>
+#include <cstdlib>
#include "gromacs/listed-forces/bonded.h"
#include "gromacs/pbcutil/rmpbc.h"
const char *topfile, t_xrama *xr, int mult)
{
t_topology *top;
- int ePBC;
real t;
top = read_top(topfile, &xr->ePBC);
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "config.h"
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/fileio/strdb.h"
#include "gromacs/math/vec.h"
int nthreads;
gmx_rng_t *trng = NULL;
#endif
- gmx_int64_t mc = 0, max;
+ gmx_int64_t mc = 0, mc_max;
gmx_rng_t rng = NULL;
/* allocate memory for pr */
rmax = norm(dist);
- pr->grn = (int)floor(rmax/pr->binwidth)+1;
- rmax = pr->grn*pr->binwidth;
+ pr->grn = static_cast<int>(std::floor(rmax/pr->binwidth)+1);
snew(pr->gr, pr->grn);
/* Special case for setting automaticaly number of mc iterations to 1% of total number of direct iterations */
if (mcover == -1)
{
- max = (gmx_int64_t)floor(0.5*0.01*isize*(isize-1));
+ mc_max = static_cast<gmx_int64_t>(std::floor(0.5*0.01*isize*(isize-1)));
}
else
{
- max = (gmx_int64_t)floor(0.5*mcover*isize*(isize-1));
+ mc_max = static_cast<gmx_int64_t>(std::floor(0.5*mcover*isize*(isize-1)));
}
rng = gmx_rng_init(seed);
#ifdef GMX_OPENMP
tid = gmx_omp_get_thread_num();
/* now starting parallel threads */
#pragma omp for
- for (mc = 0; mc < max; mc++)
+ for (mc = 0; mc < mc_max; mc++)
{
- i = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
- j = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
+ i = static_cast<int>(std::floor(gmx_rng_uniform_real(trng[tid])*isize));
+ j = static_cast<int>(std::floor(gmx_rng_uniform_real(trng[tid])*isize));
if (i != j)
{
- tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
+ tgr[tid][static_cast<int>(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]];
}
}
}
sfree(tgr);
sfree(trng);
#else
- for (mc = 0; mc < max; mc++)
+ for (mc = 0; mc < mc_max; mc++)
{
- i = (int)floor(gmx_rng_uniform_real(rng)*isize);
- j = (int)floor(gmx_rng_uniform_real(rng)*isize);
+ i = static_cast<int>(std::floor(gmx_rng_uniform_real(rng)*isize));
+ j = static_cast<int>(std::floor(gmx_rng_uniform_real(rng)*isize));
if (i != j)
{
- pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
+ pr->gr[static_cast<int>(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]];
}
}
#endif
{
for (j = 0; j < i; j++)
{
- tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
+ tgr[tid][static_cast<int>(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]];
}
}
}
{
for (j = 0; j < i; j++)
{
- pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
+ pr->gr[static_cast<int>(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]];
}
}
#endif
int i, j;
/* init data */
snew(sq, 1);
- sq->qn = (int)floor((end_q-start_q)/q_step);
+ sq->qn = static_cast<int>(std::floor((end_q-start_q)/q_step));
snew(sq->q, sq->qn);
snew(sq->s, sq->qn);
for (i = 0; i < sq->qn; i++)
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdlib.h>
+#include <cmath>
+#include <cstdlib>
+
+#include <algorithm>
#include "gromacs/fileio/matio.h"
#include "gromacs/gmxana/gstat.h"
newdata[i][j] = sd->data[i/nfac][j/nfac];
else*/
newdata[i][j] = interpolate(phi, psi, sd);
- lo = min(lo, newdata[i][j]);
- hi = max(hi, newdata[i][j]);
+ lo = std::min(lo, newdata[i][j]);
+ hi = std::max(hi, newdata[i][j]);
}
}
sprintf(buf, "%s.xpm", fn);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2014, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "princ.h"
+#include <cmath>
+
#include "gromacs/legacyheaders/txtdump.h"
#include "gromacs/legacyheaders/typedefs.h"
#include "gromacs/linearalgebra/nrjac.h"
z = inten[m][3];
n = x*x+y*y+z*z;
fprintf(stderr, "%8s %8.3f %8.3f %8.3f, norm:%8.3f, d:%8.3f, e:%8.3f\n",
- s, x, y, z, sqrt(n), d[m], e[m]);
+ s, x, y, z, std::sqrt(n), d[m], e[m]);
}
fprintf(stderr, "\n");
}
#include "sfactor.h"
+#include <cmath>
+
+#include <algorithm>
+
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/strdb.h"
#include "gromacs/fileio/trxio.h"
k_factor[YY] = 2 * M_PI / box[YY][YY];
k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];
- maxkx = (int) (end_q / k_factor[XX] + 0.5);
- maxky = (int) (end_q / k_factor[YY] + 0.5);
- maxkz = (int) (end_q / k_factor[ZZ] + 0.5);
+ maxkx = static_cast<int>(end_q / k_factor[XX] + 0.5);
+ maxky = static_cast<int>(end_q / k_factor[YY] + 0.5);
+ maxkz = static_cast<int>(end_q / k_factor[ZZ] + 0.5);
snew (counter, sf->n_angles);
fprintf(stderr, "\n");
for (i = 0; i < maxkx; i++)
{
- fprintf (stderr, "\rdone %3.1f%% ", (double)(100.0*(i+1))/maxkx);
+ fprintf (stderr, "\rdone %3.1f%% ", (100.0*(i+1))/maxkx);
kx = i * k_factor[XX];
for (j = 0; j < maxky; j++)
{
if (i != 0 || j != 0 || k != 0)
{
kz = k * k_factor[ZZ];
- krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz));
+ krr = std::sqrt (sqr (kx) + sqr (ky) + sqr (kz));
if (krr >= start_q && krr <= end_q)
{
- kr = (int) (krr/sf->ref_k + 0.5);
+ kr = static_cast<int>(krr/sf->ref_k + 0.5);
if (kr < sf->n_angles)
{
counter[kr]++; /* will be used for the copmutation
kdotx = kx * redt[p].x[XX] +
ky * redt[p].x[YY] + kz * redt[p].x[ZZ];
- tmpSF[i][j][k].re += cos (kdotx) * asf;
- tmpSF[i][j][k].im += sin (kdotx) * asf;
+ tmpSF[i][j][k].re += std::cos(kdotx) * asf;
+ tmpSF[i][j][k].im += std::sin(kdotx) * asf;
}
}
}
{
ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++)
{
- kz = k * k_factor[ZZ]; krr = sqrt (sqr (kx) + sqr (ky)
- + sqr (kz)); if (krr >= start_q && krr <= end_q)
+ kz = k * k_factor[ZZ]; krr = std::sqrt (sqr (kx) + sqr (ky)
+ + sqr (kz)); if (krr >= start_q && krr <= end_q)
{
- kr = (int) (krr / sf->ref_k + 0.5);
+ kr = static_cast<int>(krr / sf->ref_k + 0.5);
if (kr < sf->n_angles && counter[kr] != 0)
{
sf->F[group][kr] +=
structure_factor *sf;
rvec *xtop;
real **sf_table;
- int nsftable;
matrix box;
- double r_tmp;
+ real r_tmp;
gmx_structurefactors_t *gmx_sf;
real *a, *b, c;
- int success;
snew(a, 4);
snew(b, 4);
gmx_sf = gmx_structurefactors_init(fnDAT);
- success = gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c);
+ gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c);
snew (sf, 1);
sf->energy = energy;
snew (red, ng);
snew (index_atp, ng);
- r_tmp = max (box[XX][XX], box[YY][YY]);
- r_tmp = (double) max (box[ZZ][ZZ], r_tmp);
+ r_tmp = std::max(box[XX][XX], box[YY][YY]);
+ r_tmp = std::max(box[ZZ][ZZ], r_tmp);
sf->ref_k = (2.0 * M_PI) / (r_tmp);
/* ref_k will be the reference momentum unit */
- sf->n_angles = (int) (end_q / sf->ref_k + 0.5);
+ sf->n_angles = static_cast<int>(end_q / sf->ref_k + 0.5);
snew (sf->F, ng);
for (i = 0; i < ng; i++)
* sin(theta) = q/(2k) := A -> sin^2(theta) = 4*A^2 (1-A^2) ->
* -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
*/
- A = (double) (i * sf->ref_k) / (2.0 * sf->momentum);
+ A = static_cast<double>(i * sf->ref_k) / (2.0 * sf->momentum);
polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));
sf->F[g][i] *= polarization_factor;
}
* vectors. See g_sq.h for a short description of CM fit.
*/
{
- int i, success;
+ int i;
double tmp = 0.0, k2;
real *a, *b;
real c;
else
{
k2 = (sqr (sin_theta) / sqr (10.0 * lambda));
- success = gmx_structurefactors_get_sf(gsf, type, a, b, &c);
+ gmx_structurefactors_get_sf(gsf, type, a, b, &c);
tmp = c;
for (i = 0; (i < 4); i++)
{
snew (sf_table[i], n_angles);
for (j = 0; j < n_angles; j++)
{
- q = ((double) j * ref_k);
+ q = static_cast<double>(j * ref_k);
/* theta is half the angle between incoming
and scattered wavevectors. */
sin_theta = q / (2.0 * momentum);
/* \hbar \omega \lambda = hc = 1239.842 eV * nm */
- sf->momentum = ((double) (2. * 1000.0 * M_PI * sf->energy) / hc);
+ sf->momentum = (static_cast<double>(2. * 1000.0 * M_PI * sf->energy) / hc);
sf->lambda = hc / (1000.0 * sf->energy);
fprintf (stderr, "\nwavelenght = %f nm\n", sf->lambda);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff