#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"
#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
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<real> > data; /* the displacement data. First index is the group
+ number, second is frame number */
+ std::vector<real> time; /* frame time */
+ std::vector<real> mass; /* masses for mass-weighted msd */
+ matrix **datam;
+ std::vector< std::vector<gmx::RVec> > x0; /* original positions */
+ std::vector<gmx::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 */
+ std::vector<int> n_offs;
+ std::vector< std::vector<int> > 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),
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<real>()),
datam(nullptr),
- x0(nullptr),
- com(nullptr),
lsq(nullptr),
type(static_cast<msd_type>(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<int>())
{
- 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;
}
~t_corr()
{
- sfree(ndata);
- sfree(data);
for (int i = 0; i < nrestart; i++)
{
for (int j = 0; j < nmol; j++)
}
}
sfree(lsq);
- if (mass)
- {
- sfree(mass);
- }
}
};
{
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++;
}
}
-/* 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)
{
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;
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)
{
}
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;
*/
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<const int> gnx_com, int *index_com[],
real dt, real t_pdb, rvec **x_pdb, matrix box_pdb,
const gmx_output_env_t *oenv)
{
#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);
}
// 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);
{
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++)
{
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);
}
}
}
- srenew(curr->time, maxframes);
+ curr->time.resize(maxframes);
}
/* set the time */
}
/* 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);
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;
real dt, real beginfit, real endfit, const gmx_output_env_t *oenv)
{
std::unique_ptr<t_corr> msd;
- int *gnx; /* the selected groups' sizes */
- int **index; /* selected groups' indices */
+ std::vector<int> 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<real> 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)
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);
top->atoms.nr = i;
}
- DD = nullptr;
- SigmaD = nullptr;
-
if (beginfit == -1)
{
i0 = gmx::roundToInt(0.1*(msd->nframes - 1));
}
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)
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",
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[])
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff