From 45041446c5e415a67edfd986caa1c78003de8e63 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Thu, 8 Oct 2026 14:33:44 -0400 Subject: [PATCH 1/5] Add DCLV surface area and volume 3D descriptors Port RDKit's DoubleCubicLatticeVolume (Double Cubic Lattice Method surface area and volume) as seven scalar properties in a new Dclv family, one per getter: DCLVSurfaceArea, DCLVPolarSurfaceArea, DCLVVolume, DCLVVDWVolume, DCLVPolarVolume, DCLVCompactness and DCLVPackingDensity. Any selection shares one kernel launch. One warp per conformer centers the coordinates on RDKit's centre of gravity in float64, then per atom gathers neighbors in shared memory (falling back to a full scan past 128) and tests RDKit's 320 dots. The probe-expanded dots serve every property but DCLVVDWVolume and the bare van der Waals dots serve DCLVVDWVolume and DCLVPackingDensity; each set is tested only when a requested property needs it. DclvOptions exposes probeRadius, includeSandP and includeHs. Custom radii, isProtein and includeLigand are not supported. Conformers with a dummy atom (radius 0), on which RDKit crashes, and rows with non-finite coordinates produce NaN. On 484 ChEMBL conformers FULL matches RDKit to 3e-15 relative and SINGLE to 1.2e-6 relative with no surface-dot flips, so the kernel computes in the requested precision apart from the existing float64 centering. --- benchmarks/descriptors3d_bench.py | 18 + docs/api/nvmolkit.rst | 1 + nvmolkit/descriptors3d.cpp | 9 + nvmolkit/descriptors3d.py | 48 ++ nvmolkit/tests/test_descriptors3d.py | 139 ++++- src/descriptors3d.cu | 37 ++ src/descriptors3d.h | 91 +++- src/descriptors3d_dclv.cuh | 743 +++++++++++++++++++++++++++ src/descriptors3d_mol.cpp | 57 +- tests/test_descriptors3d.cpp | 80 ++- 10 files changed, 1199 insertions(+), 24 deletions(-) create mode 100644 src/descriptors3d_dclv.cuh diff --git a/benchmarks/descriptors3d_bench.py b/benchmarks/descriptors3d_bench.py index d673d6e1..baad2099 100644 --- a/benchmarks/descriptors3d_bench.py +++ b/benchmarks/descriptors3d_bench.py @@ -42,6 +42,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,), @@ -60,6 +62,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.""" @@ -89,6 +101,9 @@ def _calc_rdkit_property(mol: Chem.Mol, conf_id: int, prop: Property3D) -> float # RDKit's GETAWAY does not finish in reasonable time on multi-fragment molecules. return [math.nan] * 273 return rdMolDescriptors.CalcGETAWAY(mol, confId=conf_id) + if prop in DCLV_GETTERS: + dclv = rdMolDescriptors.DoubleCubicLatticeVolume(mol, confId=conf_id) + return getattr(dclv, DCLV_GETTERS[prop])() return RDKIT_CALCULATORS[prop](mol, confId=conf_id, useAtomicMasses=True) @@ -161,6 +176,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 diff --git a/docs/api/nvmolkit.rst b/docs/api/nvmolkit.rst index 107a7c8c..d31d6519 100644 --- a/docs/api/nvmolkit.rst +++ b/docs/api/nvmolkit.rst @@ -178,6 +178,7 @@ Conformer RMSD descriptors3d.MomentOptions descriptors3d.WhimOptions descriptors3d.GetawayOptions + descriptors3d.DclvOptions Torsion Fingerprint Deviation (TFD) ----------------------------------- diff --git a/nvmolkit/descriptors3d.cpp b/nvmolkit/descriptors3d.cpp index f1cf04dd..ba6e0480 100644 --- a/nvmolkit/descriptors3d.cpp +++ b/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 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"), diff --git a/nvmolkit/descriptors3d.py b/nvmolkit/descriptors3d.py index 235201e4..f3362f62 100644 --- a/nvmolkit/descriptors3d.py +++ b/nvmolkit/descriptors3d.py @@ -27,6 +27,13 @@ class Property3D(Enum): ``rdMolDescriptors`` module (``CalcWHIM``, ``GetUSR``, ...). Wherever a property is accepted, its string value may be used instead. + The ``DCLV`` properties are the getters of RDKit's + ``DoubleCubicLatticeVolume`` class (``DCLVSurfaceArea`` is + ``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 (dummy atoms, on which RDKit crashes) + produce NaN. + Vector properties, in RDKit's order; every other property is a scalar: - ``WHIM``: 114 values. @@ -61,6 +68,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 +113,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 +156,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 +264,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 +376,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 +449,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 +471,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, diff --git a/nvmolkit/tests/test_descriptors3d.py b/nvmolkit/tests/test_descriptors3d.py index f9f6b79b..cea1f792 100644 --- a/nvmolkit/tests/test_descriptors3d.py +++ b/nvmolkit/tests/test_descriptors3d.py @@ -13,6 +13,7 @@ from nvmolkit.descriptors3d import ( Calc3DProperties, + DclvOptions, Device3DPropertyResult, GetawayOptions, MomentOptions, @@ -33,7 +34,17 @@ Property3D.USRCAT, Property3D.GETAWAY, ) -SCALAR_PROPERTIES = tuple(prop for prop in Property3D if prop not in VECTOR_PROPERTIES) +DCLV_PROPERTIES = ( + Property3D.DCLV_SURFACE_AREA, + Property3D.DCLV_POLAR_SURFACE_AREA, + Property3D.DCLV_VOLUME, + Property3D.DCLV_VDW_VOLUME, + Property3D.DCLV_POLAR_VOLUME, + Property3D.DCLV_COMPACTNESS, + Property3D.DCLV_PACKING_DENSITY, +) +# Shape and moment scalars, compared against RDKit by _rdkit_property; the DCLV scalars have their own tests. +SCALAR_PROPERTIES = tuple(prop for prop in Property3D if prop not in VECTOR_PROPERTIES + DCLV_PROPERTIES) PAIRWISE_PROPERTIES = (Property3D.RDF, Property3D.MORSE, Property3D.AUTOCORR3D) RDKIT_PAIRWISE = { Property3D.RDF: rdMolDescriptors.CalcRDF, @@ -662,6 +673,132 @@ def test_getaway_precision_ignored_when_not_requested(precision_digits): assert result[Property3D.PMI1].torch().shape == (1,) +def _rdkit_dclv_rows(mols, options=None): + """RDKit's DoubleCubicLatticeVolume getters, one column per DCLV_PROPERTIES entry.""" + options = DclvOptions() if options is None else options + rows = [] + for mol in mols: + for conf in mol.GetConformers(): + dclv = rdMolDescriptors.DoubleCubicLatticeVolume(mol, probeRadius=options.probeRadius, confId=conf.GetId()) + polar = {"includeSandP": options.includeSandP, "includeHs": options.includeHs} + rows.append( + [ + dclv.GetSurfaceArea(), + dclv.GetPolarSurfaceArea(**polar), + dclv.GetVolume(), + dclv.GetVDWVolume(), + dclv.GetPolarVolume(**polar), + dclv.GetCompactness(), + dclv.GetPackingDensity(), + ] + ) + return np.asarray(rows) + + +def _dclv_columns(mols, precision=PrecisionMode.FULL, **kwargs): + """Calculate every DCLV property and stack them in DCLV_PROPERTIES order.""" + result = Calc3DProperties(mols, DCLV_PROPERTIES, precision=precision, **kwargs) + return np.column_stack([result[prop].numpy() for prop in DCLV_PROPERTIES]) + + +def _assert_dclv_matches(actual, expected, precision): + """Compare DCLV rows. + + The values count exposed surface dots, so a dot within rounding of a neighboring sphere could flip and move a + value by one dot's share (~1e-3 relative); none does for these geometries, and float32 stays within 1e-5. + """ + rtol = 1e-9 if precision == PrecisionMode.FULL else 1e-5 + np.testing.assert_allclose(actual, expected, rtol=rtol, atol=1e-9, equal_nan=True) + + +def _cubic_cluster(points_per_side, spacing): + grid = np.arange(points_per_side) * spacing + return np.asarray([(x, y, z) for x in grid for y in grid for z in grid]) + + +@pytest.mark.parametrize("precision", PRECISIONS) +def test_dclv_matches_rdkit(precision): + mols = [ + _embed("CC(=O)Nc1ccc(O)cc1", 3, 107), + _embed("CS(=O)(=O)Nc1ccc(P(=O)(O)O)cc1", 2, 109), + _embed("CC(=O)Nc1ccc(Br)cc1.Cl", 1, 113), + # 216 atoms 1.1 A apart: each atom has more neighbors than the per-warp neighbor cache holds. + _mol_with_conformers(".".join(["C"] * 216), [_cubic_cluster(6, 1.1)]), + _mol_with_conformers("[He]", [[(4.0, -3.0, 2.0)]]), + _mol_with_conformers("CC", [[(1.0, 1.0, 1.0), (1.0, 1.0, 1.0)]]), + ] + result = Calc3DProperties(mols, DCLV_PROPERTIES, precision=precision) + assert tuple(result) == tuple(prop.value for prop in DCLV_PROPERTIES) + assert all(result[prop].torch().shape == (9,) for prop in DCLV_PROPERTIES) + actual = np.column_stack([result[prop].numpy() for prop in DCLV_PROPERTIES]) + _assert_dclv_matches(actual, _rdkit_dclv_rows(mols), precision) + # The probe-expanded and bare-sphere dot passes run only for the properties that need them; a property's value + # must not depend on what else is requested. + for prop in DCLV_PROPERTIES: + alone = Calc3DProperties(mols, prop, precision=precision) + np.testing.assert_array_equal(alone[prop].numpy(), result[prop].numpy(), err_msg=prop.value) + + +@pytest.mark.parametrize( + "options", + [ + DclvOptions(probeRadius=0.0), + DclvOptions(probeRadius=2.0, includeSandP=True), + DclvOptions(includeHs=True), + DclvOptions(includeSandP=True, includeHs=True), + ], +) +def test_dclv_options_match_rdkit(options): + mols = [_embed("CS(=O)(=O)Nc1ccc(P(=O)(O)O)cc1", 2, 127)] + actual = _dclv_columns(mols, options=Property3DOptions(dclv=options)) + _assert_dclv_matches(actual, _rdkit_dclv_rows(mols, options), PrecisionMode.FULL) + + +@pytest.mark.parametrize("precision", PRECISIONS) +def test_dclv_zero_radius_atoms_give_nan(precision): + # RDKit crashes on dummy atoms (radius 0), so such conformers have no reference value. + coordinates = [(0.0, 0.0, 0.0), (1.5, 0.1, 0.0), (2.1, 1.4, 0.2)] + mols = [ + _mol_with_conformers("CCO.*", [[*coordinates, (0.8, 0.9, 0.1)]]), + _mol_with_conformers("CCO", [coordinates]), + ] + values = _dclv_columns(mols, precision) + assert np.isnan(values[0]).all() + _assert_dclv_matches(values[1:], _rdkit_dclv_rows(mols[1:]), precision) + + +@pytest.mark.parametrize("precision", PRECISIONS) +def test_dclv_non_finite_coordinates_give_nan(precision): + mols = [ + _mol_with_conformers("CCO", [[(0.0, 0.0, 0.0), (1.5, float("nan"), 0.0), (2.1, 1.4, 0.2)]]), + _mol_with_conformers("CCO", [[(0.0, 0.0, 0.0), (1.5, 0.1, 0.0), (2.1, 1.4, float("inf"))]]), + _embed("CCCO", 1, 47), + ] + values = _dclv_columns(mols, precision) + assert np.isnan(values[:2]).all() + _assert_dclv_matches(values[2:], _rdkit_dclv_rows(mols[2:]), precision) + + +@pytest.mark.parametrize("precision", PRECISIONS) +def test_dclv_device_coordinates_far_from_origin(precision): + mols = [_embed("CC(=O)Nc1ccc(O)cc1", 2, 131)] + coordinates = _device_result_from_molecules(mols) + near = _dclv_columns(mols, precision, coordinates=coordinates) + _assert_dclv_matches(near, _rdkit_dclv_rows(mols), precision) + coordinates.values.torch().add_(torch.tensor([1e4, -4e3, 6e3], dtype=torch.float64, device="cuda")) + far = _dclv_columns(mols, precision, coordinates=coordinates) + _assert_dclv_matches(far, near, precision) + + +def test_dclv_rejects_invalid_probe_radius_only_when_requested(): + mol = _embed("CCO", 1, 5) + for probe_radius in (-0.1, float("nan"), float("inf")): + options = Property3DOptions(dclv=DclvOptions(probeRadius=probe_radius)) + with pytest.raises(ValueError, match="probe radius"): + Calc3DProperties(mol, Property3D.DCLV_VDW_VOLUME, options=options) + assert Calc3DProperties(mol, "PMI1", options=options)["PMI1"].torch().shape == (1,) + + @pytest.mark.parametrize("precision", PRECISIONS) def test_autocorr3d_accepts_batches_without_bonds(precision): mols = [ diff --git a/src/descriptors3d.cu b/src/descriptors3d.cu index a34fa3c6..2c18477c 100644 --- a/src/descriptors3d.cu +++ b/src/descriptors3d.cu @@ -7,6 +7,7 @@ #include #include "src/descriptors3d.h" +#include "src/descriptors3d_dclv.cuh" #include "src/descriptors3d_getaway.cuh" #include "src/descriptors3d_kernel.cuh" #include "src/descriptors3d_moments.cuh" @@ -56,6 +57,20 @@ std::string_view property3DName(const Property3D property) { return "USRCAT"; case Property3D::GETAWAY: return "GETAWAY"; + case Property3D::DCLVSurfaceArea: + return "DCLVSurfaceArea"; + case Property3D::DCLVPolarSurfaceArea: + return "DCLVPolarSurfaceArea"; + case Property3D::DCLVVolume: + return "DCLVVolume"; + case Property3D::DCLVVDWVolume: + return "DCLVVDWVolume"; + case Property3D::DCLVPolarVolume: + return "DCLVPolarVolume"; + case Property3D::DCLVCompactness: + return "DCLVCompactness"; + case Property3D::DCLVPackingDensity: + return "DCLVPackingDensity"; } throw std::invalid_argument("Unknown Property3D value " + std::to_string(static_cast(property))); } @@ -75,12 +90,15 @@ using descriptors3d_detail::computeInertiaTensor; using descriptors3d_detail::computeMomentProperty; using descriptors3d_detail::computePrincipalMoments; using descriptors3d_detail::ConformerAtoms; +using descriptors3d_detail::dclvColumn; +using descriptors3d_detail::DclvOutputs; using descriptors3d_detail::kBlockSize; using descriptors3d_detail::kConformersPerBlock; using descriptors3d_detail::kGroupSize; using descriptors3d_detail::kGroupsPerWarp; using descriptors3d_detail::kWarpSize; using descriptors3d_detail::kWarpsPerBlock; +using descriptors3d_detail::launchDclvProperties; using descriptors3d_detail::launchGetawayProperties; using descriptors3d_detail::launchPairwiseProperties; using descriptors3d_detail::launchProjectionProperties; @@ -135,6 +153,13 @@ constexpr SharedStageSet directStages(const Property3D property) { case Property3D::USR: case Property3D::USRCAT: case Property3D::GETAWAY: + case Property3D::DCLVSurfaceArea: + case Property3D::DCLVPolarSurfaceArea: + case Property3D::DCLVVolume: + case Property3D::DCLVVDWVolume: + case Property3D::DCLVPolarVolume: + case Property3D::DCLVCompactness: + case Property3D::DCLVPackingDensity: return 0; } return 0; @@ -302,6 +327,7 @@ Property3DResults calc3DPropertiesGpu(const DeviceCoordView& coord Real* usrOutput = nullptr; Real* usrcatOutput = nullptr; Real* getawayOutput = nullptr; + DclvOutputs dclvOutputs; for (const Property3D property : properties) { // Bounds the work arrays, which hold one slot per known property. if (std::find(kAllProperty3D.begin(), kAllProperty3D.end(), property) == kAllProperty3D.end()) { @@ -331,6 +357,9 @@ Property3DResults calc3DPropertiesGpu(const DeviceCoordView& coord case Property3DFamily::Getaway: getawayOutput = it->second.data(); break; + case Property3DFamily::Dclv: + dclvOutputs.values[dclvColumn(property)] = it->second.data(); + break; } } @@ -341,6 +370,9 @@ Property3DResults calc3DPropertiesGpu(const DeviceCoordView& coord if (getawayOutput != nullptr && (options.getaway.precision < 1 || options.getaway.precision > 6)) { throw std::invalid_argument("GETAWAY precision must be between 1 and 6 significant digits"); } + if (dclvOutputs.any() && (!std::isfinite(options.dclv.probeRadius) || options.dclv.probeRadius < 0)) { + throw std::invalid_argument("DCLV probe radius must be finite and non-negative"); + } if (coordinates.numConformers == 0) { return results; @@ -370,6 +402,10 @@ Property3DResults calc3DPropertiesGpu(const DeviceCoordView& coord if (pairwiseOutputs.rdf != nullptr && inputs.iStateDragWeights == nullptr) { throw std::invalid_argument("I-state drag weights must not be null when RDF is requested"); } + if (dclvOutputs.any() && (inputs.vdwRadii == nullptr || inputs.dclvPolarClasses == nullptr)) { + throw std::invalid_argument( + "Van der Waals radii and polar classes must not be null when a DCLV property is requested"); + } const bool momentWeightsAreUnit = inputs.momentWeights == nullptr; const bool onlySpherocity = hasSpherocity && work.numProperties == 1; @@ -385,6 +421,7 @@ Property3DResults calc3DPropertiesGpu(const DeviceCoordView& coord launchPairwiseProperties(coordinates, inputs, pairwiseOutputs, bondTable.table, stream); launchUsrProperties(coordinates, inputs, usrOutput, usrcatOutput, stream); launchGetawayProperties(coordinates, inputs, options.getaway, bondTable.table, getawayOutput, stream); + launchDclvProperties(coordinates, inputs, options.dclv, dclvOutputs, stream); return results; } diff --git a/src/descriptors3d.h b/src/descriptors3d.h index 0c426c67..8535e92a 100644 --- a/src/descriptors3d.h +++ b/src/descriptors3d.h @@ -19,27 +19,35 @@ namespace nvMolKit { //! Per-conformer 3D properties. Names match the corresponding RDKit descriptor names. enum class Property3D : int { - PMI1 = 0, - PMI2 = 1, - PMI3 = 2, - RadiusOfGyration = 3, - NPR1 = 4, - NPR2 = 5, - InertialShapeFactor = 6, - Eccentricity = 7, - Asphericity = 8, - SpherocityIndex = 9, - PBF = 10, - WHIM = 11, - RDF = 12, - MORSE = 13, - AUTOCORR3D = 14, - USR = 15, - USRCAT = 16, - GETAWAY = 17, + PMI1 = 0, + PMI2 = 1, + PMI3 = 2, + RadiusOfGyration = 3, + NPR1 = 4, + NPR2 = 5, + InertialShapeFactor = 6, + Eccentricity = 7, + Asphericity = 8, + SpherocityIndex = 9, + PBF = 10, + WHIM = 11, + RDF = 12, + MORSE = 13, + AUTOCORR3D = 14, + USR = 15, + USRCAT = 16, + GETAWAY = 17, + //! RDKit's DoubleCubicLatticeVolume getters (GetSurfaceArea, ...). + DCLVSurfaceArea = 18, + DCLVPolarSurfaceArea = 19, + DCLVVolume = 20, + DCLVVDWVolume = 21, + DCLVPolarVolume = 22, + DCLVCompactness = 23, + DCLVPackingDensity = 24, }; -inline constexpr std::array kAllProperty3D = { +inline constexpr std::array kAllProperty3D = { Property3D::PMI1, Property3D::PMI2, Property3D::PMI3, @@ -58,6 +66,13 @@ inline constexpr std::array kAllProperty3D = { Property3D::USR, Property3D::USRCAT, Property3D::GETAWAY, + Property3D::DCLVSurfaceArea, + Property3D::DCLVPolarSurfaceArea, + Property3D::DCLVVolume, + Property3D::DCLVVDWVolume, + Property3D::DCLVPolarVolume, + Property3D::DCLVCompactness, + Property3D::DCLVPackingDensity, }; inline constexpr int kNumWhimProperties = 114; @@ -104,6 +119,7 @@ enum class Property3DFamily : int { Pairwise, //!< Sums over atom pairs weighted by atom-property pairs: RDF, MORSE and AUTOCORR3D. Usr, //!< Distance moments from four reference points: USR and USRCAT. Getaway, //!< Leverage (molecular influence) matrix descriptors: GETAWAY. + Dclv, //!< Dot-surface areas and volumes of the atom spheres: RDKit's DoubleCubicLatticeVolume (DCLV*). }; constexpr Property3DFamily property3DFamily(const Property3D property) { @@ -131,6 +147,14 @@ constexpr Property3DFamily property3DFamily(const Property3D property) { return Property3DFamily::Usr; case Property3D::GETAWAY: return Property3DFamily::Getaway; + case Property3D::DCLVSurfaceArea: + case Property3D::DCLVPolarSurfaceArea: + case Property3D::DCLVVolume: + case Property3D::DCLVVDWVolume: + case Property3D::DCLVPolarVolume: + case Property3D::DCLVCompactness: + case Property3D::DCLVPackingDensity: + return Property3DFamily::Dclv; } return Property3DFamily::Moments; } @@ -155,12 +179,31 @@ struct GetawayOptions { unsigned int precision = 2; }; +//! Property3DDeviceInputs::dclvPolarClasses bits: the atom kinds RDKit's DoubleCubicLatticeVolume counts as polar. +enum DclvPolarClass : uint8_t { + kDclvPolarNitrogenOxygen = 1u << 0, //!< Nitrogen or oxygen; always polar. + kDclvPolarSulfurPhosphorus = 1u << 1, //!< Sulfur or phosphorus; polar with DclvOptions::includeSandP. + kDclvPolarHydrogenOnNO = 1u << 2, //!< Hydrogen bonded to nitrogen or oxygen; with includeHs. + kDclvPolarHydrogenOnSulfurPhos = 1u << 3, //!< Hydrogen bonded to sulfur or phosphorus; with both options. +}; + +//! Options for the DCLV properties (RDKit's DoubleCubicLatticeVolume), with RDKit's defaults. +struct DclvOptions { + //! Radius of the solvent probe sphere in Angstrom. Must be finite and non-negative. + double probeRadius = 1.4; + //! DCLVPolarSurfaceArea and DCLVPolarVolume also count sulfur and phosphorus (RDKit's includeSandP). + bool includeSandP = false; + //! The polar DCLV properties also count hydrogens bonded to a polar atom (RDKit's includeHs). + bool includeHs = false; +}; + //! Per-family options; each family reads only its own member. PBF, RDF, MORSE, AUTOCORR3D, USR and //! USRCAT have no options. struct Property3DOptions { MomentOptions moments; WhimOptions whim; GetawayOptions getaway; + DclvOptions dclv; }; /** @@ -197,6 +240,12 @@ struct Property3DDeviceInputs { //! GETAWAY: per molecule, the coordinate row of its default (first) conformer, whose PBF decides HIC's //! dimension as in RDKit. Null, or a negative entry, uses each row's own conformer. const int32_t* defaultConformerRows = nullptr; + //! DCLV: per atom, RDKit's van der Waals radius (PeriodicTable::getRvdw). Conformers with a radius-0 atom (dummy + //! atoms, on which RDKit's DoubleCubicLatticeVolume crashes) produce NaN for every DCLV property. + const double* vdwRadii = nullptr; + //! DCLV: per atom, its DclvPolarClass bits. DclvOptions selects the classes DCLVPolarSurfaceArea and + //! DCLVPolarVolume count. + const uint8_t* dclvPolarClasses = nullptr; //! PBF, GETAWAY: per-conformer RDKit is3D flags (one per coordinate row); null treats every row as 3D. const int8_t* conformerIs3D = nullptr; //! WHIM, AUTOCORR3D, GETAWAY: largest molecule atom count in the batch (host value). Sizes WHIM's @@ -227,7 +276,9 @@ template using Property3DResults = std::unordered_map diff --git a/src/descriptors3d_dclv.cuh b/src/descriptors3d_dclv.cuh new file mode 100644 index 00000000..7b0af148 --- /dev/null +++ b/src/descriptors3d_dclv.cuh @@ -0,0 +1,743 @@ +// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +// SPDX-License-Identifier: Apache-2.0 + +#ifndef NVMOLKIT_DESCRIPTORS3D_DCLV_CUH +#define NVMOLKIT_DESCRIPTORS3D_DCLV_CUH + +#include +#include + +#include "src/descriptors3d.h" +#include "src/descriptors3d_kernel.cuh" +#include "src/utils/cuda_error_check.h" +#include "src/utils/device_vector.h" + +namespace nvMolKit::descriptors3d_detail { + +//! The DCLV properties, in RDKit's DoubleCubicLatticeVolume getter declaration order. +enum DclvColumn : int { + kDclvSurfaceArea = 0, + kDclvPolarSurfaceArea = 1, + kDclvVolume = 2, + kDclvVdwVolume = 3, + kDclvPolarVolume = 4, + kDclvCompactness = 5, + kDclvPackingDensity = 6, + kNumDclvColumns = 7, +}; + +//! Column of a Dclv-family property. +constexpr int dclvColumn(const Property3D property) { + switch (property) { + case Property3D::DCLVPolarSurfaceArea: + return kDclvPolarSurfaceArea; + case Property3D::DCLVVolume: + return kDclvVolume; + case Property3D::DCLVVDWVolume: + return kDclvVdwVolume; + case Property3D::DCLVPolarVolume: + return kDclvPolarVolume; + case Property3D::DCLVCompactness: + return kDclvCompactness; + case Property3D::DCLVPackingDensity: + return kDclvPackingDensity; + default: + return kDclvSurfaceArea; + } +} + +//! Output buffer of each requested DCLV property, indexed by DclvColumn; null when not requested. +template struct DclvOutputs { + Real* values[kNumDclvColumns] = {}; + + bool any() const { + for (const Real* output : values) { + if (output != nullptr) { + return true; + } + } + return false; + } +}; + +//! Polar classes DCLVPolarSurfaceArea and DCLVPolarVolume count under @p options, following RDKit's includeAsPolar. +inline uint8_t dclvPolarMask(const DclvOptions& options) { + uint8_t mask = kDclvPolarNitrogenOxygen; + if (options.includeSandP) { + mask |= kDclvPolarSulfurPhosphorus; + } + if (options.includeHs) { + mask |= kDclvPolarHydrogenOnNO; + if (options.includeSandP) { + mask |= kDclvPolarHydrogenOnSulfurPhos; + } + } + return mask; +} + +//! Surface dots tested per atom sphere. +constexpr int kDclvNumDots = 320; +//! RDKit's area of one dot and of the whole dot set on the unit sphere (Code/GraphMol/Descriptors/DCLV_dots.h). +constexpr double kDclvDotArea = 0.03853087; +constexpr double kDclvStandardArea = 12.3298; +//! Neighbors cached per warp in shared memory; atoms with more scan the whole conformer for every dot instead. +constexpr int kDclvNeighborCapacity = 128; +//! Widens the neighbor cutoff (Angstrom). Any atom that buries a dot lies strictly inside the cutoff, so a +//! wider list only adds atoms that bury nothing; the margin keeps rounding from dropping a borderline atom. +constexpr double kDclvNeighborMargin = 0.01; + +/** + * RDKit's dot set, 320 unit vectors "based on a 320 faced polyhedron", copied from + * Code/GraphMol/Descriptors/DCLV_dots.h. RDKit stores them as doubles initialized from these float literals, so + * converting the floats reproduces its values exactly. + */ +__device__ const float kDclvDots[kDclvNumDots][3] = { + {-0.577350f, -0.577350f, -0.577350f}, + {-0.448278f, -0.547076f, -0.706934f}, + {-0.706934f, -0.448278f, -0.547076f}, + {-0.547076f, -0.706934f, -0.448278f}, + {-0.286211f, -0.762327f, -0.580467f}, + {-0.142204f, -0.816797f, -0.559125f}, + {-0.302648f, -0.639218f, -0.706968f}, + {-0.400803f, -0.798035f, -0.449997f}, + {-0.580467f, -0.286211f, -0.762327f}, + {-0.449997f, -0.400803f, -0.798035f}, + {-0.559125f, -0.142204f, -0.816797f}, + {-0.706968f, -0.302648f, -0.639218f}, + {-0.762327f, -0.580467f, -0.286211f}, + {-0.639218f, -0.706968f, -0.302648f}, + {-0.798035f, -0.449997f, -0.400803f}, + {-0.816797f, -0.559125f, -0.142204f}, + { 0.000000f, -0.356822f, -0.934172f}, + {-0.159858f, -0.436909f, -0.885187f}, + { 0.000000f, -0.178253f, -0.983985f}, + { 0.159858f, -0.436909f, -0.885187f}, + { 0.000000f, -0.653004f, -0.757355f}, + { 0.000000f, -0.762480f, -0.647012f}, + {-0.158817f, -0.584279f, -0.795861f}, + { 0.158817f, -0.584279f, -0.795861f}, + {-0.294256f, -0.176888f, -0.939215f}, + {-0.306166f, -0.345864f, -0.886928f}, + {-0.416921f, -0.087887f, -0.904684f}, + {-0.147349f, -0.088893f, -0.985082f}, + { 0.294256f, -0.176888f, -0.939215f}, + { 0.306166f, -0.345864f, -0.886928f}, + { 0.147349f, -0.088893f, -0.985082f}, + { 0.416921f, -0.087887f, -0.904684f}, + {-0.356822f, -0.934172f, 0.000000f}, + {-0.436909f, -0.885187f, -0.159858f}, + {-0.436909f, -0.885187f, 0.159858f}, + {-0.178253f, -0.983985f, 0.000000f}, + {-0.176888f, -0.939215f, -0.294256f}, + {-0.087887f, -0.904684f, -0.416921f}, + {-0.345864f, -0.886928f, -0.306166f}, + {-0.088893f, -0.985082f, -0.147349f}, + {-0.653004f, -0.757355f, 0.000000f}, + {-0.584279f, -0.795861f, -0.158817f}, + {-0.762480f, -0.647012f, 0.000000f}, + {-0.584279f, -0.795861f, 0.158817f}, + {-0.176888f, -0.939215f, 0.294256f}, + {-0.088893f, -0.985082f, 0.147349f}, + {-0.345864f, -0.886928f, 0.306166f}, + {-0.087887f, -0.904684f, 0.416921f}, + { 0.356822f, -0.934172f, 0.000000f}, + { 0.178253f, -0.983985f, 0.000000f}, + { 0.436909f, -0.885187f, 0.159858f}, + { 0.436909f, -0.885187f, -0.159858f}, + { 0.176888f, -0.939215f, -0.294256f}, + { 0.087887f, -0.904684f, -0.416921f}, + { 0.088893f, -0.985082f, -0.147349f}, + { 0.345864f, -0.886928f, -0.306166f}, + { 0.176888f, -0.939215f, 0.294256f}, + { 0.088893f, -0.985082f, 0.147349f}, + { 0.087887f, -0.904684f, 0.416921f}, + { 0.345864f, -0.886928f, 0.306166f}, + { 0.653004f, -0.757355f, 0.000000f}, + { 0.584279f, -0.795861f, -0.158817f}, + { 0.584279f, -0.795861f, 0.158817f}, + { 0.762480f, -0.647012f, 0.000000f}, + { 0.577350f, -0.577350f, -0.577350f}, + { 0.448278f, -0.547076f, -0.706934f}, + { 0.706934f, -0.448278f, -0.547076f}, + { 0.547076f, -0.706934f, -0.448278f}, + { 0.286211f, -0.762327f, -0.580467f}, + { 0.142204f, -0.816797f, -0.559125f}, + { 0.302648f, -0.639218f, -0.706968f}, + { 0.400803f, -0.798035f, -0.449997f}, + { 0.580467f, -0.286211f, -0.762327f}, + { 0.449997f, -0.400803f, -0.798035f}, + { 0.559125f, -0.142204f, -0.816797f}, + { 0.706968f, -0.302648f, -0.639218f}, + { 0.762327f, -0.580467f, -0.286211f}, + { 0.639218f, -0.706968f, -0.302648f}, + { 0.798035f, -0.449997f, -0.400803f}, + { 0.816797f, -0.559125f, -0.142204f}, + {-0.934172f, 0.000000f, -0.356822f}, + {-0.885187f, -0.159858f, -0.436909f}, + {-0.983985f, 0.000000f, -0.178253f}, + {-0.885187f, 0.159858f, -0.436909f}, + {-0.757355f, 0.000000f, -0.653004f}, + {-0.647012f, 0.000000f, -0.762480f}, + {-0.795861f, -0.158817f, -0.584279f}, + {-0.795861f, 0.158817f, -0.584279f}, + {-0.939215f, -0.294256f, -0.176888f}, + {-0.886928f, -0.306166f, -0.345864f}, + {-0.904684f, -0.416921f, -0.087887f}, + {-0.985082f, -0.147349f, -0.088893f}, + {-0.939215f, 0.294256f, -0.176888f}, + {-0.886928f, 0.306166f, -0.345864f}, + {-0.985082f, 0.147349f, -0.088893f}, + {-0.904684f, 0.416921f, -0.087887f}, + { 0.000000f, 0.356822f, -0.934172f}, + { 0.000000f, 0.178253f, -0.983985f}, + { 0.159858f, 0.436909f, -0.885187f}, + {-0.159858f, 0.436909f, -0.885187f}, + {-0.294256f, 0.176888f, -0.939215f}, + {-0.416921f, 0.087887f, -0.904684f}, + {-0.147349f, 0.088893f, -0.985082f}, + {-0.306166f, 0.345864f, -0.886928f}, + { 0.294256f, 0.176888f, -0.939215f}, + { 0.147349f, 0.088893f, -0.985082f}, + { 0.416921f, 0.087887f, -0.904684f}, + { 0.306166f, 0.345864f, -0.886928f}, + { 0.000000f, 0.653004f, -0.757355f}, + {-0.158817f, 0.584279f, -0.795861f}, + { 0.158817f, 0.584279f, -0.795861f}, + { 0.000000f, 0.762480f, -0.647012f}, + {-0.577350f, 0.577350f, -0.577350f}, + {-0.706934f, 0.448278f, -0.547076f}, + {-0.547076f, 0.706934f, -0.448278f}, + {-0.448278f, 0.547076f, -0.706934f}, + {-0.580467f, 0.286211f, -0.762327f}, + {-0.559125f, 0.142204f, -0.816797f}, + {-0.706968f, 0.302648f, -0.639218f}, + {-0.449997f, 0.400803f, -0.798035f}, + {-0.762327f, 0.580467f, -0.286211f}, + {-0.798035f, 0.449997f, -0.400803f}, + {-0.816797f, 0.559125f, -0.142204f}, + {-0.639218f, 0.706968f, -0.302648f}, + {-0.286211f, 0.762327f, -0.580467f}, + {-0.302648f, 0.639218f, -0.706968f}, + {-0.400803f, 0.798035f, -0.449997f}, + {-0.142204f, 0.816797f, -0.559125f}, + {-0.577350f, -0.577350f, 0.577350f}, + {-0.547076f, -0.706934f, 0.448278f}, + {-0.448278f, -0.547076f, 0.706934f}, + {-0.706934f, -0.448278f, 0.547076f}, + {-0.762327f, -0.580467f, 0.286211f}, + {-0.816797f, -0.559125f, 0.142204f}, + {-0.639218f, -0.706968f, 0.302648f}, + {-0.798035f, -0.449997f, 0.400803f}, + {-0.286211f, -0.762327f, 0.580467f}, + {-0.400803f, -0.798035f, 0.449997f}, + {-0.142204f, -0.816797f, 0.559125f}, + {-0.302648f, -0.639218f, 0.706968f}, + {-0.580467f, -0.286211f, 0.762327f}, + {-0.706968f, -0.302648f, 0.639218f}, + {-0.449997f, -0.400803f, 0.798035f}, + {-0.559125f, -0.142204f, 0.816797f}, + {-0.934172f, 0.000000f, 0.356822f}, + {-0.983985f, 0.000000f, 0.178253f}, + {-0.885187f, 0.159858f, 0.436909f}, + {-0.885187f, -0.159858f, 0.436909f}, + {-0.939215f, -0.294256f, 0.176888f}, + {-0.904684f, -0.416921f, 0.087887f}, + {-0.985082f, -0.147349f, 0.088893f}, + {-0.886928f, -0.306166f, 0.345864f}, + {-0.939215f, 0.294256f, 0.176888f}, + {-0.985082f, 0.147349f, 0.088893f}, + {-0.904684f, 0.416921f, 0.087887f}, + {-0.886928f, 0.306166f, 0.345864f}, + {-0.757355f, 0.000000f, 0.653004f}, + {-0.795861f, -0.158817f, 0.584279f}, + {-0.795861f, 0.158817f, 0.584279f}, + {-0.647012f, 0.000000f, 0.762480f}, + { 0.934172f, 0.000000f, -0.356822f}, + { 0.885187f, -0.159858f, -0.436909f}, + { 0.983985f, 0.000000f, -0.178253f}, + { 0.885187f, 0.159858f, -0.436909f}, + { 0.757355f, 0.000000f, -0.653004f}, + { 0.647012f, 0.000000f, -0.762480f}, + { 0.795861f, -0.158817f, -0.584279f}, + { 0.795861f, 0.158817f, -0.584279f}, + { 0.939215f, -0.294256f, -0.176888f}, + { 0.886928f, -0.306166f, -0.345864f}, + { 0.904684f, -0.416921f, -0.087887f}, + { 0.985082f, -0.147349f, -0.088893f}, + { 0.939215f, 0.294256f, -0.176888f}, + { 0.886928f, 0.306166f, -0.345864f}, + { 0.985082f, 0.147349f, -0.088893f}, + { 0.904684f, 0.416921f, -0.087887f}, + { 0.577350f, 0.577350f, -0.577350f}, + { 0.448278f, 0.547076f, -0.706934f}, + { 0.547076f, 0.706934f, -0.448278f}, + { 0.706934f, 0.448278f, -0.547076f}, + { 0.580467f, 0.286211f, -0.762327f}, + { 0.559125f, 0.142204f, -0.816797f}, + { 0.449997f, 0.400803f, -0.798035f}, + { 0.706968f, 0.302648f, -0.639218f}, + { 0.286211f, 0.762327f, -0.580467f}, + { 0.302648f, 0.639218f, -0.706968f}, + { 0.142204f, 0.816797f, -0.559125f}, + { 0.400803f, 0.798035f, -0.449997f}, + { 0.762327f, 0.580467f, -0.286211f}, + { 0.798035f, 0.449997f, -0.400803f}, + { 0.639218f, 0.706968f, -0.302648f}, + { 0.816797f, 0.559125f, -0.142204f}, + { 0.577350f, -0.577350f, 0.577350f}, + { 0.547076f, -0.706934f, 0.448278f}, + { 0.706934f, -0.448278f, 0.547076f}, + { 0.448278f, -0.547076f, 0.706934f}, + { 0.286211f, -0.762327f, 0.580467f}, + { 0.142204f, -0.816797f, 0.559125f}, + { 0.400803f, -0.798035f, 0.449997f}, + { 0.302648f, -0.639218f, 0.706968f}, + { 0.762327f, -0.580467f, 0.286211f}, + { 0.639218f, -0.706968f, 0.302648f}, + { 0.816797f, -0.559125f, 0.142204f}, + { 0.798035f, -0.449997f, 0.400803f}, + { 0.580467f, -0.286211f, 0.762327f}, + { 0.449997f, -0.400803f, 0.798035f}, + { 0.706968f, -0.302648f, 0.639218f}, + { 0.559125f, -0.142204f, 0.816797f}, + { 0.000000f, -0.356822f, 0.934172f}, + {-0.159858f, -0.436909f, 0.885187f}, + { 0.000000f, -0.178253f, 0.983985f}, + { 0.159858f, -0.436909f, 0.885187f}, + { 0.000000f, -0.653004f, 0.757355f}, + { 0.000000f, -0.762480f, 0.647012f}, + {-0.158817f, -0.584279f, 0.795861f}, + { 0.158817f, -0.584279f, 0.795861f}, + {-0.294256f, -0.176888f, 0.939215f}, + {-0.306166f, -0.345864f, 0.886928f}, + {-0.416921f, -0.087887f, 0.904684f}, + {-0.147349f, -0.088893f, 0.985082f}, + { 0.294256f, -0.176888f, 0.939215f}, + { 0.306166f, -0.345864f, 0.886928f}, + { 0.147349f, -0.088893f, 0.985082f}, + { 0.416921f, -0.087887f, 0.904684f}, + { 0.934172f, 0.000000f, 0.356822f}, + { 0.885187f, -0.159858f, 0.436909f}, + { 0.885187f, 0.159858f, 0.436909f}, + { 0.983985f, 0.000000f, 0.178253f}, + { 0.939215f, -0.294256f, 0.176888f}, + { 0.904684f, -0.416921f, 0.087887f}, + { 0.886928f, -0.306166f, 0.345864f}, + { 0.985082f, -0.147349f, 0.088893f}, + { 0.757355f, 0.000000f, 0.653004f}, + { 0.795861f, -0.158817f, 0.584279f}, + { 0.647012f, 0.000000f, 0.762480f}, + { 0.795861f, 0.158817f, 0.584279f}, + { 0.939215f, 0.294256f, 0.176888f}, + { 0.985082f, 0.147349f, 0.088893f}, + { 0.886928f, 0.306166f, 0.345864f}, + { 0.904684f, 0.416921f, 0.087887f}, + {-0.577350f, 0.577350f, 0.577350f}, + {-0.706934f, 0.448278f, 0.547076f}, + {-0.448278f, 0.547076f, 0.706934f}, + {-0.547076f, 0.706934f, 0.448278f}, + {-0.762327f, 0.580467f, 0.286211f}, + {-0.816797f, 0.559125f, 0.142204f}, + {-0.798035f, 0.449997f, 0.400803f}, + {-0.639218f, 0.706968f, 0.302648f}, + {-0.580467f, 0.286211f, 0.762327f}, + {-0.706968f, 0.302648f, 0.639218f}, + {-0.559125f, 0.142204f, 0.816797f}, + {-0.449997f, 0.400803f, 0.798035f}, + {-0.286211f, 0.762327f, 0.580467f}, + {-0.400803f, 0.798035f, 0.449997f}, + {-0.302648f, 0.639218f, 0.706968f}, + {-0.142204f, 0.816797f, 0.559125f}, + {-0.356822f, 0.934172f, 0.000000f}, + {-0.436909f, 0.885187f, -0.159858f}, + {-0.178253f, 0.983985f, 0.000000f}, + {-0.436909f, 0.885187f, 0.159858f}, + {-0.653004f, 0.757355f, 0.000000f}, + {-0.762480f, 0.647012f, 0.000000f}, + {-0.584279f, 0.795861f, -0.158817f}, + {-0.584279f, 0.795861f, 0.158817f}, + {-0.176888f, 0.939215f, -0.294256f}, + {-0.345864f, 0.886928f, -0.306166f}, + {-0.087887f, 0.904684f, -0.416921f}, + {-0.088893f, 0.985082f, -0.147349f}, + {-0.176888f, 0.939215f, 0.294256f}, + {-0.345864f, 0.886928f, 0.306166f}, + {-0.088893f, 0.985082f, 0.147349f}, + {-0.087887f, 0.904684f, 0.416921f}, + { 0.356822f, 0.934172f, 0.000000f}, + { 0.178253f, 0.983985f, 0.000000f}, + { 0.436909f, 0.885187f, 0.159858f}, + { 0.436909f, 0.885187f, -0.159858f}, + { 0.176888f, 0.939215f, -0.294256f}, + { 0.087887f, 0.904684f, -0.416921f}, + { 0.088893f, 0.985082f, -0.147349f}, + { 0.345864f, 0.886928f, -0.306166f}, + { 0.176888f, 0.939215f, 0.294256f}, + { 0.088893f, 0.985082f, 0.147349f}, + { 0.087887f, 0.904684f, 0.416921f}, + { 0.345864f, 0.886928f, 0.306166f}, + { 0.653004f, 0.757355f, 0.000000f}, + { 0.584279f, 0.795861f, -0.158817f}, + { 0.584279f, 0.795861f, 0.158817f}, + { 0.762480f, 0.647012f, 0.000000f}, + { 0.000000f, 0.356822f, 0.934172f}, + {-0.159858f, 0.436909f, 0.885187f}, + { 0.159858f, 0.436909f, 0.885187f}, + { 0.000000f, 0.178253f, 0.983985f}, + {-0.294256f, 0.176888f, 0.939215f}, + {-0.416921f, 0.087887f, 0.904684f}, + {-0.306166f, 0.345864f, 0.886928f}, + {-0.147349f, 0.088893f, 0.985082f}, + { 0.000000f, 0.653004f, 0.757355f}, + {-0.158817f, 0.584279f, 0.795861f}, + { 0.000000f, 0.762480f, 0.647012f}, + { 0.158817f, 0.584279f, 0.795861f}, + { 0.294256f, 0.176888f, 0.939215f}, + { 0.147349f, 0.088893f, 0.985082f}, + { 0.306166f, 0.345864f, 0.886928f}, + { 0.416921f, 0.087887f, 0.904684f}, + { 0.577350f, 0.577350f, 0.577350f}, + { 0.448278f, 0.547076f, 0.706934f}, + { 0.706934f, 0.448278f, 0.547076f}, + { 0.547076f, 0.706934f, 0.448278f}, + { 0.286211f, 0.762327f, 0.580467f}, + { 0.142204f, 0.816797f, 0.559125f}, + { 0.302648f, 0.639218f, 0.706968f}, + { 0.400803f, 0.798035f, 0.449997f}, + { 0.580467f, 0.286211f, 0.762327f}, + { 0.449997f, 0.400803f, 0.798035f}, + { 0.559125f, 0.142204f, 0.816797f}, + { 0.706968f, 0.302648f, 0.639218f}, + { 0.762327f, 0.580467f, 0.286211f}, + { 0.639218f, 0.706968f, 0.302648f}, + { 0.798035f, 0.449997f, 0.400803f}, + { 0.816797f, 0.559125f, 0.142204f} +}; + +//! A cached neighbor: its position relative to the current atom and its radius. +template struct DclvNeighbor { + Real x; + Real y; + Real z; + Real radius; +}; + +//! Whether a sphere of @p reach around (@p x, @p y, @p z) contains the point (@p vx, @p vy, @p vz); strict, as in +//! RDKit's testPoint. +template +__device__ __forceinline__ bool +dclvCovers(const Real x, const Real y, const Real z, const Real reach, const Real vx, const Real vy, const Real vz) { + const Real dx = x - vx; + const Real dy = y - vy; + const Real dz = z - vz; + return dx * dx + dy * dy + dz * dz < reach * reach; +} + +//! Exposed dots of one atom sphere, for the probe-expanded sphere (surface area and volume) and the bare van der +//! Waals sphere (VDW volume): dot counts and sums of the exposed unit vectors. +template struct DclvAtomDots { + Real probeCount = 0; + Real probeX = 0; + Real probeY = 0; + Real probeZ = 0; + Real vdwCount = 0; + Real vdwX = 0; + Real vdwY = 0; + Real vdwZ = 0; +}; + +/** + * @brief Warp-collective. Tests atom @p atomIdx's dots against the other atoms and returns the warp-wide exposed + * dot counts and unit-vector sums to every lane. + * + * With @p kProbe, a dot on the probe-expanded sphere (radius r + probe) is exposed unless it lies inside another + * atom's probe-expanded sphere; with @p kVdw, a dot on the bare sphere (radius r) is exposed unless it lies inside + * another atom's bare sphere. Both use the neighbors of RDKit's findNeighbours, atoms closer than + * r_i + 2 probe + r_j, which are the only atoms that can bury either kind of dot. Disabled kinds stay zero. + */ +template +__device__ __forceinline__ DclvAtomDots dclvAtomDots(const Real* centered, + const double* radii, + const int numAtoms, + const int atomIdx, + const Real probe, + const int lane, + DclvNeighbor* neighbors) { + const Real px = centered[atomIdx * 3 + 0]; + const Real py = centered[atomIdx * 3 + 1]; + const Real pz = centered[atomIdx * 3 + 2]; + const Real radius = static_cast(radii[atomIdx]); + const Real range = radius + Real(2) * probe + static_cast(kDclvNeighborMargin); + + // Every lane has finished reading the previous atom's neighbors before they are overwritten. + __syncwarp(); + int numNeighbors = 0; + for (int base = 0; base < numAtoms; base += kWarpSize) { + const int other = base + lane; + DclvNeighbor candidate{}; + bool isNeighbor = false; + if (other < numAtoms && other != atomIdx) { + candidate.radius = static_cast(radii[other]); + candidate.x = centered[other * 3 + 0] - px; + candidate.y = centered[other * 3 + 1] - py; + candidate.z = centered[other * 3 + 2] - pz; + isNeighbor = + dclvCovers(candidate.x, candidate.y, candidate.z, range + candidate.radius, Real(0), Real(0), Real(0)); + } + const unsigned ballot = __ballot_sync(0xffffffffu, isNeighbor); + const int slot = numNeighbors + __popc(ballot & ((1u << lane) - 1u)); + if (isNeighbor && slot < kDclvNeighborCapacity) { + neighbors[slot] = candidate; + } + numNeighbors += __popc(ballot); + } + __syncwarp(); + const bool cached = numNeighbors <= kDclvNeighborCapacity; + + DclvAtomDots dots; + const Real probeScale = radius + probe; + for (int dotIdx = lane; dotIdx < kDclvNumDots; dotIdx += kWarpSize) { + const Real ux = static_cast(kDclvDots[dotIdx][0]); + const Real uy = static_cast(kDclvDots[dotIdx][1]); + const Real uz = static_cast(kDclvDots[dotIdx][2]); + // The dot on each sphere, relative to the atom, stays in registers across the neighbor loop. + const Real probeDotX = ux * probeScale; + const Real probeDotY = uy * probeScale; + const Real probeDotZ = uz * probeScale; + const Real vdwDotX = ux * radius; + const Real vdwDotY = uy * radius; + const Real vdwDotZ = uz * radius; + bool probeOpen = kProbe; + bool vdwOpen = kVdw; + auto test = [&](const DclvNeighbor& other) { + if constexpr (kProbe) { + probeOpen = + probeOpen && !dclvCovers(other.x, other.y, other.z, other.radius + probe, probeDotX, probeDotY, probeDotZ); + } + if constexpr (kVdw) { + vdwOpen = vdwOpen && !dclvCovers(other.x, other.y, other.z, other.radius, vdwDotX, vdwDotY, vdwDotZ); + } + }; + if (cached) { + for (int neighborIdx = 0; neighborIdx < numNeighbors && (probeOpen || vdwOpen); ++neighborIdx) { + test(neighbors[neighborIdx]); + } + } else { + for (int other = 0; other < numAtoms && (probeOpen || vdwOpen); ++other) { + if (other != atomIdx) { + test(DclvNeighbor{centered[other * 3 + 0] - px, + centered[other * 3 + 1] - py, + centered[other * 3 + 2] - pz, + static_cast(radii[other])}); + } + } + } + if (probeOpen) { + dots.probeCount += Real(1); + dots.probeX += ux; + dots.probeY += uy; + dots.probeZ += uz; + } + if (vdwOpen) { + dots.vdwCount += Real(1); + dots.vdwX += ux; + dots.vdwY += uy; + dots.vdwZ += uz; + } + } + if constexpr (kProbe) { + dots.probeCount = groupAllReduceSum(dots.probeCount); + dots.probeX = groupAllReduceSum(dots.probeX); + dots.probeY = groupAllReduceSum(dots.probeY); + dots.probeZ = groupAllReduceSum(dots.probeZ); + } + if constexpr (kVdw) { + dots.vdwCount = groupAllReduceSum(dots.vdwCount); + dots.vdwX = groupAllReduceSum(dots.vdwX); + dots.vdwY = groupAllReduceSum(dots.vdwY); + dots.vdwZ = groupAllReduceSum(dots.vdwZ); + } + return dots; +} + +/** + * @brief RDKit's DoubleCubicLatticeVolume for every conformer; one warp per conformer. + * + * Each warp centers its conformer on the mean atom position (RDKit's centreOfGravity) in float64 (see + * centeredPosition()), then walks the atoms: it gathers each atom's neighbors in shared memory and each lane tests + * every kWarpSize-th dot. Per atom, RDKit's getAtomSurfaceArea gives + * area = exposed * dotArea * 4 pi f^2 / standardArea and its getAtomVolume (Gauss-Ostrogradskii) gives + * volume = f^2 ((x - centre) . sum of exposed dots + f exposed) 4 pi / (3 * 320), with f = r + probe for the + * surface and volume (@p kProbe) and f = r for the VDW volume (@p kVdw). Only requested outputs are written. Rows + * with a non-finite coordinate or a zero-radius (dummy) atom produce NaN. + */ +template +__global__ void dclv3DKernel(const DeviceCoordView coordinates, + const Property3DDeviceInputs inputs, + const Real probe, + const uint8_t polarMask, + Real* __restrict__ centered, + const DclvOutputs outputs) { + __shared__ DclvNeighbor blockNeighbors[kWarpsPerBlock][kDclvNeighborCapacity]; + const int lane = static_cast(threadIdx.x) % kWarpSize; + const int warpInBlock = static_cast(threadIdx.x) / kWarpSize; + const int conformerIdx = blockIdx.x * kWarpsPerBlock + warpInBlock; + if (conformerIdx >= coordinates.numConformers) { + return; // Uniform across the warp. + } + auto writeRow = [&](const Real(&values)[kNumDclvColumns]) { + if (lane == 0) { + for (int column = 0; column < kNumDclvColumns; ++column) { + if (outputs.values[column] != nullptr) { + outputs.values[column][conformerIdx] = values[column]; + } + } + } + }; + const ConformerAtoms atoms = loadConformer(coordinates, nullptr, inputs.moleculeAtomStarts, conformerIdx); + bool valid = atoms.valid; + const double* radii = nullptr; + const uint8_t* polar = nullptr; + double centreX = 0; + double centreY = 0; + double centreZ = 0; + if (valid) { + const int moleculeStart = inputs.moleculeAtomStarts[coordinates.molIndices[conformerIdx]]; + radii = inputs.vdwRadii + moleculeStart; + polar = inputs.dclvPolarClasses + moleculeStart; + double sumX = 0; + double sumY = 0; + double sumZ = 0; + bool usable = true; + for (int atomIdx = lane; atomIdx < atoms.numAtoms; atomIdx += kWarpSize) { + const double x = atoms.positions[atomIdx * 3 + 0]; + const double y = atoms.positions[atomIdx * 3 + 1]; + const double z = atoms.positions[atomIdx * 3 + 2]; + usable = usable && isfinite(x) && isfinite(y) && isfinite(z) && radii[atomIdx] != 0.0; + sumX += x; + sumY += y; + sumZ += z; + } + valid = __all_sync(0xffffffffu, usable); + const double inverseAtoms = 1.0 / static_cast(atoms.numAtoms); + centreX = groupAllReduceSum(sumX) * inverseAtoms; + centreY = groupAllReduceSum(sumY) * inverseAtoms; + centreZ = groupAllReduceSum(sumZ) * inverseAtoms; + } + if (!valid) { + Real nanRow[kNumDclvColumns]; + for (Real& value : nanRow) { + value = static_cast(nan("")); + } + writeRow(nanRow); + return; + } + + Real* const conformerCentered = centered + (atoms.positions - coordinates.positions); + for (int atomIdx = lane; atomIdx < atoms.numAtoms; atomIdx += kWarpSize) { + centeredPosition(atoms.positions, + atomIdx, + centreX, + centreY, + centreZ, + conformerCentered[atomIdx * 3 + 0], + conformerCentered[atomIdx * 3 + 1], + conformerCentered[atomIdx * 3 + 2]); + } + __syncwarp(); + + constexpr Real kAreaScale = static_cast(kDclvDotArea * 4.0 * M_PI / kDclvStandardArea); + constexpr Real kVolumeScale = static_cast(4.0 / 3.0 * M_PI / kDclvNumDots); + Real surfaceArea = 0; + Real polarArea = 0; + Real volume = 0; + Real polarVolume = 0; + Real vdwVolume = 0; + for (int atomIdx = 0; atomIdx < atoms.numAtoms; ++atomIdx) { + const DclvAtomDots dots = dclvAtomDots(conformerCentered, + radii, + atoms.numAtoms, + atomIdx, + probe, + lane, + blockNeighbors[warpInBlock]); + const Real x = conformerCentered[atomIdx * 3 + 0]; + const Real y = conformerCentered[atomIdx * 3 + 1]; + const Real z = conformerCentered[atomIdx * 3 + 2]; + const Real radius = static_cast(radii[atomIdx]); + if constexpr (kProbe) { + const Real probeScale = radius + probe; + const Real atomArea = dots.probeCount * kAreaScale * probeScale * probeScale; + const Real atomVolume = probeScale * probeScale * + (x * dots.probeX + y * dots.probeY + z * dots.probeZ + probeScale * dots.probeCount) * + kVolumeScale; + surfaceArea += atomArea; + volume += atomVolume; + if ((polar[atomIdx] & polarMask) != 0) { + polarArea += atomArea; + polarVolume += atomVolume; + } + } + if constexpr (kVdw) { + vdwVolume += + radius * radius * (x * dots.vdwX + y * dots.vdwY + z * dots.vdwZ + radius * dots.vdwCount) * kVolumeScale; + } + } + const Real values[kNumDclvColumns] = { + surfaceArea, + polarArea, + volume, + vdwVolume, + polarVolume, + surfaceArea / cbrt(Real(36 * M_PI) * volume * volume), + vdwVolume / volume, + }; + writeRow(values); +} + +template +void launchDclvKernel(const DeviceCoordView& coordinates, + const Property3DDeviceInputs& inputs, + const DclvOptions& options, + Real* centered, + const DclvOutputs& outputs, + const cudaStream_t stream) { + const int numBlocks = (coordinates.numConformers + kWarpsPerBlock - 1) / kWarpsPerBlock; + dclv3DKernel<<>>(coordinates, + inputs, + static_cast(options.probeRadius), + dclvPolarMask(options), + centered, + outputs); + cudaCheckError(cudaGetLastError()); +} + +//! Launches the DCLV kernel when any DCLV property is requested. The probe-expanded dots serve every output but +//! DCLVVDWVolume; the bare-sphere dots serve DCLVVDWVolume and DCLVPackingDensity. Each is tested only when needed. +template +void launchDclvProperties(const DeviceCoordView& coordinates, + const Property3DDeviceInputs& inputs, + const DclvOptions& options, + const DclvOutputs& outputs, + const cudaStream_t stream) { + if (coordinates.numConformers == 0 || !outputs.any()) { + return; + } + const bool vdw = outputs.values[kDclvVdwVolume] != nullptr || outputs.values[kDclvPackingDensity] != nullptr; + bool probe = false; + for (int column = 0; column < kNumDclvColumns; ++column) { + probe |= column != kDclvVdwVolume && outputs.values[column] != nullptr; + } + AsyncDeviceVector centered(static_cast(coordinates.numAtoms) * 3, stream); + if (probe && vdw) { + launchDclvKernel(coordinates, inputs, options, centered.data(), outputs, stream); + } else if (probe) { + launchDclvKernel(coordinates, inputs, options, centered.data(), outputs, stream); + } else { + launchDclvKernel(coordinates, inputs, options, centered.data(), outputs, stream); + } +} + +} // namespace nvMolKit::descriptors3d_detail + +#endif // NVMOLKIT_DESCRIPTORS3D_DCLV_CUH diff --git a/src/descriptors3d_mol.cpp b/src/descriptors3d_mol.cpp index 2a34e885..3c5fc3b3 100644 --- a/src/descriptors3d_mol.cpp +++ b/src/descriptors3d_mol.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -34,6 +35,8 @@ struct DeviceDescriptorInputs { AsyncDeviceVector bondNeighbors; AsyncDeviceVector usrcatAtomClasses; AsyncDeviceVector heavyAtomFlags; + AsyncDeviceVector vdwRadii; + AsyncDeviceVector dclvPolarClasses; AsyncDeviceVector defaultConformerRows; AsyncDeviceVector conformerIs3D; AsyncDeviceVector moleculeAtomStarts; @@ -48,6 +51,7 @@ struct DescriptorInputNeeds { bool bondAdjacency = false; bool usrcatAtomClasses = false; bool heavyAtomFlags = false; + bool dclvAtomInputs = false; bool defaultConformerRows = false; bool conformerFlags = false; }; @@ -73,6 +77,37 @@ void writeBondAdjacency(const RDKit::ROMol& mol, } } +//! Property3DDeviceInputs::dclvPolarClasses bits of @p atom, from RDKit's includeAsPolar (DCLV.cpp). +uint8_t dclvPolarClasses(const RDKit::ROMol& mol, const RDKit::Atom& atom) { + auto heavyClass = [](const int atomicNum) -> uint8_t { + switch (atomicNum) { + case 7: + case 8: + return kDclvPolarNitrogenOxygen; + case 15: + case 16: + return kDclvPolarSulfurPhosphorus; + default: + return 0; + } + }; + if (atom.getAtomicNum() != 1) { + return heavyClass(atom.getAtomicNum()); + } + // A hydrogen is polar when a neighbor is. RDKit recurses into hydrogen neighbors, which never terminates for + // bonded hydrogens (H2); such hydrogens are not polar here. + uint8_t classes = 0; + for (const auto* neighbor : mol.atomNeighbors(&atom)) { + const uint8_t neighborClass = heavyClass(neighbor->getAtomicNum()); + if (neighborClass == kDclvPolarNitrogenOxygen) { + classes |= kDclvPolarHydrogenOnNO; + } else if (neighborClass == kDclvPolarSulfurPhosphorus) { + classes |= kDclvPolarHydrogenOnSulfurPhos; + } + } + return classes; +} + //! RDKit's USRCAT atom classes (hydrophobic, aromatic, acceptor, donor), from //! Code/GraphMol/Descriptors/USRDescriptor.cpp. const std::array, 4>& usrcatClassPatterns() { @@ -141,6 +176,8 @@ DeviceDescriptorInputs uploadDescriptorInputs(const std::vector usrcatAtomClasses(needs.usrcatAtomClasses ? static_cast(totalAtoms) : 0); std::vector heavyAtomFlags(needs.heavyAtomFlags ? static_cast(totalAtoms) : 0); + std::vector vdwRadii(needs.dclvAtomInputs ? static_cast(totalAtoms) : 0); + std::vector polarClasses(needs.dclvAtomInputs ? static_cast(totalAtoms) : 0); // Coordinate rows follow each molecule's conformers in order, so a molecule's first row is its default // conformer (RDKit's getConformer(-1)). std::vector defaultConformerRows; @@ -178,7 +215,7 @@ DeviceDescriptorInputs uploadDescriptorInputs(const std::vectorgetIdx()] = atom->getAtomicNum() > 1 ? 1 : 0; } } + if (needs.dclvAtomInputs) { + const RDKit::PeriodicTable* table = RDKit::PeriodicTable::getTable(); + for (const auto* atom : mol.atoms()) { + vdwRadii[atomStart + atom->getIdx()] = table->getRvdw(atom->getAtomicNum()); + polarClasses[atomStart + atom->getIdx()] = dclvPolarClasses(mol, *atom); + } + } if (needs.bondAdjacency) { writeBondAdjacency(mol, neighborOffsets[molIdx], bondNeighborStarts.data() + atomStart, bondNeighbors.data()); } @@ -257,6 +303,8 @@ DeviceDescriptorInputs uploadDescriptorInputs(const std::vector(bondNeighbors.size(), stream), AsyncDeviceVector(usrcatAtomClasses.size(), stream), AsyncDeviceVector(heavyAtomFlags.size(), stream), + AsyncDeviceVector(vdwRadii.size(), stream), + AsyncDeviceVector(polarClasses.size(), stream), AsyncDeviceVector(defaultConformerRows.size(), stream), AsyncDeviceVector(conformerIs3D.size(), stream), AsyncDeviceVector(atomStarts.size(), stream)}; @@ -284,6 +332,10 @@ DeviceDescriptorInputs uploadDescriptorInputs(const std::vector calc3DProperties(const std::vector calc3DProperties(const std::vector #include +#include #include #include #include @@ -13,6 +14,7 @@ #include #include +#include #include #include #include @@ -46,8 +48,23 @@ constexpr std::array kMomentProperties = { Property3D::SpherocityIndex, }; -std::unique_ptr molWithConformers(const char* smiles, const std::vector>& confs) { - std::unique_ptr mol(RDKit::SmilesToMol(smiles)); +//! Every DCLV property, in RDKit's DoubleCubicLatticeVolume getter order. +constexpr std::array kDclvProperties = { + Property3D::DCLVSurfaceArea, + Property3D::DCLVPolarSurfaceArea, + Property3D::DCLVVolume, + Property3D::DCLVVDWVolume, + Property3D::DCLVPolarVolume, + Property3D::DCLVCompactness, + Property3D::DCLVPackingDensity, +}; + +std::unique_ptr molWithConformers(const char* smiles, + const std::vector>& confs, + const bool keepHs = false) { + RDKit::SmilesParserParams params; + params.removeHs = !keepHs; + std::unique_ptr mol(RDKit::SmilesToMol(smiles, params)); for (const auto& points : confs) { auto conf = std::make_unique(mol->getNumAtoms()); for (size_t atomIdx = 0; atomIdx < points.size(); ++atomIdx) { @@ -389,6 +406,65 @@ TEST(Descriptors3DGetaway, MatchesRdkitGetaway) { } } +TEST(Descriptors3DDclv, MatchesRdkitDoubleCubicLatticeVolume) { + // Glycine with explicit hydrogens on N and O, so every polar class and option is exercised; two conformers, + // the second far from the origin. + auto mol = molWithConformers("[H]OC(=O)CN([H])[H]", + { + { {-1.9, 0.9, 0.1}, + {-1.2, 0.2, 0.0}, + {0.1, 0.6, 0.0}, + {0.4, 1.8, 0.1}, + {1.1, -0.5, -0.1}, + {2.5, -0.1, 0.0}, + {2.8, 0.5, 0.8}, + {3.0, -1.0, 0.1} }, + {{98.1, -50.9, 30.4}, + {98.8, -51.6, 30.2}, + {100.1, -51.1, 30.1}, + {100.5, -50.0, 30.5}, + {101.0, -52.2, 29.6}, + {102.3, -51.6, 29.4}, + {102.9, -52.3, 29.0}, + {102.4, -50.9, 28.6}}, + }, + true); + const std::vector mols = {mol.get()}; + for (const double probeRadius : {1.4, 0.0}) { + for (const bool includeSandP : {false, true}) { + for (const bool includeHs : {false, true}) { + nvMolKit::Property3DOptions options; + options.dclv = {probeRadius, includeSandP, includeHs}; + auto results = + nvMolKit::calc3DProperties(mols, {kDclvProperties.begin(), kDclvProperties.end()}, options, nullptr); + for (int confIdx = 0; confIdx < 2; ++confIdx) { + RDKit::Descriptors::DoubleCubicLatticeVolume dclv(*mol, false, true, probeRadius, confIdx); + const std::array expected = { + dclv.getSurfaceArea(), + dclv.getPolarSurfaceArea(includeSandP, includeHs), + dclv.getVolume(), + dclv.getVDWVolume(), + dclv.getPolarVolume(includeSandP, includeHs), + dclv.getCompactness(), + dclv.getPackingDensity(), + }; + for (size_t propertyIdx = 0; propertyIdx < kDclvProperties.size(); ++propertyIdx) { + const auto values = toHost(results.properties.at(kDclvProperties[propertyIdx])); + ASSERT_EQ(values.size(), 2u); + EXPECT_NEAR(values[confIdx], expected[propertyIdx], 1e-9 * std::max(1.0, std::abs(expected[propertyIdx]))) + << "probe " << probeRadius << ", S/P " << includeSandP << ", Hs " << includeHs << ", conformer " + << confIdx << ", " << nvMolKit::property3DName(kDclvProperties[propertyIdx]); + } + } + } + } + } + nvMolKit::Property3DOptions invalid; + invalid.dclv.probeRadius = -0.5; + EXPECT_THROW(nvMolKit::calc3DProperties(mols, {Property3D::DCLVVolume}, invalid, nullptr), + std::invalid_argument); +} + TEST(Descriptors3DPairwise, BondSearchFallbackMatchesBondDistanceTable) { // A 12-atom chain has bond distances up to 11, past the deepest lag either descriptor uses. constexpr int kChainAtoms = 12; From ec3a0cdc0d81a976a8b42e2869b3330390550ba2 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Mon, 19 Oct 2026 13:35:26 -0400 Subject: [PATCH 2/5] Validate DCLV stream lifetime and simplify device output writes --- nvmolkit/tests/test_descriptors3d.py | 23 ++++++++ src/descriptors3d_dclv.cuh | 87 ++++++++++++++++------------ 2 files changed, 73 insertions(+), 37 deletions(-) diff --git a/nvmolkit/tests/test_descriptors3d.py b/nvmolkit/tests/test_descriptors3d.py index cea1f792..32f330e5 100644 --- a/nvmolkit/tests/test_descriptors3d.py +++ b/nvmolkit/tests/test_descriptors3d.py @@ -790,6 +790,29 @@ def test_dclv_device_coordinates_far_from_origin(precision): _assert_dclv_matches(far, near, precision) +@pytest.mark.parametrize("precision", PRECISIONS) +def test_dclv_stream_ordering_empty_inputs_and_mixed_selection(precision): + mols = [_embed("CCO", 3, 137), Chem.MolFromSmiles("CC"), _embed("CCN", 2, 139)] + stream = torch.cuda.Stream() + with torch.cuda.stream(stream): + coordinates = _device_result_from_molecules(mols) + result = Calc3DProperties( + mols, (*DCLV_PROPERTIES, Property3D.PMI1), coordinates=coordinates, precision=precision, stream=stream + ) + values = [result[prop] for prop in DCLV_PROPERTIES] + moments = result[Property3D.PMI1] + labels = (result.mol_indices, result.conf_indices) + del result, coordinates + gc.collect() + _assert_dclv_matches(np.column_stack([value.numpy() for value in values]), _rdkit_dclv_rows(mols), precision) + _assert_matches_rdkit(moments.numpy(), _reference_rows(mols, (Property3D.PMI1,), True)[:, 0], precision) + np.testing.assert_array_equal(labels[0].numpy(), [0, 0, 0, 2, 2]) + np.testing.assert_array_equal(labels[1].numpy(), [0, 1, 2, 0, 1]) + for empty_mols in ([], [Chem.MolFromSmiles("CC")], [Chem.MolFromSmiles("")]): + empty = Calc3DProperties(empty_mols, DCLV_PROPERTIES, precision=precision, stream=stream) + assert all(empty[prop].torch().shape == (0,) for prop in DCLV_PROPERTIES) + + def test_dclv_rejects_invalid_probe_radius_only_when_requested(): mol = _embed("CCO", 1, 5) for probe_radius in (-0.1, float("nan"), float("inf")): diff --git a/src/descriptors3d_dclv.cuh b/src/descriptors3d_dclv.cuh index 7b0af148..ea637084 100644 --- a/src/descriptors3d_dclv.cuh +++ b/src/descriptors3d_dclv.cuh @@ -509,26 +509,33 @@ __device__ __forceinline__ DclvAtomDots dclvAtomDots(const Real* c const Real vdwDotZ = uz * radius; bool probeOpen = kProbe; bool vdwOpen = kVdw; - auto test = [&](const DclvNeighbor& other) { - if constexpr (kProbe) { - probeOpen = - probeOpen && !dclvCovers(other.x, other.y, other.z, other.radius + probe, probeDotX, probeDotY, probeDotZ); - } - if constexpr (kVdw) { - vdwOpen = vdwOpen && !dclvCovers(other.x, other.y, other.z, other.radius, vdwDotX, vdwDotY, vdwDotZ); - } - }; if (cached) { for (int neighborIdx = 0; neighborIdx < numNeighbors && (probeOpen || vdwOpen); ++neighborIdx) { - test(neighbors[neighborIdx]); + const auto& neighbor = neighbors[neighborIdx]; + if constexpr (kProbe) { + probeOpen = + probeOpen && + !dclvCovers(neighbor.x, neighbor.y, neighbor.z, neighbor.radius + probe, probeDotX, probeDotY, probeDotZ); + } + if constexpr (kVdw) { + vdwOpen = + vdwOpen && !dclvCovers(neighbor.x, neighbor.y, neighbor.z, neighbor.radius, vdwDotX, vdwDotY, vdwDotZ); + } } } else { for (int other = 0; other < numAtoms && (probeOpen || vdwOpen); ++other) { if (other != atomIdx) { - test(DclvNeighbor{centered[other * 3 + 0] - px, - centered[other * 3 + 1] - py, - centered[other * 3 + 2] - pz, - static_cast(radii[other])}); + const Real deltaX = centered[other * 3 + 0] - px; + const Real deltaY = centered[other * 3 + 1] - py; + const Real deltaZ = centered[other * 3 + 2] - pz; + const Real otherRadius = static_cast(radii[other]); + if constexpr (kProbe) { + probeOpen = + probeOpen && !dclvCovers(deltaX, deltaY, deltaZ, otherRadius + probe, probeDotX, probeDotY, probeDotZ); + } + if constexpr (kVdw) { + vdwOpen = vdwOpen && !dclvCovers(deltaX, deltaY, deltaZ, otherRadius, vdwDotX, vdwDotY, vdwDotZ); + } } } } @@ -585,15 +592,6 @@ __global__ void dclv3DKernel(const DeviceCoordView coordinates, if (conformerIdx >= coordinates.numConformers) { return; // Uniform across the warp. } - auto writeRow = [&](const Real(&values)[kNumDclvColumns]) { - if (lane == 0) { - for (int column = 0; column < kNumDclvColumns; ++column) { - if (outputs.values[column] != nullptr) { - outputs.values[column][conformerIdx] = values[column]; - } - } - } - }; const ConformerAtoms atoms = loadConformer(coordinates, nullptr, inputs.moleculeAtomStarts, conformerIdx); bool valid = atoms.valid; const double* radii = nullptr; @@ -625,11 +623,13 @@ __global__ void dclv3DKernel(const DeviceCoordView coordinates, centreZ = groupAllReduceSum(sumZ) * inverseAtoms; } if (!valid) { - Real nanRow[kNumDclvColumns]; - for (Real& value : nanRow) { - value = static_cast(nan("")); + if (lane == 0) { + for (int column = 0; column < kNumDclvColumns; ++column) { + if (outputs.values[column] != nullptr) { + outputs.values[column][conformerIdx] = static_cast(nan("")); + } + } } - writeRow(nanRow); return; } @@ -683,16 +683,29 @@ __global__ void dclv3DKernel(const DeviceCoordView coordinates, radius * radius * (x * dots.vdwX + y * dots.vdwY + z * dots.vdwZ + radius * dots.vdwCount) * kVolumeScale; } } - const Real values[kNumDclvColumns] = { - surfaceArea, - polarArea, - volume, - vdwVolume, - polarVolume, - surfaceArea / cbrt(Real(36 * M_PI) * volume * volume), - vdwVolume / volume, - }; - writeRow(values); + if (lane == 0) { + if (outputs.values[kDclvSurfaceArea] != nullptr) { + outputs.values[kDclvSurfaceArea][conformerIdx] = surfaceArea; + } + if (outputs.values[kDclvPolarSurfaceArea] != nullptr) { + outputs.values[kDclvPolarSurfaceArea][conformerIdx] = polarArea; + } + if (outputs.values[kDclvVolume] != nullptr) { + outputs.values[kDclvVolume][conformerIdx] = volume; + } + if (outputs.values[kDclvVdwVolume] != nullptr) { + outputs.values[kDclvVdwVolume][conformerIdx] = vdwVolume; + } + if (outputs.values[kDclvPolarVolume] != nullptr) { + outputs.values[kDclvPolarVolume][conformerIdx] = polarVolume; + } + if (outputs.values[kDclvCompactness] != nullptr) { + outputs.values[kDclvCompactness][conformerIdx] = surfaceArea / cbrt(Real(36 * M_PI) * volume * volume); + } + if (outputs.values[kDclvPackingDensity] != nullptr) { + outputs.values[kDclvPackingDensity][conformerIdx] = vdwVolume / volume; + } + } } template From 29638c8d3e962edc79db17d31cb3d7e7062adf5c Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Thu, 22 Oct 2026 16:16:58 -0400 Subject: [PATCH 3/5] Supply explicit DCLV reference radii for RDKit 2025.09.1 --- tests/test_descriptors3d.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/tests/test_descriptors3d.cpp b/tests/test_descriptors3d.cpp index c3e14e41..9c958954 100644 --- a/tests/test_descriptors3d.cpp +++ b/tests/test_descriptors3d.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -430,6 +431,11 @@ TEST(Descriptors3DDclv, MatchesRdkitDoubleCubicLatticeVolume) { }, true); const std::vector mols = {mol.get()}; + // RDKit 2025.09.1 requires explicit radii; empty-vector defaults were added later. + std::vector radii; + for (const auto atom : mol->atoms()) { + radii.push_back(RDKit::PeriodicTable::getTable()->getRvdw(atom->getAtomicNum())); + } for (const double probeRadius : {1.4, 0.0}) { for (const bool includeSandP : {false, true}) { for (const bool includeHs : {false, true}) { @@ -438,12 +444,7 @@ TEST(Descriptors3DDclv, MatchesRdkitDoubleCubicLatticeVolume) { auto results = nvMolKit::calc3DProperties(mols, {kDclvProperties.begin(), kDclvProperties.end()}, options, nullptr); for (int confIdx = 0; confIdx < 2; ++confIdx) { - RDKit::Descriptors::DoubleCubicLatticeVolume dclv(*mol, - std::vector{}, - false, - true, - probeRadius, - confIdx); + RDKit::Descriptors::DoubleCubicLatticeVolume dclv(*mol, radii, false, true, probeRadius, confIdx); const std::array expected = { dclv.getSurfaceArea(), dclv.getPolarSurfaceArea(includeSandP, includeHs), From 663f08495650850936283af20eec7e69349345d5 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Thu, 22 Oct 2026 16:20:51 -0400 Subject: [PATCH 4/5] Clarify DCLV property documentation --- nvmolkit/descriptors3d.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/nvmolkit/descriptors3d.py b/nvmolkit/descriptors3d.py index d104b547..3c7f0347 100644 --- a/nvmolkit/descriptors3d.py +++ b/nvmolkit/descriptors3d.py @@ -27,9 +27,9 @@ class Property3D(Enum): ``rdMolDescriptors`` module (``CalcWHIM``, ``GetUSR``, ...). Wherever a property is accepted, its string value may be used instead. - The ``DCLV`` properties are the getters of RDKit's - ``DoubleCubicLatticeVolume`` class (``DCLVSurfaceArea`` is - ``GetSurfaceArea``, ...), configured by :class:`DclvOptions`. Requesting + 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. From 4d8172bdde5aef6d0889c0a947eea71ca8f57db5 Mon Sep 17 00:00:00 2001 From: Kevin Boyd Date: Thu, 22 Oct 2026 22:26:35 +0000 Subject: [PATCH 5/5] Pass explicit DCLV radii in the Python reference and benchmark RDKit 2025.09.1's Python DoubleCubicLatticeVolume has no default-radii constructor. The Python test reference, descriptors3d_bench.py and its tests now pass per-atom van der Waals radii from the periodic table, which match the defaults of later RDKit releases exactly. Co-Authored-By: Claude Opus 5.5 --- benchmarks/descriptors3d_bench.py | 8 +++++++- benchmarks/tests/test_descriptors3d_dclv_bench.py | 14 +++++++++----- nvmolkit/tests/test_descriptors3d.py | 7 ++++++- 3 files changed, 22 insertions(+), 7 deletions(-) diff --git a/benchmarks/descriptors3d_bench.py b/benchmarks/descriptors3d_bench.py index 74328235..77fc6afb 100644 --- a/benchmarks/descriptors3d_bench.py +++ b/benchmarks/descriptors3d_bench.py @@ -139,11 +139,17 @@ def _pack_device_coordinates(mols: list[Chem.Mol]) -> Device3DResult: ) +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, confId=conf_id) + rdMolDescriptors.DoubleCubicLatticeVolume(mol, _dclv_radii(mol), confId=conf_id) if any(prop in DCLV_GETTERS for prop in properties) else None ) diff --git a/benchmarks/tests/test_descriptors3d_dclv_bench.py b/benchmarks/tests/test_descriptors3d_dclv_bench.py index 72b05f02..09a9f073 100644 --- a/benchmarks/tests/test_descriptors3d_dclv_bench.py +++ b/benchmarks/tests/test_descriptors3d_dclv_bench.py @@ -44,16 +44,20 @@ def test_dclv_reference_reuses_one_construction_per_conformer_and_timed_call(mon for conformer in molecule.GetConformers(): for prop in expected: value = ( - getattr(constructor(molecule, confId=conformer.GetId()), bench.DCLV_GETTERS[prop])() + 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, **kwargs): + def counted_constructor(molecule, radii, **kwargs): + assert radii == bench._dclv_radii(molecule) constructions.append((molecule, kwargs["confId"])) - return constructor(molecule, **kwargs) + return constructor(molecule, radii, **kwargs) monkeypatch.setattr(bench.rdMolDescriptors, "DoubleCubicLatticeVolume", counted_constructor) expected_constructions = [ @@ -129,10 +133,10 @@ def counted_property(molecule, conf_id, prop): return calculator(molecule, conf_id, prop) class CountedDclv: - def __init__(self, molecule, *, confId): + def __init__(self, molecule, radii, *, confId): constructions.append(confId) self.conf_id = confId - self.reference = constructor(molecule, confId=confId) + self.reference = constructor(molecule, radii, confId=confId) def __getattr__(self, name): def read_property(): diff --git a/nvmolkit/tests/test_descriptors3d.py b/nvmolkit/tests/test_descriptors3d.py index 32f330e5..17319dfa 100644 --- a/nvmolkit/tests/test_descriptors3d.py +++ b/nvmolkit/tests/test_descriptors3d.py @@ -676,10 +676,15 @@ def test_getaway_precision_ignored_when_not_requested(precision_digits): def _rdkit_dclv_rows(mols, options=None): """RDKit's DoubleCubicLatticeVolume getters, one column per DCLV_PROPERTIES entry.""" options = DclvOptions() if options is None else options + table = Chem.GetPeriodicTable() rows = [] for mol in mols: + # RDKit 2025.09.1 has no default-radii constructor, so pass the van der Waals radii it would use. + radii = [table.GetRvdw(atom.GetAtomicNum()) for atom in mol.GetAtoms()] for conf in mol.GetConformers(): - dclv = rdMolDescriptors.DoubleCubicLatticeVolume(mol, probeRadius=options.probeRadius, confId=conf.GetId()) + dclv = rdMolDescriptors.DoubleCubicLatticeVolume( + mol, radii, probeRadius=options.probeRadius, confId=conf.GetId() + ) polar = {"includeSandP": options.includeSandP, "includeHs": options.includeHs} rows.append( [