From: Kevin Boyd Date: Tue, 26 Oct 2021 18:13:12 +0000 (+0000) Subject: Add -maxtau option to gmx msd X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?p=alexxy%2Fgromacs.git;a=commitdiff_plain;h=d64102a1b1a599d7905ecdaa8ea47b2fcdce244a Add -maxtau option to gmx msd --- diff --git a/src/gromacs/trajectoryanalysis/modules/msd.cpp b/src/gromacs/trajectoryanalysis/modules/msd.cpp index 0a0dd25616..c227bc6c7f 100644 --- a/src/gromacs/trajectoryanalysis/modules/msd.cpp +++ b/src/gromacs/trajectoryanalysis/modules/msd.cpp @@ -431,6 +431,10 @@ private: double t0_ = 0; //! Inter-frame delta-t std::optional dt_ = std::nullopt; + //! Inter-frame maximum time delta. + double maxTau_ = std::numeric_limits::max(); + //! Tracks the index of the first valid frame. Starts at 0, increases as old frames are deleted + size_t firstValidFrame_ = 0; //! First tau value to fit from for diffusion coefficient, defaults to 0.1 * max tau real beginFit_ = -1.0; @@ -501,6 +505,14 @@ void Msd::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* se "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]", + "By default, [THISMODULE] compares all trajectory frames against every frame stored at", + "[TT]-trestart[TT] intervals, so the number of frames stored scales linearly with the", + "number of frames processed. This can lead to long analysis times and out-of-memory errors", + "for long/large trajectories, and often the data at higher time deltas lacks sufficient ", + "sampling, often manifesting as a wobbly line on the MSD plot after a straighter region at", + "lower time deltas. The [TT]-maxtau[TT] option can be used to cap the maximum time delta", + "for frame comparison, which may improve performance and can be used to avoid", + "out-of-memory issues.[PAR]" }; settings->setHelpText(desc); @@ -528,6 +540,10 @@ void Msd::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* se .description("Time between restarting points in trajectory (ps)") .defaultValue(10.0) .store(&trestart_)); + options->addOption( + DoubleOption("maxtau") + .description("Maximum time delta between frames to calculate MSDs for (ps)") + .store(&maxTau_)); options->addOption( RealOption("beginfit").description("Time point at which to start fitting.").store(&beginFit_)); options->addOption(RealOption("endfit").description("End time for fitting.").store(&endFit_)); @@ -645,8 +661,10 @@ void Msd::analyzeFrame(int gmx_unused frameNumber, // 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]); + if (const double tau = time - times_[0]; tau <= maxTau_) + { + taus_.push_back(time - times_[0]); + } for (MsdGroupData& msdData : groupData_) { @@ -655,10 +673,19 @@ void Msd::analyzeFrame(int gmx_unused frameNumber, ArrayRef 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++) + for (size_t i = firstValidFrame_; i < msdData.frames.size(); i++) { // If dt > trestart, the tau increment will be dt rather than trestart. - double tau = time - (t0_ + std::max(trestart_, *dt_) * i); + const double tau = time - (t0_ + std::max(trestart_, *dt_) * i); + if (tau > maxTau_) + { + // The (now empty) entry is no longer needed, so over time the outer vector will + // grow with extraneous empty elements persisting, but the alternative would require + // some more complicated remapping of tau to frame index. + msdData.frames[i].clear(); + firstValidFrame_ = i + 1; + continue; + } int64_t tauIndex = gmx::roundToInt64(tau / *dt_); msdData.msds[tauIndex].push_back(calcMsd_(coords, msdData.frames[i])); diff --git a/src/gromacs/trajectoryanalysis/tests/msd.cpp b/src/gromacs/trajectoryanalysis/tests/msd.cpp index 16220eafd5..08d489fc2b 100644 --- a/src/gromacs/trajectoryanalysis/tests/msd.cpp +++ b/src/gromacs/trajectoryanalysis/tests/msd.cpp @@ -225,6 +225,16 @@ TEST_F(MsdModuleTest, oneDimensionalDiffusion) runTest(CommandLine(cmdline)); } +TEST_F(MsdModuleTest, oneDimensionalDiffusionWithMaxTau) +{ + setInputFile("-f", "msd_traj.xtc"); + setInputFile("-s", "msd_coords.gro"); + setInputFile("-n", "msd.ndx"); + const char* const cmdline[] = { "-trestart", "200", "-type", "x", "-sel", "0", "-maxtau", "5" }; + 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 diff --git a/src/gromacs/trajectoryanalysis/tests/refdata/MsdModuleTest_oneDimensionalDiffusionWithMaxTau.xml b/src/gromacs/trajectoryanalysis/tests/refdata/MsdModuleTest_oneDimensionalDiffusionWithMaxTau.xml new file mode 100644 index 0000000000..c74497d964 --- /dev/null +++ b/src/gromacs/trajectoryanalysis/tests/refdata/MsdModuleTest_oneDimensionalDiffusionWithMaxTau.xml @@ -0,0 +1,49 @@ + + + + -trestart 200 -type x -sel 0 -maxtau 5 + + + + + title "Mean Squared Displacement" + xaxis label "tau (ps)" + yaxis label "MSD (nm\\S2\\N)" + TYPE xy + s0 legend "D[ particles] = 8.0000 (+/- 0.0000) (1e-5 cm^2/s)" + + + + 2 + 0.000 + 0 + + + 2 + 1.000 + 0.016 + + + 2 + 2.000 + 0.032 + + + 2 + 3.000 + 0.048 + + + 2 + 4.000 + 0.064 + + + 2 + 5.000 + 0.08 + + + + +