Generating Molecular vdW Surfaces

Background

The vdW radius of an atom is the closest approachable distance by an external object before repulsion overtakes in importance. This distance taken from all directions creates a sphere called the atomic vdW surface. To create this sphere for any atom you need the equation of a sphere’s surface area $\mathrm{4πr^2}$ substituting $\mathrm{r}$ with the atom’s vdW radius. Databases of atomic vdW radii exist with slight numerical variations . In vdw-surfgen (hence called vsg), I focus on the following 38 fairly common atoms using radii from Truhlar and co. and periodictable.com, all in Angstrom (Å i.e. $\mathrm{10^{-10} m}$).

VDW_RADII = { 
    'H': 1.20, 'C': 1.70, 'N': 1.55, 'O': 1.52, 
    'F': 1.47, 'P': 1.80, 'S': 1.80, 'Cl': 1.75,
    'Br': 1.85, 'I': 1.98, 'He': 1.40, 'Ne': 1.54,
    'Ar': 1.88, 'Kr': 2.02, 'Xe': 2.16, 'Li': 1.82, 
    'Na': 2.27, 'K': 2.75, 'Rb': 3.03, 'Cs': 3.43, 
    'Fr': 3.48, 'Be': 2.00,'Mg': 1.73, 'Ca': 2.31, 
    'Sr': 2.49, 'Ba': 2.68, 'Ra': 2.83, 'B': 1.92, 
    'Al': 1.84, 'Si': 2.10, 'Ti': 2.15, 'Fe': 2.00, 
    'Zn': 2.10, 'Cu': 1.95, 'Mn': 2.05, 'Hg': 2.05,
    'Pb': 2.02, 'U': 1.86
}

When atoms are combined in a molecule, their spheres jointly create new, complex shapes. For example:

Generating Atomic vdW Surfaces

The first challenge in vsg was to generate coordinates of points across the vdW sphere for any atom given the 3D coordinates of just the atom. This is a two part challenge:

1. Generating evenly spaced coordinates (points) to form a sphere given a variable number of points.
2. Adequately covering the vdW surface for each atom with these points, so, for example, we don’t use 1000 points for Hydrogen (H) and same number of points for Lead (Pb), instead we scale the number based on vdW radius.

The Fibonacci Sphere Algorithm

Luckily, challenge (i) has been addressed on Stack Overflow but that post assumes a level of mathematical intuition I did not explicitly possess prior to writing this. This is the code:

def fibonacci_sphere(samples):
    indices = np.arange(0, samples, dtype=float) + 0.5
    phi = np.arccos(1 - 2 * indices / samples)
    theta = np.pi * (1 + 5 ** 0.5) * indices 
    z = np.cos(phi) 
    x = np.cos(theta) * np.sin(phi) 
    y = np.sin(theta) * np.sin(phi) 
    return np.stack((x, y, z), axis=1)

The samples at the start is an integer of the number of points we intend to spread to form the sphere. The larger this number, the greater the detail we capture.

If samples= 100, np.arangefor(0,100)creates a 100 item list starting from 0 and ending at 99. The addition of 0.5 to these numbers confused me, but it turns out that shifting by 0.5 avoids sphere surface points starting at the center of the pole, which overall makes the sphere more descriptive.

The second line maps the numbers from index 0.5 to the highest zero based index +0.5, to numbers between +1 and -1. For example, if samples is 100, for index 0 (which, after shifting becomes 0.5), the bracket contains $1 - 2\left(\frac{0.5}{100}\right)$ which is ~+1 and for index 99.5, this is ~-1. This way, even if we have 2 or 200,000,000 points, they are spaced evenly and mapped between +1 and -1.

From phi, we notice the utility of this mapping as arccos(+1) = 0 and arccos(-1) =180, i.e. for any number of points we can translate the lowest of them to 0 and the highest to 180 and still maintain even spacing between each point in the list.

It becomes more useful as we go down the function, but those require a discussion of how we get the coordinates of a sphere. For that, we need to discuss the unit circle.

The Unit Circle

This is a circle centered at the origin (0,0) with a radius of 1. As a point P rotates around this circle, its coordinates can be defined solely by the angle of rotation, θ, measured counter-clockwise from the positive x-axis.

To get the x coordinate of the point on the circle, we use the trigonometric function cos(θ). We do this because, for the red-black-green triangle, $\cos(\theta) = \frac{\text{adjacent}}{\text{hypotenuse}}$. Since the adjacent side represents the length of the horizontal green line (our focus), and the hypotenuse is the radius r = 1, $\cos(\theta) = \frac{\text{x}}{\text{1}}$. With this, we can find the x-coordinate as $x = r \cos(\theta) = 1 \cdot \cos(\theta)$.

Performing this same operation using the sine function gets us the $y$ as $1 \cdot \sin(\theta)$. This means, when the angle to the positive horizontal line is θ, Point P on this circle will have its coordinates $(x,y)$ at ($\cos(\theta),$$\sin(\theta)$).

For a dynamic understanding:

Upon observation of Figs. 2 and 3, we notice that for $x$ coordinates to go from +1 to -1, you need angles 0$^{\circ}$- 180$^{\circ}$ and for $y$ coordinates, you need angles 90$^{\circ}$- 270$^{\circ}.$ Because of this 90$^{\circ}$offset, sine and cosine of an angle will always be perpendicular, directly mapping the $xy$ coordinates.

Maximizing Coordinate Space

This information of mapping coordinates to circles is modifiable to include "randomness" using the irrational Fibonacci number, $\varphi$, where $\varphi = \frac{1 + \sqrt{5}}{2}$. This number when multiplied by 2π, gives you the Fibonacci angle, 137.5$^{\circ}$. We can get a physical meaning for this number by seeing the effect of this angle on points on a circle.

In Fig. 4, we see the positioning of each new point on this circle by the rotation of the previous point on the circle by this Fibonacci angle (137.5$^{\circ}$). Notice at the end all points are placed to maximize distances from all other points, this is the effect of using this ratio for coordinates on circular structures. Using multiples of the Fibonacci angle will also give the same behavior. Let’s really consider this.

Starting from 137.5$^{\circ}$, we go to 137.5 * 2 = 275$^{\circ}$, then, 137.5 * 3 = 412.5$^{\circ}$. From this point onwards, we run into the issue of the degrees being > 360$^{\circ}$, and we begin to instead use the modulus of multiples of the golden angle, for example, $412.5^\circ \hspace{1mm} \% \hspace{1mm} 360^\circ = 52.5^\circ$, 137.5*4 = 550 and $550^\circ \hspace{1mm} \% \hspace{1mm} 360^\circ = 190^\circ$, and so on.

By the way, this same behavior also happens if you use $\frac{1}{\varphi}$ to multiply $2π$ instead of $\varphi$. That multiplication gives 222.5$^\circ$ and since 360$^\circ$ - 137.5$^\circ$ = 222.5$^\circ$ , rotating by this angle also lands you in the same spot as rotating by 137.5$^{\circ}$- giving physical meaning to fraction inversion.

Thus, we have found a way to link the angle of a point from the positive x-axis to its x and y coordinates, and another way to ensure this angle can change in a way to space out resulting coordinates while staying away from previous coordinates on the circle as much as possible. Combining these two will help us build evenly spaced points into spheres, but we need a third idea.

Spheres as Multiple Circles

If you stand above a construction site and look at a translucent sphere being built from multiple circular slices, you will see this as bigger and bigger circles being added till the middle and then, smaller circles afterwards. This even growth and contraction in radius can be connected to angles 0$^\circ$ to 180$^\circ$ - phi, the list of angles we made. For this, we'll need a function that starts at 0, grows to its peak at the middle (90 $^\circ$) and contracts back to 0 (at 180 $^\circ$). Figs. 2 and 3 have this exact function, the sine. This means that if we multiply the sine of this list by the coordinates of the unit circle, we will get a new set of coordinates of expanding and contracting circles.

With this we can see that at the start and end, the circles are so tiny that they are points and they grow towards the middle (90 $^\circ$). Adding a z-coordinate, will shift this concentric circular behavior in a third dimension (up and down) for each succeeding circle, starting from one small point at a pole to the biggest at the center and down to the point at the other pole, giving us the sphere.

The Triune Sphere

At this point, we can combine the three concepts of, (i) building a circle solely based on angles, (ii) evenly spacing coordinates using the Fibonacci angle and, (iii) building spheres from many concentric circles, to create what is called the Fibonacci sphere. This is an algorithm that uses the Fibonacci angle to spread points to maximize distance from all previous points, while also expanding and contracting the size of the circle on which these points are being rotated along a z-axis.

To do this, we divide the sphere into as many circles as there are samples, place one pair of $x,y$ coordinates on the first circle, move the angle by the Fibonacci angle, then place the next point on the coordinates of this new angle but on the next circle and keep repeating this till we exhaust thesamples variable. At the end of this process, we form a well distributed spiral down the sphere, to and from the poles.

