pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-03 21:44:30 +08:00
Code Issues Packages Projects Releases Wiki Activity

Allow creation of an empty forcefield (#7494)

Browse Source 
* add posibility to create a force field that does not have any contributions.

* move to FFConvenience.h and add C++ tests

* deallocate memory

* use range based for loop

Co-authored-by: Greg Landrum <greg.landrum@gmail.com>

* Making adjustments according to the review of @greglandrum

* undo unintentional formatting

* replace TEST_ASSERT with CHECK

---------

Co-authored-by: Greg Landrum <greg.landrum@gmail.com>
This commit is contained in:
nmaeder
2024-06-08 05:18:22 +02:00
committed by GitHub
parent 93e8a746bf
commit 685e781b54
5 changed files with 137 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
Code/GraphMol/ForceFieldHelpers/CMakeLists.txt
View File

@@ -14,6 +14,9 @@ rdkit_headers(MMFF/AtomTyper.h
rdkit_headers(CrystalFF/TorsionAngleM6.h CrystalFF/TorsionPreferences.h
DEST GraphMol/ForceFieldHelpers/CrystalFF)
rdkit_catch_test(forceFieldHelpersCatch catch_tests.cpp
LINK_LIBRARIES ForceFieldHelpers SmilesParse ForceField )
add_subdirectory(MMFF)
add_subdirectory(UFF)
add_subdirectory(CrystalFF)

42
Code/GraphMol/ForceFieldHelpers/FFConvenience.h
View File

@@ -18,10 +18,10 @@ class ROMol;
namespace ForceFieldsHelper {
namespace detail {
#ifdef RDK_BUILD_THREADSAFE_SSS
inline void OptimizeMoleculeConfsHelper_(ForceFields::ForceField ff, ROMol *mol,
std::vector<std::pair<int, double>> *res,
unsigned int threadIdx,
unsigned int numThreads, int maxIters) {
inline void OptimizeMoleculeConfsHelper_(
ForceFields::ForceField ff, ROMol *mol,
std::vector<std::pair<int, double>> *res, unsigned int threadIdx,
unsigned int numThreads, int maxIters) {
PRECONDITION(mol, "mol must not be nullptr");
PRECONDITION(res, "res must not be nullptr");
PRECONDITION(res->size() >= mol->getNumConformers(),
@@ -43,9 +43,10 @@ inline void OptimizeMoleculeConfsHelper_(ForceFields::ForceField ff, ROMol *mol,
}
}
inline void OptimizeMoleculeConfsMT(ROMol &mol, const ForceFields::ForceField &ff,
std::vector<std::pair<int, double>> &res,
int numThreads, int maxIters) {
inline void OptimizeMoleculeConfsMT(ROMol &mol,
const ForceFields::ForceField &ff,
std::vector<std::pair<int, double>> &res,
int numThreads, int maxIters) {
std::vector<std::thread> tg;
for (int ti = 0; ti < numThreads; ++ti) {
tg.emplace_back(std::thread(detail::OptimizeMoleculeConfsHelper_, ff, &mol,
@@ -60,8 +61,8 @@ inline void OptimizeMoleculeConfsMT(ROMol &mol, const ForceFields::ForceField &f
#endif
inline void OptimizeMoleculeConfsST(ROMol &mol, ForceFields::ForceField &ff,
std::vector<std::pair<int, double>> &res,
int maxIters) {
std::vector<std::pair<int, double>> &res,
int maxIters) {
PRECONDITION(res.size() >= mol.getNumConformers(),
"res.size() must be >= mol.getNumConformers()");
unsigned int i = 0;
@@ -91,7 +92,7 @@ inline void OptimizeMoleculeConfsST(ROMol &mol, ForceFields::ForceField &ff,
second: the energy
*/
inline std::pair<int, double> OptimizeMolecule(ForceFields::ForceField &ff,
int maxIters = 1000) {
int maxIters = 1000) {
ff.initialize();
int res = ff.minimize(maxIters);
double e = ff.calcEnergy();
@@ -112,8 +113,8 @@ inline std::pair<int, double> OptimizeMolecule(ForceFields::ForceField &ff,
*/
inline void OptimizeMoleculeConfs(ROMol &mol, ForceFields::ForceField &ff,
std::vector<std::pair<int, double>> &res,
int numThreads = 1, int maxIters = 1000) {
std::vector<std::pair<int, double>> &res,
int numThreads = 1, int maxIters = 1000) {
res.resize(mol.getNumConformers());
numThreads = getNumThreadsToUse(numThreads);
if (numThreads == 1) {
@@ -125,6 +126,23 @@ inline void OptimizeMoleculeConfs(ROMol &mol, ForceFields::ForceField &ff,
}
#endif
}
//! Convenience Function for generating an empty force Field with just
/// the molecules' atoms position.
/*
\param mol The molecule which positions should be added to the force field.
\param confId Id of the conformer which positions are used in the force field.
*/
inline std::unique_ptr<ForceFields::ForceField> createEmptyForceFieldForMol(
ROMol &mol, int confId = -1) {
auto res = std::make_unique<ForceFields::ForceField>();
auto &conf = mol.getConformer(confId);
for (auto &pt : conf.getPositions()) {
res->positions().push_back(&pt);
}
return res;
}
} // end of namespace ForceFieldsHelper
} // end of namespace RDKit
#endif

20
Code/GraphMol/ForceFieldHelpers/Wrap/rdForceFields.cpp
View File

@@ -86,6 +86,14 @@ python::object FFConfsHelper(ROMol &mol, ForceFields::PyForceField &ff,
return pyres;
}
ForceFields::PyForceField *CreateEmptyForceFieldForMol(ROMol &mol,
int confId = -1) {
auto ff = ForceFieldsHelper::createEmptyForceFieldForMol(mol, confId);
auto *res = new ForceFields::PyForceField(ff.release());
res->initialize();
return res;
}
ForceFields::PyForceField *UFFGetMoleculeForceField(
ROMol &mol, double vdwThresh = 10.0, int confId = -1,
bool ignoreInterfragInteractions = true) {
@@ -383,6 +391,18 @@ BOOST_PYTHON_MODULE(rdForceFieldHelpers) {
python::return_value_policy<python::manage_new_object>(),
docString.c_str());
docString =
"Get An empty Force Field, with only the positions of the atoms but no Contributions.\n\n\
\n\
ARGUMENTS :\n\n\
- mol : the molecule of interest\n\
- confId: the conformer which positions should be added to the force field.\n\
\n ";
python::def("CreateEmptyForceFieldForMol", RDKit::CreateEmptyForceFieldForMol,
(python::arg("mol"), python::arg("confId") = -1),
python::return_value_policy<python::manage_new_object>(),
docString.c_str());
docString =
"checks if MMFF parameters are available for all of a molecule's atoms\n\n\
\n\

31
Code/GraphMol/ForceFieldHelpers/Wrap/testHelpers.py
View File

@@ -378,6 +378,35 @@ M END"""
self.assertTrue(all(map(lambda i: i == 0, res)))
self.assertTrue(all(after[i] < b for i, b in enumerate(before)))
def testEmptyFF(self) -> None:
m = Chem.MolFromSmiles(
'CCCO |(-1.28533,-0.0567758,0.434662;-0.175447,0.695786,-0.299881;0.918409,-0.342619,-0.555572;1.30936,-0.801512,0.71705)|'
)
self.assertIsNotNone(m)
ff = ChemicalForceFields.CreateEmptyForceFieldForMol(m)
posa = m.GetConformer().GetAtomPosition(0)
posb = m.GetConformer().GetAtomPosition(1)
dist = (posa - posb).Length()
self.assertTrue(ff)
self.assertFalse(ff.Initialize())
self.assertEqual(ff.CalcEnergy(), 0.0)
self.assertTrue(all(v == 0.0 for v in ff.CalcGrad()))
self.assertEqual(ff.NumPoints(), m.GetNumAtoms())
self.assertEqual(len(ff.Positions()) / 3, m.GetNumAtoms())
self.assertFalse(ff.Minimize())
posa = m.GetConformer().GetAtomPosition(0)
posb = m.GetConformer().GetAtomPosition(1)
self.assertEqual((posa - posb).Length(), dist)
if __name__ == '__main__':
ff.MMFFAddDistanceConstraint(0, 1, False, 100, 100, 100)
self.assertFalse(ff.Minimize())
pos = ff.Positions()
dist = ((pos[0] - pos[3])**2 + (pos[1] - pos[4])**2 + (pos[2] - pos[5])**2)**0.5
self.assertAlmostEqual(dist, 100, delta=10e-5)
posa = m.GetConformer().GetAtomPosition(0)
posb = m.GetConformer().GetAtomPosition(1)
self.assertAlmostEqual((posa - posb).Length(), 100, delta=10e-5)
if __name__ == "__main__":
unittest.main()

54
Code/GraphMol/ForceFieldHelpers/catch_tests.cpp Normal file
View File

@@ -0,0 +1,54 @@
//
// Copyright (C) 2024 Niels Maeder and other RDKit contributors
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include <cmath>
#include <RDGeneral/test.h>
#include <catch2/catch_all.hpp>
#include <GraphMol/SmilesParse/SmilesParse.h>
#include <ForceField/MMFF/Params.h>
#include <ForceField/MMFF/BondStretch.h>
#include "FFConvenience.h"
using namespace RDKit;
TEST_CASE("Test empty force field") {
auto mol =
"CCCO |(-1.28533,-0.0567758,0.434662;-0.175447,0.695786,-0.299881;0.918409,-0.342619,-0.555572;1.30936,-0.801512,0.71705)|"_smiles;
REQUIRE(mol);
SECTION("basics") {
auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol);
REQUIRE(forceField);
forceField->initialize();
CHECK(forceField->minimize() == 0);
CHECK(forceField->calcEnergy() == 0.0);
REQUIRE(forceField->numPoints() == mol->getNumAtoms());
REQUIRE(forceField->positions().size() == mol->getNumAtoms());
auto dist = forceField->distance(0, 1);
auto pos1 = mol->getConformer().getAtomPos(0);
auto pos2 = mol->getConformer().getAtomPos(1);
CHECK(dist == (pos2 - pos1).length());
}
SECTION("add contrib and minimize") {
auto forceField = ForceFieldsHelper::createEmptyForceFieldForMol(*mol);
forceField->initialize();
auto params = ForceFields::MMFF::MMFFBond{6.0, 100.0};
auto *contrib = new ForceFields::MMFF::BondStretchContrib(forceField.get(),
0, 1, &params);
forceField->contribs().push_back(ForceFields::ContribPtr(contrib));
CHECK(forceField->minimize() == 0);
CHECK(std::round(forceField->distance(0, 1)) == 100);
auto pos1 = mol->getConformer().getAtomPos(0);
auto pos2 = mol->getConformer().getAtomPos(1);
auto dist = (pos2 - pos1).length();
CHECK(std::round(dist) == 100);
}
}