allow default radii in the DCLV calculation (#8836)

This commit is contained in:
Greg Landrum
2025-10-11 20:35:20 +02:00
committed by GitHub
parent 9fd5595f98
commit 3a4d4df17e
5 changed files with 92 additions and 57 deletions

View File

@@ -252,6 +252,14 @@ DoubleCubicLatticeVolume::DoubleCubicLatticeVolume(
object
*/
if (radii_.empty()) {
const auto *tbl = PeriodicTable::getTable();
radii_.reserve(mol.getNumAtoms());
for (const auto atom : mol.atoms()) {
radii_.push_back(tbl->getRvdw(atom->getAtomicNum()));
}
}
positions = mol.getConformer(confId).getPositions();
maxRadius = *std::max_element(radii_.begin(), radii_.end());

View File

@@ -26,12 +26,29 @@
namespace RDKit {
namespace Descriptors {
//! Class for calculation of the Shrake and Rupley surface area and volume
//! using the Double Cubic Lattice Method.
//!
//! Frank Eisenhaber, Philip Lijnzaad, Patrick Argos, Chris Sander and
//! Michael Scharf, "The Double Cubic Lattice Method: Efficient Approaches
//! to Numerical Integration of Surface Area and Volume and to Dot Surface
//! Contouring of Molecular Assemblies", Journal of Computational Chemistry,
//! Vol. 16, No. 3, pp. 273-284, 1995.
class RDKIT_DESCRIPTORS_EXPORT DoubleCubicLatticeVolume {
public:
// default params assume a small molecule and default conformer
const ROMol &mol;
std::vector<double> radii_;
bool isProtein = false;
bool includeLigand = true;
double probeRadius = 1.4;
int confId = -1;
double maxRadius = 1.7; // treat default max radius as Carbon
/*!
\param mol: input molecule or protein
\param radii: radii for atoms of input mol
\param radii: radii for atoms of input mol, empty for default values
\param isProtein: flag to calculate burried surface area of a protein ligand
complex [default=false, free ligand]
\param includeLigand: flag to trigger
@@ -42,27 +59,16 @@ class RDKIT_DESCRIPTORS_EXPORT DoubleCubicLatticeVolume {
\param confId: conformer ID to consider [default=-1]
*/
// default params assume a small molecule and default conformer
const ROMol &mol;
std::vector<double> radii_;
bool isProtein = false;
bool includeLigand = true;
double probeRadius = 1.4;
int confId = -1;
double maxRadius = 1.7; // treat default max radius as Carbon
DoubleCubicLatticeVolume(const ROMol &mol, std::vector<double> radii,
bool isProtein = false, bool includeLigand = true,
double probeRadius = 1.4, int confId = -1);
//! Class for calculation of the Shrake and Rupley surface area and volume
//! using the Double Cubic Lattice Method.
//!
//! Frank Eisenhaber, Philip Lijnzaad, Patrick Argos, Chris Sander and
//! Michael Scharf, "The Double Cubic Lattice Method: Efficient Approaches
//! to Numerical Integration of Surface Area and Volume and to Dot Surface
//! Contouring of Molecular Assemblies", Journal of Computational Chemistry,
//! Vol. 16, No. 3, pp. 273-284, 1995.
//! \overload uses default vdw radii
DoubleCubicLatticeVolume(const ROMol &mol, bool isProtein = false,
bool includeLigand = true, double probeRadius = 1.4,
int confId = -1)
: DoubleCubicLatticeVolume(mol, std::vector<double>(), isProtein,
includeLigand, probeRadius, confId) {};
// value returns

View File

@@ -1,5 +1,5 @@
//
// Copyright (C) 2007-2017 Greg Landrum
// Copyright (C) 2007-2025 Greg Landrum and other RDKit contributors
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
@@ -1712,7 +1712,11 @@ BOOST_PYTHON_MODULE(rdMolDescriptors) {
RDKit::Descriptors::DoubleCubicLatticeVolume,
boost::shared_ptr<RDKit::Descriptors::DoubleCubicLatticeVolume>>(
"DoubleCubicLatticeVolume",
"Class for the Double Cubic Lattice Volume method", python::no_init)
"Class for the Double Cubic Lattice Volume method",
python::init<const RDKit::ROMol &, bool, bool, double, int>(
(python::arg("mol"), python::arg("isProtein") = false,
python::arg("includeLigand") = true,
python::arg("probeRadius") = 1.4, python::arg("confId") = -1)))
.def("__init__",
python::make_constructor(
&getDoubleCubicLatticeVolume, python::default_call_policies(),

View File

@@ -775,14 +775,15 @@ class TestCase(unittest.TestCase):
radii2 = [Chem.GetPeriodicTable().GetRvdw(atom.GetAtomicNum()) for atom in mol2.GetAtoms()]
# test from SDF file with defaults
sdf = rdMD.DoubleCubicLatticeVolume(mol2, radii2, isProtein=False)
for sdf in (rdMD.DoubleCubicLatticeVolume(mol2, radii2, isProtein=False),
rdMD.DoubleCubicLatticeVolume(mol2, isProtein=False)):
self.assertTrue(abs(sdf.GetSurfaceArea() - 304.239) < 0.05)
self.assertTrue(abs(sdf.GetPolarSurfaceArea() - 18.7319) < 0.05)
self.assertTrue(abs(sdf.GetPolarSurfaceArea(includeSandP=True) - 64.9764) < 0.05)
self.assertTrue(abs(sdf.GetVolume() - 431.35) < 0.05)
self.assertTrue(abs(sdf.GetVDWVolume() - 119.296) < 0.05)
self.assertTrue(abs(sdf.GetPolarVolume() - 21.35) < 0.05)
self.assertTrue(abs(sdf.GetSurfaceArea() - 304.239) < 0.05)
self.assertTrue(abs(sdf.GetPolarSurfaceArea() - 18.7319) < 0.05)
self.assertTrue(abs(sdf.GetPolarSurfaceArea(includeSandP=True) - 64.9764) < 0.05)
self.assertTrue(abs(sdf.GetVolume() - 431.35) < 0.05)
self.assertTrue(abs(sdf.GetVDWVolume() - 119.296) < 0.05)
self.assertTrue(abs(sdf.GetPolarVolume() - 21.35) < 0.05)
if __name__ == '__main__':

View File

@@ -565,28 +565,40 @@ TEST_CASE("DCLV") {
const bool includeLigand = false;
const double probeRadius = 1.4;
const int confId = -1;
Descriptors::DoubleCubicLatticeVolume dclv(*m, radii, isProtein, includeLigand, probeRadius, confId);
CHECK(dclv.getSurfaceArea() == Catch::Approx(8306.62).epsilon(0.05));
CHECK(dclv.getPolarSurfaceArea(false, false) == Catch::Approx(4652.19).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) == Catch::Approx(4673.9).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getVolume() == Catch::Approx(29952.3).epsilon(0.05));
CHECK(dclv.getVDWVolume() == Catch::Approx(13541.5).epsilon(0.05));
CHECK(dclv.getPolarVolume(false, false) == Catch::Approx(17096).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getCompactness() == Catch::Approx(1.78096).epsilon(0.05));
CHECK(dclv.getPackingDensity() == Catch::Approx(0.452103).epsilon(0.05));
for (auto dclv :
{Descriptors::DoubleCubicLatticeVolume(
*m, radii, isProtein, includeLigand, probeRadius, confId),
Descriptors::DoubleCubicLatticeVolume(*m, isProtein, includeLigand,
probeRadius, confId)}) {
CHECK(dclv.getSurfaceArea() == Catch::Approx(8306.62).epsilon(0.05));
CHECK(dclv.getPolarSurfaceArea(false, false) ==
Catch::Approx(4652.19).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) ==
Catch::Approx(4673.9).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getVolume() == Catch::Approx(29952.3).epsilon(0.05));
CHECK(dclv.getVDWVolume() == Catch::Approx(13541.5).epsilon(0.05));
CHECK(dclv.getPolarVolume(false, false) ==
Catch::Approx(17096).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getCompactness() == Catch::Approx(1.78096).epsilon(0.05));
CHECK(dclv.getPackingDensity() == Catch::Approx(0.452103).epsilon(0.05));
}
}
SECTION("non-default radius") {
const bool isProtein = true;
const bool includeLigand = false;
const double probeRadius = 1.6;
const int confId = -1;
Descriptors::DoubleCubicLatticeVolume dclv(*m, radii, isProtein, includeLigand, probeRadius, confId);
Descriptors::DoubleCubicLatticeVolume dclv(
*m, radii, isProtein, includeLigand, probeRadius, confId);
CHECK(dclv.getSurfaceArea() == Catch::Approx(8112.42).epsilon(0.05));
CHECK(dclv.getPolarSurfaceArea(false, false) == Catch::Approx(4656.8).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) == Catch::Approx(4674.96).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getPolarSurfaceArea(false, false) ==
Catch::Approx(4656.8).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) ==
Catch::Approx(4674.96).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getVolume() == Catch::Approx(31591).epsilon(0.05));
CHECK(dclv.getVDWVolume() == Catch::Approx(13541.5).epsilon(0.05));
CHECK(dclv.getPolarVolume(false, false) == Catch::Approx(18402.8).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarVolume(false, false) ==
Catch::Approx(18402.8).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getCompactness() == Catch::Approx(1.67864).epsilon(0.05));
CHECK(dclv.getPackingDensity() == Catch::Approx(0.428652).epsilon(0.05));
}
@@ -595,13 +607,17 @@ TEST_CASE("DCLV") {
const bool includeLigand = true;
const double probeRadius = 1.4;
const int confId = -1;
Descriptors::DoubleCubicLatticeVolume dclv(*m, radii, isProtein, includeLigand, probeRadius, confId);
Descriptors::DoubleCubicLatticeVolume dclv(
*m, radii, isProtein, includeLigand, probeRadius, confId);
CHECK(dclv.getSurfaceArea() == Catch::Approx(8206.1).epsilon(0.05));
CHECK(dclv.getPolarSurfaceArea(false, false) == Catch::Approx(4467.92).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) == Catch::Approx(4481.19).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getPolarSurfaceArea(false, false) ==
Catch::Approx(4467.92).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) ==
Catch::Approx(4481.19).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getVolume() == Catch::Approx(30340.1).epsilon(0.05));
CHECK(dclv.getVDWVolume() == Catch::Approx(13789.7).epsilon(0.05));
CHECK(dclv.getPolarVolume(false, false) == Catch::Approx(16507).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarVolume(false, false) ==
Catch::Approx(16507).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getCompactness() == Catch::Approx(1.74438).epsilon(0.05));
CHECK(dclv.getPackingDensity() == Catch::Approx(0.454504).epsilon(0.05));
}
@@ -614,19 +630,18 @@ TEST_CASE("DCLV") {
radii.push_back(tbl->getRvdw(atom->getAtomicNum()));
}
REQUIRE(m);
const bool isProtein = false;
const bool includeLigand = true; // needs to be true for small mol
const double probeRadius = 1.4;
const int confId = -1;
Descriptors::DoubleCubicLatticeVolume dclv(*m, radii, isProtein, includeLigand, probeRadius, confId);
Descriptors::DoubleCubicLatticeVolume dclv(*m, radii);
// packing density and compactness not relevant for small mol
CHECK(dclv.getSurfaceArea() == Catch::Approx(304.239).epsilon(0.05));
CHECK(dclv.getPolarSurfaceArea(false, false) == Catch::Approx(18.7319).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) == Catch::Approx(64.9764).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getPolarSurfaceArea(false, false) ==
Catch::Approx(18.7319).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarSurfaceArea(true, false) ==
Catch::Approx(64.9764).epsilon(0.05)); // N, O, S & P, no Hs
CHECK(dclv.getVolume() == Catch::Approx(431.35).epsilon(0.05));
CHECK(dclv.getVDWVolume() == Catch::Approx(119.296).epsilon(0.05));
CHECK(dclv.getPolarVolume(false, false) == Catch::Approx(21.35).epsilon(0.05)); // N & O, no Hs
CHECK(dclv.getPolarVolume(false, false) ==
Catch::Approx(21.35).epsilon(0.05)); // N & O, no Hs
}
SECTION("Partial Test") {
std::string sdfName =
@@ -637,10 +652,11 @@ TEST_CASE("DCLV") {
std::vector<double> radii = {1.7, 1.7};
const bool isProtein = false;
const bool includeLigand = true; // needs to be true for small mol
const bool includeLigand = true; // needs to be true for small mol
const double probeRadius = 1.4;
const int confId = -1;
Descriptors::DoubleCubicLatticeVolume dclv(*m, radii, isProtein, includeLigand, probeRadius, confId);
Descriptors::DoubleCubicLatticeVolume dclv(
*m, radii, isProtein, includeLigand, probeRadius, confId);
boost::dynamic_bitset<> partialAtoms(m->getNumAtoms());
partialAtoms.set(0);
@@ -653,9 +669,9 @@ TEST_CASE("DCLV") {
partialAtoms3.set(1);
CHECK(dclv.getAtomSurfaceArea(0) == dclv.getAtomSurfaceArea(1));
CHECK(dclv.getPartialSurfaceArea(partialAtoms) == dclv.getPartialSurfaceArea(partialAtoms2));
CHECK(dclv.getPartialSurfaceArea(partialAtoms) ==
dclv.getPartialSurfaceArea(partialAtoms2));
CHECK(dclv.getPartialSurfaceArea(partialAtoms3) == dclv.getSurfaceArea());
}
}