*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
/*! \brief Shortcut macro to select modes. */
#define MODE(x) ((mode & (x)) == (x))
-typedef struct {
+typedef struct
+{
unsigned long mode;
int nrestart, nout, P, fitfn;
gmx_bool bFour, bNormalize;
} t_acf;
/*! \brief Global variable set true if initialization routines are called. */
-static gmx_bool bACFinit = FALSE;
+static gmx_bool bACFinit = FALSE;
/*! \brief Data structure for storing command line variables. */
-static t_acf acf;
+static t_acf acf;
-enum {
- enNorm, enCos, enSin
+enum
+{
+ enNorm,
+ enCos,
+ enSin
};
/*! \brief Routine to compute ACF using FFT. */
-static void low_do_four_core(int nframes, real c1[], real cfour[],
- int nCos)
+static void low_do_four_core(int nframes, real c1[], real cfour[], int nCos)
{
- int i = 0;
- std::vector<std::vector<real> > data;
+ int i = 0;
+ std::vector<std::vector<real>> data;
data.resize(1);
data[0].resize(nframes, 0);
switch (nCos)
data[0][i] = sin(c1[i]);
}
break;
- default:
- gmx_fatal(FARGS, "nCos = %d, %s %d", nCos, __FILE__, __LINE__);
+ default: gmx_fatal(FARGS, "nCos = %d, %s %d", nCos, __FILE__, __LINE__);
}
many_auto_correl(&data);
}
/*! \brief Routine to comput ACF without FFT. */
-static void do_ac_core(int nframes, int nout,
- real corr[], real c1[], int nrestart,
- unsigned long mode)
+static void do_ac_core(int nframes, int nout, real corr[], real c1[], int nrestart, unsigned long mode)
{
- int j, k, j3, jk3, m, n;
- real ccc, cth;
- rvec xj, xk;
+ int j, k, j3, jk3, m, n;
+ real ccc, cth;
+ rvec xj, xk;
if (nrestart < 1)
{
}
if (debug)
{
- fprintf(debug,
- "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n",
- nframes, nout, nrestart, mode);
+ fprintf(debug, "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n", nframes,
+ nout, nrestart, mode);
}
for (j = 0; (j < nout); j++)
/* Loop over starting points. */
for (j = 0; (j < nframes); j += nrestart)
{
- j3 = DIM*j;
+ j3 = DIM * j;
/* Loop over the correlation length for this starting point */
- for (k = 0; (k < nout) && (j+k < nframes); k++)
+ for (k = 0; (k < nout) && (j + k < nframes); k++)
{
- jk3 = DIM*(j+k);
+ jk3 = DIM * (j + k);
/* Switch over possible ACF types.
* It might be more efficient to put the loops inside the switch,
*/
if (MODE(eacNormal))
{
- corr[k] += c1[j]*c1[j+k];
+ corr[k] += c1[j] * c1[j + k];
}
else if (MODE(eacCos))
{
/* Compute the cos (phi(t)-phi(t+dt)) */
- corr[k] += std::cos(c1[j]-c1[j+k]);
+ corr[k] += std::cos(c1[j] - c1[j + k]);
}
else if (MODE(eacIden))
{
/* Check equality (phi(t)==phi(t+dt)) */
- corr[k] += (c1[j] == c1[j+k]) ? 1 : 0;
+ corr[k] += (c1[j] == c1[j + k]) ? 1 : 0;
}
else if (MODE(eacP1) || MODE(eacP2) || MODE(eacP3))
{
for (m = 0; (m < DIM); m++)
{
- xj[m] = c1[j3+m];
- xk[m] = c1[jk3+m];
+ xj[m] = c1[j3 + m];
+ xk[m] = c1[jk3 + m];
}
cth = cos_angle(xj, xk);
- if (cth-1.0 > 1.0e-15)
+ if (cth - 1.0 > 1.0e-15)
{
- printf("j: %d, k: %d, xj:(%g,%g,%g), xk:(%g,%g,%g)\n",
- j, k, xj[XX], xj[YY], xj[ZZ], xk[XX], xk[YY], xk[ZZ]);
+ printf("j: %d, k: %d, xj:(%g,%g,%g), xk:(%g,%g,%g)\n", j, k, xj[XX], xj[YY],
+ xj[ZZ], xk[XX], xk[YY], xk[ZZ]);
}
mmm = 1;
if (MODE(eacP2))
rvec xj, xk, rr;
for (m = 0; (m < DIM); m++)
{
- xj[m] = c1[j3+m];
- xk[m] = c1[jk3+m];
+ xj[m] = c1[j3 + m];
+ xk[m] = c1[jk3 + m];
}
cprod(xj, xk, rr);
{
for (m = 0; (m < DIM); m++)
{
- xj[m] = c1[j3+m];
- xk[m] = c1[jk3+m];
+ xj[m] = c1[j3 + m];
+ xk[m] = c1[jk3 + m];
}
ccc = iprod(xj, xk);
/* Correct for the number of points and copy results to the data array */
for (j = 0; (j < nout); j++)
{
- n = (nframes-j+(nrestart-1))/nrestart;
- c1[j] = corr[j]/n;
+ n = (nframes - j + (nrestart - 1)) / nrestart;
+ c1[j] = corr[j] / n;
}
}
}
else
{
- c0 = 1.0/corr[0];
+ c0 = 1.0 / corr[0];
}
for (j = 0; (j < nout); j++)
{
}
/*! \brief Routine that averages ACFs. */
-static void average_acf(gmx_bool bVerbose, int n, int nitem, real **c1)
+static void average_acf(gmx_bool bVerbose, int n, int nitem, real** c1)
{
real c0;
int i, j;
{
c0 += c1[i][j];
}
- c1[0][j] = c0/nitem;
+ c1[0][j] = c0 / nitem;
}
}
static void norm_and_scale_vectors(int nframes, real c1[], real scale)
{
int j, m;
- real *rij;
+ real* rij;
for (j = 0; (j < nframes); j++)
{
- rij = &(c1[j*DIM]);
+ rij = &(c1[j * DIM]);
unitv(rij, rij);
for (m = 0; (m < DIM); m++)
{
}
/*! \brief Debugging */
-static void dump_tmp(char *s, int n, real c[])
+static void dump_tmp(char* s, int n, real c[])
{
- FILE *fp;
+ FILE* fp;
int i;
fp = gmx_ffopen(s, "w");
}
/*! \brief High level ACF routine. */
-static void do_four_core(unsigned long mode, int nframes,
- real c1[], real csum[], real ctmp[])
+static void do_four_core(unsigned long mode, int nframes, real c1[], real csum[], real ctmp[])
{
- real *cfour;
- char buf[32];
- real fac;
- int j, m, m1;
+ real* cfour;
+ char buf[32];
+ real fac;
+ int j, m, m1;
snew(cfour, nframes);
low_do_four_core(nframes, ctmp, cfour, enCos);
for (j = 0; (j < nframes); j++)
{
- c1[j] = cfour[j];
+ c1[j] = cfour[j];
}
/* Sine term of AC function */
low_do_four_core(nframes, ctmp, cfour, enSin);
for (j = 0; (j < nframes); j++)
{
- c1[j] += cfour[j];
+ c1[j] += cfour[j];
csum[j] = c1[j];
}
}
*/
for (j = 0; (j < nframes); j++)
{
- csum[j] = -0.5*(nframes-j);
+ csum[j] = -0.5 * (nframes - j);
}
/***** DIAGONAL ELEMENTS ************/
/* Copy the vector data in a linear array */
for (j = 0; (j < nframes); j++)
{
- ctmp[j] = gmx::square(c1[DIM*j+m]);
+ ctmp[j] = gmx::square(c1[DIM * j + m]);
}
if (debug)
{
fac = 1.5;
for (j = 0; (j < nframes); j++)
{
- csum[j] += fac*(cfour[j]);
+ csum[j] += fac * (cfour[j]);
}
}
/******* OFF-DIAGONAL ELEMENTS **********/
for (m = 0; (m < DIM); m++)
{
/* Copy the vector data in a linear array */
- m1 = (m+1) % DIM;
+ m1 = (m + 1) % DIM;
for (j = 0; (j < nframes); j++)
{
- ctmp[j] = c1[DIM*j+m]*c1[DIM*j+m1];
+ ctmp[j] = c1[DIM * j + m] * c1[DIM * j + m1];
}
if (debug)
fac = 3.0;
for (j = 0; (j < nframes); j++)
{
- csum[j] += fac*cfour[j];
+ csum[j] += fac * cfour[j];
}
}
}
/* Copy the vector data in a linear array */
for (j = 0; (j < nframes); j++)
{
- ctmp[j] = c1[DIM*j+m];
+ ctmp[j] = c1[DIM * j + m];
}
low_do_four_core(nframes, ctmp, cfour, enNorm);
for (j = 0; (j < nframes); j++)
sfree(cfour);
for (j = 0; (j < nframes); j++)
{
- c1[j] = csum[j]/static_cast<real>(nframes-j);
+ c1[j] = csum[j] / static_cast<real>(nframes - j);
}
}
-void low_do_autocorr(const char *fn, const gmx_output_env_t *oenv, const char *title,
- int nframes, int nitem, int nout, real **c1,
- real dt, unsigned long mode, int nrestart,
- gmx_bool bAver, gmx_bool bNormalize,
- gmx_bool bVerbose, real tbeginfit, real tendfit,
- int eFitFn)
+void low_do_autocorr(const char* fn,
+ const gmx_output_env_t* oenv,
+ const char* title,
+ int nframes,
+ int nitem,
+ int nout,
+ real** c1,
+ real dt,
+ unsigned long mode,
+ int nrestart,
+ gmx_bool bAver,
+ gmx_bool bNormalize,
+ gmx_bool bVerbose,
+ real tbeginfit,
+ real tendfit,
+ int eFitFn)
{
- FILE *fp, *gp = nullptr;
- int i;
- real *csum;
- real *ctmp, *fit;
- real sum, Ct2av, Ctav;
- gmx_bool bFour = acf.bFour;
+ FILE * fp, *gp = nullptr;
+ int i;
+ real* csum;
+ real * ctmp, *fit;
+ real sum, Ct2av, Ctav;
+ gmx_bool bFour = acf.bFour;
/* Check flags and parameters */
nout = get_acfnout();
if (nout == -1)
{
- nout = acf.nout = (nframes+1)/2;
+ nout = acf.nout = (nframes + 1) / 2;
}
else if (nout > nframes)
{
if (MODE(eacCos) && MODE(eacVector))
{
- gmx_fatal(FARGS, "Incompatible options bCos && bVector (%s, %d)",
- __FILE__, __LINE__);
+ gmx_fatal(FARGS, "Incompatible options bCos && bVector (%s, %d)", __FILE__, __LINE__);
}
if ((MODE(eacP3) || MODE(eacRcross)) && bFour)
{
{
printf("Will calculate %s of %d thingies for %d frames\n",
title ? title : "autocorrelation", nitem, nframes);
- printf("bAver = %s, bFour = %s bNormalize= %s\n",
- gmx::boolToString(bAver), gmx::boolToString(bFour),
- gmx::boolToString(bNormalize));
+ printf("bAver = %s, bFour = %s bNormalize= %s\n", gmx::boolToString(bAver),
+ gmx::boolToString(bFour), gmx::boolToString(bNormalize));
printf("mode = %lu, dt = %g, nrestart = %d\n", mode, dt, nrestart);
}
/* Allocate temp arrays */
*/
for (int i = 0; i < nitem; i++)
{
- if (bVerbose && (((i % 100) == 0) || (i == nitem-1)))
+ if (bVerbose && (((i % 100) == 0) || (i == nitem - 1)))
{
- fprintf(stderr, "\rThingie %d", i+1);
+ fprintf(stderr, "\rThingie %d", i + 1);
fflush(stderr);
}
sum = print_and_integrate(fp, nout, dt, c1[i], nullptr, 1);
if (debug)
{
- fprintf(debug,
- "CORRelation time (integral over corrfn %d): %g (ps)\n",
- i, sum);
+ fprintf(debug, "CORRelation time (integral over corrfn %d): %g (ps)\n", i, sum);
}
}
- Ctav += sum;
- Ct2av += sum*sum;
+ Ctav += sum;
+ Ct2av += sum * sum;
if (debug)
{
fprintf(gp, "%5d %.3f\n", i, sum);
}
if (nitem > 1)
{
- Ctav /= nitem;
+ Ctav /= nitem;
Ct2av /= nitem;
- printf("Average correlation time %.3f Std. Dev. %.3f Error %.3f (ps)\n",
- Ctav, std::sqrt((Ct2av - gmx::square(Ctav))),
- std::sqrt((Ct2av - gmx::square(Ctav))/(nitem-1)));
+ printf("Average correlation time %.3f Std. Dev. %.3f Error %.3f (ps)\n", Ctav,
+ std::sqrt((Ct2av - gmx::square(Ctav))),
+ std::sqrt((Ct2av - gmx::square(Ctav)) / (nitem - 1)));
}
}
if (fp)
}
/*! \brief Legend for selecting Legendre polynomials. */
-static const char *Leg[] = { nullptr, "0", "1", "2", "3", nullptr };
+static const char* Leg[] = { nullptr, "0", "1", "2", "3", nullptr };
-t_pargs *add_acf_pargs(int *npargs, t_pargs *pa)
+t_pargs* add_acf_pargs(int* npargs, t_pargs* pa)
{
- t_pargs acfpa[] = {
- { "-acflen", FALSE, etINT, {&acf.nout},
+ t_pargs acfpa[] = {
+ { "-acflen",
+ FALSE,
+ etINT,
+ { &acf.nout },
"Length of the ACF, default is half the number of frames" },
- { "-normalize", FALSE, etBOOL, {&acf.bNormalize},
- "Normalize ACF" },
- { "-fftcorr", FALSE, etBOOL, {&acf.bFour},
+ { "-normalize", FALSE, etBOOL, { &acf.bNormalize }, "Normalize ACF" },
+ { "-fftcorr",
+ FALSE,
+ etBOOL,
+ { &acf.bFour },
"HIDDENUse fast fourier transform for correlation function" },
- { "-nrestart", FALSE, etINT, {&acf.nrestart},
+ { "-nrestart",
+ FALSE,
+ etINT,
+ { &acf.nrestart },
"HIDDENNumber of frames between time origins for ACF when no FFT is used" },
- { "-P", FALSE, etENUM, {Leg},
- "Order of Legendre polynomial for ACF (0 indicates none)" },
- { "-fitfn", FALSE, etENUM, {s_ffn},
- "Fit function" },
- { "-beginfit", FALSE, etREAL, {&acf.tbeginfit},
+ { "-P", FALSE, etENUM, { Leg }, "Order of Legendre polynomial for ACF (0 indicates none)" },
+ { "-fitfn", FALSE, etENUM, { s_ffn }, "Fit function" },
+ { "-beginfit",
+ FALSE,
+ etREAL,
+ { &acf.tbeginfit },
"Time where to begin the exponential fit of the correlation function" },
- { "-endfit", FALSE, etREAL, {&acf.tendfit},
- "Time where to end the exponential fit of the correlation function, -1 is until the end" },
+ { "-endfit",
+ FALSE,
+ etREAL,
+ { &acf.tendfit },
+ "Time where to end the exponential fit of the correlation function, -1 is until the "
+ "end" },
};
- t_pargs *ppa;
+ t_pargs* ppa;
int i, npa;
npa = asize(acfpa);
- snew(ppa, *npargs+npa);
+ snew(ppa, *npargs + npa);
for (i = 0; (i < *npargs); i++)
{
ppa[i] = pa[i];
}
for (i = 0; (i < npa); i++)
{
- ppa[*npargs+i] = acfpa[i];
+ ppa[*npargs + i] = acfpa[i];
}
(*npargs) += npa;
return ppa;
}
-void do_autocorr(const char *fn, const gmx_output_env_t *oenv, const char *title,
- int nframes, int nitem, real **c1,
- real dt, unsigned long mode, gmx_bool bAver)
+void do_autocorr(const char* fn,
+ const gmx_output_env_t* oenv,
+ const char* title,
+ int nframes,
+ int nitem,
+ real** c1,
+ real dt,
+ unsigned long mode,
+ gmx_bool bAver)
{
if (!bACFinit)
{
switch (acf.P)
{
- case 1:
- mode = mode | eacP1;
- break;
- case 2:
- mode = mode | eacP2;
- break;
- case 3:
- mode = mode | eacP3;
- break;
- default:
- break;
+ case 1: mode = mode | eacP1; break;
+ case 2: mode = mode | eacP2; break;
+ case 3: mode = mode | eacP3; break;
+ default: break;
}
- low_do_autocorr(fn, oenv, title, nframes, nitem, acf.nout, c1, dt, mode,
- acf.nrestart, bAver, acf.bNormalize,
- bDebugMode(), acf.tbeginfit, acf.tendfit,
- acf.fitfn);
+ low_do_autocorr(fn, oenv, title, nframes, nitem, acf.nout, c1, dt, mode, acf.nrestart, bAver,
+ acf.bNormalize, bDebugMode(), acf.tbeginfit, acf.tendfit, acf.fitfn);
}
int get_acfnout()
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff