X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_hbond.cpp;fp=src%2Fgromacs%2Fgmxana%2Fgmx_hbond.c;h=a33a1e6bdb82967c47212582a1d572721e224164;hb=3a5934fccd489925ab5e80aeb4bf703040137071;hp=8bde9f1ac0724d8889a63ae5ef96dc37036d1acc;hpb=6ecd5d46bf94b6ef898e8fac4a7770bd7b2ee5d5;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_hbond.c b/src/gromacs/gmxana/gmx_hbond.cpp similarity index 96% rename from src/gromacs/gmxana/gmx_hbond.c rename to src/gromacs/gmxana/gmx_hbond.cpp index 8bde9f1ac0..a33a1e6bdb 100644 --- a/src/gromacs/gmxana/gmx_hbond.c +++ b/src/gromacs/gmxana/gmx_hbond.cpp @@ -38,7 +38,11 @@ #include "config.h" -#include +#include +#include +#include + +#include #include "gromacs/commandline/pargs.h" #include "gromacs/correlationfunctions/autocorr.h" @@ -50,6 +54,7 @@ #include "gromacs/fileio/trxio.h" #include "gromacs/fileio/xvgr.h" #include "gromacs/gmxana/gmx_ana.h" +#include "gromacs/gmxana/gstat.h" #include "gromacs/legacyheaders/copyrite.h" #include "gromacs/legacyheaders/macros.h" #include "gromacs/legacyheaders/txtdump.h" @@ -222,8 +227,6 @@ static void mk_hbmap(t_hbdata *hb) static void add_frames(t_hbdata *hb, int nframes) { - int i, j, k, l; - if (nframes >= hb->max_frames) { hb->max_frames += 4096; @@ -284,7 +287,7 @@ static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb) { int i, j, n; t_hbond *hb = hbd->hbmap[id][ia]; - int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]); + int maxhydro = std::min(hbd->maxhydro, hbd->d.nhydro[id]); int wlen = hbd->wordlen; int delta = 32*wlen; @@ -723,11 +726,10 @@ static void search_donors(t_topology *top, int isize, atom_id *index, t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable) { - int i, j, nra, n; + int i, j; t_functype func_type; t_ilist *interaction; atom_id nr1, nr2, nr3; - gmx_bool stop; if (!ddd->dptr) { @@ -865,7 +867,7 @@ static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid) { for (i = 0; i < DIM; i++) { - ngrid[i] = (box[i][i]/(1.2*rcut)); + ngrid[i] = static_cast(box[i][i]/(1.2*rcut)); } } @@ -917,7 +919,7 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, int i, m, gr, xi, yi, zi, nr; atom_id *ad; ivec grididx; - rvec invdelta, dshell, xtemp = {0, 0, 0}; + rvec invdelta, dshell; t_ncell *newgrid; gmx_bool bDoRshell, bInShell, bAcc; real rshell2 = 0; @@ -1035,7 +1037,7 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, for (m = DIM-1; m >= 0; m--) { /* determine grid index of atom */ - grididx[m] = x[ad[i]][m]*invdelta[m]; + grididx[m] = static_cast(x[ad[i]][m]*invdelta[m]); grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m]; } } @@ -1236,9 +1238,8 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, real *d_ha, real *ang, gmx_bool bDA, int *hhh, gmx_bool bContact, gmx_bool bMerge) { - int h, hh, id, ja, ihb; - rvec r_da, r_ha, r_dh, r = {0, 0, 0}; - ivec ri; + int h, hh, id; + rvec r_da, r_ha, r_dh; real rc2, r2c2, rda2, rha2, ca; gmx_bool HAinrange = FALSE; /* If !bDA. Needed for returning hbDist in a correct way. */ gmx_bool daSwap = FALSE; @@ -1249,7 +1250,7 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, } if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) || - ((ja = acceptor_index(&hb->a, grpa, a)) == NOTSET)) + (acceptor_index(&hb->a, grpa, a) == NOTSET)) { return hbNo; } @@ -1332,8 +1333,8 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, if (ca >= ccut) { *hhh = hh; - *d_ha = sqrt(bDA ? rda2 : rha2); - *ang = acos(ca); + *d_ha = std::sqrt(bDA ? rda2 : rha2); + *ang = std::acos(ca); return hbHB; } } @@ -1363,8 +1364,8 @@ static void do_merge(t_hbdata *hb, int ntmp, /* Decide where to start from when merging */ n00 = hb0->n0; n01 = hb1->n0; - nn0 = min(n00, n01); - nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0; + nn0 = std::min(n00, n01); + nnframes = std::max(n00 + hb0->nframes, n01 + hb1->nframes) - nn0; /* Initiate tmp arrays */ for (m = 0; (m < ntmp); m++) { @@ -1414,7 +1415,7 @@ static void do_merge(t_hbdata *hb, int ntmp, static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) { - int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp; + int i, inrnew, indnew, j, ii, jj, id, ia, ntmp; unsigned int *htmp, *gtmp; t_hbond *hb0, *hb1; @@ -1482,8 +1483,6 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t) { int i, j, k, n_bound[MAXHH], nbtot; - h_id nhb; - /* Set array to 0 */ for (k = 0; (k < MAXHH); k++) @@ -1692,7 +1691,7 @@ static real calc_dg(real tau, real temp) } else { - return kbt*log(kbt*tau/PLANCK); + return kbt*std::log(kbt*tau/PLANCK); } } @@ -1743,15 +1742,14 @@ static real compute_weighted_rates(int n, real t[], real ct[], real nt[], tl.ndelta = NK; for (j = 0; (j < NK); j++) { - /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */ kkk += tl.kkk[0]; kkp += tl.kkk[1]; kk2 += sqr(tl.kkk[0]); kp2 += sqr(tl.kkk[1]); tl.n0++; } - *sigma_k = sqrt(kk2/NK - sqr(kkk/NK)); - *sigma_kp = sqrt(kp2/NK - sqr(kkp/NK)); + *sigma_k = std::sqrt(kk2/NK - sqr(kkk/NK)); + *sigma_kp = std::sqrt(kp2/NK - sqr(kkp/NK)); return chi2; } @@ -1762,7 +1760,7 @@ void analyse_corr(int n, real t[], real ct[], real nt[], real kt[], { int i0, i; real k = 1, kp = 1, kow = 1; - real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg; + real Q = 0, chi2, tau_hb, dtau, tau_rlx, e_1, sigma_k, sigma_kp, ddg; double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0; gmx_bool bError = (sigma_ct != NULL) && (sigma_nt != NULL) && (sigma_kt != NULL); @@ -1794,9 +1792,9 @@ void analyse_corr(int n, real t[], real ct[], real nt[], real kt[], { chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]); } - chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt, - sigma_kt, &k, &kp, - &sigma_k, &sigma_kp, fit_start); + compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt, + sigma_kt, &k, &kp, + &sigma_k, &sigma_kp, fit_start); Q = 0; /* quality_of_fit(chi2, 2);*/ ddg = BOLTZ*temp*sigma_k/k; printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n", @@ -1843,7 +1841,7 @@ void analyse_corr(int n, real t[], real ct[], real nt[], real kt[], tau_hb = evaluate_integral(n, t, ct, NULL, (t[n-1]-t[0])/2, &dtau); printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb, bError ? " " : "", tau_hb, calc_dg(tau_hb, temp)); - e_1 = exp(-1); + e_1 = std::exp(-1.0); for (i = 0; (i < n-2); i++) { if ((ct[i] > e_1) && (ct[i+1] <= e_1)) @@ -1902,7 +1900,7 @@ static void normalizeACF(real *ct, real *gt, int nhb, int len) ct_fac = 1.0/ct[0]; if (nhb != 0) { - gt_fac = 1.0/(real)nhb; + gt_fac = 1.0/nhb; } printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac); @@ -1922,12 +1920,7 @@ static void do_hbac(const char *fn, t_hbdata *hb, int nThreads) { FILE *fp; - int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets; - const char *legNN[] = { - "Ac(t)", - "Ac'(t)" - }; - static char **legGem; + int i, j, k, m, ihb, idist, n2, nn; const char *legLuzar[] = { "Ac\\sfin sys\\v{}\\z{}(t)", @@ -1937,20 +1930,15 @@ static void do_hbac(const char *fn, t_hbdata *hb, }; gmx_bool bNorm = FALSE, bOMP = FALSE; double nhb = 0; - int nhbi = 0; real *rhbex = NULL, *ht, *gt, *ght, *dght, *kt; - real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct; + real *ct, tail, tail2, dtail, *cct; const real tol = 1e-3; - int nframes = hb->nframes, nf; + int nframes = hb->nframes; unsigned int **h = NULL, **g = NULL; - int nh, nhbonds, nhydro, ngh; + int nh, nhbonds, nhydro; t_hbond *hbh; - int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN; - gmx_bool c; int acType; - double *ctdouble, *timedouble, *fittedct; - double fittolerance = 0.1; - int *dondata = NULL, thisThread; + int *dondata = NULL; enum { AC_NONE, AC_NN, AC_GEM, AC_LUZAR @@ -1989,12 +1977,10 @@ static void do_hbac(const char *fn, t_hbdata *hb, /* Total number of hbonds analyzed here */ nhbonds = 0; - ngh = 0; - anhb = 0; if (acType != AC_LUZAR && bOMP) { - nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads()); + nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads()); gmx_omp_set_num_threads(nThreads); snew(dondata, nThreads); @@ -2060,7 +2046,7 @@ static void do_hbac(const char *fn, t_hbdata *hb, } } - nf = hbh->nframes; + int nf = hbh->nframes; for (nh = 0; (nh < nhydro); nh++) { int nrint = bContact ? hb->nrdist : hb->nrhb; @@ -2120,7 +2106,7 @@ static void do_hbac(const char *fn, t_hbdata *hb, fprintf(stderr, "\n"); sfree(h); sfree(g); - normalizeACF(ct, ght, nhb, nn); + normalizeACF(ct, ght, static_cast(nhb), nn); /* Determine tail value for statistics */ tail = 0; @@ -2132,7 +2118,7 @@ static void do_hbac(const char *fn, t_hbdata *hb, } tail /= (nn - nn/2); tail2 /= (nn - nn/2); - dtail = sqrt(tail2-tail*tail); + dtail = std::sqrt(tail2-tail*tail); /* Check whether the ACF is long enough */ if (dtail > tol) @@ -2189,7 +2175,7 @@ static void do_hbac(const char *fn, t_hbdata *hb, static void init_hbframe(t_hbdata *hb, int nframes, real t) { - int i, j, m; + int i; hb->time[nframes] = t; hb->nhb[nframes] = 0; @@ -2362,7 +2348,6 @@ static void dump_hbmap(t_hbdata *hb, * It mimics add_frames() and init_frame() to some extent. */ static void sync_hbdata(t_hbdata *p_hb, int nframes) { - int i; if (nframes >= p_hb->max_frames) { p_hb->max_frames += 4096; @@ -2374,22 +2359,15 @@ static void sync_hbdata(t_hbdata *p_hb, int nframes) { srenew(p_hb->danr, p_hb->max_frames); } - memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); - memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); + std::memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); + std::memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); p_hb->nhb[nframes] = 0; p_hb->ndist[nframes] = 0; } p_hb->nframes = nframes; -/* for (i=0;) */ -/* { */ -/* p_hb->nhx[nframes][i] */ -/* } */ - memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */ - - /* hb->per will remain constant througout the frame loop, - * even though the data its members point to will change, - * hence no need for re-syncing. */ + + std::memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */ } int gmx_hbond(int argc, char *argv[]) @@ -2465,14 +2443,12 @@ int gmx_hbond(int argc, char *argv[]) }; static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1; - static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, D = -1; + static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15; static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE; - static int nDump = 0, nFitPoints = 100; - static int nThreads = 0, nBalExp = 4; + static int nDump = 0; + static int nThreads = 0; static gmx_bool bContact = FALSE; - static real logAfterTime = 10, gemBallistic = 0.2; /* ps */ - static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL}; /* options */ t_pargs pa [] = { @@ -2553,23 +2529,21 @@ int gmx_hbond(int argc, char *argv[]) matrix box; real t, ccut, dist = 0.0, ang = 0.0; double max_nhb, aver_nhb, aver_dist; - int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel; + int h = 0, i = 0, j, k = 0, ogrp, nsel; int xi, yi, zi, ai; int xj, yj, zj, aj, xjj, yjj, zjj; - int xk, yk, zk, ak, xkk, ykk, zkk; - gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric; - int *adist, *rdist, *aptr, *rprt; - int grp, nabin, nrbin, bin, resdist, ihb; + gmx_bool bSelected, bHBmap, bStop, bTwo, bBox, bTric; + int *adist, *rdist; + int grp, nabin, nrbin, resdist, ihb; char **leg; - t_hbdata *hb, *hbptr; - FILE *fp, *fpins = NULL, *fpnhb = NULL, *donor_properties = NULL; + t_hbdata *hb; + FILE *fp, *fpnhb = NULL, *donor_properties = NULL; t_gridcell ***grid; - t_ncell *icell, *jcell, *kcell; + t_ncell *icell, *jcell; ivec ngrid; unsigned char *datable; output_env_t oenv; - int NN; - int ii, jj, hh, actual_nThreads; + int ii, hh, actual_nThreads; int threadNr = 0; gmx_bool bParallel; gmx_bool bEdge_yjj, bEdge_xjj, bOMP; @@ -2594,7 +2568,7 @@ int gmx_hbond(int argc, char *argv[]) /* process input */ bSelected = FALSE; - ccut = cos(acut*DEG2RAD); + ccut = std::cos(acut*DEG2RAD); if (bContact) { @@ -2812,13 +2786,11 @@ int gmx_hbond(int argc, char *argv[]) bBox = ir.ePBC != epbcNONE; grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid); - nabin = acut/abin; - nrbin = rcut/rbin; + nabin = static_cast(acut/abin); + nrbin = static_cast(rcut/rbin); snew(adist, nabin+1); snew(rdist, nrbin+1); - bParallel = FALSE; - #ifndef GMX_OPENMP #define __ADIST adist #define __RDIST rdist @@ -2838,7 +2810,7 @@ int gmx_hbond(int argc, char *argv[]) if (bParallel) { - actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads()); + actual_nThreads = std::min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads()); gmx_omp_set_num_threads(actual_nThreads); printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads); @@ -2886,12 +2858,12 @@ int gmx_hbond(int argc, char *argv[]) #pragma omp parallel \ firstprivate(i) \ - private(j, h, ii, jj, hh, \ + private(j, h, ii, hh, \ xi, yi, zi, xj, yj, zj, threadNr, \ dist, ang, icell, jcell, \ grp, ogrp, ai, aj, xjj, yjj, zjj, \ - xk, yk, zk, ihb, id, resdist, \ - xkk, ykk, zkk, kcell, ak, k, bTric, \ + ihb, resdist, \ + k, bTric, \ bEdge_xjj, bEdge_yjj) \ default(shared) { /* Start of parallel region */ @@ -3033,21 +3005,19 @@ int gmx_hbond(int argc, char *argv[]) gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist); } ang *= RAD2DEG; - __ADIST[(int)( ang/abin)]++; - __RDIST[(int)(dist/rbin)]++; + __ADIST[static_cast( ang/abin)]++; + __RDIST[static_cast(dist/rbin)]++; if (!bTwo) { - int id, ia; - if ((id = donor_index(&hb->d, grp, i)) == NOTSET) + if (donor_index(&hb->d, grp, i) == NOTSET) { gmx_fatal(FARGS, "Invalid donor %d", i); } - if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET) + if (acceptor_index(&hb->a, ogrp, j) == NOTSET) { gmx_fatal(FARGS, "Invalid acceptor %d", j); } - resdist = abs(top.atoms.atom[i].resind- - top.atoms.atom[j].resind); + resdist = std::abs(top.atoms.atom[i].resind-top.atoms.atom[j].resind); if (resdist >= max_hx) { resdist = max_hx-1; @@ -3238,7 +3208,7 @@ int gmx_hbond(int argc, char *argv[]) /* Print HB distance distribution */ if (opt2bSet("-dist", NFILE, fnm)) { - long sum; + int sum; sum = 0; for (i = 0; i < nrbin; i++) @@ -3253,7 +3223,7 @@ int gmx_hbond(int argc, char *argv[]) "Hydrogen - Acceptor Distance (nm)", "", oenv); for (i = 0; i < nrbin; i++) { - fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum)); + fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*sum)); } xvgrclose(fp); } @@ -3274,7 +3244,7 @@ int gmx_hbond(int argc, char *argv[]) "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv); for (i = 0; i < nabin; i++) { - fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum)); + fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*sum)); } xvgrclose(fp); } @@ -3308,16 +3278,12 @@ int gmx_hbond(int argc, char *argv[]) Added support for -contact in ac and hbm calculations below. - Erik Marklund, May 29, 2006 */ - ivec itmp; - rvec rtmp; if (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)) { please_cite(stdout, "Spoel2006b"); } if (opt2bSet("-ac", NFILE, fnm)) { - char *gemstring = NULL; - do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump, bMerge, bContact, fit_start, temp, r2cut > 0, oenv, nThreads);