X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_anaeig.cpp;fp=src%2Fgromacs%2Fgmxana%2Fgmx_anaeig.c;h=2421ca402d65697507db74297cbd40705fe22298;hb=c3f2d46e4047f0c465f7234b3784a2fa6f02a065;hp=08441f75603d56c40d577ba27323d4bb782a83b9;hpb=0595b4a4c763a0bc574658992081abf8b0abc3fe;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_anaeig.c b/src/gromacs/gmxana/gmx_anaeig.cpp similarity index 91% rename from src/gromacs/gmxana/gmx_anaeig.c rename to src/gromacs/gmxana/gmx_anaeig.cpp index 08441f7560..2421ca402d 100644 --- a/src/gromacs/gmxana/gmx_anaeig.c +++ b/src/gromacs/gmxana/gmx_anaeig.cpp @@ -36,9 +36,11 @@ */ #include "gmxpre.h" -#include -#include -#include +#include +#include +#include + +#include #include "gromacs/commandline/pargs.h" #include "gromacs/fileio/confio.h" @@ -61,12 +63,13 @@ #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); @@ -74,10 +77,11 @@ static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip { 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) { @@ -88,7 +92,6 @@ static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip 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", @@ -99,14 +102,12 @@ static void calc_entropy_schlitter(FILE *fp, int n, int nskip, 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); @@ -116,7 +117,7 @@ static void calc_entropy_schlitter(FILE *fp, int n, int nskip, for (i = 0; (i < n-nskip); i++) { dd = 1+kteh*eigval[i]; - deter += log(dd); + deter += std::log(dd); } S = 0.5*RGAS*deter; @@ -134,7 +135,7 @@ static real tick_spacing(real range, int minticks) 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; @@ -152,7 +153,7 @@ static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph, { 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) @@ -163,50 +164,50 @@ static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph, { 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); @@ -230,8 +231,8 @@ static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph, { 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"); @@ -241,12 +242,12 @@ static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph, 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"); @@ -256,7 +257,7 @@ static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph, { 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) ? "&" : ""); } @@ -273,13 +274,13 @@ static void 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; @@ -290,7 +291,7 @@ compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2, 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++) @@ -305,9 +306,18 @@ compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2, 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]; @@ -317,10 +327,10 @@ compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2, 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(100*sum1/trace1+0.5), static_cast(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++) @@ -345,13 +355,13 @@ compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2, } 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)); } @@ -364,7 +374,7 @@ static void inprod_matrix(const char *matfile, int natoms, 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); @@ -396,7 +406,7 @@ static void inprod_matrix(const char *matfile, int natoms, snew(mat, nx); snew(t_x, nx); - max = 0; + maxval = 0; for (x1 = 0; x1 < nx; x1++) { snew(mat[x1], ny); @@ -429,10 +439,10 @@ static void inprod_matrix(const char *matfile, int natoms, { 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]; } } } @@ -442,7 +452,7 @@ static void inprod_matrix(const char *matfile, int natoms, 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); } @@ -508,7 +518,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb 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; @@ -642,6 +652,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb 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++) { @@ -665,7 +676,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb 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) ? "&" : ""); } @@ -729,7 +740,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb 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(std::ceil(i*fact)); atoms.resinfo[i].ic = ' '; x[i][XX] = inprod[0][i]; x[i][YY] = inprod[1][i]; @@ -741,6 +752,8 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb } 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) @@ -750,7 +763,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb 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; @@ -779,6 +792,7 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb 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"); @@ -810,16 +824,16 @@ static void project(const char *trajfile, t_topology *top, int ePBC, matrix topb 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 { @@ -906,7 +920,7 @@ static void rmsf(const char *outfile, int natoms, real *sqrtm, 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; @@ -939,7 +953,7 @@ static void rmsf(const char *outfile, int natoms, real *sqrtm, 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, @@ -1053,24 +1067,22 @@ int gmx_anaeig(int argc, char *argv[]) }; #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; @@ -1129,7 +1141,6 @@ int gmx_anaeig(int argc, char *argv[]) OverlapFile = opt2fn_null("-over", NFILE, fnm); InpMatFile = ftp2fn_null(efXPM, NFILE, fnm); - bTop = fn2bTPX(topfile); bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile || FilterFile || ExtremeFile; bFirstLastSet = @@ -1204,6 +1215,11 @@ int gmx_anaeig(int argc, char *argv[]) gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match"); } } + else + { + nvec2 = 0; + neig2 = 0; + } if (Eig2File != NULL) { @@ -1314,7 +1330,7 @@ int gmx_anaeig(int argc, char *argv[]) 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 @@ -1339,7 +1355,7 @@ int gmx_anaeig(int argc, char *argv[]) } } 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) @@ -1362,7 +1378,7 @@ int gmx_anaeig(int argc, char *argv[]) { /* 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;