Surface charge method for molecular surfaces with curved areal elements I. Spherical triangles

Parametrizing a curved surface with flat triangles in electrostatics problems creates a diverging electric field. One way to avoid this is to have curved areal elements. However, charge density integration over curved patches appears difficult. This paper, dealing with spherical triangles, is the first in a series aiming to solve this problem. Here, we lay the ground work for employing curved patches for applying the surface charge method to electrostatics. We show analytically how one may control the accuracy by expanding in powers of the the arc length (multiplied by the curvature). To accommodate not extremely small curved areal elements, we have provided enough details to include higher order corrections that are needed for better accuracy when slightly larger surface elements are used.


Introduction
Not only of general importance at the molecular level, electrostatic interactions are crucial for biomolecular systems which consist of many charged molecules embedded in the polar solvent water. It is, however, impractical to apply precise quantum methods to even a moderately large biomolecular system. Therefore, numerous efforts were invested in classical approaches [1], which are often categorized according to how the solvent water is treated. Explicit solvent methods treat each water molecule at the level of atomic detail, while implicit solvent methods replace the individual water molecules with some type of smoothed out version.
Explicit solvent methods, such as TIPnP [2], allow description of biomolecular systems at finer detail, but their lack of incorporating mutual polarization of molecules can be a critical deficiency sometimes. The implicit solvent methods [3,4] are in principle less computationally intensive when larger systems are considered; however, their application is limited to systems where fine details of solute-solvent interactions do not play a major role.
We have earlier suggested a classical model of dielectric spheres [5] and used it to approximate atoms and molecules by dielectric polarizable bodies [6]. If the bodies are embedded in a dielectric medium, this model becomes an implicit solvent model and the bodies represent the biomolecules while the solvent is modeled by the dielectric continuum. If, however, the bodies are considered to be in the vacuum, the model becomes an explicit solvent model in which biomolecules and waters are modeled explicitly, with their own parameters. This model readily includes polarization effect and was previously used to investigate the accuracy of our classical formalism in describing interactions between atomic-sized objects [6]. We find that the dielectric spheres model is surprisingly accurate down to distances where chemical bonds start to form. The formalism introduced in [5] allows one to find a solution with arbitrary accuracy for an arbitrary number of interacting dielectric spheres with point charges at their centers.
To extend this formalism to model molecules, one need to account for the non-spherical nature of chemical bonds by either introducing permanent and inducible higher multipoles [7] or expanding the scope of dielectric bodies to include non-spherical shapes. The former extension [7] is appropriate when the separation distances among the molecules are large and details of geometry play a lesser role; the latter extension may better represent regions of significant electronic density and thus might be suitable for shorter separations. For the latter, it is necessary to define the surface of the non-spherical molecule. A widely employed strategy for doing this is to create a grid of points at the outer atoms of the molecule. These points are then used to generate a set of flat triangular patches. The assumption is that one need only use a sufficient number of patches and then an accurate calculation will emerge. Indeed, the approximate surface certainly gets closer and closer to a smooth surface corresponding to the shape of the molecule. However, it does so at the expense of introducing more and more vertices, where the approximate surface has no well-defined normal direction.
For the purpose of finding the solvent accessible surface area or similar quantities, these points of undefined surface normal pose no problem (even though using curved surface element were shown to improve the results [8]). However, it is well known that at such points the electric field diverges. Therefore, when calculating the electrostatic interaction, the vertices of the flat patches create inherently uncontrollable sources of inaccuracy. Developing a better strategy for representing surfaces that are not completely spherical is thus a necessary task.
Two primary difficulties exist. First, with a constraint that the patch must include three specified spatial points, there exists no unique solution when curved surfaces are allowed. The parametrization of the patch becomes a problem. Second, even if one can parametrize a surface patch, one still needs a small expansion parameter to control the accuracy of the results obtained. In the dielectric sphere case, the small expansion parameter is naturally the radius of the sphere divided by the separation distances among the spheres. To overcome the first problem, we need a crisp rule for defining the molecular surface. For this, we use the molecular surface definition of Richards [9] and Connolly [10]. With this definition, only spherical (whole or partial) and partial toroidal surfaces come into play. When limited to those surfaces, however, a small intra-molecular expansion parameter emerges. It is the arc length on the surface multiplied by the local curvature. This paper, focusing on spherical triangles, is the first in a series to investigate the use of curved surface elements for applying the surface charge method (SCM) [5,11] in electrostatics. When a set of points on a sphere is chosen, the set of corresponding charge density values, one for each point/vertex, becomes the set of variables to be solved for. The basic idea of SCM is to write the electric potential and field near the dielectric interfaces in terms of the surface charge density values and then determine those values by applying the appropriate boundary condition.
In the sections below, we begin by describing the surface types that we shall encounter, followed by mentioning the fundamental equations needed for specifying a spherical triangle. We then describe, given the (not yet known) density values at the selected vertices, how may one smooth out the charge distribution completely independent of the coordinate choices. After illustrating explicitly the electric field discontinuity due to nonzero surface charge density, we proceed to calculate the electric field contribution right above/below a vertex of a spherical triangle from the surface charge density of that spherical triangle. We show there the emergence of a small expansion parameter. We then calculate the electric field produced by the charge density of a spherical triangle at a vertex of a non-neighboring spherical triangle. Here we also mention cases when numeric integration is the only way to proceed. We then consider, at a location outside a sphere i, the potential and the field contribution due to the charge density of a triangle on the sphere i. Some relevant mathematical formulas/derivations are provided in the appendices.

Surface types
In the context of implicit solvent, one may use a probe sphere to represent the solvent molecule and represent each atom by a sphere of finite radius. This leads to two surfaces: the solvent accessible surface and the molecular surface. The former is defined by the surface traced by the center of the probe sphere while rolling around the atoms of a molecule; the latter, the molecular surface [9,10,12], is defined by tracing the inwardfacing surface of the probe. In the context of explicit solvent model, molecules of interest and solvent molecules are modeled separately. Here a probe sphere is introduced simply to smooth out the surfaces of dielectric objects so that the surface normal is well defined everywhere in order to properly apply the boundary condition; the size of the probe has no direct physical interpretation.
When using the molecular surface definition [9,10,12] above, there will be three types of surface regions for each generic molecular surface. The type I surface are inner toroidal surfaces, made of the inward-facing surface of the probe, that occur when rolling the probe sphere between two atoms. The type II are the partial inward-facing surfaces of the probe sphere right at the stopping point while its rolling between two atoms is terminated by other atoms in the way. The type III are partial outward-facing spherical surfaces associated with different atoms. It is important to note that the adjoining partial surfaces can only be either the jointing of type I and type II or the jointing of type I and type III. Type II and type III surfaces do not join together. Furthermore, the jointing line between different surface types, when viewed from the perspective of the associated spheres, can be parameterized in such a way that it is of constant polar angle (i.e. only the azimuthal angle varies).
When treating a molecule as an dielectric object with a prescribed dielectric constant and surface delineated by the method of Richards [9] and Connolly [10,12], we know from basic electrostatics that the induced charge only exists on the dielectric boundaries. Consequently, in terms of numerically solving the associated electrostatic problem, we need to assign a charge density value for each surface elements on the dielectric boundary. As we mentioned in the Introduction, using flat triangulation will lead to the undesirable consequence of a diverging electric field. Our aim here is to replace the flat triangles with other curved surface elements to eliminate the unphysical diverging electric field. For surface elements, we choose to use spherical triangles for type II and type III partial surfaces and curved rectangles for type I partial surfaces. This last point will be elaborated later in a subsequent manuscript.
One important point to note is that the curves connecting surface type I and surface type III (or II) can be viewed as a curve of fixed latitude on surface type III (or II). This introduces a special need for treating the surface element as we will detail in the subsequent manuscript. Our main objective is to calculate the electric field produced by the curved surface elements to minimize the amount of brute force numerical integration. In order to have a smooth surface charge density profile, it is also important to develop a method to interpolate the charge density values from the charge densities specified at discrete points. We shall show how we may achieve this goal in a coordinate free fashion.

Spherical triangles
When we partition (a part of) a spherical surface, we may first choose individual points, and then connect between points with the big arc. This uniquely specifies a spherical triangle when given its three vertices. For two points along the same longitude, the great arc is simply the segment of the longitude line connecting these two points; the great arc connecting two points of different longitude is given by the geodesic equation: cot θ = a cos(ϕ − ϕ 0 ) . (1) (For completeness, a derivation is provided in appendix A.)

Parameters for geodesic equation
Without losing generality, let us denote the first of the three vertices of the spherical triangle as the north pole of the sphere. Let the second vertex have spherical coordinates (θ = β 3 , ϕ = 0) and the third vertex have coordinates (θ = β 2 , ϕ = α). Obviously, the great arcs connecting the north pole to both (β 3 , 0) and (β 2 , α) are longitudinal lines. The reason we let the second vertex have spherical coordinates (β 3 , 0) is because β 3 is the length of the arc opposite to cot β 2 = a cos(α − ϕ 0 ) = a(cos ϕ 0 cos α + sin ϕ 0 sin α) (3) we may learn the parameters of the geodesic. Note that we must have 0 < α < π. We then express ϕ 0 and a in terms of the angle α. Multiplying equation (2) by cos α and subtracting it from both sides of equation (3) we obtain cot β 2 − cos α cot β 3 = a sin α ϕ 0 .
We thus have (using the ratio of the two equations above) With similar algebraic operations, we also have From equation (4), we also have ϕ 0 = tan −1 tan β 3 − cos α tan β 2 sin α tan β 2 .
One interesting exercise is to find the length of the great arc connecting point 2 and 3. Given equation (1), we know And the great arc angle is thus given by Using equation (7) Hence 1 + a 2 = sin β 1 ∕ sin α sin β 2 sin β 3 .

Smooth charge density parametrization
On the three vertices of a spherical triangle, one may assign three charge density values. Physically speaking, the charge distribution should be smooth and thus we seek to interpolate the charge density values from those specified on the three vertices. Along a great arc side connecting two vertices of the spherical triangle, we expect the charge density to vary smoothly along the side from one vertex to the other. The charge density profile also should be smooth within each spherical triangle. To achieve this goal, we adapt the following interpolation scheme starting with a smooth charge distribution on a flat triangle and project it to the spherical triangle.
Picking any two of the three vertices of the sphereical triangle, one finds that the vectors pointing from the sphere center to these two vertices form a plane that contains the great arc connecting these two vertices. For example, let us focus on the plane containing the sphere center, vertex 1 (the north pole), and vertex 2 (the point with coordinates (θ = β 3 , ϕ = 0)).
Vectors from the center of the sphere to any point on the line can be extended to the surface of the sphere, the collection of which forms the great arc connecting vertex 1 and vertex 2 on the spherical surface, see figure 1 as an illustration. Any point in the interior of the flat triangle is thus mapped to a point inside the spherical triangle via extending vectors originating from the sphere center to points on the flat triangle further outward to the sphere surface. Evidently, this mapping is one to one and continuous, meaning that close points on the flat triangle are mapped to close points on the sphere.
We next build a smooth function on the flat triangle and then extend it to the spherical triangle to achieve a continuous charge distribution over all spherical triangles. Assume that we assign a density value σ i to vertex i. A simple way to build a smooth function is to interpolate linearly on the flat triangle. For any point p inside the flat triangle, let us extend the line from vertex 1 to p further till it reaches the edge L 1 connecting vertices 2 and 3, and let us call the intersection point p′. See figure 2 for a graphical illustration. Let the line length from vertex 1 to p′ be ℓ and the line length from vertex 1 to p be tℓ (with 0 ⩽ t ⩽ 1).
Let us also use L 1 to represent the length of edge L 1 . Now assume the distance between vertex 2 and p′ to be L 1 λ (with 0 ⩽ λ ⩽ 1). It is natural to assign the charge density (1 − λ) σ 2 + λσ 3 to point p′. Similarly, for point p, we may assign the charge density to be (1 − t) ] Note that at the center of mass of the flat triangle where λ = 1/2 and t = 2/3, we find that the interpolated density takes the desirable value (σ 1 + σ 2 + σ 3 )/3.
With all the needed components developed, we may now write down the interpolated charge density profile on the spherical triangle as where t(θ) is given by equation (11) and ⩽(ϕ) is given by equation (12). Note that for a given ϕ angle, θ ranges from 0 to β(ϕ), which by equation (1) is with ϕ 0 given by equation (6) and a by equation (7).
Given the coordinate choice of the three vertices (with one being the north pole of the sphere) of the spherical triangle and the charge interpolation scheme above, we may now compute the electric field produced by the spherical triangle at any arbitrary point. We separate three cases: first, the potentials right above and below the north pole; second, the electric field at any point on the same sphere but outside the spherical triangle; third, the electric field at any point outside the sphere.

Isolating electric field discontinuity right above and below the north pole
For the first case mentioned above, the surface charge density right at the north pole produces electric field of opposite direction right above and below the north pole. This suggests that one may calculate directly the electric field right at the north pole and incorporate the discontinuity of electric field afterwards. To illustrate this idea, let us consider the electric field produced by a spherical cap (parametrized by 0 ⩽ θ ⩽ β, 0 ⩽ ϕ ⩽ 2π) with uniform surface charge density σ. In the Cartesian coordinate system, a point near the north pole has coordinates (0, 0, R + ϵ), with ϵ > 0 (ϵ < 0) indicating point right above (below) the north pole. Of course, the north pole location is at ϵ = 0.
Consider the electric field produced by the spherical cap. Due to azimuthal symmetry, the net field can only be along the z direction: with s ≡ R/(R + ϵ). Since we are interested in the limit |ϵ| → 0, we consider the limit s → 1. When s → 1 from the s < 1 end, it means ϵ > 0, corresponding to the case when the point is right above the north pole. When s → 1 from s > 1 end, it means ϵ < 0, corresponding to the case when the point is right below the north pole. Continuing the calculation, we find Note that when ϵ = 0, the integral (16) superficially reduces to 2πσ (1 − cos β) ∕ 2, which is exactly the first contribution in the equation above. The second (nontrivial) contribution in the equation above arises from the difference between the electric fields right above and below the north pole. It contributes 2πσ (−2πσ) for the point immediately above (below) the north pole. This indicates that we may obtain the correct field strength by naïvely computing the field right at the north pole and add the discontinuous part afterwards.

Electric field near the north pole (a vertex of the charged spherical triangle)
Without loss of generality, we can rotate the coordinates such that the vertex of interest (near which the electric field should be calculated) becomes the north pole, and the other two vertices have coordinates (θ = β 3 , ϕ = 0) and (θ = β 2 , ϕ = α). Obviously, the great arcs connecting the north pole to both (β 3 , 0) and (β 2 , α) are longitudinal lines. The great-circle arc connecting (β 3 , 0) and (β 2 , α), on the other hand, follows the parametric form (1) with parameters specified by equations (6) and (7) in section 3.1.
Using the notation and development in section 3.2, we may write the continuous portion of the electric field along the z direction at the north pole (vertex 1) as with t(θ) and λ(ϕ) given by equations (11) and (12) respectively. The θ integral in When the β 2 and β 3 angles are small, the parameter a in equation (1) is large, implying that β(ϕ) is also small. Consequently, one may expand 1 − cos β(ϕ) using and In fact, the above expansion can be evaluated to arbitrary accuracy since 1/a is of order β (see equation (7))and the denominator of every term contains odd powers of a cos(ϕ − ϕ 0 ) whose integral can be calculated with the understanding that both (−1)!! and 0! are defined to be 1.
We now turn to the integral containing t( because both 1 − cos β and δ/R are (β 2 ). The θ integral in equation (18) can now be evaluated: The first part of integral (22) can be rewritten as (with η 2 ≡ (C 2 + D 2 )/C 2 ) Note that η cos β remains greater than one for small β and δ/R. This point can be seen as follows. Testing if η 2 cos 2 β > 1 leads to testing the positivity of which is easily satisfied for small β and δ/R. One then continues with the remaining integral Therefore, the first part of integral (22) (the total contribution from integral (23)) can be written as The second part of the integral (22) can be done in a similar fashion Combining integrals (24) and (25), we have the full results for the θ integral (22) By multiplying the result above by [(σ 2 − σ 1 ) + λ(ϕ)(σ 3 − σ 2 )] and integrating over the ϕ angle from 0 to α, we obtain the second half of the integral (18). Prior to doing so, let us further analyze δ(ϕ) and η to obtain the expansion scheme. We already know that for a small spherical triangle (hence small β) δ/R is of order β 2 , η is of order 1/β, C is of order β 2 and D is of order β. Therefore, η 2 C is of order 1. We expand expression (26) as follows.
The first and the second terms inside the curly brackets can be expanded as Similarly, the terms inside the square brackets can be expanded as follows: Combining results from (26)-(28), we may write In order to perform the ϕ integral, we need to express sin β 2 and cos β 2 in terms of ϕ. Using cot β = a cos(ϕ − ϕ 0 ), we have This confirms that 1/(a cos(ϕ − ϕ 0 )) is indeed of order β. Similarly, There are a few more details we need to attend to. If we examine equation (29), we notice that there is an order 1 combination η 2 C that needs to be worked out and expanded. Also, terms inside the square brackets of equation (29) properly. Since the order β 2 small quantity δ/R appears in C(ϕ) and D(ϕ), we also need to express adequately the leading behavior of δ/R to facilitate the proper expansions needed. We shall name the order β 2 small quantity 1 − r 2 ⋅ r 3 as ζ 23 . To simplify our expansion, without loss of generality, we shall introduce an order 1 ϕ-dependent quantity We now further examine the coefficients C(ϕ) and D(ϕ). and leading to We thus have (1 + 2h ) 2 1 + 4h 2 − 8h + 3 4a 2 cos 2 (ϕ − ϕ 0 ) + 28h 2 − 20h − 1 16a 4 cos 4 (ϕ − ϕ 0 ) + … .
(33) Therefore, we have Also, the order 1 quantity 2/η 2 C may also be expanded as and finally the expression inside the curly brackets of equation (29), after some tedious algebra, may be expanded in powers of 1/a cos(ϕ − ϕ 0 ) as We may now rewrite the expression (29) as ∫ 0 Note that the order 1 variable h depends on the azimuthal angle ϕ. The expression above will be multiplied by [(σ 2 − σ 1 ) + λ(ϕ)(σ 3 − σ 2 )] and then integrated over ϕ from 0 to α. In equation (34), we have kept terms all the way up to β 7 . However, as a proof of principle we will only keep terms up to β 3 in the remainder of the manuscript. (Although the systematic approach to reach higher order terms will be described, we note that the ϕ integrals can be quite tedious for higher order terms.) Given h ≡ ζ 23 λ(ϕ)(1 − λ(ϕ))a 2 cos 2 (ϕ − ϕ 0 ), integrals of the following sorts (ϕ dependence of λ suppressed) must appear with d = 0 or 1, p an odd positive integer, and the nonnegative integer ℓ < p/2. The case of ℓ = d = 0 has been worked out explicitly in (21). Therefore, we concentrate here on ℓ + d ⩾ 1, meaning the integral exponent of λ is positive.
From equation (12), we see that both λ and 1 − λ are rational function of sin ϕ and cos ϕ. If one substitutes tan and consequently, We also write cos(ϕ − ϕ 0 ) = cos ϕ 0 cos ϕ + sin ϕ 0 sin ϕ = cos ϕ 0 It turns out that tan ϕ 0 and (A + B 2 )/B 1 are very close when β ≪ 1. From equation (4) tan Below we will illustrate the basic idea how the order of β are incorporated in ℱ(ℓ, d, p) when By writing τ′ = τ + (τ′ − τ), we may expand part of the integrand of the above integral as where the negative binomial coefficient C k −n (with n ⩾ 1) is a short hand for Consequently, With (τ′ − τ) being of order (β 2 ), the above expansion of ℱ provides a systematic inclusion of higher order β contributions. We should note that the factor 1 − x 2 + 2τx is of order 1 and is always large compared to β in the range 0 ⩽ ϕ ⩽ α. Given that f(   However, we also note that Because the integrands are all rational functions of x, they can be integrated exactly using, for example, the Ostrogradskiy-Hermite method. In appendix C, we list some results of these integrals. Note that σ i − σ j (when i and j are neighboring vertices and when the spherical triangles are small) should be of order β based on the charge distribution continuity. (We may imagine Taylor expanding the charge distribution around a fixed point (say vertex i of the triangle considered.) Using equations (19) and (34), the electric field contribution in equation (18) accurate up to order (β 5 ) can be written as Looking up the integral results from appendix C, one sees that the expression becomes complicated as one goes to higher order terms. However, we have made an effort to keep most of the results expanded in powers of cos ϕ 0 and cos(α λ ϕ 0 ) although the expressions may still contain a single power of sin ϕ 0 and/or sin(α − ϕ 0 ). This choice is for computational convenience as one only need to obtain cos ϕ 0 (or cos(α − ϕ 0 )) once and then store it and raise it to any power as needed. Even though one may in principle keep as many higher order terms as one wishes, for illustration purpose, we will from this point on only keep to order β 3 accuracy as previously mentioned. By looking at the development up to this point, the readers would have gotten the idea of how to keep higher order terms when needed. Remember that (ℓ, d, k, p) is of order β p+2k and that (σ 2 − σ 1 )/σ 1(2) is of order β, we therefore have (with accuracy up to order β 3 ) We can thus compute the field contribution of all spherical triangles containing the north pole as the one of the vertices. Once contributions of all those spherical triangles are summed, we still need to remember to add 2πσ 1 to the field right above the north pole and −2πσ 1 to the field right below the north pole.

Radial electric field on the same sphere but outside the charged spherical triangle
Now consider the radial electric field, produced by the spherical triangle discussed earlier, at a point outside the solid angle region of the spherical triangle. Let us denote the radial direction by the polar angles (θ, ϕ). For α < ϕ < 2π, θ can take the full range 0 ⩽ θ ⩽ π. However, when 0 ⩽ ϕ ⩽ α, the range of θ becomes where r < = min(r , r), r ≡ r ∕ r, and in spherical polar coordinates r ≡ (θ, ϕ).
Electric potential produced by the spherical triangle (with 0 ⩽ ϕ ⩽ α and 0 ⩽ θ ⩽ tan −1 The electric field at the spherical surface will be needed when enforcing the boundary condition. If we take the derivative by assuming first r > R, we obtain If we take the derivative by assuming first r < R, we obtain These two expressions suggest that This is correct since direction of r is outside the spherical triangle while the direction of r is within the spherical triangle. Therefore, r ≠ r always. This indicates that we can actually average the two expressions for E r above to arrive at The correctness of the expression above can be verified by direct evaluation. Let us consider the radial field produced at location r by a charge q t at r: q t 2R 2 1 2 − 2r ⋅ r ni exact agreement with (43). This is a great simplification as the radial electric field expression (43) happens to be 1 2R times the electric potential.
In the previous subsection, the electric field at the north pole was computed using β (the range of polar θ angle) as the smallness parameter. In fact, the leading term in electric field is of order β and we worked out its leading correction explicitly to (β 3 ), that is β 2 smaller than the leading term. Our case here is slightly different. Previously, we have r z and thus presents some difficulty for integration over θ and ϕ. Specifically, there is no clean way to perform an expansion. For example, although we know for sure −1 < r ⋅ r < 1, expanding 1 ∕ 1 − r ⋅ r using r ⋅ r as the expansion parameter can be extremely slow converging when θ is small because the integration domain is an area close to and include the north pole. To draw a close analog to the previous case, we write  To proceed with the expansion, we need to first find the average r ave = (sin θ cos ϕ, sin θ sin ϕ, cos θ) ≡ (t x , t y , t z ) over the integration domain ∫ 0 Before we perform the averages one by one, let us first find out the area of the domain of integration, which we will need when computing the average. (Let y = sin(ϕ − ϕ 0 ) be the substitution variable below.) We now calculate the average below sin θ cos ϕ = ∫ 0 α cos ϕ dϕ ∫ 0 where equations (8) and (9) have been applied to simplify the expression above. Similarly, Let us define the following notations: t x ≡ sin θ cos ϕ, t y ≡ sin θ sin ϕ, and t z ≡ cos θ. Since , r ⋅ [r − r ave ] = sin θ cos ϕ(sin θ cos ϕ − t x ) + sin θ sin ϕ(sin θ sin ϕ − t y ) + cos θ (cos θ − t z ) ≡ sin θ cos ϕ(r x − t x ) + sin θ sin ϕ(r y − t y ) + cos θ (r z − t z ) .
If we keep up to the third term (quadratic in r ⋅ [r − r ave ]), we will need all the following integrals evaluated Here, we can in principle expand the expression to higher order β to gain higher order accuracy. However, as we mentioned before, we only keep correction terms of relative order β 2 to the leading term, which itself is of order β 2 for q 0 . (This is because the area of the spherical triangle is proportional to β 2 .) Expression (C.1) containing higher order terms are given in appendix C. We thus expand (50) to ∫ 0 α dϕ ∫ 0 β(ϕ) sin θdθ (sin θ cos ϕ, sin θ sin ϕ, cos θ) σ(θ, ϕ) .
So in terms of the θ integral, we need the following integrals Both integrals above again can be expanded to higher orders (see appendix C), but we again only keep correction terms of relative order of β 2 to the leading terms. Integral (52) can be evaluated to yield  The integrals needed are provided in appendix B. The integral (53) can also be evaluated to yield which upon evaluating ∫ 0 α dϕ yields 〈r z 〉 q 0 The integrals needed are provided in appendix B.
The calculations shown in the previous section and their extensions provide the formulas needed and we thus will not repeat them here.

Discussions and conclusions
The reader might have noticed that we did not discuss the possible singularities [13,14] that may arise from the molecular surface definition of Richards and Connolly. This is because we are mainly interested in applying the classical formalism here to the explicit solvent model where we are allowed to make the probe sphere (whose sole functionality is to provide a smooth surface for applying boundary condition) small enough to avoid these potential singularities.
In this manuscript, we lay the first part of the ground work for employing curved patches for applying surface charge method in electrostatics. We have analytically shown how one may control the accuracy by expanding in powers of the arc length (multiplied by the curvature). Before we conclude this paper, let us illustrate with the computation of the area of a spherical triangle using the series expansion method outlined here versus the exact answer. Let us consider a spherical triangle with one vertex on the north pole, (0, 0, 1), the second vertex having coordinates (sin β, 0, cos β) and the third vertex (sin β cos α, sin β sin α, cos β) with β ≪ 1 and α = π/3. We then have ϕ 0 = π/6 and a = 2 3 cot β. We have evaluated the area exactly in equation (46)  In the forthcoming manuscript, we shall cover the other type of surface, the toroid, along with a proper way to deal with triangles with one side not belonging to a great arc, which occur at the jointing region where spherical surface and toroidal surface meet. In addition, we will provide a discussion on the disappearance of divergence of the electric field when curved surface patches, developed in this paper and the forthcoming paper, are employed. We will also, in the third paper of the series, apply the fundamental work developed here and in the forthcoming paper to investigate molecular interaction at short distance. In this case, because the electric field experienced by a molecule may vary significantly from one of its atoms to the others, approximating the molecule as a single dielectric sphere with multipole moments may be insufficient; using the more general molecular shape, composed of partial toroids and spheres, may be better suited to deal with cases where external electric field has non-negligible variation across a molecule.    The following are some identities we used to somewhat symmetrize cos ϕ 0 and cos(α − ϕ 0 ) in the integrals above.