X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_current.cpp;fp=src%2Fgromacs%2Fgmxana%2Fgmx_current.c;h=833e5a5ef306b27f89c145ac744b531a9532cf76;hb=3a5934fccd489925ab5e80aeb4bf703040137071;hp=107d87223055f1bb335e96f62d8c06d7623d18bc;hpb=6ecd5d46bf94b6ef898e8fac4a7770bd7b2ee5d5;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_current.c b/src/gromacs/gmxana/gmx_current.cpp similarity index 94% rename from src/gromacs/gmxana/gmx_current.c rename to src/gromacs/gmxana/gmx_current.cpp index 107d872230..833e5a5ef3 100644 --- a/src/gromacs/gmxana/gmx_current.c +++ b/src/gromacs/gmxana/gmx_current.cpp @@ -34,8 +34,8 @@ */ #include "gmxpre.h" -#include -#include +#include +#include #include "gromacs/commandline/pargs.h" #include "gromacs/fileio/confio.h" @@ -52,9 +52,9 @@ #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) @@ -139,7 +139,7 @@ static gmx_bool precalc(t_topology top, real mass2[], real qmol[]) } } - 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; @@ -275,12 +275,10 @@ static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int { 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) { @@ -288,8 +286,9 @@ static real calc_cacf(FILE *fcacf, real prefactor, real cacf[], real time[], int while (i < nfr) { + real corint; - rfr = (real) (nfr/nshift-i/nshift); + rfr = static_cast(nfr+i)/nshift; cacf[i] /= rfr; if (time[vfr[i]] <= time[vfr[ei]]) @@ -327,14 +326,9 @@ static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], i { 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++) { @@ -353,18 +347,17 @@ static void calc_mjdsp(FILE *fmjdsp, real prefactor, real dsp2[], real time[], 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; @@ -374,7 +367,6 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor, 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; @@ -392,12 +384,9 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor, 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 @@ -449,7 +438,7 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor, do { - refr = (real)(nfr+1); + refr = (nfr+1); if (nfr >= nalloc) { @@ -470,8 +459,7 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor, xshfr[i] = 0.0; } } - assert(time != NULL); - + GMX_RELEASE_ASSERT(time != NULL, "Memory not allocated correctly - time array is NULL"); if (nfr == 0) { @@ -652,9 +640,6 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor, if (v0 != NULL) { - trust *= (real) nvfr; - itrust = (int) trust; - if (bINT) { @@ -678,15 +663,15 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor, 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; @@ -729,7 +714,7 @@ static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor, 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) @@ -792,7 +777,6 @@ int gmx_current(int argc, char *argv[]) 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; @@ -814,8 +798,6 @@ int gmx_current(int argc, char *argv[]) "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."} }; @@ -834,21 +816,13 @@ int gmx_current(int argc, char *argv[]) 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; @@ -915,7 +889,7 @@ int gmx_current(int argc, char *argv[]) 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); @@ -933,7 +907,7 @@ int gmx_current(int argc, char *argv[]) snew(mass2, top.atoms.nr); snew(qmol, top.atoms.nr); - bNEU = precalc(top, mass2, qmol); + precalc(top, mass2, qmol); snew(indexm, isize); @@ -990,7 +964,7 @@ int gmx_current(int argc, char *argv[]) * 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);