Resolves part of #3581.
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
+To align input structure and density data, a transformation matrix
+:math:`\mathbf{A}` and shift vector :math:`\mathbf{v_{\mathrm{shift}}}` may be
+defined that transform 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{.}
+.. math:: U = U_{\mathrm{forcefield}}(\mathbf{r}) - k S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{A 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{.}
+.. math:: \mathbf{F}_{\mathrm{density}} = k \nabla_{\mathbf{r}} S[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}\!(\mathbf{A r+v_{\mathrm{shift}}})]\,\mathrm{.}
:label:eqndensten
Future developments
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
-"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+Density-guided simulations can apply matrix multiplication and 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.
+The mdp option "density-guided-simulation-transformation-matrix" allows to
+define a matrix with which to multiply the structure coordinates, before the shift
+vector is applied. This allows arbitrary rotation, skewing and scaling of input
+structures with respect to the input densities.
A typical use case are membrane-embedded proteins which cannot easily be
-shifted.
+shifted and rotated within membranes.
Lower energy drift due to SETTLE
""""""""""""""""""""""""""""""""
and energies. Corresponds to a shift of the input density in the opposite
direction by (-1) * density-guided-simulation-shift-vector.
+.. mdp:: density-guided-simulation-transformation-matrix
+
+ (1,0,0,0,1,0,0,0,1) Multiply all atoms with this matrix 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 transformation of the input density in the
+ opposite by the inverse of this matrix. The matrix is given in row-major order.
+
User defined thingies
^^^^^^^^^^^^^^^^^^^^^
std::optional<ExponentialMovingAverage> expAverageSimilarity_;
//! Optionally translate the structure
- std::optional<TranslateAndScale> translate_;
+ std::optional<AffineTransformation> affineTransformation_;
};
DensityFittingForceProvider::Impl::~Impl() = default;
// 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)
+ // set up optional coordinate transformation if the transformation string contains data
+ const std::optional<std::array<real, 9>> transformationMatrixParametersAsArray =
+ parsedArrayFromInputString<real, 9>(parameters_.transformationMatrixString_);
+ if (translationParametersAsArray || transformationMatrixParametersAsArray)
{
- translate_.emplace(RVec(1, 1, 1), RVec((*translationParametersAsArray)[XX],
- (*translationParametersAsArray)[YY],
- (*translationParametersAsArray)[ZZ]));
+ Matrix3x3 translationMatrix = transformationMatrixParametersAsArray.has_value()
+ ? *transformationMatrixParametersAsArray
+ : identityMatrix<real, 3>();
+ RVec translationVector = translationParametersAsArray.has_value()
+ ? RVec((*translationParametersAsArray)[XX],
+ (*translationParametersAsArray)[YY],
+ (*translationParametersAsArray)[ZZ])
+ : RVec(0, 0, 0);
+ affineTransformation_.emplace(translationMatrix.asConstView(), translationVector);
}
referenceDensityCenter_ = { real(referenceDensity.extent(XX)) / 2,
std::begin(transformedCoordinates_),
[&forceProviderInput](int index) { return forceProviderInput.x_[index]; });
- // optionally apply additional structure transformations
- if (translate_)
+ // apply additional structure transformations
+ if (affineTransformation_)
{
- (*translate_)(transformedCoordinates_);
+ (*affineTransformation_)(transformedCoordinates_);
}
// pick periodic image that is closest to the center of the reference density
};
densityfittingMdpTransformFromString<std::string>(rules, stringRVecToStringRVecWithCheck,
c_translationTag_);
+
+ const auto& stringMatrixToStringMatrixWithCheck = [](const std::string& str) {
+ return stringIdentityTransformWithArrayCheck<real, 9>(
+ str, "Reading nine real values as vector while parsing the .mdp input failed in "
+ + DensityFittingModuleInfo::name_ + ".");
+ };
+ densityfittingMdpTransformFromString<std::string>(rules, stringMatrixToStringMatrixWithCheck,
+ c_transformationMatrixTag_);
}
//! Name the methods that may be used to evaluate similarity between densities
section.addOption(RealOption(c_adaptiveForceScalingTimeConstantTag_.c_str())
.store(¶meters_.adaptiveForceScalingTimeConstant_));
section.addOption(StringOption(c_translationTag_.c_str()).store(¶meters_.translationString_));
+ section.addOption(
+ StringOption(c_transformationMatrixTag_.c_str()).store(¶meters_.transformationMatrixString_));
}
bool DensityFittingOptions::active() const
const std::string c_translationTag_ = "shift-vector";
+ const std::string c_transformationMatrixTag_ = "transformation-matrix";
+
DensityFittingParameters parameters_;
};
real adaptiveForceScalingTimeConstant_ = 4;
//! Translation of the structure, so that the coordinates that are fitted are x+translation
std::string translationString_ = "";
+ //! Linear transformat of the structure, so that the coordinates that are fitted are Matrix * x
+ std::string transformationMatrixString_ = "";
};
/*!\brief Check if two structs holding density fitting parameters are equal.
const std::string mdpSkipDensityfittingEveryOtherStep_ = formatString(
"nstenergy = 2\n"
"density-guided-simulation-nst = 2\n");
+ //! A properly set shift vector
const std::string mdpTranslationSet_ =
formatString("density-guided-simulation-shift-vector = 0.1 -0.2 0.3\n");
+ //! A shift vector that is lacking an entry
const std::string mdpTranslationSetWrongValues_ =
formatString("density-guided-simulation-shift-vector = 0.1 -0.2\n");
+ //! A 45 degree rotation around the y axis expressed as matrix transformation
+ const std::string mdpTransformationMatrix1degAroundY_ = formatString(
+ "density-guided-simulation-transformation-matrix = 0.9998477 0.0000000 0.0174524 "
+ "0.0000000 1.0000000 0.0000000 -0.0174524 0.0000000 0.9998477 \n");
+ //! The identity matrix as transformation matrix
+ const std::string mdpTransformationMatrixIdentity_ = formatString(
+ "density-guided-simulation-transformation-matrix = 1 0 0 "
+ "0 1 0 0 0 1 \n");
+ //! A transformation matrix string where only eight values are given
+ const std::string mdpTransformationMatrixWrongValues_ = formatString(
+ "density-guided-simulation-transformation-matrix = 0.7071068 0.0000000 0.7071068 "
+ "0.0000000 0.0000000 -0.7071068 0.0000000 0.7071068 \n");
+
//! The command line to call mdrun
CommandLine commandLineForMdrun_;
};
GMX_EXPECT_DEATH_IF_SUPPORTED(runner_.callGrompp(), ".*Reading three real values.*");
}
+/* Fit a subset of three of twelve argon atoms into a reference density
+ * that are rotated around the simulation box origin by a matrix multiplication.
+ *
+ * All density fitting mdp parameters are set to defaults
+ */
+TEST_F(DensityFittingTest, EnergyMinimizationEnergyCorrectInnerProductTranslationAndTransformationMatrix)
+{
+ runner_.useStringAsMdpFile(mdpEminDensfitYesUnsetValues + mdpTranslationSet_
+ + mdpTransformationMatrix1degAroundY_);
+
+ 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, EnergyMinimizationEnergyMatrixTransfromationOff)
+{
+ runner_.useStringAsMdpFile(mdpEminDensfitYesUnsetValues + mdpTransformationMatrixWrongValues_);
+
+ GMX_EXPECT_DEATH_IF_SUPPORTED(runner_.callGrompp(), ".*Reading nine real values.*");
+}
+
+/* Fit a subset of three of twelve argon atoms into a reference density
+ * where the given matrix transformation is the identity transformation.
+ *
+ * All density fitting mdp parameters are set to defaults
+ */
+TEST_F(DensityFittingTest, EnergyMinimizationEnergyCorrectInnerProductIdentityMatrix)
+{
+ runner_.useStringAsMdpFile(mdpEminDensfitYesUnsetValues + mdpTransformationMatrixIdentity_);
+
+ ASSERT_EQ(0, runner_.callGrompp());
+ ASSERT_EQ(0, runner_.callMdrun(commandLineForMdrun_));
+
+ const real expectedEnergyTermMagnitude = -8991;
+ checkMdrun(expectedEnergyTermMagnitude);
+}
+
/* Like above, but with as many parameters reversed as possible
*/
TEST_F(DensityFittingTest, EnergyMinimizationEnergyCorrectForRelativeEntropy)
--- /dev/null
+<?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">-3379.9202</Real>
+ <Real Name="Time 1.000000 Step 1 in frame 1">-3591.7776</Real>
+ <Real Name="Time 2.000000 Step 2 in frame 2">-3856.54</Real>
+ </Energy>
+ <Energy Name="Density fitting">
+ <Real Name="Time 0.000000 Step 0 in frame 0">-3378.9026</Real>
+ <Real Name="Time 1.000000 Step 1 in frame 1">-3590.7383</Real>
+ <Real Name="Time 2.000000 Step 2 in frame 2">-3855.5586</Real>
+ </Energy>
+</ReferenceData>
--- /dev/null
+<?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">-8933.2021</Real>
+ <Real Name="Time 1.000000 Step 1 in frame 1">-9308.2793</Real>
+ <Real Name="Time 2.000000 Step 2 in frame 2">-9760.4385</Real>
+ </Energy>
+ <Energy Name="Density fitting">
+ <Real Name="Time 0.000000 Step 0 in frame 0">-8932.1846</Real>
+ <Real Name="Time 1.000000 Step 1 in frame 1">-9307.2402</Real>
+ <Real Name="Time 2.000000 Step 2 in frame 2">-9759.457</Real>
+ </Energy>
+</ReferenceData>
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff