X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_analyze.cpp;fp=src%2Fgromacs%2Fgmxana%2Fgmx_analyze.c;h=c38228272764f6d0e18daf39be9e4767c2fa1b8f;hb=c3f2d46e4047f0c465f7234b3784a2fa6f02a065;hp=5d1c8c4a59d4a91a5955140c6e329d85a4415bde;hpb=0595b4a4c763a0bc574658992081abf8b0abc3fe;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_analyze.c b/src/gromacs/gmxana/gmx_analyze.cpp similarity index 94% rename from src/gromacs/gmxana/gmx_analyze.c rename to src/gromacs/gmxana/gmx_analyze.cpp index 5d1c8c4a59..c382282727 100644 --- a/src/gromacs/gmxana/gmx_analyze.c +++ b/src/gromacs/gmxana/gmx_analyze.cpp @@ -36,9 +36,9 @@ */ #include "gmxpre.h" -#include -#include -#include +#include +#include +#include #include "gromacs/commandline/pargs.h" #include "gromacs/correlationfunctions/autocorr.h" @@ -80,7 +80,7 @@ static void power_fit(int n, int nset, real **val, real *t) { if (t[0] > 0) { - x[i] = log(t[i]); + x[i] = std::log(t[i]); } } } @@ -89,16 +89,15 @@ static void power_fit(int n, int nset, real **val, real *t) 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(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) { @@ -106,7 +105,7 @@ static void power_fit(int n, int nset, real **val, real *t) } 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); @@ -130,7 +129,7 @@ static real cosine_content(int nhp, int n, real *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]; } @@ -244,36 +243,36 @@ void histogram(const char *distfile, real binwidth, int n, int nset, real **val, { 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(((maxval - minval)/binwidth + 0.5) + 1); fprintf(stderr, "Making distributions with %d bins\n", nbin); snew(histo, nbin); fp = xvgropen(distfile, "Distribution", "", "", oenv); @@ -285,11 +284,11 @@ void histogram(const char *distfile, real binwidth, int n, int nset, real **val, } for (i = 0; i < n; i++) { - histo[(int)((val[s][i] - min)/binwidth + 0.5)]++; + histo[static_cast((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(histo[i])/(n*binwidth)); } if (s < nset-1) { @@ -336,9 +335,9 @@ static void average(const char *avfile, int avbar_opt, { snew(tmp, nset); fprintf(fp, "@TYPE xydydy\n"); - edge = (int)(nset*0.05+0.5); + edge = static_cast(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(100*(nset-2*edge)/nset+0.5)); } else { @@ -375,11 +374,11 @@ static void average(const char *avfile, int avbar_opt, } 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); } @@ -408,7 +407,7 @@ static real optimal_error_estimate(double sigma, double fitparm[], real tTotal) 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, @@ -445,7 +444,7 @@ 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); @@ -456,7 +455,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, nbr = nb_min; while (nbr <= n) { - bs = n/(int)nbr; + bs = n/static_cast(nbr); if (bs != prev_bs) { nb = n/bs; @@ -482,7 +481,6 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, nbs++; } nbr *= spacing; - nb = (int)(nbr+0.5); prev_bs = bs; } if (sig[s] == 0) @@ -507,7 +505,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, * 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 { @@ -553,7 +551,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, { 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) @@ -563,7 +561,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, /* 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); @@ -637,8 +635,8 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, } 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) @@ -653,7 +651,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, ac[i] = val[s][i] - av[s]; if (i > 0) { - fitsig[i] = sqrt(i); + fitsig[i] = std::sqrt(static_cast(i)); } else { @@ -677,7 +675,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, /* 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; @@ -695,7 +693,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, 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); @@ -714,9 +712,8 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, 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"); @@ -737,7 +734,7 @@ static void luzar_correl(int nn, real *time, int nset, real **val, real temp, { 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); @@ -760,13 +757,13 @@ static void filter(real flen, int n, int nset, real **val, real dt) int f, s, i, j; double *filt, sum, vf, fluc, fluctot; - f = (int)(flen/(2*dt)); + f = static_cast(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++) @@ -790,9 +787,9 @@ static void filter(real flen, int n, int nset, real **val, real dt) } 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); @@ -926,10 +923,10 @@ static void print_fitted_function(const char *fitfile, 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++) { @@ -1039,9 +1036,9 @@ int gmx_analyze(int argc, char *argv[]) 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] = { @@ -1106,7 +1103,7 @@ int gmx_analyze(int argc, char *argv[]) }; #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; @@ -1232,10 +1229,10 @@ int gmx_analyze(int argc, char *argv[]) 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 { @@ -1243,7 +1240,7 @@ int gmx_analyze(int argc, char *argv[]) } 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"); @@ -1257,7 +1254,7 @@ int gmx_analyze(int argc, char *argv[]) { out = xvgropen(msdfile, "Mean square displacement", "time", "MSD (nm\\S2\\N)", oenv); - nlast = (int)(n*frac); + nlast = static_cast(n*frac); for (s = 0; s < nset; s++) { for (j = 0; j <= nlast; j++) @@ -1271,7 +1268,7 @@ int gmx_analyze(int argc, char *argv[]) { 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)