From Fig. 8, my early assertion that the greater the number ofsamples, the greater the detail captured, makes more sense. Anyway, at the end of this process we get coordinates for points that build the vdW sphere around a single coordinate - in this case the atomic coordinate.

After this, we can move on to building the molecular surfaces.

Linear Combination of Atomic Spheres

To create the molecular vdW surface, we combine the spheres from the atoms. We do this in a stepwise fashion, viz:

vdW Scaling

For each element in the molecule, we get the vdW radius and scale it by whatever scalevariable is provided. For example, for Hydrogen, vdW radius is 1.2. If we want to focus on the 2 vdW distance surface, we multiply 1.2 by 2 and use 2.4Å as the radius. This is done per element and their the vdW sphere surface areas are calculated. For example:

Detailed Sphere Generation

Here we handle challenge (ii) of the Generating Atomic vdW Surfaces section. We'll need to cover this area with the density of points provided. By default I use a density of 1000, meaning if the sphere’s area is 1$\text{\AA}^2$ we use 1000 points. For H atom in our 2.4Å radius illustration, we use 72,400 points as the samples variable. After this, we create a unit sphere with this number of points evenly spaced, using the Fibonacci sphere algorithm.

Atomic Sphere Scaling and Centering

Each sphere constructed is centered at the origin (0,0,0), but since each atom exists at it’s own coordinate, adding the atom’s coordinates to the points that make the sphere moves them to make the atom the new origin. At this point, each atom’s sphere has a radius of 1; to get the chemically correct sphere of interest per atom, we multiply the component coordinates of each point of the unit sphere by the scaled vdW radius - in our Hydrogen example, 2.4 and you get the coordinates on the atom’s scaled vdW sphere surface.

For example, let a point on the unit sphere be:

and let the Hydrogen atom be centered at:

If the scaled van der Waals radius is $r_{\text{vdW, scaled}} = 2.4\,\text{\AA}$, then each surface point on Hydrogen’s vdW sphere is:

Substituting the values, $\mathbf{r}_{\text{vdW, H}}$ the new atom surface specific coordinates for this point

Doing this for all atoms and surface coordinates gives us the coordinates of points that cover each of their atomic vdW spheres.

Caveat

However, in a molecule with multiple atoms, the new vdW surface does not include the internal sections of the atoms which crossover, e.g. in Fig. 1, the surface only covers the outer parts of the atoms, even though the atoms on their own would have had spheres.

Thus, for chemical meaning, we need to correctly exclude those atomic surface points which will be contained in a nearby atom’s vdW sphere when they are both in a molecule, these are hereby called occluded points.

Removing Occluded Points

To do this, we create aTruearray of same length as the total number of points for all atoms’ spheres. Then for each atom, we run this loop:

for j in range(n_atoms):
     if j == i:
       continue
       other_pos = coords_np[j]
       other_r = vdw_radii[j]
       dists = np.linalg.norm(points - other_pos, axis=1)
       keep_mask &= (dists > other_r)
   filtered_points = points[keep_mask]
   surface_points.append(filtered_points)

This gets the vdW radii and coordinates of all other atoms, then uses the coordinates to get the distances from the nearby atom to each surface point of the current atom.

E.g. If there's a Na atom at $[1.0, 1.0, 0.0]$ and we have surface points at:

points - other_pos gives:

These are difference vectors, i.e distances in each dimension. Euclidean normalization of this then gives us the nearby atom's total distance magnitude from each nearby surface point:

For each surface index in the startingkeep_mask, if the corresponding difference magnitude from the nearby atom of current focus is not > the nearby atom's vDW radius, we change corresponding index’s value in thekeep_maskto False. Therefore, if the vdw radius for the atom of current focus is 1.27, the resultantkeep_maskfor the above list, becomes,keep_mask = [ True, False, False ].

After removing the occluded points between the atom i ’s surface points for an atomj, the code moves onto the nextjand reiterates until only vDW surface points not covered or included in any other vDW regions are kept as thefiltered_pointsfor that atomi.These coordinates are then combined with asurface_pointsand returned.

Thissurface_pointslist has the final coordinates of all the spherical surface points for all the atoms, which are not contained in the internal mix of the molecular structure. While this is good enough, it is not perfect as we run into the challeng of clusters around certain regions.

This clustering is because we exclude points based on being contained in other atoms’ vdw region and not based on distance to other surface points. This meant some points technically covered same surface area, usually at the areas where neighboring sphere surfaces meet.

To avoid this, we need to rearrange the final surface points.

Optimizing Surface Points

My approach was to take the molecular surface results and space the points out while removing extraneous points.

At this point of the process, we have the list of coordinates for the points on the molecular surface and our desired density - default is 1 point per $\text\AA ^2$. Our goal is therefore to arrange the points to approximate this density while avoiding the clustering at certain regions.

This is acheived with the following greedy sampling algorithm. ^[1]

def optimize_surface_density(surface_points, density):
    if len(surface_points) == 0:
        return surface_points

    spacing = 1.0 / np.sqrt(density)
    tree = cKDTree(surface_points)
    indices = np.random.permutation(len(surface_points))
    kept = []
    remaining = set(range(len(surface_points)))

    current = indices[0]
    with tqdm(total=len(surface_points), desc="🔄 Optimizing density", unit="pts", dynamic_ncols=True, leave=False) as pbar:
        while remaining:
            kept.append(current)
            neighbors = tree.query_ball_point(surface_points[current], spacing)
            remaining -= set(neighbors)
            pbar.update(len(neighbors))
            if not remaining:
                break
            candidates = list(remaining)
            dists, _ = tree.query(surface_points[candidates], k=min(len(kept), 5))
            if len(kept) == 1:
                min_dists = dists
            else:
                min_dists = np.min(dists, axis=1) if dists.ndim > 1 else dists
            current = candidates[np.argmax(min_dists)]

    return surface_points[kept]

To understand this, let’s first consider what we mean by density, $\rho$ and translate it to a minimum distance between each potentially organized point.

If ${\rho} = \frac{N_{\text{points}}}{A_{\text{total surface}}}$, ${A_{\text{total surface}}}= \frac{N_{\text{points}}}{\rho}$ and for a single point, ${A_{\text{point covered surface}}}= \frac{1}{\rho}$. If the points are in a square, ${A_{\text{point covered surface}}}= \frac{1}{\rho}=\text{L}^2$ and from this, $\text{L} = \frac{1}{\sqrt{\rho}}$. This $\text{L}$ is used as the minimum distance between points and is unique for each $\rho$ and ${A_{\text{total surface}}}$ pair.

Armed with this minimum distance filter, we build a kDtree map of the surface coordinates we have and begin our neighbor search and filtering.

To avoid searching for neighbors from only one side of the molecule, we randomise thesurface pointlist’s indices and save in a listindices . This is then converted into aremainingset for fast removal of multiple indices.

Below, I detail the core loop in English:

Initialize a`kept`list.
		Starting from the first item in this randomized `indices` list:
		While there are points available in `remaining`:
			 Add the current point to `kept`.
			 
			 Ask the tree, “Which points are within a distance `spacing` of
			 the location `surface_points[current]`?” 
			
			 This returns indices of the neigboring points in the 
			 `surface_points`list.
			
			 Remove these indices from the`remaining`set.
			
			 If remaining is empty: exit the loop. 
			
			 Because sets are not indexed in Python, convert the current 
			 existing `remaining` set into a list of`candidates`.

			 This `candidates`is a list of remaining well distanced points from 
		   the current point in the loop. While these are adequately spaced 
		   from the current focus point, we need to check that they are also 
		   spaced from all points in `surface_points` and then select the 
		   most spaced away of the `candidates`.

We do this usingtree.query(surface_points[candidates], k=min(len(kept), 5))^[2],getting the distances of X number of points nearby, where X is whichever is lower between the length ofkeptor 5.

At the end of this, if:		 
			 
			 There's only 1 item in `kept`, `dists` will have shape (N_candidates,)
			 containine the distances for each candidate to the nearest point in
			 surface points and this will be the minimum distances list.
			 
			 However, if there are multiple items in `kept`,
			 `dists` will have shape (N_candidates,N_kept).
			 
				E.g [[1, 3, 4],   # candidate 1 → distances to 3 kept points
						 [10, 5, 7]]    # candidate 2 → distances to 3 kept points
						 
				Applying min_dists = np.min(dists, axis=1), on this picks 
				the max per row and for the above sample, we get a minimum
				distance list of [1,5], meaning the min of row 0 of `dists` 
				is the item in index 0 of `min_dists`. 
				
				Given these minimum distance lists from either instance, 
				`current = candidates[np.argmax(min_dists)]`
				get the max distance in min_dists list, uses its index 
				(since it is carried over from the candidates list) to get the 
				coordinates of that item.
				
				Saves current to `kept` and restarts the loop. 

This algorithm gives us the following vdW surface which contains less points and no clusters, a satisfactory conclusion.

That was a lot. Haha.