# Computing the surface area and volume of molecules

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 Chartingget3DPlotRange[] 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"}

• A couple of questions: (1) Is the data MMA uses to create those space filling models intended to accurately represent their VDW volumes, as opposed to merely being sufficient to create reasonable visualizations? (2) Have you compared the values generated by your code to the best available estimates of VDW volumes? This paper is a bit old, but might provide a starting point. [It provides a table showing VDW volumes calculated with their fast, simplified method (col D), and with the TSAR molecular modeling software (col E).] – theorist Sep 16 '19 at 6:25
• 1. Indeed, I am making the tacit assumption that the coordinates are "sensible" (e.g. they are the coordinates of an appropriate molecular conformer), and that the van der Waals radii are correct (cursory checking seems to indicate that they are using the classical Bondi values), so that in the abstract, this is a problem of finding the surface area and volume of a collection of balls. – J. M.'s torpor Sep 16 '19 at 12:36
• 2. I've seen that paper you refer to, but the problem is now shifted to detecting/counting rings and aromaticity in a given molecule, which does not look straightforward to do especially in older Mathematica versions. – J. M.'s torpor Sep 16 '19 at 12:36

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,
"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.

• This looks nice; I'll try this out once I get my computer back. – J. M.'s torpor Sep 16 '19 at 12:30

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[
];
];
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]}

• This is certainly a great solution for version 12; I do wonder how hard this would be to do in older versions... – J. M.'s torpor Sep 21 '19 at 5:04
• Before the version 12 functions were added, I used the CDK java library quite a bit. Here is an example using it from through JLink. And looking here I see that they do have "nAromRings" and "nSmallRings" descriptors that you could use to construct the above function. – Jason B. Sep 21 '19 at 16:08

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 =
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"}}];`.)