*/
#include "gmxpre.h"
-#include <assert.h>
-#include <stdlib.h>
+#include <cmath>
+#include <cstdlib>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/topology/index.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
-#define SQR(x) (pow(x, 2.0))
#define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
static void index_atom2mol(int *n, int *index, t_block *mols)
}
}
- if (fabs(qall) > 0.01)
+ if (std::abs(qall) > 0.01)
{
printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n", qall);
bNEU = FALSE;
{
int i;
- real corint;
real deltat = 0.0;
real rfr;
real sigma = 0.0;
real sigma_ret = 0.0;
- corint = 0.0;
if (nfr > 1)
{
while (i < nfr)
{
+ real corint;
- rfr = (real) (nfr/nshift-i/nshift);
+ rfr = static_cast<real>(nfr+i)/nshift;
cacf[i] /= rfr;
if (time[vfr[i]] <= time[vfr[ei]])
{
int i;
- real rtmp;
- real rfr;
-
fprintf(fmjdsp, "#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n", prefactor);
-
-
for (i = 0; i < nfr; i++)
{
static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,
FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,
int ePBC, t_topology top, t_trxframe fr, real temp,
- real trust, real bfit, real efit, real bvit, real evit,
+ real bfit, real efit, real bvit, real evit,
t_trxstatus *status, int isize, int nmols, int nshift,
atom_id *index0, int indexm[], real mass2[],
real qmol[], real eps_rf, const output_env_t oenv)
{
- int i, j, k, l, f;
- int valloc, nalloc, nfr, nvfr, m, itrust = 0;
+ int i, j;
+ int valloc, nalloc, nfr, nvfr;
int vshfr;
real *xshfr = NULL;
int *vfr = NULL;
real refr = 0.0;
- real rfr = 0.0;
real *cacf = NULL;
real *time = NULL;
real *djc = NULL;
real volume;
real volume_av = 0.0;
real dk_s, dk_d, dk_f;
- real dm_s, dm_d;
real mj = 0.0;
real mj2 = 0.0;
real mjd = 0.0;
real *dsp2 = NULL;
real t0;
real rtmp;
- real qtmp;
-
-
- rvec tmp;
- rvec *mtrans = NULL;
+ rvec tmp;
+ rvec *mtrans = NULL;
/*
* Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
do
{
- refr = (real)(nfr+1);
+ refr = (nfr+1);
if (nfr >= nalloc)
{
xshfr[i] = 0.0;
}
}
- assert(time != NULL);
-
+ GMX_RELEASE_ASSERT(time != NULL, "Memory not allocated correctly - time array is NULL");
if (nfr == 0)
{
if (v0 != NULL)
{
- trust *= (real) nvfr;
- itrust = (int) trust;
-
if (bINT)
{
for (i = ii; i <= ie; i++)
{
- xfit[i-ii] = log(time[vfr[i]]);
- rtmp = fabs(cacf[i]);
- yfit[i-ii] = log(rtmp);
+ xfit[i-ii] = std::log(time[vfr[i]]);
+ rtmp = std::abs(cacf[i]);
+ yfit[i-ii] = std::log(rtmp);
}
lsq_y_ax_b(ie-ii, xfit, yfit, &sigma, &malt, &err, &rest);
- malt = exp(malt);
+ malt = std::exp(malt);
sigma += 1.0;
if (bACF && (ii < nvfr))
{
fprintf(stderr, "Integral and integrated fit to the current acf yields at t=%f:\n", time[vfr[ii]]);
- fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*pow(time[vfr[ii]], sigma), sgk);
+ fprintf(stderr, "sigma=%8.3f (pure integral: %.3f)\n", sgk-malt*std::pow(time[vfr[ii]], sigma), sgk);
}
if (ei > bi)
static int nshift = 1000;
static real temp = 300.0;
- static real trust = 0.25;
static real eps_rf = 0.0;
static gmx_bool bNoJump = TRUE;
static real bfit = 100.0;
"Begin of the fit of the current autocorrelation function to a*t^b."},
{ "-evit", FALSE, etREAL, {&evit},
"End of the fit of the current autocorrelation function to a*t^b."},
- { "-tr", FALSE, etREAL, {&trust},
- "Fraction of the trajectory taken into account for the integral."},
{ "-temp", FALSE, etREAL, {&temp},
"Temperature for calculating epsilon."}
};
int isize;
t_trxstatus *status;
int flags = 0;
- gmx_bool bTop;
- gmx_bool bNEU;
gmx_bool bACF;
gmx_bool bINT;
int ePBC = -1;
- int natoms;
int nmols;
- int i, j, k = 0, l;
- int step;
- real t;
- real lambda;
+ int i;
real *qmol;
FILE *outf = NULL;
- FILE *outi = NULL;
- FILE *tfile = NULL;
FILE *mcor = NULL;
FILE *fmj = NULL;
FILE *fmd = NULL;
bACF = opt2bSet("-caf", NFILE, fnm);
bINT = opt2bSet("-mc", NFILE, fnm);
- bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
+ read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, &vtop, box, TRUE);
snew(mass2, top.atoms.nr);
snew(qmol, top.atoms.nr);
- bNEU = precalc(top, mass2, qmol);
+ precalc(top, mass2, qmol);
snew(indexm, isize);
* and calculates the requested quantities */
dielectric(fmj, fmd, outf, fcur, mcor, fmjdsp, bNoJump, bACF, bINT, ePBC, top, fr,
- temp, trust, bfit, efit, bvit, evit, status, isize, nmols, nshift,
+ temp, bfit, efit, bvit, evit, status, isize, nmols, nshift,
index0, indexm, mass2, qmol, eps_rf, oenv);
xvgrclose(fmj);