X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_energy.cpp;fp=src%2Fgromacs%2Fgmxana%2Fgmx_energy.c;h=b95200a5cf9b5eb0dac0a2485f1781ab043fc8e1;hb=c3f2d46e4047f0c465f7234b3784a2fa6f02a065;hp=fc2b29bbf882e5de776018d9b4364b41f0bb6438;hpb=0595b4a4c763a0bc574658992081abf8b0abc3fe;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_energy.c b/src/gromacs/gmxana/gmx_energy.cpp similarity index 95% rename from src/gromacs/gmxana/gmx_energy.c rename to src/gromacs/gmxana/gmx_energy.cpp index fc2b29bbf8..b95200a5cf 100644 --- a/src/gromacs/gmxana/gmx_energy.c +++ b/src/gromacs/gmxana/gmx_energy.cpp @@ -36,9 +36,11 @@ */ #include "gmxpre.h" -#include -#include -#include +#include +#include +#include + +#include #include "gromacs/commandline/pargs.h" #include "gromacs/correlationfunctions/autocorr.h" @@ -60,6 +62,7 @@ #include "gromacs/topology/mtop_util.h" #include "gromacs/utility/cstringutil.h" #include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/gmxassert.h" #include "gromacs/utility/smalloc.h" static real minthird = -1.0/3.0, minsixth = -1.0/6.0; @@ -94,7 +97,7 @@ static double mypow(double x, double y) { if (x > 0) { - return pow(x, y); + return std::pow(x, y); } else { @@ -159,7 +162,7 @@ static int *select_it(int nre, char *nm[], int *nset) static void chomp(char *buf) { - int len = strlen(buf); + int len = std::strlen(buf); while ((len > 0) && (buf[len-1] == '\n')) { @@ -171,7 +174,7 @@ static void chomp(char *buf) static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) { gmx_bool *bE; - int n, k, kk, j, i, nmatch, nind, nss; + int k, kk, j, i, nmatch, nind, nss; int *set; gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE; char *ptr, buf[STRLEN]; @@ -196,7 +199,7 @@ static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) { newnm[k] = gmx_strdup(nm[k].name); /* Insert dashes in all the names */ - while ((ptr = strchr(newnm[k], ' ')) != NULL) + while ((ptr = std::strchr(newnm[k], ' ')) != NULL) { *ptr = '-'; } @@ -211,7 +214,7 @@ static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) bLong = FALSE; for (kk = k; kk < k+4; kk++) { - if (kk < nre && strlen(nm[kk].name) > 14) + if (kk < nre && std::strlen(nm[kk].name) > 14) { bLong = TRUE; } @@ -258,7 +261,7 @@ static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) trim(buf); /* Empty line means end of input */ - bEOF = (strlen(buf) == 0); + bEOF = (std::strlen(buf) == 0); if (!bEOF) { ptr = buf; @@ -287,7 +290,6 @@ static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) else { /* Now try to read a string */ - i = strlen(ptr); nmatch = 0; for (nind = 0; nind < nre; nind++) { @@ -299,7 +301,7 @@ static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) } if (nmatch == 0) { - i = strlen(ptr); + i = std::strlen(ptr); nmatch = 0; for (nind = 0; nind < nre; nind++) { @@ -317,12 +319,12 @@ static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset) } } /* Look for the first space, and remove spaces from there */ - if ((ptr = strchr(ptr, ' ')) != NULL) + if ((ptr = std::strchr(ptr, ' ')) != NULL) { trim(ptr); } } - while (!bEOF && (ptr && (strlen(ptr) > 0))); + while (!bEOF && (ptr && (std::strlen(ptr) > 0))); } } @@ -355,7 +357,6 @@ static void get_dhdl_parms(const char *topnm, t_inputrec *ir) { gmx_mtop_t mtop; int natoms; - t_iatom *iatom; matrix box; /* all we need is the ir to be able to write the label */ @@ -508,8 +509,8 @@ static void calc_violations(real rt[], real rav3[], int nb, int index[], rav += sqr(rav3[j]); rsum += mypow(rt[j], -6); } - rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]); - rav = max(0.0, mypow(rav, -sixth)-bounds[i]); + rsum = std::max(0.0, mypow(rsum, -sixth)-bounds[i]); + rav = std::max(0.0, mypow(rav, -sixth)-bounds[i]); sumt += rsum; sumaver += rav; @@ -549,10 +550,10 @@ static void analyse_disre(const char *voutfn, int nframes, { sumaver += sqr(violaver[j]/nframes); } - sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]); + sumaver = std::max(0.0, mypow(sumaver, minsixth)-bounds[i]); sumt += sumaver; - sum = max(sum, sumaver); + sum = std::max(sum, sumaver); fprintf(vout, "%10d %10.5e\n", index[i], sumaver); } #ifdef DEBUG @@ -692,7 +693,7 @@ static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax) { int nb, i, f, nee; double sum, sum2, sump, see2; - gmx_int64_t steps, np, p, bound_nb; + gmx_int64_t np, p, bound_nb; enerdat_t *ed; exactsum_t *es; gmx_bool bAllZero; @@ -819,11 +820,11 @@ static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax) edat->s[i].av = sum/np; if (ed->bExactStat) { - edat->s[i].rmsd = sqrt(sum2/np); + edat->s[i].rmsd = std::sqrt(sum2/np); } else { - edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av)); + edat->s[i].rmsd = std::sqrt(sum2/np - dsqr(edat->s[i].av)); } if (edat->nframes > 1) @@ -858,7 +859,7 @@ static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax) } if (nee > 0) { - edat->s[i].ee = sqrt(see2/nee); + edat->s[i].ee = std::sqrt(see2/nee); } else { @@ -940,7 +941,7 @@ static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *ed /* Remove the drift by performing a fit to y = ax+b. Uses 5 iterations. */ int i, j, k; - double delta, d, sd, sd2; + double delta; edat->npoints = edat->nframes; edat->nsteps = edat->nframes; @@ -974,7 +975,7 @@ static void calc_fluctuation_props(FILE *fp, int nbmin, int nbmax) { int i, j; - double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet; + double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, varet; double NANO3; enum { eVol, eEnth, eTemp, eEtot, eNR @@ -1029,13 +1030,11 @@ static void calc_fluctuation_props(FILE *fp, cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt)); } /* Total energy */ - et = varet = NOTSET; if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET)) { /* Only compute cv in constant volume runs, which we can test by checking whether the enthalpy was computed. */ - et = edat->s[ii[eEtot]].av; varet = sqr(edat->s[ii[eEtot]].rmsd); cv = KILO*((varet/nmol)/(BOLTZ*tt*tt)); } @@ -1142,14 +1141,12 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, /* Check out the printed manual for equations! */ double Dt, aver, stddev, errest, delta_t, totaldrift; enerdata_t *esum = NULL; - real xxx, integral, intBulk, Temp = 0, Pres = 0; - real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest; + real integral, intBulk, Temp = 0, Pres = 0; + real pr_aver, pr_stddev, pr_errest; double beta = 0, expE, expEtot, *fee = NULL; gmx_int64_t nsteps; int nexact, nnotexact; - double x1m, x1mk; - int i, j, k, nout; - real chi2; + int i, j, nout; char buf[256], eebuf[100]; nsteps = step - start_step + 1; @@ -1252,24 +1249,24 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, expE = 0; for (j = 0; (j < edat->nframes); j++) { - expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol); + expE += std::exp(beta*(edat->s[i].ener[j] - aver)/nmol); } if (bSum) { expEtot += expE/edat->nframes; } - fee[i] = log(expE/edat->nframes)/beta + aver/nmol; + fee[i] = std::log(expE/edat->nframes)/beta + aver/nmol; } - if (strstr(leg[i], "empera") != NULL) + if (std::strstr(leg[i], "empera") != NULL) { Temp = aver; } - else if (strstr(leg[i], "olum") != NULL) + else if (std::strstr(leg[i], "olum") != NULL) { Vaver = aver; } - else if (strstr(leg[i], "essure") != NULL) + else if (std::strstr(leg[i], "essure") != NULL) { Pres = aver; } @@ -1320,7 +1317,7 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, if (bFee) { fprintf(stdout, " %10g %10g\n", - log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta); + std::log(expEtot)/beta + esum->s[0].av/nmol, std::log(expEtot)/beta); } else { @@ -1393,13 +1390,13 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/ /* Do it for shear viscosity */ - strcpy(buf, "Shear Viscosity"); + std::strcpy(buf, "Shear Viscosity"); low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3, (edat->nframes+1)/2, eneset, Dt, eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0); /* Now for bulk viscosity */ - strcpy(buf, "Bulk Viscosity"); + std::strcpy(buf, "Bulk Viscosity"); low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1, (edat->nframes+1)/2, &(eneset[11]), Dt, eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0); @@ -1434,11 +1431,11 @@ static void analyse_ener(gmx_bool bCorr, const char *corrfn, { if (bFluct) { - strcpy(buf, "Autocorrelation of Energy Fluctuations"); + std::strcpy(buf, "Autocorrelation of Energy Fluctuations"); } else { - strcpy(buf, "Energy Autocorrelation"); + std::strcpy(buf, "Energy Autocorrelation"); } #if 0 do_autocorr(corrfn, oenv, buf, edat->nframes, @@ -1478,7 +1475,7 @@ static void fec(const char *ene2fn, const char *runavgfn, }; FILE *fp; ener_file_t enx; - int nre, timecheck, step, nenergy, nenergy2, maxenergy; + int timecheck, nenergy, nenergy2, maxenergy; int i, j; gmx_bool bCont; real aver, beta; @@ -1528,6 +1525,8 @@ static void fec(const char *ene2fn, const char *runavgfn, srenew(eneset2[i], maxenergy); } } + GMX_RELEASE_ASSERT(time != NULL, "trying to dereference NULL time pointer"); + if (fr->t != time[nenergy2]) { fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n", @@ -1549,7 +1548,7 @@ static void fec(const char *ene2fn, const char *runavgfn, fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n", edat->nframes, nenergy2); } - nenergy = min(edat->nframes, nenergy2); + nenergy = std::min(edat->nframes, nenergy2); /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */ fp = NULL; @@ -1573,14 +1572,14 @@ static void fec(const char *ene2fn, const char *runavgfn, for (j = 0; j < nenergy; j++) { dE = eneset2[i][j] - edat->s[i].ener[j]; - sum += exp(-dE*beta); + sum += std::exp(-dE*beta); if (fp) { fprintf(fp, "%10g %10g %10g\n", - time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) ); + time[j], dE, -BOLTZ*reftemp*std::log(sum/(j+1)) ); } } - aver = -BOLTZ*reftemp*log(sum/nenergy); + aver = -BOLTZ*reftemp*std::log(sum/nenergy); fprintf(stdout, "%-24s %10g\n", leg[i], aver); } if (fp) @@ -1599,17 +1598,14 @@ static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda"; char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN]; char buf[STRLEN]; - gmx_bool first = FALSE; int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0; int i, j, k; /* coll data */ - double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0; + double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0; static int setnr = 0; double *native_lambda_vec = NULL; const char **lambda_components = NULL; int n_lambda_vec = 0; - gmx_bool changing_lambda = FALSE; - int lambda_fep_state; /* now count the blocks & handle the global dh data */ for (i = 0; i < fr->nblock; i++) @@ -1633,15 +1629,12 @@ static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, } /* read the data from the DHCOLL block */ - temp = fr->block[i].sub[0].dval[0]; - start_time = fr->block[i].sub[0].dval[1]; - delta_time = fr->block[i].sub[0].dval[2]; + temp = fr->block[i].sub[0].dval[0]; + start_time = fr->block[i].sub[0].dval[1]; + delta_time = fr->block[i].sub[0].dval[2]; start_lambda = fr->block[i].sub[0].dval[3]; - delta_lambda = fr->block[i].sub[0].dval[4]; - changing_lambda = (delta_lambda != 0); if (fr->block[i].nsub > 1) { - lambda_fep_state = fr->block[i].sub[1].ival[0]; if (n_lambda_vec == 0) { n_lambda_vec = fr->block[i].sub[1].ival[1]; @@ -1775,14 +1768,12 @@ static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, } } } - (*samples) += (int)(sum/nblock_hist); + (*samples) += static_cast(sum/nblock_hist); } else { /* raw dh */ int len = 0; - char **setnames = NULL; - int nnames = nblock_dh; for (i = 0; i < fr->nblock; i++) { @@ -1828,7 +1819,7 @@ static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, if (j == 1 && ir->bExpanded) { - fprintf(*fp_dhdl, "%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */ + fprintf(*fp_dhdl, "%4d", static_cast(value)); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */ } else { @@ -1997,13 +1988,11 @@ int gmx_energy(int argc, char *argv[]) FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL; FILE *fp_dhdl = NULL; - FILE **drout; ener_file_t fp; int timecheck = 0; gmx_mtop_t mtop; gmx_localtop_t *top = NULL; t_inputrec ir; - t_energy **ee; enerdata_t edat; gmx_enxnm_t *enm = NULL; t_enxframe *frame, *fr = NULL; @@ -2017,15 +2006,14 @@ int gmx_energy(int argc, char *argv[]) int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL; int nbounds = 0, npairs; gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL; - gmx_bool bFoundStart, bCont, bEDR, bVisco; - double sum, sumaver, sumt, ener, dbl; + gmx_bool bFoundStart, bCont, bVisco; + double sum, sumaver, sumt, dbl; double *time = NULL; real Vaver; int *set = NULL, i, j, k, nset, sss; gmx_bool *bIsEner = NULL; char **pairleg, **odtleg, **otenleg; char **leg = NULL; - char **nms; char *anm_j, *anm_k, *resnm_j, *resnm_k; int resnr_j, resnr_k; const char *orinst_sub = "@ subtitle \"instantaneous\"\n"; @@ -2099,7 +2087,7 @@ int gmx_energy(int argc, char *argv[]) { for (i = 0; i < nre; i++) { - if (strstr(enm[i].name, setnm[j])) + if (std::strstr(enm[i].name, setnm[j])) { set[j] = i; break; @@ -2134,16 +2122,16 @@ int gmx_energy(int argc, char *argv[]) { for (j = 0; j < i; j++) { - if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0) + if (std::strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0) { break; } } if (j == i) { - strcat(buf, ", ("); - strcat(buf, enm[set[i]].unit); - strcat(buf, ")"); + std::strcat(buf, ", ("); + std::strcat(buf, enm[set[i]].unit); + std::strcat(buf, ")"); } } out = xvgropen(opt2fn("-o", NFILE, fnm), "GROMACS Energies", "Time (ps)", buf, @@ -2375,20 +2363,20 @@ int gmx_energy(int argc, char *argv[]) if (edat.nframes % 1000 == 0) { srenew(edat.step, edat.nframes+1000); - memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0])); + std::memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0])); srenew(edat.steps, edat.nframes+1000); - memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0])); + std::memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0])); srenew(edat.points, edat.nframes+1000); - memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0])); + std::memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0])); for (i = 0; i < nset; i++) { srenew(edat.s[i].ener, edat.nframes+1000); - memset(&(edat.s[i].ener[edat.nframes]), 0, - 1000*sizeof(edat.s[i].ener[0])); + std::memset(&(edat.s[i].ener[edat.nframes]), 0, + 1000*sizeof(edat.s[i].ener[0])); srenew(edat.s[i].es, edat.nframes+1000); - memset(&(edat.s[i].es[edat.nframes]), 0, - 1000*sizeof(edat.s[i].es[0])); + std::memset(&(edat.s[i].es[edat.nframes]), 0, + 1000*sizeof(edat.s[i].es[0])); } } @@ -2531,6 +2519,7 @@ int gmx_energy(int argc, char *argv[]) *******************************************/ if (ndisre > 0) { + GMX_RELEASE_ASSERT(blk_disre != NULL, "Trying to dereference NULL blk_disre pointer"); #ifndef GMX_DOUBLE float *disre_rt = blk_disre->sub[0].fval; float *disre_rm3tav = blk_disre->sub[1].fval; @@ -2586,7 +2575,7 @@ int gmx_energy(int argc, char *argv[]) print_time(out, fr->t); print1(out, bDp, fr->ener[set[0]].e); print1(out, bDp, fr->ener[set[0]].esum/fr->nsum); - print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum)); + print1(out, bDp, std::sqrt(fr->ener[set[0]].eav/fr->nsum)); fprintf(out, "\n"); } } @@ -2639,7 +2628,7 @@ int gmx_energy(int argc, char *argv[]) vals = blk->sub[0].dval; #endif - if (blk->sub[0].nr != (size_t)nor) + if (blk->sub[0].nr != nor) { gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr); } @@ -2697,7 +2686,7 @@ int gmx_energy(int argc, char *argv[]) vals = blk->sub[0].dval; #endif - if (blk->sub[0].nr != (size_t)(nex*12)) + if (blk->sub[0].nr != nex*12) { gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)", blk->sub[0].nr/12, nex); @@ -2780,7 +2769,7 @@ int gmx_energy(int argc, char *argv[]) } for (i = 0; i < nor; i++) { - fprintf(out, "%5d %g\n", or_label[i], sqrt(odrms[i]/norfr)); + fprintf(out, "%5d %g\n", or_label[i], std::sqrt(odrms[i]/norfr)); } xvgrclose(out); }