*/
#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;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff