Computing the surface area and volume of molecules

Question

Here is a somewhat difficult problem I've had for a long while.

Space-filling models of molecules can be easily generated by ChemicalData[], Import[] with an appropriate file, or the new MoleculePlot3D[] functionality in version 12. For instance,

m1 = MoleculePlot3D[Entity["Chemical", "Phenacetin"], PlotTheme -> "Spacefilling"]

or in older versions,

m2 = ChemicalData["Phenacetin", "SpaceFillingMoleculePlot"]

(Before I proceed further, let me make a warning for people doing comparisons on older systems: the new system uses ångström units, while the old method uses picometer units, with the conversion $1 \ \mathrm{\unicode{x212B}} = 100 \ \mathrm{pm}$. You can see this by evaluating Charting`get3DPlotRange[] on both plots.)

Now, apart from display purposes, the space-filling model can also be used to compute the so-called van der Waals area and van der Waals volume of the molecule, where the molecule is treated as a collection of rigid spheres. A simple way to generate the associated van der Waals region would be something like

vdwmol = RegionUnion @@ (Ball @@@ Cases[Normal[m1], _Sphere, ∞]);

In version 12 (whose online version I am using), one can then use

{Quantity[Area[RegionBoundary[vdwmol]], "Angstroms"^2],
 Quantity[Volume[vdwmol], "Angstroms"^3]}
   {Quantity[209.076, ("Angstroms")^2], Quantity[170.445, ("Angstroms")^3]}

(replace "Angstroms" with "Picometers" if you used m2 instead) to compute the surface area and volume, but this facility is not entirely available in older versions. For example, the area computation will not work in version 11, and none of these methods work in the much older versions that do not have region/mesh functionality. Apart from this, having alternative methods available would be a useful way to check if one is getting reasonable results.

What are (relatively) quick and (reasonably) accurate methods to compute the van der Waals area and volume of a given molecule?

Here, "quick" can be e.g. "won't take more than $2$ minutes", and "accurate" would be $\approx 4$ digits of accuracy.

One possibility would be to perform a preliminary discretization via (Boundary)DiscretizeRegion[] before using Area[]/Volume[], but this moves the difficulty to the discretization step itself, and thus one must come up with a quick way to discretize the molecule. (Nevertheless, the (Boundary)MeshRegion[] object thus produced would be a useful object on its own, so methods for quickly producing van der Waals meshes are quite welcome.)

Otherwise, one has to use NIntegrate[]. For example, you can compute the volume like so:

NIntegrate[Boole[Or @@ Cases[Normal[m1],
                             Sphere[c_, r_] :> Norm[{x, y, z} - c] < r, ∞]],
           {x, -∞, ∞}, {y, -∞, ∞}, {z, -∞, ∞}]
   170.29

but computing the surface area with NIntegrate[] is not terribly straightforward.

Of course, other methods for computing the van der Waals area and volume based on the output of MoleculePlot3D[] (or equivalent functions) would be welcome.

If one wants to test a proposed method on molecules of varying complexity, I would suggest the following ChemicalData[] entries: {"Ethylene", "Adamantane", "Mirex", "Phenacetin", "Caffeine", "Estrone", "Erythromycin"}

Answer 1

I can offer a solution by discretization to ElementMesh with external meshing software Gmsh. Currently it produces better quality mesh on curved surfaces that built-in methods. We need packages GmshLink, ImportMesh and MeshTools.

Needs["GmshLink`"]
(*Set your path to Gmsh executable.*)
$GmshDirectory = "my_path_to_current_release\\gmsh-4.4.1-Windows64";

OP's example for molecule and extracted geometric region (RegionUnion of Balls).

m1 = MoleculePlot3D[Entity["Chemical", "Phenacetin"], PlotTheme -> "Spacefilling"]
vdwmol = RegionUnion @@ (Ball @@@ Cases[Normal[m1], _Sphere, ∞]);

Create ElementMesh object from symbolic region. One needs to play around with different options of GmshGenerator, especially MaxCellMeasure to control density of mesh and therefore accuracy of computed volume and surface. Disabling "OptimizeQuality" makes discretization procedure a little faster.

mesh = ToElementMesh[
  vdwmol,
  "ElementMeshGenerator" -> {GmshGenerator,
    "OptimizeQuality" -> True,
    "GmshAlgorithm" -> "MeshAdapt",
    "MeshOrder" -> 2,
    MaxCellMeasure -> 0.25
    }
  ]

(* ElementMesh[{{-6.29722, 5.70905}, {-3.2177, 3.75228}, {-3.6796, 
   4.2785}}, {TetrahedronElement["<" 47739 ">"]}] *)

mesh["Wireframe"[
  "MeshElement" -> "MeshElements",
  "MeshElementStyle" -> FaceForm@LightBlue
  ]]

Histogram[Flatten@mesh["Quality"], PlotRange -> {{0, 1}, Automatic}]

MeshTools package contains functions for calculation of ElementMesh volume and surface.

Needs["MeshTools`"]

Total@Flatten@MeshElementMeasure[mesh]
(* 171.261 *)

Total@Flatten@BoundaryElementMeasure[mesh]
(* 215.234 *)

This example took about 30 seconds on my computer and discretization time is heavily affected by MaxCellMeasure. One could also use other meshing software and import is result as ElementMesh object.

Answer 2

The paper Fast Calculation of van der Waals Volume as a Sum of Atomic and Bond Contributions and Its Application to Drug Compounds, pointed out by theorist gives an approximate formula for the molecular volume as

$$ V_{\mathrm{vdW}} = \sum \mathrm{all\: atom\: contributions} - 5.92 N_\mathrm{B} - 14.7 R_\mathrm{A} - 3.8 R_\mathrm{NA} $$

where $N_\mathrm{B}$ is the number of bonds, $R_\mathrm{A}$ the number of aromatic rings, and $R_\mathrm{NA}$ the number of nonaromatic rings.

vdwSurfaceArea[mol_ ? MoleculeQ] := Module[
    {atomRadii, atomContributions, nb, ra, rna},
    atomRadii = AtomList[
        mol, All, EntityProperty["Element", "VanDerWaalsRadius"]
    ];
    If[MemberQ[atomRadii, _Missing],
        Return[Missing @ "VDWRadiusMissing", Module]
    ];
    atomRadii = UnitConvert[atomRadii, "Angstroms"];
    atomContributions = N @ Total[(4 / 3) * Pi * atomRadii^3];
    nb = BondCount[mol];
    ra = mol["AromaticRingCount"];
    rna = mol["AliphaticRingCount"];
    atomContributions - Quantity[5.92 nb + 14.7 ra + 3.8 rna, "Angstroms" ^ 3];
]

The results match reasonably closely with those in the supporting info for the linked paper. For the molecule in Pinti's answer you get something very close to the mesh volume:

In[73]:= AbsoluteTiming[
 vdwSurfaceArea[Molecule["N-(4-ethoxyphenyl)acetamide"]]]

Out[73]= {0.01516, Quantity[176.291, ("Angstroms")^3]}

Answer 3

The method used in Jason's answer, as noted, requires version 12. For those in earlier versions who are unable or unwilling to install cheminformatics toolboxes like RDKit, ChEMBL Beaker (see also the associated paper) provides an API that can be used to compute the ring counts required by the underlying formula, which is used in the code that follows:

vdwVolume[s_] := Module[{mol, molString, raw},
   mol = Table[ChemicalData[s, p],
               {p, {"VertexTypes", "EdgeRules", "EdgeTypes", "AtomPositions"}}]; 
   molString =
   ExportString[Thread[{"VertexTypes", "EdgeRules", "EdgeTypes", "VertexCoordinates"} -> 
                       MapAt[QuantityMagnitude, mol, 4]], {"MOL", "Rules"}];
   raw = URLExecute["https://www.ebi.ac.uk/chembl/api/utils/descriptors", {}, 
                    "RawJSON", Method -> "POST", "Headers" -> {"Expect" -> ""}, 
                    "Body" -> molString];
   UnitConvert[Total[Times @@@ MapAt[(4/3 π ElementData[#, "VanDerWaalsRadius"]^3) &, 
                                     Tally[mol[[1]]], {All, 1}]] - 
               Quantity[{3.8, 14.7} .
                        First[Lookup[raw, {"NumAliphaticRings", "NumAromaticRings"}]] + 
                        5.92 Length[mol[[2]]], "Angstroms"^3], "Angstroms"^3]]

I have omitted error checking/trapping from this code to keep it short; for actual applications, you should add in the necessary safeguards. Additionally, in version 11, one should instead use HTTPRequest[] in URLExecute[].

For example,

QuantityMagnitude[vdwVolume /@ {"Ethylene", "Adamantane", "Mirex", "Phenacetin",
                                "Caffeine", "Estrone", "Erythromycin"}]
   {40.512, 140.647, 298.467, 176.291, 151.433, 266.593, 731.435}

(If you want to use this method directly on MOL files, you can replace the first two assignments with something like molString = Import["fileName.mol", "Text"]; mol = ImportString[molString, {"MOL", {"VertexTypes", "VertexCoordinates"}}];.)