#include "config.h"
-#include <math.h>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/gmxana/gstat.h"
#include "gromacs/legacyheaders/copyrite.h"
#include "gromacs/legacyheaders/macros.h"
#include "gromacs/legacyheaders/txtdump.h"
static void add_frames(t_hbdata *hb, int nframes)
{
- int i, j, k, l;
-
if (nframes >= hb->max_frames)
{
hb->max_frames += 4096;
{
int i, j, n;
t_hbond *hb = hbd->hbmap[id][ia];
- int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]);
+ int maxhydro = std::min(hbd->maxhydro, hbd->d.nhydro[id]);
int wlen = hbd->wordlen;
int delta = 32*wlen;
t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
unsigned char *datable)
{
- int i, j, nra, n;
+ int i, j;
t_functype func_type;
t_ilist *interaction;
atom_id nr1, nr2, nr3;
- gmx_bool stop;
if (!ddd->dptr)
{
{
for (i = 0; i < DIM; i++)
{
- ngrid[i] = (box[i][i]/(1.2*rcut));
+ ngrid[i] = static_cast<int>(box[i][i]/(1.2*rcut));
}
}
int i, m, gr, xi, yi, zi, nr;
atom_id *ad;
ivec grididx;
- rvec invdelta, dshell, xtemp = {0, 0, 0};
+ rvec invdelta, dshell;
t_ncell *newgrid;
gmx_bool bDoRshell, bInShell, bAcc;
real rshell2 = 0;
for (m = DIM-1; m >= 0; m--)
{ /* determine grid index of atom */
- grididx[m] = x[ad[i]][m]*invdelta[m];
+ grididx[m] = static_cast<int>(x[ad[i]][m]*invdelta[m]);
grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
}
}
real *d_ha, real *ang, gmx_bool bDA, int *hhh,
gmx_bool bContact, gmx_bool bMerge)
{
- int h, hh, id, ja, ihb;
- rvec r_da, r_ha, r_dh, r = {0, 0, 0};
- ivec ri;
+ int h, hh, id;
+ rvec r_da, r_ha, r_dh;
real rc2, r2c2, rda2, rha2, ca;
gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */
gmx_bool daSwap = FALSE;
}
if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
- ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET))
+ (acceptor_index(&hb->a, grpa, a) == NOTSET))
{
return hbNo;
}
if (ca >= ccut)
{
*hhh = hh;
- *d_ha = sqrt(bDA ? rda2 : rha2);
- *ang = acos(ca);
+ *d_ha = std::sqrt(bDA ? rda2 : rha2);
+ *ang = std::acos(ca);
return hbHB;
}
}
/* Decide where to start from when merging */
n00 = hb0->n0;
n01 = hb1->n0;
- nn0 = min(n00, n01);
- nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
+ nn0 = std::min(n00, n01);
+ nnframes = std::max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0;
/* Initiate tmp arrays */
for (m = 0; (m < ntmp); m++)
{
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;
+ int i, inrnew, indnew, j, ii, jj, id, ia, ntmp;
unsigned int *htmp, *gtmp;
t_hbond *hb0, *hb1;
static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
{
int i, j, k, n_bound[MAXHH], nbtot;
- h_id nhb;
-
/* Set array to 0 */
for (k = 0; (k < MAXHH); k++)
}
else
{
- return kbt*log(kbt*tau/PLANCK);
+ return kbt*std::log(kbt*tau/PLANCK);
}
}
tl.ndelta = NK;
for (j = 0; (j < NK); j++)
{
- /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
kkk += tl.kkk[0];
kkp += tl.kkk[1];
kk2 += sqr(tl.kkk[0]);
kp2 += sqr(tl.kkk[1]);
tl.n0++;
}
- *sigma_k = sqrt(kk2/NK - sqr(kkk/NK));
- *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK));
+ *sigma_k = std::sqrt(kk2/NK - sqr(kkk/NK));
+ *sigma_kp = std::sqrt(kp2/NK - sqr(kkp/NK));
return chi2;
}
{
int i0, i;
real k = 1, kp = 1, kow = 1;
- real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg;
+ real Q = 0, chi2, tau_hb, dtau, tau_rlx, e_1, sigma_k, sigma_kp, ddg;
double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL);
{
chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]);
}
- chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
- sigma_kt, &k, &kp,
- &sigma_k, &sigma_kp, fit_start);
+ compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
+ sigma_kt, &k, &kp,
+ &sigma_k, &sigma_kp, fit_start);
Q = 0; /* quality_of_fit(chi2, 2);*/
ddg = BOLTZ*temp*sigma_k/k;
printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau);
printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
- e_1 = exp(-1);
+ e_1 = std::exp(-1.0);
for (i = 0; (i < n-2); i++)
{
if ((ct[i] > e_1) && (ct[i+1] <= e_1))
ct_fac = 1.0/ct[0];
if (nhb != 0)
{
- gt_fac = 1.0/(real)nhb;
+ gt_fac = 1.0/nhb;
}
printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
int nThreads)
{
FILE *fp;
- int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets;
- const char *legNN[] = {
- "Ac(t)",
- "Ac'(t)"
- };
- static char **legGem;
+ int i, j, k, m, ihb, idist, n2, nn;
const char *legLuzar[] = {
"Ac\\sfin sys\\v{}\\z{}(t)",
};
gmx_bool bNorm = FALSE, bOMP = FALSE;
double nhb = 0;
- int nhbi = 0;
real *rhbex = NULL, *ht, *gt, *ght, *dght, *kt;
- real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct;
+ real *ct, tail, tail2, dtail, *cct;
const real tol = 1e-3;
- int nframes = hb->nframes, nf;
+ int nframes = hb->nframes;
unsigned int **h = NULL, **g = NULL;
- int nh, nhbonds, nhydro, ngh;
+ int nh, nhbonds, nhydro;
t_hbond *hbh;
- int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN;
- gmx_bool c;
int acType;
- double *ctdouble, *timedouble, *fittedct;
- double fittolerance = 0.1;
- int *dondata = NULL, thisThread;
+ int *dondata = NULL;
enum {
AC_NONE, AC_NN, AC_GEM, AC_LUZAR
/* Total number of hbonds analyzed here */
nhbonds = 0;
- ngh = 0;
- anhb = 0;
if (acType != AC_LUZAR && bOMP)
{
- nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
+ nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
gmx_omp_set_num_threads(nThreads);
snew(dondata, nThreads);
}
}
- nf = hbh->nframes;
+ int nf = hbh->nframes;
for (nh = 0; (nh < nhydro); nh++)
{
int nrint = bContact ? hb->nrdist : hb->nrhb;
fprintf(stderr, "\n");
sfree(h);
sfree(g);
- normalizeACF(ct, ght, nhb, nn);
+ normalizeACF(ct, ght, static_cast<int>(nhb), nn);
/* Determine tail value for statistics */
tail = 0;
}
tail /= (nn - nn/2);
tail2 /= (nn - nn/2);
- dtail = sqrt(tail2-tail*tail);
+ dtail = std::sqrt(tail2-tail*tail);
/* Check whether the ACF is long enough */
if (dtail > tol)
static void init_hbframe(t_hbdata *hb, int nframes, real t)
{
- int i, j, m;
+ int i;
hb->time[nframes] = t;
hb->nhb[nframes] = 0;
* It mimics add_frames() and init_frame() to some extent. */
static void sync_hbdata(t_hbdata *p_hb, int nframes)
{
- int i;
if (nframes >= p_hb->max_frames)
{
p_hb->max_frames += 4096;
{
srenew(p_hb->danr, p_hb->max_frames);
}
- memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
- memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
+ std::memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
+ std::memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
p_hb->nhb[nframes] = 0;
p_hb->ndist[nframes] = 0;
}
p_hb->nframes = nframes;
-/* for (i=0;) */
-/* { */
-/* p_hb->nhx[nframes][i] */
-/* } */
- memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
-
- /* hb->per will remain constant througout the frame loop,
- * even though the data its members point to will change,
- * hence no need for re-syncing. */
+
+ std::memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
}
int gmx_hbond(int argc, char *argv[])
};
static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
- static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, D = -1;
+ static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15;
static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
- static int nDump = 0, nFitPoints = 100;
- static int nThreads = 0, nBalExp = 4;
+ static int nDump = 0;
+ static int nThreads = 0;
static gmx_bool bContact = FALSE;
- static real logAfterTime = 10, gemBallistic = 0.2; /* ps */
- static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL};
/* options */
t_pargs pa [] = {
matrix box;
real t, ccut, dist = 0.0, ang = 0.0;
double max_nhb, aver_nhb, aver_dist;
- int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel;
+ int h = 0, i = 0, j, k = 0, ogrp, nsel;
int xi, yi, zi, ai;
int xj, yj, zj, aj, xjj, yjj, zjj;
- int xk, yk, zk, ak, xkk, ykk, zkk;
- gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric;
- int *adist, *rdist, *aptr, *rprt;
- int grp, nabin, nrbin, bin, resdist, ihb;
+ gmx_bool bSelected, bHBmap, bStop, bTwo, bBox, bTric;
+ int *adist, *rdist;
+ int grp, nabin, nrbin, resdist, ihb;
char **leg;
- t_hbdata *hb, *hbptr;
- FILE *fp, *fpins = NULL, *fpnhb = NULL, *donor_properties = NULL;
+ t_hbdata *hb;
+ FILE *fp, *fpnhb = NULL, *donor_properties = NULL;
t_gridcell ***grid;
- t_ncell *icell, *jcell, *kcell;
+ t_ncell *icell, *jcell;
ivec ngrid;
unsigned char *datable;
output_env_t oenv;
- int NN;
- int ii, jj, hh, actual_nThreads;
+ int ii, hh, actual_nThreads;
int threadNr = 0;
gmx_bool bParallel;
gmx_bool bEdge_yjj, bEdge_xjj, bOMP;
/* process input */
bSelected = FALSE;
- ccut = cos(acut*DEG2RAD);
+ ccut = std::cos(acut*DEG2RAD);
if (bContact)
{
bBox = ir.ePBC != epbcNONE;
grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
- nabin = acut/abin;
- nrbin = rcut/rbin;
+ nabin = static_cast<int>(acut/abin);
+ nrbin = static_cast<int>(rcut/rbin);
snew(adist, nabin+1);
snew(rdist, nrbin+1);
- bParallel = FALSE;
-
#ifndef GMX_OPENMP
#define __ADIST adist
#define __RDIST rdist
if (bParallel)
{
- actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
+ actual_nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
gmx_omp_set_num_threads(actual_nThreads);
printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
#pragma omp parallel \
firstprivate(i) \
- private(j, h, ii, jj, hh, \
+ private(j, h, ii, hh, \
xi, yi, zi, xj, yj, zj, threadNr, \
dist, ang, icell, jcell, \
grp, ogrp, ai, aj, xjj, yjj, zjj, \
- xk, yk, zk, ihb, id, resdist, \
- xkk, ykk, zkk, kcell, ak, k, bTric, \
+ ihb, resdist, \
+ k, bTric, \
bEdge_xjj, bEdge_yjj) \
default(shared)
{ /* Start of parallel region */
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)]++;
+ __ADIST[static_cast<int>( ang/abin)]++;
+ __RDIST[static_cast<int>(dist/rbin)]++;
if (!bTwo)
{
- int id, ia;
- if ((id = donor_index(&hb->d, grp, i)) == NOTSET)
+ if (donor_index(&hb->d, grp, i) == NOTSET)
{
gmx_fatal(FARGS, "Invalid donor %d", i);
}
- if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET)
+ if (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);
+ resdist = std::abs(top.atoms.atom[i].resind-top.atoms.atom[j].resind);
if (resdist >= max_hx)
{
resdist = max_hx-1;
/* Print HB distance distribution */
if (opt2bSet("-dist", NFILE, fnm))
{
- long sum;
+ int sum;
sum = 0;
for (i = 0; i < nrbin; i++)
"Hydrogen - Acceptor Distance (nm)", "", oenv);
for (i = 0; i < nrbin; i++)
{
- fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
+ fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*sum));
}
xvgrclose(fp);
}
"Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
for (i = 0; i < nabin; i++)
{
- fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
+ fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*sum));
}
xvgrclose(fp);
}
Added support for -contact in ac and hbm calculations below.
- Erik Marklund, May 29, 2006
*/
- ivec itmp;
- rvec rtmp;
if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm))
{
please_cite(stdout, "Spoel2006b");
}
if (opt2bSet("-ac", NFILE, fnm))
{
- char *gemstring = NULL;
-
do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
bMerge, bContact, fit_start, temp, r2cut > 0, oenv,
nThreads);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff