Computing the Volume, Surface Area, Mean, and Gaussian Curvatures of Molecules and Their Derivatives

Geometry is crucial in our efforts to comprehend the structures and dynamics of biomolecules. For example, volume, surface area, and integrated mean and Gaussian curvature of the union of balls representing a molecule are used to quantify its interactions with the water surrounding it in the morphometric implicit solvent models. The Alpha Shape theory provides an accurate and reliable method for computing these geometric measures. In this paper, we derive homogeneous formulas for the expressions of these measures and their derivatives with respect to the atomic coordinates, and we provide algorithms that implement them into a new software package, AlphaMol. The only variables in these formulas are the interatomic distances, making them insensitive to translations and rotations. AlphaMol includes a sequential algorithm and a parallel algorithm. In the parallel version, we partition the atoms of the molecule of interest into 3D rectangular blocks, using a kd-tree algorithm. We then apply the sequential algorithm of AlphaMol to each block, augmented by a buffer zone to account for atoms whose ball representations may partially cover the block. The current parallel version of AlphaMol leads to a 20-fold speed-up compared to an independent serial implementation when using 32 processors. For instance, it takes 31 s to compute the geometric measures and derivatives of each atom in a viral capsid with more than 26 million atoms on 32 Intel processors running at 2.7 GHz. The presence of the buffer zones, however, leads to redundant computations, which ultimately limit the impact of using multiple processors. AlphaMol is available as an OpenSource software.

A. Surface Areas and Volumes of Ball Intersections, and Derivatives Several formulas have been presented for the volume and surface areas of the intersection of two, three, and four spheres with unequal radii; see for example. [1][2][3] Here we describe versions of these formulas that depend only on the radii of the spheres and the distance between their centers.
Derivations of these versions were originally given in Mach and Koehl, 2011;4 they are provided here for sake of completeness, as well as a support to provide the derivatives of these geometric measures with respect to the distance between the sphere centers.
In the following, we write B i , B j , B k , B l be four balls with bounding spheres S i , S j , S k , S l , centers z i , z j , z k , z l , and radii r i , r j , r k , r l .

Intersection of two balls
The intersection between the two balls B i and B j is the union of two caps, as illustrated in Figure   3 in the main paper, which we redraw here as it will be referred to many times below: The inter- section of the two caps is a disk with center y ij and radius r ij = z i − z j . We also define the height of cap from S i , which we denote h i;j . As noted in Bryant et al. 2004, 5 the signed distance between z i and the plane that contains the disk is: Hence, Setting Proposition S.1: The intersection between two balls is illustrated in Figure 3 in the main paper.
We have: ∂A i;j ∂r ij = 2πr i λ i;j , (S.8) Intersection of three balls. B. The tetrahedron, T , used to measure the intersection of the three balls. z i , z j , z k are the centers of the spheres, and P i;jk is one of the two points common to all three spheres.
The contribution of B i to the volume and surface area of the intersection of B i , B j , B k is defined by the intersection of the caps C i;j and C i;k , illustrated in orange and blue in panel A of Figure S.2.
The three spheres, S i , S j , S k , intersect in two points, P i;jk and P i;kj . We consider the tetrahedron, The dihedral angles at the edges z i z j and z i z k are denoted θ ij;k and θ ik;j , respectively, while ψ i is the dihedral angle at the edge z i P ijk .
Proposition S.2: The contributions of S i and B i to the surface area and volume of the triple intersection are where the dihedral angles are computed from the edge lengths of the tetrahedron T ; see section D below. Formulas for the contributions of B j and B k to the intersection are easily deduced by index permutation on these equations. The proof of Proposition A.2 can be found in Appendix A of Mach and Koehl, 2011 4 (Lemma 2).
Formulas for the derivatives with respect to edge lengths of the terms A i;jk and V i;jk are straightforward from their analytical expressions. We express them here for A i;jk : Proposition S.3: with similar expressions for the derivatives of V i;jk . Note that those derivatives require that that the corresponding derivatives of the dihedral angles they include are known. Computations of those derivatives are provided in section D below.

Intersection of four balls
The formula for the weighted area and volume of the contribution of a ball to the union of balls removes all contributions of the intersections of four balls, with the exception of a term vol F i;jkl , which we describe here. Let B i , B j , B k , B l have a non-empty common intersection. Their centers define a tetrahedron, T , with faces T i , T j , T k , T l , defined such that z a / ∈ T a for all a = i, j, k, l. We denote the dihedral angle of T between the two faces that share the edge z i z j by φ ij;kl . Let F i;jkl be the region delimited by the tetrahedron T and the three Voronoi planes that separate B i from Proposition S.4: The volume of F i;jkl is vol F i;jkl = 1 6 (r i − h i;j )r 2 ij 2 cos θ ij;k cos θ ij;l − (cos 2 θ ij;k + cos 2 θ ij;l ) cos φ ij;kl sin φ ij;kl + 1 6 (r i − h i;k )r 2 ik 2 cos θ ik;j cos θ ik;l − (cos 2 θ ik;j + cos 2 θ ik;l ) cos φ ik;jl sin φ ik;jl Derivatives of γ i , γ ij , and γ ijk The angular coefficient, γ i , of a vertex z i is computed over all tetrahedra of K that contain i. If z i is such that it belongs to at least one tetrahedron of K that also contains z a and z b , then In all other cases, ∂γ i ∂r ab = 0. Similarly, if τ ijab ∈ K, and 0 otherwise. The derivatives of the dihedral angles of a tetrahedron with respect to its edge lengths are given in section D below. The derivatives of γ ijk are piecewise zero because the γ ijk are piecewise constant. Their values change at non-generic states, where their derivatives are not defined. 5,6 Recall that the mean curvature is the sum of two terms: the contribution of the spherical patches and of the accessible circular arcs at the intersections of pairs of spheres. Along these arcs, the mean curvature is partitioned equally between the two spheres involved. The contribution of a ball i to the total mean curvature of a union of balls is then in which A i is the contribution of B i to the total surface area of the union of balls. Expressions for A i and its derivatives have been given above. Here we are concerned with the three additional terms: r i;j , α ij , and σ ij . Recall that r ij is the radius of the circle S ij at which S i and S j intersect; see

Derivative of α ij
The angle between the unit normals of S i and S j at any point P on the circle S ij is denoted α ij .
Using the law of cosine in the triangle ∆z i z j P , we get A straightforward differentiation of Equation (S.20) gives the derivatives of α ij with respect to the distance, r ij , between the centers z i and z j of S i and S j : Let S i and S j be two intersecting spheres, and let S ij at which they intersect. σ ij is the fraction of the length of S ij that is accessible, i.e. not covered by any other spheres. In references, 5,7 Edelsbrunner and co-workers developed extrinsic formulas for σ ij and its derivatives, which are based on the Cartesian coordinates of the centers of the spheres that intersect S ij . Here we revisit this problems and derive intrinsic formulas based on the distances between the sphere centers. We first establish Proposition S.5: The fraction of S ij that is not covered is given by in which the angles θ ij;k and φ ij;kl are the dihedral angles at the edge z i z j in the tetrahedra T 3 = z i z j z k P ijk and T 4 = z i z j z k z l (see section above for details), and γ ijk is defined as The circle of intersection of two spheres. We consider two spheres, S i and S j , with nonempty intersection. The Voronoi plane between their centers cuts the two spheres in a circle with center y ij and radius r ij . In A, this circle is partially covered by a third sphere, S k , that intersects the circle at two points P ijk and P ijkl . In B, the circle is partially covered by two sphere S k and S l . We assume that the corresponding caps between P ijk and P ikj for sphere S k , and between P ijl and P ilj for sphere S l , are not nested. This is always the case if the spheres arise from triangles in the dual complex.
The proof follows from the inclusion-exclusion principle. We illustrate it here for two, three, and four spheres. In the simplest case of two intersecting spheres, the entire length of the circle is accessible, and σ ij = 1. When three spheres, S i , S j , S k , intersect in an accessible point, then they form a triangle in the dual complex, K, and S ij is partially covered by S k ; see Figure S.3A. The arc S ij that is covered connects the points P ijk and P ikj , and its length is r ij · 2θ ij;k . In this case, When four spheres, S i , S j , S k , S l intersect in such a way that they form a tetrahedron in the dual complex, K, then S ij is partially covered by both S k and S l ; see Figure S.3B. The length of S ij that is covered is the sum of arc P ijk P ikj , whose length is r ij · 2θ ij;k , and the arc P ijl P ilj , whose length is r ij · 2θ ij;l , minus the sub arc common to those two arcs, whose length is r ij · x. Notice that x = θ ij;k + θ ij;l − φ ij;kl . Therefore The extensions of the three cases to include all simplices of the dual complex that contain the edge z i z j leads to Equation (S.22). Formulas for the derivatives of σ ij with respect to edge lengths are then straightforward: in which the derivatives of the dihedral angles of a tetrahedron as a function of edge lengths are provided in section D below.

C. Integrated Gaussian Curvature, and Derivative
Recall that the Gaussian curvature is the sum of three terms that account for the spherical patches, the circular arcs between spheres, and corners on the boundary of the union of balls: All variables have been defined above, with the exception of λ ij and σ i;jk , which we describe below.

Derivative of λ
The derivative of λ ij with respect to r ij , the distance between z i and z j , is therefore Derivative of σ i;jk σ i;jk is the corner term that is specific to the Gaussian curvature. Such a corner exists for each accessible triangle of the dual complex, and there are two corners for each triangle that is accessible from both sides. Let us consider one such triangle, with vertices z i , z j , z k , which are the centers of the spheres S i , S j , S k . These three spheres intersect at two corners, P ijk and P ikj (see We consider three intersecting spheres, S i , S j , S k , and one of the point, P ijk , common to the three spheres. A) The outward unit normals, n i , n j , n k , at this point form a spherical triangle on the unit sphere whose area is the contribution of P ijk to the Gaussian curvature. To divide this contribution, we use the center, c, of the unique circle that passes through n i , n j , n k . B) The spherical triangle is divided into three quadrangles by connecting c to the midpoints of the three edges. The area σ i;jk is the contribution to the Gaussian curvature we associate with S i . C) Note that c may be outside of the spherical triangle n i n j n k . In the example shown, the oriented area, A i , of n j cn k is negative. D) In the limiting case, c lies on an edge of the spherical triangle, here n j n k , which implies A i = 0. This special case leads to a singularity that needs to be handled explicitly (see text for details).
all these derivatives will be expressed as functions of the inter-vertex distances, which in this case are r ij , r jk , r ik . These formulas have been derived in Akopyan and Edelsbrunner, 2020. 8 Fix P ijk and let n i , n j , n k be the unit outward normals of the three spheres at P ijk . The total contribution of this corner to the Gaussian curvature is equal to the area σ ijk of the spherical triangle with vertices n i , n j , n k on the unit sphere; see Figure S.4. The geodesic lengths of the sides are α ij , α jk , α ik , as defined in Figure S.1. We begin with the following Proposition S.6: Let T be a spherical triangle with side lengths α, β, γ, and set a = cos 2 (α/2), b = cos 2 (β/2), c = cos 2 (γ ik /2). The area of this triangle is In addition, if R(a, b, c) is the radius of the circumcircle and r(a, b, c) = cos 2 (R(a, b, c)/2), then  S(a, b, c), is the total contribution of the corner P ijk to the Gaussian curvature. Since S i , S j , S k have different weights, we break down this contribution to individual contributions of the spheres: σ ijk = σ i;jk + σ + j : ki + σ k;ij , in which σ i;jk = ω i σ ijk , etc. The division is based on the position of the spherical circumcenter, c, of n i , n j , n k ; see S.4A for details. Note that the coefficients ω i , ω j , ω k can be seen as spherical barycentric coordinates of c with respect to the spherical triangle n i n j n k . To compute these coordinates, we split n i n j n k in two different ways. First, the triangle is split into three triangles, cn j n k , n i cn k , n i cn k , with surface areas A i , A j , A k . Let R(a, b, c) is the spherical radius of the circumcircle of n i , n j , n k , and with r(a, b, c) = cos 2 (R(a, b, c)/2) we get Second, we let m i , m j , m k be the midpoints of z j z k , z i z k , z i z j , and we subdivide n i n j n k into three quadrangles, n i m k cm j , n j m k cm i , n k m j cm i , with areas σ i;jk , σ j;ik , and σ k;ij , respectively. To establish the correspondence between the areas A and the areas σ, we need to take into account the possibility that c falls outside of the triangle n i n j n k . In the case illustrated in Figure  when sin 2 (α ij /2) + sin 2 (α ik /2) = sin 2 (α jk /2) or, equivalently, when a + c = 1 + b (see Akopyan and Edelsbrunner, 2020 8 for details). When a + c ≤ 1 + b, n i and c lie on the same side of n j n k .
We define

A singularity when computing the derivatives
To compute the derivatives of the different surface areas, σ i;jk , σ j;ik , σ k;ij , we need the derivatives of the areas A i , A j , A k . These terms are all computed as surface areas of spherical triangles, given by Proposition S.6. We recall Equation 21 in Akopyan and Edelsbrunner, 2020: 8 Using Proposition S.6, we rewrite it as Difficulties arise when the center, c, of the circle that passes through n i , n j , n k lies on one of the edges of the spherical triangle; see Figure S.4D. Let us assume, for example, that c lies on the edge n j n k . Then A i = 0. According to Equation (S.42), however, with similar expressions for the derivatives with respect to a and b.

D. The Geometry of a Tetrahedron
Let us consider the tetrahedron, T , with vertices P 1 , P 2 , P 3 , P 4 . The four faces of T are T 1 = P 2 P 3 P 4 , T 2 = P 1 P 3 P 4 , T 3 = P 1 P 2 P 4 , and T 4 = P 1 P 2 P 3 , and we write s 1 , s 2 , s 3 , s 4 for their areas, respectively. We denote the dihedral angle between T i and T j as θ ij and the length of the edge connecting P i and P j as r ij .

Volume and Surface area
The Cayley-Menger matrix, M , associated with the tetrahedron, T , is

Dihedral angles
We use the law of cotangents of dihedral angles 4 to express the relationship between the above two determinants and the dihedral angles of the tetrahedron:

Derivative of the dihedral angles of a tetrahedron
Deriving Equation (S.50) with respect to the length, r ab , of the edge connecting the vertices P a and P b , we get −(1 + cot 2 θ ij ) ∂θ ij ∂r ab = −δ ij;ab cot θ ij r ij + 1 24