From: Kevin Boyd Date: Sun, 7 Oct 2018 16:05:32 +0000 (-0400) Subject: Partial conversion of gmx msd to std::vector X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=a2833204d1227efa63eed022c1ae8d8ca6db4c62;p=alexxy%2Fgromacs.git Partial conversion of gmx msd to std::vector Changed raw arrays in the primary struct to vector, with the exception of Matrix - an STL friendly version of this is still in the gerrit pipeline refs #2368 Change-Id: Ic41e6a1f78d53c2d5134e5831536a714b0a38a5b --- diff --git a/src/gromacs/gmxana/gmx_msd.cpp b/src/gromacs/gmxana/gmx_msd.cpp index d81ce39ef9..8cbd28a23e 100644 --- a/src/gromacs/gmxana/gmx_msd.cpp +++ b/src/gromacs/gmxana/gmx_msd.cpp @@ -50,6 +50,7 @@ #include "gromacs/math/functions.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" +#include "gromacs/math/vectypes.h" #include "gromacs/pbcutil/rmpbc.h" #include "gromacs/statistics/statistics.h" #include "gromacs/topology/index.h" @@ -60,7 +61,7 @@ #include "gromacs/utility/gmxassert.h" #include "gromacs/utility/smalloc.h" -#define FACTOR 1000.0 /* Convert nm^2/ps to 10e-5 cm^2/s */ +static constexpr double diffusionConversionFactor = 1000.0; /* Convert nm^2/ps to 10e-5 cm^2/s */ /* NORMAL = total diffusion coefficient (default). X,Y,Z is diffusion coefficient in X,Y,Z direction. LATERAL is diffusion coefficient in plane perpendicular to axis @@ -69,33 +70,34 @@ typedef enum { NOT_USED, NORMAL, X, Y, Z, LATERAL } msd_type; +// TODO : Group related fields into a struct struct t_corr { - real t0; /* start time and time increment between */ - real delta_t; /* time between restart points */ - real beginfit, /* the begin/end time for fits as reals between */ - endfit; /* 0 and 1 */ - real dim_factor; /* the dimensionality factor for the diffusion - constant */ - real **data; /* the displacement data. First index is the group - number, second is frame number */ - real *time; /* frame time */ - real *mass; /* masses for mass-weighted msd */ - matrix **datam; - rvec **x0; /* original positions */ - rvec *com; /* center of mass correction for each frame */ - gmx_stats_t **lsq; /* fitting stats for individual molecule msds */ - msd_type type; /* the type of msd to calculate (lateral, etc.)*/ - int axis; /* the axis along which to calculate */ - int ncoords; - int nrestart; /* number of restart points */ - int nmol; /* number of molecules (for bMol) */ - int nframes; /* number of frames */ - int nlast; - int ngrp; /* number of groups to use for msd calculation */ - int *n_offs; - int **ndata; /* the number of msds (particles/mols) per data - point. */ - t_corr(int nrgrp, int type, int axis, real dim_factor, int nmol, + real t0; /* start time and time increment between */ + real delta_t; /* time between restart points */ + real beginfit, /* the begin/end time for fits as reals between */ + endfit; /* 0 and 1 */ + real dim_factor; /* the dimensionality factor for the diffusion + constant */ + std::vector< std::vector > data; /* the displacement data. First index is the group + number, second is frame number */ + std::vector time; /* frame time */ + std::vector mass; /* masses for mass-weighted msd */ + matrix **datam; + std::vector< std::vector > x0; /* original positions */ + std::vector com; /* center of mass correction for each frame */ + gmx_stats_t **lsq; /* fitting stats for individual molecule msds */ + msd_type type; /* the type of msd to calculate (lateral, etc.)*/ + int axis; /* the axis along which to calculate */ + int ncoords; + int nrestart; /* number of restart points */ + int nmol; /* number of molecules (for bMol) */ + int nframes; /* number of frames */ + int nlast; + int ngrp; /* number of groups to use for msd calculation */ + std::vector n_offs; + std::vector< std::vector > ndata; /* the number of msds (particles/mols) per data + point. */ + t_corr(int nrgrp, int type, int axis, real dim_factor, int nrmol, gmx_bool bTen, gmx_bool bMass, real dt, const t_topology *top, real beginfit, real endfit) : t0(0), @@ -103,53 +105,39 @@ struct t_corr { beginfit((1 - 2*GMX_REAL_EPS)*beginfit), endfit((1 + 2*GMX_REAL_EPS)*endfit), dim_factor(dim_factor), - data(nullptr), - time(nullptr), - mass(nullptr), + data(nrgrp, std::vector()), datam(nullptr), - x0(nullptr), - com(nullptr), lsq(nullptr), type(static_cast(type)), axis(axis), ncoords(0), nrestart(0), - nmol(nmol), + nmol(nrmol), nframes(0), nlast(0), ngrp(nrgrp), - n_offs(nullptr), - ndata(nullptr) + ndata(nrgrp, std::vector()) { - snew(ndata, nrgrp); - snew(data, nrgrp); + if (bTen) { snew(datam, nrgrp); - } - for (int i = 0; (i < nrgrp); i++) - { - ndata[i] = nullptr; - data[i] = nullptr; - if (bTen) + for (int i = 0; i < nrgrp; i++) { datam[i] = nullptr; } } + if (nmol > 0) { - snew(mass, nmol); - for (int i = 0; i < nmol; i++) - { - mass[i] = 1; - } + mass.resize(nmol, 1); } else { if (bMass) { const t_atoms *atoms = &top->atoms; - snew(mass, atoms->nr); + mass.resize(atoms->nr); for (int i = 0; (i < atoms->nr); i++) { mass[i] = atoms->atom[i].m; @@ -159,8 +147,6 @@ struct t_corr { } ~t_corr() { - sfree(ndata); - sfree(data); for (int i = 0; i < nrestart; i++) { for (int j = 0; j < nmol; j++) @@ -169,10 +155,6 @@ struct t_corr { } } sfree(lsq); - if (mass) - { - sfree(mass); - } } }; @@ -247,7 +229,7 @@ static void calc_corr(t_corr *curr, int nr, int nx, int index[], rvec xc[], { if ((thistime(curr) >= (curr->nlast*curr->delta_t)) && (nr == 0)) { - std::memcpy(curr->x0[curr->nlast], xc, curr->ncoords*sizeof(xc[0])); + std::memcpy(curr->x0[curr->nlast].data()->as_vec(), xc, curr->ncoords*sizeof(xc[0])); curr->n_offs[curr->nlast] = curr->nframes; copy_rvec(com, curr->com[curr->nlast]); curr->nlast++; @@ -281,7 +263,7 @@ static void calc_corr(t_corr *curr, int nr, int nx, int index[], rvec xc[], } } -/* the non-mass-weighted mean-squared displacement calcuation */ +/* the non-mass-weighted mean-squared displacement calculation */ static real calc1_norm(t_corr *curr, int nx, const int index[], int nx0, rvec xc[], const rvec dcom, gmx_bool bTen, matrix mat) { @@ -565,7 +547,6 @@ static void printmol(t_corr *curr, const char *fn, const char *fn_pdb, const int *molindex, const t_topology *top, rvec *x, int ePBC, matrix box, const gmx_output_env_t *oenv) { -#define NDIST 100 FILE *out; gmx_stats_t lsq1; int i, j; @@ -573,7 +554,7 @@ static void printmol(t_corr *curr, const char *fn, t_pdbinfo *pdbinfo = nullptr; const int *mol2a = nullptr; - out = xvgropen(fn, "Diffusion Coefficients / Molecule", "Molecule", "D", oenv); + out = xvgropen(fn, "Diffusion Coefficients / Molecule", "Molecule", "D (1e-5 cm^2/s)", oenv); if (fn_pdb) { @@ -597,7 +578,7 @@ static void printmol(t_corr *curr, const char *fn, } gmx_stats_get_ab(lsq1, elsqWEIGHT_NONE, &a, &b, nullptr, nullptr, nullptr, nullptr); gmx_stats_free(lsq1); - D = a*FACTOR/curr->dim_factor; + D = a*diffusionConversionFactor/curr->dim_factor; if (D < 0) { D = 0; @@ -654,7 +635,7 @@ static void printmol(t_corr *curr, const char *fn, */ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int ePBC, gmx_bool bMol, int gnx[], int *index[], - t_calc_func *calc1, gmx_bool bTen, int *gnx_com, int *index_com[], + t_calc_func *calc1, gmx_bool bTen, gmx::ArrayRef gnx_com, int *index_com[], real dt, real t_pdb, rvec **x_pdb, matrix box_pdb, const gmx_output_env_t *oenv) { @@ -673,7 +654,7 @@ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int eP #ifdef DEBUG fprintf(stderr, "Read %d atoms for first frame\n", natoms); #endif - if ((gnx_com != nullptr) && natoms < top->atoms.nr) + if ((!gnx_com.empty()) && natoms < top->atoms.nr) { fprintf(stderr, "WARNING: The trajectory only contains part of the system (%d of %d atoms) and therefore the COM motion of only this part of the system will be removed\n", natoms, top->atoms.nr); } @@ -682,7 +663,7 @@ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int eP // if com is requested, the data structure needs to be large enough to do this // to prevent overflow - if (bMol && !gnx_com) + if (bMol && gnx_com.empty()) { curr->ncoords = curr->nmol; snew(xa[0], curr->ncoords); @@ -732,10 +713,10 @@ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int eP { curr->nrestart++; - srenew(curr->x0, curr->nrestart); - snew(curr->x0[curr->nrestart-1], curr->ncoords); - srenew(curr->com, curr->nrestart); - srenew(curr->n_offs, curr->nrestart); + curr->x0.resize(curr->nrestart); + curr->x0[curr->nrestart-1].resize(curr->ncoords); + curr->com.resize(curr->nrestart); + curr->n_offs.resize(curr->nrestart); srenew(curr->lsq, curr->nrestart); snew(curr->lsq[curr->nrestart-1], curr->nmol); for (i = 0; i < curr->nmol; i++) @@ -756,20 +737,17 @@ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int eP { for (i = 0; (i < curr->ngrp); i++) { - curr->ndata[i] = nullptr; - curr->data[i] = nullptr; if (bTen) { curr->datam[i] = nullptr; } } - curr->time = nullptr; } maxframes += 10; for (i = 0; (i < curr->ngrp); i++) { - srenew(curr->ndata[i], maxframes); - srenew(curr->data[i], maxframes); + curr->ndata[i].resize(maxframes); + curr->data[i].resize(maxframes); if (bTen) { srenew(curr->datam[i], maxframes); @@ -784,7 +762,7 @@ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int eP } } } - srenew(curr->time, maxframes); + curr->time.resize(maxframes); } /* set the time */ @@ -819,7 +797,7 @@ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int eP } /* calculate the center of mass */ - if (gnx_com) + if (!gnx_com.empty()) { calc_com(bMol, gnx_com[0], index_com[0], xa[cur], xa[prev], box, &top->atoms, com); @@ -829,7 +807,7 @@ static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int eP for (i = 0; (i < curr->ngrp); i++) { /* calculate something useful, like mean square displacements */ - calc_corr(curr, i, gnx[i], index[i], xa[cur], (gnx_com != nullptr), com, + calc_corr(curr, i, gnx[i], index[i], xa[cur], (!gnx_com.empty()), com, calc1, bTen); } cur = prev; @@ -896,32 +874,32 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_ real dt, real beginfit, real endfit, const gmx_output_env_t *oenv) { std::unique_ptr msd; - int *gnx; /* the selected groups' sizes */ - int **index; /* selected groups' indices */ + std::vector gnx, gnx_com; /* the selected groups' sizes */ + int **index; /* selected groups' indices */ char **grpname; int i, i0, i1, j, N, nat_trx; - real *DD, *SigmaD, a, a2, b, r, chi2; + std::vector SigmaD, DD; + real a, a2, b, r, chi2; rvec *x = nullptr; matrix box; - int *gnx_com = nullptr; /* the COM removal group size */ int **index_com = nullptr; /* the COM removal group atom indices */ char **grpname_com = nullptr; /* the COM removal group name */ - snew(gnx, nrgrp); + gnx.resize(nrgrp); snew(index, nrgrp); snew(grpname, nrgrp); fprintf(stderr, "\nSelect a group to calculate mean squared displacement for:\n"); - get_index(&top->atoms, ndx_file, nrgrp, gnx, index, grpname); + get_index(&top->atoms, ndx_file, nrgrp, gnx.data(), index, grpname); if (bRmCOMM) { - snew(gnx_com, 1); + gnx_com.resize(1); snew(index_com, 1); snew(grpname_com, 1); fprintf(stderr, "\nNow select a group for center of mass removal:\n"); - get_index(&top->atoms, ndx_file, 1, gnx_com, index_com, grpname_com); + get_index(&top->atoms, ndx_file, 1, gnx_com.data(), index_com, grpname_com); } if (mol_file) @@ -934,7 +912,7 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_ bTen, bMW, dt, top, beginfit, endfit); nat_trx = - corr_loop(msd.get(), trx_file, top, ePBC, mol_file ? gnx[0] != 0 : false, gnx, index, + corr_loop(msd.get(), trx_file, top, ePBC, mol_file ? gnx[0] != 0 : false, gnx.data(), index, (mol_file != nullptr) ? calc1_mol : (bMW ? calc1_mw : calc1_norm), bTen, gnx_com, index_com, dt, t_pdb, pdb_file ? &x : nullptr, box, oenv); @@ -969,9 +947,6 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_ top->atoms.nr = i; } - DD = nullptr; - SigmaD = nullptr; - if (beginfit == -1) { i0 = gmx::roundToInt(0.1*(msd->nframes - 1)); @@ -1008,8 +983,8 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_ } else { - snew(DD, msd->ngrp); - snew(SigmaD, msd->ngrp); + DD.resize(msd->ngrp); + SigmaD.resize(msd->ngrp); for (j = 0; j < msd->ngrp; j++) { if (N >= 4) @@ -1023,8 +998,8 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_ SigmaD[j] = 0; } lsq_y_ax_b(N, &(msd->time[i0]), &(msd->data[j][i0]), &(DD[j]), &b, &r, &chi2); - DD[j] *= FACTOR/msd->dim_factor; - SigmaD[j] *= FACTOR/msd->dim_factor; + DD[j] *= diffusionConversionFactor/msd->dim_factor; + SigmaD[j] *= diffusionConversionFactor/msd->dim_factor; if (DD[j] > 0.01 && DD[j] < 1e4) { fprintf(stdout, "D[%10s] %.4f (+/- %.4f) 1e-5 cm^2/s\n", @@ -1041,7 +1016,7 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_ corr_print(msd.get(), bTen, msd_file, "Mean Square Displacement", "MSD (nm\\S2\\N)", - msd->time[msd->nframes-1], beginfit, endfit, DD, SigmaD, grpname, oenv); + msd->time[msd->nframes-1], beginfit, endfit, DD.data(), SigmaD.data(), grpname, oenv); } int gmx_msd(int argc, char *argv[]) diff --git a/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolMassWeighted.xml b/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolMassWeighted.xml index dfbab14e37..9f4b165984 100644 --- a/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolMassWeighted.xml +++ b/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolMassWeighted.xml @@ -7,7 +7,7 @@ diff --git a/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolNonMassWeighted.xml b/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolNonMassWeighted.xml index dfbab14e37..9f4b165984 100644 --- a/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolNonMassWeighted.xml +++ b/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolNonMassWeighted.xml @@ -7,7 +7,7 @@ diff --git a/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolSelected.xml b/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolSelected.xml index d315daf378..4c05484b73 100644 --- a/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolSelected.xml +++ b/src/gromacs/gmxana/tests/refdata/MsdMolTest_diffMolSelected.xml @@ -7,7 +7,7 @@