}
fprintf(log, "Min data: %10g Max data: %10g\n", minx, maxx);
}
- dx = (double)npoints/(maxx-minx);
+ dx = npoints/(maxx-minx);
if (debug)
{
fprintf(debug,
for (i = 0; (i < nangles); i++)
{
real dd = angles[cur][i];
- angles[cur][i] = std::atan2(sin(dd), cos(dd));
+ angles[cur][i] = std::atan2(std::sin(dd), std::cos(dd));
}
}
else
{
for (j = i+1; j < nn; j++)
{
- x = (int)(fac*mat[i][j]+0.5);
+ x = static_cast<int>(fac*mat[i][j]+0.5);
if (x <= 100)
{
histo[x]++;
/* end fixing aromatics */
/* Special case for Pro, has no H */
- if (strcmp(thisres, "PRO") == 0)
+ if (std::strcmp(thisres, "PRO") == 0)
{
atm.H = atm.Cn[4];
}
*/
#include "gmxpre.h"
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/legacyheaders/copyrite.h"
#include "gromacs/legacyheaders/macros.h"
#include "gromacs/legacyheaders/typedefs.h"
ptr = fgets2(buf, 255, fp);
if (ptr)
{
- if (strstr(buf, "Intermolecular") != NULL)
+ if (std::strstr(buf, "Intermolecular") != NULL)
{
- ptr = strchr(buf, '=');
+ ptr = std::strchr(buf, '=');
sscanf(ptr+1, "%lf", &e);
pdbf->edocked = e;
}
- else if (strstr(buf, "Estimated Free") != NULL)
+ else if (std::strstr(buf, "Estimated Free") != NULL)
{
- ptr = strchr(buf, '=');
+ ptr = std::strchr(buf, '=');
sscanf(ptr+1, "%lf", &e);
pdbf->efree = e;
}
char buf[256], name[256];
gmx_bool bExist;
- strcpy(buf, fn);
- buf[strlen(buf)-4] = '\0';
- maxpdbf = 100;
+ std::strcpy(buf, fn);
+ buf[std::strlen(buf)-4] = '\0';
+ maxpdbf = 100;
snew(pdbf, maxpdbf);
i = 0;
do
static real rmsd_dist(t_pdbfile *pa, t_pdbfile *pb, gmx_bool bRMSD)
{
int i;
- real rmsd, dist;
+ real rmsd;
rvec acm, bcm, dx;
if (bRMSD)
rvec_sub(pa->x[i], pb->x[i], dx);
rmsd += iprod(dx, dx);
}
- rmsd = sqrt(rmsd/pa->atoms.nr);
+ rmsd = std::sqrt(rmsd/pa->atoms.nr);
}
else
{
- dist = 0;
clear_rvec(acm);
clear_rvec(bcm);
for (i = 0; (i < pa->atoms.nr); i++)
*/
#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/fileio/confio.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
{
int i;
- double hwkT, w, dS, S = 0;
+ double hwkT, dS, S = 0;
double hbar, lambda;
hbar = PLANCK1/(2*M_PI);
{
if (eigval[i] > 0)
{
+ double w;
lambda = eigval[i]*AMU;
- w = sqrt(BOLTZMANN*temp/lambda)/NANO;
+ w = std::sqrt(BOLTZMANN*temp/lambda)/NANO;
hwkT = (hbar*w)/(BOLTZMANN*temp);
- dS = (hwkT/gmx_expm1(hwkT) - gmx_log1p(-exp(-hwkT)));
+ dS = (hwkT/gmx_expm1(hwkT) - gmx_log1p(-std::exp(-hwkT)));
S += dS;
if (debug)
{
else
{
fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
- w = 0;
}
}
fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
real *eigval, real temp)
{
double dd, deter;
- int *indx;
- int i, j, k, m;
- char buf[256];
+ int i;
double hbar, kt, kteh, S;
hbar = PLANCK1/(2*M_PI);
kt = BOLTZMANN*temp;
- kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
+ kteh = kt*std::exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
if (debug)
{
fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
for (i = 0; (i < n-nskip); i++)
{
dd = 1+kteh*eigval[i];
- deter += log(dd);
+ deter += std::log(dd);
}
S = 0.5*RGAS*deter;
return 1.0;
}
- sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
+ sp = 0.2*std::exp(std::log(10.0)*std::ceil(std::log(range)/std::log(10.0)));
while (range/sp < minticks-1)
{
sp = sp/2;
{
FILE *out;
int g, s, i;
- real min, max, xsp, ysp;
+ real ymin, ymax, xsp, ysp;
out = gmx_ffopen(file, "w");
if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
{
if (y)
{
- min = y[g][0];
- max = y[g][0];
+ ymin = y[g][0];
+ ymax = y[g][0];
for (i = 0; i < n; i++)
{
- if (y[g][i] < min)
+ if (y[g][i] < ymin)
{
- min = y[g][i];
+ ymin = y[g][i];
}
- if (y[g][i] > max)
+ if (y[g][i] > ymax)
{
- max = y[g][i];
+ ymax = y[g][i];
}
}
}
else
{
- min = sy[g][0][0];
- max = sy[g][0][0];
+ ymin = sy[g][0][0];
+ ymax = sy[g][0][0];
for (s = 0; s < nsetspergraph; s++)
{
for (i = 0; i < n; i++)
{
- if (sy[g][s][i] < min)
+ if (sy[g][s][i] < ymin)
{
- min = sy[g][s][i];
+ ymin = sy[g][s][i];
}
- if (sy[g][s][i] > max)
+ if (sy[g][s][i] > ymax)
{
- max = sy[g][s][i];
+ ymax = sy[g][s][i];
}
}
}
}
if (bZero)
{
- min = 0;
+ ymin = 0;
}
else
{
- min = min-0.1*(max-min);
+ ymin = ymin-0.1*(ymax-ymin);
}
- max = max+0.1*(max-min);
- xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
- ysp = tick_spacing(max-min, 3);
+ ymax = ymax+0.1*(ymax-ymin);
+ xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
+ ysp = tick_spacing(ymax-ymin, 3);
if (output_env_get_print_xvgr_codes(oenv))
{
fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
{
fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
- fprintf(out, "@ world ymin %g\n", min);
- fprintf(out, "@ world ymax %g\n", max);
+ fprintf(out, "@ world ymin %g\n", ymin);
+ fprintf(out, "@ world ymax %g\n", ymax);
}
fprintf(out, "@ view xmin 0.15\n");
fprintf(out, "@ view xmax 0.85\n");
fprintf(out, "@ xaxis tick major %g\n", xsp);
fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
fprintf(out, "@ xaxis ticklabel start type spec\n");
- fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp);
+ fprintf(out, "@ xaxis ticklabel start %g\n", std::ceil(ymin/xsp)*xsp);
fprintf(out, "@ yaxis tick major %g\n", ysp);
fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
fprintf(out, "@ yaxis ticklabel start type spec\n");
- fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp);
- if ((min < 0) && (max > 0))
+ fprintf(out, "@ yaxis ticklabel start %g\n", std::ceil(ymin/ysp)*ysp);
+ if ((ymin < 0) && (ymax > 0))
{
fprintf(out, "@ zeroxaxis bar on\n");
fprintf(out, "@ zeroxaxis bar linestyle 3\n");
{
for (i = 0; i < n; i++)
{
- if (bSplit && i > 0 && fabs(x[i]) < 1e-5)
+ if (bSplit && i > 0 && std::abs(x[i]) < 1e-5)
{
fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
real *eigval1, int neig1, real *eigval2, int neig2)
{
- int n, nrow;
+ int n;
int i, j, k;
double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
- n = min(n1, n2);
+ n = std::min(n1, n2);
- n = min(n, min(neig1, neig2));
+ n = std::min(n, std::min(neig1, neig2));
fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
sum1 = 0;
eigval1[i] = 0;
}
sum1 += eigval1[i];
- eigval1[i] = sqrt(eigval1[i]);
+ eigval1[i] = std::sqrt(eigval1[i]);
}
trace1 = sum1;
for (i = n; i < neig1; i++)
eigval2[i] = 0;
}
sum2 += eigval2[i];
- eigval2[i] = sqrt(eigval2[i]);
+ eigval2[i] = std::sqrt(eigval2[i]);
}
trace2 = sum2;
+
+
+ // The static analyzer appears to be confused by the fact that the loop below
+ // starts from n instead of 0. However, given all the complex code it's
+ // better to be safe than sorry, so we check it with an assert.
+ // If we are in this comparison routine in the first place, neig2 should not be 0,
+ // so eigval2 should always be a valid pointer.
+ GMX_RELEASE_ASSERT(eigval2 != NULL, "NULL pointer provided for eigval2");
+
for (i = n; i < neig2; i++)
{
trace2 += eigval2[i];
if (neig1 != n || neig2 != n)
{
fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
- (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
+ static_cast<int>(100*sum1/trace1+0.5), static_cast<int>(100*sum2/trace2+0.5));
}
fprintf(stdout, "Square root of the traces: %g and %g\n",
- sqrt(sum1), sqrt(sum2));
+ std::sqrt(sum1), std::sqrt(sum2));
sab = 0;
for (i = 0; i < n; i++)
}
fprintf(stdout, "The overlap of the covariance matrices:\n");
- fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
- tmp = 1-sab/sqrt(sum1*sum2);
+ fprintf(stdout, " normalized: %.3f\n", 1-std::sqrt(samsb2/(sum1+sum2)));
+ tmp = 1-sab/std::sqrt(sum1*sum2);
if (tmp < 0)
{
tmp = 0;
}
- fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
+ fprintf(stdout, " shape: %.3f\n", 1-std::sqrt(tmp));
}
real **mat;
int i, x1, y1, x, y, nlevels;
int nx, ny;
- real inp, *t_x, *t_y, max;
+ real inp, *t_x, *t_y, maxval;
t_rgb rlo, rhi;
snew(t_y, nvec2);
snew(mat, nx);
snew(t_x, nx);
- max = 0;
+ maxval = 0;
for (x1 = 0; x1 < nx; x1++)
{
snew(mat[x1], ny);
{
inp += iprod(eigvec1[x][i], eigvec2[y][i]);
}
- mat[x1][y1] = fabs(inp);
- if (mat[x1][y1] > max)
+ mat[x1][y1] = std::abs(inp);
+ if (mat[x1][y1] > maxval)
{
- max = mat[x1][y1];
+ maxval = mat[x1][y1];
}
}
}
nlevels = 41;
out = gmx_ffopen(matfile, "w");
write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
- nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
+ nx, ny, t_x, t_y, mat, 0.0, maxval, rlo, rhi, &nlevels);
gmx_ffclose(out);
}
atom_id *all_at;
matrix box;
rvec *xread, *x;
- real t, inp, **inprod = NULL, min = 0, max = 0;
+ real t, inp, **inprod = NULL;
char str[STRLEN], str2[STRLEN], **ylabel, *c;
real fact;
gmx_rmpbc_t gpbc = NULL;
if (projfile)
{
+ GMX_RELEASE_ASSERT(inprod != NULL, "inprod must be non-NULL if projfile is non-NULL");
snew(ylabel, noutvec);
for (v = 0; v < noutvec; v++)
{
oenv);
for (i = 0; i < nframes; i++)
{
- if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
+ if (bSplit && i > 0 && std::abs(inprod[noutvec][i]) < 1e-5)
{
fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
}
atoms.atomname[i] = &atnm;
atoms.atom[i].resind = i;
atoms.resinfo[i].name = &resnm;
- atoms.resinfo[i].nr = ceil(i*fact);
+ atoms.resinfo[i].nr = static_cast<int>(std::ceil(i*fact));
atoms.resinfo[i].ic = ' ';
x[i][XX] = inprod[0][i];
x[i][YY] = inprod[1][i];
}
if ( ( b4D || bSplit ) && bPDB)
{
+ GMX_RELEASE_ASSERT(inprod != NULL, "inprod must be non-NULL with 4D or split PDB output options");
+
out = gmx_ffopen(threedplotfile, "w");
fprintf(out, "HEADER %s\n", str);
if (b4D)
j = 0;
for (i = 0; i < atoms.nr; i++)
{
- if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
+ if (j > 0 && bSplit && std::abs(inprod[noutvec][i]) < 1e-5)
{
fprintf(out, "TER\n");
j = 0;
snew(pmax, noutvec_extr);
if (extreme == 0)
{
+ GMX_RELEASE_ASSERT(inprod != NULL, "inprod must be non-NULL");
fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
fprintf(stderr,
"%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
pmax[0] = extreme;
}
/* build format string for filename: */
- strcpy(str, extremefile); /* copy filename */
- c = strrchr(str, '.'); /* find where extention begins */
- strcpy(str2, c); /* get extention */
- sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
+ std::strcpy(str, extremefile); /* copy filename */
+ c = std::strrchr(str, '.'); /* find where extention begins */
+ std::strcpy(str2, c); /* get extention */
+ sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
for (v = 0; v < noutvec_extr; v++)
{
/* make filename using format string */
if (noutvec_extr == 1)
{
- strcpy(str2, extremefile);
+ std::strcpy(str2, extremefile);
}
else
{
real *eigval, int neig,
const output_env_t oenv)
{
- int nrow, g, v, i;
+ int g, v, i;
real *x, **y;
char str[STRLEN], **ylabel;
snew(y[g], natoms);
for (i = 0; i < natoms; i++)
{
- y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
+ y[g][i] = std::sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
}
}
write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
};
#define NPA asize(pa)
- FILE *out;
- int status, trjout;
t_topology top;
int ePBC = -1;
t_atoms *atoms = NULL;
rvec *xtop, *xref1, *xref2, *xrefp = NULL;
gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
- rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
+ rvec *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
matrix topbox;
- real xid, totmass, *sqrtm, *w_rls, t, lambda;
- int natoms, step;
+ real totmass, *sqrtm, *w_rls, t;
+ int natoms;
char *grpname;
const char *indexfile;
char title[STRLEN];
int i, j, d;
int nout, *iout, noutvec, *outvec, nfit;
- atom_id *index, *ifit;
+ atom_id *index = NULL, *ifit = NULL;
const char *VecFile, *Vec2File, *topfile;
const char *EigFile, *Eig2File;
const char *CompFile, *RmsfFile, *ProjOnVecFile;
OverlapFile = opt2fn_null("-over", NFILE, fnm);
InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
- bTop = fn2bTPX(topfile);
bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
|| FilterFile || ExtremeFile;
bFirstLastSet =
gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
}
}
+ else
+ {
+ nvec2 = 0;
+ neig2 = 0;
+ }
if (Eig2File != NULL)
{
proj_unit = "u\\S1/2\\Nnm";
for (i = 0; (i < natoms); i++)
{
- sqrtm[i] = sqrt(atoms->atom[index[i]].m);
+ sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
}
}
else
}
}
fprintf(stdout, "RMSD (without fit) between the two average structures:"
- " %.3f (nm)\n\n", sqrt(t/totmass));
+ " %.3f (nm)\n\n", std::sqrt(t/totmass));
}
if (last == -1)
{
/* make an index of first+(0,1,2) and last */
nout = bPDB3D ? 4 : 3;
- nout = min(last-first+1, nout);
+ nout = std::min(last-first+1, nout);
snew(iout, nout);
iout[0] = first-1;
iout[1] = first;
*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
{
if (t[0] > 0)
{
- x[i] = log(t[i]);
+ x[i] = std::log(t[i]);
}
}
}
fprintf(stdout, "First time is not larger than 0, using index number as time for power fit\n");
for (i = 0; i < n; i++)
{
- x[i] = gmx_log1p(i);
+ x[i] = gmx_log1p(static_cast<real>(i));
}
}
for (s = 0; s < nset; s++)
{
- i = 0;
for (i = 0; i < n && val[s][i] >= 0; i++)
{
- y[i] = log(val[s][i]);
+ y[i] = std::log(val[s][i]);
}
if (i < n)
{
}
lsq_y_ax_b(i, x, y, &a, &b, &r, &quality);
fprintf(stdout, "Power fit set %3d: error %.3f a %g b %g\n",
- s+1, quality, a, exp(b));
+ s+1, quality, a, std::exp(b));
}
sfree(y);
yyint = 0;
for (i = 0; i < n; i++)
{
- cosyint += cos(fac*i)*y[i];
+ cosyint += std::cos(fac*i)*y[i];
yyint += y[i]*y[i];
}
{
FILE *fp;
int i, s;
- double min, max;
+ double minval, maxval;
int nbin;
gmx_int64_t *histo;
- min = val[0][0];
- max = val[0][0];
+ minval = val[0][0];
+ maxval = val[0][0];
for (s = 0; s < nset; s++)
{
for (i = 0; i < n; i++)
{
- if (val[s][i] < min)
+ if (val[s][i] < minval)
{
- min = val[s][i];
+ minval = val[s][i];
}
- else if (val[s][i] > max)
+ else if (val[s][i] > maxval)
{
- max = val[s][i];
+ maxval = val[s][i];
}
}
}
- min = binwidth*floor(min/binwidth);
- max = binwidth*ceil(max/binwidth);
- if (min != 0)
+ minval = binwidth*std::floor(minval/binwidth);
+ maxval = binwidth*std::ceil(maxval/binwidth);
+ if (minval != 0)
{
- min -= binwidth;
+ minval -= binwidth;
}
- max += binwidth;
+ maxval += binwidth;
- nbin = (int)((max - min)/binwidth + 0.5) + 1;
+ nbin = static_cast<int>(((maxval - minval)/binwidth + 0.5) + 1);
fprintf(stderr, "Making distributions with %d bins\n", nbin);
snew(histo, nbin);
fp = xvgropen(distfile, "Distribution", "", "", oenv);
}
for (i = 0; i < n; i++)
{
- histo[(int)((val[s][i] - min)/binwidth + 0.5)]++;
+ histo[static_cast<int>((val[s][i] - minval)/binwidth + 0.5)]++;
}
for (i = 0; i < nbin; i++)
{
- fprintf(fp, " %g %g\n", min+i*binwidth, (double)histo[i]/(n*binwidth));
+ fprintf(fp, " %g %g\n", minval+i*binwidth, static_cast<double>(histo[i])/(n*binwidth));
}
if (s < nset-1)
{
{
snew(tmp, nset);
fprintf(fp, "@TYPE xydydy\n");
- edge = (int)(nset*0.05+0.5);
+ edge = static_cast<int>(nset*0.05+0.5);
fprintf(stdout, "Errorbars: discarding %d points on both sides: %d%%"
- " interval\n", edge, (int)(100*(nset-2*edge)/nset+0.5));
+ " interval\n", edge, static_cast<int>(100*(nset-2*edge)/nset+0.5));
}
else
{
}
if (avbar_opt == avbarSTDDEV)
{
- err = sqrt(var/nset);
+ err = std::sqrt(var/nset);
}
else
{
- err = sqrt(var/(nset*(nset-1)));
+ err = std::sqrt(var/(nset*(nset-1)));
}
fprintf(fp, " %g", err);
}
return 0;
}
- return sigma*sqrt(2*ss/tTotal);
+ return sigma*std::sqrt(2*ss/tTotal);
}
static void estimate_error(const char *eefile, int nb_min, int resol, int n,
xvgr_legend(fp, 2*nset, (const char**)leg, oenv);
sfree(leg);
- spacing = pow(2, 1.0/resol);
+ spacing = std::pow(2.0, 1.0/resol);
snew(tbs, n);
snew(ybs, n);
snew(fitsig, n);
nbr = nb_min;
while (nbr <= n)
{
- bs = n/(int)nbr;
+ bs = n/static_cast<int>(nbr);
if (bs != prev_bs)
{
nb = n/bs;
nbs++;
}
nbr *= spacing;
- nb = (int)(nbr+0.5);
prev_bs = bs;
}
if (sig[s] == 0)
* For a single exponential autocorrelation: ybs(tau1) = 2/e tau1
* From this we take our initial guess for tau1.
*/
- twooe = 2/exp(1);
+ twooe = 2.0/std::exp(1.0);
i = -1;
do
{
{
dens = 0.5*(tbs[i+1]/tbs[i-1] - 1);
}
- fitsig[i] = sqrt((tau_sig + tbs[i])/dens);
+ fitsig[i] = std::sqrt((tau_sig + tbs[i])/dens);
}
if (!bSingleExpFit)
/* We set the initial guess for tau2
* to halfway between tau1_est and the total time (on log scale).
*/
- fitparm[2] = sqrt(tau1_est*(n-1)*dt);
+ fitparm[2] = std::sqrt(tau1_est*(n-1)*dt);
do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv,
bDebugMode(), effnERREST, fitparm, 0,
NULL);
}
for (i = 0; i < nbs; i++)
{
- fprintf(fp, "%g %g %g\n", tbs[i], sig[s]*sqrt(ybs[i]/(n*dt)),
- sig[s]*sqrt(fit_function(effnERREST, fitparm, tbs[i])/(n*dt)));
+ fprintf(fp, "%g %g %g\n", tbs[i], sig[s]*std::sqrt(ybs[i]/(n*dt)),
+ sig[s]*std::sqrt(fit_function(effnERREST, fitparm, tbs[i])/(n*dt)));
}
if (bFitAc)
ac[i] = val[s][i] - av[s];
if (i > 0)
{
- fitsig[i] = sqrt(i);
+ fitsig[i] = std::sqrt(static_cast<real>(i));
}
else
{
/* Generate more or less appropriate sigma's */
for (i = 0; i <= fitlen; i++)
{
- fitsig[i] = sqrt(acint + dt*i);
+ fitsig[i] = std::sqrt(acint + dt*i);
}
ac_fit[0] = 0.5*acint;
for (i = 0; i < nbs; i++)
{
fprintf(fp, "%g %g\n", tbs[i],
- sig[s]*sqrt(fit_function(effnERREST, ac_fit, tbs[i]))/(n*dt));
+ sig[s]*std::sqrt(fit_function(effnERREST, ac_fit, tbs[i]))/(n*dt));
}
sfree(ac);
static void luzar_correl(int nn, real *time, int nset, real **val, real temp,
gmx_bool bError, real fit_start)
{
- const real tol = 1e-8;
real *kt;
- real weight, d2;
+ real d2;
int j;
please_cite(stdout, "Spoel2006b");
{
d2 += sqr(kt[j] - val[3][j]);
}
- fprintf(debug, "RMS difference in derivatives is %g\n", sqrt(d2/nn));
+ fprintf(debug, "RMS difference in derivatives is %g\n", std::sqrt(d2/nn));
}
analyse_corr(nn, time, val[0], val[2], kt, NULL, NULL, NULL, fit_start,
temp);
int f, s, i, j;
double *filt, sum, vf, fluc, fluctot;
- f = (int)(flen/(2*dt));
+ f = static_cast<int>(flen/(2*dt));
snew(filt, f+1);
filt[0] = 1;
sum = 1;
for (i = 1; i <= f; i++)
{
- filt[i] = cos(M_PI*dt*i/flen);
+ filt[i] = std::cos(M_PI*dt*i/flen);
sum += 2*filt[i];
}
for (i = 0; i <= f; i++)
}
fluc /= n - 2*f;
fluctot += fluc;
- fprintf(stdout, "Set %3d filtered fluctuation: %12.6e\n", s+1, sqrt(fluc));
+ fprintf(stdout, "Set %3d filtered fluctuation: %12.6e\n", s+1, std::sqrt(fluc));
}
- fprintf(stdout, "Overall filtered fluctuation: %12.6e\n", sqrt(fluctot/nset));
+ fprintf(stdout, "Overall filtered fluctuation: %12.6e\n", std::sqrt(fluctot/nset));
fprintf(stdout, "\n");
sfree(filt);
int s, buflen = 0;
if (NULL != fn_fitted)
{
- buflen = strlen(fn_fitted)+32;
+ buflen = std::strlen(fn_fitted)+32;
snew(buf2, buflen);
- strncpy(buf2, fn_fitted, buflen);
- buf2[strlen(buf2)-4] = '\0';
+ std::strncpy(buf2, fn_fitted, buflen);
+ buf2[std::strlen(buf2)-4] = '\0';
}
for (s = 0; s < nset; s++)
{
static real tb = -1, te = -1, frac = 0.5, filtlen = 0, binwidth = 0.1, aver_start = 0;
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, nFitPoints = 100;
- static real temp = 298.15, fit_start = 1, fit_end = 60, diffusion = 5e-5, rcut = 0.35;
+ static gmx_bool bIntegrate = FALSE, bRegression = FALSE, bLuzar = FALSE;
+ static int nsets_in = 1, d = 1, nb_min = 4, resol = 10;
+ static real temp = 298.15, fit_start = 1, fit_end = 60;
/* must correspond to enum avbar* declared at beginning of file */
static const char *avbar_opt[avbarNR+1] = {
};
#define NPA asize(pa)
- FILE *out, *out_fit;
+ FILE *out;
int n, nlast, s, nset, i, j = 0;
real **val, *t, dt, tot, error;
double *av, *sig, cum1, cum2, cum3, cum4, db;
cum3 /= n;
cum4 /= n;
av[s] = cum1;
- sig[s] = sqrt(cum2);
+ sig[s] = std::sqrt(cum2);
if (n > 1)
{
- error = sqrt(cum2/(n-1));
+ error = std::sqrt(cum2/(n-1));
}
else
{
}
printf("SS%d %13.6e %12.6e %12.6e %6.3f %6.3f\n",
s+1, av[s], sig[s], error,
- sig[s] ? cum3/(sig[s]*sig[s]*sig[s]*sqrt(8/M_PI)) : 0,
+ sig[s] ? cum3/(sig[s]*sig[s]*sig[s]*std::sqrt(8/M_PI)) : 0,
sig[s] ? cum4/(sig[s]*sig[s]*sig[s]*sig[s]*3)-1 : 0);
}
printf("\n");
{
out = xvgropen(msdfile, "Mean square displacement",
"time", "MSD (nm\\S2\\N)", oenv);
- nlast = (int)(n*frac);
+ nlast = static_cast<int>(n*frac);
for (s = 0; s < nset; s++)
{
for (j = 0; j <= nlast; j++)
{
tot += sqr(val[s][i]-val[s][i+j]);
}
- tot /= (real)(n-j);
+ tot /= (n-j);
fprintf(out, " %g %8g\n", dt*j, tot);
}
if (s < nset-1)
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
#include "gromacs/topology/index.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
static void dump_dih_trr(int nframes, int nangles, real **dih, const char *fn,
{
for (m = 0; (m < 2); m++)
{
- x[k][l] = (m == 0) ? cos(dih[j][i]) : sin(dih[j][i]);
+ // This is just because the compler and static-analyzer cannot
+ // know that dih[j][i] is always valid. Since it occurs in the innermost
+ // loop over angles and will only trigger on coding errors, we
+ // only enable it for debug builds.
+ GMX_ASSERT(dih != NULL && dih[j] != NULL, "Incorrect dihedral array data");
+ x[k][l] = (m == 0) ? std::cos(dih[j][i]) : std::sin(dih[j][i]);
l++;
if (l == DIM)
{
};
FILE *out;
- real tmp, dt;
- int status, isize;
+ real dt;
+ int isize;
atom_id *index;
char *grpname;
- real maxang, Jc, S2, norm_fac, maxstat;
+ real maxang, S2, norm_fac, maxstat;
unsigned long mode;
int nframes, maxangstat, mult, *angstat;
- int i, j, total, nangles, natoms, nat2, first, last, angind;
+ int i, j, nangles, first, last;
gmx_bool bAver, bRb, bPeriodic,
- bFrac, /* calculate fraction too? */
- bTrans, /* worry about transtions too? */
- bCorr; /* correlation function ? */
- real t, aa, aver, aver2, aversig, fraction; /* fraction trans dihedrals */
+ bFrac, /* calculate fraction too? */
+ bTrans, /* worry about transtions too? */
+ bCorr; /* correlation function ? */
+ real aver, aver2, aversig; /* fraction trans dihedrals */
double tfrac = 0;
char title[256];
real **dih = NULL; /* mega array with all dih. angles at all times*/
- char buf[80];
real *time, *trans_frac, *aver_angle;
t_filenm fnm[] = {
{ efTRX, "-f", NULL, ffREAD },
mult = 4;
maxang = 360.0;
bRb = FALSE;
+
+ GMX_RELEASE_ASSERT(opt[0] != NULL, "Internal option inconsistency; opt[0]==NULL after processing");
+
switch (opt[0][0])
{
case 'a':
}
/* Calculate bin size */
- maxangstat = (int)(maxang/binwidth+0.5);
+ maxangstat = static_cast<int>(maxang/binwidth+0.5);
binwidth = maxang/maxangstat;
rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
if (bPBC)
{
real dd = dih[j][i];
- fprintf(out, " %8.3f", atan2(sin(dd), cos(dd))*RAD2DEG);
+ fprintf(out, " %8.3f", std::atan2(std::sin(dd), std::cos(dd))*RAD2DEG);
}
else
{
aver += RAD2DEG*aver_angle[i];
aver2 += sqr(RAD2DEG*aver_angle[i]);
}
- aver /= (real) nframes;
- aver2 /= (real) nframes;
- aversig = sqrt(aver2-sqr(aver));
+ aver /= nframes;
+ aver2 /= nframes;
+ aversig = std::sqrt(aver2-sqr(aver));
printf("Found points in the range from %d to %d (max %d)\n",
first, last, maxangstat);
printf(" < angle > = %g\n", aver);
maxstat = 0;
for (i = first; (i <= last); i++)
{
- maxstat = max(maxstat, angstat[i]*norm_fac);
+ maxstat = std::max(maxstat, angstat[i]*norm_fac);
}
if (output_env_get_print_xvgr_codes(oenv))
{
*/
#include "gmxpre.h"
-#include <ctype.h>
-#include <float.h>
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cctype>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
+#include <limits>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/enxio.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/dir_separator.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/snprintf.h"
srenew(lc->names, lc->Nalloc);
}
snew(lc->names[lc->N], name_length+1);
- strncpy(lc->names[lc->N], name, name_length);
+ std::strncpy(lc->names[lc->N], name, name_length);
lc->N++;
}
{
return FALSE;
}
- len = strlen(lc->names[index]);
+ len = std::strlen(lc->names[index]);
if (len != name_length)
{
return FALSE;
}
- if (strncmp(lc->names[index], name, name_length) == 0)
+ if (std::strncmp(lc->names[index], name, name_length) == 0)
{
return TRUE;
}
for (i = 0; i < lc->N; i++)
{
- if (strncmp(lc->names[i], name, name_length) == 0)
+ if (std::strncmp(lc->names[i], name, name_length) == 0)
{
return i;
}
static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
{
int i;
- size_t np;
str[0] = 0; /* reset the string */
if (lv->dhdl < 0)
}
if (lv->lc->N > 1)
{
- str += sprintf(str, ")");
+ sprintf(str, ")");
}
}
else
{
/* this lambda vector describes a derivative */
str += sprintf(str, "dH/dl");
- if (strlen(lv->lc->names[lv->dhdl]) > 0)
+ if (std::strlen(lv->lc->names[lv->dhdl]) > 0)
{
sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
}
/* write a shortened version of the lambda vec to a preallocated string */
static void lambda_vec_print_short(const lambda_vec_t *lv, char *str)
{
- int i;
- size_t np;
-
if (lv->index >= 0)
{
sprintf(str, "%6d", lv->index);
static void lambda_vec_print_intermediate(const lambda_vec_t *a,
const lambda_vec_t *b, char *str)
{
- int i;
- size_t np;
-
str[0] = 0;
if ( (a->index >= 0) && (b->index >= 0) )
{
- sprintf(str, "%6.3f", ((double)a->index+(double)b->index)/2.);
+ sprintf(str, "%6.3f", (a->index+b->index)/2.0);
}
else
{
if ( (a->dhdl < 0) && (b->dhdl < 0) )
{
- sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.);
+ sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.0);
}
}
}
double df = a->val[i] - b->val[i];
ret += df*df;
}
- return sqrt(ret);
+ return std::sqrt(ret);
}
static gmx_bool lambda_vec_cmp_native(const lambda_vec_t *a,
const lambda_vec_t *b)
{
- int i;
-
if (a->lc != b->lc)
{
gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
/* check if there's room */
if ( (sc->nsamples + 1) > sc->nsamples_alloc)
{
- sc->nsamples_alloc = max(2*sc->nsamples_alloc, 2);
+ sc->nsamples_alloc = std::max(2*sc->nsamples_alloc, 2);
srenew(sc->s, sc->nsamples_alloc);
srenew(sc->r, sc->nsamples_alloc);
}
}
else
{
- *nbin = (xmax-(*xmin))/(*dx);
+ *nbin = static_cast<int>((xmax-(*xmin))/(*dx));
}
if (*nbin > *nbin_alloc)
/* calculate the bin corresponding to the middle of the
original bin */
double x = hdx*(j+0.5) + xmin_hist;
- int binnr = (int)((x-(*xmin))/(*dx));
+ int binnr = static_cast<int>((x-(*xmin))/(*dx));
if (binnr >= *nbin || binnr < 0)
{
int endi = sc->r[i].end;
for (j = starti; j < endi; j++)
{
- int binnr = (int)((sc->s[i]->du[j]-(*xmin))/(*dx));
+ int binnr = static_cast<int>((sc->s[i]->du[j]-(*xmin))/(*dx));
if (binnr >= *nbin || binnr < 0)
{
binnr = (*nbin)-1;
int nbin = 0;
int nbin_alloc = 0;
double dx = 0;
- double min = 0;
+ double minval = 0;
int i;
lambda_data_t *bl_head = sd->lb;
xvgr_new_dataset(fp, 0, 0, NULL, oenv);
}
- sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &min,
+ sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &minval,
nbin_default);
for (i = 0; i < nbin; i++)
{
- double xmin = i*dx + min;
- double xmax = (i+1)*dx + min;
+ double xmin = i*dx + minval;
+ double xmax = (i+1)*dx + minval;
fprintf(fp, "%g %d\n%g %d\n", xmin, hist[i], xmax, hist[i]);
}
lambda_data_t *bl;
int nlambda = 0;
barres_t *res;
- int i;
gmx_bool dhdl = FALSE;
gmx_bool first = TRUE;
lambda_data_t *bl_head = sd->lb;
Wfac = delta_lambda;
}
- disc_err = max(disc_err, Wfac*br->a->s[j]->hist->dx[0]);
+ disc_err = std::max(disc_err, Wfac*br->a->s[j]->hist->dx[0]);
}
}
for (j = 0; j < br->b->nsamples; j++)
{
Wfac = delta_lambda;
}
- disc_err = max(disc_err, Wfac*br->b->s[j]->hist->dx[0]);
+ disc_err = std::max(disc_err, Wfac*br->b->s[j]->hist->dx[0]);
}
}
}
double end_time;
if (s->start_time < begin_t)
{
- r->start = (int)((begin_t - s->start_time)/s->delta_time);
+ r->start = static_cast<int>((begin_t - s->start_time)/s->delta_time);
}
end_time = s->delta_time*s->ndu + s->start_time;
if (end_time > end_t)
{
- r->end = (int)((end_t - s->start_time)/s->delta_time);
+ r->end = static_cast<int>((end_t - s->start_time)/s->delta_time);
}
}
else
int i, int ni)
{
int j;
- int hist_start, hist_end;
gmx_int64_t ntot_start;
gmx_int64_t ntot_end;
/* now fix start and end fields */
/* the casts avoid possible overflows */
- ntot_start = (gmx_int64_t)(sc_orig->ntot*(double)i/(double)ni);
- ntot_end = (gmx_int64_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
+ ntot_start = static_cast<gmx_int64_t>(sc_orig->ntot*static_cast<double>(i)/static_cast<double>(ni));
+ ntot_end = static_cast<gmx_int64_t>(sc_orig->ntot*static_cast<double>(i+1)/static_cast<double>(ni));
ntot_so_far = 0;
for (j = 0; j < sc->nsamples; j++)
{
new_end = 0;
}
/* and write the new range */
- sc->r[j].start = (int)new_start;
- sc->r[j].end = (int)new_end;
+ GMX_RELEASE_ASSERT(new_start <= std::numeric_limits<int>::max(), "Value of 'new_start' too large for int converstion");
+ GMX_RELEASE_ASSERT(new_end <= std::numeric_limits<int>::max(), "Value of 'new_end' too large for int converstion");
+ sc->r[j].start = static_cast<int>(new_start);
+ sc->r[j].end = static_cast<int>(new_end);
}
else
{
/* first calculate normalized bounds
(where 0 is the start of the hist range, and 1 the end) */
- ntot_start_norm = (ntot_start-ntot_so_far)/(double)ntot_add;
- ntot_end_norm = (ntot_end-ntot_so_far)/(double)ntot_add;
+ ntot_start_norm = (ntot_start-ntot_so_far)/static_cast<double>(ntot_add);
+ ntot_end_norm = (ntot_end-ntot_so_far)/static_cast<double>(ntot_add);
/* now fix the boundaries */
- ntot_start_norm = min(1, max(0., ntot_start_norm));
- ntot_end_norm = max(0, min(1., ntot_end_norm));
+ ntot_start_norm = std::min(1.0, std::max(0.0, ntot_start_norm));
+ ntot_end_norm = std::max(0.0, std::min(1.0, ntot_end_norm));
/* and calculate the overlap */
overlap = ntot_end_norm - ntot_start_norm;
{
int i, j;
- *Wmin = FLT_MAX;
- *Wmax = -FLT_MAX;
+ *Wmin = std::numeric_limits<float>::max();
+ *Wmax = -std::numeric_limits<float>::max();
for (i = 0; i < sc->nsamples; i++)
{
{
for (j = r->start; j < r->end; j++)
{
- *Wmin = min(*Wmin, s->du[j]*Wfac);
- *Wmax = max(*Wmax, s->du[j]*Wfac);
+ *Wmin = std::min(*Wmin, s->du[j]*Wfac);
+ *Wmax = std::max(*Wmax, s->du[j]*Wfac);
}
}
else
for (j = s->hist->nbin[hd]-1; j >= 0; j--)
{
- *Wmin = min(*Wmin, Wfac*(s->hist->x0[hd])*dx);
- *Wmax = max(*Wmax, Wfac*(s->hist->x0[hd])*dx);
+ *Wmin = std::min(*Wmin, Wfac*(s->hist->x0[hd])*dx);
+ *Wmax = std::max(*Wmax, Wfac*(s->hist->x0[hd])*dx);
/* look for the highest value bin with values */
if (s->hist->bin[hd][j] > 0)
{
- *Wmin = min(*Wmin, Wfac*(j+s->hist->x0[hd]+1)*dx);
- *Wmax = max(*Wmax, Wfac*(j+s->hist->x0[hd]+1)*dx);
+ *Wmin = std::min(*Wmin, Wfac*(j+s->hist->x0[hd]+1)*dx);
+ *Wmax = std::max(*Wmax, Wfac*(j+s->hist->x0[hd]+1)*dx);
break;
}
}
for (i = 0; i < n; i++)
{
- sum += 1./(1. + exp(Wfac*W[i] + sbMmDG));
+ sum += 1./(1. + std::exp(Wfac*W[i] + sbMmDG));
}
return sum;
{
double sum = 0.;
int i;
- int max;
+ int maxbin;
/* normalization factor multiplied with bin width and
number of samples (we normalize through M): */
double normdx = 1.;
{
hd = 1;
}
- dx = hist->dx[hd];
- max = hist->nbin[hd]-1;
+ dx = hist->dx[hd];
+ maxbin = hist->nbin[hd]-1;
if (type == 1)
{
- max = hist->nbin[hd]; /* we also add whatever was out of range */
+ maxbin = hist->nbin[hd]; /* we also add whatever was out of range */
}
- for (i = 0; i < max; i++)
+ for (i = 0; i < maxbin; i++)
{
double x = Wfac*((i+hist->x0[hd])+0.5)*dx; /* bin middle */
double pxdx = hist->bin[0][i]*normdx; /* p(x)dx */
- sum += pxdx/(1. + exp(x + sbMmDG));
+ sum += pxdx/(1. + std::exp(x + sbMmDG));
}
return sum;
double temp, double tol, int type)
{
double kT, beta, M;
- double DG;
- int i, j;
+ int i;
double Wfac1, Wfac2, Wmin, Wmax;
double DG0, DG1, DG2, dDG1;
- double sum1, sum2;
double n1, n2; /* numbers of samples as doubles */
kT = BOLTZ*temp;
n1 = ca->ntot;
n2 = cb->ntot;
- M = log(n1/n2);
+ M = std::log(n1/n2);
/*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
if (ca->foreign_lambda->dhdl < 0)
sample_coll_min_max(ca, Wfac1, &Wmin1, &Wmax1);
sample_coll_min_max(cb, Wfac2, &Wmin2, &Wmax2);
- Wmin = min(Wmin1, Wmin2);
- Wmax = max(Wmax1, Wmax2);
+ Wmin = std::min(Wmin1, Wmin2);
+ Wmax = std::max(Wmax1, Wmax2);
}
DG0 = Wmin;
Wfac2 = -beta*delta_lambda;
}
- M = log(n1/n2);
+ M = std::log(n1/n2);
/* calculate average in both directions */
{
for (j = r->start; j < r->end; j++)
{
- sigmafact += 1./(2. + 2.*cosh((M + Wfac1*s->du[j] - dg)));
+ sigmafact += 1./(2. + 2.*std::cosh((M + Wfac1*s->du[j] - dg)));
}
}
else
double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
- sigmafact += pxdx/(2. + 2.*cosh((M + x - dg)));
+ sigmafact += pxdx/(2. + 2.*std::cosh((M + x - dg)));
}
}
}
{
for (j = r->start; j < r->end; j++)
{
- sigmafact += 1./(2. + 2.*cosh((M - Wfac2*s->du[j] - dg)));
+ sigmafact += 1./(2. + 2.*std::cosh((M - Wfac2*s->du[j] - dg)));
}
}
else
double x = Wfac2*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
- sigmafact += pxdx/(2. + 2.*cosh((M - x - dg)));
+ sigmafact += pxdx/(2. + 2.*std::cosh((M - x - dg)));
}
}
}
/* Eq. 10 from
Shirts, Bair, Hooker & Pande, Phys. Rev. Lett 91, 140601 (2003): */
- *stddev = sqrt(((1./sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
+ *stddev = std::sqrt(((1.0/sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
}
int npee, p;
double dg_sig2, sa_sig2, sb_sig2, stddev_sig2; /* intermediate variance values
for calculated quantities */
- int nsample1, nsample2;
double temp = br->a->temp;
- int i, j;
+ int i;
double dg_min, dg_max;
gmx_bool have_hist = FALSE;
dg_min = calc_bar_lowlevel(br->a, br->b, temp, tol, -1);
dg_max = calc_bar_lowlevel(br->a, br->b, temp, tol, 1);
- if (fabs(dg_max - dg_min) > GMX_REAL_EPS*10)
+ if (std::abs(dg_max - dg_min) > GMX_REAL_EPS*10)
{
/* the histogram range error is the biggest of the differences
between the best estimate and the extremes */
- br->dg_histrange_err = fabs(dg_max - dg_min);
+ br->dg_histrange_err = std::abs(dg_max - dg_min);
}
br->dg_disc_err = 0.;
for (i = 0; i < br->a->nsamples; i++)
{
if (br->a->s[i]->hist)
{
- br->dg_disc_err = max(br->dg_disc_err, br->a->s[i]->hist->dx[0]);
+ br->dg_disc_err = std::max(br->dg_disc_err, br->a->s[i]->hist->dx[0]);
}
}
for (i = 0; i < br->b->nsamples; i++)
{
if (br->b->s[i]->hist)
{
- br->dg_disc_err = max(br->dg_disc_err, br->b->s[i]->hist->dx[0]);
+ br->dg_disc_err = std::max(br->dg_disc_err, br->b->s[i]->hist->dx[0]);
}
}
}
dstddev2 /= npee;
stddev_sig2 += (dstddev2-dstddev*dstddev)/(npee-1);
}
- br->dg_err = sqrt(dg_sig2/(npee_max - npee_min + 1));
- br->sa_err = sqrt(sa_sig2/(npee_max - npee_min + 1));
- br->sb_err = sqrt(sb_sig2/(npee_max - npee_min + 1));
- br->dg_stddev_err = sqrt(stddev_sig2/(npee_max - npee_min + 1));
+ br->dg_err = std::sqrt(dg_sig2/(npee_max - npee_min + 1));
+ br->sa_err = std::sqrt(sa_sig2/(npee_max - npee_min + 1));
+ br->sb_err = std::sqrt(sb_sig2/(npee_max - npee_min + 1));
+ br->dg_stddev_err = std::sqrt(stddev_sig2/(npee_max - npee_min + 1));
}
}
svar += (s2 - s*s)/(nb - 1);
}
- return sqrt(svar/(nbmax + 1 - nbmin));
+ return std::sqrt(svar/(nbmax + 1 - nbmin));
}
/* first find the end of the name */
if (!name_end_found)
{
- if (isspace(*str) || (*str == '=') )
+ if (std::isspace(*str) || (*str == '=') )
{
name_end_found = TRUE;
}
}
else
{
- if (!( isspace(*str) || (*str == '=') ))
+ if (!( std::isspace(*str) || (*str == '=') ))
{
return str;
}
}
-
/* read a vector-notation description of a lambda vector */
static gmx_bool read_lambda_compvec(const char *str,
lambda_vec_t *lv,
{
if (!start_reached)
{
- if (isalnum(*str))
+ if (std::isalnum(*str))
{
vector = FALSE;
start_reached = TRUE;
vector = TRUE;
start_reached = TRUE;
}
- else if (!isspace(*str))
+ else if (!std::isspace(*str))
{
gmx_fatal(FARGS, "Error in lambda components in %s", fn);
}
{
if (val_start)
{
- if (isspace(*str) || *str == ')' || *str == ',' || *str == '\0')
+ if (std::isspace(*str) || *str == ')' || *str == ',' || *str == '\0')
{
/* end of value */
if (lv == NULL)
}
}
}
- else if (isalnum(*str))
+ else if (std::isalnum(*str))
{
val_start = str;
}
}
else
{
+ GMX_RELEASE_ASSERT(lc_in != NULL, "Internal inconsistency? lc_in==NULL");
if (n == lc_in->N)
{
return OK;
const char *legend,
lambda_vec_t *lam)
{
- double lambda = 0;
const char *ptr = NULL, *ptr2 = NULL;
gmx_bool ok = FALSE;
gmx_bool bdhdl = FALSE;
ptr2 = legend;
do
{
- ptr2 = strstr(ptr2, tostr);
+ ptr2 = std::strstr(ptr2, tostr);
if (ptr2 != NULL)
{
ptr = ptr2;
if (ptr)
{
- ptr += strlen(tostr)-1; /* and advance past that 'to' */
+ ptr += std::strlen(tostr)-1; /* and advance past that 'to' */
}
else
{
/* look for the = sign */
- ptr = strrchr(legend, '=');
+ ptr = std::strrchr(legend, '=');
if (!ptr)
{
/* otherwise look for the last space */
- ptr = strrchr(legend, ' ');
+ ptr = std::strrchr(legend, ' ');
}
}
- if (strstr(legend, "dH"))
+ if (std::strstr(legend, "dH"))
{
ok = TRUE;
bdhdl = TRUE;
}
- else if (strchr(legend, 'D') != NULL && strchr(legend, 'H') != NULL)
+ else if (std::strchr(legend, 'D') != NULL && std::strchr(legend, 'H') != NULL)
{
ok = TRUE;
bdhdl = FALSE;
}
- else /*if (strstr(legend, "pV"))*/
+ else /*if (std::strstr(legend, "pV"))*/
{
return FALSE;
}
{
int dhdl_index;
const char *end;
- char buf[STRLEN];
- ptr = strrchr(legend, '=');
+ ptr = std::strrchr(legend, '=');
end = ptr;
if (ptr)
{
gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
}
}
- while (!isspace(*ptr))
+ while (!std::isspace(*ptr))
{
ptr--;
if (ptr < legend)
}
}
ptr++;
- strncpy(buf, ptr, (end-ptr));
- buf[(end-ptr)] = '\0';
dhdl_index = lambda_components_find(lam->lc, ptr, (end-ptr));
if (dhdl_index < 0)
{
char buf[STRLEN];
- strncpy(buf, ptr, (end-ptr));
+ std::strncpy(buf, ptr, (end-ptr));
buf[(end-ptr)] = '\0';
gmx_fatal(FARGS,
"Did not find lambda component for '%s' in %s",
bFound = FALSE;
/* first check for a state string */
- ptr = strstr(subtitle, "state");
+ ptr = std::strstr(subtitle, "state");
if (ptr)
{
int index = -1;
ptr = find_value(ptr);
if (ptr)
{
- index = strtol(ptr, &end, 10);
+ index = std::strtol(ptr, &end, 10);
if (ptr == end)
{
gmx_fatal(FARGS, "Incomplete state data in %s", fn);
return FALSE;
}
/* now find the lambda vector component names */
- while (*ptr != '(' && !isalnum(*ptr))
+ while (*ptr != '(' && !std::isalnum(*ptr))
{
ptr++;
if (*ptr == '\0')
{
/* compatibility mode: check for lambda in other ways. */
/* plain text lambda string */
- ptr = strstr(subtitle, "lambda");
+ ptr = std::strstr(subtitle, "lambda");
if (ptr == NULL)
{
/* xmgrace formatted lambda string */
- ptr = strstr(subtitle, "\\xl\\f{}");
+ ptr = std::strstr(subtitle, "\\xl\\f{}");
}
if (ptr == NULL)
{
/* xmgr formatted lambda string */
- ptr = strstr(subtitle, "\\8l\\4");
+ ptr = std::strstr(subtitle, "\\8l\\4");
}
if (ptr != NULL)
{
- ptr = strstr(ptr, "=");
+ ptr = std::strstr(ptr, "=");
}
if (ptr != NULL)
{
return bFound;
}
-static void filename2lambda(const char *fn)
+static double filename2lambda(const char *fn)
{
double lambda;
const char *ptr, *digitptr;
}
/* save the last position of a digit between the last two
separators = in the last dirname */
- if (dirsep > 0 && isdigit(*ptr))
+ if (dirsep > 0 && std::isdigit(*ptr))
{
digitptr = ptr;
}
{
gmx_fatal(FARGS, "Malformed number in file path '%s'", fn);
}
+ return lambda;
}
static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
int np;
gmx_bool native_lambda_read = FALSE;
char buf[STRLEN];
- lambda_vec_t lv;
xvg_init(ba);
if (subtitle != NULL)
{
/* try to extract temperature */
- ptr = strstr(subtitle, "T =");
+ ptr = std::strstr(subtitle, "T =");
if (ptr != NULL)
{
ptr += 3;
{
if (!native_lambda_read)
{
- /* Deduce lambda from the file name */
+ // Deduce lambda from the file name.
+ // Updated the routine filename2lambda to actually return lambda
+ // to avoid C++ warnings, but this usage does not seem to need it?
+ // EL 2015-07-10
filename2lambda(fn);
native_lambda_read = TRUE;
}
xvg_t *barsim;
samples_t *s;
int i;
- double *lambda;
snew(barsim, 1);
double *last_t, const char *filename)
{
int i, j;
- gmx_bool allocated;
- double old_foreign_lambda;
lambda_vec_t *foreign_lambda;
int type;
samples_t *s; /* convenience pointer */
int i, j;
samples_t *s;
int nhist;
- double old_foreign_lambda;
lambda_vec_t *foreign_lambda;
int type;
int nbins[2];
snew(foreign_lambda, 1);
lambda_vec_init(foreign_lambda, native_lambda->lc);
lambda_vec_copy(foreign_lambda, native_lambda);
- type = (int)(blk->sub[1].lval[1]);
+ type = static_cast<int>(blk->sub[1].lval[1]);
if (type == dhbtDH)
{
double old_foreign_lambda;
for (i = 0; i < nhist; i++)
{
- int nbin;
gmx_int64_t sum = 0;
for (j = 0; j < s->hist->nbin[i]; j++)
{
- int binv = (int)(blk->sub[i+2].ival[j]);
+ int binv = static_cast<int>(blk->sub[i+2].ival[j]);
s->hist->bin[i][j] = binv;
sum += binv;
int *npts = NULL; /* array to keep count & print at end */
lambda_vec_t **lambdas = NULL; /* array to keep count & print at end */
lambda_vec_t *native_lambda;
- double end_time; /* the end time of the last batch of samples */
int nsamples = 0;
lambda_vec_t start_lambda;
}
/* read the data from the DHCOLL block */
- rtemp = 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];
+ rtemp = 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];
old_start_lambda = fr->block[i].sub[0].dval[3];
delta_lambda = fr->block[i].sub[0].dval[4];
{
/* check the components */
lambda_components_check(&(sd->lc), j, name,
- strlen(name));
+ std::strlen(name));
}
else
{
lambda_components_add(&(sd->lc), name,
- strlen(name));
+ std::strlen(name));
}
}
lambda_vec_init(&start_lambda, &(sd->lc));
gmx_fatal(FARGS, "Can't handle both raw delta U data and histograms in the same file %s", fn);
}
- if (nsamples > 0)
+ if (nsamples == 0)
{
+ /* this is the first round; allocate the associated data
+ structures */
+ /*native_lambda=start_lambda;*/
+ lambda_vec_init(native_lambda, &(sd->lc));
+ lambda_vec_copy(native_lambda, &start_lambda);
+ nsamples = nblocks_raw+nblocks_hist;
+ snew(nhists, nsamples);
+ snew(npts, nsamples);
+ snew(lambdas, nsamples);
+ snew(samples_rawdh, nsamples);
+ for (i = 0; i < nsamples; i++)
+ {
+ nhists[i] = 0;
+ npts[i] = 0;
+ lambdas[i] = NULL;
+ samples_rawdh[i] = NULL; /* init to NULL so we know which
+ ones contain values */
+ }
+ }
+ else
+ {
+ // nsamples > 0 means this is NOT the first iteration
+
/* check the native lambda */
if (!lambda_vec_same(&start_lambda, native_lambda) )
{
}
/* check whether last iterations's end time matches with
the currrent start time */
- if ( (fabs(last_t - start_time) > 2*delta_time) && last_t >= 0)
+ if ( (std::abs(last_t - start_time) > 2*delta_time) && last_t >= 0)
{
/* it didn't. We need to store our samples and reallocate */
+
for (i = 0; i < nsamples; i++)
{
+ // nsamples is always >0 here, so samples_rawdh must be a valid pointer. Unfortunately
+ // cppcheck does not understand the logic unless the assert is inside the loop, but
+ // this is not performance-sensitive code.
+ GMX_RELEASE_ASSERT(samples_rawdh != NULL, "samples_rawdh==NULL with nsamples>0");
if (samples_rawdh[i])
{
/* insert it into the existing list */
}
}
}
- else
- {
- /* this is the first round; allocate the associated data
- structures */
- /*native_lambda=start_lambda;*/
- lambda_vec_init(native_lambda, &(sd->lc));
- lambda_vec_copy(native_lambda, &start_lambda);
- nsamples = nblocks_raw+nblocks_hist;
- snew(nhists, nsamples);
- snew(npts, nsamples);
- snew(lambdas, nsamples);
- snew(samples_rawdh, nsamples);
- for (i = 0; i < nsamples; i++)
- {
- nhists[i] = 0;
- npts[i] = 0;
- lambdas[i] = NULL;
- samples_rawdh[i] = NULL; /* init to NULL so we know which
- ones contain values */
- }
- }
/* and read them */
k = 0; /* counter for the lambdas, etc. arrays */
}
else if (fr->block[i].id == enxDHHIST)
{
- int type = (int)(fr->block[i].sub[1].lval[1]);
+ int type = static_cast<int>(fr->block[i].sub[1].lval[1]);
if (type == dhbtDH || type == dhbtDHDL)
{
int j;
int nd = 2, nbmin = 5, nbmax = 5;
int nbin = 100;
gmx_bool use_dhdl = FALSE;
- gmx_bool calc_s, calc_v;
- t_pargs pa[] = {
+ t_pargs pa[] = {
{ "-b", FALSE, etREAL, {&begin}, "Begin time for BAR" },
{ "-e", FALSE, etREAL, {&end}, "End time for BAR" },
{ "-temp", FALSE, etREAL, {&temp}, "Temperature (K)" },
};
#define NFILE asize(fnm)
- int f, i, j;
+ int f;
int nf = 0; /* file counter */
- int nbs;
int nfile_tot; /* total number of input files */
int nxvgfile = 0;
int nedrfile = 0;
int nresults; /* number of results in results array */
double *partsum;
- double prec, dg_tot, dg, sig, dg_tot_max, dg_tot_min;
+ double prec, dg_tot;
FILE *fpb, *fpi;
char dgformat[20], xvg2format[STRLEN], xvg3format[STRLEN];
char buf[STRLEN], buf2[STRLEN];
char ktformat[STRLEN], sktformat[STRLEN];
char kteformat[STRLEN], skteformat[STRLEN];
output_env_t oenv;
- double kT, beta;
+ double kT;
gmx_bool result_OK = TRUE, bEE = TRUE;
gmx_bool disc_err = FALSE;
{
gmx_fatal(FARGS, "Can not have negative number of digits");
}
- prec = pow(10, -nd);
+ prec = std::pow(10.0, static_cast<double>(-nd));
snew(partsum, (nbmax+1)*(nbmax+1));
nf = 0;
if (sum_disc_err > prec)
{
prec = sum_disc_err;
- nd = ceil(-log10(prec));
+ nd = static_cast<int>(std::ceil(-std::log10(prec)));
printf("WARNING: setting the precision to %g because that is the minimum\n reasonable number, given the expected discretization error.\n", prec);
}
/* print results in kT */
kT = BOLTZ*temp;
- beta = 1/kT;
printf("\nTemperature: %g K\n", temp);
{
stat_err = bar_err(nbmin, nbmax, partsum)*kT;
printf(" +/- ");
- printf(dgformat, max(max(stat_err, sum_disc_err), sum_histrange_err));
+ printf(dgformat, std::max(std::max(stat_err, sum_disc_err), sum_histrange_err));
}
printf("\n");
if (disc_err)
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
snew(mtot, bun->n);
+ clear_rvec(axis[0]);
+ clear_rvec(axis[1]);
+
for (end = 0; end < bun->nend; end++)
{
for (i = 0; i < bun->n; i++)
{
t_trxframe frout;
static rvec *xout = NULL;
- int i, end;
+ int i;
if (xout == NULL)
{
{ "-z", FALSE, etBOOL, {&bZ},
"Use the [IT]z[it]-axis as reference instead of the average axis" }
};
- FILE *out, *flen, *fdist, *fz, *ftilt, *ftiltr, *ftiltl;
+ FILE *flen, *fdist, *fz, *ftilt, *ftiltr, *ftiltl;
FILE *fkink = NULL, *fkinkr = NULL, *fkinkl = NULL;
t_trxstatus *status;
t_trxstatus *fpdb;
t_trxframe fr;
t_atoms outatoms;
real t, comp;
- int natoms;
char *grpname[MAX_ENDS], title[256];
/* FIXME: The constness should not be cast away */
char *anm = (char *)"CA", *rnm = (char *)"GLY";
- int i, j, gnx[MAX_ENDS];
+ int i, gnx[MAX_ENDS];
atom_id *index[MAX_ENDS];
t_bundle bun;
gmx_bool bKink;
fprintf(ftilt, " %6g", RAD2DEG*acos(bun.dir[i][ZZ]));
comp = bun.mid[i][XX]*bun.dir[i][XX]+bun.mid[i][YY]*bun.dir[i][YY];
fprintf(ftiltr, " %6g", RAD2DEG*
- asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
+ std::asin(comp/std::sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
comp = bun.mid[i][YY]*bun.dir[i][XX]-bun.mid[i][XX]*bun.dir[i][YY];
fprintf(ftiltl, " %6g", RAD2DEG*
- asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
+ std::asin(comp/std::sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
if (bKink)
{
rvec_sub(bun.end[0][i], bun.end[2][i], va);
copy_rvec(bun.mid[i], vr);
vr[ZZ] = 0;
unitv_no_table(vr, vr);
- fprintf(fkinkr, " %6g", RAD2DEG*asin(iprod(vc, vr)));
+ fprintf(fkinkr, " %6g", RAD2DEG*std::asin(iprod(vc, vr)));
vl[XX] = vr[YY];
vl[YY] = -vr[XX];
vl[ZZ] = 0;
- fprintf(fkinkl, " %6g", RAD2DEG*asin(iprod(vc, vl)));
+ fprintf(fkinkl, " %6g", RAD2DEG*std::asin(iprod(vc, vl)));
}
}
fprintf(flen, "\n");
*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
+#include <cmath>
+#include <cstdio>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
#define NPP asize(map)
int x, y;
-#define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
+#define INDEX(ppp) (((static_cast<int> (360+ppp*RAD2DEG)) % 360)/6)
x = INDEX(phi);
y = INDEX(psi);
#undef INDEX
- return (gmx_bool) map[x][y];
+
+ return (map[x][y] == '1') ? TRUE : FALSE;
}
atom_id *make_chi_ind(int nl, t_dlist dl[], int *ndih)
{
mult = 3;
- return (int) (chi*mult/360.0);
+ return static_cast<int>(chi*mult/360.0);
}
snew(data, nf);
if (bRAD)
{
- strcpy(ystr, "Angle (rad)");
+ std::strcpy(ystr, "Angle (rad)");
}
else
{
- strcpy(ystr, "Angle (degrees)");
+ std::strcpy(ystr, "Angle (degrees)");
}
/* Dump em all */
real *normhisto;
real **Jc, **Jcsig;
int ****his_aa_ss = NULL;
- int ***his_aa, **his_aa1, *histmp;
+ int ***his_aa, *histmp;
int i, j, k, m, n, nn, Dih, nres, hindex, angle;
gmx_bool bBfac, bOccup;
char hisfile[256], hhisfile[256], sshisfile[256], title[256], *ss_str = NULL;
}
if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac)))
{
- hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
+ hindex = static_cast<int>(((dih[j][0]+M_PI)*nbin)/(2*M_PI));
range_check(hindex, 0, nbin);
/* Assign dihedral to either of the structure determined
sprintf(title, "\\xc\\f{}\\s%d\\N Distribution for %s",
Dih-NONCHI+1, residue_name);
}
- strcpy(hhisfile, hisfile);
- strcat(hhisfile, ".xvg");
+ std::strcpy(hhisfile, hisfile);
+ std::strcat(hhisfile, ".xvg");
fp = xvgropen(hhisfile, title, "Degrees", "", oenv);
if (output_env_get_print_xvgr_codes(oenv))
{
fprintf(fp, "%10g %10g\n", phi, psi);
if (bViol)
{
- fprintf(gp, "%d\n", !bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]));
+ fprintf(gp, "%d\n", (bAllowed(dih[Phi][j], RAD2DEG*dih[Psi][j]) == FALSE) );
}
if (bOm)
{
omega = RAD2DEG*dih[Om][j];
- mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
+ mat[static_cast<int>(((phi*NMAT)/360)+NMAT/2)][static_cast<int>(((psi*NMAT)/360)+NMAT/2)]
+= omega;
}
}
for (k = 0; (k < NMAT); k++)
{
mat[j][k] /= nf;
- lo = min(mat[j][k], lo);
- hi = max(mat[j][k], hi);
+ lo = std::min(mat[j][k], lo);
+ hi = std::max(mat[j][k], hi);
}
}
/* Symmetrise */
- if (fabs(lo) > fabs(hi))
+ if (std::abs(lo) > std::abs(hi))
{
hi = -lo;
}
{
/* based on order_params below */
FILE *fp;
- int nh[edMax];
int i, Dih, Xi;
/* must correspond with enum in pp2shift.h:38 */
oenv);
xvgr_legend(fp, NONCHI+maxchi, (const char**)leg, oenv);
- for (Dih = 0; (Dih < edMax); Dih++)
- {
- nh[Dih] = 0;
- }
-
fprintf(fp, "%5s ", "#Res.");
fprintf(fp, "%10s %10s %10s ", leg[edPhi], leg[edPsi], leg[edOmega]);
for (Xi = 0; Xi < maxchi; Xi++)
x0 = y0 = z0 = 1000.0;
for (i = 0; (i < atoms->nr); i++)
{
- x0 = min(x0, x[i][XX]);
- y0 = min(y0, x[i][YY]);
- z0 = min(z0, x[i][ZZ]);
+ x0 = std::min(x0, x[i][XX]);
+ y0 = std::min(y0, x[i][YY]);
+ z0 = std::min(z0, x[i][ZZ]);
}
x0 *= 10.0; /* nm -> angstrom */
y0 *= 10.0; /* nm -> angstrom */
atom_id isize, *index;
int ndih, nactdih, nf;
real **dih, *trans_frac, *aver_angle, *time;
- int i, j, **chi_lookup, *multiplicity;
+ int i, **chi_lookup, *multiplicity;
t_filenm fnm[] = {
{ efSTX, "-s", NULL, ffREAD },
snew(multiplicity, ndih);
mk_multiplicity_lookup(multiplicity, maxchi, nlist, dlist, ndih);
- strcpy(grpname, "All residues, ");
+ std::strcpy(grpname, "All residues, ");
if (bPhi)
{
- strcat(grpname, "Phi ");
+ std::strcat(grpname, "Phi ");
}
if (bPsi)
{
- strcat(grpname, "Psi ");
+ std::strcat(grpname, "Psi ");
}
if (bOmega)
{
- strcat(grpname, "Omega ");
+ std::strcat(grpname, "Omega ");
}
if (bChi)
{
- strcat(grpname, "Chi 1-");
- sprintf(grpname + strlen(grpname), "%i", maxchi);
+ std::strcat(grpname, "Chi 1-");
+ sprintf(grpname + std::strlen(grpname), "%i", maxchi);
}
*/
#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/fileio/confio.h"
/* Generate new swapping candidates */
do
{
- iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
- jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
+ iswap = static_cast<int>(1+(nn-2)*gmx_rng_uniform_real(rng));
+ jswap = static_cast<int>(1+(nn-2)*gmx_rng_uniform_real(rng));
}
while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
else if (kT > 0)
{
/* Try Monte Carlo step */
- prob = exp(-(enext-ecur)/(enorm*kT));
+ prob = std::exp(-(enext-ecur)/(enorm*kT));
}
if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
}
r2 /= (isize*(isize-1))/2;
- return sqrt(r2);
+ return std::sqrt(r2);
}
static int rms_dist_comp(const void *a, const void *b)
t_dist *row;
t_clustid *c;
int **nnb;
- int i, j, k, cid, diff, max;
+ int i, j, k, cid, diff, maxval;
gmx_bool bChange;
real **mcpy = NULL;
else
{
/* Put all neighbors nearer than rmsdcut in the list */
- max = 0;
- k = 0;
+ maxval = 0;
+ k = 0;
for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
{
if (row[j].j != i)
{
- if (k >= max)
+ if (k >= maxval)
{
- max += 10;
- srenew(nnb[i], max);
+ maxval += 10;
+ srenew(nnb[i], maxval);
}
nnb[i][k] = row[j].j;
k++;
}
}
- if (k == max)
+ if (k == maxval)
{
- srenew(nnb[i], max+1);
+ srenew(nnb[i], maxval+1);
}
nnb[i][k] = -1;
}
}
}
-/* Again, I don't see the point in this... (AF) */
-/* for(i=0; (i<n1); i++) { */
-/* for(j=0; (j<n1); j++) */
-/* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
-/* } */
-/* for(i=0; (i<n1); i++) { */
-/* for(j=0; (j<n1); j++) */
-/* mat[i][j] = mcpy[i][j]; */
-/* } */
done_matrix(n1, &mcpy);
sfree(c);
{
t_dist *row;
t_nnb *nnb;
- int i, j, k, j1, max;
+ int i, j, k, j1, maxval;
/* Put all neighbors nearer than rmsdcut in the list */
fprintf(stderr, "Making list of neighbors within cutoff ");
snew(row, n1);
for (i = 0; (i < n1); i++)
{
- max = 0;
- k = 0;
+ maxval = 0;
+ k = 0;
/* put all neighbors within cut-off in list */
for (j = 0; j < n1; j++)
{
if (mat[i][j] < rmsdcut)
{
- if (k >= max)
+ if (k >= maxval)
{
- max += 10;
- srenew(nnb[i].nb, max);
+ maxval += 10;
+ srenew(nnb[i].nb, maxval);
}
nnb[i].nb[k] = j;
k++;
static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
{
- int i, j, v;
+ int i, j;
for (i = 0; i < nf; i++)
{
const char *ext;
char buf[STRLEN];
- if (strchr(fn, '%'))
+ if (std::strchr(fn, '%'))
{
gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
}
/* number of digits needed in numbering */
- i = (int)(log(maxnr)/log(10)) + 1;
+ i = static_cast<int>((std::log(static_cast<real>(maxnr))/std::log(10.0)) + 1);
/* split fn and ext */
- ext = strrchr(fn, '.');
+ ext = std::strrchr(fn, '.');
if (!ext)
{
gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
ext++;
/* insert e.g. '%03d' between fn and ext */
sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
- snew(fnout, strlen(buf)+1);
- strcpy(fnout, buf);
+ snew(fnout, std::strlen(buf)+1);
+ std::strcpy(fnout, buf);
return fnout;
}
ntrans[clust->cl[i-1]-1]++;
ntrans[clust->cl[i]-1]++;
trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
- maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
+ maxtrans = static_cast<int>(std::max(static_cast<real>(maxtrans), trans[clust->cl[i]-1][clust->cl[i-1]-1]));
}
}
ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
if (transfn)
{
fp = gmx_ffopen(transfn, "w");
- i = min(maxtrans+1, 80);
+ i = std::min(maxtrans+1, 80);
write_xpm(fp, 0, "Cluster Transitions", "# transitions",
"from cluster", "to cluster",
clust->ncl, clust->ncl, axis, axis, trans,
/* do we write all structures? */
if (write_ncl)
{
- trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
+ trxsfn = parse_filename(trxfn, std::max(write_ncl, clust->ncl));
snew(bWrite, nf);
}
ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
{
do_fit(natom, mass, xtps, xav);
}
- r = cl;
write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
}
}
for (j = i; j < mat->nx; j++)
{
rms->sumrms += rms->mat[i][j];
- rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
+ rms->maxrms = std::max(rms->maxrms, rms->mat[i][j]);
if (j != i)
{
- rms->minrms = min(rms->minrms, rms->mat[i][j]);
+ rms->minrms = std::min(rms->minrms, rms->mat[i][j]);
}
}
}
real *eigenvectors;
int isize = 0, ifsize = 0, iosize = 0;
- atom_id *index = NULL, *fitidx, *outidx;
+ atom_id *index = NULL, *fitidx = NULL, *outidx = NULL;
char *grpname;
real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
char buf[STRLEN], buf1[80], title[STRLEN];
else /* !bReadMat */
{
rms = init_mat(nf, method == m_diagonalize);
- nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
+ nrms = (static_cast<gmx_int64_t>(nf)*static_cast<gmx_int64_t>(nf-1))/2;
if (!bRMSdist)
{
fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
rmsd = rmsdev(isize, mass, xx[i2], x1);
set_mat_entry(rms, i1, i2, rmsd);
}
- nrms -= (gmx_int64_t) (nf-i1-1);
- fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
+ nrms -= nf-i1-1;
+ fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 " ", nrms);
}
sfree(x1);
}
calc_dist(isize, xx[i2], d2);
set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
}
- nrms -= (nf-i1-1);
- fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
+ nrms -= nf-i1-1;
+ fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 " ", nrms);
}
/* Clean up work arrays */
for (i = 0; (i < isize); i++)
/* Do a diagonalization */
snew(eigenvalues, nf);
snew(eigenvectors, nf*nf);
- memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
+ std::memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
sfree(eigenvectors);
{
useatoms.atomname[i] = top.atoms.atomname[index[i]];
useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
- useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
+ useatoms.nres = std::max(useatoms.nres, useatoms.atom[i].resind+1);
copy_rvec(xtps[index[i]], usextps[i]);
}
useatoms.nr = isize;
{
write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
- rlo_top, rhi_top, 0.0, (real) ncluster,
+ rlo_top, rhi_top, 0.0, ncluster,
&ncluster, TRUE, rlo_bot, rhi_bot);
}
else
*/
#include "gmxpre.h"
-#include <math.h>
+#include <cmath>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/matio.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
static void clust_size(const char *ndx, const char *trx, const char *xpm,
t_block *mols = NULL;
gmx_mtop_atomlookup_t alook;
t_atom *atom;
- int version, generation, ii, jj, nsame;
+ int version, generation, ii, jj;
real temp, tfac;
/* Cluster size distribution (matrix) */
real **cs_dist = NULL;
real tf, dx2, cut2, *t_x = NULL, *t_y, cmid, cmax, cav, ekin;
- int i, j, k, ai, aj, ak, ci, cj, nframe, nclust, n_x, n_y, max_size = 0;
+ int i, j, k, ai, aj, ci, cj, nframe, nclust, n_x, max_size = 0;
int *clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto;
t_rgb rlo = { 1.0, 1.0, 1.0 };
printf("Using molecules rather than atoms. Not reading index file %s\n",
ndx);
}
+ GMX_RELEASE_ASSERT(mtop != NULL, "Trying to access mtop->mols from NULL mtop pointer");
mols = &(mtop->mols);
/* Make dummy index */
/* Compute distance */
if (bMol)
{
+ GMX_RELEASE_ASSERT(mols != NULL, "Cannot access index[] from NULL mols pointer");
bSame = FALSE;
for (ii = mols->index[ai]; !bSame && (ii < mols->index[ai+1]); ii++)
{
{
nclust++;
cs_dist[n_x-1][ci-1] += 1.0;
- max_size = max(max_size, ci);
+ max_size = std::max(max_size, ci);
if (ci > 1)
{
cav += ci;
{
if (bMol)
{
+ GMX_RELEASE_ASSERT(mols != NULL, "Cannot access index[] from NULL mols pointer");
for (j = mols->index[i]; (j < mols->index[i+1]); j++)
{
fprintf(fp, "%d\n", j+1);
nelem += cs_dist[i][j];
}
fprintf(fp, "%5d %8.3f\n", j+1, nelem/n_x);
- nhisto += (int)((j+1)*nelem/n_x);
+ nhisto += static_cast<int>((j+1)*nelem/n_x);
}
fprintf(fp, "%5d %8.3f\n", j+1, 0.0);
xvgrclose(fp);
{
cmid = cs_dist[i][j];
}
- cmax = max(cs_dist[i][j], cmax);
+ cmax = std::max(cs_dist[i][j], cmax);
}
}
fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
{
cmid = cs_dist[i][j];
}
- cmax = max(cs_dist[i][j], cmax);
+ cmax = std::max(cs_dist[i][j], cmax);
}
}
fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
};
#define NPA asize(pa)
const char *fnNDX, *fnTPR;
- gmx_bool bSQ, bRDF;
t_rgb rgblo, rgbhi;
t_filenm fnm[] = {
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
{
fprintf(debug, " %s-%s", s1, s2);
}
- return strcmp(s1, s2);
+ return std::strcmp(s1, s2);
}
int find_next_match_atoms_in_res(int *i1, atom_id index1[],
int dx, dy, dmax, cmp;
gmx_bool bFW = FALSE;
- dx = dy = 0;
cmp = NOTSET;
- dmax = max(m1-*i1, m2-*i2);
- for (dx = 0; dx < dmax && cmp != 0; dx++)
+ dmax = std::max(m1-*i1, m2-*i2);
+ for (dx = 0, dy = 0; dx < dmax && cmp != 0; dx++)
{
for (dy = dx; dy < dmax && cmp != 0; dy++)
{
int dx, dy, dmax, cmp, rr1, rr2;
gmx_bool bFW = FALSE, bFF = FALSE;
- dx = dy = 0;
rr1 = 0;
while (rr1 < isize1 && *rnr1 != index1[rr1])
{
}
cmp = NOTSET;
- dmax = max(isize1-rr1, isize2-rr2);
+ dmax = std::max(isize1-rr1, isize2-rr2);
if (debug)
{
fprintf(debug, " R:%d-%d:%d-%d:%d ",
void find_matching_names(int *isize1, atom_id index1[], t_atoms *atoms1,
int *isize2, atom_id index2[], t_atoms *atoms2)
{
- int i, i1, i2, ii1, ii2, m1, m2;
+ int i1, i2, ii1, ii2, m1, m2;
int atcmp, rescmp;
- int r, rnr1, rnr2, prnr1, prnr2;
+ int rnr1, rnr2, prnr1, prnr2;
int rsize1, rsize2;
int *rindex1, *rindex2;
- char *resnm1, *resnm2, *atnm1, *atnm2;
char ***atnms1, ***atnms2;
t_resinfo *resinfo1, *resinfo2;
fprintf(debug, "R: %s%d %s%d\n",
*resinfo1[rnr1].name, rnr1, *resinfo2[rnr2].name, rnr2);
}
- rescmp = strcmp(*resinfo1[rnr1].name, *resinfo2[rnr2].name);
+ rescmp = std::strcmp(*resinfo1[rnr1].name, *resinfo2[rnr2].name);
}
if (debug)
{
output_env_t oenv;
/* counters */
- int i, j, m;
+ int i, m;
/* center of mass calculation */
- real tmas1, tmas2;
rvec xcm1, xcm2;
/* variables for fit */
{
name1 = *atoms1->atomname[index1[i]];
name2 = *atoms2->atomname[index2[i]];
- if (strcmp(name1, name2))
+ if (std::strcmp(name1, name2))
{
if (warn < 20)
{
rms += msd*mass;
msds[at] += msd;
}
- maxmsd = max(maxmsd, msds[at]);
- minmsd = min(minmsd, msds[at]);
+ maxmsd = std::max(maxmsd, msds[at]);
+ minmsd = std::min(minmsd, msds[at]);
totmass += mass;
}
- rms = sqrt(rms/totmass);
+ rms = std::sqrt(rms/totmass);
printf("Root mean square deviation after lsq fit = %g nm\n", rms);
if (bBfac)
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
};
FILE *out = NULL; /* initialization makes all compilers happy */
t_trxstatus *status;
- t_trxstatus *trjout;
t_topology top;
int ePBC;
t_atoms *atoms;
real t, tstart, tend, **mat2;
real xj, *w_rls = NULL;
real min, max, *axis;
- int ntopatoms, step;
- int natoms, nat, count, nframes0, nframes, nlevels;
+ int natoms, nat, nframes0, nframes, nlevels;
gmx_int64_t ndim, i, j, k, l;
int WriteXref;
const char *fitfile, *trxfile, *ndxfile;
const char *eigvalfile, *eigvecfile, *averfile, *logfile;
const char *asciifile, *xpmfile, *xpmafile;
- char str[STRLEN], *fitname, *ananame, *pcwd;
+ char str[STRLEN], *fitname, *ananame;
int d, dj, nfit;
atom_id *index, *ifit;
gmx_bool bDiffMass1, bDiffMass2;
{
if (bM)
{
- sqrtm[i] = sqrt(atoms->atom[index[i]].m);
+ sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
if (i)
{
bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i-1]);
if (bFit && bDiffMass1 && !bDiffMass2)
{
bDiffMass1 = natoms != nfit;
- i = 0;
for (i = 0; (i < natoms) && !bDiffMass1; i++)
{
bDiffMass1 = index[i] != ifit[i];
snew(x, natoms);
snew(xav, natoms);
ndim = natoms*DIM;
- if (sqrt(GMX_INT64_MAX) < ndim)
+ if (std::sqrt(static_cast<real>(GMX_INT64_MAX)) < static_cast<real>(ndim))
{
gmx_fatal(FARGS, "Number of degrees of freedoms to large for matrix.\n");
}
atoms, xread, NULL, epbcNONE, zerobox, natoms, index);
sfree(xread);
- fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", (int)ndim, (int)ndim);
+ fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", static_cast<int>(ndim), static_cast<int>(ndim));
nframes = 0;
nat = read_first_x(oenv, &status, trxfile, &t, &xread, box);
tstart = t;
snew(eigenvalues, ndim);
snew(eigenvectors, ndim*ndim);
- memcpy(eigenvectors, mat, ndim*ndim*sizeof(real));
+ std::memcpy(eigenvectors, mat, ndim*ndim*sizeof(real));
fprintf(stderr, "\nDiagonalizing ...\n");
fflush(stderr);
eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
}
fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n",
sum, bM ? "u " : "");
- if (fabs(trace-sum) > 0.01*trace)
+ if (std::abs(trace-sum) > 0.01*trace)
{
fprintf(stderr, "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
}
end = nframes-1;
fprintf(stderr, "\nWARNING: there are fewer frames in your trajectory than there are\n");
fprintf(stderr, "degrees of freedom in your system. Only generating the first\n");
- fprintf(stderr, "%d out of %d eigenvectors and eigenvalues.\n", end, (int)ndim);
+ fprintf(stderr, "%d out of %d eigenvectors and eigenvalues.\n", end, static_cast<int>(ndim));
}
else
{
"Eigenvector index", str, oenv);
for (i = 0; (i < end); i++)
{
- fprintf (out, "%10d %g\n", (int)i+1, eigenvalues[ndim-1-i]);
+ fprintf (out, "%10d %g\n", static_cast<int>(i+1), eigenvalues[ndim-1-i]);
}
xvgrclose(out);
{
fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
}
- fprintf(out, "Diagonalized the %dx%d covariance matrix\n", (int)ndim, (int)ndim);
+ fprintf(out, "Diagonalized the %dx%d covariance matrix\n", static_cast<int>(ndim), static_cast<int>(ndim));
fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n",
trace);
fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n",
sum);
- fprintf(out, "Wrote %d eigenvalues to %s\n", (int)end, eigvalfile);
+ fprintf(out, "Wrote %d eigenvalues to %s\n", static_cast<int>(end), eigvalfile);
if (WriteXref == eWXR_YES)
{
fprintf(out, "Wrote reference structure to %s\n", eigvecfile);
*/
#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);
*/
#include "gmxpre.h"
-#include <ctype.h>
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cctype>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
typedef struct {
t_electron *tmp1, *tmp2;
tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
- return strcmp(tmp1->atomname, tmp2->atomname);
+ return std::strcmp(tmp1->atomname, tmp2->atomname);
}
int get_electrons(t_electron **eltab, const char *fn)
if (!*nslices)
{
- *nslices = (int)(box[axis][axis] * 10); /* default value */
+ *nslices = static_cast<int>(box[axis][axis] * 10); /* default value */
fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
}
/* determine which slice atom is in */
if (bCenter)
{
- slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
+ slice = static_cast<int>(std::floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2);
}
else
{
- slice = (z / (*slWidth));
+ slice = static_cast<int>(z / (*slWidth));
}
sought.nr_el = 0;
sought.atomname = gmx_strdup(*(top->atoms.atomname[index[n][i]]));
double invvol;
int natoms; /* nr. atoms in trj */
t_trxstatus *status;
- int **slCount, /* nr. of atoms in one slice for a group */
- i, j, n, /* loop indices */
- ax1 = 0, ax2 = 0,
+ int i, n, /* loop indices */
nr_frames = 0, /* number of frames */
slice; /* current slice */
real t,
z;
real boxSz, aveBox;
- char *buf; /* for tmp. keeping atomname */
gmx_rmpbc_t gpbc = NULL;
if (axis < 0 || axis >= DIM)
if (!*nslices)
{
- *nslices = (int)(box[axis][axis] * 10); /* default value */
+ *nslices = static_cast<int>(box[axis][axis] * 10); /* default value */
fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
}
/* determine which slice atom is in */
if (bCenter)
{
- slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
+ slice = static_cast<int>(std::floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2);
}
else
{
- slice = floor(z / (*slWidth));
+ slice = static_cast<int>(std::floor(z / (*slWidth)));
}
/* Slice should already be 0<=slice<nslices, but we just make
xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
}
+ GMX_RELEASE_ASSERT(dens_opt[0] != NULL, "Option setting inconsistency; dens_opt[0] is NULL");
+
switch (dens_opt[0][0])
{
case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
}
else
{
- ncenter = 0;
+ ncenter = 0;
+ index_center = NULL;
}
fprintf(stderr, "\nSelect %d group%s to calculate density for:\n", ngrps, (ngrps > 1) ? "s" : "");
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
int gmx_densmap(int argc, char *argv[])
matrix box;
real t, m, mtot;
t_pbc pbc;
- int cav = 0, c1 = 0, c2 = 0, natoms;
+ int cav = 0, c1 = 0, c2 = 0;
char **grpname, title[256], buf[STRLEN];
const char *unit;
int i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
}
}
- if (strcmp(eunit[0], "nm-3") == 0)
+ GMX_RELEASE_ASSERT(eunit[0] != NULL, "Option setting inconsistency; eunit[0] is NULL");
+
+ if (std::strcmp(eunit[0], "nm-3") == 0)
{
nmpower = -3;
unit = "(nm^-3)";
}
- else if (strcmp(eunit[0], "nm-2") == 0)
+ else if (std::strcmp(eunit[0], "nm-2") == 0)
{
nmpower = -2;
unit = "(nm^-2)";
}
}
+ GMX_RELEASE_ASSERT(eaver[0] != NULL, "Option setting inconsistency; eaver[0] is NULL");
+
switch (eaver[0][0])
{
case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
}
- natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
+ read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
if (!bRadial)
{
if (n1 == 0)
{
- n1 = (int)(box[c1][c1]/bin + 0.5);
+ n1 = static_cast<int>(box[c1][c1]/bin + 0.5);
}
if (n2 == 0)
{
- n2 = (int)(box[c2][c2]/bin + 0.5);
+ n2 = static_cast<int>(box[c2][c2]/bin + 0.5);
}
}
else
{
- n1 = (int)(2*amax/bin + 0.5);
- nradial = (int)(rmax/bin + 0.5);
+ n1 = static_cast<int>(2*amax/bin + 0.5);
+ nradial = static_cast<int>(rmax/bin + 0.5);
/* cppcheck-suppress zerodiv fixed in 1.68-dev */
invspa = n1/(2*amax);
/* cppcheck-suppress zerodiv fixed in 1.68-dev */
{
m2 += 1;
}
- grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
+ grid[static_cast<int>(m1*n1)][static_cast<int>(m2*n2)] += invcellvol;
}
}
}
j = index[i];
pbc_dx(&pbc, x[j], center, dx);
axial = iprod(dx, direction);
- r = sqrt(norm2(dx) - axial*axial);
+ r = std::sqrt(norm2(dx) - axial*axial);
if (axial >= -amax && axial < amax && r < rmax)
{
if (bMirror)
{
r += rmax;
}
- grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
+ grid[static_cast<int>((axial + amax)*invspa)][static_cast<int>(r*invspz)] += 1;
}
}
}
{
if (!bXmax)
{
- sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
+ sprintf(buf+std::strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
}
else if (!bXmin)
{
- sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
+ sprintf(buf+std::strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
}
else
{
- sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
+ sprintf(buf+std::strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
}
}
if (ftp2bSet(efDAT, NFILE, fnm))
*/
#include "gmxpre.h"
-#include <ctype.h>
-#include <math.h>
-#include <string.h>
+#include <cctype>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/smalloc.h"
-#ifdef GMX_DOUBLE
-#define FLOOR(x) ((int) floor(x))
-#else
-#define FLOOR(x) ((int) floorf(x))
-#endif
enum {
methSEL, methBISECT, methFUNCFIT, methNR
}
- *zslices = 1+FLOOR(box[axis][axis]/bwz);
- *yslices = 1+FLOOR(box[ax2][ax2]/bw);
- *xslices = 1+FLOOR(box[ax1][ax1]/bw);
+ *zslices = 1+static_cast<int>(std::floor(box[axis][axis]/bwz));
+ *yslices = 1+static_cast<int>(std::floor(box[ax2][ax2]/bw));
+ *xslices = 1+static_cast<int>(std::floor(box[ax1][ax1]/bw));
if (bps1d)
{
if (*xslices < *yslices)
z -= box[axis][axis];
}
- slicex = ((int) (x/bbww[XX])) % *xslices;
- slicey = ((int) (y/bbww[YY])) % *yslices;
- slicez = ((int) (z/bbww[ZZ])) % *zslices;
+ slicex = static_cast<int>(x/bbww[XX]) % *xslices;
+ slicey = static_cast<int>(y/bbww[YY]) % *yslices;
+ slicez = static_cast<int>(z/bbww[ZZ]) % *zslices;
Densslice[slicex][slicey][slicez] += (top->atoms.atom[index[0][j]].m*dscale);
}
real std, var;
int i, j, n, order;
order = ftsize/2;
- std = ((real)order/2.0);
+ std = order/2.0;
var = std*std;
snew(kernel, ftsize);
gausskernel(kernel, ftsize, var);
*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/correlationfunctions/expfit.h"
#include "gromacs/correlationfunctions/integrate.h"
int calc_nbegin(int nx, real x[], real tbegin)
{
int nbegin;
- real dt, dtt;
/* Assume input x is sorted */
for (nbegin = 0; (nbegin < nx) && (x[nbegin] <= tbegin); nbegin++)
}
/* Take the one closest to tbegin */
- if (fabs(x[nbegin]-tbegin) > fabs(x[nbegin-1]-tbegin))
+ if (std::abs(x[nbegin]-tbegin) > std::abs(x[nbegin-1]-tbegin))
{
nbegin--;
}
}
else
{
- i0 = max(0, nbegin);
- i1 = min(nx-1, nbegin+nsmooth);
+ i0 = std::max(0, nbegin);
+ i1 = std::min(nx-1, nbegin+nsmooth);
printf("Making smooth transition from %d through %d\n", i0, i1);
for (i = 0; (i < i0); i++)
{
}
for (i = i0; (i <= i1); i++)
{
- fx = (i1-i)/(real)(i1-i0);
- fy = (i-i0)/(real)(i1-i0);
+ fx = static_cast<real>(i1-i)/(i1-i0);
+ fy = static_cast<real>(i-i0)/(i1-i0);
if (debug)
{
fprintf(debug, "x: %g factors for smoothing: %10g %10g\n", x[i], fx, fy);
FILE *fp, *cp;
t_complex *tmp, gw, hw, kw;
int i, nnx, nxsav;
- real fac, nu, dt, *ptr, maxeps, numax;
+ real fac, nu, dt, maxeps, numax;
gmx_fft_t fft;
int fftcode;
nx = read_xvg(opt2fn("-f", NFILE, fnm), &yd, &ny);
dt = yd[0][1] - yd[0][0];
- nxtail = min(tail/dt, nx);
+ nxtail = std::min(static_cast<int>(tail/dt), nx);
printf("Read data set containing %d colums and %d rows\n", ny, nx);
printf("Assuming (from data) that timestep is %g, nxtail = %d\n",
snew(y, 6);
for (i = 0; (i < ny); i++)
{
- snew(y[i], max(nx, nxtail));
+ snew(y[i], std::max(nx, nxtail));
}
for (i = 0; (i < nx); i++)
{
srenew(y, 3);
snew(y[2], nx);
- fac = 1.0/((real)nx);
+ fac = 1.0/nx;
for (i = 0; (i < nx); i++)
{
y[2][i] = fac;
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include <algorithm>
}
if (index < gb->nx)
{
- alpha = acos(cosa);
+ alpha = std::acos(cosa);
if (gb->bPhi)
{
cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
static void rvec2sprvec(rvec dipcart, rvec dipsp)
{
clear_rvec(dipsp);
- dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
- dipsp[1] = atan2(dipcart[YY], dipcart[XX]); /* Theta */
- dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
+ dipsp[0] = std::sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
+ dipsp[1] = std::atan2(dipcart[YY], dipcart[XX]); /* Theta */
+ dipsp[2] = std::atan2(std::sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]); /* Phi */
}
qtot = 0;
for (j = j0; j < j1; j++)
{
- q = fabs(atom[j].q);
+ q = std::abs(atom[j].q);
qtot += q;
for (k = 0; k < DIM; k++)
{
phi = dih_angle(xi, xj, xk, xl, &pbc,
r_ij, r_kj, r_kl, mm, nn, /* out */
&sign, &t1, &t2, &t3);
- cosa = cos(phi);
+ cosa = std::cos(phi);
}
else
{
* Multiply by 2 because we only take half the matrix of interactions
* into account.
*/
- fac = 2.0/((double) ngrp * (double) nframes);
+ fac = 2.0/(ngrp * nframes);
x0 = 0;
for (i = 0; i < last; i++)
qtot += atom[a].q;
}
/* This check is only for the count print */
- if (fabs(qtot) > 0.01)
+ if (std::abs(qtot) > 0.01)
{
ncharged++;
}
xdim += x[k][idim];
}
xdim /= (k1-k0);
- k = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
+ k = (static_cast<int>(xdim*nslice/box[idim][idim] + nslice)) % nslice;
rvec_inc(slab_dipole[k], mu);
}
ddc1 = ddc2 = ddc3 = 0;
for (i = k = 0; (i < n); i++)
{
- ddc1 += fabs(cos_angle(dip[i], xxx));
- ddc2 += fabs(cos_angle(dip[i], yyy));
- ddc3 += fabs(cos_angle(dip[i], zzz));
+ ddc1 += std::abs(cos_angle(dip[i], xxx));
+ ddc2 += std::abs(cos_angle(dip[i], yyy));
+ ddc3 += std::abs(cos_angle(dip[i], zzz));
if (bPairs)
{
for (j = i+1; (j < n); j++, k++)
{
dc = cos_angle(dip[i], dip[j]);
- d += fabs(dc);
+ d += std::abs(dc);
}
}
}
/* Determine the indexes of the energy grps we need */
for (i = 0; (i < nre); i++)
{
- if (strstr(enm[i].name, "Volume"))
+ if (std::strstr(enm[i].name, "Volume"))
{
iVol = i;
}
- else if (strstr(enm[i].name, "Mu-X"))
+ else if (std::strstr(enm[i].name, "Mu-X"))
{
iMu[XX] = i;
}
- else if (strstr(enm[i].name, "Mu-Y"))
+ else if (std::strstr(enm[i].name, "Mu-Y"))
{
iMu[YY] = i;
}
- else if (strstr(enm[i].name, "Mu-Z"))
+ else if (std::strstr(enm[i].name, "Mu-Z"))
{
iMu[ZZ] = i;
}
{
/* Use 0.7 iso 0.5 to account for pressure scaling */
/* rcut = 0.7*sqrt(max_cutoff2(box)); */
- rcut = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
+ rcut = 0.7*std::sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
}
M_av[m] += dipole[i][m]; /* M per frame */
mu_mol += dipole[i][m]*dipole[i][m]; /* calc. mu for distribution */
}
- mu_mol = sqrt(mu_mol);
+ mu_mol = std::sqrt(mu_mol);
mu_ave += mu_mol; /* calc. the average mu */
/* Update the dipole distribution */
- ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
+ ibin = static_cast<int>(ndipbin*mu_mol/mu_max + 0.5);
if (ibin < ndipbin)
{
dipole_bin[ibin]++;
if (cosaver)
{
compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
- rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
- sqr(dipaxis[YY]-0.5)+
- sqr(dipaxis[ZZ]-0.5));
+ rms_cos = std::sqrt(sqr(dipaxis[XX]-0.5)+
+ sqr(dipaxis[YY]-0.5)+
+ sqr(dipaxis[ZZ]-0.5));
if (bPairs)
{
fprintf(caver, "%10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
{
fprintf(outmtot, "%10g %12.8e %12.8e %12.8e %12.8e\n",
t, M_av[XX], M_av[YY], M_av[ZZ],
- sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
+ std::sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
}
for (m = 0; (m < DIM); m++)
{
do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
teller, gnx_tot, muall, dt,
- mode, strcmp(corrtype, "molsep"));
+ mode, std::strcmp(corrtype, "molsep"));
}
}
}
for (i = 0; (i < ndipbin); i++)
{
fprintf(outdd, "%10g %10f\n",
- (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
+ (i*mu_max)/ndipbin, static_cast<real>(dipole_bin[i])/teller);
}
xvgrclose(outdd);
sfree(dipole_bin);
*/
#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/fileio/confio.h"
tpa = (t_toppop *)a;
tpb = (t_toppop *)b;
- return (1e7*(tpb->v-tpa->v));
+ return static_cast<int>(1e7*(tpb->v-tpa->v));
}
static void add5(int ndr, real viol)
{
t_iatom *forceatoms;
int i, j, nat, n, type, nviol, ndr, label;
- real ener, rt, mviol, tviol, viol, lam, dvdl, drt;
+ real rt, mviol, tviol, viol, lam, dvdl, drt;
rvec *fshift;
static gmx_bool bFirst = TRUE;
gmx_fatal(FARGS, "ndr = %d, rt_6 = %f", ndr, fcd->disres.Rt_6[0]);
}
- rt = pow(fcd->disres.Rt_6[0], -1.0/6.0);
+ rt = std::pow(fcd->disres.Rt_6[0], static_cast<real>(-1.0/6.0));
dr[clust_id].aver1[ndr] += rt;
dr[clust_id].aver2[ndr] += sqr(rt);
- drt = pow(rt, -3.0);
+ drt = std::pow(rt, static_cast<real>(-3.0));
dr[clust_id].aver_3[ndr] += drt;
dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[0];
snew(fshift, SHIFTS);
- ener = interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
- (const rvec*)x, f, fshift,
- pbc, g, lam, &dvdl, NULL, fcd, NULL);
+ interaction_function[F_DISRES].ifunc(n, &forceatoms[i], forceparams,
+ (const rvec*)x, f, fshift,
+ pbc, g, lam, &dvdl, NULL, fcd, NULL);
sfree(fshift);
viol = fcd->disres.sumviol;
{
if (index[j] == forceparams[type].disres.label)
{
- vvindex[j] = pow(fcd->disres.Rt_6[0], -1.0/6.0);
+ vvindex[j] = std::pow(fcd->disres.Rt_6[0], static_cast<real>(-1.0/6.0));
}
}
}
default:
gmx_incons("Dumping violations");
}
- viol_max = max(viol_max, viol);
+ viol_max = std::max(viol_max, viol);
if (viol > 0)
{
nviol++;
drs[i].bCore = is_core(i, isize, index);
drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
drs[i].r = dr->aver1[i]/nsteps;
- drs[i].rT3 = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
- drs[i].rT6 = pow(dr->aver_6[i]/nsteps, -1.0/6.0);
- drs[i].viol = max(0, drs[i].r-drs[i].up1);
- drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
- drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
+ drs[i].rT3 = std::pow(dr->aver_3[i]/nsteps, static_cast<real>(-1.0/3.0));
+ drs[i].rT6 = std::pow(dr->aver_6[i]/nsteps, static_cast<real>(-1.0/6.0));
+ drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
+ drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
+ drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
if (atoms)
{
int j1 = disres->iatoms[j+1];
{
gmx_fatal(FARGS, "dr[%d].aver_3[%d] = %f", k, i, dr[k].aver_3[i]);
}
- drs[i].rT3 = pow(dr[k].aver_3[i]/dr[k].nframes, -1.0/3.0);
- drs[i].rT6 = pow(dr[k].aver_6[i]/dr[k].nframes, -1.0/6.0);
- drs[i].viol = max(0, drs[i].r-drs[i].up1);
- drs[i].violT3 = max(0, drs[i].rT3-drs[i].up1);
- drs[i].violT6 = max(0, drs[i].rT6-drs[i].up1);
+ drs[i].rT3 = std::pow(dr[k].aver_3[i]/dr[k].nframes, static_cast<real>(-1.0/3.0));
+ drs[i].rT6 = std::pow(dr[k].aver_6[i]/dr[k].nframes, static_cast<real>(-1.0/6.0));
+ drs[i].viol = std::max(0.0, static_cast<double>(drs[i].r-drs[i].up1));
+ drs[i].violT3 = std::max(0.0, static_cast<double>(drs[i].rT3-drs[i].up1));
+ drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6-drs[i].up1));
sumV += drs[i].viol;
sumVT3 += drs[i].violT3;
sumVT6 += drs[i].violT6;
- maxV = max(maxV, drs[i].viol);
- maxVT3 = max(maxVT3, drs[i].violT3);
- maxVT6 = max(maxVT6, drs[i].violT6);
+ maxV = std::max(maxV, static_cast<double>(drs[i].viol));
+ maxVT3 = std::max(maxVT3, static_cast<double>(drs[i].violT3));
+ maxVT6 = std::max(maxVT6, static_cast<double>(drs[i].violT6));
}
- if (strcmp(clust_name[k], "1000") == 0)
+ if (std::strcmp(clust_name[k], "1000") == 0)
{
mmm++;
}
int *resnr;
int n_res, a_offset, mb, mol, a;
t_atoms *atoms;
- int iii, i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
+ int i, j, nra, nratoms, tp, ri, rj, index, nlabel, label;
atom_id ai, aj, *ptr;
real **matrix, *t_res, hi, *w_dr, rav, rviol;
t_rgb rlo = { 1, 1, 1 };
rj = resnr[aj];
if (bThird)
{
- rav = pow(dr->aver_3[i]/nsteps, -1.0/3.0);
+ rav = std::pow(dr->aver_3[i]/nsteps, static_cast<real>(-1.0/3.0));
}
else
{
{
fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
}
- rviol = max(0, rav-idef->iparams[tp].disres.up1);
+ rviol = std::max(static_cast<real>(0.0), rav-idef->iparams[tp].disres.up1);
matrix[ri][rj] += w_dr[i]*rviol;
matrix[rj][ri] += w_dr[i]*rviol;
- hi = max(hi, matrix[ri][rj]);
- hi = max(hi, matrix[rj][ri]);
+ hi = std::max(hi, matrix[ri][rj]);
+ hi = std::max(hi, matrix[rj][ri]);
}
}
*/
#include "gmxpre.h"
-#include <stdlib.h>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/strdb.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/macros.h"
#include "gromacs/legacyheaders/typedefs.h"
{
fgets2(buf, STRLEN, tapeout);
}
- while (strstr(buf, "KAPPA") == NULL);
+ while (std::strstr(buf, "KAPPA") == NULL);
if (bFirst)
{
/* Since we also have empty lines in the dssp output (temp) file,
for (i = 0; i < nr; i++)
{
c.c1 = ssbuf[i];
- mat->matrix[frame][i] = max(0, searchcmap(mat->nmap, mat->map, c));
+ mat->matrix[frame][i] = std::max(static_cast<t_matelmt>(0), searchcmap(mat->nmap, mat->map, c));
}
frame++;
mat->nx = frame;
for (i = 0; (i < atoms->nr); i++)
{
- if (strcmp(*(atoms->atomname[i]), "OXT") == 0)
+ if (std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
{
*atoms->atomname[i] = OOO;
}
- else if (strcmp(*(atoms->atomname[i]), "O1") == 0)
+ else if (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
{
*atoms->atomname[i] = OOO;
}
- else if (strcmp(*(atoms->atomname[i]), "OC1") == 0)
+ else if (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
{
*atoms->atomname[i] = OOO;
}
{
for (j = 0; j < nres; j++)
{
- lo = min(lo, accr[i][j]);
- hi = max(hi, accr[i][j]);
+ lo = std::min(lo, accr[i][j]);
+ hi = std::max(hi, accr[i][j]);
}
}
fp = gmx_ffopen(fn, "w");
- nlev = hi-lo+1;
+ nlev = static_cast<int>(hi-lo+1);
write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)",
"Time", "Residue Index", nframe, nres,
mat->axis_x, mat->axis_y, accr, lo, hi, rlo, rhi, &nlev);
{
fprintf(fp, "@ subtitle \"Structure = ");
}
- for (s = 0; s < strlen(ss_string); s++)
+ for (s = 0; s < std::strlen(ss_string); s++)
{
if (s > 0)
{
for (f = 0; f < mat->nx; f++)
{
ss_count = 0;
- for (s = 0; s < (size_t)mat->nmap; s++)
+ for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
{
count[s] = 0;
}
count[mat->matrix[f][r]]++;
total[mat->matrix[f][r]]++;
}
- for (s = 0; s < (size_t)mat->nmap; s++)
+ for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
{
- if (strchr(ss_string, map[s].code.c1))
+ if (std::strchr(ss_string, map[s].code.c1))
{
ss_count += count[s];
total_count += count[s];
}
}
fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
- for (s = 0; s < (size_t)mat->nmap; s++)
+ for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
{
fprintf(fp, " %5d", count[s]);
}
}
/* now print column totals */
fprintf(fp, "%-8s %5d", "# Totals", total_count);
- for (s = 0; s < (size_t)mat->nmap; s++)
+ for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
{
fprintf(fp, " %5d", total[s]);
}
fprintf(fp, "\n");
/* now print percentages */
- fprintf(fp, "%-8s %5.2f", "# SS %", total_count / (real) (mat->nx * mat->ny));
- for (s = 0; s < (size_t)mat->nmap; s++)
+ fprintf(fp, "%-8s %5.2f", "# SS %", total_count / static_cast<real>(mat->nx * mat->ny));
+ for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
{
- fprintf(fp, " %5.2f", total[s] / (real) (mat->nx * mat->ny));
+ fprintf(fp, " %5.2f", total[s] / static_cast<real>(mat->nx * mat->ny));
}
fprintf(fp, "\n");
}
fprintf(stderr, "There are %d residues in your selected group\n", nres);
- strcpy(pdbfile, "ddXXXXXX");
+ std::strcpy(pdbfile, "ddXXXXXX");
gmx_tmpnam(pdbfile);
if ((tmpf = fopen(pdbfile, "w")) == NULL)
{
fclose(tmpf);
}
- strcpy(tmpfile, "ddXXXXXX");
+ std::strcpy(tmpfile, "ddXXXXXX");
gmx_tmpnam(tmpfile);
if ((tmpf = fopen(tmpfile, "w")) == NULL)
{
for (i = 0; i < atoms->nres; i++)
{
- av_area[i] = (average_area[i] / (real)nframe);
+ av_area[i] = (average_area[i] / static_cast<real>(nframe));
}
norm_acc(atoms, nres, av_area, norm_av_area);
*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
static double FD(double Delta, double f)
{
- return (2*pow(Delta, -4.5)*pow(f, 7.5) -
- 6*pow(Delta, -3)*pow(f, 5) -
- pow(Delta, -1.5)*pow(f, 3.5) +
- 6*pow(Delta, -1.5)*pow(f, 2.5) +
+ return (2*std::pow(Delta, -4.5)*std::pow(f, 7.5) -
+ 6*std::pow(Delta, -3.0)*std::pow(f, 5.0) -
+ std::pow(Delta, -1.5)*std::pow(f, 3.5) +
+ 6*std::pow(Delta, -1.5)*std::pow(f, 2.5) +
2*f - 2);
}
static double YYY(double f, double y)
{
- return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
+ return (2*std::pow(y*f, 3.0) - sqr(f)*y*(1+6*y) +
(2+6*y)*f - 2);
}
{
return 0;
}
- return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
+ return ((1+y+sqr(y)-std::pow(y, 3.0))/(std::pow(1-y, 3.0)));
}
static double bisector(double Delta, double tol,
double ff0, double ff1,
double ff(double, double))
{
- double fd0, fd, fd1, f, f0, f1;
+ double fd, f, f0, f1;
double tolmin = 1e-8;
f0 = ff0;
do
{
- fd0 = ff(Delta, f0);
- fd1 = ff(Delta, f1);
f = (f0+f1)*0.5;
fd = ff(Delta, f);
if (fd < 0)
{
double y1, y2;
- y1 = pow(f/Delta, 1.5);
+ y1 = std::pow(f/Delta, 1.5);
y2 = bisector(f, toler, 0, 10000, YYY);
- if (fabs((y1-y2)/(y1+y2)) > 100*toler)
+ if (std::abs((y1-y2)/(y1+y2)) > 100*toler)
{
fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
y1, y2);
{
double fy = f*y;
- return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
+ return BOLTZ*(std::log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
}
static real wCsolid(real nu, real beta)
}
else
{
- ebn = exp(bhn);
+ ebn = std::exp(bhn);
koko = sqr(1-ebn);
return sqr(bhn)*ebn/koko;
}
}
else
{
- return bhn/gmx_expm1(bhn) - gmx_log1p(-exp(-bhn));
+ return bhn/gmx_expm1(bhn) - gmx_log1p(-std::exp(-bhn));
}
}
}
else
{
- return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
+ return std::log((1-std::exp(-bhn))/(std::exp(-bhn/2))) - std::log(bhn);
}
}
matrix box;
int gnx;
char title[256];
- real t0, t1, m;
+ real t0, t1;
t_trxstatus *status;
- int nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
- double rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
+ int nV, nframes, n_alloc, i, j, fftcode, Nmol, Natom;
+ double rho, dt, Vsum, V, tmass, dostot, dos2;
real **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
output_env_t oenv;
gmx_fft_t fft;
- double cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
+ double cP, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
double wCdiff, wSdiff, wAdiff, wEdiff;
int grpNatoms;
atom_id *index;
n_alloc = 0;
nframes = 0;
- Vsum = V2sum = 0;
+ Vsum = 0;
nV = 0;
do
{
if (fr.bBox)
{
V = det(fr.box);
- V2sum += V*V;
Vsum += V;
nV++;
}
dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
if (bAbsolute)
{
- dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
+ dos[DOS][j] = bfac*std::sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
}
else
{
DoS0 = dos[DOS][0];
/* Note this eqn. is incorrect in Pascal2011a! */
- Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
- pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
+ Delta = ((2*DoS0/(9*Natom))*std::sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
+ std::pow((Natom/V), 1.0/3.0)*std::pow(6.0/M_PI, 2.0/3.0));
f = calc_fluidicity(Delta, toler);
y = calc_y(f, Delta, toler);
z = calc_compress(y);
- Sig = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
+ Sig = BOLTZ*(5.0/2.0+std::log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
Shs = Sig+calc_Shs(f, y);
rho = (tmass*AMU)/(V*NANO*NANO*NANO);
- sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
+ sigHS = std::pow(6*y*V/(M_PI*Natom), 1.0/3.0);
fprintf(fplog, "System = \"%s\"\n", title);
fprintf(fplog, "Nmol = %d\n", Nmol);
#include "gromacs/fileio/trx.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/legacyheaders/copyrite.h"
#include "gromacs/legacyheaders/macros.h"
#include "gromacs/math/vec.h"
#define NFILE asize(fnm)
- const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL, *out_datfile = NULL;
+ const char *in_trajfile, *out_xvgrkfile = NULL, *out_xvginstefffile = NULL, *out_xvgrhistfile = NULL, *out_xvgkhistfile = NULL, *out_datfile = NULL;
gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE;
int ndon, nacc;
atom_id *donindex, *accindex;
/* Check command line options for filenames and set bool flags when switch used*/
in_trajfile = opt2fn("-f", NFILE, fnm);
- in_ndxfile = opt2fn("-n", NFILE, fnm);
out_xvgrkfile = opt2fn("-ot", NFILE, fnm);
out_xvgrhistfile = opt2fn("-rhist", NFILE, fnm);
out_xvgkhistfile = opt2fn("-khist", NFILE, fnm);
unitv(donvec, donvec);
unitv(accvec, accvec);
- svmul((real) 1. / ndon, donpos, donpos);
- svmul((real) 1. / nacc, accpos, accpos);
+ svmul(1.0 / ndon, donpos, donpos);
+ svmul(1.0 / nacc, accpos, accpos);
if (bPBCdist)
{
for (i = 1; i < rkcount; i++)
{
- bin = (int) ((rvalues[i] - rmin) / rincr);
+ bin = static_cast<int>((rvalues[i] - rmin) / rincr);
rhist[bin] += 1;
}
if (bNormHist)
for (i = 1; i < rkcount; i++)
{
- bin = (int) ((kappa2values[i] - kmin) / kincr);
+ bin = static_cast<int>((kappa2values[i] - kmin) / kincr);
khist[bin] += 1;
}
if (bNormHist)
*/
#include "gmxpre.h"
+#include <cmath>
+
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/trxio.h"
mat4 Rx, Ry, Rz, Rinvy, Rinvz, Mtot;
mat4 temp1, temp2, temp3;
vec4 xv;
- int i, j, ai;
+ int i, ai;
rvec_sub(tail, head, arrow);
arrow_len = norm(arrow);
}
/* Compute theta and phi that describe the arrow */
- theta = acos(arrow[ZZ]/arrow_len);
- phi = atan2(arrow[YY]/arrow_len, arrow[XX]/arrow_len);
+ theta = std::acos(arrow[ZZ]/arrow_len);
+ phi = std::atan2(arrow[YY]/arrow_len, arrow[XX]/arrow_len);
if (debug)
{
fprintf(debug, "Phi = %.1f, Theta = %.1f\n", RAD2DEG*phi, RAD2DEG*theta);
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
-typedef struct
-{
- char sanm[12];
- int natm;
- int nw;
- char anm[6][12];
-} t_simat;
-
-typedef struct
-{
- char reso[12];
- char resn[12];
- int nsatm;
- t_simat sat[3];
-} t_simlist;
real calc_mass(t_atoms *atoms, gmx_bool bGetMass, gmx_atomprop_t aps)
{
return tmass;
}
-real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
- rvec max, gmx_bool bDiam)
+real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec minval,
+ rvec maxval, gmx_bool bDiam)
{
real diam2, d;
- char *grpnames;
int ii, i, j;
clear_rvec(geom_center);
diam2 = 0;
if (isize == 0)
{
- clear_rvec(min);
- clear_rvec(max);
+ clear_rvec(minval);
+ clear_rvec(maxval);
}
else
{
}
for (j = 0; j < DIM; j++)
{
- min[j] = max[j] = x[ii][j];
+ minval[j] = maxval[j] = x[ii][j];
}
for (i = 0; i < isize; i++)
{
rvec_inc(geom_center, x[ii]);
for (j = 0; j < DIM; j++)
{
- if (x[ii][j] < min[j])
+ if (x[ii][j] < minval[j])
{
- min[j] = x[ii][j];
+ minval[j] = x[ii][j];
}
- if (x[ii][j] > max[j])
+ if (x[ii][j] > maxval[j])
{
- max[j] = x[ii][j];
+ maxval[j] = x[ii][j];
}
}
if (bDiam)
for (j = i + 1; j < isize; j++)
{
d = distance2(x[ii], x[index[j]]);
- diam2 = max(d, diam2);
+ diam2 = std::max(d, diam2);
}
}
else
for (j = i + 1; j < isize; j++)
{
d = distance2(x[i], x[j]);
- diam2 = max(d, diam2);
+ diam2 = std::max(d, diam2);
}
}
}
svmul(1.0 / isize, geom_center, geom_center);
}
- return sqrt(diam2);
+ return std::sqrt(diam2);
}
void center_conf(int natom, rvec *x, rvec center, rvec geom_cent)
void set_pdb_conf_bfac(int natoms, int nres, t_atoms *atoms, int n_bfac,
double *bfac, int *bfac_nr, gmx_bool peratom)
{
- FILE *out;
real bfac_min, bfac_max;
int i, n;
gmx_bool found;
bfac_max /= 10;
bfac_min /= 10;
}
- while ((fabs(bfac_max) < 0.5) && (fabs(bfac_min) < 0.5))
+ while ((std::abs(bfac_max) < 0.5) && (std::abs(bfac_min) < 0.5))
{
fprintf(stderr,
"Range of values for B-factors too small (min %g, max %g) "
zmin = 1e10;
for (i = 0; (i < natoms); i++)
{
- xmin = min(xmin, x[i][XX]);
- ymin = min(ymin, x[i][YY]);
- zmin = min(zmin, x[i][ZZ]);
- bfac_min = min(bfac_min, atoms->pdbinfo[i].bfac);
- bfac_max = max(bfac_max, atoms->pdbinfo[i].bfac);
+ xmin = std::min(xmin, x[i][XX]);
+ ymin = std::min(ymin, x[i][YY]);
+ zmin = std::min(zmin, x[i][ZZ]);
+ bfac_min = std::min(bfac_min, atoms->pdbinfo[i].bfac);
+ bfac_max = std::max(bfac_max, atoms->pdbinfo[i].bfac);
}
fprintf(stderr, "B-factors range from %g to %g\n", bfac_min, bfac_max);
for (i = 1; (i < 12); i++)
a0++;
r0++;
- nx = (int) (gridsize[XX] + 0.5);
- ny = (int) (gridsize[YY] + 0.5);
- nz = (int) (gridsize[ZZ] + 0.5);
+ nx = static_cast<int>(gridsize[XX] + 0.5);
+ ny = static_cast<int>(gridsize[YY] + 0.5);
+ nz = static_cast<int>(gridsize[ZZ] + 0.5);
nbox = nx * ny * nz;
if (TRICLINIC(box))
{
}
for (i = 0; i < 24; i += 2)
{
- fprintf(out, "CONECT%5d%5d\n", a0 + rectedge[i], a0 + rectedge[i
- + 1]);
+ fprintf(out, "CONECT%5d%5d\n", a0 + rectedge[i], a0 + rectedge[i + 1]);
}
}
}
real ux, uy, uz, costheta, sintheta;
costheta = cos_angle(principal_axis, targetvec);
- sintheta = sqrt(1.0-costheta*costheta); /* sign is always positive since 0<theta<pi */
+ sintheta = s
td::sqrt(1.0-costheta*costheta); /* sign is always positive since 0<theta<pi */
/* Determine rotation from cross product with target vector */
cprod(principal_axis, targetvec, rotvec);
"For complex molecules, the periodicity removal routine may break down, "
"in that case you can use [gmx-trjconv]."
};
- static real dist = 0.0, rbox = 0.0, to_diam = 0.0;
+ static real dist = 0.0;
static gmx_bool bNDEF = FALSE, bRMPBC = FALSE, bCenter = FALSE, bReadVDW =
FALSE, bCONECT = FALSE;
static gmx_bool peratom = FALSE, bLegend = FALSE, bOrient = FALSE, bMead =
t_topology *top = NULL;
t_atoms atoms;
char *grpname, *sgrpname, *agrpname;
- int isize, ssize, tsize, asize;
- atom_id *index, *sindex, *tindex, *aindex;
- rvec *x, *v, gc, min, max, size;
+ int isize, ssize, asize;
+ atom_id *index, *sindex, *aindex;
+ rvec *x, *v, gc, rmin, rmax, size;
int ePBC;
matrix box, rotmatrix, trans;
rvec princd, tmpvec;
- gmx_bool bIndex, bSetSize, bSetAng, bCubic, bDist, bSetCenter, bAlign;
+ gmx_bool bIndex, bSetSize, bSetAng, bDist, bSetCenter, bAlign;
gmx_bool bHaveV, bScale, bRho, bTranslate, bRotate, bCalcGeom, bCalcDiam;
- real xs, ys, zs, xcent, ycent, zcent, diam = 0, mass = 0, d, vdw;
+ real diam = 0, mass = 0, d, vdw;
gmx_atomprop_t aps;
gmx_conect conect;
output_env_t oenv;
}
bScale = bScale || bRho;
bCalcGeom = bCenter || bRotate || bOrient || bScale;
- bCalcDiam = btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o';
+
+ GMX_RELEASE_ASSERT(btype[0] != NULL, "Option setting inconsistency; btype[0] is NULL");
+
+ bCalcDiam = (btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o');
infile = ftp2fn(efSTX, NFILE, fnm);
if (bMead)
{
real vol = det(box);
printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
- vol, 100*((int)(vol*4.5)));
+ vol, 100*(static_cast<int>(vol*4.5)));
}
if (bMead || bGrasp || bCONECT)
{
sig6 = c12/c6;
}
- vdw = 0.5*pow(sig6, 1.0/6.0);
+ vdw = 0.5*std::pow(sig6, static_cast<real>(1.0/6.0));
}
else
{
ssize = atoms.nr;
sindex = NULL;
}
- diam = calc_geom(ssize, sindex, x, gc, min, max, bCalcDiam);
- rvec_sub(max, min, size);
+ diam = calc_geom(ssize, sindex, x, gc, rmin, rmax, bCalcDiam);
+ rvec_sub(rmax, rmin, size);
printf(" system size :%7.3f%7.3f%7.3f (nm)\n",
size[XX], size[YY], size[ZZ]);
if (bCalcDiam)
norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
printf(" box angles :%7.2f%7.2f%7.2f (degrees)\n",
norm2(box[ZZ]) == 0 ? 0 :
- RAD2DEG*acos(cos_angle_no_table(box[YY], box[ZZ])),
+ RAD2DEG*std::acos(cos_angle_no_table(box[YY], box[ZZ])),
norm2(box[ZZ]) == 0 ? 0 :
- RAD2DEG*acos(cos_angle_no_table(box[XX], box[ZZ])),
+ RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[ZZ])),
norm2(box[YY]) == 0 ? 0 :
- RAD2DEG*acos(cos_angle_no_table(box[XX], box[YY])));
+ RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[YY])));
printf(" box volume :%7.2f (nm^3)\n", det(box));
}
"zero mass (%g) or volume (%g)\n", mass, vol);
}
- scale[XX] = scale[YY] = scale[ZZ] = pow(dens/rho, 1.0/3.0);
+ scale[XX] = scale[YY] = scale[ZZ] = std::pow(dens/rho, static_cast<real>(1.0/3.0));
fprintf(stderr, "Scaling all box vectors by %g\n", scale[XX]);
}
scale_conf(atoms.nr, x, box, scale);
if (bCalcGeom)
{
/* recalc geometrical center and max and min coordinates and size */
- calc_geom(ssize, sindex, x, gc, min, max, FALSE);
- rvec_sub(max, min, size);
+ calc_geom(ssize, sindex, x, gc, rmin, rmax, FALSE);
+ rvec_sub(rmax, rmin, size);
if (bScale || bOrient || bRotate)
{
printf("new system size : %6.3f %6.3f %6.3f\n",
}
}
- if (bSetSize || bDist || (btype[0][0] == 't' && bSetAng))
+ if ((btype[0] != NULL) && (bSetSize || bDist || (btype[0][0] == 't' && bSetAng)))
{
ePBC = epbcXYZ;
if (!(bSetSize || bDist))
box[YY][YY] = d;
box[ZZ][XX] = d/2;
box[ZZ][YY] = d/2;
- box[ZZ][ZZ] = d*sqrt(2)/2;
+ box[ZZ][ZZ] = d*std::sqrt(2.0)/2.0;
}
else
{
box[XX][XX] = d;
box[YY][XX] = d/3;
- box[YY][YY] = d*sqrt(2)*2/3;
+ box[YY][YY] = d*std::sqrt(2.0)*2.0/3.0;
box[ZZ][XX] = -d/3;
- box[ZZ][YY] = d*sqrt(2)/3;
- box[ZZ][ZZ] = d*sqrt(6)/3;
+ box[ZZ][YY] = d*std::sqrt(2.0)/3.0;
+ box[ZZ][ZZ] = d*std::sqrt(6.0)/3.0;
}
break;
}
/* print some */
if (bCalcGeom)
{
- calc_geom(ssize, sindex, x, gc, min, max, FALSE);
+ calc_geom(ssize, sindex, x, gc, rmin, rmax, FALSE);
printf("new center :%7.3f%7.3f%7.3f (nm)\n", gc[XX], gc[YY], gc[ZZ]);
}
if (bOrient || bScale || bDist || bSetSize)
norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
printf("new box angles :%7.2f%7.2f%7.2f (degrees)\n",
norm2(box[ZZ]) == 0 ? 0 :
- RAD2DEG*acos(cos_angle_no_table(box[YY], box[ZZ])),
+ RAD2DEG*std::acos(cos_angle_no_table(box[YY], box[ZZ])),
norm2(box[ZZ]) == 0 ? 0 :
- RAD2DEG*acos(cos_angle_no_table(box[XX], box[ZZ])),
+ RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[ZZ])),
norm2(box[YY]) == 0 ? 0 :
- RAD2DEG*acos(cos_angle_no_table(box[XX], box[YY])));
+ RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[YY])));
printf("new box volume :%7.2f (nm^3)\n", det(box));
}
*/
#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/fileio/enxio.h"
real *readtime, real *timestep, int *nremax)
{
/* Check number of energy terms and start time of all files */
- int f, i, nre, nremin = 0, nresav = 0;
+ int f, nre, nremin = 0, nresav = 0;
ener_file_t in;
real t1, t2;
char inputstring[STRLEN];
}
else
{
- nremin = min(nremin, fr->nre);
- *nremax = max(*nremax, fr->nre);
+ nremin = std::min(nremin, fr->nre);
+ *nremax = std::max(*nremax, fr->nre);
if (nre != nresav)
{
fprintf(stderr,
{
gmx_fatal(FARGS, "Error reading user input");
}
- inputstring[strlen(inputstring)-1] = 0;
+ inputstring[std::strlen(inputstring)-1] = 0;
if (inputstring[0] == 'c' || inputstring[0] == 'C')
{
fr->ener[i].eav +
dsqr(ee_sum[i].esum/nsum
- (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
- nsum*(nsum + fr->nsum)/(double)fr->nsum;
+ nsum*(nsum + fr->nsum)/static_cast<double>(fr->nsum);
ee_sum[i].esum += fr->ener[i].esum;
}
}
t_enxframe *fr, *fro;
gmx_int64_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
t_energy *ee_sum;
- gmx_int64_t lastfilestep, laststep, startstep, startstep_file = 0;
+ gmx_int64_t lastfilestep, laststep, startstep_file = 0;
int noutfr;
- int nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL;
+ int nre, nremax, this_nre, nfile, f, i, kkk, nset, *set = NULL;
double last_t;
char **fnms;
- real *readtime, *settime, timestep, t1, tadjust;
- char inputstring[STRLEN], *chptr, buf[22], buf2[22], buf3[22];
- gmx_bool ok;
+ real *readtime, *settime, timestep, tadjust;
+ char buf[22], buf2[22];
int *cont_type;
gmx_bool bNewFile, bFirst, bNewOutput;
output_env_t oenv;
timestep = 0.0;
snew(fnms, argc);
lastfilestep = 0;
- laststep = startstep = 0;
+ laststep = 0;
nfile = opt2fns(&fnms, "-f", NFILE, fnm);
if (bFirst)
{
bFirst = FALSE;
- startstep = fr->step;
}
if (bWrite)
{
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/enxio.h"
static int search_str2(int nstr, char **str, char *key)
{
int i, n;
- int keylen = strlen(key);
+ int keylen = std::strlen(key);
/* Linear search */
n = 0;
while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
fprintf(stderr, "Will read groupnames from inputfile\n");
ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
fprintf(stderr, "Read %d groups\n", ngroups);
- snew(set, sqr(ngroups)*egNR/2);
+ snew(set, static_cast<size_t>(sqr(ngroups)*egNR/2));
n = 0;
prevk = 0;
for (i = 0; (i < ngroups); i++)
#endif
for (k = prevk; (k < prevk+nre); k++)
{
- if (strcmp(enm[k%nre].name, groupname) == 0)
+ if (std::strcmp(enm[k%nre].name, groupname) == 0)
{
set[n++] = k;
break;
expE = 0;
for (k = 0; (k < nenergy); k++)
{
- expE += exp(beta*(e[i][k]-eaver[i]));
+ expE += std::exp(beta*(e[i][k]-eaver[i]));
}
- efree[i] = log(expE/nenergy)/beta + eaver[i];
+ efree[i] = std::log(expE/nenergy)/beta + eaver[i];
if (bRef)
{
n = search_str2(neref, erefres, groups[i]);
*/
#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);
}
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
sum = 0;
for (i = 0; i < nffr; i++)
{
- filt[i] = cos(2*M_PI*(i - nf + 1)/(real)flen) + 1;
+ filt[i] = std::cos(2*M_PI*(i - nf + 1)/static_cast<real>(flen)) + 1;
sum += filt[i];
}
fprintf(stdout, "filter weights:");
*/
#include "gmxpre.h"
-#include <ctype.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cctype>
+#include <cstdlib>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
rvec dx;
gmx_int64_t maxrand;
- ei = -1;
nw = *nwater;
maxrand = nw;
maxrand *= 1000;
do
{
- ei = nw*gmx_rng_uniform_real(rng);
+ ei = static_cast<int>(nw*gmx_rng_uniform_real(rng));
maxrand--;
}
while (bSet[ei] && (maxrand > 0));
int i;
str = gmx_strdup(mname);
- i = strlen(str)-1;
- while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
+ i = std::strlen(str)-1;
+ while (i > 1 && (std::isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
{
str[i] = '\0';
i--;
{
FILE *fpin, *fpout;
char buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
- int line, i, nsol, nmol_line, sol_line, nsol_last;
+ int line, i, nmol_line, sol_line, nsol_last;
gmx_bool bMolecules;
char temporary_filename[STRLEN];
printf("\nProcessing topology\n");
fpin = gmx_ffopen(topinout, "r");
- strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
+ std::strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
gmx_tmpnam(temporary_filename);
fpout = gmx_ffopen(temporary_filename, "w");
while (fgets(buf, STRLEN, fpin))
{
line++;
- strcpy(buf2, buf);
- if ((temp = strchr(buf2, '\n')) != NULL)
+ std::strcpy(buf2, buf);
+ if ((temp = std::strchr(buf2, '\n')) != NULL)
{
temp[0] = '\0';
}
if (buf2[0] == '[')
{
buf2[0] = ' ';
- if ((temp = strchr(buf2, '\n')) != NULL)
+ if ((temp = std::strchr(buf2, '\n')) != NULL)
{
temp[0] = '\0';
}
rtrim(buf2);
- if (buf2[strlen(buf2)-1] == ']')
+ if (buf2[std::strlen(buf2)-1] == ']')
{
- buf2[strlen(buf2)-1] = '\0';
+ buf2[std::strlen(buf2)-1] = '\0';
ltrim(buf2);
rtrim(buf2);
bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
+
+#include <algorithm>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
{
if (disre_frac > 0)
{
- dd = min(disre_dist, disre_frac*d);
+ dd = std::min(disre_dist, disre_frac*d);
}
else
{
dd = disre_dist;
}
- lo = max(0, d-dd);
+ lo = std::max(static_cast<real>(0.0), d-dd);
hi = d+dd;
fprintf(out, "%5d %5d %1d %5d %10d %10g %10g %10g %10g\n",
ind_grp[i]+1, ind_grp[j]+1, 1, k, 1,
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
}
for (m = 0; (m < DIM); m++)
{
- d[m] = sqrt(d[m]/tm);
+ d[m] = std::sqrt(d[m]/tm);
}
#ifdef DEBUG
pr_rvecs(stderr, 0, "trans", trans, DIM);
ii = index[i];
if (bQ)
{
- m0 = fabs(atom[ii].q);
+ m0 = std::abs(atom[ii].q);
}
else
{
for (m = 0; (m < DIM); m++)
{
- gvec[m] = sqrt((gyro-comp[m])/tm);
+ gvec[m] = std::sqrt((gyro-comp[m])/tm);
}
- return sqrt(gyro/tm);
+ return std::sqrt(gyro/tm);
}
void calc_gyro_z(rvec x[], matrix box,
}
for (j = 0; j < 2; j++)
{
- zi = zf + j;
+ zi = static_cast<int>(zf + j);
if (zi == nz)
{
zi = 0;
}
- w = atom[ii].m*(1 + cos(M_PI*(zf - zi)));
+ w = atom[ii].m*(1 + std::cos(M_PI*(zf - zi)));
inertia[zi][0] += w*sqr(x[ii][YY]);
inertia[zi][1] += w*sqr(x[ii][XX]);
inertia[zi][2] -= w*x[ii][XX]*x[ii][YY];
{
inertia[j][i] /= tm[j];
}
- sdet = sqrt(sqr(inertia[j][0] - inertia[j][1]) + 4*sqr(inertia[j][2]));
- e1 = sqrt(0.5*(inertia[j][0] + inertia[j][1] + sdet));
- e2 = sqrt(0.5*(inertia[j][0] + inertia[j][1] - sdet));
+ sdet = std::sqrt(sqr(inertia[j][0] - inertia[j][1]) + 4*sqr(inertia[j][2]));
+ e1 = std::sqrt(0.5*(inertia[j][0] + inertia[j][1] + sdet));
+ e2 = std::sqrt(0.5*(inertia[j][0] + inertia[j][1] - sdet));
fprintf(out, " %5.3f %5.3f", e1, e2);
}
fprintf(out, "\n");
real t, t0, tm, gyro;
int natoms;
char *grpname, title[256];
- int i, j, m, gnx, nam, mol;
+ int j, m, gnx, nam, mol;
atom_id *index;
output_env_t oenv;
gmx_rmpbc_t gpbc = NULL;
}
gyro = 0;
clear_rvec(gvec);
+ clear_rvec(gvec1);
clear_rvec(d);
+ clear_rvec(d1);
for (mol = 0; mol < nmol; mol++)
{
tm = sub_xcm(nz == 0 ? x_s : x, nam, index+mol*nam, top.atoms.atom, xcm, bQ);
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/trxio.h"
if (!*nslices)
{
- *nslices = (int)(box[axis][axis] * 10); /* default value */
+ *nslices = static_cast<int>(box[axis][axis] * 10); /* default value */
}
*/
if (bMicel)
- { /* this is for spherical interfaces */
- rvec_sub(com, x0[index[3*i]], normal); /* vector from Oxygen to COM */
- slice = norm(normal)/(*slWidth); /* spherical slice */
+ { /* this is for spherical interfaces */
+ rvec_sub(com, x0[index[3*i]], normal); /* vector from Oxygen to COM */
+ slice = static_cast<int>(norm(normal)/(*slWidth)); /* spherical slice */
sum[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
frame[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
/* this is for flat interfaces */
/* determine which slice atom is in */
- slice = (x0[index[3*i]][axis] / (*slWidth));
+ slice = static_cast<int>(x0[index[3*i]][axis] / (*slWidth));
if (slice < 0 || slice >= *nslices)
{
fprintf(stderr, "Coordinate: %f ", x0[index[3*i]][axis]);
#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);
*/
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
output_env_t oenv;
char buf[54];
t_trxstatus *status;
- int natoms, nre, nres;
+ int natoms, nres;
t_bb *bb;
- int i, j, m, nall, nbb, nca, teller, nSel = 0;
+ int i, j, nall, nbb, nca, teller;
atom_id *bbindex, *caindex, *allindex;
t_topology *top;
int ePBC;
{
fprintf(xf[efhRMSA].fp, "%10d %10g\n", r0+i, bb[i].rmsa/bb[i].nrms);
}
- fprintf(xf[efhAHX].fp, "%10d %10g\n", r0+i, (bb[i].nhx*100.0)/(real )teller);
+ fprintf(xf[efhAHX].fp, "%10d %10g\n", r0+i, (bb[i].nhx*100.0)/static_cast<real>(teller));
fprintf(xf[efhJCA].fp, "%10d %10g\n",
- r0+i, 140.3+(bb[i].jcaha/(double)teller));
+ r0+i, 140.3+(bb[i].jcaha/static_cast<double>(teller)));
}
for (i = 0; (i < efhNR); i++)
*/
#include "gmxpre.h"
-#include <math.h>
+#include <cmath>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/trxio.h"
t_topology *top = NULL;
real t;
- rvec *x = NULL, dx;
+ rvec *x = NULL;
matrix box;
t_trxstatus *status;
int natoms;
real theta1, theta2, theta3;
- int d, i, j, teller = 0;
+ int i, j, teller = 0;
int iCA, iSC;
atom_id *ind_CA;
atom_id *ind_SC;
char *gn_CA;
char *gn_SC;
- rvec averageaxis;
- rvec v1, v2, p1, p2, vtmp, vproj;
+ rvec v1, v2;
rvec *x_CA, *x_SC;
rvec *r12;
rvec *r23;
rvec *residuevector;
rvec *sidechainvector;
- rvec axes_t0[3];
- rvec axes[3];
rvec *residuehelixaxis_t0;
rvec *residuevector_t0;
rvec *axis3_t0;
real *rise, *residuerise;
real *residuebending;
- real tmp, rotangle;
+ real tmp;
real weight[3];
t_pbc pbc;
matrix A;
svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);
tmp = cos_angle(diff13[i], diff24[i]);
- twist[i] = 180.0/M_PI * acos( tmp );
- radius[i] = sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
- rise[i] = fabs(iprod(r23[i], helixaxis[i]));
+ twist[i] = 180.0/M_PI * std::acos( tmp );
+ radius[i] = std::sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
+ rise[i] = std::abs(iprod(r23[i], helixaxis[i]));
svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
svmul(radius[i]/norm(diff24[i]), diff24[i], v2);
residueradius[i] = 0.5*(radius[i-2]+radius[i-1]);
residuetwist[i] = 0.5*(twist[i-2]+twist[i-1]);
residuerise[i] = 0.5*(rise[i-2]+rise[i-1]);
- residuebending[i] = 180.0/M_PI*acos( cos_angle(helixaxis[i-2], helixaxis[i-1]) );
+ residuebending[i] = 180.0/M_PI*std::acos( cos_angle(helixaxis[i-2], helixaxis[i-1]) );
}
residueradius[iCA-2] = radius[iCA-4];
residuetwist[iCA-2] = twist[iCA-4];
fprintf(fpradius, "%15.12g ", residueradius[i]);
fprintf(fptwist, "%15.12g ", residuetwist[i]);
fprintf(fpbending, "%15.12g ", residuebending[i]);
-
- /* angle with local vector? */
- /*
- printf("res[%2d]: axis: %g %g %g origin: %g %g %g vector: %g %g %g angle: %g\n",i,
- residuehelixaxis[i][0],
- residuehelixaxis[i][1],
- residuehelixaxis[i][2],
- residueorigin[i][0],
- residueorigin[i][1],
- residueorigin[i][2],
- residuevector[i][0],
- residuevector[i][1],
- residuevector[i][2],
- 180.0/M_PI*acos( cos_angle(residuevector[i],residuehelixaxis[i]) ));
- */
- /* fprintf(fp,"%15.12g %15.12g %15.12g %15.12g %15.12g %15.12g\n",
- residuehelixaxis[i][0],
- residuehelixaxis[i][1],
- residuehelixaxis[i][2],
- residuevector[i][0],
- residuevector[i][1],
- residuevector[i][2]);
- */
}
}
fprintf(fprise, "\n");
copy_rvec(residuevector[i], newaxes[1]);
copy_rvec(axis3[i], newaxes[2]);
- /*
- printf("frame %d, i=%d:\n old: %g %g %g , %g %g %g , %g %g %g\n new: %g %g %g , %g %g %g , %g %g %g\n",
- teller,i,
- refaxes[0][0],refaxes[0][1],refaxes[0][2],
- refaxes[1][0],refaxes[1][1],refaxes[1][2],
- refaxes[2][0],refaxes[2][1],refaxes[2][2],
- newaxes[0][0],newaxes[0][1],newaxes[0][2],
- newaxes[1][0],newaxes[1][1],newaxes[1][2],
- newaxes[2][0],newaxes[2][1],newaxes[2][2]);
- */
-
/* rotate reference frame onto unit axes */
calc_fit_R(3, 3, weight, unitaxes, refaxes, A);
for (j = 0; j < 3; j++)
* A contains rotation column vectors.
*/
- /*
- printf("frame %d, i=%d, A: %g %g %g , %g %g %g , %g %g %g\n",
- teller,i,A[0][0],A[0][1],A[0][2],A[1][0],A[1][1],A[1][2],A[2][0],A[2][1],A[2][2]);
- */
-
- theta1 = 180.0/M_PI*atan2(A[0][2], A[0][0]);
- theta2 = 180.0/M_PI*asin(-A[0][1]);
- theta3 = 180.0/M_PI*atan2(A[2][1], A[1][1]);
+ theta1 = 180.0/M_PI*std::atan2(A[0][2], A[0][0]);
+ theta2 = 180.0/M_PI*std::asin(-A[0][1]);
+ theta3 = 180.0/M_PI*std::atan2(A[2][1], A[1][1]);
- tilt = sqrt(theta1*theta1+theta2*theta2);
+ tilt = std::sqrt(theta1*theta1+theta2*theta2);
rotation = theta3;
fprintf(fptheta1, "%15.12g ", theta1);
fprintf(fptheta2, "%15.12g ", theta2);
#include "gmxpre.h"
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxana/binsearch.h"
+#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/gmxana/gstat.h"
#include "gromacs/gmxana/powerspect.h"
#include "gromacs/legacyheaders/macros.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
/* Print name of first atom in all groups in index file */
rmean = 0;
for (j = 0; (j < 4); j++)
{
- r_nn[j][i] = sqrt(r_nn[j][i]);
+ r_nn[j][i] = std::sqrt(r_nn[j][i]);
rmean += r_nn[j][i];
}
rmean /= 4;
output_env_t oenv)
{
FILE *fpsg = NULL, *fpsk = NULL;
- char *sgslfn = "sg_ang_mesh"; /* Hardcoded filenames for debugging*/
- char *skslfn = "sk_dist_mesh";
t_topology top;
int ePBC;
- char title[STRLEN], subtitle[STRLEN];
+ char title[STRLEN];
t_trxstatus *status;
int natoms;
real t;
rvec *xtop, *x;
matrix box;
- real sg, sk, sgintf, pos;
+ real sg, sk, sgintf;
atom_id **index = NULL;
char **grpname = NULL;
int i, j, k, n, *isize, ng, nslicez, framenr;
read_tps_conf(fnTPS, title, &top, &ePBC, &xtop, NULL, box, FALSE);
- *nslicex = (int)(box[XX][XX]/binw + onehalf); /*Calculate slicenr from binwidth*/
- *nslicey = (int)(box[YY][YY]/binw + onehalf);
- nslicez = (int)(box[ZZ][ZZ]/binw + onehalf);
+ *nslicex = static_cast<int>(box[XX][XX]/binw + onehalf); /*Calculate slicenr from binwidth*/
+ *nslicey = static_cast<int>(box[YY][YY]/binw + onehalf);
+ nslicez = static_cast<int>(box[ZZ][ZZ]/binw + onehalf);
find_tetra_order_grid(top, ePBC, natoms, box, x, isize[0], index[0],
&sg, &sk, *nslicex, *nslicey, nslicez, sg_grid, sk_grid);
+ GMX_RELEASE_ASSERT(sk_fravg != NULL, "Trying to dereference NULL sk_fravg pointer");
for (i = 0; i < *nslicex; i++)
{
for (j = 0; j < *nslicey; j++)
if (framenr%tblock == 0)
{
+ GMX_RELEASE_ASSERT(sk_4d != NULL, "Trying to dereference NULL sk_4d pointer");
sk_4d[*nframes] = sk_fravg;
sg_4d[*nframes] = sg_fravg;
(*nframes)++;
/*Debugging for printing out the entire order parameter meshes.*/
if (debug)
{
- fpsg = xvgropen(sgslfn, "S\\sg\\N Angle Order Parameter / Meshpoint", "(nm)", "S\\sg\\N", oenv);
- fpsk = xvgropen(skslfn, "S\\sk\\N Distance Order Parameter / Meshpoint", "(nm)", "S\\sk\\N", oenv);
+ fpsg = xvgropen("sg_ang_mesh", "S\\sg\\N Angle Order Parameter / Meshpoint", "(nm)", "S\\sg\\N", oenv);
+ fpsk = xvgropen("sk_dist_mesh", "S\\sk\\N Distance Order Parameter / Meshpoint", "(nm)", "S\\sk\\N", oenv);
for (n = 0; n < (*nframes); n++)
{
fprintf(fpsg, "%i\n", n);
}
- /*sfree(perm);*/
sfree(sk_4d);
sfree(sg_4d);
- /*sfree(sg_grid);*/
- /*sfree(sk_grid);*/
-
-
}
static void writesurftoxpms(real ***surf, int tblocks, int xbins, int ybins, real bw, char **outfiles, int maplevels )
static gmx_bool bRawOut = FALSE;
int frames, xslices, yslices; /* Dimensions of interface arrays*/
real ***intfpos; /* Interface arrays (intfnr,t,xy) -potentially large */
- static char *normal_axis[] = { NULL, "z", "x", "y", NULL };
+ static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
t_pargs pa[] = {
{ "-d", FALSE, etENUM, {normal_axis},
trxfnm = ftp2fn(efTRX, NFILE, fnm);
/* Calculate axis */
- if (strcmp(normal_axis[0], "x") == 0)
+ GMX_RELEASE_ASSERT(normal_axis[0] != NULL, "Option setting inconsistency; normal_axis[0] is NULL");
+ if (std::strcmp(normal_axis[0], "x") == 0)
{
axis = XX;
}
- else if (strcmp(normal_axis[0], "y") == 0)
+ else if (std::strcmp(normal_axis[0], "y") == 0)
{
axis = YY;
}
- else if (strcmp(normal_axis[0], "z") == 0)
+ else if (std::strcmp(normal_axis[0], "z") == 0)
{
axis = ZZ;
}
*/
#include "gmxpre.h"
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/enxio.h"
/* Skip until we come to pressure */
for (i = 0; (i < F_NRE); i++)
{
- if (strcmp(names[i].name, interaction_function[F_PRES].longname) == 0)
+ if (std::strcmp(names[i].name, interaction_function[F_PRES].longname) == 0)
{
break;
}
snew(ld, 1);
for (; (i < nre); i++)
{
- if ((strstr(names[i].name, ligand) != NULL) &&
- (strstr(names[i].name, self) == NULL))
+ if ((std::strstr(names[i].name, ligand) != NULL) &&
+ (std::strstr(names[i].name, self) == NULL))
{
- if (strstr(names[i].name, "LJ") != NULL)
+ if (std::strstr(names[i].name, "LJ") != NULL)
{
ld->nlj++;
srenew(ld->lj, ld->nlj);
ld->lj[ld->nlj-1] = i;
}
- else if (strstr(names[i].name, "Coul") != NULL)
+ else if (std::strstr(names[i].name, "Coul") != NULL)
{
ld->nqq++;
srenew(ld->qq, ld->nqq);
FILE *out;
int nre, nframes = 0, ct = 0;
ener_file_t fp;
- gmx_bool bCont;
t_liedata *ld;
gmx_enxnm_t *enm = NULL;
t_enxframe *fr;
if (nframes > 0)
{
printf("DGbind = %.3f (%.3f)\n", lieaver/nframes,
- sqrt(lieav2/nframes-sqr(lieaver/nframes)));
+ std::sqrt(lieav2/nframes-sqr(lieaver/nframes)));
}
do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
#include "config.h"
+#include <cmath>
+
#include <algorithm>
#include "gromacs/commandline/pargs.h"
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += pow( tmp, -n );
+ nom += std::pow( tmp, -n );
}
for (i = SUMORDER; i > 0; i--)
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += pow( tmp, -n );
+ nom += std::pow( tmp, -n );
}
tmp = m / K;
tmp *= 2.0*M_PI;
- denom = pow( tmp, -n )+nom;
+ denom = std::pow( tmp, -n )+nom;
return -nom/denom;
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += pow( tmp, -2*n );
+ nom += std::pow( tmp, -2*n );
}
for (i = SUMORDER; i > 0; i--)
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += pow( tmp, -2*n );
+ nom += std::pow( tmp, -2*n );
}
for (i = -SUMORDER; i < SUMORDER+1; i++)
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- denom += pow( tmp, -n );
+ denom += std::pow( tmp, -n );
}
tmp = eps_poly1(m, K, n);
return nom / denom / denom + tmp*tmp;
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += i * pow( tmp, -2*n );
+ nom += i * std::pow( tmp, -2*n );
}
for (i = SUMORDER; i > 0; i--)
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += i * pow( tmp, -2*n );
+ nom += i * std::pow( tmp, -2*n );
}
for (i = -SUMORDER; i < SUMORDER+1; i++)
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- denom += pow( tmp, -n );
+ denom += std::pow( tmp, -n );
}
return 2.0 * M_PI * nom / denom / denom;
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += i * i * pow( tmp, -2*n );
+ nom += i * i * std::pow( tmp, -2*n );
}
for (i = SUMORDER; i > 0; i--)
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- nom += i * i * pow( tmp, -2*n );
+ nom += i * i * std::pow( tmp, -2*n );
}
for (i = -SUMORDER; i < SUMORDER+1; i++)
{
tmp = m / K + i;
tmp *= 2.0*M_PI;
- denom += pow( tmp, -n );
+ denom += std::pow( tmp, -n );
}
return 4.0 * M_PI * M_PI * nom / denom / denom;
for (i = -SUMORDER; i < 0; i++)
{
- tmp = -sin(2.0 * M_PI * i * K * rcoord);
+ tmp = -std::sin(2.0 * M_PI * i * K * rcoord);
tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
- tmp2 = pow(tmp1, -n);
+ tmp2 = std::pow(tmp1, -n);
nom += tmp * tmp2 * i;
denom += tmp2;
}
for (i = SUMORDER; i > 0; i--)
{
- tmp = -sin(2.0 * M_PI * i * K * rcoord);
+ tmp = -std::sin(2.0 * M_PI * i * K * rcoord);
tmp1 = 2.0 * M_PI * m / K + 2.0 * M_PI * i;
- tmp2 = pow(tmp1, -n);
+ tmp2 = std::pow(tmp1, -n);
nom += tmp * tmp2 * i;
denom += tmp2;
}
xtot = stopglobal*2+1;
if (PAR(cr))
{
- x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
+ x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
startlocal = startglobal + x_per_core*cr->nodeid;
stoplocal = startlocal + x_per_core -1;
if (stoplocal > stopglobal)
svmul(nz, info->recipbox[ZZ], tmpvec);
rvec_add(gridpxy, tmpvec, gridp);
tmp = norm2(gridp);
- coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
+ coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
coeff /= 2.0 * M_PI * info->volume * tmp;
coeff2 = tmp;
/* Here xtot is the number of samples taken for the Monte Carlo calculation
* of the average of term IV of equation 35 in Wang2010. Round up to a
* number of samples that is divisible by the number of nodes */
- x_per_core = static_cast<int>(ceil(info->fracself * nr / cr->nnodes));
+ x_per_core = static_cast<int>(std::ceil(info->fracself * nr / cr->nnodes));
xtot = x_per_core * cr->nnodes;
}
else
{
/* In this case we use all nr particle positions */
xtot = nr;
- x_per_core = static_cast<int>(ceil(static_cast<real>(xtot) / cr->nnodes));
+ x_per_core = static_cast<int>(std::ceil(static_cast<real>(xtot) / cr->nnodes));
}
startlocal = x_per_core * cr->nodeid;
{
for (i = 0; i < xtot; i++)
{
- numbers[i] = static_cast<int>(floor(gmx_rng_uniform_real(rng) * nr));
+ numbers[i] = static_cast<int>(std::floor(gmx_rng_uniform_real(rng) * nr));
}
}
/* Broadcast the random number array to the other nodes */
svmul(nz, info->recipbox[ZZ], tmpvec);
rvec_add(gridpxy, tmpvec, gridp);
tmp = norm2(gridp);
- coeff = exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
+ coeff = std::exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
coeff /= tmp;
e_rec3x += coeff*eps_self(nx, info->nkx[0], info->recipbox[XX], info->pme_order[0], x[ci]);
e_rec3y += coeff*eps_self(ny, info->nky[0], info->recipbox[YY], info->pme_order[0], x[ci]);
e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
*/
- e_rec = sqrt(e_rec1+e_rec2+e_rec3);
+ e_rec = std::sqrt(e_rec1+e_rec2+e_rec3);
return ONE_4PI_EPS0 * e_rec;
edir = info->e_dir[0];
erec = info->e_rec[0];
derr = edir-erec;
- while (fabs(derr/std::min(erec, edir)) > 1e-4)
+ while (std::abs(derr/std::min(erec, edir)) > 1e-4)
{
beta = info->ewald_beta[0];
if (MASTER(cr))
{
i++;
- fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, fabs(derr));
+ fprintf(stderr, "difference between real and rec. space error (step %d): %g\n", i, std::abs(derr));
fprintf(stderr, "old beta: %f\n", beta0);
fprintf(stderr, "new beta: %f\n", beta);
}
#include "config.h"
-#include <ctype.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cctype>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
#include <algorithm>
#include <sstream>
}
for (i = 0; i < nl-1; i++)
{
- if (fabs(y[0][i+1]-y[0][i]-opt->tabDz) > opt->tabDz/1e6)
+ if (std::abs(y[0][i+1]-y[0][i]-opt->tabDz) > opt->tabDz/1e6)
{
gmx_fatal(FARGS, "z-values in %s are not equally spaced.\n", ny, fn);
}
}
ist.str(line);
ist >> Buffer0 >> Buffer1 >> Buffer2;
- if (strcmp(Buffer1, "UMBRELLA"))
+ if (std::strcmp(Buffer1, "UMBRELLA"))
{
gmx_fatal(FARGS, "This does not appear to be a valid pdo file. Found %s, expected %s\n"
"(Found in first line: `%s')\n",
Buffer1, "UMBRELLA", line);
}
- if (strcmp(Buffer2, "3.0"))
+ if (std::strcmp(Buffer2, "3.0"))
{
gmx_fatal(FARGS, "This does not appear to be a version 3.0 pdo file");
}
}
ist.str(line);
ist >> Buffer3;
- if (strcmp(Buffer3, "#####") != 0)
+ if (std::strcmp(Buffer3, "#####") != 0)
{
gmx_fatal(FARGS, "Expected '#####', found %s. Hick.\n", Buffer3);
}
return NULL;
}
p = ptr;
- while ((strchr(ptr, '\n') == NULL) && (!feof(fp)))
+ while ((std::strchr(ptr, '\n') == NULL) && (!feof(fp)))
{
/* This line is longer than len characters, let's increase len! */
*len += STRLEN;
break;
}
}
- slen = strlen(ptr);
+ slen = std::strlen(ptr);
if (ptr[slen-1] == '\n')
{
ptr[slen-1] = '\0';
gmx_bool bGetMinMax, real *mintmp, real *maxtmp)
{
int i, inttemp, bins, count, ntot;
- real min, max, minfound = 1e20, maxfound = -1e20;
+ real minval, maxval, minfound = 1e20, maxfound = -1e20;
double temp, time, time0 = 0, dt;
char *ptr = 0;
t_UmbrellaWindow * window = 0;
if (!bGetMinMax)
{
- bins = opt->bins;
- min = opt->min;
- max = opt->max;
+ bins = opt->bins;
+ minval = opt->min;
+ maxval = opt->max;
window = win+fileno;
/* Need to alocate memory and set up structure */
{
minfound = 1e20;
maxfound = -1e20;
- min = max = bins = 0; /* Get rid of warnings */
+ minval = maxval = bins = 0; /* Get rid of warnings */
}
count = 0;
{
trim(ptr);
- if (ptr[0] == '#' || strlen(ptr) < 2)
+ if (ptr[0] == '#' || std::strlen(ptr) < 2)
{
continue;
}
/* Initiate format string */
fmtign[0] = '\0';
- strcat(fmtign, "%*s");
+ std::strcat(fmtign, "%*s");
sscanf(ptr, fmtlf, &time); /* printf("Time %f\n",time); */
/* Round time to fs */
{
for (i = 0; i < header->nPull; ++i)
{
- strcpy(fmt, fmtign);
- strcat(fmt, "%lf"); /* Creating a format stings such as "%*s...%*s%lf" */
- strcat(fmtign, "%*s"); /* ignoring one more entry in the next loop */
+ std::strcpy(fmt, fmtign);
+ std::strcat(fmt, "%lf"); /* Creating a format stings such as "%*s...%*s%lf" */
+ std::strcat(fmtign, "%*s"); /* ignoring one more entry in the next loop */
if (sscanf(ptr, fmt, &temp))
{
temp += header->UmbPos[i][0];
window->ztime[i][ntot] = temp;
}
- temp -= min;
- temp /= (max-min);
+ temp -= minval;
+ temp /= (maxval-minval);
temp *= bins;
- temp = floor(temp);
+ temp = std::floor(temp);
inttemp = static_cast<int> (temp);
if (opt->bCycl)
int jl, ju;
double pl, pu, dz, dp;
- jl = static_cast<int> (floor((dist-opt->tabMin)/opt->tabDz));
+ jl = static_cast<int> (std::floor((dist-opt->tabMin)/opt->tabDz));
ju = jl+1;
if (jl < 0 || ju >= opt->tabNbins)
{
{
U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
- contrib1 = profile[k]*exp(-U/(8.314e-3*opt->Temperature));
- contrib2 = window[i].N[j]*exp(-U/(8.314e-3*opt->Temperature) + window[i].z[j]);
+ contrib1 = profile[k]*std::exp(-U/(8.314e-3*opt->Temperature));
+ contrib2 = window[i].N[j]*std::exp(-U/(8.314e-3*opt->Temperature) + window[i].z[j]);
window[i].bContrib[j][k] = (contrib1 > wham_contrib_lim || contrib2 > wham_contrib_lim);
bAnyContrib = (bAnyContrib | window[i].bContrib[j][k]);
if (window[i].bContrib[j][k])
{
U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
- denom += invg*window[j].N[k]*exp(-U/(8.314e-3*opt->Temperature) + window[j].z[k]);
+ denom += invg*window[j].N[k]*std::exp(-U/(8.314e-3*opt->Temperature) + window[j].z[k]);
}
}
profile[i] = num/denom;
{
U = tabulated_pot(distance, opt); /* Use tabulated potential */
}
- total += profile[k]*exp(-U/(8.314e-3*opt->Temperature));
+ total += profile[k]*std::exp(-U/(8.314e-3*opt->Temperature));
}
/* Avoid floating point exception if window is far outside min and max */
if (total != 0.0)
{
- total = -log(total);
+ total = -std::log(total);
}
else
{
total = 1000.0;
}
- temp = fabs(total - window[i].z[j]);
+ temp = std::abs(total - window[i].z[j]);
if (temp > maxloc)
{
maxloc = temp;
z = min+(i+0.5)*dz;
zsym = -z;
/* bin left of zsym */
- j = static_cast<int> (floor((zsym-min)/dz-0.5));
+ j = static_cast<int> (std::floor((zsym-min)/dz-0.5));
if (j >= 0 && (j+1) < bins)
{
/* interpolate profile linearly between bins j and j+1 */
}
}
- memcpy(profile, prof2, bins*sizeof(double));
+ std::memcpy(profile, prof2, bins*sizeof(double));
sfree(prof2);
}
{
if (profile[i] > 0.0)
{
- profile[i] = -log(profile[i])*unit_factor;
+ profile[i] = -std::log(profile[i])*unit_factor;
}
}
if (opt->bs_verbose)
{
- snew(fn, strlen(fnhist)+10);
- snew(buf, strlen(fnhist)+10);
- sprintf(fn, "%s_cumul.xvg", strncpy(buf, fnhist, strlen(fnhist)-4));
+ snew(fn, std::strlen(fnhist)+10);
+ snew(buf, std::strlen(fnhist)+10);
+ sprintf(fn, "%s_cumul.xvg", std::strncpy(buf, fnhist, std::strlen(fnhist)-4));
fp = xvgropen(fn, "CDFs of umbrella windows", xlabel, "CDF", opt->oenv);
}
"cannot be predicted. You have 3 options:\n"
"1) Make g_wham estimate the ACTs (options -ac and -acsig).\n"
"2) Calculate the ACTs by yourself (e.g. with g_analyze) and provide them\n");
- strcat(errstr,
- " with option -iiact for all umbrella windows.\n"
- "3) If all ACTs are identical and know, you can define them with -bs-tau.\n"
- " Use option (3) only if you are sure what you're doing, you may severely\n"
- " underestimate the error if a too small ACT is given.\n");
+ std::strcat(errstr,
+ " with option -iiact for all umbrella windows.\n"
+ "3) If all ACTs are identical and know, you can define them with -bs-tau.\n"
+ " Use option (3) only if you are sure what you're doing, you may severely\n"
+ " underestimate the error if a too small ACT is given.\n");
gmx_fatal(FARGS, errstr);
}
Note: The ACT of the flat distribution and of the generated histogram is not
100% exactly tau, but near tau (my test was 3.8 instead of 4).
*/
- a = exp(-1.0/tausteps);
- ap = sqrt(1-a*a);
- invsqrt2 = 1./sqrt(2.0);
+ a = std::exp(-1.0/tausteps);
+ ap = std::sqrt(1.0-a*a);
+ invsqrt2 = 1.0/std::sqrt(2.0);
/* init random sequence */
x = gmx_rng_gaussian_table(opt->rng);
y = gmx_rng_gaussian_table(opt->rng);
x = a*x+ap*y;
z = x*sig+mu;
- ibin = static_cast<int> (floor((z-opt->min)/opt->dz));
+ ibin = static_cast<int> (std::floor((z-opt->min)/opt->dz));
if (opt->bCycl)
{
if (ibin < 0)
if (bs_index >= 0)
{
- snew(fn, strlen(fnhist)+10);
- snew(buf, strlen(fnhist)+1);
- sprintf(fn, "%s_bs%d.xvg", strncpy(buf, fnhist, strlen(fnhist)-4), bs_index);
+ snew(fn, std::strlen(fnhist)+10);
+ snew(buf, std::strlen(fnhist)+1);
+ sprintf(fn, "%s_bs%d.xvg", std::strncpy(buf, fnhist, std::strlen(fnhist)-4), bs_index);
sprintf(title, "Umbrella histograms. Bootstrap #%d", bs_index);
}
else
{
fn = gmx_strdup(fnhist);
- strcpy(title, "Umbrella histograms");
+ std::strcpy(title, "Umbrella histograms");
}
fp = xvgropen(fn, title, xlabel, "count", opt->oenv);
/* Write histograms */
for (l = 0; l < bins; ++l)
{
- fprintf(fp, "%e\t", (double)(l+0.5)*opt->dz+opt->min);
+ fprintf(fp, "%e\t", (l+0.5)*opt->dz+opt->min);
for (i = 0; i < nWindows; ++i)
{
for (j = 0; j < window[i].nPull; ++j)
i = 0;
bExact = FALSE;
maxchange = 1e20;
- memcpy(bsProfile, profile, opt->bins*sizeof(double)); /* use profile as guess */
+ std::memcpy(bsProfile, profile, opt->bins*sizeof(double)); /* use profile as guess */
do
{
if ( (i%opt->stepUpdateContrib) == 0)
bsProfiles_av [i] /= opt->nBootStrap;
bsProfiles_av2[i] /= opt->nBootStrap;
tmp = bsProfiles_av2[i]-sqr(bsProfiles_av[i]);
- stddev = (tmp >= 0.) ? sqrt(tmp) : 0.; /* Catch rouding errors */
+ stddev = (tmp >= 0.) ? std::sqrt(tmp) : 0.; /* Catch rouding errors */
fprintf(fp, "%e\t%e\t%e\n", (i+0.5)*opt->dz+opt->min, bsProfiles_av [i], stddev);
}
xvgrclose(fp);
int whaminFileType(char *fn)
{
int len;
- len = strlen(fn);
- if (strcmp(fn+len-3, "tpr") == 0)
+ len = std::strlen(fn);
+ if (std::strcmp(fn+len-3, "tpr") == 0)
{
return whamin_tpr;
}
- else if (strcmp(fn+len-3, "xvg") == 0 || strcmp(fn+len-6, "xvg.gz") == 0)
+ else if (std::strcmp(fn+len-3, "xvg") == 0 || std::strcmp(fn+len-6, "xvg.gz") == 0)
{
return whamin_pullxf;
}
- else if (strcmp(fn+len-3, "pdo") == 0 || strcmp(fn+len-6, "pdo.gz") == 0)
+ else if (std::strcmp(fn+len-3, "pdo") == 0 || std::strcmp(fn+len-6, "pdo.gz") == 0)
{
return whamin_pdo;
}
sizenow = 0;
while (fgets(tmp, sizeof(tmp), fp) != NULL)
{
- if (strlen(tmp) >= WHAM_MAXFILELEN)
+ if (std::strlen(tmp) >= WHAM_MAXFILELEN)
{
gmx_fatal(FARGS, "Filename too long in %s. Only %d characters allowed.\n", fn, WHAM_MAXFILELEN);
}
}
}
/* remove newline if there is one */
- if (tmp[strlen(tmp)-1] == '\n')
+ if (tmp[std::strlen(tmp)-1] == '\n')
{
- tmp[strlen(tmp)-1] = '\0';
+ tmp[std::strlen(tmp)-1] = '\0';
}
- strcpy(filename[nread], tmp);
+ std::strcpy(filename[nread], tmp);
if (opt->verbose)
{
printf("Found file %s in %s\n", filename[nread], fn);
static gmx_bool bFirst = 1;
/* gzipped pdo file? */
- if ((strcmp(fn+strlen(fn)-3, ".gz") == 0))
+ if ((std::strcmp(fn+std::strlen(fn)-3, ".gz") == 0))
{
/* search gunzip executable */
if (!(Path = getenv("GMX_PATH_GZIP")))
{
r2 += sqr(dx[i][line]);
}
- return sqrt(r2);
+ return std::sqrt(r2);
}
//! Read pullx.xvg or pullf.xvg
window->ztime[gUsed][ntot] = pos;
}
- ibin = static_cast<int> (floor((pos-min)/(max-min)*bins));
+ ibin = static_cast<int> (std::floor((pos-min)/(max-min)*bins));
if (opt->bCycl)
{
if (ibin < 0)
siglim = 3.0*opt->sigSmoothIact;
siglim2 = dsqr(siglim);
/* pre-factor of Gaussian */
- gaufact = 1.0/(sqrt(2*M_PI)*opt->sigSmoothIact);
+ gaufact = 1.0/(std::sqrt(2*M_PI)*opt->sigSmoothIact);
invtwosig2 = 0.5/dsqr(opt->sigSmoothIact);
for (i = 0; i < nwins; i++)
dpos2 = dsqr(window[j].pos[jg]-pos);
if (dpos2 < siglim2)
{
- w = gaufact*exp(-dpos2*invtwosig2);
+ w = gaufact*std::exp(-dpos2*invtwosig2);
weight += w;
tausm += w*window[j].tau[jg];
/*printf("Weight %g dpos2=%g pos=%g gaufact=%g invtwosig2=%g\n",
diff = ztime[k]-av;
sum2 += diff*diff;
}
- sig = sqrt(sum2/ntot);
+ sig = std::sqrt(sum2/ntot);
window[i].aver[ig] = av;
/* Note: This estimate for sigma is biased from the limited sampling.
for (ig = 0; ig < window[i].nPull; ig++)
{
hispos = window[i].pos[ig];
- dist = fabs(hispos-pos);
+ dist = std::abs(hispos-pos);
/* average force within bin */
if (dist < dz/2)
{
*/
for (j = 0; j < opt->bins; ++j)
{
- pot[j] = exp(-pot[j]/(8.314e-3*opt->Temperature));
+ pot[j] = std::exp(-pot[j]/(8.314e-3*opt->Temperature));
}
calc_z(pot, window, nWindows, opt, TRUE);
int i, n, is[2];
int cur = 0;
- if (strlen(ptr) == 0)
+ if (std::strlen(ptr) == 0)
{
return 0;
}
fmtign[0] = '\0';
for (i = 0; i < n; i++)
{
- strcpy(fmt, fmtign);
- strcat(fmt, "%d");
+ std::strcpy(fmt, fmtign);
+ std::strcat(fmt, "%d");
if (sscanf(ptr, fmt, &temp))
{
opt->groupsel[iline].bUse[i] = (temp > 0);
opt->groupsel[iline].nUse++;
}
}
- strcat(fmtign, "%*s");
+ std::strcat(fmtign, "%*s");
}
iline++;
}
xlabel, "count", opt.oenv);
for (l = 0; l < opt.bins; ++l)
{
- fprintf(histout, "%e\t", (double)(l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
+ fprintf(histout, "%e\t", (l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
for (i = 0; i < nwins; ++i)
{
for (j = 0; j < window[i].nPull; ++j)
if (opt.bLog)
{
prof_normalization_and_unit(profile, &opt);
- strcpy(ylabel, en_unit_label[opt.unit]);
- strcpy(title, "Umbrella potential");
+ std::strcpy(ylabel, en_unit_label[opt.unit]);
+ std::strcpy(title, "Umbrella potential");
}
else
{
- strcpy(ylabel, "Density of states");
- strcpy(title, "Density of states");
+ std::strcpy(ylabel, "Density of states");
+ std::strcpy(title, "Density of states");
}
/* symmetrize profile around z=0? */
profout = xvgropen(opt2fn("-o", NFILE, fnm), title, xlabel, ylabel, opt.oenv);
for (i = 0; i < opt.bins; ++i)
{
- fprintf(profout, "%e\t%e\n", (double)(i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
+ fprintf(profout, "%e\t%e\n", (i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
}
xvgrclose(profout);
printf("Wrote %s\n", opt2fn("-o", NFILE, fnm));
#include "nrama.h"
#include <cstdlib>
+#include <cstring>
#include "gromacs/listed-forces/bonded.h"
#include "gromacs/pbcutil/rmpbc.h"
for (i = start; (i < nr); i++)
{
- if (strcmp(find, *names[i]) == 0)
+ if (std::strcmp(find, *names[i]) == 0)
{
return i;
}
{
for (j = 0; j < pr->grn; j++)
{
- sq->s[i] += (pr->gr[j]/pr->r[j])*sin(sq->q[i]*pr->r[j]);
+ sq->s[i] += (pr->gr[j]/pr->r[j])*std::sin(sq->q[i]*pr->r[j]);
}
sq->s[i] /= sq->q[i];
}
{
for (j = 0; j < pr->grn; j++)
{
- sq->s[i] += (pr->gr[j]/pr->r[j])*sin(sq->q[i]*pr->r[j]);
+ sq->s[i] += (pr->gr[j]/pr->r[j])*std::sin(sq->q[i]*pr->r[j]);
}
sq->s[i] /= sq->q[i];
}
fphi = phi*sd->dx;
fpsi = psi*sd->dy;
- iphi = (int)fphi;
- ipsi = (int)fpsi;
+ iphi = static_cast<int>(fphi);
+ ipsi = static_cast<int>(fpsi);
fphi -= iphi; /* Fraction (offset from gridpoint) */
fpsi -= ipsi;
/* Sort eigenvalues in ascending order */
#define SWAPPER(i) \
- if (fabs(dd[i+1]) < fabs(dd[i])) { \
+ if (std::abs(dd[i+1]) < std::abs(dd[i])) { \
temp = dd[i]; \
for (j = 0; (j < NDIM); j++) { tvec[j] = ev[j][i]; } \
dd[i] = dd[i+1]; \
{
if (bQ)
{
- m0 = fabs(atom[ii].q);
+ m0 = std::abs(atom[ii].q);
}
else
{
#include "sfactor.h"
#include <cmath>
+#include <cstring>
#include <algorithm>
for (i = 0; (i < asize(uh)); i++)
{
- if (strcmp(name, uh[i].name) == 0)
+ if (std::strcmp(name, uh[i].name) == 0)
{
return NCMT-1+uh[i].nh;
}
for (i = 0; (i < NCMT); i++)
{
- if (strncmp (name, gsft->atomnm[i], strlen(gsft->atomnm[i])) == 0)
+ if (std::strncmp (name, gsft->atomnm[i], std::strlen(gsft->atomnm[i])) == 0)
{
tndx[cnt] = i;
cnt++;
nrc = 0;
for (i = 0; i < cnt; i++)
{
- if (strlen(gsft->atomnm[tndx[i]]) > (size_t)nrc)
+ if (std::strlen(gsft->atomnm[tndx[i]]) > (size_t)nrc)
{
- nrc = strlen(gsft->atomnm[tndx[i]]);
+ nrc = std::strlen(gsft->atomnm[tndx[i]]);
fndx = tndx[i];
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff