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"}