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;
435 //! First tau value to fit from for diffusion coefficient, defaults to 0.1 * max tau
436 real beginFit_ = -1.0;
437 //! Final tau value to fit to for diffusion coefficient, defaults to 0.9 * max tau
440 //! All selection group-specific data stored here.
441 std::vector<MsdGroupData> groupData_;
443 //! Time of all frames.
444 std::vector<double> times_;
445 //! Taus for output - won't know the size until the end.
446 std::vector<double> taus_;
447 //! Tau indices for fitting.
448 size_t beginFitIndex_ = 0;
449 size_t endFitIndex_ = 0;
451 // MSD per-molecule stuff
452 //! Are we doing molecule COM-based MSDs?
453 bool molSelected_ = false;
454 //! Per molecule topology information and MSD accumulators.
455 std::vector<MoleculeData> molecules_;
456 //! Atom index -> mol index map
457 std::vector<int> moleculeIndexMappings_;
460 AnalysisData msdPlotData_;
461 AnalysisData msdMoleculePlotData_;
463 AnalysisDataPlotSettings plotSettings_;
464 //! Per-tau MSDs for each selected group
466 //! Per molecule diffusion coefficients if -mol is selected.
467 std::string moleculeOutput_;
471 Msd::Msd() = default;
473 void Msd::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
475 static const char* const desc[] = {
476 "[THISMODULE] computes the mean square displacement (MSD) of atoms from",
477 "a set of initial positions. This provides an easy way to compute",
478 "the diffusion constant using the Einstein relation.",
479 "The time between the reference points for the MSD calculation",
480 "is set with [TT]-trestart[tt].",
481 "The diffusion constant is calculated by least squares fitting a",
482 "straight line (D*t + c) through the MSD(t) from [TT]-beginfit[tt] to",
483 "[TT]-endfit[tt] (note that t is time from the reference positions,",
484 "not simulation time). An error estimate given, which is the difference",
485 "of the diffusion coefficients obtained from fits over the two halves",
486 "of the fit interval.[PAR]",
487 "There are three, mutually exclusive, options to determine different",
488 "types of mean square displacement: [TT]-type[tt], [TT]-lateral[tt]",
489 "and [TT]-ten[tt]. Option [TT]-ten[tt] writes the full MSD tensor for",
490 "each group, the order in the output is: trace xx yy zz yx zx zy.[PAR]",
491 "If [TT]-mol[tt] is set, [THISMODULE] plots the MSD for individual molecules",
492 "(including making molecules whole across periodic boundaries): ",
493 "for each individual molecule a diffusion constant is computed for ",
494 "its center of mass. The chosen index group will be split into ",
495 "molecules. With -mol, only one index group can be selected.[PAR]",
496 "The diffusion coefficient is determined by linear regression of the MSD.",
497 "When [TT]-beginfit[tt] is -1, fitting starts at 10%",
498 "and when [TT]-endfit[tt] is -1, fitting goes to 90%.",
499 "Using this option one also gets an accurate error estimate",
500 "based on the statistics between individual molecules.",
501 "Note that this diffusion coefficient and error estimate are only",
502 "accurate when the MSD is completely linear between",
503 "[TT]-beginfit[tt] and [TT]-endfit[tt].[PAR]",
505 settings->setHelpText(desc);
508 options->addOption(SelectionOption("sel")
509 .storeVector(&selections_)
513 .description("Selections to compute MSDs for from the reference"));
515 // Select MSD type - defaults to 3D if neither option is selected.
516 EnumerationArray<SingleDimDiffType, const char*> enumTypeNames = { "x", "y", "z", "unused" };
517 EnumerationArray<TwoDimDiffType, const char*> enumLateralNames = { "x", "y", "z", "unused" };
518 options->addOption(EnumOption<SingleDimDiffType>("type")
519 .enumValue(enumTypeNames)
520 .store(&singleDimType_)
521 .defaultValue(SingleDimDiffType::Unused));
522 options->addOption(EnumOption<TwoDimDiffType>("lateral")
523 .enumValue(enumLateralNames)
525 .defaultValue(TwoDimDiffType::Unused));
527 options->addOption(DoubleOption("trestart")
528 .description("Time between restarting points in trajectory (ps)")
532 RealOption("beginfit").description("Time point at which to start fitting.").store(&beginFit_));
533 options->addOption(RealOption("endfit").description("End time for fitting.").store(&endFit_));
536 options->addOption(FileNameOption("o")
537 .filetype(OptionFileType::Plot)
540 .defaultBasename("msdout")
541 .description("MSD output"));
543 FileNameOption("mol")
544 .filetype(OptionFileType::Plot)
546 .store(&moleculeOutput_)
547 .storeIsSet(&molSelected_)
548 .defaultBasename("diff_mol")
549 .description("Report diffusion coefficients for each molecule in selection"));
552 void Msd::optionsFinished(TrajectoryAnalysisSettings gmx_unused* settings)
554 if (singleDimType_ != SingleDimDiffType::Unused && twoDimType_ != TwoDimDiffType::Unused)
556 std::string errorMessage =
557 "Options -type and -lateral are mutually exclusive. Choose one or neither (for 3D "
559 GMX_THROW(InconsistentInputError(errorMessage.c_str()));
561 if (selections_.size() > 1 && molSelected_)
563 std::string errorMessage =
564 "Cannot have multiple groups selected with -sel when using -mol.";
565 GMX_THROW(InconsistentInputError(errorMessage.c_str()));
570 void Msd::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
572 plotSettings_ = settings.plotSettings();
574 // Enumeration helpers for dispatching the right MSD calculation type.
575 const EnumerationArray<SingleDimDiffType, decltype(&calcAverageDisplacement<true, true, true>)>
576 oneDimensionalMsdFunctions = { &calcAverageDisplacement<true, false, false>,
577 &calcAverageDisplacement<false, true, false>,
578 &calcAverageDisplacement<false, false, true> };
579 const EnumerationArray<TwoDimDiffType, decltype(&calcAverageDisplacement<true, true, true>)>
580 twoDimensionalMsdFunctions = { calcAverageDisplacement<false, true, true>,
581 calcAverageDisplacement<true, false, true>,
582 calcAverageDisplacement<true, true, false> };
584 // Parse dimensionality and assign the MSD calculating function.
585 // Note if we don't hit either of these cases, we're computing 3D MSDs.
586 if (singleDimType_ != SingleDimDiffType::Unused)
588 calcMsd_ = oneDimensionalMsdFunctions[singleDimType_];
589 diffusionCoefficientDimensionFactor_ = c_1DdiffusionDimensionFactor;
591 else if (twoDimType_ != TwoDimDiffType::Unused)
593 calcMsd_ = twoDimensionalMsdFunctions[twoDimType_];
594 diffusionCoefficientDimensionFactor_ = c_2DdiffusionDimensionFactor;
597 // TODO validate that we have mol info and not atom only - and masses, and topology.
600 Selection& sel = selections_[0];
601 const int nMol = sel.initOriginalIdsToGroup(top.mtop(), INDEX_MOL);
603 gmx::ArrayRef<const int> mappedIds = selections_[0].mappedIds();
604 moleculeIndexMappings_.resize(selections_[0].posCount());
605 std::copy(mappedIds.begin(), mappedIds.end(), moleculeIndexMappings_.begin());
607 // Precalculate each molecules mass for speeding up COM calculations.
608 ArrayRef<const real> masses = sel.masses();
610 molecules_.resize(nMol);
611 for (int i = 0; i < sel.posCount(); i++)
613 molecules_[mappedIds[i]].atomCount++;
614 molecules_[mappedIds[i]].mass += masses[i];
618 // Accumulated frames and results
619 for (const Selection& sel : selections_)
621 groupData_.emplace_back(sel, molecules_, moleculeIndexMappings_);
625 void Msd::initAfterFirstFrame(const TrajectoryAnalysisSettings gmx_unused& settings, const t_trxframe& fr)
627 t0_ = std::round(fr.time);
631 void Msd::analyzeFrame(int gmx_unused frameNumber,
632 const t_trxframe& frame,
634 TrajectoryAnalysisModuleData* pdata)
636 const real time = std::round(frame.time);
637 // Need to populate dt on frame 2;
638 if (!dt_.has_value() && !times_.empty())
640 dt_ = time - times_[0];
643 // Each frame gets an entry in times, but frameTimes only updates if we're at a restart.
644 times_.push_back(time);
646 // Each frame will get a tau between it and frame 0, and all other frame combos should be
648 // \todo this will no longer hold exactly when maxtau is added
649 taus_.push_back(time - times_[0]);
651 for (MsdGroupData& msdData : groupData_)
653 const Selection& sel = pdata->parallelSelection(msdData.sel);
655 ArrayRef<const RVec> coords = msdData.coordinateManager_.buildCoordinates(sel, pbc);
657 // For each preceding frame, calculate tau and do comparison.
658 for (size_t i = 0; i < msdData.frames.size(); i++)
660 // If dt > trestart, the tau increment will be dt rather than trestart.
661 double tau = time - (t0_ + std::max<double>(trestart_, *dt_) * i);
662 int64_t tauIndex = gmx::roundToInt64(tau / *dt_);
663 msdData.msds[tauIndex].push_back(calcMsd_(coords, msdData.frames[i]));
665 for (size_t molInd = 0; molInd < molecules_.size(); molInd++)
667 molecules_[molInd].msdData[tauIndex].push_back(
668 calcMsd_(arrayRefFromArray(&coords[molInd], 1),
669 arrayRefFromArray(&msdData.frames[i][molInd], 1)));
674 // We only store the frame for the future if it's a restart per -trestart.
675 if (bRmod(time, t0_, trestart_))
677 msdData.frames.emplace_back(coords.begin(), coords.end());
682 //! Calculate the tau index for fitting. If userFitTau < 0, uses the default fraction of max tau.
683 static size_t calculateFitIndex(const int userFitTau,
684 const double defaultTauFraction,
690 return gmx::roundToInt((numTaus - 1) * defaultTauFraction);
692 return std::min<size_t>(numTaus - 1, gmx::roundToInt(static_cast<double>(userFitTau) / dt));
696 void Msd::finishAnalysis(int gmx_unused nframes)
698 static constexpr double c_defaultStartFitIndexFraction = 0.1;
699 static constexpr double c_defaultEndFitIndexFraction = 0.9;
700 beginFitIndex_ = calculateFitIndex(beginFit_, c_defaultStartFitIndexFraction, taus_.size(), *dt_);
701 endFitIndex_ = calculateFitIndex(endFit_, c_defaultEndFitIndexFraction, taus_.size(), *dt_);
702 const int numTausForFit = 1 + endFitIndex_ - beginFitIndex_;
704 // These aren't used, except for correlationCoefficient, which is used to estimate error if
705 // enough points are available.
706 real b = 0.0, correlationCoefficient = 0.0, chiSquared = 0.0;
708 for (MsdGroupData& msdData : groupData_)
710 msdData.msdSums = msdData.msds.averageMsds();
712 if (numTausForFit >= 4)
714 const int halfNumTaus = numTausForFit / 2;
715 const int secondaryStartIndex = beginFitIndex_ + halfNumTaus;
716 // Split the fit in 2, and compare the results of each fit;
717 real a = 0.0, a2 = 0.0;
718 lsq_y_ax_b_xdouble(halfNumTaus,
719 &taus_[beginFitIndex_],
720 &msdData.msdSums[beginFitIndex_],
723 &correlationCoefficient,
725 lsq_y_ax_b_xdouble(halfNumTaus,
726 &taus_[secondaryStartIndex],
727 &msdData.msdSums[secondaryStartIndex],
730 &correlationCoefficient,
732 msdData.sigma = std::abs(a - a2);
734 lsq_y_ax_b_xdouble(numTausForFit,
735 &taus_[beginFitIndex_],
736 &msdData.msdSums[beginFitIndex_],
737 &msdData.diffusionCoefficient,
739 &correlationCoefficient,
741 msdData.diffusionCoefficient *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
742 msdData.sigma *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
745 for (MoleculeData& molecule : molecules_)
747 std::vector<real> msds = molecule.msdData.averageMsds();
748 lsq_y_ax_b_xdouble(numTausForFit,
749 &taus_[beginFitIndex_],
750 &msds[beginFitIndex_],
751 &molecule.diffusionCoefficient,
753 &correlationCoefficient,
755 molecule.diffusionCoefficient *= c_diffusionConversionFactor / diffusionCoefficientDimensionFactor_;
759 void Msd::writeOutput()
761 // AnalysisData currently doesn't support changing column counts after analysis has started.
762 // We can't determine the number of tau values until the trajectory is fully read, so analysis
763 // data construction and plotting are done here.
764 AnalysisDataPlotModulePointer msdPlotModule(new AnalysisDataPlotModule(plotSettings_));
765 msdPlotModule->setFileName(output_);
766 msdPlotModule->setTitle("Mean Squared Displacement");
767 msdPlotModule->setXLabel("tau (ps)");
768 msdPlotModule->setYLabel(R"(MSD (nm\\S2\\N))");
769 msdPlotModule->setYFormat(10, 6, 'g');
770 for (const auto& group : groupData_)
772 const real D = group.diffusionCoefficient;
773 if (D > 0.01 && D < 1e4)
775 msdPlotModule->appendLegend(formatString(
776 "D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)", group.sel.name(), D, group.sigma));
780 msdPlotModule->appendLegend(formatString(
781 "D[%10s] = %.4g (+/- %.4f) (1e-5 cm^2/s)", group.sel.name(), D, group.sigma));
784 msdPlotData_.addModule(msdPlotModule);
785 msdPlotData_.setDataSetCount(groupData_.size());
786 for (size_t i = 0; i < groupData_.size(); i++)
788 msdPlotData_.setColumnCount(i, 1);
790 AnalysisDataHandle dh = msdPlotData_.startData({});
791 for (size_t tauIndex = 0; tauIndex < taus_.size(); tauIndex++)
793 dh.startFrame(tauIndex, taus_[tauIndex]);
794 for (size_t dataSetIndex = 0; dataSetIndex < groupData_.size(); dataSetIndex++)
796 dh.selectDataSet(dataSetIndex);
797 dh.setPoint(0, groupData_[dataSetIndex].msdSums[tauIndex]);
805 AnalysisDataPlotModulePointer molPlotModule(new AnalysisDataPlotModule(plotSettings_));
806 molPlotModule->setFileName(moleculeOutput_);
807 molPlotModule->setTitle("Mean Squared Displacement / Molecule");
808 molPlotModule->setXLabel("Molecule");
809 molPlotModule->setYLabel("D(1e-5 cm^2/s)");
810 molPlotModule->setYFormat(10, 0, 'g');
811 msdMoleculePlotData_.addModule(molPlotModule);
812 msdMoleculePlotData_.setDataSetCount(1);
813 msdMoleculePlotData_.setColumnCount(0, 1);
814 AnalysisDataHandle molDh = msdMoleculePlotData_.startData({});
815 for (size_t moleculeIndex = 0; moleculeIndex < molecules_.size(); moleculeIndex++)
817 molDh.startFrame(moleculeIndex, moleculeIndex);
818 molDh.setPoint(0, molecules_[moleculeIndex].diffusionCoefficient);
826 const char MsdInfo::name[] = "msd";
827 const char MsdInfo::shortDescription[] = "Compute mean squared displacements";
828 TrajectoryAnalysisModulePointer MsdInfo::create()
830 return TrajectoryAnalysisModulePointer(new Msd);
834 } // namespace analysismodules