Also, please use the syntax :issue:`number` to reference issues on GitLab, without the
a space between the colon and number!
+gmx msd has been migrated to the trajectoryanalysis framework
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The tool now uses the |Gromacs| selection syntax. Rather than piping selections via stdin,
+selections are now made using the "-sel" option.
+
+This migration comes with about a 20% speedup in execution time.
+
+TODO: Modify/Delete this segment as features are added back in.
+Some rarely used features have yet to be migrated, including:
+
+- Mass weighting of MSDs cannot currently be turned on or off. It is set to on when -mol is set, otherwise off.
+- The -tensor option is not yet implemented.
+- System COM removal with -rmcomm has not yet been implemented.
+- B-factor writing using the -pdb option is not yet supported.
+
+:issue:`2368`
\ No newline at end of file
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2008, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
int gmx_mindist(int argc, char* argv[]);
-int gmx_msd(int argc, char* argv[]);
-
int gmx_nmeig(int argc, char* argv[]);
int gmx_nmens(int argc, char* argv[]);
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include <cmath>
-#include <cstring>
-
-#include <memory>
-
-#include "gromacs/commandline/pargs.h"
-#include "gromacs/commandline/viewit.h"
-#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/trxio.h"
-#include "gromacs/fileio/xvgr.h"
-#include "gromacs/gmxana/gmx_ana.h"
-#include "gromacs/gmxana/gstat.h"
-#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/topology/topology.h"
-#include "gromacs/utility/arraysize.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/futil.h"
-#include "gromacs/utility/gmxassert.h"
-#include "gromacs/utility/smalloc.h"
-
-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
- */
-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 */
- 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),
- delta_t(dt),
- beginfit((1 - 2 * GMX_REAL_EPS) * beginfit),
- endfit((1 + 2 * GMX_REAL_EPS) * endfit),
- dim_factor(dim_factor),
- data(nrgrp, std::vector<real>()),
- datam(nullptr),
- lsq(nullptr),
- type(static_cast<msd_type>(type)),
- axis(axis),
- ncoords(0),
- nrestart(0),
- nmol(nrmol),
- nframes(0),
- nlast(0),
- ngrp(nrgrp),
- ndata(nrgrp, std::vector<int>())
- {
-
- if (bTen)
- {
- snew(datam, nrgrp);
- for (int i = 0; i < nrgrp; i++)
- {
- datam[i] = nullptr;
- }
- }
-
- if (nmol > 0)
- {
- mass.resize(nmol, 1);
- }
- else
- {
- if (bMass)
- {
- const t_atoms* atoms = &top->atoms;
- mass.resize(atoms->nr);
- for (int i = 0; (i < atoms->nr); i++)
- {
- mass[i] = atoms->atom[i].m;
- }
- }
- }
- }
- ~t_corr()
- {
- for (int i = 0; i < nrestart; i++)
- {
- for (int j = 0; j < nmol; j++)
- {
- gmx_stats_free(lsq[i][j]);
- }
- }
- sfree(lsq);
- }
-};
-
-typedef real
-t_calc_func(t_corr* curr, int nx, const int index[], int nx0, rvec xc[], const rvec dcom, gmx_bool bTen, matrix mat);
-
-static real thistime(t_corr* curr)
-{
- return curr->time[curr->nframes];
-}
-
-static int in_data(t_corr* curr, int nx00)
-{
- return curr->nframes - curr->n_offs[nx00];
-}
-
-static void corr_print(t_corr* curr,
- gmx_bool bTen,
- const char* fn,
- const char* title,
- const char* yaxis,
- real msdtime,
- real beginfit,
- real endfit,
- real* DD,
- real* SigmaD,
- char* grpname[],
- const gmx_output_env_t* oenv)
-{
- FILE* out;
- int i, j;
-
- out = xvgropen(fn, title, output_env_get_xvgr_tlabel(oenv), yaxis, oenv);
- if (DD)
- {
- fprintf(out,
- "# MSD gathered over %g %s with %d restarts\n",
- msdtime,
- output_env_get_time_unit(oenv).c_str(),
- curr->nrestart);
- fprintf(out,
- "# Diffusion constants fitted from time %g to %g %s\n",
- beginfit,
- endfit,
- output_env_get_time_unit(oenv).c_str());
- for (i = 0; i < curr->ngrp; i++)
- {
- fprintf(out, "# D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)\n", grpname[i], DD[i], SigmaD[i]);
- }
- }
- for (i = 0; i < curr->nframes; i++)
- {
- fprintf(out, "%10g", output_env_conv_time(oenv, curr->time[i]));
- for (j = 0; j < curr->ngrp; j++)
- {
- fprintf(out, " %10g", curr->data[j][i]);
- if (bTen)
- {
- fprintf(out,
- " %10g %10g %10g %10g %10g %10g",
- curr->datam[j][i][XX][XX],
- curr->datam[j][i][YY][YY],
- curr->datam[j][i][ZZ][ZZ],
- curr->datam[j][i][YY][XX],
- curr->datam[j][i][ZZ][XX],
- curr->datam[j][i][ZZ][YY]);
- }
- }
- fprintf(out, "\n");
- }
- xvgrclose(out);
-}
-
-/* called from corr_loop, to do the main calculations */
-static void
-calc_corr(t_corr* curr, int nr, int nx, int index[], rvec xc[], gmx_bool bRmCOMM, rvec com, t_calc_func* calc1, gmx_bool bTen)
-{
- int nx0;
- real g;
- matrix mat;
- rvec dcom;
-
- /* Check for new starting point */
- if (curr->nlast < curr->nrestart)
- {
- if ((thistime(curr) >= (curr->nlast * curr->delta_t)) && (nr == 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++;
- }
- }
-
- /* nx0 appears to be the number of new starting points,
- * so for all starting points, call calc1.
- */
- for (nx0 = 0; (nx0 < curr->nlast); nx0++)
- {
- if (bRmCOMM)
- {
- rvec_sub(com, curr->com[nx0], dcom);
- }
- else
- {
- clear_rvec(dcom);
- }
- g = calc1(curr, nx, index, nx0, xc, dcom, bTen, mat);
-#ifdef DEBUG2
- printf("g[%d]=%g\n", nx0, g);
-#endif
- curr->data[nr][in_data(curr, nx0)] += g;
- if (bTen)
- {
- m_add(curr->datam[nr][in_data(curr, nx0)], mat, curr->datam[nr][in_data(curr, nx0)]);
- }
- curr->ndata[nr][in_data(curr, nx0)]++;
- }
-}
-
-/* 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)
-{
- int i, ix, m, m2;
- real g, r, r2;
- rvec rv;
-
- g = 0.0;
- clear_mat(mat);
-
- for (i = 0; (i < nx); i++)
- {
- ix = index[i];
- r2 = 0.0;
- switch (curr->type)
- {
- case NORMAL:
- for (m = 0; (m < DIM); m++)
- {
- rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
- r2 += rv[m] * rv[m];
- if (bTen)
- {
- for (m2 = 0; m2 <= m; m2++)
- {
- mat[m][m2] += rv[m] * rv[m2];
- }
- }
- }
- break;
- case X:
- case Y:
- case Z:
- r = xc[ix][curr->type - X] - curr->x0[nx0][ix][curr->type - X] - dcom[curr->type - X];
- r2 += r * r;
- break;
- case LATERAL:
- for (m = 0; (m < DIM); m++)
- {
- if (m != curr->axis)
- {
- r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
- r2 += r * r;
- }
- }
- break;
- default: gmx_fatal(FARGS, "Error: did not expect option value %d", curr->type);
- }
- g += r2;
- }
- g /= nx;
- msmul(mat, 1.0 / nx, mat);
-
- return g;
-}
-
-/* calculate the com of molecules in x and put it into xa */
-static void
-calc_mol_com(int nmol, const int* molindex, const t_block* mols, const t_atoms* atoms, rvec* x, rvec* xa)
-{
- int m, mol, i, d;
- rvec xm;
- real mass, mtot;
-
- for (m = 0; m < nmol; m++)
- {
- mol = molindex[m];
- clear_rvec(xm);
- mtot = 0;
- for (i = mols->index[mol]; i < mols->index[mol + 1]; i++)
- {
- mass = atoms->atom[i].m;
- for (d = 0; d < DIM; d++)
- {
- xm[d] += mass * x[i][d];
- }
- mtot += mass;
- }
- svmul(1 / mtot, xm, xa[m]);
- }
-}
-
-static real calc_one_mw(t_corr* curr, int ix, int nx0, rvec xc[], real* tm, const rvec dcom, gmx_bool bTen, matrix mat)
-{
- real r2, r, mm;
- rvec rv;
- int m, m2;
-
- mm = curr->mass[ix];
- if (mm == 0)
- {
- return 0;
- }
- (*tm) += mm;
- r2 = 0.0;
- switch (curr->type)
- {
- case NORMAL:
- for (m = 0; (m < DIM); m++)
- {
- rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
- r2 += mm * rv[m] * rv[m];
- if (bTen)
- {
- for (m2 = 0; m2 <= m; m2++)
- {
- mat[m][m2] += mm * rv[m] * rv[m2];
- }
- }
- }
- break;
- case X:
- case Y:
- case Z:
- r = xc[ix][curr->type - X] - curr->x0[nx0][ix][curr->type - X] - dcom[curr->type - X];
- r2 = mm * r * r;
- break;
- case LATERAL:
- for (m = 0; (m < DIM); m++)
- {
- if (m != curr->axis)
- {
- r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
- r2 += mm * r * r;
- }
- }
- break;
- default: gmx_fatal(FARGS, "Options got screwed. Did not expect value %d\n", curr->type);
- } /* end switch */
- return r2;
-}
-
-/* the normal, mass-weighted mean-squared displacement calcuation */
-static real calc1_mw(t_corr* curr, int nx, const int index[], int nx0, rvec xc[], const rvec dcom, gmx_bool bTen, matrix mat)
-{
- int i;
- real g, tm;
-
- g = tm = 0.0;
- clear_mat(mat);
- for (i = 0; (i < nx); i++)
- {
- g += calc_one_mw(curr, index[i], nx0, xc, &tm, dcom, bTen, mat);
- }
-
- g /= tm;
- if (bTen)
- {
- msmul(mat, 1 / tm, mat);
- }
-
- return g;
-}
-
-/* prepare the coordinates by removing periodic boundary crossings.
- gnx = the number of atoms/molecules
- index = the indices
- xcur = the current coordinates
- xprev = the previous coordinates
- box = the box matrix */
-static void prep_data(gmx_bool bMol, int gnx, const int index[], rvec xcur[], rvec xprev[], matrix box)
-{
- int i, m, ind;
- rvec hbox;
-
- /* Remove periodicity */
- for (m = 0; (m < DIM); m++)
- {
- hbox[m] = 0.5 * box[m][m];
- }
-
- for (i = 0; (i < gnx); i++)
- {
- if (bMol)
- {
- ind = i;
- }
- else
- {
- ind = index[i];
- }
-
- for (m = DIM - 1; m >= 0; m--)
- {
- if (hbox[m] == 0)
- {
- continue;
- }
- while (xcur[ind][m] - xprev[ind][m] <= -hbox[m])
- {
- rvec_inc(xcur[ind], box[m]);
- }
- while (xcur[ind][m] - xprev[ind][m] > hbox[m])
- {
- rvec_dec(xcur[ind], box[m]);
- }
- }
- }
-}
-
-/* calculate the center of mass for a group
- gnx = the number of atoms/molecules
- index = the indices
- xcur = the current coordinates
- xprev = the previous coordinates
- box = the box matrix
- atoms = atom data (for mass)
- com(output) = center of mass */
-static void
-calc_com(gmx_bool bMol, int gnx, int index[], rvec xcur[], rvec xprev[], matrix box, const t_atoms* atoms, rvec com)
-{
- int i, m, ind;
- real mass;
- double tmass;
- dvec sx;
-
- clear_dvec(sx);
- tmass = 0;
-
- prep_data(bMol, gnx, index, xcur, xprev, box);
- for (i = 0; (i < gnx); i++)
- {
- if (bMol)
- {
- ind = i;
- }
- else
- {
- ind = index[i];
- }
-
-
- mass = atoms->atom[ind].m;
- for (m = 0; m < DIM; m++)
- {
- sx[m] += mass * xcur[ind][m];
- }
- tmass += mass;
- }
- for (m = 0; m < DIM; m++)
- {
- com[m] = sx[m] / tmass;
- }
-}
-
-
-static real calc1_mol(t_corr* curr,
- int nx,
- const int gmx_unused index[],
- int nx0,
- rvec xc[],
- const rvec dcom,
- gmx_bool bTen,
- matrix mat)
-{
- int i;
- real g, tm, gtot, tt;
-
- tt = curr->time[in_data(curr, nx0)];
- gtot = 0;
- tm = 0;
- clear_mat(mat);
- for (i = 0; (i < nx); i++)
- {
- g = calc_one_mw(curr, i, nx0, xc, &tm, dcom, bTen, mat);
- /* We don't need to normalize as the mass was set to 1 */
- gtot += g;
- if (tt >= curr->beginfit && (curr->endfit < 0 || tt <= curr->endfit))
- {
- gmx_stats_add_point(curr->lsq[nx0][i], tt, g, 0, 0);
- }
- }
- msmul(mat, 1.0 / nx, mat);
-
- return gtot / nx;
-}
-
-static void printmol(t_corr* curr,
- const char* fn,
- const char* fn_pdb,
- const int* molindex,
- const t_topology* top,
- rvec* x,
- PbcType pbcType,
- matrix box,
- const gmx_output_env_t* oenv)
-{
- FILE* out;
- gmx_stats_t lsq1;
- int i, j;
- real a, b, D, Dav, D2av, VarD, sqrtD, sqrtD_max, scale;
- t_pdbinfo* pdbinfo = nullptr;
- const int* mol2a = nullptr;
-
- out = xvgropen(fn, "Diffusion Coefficients / Molecule", "Molecule", "D (1e-5 cm^2/s)", oenv);
-
- if (fn_pdb)
- {
- pdbinfo = top->atoms.pdbinfo;
- mol2a = top->mols.index;
- }
-
- Dav = D2av = 0;
- sqrtD_max = 0;
- for (i = 0; (i < curr->nmol); i++)
- {
- lsq1 = gmx_stats_init();
- for (j = 0; (j < curr->nrestart); j++)
- {
- real xx, yy, dx, dy;
-
- while (gmx_stats_get_point(curr->lsq[j][i], &xx, &yy, &dx, &dy, 0) == estatsOK)
- {
- gmx_stats_add_point(lsq1, xx, yy, dx, dy);
- }
- }
- gmx_stats_get_ab(lsq1, elsqWEIGHT_NONE, &a, &b, nullptr, nullptr, nullptr, nullptr);
- gmx_stats_free(lsq1);
- D = a * diffusionConversionFactor / curr->dim_factor;
- if (D < 0)
- {
- D = 0;
- }
- Dav += D;
- D2av += gmx::square(D);
- fprintf(out, "%10d %10g\n", i, D);
- if (pdbinfo)
- {
- sqrtD = std::sqrt(D);
- if (sqrtD > sqrtD_max)
- {
- sqrtD_max = sqrtD;
- }
- for (j = mol2a[molindex[i]]; j < mol2a[molindex[i] + 1]; j++)
- {
- pdbinfo[j].bfac = sqrtD;
- }
- }
- }
- xvgrclose(out);
- fprintf(stdout, "Wrote per-molecule output to %s\n", fn);
- do_view(oenv, fn, "-graphtype bar");
-
- /* Compute variance, stddev and error */
- Dav /= curr->nmol;
- D2av /= curr->nmol;
- VarD = D2av - gmx::square(Dav);
- printf("<D> = %.4f Std. Dev. = %.4f Error = %.4f\n", Dav, std::sqrt(VarD), std::sqrt(VarD / curr->nmol));
-
- if (fn_pdb && x)
- {
- scale = 1;
- while (scale * sqrtD_max > 10)
- {
- scale *= 0.1;
- }
- while (scale * sqrtD_max < 0.1)
- {
- scale *= 10;
- }
- GMX_RELEASE_ASSERT(pdbinfo != nullptr, "Internal error - pdbinfo not set for PDB input");
- for (i = 0; i < top->atoms.nr; i++)
- {
- pdbinfo[i].bfac *= scale;
- }
- write_sto_conf(fn_pdb, "molecular MSD", &top->atoms, x, nullptr, pbcType, box);
- fprintf(stdout, "Wrote frame for -tpdb to %s\n", fn_pdb);
- }
-}
-
-/* this is the main loop for the correlation type functions
- * fx and nx are file pointers to things like read_first_x and
- * read_next_x
- */
-static int corr_loop(t_corr* curr,
- const char* fn,
- const t_topology* top,
- PbcType pbcType,
- gmx_bool bMol,
- int gnx[],
- int* index[],
- 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)
-{
- rvec* x[2]; /* the coordinates to read */
- rvec* xa[2]; /* the coordinates to calculate displacements for */
- rvec com = { 0 };
- real t, t_prev = 0;
- int natoms, i, j, cur = 0, maxframes = 0;
- t_trxstatus* status;
-#define prev (1 - cur)
- matrix box;
- gmx_bool bFirst;
- gmx_rmpbc_t gpbc = nullptr;
-
- natoms = read_first_x(oenv, &status, fn, &curr->t0, &(x[cur]), box);
-#ifdef DEBUG
- fprintf(stderr, "Read %d atoms for first frame\n", natoms);
-#endif
- 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);
- }
-
- snew(x[prev], natoms);
-
- // if com is requested, the data structure needs to be large enough to do this
- // to prevent overflow
- if (bMol && gnx_com.empty())
- {
- curr->ncoords = curr->nmol;
- snew(xa[0], curr->ncoords);
- snew(xa[1], curr->ncoords);
- }
- else
- {
- curr->ncoords = natoms;
- xa[0] = x[0];
- xa[1] = x[1];
- }
-
- bFirst = TRUE;
- t = curr->t0;
- if (x_pdb)
- {
- *x_pdb = nullptr;
- }
-
- if (bMol)
- {
- gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
- }
-
- /* the loop over all frames */
- do
- {
- if (x_pdb
- && ((bFirst && t_pdb < t)
- || (!bFirst && t_pdb > t - 0.5 * (t - t_prev) && t_pdb < t + 0.5 * (t - t_prev))))
- {
- if (*x_pdb == nullptr)
- {
- snew(*x_pdb, natoms);
- }
- for (i = 0; i < natoms; i++)
- {
- copy_rvec(x[cur][i], (*x_pdb)[i]);
- }
- copy_mat(box, box_pdb);
- }
-
-
- /* check whether we've reached a restart point */
- if (bRmod(t, curr->t0, dt))
- {
- 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++)
- {
- curr->lsq[curr->nrestart - 1][i] = gmx_stats_init();
- }
-
- if (debug)
- {
- fprintf(debug, "Extended data structures because of new restart %d\n", curr->nrestart);
- }
- }
- /* create or extend the frame-based arrays */
- if (curr->nframes >= maxframes - 1)
- {
- if (maxframes == 0)
- {
- for (i = 0; (i < curr->ngrp); i++)
- {
- if (bTen)
- {
- curr->datam[i] = nullptr;
- }
- }
- }
- maxframes += 10;
- for (i = 0; (i < curr->ngrp); i++)
- {
- curr->ndata[i].resize(maxframes);
- curr->data[i].resize(maxframes);
- if (bTen)
- {
- srenew(curr->datam[i], maxframes);
- }
- for (j = maxframes - 10; j < maxframes; j++)
- {
- curr->ndata[i][j] = 0;
- curr->data[i][j] = 0;
- if (bTen)
- {
- clear_mat(curr->datam[i][j]);
- }
- }
- }
- curr->time.resize(maxframes);
- }
-
- /* set the time */
- curr->time[curr->nframes] = t - curr->t0;
-
- /* make the molecules whole */
- if (bMol)
- {
- gmx_rmpbc(gpbc, natoms, box, x[cur]);
- }
-
- /* calculate the molecules' centers of masses and put them into xa */
- // NOTE and WARNING! If above both COM removal and individual molecules have been
- // requested, x and xa point to the same memory, and the coordinate
- // data becomes overwritten by the molecule data.
- if (bMol)
- {
- calc_mol_com(gnx[0], index[0], &top->mols, &top->atoms, x[cur], xa[cur]);
- }
-
- /* for the first frame, the previous frame is a copy of the first frame */
- if (bFirst)
- {
- std::memcpy(xa[prev], xa[cur], curr->ncoords * sizeof(xa[prev][0]));
- bFirst = FALSE;
- }
-
- /* first remove the periodic boundary condition crossings */
- for (i = 0; i < curr->ngrp; i++)
- {
- prep_data(bMol, gnx[i], index[i], xa[cur], xa[prev], box);
- }
-
- /* calculate the center of mass */
- if (!gnx_com.empty())
- {
- GMX_RELEASE_ASSERT(index_com != nullptr,
- "Center-of-mass removal must have valid index group");
- calc_com(bMol, gnx_com[0], index_com[0], xa[cur], xa[prev], box, &top->atoms, com);
- }
-
- /* loop over all groups in index file */
- 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.empty()), com, calc1, bTen);
- }
- cur = prev;
- t_prev = t;
-
- curr->nframes++;
- } while (read_next_x(oenv, status, &t, x[cur], box));
- fprintf(stderr,
- "\nUsed %d restart points spaced %g %s over %g %s\n\n",
- curr->nrestart,
- output_env_conv_time(oenv, dt),
- output_env_get_time_unit(oenv).c_str(),
- output_env_conv_time(oenv, curr->time[curr->nframes - 1]),
- output_env_get_time_unit(oenv).c_str());
-
- if (bMol)
- {
- gmx_rmpbc_done(gpbc);
- }
-
- close_trx(status);
-
- return natoms;
-}
-
-static void index_atom2mol(int* n, int* index, const t_block* mols)
-{
- int nat, i, nmol, mol, j;
-
- nat = *n;
- i = 0;
- nmol = 0;
- mol = 0;
- while (i < nat)
- {
- while (index[i] > mols->index[mol])
- {
- mol++;
- if (mol >= mols->nr)
- {
- gmx_fatal(FARGS, "Atom index out of range: %d", index[i] + 1);
- }
- }
- for (j = mols->index[mol]; j < mols->index[mol + 1]; j++)
- {
- if (i >= nat || index[i] != j)
- {
- gmx_fatal(FARGS, "The index group does not consist of whole molecules");
- }
- i++;
- }
- index[nmol++] = mol;
- }
-
- fprintf(stderr, "Split group of %d atoms into %d molecules\n", nat, nmol);
-
- *n = nmol;
-}
-
-static void do_corr(const char* trx_file,
- const char* ndx_file,
- const char* msd_file,
- const char* mol_file,
- const char* pdb_file,
- real t_pdb,
- int nrgrp,
- t_topology* top,
- PbcType pbcType,
- gmx_bool bTen,
- gmx_bool bMW,
- gmx_bool bRmCOMM,
- int type,
- real dim_factor,
- int axis,
- real dt,
- real beginfit,
- real endfit,
- const gmx_output_env_t* oenv)
-{
- std::unique_ptr<t_corr> msd;
- 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;
- std::vector<real> SigmaD, DD;
- real a, a2, b, r, chi2;
- rvec* x = nullptr;
- matrix box;
- int** index_com = nullptr; /* the COM removal group atom indices */
- char** grpname_com = nullptr; /* the COM removal group name */
-
- 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.data(), index, grpname);
-
- if (bRmCOMM)
- {
- 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.data(), index_com, grpname_com);
- }
-
- if (mol_file)
- {
- index_atom2mol(&gnx[0], index[0], &top->mols);
- }
-
- msd = std::make_unique<t_corr>(
- nrgrp, type, axis, dim_factor, mol_file == nullptr ? 0 : gnx[0], bTen, bMW, dt, top, beginfit, endfit);
-
- nat_trx = corr_loop(msd.get(),
- trx_file,
- top,
- pbcType,
- 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);
-
- /* Correct for the number of points */
- for (j = 0; (j < msd->ngrp); j++)
- {
- for (i = 0; (i < msd->nframes); i++)
- {
- msd->data[j][i] /= msd->ndata[j][i];
- if (bTen)
- {
- msmul(msd->datam[j][i], 1.0 / msd->ndata[j][i], msd->datam[j][i]);
- }
- }
- }
-
- if (mol_file)
- {
- if (pdb_file && x == nullptr)
- {
- fprintf(stderr,
- "\nNo frame found need time tpdb = %g ps\n"
- "Can not write %s\n\n",
- t_pdb,
- pdb_file);
- }
- i = top->atoms.nr;
- top->atoms.nr = nat_trx;
- if (pdb_file && top->atoms.pdbinfo == nullptr)
- {
- snew(top->atoms.pdbinfo, top->atoms.nr);
- }
- printmol(msd.get(), mol_file, pdb_file, index[0], top, x, pbcType, box, oenv);
- top->atoms.nr = i;
- }
-
- if (beginfit == -1)
- {
- i0 = gmx::roundToInt(0.1 * (msd->nframes - 1));
- beginfit = msd->time[i0];
- }
- else
- {
- for (i0 = 0; i0 < msd->nframes && msd->time[i0] < beginfit; i0++) {}
- }
-
- if (endfit == -1)
- {
- i1 = gmx::roundToInt(0.9 * (msd->nframes - 1)) + 1;
- endfit = msd->time[i1 - 1];
- }
- else
- {
- for (i1 = i0; i1 < msd->nframes && msd->time[i1] <= endfit; i1++) {}
- }
- fprintf(stdout, "Fitting from %g to %g %s\n\n", beginfit, endfit, output_env_get_time_unit(oenv).c_str());
-
- N = i1 - i0;
- if (N <= 2)
- {
- fprintf(stdout,
- "Not enough points for fitting (%d).\n"
- "Can not determine the diffusion constant.\n",
- N);
- }
- else
- {
- DD.resize(msd->ngrp);
- SigmaD.resize(msd->ngrp);
- for (j = 0; j < msd->ngrp; j++)
- {
- if (N >= 4)
- {
- lsq_y_ax_b(N / 2, &(msd->time[i0]), &(msd->data[j][i0]), &a, &b, &r, &chi2);
- lsq_y_ax_b(N / 2, &(msd->time[i0 + N / 2]), &(msd->data[j][i0 + N / 2]), &a2, &b, &r, &chi2);
- SigmaD[j] = std::abs(a - a2);
- }
- else
- {
- SigmaD[j] = 0;
- }
- lsq_y_ax_b(N, &(msd->time[i0]), &(msd->data[j][i0]), &(DD[j]), &b, &r, &chi2);
- 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", grpname[j], DD[j], SigmaD[j]);
- }
- else
- {
- fprintf(stdout, "D[%10s] %.4g (+/- %.4g) 1e-5 cm^2/s\n", grpname[j], DD[j], SigmaD[j]);
- }
- }
- }
- /* Print mean square displacement */
- corr_print(msd.get(),
- bTen,
- msd_file,
- "Mean Square Displacement",
- "MSD (nm\\S2\\N)",
- msd->time[msd->nframes - 1],
- beginfit,
- endfit,
- DD.data(),
- SigmaD.data(),
- grpname,
- oenv);
-}
-
-int gmx_msd(int argc, char* argv[])
-{
- const char* desc[] = {
- "[THISMODULE] computes the mean square displacement (MSD) of atoms from",
- "a set of initial positions. This provides an easy way to compute",
- "the diffusion constant using the Einstein relation.",
- "The time between the reference points for the MSD calculation",
- "is set with [TT]-trestart[tt].",
- "The diffusion constant is calculated by least squares fitting a",
- "straight line (D*t + c) through the MSD(t) from [TT]-beginfit[tt] to",
- "[TT]-endfit[tt] (note that t is time from the reference positions,",
- "not simulation time). An error estimate given, which is the difference",
- "of the diffusion coefficients obtained from fits over the two halves",
- "of the fit interval.[PAR]",
- "There are three, mutually exclusive, options to determine different",
- "types of mean square displacement: [TT]-type[tt], [TT]-lateral[tt]",
- "and [TT]-ten[tt]. Option [TT]-ten[tt] writes the full MSD tensor for",
- "each group, the order in the output is: trace xx yy zz yx zx zy.[PAR]",
- "If [TT]-mol[tt] is set, [THISMODULE] plots the MSD for individual molecules",
- "(including making molecules whole across periodic boundaries): ",
- "for each individual molecule a diffusion constant is computed for ",
- "its center of mass. The chosen index group will be split into ",
- "molecules.[PAR]",
- "The default way to calculate a MSD is by using mass-weighted averages.",
- "This can be turned off with [TT]-nomw[tt].[PAR]",
- "With the option [TT]-rmcomm[tt], the center of mass motion of a ",
- "specific group can be removed. For trajectories produced with ",
- "GROMACS this is usually not necessary, ",
- "as [gmx-mdrun] usually already removes the center of mass motion.",
- "When you use this option be sure that the whole system is stored",
- "in the trajectory file.[PAR]",
- "The diffusion coefficient is determined by linear regression of the MSD.",
- "When [TT]-beginfit[tt] is -1, fitting starts at 10%",
- "and when [TT]-endfit[tt] is -1, fitting goes to 90%.",
- "Using this option one also gets an accurate error estimate",
- "based on the statistics between individual molecules.",
- "Note that this diffusion coefficient and error estimate are only",
- "accurate when the MSD is completely linear between",
- "[TT]-beginfit[tt] and [TT]-endfit[tt].[PAR]",
- "Option [TT]-pdb[tt] writes a [REF].pdb[ref] file with the coordinates of the frame",
- "at time [TT]-tpdb[tt] with in the B-factor field the square root of",
- "the diffusion coefficient of the molecule.",
- "This option implies option [TT]-mol[tt]."
- };
- static const char* normtype[] = { nullptr, "no", "x", "y", "z", nullptr };
- static const char* axtitle[] = { nullptr, "no", "x", "y", "z", nullptr };
- static int ngroup = 1;
- static real dt = 10;
- static real t_pdb = 0;
- static real beginfit = -1;
- static real endfit = -1;
- static gmx_bool bTen = FALSE;
- static gmx_bool bMW = TRUE;
- static gmx_bool bRmCOMM = FALSE;
- t_pargs pa[] = {
- { "-type", FALSE, etENUM, { normtype }, "Compute diffusion coefficient in one direction" },
- { "-lateral",
- FALSE,
- etENUM,
- { axtitle },
- "Calculate the lateral diffusion in a plane perpendicular to" },
- { "-ten", FALSE, etBOOL, { &bTen }, "Calculate the full tensor" },
- { "-ngroup", FALSE, etINT, { &ngroup }, "Number of groups to calculate MSD for" },
- { "-mw", FALSE, etBOOL, { &bMW }, "Mass weighted MSD" },
- { "-rmcomm", FALSE, etBOOL, { &bRmCOMM }, "Remove center of mass motion" },
- { "-tpdb", FALSE, etTIME, { &t_pdb }, "The frame to use for option [TT]-pdb[tt] (%t)" },
- { "-trestart", FALSE, etTIME, { &dt }, "Time between restarting points in trajectory (%t)" },
- { "-beginfit",
- FALSE,
- etTIME,
- { &beginfit },
- "Start time for fitting the MSD (%t), -1 is 10%" },
- { "-endfit", FALSE, etTIME, { &endfit }, "End time for fitting the MSD (%t), -1 is 90%" }
- };
-
- t_filenm fnm[] = {
- { efTRX, nullptr, nullptr, ffREAD }, { efTPS, nullptr, nullptr, ffREAD },
- { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, nullptr, "msd", ffWRITE },
- { efXVG, "-mol", "diff_mol", ffOPTWR }, { efPDB, "-pdb", "diff_mol", ffOPTWR }
- };
-#define NFILE asize(fnm)
-
- t_topology top;
- PbcType pbcType;
- matrix box;
- const char * trx_file, *tps_file, *ndx_file, *msd_file, *mol_file, *pdb_file;
- rvec* xdum;
- gmx_bool bTop;
- int axis, type;
- real dim_factor;
- gmx_output_env_t* oenv;
-
- if (!parse_common_args(&argc,
- argv,
- PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_TIME_UNIT,
- NFILE,
- fnm,
- asize(pa),
- pa,
- asize(desc),
- desc,
- 0,
- nullptr,
- &oenv))
- {
- return 0;
- }
- trx_file = ftp2fn_null(efTRX, NFILE, fnm);
- tps_file = ftp2fn_null(efTPS, NFILE, fnm);
- ndx_file = ftp2fn_null(efNDX, NFILE, fnm);
- msd_file = ftp2fn_null(efXVG, NFILE, fnm);
- pdb_file = opt2fn_null("-pdb", NFILE, fnm);
- if (pdb_file)
- {
- mol_file = opt2fn("-mol", NFILE, fnm);
- }
- else
- {
- mol_file = opt2fn_null("-mol", NFILE, fnm);
- }
-
- if (ngroup < 1)
- {
- gmx_fatal(FARGS, "Must have at least 1 group (now %d)", ngroup);
- }
- if (mol_file && ngroup > 1)
- {
- gmx_fatal(FARGS, "With molecular msd can only have 1 group (now %d)", ngroup);
- }
-
-
- if (mol_file)
- {
- bMW = TRUE;
- fprintf(stderr, "Calculating diffusion coefficients for molecules.\n");
- }
-
- GMX_RELEASE_ASSERT(normtype[0] != nullptr, "Options inconsistency; normtype[0]Â is NULL");
- GMX_RELEASE_ASSERT(axtitle[0] != nullptr, "Options inconsistency; axtitle[0]Â is NULL");
-
- if (normtype[0][0] != 'n')
- {
- type = normtype[0][0] - 'x' + X; /* See defines above */
- dim_factor = 2.0;
- }
- else
- {
- type = NORMAL;
- dim_factor = 6.0;
- }
- if ((type == NORMAL) && (axtitle[0][0] != 'n'))
- {
- type = LATERAL;
- dim_factor = 4.0;
- axis = (axtitle[0][0] - 'x'); /* See defines above */
- }
- else
- {
- axis = 0;
- }
-
- if (bTen && type != NORMAL)
- {
- gmx_fatal(FARGS, "Can only calculate the full tensor for 3D msd");
- }
-
- bTop = read_tps_conf(tps_file, &top, &pbcType, &xdum, nullptr, box, bMW || bRmCOMM);
- if (mol_file && !bTop)
- {
- gmx_fatal(FARGS, "Could not read a topology from %s. Try a tpr file instead.", tps_file);
- }
-
- do_corr(trx_file,
- ndx_file,
- msd_file,
- mol_file,
- pdb_file,
- t_pdb,
- ngroup,
- &top,
- pbcType,
- bTen,
- bMW,
- bRmCOMM,
- type,
- dim_factor,
- axis,
- dt,
- beginfit,
- endfit,
- oenv);
-
- done_top(&top);
- view_all(oenv, NFILE, fnm);
-
- return 0;
-}
# This file is part of the GROMACS molecular simulation package.
#
# Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
-# Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+# Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
entropy.cpp
gmx_traj.cpp
gmx_mindist.cpp
- gmx_msd.cpp
)
gmx_register_gtest_test(GmxAnaTest ${exename} INTEGRATION_TEST IGNORE_LEAKS)
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-/*! \internal \file
- * \brief
- * Tests for gmx msd.
- *
- * \author Kevin Boyd <kevin.boyd@uconn.edu>
- */
-
-#include "gmxpre.h"
-
-#include <cstdio>
-#include <cstdlib>
-
-#include "gromacs/gmxana/gmx_ana.h"
-#include "gromacs/gmxpreprocess/grompp.h"
-#include "gromacs/utility/futil.h"
-#include "gromacs/utility/path.h"
-#include "gromacs/utility/textreader.h"
-
-#include "testutils/cmdlinetest.h"
-#include "testutils/refdata.h"
-#include "testutils/testfilemanager.h"
-#include "testutils/textblockmatchers.h"
-#include "testutils/xvgtest.h"
-
-namespace
-{
-
-using gmx::test::CommandLine;
-using gmx::test::XvgMatch;
-
-class MsdTest : public gmx::test::CommandLineTestBase
-{
-public:
- MsdTest()
- {
- setOutputFile("-o", "msd.xvg", XvgMatch());
- setInputFile("-f", "msd_traj.xtc");
- setInputFile("-s", "msd_coords.gro");
- setInputFile("-n", "msd.ndx");
- }
-
- void runTest(const CommandLine& args)
- {
- CommandLine& cmdline = commandLine();
- cmdline.merge(args);
- ASSERT_EQ(0, gmx_msd(cmdline.argc(), cmdline.argv()));
- checkOutputFiles();
- }
-};
-
-class MsdMolTest : public gmx::test::CommandLineTestBase
-{
-public:
- MsdMolTest()
- {
- double tolerance = 1e-5;
- XvgMatch xvg;
- XvgMatch& toler = xvg.tolerance(gmx::test::relativeToleranceAsFloatingPoint(1, tolerance));
- setOutputFile("-mol", "msdmol.xvg", toler);
- }
-
- void runTest(const CommandLine& args, const char* ndxfile, const std::string& simulationName)
- {
- setInputFile("-f", simulationName + ".pdb");
- std::string tpr = fileManager().getTemporaryFilePath(".tpr");
- std::string mdp = fileManager().getTemporaryFilePath(".mdp");
- FILE* fp = fopen(mdp.c_str(), "w");
- fprintf(fp, "cutoff-scheme = verlet\n");
- fprintf(fp, "rcoulomb = 0.85\n");
- fprintf(fp, "rvdw = 0.85\n");
- fprintf(fp, "rlist = 0.85\n");
- fclose(fp);
-
- // Prepare a .tpr file
- {
- CommandLine caller;
- auto simDB = gmx::test::TestFileManager::getTestSimulationDatabaseDirectory();
- auto base = gmx::Path::join(simDB, simulationName);
- caller.append("grompp");
- caller.addOption("-maxwarn", 0);
- caller.addOption("-f", mdp.c_str());
- std::string gro = (base + ".pdb");
- caller.addOption("-c", gro.c_str());
- std::string top = (base + ".top");
- caller.addOption("-p", top.c_str());
- std::string ndx = (base + ".ndx");
- caller.addOption("-n", ndx.c_str());
- caller.addOption("-o", tpr.c_str());
- ASSERT_EQ(0, gmx_grompp(caller.argc(), caller.argv()));
- }
- // Run the MSD analysis
- {
- setInputFile("-n", ndxfile);
- CommandLine& cmdline = commandLine();
- cmdline.merge(args);
- cmdline.addOption("-s", tpr.c_str());
- ASSERT_EQ(0, gmx_msd(cmdline.argc(), cmdline.argv()));
- checkOutputFiles();
- }
- }
-};
-
-/* msd_traj.xtc contains a 10 frame (1 ps per frame) simulation
- * containing 3 atoms, with different starting positions but identical
- * displacements. The displacements are calculated to yield the following
- * diffusion coefficients when lag is calculated ONLY FROM TIME 0
- * D_x = 8 * 10 ^ -5 cm^2 /s, D_y = 4 * 10^ -5 cm^2 /s , D_z = 0
- *
- * To test for these results, -trestart is set to a larger value than the
- * total simulation length, so that only lag 0 is calculated
- */
-
-// for 3D, (8 + 4 + 0) / 3 should yield 4 cm^2 / s
-TEST_F(MsdTest, threeDimensionalDiffusion)
-{
- const char* const cmdline[] = {
- "msd", "-mw", "no", "-trestart", "200",
- };
- runTest(CommandLine(cmdline));
-}
-
-// for lateral z, (8 + 4) / 2 should yield 6 cm^2 /s
-TEST_F(MsdTest, twoDimensionalDiffusion)
-{
- const char* const cmdline[] = { "msd", "-mw", "no", "-trestart", "200", "-lateral", "z" };
- runTest(CommandLine(cmdline));
-}
-
-// for type x, should yield 8 cm^2 / s
-TEST_F(MsdTest, oneDimensionalDiffusion)
-{
- const char* const cmdline[] = { "msd", "-mw", "no", "-trestart", "200", "-type", "x" };
- runTest(CommandLine(cmdline));
-}
-
-// Test the diffusion per molecule output, mass weighted
-TEST_F(MsdMolTest, diffMolMassWeighted)
-{
- const char* const cmdline[] = { "msd", "-trestart", "200" };
- runTest(CommandLine(cmdline), "spc5.ndx", "spc5");
-}
-
-// Test the diffusion per molecule output, non-mass weighted
-TEST_F(MsdMolTest, diffMolNonMassWeighted)
-{
- const char* const cmdline[] = { "msd", "-trestart", "200", "-mw", "no" };
- runTest(CommandLine(cmdline), "spc5.ndx", "spc5");
-}
-
-// Test the diffusion per molecule output, with selection
-TEST_F(MsdMolTest, diffMolSelected)
-{
- const char* const cmdline[] = { "msd", "-trestart", "200" };
- runTest(CommandLine(cmdline), "spc5_3.ndx", "spc5");
-}
-
-} // namespace
+++ /dev/null
-<?xml version="1.0"?>
-<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
-<ReferenceData>
- <OutputFiles Name="Files">
- <File Name="-mol">
- <XvgLegend Name="Legend">
- <String Name="XvgLegend"><![CDATA[
-title "Diffusion Coefficients / Molecule"
-xaxis label "Molecule"
-yaxis label "D (1e-5 cm^2/s)"
-TYPE xy
-]]></String>
- </XvgLegend>
- <XvgData Name="Data">
- <Sequence Name="Row0">
- <Int Name="Length">2</Int>
- <Real>0</Real>
- <Real>0.747088</Real>
- </Sequence>
- <Sequence Name="Row1">
- <Int Name="Length">2</Int>
- <Real>1</Real>
- <Real>0.129497</Real>
- </Sequence>
- <Sequence Name="Row2">
- <Int Name="Length">2</Int>
- <Real>2</Real>
- <Real>1.00167</Real>
- </Sequence>
- <Sequence Name="Row3">
- <Int Name="Length">2</Int>
- <Real>3</Real>
- <Real>21.2013</Real>
- </Sequence>
- <Sequence Name="Row4">
- <Int Name="Length">2</Int>
- <Real>4</Real>
- <Real>9.28605</Real>
- </Sequence>
- </XvgData>
- </File>
- </OutputFiles>
-</ReferenceData>
+++ /dev/null
-<?xml version="1.0"?>
-<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
-<ReferenceData>
- <OutputFiles Name="Files">
- <File Name="-mol">
- <XvgLegend Name="Legend">
- <String Name="XvgLegend"><![CDATA[
-title "Diffusion Coefficients / Molecule"
-xaxis label "Molecule"
-yaxis label "D (1e-5 cm^2/s)"
-TYPE xy
-]]></String>
- </XvgLegend>
- <XvgData Name="Data">
- <Sequence Name="Row0">
- <Int Name="Length">2</Int>
- <Real>0</Real>
- <Real>0.747088</Real>
- </Sequence>
- <Sequence Name="Row1">
- <Int Name="Length">2</Int>
- <Real>1</Real>
- <Real>0.129497</Real>
- </Sequence>
- <Sequence Name="Row2">
- <Int Name="Length">2</Int>
- <Real>2</Real>
- <Real>1.00167</Real>
- </Sequence>
- <Sequence Name="Row3">
- <Int Name="Length">2</Int>
- <Real>3</Real>
- <Real>21.2013</Real>
- </Sequence>
- <Sequence Name="Row4">
- <Int Name="Length">2</Int>
- <Real>4</Real>
- <Real>9.28605</Real>
- </Sequence>
- </XvgData>
- </File>
- </OutputFiles>
-</ReferenceData>
+++ /dev/null
-<?xml version="1.0"?>
-<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
-<ReferenceData>
- <OutputFiles Name="Files">
- <File Name="-mol">
- <XvgLegend Name="Legend">
- <String Name="XvgLegend"><![CDATA[
-title "Diffusion Coefficients / Molecule"
-xaxis label "Molecule"
-yaxis label "D (1e-5 cm^2/s)"
-TYPE xy
-]]></String>
- </XvgLegend>
- <XvgData Name="Data">
- <Sequence Name="Row0">
- <Int Name="Length">2</Int>
- <Real>0</Real>
- <Real>0.747088</Real>
- </Sequence>
- <Sequence Name="Row1">
- <Int Name="Length">2</Int>
- <Real>1</Real>
- <Real>0.129497</Real>
- </Sequence>
- <Sequence Name="Row2">
- <Int Name="Length">2</Int>
- <Real>2</Real>
- <Real>21.2013</Real>
- </Sequence>
- </XvgData>
- </File>
- </OutputFiles>
-</ReferenceData>
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2010,2011,2012,2013,2014 by the GROMACS development team.
- * Copyright (c) 2016,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2016,2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "modules/distance.h"
#include "modules/extract_cluster.h"
#include "modules/freevolume.h"
+#include "modules/msd.h"
#include "modules/pairdist.h"
#include "modules/rdf.h"
#include "modules/sasa.h"
registerModule<DistanceInfo>(manager, group);
registerModule<ExtractClusterInfo>(manager, group);
registerModule<FreeVolumeInfo>(manager, group);
+ registerModule<MsdInfo>(manager, group);
registerModule<PairDistanceInfo>(manager, group);
registerModule<RdfInfo>(manager, group);
registerModule<SasaInfo>(manager, group);
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Defines the trajectory analysis module for mean squared displacement calculations.
+ *
+ * \author Kevin Boyd <kevin44boyd@gmail.com>
+ * \ingroup module_trajectoryanalysis
+ */
+
+#include "gmxpre.h"
+
+#include "msd.h"
+
+#include <numeric>
+
+#include "gromacs/analysisdata/analysisdata.h"
+#include "gromacs/analysisdata/modules/average.h"
+#include "gromacs/analysisdata/modules/plot.h"
+#include "gromacs/analysisdata/paralleloptions.h"
+#include "gromacs/fileio/oenv.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/options/basicoptions.h"
+#include "gromacs/options/filenameoption.h"
+#include "gromacs/options/ioptionscontainer.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/selection/selectionoption.h"
+#include "gromacs/statistics/statistics.h"
+#include "gromacs/trajectoryanalysis/analysissettings.h"
+#include "gromacs/trajectoryanalysis/topologyinformation.h"
+#include "gromacs/trajectory/trajectoryframe.h"
+#include "gromacs/utility/stringutil.h"
+#include "gromacs/utility.h"
+
+namespace gmx
+{
+namespace analysismodules
+{
+
+namespace
+{
+
+//! Convert nm^2/ps to 10e-5 cm^2/s
+constexpr double c_diffusionConversionFactor = 1000.0;
+//! Three dimensional diffusion coefficient multiplication constant.
+constexpr double c_3DdiffusionDimensionFactor = 6.0;
+//! Two dimensional diffusion coefficient multiplication constant.
+constexpr double c_2DdiffusionDimensionFactor = 4.0;
+//! One dimensional diffusion coefficient multiplication constant.
+constexpr double c_1DdiffusionDimensionFactor = 2.0;
+
+
+/*! \brief Mean Squared Displacement data accumulator
+ *
+ * This class is used to accumulate individual MSD data points
+ * and emit tau-averaged results once data is finished collecting. Displacements at each observed
+ * time difference (tau) are recorded from the trajectory. Because it is not known in advance which
+ * time differences will be observed from the trajectory, this data structure is built adaptively.
+ * New columns corresponding to observed time differences are added as needed, and additional
+ * observations at formerly observed time differences are added to those columns. Separate time lags
+ * will likely have differing total data points.
+ *
+ * Data columns per tau are accessed via operator[], which always guarantees
+ * a column is initialized and returns an MsdColumProxy to the column that can push data.
+ */
+class MsdData
+{
+public:
+ //! Proxy to a MsdData tau column vector. Supports only push_back.
+ class MsdColumnProxy
+ {
+ public:
+ MsdColumnProxy(std::vector<double>& column) : column_(column) {}
+
+ void push_back(double value) { column_.push_back(value); }
+
+ private:
+ std::vector<double>& column_;
+ };
+ //! Returns a proxy to the column for the given tau index. Guarantees that the column is initialized.
+ MsdColumnProxy operator[](size_t index)
+ {
+ if (msds_.size() <= index)
+ {
+ msds_.resize(index + 1);
+ }
+ return MsdColumnProxy(msds_[index]);
+ }
+ /*! \brief Compute per-tau MSDs averaged over all added points.
+ *
+ * The resulting vector is size(max tau index). Any indices
+ * that have no data points have MSD set to 0.
+ *
+ * \return Average MSD per tau
+ */
+ [[nodiscard]] std::vector<real> averageMsds() const;
+
+private:
+ //! Results - first indexed by tau, then data points
+ std::vector<std::vector<double>> msds_;
+};
+
+
+std::vector<real> MsdData::averageMsds() const
+{
+ std::vector<real> msdSums;
+ msdSums.reserve(msds_.size());
+ for (gmx::ArrayRef<const double> msdValues : msds_)
+ {
+ if (msdValues.empty())
+ {
+ msdSums.push_back(0.0);
+ continue;
+ }
+ msdSums.push_back(std::accumulate(msdValues.begin(), msdValues.end(), 0.0, std::plus<>())
+ / msdValues.size());
+ }
+ return msdSums;
+}
+
+/*! \brief Calculates 1,2, or 3D distance for two vectors.
+ *
+ * \todo Remove NOLINTs once clang-tidy is updated to v11, it should be able to handle constexpr.
+ *
+ * \tparam x If true, calculate x dimension of displacement
+ * \tparam y If true, calculate y dimension of displacement
+ * \tparam z If true, calculate z dimension of displacement
+ * \param c1 First point
+ * \param c2 Second point
+ * \return Euclidian distance for the given dimension.
+ */
+template<bool x, bool y, bool z>
+inline double calcSingleSquaredDistance(RVec c1, const RVec c2)
+{
+ static_assert(x || y || z, "zero-dimensional MSD selected");
+ const DVec firstCoords = c1.toDVec();
+ const DVec secondCoords = c2.toDVec();
+ double result = 0;
+ if constexpr (x)
+ {
+ result += (firstCoords[XX] - secondCoords[XX]) * (firstCoords[XX] - secondCoords[XX]);
+ }
+ if constexpr (y) // NOLINT(readability-misleading-indentation): clang-tidy-9 can't handle if constexpr (https://bugs.llvm.org/show_bug.cgi?id=32203)
+ {
+ result += (firstCoords[YY] - secondCoords[YY]) * (firstCoords[YY] - secondCoords[YY]);
+ }
+ if constexpr (z) // NOLINT(readability-misleading-indentation)
+ {
+ result += (firstCoords[ZZ] - secondCoords[ZZ]) * (firstCoords[ZZ] - secondCoords[ZZ]);
+ }
+ return result; // NOLINT(readability-misleading-indentation)
+}
+
+/*! \brief Calculate average displacement between sets of points
+ *
+ * Each displacement c1[i] - c2[i] is calculated and the distances
+ * are averaged.
+ *
+ * \tparam x If true, calculate x dimension of displacement
+ * \tparam y If true, calculate y dimension of displacement
+ * \tparam z If true, calculate z dimension of displacement
+ * \param c1 First vector
+ * \param c2 Second vector
+ * \param numberOfValues
+ * \return Per-particle averaged distance
+ */
+template<bool x, bool y, bool z>
+double calcAverageDisplacement(ArrayRef<const RVec> c1, ArrayRef<const RVec> c2)
+{
+ double result = 0;
+ for (size_t i = 0; i < c1.size(); i++)
+ {
+ result += calcSingleSquaredDistance<x, y, z>(c1[i], c2[i]);
+ }
+ return result / c1.size();
+}
+
+
+//! Describes 1D MSDs, in the given dimension.
+enum class SingleDimDiffType : int
+{
+ X = 0,
+ Y,
+ Z,
+ Unused,
+ Count,
+};
+
+//! Describes 2D MSDs, in the plane normal to the given dimension.
+enum class TwoDimDiffType : int
+{
+ NormalToX = 0,
+ NormalToY,
+ NormalToZ,
+ Unused,
+ Count,
+};
+
+/*! \brief Removes jumps across periodic boundaries for currentFrame, based on the positions in
+ * previousFrame. Updates currentCoords in place.
+ */
+void removePbcJumps(ArrayRef<RVec> currentCoords, ArrayRef<const RVec> previousCoords, t_pbc* pbc)
+{
+ // There are two types of "pbc removal" in gmx msd. The first happens in the trajectoryanalysis
+ // framework, which makes molecules whole across periodic boundaries and is done
+ // automatically where the inputs support it. This lambda performs the second PBC correction, where
+ // any "jump" across periodic boundaries BETWEEN FRAMES is put back. The order of these
+ // operations is important - since the first transformation may only apply to part of a
+ // molecule (e.g., one half in/out of the box is put on one side of the box), the
+ // subsequent step needs to be applied to the molecule COM rather than individual atoms, or
+ // we'd have a clash where the per-mol PBC removal moves an atom that gets put back into
+ // it's original position by the second transformation. Therefore, this second transformation
+ // is applied *after* per molecule coordinates have been consolidated into COMs.
+ auto pbcRemover = [pbc](RVec in, RVec prev) {
+ rvec dx;
+ pbc_dx(pbc, in, prev, dx);
+ return prev + dx;
+ };
+ std::transform(
+ currentCoords.begin(), currentCoords.end(), previousCoords.begin(), currentCoords.begin(), pbcRemover);
+}
+
+//! Holds data needed for MSD calculations for a single molecule, if requested.
+struct MoleculeData
+{
+ //! Number of atoms in the molecule.
+ int atomCount = 0;
+ //! Total mass.
+ double mass = 0;
+ //! MSD accumulator and calculator for the molecule
+ MsdData msdData;
+ //! Calculated diffusion coefficient
+ real diffusionCoefficient = 0;
+};
+
+/*! \brief Handles coordinate operations for MSD calculations.
+ *
+ * Can be used to hold coordinates for individual atoms as well as molecules COMs. Handles PBC
+ * jump removal between consecutive frames.
+ *
+ * Only previous_ contains valid data between calls to buildCoordinates(), although both vectors
+ * are maintained at the correct size for the number of coordinates needed.
+ */
+class MsdCoordinateManager
+{
+public:
+ MsdCoordinateManager(const int numAtoms,
+ ArrayRef<const MoleculeData> molecules,
+ ArrayRef<const int> moleculeIndexMapping) :
+ current_(numAtoms),
+ previous_(numAtoms),
+ molecules_(molecules),
+ moleculeIndexMapping_(moleculeIndexMapping)
+ {
+ }
+ /*! \brief Prepares coordinates for the current frame.
+ *
+ * Reads in selection data, and returns an ArrayRef of the particle positions or molecule
+ * centers of mass (if the molecules input is not empty). Removes jumps across periodic
+ * boundaries based on the previous frame coordinates, except for the first frame built with
+ * builtCoordinates(), which has no previous frame as a reference.
+ *
+ * @param sel The selection object which holds coordinates
+ * @param molecules Per-molecule description of atom counts and masses
+ * @param moleculeIndexMapping Mapping of atom index to molecule index.
+ * @return The current frames coordinates in proper format.
+ */
+ ArrayRef<const RVec> buildCoordinates(const Selection& sel, t_pbc* pbc);
+
+private:
+ //! The current coordinates.
+ std::vector<RVec> current_;
+ //! The previous frame's coordinates.
+ std::vector<RVec> previous_;
+ //! Molecule data.
+ ArrayRef<const MoleculeData> molecules_;
+ //! Mapping of atom indices to molecle indices;
+ ArrayRef<const int> moleculeIndexMapping_;
+ //! Tracks if a previous frame exists to compare with for PBC handling.
+ bool containsPreviousFrame_ = false;
+};
+
+ArrayRef<const RVec> MsdCoordinateManager::buildCoordinates(const Selection& sel, t_pbc* pbc)
+{
+
+ if (molecules_.empty())
+ {
+ std::copy(sel.coordinates().begin(), sel.coordinates().end(), current_.begin());
+ }
+ else
+ {
+ // Prepare for molecule COM calculation, then sum up all positions per molecule.
+
+ std::fill(current_.begin(), current_.end(), RVec(0, 0, 0));
+ gmx::ArrayRef<const real> masses = sel.masses();
+ for (int i = 0; i < sel.posCount(); i++)
+ {
+ const int moleculeIndex = moleculeIndexMapping_[i];
+ // accumulate ri * mi, and do division at the end to minimize number of divisions.
+ current_[moleculeIndex] += RVec(sel.position(i).x()) * masses[i];
+ }
+ // Divide accumulated mass * positions to get COM.
+ std::transform(current_.begin(),
+ current_.end(),
+ molecules_.begin(),
+ current_.begin(),
+ [](const RVec& position, const MoleculeData& molecule) -> RVec {
+ return position / molecule.mass;
+ });
+ }
+
+ if (containsPreviousFrame_)
+ {
+ removePbcJumps(current_, previous_, pbc);
+ }
+ else
+ {
+ containsPreviousFrame_ = true;
+ }
+
+ // Previous is no longer needed, swap with current and return "current" coordinates which
+ // now reside in previous.
+ current_.swap(previous_);
+ return previous_;
+}
+
+
+//! Holds per-group coordinates, analysis, and results.
+struct MsdGroupData
+{
+ explicit MsdGroupData(const Selection& inputSel,
+ ArrayRef<const MoleculeData> molecules,
+ ArrayRef<const int> moleculeAtomMapping) :
+ sel(inputSel),
+ coordinateManager_(molecules.empty() ? sel.posCount() : molecules.size(), molecules, moleculeAtomMapping)
+ {
+ }
+
+ //! Selection associated with this group.
+ const Selection& sel;
+
+ //! Stored coordinates, indexed by frame then atom number.
+ std::vector<std::vector<RVec>> frames;
+
+ //! MSD result accumulator
+ MsdData msds;
+ //! Coordinate handler for the group.
+ MsdCoordinateManager coordinateManager_;
+ //! Collector for processed MSD averages per tau
+ std::vector<real> msdSums;
+ //! Fitted diffusion coefficient
+ real diffusionCoefficient = 0.0;
+ //! Uncertainty of diffusion coefficient
+ double sigma = 0.0;
+};
+
+} // namespace
+
+/*! \brief Implements the gmx msd module
+ *
+ * \todo Implement -(no)mw. Right now, all calculations are mass-weighted with -mol, and not otherwise
+ * \todo Implement -tensor for full MSD tensor calculation
+ * \todo Implement -rmcomm for total-frame COM removal
+ * \todo Implement -pdb for molecule B factors
+ * \todo Implement -maxtau option proposed at https://gitlab.com/gromacs/gromacs/-/issues/3870
+ * \todo Update help text as options are added and clarifications decided on at https://gitlab.com/gromacs/gromacs/-/issues/3869
+ */
+class Msd : public TrajectoryAnalysisModule
+{
+public:
+ Msd();
+
+ void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
+ void optionsFinished(TrajectoryAnalysisSettings* settings) override;
+ void initAfterFirstFrame(const TrajectoryAnalysisSettings& settings, const t_trxframe& fr) override;
+ void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
+ void analyzeFrame(int frameNumber,
+ const t_trxframe& frame,
+ t_pbc* pbc,
+ TrajectoryAnalysisModuleData* pdata) override;
+ void finishAnalysis(int nframes) override;
+ void writeOutput() override;
+
+private:
+ //! Selections for MSD output
+ SelectionList selections_;
+
+ //! MSD type information, for -type {x,y,z}
+ SingleDimDiffType singleDimType_ = SingleDimDiffType::Unused;
+ //! MSD type information, for -lateral {x,y,z}
+ TwoDimDiffType twoDimType_ = TwoDimDiffType::Unused;
+ //! Diffusion coefficient conversion factor
+ double diffusionCoefficientDimensionFactor_ = c_3DdiffusionDimensionFactor;
+ //! Method used to calculate MSD - changes based on dimensonality.
+ std::function<double(ArrayRef<const RVec>, ArrayRef<const RVec>)> calcMsd_ =
+ calcAverageDisplacement<true, true, true>;
+
+ //! Picoseconds between restarts
+ double trestart_ = 10.0;
+ //! Initial time
+ double t0_ = 0;
+ //! Inter-frame delta-t
+ std::optional<double> dt_ = std::nullopt;
+
+ //! First tau value to fit from for diffusion coefficient, defaults to 0.1 * max tau
+ real beginFit_ = -1.0;
+ //! Final tau value to fit to for diffusion coefficient, defaults to 0.9 * max tau
+ real endFit_ = -1.0;
+
+ //! All selection group-specific data stored here.
+ std::vector<MsdGroupData> groupData_;
+
+ //! Time of all frames.
+ std::vector<double> times_;
+ //! Taus for output - won't know the size until the end.
+ std::vector<double> taus_;
+ //! Tau indices for fitting.
+ size_t beginFitIndex_ = 0;
+ size_t endFitIndex_ = 0;
+
+ // MSD per-molecule stuff
+ //! Are we doing molecule COM-based MSDs?
+ bool molSelected_ = false;
+ //! Per molecule topology information and MSD accumulators.
+ std::vector<MoleculeData> molecules_;
+ //! Atom index -> mol index map
+ std::vector<int> moleculeIndexMappings_;
+
+ // Output stuff
+ AnalysisData msdPlotData_;
+ AnalysisData msdMoleculePlotData_;
+
+ AnalysisDataPlotSettings plotSettings_;
+ //! Per-tau MSDs for each selected group
+ std::string output_;
+ //! Per molecule diffusion coefficients if -mol is selected.
+ std::string moleculeOutput_;
+};
+
+
+Msd::Msd() = default;
+
+void Msd::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
+{
+ static const char* const desc[] = {
+ "[THISMODULE] computes the mean square displacement (MSD) of atoms from",
+ "a set of initial positions. This provides an easy way to compute",
+ "the diffusion constant using the Einstein relation.",
+ "The time between the reference points for the MSD calculation",
+ "is set with [TT]-trestart[tt].",
+ "The diffusion constant is calculated by least squares fitting a",
+ "straight line (D*t + c) through the MSD(t) from [TT]-beginfit[tt] to",
+ "[TT]-endfit[tt] (note that t is time from the reference positions,",
+ "not simulation time). An error estimate given, which is the difference",
+ "of the diffusion coefficients obtained from fits over the two halves",
+ "of the fit interval.[PAR]",
+ "There are three, mutually exclusive, options to determine different",
+ "types of mean square displacement: [TT]-type[tt], [TT]-lateral[tt]",
+ "and [TT]-ten[tt]. Option [TT]-ten[tt] writes the full MSD tensor for",
+ "each group, the order in the output is: trace xx yy zz yx zx zy.[PAR]",
+ "If [TT]-mol[tt] is set, [THISMODULE] plots the MSD for individual molecules",
+ "(including making molecules whole across periodic boundaries): ",
+ "for each individual molecule a diffusion constant is computed for ",
+ "its center of mass. The chosen index group will be split into ",
+ "molecules. With -mol, only one index group can be selected.[PAR]",
+ "The diffusion coefficient is determined by linear regression of the MSD.",
+ "When [TT]-beginfit[tt] is -1, fitting starts at 10%",
+ "and when [TT]-endfit[tt] is -1, fitting goes to 90%.",
+ "Using this option one also gets an accurate error estimate",
+ "based on the statistics between individual molecules.",
+ "Note that this diffusion coefficient and error estimate are only",
+ "accurate when the MSD is completely linear between",
+ "[TT]-beginfit[tt] and [TT]-endfit[tt].[PAR]",
+ };
+ settings->setHelpText(desc);
+
+ // Selections
+ options->addOption(SelectionOption("sel")
+ .storeVector(&selections_)
+ .required()
+ .onlyStatic()
+ .multiValue()
+ .description("Selections to compute MSDs for from the reference"));
+
+ // Select MSD type - defaults to 3D if neither option is selected.
+ EnumerationArray<SingleDimDiffType, const char*> enumTypeNames = { "x", "y", "z", "unused" };
+ EnumerationArray<TwoDimDiffType, const char*> enumLateralNames = { "x", "y", "z", "unused" };
+ options->addOption(EnumOption<SingleDimDiffType>("type")
+ .enumValue(enumTypeNames)
+ .store(&singleDimType_)
+ .defaultValue(SingleDimDiffType::Unused));
+ options->addOption(EnumOption<TwoDimDiffType>("lateral")
+ .enumValue(enumLateralNames)
+ .store(&twoDimType_)
+ .defaultValue(TwoDimDiffType::Unused));
+
+ options->addOption(DoubleOption("trestart")
+ .description("Time between restarting points in trajectory (ps)")
+ .defaultValue(10.0)
+ .store(&trestart_));
+ options->addOption(RealOption("beginfit").description("").store(&beginFit_));
+ options->addOption(RealOption("endfit").description("").store(&endFit_));
+
+ // Output options
+ options->addOption(FileNameOption("o")
+ .filetype(eftPlot)
+ .outputFile()
+ .store(&output_)
+ .defaultBasename("msdout")
+ .description("MSD output"));
+ options->addOption(
+ FileNameOption("mol")
+ .filetype(eftPlot)
+ .outputFile()
+ .store(&moleculeOutput_)
+ .storeIsSet(&molSelected_)
+ .defaultBasename("diff_mol")
+ .description("Report diffusion coefficients for each molecule in selection"));
+}
+
+void Msd::optionsFinished(TrajectoryAnalysisSettings gmx_unused* settings)
+{
+ if (singleDimType_ != SingleDimDiffType::Unused && twoDimType_ != TwoDimDiffType::Unused)
+ {
+ std::string errorMessage =
+ "Options -type and -lateral are mutually exclusive. Choose one or neither (for 3D "
+ "MSDs).";
+ GMX_THROW(InconsistentInputError(errorMessage.c_str()));
+ }
+ if (selections_.size() > 1 && molSelected_)
+ {
+ std::string errorMessage =
+ "Cannot have multiple groups selected with -sel when using -mol.";
+ GMX_THROW(InconsistentInputError(errorMessage.c_str()));
+ }
+}
+
+
+void Msd::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
+{
+ plotSettings_ = settings.plotSettings();
+
+ // Enumeration helpers for dispatching the right MSD calculation type.
+ const EnumerationArray<SingleDimDiffType, decltype(&calcAverageDisplacement<true, true, true>)>
+ oneDimensionalMsdFunctions = { &calcAverageDisplacement<true, false, false>,
+ &calcAverageDisplacement<false, true, false>,
+ &calcAverageDisplacement<false, false, true> };
+ const EnumerationArray<TwoDimDiffType, decltype(&calcAverageDisplacement<true, true, true>)>
+ twoDimensionalMsdFunctions = { calcAverageDisplacement<false, true, true>,
+ calcAverageDisplacement<true, false, true>,
+ calcAverageDisplacement<true, true, false> };
+
+ // Parse dimensionality and assign the MSD calculating function.
+ // Note if we don't hit either of these cases, we're computing 3D MSDs.
+ if (singleDimType_ != SingleDimDiffType::Unused)
+ {
+ calcMsd_ = oneDimensionalMsdFunctions[singleDimType_];
+ diffusionCoefficientDimensionFactor_ = c_1DdiffusionDimensionFactor;
+ }
+ else if (twoDimType_ != TwoDimDiffType::Unused)
+ {
+ calcMsd_ = twoDimensionalMsdFunctions[twoDimType_];
+ diffusionCoefficientDimensionFactor_ = c_2DdiffusionDimensionFactor;
+ }
+
+ // TODO validate that we have mol info and not atom only - and masses, and topology.
+ if (molSelected_)
+ {
+ Selection& sel = selections_[0];
+ const int nMol = sel.initOriginalIdsToGroup(top.mtop(), INDEX_MOL);
+
+ gmx::ArrayRef<const int> mappedIds = selections_[0].mappedIds();
+ moleculeIndexMappings_.resize(selections_[0].posCount());
+ std::copy(mappedIds.begin(), mappedIds.end(), moleculeIndexMappings_.begin());
+
+ // Precalculate each molecules mass for speeding up COM calculations.
+ ArrayRef<const real> masses = sel.masses();
+
+ molecules_.resize(nMol);
+ for (int i = 0; i < sel.posCount(); i++)
+ {
+ molecules_[mappedIds[i]].atomCount++;
+ molecules_[mappedIds[i]].mass += masses[i];
+ }
+ }
+
+ // Accumulated frames and results
+ for (const Selection& sel : selections_)
+ {
+ groupData_.emplace_back(sel, molecules_, moleculeIndexMappings_);
+ }
+}
+
+void Msd::initAfterFirstFrame(const TrajectoryAnalysisSettings gmx_unused& settings, const t_trxframe& fr)
+{
+ t0_ = std::round(fr.time);
+}
+
+
+void Msd::analyzeFrame(int gmx_unused frameNumber,
+ const t_trxframe& frame,
+ t_pbc* pbc,
+ TrajectoryAnalysisModuleData* pdata)
+{
+ const real time = std::round(frame.time);
+ // Need to populate dt on frame 2;
+ if (!dt_.has_value() && !times_.empty())
+ {
+ dt_ = time - times_[0];
+ }
+
+ // Each frame gets an entry in times, but frameTimes only updates if we're at a restart.
+ times_.push_back(time);
+
+ // Each frame will get a tau between it and frame 0, and all other frame combos should be
+ // covered by this.
+ // \todo this will no longer hold exactly when maxtau is added
+ taus_.push_back(time - times_[0]);
+
+ for (MsdGroupData& msdData : groupData_)
+ {
+ //NOLINTNEXTLINE(readability-static-accessed-through-instance)
+ const Selection& sel = pdata->parallelSelection(msdData.sel);
+
+ ArrayRef<const RVec> coords = msdData.coordinateManager_.buildCoordinates(sel, pbc);
+
+ // For each preceding frame, calculate tau and do comparison.
+ for (size_t i = 0; i < msdData.frames.size(); i++)
+ {
+ // If dt > trestart, the tau increment will be dt rather than trestart.
+ double tau = time - (t0_ + std::max<double>(trestart_, *dt_) * i);
+ int64_t tauIndex = gmx::roundToInt64(tau / *dt_);
+ msdData.msds[tauIndex].push_back(calcMsd_(coords, msdData.frames[i]));
+
+ for (size_t molInd = 0; molInd < molecules_.size(); molInd++)
+ {
+ molecules_[molInd].msdData[tauIndex].push_back(
+ calcMsd_(arrayRefFromArray(&coords[molInd], 1),
+ arrayRefFromArray(&msdData.frames[i][molInd], 1)));
+ }
+ }
+
+
+ // We only store the frame for the future if it's a restart per -trestart.
+ if (bRmod(time, t0_, trestart_))
+ {
+ msdData.frames.emplace_back(coords.begin(), coords.end());
+ }
+ }
+}
+
+//! Calculate the tau index for fitting. If userFitTau < 0, uses the default fraction of max tau.
+static size_t calculateFitIndex(const int userFitTau,
+ const double defaultTauFraction,
+ const int numTaus,
+ const double dt)
+{
+ if (userFitTau < 0)
+ {
+ return gmx::roundToInt((numTaus - 1) * defaultTauFraction);
+ }
+ return std::min<size_t>(numTaus - 1, gmx::roundToInt(static_cast<double>(userFitTau) / dt));
+}
+
+
+void Msd::finishAnalysis(int gmx_unused nframes)
+{
+ static constexpr double c_defaultStartFitIndexFraction = 0.1;
+ static constexpr double c_defaultEndFitIndexFraction = 0.9;
+ beginFitIndex_ = calculateFitIndex(beginFit_, c_defaultStartFitIndexFraction, taus_.size(), *dt_);
+ endFitIndex_ = calculateFitIndex(endFit_, c_defaultEndFitIndexFraction, taus_.size(), *dt_);
+ const int numTausForFit = 1 + endFitIndex_ - beginFitIndex_;
+
+ // These aren't used, except for correlationCoefficient, which is used to estimate error if
+ // enough points are available.
+ real b = 0.0, correlationCoefficient = 0.0, chiSquared = 0.0;
+
+ for (MsdGroupData& msdData : groupData_)
+ {
+ msdData.msdSums = msdData.msds.averageMsds();
+
+ if (numTausForFit >= 4)
+ {
+ const int halfNumTaus = numTausForFit / 2;
+ const int secondaryStartIndex = beginFitIndex_ + halfNumTaus;
+ // Split the fit in 2, and compare the results of each fit;
+ real a = 0.0, a2 = 0.0;
+ lsq_y_ax_b_xdouble(halfNumTaus,
+ &taus_[beginFitIndex_],
+ &msdData.msdSums[beginFitIndex_],
+ &a,
+ &b,
+ &correlationCoefficient,
+ &chiSquared);
+ lsq_y_ax_b_xdouble(halfNumTaus,
+ &taus_[secondaryStartIndex],
+ &msdData.msdSums[secondaryStartIndex],
+ &a2,
+ &b,
+ &correlationCoefficient,
+ &chiSquared);
+ msdData.sigma = std::abs(a - a2);
+ }
+ lsq_y_ax_b_xdouble(numTausForFit,
+ &taus_[beginFitIndex_],
+ &msdData.msdSums[beginFitIndex_],
+ &msdData.diffusionCoefficient,
+ &b,
+ &correlationCoefficient,
+ &chiSquared);
+ msdData.diffusionCoefficient *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
+ msdData.sigma *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
+ }
+
+ for (MoleculeData& molecule : molecules_)
+ {
+ std::vector<real> msds = molecule.msdData.averageMsds();
+ lsq_y_ax_b_xdouble(numTausForFit,
+ &taus_[beginFitIndex_],
+ &msds[beginFitIndex_],
+ &molecule.diffusionCoefficient,
+ &b,
+ &correlationCoefficient,
+ &chiSquared);
+ molecule.diffusionCoefficient *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
+ }
+}
+
+void Msd::writeOutput()
+{
+ // AnalysisData currently doesn't support changing column counts after analysis has started.
+ // We can't determine the number of tau values until the trajectory is fully read, so analysis
+ // data construction and plotting are done here.
+ AnalysisDataPlotModulePointer msdPlotModule(new AnalysisDataPlotModule(plotSettings_));
+ msdPlotModule->setFileName(output_);
+ msdPlotModule->setTitle("Mean Squared Displacement");
+ msdPlotModule->setXLabel("tau (ps)");
+ msdPlotModule->setYLabel(R"(MSD (nm\\S2\\N))");
+ msdPlotModule->setYFormat(10, 6, 'g');
+ for (const auto& group : groupData_)
+ {
+ const real D = group.diffusionCoefficient;
+ if (D > 0.01 && D < 1e4)
+ {
+ msdPlotModule->appendLegend(formatString(
+ "D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)", group.sel.name(), D, group.sigma));
+ }
+ else
+ {
+ msdPlotModule->appendLegend(formatString(
+ "D[%10s] = %.4g (+/- %.4f) (1e-5 cm^2/s)", group.sel.name(), D, group.sigma));
+ }
+ }
+ msdPlotData_.addModule(msdPlotModule);
+ msdPlotData_.setDataSetCount(groupData_.size());
+ for (size_t i = 0; i < groupData_.size(); i++)
+ {
+ msdPlotData_.setColumnCount(i, 1);
+ }
+ AnalysisDataHandle dh = msdPlotData_.startData({});
+ for (size_t tauIndex = 0; tauIndex < taus_.size(); tauIndex++)
+ {
+ dh.startFrame(tauIndex, taus_[tauIndex]);
+ for (size_t dataSetIndex = 0; dataSetIndex < groupData_.size(); dataSetIndex++)
+ {
+ dh.selectDataSet(dataSetIndex);
+ dh.setPoint(0, groupData_[dataSetIndex].msdSums[tauIndex]);
+ }
+ dh.finishFrame();
+ }
+ dh.finishData();
+
+ if (molSelected_)
+ {
+ AnalysisDataPlotModulePointer molPlotModule(new AnalysisDataPlotModule(plotSettings_));
+ molPlotModule->setFileName(moleculeOutput_);
+ molPlotModule->setTitle("Mean Squared Displacement / Molecule");
+ molPlotModule->setXLabel("Molecule");
+ molPlotModule->setYLabel("D(1e-5 cm^2/s)");
+ molPlotModule->setYFormat(10, 0, 'g');
+ msdMoleculePlotData_.addModule(molPlotModule);
+ msdMoleculePlotData_.setDataSetCount(1);
+ msdMoleculePlotData_.setColumnCount(0, 1);
+ AnalysisDataHandle molDh = msdMoleculePlotData_.startData({});
+ for (size_t moleculeIndex = 0; moleculeIndex < molecules_.size(); moleculeIndex++)
+ {
+ molDh.startFrame(moleculeIndex, moleculeIndex);
+ molDh.setPoint(0, molecules_[moleculeIndex].diffusionCoefficient);
+ molDh.finishFrame();
+ }
+ molDh.finishData();
+ }
+}
+
+
+const char MsdInfo::name[] = "msd";
+const char MsdInfo::shortDescription[] = "Compute mean squared displacements";
+TrajectoryAnalysisModulePointer MsdInfo::create()
+{
+ return TrajectoryAnalysisModulePointer(new Msd);
+}
+
+
+} // namespace analysismodules
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020,2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Declares trajectory analysis module for Mean Squared Displacement calculations.
+ *
+ * \author Kevin Boyd <kevin44boyd@gmail.com>
+ * \ingroup module_trajectoryanalysis
+ */
+#ifndef GMX_TRAJECTORYANALYSIS_MODULES_MSD_H
+#define GMX_TRAJECTORYANALYSIS_MODULES_MSD_H
+
+#include "gromacs/trajectoryanalysis/analysismodule.h"
+
+namespace gmx::analysismodules
+{
+
+class MsdInfo
+{
+public:
+ static const char name[];
+ static const char shortDescription[];
+ static TrajectoryAnalysisModulePointer create();
+};
+
+} // namespace gmx::analysismodules
+
+#endif // GMX_TRAJECTORYANALYSIS_MODULES_MSD_H
# This file is part of the GROMACS molecular simulation package.
#
# Copyright (c) 2010,2012,2013,2014,2015 by the GROMACS development team.
-# Copyright (c) 2016,2017,2018,2019,2020, by the GROMACS development team, led by
+# Copyright (c) 2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
distance.cpp
extract_cluster.cpp
freevolume.cpp
+ msd.cpp
pairdist.cpp
rdf.cpp
sasa.cpp
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020,2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests for functionality of the "msd" trajectory analysis module.
+ *
+ * \author Kevin Boyd <kevin44boyd@gmail.com>
+ * \ingroup module_trajectoryanalysis
+ */
+#include "gmxpre.h"
+
+#include "gromacs/trajectoryanalysis/modules/msd.h"
+
+#include <gtest/gtest.h>
+
+#include "gromacs/gmxpreprocess/grompp.h"
+#include "gromacs/utility/path.h"
+#include "gromacs/utility/stringutil.h"
+#include "gromacs/utility/textstream.h"
+
+#include "testutils/cmdlinetest.h"
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+#include "testutils/testfilemanager.h"
+#include "testutils/textblockmatchers.h"
+
+#include "moduletest.h"
+
+namespace gmx::test
+{
+
+namespace
+{
+
+using gmx::test::CommandLine;
+
+/*! \brief Returns whether or not we care about a header line in an xvg file, for matching purposes.
+ *
+ * \todo This is mostly taken from xvgtest.cpp. We could probably create some modular checker
+ * functionality where each line is compared against a set of subcheckers automatically. Then we
+ * could build matchers like this out of the modular components.
+ */
+bool isRelevantXvgHeader(const std::string& line)
+{
+ return startsWith(line, "@")
+ && (contains(line, " title ") || contains(line, " subtitle ") || contains(line, " label ")
+ || contains(line, "@TYPE ") || contains(line, " legend \""));
+}
+
+//! Helper function to check a single xvg value in a sequence.
+void checkXvgDataPoint(TestReferenceChecker* checker, const std::string& value)
+{
+ checker->checkRealFromString(value, nullptr);
+}
+
+/*! \brief MsdMatcher is effectively an extension of XvgMatcher for gmx msd results.
+ *
+ * In addition to the usual fields XvgMatcher checks, MsdMatcher checks for properly reported
+ * diffusion coefficients.
+ */
+class MsdMatcher : public ITextBlockMatcher
+{
+public:
+ MsdMatcher() = default;
+
+ void checkStream(TextInputStream* stream, TestReferenceChecker* checker) override
+ {
+ checker->setDefaultTolerance(
+ gmx::test::FloatingPointTolerance(gmx::test::absoluteTolerance(1.0e-5)));
+ TestReferenceChecker dCoefficientChecker(
+ checker->checkCompound("XvgLegend", "DiffusionCoefficient"));
+ TestReferenceChecker legendChecker(checker->checkCompound("XvgLegend", "Legend"));
+ TestReferenceChecker dataChecker(checker->checkCompound("XvgData", "Data"));
+ std::string line;
+ int rowCount = 0;
+ while (stream->readLine(&line))
+ {
+ // Legend and titles.
+ if (isRelevantXvgHeader(line))
+ {
+ legendChecker.checkString(stripString(line.substr(1)), nullptr);
+ continue;
+ }
+ // All other comment lines we don't care about.
+ if (startsWith(line, "#") || startsWith(line, "@"))
+ {
+ continue;
+ }
+ // Actual data.
+ const std::vector<std::string> columns = splitString(line);
+ const std::string id = formatString("Row%d", rowCount);
+ dataChecker.checkSequence(columns.begin(), columns.end(), id.c_str(), &checkXvgDataPoint);
+ rowCount++;
+ }
+ }
+};
+
+class MsdMatch : public ITextBlockMatcherSettings
+{
+public:
+ [[nodiscard]] TextBlockMatcherPointer createMatcher() const override
+ {
+ return std::make_unique<MsdMatcher>();
+ }
+};
+
+class MsdModuleTest : public gmx::test::TrajectoryAnalysisModuleTestFixture<gmx::analysismodules::MsdInfo>
+{
+public:
+ MsdModuleTest() { setOutputFile("-o", "msd.xvg", MsdMatch()); }
+ // Creates a TPR for the given starting structure and topology. Builds an mdp in place prior
+ // to calling grompp. sets the -s input to the generated tpr
+ void createTpr(const std::string& structure, const std::string& topology, const std::string& index)
+ {
+ std::string tpr = fileManager().getTemporaryFilePath(".tpr");
+ std::string mdp = fileManager().getTemporaryFilePath(".mdp");
+ FILE* fp = fopen(mdp.c_str(), "w");
+ fprintf(fp, "cutoff-scheme = verlet\n");
+ fprintf(fp, "rcoulomb = 0.85\n");
+ fprintf(fp, "rvdw = 0.85\n");
+ fprintf(fp, "rlist = 0.85\n");
+ fclose(fp);
+
+ // Prepare a .tpr file
+ CommandLine caller;
+ const auto* simDB = gmx::test::TestFileManager::getTestSimulationDatabaseDirectory();
+ caller.append("grompp");
+ caller.addOption("-maxwarn", 0);
+ caller.addOption("-f", mdp.c_str());
+ std::string gro = gmx::Path::join(simDB, structure);
+ caller.addOption("-c", gro.c_str());
+ std::string top = gmx::Path::join(simDB, topology);
+ caller.addOption("-p", top.c_str());
+ std::string ndx = gmx::Path::join(simDB, index);
+ caller.addOption("-n", ndx.c_str());
+ caller.addOption("-o", tpr.c_str());
+ ASSERT_EQ(0, gmx_grompp(caller.argc(), caller.argv()));
+
+ // setInputFile() doesn't like the temporary tpr path.
+ CommandLine& cmdline = commandLine();
+ cmdline.addOption("-s", tpr.c_str());
+ }
+
+ // Convenience function to set input trajectory, tpr, and index, if all of the input files
+ // share a common prefix.
+ void setAllInputs(const std::string& prefix)
+ {
+ setInputFile("-f", prefix + ".xtc");
+ setInputFile("-n", prefix + ".ndx");
+ createTpr(prefix + ".gro", prefix + ".top", prefix + ".ndx");
+ }
+};
+
+// ----------------------------------------------
+// These tests check the basic MSD and diffusion coefficient reporting capabilities on toy systems.
+// Trestart is set to larger than the size of the trajectory so that all frames are only compared
+// against the first frame and diffusion coefficients can be predicted exactly.
+// ----------------------------------------------
+
+// for 3D, (8 + 4 + 0) / 3 should yield 4 cm^2 / s
+TEST_F(MsdModuleTest, threeDimensionalDiffusion)
+{
+ setInputFile("-f", "msd_traj.xtc");
+ setInputFile("-s", "msd_coords.gro");
+ setInputFile("-n", "msd.ndx");
+ const char* const cmdline[] = {
+ "-trestart",
+ "200",
+ "-sel",
+ "0",
+ };
+ runTest(CommandLine(cmdline));
+}
+
+// for lateral z, (8 + 4) / 2 should yield 6 cm^2 /s
+TEST_F(MsdModuleTest, twoDimensionalDiffusion)
+{
+ setInputFile("-f", "msd_traj.xtc");
+ setInputFile("-s", "msd_coords.gro");
+ setInputFile("-n", "msd.ndx");
+ const char* const cmdline[] = { "-trestart", "200", "-lateral", "z", "-sel", "0" };
+ runTest(CommandLine(cmdline));
+}
+
+// for type x, should yield 8 cm^2 / s
+TEST_F(MsdModuleTest, oneDimensionalDiffusion)
+{
+ setInputFile("-f", "msd_traj.xtc");
+ setInputFile("-s", "msd_coords.gro");
+ setInputFile("-n", "msd.ndx");
+ const char* const cmdline[] = { "-trestart", "200", "-type", "x", "-sel", "0" };
+ runTest(CommandLine(cmdline));
+}
+
+// -------------------------------------------------------------------------
+// These tests operate on a more realistic trajectory, with a solvated protein,
+// with 20 frames at a 2 ps dt. Note that this box is also non-square, so we're validating
+// non-trivial pbc removal.
+
+// Group 1 = protein, 2 = all waters, 3 = a subset of 6 water molecules
+// ------------------------------------------------------------------------
+TEST_F(MsdModuleTest, multipleGroupsWork)
+{
+ setAllInputs("alanine_vsite_solvated");
+ // Restart every frame, select protein and water separately. Note that the reported diffusion
+ // coefficient for protein is not well-sampled and doesn't correspond to anything physical.
+ const char* const cmdline[] = { "-trestart", "2", "-sel", "1;2" };
+ runTest(CommandLine(cmdline));
+}
+
+TEST_F(MsdModuleTest, trestartLessThanDt)
+{
+ setAllInputs("alanine_vsite_solvated");
+ const char* const cmdline[] = { "-trestart", "1", "-sel", "2" };
+ runTest(CommandLine(cmdline));
+}
+
+TEST_F(MsdModuleTest, trestartGreaterThanDt)
+{
+ setAllInputs("alanine_vsite_solvated");
+ const char* const cmdline[] = { "-trestart", "10", "-sel", "2" };
+ runTest(CommandLine(cmdline));
+}
+
+TEST_F(MsdModuleTest, molTest)
+{
+ setAllInputs("alanine_vsite_solvated");
+ setOutputFile("-mol", "diff_mol.xvg", MsdMatch());
+ const char* const cmdline[] = { "-trestart", "10", "-sel", "3" };
+ runTest(CommandLine(cmdline));
+}
+
+TEST_F(MsdModuleTest, beginFit)
+{
+ setAllInputs("alanine_vsite_solvated");
+ const char* const cmdline[] = { "-trestart", "2", "-sel", "3", "-beginfit", "20" };
+ runTest(CommandLine(cmdline));
+}
+
+TEST_F(MsdModuleTest, endFit)
+{
+ setAllInputs("alanine_vsite_solvated");
+ const char* const cmdline[] = { "-trestart", "2", "-sel", "3", "-endfit", "25" };
+ runTest(CommandLine(cmdline));
+}
+
+TEST_F(MsdModuleTest, notEnoughPointsForFitErrorEstimate)
+{
+ setAllInputs("alanine_vsite_solvated");
+ const char* const cmdline[] = { "-trestart", "2", "-beginfit", "5", "-endfit",
+ "9", "-lateral", "x", "-sel", "all" };
+ runTest(CommandLine(cmdline));
+}
+
+} // namespace
+
+} // namespace gmx::test
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <String Name="CommandLine">-trestart 2 -sel 3 -beginfit 20</String>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[some_water_subset] = 4.2232 (+/- 0.4722) (1e-5 cm^2/s)"</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">2</Int>
+ <Real>0.000</Real>
+ <Real>0</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">2</Int>
+ <Real>2.000</Real>
+ <Real>0.0680311</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">2</Int>
+ <Real>4.000</Real>
+ <Real>0.126496</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">2</Int>
+ <Real>6.000</Real>
+ <Real>0.185337</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">2</Int>
+ <Real>8.000</Real>
+ <Real>0.239227</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">2</Int>
+ <Real>10.000</Real>
+ <Real>0.278451</Real>
+ </Sequence>
+ <Sequence Name="Row6">
+ <Int Name="Length">2</Int>
+ <Real>12.000</Real>
+ <Real>0.324372</Real>
+ </Sequence>
+ <Sequence Name="Row7">
+ <Int Name="Length">2</Int>
+ <Real>14.000</Real>
+ <Real>0.370175</Real>
+ </Sequence>
+ <Sequence Name="Row8">
+ <Int Name="Length">2</Int>
+ <Real>16.000</Real>
+ <Real>0.43655</Real>
+ </Sequence>
+ <Sequence Name="Row9">
+ <Int Name="Length">2</Int>
+ <Real>18.000</Real>
+ <Real>0.499879</Real>
+ </Sequence>
+ <Sequence Name="Row10">
+ <Int Name="Length">2</Int>
+ <Real>20.000</Real>
+ <Real>0.553128</Real>
+ </Sequence>
+ <Sequence Name="Row11">
+ <Int Name="Length">2</Int>
+ <Real>22.000</Real>
+ <Real>0.582381</Real>
+ </Sequence>
+ <Sequence Name="Row12">
+ <Int Name="Length">2</Int>
+ <Real>24.000</Real>
+ <Real>0.627822</Real>
+ </Sequence>
+ <Sequence Name="Row13">
+ <Int Name="Length">2</Int>
+ <Real>26.000</Real>
+ <Real>0.681663</Real>
+ </Sequence>
+ <Sequence Name="Row14">
+ <Int Name="Length">2</Int>
+ <Real>28.000</Real>
+ <Real>0.770022</Real>
+ </Sequence>
+ <Sequence Name="Row15">
+ <Int Name="Length">2</Int>
+ <Real>30.000</Real>
+ <Real>0.835282</Real>
+ </Sequence>
+ <Sequence Name="Row16">
+ <Int Name="Length">2</Int>
+ <Real>32.000</Real>
+ <Real>0.908245</Real>
+ </Sequence>
+ <Sequence Name="Row17">
+ <Int Name="Length">2</Int>
+ <Real>34.000</Real>
+ <Real>0.90827</Real>
+ </Sequence>
+ <Sequence Name="Row18">
+ <Int Name="Length">2</Int>
+ <Real>36.000</Real>
+ <Real>0.890272</Real>
+ </Sequence>
+ <Sequence Name="Row19">
+ <Int Name="Length">2</Int>
+ <Real>38.000</Real>
+ <Real>0.856884</Real>
+ </Sequence>
+ <Sequence Name="Row20">
+ <Int Name="Length">2</Int>
+ <Real>40.000</Real>
+ <Real>0.913993</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ </OutputFiles>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <String Name="CommandLine">-trestart 2 -sel 3 -endfit 25</String>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[some_water_subset] = 4.2360 (+/- 0.4344) (1e-5 cm^2/s)"</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">2</Int>
+ <Real>0.000</Real>
+ <Real>0</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">2</Int>
+ <Real>2.000</Real>
+ <Real>0.0680311</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">2</Int>
+ <Real>4.000</Real>
+ <Real>0.126496</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">2</Int>
+ <Real>6.000</Real>
+ <Real>0.185337</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">2</Int>
+ <Real>8.000</Real>
+ <Real>0.239227</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">2</Int>
+ <Real>10.000</Real>
+ <Real>0.278451</Real>
+ </Sequence>
+ <Sequence Name="Row6">
+ <Int Name="Length">2</Int>
+ <Real>12.000</Real>
+ <Real>0.324372</Real>
+ </Sequence>
+ <Sequence Name="Row7">
+ <Int Name="Length">2</Int>
+ <Real>14.000</Real>
+ <Real>0.370175</Real>
+ </Sequence>
+ <Sequence Name="Row8">
+ <Int Name="Length">2</Int>
+ <Real>16.000</Real>
+ <Real>0.43655</Real>
+ </Sequence>
+ <Sequence Name="Row9">
+ <Int Name="Length">2</Int>
+ <Real>18.000</Real>
+ <Real>0.499879</Real>
+ </Sequence>
+ <Sequence Name="Row10">
+ <Int Name="Length">2</Int>
+ <Real>20.000</Real>
+ <Real>0.553128</Real>
+ </Sequence>
+ <Sequence Name="Row11">
+ <Int Name="Length">2</Int>
+ <Real>22.000</Real>
+ <Real>0.582381</Real>
+ </Sequence>
+ <Sequence Name="Row12">
+ <Int Name="Length">2</Int>
+ <Real>24.000</Real>
+ <Real>0.627822</Real>
+ </Sequence>
+ <Sequence Name="Row13">
+ <Int Name="Length">2</Int>
+ <Real>26.000</Real>
+ <Real>0.681663</Real>
+ </Sequence>
+ <Sequence Name="Row14">
+ <Int Name="Length">2</Int>
+ <Real>28.000</Real>
+ <Real>0.770022</Real>
+ </Sequence>
+ <Sequence Name="Row15">
+ <Int Name="Length">2</Int>
+ <Real>30.000</Real>
+ <Real>0.835282</Real>
+ </Sequence>
+ <Sequence Name="Row16">
+ <Int Name="Length">2</Int>
+ <Real>32.000</Real>
+ <Real>0.908245</Real>
+ </Sequence>
+ <Sequence Name="Row17">
+ <Int Name="Length">2</Int>
+ <Real>34.000</Real>
+ <Real>0.90827</Real>
+ </Sequence>
+ <Sequence Name="Row18">
+ <Int Name="Length">2</Int>
+ <Real>36.000</Real>
+ <Real>0.890272</Real>
+ </Sequence>
+ <Sequence Name="Row19">
+ <Int Name="Length">2</Int>
+ <Real>38.000</Real>
+ <Real>0.856884</Real>
+ </Sequence>
+ <Sequence Name="Row20">
+ <Int Name="Length">2</Int>
+ <Real>40.000</Real>
+ <Real>0.913993</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ </OutputFiles>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <String Name="CommandLine">-trestart 10 -sel 3</String>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[some_water_subset] = 3.9792 (+/- 0.4600) (1e-5 cm^2/s)"</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">2</Int>
+ <Real>0.000</Real>
+ <Real>0</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">2</Int>
+ <Real>2.000</Real>
+ <Real>0.0708522</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">2</Int>
+ <Real>4.000</Real>
+ <Real>0.0947372</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">2</Int>
+ <Real>6.000</Real>
+ <Real>0.173871</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">2</Int>
+ <Real>8.000</Real>
+ <Real>0.224929</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">2</Int>
+ <Real>10.000</Real>
+ <Real>0.274871</Real>
+ </Sequence>
+ <Sequence Name="Row6">
+ <Int Name="Length">2</Int>
+ <Real>12.000</Real>
+ <Real>0.311595</Real>
+ </Sequence>
+ <Sequence Name="Row7">
+ <Int Name="Length">2</Int>
+ <Real>14.000</Real>
+ <Real>0.356893</Real>
+ </Sequence>
+ <Sequence Name="Row8">
+ <Int Name="Length">2</Int>
+ <Real>16.000</Real>
+ <Real>0.459085</Real>
+ </Sequence>
+ <Sequence Name="Row9">
+ <Int Name="Length">2</Int>
+ <Real>18.000</Real>
+ <Real>0.46479</Real>
+ </Sequence>
+ <Sequence Name="Row10">
+ <Int Name="Length">2</Int>
+ <Real>20.000</Real>
+ <Real>0.545767</Real>
+ </Sequence>
+ <Sequence Name="Row11">
+ <Int Name="Length">2</Int>
+ <Real>22.000</Real>
+ <Real>0.505186</Real>
+ </Sequence>
+ <Sequence Name="Row12">
+ <Int Name="Length">2</Int>
+ <Real>24.000</Real>
+ <Real>0.574966</Real>
+ </Sequence>
+ <Sequence Name="Row13">
+ <Int Name="Length">2</Int>
+ <Real>26.000</Real>
+ <Real>0.66103</Real>
+ </Sequence>
+ <Sequence Name="Row14">
+ <Int Name="Length">2</Int>
+ <Real>28.000</Real>
+ <Real>0.64775</Real>
+ </Sequence>
+ <Sequence Name="Row15">
+ <Int Name="Length">2</Int>
+ <Real>30.000</Real>
+ <Real>0.749603</Real>
+ </Sequence>
+ <Sequence Name="Row16">
+ <Int Name="Length">2</Int>
+ <Real>32.000</Real>
+ <Real>0.804233</Real>
+ </Sequence>
+ <Sequence Name="Row17">
+ <Int Name="Length">2</Int>
+ <Real>34.000</Real>
+ <Real>0.828547</Real>
+ </Sequence>
+ <Sequence Name="Row18">
+ <Int Name="Length">2</Int>
+ <Real>36.000</Real>
+ <Real>0.909855</Real>
+ </Sequence>
+ <Sequence Name="Row19">
+ <Int Name="Length">2</Int>
+ <Real>38.000</Real>
+ <Real>0.721128</Real>
+ </Sequence>
+ <Sequence Name="Row20">
+ <Int Name="Length">2</Int>
+ <Real>40.000</Real>
+ <Real>0.899906</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ <File Name="-mol">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement / Molecule"</String>
+ <String>xaxis label "Molecule"</String>
+ <String>yaxis label "D(1e-5 cm^2/s)"</String>
+ <String>TYPE xy</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">2</Int>
+ <Real>0.000</Real>
+ <Real>4</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">2</Int>
+ <Real>1.000</Real>
+ <Real>2</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">2</Int>
+ <Real>2.000</Real>
+ <Real>3</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">2</Int>
+ <Real>3.000</Real>
+ <Real>3</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">2</Int>
+ <Real>4.000</Real>
+ <Real>1e+01</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">2</Int>
+ <Real>5.000</Real>
+ <Real>1</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ </OutputFiles>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <String Name="CommandLine">-trestart 2 -sel 1;2</String>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[ Protein] = 0.0005734 (+/- 1.2675) (1e-5 cm^2/s)"</String>
+ <String>s1 legend "D[ Water] = 5.4502 (+/- 0.2739) (1e-5 cm^2/s)"</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">3</Int>
+ <Real>0.000</Real>
+ <Real>0</Real>
+ <Real>0</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">3</Int>
+ <Real>2.000</Real>
+ <Real>0.0218917</Real>
+ <Real>0.0722943</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">3</Int>
+ <Real>4.000</Real>
+ <Real>0.0376418</Real>
+ <Real>0.13585</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">3</Int>
+ <Real>6.000</Real>
+ <Real>0.0505145</Real>
+ <Real>0.19875</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">3</Int>
+ <Real>8.000</Real>
+ <Real>0.0618931</Real>
+ <Real>0.262333</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">3</Int>
+ <Real>10.000</Real>
+ <Real>0.0757011</Real>
+ <Real>0.324822</Real>
+ </Sequence>
+ <Sequence Name="Row6">
+ <Int Name="Length">3</Int>
+ <Real>12.000</Real>
+ <Real>0.0866048</Real>
+ <Real>0.388126</Real>
+ </Sequence>
+ <Sequence Name="Row7">
+ <Int Name="Length">3</Int>
+ <Real>14.000</Real>
+ <Real>0.0937544</Real>
+ <Real>0.45144</Real>
+ </Sequence>
+ <Sequence Name="Row8">
+ <Int Name="Length">3</Int>
+ <Real>16.000</Real>
+ <Real>0.0982917</Real>
+ <Real>0.515772</Real>
+ </Sequence>
+ <Sequence Name="Row9">
+ <Int Name="Length">3</Int>
+ <Real>18.000</Real>
+ <Real>0.0961307</Real>
+ <Real>0.581858</Real>
+ </Sequence>
+ <Sequence Name="Row10">
+ <Int Name="Length">3</Int>
+ <Real>20.000</Real>
+ <Real>0.0930816</Real>
+ <Real>0.645979</Real>
+ </Sequence>
+ <Sequence Name="Row11">
+ <Int Name="Length">3</Int>
+ <Real>22.000</Real>
+ <Real>0.0958314</Real>
+ <Real>0.708053</Real>
+ </Sequence>
+ <Sequence Name="Row12">
+ <Int Name="Length">3</Int>
+ <Real>24.000</Real>
+ <Real>0.102653</Real>
+ <Real>0.772597</Real>
+ </Sequence>
+ <Sequence Name="Row13">
+ <Int Name="Length">3</Int>
+ <Real>26.000</Real>
+ <Real>0.0928015</Real>
+ <Real>0.84093</Real>
+ </Sequence>
+ <Sequence Name="Row14">
+ <Int Name="Length">3</Int>
+ <Real>28.000</Real>
+ <Real>0.0763908</Real>
+ <Real>0.908325</Real>
+ </Sequence>
+ <Sequence Name="Row15">
+ <Int Name="Length">3</Int>
+ <Real>30.000</Real>
+ <Real>0.0721617</Real>
+ <Real>0.975539</Real>
+ </Sequence>
+ <Sequence Name="Row16">
+ <Int Name="Length">3</Int>
+ <Real>32.000</Real>
+ <Real>0.0722768</Real>
+ <Real>1.04299</Real>
+ </Sequence>
+ <Sequence Name="Row17">
+ <Int Name="Length">3</Int>
+ <Real>34.000</Real>
+ <Real>0.0506321</Real>
+ <Real>1.11211</Real>
+ </Sequence>
+ <Sequence Name="Row18">
+ <Int Name="Length">3</Int>
+ <Real>36.000</Real>
+ <Real>0.0367256</Real>
+ <Real>1.19384</Real>
+ </Sequence>
+ <Sequence Name="Row19">
+ <Int Name="Length">3</Int>
+ <Real>38.000</Real>
+ <Real>0.0418762</Real>
+ <Real>1.27086</Real>
+ </Sequence>
+ <Sequence Name="Row20">
+ <Int Name="Length">3</Int>
+ <Real>40.000</Real>
+ <Real>0.052471</Real>
+ <Real>1.34173</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ </OutputFiles>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <String Name="CommandLine">-trestart 2 -beginfit 5 -endfit 9 -lateral x -sel all</String>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[ all] = 5.0364 (+/- 0.0000) (1e-5 cm^2/s)"</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">2</Int>
+ <Real>0.000</Real>
+ <Real>0</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">2</Int>
+ <Real>2.000</Real>
+ <Real>0.0469217</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">2</Int>
+ <Real>4.000</Real>
+ <Real>0.0883056</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">2</Int>
+ <Real>6.000</Real>
+ <Real>0.128427</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">2</Int>
+ <Real>8.000</Real>
+ <Real>0.168888</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">2</Int>
+ <Real>10.000</Real>
+ <Real>0.209869</Real>
+ </Sequence>
+ <Sequence Name="Row6">
+ <Int Name="Length">2</Int>
+ <Real>12.000</Real>
+ <Real>0.251233</Real>
+ </Sequence>
+ <Sequence Name="Row7">
+ <Int Name="Length">2</Int>
+ <Real>14.000</Real>
+ <Real>0.292567</Real>
+ </Sequence>
+ <Sequence Name="Row8">
+ <Int Name="Length">2</Int>
+ <Real>16.000</Real>
+ <Real>0.334993</Real>
+ </Sequence>
+ <Sequence Name="Row9">
+ <Int Name="Length">2</Int>
+ <Real>18.000</Real>
+ <Real>0.378822</Real>
+ </Sequence>
+ <Sequence Name="Row10">
+ <Int Name="Length">2</Int>
+ <Real>20.000</Real>
+ <Real>0.420988</Real>
+ </Sequence>
+ <Sequence Name="Row11">
+ <Int Name="Length">2</Int>
+ <Real>22.000</Real>
+ <Real>0.461147</Real>
+ </Sequence>
+ <Sequence Name="Row12">
+ <Int Name="Length">2</Int>
+ <Real>24.000</Real>
+ <Real>0.503649</Real>
+ </Sequence>
+ <Sequence Name="Row13">
+ <Int Name="Length">2</Int>
+ <Real>26.000</Real>
+ <Real>0.548967</Real>
+ </Sequence>
+ <Sequence Name="Row14">
+ <Int Name="Length">2</Int>
+ <Real>28.000</Real>
+ <Real>0.594155</Real>
+ </Sequence>
+ <Sequence Name="Row15">
+ <Int Name="Length">2</Int>
+ <Real>30.000</Real>
+ <Real>0.640334</Real>
+ </Sequence>
+ <Sequence Name="Row16">
+ <Int Name="Length">2</Int>
+ <Real>32.000</Real>
+ <Real>0.684822</Real>
+ </Sequence>
+ <Sequence Name="Row17">
+ <Int Name="Length">2</Int>
+ <Real>34.000</Real>
+ <Real>0.729676</Real>
+ </Sequence>
+ <Sequence Name="Row18">
+ <Int Name="Length">2</Int>
+ <Real>36.000</Real>
+ <Real>0.784167</Real>
+ </Sequence>
+ <Sequence Name="Row19">
+ <Int Name="Length">2</Int>
+ <Real>38.000</Real>
+ <Real>0.843655</Real>
+ </Sequence>
+ <Sequence Name="Row20">
+ <Int Name="Length">2</Int>
+ <Real>40.000</Real>
+ <Real>0.886678</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ </OutputFiles>
+</ReferenceData>
<?xml version="1.0"?>
<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
<ReferenceData>
+ <String Name="CommandLine">-trestart 200 -type x -sel 0</String>
<OutputFiles Name="Files">
<File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
<XvgLegend Name="Legend">
- <String Name="XvgLegend"><![CDATA[
-title "Mean Square Displacement"
-xaxis label "Time (ps)"
-yaxis label "MSD (nm\S2\N)"
-TYPE xy
-]]></String>
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[ particles] = 8.0000 (+/- 0.0000) (1e-5 cm^2/s)"</String>
</XvgLegend>
<XvgData Name="Data">
<Sequence Name="Row0">
<Int Name="Length">2</Int>
- <Real>0</Real>
+ <Real>0.000</Real>
<Real>0</Real>
</Sequence>
<Sequence Name="Row1">
<Int Name="Length">2</Int>
- <Real>1</Real>
+ <Real>1.000</Real>
<Real>0.016</Real>
</Sequence>
<Sequence Name="Row2">
<Int Name="Length">2</Int>
- <Real>2</Real>
+ <Real>2.000</Real>
<Real>0.032</Real>
</Sequence>
<Sequence Name="Row3">
<Int Name="Length">2</Int>
- <Real>3</Real>
+ <Real>3.000</Real>
<Real>0.048</Real>
</Sequence>
<Sequence Name="Row4">
<Int Name="Length">2</Int>
- <Real>4</Real>
+ <Real>4.000</Real>
<Real>0.064</Real>
</Sequence>
<Sequence Name="Row5">
<Int Name="Length">2</Int>
- <Real>5</Real>
+ <Real>5.000</Real>
<Real>0.08</Real>
</Sequence>
<Sequence Name="Row6">
<Int Name="Length">2</Int>
- <Real>6</Real>
+ <Real>6.000</Real>
<Real>0.096</Real>
</Sequence>
<Sequence Name="Row7">
<Int Name="Length">2</Int>
- <Real>7</Real>
+ <Real>7.000</Real>
<Real>0.112</Real>
</Sequence>
<Sequence Name="Row8">
<Int Name="Length">2</Int>
- <Real>8</Real>
+ <Real>8.000</Real>
<Real>0.128</Real>
</Sequence>
<Sequence Name="Row9">
<Int Name="Length">2</Int>
- <Real>9</Real>
+ <Real>9.000</Real>
<Real>0.144</Real>
</Sequence>
</XvgData>
<?xml version="1.0"?>
<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
<ReferenceData>
+ <String Name="CommandLine">-trestart 200 -sel 0</String>
<OutputFiles Name="Files">
<File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
<XvgLegend Name="Legend">
- <String Name="XvgLegend"><![CDATA[
-title "Mean Square Displacement"
-xaxis label "Time (ps)"
-yaxis label "MSD (nm\S2\N)"
-TYPE xy
-]]></String>
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[ particles] = 4.0000 (+/- 0.0000) (1e-5 cm^2/s)"</String>
</XvgLegend>
<XvgData Name="Data">
<Sequence Name="Row0">
<Int Name="Length">2</Int>
- <Real>0</Real>
+ <Real>0.000</Real>
<Real>0</Real>
</Sequence>
<Sequence Name="Row1">
<Int Name="Length">2</Int>
- <Real>1</Real>
+ <Real>1.000</Real>
<Real>0.024</Real>
</Sequence>
<Sequence Name="Row2">
<Int Name="Length">2</Int>
- <Real>2</Real>
+ <Real>2.000</Real>
<Real>0.048</Real>
</Sequence>
<Sequence Name="Row3">
<Int Name="Length">2</Int>
- <Real>3</Real>
+ <Real>3.000</Real>
<Real>0.072</Real>
</Sequence>
<Sequence Name="Row4">
<Int Name="Length">2</Int>
- <Real>4</Real>
+ <Real>4.000</Real>
<Real>0.096</Real>
</Sequence>
<Sequence Name="Row5">
<Int Name="Length">2</Int>
- <Real>5</Real>
+ <Real>5.000</Real>
<Real>0.12</Real>
</Sequence>
<Sequence Name="Row6">
<Int Name="Length">2</Int>
- <Real>6</Real>
+ <Real>6.000</Real>
<Real>0.144</Real>
</Sequence>
<Sequence Name="Row7">
<Int Name="Length">2</Int>
- <Real>7</Real>
+ <Real>7.000</Real>
<Real>0.168</Real>
</Sequence>
<Sequence Name="Row8">
<Int Name="Length">2</Int>
- <Real>8</Real>
+ <Real>8.000</Real>
<Real>0.192</Real>
</Sequence>
<Sequence Name="Row9">
<Int Name="Length">2</Int>
- <Real>9</Real>
+ <Real>9.000</Real>
<Real>0.216</Real>
</Sequence>
</XvgData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <String Name="CommandLine">-trestart 10 -sel 2</String>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[ Water] = 5.3526 (+/- 0.0465) (1e-5 cm^2/s)"</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">2</Int>
+ <Real>0.000</Real>
+ <Real>0</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">2</Int>
+ <Real>2.000</Real>
+ <Real>0.0713805</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">2</Int>
+ <Real>4.000</Real>
+ <Real>0.136132</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">2</Int>
+ <Real>6.000</Real>
+ <Real>0.202192</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">2</Int>
+ <Real>8.000</Real>
+ <Real>0.269467</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">2</Int>
+ <Real>10.000</Real>
+ <Real>0.333361</Real>
+ </Sequence>
+ <Sequence Name="Row6">
+ <Int Name="Length">2</Int>
+ <Real>12.000</Real>
+ <Real>0.389284</Real>
+ </Sequence>
+ <Sequence Name="Row7">
+ <Int Name="Length">2</Int>
+ <Real>14.000</Real>
+ <Real>0.449046</Real>
+ </Sequence>
+ <Sequence Name="Row8">
+ <Int Name="Length">2</Int>
+ <Real>16.000</Real>
+ <Real>0.514649</Real>
+ </Sequence>
+ <Sequence Name="Row9">
+ <Int Name="Length">2</Int>
+ <Real>18.000</Real>
+ <Real>0.59199</Real>
+ </Sequence>
+ <Sequence Name="Row10">
+ <Int Name="Length">2</Int>
+ <Real>20.000</Real>
+ <Real>0.665065</Real>
+ </Sequence>
+ <Sequence Name="Row11">
+ <Int Name="Length">2</Int>
+ <Real>22.000</Real>
+ <Real>0.699726</Real>
+ </Sequence>
+ <Sequence Name="Row12">
+ <Int Name="Length">2</Int>
+ <Real>24.000</Real>
+ <Real>0.74971</Real>
+ </Sequence>
+ <Sequence Name="Row13">
+ <Int Name="Length">2</Int>
+ <Real>26.000</Real>
+ <Real>0.830236</Real>
+ </Sequence>
+ <Sequence Name="Row14">
+ <Int Name="Length">2</Int>
+ <Real>28.000</Real>
+ <Real>0.912884</Real>
+ </Sequence>
+ <Sequence Name="Row15">
+ <Int Name="Length">2</Int>
+ <Real>30.000</Real>
+ <Real>0.984772</Real>
+ </Sequence>
+ <Sequence Name="Row16">
+ <Int Name="Length">2</Int>
+ <Real>32.000</Real>
+ <Real>1.03085</Real>
+ </Sequence>
+ <Sequence Name="Row17">
+ <Int Name="Length">2</Int>
+ <Real>34.000</Real>
+ <Real>1.08669</Real>
+ </Sequence>
+ <Sequence Name="Row18">
+ <Int Name="Length">2</Int>
+ <Real>36.000</Real>
+ <Real>1.18283</Real>
+ </Sequence>
+ <Sequence Name="Row19">
+ <Int Name="Length">2</Int>
+ <Real>38.000</Real>
+ <Real>1.26495</Real>
+ </Sequence>
+ <Sequence Name="Row20">
+ <Int Name="Length">2</Int>
+ <Real>40.000</Real>
+ <Real>1.34173</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ </OutputFiles>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <String Name="CommandLine">-trestart 1 -sel 2</String>
+ <OutputFiles Name="Files">
+ <File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
+ <XvgLegend Name="Legend">
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[ Water] = 5.4502 (+/- 0.2739) (1e-5 cm^2/s)"</String>
+ </XvgLegend>
+ <XvgData Name="Data">
+ <Sequence Name="Row0">
+ <Int Name="Length">2</Int>
+ <Real>0.000</Real>
+ <Real>0</Real>
+ </Sequence>
+ <Sequence Name="Row1">
+ <Int Name="Length">2</Int>
+ <Real>2.000</Real>
+ <Real>0.0722943</Real>
+ </Sequence>
+ <Sequence Name="Row2">
+ <Int Name="Length">2</Int>
+ <Real>4.000</Real>
+ <Real>0.13585</Real>
+ </Sequence>
+ <Sequence Name="Row3">
+ <Int Name="Length">2</Int>
+ <Real>6.000</Real>
+ <Real>0.19875</Real>
+ </Sequence>
+ <Sequence Name="Row4">
+ <Int Name="Length">2</Int>
+ <Real>8.000</Real>
+ <Real>0.262333</Real>
+ </Sequence>
+ <Sequence Name="Row5">
+ <Int Name="Length">2</Int>
+ <Real>10.000</Real>
+ <Real>0.324822</Real>
+ </Sequence>
+ <Sequence Name="Row6">
+ <Int Name="Length">2</Int>
+ <Real>12.000</Real>
+ <Real>0.388126</Real>
+ </Sequence>
+ <Sequence Name="Row7">
+ <Int Name="Length">2</Int>
+ <Real>14.000</Real>
+ <Real>0.45144</Real>
+ </Sequence>
+ <Sequence Name="Row8">
+ <Int Name="Length">2</Int>
+ <Real>16.000</Real>
+ <Real>0.515772</Real>
+ </Sequence>
+ <Sequence Name="Row9">
+ <Int Name="Length">2</Int>
+ <Real>18.000</Real>
+ <Real>0.581858</Real>
+ </Sequence>
+ <Sequence Name="Row10">
+ <Int Name="Length">2</Int>
+ <Real>20.000</Real>
+ <Real>0.645979</Real>
+ </Sequence>
+ <Sequence Name="Row11">
+ <Int Name="Length">2</Int>
+ <Real>22.000</Real>
+ <Real>0.708053</Real>
+ </Sequence>
+ <Sequence Name="Row12">
+ <Int Name="Length">2</Int>
+ <Real>24.000</Real>
+ <Real>0.772597</Real>
+ </Sequence>
+ <Sequence Name="Row13">
+ <Int Name="Length">2</Int>
+ <Real>26.000</Real>
+ <Real>0.84093</Real>
+ </Sequence>
+ <Sequence Name="Row14">
+ <Int Name="Length">2</Int>
+ <Real>28.000</Real>
+ <Real>0.908325</Real>
+ </Sequence>
+ <Sequence Name="Row15">
+ <Int Name="Length">2</Int>
+ <Real>30.000</Real>
+ <Real>0.975539</Real>
+ </Sequence>
+ <Sequence Name="Row16">
+ <Int Name="Length">2</Int>
+ <Real>32.000</Real>
+ <Real>1.04299</Real>
+ </Sequence>
+ <Sequence Name="Row17">
+ <Int Name="Length">2</Int>
+ <Real>34.000</Real>
+ <Real>1.11211</Real>
+ </Sequence>
+ <Sequence Name="Row18">
+ <Int Name="Length">2</Int>
+ <Real>36.000</Real>
+ <Real>1.19384</Real>
+ </Sequence>
+ <Sequence Name="Row19">
+ <Int Name="Length">2</Int>
+ <Real>38.000</Real>
+ <Real>1.27086</Real>
+ </Sequence>
+ <Sequence Name="Row20">
+ <Int Name="Length">2</Int>
+ <Real>40.000</Real>
+ <Real>1.34173</Real>
+ </Sequence>
+ </XvgData>
+ </File>
+ </OutputFiles>
+</ReferenceData>
<?xml version="1.0"?>
<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
<ReferenceData>
+ <String Name="CommandLine">-trestart 200 -lateral z -sel 0</String>
<OutputFiles Name="Files">
<File Name="-o">
+ <XvgLegend Name="DiffusionCoefficient"></XvgLegend>
<XvgLegend Name="Legend">
- <String Name="XvgLegend"><![CDATA[
-title "Mean Square Displacement"
-xaxis label "Time (ps)"
-yaxis label "MSD (nm\S2\N)"
-TYPE xy
-]]></String>
+ <String>title "Mean Squared Displacement"</String>
+ <String>xaxis label "tau (ps)"</String>
+ <String>yaxis label "MSD (nm\\S2\\N)"</String>
+ <String>TYPE xy</String>
+ <String>s0 legend "D[ particles] = 6.0000 (+/- 0.0000) (1e-5 cm^2/s)"</String>
</XvgLegend>
<XvgData Name="Data">
<Sequence Name="Row0">
<Int Name="Length">2</Int>
- <Real>0</Real>
+ <Real>0.000</Real>
<Real>0</Real>
</Sequence>
<Sequence Name="Row1">
<Int Name="Length">2</Int>
- <Real>1</Real>
+ <Real>1.000</Real>
<Real>0.024</Real>
</Sequence>
<Sequence Name="Row2">
<Int Name="Length">2</Int>
- <Real>2</Real>
+ <Real>2.000</Real>
<Real>0.048</Real>
</Sequence>
<Sequence Name="Row3">
<Int Name="Length">2</Int>
- <Real>3</Real>
+ <Real>3.000</Real>
<Real>0.072</Real>
</Sequence>
<Sequence Name="Row4">
<Int Name="Length">2</Int>
- <Real>4</Real>
+ <Real>4.000</Real>
<Real>0.096</Real>
</Sequence>
<Sequence Name="Row5">
<Int Name="Length">2</Int>
- <Real>5</Real>
+ <Real>5.000</Real>
<Real>0.12</Real>
</Sequence>
<Sequence Name="Row6">
<Int Name="Length">2</Int>
- <Real>6</Real>
+ <Real>6.000</Real>
<Real>0.144</Real>
</Sequence>
<Sequence Name="Row7">
<Int Name="Length">2</Int>
- <Real>7</Real>
+ <Real>7.000</Real>
<Real>0.168</Real>
</Sequence>
<Sequence Name="Row8">
<Int Name="Length">2</Int>
- <Real>8</Real>
+ <Real>8.000</Real>
<Real>0.192</Real>
</Sequence>
<Sequence Name="Row9">
<Int Name="Length">2</Int>
- <Real>9</Real>
+ <Real>9.000</Real>
<Real>0.216</Real>
</Sequence>
</XvgData>
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
- * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
registerModule(manager, &gmx_mdmat, "mdmat", "Calculate residue contact maps");
registerModule(
manager, &gmx_mindist, "mindist", "Calculate the minimum distance between two groups");
- registerModule(manager, &gmx_msd, "msd", "Calculates mean square displacements");
registerModule(manager, &gmx_nmeig, "nmeig", "Diagonalize the Hessian for normal mode analysis");
registerModule(manager,
&gmx_nmens,
gmx::CommandLineModuleGroup group =
manager->addModuleGroup("Mass distribution properties over time");
group.addModule("gyrate");
- group.addModule("msd");
group.addModule("polystat");
group.addModule("rdf");
group.addModule("rotacf");
[ System ]
- 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
- 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
- 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
- 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
- 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
- 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
- 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
- 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
- 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
- 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
- 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
- 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
- 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
- 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
- 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
- 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
- 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
- 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
- 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
- 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
- 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315
- 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
- 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
- 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
- 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375
- 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390
- 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
- 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
- 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435
- 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450
- 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465
- 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
- 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495
- 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
- 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525
- 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540
- 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555
- 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
- 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585
- 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600
- 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615
- 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630
- 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645
- 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
- 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675
- 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690
- 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705
- 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720
- 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735
- 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750
- 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765
- 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780
- 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795
- 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810
- 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825
- 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840
- 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855
- 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870
- 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885
- 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900
- 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915
- 916 917 918 919 920 921 922 923
+ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
+ 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
+ 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
+ 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
+ 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
+ 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
+ 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
+ 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
+ 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
+ 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
+ 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
+ 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
+ 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
+ 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
+ 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
+ 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
+ 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
+ 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
+ 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
+ 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
+ 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315
+ 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
+ 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
+ 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
+ 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375
+ 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390
+ 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
+ 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
+ 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435
+ 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450
+ 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465
+ 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
+ 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495
+ 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
+ 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525
+ 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540
+ 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555
+ 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570
+ 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585
+ 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600
+ 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615
+ 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630
+ 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645
+ 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
+ 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675
+ 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690
+ 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705
+ 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720
+ 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735
+ 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750
+ 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765
+ 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780
+ 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795
+ 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810
+ 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825
+ 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840
+ 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855
+ 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870
+ 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885
+ 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900
+ 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915
+ 916 917 918 919 920 921 922 923
+[ Protein ]
+ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
+ 16 17 18 19 20 21 22 23 24 25 26 27 28 29
+[ Water ]
+ 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
+ 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
+ 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
+ 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
+ 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
+ 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
+ 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
+ 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
+ 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
+ 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
+ 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194
+ 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
+ 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
+ 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
+ 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
+ 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
+ 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284
+ 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
+ 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314
+ 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329
+ 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344
+ 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
+ 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374
+ 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389
+ 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404
+ 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419
+ 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434
+ 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449
+ 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464
+ 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479
+ 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494
+ 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509
+ 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
+ 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539
+ 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
+ 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569
+ 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584
+ 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599
+ 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614
+ 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629
+ 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644
+ 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659
+ 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674
+ 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689
+ 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704
+ 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719
+ 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734
+ 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749
+ 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764
+ 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779
+ 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794
+ 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809
+ 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824
+ 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839
+ 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854
+ 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869
+ 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884
+ 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899
+ 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914
+ 915 916 917 918 919 920 921 922 923
+[ some_water_subset ]
+ 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
+ 45 46 47
--- /dev/null
+[ particles ]
+1 2 3
--- /dev/null
+msd_beads
+ 3
+ 1A A 1 2.167 1.833 1.500
+ 2A A 2 0.000 0.000 0.000
+ 3A A 3 3.167 3.833 4.500
+ 5.00000 5.00000 5.00000
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff