Add DCLV 3D descriptors

benchmarks/descriptors3d_bench.py

@@ -46,6 +46,8 @@
 # Absolute tolerances for USR skews (whole molecule, USRCAT atom classes): the cube root of a near-zero third
 # moment turns rounding noise into visible values, in RDKit as well (a two-atom class has an exact skew of 0).
 USR_SKEW_TOLERANCES = {PrecisionMode.SINGLE: (5e-3, 0.1), PrecisionMode.FULL: (1e-4, 1e-3)}
+# Relative tolerances for the DCLV properties, which sum float32 dot contributions over every atom.
+DCLV_RELATIVE_TOLERANCES = {PrecisionMode.SINGLE: 1e-5, PrecisionMode.FULL: 1e-9}
 
 PROPERTY_SETS = {
     "single": (Property3D.RADIUS_OF_GYRATION,),
@@ -64,6 +66,16 @@
     Property3D.ASPHERICITY: rdMolDescriptors.CalcAsphericity,
 }
 
+DCLV_GETTERS = {
+    Property3D.DCLV_SURFACE_AREA: "GetSurfaceArea",
+    Property3D.DCLV_POLAR_SURFACE_AREA: "GetPolarSurfaceArea",
+    Property3D.DCLV_VOLUME: "GetVolume",
+    Property3D.DCLV_VDW_VOLUME: "GetVDWVolume",
+    Property3D.DCLV_POLAR_VOLUME: "GetPolarVolume",
+    Property3D.DCLV_COMPACTNESS: "GetCompactness",
+    Property3D.DCLV_PACKING_DENSITY: "GetPackingDensity",
+}
+
 
 def _calc_rdkit_property(mol: Chem.Mol, conf_id: int, prop: Property3D) -> float | list[float]:
     """Calculate one property using the options mirrored by the benchmark."""
@@ -127,17 +139,36 @@ def _pack_device_coordinates(mols: list[Chem.Mol]) -> Device3DResult:
     )
 
 
-def _calc_rdkit(mols: list[Chem.Mol], properties: tuple[Property3D, ...]) -> dict[Property3D, np.ndarray]:
-    """Calculate reference arrays in molecule/conformer order."""
+def _dclv_radii(mol: Chem.Mol) -> list[float]:
+    """Van der Waals radii per atom; RDKit 2025.09.1 has no default-radii DCLV constructor."""
+    table = Chem.GetPeriodicTable()
+    return [table.GetRvdw(atom.GetAtomicNum()) for atom in mol.GetAtoms()]
+
+
+def _calc_rdkit_conformer(
+    mol: Chem.Mol, conf_id: int, properties: tuple[Property3D, ...]
+) -> dict[Property3D, float | list[float]]:
+    dclv = (
+        rdMolDescriptors.DoubleCubicLatticeVolume(mol, _dclv_radii(mol), confId=conf_id)
+        if any(prop in DCLV_GETTERS for prop in properties)
+        else None
+    )
     return {
-        prop: np.asarray(
-            [_calc_rdkit_property(mol, conf.GetId(), prop) for mol in mols for conf in mol.GetConformers()],
-            dtype=np.float64,
-        )
+        prop: getattr(dclv, DCLV_GETTERS[prop])() if prop in DCLV_GETTERS else _calc_rdkit_property(mol, conf_id, prop)
         for prop in properties
     }
 
 
+def _calc_rdkit(mols: list[Chem.Mol], properties: tuple[Property3D, ...]) -> dict[Property3D, np.ndarray]:
+    """Calculate reference arrays in molecule/conformer order."""
+    values = {prop: [] for prop in properties}
+    for mol in mols:
+        for conf in mol.GetConformers():
+            for prop, value in _calc_rdkit_conformer(mol, conf.GetId(), properties).items():
+                values[prop].append(value)
+    return {prop: np.asarray(rows, dtype=np.float64) for prop, rows in values.items()}
+
+
 def _validate(
     mols: list[Chem.Mol],
     properties: tuple[Property3D, ...],
@@ -165,6 +196,9 @@ def _validate(
                 computed = ~np.isnan(expected[prop]).all(axis=1)
                 np.testing.assert_allclose(actual[computed], expected[prop][computed], rtol=0, atol=max(atol, 1.1e-3))
                 continue
+            if prop in DCLV_GETTERS:
+                np.testing.assert_allclose(actual, expected[prop], rtol=DCLV_RELATIVE_TOLERANCES[precision], atol=atol)
+                continue
             if prop not in (Property3D.USR, Property3D.USRCAT):
                 np.testing.assert_allclose(actual, expected[prop], rtol=rtol, atol=atol)
                 continue
@@ -196,8 +230,7 @@ def _rdkit_worker_range(bounds: tuple[int, int]) -> int:
     done = 0
     for mol in _WORKER_MOLS[start:stop]:
         for conf in mol.GetConformers():
-            for prop in _WORKER_PROPERTIES:
-                _calc_rdkit_property(mol, conf.GetId(), prop)
+            _calc_rdkit_conformer(mol, conf.GetId(), _WORKER_PROPERTIES)
             done += 1
     return done
 

benchmarks/tests/test_descriptors3d_dclv_bench.py

@@ -0,0 +1,161 @@
+# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+# SPDX-License-Identifier: Apache-2.0
+
+import descriptors3d_bench as bench
+import numpy as np
+import pytest
+from rdkit import Chem
+
+
+@pytest.fixture
+def molecules():
+    result = []
+    for smiles, conformer_ids in [("CCO", [7, 11]), ("CC", []), ("CO", [23])]:
+        molecule = Chem.MolFromSmiles(smiles)
+        for conformer_id in conformer_ids:
+            conformer = Chem.Conformer(molecule.GetNumAtoms())
+            conformer.SetId(conformer_id)
+            for atom in range(molecule.GetNumAtoms()):
+                conformer.SetAtomPosition(atom, (1.2 * atom, 0.2 * atom * atom, 0.01 * conformer_id * atom))
+            molecule.AddConformer(conformer, assignId=False)
+        result.append(molecule)
+    return result
+
+
+@pytest.mark.parametrize(
+    "properties",
+    [
+        (bench.Property3D.DCLV_VOLUME,),
+        tuple(bench.DCLV_GETTERS),
+        tuple(bench.Property3D),
+        (bench.Property3D.DCLV_VOLUME, bench.Property3D.PMI1, bench.Property3D.DCLV_VOLUME),
+        (
+            bench.Property3D.DCLV_VOLUME,
+            bench.Property3D.PMI1,
+            bench.Property3D.DCLV_SURFACE_AREA,
+            bench.Property3D.USR,
+        ),
+    ],
+)
+def test_dclv_reference_reuses_one_construction_per_conformer_and_timed_call(monkeypatch, molecules, properties):
+    constructor = bench.rdMolDescriptors.DoubleCubicLatticeVolume
+    expected = {prop: [] for prop in properties}
+    for molecule in molecules:
+        for conformer in molecule.GetConformers():
+            for prop in expected:
+                value = (
+                    getattr(
+                        constructor(molecule, bench._dclv_radii(molecule), confId=conformer.GetId()),
+                        bench.DCLV_GETTERS[prop],
+                    )()
+                    if prop in bench.DCLV_GETTERS
+                    else bench._calc_rdkit_property(molecule, conformer.GetId(), prop)
+                )
+                expected[prop].append(value)
+    constructions = []
+
+    def counted_constructor(molecule, radii, **kwargs):
+        assert radii == bench._dclv_radii(molecule)
+        constructions.append((molecule, kwargs["confId"]))
+        return constructor(molecule, radii, **kwargs)
+
+    monkeypatch.setattr(bench.rdMolDescriptors, "DoubleCubicLatticeVolume", counted_constructor)
+    expected_constructions = [
+        (molecule, conformer.GetId()) for molecule in molecules for conformer in molecule.GetConformers()
+    ]
+    for repeat in range(2):
+        actual = bench._calc_rdkit(molecules, properties)
+        assert constructions == expected_constructions * (repeat + 1)
+        assert tuple(actual) == tuple(expected)
+        for prop in expected:
+            assert actual[prop].dtype == np.float64
+            np.testing.assert_allclose(actual[prop], expected[prop], rtol=0, atol=0, equal_nan=True)
+
+
+@pytest.mark.parametrize("properties", [(), (bench.Property3D.PMI1,), (bench.Property3D.USR,)])
+def test_reference_without_dclv_does_not_construct_it(monkeypatch, molecules, properties):
+    def unexpected_constructor(*args, **kwargs):
+        pytest.fail("DCLV constructed when no DCLV property was requested")
+
+    monkeypatch.setattr(bench.rdMolDescriptors, "DoubleCubicLatticeVolume", unexpected_constructor)
+    actual = bench._calc_rdkit(molecules, properties)
+    assert tuple(actual) == properties
+    for prop in properties:
+        expected = [
+            bench._calc_rdkit_property(molecule, conformer.GetId(), prop)
+            for molecule in molecules
+            for conformer in molecule.GetConformers()
+        ]
+        np.testing.assert_allclose(actual[prop], expected, rtol=0, atol=0, equal_nan=True)
+
+
+@pytest.mark.parametrize("empty_molecules", [[], [Chem.MolFromSmiles("CCO")]])
+def test_empty_reference_has_empty_float64_columns(monkeypatch, empty_molecules):
+    def unexpected_constructor(*args, **kwargs):
+        pytest.fail("DCLV constructed without a conformer")
+
+    monkeypatch.setattr(bench.rdMolDescriptors, "DoubleCubicLatticeVolume", unexpected_constructor)
+    properties = tuple(bench.DCLV_GETTERS) + (bench.Property3D.PMI1,)
+    actual = bench._calc_rdkit(empty_molecules, properties)
+    assert tuple(actual) == properties
+    assert all(values.shape == (0,) and values.dtype == np.float64 for values in actual.values())
+
+
+def test_dclv_reference_propagates_constructor_errors(monkeypatch, molecules):
+    def failed_constructor(*args, **kwargs):
+        raise ValueError("invalid reference geometry")
+
+    monkeypatch.setattr(bench.rdMolDescriptors, "DoubleCubicLatticeVolume", failed_constructor)
+    with pytest.raises(ValueError, match="invalid reference geometry"):
+        bench._calc_rdkit(molecules, (bench.Property3D.DCLV_VOLUME,))
+
+
+@pytest.mark.parametrize("bounds,expected_ids", [((0, 3), [7, 11, 23]), ((1, 3), [23]), ((0, 0), [])])
+@pytest.mark.parametrize(
+    "properties",
+    [
+        tuple(bench.DCLV_GETTERS),
+        (bench.Property3D.PMI1,),
+        (bench.Property3D.DCLV_VOLUME, bench.Property3D.PMI1, bench.Property3D.USR),
+    ],
+)
+def test_timed_worker_constructs_dclv_once_and_reads_each_requested_getter(
+    monkeypatch, molecules, bounds, expected_ids, properties
+):
+    constructor = bench.rdMolDescriptors.DoubleCubicLatticeVolume
+    calculator = bench._calc_rdkit_property
+    constructions = []
+    getters = []
+    other_properties = []
+
+    def counted_property(molecule, conf_id, prop):
+        other_properties.append((conf_id, prop))
+        return calculator(molecule, conf_id, prop)
+
+    class CountedDclv:
+        def __init__(self, molecule, radii, *, confId):
+            constructions.append(confId)
+            self.conf_id = confId
+            self.reference = constructor(molecule, radii, confId=confId)
+
+        def __getattr__(self, name):
+            def read_property():
+                getters.append((self.conf_id, name))
+                return getattr(self.reference, name)()
+
+            return read_property
+
+    monkeypatch.setattr(bench.rdMolDescriptors, "DoubleCubicLatticeVolume", CountedDclv)
+    monkeypatch.setattr(bench, "_calc_rdkit_property", counted_property)
+    monkeypatch.setattr(bench, "_WORKER_MOLS", molecules)
+    monkeypatch.setattr(bench, "_WORKER_PROPERTIES", properties)
+    dclv_properties = [prop for prop in properties if prop in bench.DCLV_GETTERS]
+    for repeat in range(2):
+        assert bench._rdkit_worker_range(bounds) == len(expected_ids)
+        assert constructions == (expected_ids * (repeat + 1) if dclv_properties else [])
+        assert getters == [
+            (conf_id, bench.DCLV_GETTERS[prop]) for conf_id in expected_ids for prop in dclv_properties
+        ] * (repeat + 1)
+        assert other_properties == [
+            (conf_id, prop) for conf_id in expected_ids for prop in properties if prop not in bench.DCLV_GETTERS
+        ] * (repeat + 1)

docs/api/nvmolkit.rst

@@ -178,6 +178,7 @@ Conformer RMSD
    descriptors3d.MomentOptions
    descriptors3d.WhimOptions
    descriptors3d.GetawayOptions
+   descriptors3d.DclvOptions
 
 Torsion Fingerprint Deviation (TFD)
 -----------------------------------

nvmolkit/descriptors3d.cpp

@@ -90,6 +90,9 @@ bp::object calc3DProperties(const bp::list&                mols,
                             const bool                     useAtomicMasses,
                             const double                   whimThreshold,
                             const unsigned int             getawayPrecision,
+                            const double                   dclvProbeRadius,
+                            const bool                     dclvIncludeSandP,
+                            const bool                     dclvIncludeHs,
                             const bp::object&              coordinates,
                             const nvMolKit::PrecisionMode& precision,
                             const int                      preprocessingThreads,
@@ -102,6 +105,9 @@ bp::object calc3DProperties(const bp::list&                mols,
   options.moments.useAtomicMasses = useAtomicMasses;
   options.whim.threshold          = whimThreshold;
   options.getaway.precision       = getawayPrecision;
+  options.dclv.probeRadius        = dclvProbeRadius;
+  options.dclv.includeSandP       = dclvIncludeSandP;
+  options.dclv.includeHs          = dclvIncludeHs;
 
   std::vector<nvMolKit::Property3D> properties;
   for (int i = 0; i < bp::len(propertyNames); ++i) {
@@ -134,6 +140,9 @@ BOOST_PYTHON_MODULE(_descriptors3d) {
            bp::arg("useAtomicMasses"),
            bp::arg("whimThreshold"),
            bp::arg("getawayPrecision"),
+           bp::arg("dclvProbeRadius"),
+           bp::arg("dclvIncludeSandP"),
+           bp::arg("dclvIncludeHs"),
            bp::arg("coordinates"),
            bp::arg("precision"),
            bp::arg("preprocessingThreads"),

nvmolkit/descriptors3d.py

@@ -27,6 +27,12 @@ class Property3D(Enum):
     ``rdMolDescriptors`` module (``CalcWHIM``, ``GetUSR``, ...). Wherever a
     property is accepted, its string value may be used instead.
 
+    The ``DCLV`` properties correspond to the getter methods of RDKit's
+    ``DoubleCubicLatticeVolume`` class (for example, ``DCLVSurfaceArea``
+    corresponds to ``GetSurfaceArea``), configured by :class:`DclvOptions`. Requesting
+    several costs one pass over the surface dots. Conformers containing an atom
+    with a zero van der Waals radius produce NaN.
+
     Vector properties, in RDKit's order; every other property is a scalar:
 
     - ``WHIM``: 114 values.
@@ -61,6 +67,13 @@ class Property3D(Enum):
     USR = "USR"
     USRCAT = "USRCAT"
     GETAWAY = "GETAWAY"
+    DCLV_SURFACE_AREA = "DCLVSurfaceArea"
+    DCLV_POLAR_SURFACE_AREA = "DCLVPolarSurfaceArea"
+    DCLV_VOLUME = "DCLVVolume"
+    DCLV_VDW_VOLUME = "DCLVVDWVolume"
+    DCLV_POLAR_VOLUME = "DCLVPolarVolume"
+    DCLV_COMPACTNESS = "DCLVCompactness"
+    DCLV_PACKING_DENSITY = "DCLVPackingDensity"
 
 
 @dataclass(frozen=True)
@@ -99,6 +112,27 @@ class GetawayOptions:
     precision: int = 2
 
 
+@dataclass(frozen=True)
+class DclvOptions:
+    """Options for the ``DCLV`` properties (RDKit's ``DoubleCubicLatticeVolume``), with RDKit's defaults.
+
+    Atoms take RDKit's default van der Waals radii. RDKit's ``radii``, ``isProtein`` and
+    ``includeLigand`` arguments are not supported (``isProtein=False``).
+
+    Attributes:
+        probeRadius: Radius of the solvent probe sphere in Angstrom, RDKit's default ``1.4``. Must be
+            finite and non-negative.
+        includeSandP: ``DCLVPolarSurfaceArea`` and ``DCLVPolarVolume`` also count sulfur and phosphorus, as
+            RDKit's ``includeSandP`` argument of ``GetPolarSurfaceArea`` and ``GetPolarVolume``.
+        includeHs: The polar values also count hydrogens bonded to a polar atom, as RDKit's
+            ``includeHs`` argument.
+    """
+
+    probeRadius: float = 1.4
+    includeSandP: bool = False
+    includeHs: bool = False
+
+
 @dataclass(frozen=True)
 class Property3DOptions:
     """Per-family options for :func:`Calc3DProperties`; each family reads only its own member.
@@ -121,11 +155,13 @@ class Property3DOptions:
         moments: Options for the moment-based properties.
         whim: Options for ``WHIM``.
         getaway: Options for ``GETAWAY``.
+        dclv: Options for the ``DCLV`` properties.
     """
 
     moments: MomentOptions = MomentOptions()
     whim: WhimOptions = WhimOptions()
     getaway: GetawayOptions = GetawayOptions()
+    dclv: DclvOptions = DclvOptions()
 
 
 @dataclass(frozen=True)
@@ -227,6 +263,9 @@ def dense(self, pad_value: float = float("nan")) -> Dense3DPropertyResult:
         return Dense3DPropertyResult(values=values, conf_mask=conf_mask)
 
 
+_DCLV_PROPERTIES = frozenset(prop for prop in Property3D if prop.value.startswith("DCLV"))
+
+
 def _normalize_property(property_name: Property3D | str) -> Property3D:
     if isinstance(property_name, Property3D):
         return property_name
@@ -336,6 +375,9 @@ def Calc3DProperties(
             MORSE's first scattering value (a large sum of positive terms),
             USR's reference-atom choice (near-equidistant atoms) and GETAWAY's
             leverages (rounded and clustered for ITH and ISH).
+            The ``DCLV`` properties count exposed surface dots, and a dot
+            within rounding of a neighboring sphere can flip between precisions
+            (or against RDKit), which moves the values by one dot's share.
         hardwareOptions: Only ``preprocessingThreads`` applies: the CPU
             threads used to extract coordinates and atom weights from the
             molecules (default ``-1``, all threads).
@@ -406,6 +448,8 @@ def Calc3DProperties(
             raise TypeError(f"GETAWAY precision must be an int, got {type(getaway_precision).__name__}")
         if not 1 <= getaway_precision <= 6:
             raise ValueError(f"GETAWAY precision must be between 1 and 6 significant digits, got {getaway_precision}")
+    # Options of unrequested families are ignored, even when invalid; the native code validates the probe radius.
+    dclv = options.dclv if any(prop in _DCLV_PROPERTIES for prop in normalized_properties) else DclvOptions()
     if hardwareOptions is None:
         hardwareOptions = HardwareOptions()
     elif not isinstance(hardwareOptions, HardwareOptions):
@@ -426,6 +470,9 @@ def Calc3DProperties(
         options.moments.useAtomicMasses,
         options.whim.threshold,
         getaway_precision,
+        float(dclv.probeRadius),
+        bool(dclv.includeSandP),
+        bool(dclv.includeHs),
         coordinate_interfaces,
         precision,
         hardwareOptions.preprocessingThreads,