biod.pnpi.spb.ru
/
alexxy
/
gromacs.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Merge branch 'master' into pygromacs
[alexxy/gromacs.git]
/
src
/
gromacs
/
gmxana
/
gmx_analyze.cpp
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 5d1c8c4a59d4a91a5955140c6e329d85a4415bde..c38228272764f6d0e18daf39be9e4767c2fa1b8f 100644
(file)
--- a/
src/gromacs/gmxana/gmx_analyze.c
+++ b/
src/gromacs/gmxana/gmx_analyze.cpp
@@
-36,9
+36,9
@@
*/
#include "gmxpre.h"
*/
#include "gmxpre.h"
-#include <
math.
h>
-#include <
stdlib.h
>
-#include <
string.h
>
+#include <
cmat
h>
+#include <
cstdlib
>
+#include <
cstring
>
#include "gromacs/commandline/pargs.h"
#include "gromacs/correlationfunctions/autocorr.h"
#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)
{
{
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++)
{
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++)
{
}
}
for (s = 0; s < nset; s++)
{
- i = 0;
for (i = 0; i < n && val[s][i] >= 0; i++)
{
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)
{
}
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",
}
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);
}
sfree(y);
@@
-130,7
+129,7
@@
static real cosine_content(int nhp, int n, real *y)
yyint = 0;
for (i = 0; i < n; i++)
{
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];
}
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;
{
FILE *fp;
int i, s;
- double min
, max
;
+ double min
val, maxval
;
int nbin;
gmx_int64_t *histo;
int nbin;
gmx_int64_t *histo;
- min = val[0][0];
- max = val[0][0];
+ min
val
= val[0][0];
+ max
val
= val[0][0];
for (s = 0; s < nset; s++)
{
for (i = 0; i < n; i++)
{
for (s = 0; s < nset; s++)
{
for (i = 0; i < n; i++)
{
- if (val[s][i] < min)
+ if (val[s][i] < min
val
)
{
{
- min = val[s][i];
+ min
val
= val[s][i];
}
}
- else if (val[s][i] > max)
+ else if (val[s][i] > max
val
)
{
{
- max = val[s][i];
+ max
val
= val[s][i];
}
}
}
}
}
}
- min
= binwidth*floor(min
/binwidth);
- max
= binwidth*ceil(max
/binwidth);
- if (min != 0)
+ min
val = binwidth*std::floor(minval
/binwidth);
+ max
val = binwidth*std::ceil(maxval
/binwidth);
+ if (min
val
!= 0)
{
{
- min -= binwidth;
+ min
val
-= binwidth;
}
}
- max += binwidth;
+ max
val
+= 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);
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++)
{
}
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++)
{
}
for (i = 0; i < nbin; i++)
{
- fprintf(fp, " %g %g\n", min
+i*binwidth, (double)histo[i]
/(n*binwidth));
+ fprintf(fp, " %g %g\n", min
val+i*binwidth, static_cast<double>(histo[i])
/(n*binwidth));
}
if (s < nset-1)
{
}
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");
{
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%%"
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
{
}
else
{
@@
-375,11
+374,11
@@
static void average(const char *avfile, int avbar_opt,
}
if (avbar_opt == avbarSTDDEV)
{
}
if (avbar_opt == avbarSTDDEV)
{
- err = sqrt(var/nset);
+ err = s
td::s
qrt(var/nset);
}
else
{
}
else
{
- err = sqrt(var/(nset*(nset-1)));
+ err = s
td::s
qrt(var/(nset*(nset-1)));
}
fprintf(fp, " %g", err);
}
}
fprintf(fp, " %g", err);
}
@@
-408,7
+407,7
@@
static real optimal_error_estimate(double sigma, double fitparm[], real tTotal)
return 0;
}
return 0;
}
- return sigma*sqrt(2*ss/tTotal);
+ return sigma*s
td::s
qrt(2*ss/tTotal);
}
static void estimate_error(const char *eefile, int nb_min, int resol, int n,
}
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);
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);
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)
{
nbr = nb_min;
while (nbr <= n)
{
- bs = n/
(int)nbr
;
+ bs = n/
static_cast<int>(nbr)
;
if (bs != prev_bs)
{
nb = n/bs;
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;
nbs++;
}
nbr *= spacing;
- nb = (int)(nbr+0.5);
prev_bs = bs;
}
if (sig[s] == 0)
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.
*/
* 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
{
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);
}
{
dens = 0.5*(tbs[i+1]/tbs[i-1] - 1);
}
- fitsig[i] = sqrt((tau_sig + tbs[i])/dens);
+ fitsig[i] = s
td::s
qrt((tau_sig + tbs[i])/dens);
}
if (!bSingleExpFit)
}
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).
*/
/* 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] = s
td::s
qrt(tau1_est*(n-1)*dt);
do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv,
bDebugMode(), effnERREST, fitparm, 0,
NULL);
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++)
{
}
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]*s
td::s
qrt(ybs[i]/(n*dt)),
+ sig[s]*s
td::s
qrt(fit_function(effnERREST, fitparm, tbs[i])/(n*dt)));
}
if (bFitAc)
}
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)
{
ac[i] = val[s][i] - av[s];
if (i > 0)
{
- fitsig[i] = s
qrt(i
);
+ fitsig[i] = s
td::sqrt(static_cast<real>(i)
);
}
else
{
}
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++)
{
/* Generate more or less appropriate sigma's */
for (i = 0; i <= fitlen; i++)
{
- fitsig[i] = sqrt(acint + dt*i);
+ fitsig[i] = s
td::s
qrt(acint + dt*i);
}
ac_fit[0] = 0.5*acint;
}
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],
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]*s
td::s
qrt(fit_function(effnERREST, ac_fit, tbs[i]))/(n*dt));
}
sfree(ac);
}
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)
{
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 *kt;
- real
weight,
d2;
+ real d2;
int j;
please_cite(stdout, "Spoel2006b");
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]);
}
{
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", s
td::s
qrt(d2/nn));
}
analyse_corr(nn, time, val[0], val[2], kt, NULL, NULL, NULL, fit_start,
temp);
}
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;
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++)
{
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++)
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;
}
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, s
td::s
qrt(fluc));
}
}
- fprintf(stdout, "Overall filtered fluctuation: %12.6e\n", sqrt(fluctot/nset));
+ fprintf(stdout, "Overall filtered fluctuation: %12.6e\n", s
td::s
qrt(fluctot/nset));
fprintf(stdout, "\n");
sfree(filt);
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)
{
int s, buflen = 0;
if (NULL != fn_fitted)
{
- buflen = strlen(fn_fitted)+32;
+ buflen = st
d::st
rlen(fn_fitted)+32;
snew(buf2, buflen);
snew(buf2, buflen);
- strncpy(buf2, fn_fitted, buflen);
- buf2[strlen(buf2)-4] = '\0';
+ st
d::st
rncpy(buf2, fn_fitted, buflen);
+ buf2[st
d::st
rlen(buf2)-4] = '\0';
}
for (s = 0; s < nset; s++)
{
}
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 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] = {
/* 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)
};
#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;
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;
cum3 /= n;
cum4 /= n;
av[s] = cum1;
- sig[s] = sqrt(cum2);
+ sig[s] = s
td::s
qrt(cum2);
if (n > 1)
{
if (n > 1)
{
- error = sqrt(cum2/(n-1));
+ error = s
td::s
qrt(cum2/(n-1));
}
else
{
}
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,
}
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]*s
td::s
qrt(8/M_PI)) : 0,
sig[s] ? cum4/(sig[s]*sig[s]*sig[s]*sig[s]*3)-1 : 0);
}
printf("\n");
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);
{
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++)
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 += 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)
fprintf(out, " %g %8g\n", dt*j, tot);
}
if (s < nset-1)