The Weighted Gaussian Curvature Derivative of a Space-Filling Diagram

TL;DR

This paper derives a formula for the derivative of the weighted Gaussian curvature within the morphometric approach to solvation free energy, enabling comprehensive differentiation of the energy expression based on geometric measures.

Contribution

It provides the first explicit formula for the derivative of weighted Gaussian curvature, completing the set of derivatives needed for the morphometric solvation energy model.

Findings

Derived the formula for the derivative of weighted Gaussian curvature.

Enabled full differentiation of the morphometric solvation energy expression.

Facilitated more accurate and efficient calculations in solvation modeling.

Abstract

The morphometric approach [HRC13,RHK06] writes the solvation free energy as a linear combination of weighted versions of the volume, area, mean curvature, and Gaussian curvature of the space-filling diagram. We give a formula for the derivative of the weighted Gaussian curvature. Together with the derivatives of the weighted volume in [EdKo03], the weighted area in [BEKL04], and the weighted mean curvature in [AkEd19], this yields the derivative of the morphometric expression of solvation free energy.

Tables

Table 1. Table 1: Notation for concepts, sets, functions, vectors, variables.

$S_i=\mathrm{bd}{B}_i$	sphere bounding ball
$S_{ij}=\mathrm{bd}{B}_{ij}$	circle bounding disk
$S_{ijk}=\mathrm{bd}{B}_{ijk}$	pair of points bounding line segment
$x_i,x_{ij},x_{ijk}$	centers of $S_i$, $S_{ij}$, $S_{ijk}$
$r_i,r_{ij},r_{ijk}$	radii of $S_i$, $S_{ij}$, $S_{ijk}$
${𝐮}_{ij}=\frac{x_i-x_j}{\Vert x_i-x_j\Vert }$	unit vector between centers
${𝐮}_{ijk}=\frac{{𝐮}_{ik}-⟨{𝐮}_{ik},{𝐮}_{ij}⟩{𝐮}_{ij}}{\Vert {𝐮}_{ik}-⟨{𝐮}_{ik},{𝐮}_{ij}⟩{𝐮}_{ij}\Vert }$	unit normal to ${𝐮}_{ij}$ with positive component in direction ${𝐮}_{ik}$
${\nu }_i=\frac{\mathrm{𝑉𝑜𝑙𝑢𝑚𝑒}\left(B_i\cap V_i\right)}{\mathrm{𝑉𝑜𝑙𝑢𝑚𝑒}\left(B_i\right)}$	volume fraction of ball
${\nu }_{ij}=\frac{\mathrm{𝐴𝑟𝑒𝑎}\left(B_{ij}\cap V_{ij}\right)}{\mathrm{𝐴𝑟𝑒𝑎}\left(B_{ij}\right)}$	area fraction of disk
${\nu }_{ijk}=\frac{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(B_{ijk}\cap V_{ijk}\right)}{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(B_{ijk}\right)}$	length fraction of line segment
${\nu }_{ijk\ell }=\frac{\#\left(B_{ijk\ell }\cap V_{ijk\ell }\right)}{\#\left(B_{ijk\ell }\right)}$	$0$ or $1$
${\sigma }_i=\frac{\mathrm{𝐴𝑟𝑒𝑎}\left(S_i\cap V_i\right)}{\mathrm{𝐴𝑟𝑒𝑎}\left(B_i\right)}$	area fraction of sphere
${\sigma }_{ij}=\frac{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(S_{ij}\cap V_{ij}\right)}{\mathrm{𝐿𝑒𝑛𝑔𝑡ℎ}\left(S_{ij}\right)}$	length fraction of circle
${\sigma }_{ijk}=\frac{\#\left(S_{ijk}\cap V_{ijk}\right)}{\#\left(S_{ijk}\right)}$	$0$, $\frac{1}{2}$, or $1$
$w_i;{\varphi }_{ij};{\lambda }_{ij}$	real weight; angle between normals; combined length of projection
${\varphi }_{i,jk}={\alpha }_i{\varphi }_{ijk}$	area of spherical quadrangle; area fraction, area of spherical triangle
$X,\bigcup X$	set of balls, space-filling diagram
$𝔽,{𝕊}^2,{𝔹}^3,{ℝ}^3$	surface, unit sphere, unit ball, Euclidean space
$𝐱,𝐭;𝐯,𝐚,𝐦,𝐠$	state, momentum; gradients
$F:{ℝ}^{3n}\to ℝ$	intrinsic volume function
${\kappa }_1\left(x\right),{\kappa }_2\left(x\right);\chi \left[𝔽\right]$	principal curvatures; Euler characteristic
$P,{\alpha }^P;{\xi }_i$	intersection point, parametrizing angle; signed distance
${𝐧}_i,{𝐦}_i,𝐳;{\alpha }_i$	unit normal, midpoint, center; area fraction
$a,b,c;s$	squared cosines; half-perimeter
$S,A,B,C$	area functions
${\mathrm{sgm}}_{i,jk}$	orientation of isosceles triangle
$R_{ijk};r$	radius; squared cosine of half-radius
$u,v,w,t$	generic functions
${\mathrm{𝑔𝑎𝑢𝑠𝑠}}^{\prime }=d^{\prime }+e^{\prime }+f^{\prime }+h^{\prime }$	derivatives
$𝐝,𝐞,𝐟,𝐡;{𝐝}_i,{𝐞}_i,{𝐟}_i,{𝐡}_i$	gradients; as they apply to $x_i$
$p_{ij},q_{ij},s_{ij}$	coefficients
$d_{ij},d_{ijk},e_{ijk},{\overline{e}}_{ijk},f_{ij},h_{i,jk}$	coefficients
$D;{𝐚}_{ijk},{𝐛}_{ijk},{𝐜}_{ijk},{𝐝}_{ijk}$	auxiliary constant; auxiliary vectors
${𝕄}_I,{𝕄}_{\mathrm{II}}$	subsets of ${ℝ}^{3n}$

Equations

$$\displaystyle{\mbox{$\it Gauss$}}({{\mathbb{F}}})$$

$$\mathbf{x})$$

$$\displaystyle{\mbox{$\it area$}}({\bf x})$$

$$\displaystyle{\mbox{$\it mean$}}({\bf x})$$

$$\displaystyle{\mbox{$\it gauss$}}({\bf x})$$

$$\displaystyle{\mbox{$\it Area$}}({\bf n}_{i}{\bf m}_{k}{\bf z}{\bf m}_{j})$$

$${\xi }_i$$

$${\lambda }_{ij}$$

$$\frac{\mathrm{d}{\lambda }_{ij}}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$${\lambda }_{ij}^{\prime }$$

$$\sin \frac{{\varphi }_{ijk}}{2}$$

$$\sin s\sin (s-{\varphi }_{ij})\sin (s-{\varphi }_{jk})\sin (s-{\varphi }_{ki})$$

$$S(a,b,c)$$

$${\varphi }_{i,jk}$$

$${\varphi }_{j,ki}$$

$${\varphi }_{k,ij}$$

$$\displaystyle{\mbox{$\rm sgm$}}_{i,jk}$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}$$

$$\frac{\mathrm{d}t}{\mathrm{d}w}$$

$$\frac{\mathrm{d}w}{\mathrm{d}a}$$

$$\frac{\mathrm{d}S}{\mathrm{d}a}$$

$$\tan R_{ijk}$$

$$\tan R(a,b,c)$$

$$r(a,b,c)$$

$$=\frac{1}{2}+\frac{1}{2}\sqrt{\frac{4abc-(a+b+c-1{)}^2}{(a-1{)}^2+(b-1{)}^2+(c-1{)}^2-(a-b{)}^2-(a-c{)}^2-(b-c{)}^2}}.$$

$$u(a,b,c)$$

$$v(a,b,c)$$

$$\frac{\mathrm{d}r}{\mathrm{d}a}$$

$$\frac{\mathrm{d}{\varphi }_{ij}}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$$\frac{\mathrm{d}A}{\mathrm{d}\Vert x_i-x_j\Vert }$$

$$\frac{\mathrm{d}A}{\mathrm{d}\Vert x_j-x_k\Vert }$$

$$\frac{\mathrm{d}A}{\mathrm{d}\Vert x_k-x_i\Vert }$$

$${\varphi }_{i,jk}^{\prime }$$

$$p_{ij}$$

$$p_{jk}$$

$$p_{ki}$$

$${\varphi }_{j,ki}^{\prime }$$

Full text

††volume: (

1 Introduction

This is a companion paper to [8] and [5] on the derivatives of the weighted volume and of the weighed area of a space-filling diagram, and particularly of [2] on the derivative of the weighted mean curvature, with which this paper shares the notation. Together with these three papers, we complete the explicit description of the derivatives of the weighted intrinsic volumes in ${{\mathbb{R}}}^{3}$. All four derivatives are needed to realize the morphometric approach of Mecke and Roth to the solvation free energy; see [11, 12, 13, 14, 15]. Getting its inspiration from the theory of intrinsic volumes, this approach writes the solvation free energy as a linear combination of the weighted volume, the weighted area, the weighted mean curvature, and the weighted Gaussian curvature. Compare this with Hadwiger’s characterization theorem, which asserts that any measure of convex bodies in ${{\mathbb{R}}}^{3}$ that is invariant under rigid motion, continuous, and additive is a linear combination of the four (unweighted) intrinsic volumes [10].

Among the four intrinsic volumes, the Gaussian curvature is distinguished by the Gauss–Bonnet theorem, which asserts that it is $2\pi$ times the Euler characteristic of the surface [4]. Translated to a set of spheres in motion, this implies that the curvature function is piecewise constant and thus lacks any interesting derivative. The situation changes drastically in the weighted case, in which the contribution of a sphere to the Gaussian curvature is scaled by a real weight that represents a physical property of the corresponding atom, such as its hydrophobicity. The main result of this paper is an explicit formula for the derivative of the weighted Gaussian curvature of a space-filling diagram, and an analysis of the cases in which this derivative either does not exist or is not continuous. In the special case in which all weights are equal, we still get zero derivative almost everywhere and jumps with at places of undefined derivative when the space-filling diagram undergoes topological changes.

Outline. Section 2 reviews the background relevant to this paper. Section 3 discusses the weighted Gaussian curvature of a space-filling diagram. Section 4 derives the constituents of the weighted Gaussian curvature. Section 5 states the derivative in terms of the gradient. Section 6 analyzes the configurations at which the gradient is not continuous. Section 7 concludes this paper.

2 Background

This section reviews the space-filling diagram and its dual alpha shape, both geometric models of molecules.

Space-filling diagrams. Throughout this paper, we let $X$ be a collection of $n$ closed balls in ${{\mathbb{R}}}^{3}$, and to remind us of the connection to chemistry and biology, we call the union of these balls a space-filling diagram. Writing $x_{i}\in{{\mathbb{R}}}^{3}$ for the center and $r_{i}>0$ for the radius, the corresponding ball is $B_{i}=\{a\in{{\mathbb{R}}}^{3}\mid{\|{a}-{x_{i}}\|}\leq r_{i}\}$, and the space-filling diagram is $\bigcup X=\bigcup_{i=0}^{n-1}B_{i}$; see the $2$-dimensional case illustrated in Figure 1, in which $\bigcup X$ is a union of disks. Its boundary consists of (spherical) patches of the form $S_{i}\setminus\bigcup_{j\neq i}{\rm int\,}{B_{j}}$, in which $S_{i}={\rm bd\,}{B_{i}}$. These patches meet along closed (circular) arcs of the form $S_{ij}\setminus\bigcup_{k\neq i,j}{\rm int\,}{B_{k}}$, in which $S_{ij}=S_{i}\cap S_{j}$. Finally, these arcs meet at corners of the form $S_{ijk}\setminus\bigcup_{\ell\neq i,j,k}{\rm int\,}{B_{\ell}}$, in which $S_{ijk}=S_{i}\cap S_{j}\cap S_{k}$. For ease of reference, we collect and extend the introduced notation in Table 1, which is given in the appendix.

We are interested in the volume of the space-filling diagram and the area, mean curvature, and Gaussian curvature of its boundary. Since this boundary is not necessarily smoothly embedded in ${{\mathbb{R}}}^{3}$, we use discrete formulas for the curvatures that agree with the smooth formulas in the limit. For the mean curvature, we have non-zero contributions from the patches and the arcs; see Equation (4), and for the Gaussian curvature, we have non-zero contributions from the patches, the arcs, and the corners; see Equation (5). However, according to the Gauss–Bonnet theorem [4], the Gaussian curvature is invariant under deformations that preserve its topology, and it changes by integer multiples of $2\pi$ whenever the topology is not preserved. The derivative under motion of the balls is therefore not very informative, namely generically zero and occasionally undefined. This is in sharp contrast of the weighted case.

Voronoi decomposition and dual alpha complex. We need a decomposition of the space-filling diagram and of its boundary to define what we mean by its weighted measures. Recall that $B_{i}$ is the closed ball with center $x_{i}$ and radius $r_{i}$. The power distance of a point $a\in{{\mathbb{R}}}^{3}$ from $B_{i}$ is $\pi_{i}(a)={\|{a}-{x_{i}}\|}^{2}-r_{i}^{2}$. Note that $B_{i}=\pi_{i}^{-1}(-\infty,0]$. The Voronoi domain of $B_{i}$ is the set $V_{i}$ of points $a\in{{\mathbb{R}}}^{3}$ for which $\pi_{i}(a)\leq\pi_{j}(a)$ for $0\leq j<n$. It is a possibly unbounded and always closed and convex $3$-dimensional polyhedron; see Figure 1 for the $2$-dimensional case in which the Voronoi domains are convex polygons. Generically, $V_{ij}=V_{i}\cap V_{j}$ is either empty or $2$-dimensional, $V_{ijk}=V_{i}\cap V_{j}\cap V_{k}$ is either empty or $1$-dimensional, and $V_{ijk\ell}=V_{i}\cap V_{j}\cap V_{k}\cap V_{\ell}$ is either empty or a point. The Voronoi domains decompose the space-filling diagram into convex sets of the form $B_{i}\cap V_{i}$. The boundary of such a set consist of a spherical patch (which also belong to the boundary of the space-filling diagram) and a collection of flat patches (along which the sets meet). The spherical patch meets the flat patches along circular arcs on the boundary of $\bigcup X$, and the flat patches meet along straight line segments inside $\bigcup X$.

We use the notation in Table 1 for the fractional measures of the pieces in the decomposition. They are conveniently computed using the alpha complex of $X$, which is the nerve of the collection of $B_{i}\cap V_{i}$ [9]. We recall that the Delaunay mosaic is the nerve of the collection of $V_{i}$, which implies that the alpha complex is a subcomplex of the Delaunay mosaic. Just like the Delaunay mosaic, the alpha complex is geometrically realized by mapping its vertices to the corresponding points in ${{\mathbb{R}}}^{3}$. The alpha shape is the underlying space of the alpha complex, and we stress that this is a set of points while the alpha complex is a set of simplices. The boundary of the alpha shape consists of all points for which no open neighborhood is contained in the set. It is the underlying space of a subcomplex of the alpha complex, and we refer to its members as boundary simplices. The alpha complex contains a simplex for every non-empty common intersection of Voronoi domains, which makes it a convenient book-keeping device for the intended computations. To appreciate this facility, it is useful to familiarize ourselves with the details of the correspondence between the boundaries of the space-filling diagram and the alpha shape. Whenever $S_{i}$ contributes a non-zero number of patches, its center, $x_{i}$, is a boundary vertex of the alpha shape. The number and structure of the patches can be extracted from the structure of simplices in the alpha complex that share $x_{i}$. Whenever $S_{ij}$ contributes a non-zero number of arcs, the edge connecting $x_{i}$ to $x_{j}$ belongs to the boundary of the alpha shape. The triangles and tetrahedra that share the edge are arranged cyclically around the edge, and the number of gaps between contiguous triangles in the cyclic order is equal to the number of contributed arcs. Finally, whenever $S_{ijk}$ contributes a non-zero number of corners—namely one or two—the triangle with vertices $x_{i},x_{j},x_{k}$ belongs to the boundary of the alpha shape. The number of contributed corners is equal to the number of sides of the triangle that face the outside of the alpha shape, which is again one or two. Similar but simpler rules of correspondence between the boundary of the union of disks and the dual alpha shape in two dimensions can be seen in Figure 1.

3 Weighted Gaussian Curvature

According to the Gauss–Bonnet theorem [4], the (integrated) Gaussian curvature of a closed surface, ${{\mathbb{F}}}$, that is smoothly embedded in ${{\mathbb{R}}}^{3}$ is a fixed multiple of the Euler characteristic of the surface:

[TABLE]

in which $\kappa_{1}(x)$ and $\kappa_{2}(x)$ are the two principal curvatures at $x$. The boundary of a space-filling diagram is generally not smoothly embedded, but the extension of the Gaussian curvature used in this paper still satisfies (1). This implies that the (unweighted) Gaussian curvature is constant while the spheres move, until a time at which the surface changes topologically and the Gaussian curvature jumps to a possibly different value. In other words, the derivative of the unweighted Gaussian curvature is generally zero, and undefined whenever the topology changes. The weighted case is radically different because the contributions of the individual spheres to the weighted Gaussian curvature do not cancel each others changes. In contrast to the unweighted case, the computation of the derivative of the weighted Gaussian curvature of a space-filling diagram is a meaningful as well as challenging undertaking.

Weighted intrinsic volume. Since the boundary of a space-filling diagram is only piecewise smooth, we compute its weighted Gaussian curvature separately for the sphere patches, the circular arcs, and the corners. To facilitate a comparison, it is convenient to give such formulas for the weighted versions of all four intrinsic volumes in ${{\mathbb{R}}}^{3}$: the weighted volume, the weighted area, the weighted mean curvature, and the weighted Gaussian curvature. Recall that $\bigcup X$ is the union of the balls $B_{i}$, for $0\leq i<n$. Its state, ${\bf x}\in{{\mathbb{R}}}^{3n}$, is the concatenation of the center vectors. In other words, for $1\leq\ell\leq 3$, the $(3i+\ell)$-th coordinate of ${\bf x}$ is the $\ell$-th coordinate of $x_{i}$. Decomposing the space-filling diagram with the Voronoi tessellation, we write ${{w_{i}}}\in{{\mathbb{R}}}$ for the weight of $B_{i}$ and use the notation introduced in Table 1 to give formulas for the weighted intrinsic volumes.

Proposition 3.1** (Weighted Intrinsic Volumes).**

The weighted volume, surface area, mean curvature, and Gaussian curvature of the space-filling diagram, $\bigcup X$, are

[TABLE]

The weighted Gaussian curvature function decomposes into three sums, which account for the sphere patches, the circular arcs, and the corners, respectively. The sums are over ordered index sequences. This implies that the same three indices appear six times, which explains the coefficient of $\tfrac{1}{3}$ in front of this sum. We will see that $\phi_{i,jk}=\alpha_{i}\phi_{ijk}$ is the area of a spherical quadrangle, so we can simplify the last sum. Formulas for the derivatives of $\sigma_{i}$ and of $\sigma_{ij}$ can be found in [2]. We also need derivatives of $\lambda_{ij}$ (the combined length of the projection of two normals; see Figure 3), and of $\phi_{i,jk}$ (the area of a spherical quadrangle), which will be computed in Sections 4.1 and 4.2. Similar to the state, ${\bf x}$, we write ${\bf t}\in{{\mathbb{R}}}^{3n}$ for the momentum, which is obtained by concatenating the velocity vectors of the $n$ balls in the order of their indices. We write $\lambda_{ij}(\tau)$ for $\lambda_{ij}$ at state ${\bf x}+\tau{\bf t}$ and $\lambda_{ij}^{\prime}$ for the derivative of $\lambda_{ij}$ at $\tau=0$, and similarly for $\sigma_{ijk}$ and $\phi_{i,jk}$. Derivatives with respect to parameters other than $\tau$ are explicitly stated as such.

How to split a corner. Every corner of the space-filling diagram has a non-zero contribution to the Gaussian curvature, and since the spheres have weights, we need to decide how to distribute this contribution among the generically three spheres that meet at the corner. We do this in a way that is consistent with splitting the curvature concentrated along an arc in equal halves, which is motivated by the Appolonius diagram of the spheres; see [2].

Consider spheres $S_{i},S_{j},S_{k}$, and suppose that $P$ is one of the two points in which the spheres generically intersect. Write ${\bf n}_{i},{\bf n}_{j},{\bf n}_{k}$ for the unit outward normals of the spheres at $P$, and recall that $\phi_{ijk}$ is the area of the spherical triangle with vertices ${\bf n}_{i},{\bf n}_{j},{\bf n}_{k}$ on the unit sphere; see Figure 2.

We split the contribution of the corner using $\alpha_{i}+\alpha_{j}+\alpha_{k}=1$, in which the three coefficients are the fractions of the areas of the spherical quadrangles obtained by decomposing the spherical triangle with the great-circle arcs that connect the center of the circumcircle, denoted ${\bf z}$, with the midpoints of the sides, denoted ${\bf m}_{i},{\bf m}_{j},{\bf m}_{k}$. By connecting ${\bf z}$ with the three vertices, we further decompose each quadrangle into two spherical triangles, which we refer to as cones of ${\bf z}$ over the half-sides; see again Figure 2. We remark that ${\bf z}$ can lie outside the spherical triangle, in which case the areas of two of the cones are considered to be negative. The cones come in pairs of equal area, and each pair glues together to form an isosceles spherical triangle. As remarked, exactly one of the three isosceles triangles has negative area if $z$ is not contained in the original spherical triangle. Since each of the three quadrangles is the union of two cones, we can write its area as half the sum of the area of two isosceles triangles:

[TABLE]

and similarly for ${\it Area}({\bf n}_{j}{\bf m}_{i}{\bf z}{\bf m}_{k})$ and ${\it Area}({\bf n}_{k}{\bf m}_{j}{\bf z}{\bf m}_{i})$. Equivalently, the area of the quadrangle is half of ${\it Area}({\bf n}_{i}{\bf n}_{j}{\bf n}_{k})-{\it Area}({\bf n}_{j}{\bf n}_{k}{\bf z})$, but we prefer (6) for our computations.

4 Derivatives

We reuse the expressions for $\sigma_{i}^{\prime}$ and $\sigma_{ij}^{\prime}$ in [2], and we develop formulas for $\lambda_{ij}^{\prime}$, $\sigma_{ijk}^{\prime}$, $\phi_{i,jk}^{\prime}$. The middle of the three, $\sigma_{ijk}^{\prime}$, is piecewise zero because $\sigma_{ijk}$ is piecewise constant. As will be discussed in Section 6, $\sigma_{ijk}$ changes at non-generic states, where its derivative is not defined. We first consider the derivative of $\lambda_{ij}$ and second that of $\phi_{i,jk}$.

4.1 Derivative of $\lambda_{ij}$

Given spheres $S_{i}$ and $S_{j}$, we define $\lambda_{ij}$ as the combined length of the unit normal vectors at a point $P\in S_{i}\cap S_{j}$ after projection to the line that passes through the centers, $x_{i}$ and $x_{j}$; see Figure 3. If the two vectors point in the same direction—after projection—then we take one length negative, so more precisely, $\lambda_{ij}$ is the distance between the projections of the endpoints of the unit vectors.

To compute the derivative of $\lambda_{ij}$, it will be useful to recall that the signed distance of $x_{i}$ from the Voronoi plane of the two spheres is

[TABLE]

see [2]. The length of the unit normal vector of $S_{i}$ at $P$, after projection, is $\xi_{i}/r_{i}$, and since $\xi_{i}$ is signed so is this length, as required. Adding the two signed lengths, we get

[TABLE]

The derivative with respect to the distance between the two centers is therefore

[TABLE]

To translate this result to the derivative with respect to the motion vectors ${\bf t}_{i}$ of $x_{i}$ and ${\bf t}_{j}$ of $x_{j}$, we note that ${\bf u}_{ij}=(x_{i}-x_{j})/{\|{x_{i}}-{x_{j}}\|}$ is the unit vector that points from $x_{j}$ to $x_{i}$.

Lemma 4.1** (Derivative of $\lambda_{ij}$).**

The derivative of the combined length of the unit normal vectors at an intersection point of spheres $S_{i}$ and $S_{j}$—after projection to the line passing through the centers—at state ${\bf x}$ with momentum ${\bf t}$ is

[TABLE]

in which the coefficient is given in (9).

4.2 Derivative of $\phi_{i,jk}$

Instead of deriving $\alpha_{i}$ and $\phi_{ijk}$ separately, we recall that $\phi_{i,jk}=\alpha_{i}\phi_{ijk}$ is the area of a spherical quadrangle, which we can derive directly. We begin with intrinsic formulas for the area and the radius of a spherical triangle.

Area of spherical triangle. Recall that $\phi_{ijk}$ is the area and $\phi_{ij},\phi_{jk},\phi_{ki}$ are the lengths of the sides of the spherical triangle with vertices ${\bf n}_{i},{\bf n}_{j},{\bf n}_{k}$. Writing $s=\tfrac{1}{2}(\phi_{ij}+\phi_{jk}+\phi_{ki})$ for half the perimeter, we have the following formula for the sine of half the area:

[TABLE]

see [6, page 230, equation (225)]. To rewrite this formula, we introduce the squared cosines of the half-lengths as parameters. Geometrically, $\cos^{2}\tfrac{x}{2}$ is the distance from the center to the midpoint of the straight edge that connects the endpoints of an arc of length $x$ on the unit sphere. The square of this distance conveniently relates to Pythagoras’ theorem for right-angled triangles. In a first step, we rewrite the product of the four sines. {formula}[Product of Sines]

Writing $a=\cos^{2}\frac{\phi_{ij}}{2}$, $b=\cos^{2}\frac{\phi_{jk}}{2}$, $c=\cos^{2}\frac{\phi_{ki}}{2}$, the expression under the square root in (11) satisfies

[TABLE]

We refer to Appendix A for a proof. We can thus rewrite (11) to give the area of any spherical triangle given the squared cosines of its half-lengths:

[TABLE]

By (11) and (12), we have $\phi_{ijk}=S(a,b,c)$ with parameters as given in Formula 4.2. We introduce notation for the areas of the isosceles triangles obtained by gluing equal-area cones: $A(a,r)=S(a,r,r)$, $B(b,r)=S(r,b,r)$, $C(c,r)=S(r,r,c)$, with $r=\cos^{2}\frac{R_{ijk}}{2}$, in which $R_{ijk}$ is the radius of the cap bounded  by the circle that passes through ${\bf n}_{i},{\bf n}_{j},{\bf n}_{k}$. By symmetry between the arguments of $S$, all three functions are the same, but it is convenient to distinguish them so we do not lose the connection to the sides of the spherical triangle. Using (6), we get the areas of the spherical quadrangles:

[TABLE]

in which the signs depend on whether ${\bf z}$ lies inside or outside the triangle spanned by ${\bf n}_{i},{\bf n}_{j},{\bf n}_{k}$. Specifically, ${\rm sgm}_{i,jk}=+1$ if ${\bf z}$ and ${\bf n}_{i}$ lie on the same side of the great-circle that passes through ${\bf n}_{j}$ and ${\bf n}_{k}$, and ${\rm sgm}_{i,jk}=-1$ if the two points lie on opposite sides of this great-circle. In the boundary case, when ${\bf z}$ lies on the great-circle, ${\bf n}_{i}$, ${\bf n}_{j}$, ${\bf n}_{k}$ form a Euclidean right-angled triangle. Pythagoras’ theorem for its edges is equivalent to $\sin^{2}\frac{\phi_{ij}}{2}+\sin^{2}\frac{\phi_{ki}}{2}=\sin^{2}\frac{\phi_{jk}}{2}$. If ${\bf z}$ and ${\bf n}_{i}$ lie on the same side of the great-circle, then the left-hand side is strictly larger than the right-hand side. Substituting $\sin^{2}x=1-\cos^{2}x$, we get

[TABLE]

Since $S(a,b,c)$ is symmetric in its three arguments, it suffices to compute the derivative with respect to one of them. To this end, we recall a few basic rules of differentiation, namely ${\rm d}\arcsin x/{\rm d}x=1/\sqrt{1-x^{2}}$, ${\rm d}\sqrt{x}/{\rm d}x=1/\sqrt{4x}$, and ${\rm d}(\tfrac{u}{v})/{\rm d}x=(v{\rm d}u/{\rm d}x-u{\rm d}v/{\rm d}x)/v^{2}$ for the ratio of two functions in $x$. Setting $u(a,b,c)=4abc-(a+b+c-1)^{2}$ and $v(a,b,c)=4abc$, we get ${\rm d}u/{\rm d}a=4bc-2(a+b+c-1)$ and ${\rm d}v/{\rm d}a=4bc$. Furthermore setting $w=\tfrac{u}{v}$ and $t=\sqrt{w}$, we have $S(a,b,c)=2\arcsin t$ and therefore

[TABLE]

Multiplying the three right-hand sides and simplifying the result, we get the derivative of the area with respect to the squared cosine of one half-length:

[TABLE]

Permuting the three arguments, we get symmetric expressions for ${\rm d}S/{\rm d}b$ and for ${\rm d}S/{\rm d}c$.

Radius of spherical triangle. We need a formula for the radius of the cap, $R_{ijk}$, in the same parametrization:

[TABLE]

see [6, page 246, equation (307)]. Writing this in terms of the squared cosines, we get

[TABLE]

To turn this into a formula for the squared cosine of the half-radius, we recall two basic trigonometric identities: $\cos^{2}\tfrac{x}{2}=\tfrac{1}{2}[1+\cos x]=\tfrac{1}{2}[1+1/\sqrt{\tan^{2}x+1}]$. Writing $r(a,b,c)=\cos^{2}\frac{R(a,b,c)}{2}$ and using (25), this gives

[TABLE]

Since $r(a,b,c)$ is symmetric in the three arguments, it suffices to compute the derivative with respect to one of them. Referring to (27), we set

[TABLE]

which give ${{\rm d}u}/{{\rm d}a}=4bc-2(a+b+c-1)$ and ${{\rm d}v}/{{\rm d}a}=2(-a+b+c-1)$. Observing that $r(a,b,c)=\tfrac{1}{2}+\tfrac{1}{2}\sqrt{u(a,b,c)/v(a,b,c)}$, we see that its derivative with respect to the first argument is $(v{\rm d}u/{\rm d}a-u{\rm d}v/{\rm d}a)/\left(4v\sqrt{v}\sqrt{u}\right)$, which implies

[TABLE]

Symmetric expressions hold for ${\rm d}r/{\rm d}b$ and for ${\rm d}r/{\rm d}c$.

Using the chain rule. We get the derivatives of the areas of the spherical quadrangles by combining the results with the chain rule. Beyond the relations above, we need the derivative of the length of a side, which we find in [2, equation (30)]:

[TABLE]

We recall that ${\rm d}\cos^{2}\frac{x}{2}/{\rm d}x=-\sin\frac{x}{2}\cos\frac{x}{2}$, and we use (14) to express the area of the spherical quadrangle: $\phi_{i,jk}=\tfrac{1}{2}[{\rm sgm}_{k,ij}A(a,r)+{\rm sgm}_{j,ki}C(c,r)]$, in which $a=\cos^{2}\frac{\phi_{ij}}{2}$, $b=\cos^{2}\frac{\phi_{jk}}{2}$, $c=\cos^{2}\frac{\phi_{ki}}{2}$, and $r=\cos^{2}\frac{R_{ijk}}{2}$ as given in (27). The derivatives of $A$ with respect to its arguments are ${\rm d}A/{\rm d}a={\rm d}S/{\rm d}a$ and ${\rm d}A/{\rm d}r={\rm d}S/{\rm d}b+{\rm d}S/{\rm d}c$, both evaluated at $b=c=r$. The derivatives with respect to the distances between the centers of the spheres are therefore

[TABLE]

(33) and (34) are symmetric because the same two sides of the isosceles triangle depend on ${\|{x_{j}}-{x_{k}}\|}$ and on ${\|{x_{k}}-{x_{i}}\|}$, while (32) is different because all three sides depend on ${\|{x_{i}}-{x_{j}}\|}$. We can now summarize and state the derivative of $\phi_{i,jk}$ with respect to the motion vector.

Lemma 4.2** (Derivative of $\phi_{i,jk}$).**

The derivative of the area of the spherical quadrangle with vertices ${\bf n}_{i},{\bf m}_{k},{\bf z},{\bf m}_{j}$ on the unit sphere is

[TABLE]

in which the derivatives on the right-hand side are given in or symmetric to (32), (33), (34).

Starting with (15), (16) instead of (14), we get

[TABLE]

with symmetrically defined coefficients.

5 Gradients

For computational purposes, it is convenient to write the derivative of the weighted Gaussian curvature function, ${\it gauss}\colon{{\mathbb{R}}}^{3n}\to{{\mathbb{R}}}$, in terms of the gradient of ${\it gauss}$ at ${\bf x}\in{{\mathbb{R}}}^{3n}$, denoted ${\bf g}=\nabla_{\!{\bf x}}{\it gauss}$. Recalling (5), this derivative is ${\it gauss}^{\prime}=d^{\prime}+e^{\prime}+f^{\prime}+h^{\prime}$, with

[TABLE]

Writing ${\bf g}=[g_{1},g_{2},\ldots,g_{3n}]^{T}$, we recall that ${\bf g}_{i}=[g_{3i+1},g_{3i+2},g_{3i+3}]^{T}$ is the $3$-dimensional gradient that applies to $x_{i}$. Using boldface letters for the gradients of $d$, $e$, $f$, $h$, and similar conventions for the $3$-dimensional sub-vectors, we have ${\bf g}={\bf d}+{\bf e}+{\bf f}+{\bf h}$ and ${\bf g}_{i}={\bf d}_{i}+{\bf e}_{i}+{\bf f}_{i}+{\bf h}_{i}$ for $0\leq i<n$. We get the gradients by redistributing the relevant derivatives given in [2] and stated in Lemmas 4.1 and 4.2.

First term. The derivative of $\sigma_{i}$ is given in [2, equations (19) and (21)] as

[TABLE]

in which the first sum is over the boundary edges of the alpha complex incident to $x_{i}$, and the second sum is over the triangles incident to these edges. We can therefore rewrite (41) as

[TABLE]

in which the first sum ranges over all (directed) boundary edges of the alpha shape, and second sum ranges over all triangles incident to these edges. We get ${\bf d}_{i}$ by accumulating the relevant contributions. Using ${\bf u}_{ji}=-{\bf u}_{ij}$ and ${\bf u}_{jik}={\bf u}_{ijk}$, we get

[TABLE]

in which the first sum is over all boundary edges incident to $x_{i}$, and the second sum is over all triangles incident to these edges. Compare (49) with the weighted area gradient given in [5] and with the first term of the weighted mean curvature gradient given in [2].

Second term. This case is complicated by the change of motion needed to compute the derivative of $\sigma_{ij}$. Following the exposition in [2, Section 4], we introduce some notation before we can rewrite (42). Setting $D=\tfrac{1}{2}[1-(r_{i}^{2}-r_{j}^{2})/{\|{x_{i}}-{x_{j}}\|}^{2}]$ and writing ${\bf d}_{ijk}=(x_{k}-Dx_{j}+(D-1)x_{i})/{\|{x_{i}}-{x_{j}}\|}$, we introduce three vectors so that the scalar product of the new motion vector with $x_{k}-P$ can be written as a sum of three scalar products, each involving only one of the original motion vectors:

[TABLE]

compare with [2, equations (54) to (59)]. In addition, we recall the derivatives of $\sigma_{ij}$, $r_{ij}$, and $\alpha^{P}$ from [2]:

[TABLE]

in which ${\alpha}^{P}$ is the angle parametrizing the motion of the corner, $P$, along the circle, $S_{ij}$. With this, we rewrite (42) as

[TABLE]

in which the sum is over the three ordered versions of all oriented boundary triangles of the alpha shape, and $P$ is the corresponding corner of the space-filling diagram. Redistributing the terms, we get

[TABLE]

in which the sum is again over the three ordered versions of all oriented boundary triangles of the alpha shape that share $x_{i}$.

Third term. Using Lemma 4.1, we rewrite (43) as

[TABLE]

in which the sum ranges over all (directed) boundary edges of the alpha shape. We get ${\bf f}_{i}$ by accumulating the relevant contributions, as before:

[TABLE]

in which the sum is over boundary edges incident to $x_{i}$.

Fourth term. Using Lemma 4.2, we rewrite (44) as

[TABLE]

in which the coefficients are defined in and after Lemma 4.2, and the sum is over all (ordered) triangles in the boundary of the alpha shape. We get ${\bf h}_{i}$ by accumulating the relevant contributions:

[TABLE]

in which the sum is over all ordered boundary triangles of the alpha shape that are incident to $x_{i}$. To see this, we observe that the three ordered versions of the same orientation, $ijk,jki,kij$, all contribute the same term. By symmetry of the spherical quadrangles defined for the opposite orientation, the corresponding three ordered versions, $ikj,kji,jik$, also contribute the same term.

Summary. We finally get the gradient of the weighted Gaussian curvature function by adding the gradients of the four component functions:

Theorem 5.1** (Gradient of Weighted Gaussian Curvature).**

The derivative of the weighted Gaussian curvature of the space-filling diagram at state ${\bf x}$ with momentum ${\bf t}$ is ${\rm D}{\it gauss}_{{\bf x}}({\bf t})={\langle{\bf g},{\bf t}\rangle}$, in which ${\bf g}_{i}={\bf d}_{i}+{\bf e}_{i}+{\bf f}_{i}+{\bf h}_{i}$ as given in (49), (60), (63), (68), for $0\leq i<n$.

6 Continuity

Like the derivatives of the other intrinsic volumes of a space-filling diagram, the derivative of the weighted Gaussian curvature is almost everywhere but not everywhere continuous. Indeed, the situation is sufficiently similar to the mean curvature case that we can follow the analysis in [2] and be brief.

General position. The derivative of ${\it gauss}\colon{{\mathbb{R}}}^{3n}\to{{\mathbb{R}}}$ is continuous everywhere except at non-generic states. To define what this means, we recall that ${\bf x}\in{{\mathbb{R}}}^{3n}$ encodes a collection, $X$, of $n$ closed balls in ${{\mathbb{R}}}^{3}$. We say $X$ is in general position and, equivalently, ${\bf x}$ is generic, if the following two conditions are satisfied:

1. I.

the common intersection of $p+1$ Voronoi domains is either empty or a convex polyhedron of dimension $3-p$; 2. II.

the common intersection of $p+1$ spheres bounding balls in $X$ is either empty or a sphere of dimension $2-p$.

Condition I implies that the Delaunay mosaic of $X$ is a simplicial complex in ${{\mathbb{R}}}^{3}$. When the balls move, the mosaic undergoes combinatorial changes, and these happen at states that violate Condition I. In addition, Condition II captures states at which simplices in the mosaic are added to or removed from the alpha complex.

Each violation of the two conditions corresponds to a submanifold of dimension at most $3n-1$ in ${{\mathbb{R}}}^{3n}$. Write ${{\mathbb{M}}_{\rm{I}}}$ and ${{\mathbb{M}}_{\rm{II}}}$ for the unions of submanifolds that correspond to violations of Condition I and of Condition II, respectively. Since there are only finitely many submanifolds, and at least some of each type have dimension $3n-1$, ${{\mathbb{M}}_{\rm{I}}}$ and ${{\mathbb{M}}_{\rm{II}}}$ have dimension $3n-1$ each.

Combinatorial changes. The submanifolds decompose ${{\mathbb{R}}}^{3n}$ into cells, and we refer to the connected components of ${{\mathbb{R}}}^{3n}\setminus({{\mathbb{M}}_{\rm{I}}}\cup{{\mathbb{M}}_{\rm{II}}})$ as chambers. As long as the trajectory of states stays inside a chamber, the Delaunay mosaic and the alpha complex are invariant, combinatorially, and so are the formulas for the weighted Gaussian curvature and its derivative. The functions defined by the formulas are continuous, which implies that ${\it gauss}$ and $\nabla{\it gauss}$ restricted to a single chamber are continuous.

When the trajectory passes from one chamber to another, it necessarily crosses a submanifold, and we let ${\bf x}$ be the corresponding non-generic state. If ${\bf x}\in{{\mathbb{M}}_{\rm{I}}}$, then the associated combinatorial change is a flip in the Delaunay mosaic; see [7] for details. As argued in [2], ${\it mean}$ and $\nabla{\it mean}$ are both continuous at ${\bf x}$, and the same argument implies that also ${\it gauss}$ and $\nabla{\it gauss}$ are continuous at ${\bf x}$. If ${\bf x}\in{{\mathbb{M}}_{\rm{II}}}$, then the associated combinatorial change is the addition of an interval of simplices in the Delaunay mosaic to the alpha complex or, symmetrically, the removal of such an interval from the alpha complex. Following [2], we distinguish between a singular interval, which contains only one simplex, and a non-singular interval, which contains two or more simplices; see [3] for background on the discrete Morse theory of Delaunay mosaics. A singular interval corresponds to two, three, or four spheres whose common intersection is a single point. If this point happens to belong to the boundary of the space-filling diagram, then the operation changes the topology type of the surface. By the Gauss–Bonnet theorem, this implies that in the unweighted case, ${\it gauss}$ changes by a non-zero integer multiple of $2\pi$, and $\nabla{\it gauss}$ is undefined. Similarly, each non-critical interval corresponds to two, three, or four spheres whose common intersection is a single point, with the difference that now one of the spheres is contained in the union of the balls bounded by the other one, two, or three spheres. There are six basic cases, each with its own implications on the continuity of $\nabla{\it gauss}$. The details are not important, so we refer to the analogous analysis in [2].

Theorem 6.1** (Continuity of Gradient).**

The gradient of the weighted Gaussian curvature of a space-filling diagram of $n$ closed balls in ${{\mathbb{R}}}^{3}$ is continuous provided the state ${\bf x}\in{{\mathbb{R}}}^{3n}$ of the diagram does not belong to ${{\mathbb{M}}_{\rm{II}}}$, which is a $(3n-1)$-dimensional subset of ${{\mathbb{R}}}^{3n}$.

Not all states in ${{\mathbb{M}}_{\rm{II}}}$ are problematic during molecular dynamics computations. Some are excluded by the forces, such as two identical balls, and others do not cause discontinuities in the derivative, such as non-generic configurations in the interior of the space-filling diagram. The relevant non-generic configurations are when two, three, or four balls touch in a single point at the surface of the space-filling diagram. For two balls, this event either merges/splits components, or it closes/breaks a loop. For three balls, the event either fills/opens a tunnel, or it completes/punctures a shell surrounding a void. For four balls, the event starts/drowns a void.

7 Discussion

The main contribution of this paper is the analysis of the derivative of the weighted Gaussian curvature of a space-filling diagram. Specifically, we give an explicit description of the gradient of the weighted Gaussian curvature function, which for $n$ spheres is a map from ${{\mathbb{R}}}^{3n}$ to ${{\mathbb{R}}}$. In addition, we study the subset of ${{\mathbb{R}}}^{3n}$ at which the derivative violates continuity. In summary, this is sufficient information for an efficient and robust implementation of the weighted Gaussian curvature derivative, one that can be added to the inner loop of a molecular dynamics simulation of a physical system.

Acknowledgment

The authors of this paper thank Roland Roth for suggesting the analysis of the weighted curvature derivatives for the purpose of improving molecular dynamics simulations. They also thank Patrice Koehl for the implementation of the formulas and for his encouragement and advise along the road. Finally, they thank two anonymous reviewers for their constructive criticism.

Appendix A Proof of Formula 4.2

We present a proof of the formula from Section 4 after restating it for convenience.

Formula 3 (Product of Sines). Writing $a=\cos^{2}\frac{\phi_{ij}}{2}$, $b=\cos^{2}\frac{\phi_{jk}}{2}$, and $c=\cos^{2}\frac{\phi_{ki}}{2}$, the expression under the square root in (11) satisfies

[TABLE]

Proof A.1**.**

Recall the trigonometric identity $\sin(\alpha+\beta+\gamma)=\cos\alpha\cos\beta\sin\gamma+\cos\alpha\sin\beta\cos\gamma+\sin\alpha\cos\beta\cos\gamma-\sin\alpha\sin\beta\sin\gamma$. Setting

[TABLE]

and observing that $\cos(-\alpha)=\cos\alpha$ as well as $\sin(-\alpha)=-\sin\alpha$, we get

[TABLE]

Multiplying the four right-hand sides, we get $-(A^{4}+B^{4}+C^{4}+D^{4})+8ABCD+2(A^{2}B^{2}+A^{2}C^{2}+A^{2}D^{2}+B^{2}C^{2}+B^{2}D^{2}+C^{2}D^{2})$. To rewrite this expression, we recall that

[TABLE]

Developing all the right-hand sides and canceling the terms of degree 4, 5, 6 and more, we get the claimed relation.

Appendix B Notation