*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
{
if (x > 0)
{
- return pow(x, y);
+ return std::pow(x, y);
}
else
{
static void chomp(char *buf)
{
- int len = strlen(buf);
+ int len = std::strlen(buf);
while ((len > 0) && (buf[len-1] == '\n'))
{
static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
{
gmx_bool *bE;
- int n, k, kk, j, i, nmatch, nind, nss;
+ int k, kk, j, i, nmatch, nind, nss;
int *set;
gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
char *ptr, buf[STRLEN];
{
newnm[k] = gmx_strdup(nm[k].name);
/* Insert dashes in all the names */
- while ((ptr = strchr(newnm[k], ' ')) != NULL)
+ while ((ptr = std::strchr(newnm[k], ' ')) != NULL)
{
*ptr = '-';
}
bLong = FALSE;
for (kk = k; kk < k+4; kk++)
{
- if (kk < nre && strlen(nm[kk].name) > 14)
+ if (kk < nre && std::strlen(nm[kk].name) > 14)
{
bLong = TRUE;
}
trim(buf);
/* Empty line means end of input */
- bEOF = (strlen(buf) == 0);
+ bEOF = (std::strlen(buf) == 0);
if (!bEOF)
{
ptr = buf;
else
{
/* Now try to read a string */
- i = strlen(ptr);
nmatch = 0;
for (nind = 0; nind < nre; nind++)
{
}
if (nmatch == 0)
{
- i = strlen(ptr);
+ i = std::strlen(ptr);
nmatch = 0;
for (nind = 0; nind < nre; nind++)
{
}
}
/* Look for the first space, and remove spaces from there */
- if ((ptr = strchr(ptr, ' ')) != NULL)
+ if ((ptr = std::strchr(ptr, ' ')) != NULL)
{
trim(ptr);
}
}
- while (!bEOF && (ptr && (strlen(ptr) > 0)));
+ while (!bEOF && (ptr && (std::strlen(ptr) > 0)));
}
}
{
gmx_mtop_t mtop;
int natoms;
- t_iatom *iatom;
matrix box;
/* all we need is the ir to be able to write the label */
rav += sqr(rav3[j]);
rsum += mypow(rt[j], -6);
}
- rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]);
- rav = max(0.0, mypow(rav, -sixth)-bounds[i]);
+ rsum = std::max(0.0, mypow(rsum, -sixth)-bounds[i]);
+ rav = std::max(0.0, mypow(rav, -sixth)-bounds[i]);
sumt += rsum;
sumaver += rav;
{
sumaver += sqr(violaver[j]/nframes);
}
- sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
+ sumaver = std::max(0.0, mypow(sumaver, minsixth)-bounds[i]);
sumt += sumaver;
- sum = max(sum, sumaver);
+ sum = std::max(sum, sumaver);
fprintf(vout, "%10d %10.5e\n", index[i], sumaver);
}
#ifdef DEBUG
{
int nb, i, f, nee;
double sum, sum2, sump, see2;
- gmx_int64_t steps, np, p, bound_nb;
+ gmx_int64_t np, p, bound_nb;
enerdat_t *ed;
exactsum_t *es;
gmx_bool bAllZero;
edat->s[i].av = sum/np;
if (ed->bExactStat)
{
- edat->s[i].rmsd = sqrt(sum2/np);
+ edat->s[i].rmsd = std::sqrt(sum2/np);
}
else
{
- edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
+ edat->s[i].rmsd = std::sqrt(sum2/np - dsqr(edat->s[i].av));
}
if (edat->nframes > 1)
}
if (nee > 0)
{
- edat->s[i].ee = sqrt(see2/nee);
+ edat->s[i].ee = std::sqrt(see2/nee);
}
else
{
/* Remove the drift by performing a fit to y = ax+b.
Uses 5 iterations. */
int i, j, k;
- double delta, d, sd, sd2;
+ double delta;
edat->npoints = edat->nframes;
edat->nsteps = edat->nframes;
int nbmin, int nbmax)
{
int i, j;
- double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
+ double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, varet;
double NANO3;
enum {
eVol, eEnth, eTemp, eEtot, eNR
cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
}
/* Total energy */
- et = varet = NOTSET;
if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
{
/* Only compute cv in constant volume runs, which we can test
by checking whether the enthalpy was computed.
*/
- et = edat->s[ii[eEtot]].av;
varet = sqr(edat->s[ii[eEtot]].rmsd);
cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
}
/* Check out the printed manual for equations! */
double Dt, aver, stddev, errest, delta_t, totaldrift;
enerdata_t *esum = NULL;
- real xxx, integral, intBulk, Temp = 0, Pres = 0;
- real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
+ real integral, intBulk, Temp = 0, Pres = 0;
+ real pr_aver, pr_stddev, pr_errest;
double beta = 0, expE, expEtot, *fee = NULL;
gmx_int64_t nsteps;
int nexact, nnotexact;
- double x1m, x1mk;
- int i, j, k, nout;
- real chi2;
+ int i, j, nout;
char buf[256], eebuf[100];
nsteps = step - start_step + 1;
expE = 0;
for (j = 0; (j < edat->nframes); j++)
{
- expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
+ expE += std::exp(beta*(edat->s[i].ener[j] - aver)/nmol);
}
if (bSum)
{
expEtot += expE/edat->nframes;
}
- fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
+ fee[i] = std::log(expE/edat->nframes)/beta + aver/nmol;
}
- if (strstr(leg[i], "empera") != NULL)
+ if (std::strstr(leg[i], "empera") != NULL)
{
Temp = aver;
}
- else if (strstr(leg[i], "olum") != NULL)
+ else if (std::strstr(leg[i], "olum") != NULL)
{
Vaver = aver;
}
- else if (strstr(leg[i], "essure") != NULL)
+ else if (std::strstr(leg[i], "essure") != NULL)
{
Pres = aver;
}
if (bFee)
{
fprintf(stdout, " %10g %10g\n",
- log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
+ std::log(expEtot)/beta + esum->s[0].av/nmol, std::log(expEtot)/beta);
}
else
{
/*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
/* Do it for shear viscosity */
- strcpy(buf, "Shear Viscosity");
+ std::strcpy(buf, "Shear Viscosity");
low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
(edat->nframes+1)/2, eneset, Dt,
eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
/* Now for bulk viscosity */
- strcpy(buf, "Bulk Viscosity");
+ std::strcpy(buf, "Bulk Viscosity");
low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
(edat->nframes+1)/2, &(eneset[11]), Dt,
eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
{
if (bFluct)
{
- strcpy(buf, "Autocorrelation of Energy Fluctuations");
+ std::strcpy(buf, "Autocorrelation of Energy Fluctuations");
}
else
{
- strcpy(buf, "Energy Autocorrelation");
+ std::strcpy(buf, "Energy Autocorrelation");
}
#if 0
do_autocorr(corrfn, oenv, buf, edat->nframes,
};
FILE *fp;
ener_file_t enx;
- int nre, timecheck, step, nenergy, nenergy2, maxenergy;
+ int timecheck, nenergy, nenergy2, maxenergy;
int i, j;
gmx_bool bCont;
real aver, beta;
srenew(eneset2[i], maxenergy);
}
}
+ GMX_RELEASE_ASSERT(time != NULL, "trying to dereference NULL time pointer");
+
if (fr->t != time[nenergy2])
{
fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
edat->nframes, nenergy2);
}
- nenergy = min(edat->nframes, nenergy2);
+ nenergy = std::min(edat->nframes, nenergy2);
/* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
fp = NULL;
for (j = 0; j < nenergy; j++)
{
dE = eneset2[i][j] - edat->s[i].ener[j];
- sum += exp(-dE*beta);
+ sum += std::exp(-dE*beta);
if (fp)
{
fprintf(fp, "%10g %10g %10g\n",
- time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
+ time[j], dE, -BOLTZ*reftemp*std::log(sum/(j+1)) );
}
}
- aver = -BOLTZ*reftemp*log(sum/nenergy);
+ aver = -BOLTZ*reftemp*std::log(sum/nenergy);
fprintf(stdout, "%-24s %10g\n", leg[i], aver);
}
if (fp)
const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
char buf[STRLEN];
- gmx_bool first = FALSE;
int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
int i, j, k;
/* coll data */
- double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
+ double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0;
static int setnr = 0;
double *native_lambda_vec = NULL;
const char **lambda_components = NULL;
int n_lambda_vec = 0;
- gmx_bool changing_lambda = FALSE;
- int lambda_fep_state;
/* now count the blocks & handle the global dh data */
for (i = 0; i < fr->nblock; i++)
}
/* read the data from the DHCOLL block */
- temp = fr->block[i].sub[0].dval[0];
- start_time = fr->block[i].sub[0].dval[1];
- delta_time = fr->block[i].sub[0].dval[2];
+ temp = fr->block[i].sub[0].dval[0];
+ start_time = fr->block[i].sub[0].dval[1];
+ delta_time = fr->block[i].sub[0].dval[2];
start_lambda = fr->block[i].sub[0].dval[3];
- delta_lambda = fr->block[i].sub[0].dval[4];
- changing_lambda = (delta_lambda != 0);
if (fr->block[i].nsub > 1)
{
- lambda_fep_state = fr->block[i].sub[1].ival[0];
if (n_lambda_vec == 0)
{
n_lambda_vec = fr->block[i].sub[1].ival[1];
}
}
}
- (*samples) += (int)(sum/nblock_hist);
+ (*samples) += static_cast<int>(sum/nblock_hist);
}
else
{
/* raw dh */
int len = 0;
- char **setnames = NULL;
- int nnames = nblock_dh;
for (i = 0; i < fr->nblock; i++)
{
if (j == 1 && ir->bExpanded)
{
- fprintf(*fp_dhdl, "%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
+ fprintf(*fp_dhdl, "%4d", static_cast<int>(value)); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
}
else
{
FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
FILE *fp_dhdl = NULL;
- FILE **drout;
ener_file_t fp;
int timecheck = 0;
gmx_mtop_t mtop;
gmx_localtop_t *top = NULL;
t_inputrec ir;
- t_energy **ee;
enerdata_t edat;
gmx_enxnm_t *enm = NULL;
t_enxframe *frame, *fr = NULL;
int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
int nbounds = 0, npairs;
gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
- gmx_bool bFoundStart, bCont, bEDR, bVisco;
- double sum, sumaver, sumt, ener, dbl;
+ gmx_bool bFoundStart, bCont, bVisco;
+ double sum, sumaver, sumt, dbl;
double *time = NULL;
real Vaver;
int *set = NULL, i, j, k, nset, sss;
gmx_bool *bIsEner = NULL;
char **pairleg, **odtleg, **otenleg;
char **leg = NULL;
- char **nms;
char *anm_j, *anm_k, *resnm_j, *resnm_k;
int resnr_j, resnr_k;
const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
{
for (i = 0; i < nre; i++)
{
- if (strstr(enm[i].name, setnm[j]))
+ if (std::strstr(enm[i].name, setnm[j]))
{
set[j] = i;
break;
{
for (j = 0; j < i; j++)
{
- if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
+ if (std::strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
{
break;
}
}
if (j == i)
{
- strcat(buf, ", (");
- strcat(buf, enm[set[i]].unit);
- strcat(buf, ")");
+ std::strcat(buf, ", (");
+ std::strcat(buf, enm[set[i]].unit);
+ std::strcat(buf, ")");
}
}
out = xvgropen(opt2fn("-o", NFILE, fnm), "GROMACS Energies", "Time (ps)", buf,
if (edat.nframes % 1000 == 0)
{
srenew(edat.step, edat.nframes+1000);
- memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
+ std::memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
srenew(edat.steps, edat.nframes+1000);
- memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
+ std::memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
srenew(edat.points, edat.nframes+1000);
- memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
+ std::memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
for (i = 0; i < nset; i++)
{
srenew(edat.s[i].ener, edat.nframes+1000);
- memset(&(edat.s[i].ener[edat.nframes]), 0,
- 1000*sizeof(edat.s[i].ener[0]));
+ std::memset(&(edat.s[i].ener[edat.nframes]), 0,
+ 1000*sizeof(edat.s[i].ener[0]));
srenew(edat.s[i].es, edat.nframes+1000);
- memset(&(edat.s[i].es[edat.nframes]), 0,
- 1000*sizeof(edat.s[i].es[0]));
+ std::memset(&(edat.s[i].es[edat.nframes]), 0,
+ 1000*sizeof(edat.s[i].es[0]));
}
}
*******************************************/
if (ndisre > 0)
{
+ GMX_RELEASE_ASSERT(blk_disre != NULL, "Trying to dereference NULL blk_disre pointer");
#ifndef GMX_DOUBLE
float *disre_rt = blk_disre->sub[0].fval;
float *disre_rm3tav = blk_disre->sub[1].fval;
print_time(out, fr->t);
print1(out, bDp, fr->ener[set[0]].e);
print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
- print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
+ print1(out, bDp, std::sqrt(fr->ener[set[0]].eav/fr->nsum));
fprintf(out, "\n");
}
}
vals = blk->sub[0].dval;
#endif
- if (blk->sub[0].nr != (size_t)nor)
+ if (blk->sub[0].nr != nor)
{
gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
}
vals = blk->sub[0].dval;
#endif
- if (blk->sub[0].nr != (size_t)(nex*12))
+ if (blk->sub[0].nr != nex*12)
{
gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
blk->sub[0].nr/12, nex);
}
for (i = 0; i < nor; i++)
{
- fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr));
+ fprintf(out, "%5d %g\n", or_label[i], std::sqrt(odrms[i]/norfr));
}
xvgrclose(out);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff