2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2021, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Defines the trajectory analysis module for mean squared displacement calculations.
39 * \author Kevin Boyd <kevin44boyd@gmail.com>
40 * \ingroup module_trajectoryanalysis
50 #include "gromacs/analysisdata/analysisdata.h"
51 #include "gromacs/analysisdata/modules/average.h"
52 #include "gromacs/analysisdata/modules/plot.h"
53 #include "gromacs/analysisdata/paralleloptions.h"
54 #include "gromacs/fileio/oenv.h"
55 #include "gromacs/fileio/trxio.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/vectypes.h"
58 #include "gromacs/options/basicoptions.h"
59 #include "gromacs/options/filenameoption.h"
60 #include "gromacs/options/ioptionscontainer.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/selection/selectionoption.h"
63 #include "gromacs/statistics/statistics.h"
64 #include "gromacs/trajectoryanalysis/analysissettings.h"
65 #include "gromacs/trajectoryanalysis/topologyinformation.h"
66 #include "gromacs/trajectory/trajectoryframe.h"
67 #include "gromacs/utility/stringutil.h"
68 #include "gromacs/utility.h"
72 namespace analysismodules
78 //! Convert nm^2/ps to 10e-5 cm^2/s
79 constexpr double c_diffusionConversionFactor = 1000.0;
80 //! Three dimensional diffusion coefficient multiplication constant.
81 constexpr double c_3DdiffusionDimensionFactor = 6.0;
82 //! Two dimensional diffusion coefficient multiplication constant.
83 constexpr double c_2DdiffusionDimensionFactor = 4.0;
84 //! One dimensional diffusion coefficient multiplication constant.
85 constexpr double c_1DdiffusionDimensionFactor = 2.0;
88 /*! \brief Mean Squared Displacement data accumulator
90 * This class is used to accumulate individual MSD data points
91 * and emit tau-averaged results once data is finished collecting. Displacements at each observed
92 * time difference (tau) are recorded from the trajectory. Because it is not known in advance which
93 * time differences will be observed from the trajectory, this data structure is built adaptively.
94 * New columns corresponding to observed time differences are added as needed, and additional
95 * observations at formerly observed time differences are added to those columns. Separate time lags
96 * will likely have differing total data points.
98 * Data columns per tau are accessed via operator[], which always guarantees
99 * a column is initialized and returns an MsdColumProxy to the column that can push data.
104 //! Proxy to a MsdData tau column vector. Supports only push_back.
108 MsdColumnProxy(std::vector<double>* column) : column_(column) {}
110 void push_back(double value) { column_->push_back(value); }
113 std::vector<double>* column_;
115 //! Returns a proxy to the column for the given tau index. Guarantees that the column is initialized.
116 MsdColumnProxy operator[](size_t index)
118 if (msds_.size() <= index)
120 msds_.resize(index + 1);
122 return MsdColumnProxy(&msds_[index]);
124 /*! \brief Compute per-tau MSDs averaged over all added points.
126 * The resulting vector is size(max tau index). Any indices
127 * that have no data points have MSD set to 0.
129 * \return Average MSD per tau
131 [[nodiscard]] std::vector<real> averageMsds() const;
134 //! Results - first indexed by tau, then data points
135 std::vector<std::vector<double>> msds_;
139 std::vector<real> MsdData::averageMsds() const
141 std::vector<real> msdSums;
142 msdSums.reserve(msds_.size());
143 for (gmx::ArrayRef<const double> msdValues : msds_)
145 if (msdValues.empty())
147 msdSums.push_back(0.0);
150 msdSums.push_back(std::accumulate(msdValues.begin(), msdValues.end(), 0.0, std::plus<>())
156 /*! \brief Calculates 1,2, or 3D distance for two vectors.
158 * \todo Remove NOLINTs once clang-tidy is updated to v11, it should be able to handle constexpr.
160 * \tparam x If true, calculate x dimension of displacement
161 * \tparam y If true, calculate y dimension of displacement
162 * \tparam z If true, calculate z dimension of displacement
163 * \param[in] c1 First point
164 * \param[in] c2 Second point
165 * \return Euclidian distance for the given dimension.
167 template<bool x, bool y, bool z>
168 inline double calcSingleSquaredDistance(const RVec c1, const RVec c2)
170 static_assert(x || y || z, "zero-dimensional MSD selected");
171 const DVec firstCoords = c1.toDVec();
172 const DVec secondCoords = c2.toDVec();
176 result += (firstCoords[XX] - secondCoords[XX]) * (firstCoords[XX] - secondCoords[XX]);
178 // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
181 result += (firstCoords[YY] - secondCoords[YY]) * (firstCoords[YY] - secondCoords[YY]);
183 // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
186 result += (firstCoords[ZZ] - secondCoords[ZZ]) * (firstCoords[ZZ] - secondCoords[ZZ]);
188 // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
192 /*! \brief Calculate average displacement between sets of points
194 * Each displacement c1[i] - c2[i] is calculated and the distances
197 * \tparam x If true, calculate x dimension of displacement
198 * \tparam y If true, calculate y dimension of displacement
199 * \tparam z If true, calculate z dimension of displacement
200 * \param[in] c1 First vector
201 * \param[in] c2 Second vector
202 * \return Per-particle averaged distance
204 template<bool x, bool y, bool z>
205 double calcAverageDisplacement(ArrayRef<const RVec> c1, ArrayRef<const RVec> c2)
208 for (size_t i = 0; i < c1.size(); i++)
210 result += calcSingleSquaredDistance<x, y, z>(c1[i], c2[i]);
212 return result / c1.size();
216 //! Describes 1D MSDs, in the given dimension.
217 enum class SingleDimDiffType : int
226 //! Describes 2D MSDs, in the plane normal to the given dimension.
227 enum class TwoDimDiffType : int
236 /*! \brief Removes jumps across periodic boundaries for currentFrame, based on the positions in
237 * previousFrame. Updates currentCoords in place.
239 void removePbcJumps(ArrayRef<RVec> currentCoords, ArrayRef<const RVec> previousCoords, t_pbc* pbc)
241 // There are two types of "pbc removal" in gmx msd. The first happens in the trajectoryanalysis
242 // framework, which makes molecules whole across periodic boundaries and is done
243 // automatically where the inputs support it. This lambda performs the second PBC correction, where
244 // any "jump" across periodic boundaries BETWEEN FRAMES is put back. The order of these
245 // operations is important - since the first transformation may only apply to part of a
246 // molecule (e.g., one half in/out of the box is put on one side of the box), the
247 // subsequent step needs to be applied to the molecule COM rather than individual atoms, or
248 // we'd have a clash where the per-mol PBC removal moves an atom that gets put back into
249 // it's original position by the second transformation. Therefore, this second transformation
250 // is applied *after* per molecule coordinates have been consolidated into COMs.
251 auto pbcRemover = [pbc](RVec in, RVec prev) {
253 pbc_dx(pbc, in, prev, dx);
257 currentCoords.begin(), currentCoords.end(), previousCoords.begin(), currentCoords.begin(), pbcRemover);
260 //! Holds data needed for MSD calculations for a single molecule, if requested.
263 //! Number of atoms in the molecule.
267 //! MSD accumulator and calculator for the molecule
269 //! Calculated diffusion coefficient
270 real diffusionCoefficient = 0;
273 /*! \brief Handles coordinate operations for MSD calculations.
275 * Can be used to hold coordinates for individual atoms as well as molecules COMs. Handles PBC
276 * jump removal between consecutive frames.
278 * Only previous_ contains valid data between calls to buildCoordinates(), although both vectors
279 * are maintained at the correct size for the number of coordinates needed.
281 class MsdCoordinateManager
284 MsdCoordinateManager(const int numAtoms,
285 ArrayRef<const MoleculeData> molecules,
286 ArrayRef<const int> moleculeIndexMapping) :
287 current_(numAtoms), previous_(numAtoms), molecules_(molecules), moleculeIndexMapping_(moleculeIndexMapping)
290 /*! \brief Prepares coordinates for the current frame.
292 * Reads in selection data, and returns an ArrayRef of the particle positions or molecule
293 * centers of mass (if the molecules input is not empty). Removes jumps across periodic
294 * boundaries based on the previous frame coordinates, except for the first frame built with
295 * builtCoordinates(), which has no previous frame as a reference.
297 * \param[in] sel The selection object which holds coordinates
298 * \param[in] pbc Information about periodicity.
299 * \returns The current frames coordinates in proper format.
301 ArrayRef<const RVec> buildCoordinates(const Selection& sel, t_pbc* pbc);
304 //! The current coordinates.
305 std::vector<RVec> current_;
306 //! The previous frame's coordinates.
307 std::vector<RVec> previous_;
309 ArrayRef<const MoleculeData> molecules_;
310 //! Mapping of atom indices to molecle indices;
311 ArrayRef<const int> moleculeIndexMapping_;
312 //! Tracks if a previous frame exists to compare with for PBC handling.
313 bool containsPreviousFrame_ = false;
316 ArrayRef<const RVec> MsdCoordinateManager::buildCoordinates(const Selection& sel, t_pbc* pbc)
319 if (molecules_.empty())
321 std::copy(sel.coordinates().begin(), sel.coordinates().end(), current_.begin());
325 // Prepare for molecule COM calculation, then sum up all positions per molecule.
327 std::fill(current_.begin(), current_.end(), RVec(0, 0, 0));
328 gmx::ArrayRef<const real> masses = sel.masses();
329 for (int i = 0; i < sel.posCount(); i++)
331 const int moleculeIndex = moleculeIndexMapping_[i];
332 // accumulate ri * mi, and do division at the end to minimize number of divisions.
333 current_[moleculeIndex] += RVec(sel.position(i).x()) * masses[i];
335 // Divide accumulated mass * positions to get COM.
336 std::transform(current_.begin(),
340 [](const RVec& position, const MoleculeData& molecule) -> RVec {
341 return position / molecule.mass;
345 if (containsPreviousFrame_)
347 removePbcJumps(current_, previous_, pbc);
351 containsPreviousFrame_ = true;
354 // Previous is no longer needed, swap with current and return "current" coordinates which
355 // now reside in previous.
356 current_.swap(previous_);
361 //! Holds per-group coordinates, analysis, and results.
364 explicit MsdGroupData(const Selection& inputSel,
365 ArrayRef<const MoleculeData> molecules,
366 ArrayRef<const int> moleculeAtomMapping) :
368 coordinateManager_(molecules.empty() ? sel.posCount() : molecules.size(), molecules, moleculeAtomMapping)
372 //! Selection associated with this group.
373 const Selection& sel;
375 //! Stored coordinates, indexed by frame then atom number.
376 std::vector<std::vector<RVec>> frames;
378 //! MSD result accumulator
380 //! Coordinate handler for the group.
381 MsdCoordinateManager coordinateManager_;
382 //! Collector for processed MSD averages per tau
383 std::vector<real> msdSums;
384 //! Fitted diffusion coefficient
385 real diffusionCoefficient = 0.0;
386 //! Uncertainty of diffusion coefficient
392 /*! \brief Implements the gmx msd module
394 * \todo Implement -tensor for full MSD tensor calculation
395 * \todo Implement -rmcomm for total-frame COM removal
396 * \todo Implement -pdb for molecule B factors
398 class Msd : public TrajectoryAnalysisModule
403 void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
404 void optionsFinished(TrajectoryAnalysisSettings* settings) override;
405 void initAfterFirstFrame(const TrajectoryAnalysisSettings& settings, const t_trxframe& fr) override;
406 void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
407 void analyzeFrame(int frameNumber,
408 const t_trxframe& frame,
410 TrajectoryAnalysisModuleData* pdata) override;
411 void finishAnalysis(int nframes) override;
412 void writeOutput() override;
415 //! Selections for MSD output
416 SelectionList selections_;
418 //! MSD type information, for -type {x,y,z}
419 SingleDimDiffType singleDimType_ = SingleDimDiffType::Unused;
420 //! MSD type information, for -lateral {x,y,z}
421 TwoDimDiffType twoDimType_ = TwoDimDiffType::Unused;
422 //! Diffusion coefficient conversion factor
423 double diffusionCoefficientDimensionFactor_ = c_3DdiffusionDimensionFactor;
424 //! Method used to calculate MSD - changes based on dimensonality.
425 std::function<double(ArrayRef<const RVec>, ArrayRef<const RVec>)> calcMsd_ =
426 calcAverageDisplacement<true, true, true>;
428 //! Picoseconds between restarts
429 double trestart_ = 10.0;
432 //! Inter-frame delta-t
433 std::optional<double> dt_ = std::nullopt;
434 //! Inter-frame maximum time delta.
435 double maxTau_ = std::numeric_limits<double>::max();
436 //! Tracks the index of the first valid frame. Starts at 0, increases as old frames are deleted
437 size_t firstValidFrame_ = 0;
439 //! First tau value to fit from for diffusion coefficient, defaults to 0.1 * max tau
440 real beginFit_ = -1.0;
441 //! Final tau value to fit to for diffusion coefficient, defaults to 0.9 * max tau
444 //! All selection group-specific data stored here.
445 std::vector<MsdGroupData> groupData_;
447 //! Time of all frames.
448 std::vector<double> times_;
449 //! Taus for output - won't know the size until the end.
450 std::vector<double> taus_;
451 //! Tau indices for fitting.
452 size_t beginFitIndex_ = 0;
453 size_t endFitIndex_ = 0;
455 // MSD per-molecule stuff
456 //! Are we doing molecule COM-based MSDs?
457 bool molSelected_ = false;
458 //! Per molecule topology information and MSD accumulators.
459 std::vector<MoleculeData> molecules_;
460 //! Atom index -> mol index map
461 std::vector<int> moleculeIndexMappings_;
464 AnalysisData msdPlotData_;
465 AnalysisData msdMoleculePlotData_;
467 AnalysisDataPlotSettings plotSettings_;
468 //! Per-tau MSDs for each selected group
470 //! Per molecule diffusion coefficients if -mol is selected.
471 std::string moleculeOutput_;
475 Msd::Msd() = default;
477 void Msd::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
479 static const char* const desc[] = {
480 "[THISMODULE] computes the mean square displacement (MSD) of atoms from",
481 "a set of initial positions. This provides an easy way to compute",
482 "the diffusion constant using the Einstein relation.",
483 "The time between the reference points for the MSD calculation",
484 "is set with [TT]-trestart[tt].",
485 "The diffusion constant is calculated by least squares fitting a",
486 "straight line (D*t + c) through the MSD(t) from [TT]-beginfit[tt] to",
487 "[TT]-endfit[tt] (note that t is time from the reference positions,",
488 "not simulation time). An error estimate given, which is the difference",
489 "of the diffusion coefficients obtained from fits over the two halves",
490 "of the fit interval.[PAR]",
491 "There are three, mutually exclusive, options to determine different",
492 "types of mean square displacement: [TT]-type[tt], [TT]-lateral[tt]",
493 "and [TT]-ten[tt]. Option [TT]-ten[tt] writes the full MSD tensor for",
494 "each group, the order in the output is: trace xx yy zz yx zx zy.[PAR]",
495 "If [TT]-mol[tt] is set, [THISMODULE] plots the MSD for individual molecules",
496 "(including making molecules whole across periodic boundaries): ",
497 "for each individual molecule a diffusion constant is computed for ",
498 "its center of mass. The chosen index group will be split into ",
499 "molecules. With -mol, only one index group can be selected.[PAR]",
500 "The diffusion coefficient is determined by linear regression of the MSD.",
501 "When [TT]-beginfit[tt] is -1, fitting starts at 10%",
502 "and when [TT]-endfit[tt] is -1, fitting goes to 90%.",
503 "Using this option one also gets an accurate error estimate",
504 "based on the statistics between individual molecules.",
505 "Note that this diffusion coefficient and error estimate are only",
506 "accurate when the MSD is completely linear between",
507 "[TT]-beginfit[tt] and [TT]-endfit[tt].[PAR]",
508 "By default, [THISMODULE] compares all trajectory frames against every frame stored at",
509 "[TT]-trestart[TT] intervals, so the number of frames stored scales linearly with the",
510 "number of frames processed. This can lead to long analysis times and out-of-memory errors",
511 "for long/large trajectories, and often the data at higher time deltas lacks sufficient ",
512 "sampling, often manifesting as a wobbly line on the MSD plot after a straighter region at",
513 "lower time deltas. The [TT]-maxtau[TT] option can be used to cap the maximum time delta",
514 "for frame comparison, which may improve performance and can be used to avoid",
515 "out-of-memory issues.[PAR]"
517 settings->setHelpText(desc);
520 options->addOption(SelectionOption("sel")
521 .storeVector(&selections_)
525 .description("Selections to compute MSDs for from the reference"));
527 // Select MSD type - defaults to 3D if neither option is selected.
528 EnumerationArray<SingleDimDiffType, const char*> enumTypeNames = { "x", "y", "z", "unused" };
529 EnumerationArray<TwoDimDiffType, const char*> enumLateralNames = { "x", "y", "z", "unused" };
530 options->addOption(EnumOption<SingleDimDiffType>("type")
531 .enumValue(enumTypeNames)
532 .store(&singleDimType_)
533 .defaultValue(SingleDimDiffType::Unused));
534 options->addOption(EnumOption<TwoDimDiffType>("lateral")
535 .enumValue(enumLateralNames)
537 .defaultValue(TwoDimDiffType::Unused));
539 options->addOption(DoubleOption("trestart")
540 .description("Time between restarting points in trajectory (ps)")
544 DoubleOption("maxtau")
545 .description("Maximum time delta between frames to calculate MSDs for (ps)")
548 RealOption("beginfit").description("Time point at which to start fitting.").store(&beginFit_));
549 options->addOption(RealOption("endfit").description("End time for fitting.").store(&endFit_));
552 options->addOption(FileNameOption("o")
553 .filetype(OptionFileType::Plot)
556 .defaultBasename("msdout")
557 .description("MSD output"));
559 FileNameOption("mol")
560 .filetype(OptionFileType::Plot)
562 .store(&moleculeOutput_)
563 .storeIsSet(&molSelected_)
564 .defaultBasename("diff_mol")
565 .description("Report diffusion coefficients for each molecule in selection"));
568 void Msd::optionsFinished(TrajectoryAnalysisSettings gmx_unused* settings)
570 if (singleDimType_ != SingleDimDiffType::Unused && twoDimType_ != TwoDimDiffType::Unused)
572 std::string errorMessage =
573 "Options -type and -lateral are mutually exclusive. Choose one or neither (for 3D "
575 GMX_THROW(InconsistentInputError(errorMessage.c_str()));
577 if (selections_.size() > 1 && molSelected_)
579 std::string errorMessage =
580 "Cannot have multiple groups selected with -sel when using -mol.";
581 GMX_THROW(InconsistentInputError(errorMessage.c_str()));
586 void Msd::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
588 plotSettings_ = settings.plotSettings();
590 // Enumeration helpers for dispatching the right MSD calculation type.
591 const EnumerationArray<SingleDimDiffType, decltype(&calcAverageDisplacement<true, true, true>)>
592 oneDimensionalMsdFunctions = { &calcAverageDisplacement<true, false, false>,
593 &calcAverageDisplacement<false, true, false>,
594 &calcAverageDisplacement<false, false, true> };
595 const EnumerationArray<TwoDimDiffType, decltype(&calcAverageDisplacement<true, true, true>)>
596 twoDimensionalMsdFunctions = { calcAverageDisplacement<false, true, true>,
597 calcAverageDisplacement<true, false, true>,
598 calcAverageDisplacement<true, true, false> };
600 // Parse dimensionality and assign the MSD calculating function.
601 // Note if we don't hit either of these cases, we're computing 3D MSDs.
602 if (singleDimType_ != SingleDimDiffType::Unused)
604 calcMsd_ = oneDimensionalMsdFunctions[singleDimType_];
605 diffusionCoefficientDimensionFactor_ = c_1DdiffusionDimensionFactor;
607 else if (twoDimType_ != TwoDimDiffType::Unused)
609 calcMsd_ = twoDimensionalMsdFunctions[twoDimType_];
610 diffusionCoefficientDimensionFactor_ = c_2DdiffusionDimensionFactor;
613 // TODO validate that we have mol info and not atom only - and masses, and topology.
616 Selection& sel = selections_[0];
617 const int nMol = sel.initOriginalIdsToGroup(top.mtop(), INDEX_MOL);
619 gmx::ArrayRef<const int> mappedIds = selections_[0].mappedIds();
620 moleculeIndexMappings_.resize(selections_[0].posCount());
621 std::copy(mappedIds.begin(), mappedIds.end(), moleculeIndexMappings_.begin());
623 // Precalculate each molecules mass for speeding up COM calculations.
624 ArrayRef<const real> masses = sel.masses();
626 molecules_.resize(nMol);
627 for (int i = 0; i < sel.posCount(); i++)
629 molecules_[mappedIds[i]].atomCount++;
630 molecules_[mappedIds[i]].mass += masses[i];
634 // Accumulated frames and results
635 for (const Selection& sel : selections_)
637 groupData_.emplace_back(sel, molecules_, moleculeIndexMappings_);
641 void Msd::initAfterFirstFrame(const TrajectoryAnalysisSettings gmx_unused& settings, const t_trxframe& fr)
643 t0_ = std::round(fr.time);
647 void Msd::analyzeFrame(int gmx_unused frameNumber,
648 const t_trxframe& frame,
650 TrajectoryAnalysisModuleData* pdata)
652 const real time = std::round(frame.time);
653 // Need to populate dt on frame 2;
654 if (!dt_.has_value() && !times_.empty())
656 dt_ = time - times_[0];
659 // Each frame gets an entry in times, but frameTimes only updates if we're at a restart.
660 times_.push_back(time);
662 // Each frame will get a tau between it and frame 0, and all other frame combos should be
664 if (const double tau = time - times_[0]; tau <= maxTau_)
666 taus_.push_back(time - times_[0]);
669 for (MsdGroupData& msdData : groupData_)
671 const Selection& sel = pdata->parallelSelection(msdData.sel);
673 ArrayRef<const RVec> coords = msdData.coordinateManager_.buildCoordinates(sel, pbc);
675 // For each preceding frame, calculate tau and do comparison.
676 for (size_t i = firstValidFrame_; i < msdData.frames.size(); i++)
678 // If dt > trestart, the tau increment will be dt rather than trestart.
679 const double tau = time - (t0_ + std::max<double>(trestart_, *dt_) * i);
682 // The (now empty) entry is no longer needed, so over time the outer vector will
683 // grow with extraneous empty elements persisting, but the alternative would require
684 // some more complicated remapping of tau to frame index.
685 msdData.frames[i].clear();
686 firstValidFrame_ = i + 1;
689 int64_t tauIndex = gmx::roundToInt64(tau / *dt_);
690 msdData.msds[tauIndex].push_back(calcMsd_(coords, msdData.frames[i]));
692 for (size_t molInd = 0; molInd < molecules_.size(); molInd++)
694 molecules_[molInd].msdData[tauIndex].push_back(
695 calcMsd_(arrayRefFromArray(&coords[molInd], 1),
696 arrayRefFromArray(&msdData.frames[i][molInd], 1)));
701 // We only store the frame for the future if it's a restart per -trestart.
702 if (bRmod(time, t0_, trestart_))
704 msdData.frames.emplace_back(coords.begin(), coords.end());
709 //! Calculate the tau index for fitting. If userFitTau < 0, uses the default fraction of max tau.
710 static size_t calculateFitIndex(const int userFitTau,
711 const double defaultTauFraction,
717 return gmx::roundToInt((numTaus - 1) * defaultTauFraction);
719 return std::min<size_t>(numTaus - 1, gmx::roundToInt(static_cast<double>(userFitTau) / dt));
723 void Msd::finishAnalysis(int gmx_unused nframes)
725 static constexpr double c_defaultStartFitIndexFraction = 0.1;
726 static constexpr double c_defaultEndFitIndexFraction = 0.9;
727 beginFitIndex_ = calculateFitIndex(beginFit_, c_defaultStartFitIndexFraction, taus_.size(), *dt_);
728 endFitIndex_ = calculateFitIndex(endFit_, c_defaultEndFitIndexFraction, taus_.size(), *dt_);
729 const int numTausForFit = 1 + endFitIndex_ - beginFitIndex_;
731 // These aren't used, except for correlationCoefficient, which is used to estimate error if
732 // enough points are available.
733 real b = 0.0, correlationCoefficient = 0.0, chiSquared = 0.0;
735 for (MsdGroupData& msdData : groupData_)
737 msdData.msdSums = msdData.msds.averageMsds();
739 if (numTausForFit >= 4)
741 const int halfNumTaus = numTausForFit / 2;
742 const int secondaryStartIndex = beginFitIndex_ + halfNumTaus;
743 // Split the fit in 2, and compare the results of each fit;
744 real a = 0.0, a2 = 0.0;
745 lsq_y_ax_b_xdouble(halfNumTaus,
746 &taus_[beginFitIndex_],
747 &msdData.msdSums[beginFitIndex_],
750 &correlationCoefficient,
752 lsq_y_ax_b_xdouble(halfNumTaus,
753 &taus_[secondaryStartIndex],
754 &msdData.msdSums[secondaryStartIndex],
757 &correlationCoefficient,
759 msdData.sigma = std::abs(a - a2);
761 lsq_y_ax_b_xdouble(numTausForFit,
762 &taus_[beginFitIndex_],
763 &msdData.msdSums[beginFitIndex_],
764 &msdData.diffusionCoefficient,
766 &correlationCoefficient,
768 msdData.diffusionCoefficient *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
769 msdData.sigma *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
772 for (MoleculeData& molecule : molecules_)
774 std::vector<real> msds = molecule.msdData.averageMsds();
775 lsq_y_ax_b_xdouble(numTausForFit,
776 &taus_[beginFitIndex_],
777 &msds[beginFitIndex_],
778 &molecule.diffusionCoefficient,
780 &correlationCoefficient,
782 molecule.diffusionCoefficient *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
786 void Msd::writeOutput()
788 // AnalysisData currently doesn't support changing column counts after analysis has started.
789 // We can't determine the number of tau values until the trajectory is fully read, so analysis
790 // data construction and plotting are done here.
791 AnalysisDataPlotModulePointer msdPlotModule(new AnalysisDataPlotModule(plotSettings_));
792 msdPlotModule->setFileName(output_);
793 msdPlotModule->setTitle("Mean Squared Displacement");
794 msdPlotModule->setXLabel("tau (ps)");
795 msdPlotModule->setYLabel(R"(MSD (nm\\S2\\N))");
796 msdPlotModule->setYFormat(10, 6, 'g');
797 for (const auto& group : groupData_)
799 const real D = group.diffusionCoefficient;
800 if (D > 0.01 && D < 1e4)
802 msdPlotModule->appendLegend(formatString(
803 "D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)", group.sel.name(), D, group.sigma));
807 msdPlotModule->appendLegend(formatString(
808 "D[%10s] = %.4g (+/- %.4f) (1e-5 cm^2/s)", group.sel.name(), D, group.sigma));
811 msdPlotData_.addModule(msdPlotModule);
812 msdPlotData_.setDataSetCount(groupData_.size());
813 for (size_t i = 0; i < groupData_.size(); i++)
815 msdPlotData_.setColumnCount(i, 1);
817 AnalysisDataHandle dh = msdPlotData_.startData({});
818 for (size_t tauIndex = 0; tauIndex < taus_.size(); tauIndex++)
820 dh.startFrame(tauIndex, taus_[tauIndex]);
821 for (size_t dataSetIndex = 0; dataSetIndex < groupData_.size(); dataSetIndex++)
823 dh.selectDataSet(dataSetIndex);
824 dh.setPoint(0, groupData_[dataSetIndex].msdSums[tauIndex]);
832 AnalysisDataPlotModulePointer molPlotModule(new AnalysisDataPlotModule(plotSettings_));
833 molPlotModule->setFileName(moleculeOutput_);
834 molPlotModule->setTitle("Mean Squared Displacement / Molecule");
835 molPlotModule->setXLabel("Molecule");
836 molPlotModule->setYLabel("D(1e-5 cm^2/s)");
837 molPlotModule->setYFormat(10, 0, 'g');
838 msdMoleculePlotData_.addModule(molPlotModule);
839 msdMoleculePlotData_.setDataSetCount(1);
840 msdMoleculePlotData_.setColumnCount(0, 1);
841 AnalysisDataHandle molDh = msdMoleculePlotData_.startData({});
842 for (size_t moleculeIndex = 0; moleculeIndex < molecules_.size(); moleculeIndex++)
844 molDh.startFrame(moleculeIndex, moleculeIndex);
845 molDh.setPoint(0, molecules_[moleculeIndex].diffusionCoefficient);
853 const char MsdInfo::name[] = "msd";
854 const char MsdInfo::shortDescription[] = "Compute mean squared displacements";
855 TrajectoryAnalysisModulePointer MsdInfo::create()
857 return TrajectoryAnalysisModulePointer(new Msd);
861 } // namespace analysismodules