Allow density-guided simulations to translate structure

Density-guided simulations acquire a new feature that allows a
user-defined shift on the structure before forces are calculated.

First step to address #3581

diff --git a/docs/reference-manual/special/density-guided-simulation.rst b/docs/reference-manual/special/density-guided-simulation.rst
index e442dfad00bf89bf3211d02ac1cfc8b0e8c53955..cb5315bf7b17989c7b0a716cdc7e0375a7335a80 100644 (file)
--- a/docs/reference-manual/special/density-guided-simulation.rst
+++ b/docs/reference-manual/special/density-guided-simulation.rst
@@ -202,6 +202,19 @@ the force constant is increased by multiplying by :math:`1+2\delta t_{\mathrm{de
 Note that adaptive force scaling does not conserve energy and will ultimately lead to very high
 forces when similarity cannot be increased further.
 
+Aligning input structure and density
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+To align input structure and density data, a shift vector 
+:math:`v_{\mathrm{shift}}` may be defined that translates the input structure
+atom coordinates before evaluating density-guided-simulation energies and forces,
+so that 
+
+.. math:: U = U_{\mathrm{forcefield}}(\mathbf{r}) - k S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r+v_{\mathrm{shift}}})]\,\mathrm{.}
+          :label:eqndensnine
+
+.. math:: \mathbf{F}_{\mathrm{density}} = k \nabla_{\mathbf{r}} S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{r+v_{\mathrm{shift}}})]\,\mathrm{.}
+          :label:eqndensten
+
 Future developments
 ^^^^^^^^^^^^^^^^^^^
 

diff --git a/docs/release-notes/2021/major/features.rst b/docs/release-notes/2021/major/features.rst
index ac0c3462de255fd9fdff3cb4614b9ce7eb06d6da..27a93ec0434aaaad790932ca924db04a4d4ea3a7 100644 (file)
--- a/docs/release-notes/2021/major/features.rst
+++ b/docs/release-notes/2021/major/features.rst
@@ -13,6 +13,17 @@ Virtual site with single constructing atom
 Added a virtual site that is constructed on top if its single constructing
 atom. This can be useful for free-energy calculations.
 
+Density-guided simulations can apply shift vector to structures
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The new mdp option "density-guided-simulation-shift-vector" defines a
+shift vector that shifts the density-guided simulation group before the 
+density forces are evaluated. With a known shift vector that aligns structure
+and input density, this feature enables structure refinement to non-aligned
+densities without the need to manipulate the input density data or structure.
+A typical use case are membrane-embedded proteins which cannot easily be
+shifted.
+
 Lower energy drift due to SETTLE
 """"""""""""""""""""""""""""""""
 

diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst
index 90addbe69614216fd87c2e4b774a4b921b56abde..cd78e3151ec8c5cb493b2662ebf853ac53e23543 100644 (file)
--- a/docs/user-guide/mdp-options.rst
+++ b/docs/user-guide/mdp-options.rst
@@ -3228,6 +3228,14 @@ electron-microscopy experiments. (See the `reference manual`_ for details)
    (4) [ps] Couple force constant to increase in similarity with reference density
    with this time constant. Larger times result in looser coupling.
 
+.. mdp:: density-guided-simulation-shift-vector
+
+   (0,0,0) [nm] Add this vector to all atoms in the 
+   density-guided-simulation-group before calculating forces and energies for
+   density-guided-simulations. Affects only the density-guided-simulation forces
+   and energies. Corresponds to a shift of the input density in the opposite
+   direction by (-1) * density-guided-simulation-shift-vector.
+
 User defined thingies
 ^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/applied_forces/densityfitting/densityfittingforceprovider.cpp b/src/gromacs/applied_forces/densityfitting/densityfittingforceprovider.cpp
index f98eb4c092cea7f1f4af37136116690bb7bb4820..b2963f6e03b58d8e786c8f1c516d48f277030392 100644 (file)
--- a/src/gromacs/applied_forces/densityfitting/densityfittingforceprovider.cpp
+++ b/src/gromacs/applied_forces/densityfitting/densityfittingforceprovider.cpp
@@ -57,6 +57,7 @@
 #include "gromacs/mdtypes/enerdata.h"
 #include "gromacs/mdtypes/forceoutput.h"
 #include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/strconvert.h"
 
 #include "densityfittingamplitudelookup.h"
 #include "densityfittingparameters.h"
@@ -179,6 +180,9 @@ private:
 
     //! Optionally scale the force according to a moving average of the similarity
     std::optional<ExponentialMovingAverage> expAverageSimilarity_;
+
+    //! Optionally translate the structure
+    std::optional<TranslateAndScale> translate_;
 };
 
 DensityFittingForceProvider::Impl::~Impl() = default;
@@ -214,6 +218,17 @@ DensityFittingForceProvider::Impl::Impl(const DensityFittingParameters&
                         / (simulationTimeStep * parameters_.calculationIntervalInSteps_),
                 state.exponentialMovingAverageState_));
     }
+
+    // set up optional coordinate translation if the translation string contains a vector
+    const std::optional<std::array<real, 3>> translationParametersAsArray =
+            parsedArrayFromInputString<real, 3>(parameters_.translationString_);
+    if (translationParametersAsArray)
+    {
+        translate_.emplace(RVec(1, 1, 1), RVec((*translationParametersAsArray)[XX],
+                                               (*translationParametersAsArray)[YY],
+                                               (*translationParametersAsArray)[ZZ]));
+    }
+
     referenceDensityCenter_ = { real(referenceDensity.extent(XX)) / 2,
                                 real(referenceDensity.extent(YY)) / 2,
                                 real(referenceDensity.extent(ZZ)) / 2 };
@@ -246,6 +261,12 @@ void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput
                    std::begin(transformedCoordinates_),
                    [&forceProviderInput](int index) { return forceProviderInput.x_[index]; });
 
+    // optionally apply additional structure transformations
+    if (translate_)
+    {
+        (*translate_)(transformedCoordinates_);
+    }
+
     // pick periodic image that is closest to the center of the reference density
     {
         t_pbc pbc;

diff --git a/src/gromacs/applied_forces/densityfitting/densityfittingoptions.cpp b/src/gromacs/applied_forces/densityfitting/densityfittingoptions.cpp
index 6be0dbf1bad1ebd48ebd3d1a2f19351cae898bdf..298d6913239bd4b636e8851996c4649e77ba339c 100644 (file)
--- a/src/gromacs/applied_forces/densityfitting/densityfittingoptions.cpp
+++ b/src/gromacs/applied_forces/densityfitting/densityfittingoptions.cpp
@@ -146,6 +146,13 @@ void DensityFittingOptions::initMdpTransform(IKeyValueTreeTransformRules* rules)
     densityfittingMdpTransformFromString<bool>(rules, &fromStdString<bool>, c_adaptiveForceScalingTag_);
     densityfittingMdpTransformFromString<real>(rules, &fromStdString<real>,
                                                c_adaptiveForceScalingTimeConstantTag_);
+    const auto& stringRVecToStringRVecWithCheck = [](const std::string& str) {
+        return stringIdentityTransformWithArrayCheck<real, 3>(
+                str, "Reading three real values as vector while parsing the .mdp input failed in "
+                             + DensityFittingModuleInfo::name_ + ".");
+    };
+    densityfittingMdpTransformFromString<std::string>(rules, stringRVecToStringRVecWithCheck,
+                                                      c_translationTag_);
 }
 
 //! Name the methods that may be used to evaluate similarity between densities
@@ -237,6 +244,7 @@ void DensityFittingOptions::initMdpOptions(IOptionsContainerWithSections* option
             BooleanOption(c_adaptiveForceScalingTag_.c_str()).store(&parameters_.adaptiveForceScaling_));
     section.addOption(RealOption(c_adaptiveForceScalingTimeConstantTag_.c_str())
                               .store(&parameters_.adaptiveForceScalingTimeConstant_));
+    section.addOption(StringOption(c_translationTag_.c_str()).store(&parameters_.translationString_));
 }
 
 bool DensityFittingOptions::active() const

diff --git a/src/gromacs/applied_forces/densityfitting/densityfittingoptions.h b/src/gromacs/applied_forces/densityfitting/densityfittingoptions.h
index be2f65b61974a26672be6d7eb602000bfd135aed..b361eea5118debfd6de7859a3f3df8fe59170b7d 100644 (file)
--- a/src/gromacs/applied_forces/densityfitting/densityfittingoptions.h
+++ b/src/gromacs/applied_forces/densityfitting/densityfittingoptions.h
@@ -136,6 +136,7 @@ private:
     const std::string c_adaptiveForceScalingTimeConstantTag_ =
             "adaptive-force-scaling-time-constant";
 
+    const std::string c_translationTag_ = "shift-vector";
 
     DensityFittingParameters parameters_;
 };

diff --git a/src/gromacs/applied_forces/densityfitting/densityfittingparameters.h b/src/gromacs/applied_forces/densityfitting/densityfittingparameters.h
index 17a70545f0836ec2105265778e0f8e1cb0ea2d42..5cfef20dc0f582c20797a0885256fbee8bb5098f 100644 (file)
--- a/src/gromacs/applied_forces/densityfitting/densityfittingparameters.h
+++ b/src/gromacs/applied_forces/densityfitting/densityfittingparameters.h
@@ -81,6 +81,8 @@ struct DensityFittingParameters
     bool adaptiveForceScaling_ = false;
     //! The time constant for the adaptive force scaling in ps
     real adaptiveForceScalingTimeConstant_ = 4;
+    //! Translation of the structure, so that the coordinates that are fitted are x+translation
+    std::string translationString_ = "";
 };
 
 /*!\brief Check if two structs holding density fitting parameters are equal.

diff --git a/src/programs/mdrun/tests/densityfittingmodule.cpp b/src/programs/mdrun/tests/densityfittingmodule.cpp
index ab7cc8798f04c283d435349198830bcaeef0314f..91c74f264816957578701f10cc620c3d80856c03 100644 (file)
--- a/src/programs/mdrun/tests/densityfittingmodule.cpp
+++ b/src/programs/mdrun/tests/densityfittingmodule.cpp
@@ -131,6 +131,10 @@ public:
     const std::string mdpSkipDensityfittingEveryOtherStep_ = formatString(
             "nstenergy = 2\n"
             "density-guided-simulation-nst = 2\n");
+    const std::string mdpTranslationSet_ =
+            formatString("density-guided-simulation-shift-vector = 0.1 -0.2 0.3\n");
+    const std::string mdpTranslationSetWrongValues_ =
+            formatString("density-guided-simulation-shift-vector = 0.1 -0.2\n");
     //! The command line to call mdrun
     CommandLine commandLineForMdrun_;
 };
@@ -151,6 +155,34 @@ TEST_F(DensityFittingTest, EnergyMinimizationEnergyCorrectInnerProduct)
     checkMdrun(expectedEnergyTermMagnitude);
 }
 
+/* Fit a subset of three of twelve argon atoms into a reference density
+ * whose origin is offset from the simulation box origin.
+ *
+ * All density fitting mdp parameters are set to defaults
+ */
+TEST_F(DensityFittingTest, EnergyMinimizationEnergyCorrectInnerProductTranslation)
+{
+    runner_.useStringAsMdpFile(mdpEminDensfitYesUnsetValues + mdpTranslationSet_);
+
+    ASSERT_EQ(0, runner_.callGrompp());
+    ASSERT_EQ(0, runner_.callMdrun(commandLineForMdrun_));
+
+    const real expectedEnergyTermMagnitude = -8991;
+    checkMdrun(expectedEnergyTermMagnitude);
+}
+
+/* Fit a subset of three of twelve argon atoms into a reference density
+ * whose origin is offset from the simulation box origin.
+ *
+ * All density fitting mdp parameters are set to defaults
+ */
+TEST_F(DensityFittingTest, EnergyMinimizationEnergyTranslationParametersOff)
+{
+    runner_.useStringAsMdpFile(mdpEminDensfitYesUnsetValues + mdpTranslationSetWrongValues_);
+
+    GMX_EXPECT_DEATH_IF_SUPPORTED(runner_.callGrompp(), ".*Reading three real values.*");
+}
+
 /* Like above, but with as many parameters reversed as possible
  */
 TEST_F(DensityFittingTest, EnergyMinimizationEnergyCorrectForRelativeEntropy)

diff --git a/src/programs/mdrun/tests/refdata/DensityFittingTest_EnergyMinimizationEnergyCorrectInnerProductTranslation.xml b/src/programs/mdrun/tests/refdata/DensityFittingTest_EnergyMinimizationEnergyCorrectInnerProductTranslation.xml
new file mode 100644 (file)
index 0000000..4d3ec70
--- /dev/null
+++ b/src/programs/mdrun/tests/refdata/DensityFittingTest_EnergyMinimizationEnergyCorrectInnerProductTranslation.xml
@@ -0,0 +1,14 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <Energy Name="Potential">
+    <Real Name="Time 0.000000 Step 0 in frame 0">-8991.123</Real>
+    <Real Name="Time 1.000000 Step 1 in frame 1">-9367.3584</Real>
+    <Real Name="Time 2.000000 Step 2 in frame 2">-9820.8447</Real>
+  </Energy>
+  <Energy Name="Density fitting">
+    <Real Name="Time 0.000000 Step 0 in frame 0">-8990.1055</Real>
+    <Real Name="Time 1.000000 Step 1 in frame 1">-9366.3193</Real>
+    <Real Name="Time 2.000000 Step 2 in frame 2">-9819.8525</Real>
+  </Energy>
+</ReferenceData>