Potential applications of solute-solvent interactions include synthetic medical design, charge transfer, simulation of membranes and nanotubes as well as studies of biological process. The molecular surface is surrounded by the solvent media consisting of mobile ions admitting a specified permittivity. The domain enclosed within the molecular surface consists of the solute composed of fixed atoms admitting their respective charges   . We focus on the Boundary Element Method (BEM)     using wavelets on NURBS (Non-Uniform Rational B-Splines) patches     to solve the Polarizable Continuum
Model (PCM) which treats a particular case of the Poisson-Boltzmann equation when the coefficient related to the modified Debye-Hückel parameter vanishes. Approaches based on BEM    become a preferred tool for solving linearized solvent problems because it is inconvenient to apply 3D-solvers such as FEM (Finite Element Method)    for several reasons. For one, the related transmission problem is solved by FEM in the entire 3D-space    while only the solution on the molecular surface is required. In addition, most FEM programs and analysis assume that the right hand side (RHS) is a sufficiently smooth function whereas one considers in the solvent problem a RHS which is a sum of nonsmooth Diracs defined only in the sense of distribution  leading to very fine adaptive mesh refinements. In the opposite, FEM appears to be the only recourse to treat nonlinear Poisson-Boltzmann   . The principal BEM unknown is the apparent surface charge that is not a physical entity but a fictive value which behaves very much like a density of charge distributed on the molecular surface. In the application, the entire solution to the original transmission problem of the PCM is not needed as only the apparent surface charge suffices to deduce the electrostatic reaction potential. The initial transmission partial differential equations are recast as a PCM integral equation involving the single layer and the double layer operators. That will be decoupled into a pair of equations which are solved separately. First, one solves a second kind integral equation governed by the double layer operator only. The solution to that first problem serves subsequently as a RHS to another equation involving the single layer operator. The linear systems from BEM     are fully populated in the case that the classical polynomial basis is utilized. Wavelet basis   serves well as an efficient technique to obtain quasi-sparse matrices  . The construction of the wavelet basis    starts from univariate basis on the unit interval. Taking the tensor product in the unit square and mapping onto the patches result in the wavelet basis on the whole manifold. We assume the absence of anisotropic patches where some curved edges are relatively much longer in comparison to other ones. In addition, the surface areas of the NURBS patches ought not be considerably dissimilar. Accomplishing such a geometric property is not straightforward if the only available information consists of the nuclei coordinates and the Van-der-Waals radii of the atoms. Formerly, we have made efforts to achieve a geometric structure resulting in patches admitting comparable surface areas as well as isotropically shaped structure. The inputs are usually acquired from the format PDB or PQR for which some conversion techniques exist  . The system for the double layer operator produces a bounded condition number and no preconditioner is required. In the opposite, the single layer potential is severely badly conditioned and a domain decomposition technique  is used as a preconditioner for the reduction of the condition number which substantially decreases the required number of iterations to solve the linear system iteratively. Before proceeding further, we briefly survey our previous papers as well as computer implementations. A splitting method for CAD surfaces has been
proposed in   for BEM simulation. Additionally, methods for checking the regularity of the mappings have been proved in  . The surface structure, which is required by the wavelet-BEM, is unfortunately very complicated to construct in contrast to the standard mesh generation  . While approximations are required to obtain global continuity in  for CAD objects, it can be achieved exactly for molecular surfaces in  . Furthermore, some quantum simulation by using wavelet BEM on single processor is described in  . A wavelet BEM simulation using domain decomposition techniques was described in  which treats the case of ASM (Additive Schwarz Method). By using Sobolev norms with negative orders  , it was utilized as an efficient preconditioner for the wavelet single layer potential which is badly conditioned. Recently, we gained experience  in elasticity nanosimulation where one has used nanotube fibers immersed in polymer matrices as quantum models. Separate studies dedicated exclusively for the wavelet double layer potential are documented in  which describes as well the modeling techniques to achieves the molecular NURBS patching. In this current paper, we mainly compare computational results with analytical solutions. We display numerical outcomes pertaining to the single and double layer operators separately for complete molecules. The only well-known molecular model that has explicit analytical formula for the solvation energy appears to be the single atom, which is also known as Born ion. More precisely, it consists of an atom admitting specified charge and radius which is surrounded by the solvent media with given permittivity. More details pertinent to chemical solvation energy for more complex molecules will be deferred to a subsequent paper.
2. Transmission Problem on NURBS Patches
We present in this section the required geometric structure needed by the wavelet formulation of the PCM. The problem setting of the transmission problem which is related to the PCM will also be developed. The pertaining integral equation formulation is introduced because we need only the values on the interface . We consider atoms which are located in the solute region and we suppose the geometry satisfies the following conditions:
・ We have a covering of the surface by four-sided patches ,
・ The intersection of two different patches and is supposed to be either empty, a common curvilinear edge or a common vertex,
・ Each patch where is the image by which is described by a bivariate function that is bijective, sufficiently smooth and admitting bounded Jacobians,
・ The patch decomposition has a global continuity: for each pair of patches , sharing a curvilinear edge, the parametric representation is subject to a matching condition. That is, a bijective affine mapping exists such that for all on the common curvilinear edge, one has . In other words, the images of the functions and agree pointwise at common edges after some reorientation,
・ The manifold is orientable and the normal vector is consistently pointing outward for any ,
・ The patches admit similar surface areas.
An illustration of the above surface structure is depicted in Figure 1. The CAD representation of the former mappings uses the concept of B-spline and NURBS    . Consider two integers such that . The interval [0, 1] is subdivided by a knot sequence such that for and such that the initial and the final entries of the knot sequence are clamped and . One defines the B-splines    basis functions as
where we employ the divided difference in which we use the truncated power functions given by if , while it is zero otherwise. The integer k controls the polynomial degree of the B-spline which admits an overall smoothness of while the integer n controls the number of B-spline functions for which each B-spline basis is supported by . The NURBS patch admitting the control points and weights is expressed as
which describes thus a parametric function from the unit square onto the four-sided patch . The similarity of the surface areas are needed in practice for the wavelet threshold to function in a unified manner. We
Figure 1. Patch representation of a Water Cluster with 1089 NURBS.
will the wavelet threshold to function in a unified manner. We will consider only geometries which are globally smooth and which admit moderate curvature. For each patch , the Gram determinant is denoted by
which quantifies the norm of the normal vector at the image of any point belonging to the unit square by the parametric function . The Gram determinant and its partial derivatives are assumed to be bounded
for where for sufficiently large. We consider the transmission problem for a parameter related to the dielectric coefficient:
where denotes the electric charge of the k-th atom while is the Dirac distribution centered at the coordinates of the nucleus . It consists of a special case of the Poisson-Boltzmann equation in the situation that the effect of the Debye-Hückel parameter is neglected. We are not solving those equations directly, rather we solve only the pertaining integral equations which are located on the interface surface . Introduce the single layer , the double layer and the adjoint double layer operators for as follows
Consider the components of the solution u inside the solute and the solvent respectively. We recall very briefly the transformation of the transmission problem (5)-(9) into integral equations which follow the procedure of   where one replaces by such that
The apparent surface charge is defined as so that one obtains
which shows that behaves as a charge distribution over while the second term is a sum over the actual charges . The representation formula  yields
where represents the conormal derivative. By combining them with the transmission jump conditions and , one obtains
of which the full detail is described in   . The fact that , and the Newton potential lead to
Apply the representation formula to the Newton potential which is harmonic to obtain which results in
In practice, we are not interested in the entire solution of the transmission problem but only in the electrostatic reaction potential
Thus, our objective is to seek for the apparent surface charge satisfying
where and are constant parameters while m and g are given functions which are sufficiently smooth. For the practical applications related to (5)-(9), the coefficients and and the functions m and g are as follows
In order to simplify the presentation, we assume henceforth the parameters and are unity. Everything remains valid with little modifications for other constant parameters. The above problem (19) is decoupled as two integral equations:
That is, the solution to (22) is used as a RHS in (21). To numerically treat those last problems by means of the wavelet technique, several approximations are involved:
・ The Galerkin error related to the single layer by using a finite dimensional space ,
・ Similar Galerkin error but for the double layer potential ,
・ Replacing the RHS for by the approximate solution from the double layer equation,
・ Replacing the operator by its truncated version ,
・ Similar truncation but for substituting by .
We will consider the space of piecewise constants whose construction will be specified later on. Concerning the discrete Galerkin variational formulation in for discretizing the integral Equations ((21) and (22)), one searches for such that for all
where and are respectively the kernel functions for the single layer and double layer . Most of the following theoretical works are inspired from different sources including          .
3. A-Priori Estimate for Multi-Wavelet PCM
We consider in this section the multi-wavelet Galerkin formulation of the PCM problem. We recall several important properties which are used in the deduction of the a-priori error estimate when no wavelet. For the unidimensional basis functions, we introduce the next knot sequence on level ,
The internal knots on the next level are obtained by inserting a new knot in the middle of two consecutive knots on the lower level . Introduce the piecewise constant linear space in the unit interval on level :
where designates the characteristic function with respect to . By using the two scale relation
and the inclusion , the spaces form a nested sequence of subspaces:
On account of the nestedness (28), the space is expressed as an orthogonal sum
with respect to the -scalar product where is the complementary wavelet space
For the explicit expression of the wavelet functions , we use the Haar wavelet defined on by
whose relation with the single scale basis is . By using dilation and shift, one defines for and
The wavelet functions constitute an orthonormal basis
where the first Dirac comes from the inter-level orthogonality while the second Dirac is justified by the non-overlapping of and on the same level . By applying the decomposition (29) recursively, one obtains on the maximal level
so that we have the dimensionalities
Thus, a function has two representations: in the single-scale basis and in the multiscale basis, we have respectively
From here onward, we use the notation if there is a constant c such that in which is independent on the level and the maximal level. Besides, the notation amounts to . Moreover, we denote
when the expression of is complicated. It has the properties:
The next norm equivalences related to the coefficients are valid  
Due to the property (33) and , the orthogonal projection of any onto verifies
The 2D-wavelet space on the unit square , is defined for any level in term of tensor product as follows
On each patch ( ) and on the whole surface , we define for the level
The space corresponds to the piecewise constant functions on a tensor product mesh admitting a step size of . It is deduced from the above construction that the next nested inclusions are valid
Some immediate properties for all wavelet functions are as follows:
They are easily derived from the 1D definitions where and the boundedness of the mappings .
The space of square integrable functions is
which is equipped with the following scalar product and norm after transformation onto ,
The Sobolev space on for a non-negative integer is
where the differentiation is interpreted in the sense of distribution  such that for all compactly supported smooth functions . The Sobolev space is endowed with the norm
Concerning a real positive order such that is an integer and , the Sobolev space consists of the functions such that their norms with respect to the next Slobodeckij norm are finite
For negative orders, one employs the dual spaces equipped with the dual norm
We will denote the orthogonal projection with respect to the -scalar product onto by such that
Theorem 1. On the 2D level , consider the discrete Galerkin PCM problems for seeking such that
Then, the accuracy between the apparent surface charge from (21) and from (63) for and satisfies
In addition, the accuracy of the reaction potential from (18) and
Proof. The single layer operator can be decomposed  as where the principal part is an elliptic operator and is a compact operator. Therefore, it is a Fredholm operator of zero index which is in addition injective  . As a consequence, Equations ((21) and (63)) are uniquely solvable. On account of the Strang lemma  with perturbed right hand side, the orthogonal projection leads to
On the 2D-level for piecewise constant setup, one has
In addition, apply the inverse estimate to obtain
A combination of the above relations yields therefore
On the other hand, for the estimation of , one uses the same reasoning as above where is now the compact operator. Hence, one has the estimation
The a-priori error estimate for the apparent surface charge follows
Concerning the reaction potential, one obtains on account of the Cauchy- Schwarz inequality
We introduce the minimal distance between the nuclei coordinates and the molecular surface ,
By using the orthogonal projection and an inverse inequality, one obtains
As a consequence, one concludes the following accuracy by applying (74),
4. Multiscale for Solvation Energy
This section will be occupied by the treatment of the PCM by means of the wavelet technique. The matrix entries where the distance between the supports is beyond a level dependent threshold are annihilated. For a matrix , define
For two levels , a threshold whose expression will be specified subsequently is used. We consider the matrix coefficients
They correspond to the matrix entries in (23) and (24) for the single layer kernel and the double layer kernel after transformation onto . The following level depended truncation procedure for each pair of levels is applied:
The block matrices and have respectively the entries and such that we have blockwise:
For the single layer and the double layer, define respectively
For the involved kernels, one has on account of the Calderon-Zygmund inequality:
Lemma 1. Fix any real parameters and consider two levels . Suppose (see (39), (40),(41))
where and denotes the diameter of the supports of any wavelet on level . Then, the following error estimates in matrix norms are fulfilled:
By using the thresholds for all as in (88) such that the parameters satisfy
one has for and in the case of the single layer operator
while the double layer operator yields
Proof. The value of from (83) will be examined on a pair of patches . The indices are dropped to simplify the notation. Consider a tensor product wavelet basis and fix as the midpoints of and . By considering any , the Taylor expansion leads to
for some . For the first term, by taking the integration over , one obtains
which is zero due to (33). We will consider now the value of in the case . For the summation over from (95), supposing that the Jacobians of the mappings and are bounded, one has for and
Since , one has the next inequalities for and :
which, when inserted into (97), leads to
A combination of (95), (96) and (101) leads to
On account of the properties (51)-(54), one deduces
The function is continuous on the compact . Hence, the value
is bounded independently of . As a consequence, use the expression of to obtain
As a consequence, one deduces the norm estimate
For any , we have
where the last equality uses a local mapping and polar coordinates. By using the measure property (51), obtain . The 1-norm is processed in an analogous manner and therefore
Consider the orthogonal projection onto and introduce
Blockwise according to the levels
By utilizing the next estimates
one obtains, based on the Cauchy-Schwarz inequality and (91),
The use of the sum of geometric sequence leads to
For piecewise constant setup, obtain for the constant parameters
For the single layer , set and use to obtain from (86)
Similarly for the double layer, setting results in
Theorem 2. Consider the truncated PCM equations on level :
where the operators and are obtained from and by using the threshold and respectively.
Suppose the constant parameters are chosen such that
For patches, the numbers of nonzero matrix coefficients are then reduced from to and .
In addition, the accuracy between in (21) and in (117) for and verifies
Further, the error between the reaction potential from (18) and
Proof. We consider first the number of nonzero entries of the compressed matrices and . For computing for each fixed level pair , consider a basis function on level . According to  , the number of basis functions on level whose support is within a distance of from the support is in the worst case . For that, Lemma 3.2. of  uses an argument using a sphere of radius centered at a point located upon the support . There are basis functions on level for each patch where . The number of the interlevel nonzero entries on all patches is thus
By summing over all levels
The first sum on all patches verifies
For the second sum , the expression of the threshold leads to
By using sum of geometric sequence, deduce from
In a similar manner, one obtains
Hence, obtain for the second sum by using
Hence, for all patches. Concerning the first Equation (117), the next Strang relation  holds
By inserting (115) in the last estimate, the next inequality is fulfilled
In the same manner as (70), one has
On the other hand, for the second Equation (118) by applying the Strang lemma  for the second time
By using (116), deduce
Based on the last estimate and (127), deduce
Concerning the accuracy with respect to the electrostatic reaction potential , one proceeds as in (75) and (76) in order to obtain
5. Practical Implementation and Numerical Results
We want to present in this section some results of the former method. We will describe first the process of obtaining the computational results. We distinguish
Figure 2. Relative errors for several RHS in function of levels: (a) double layer, (b) single layer.
separate tests for the single layer operator, the double layer operator and the electrostatic reaction potential. The exact solutions which have zero Laplacians are specified by the user in the case of the double/single layer potentials. The corresponding right hand sides are then obtained by applying the Lapacian operator to the exact solutions as illustrated in Figure 2(a) and Figure 2(b). Transforming Laplace problems into integral equations using double/single layers is found in standard textbooks of integral equations   . The assembly of the basis functions consists in constructing piecewise constant wavelets on the unit interval. After taking the tensor product, one obtains basis functions on the unit square which are mapped by the NURBS functions onto the patches that are embedded in the 3D-space. The assembly of the matrix entries uses hierarchical integrations which are described in our former paper  . That procedure can be achieved on arbitrary geometries. The linear system obtained from the double layer operator does not require any preconditioner as a direct GMRes linear solver suffices to solve the system. For the single layer operator however, we use a domain decomposition preconditioner to reduce the number of iterations. The error computation consists in comparing the chosen exact solution with the outputs which are acquired by applying the BEM-simulation. The reason for displaying separate results pertaining to the double layer potential and the single layer potential is that if they are solved individually, analytical results for arbitrary boundary exist for comparisons. As for the computation of the electrostatic reaction potential, we make use of the Born ion. That is, the only well-known analytical result for the reaction potential is a simple geometry composed of a single atom. For that case, the process consists in fixing some configurations by setting the values of the radius, the electric charges as well as the dielectric parameter. The actual process consists therefore in applying the PCM simulation on a resolution which is specified by the wavelet level. The electrostatic potential is obtained by traversing the patches and by computing some integrals related to the apparent surface charge as specified in (18). For the numerical computation of those integrals, a standard Gaussian quadrature for smooth integrand suffices because the nuclei are not located on the NURBS patches . Different parameters are varied in order to validate the accuracy of the PCM-simulation.
5.1. Molecular Models
We employ different sorts of quantum models including molecules which are acquired from PDB/PQR files  as well as water clusters which are obtained from a former MD (Molecular Dynamics) simulation. When the MD-iteration attains its equilibrium state where the total energy becomes stable, a water cluster is obtained by extracting the water molecules which are contained in some given large sphere whose radius controls the final size of the water cluster. The formerly proposed method assumes that the sizes of the patches are proportional. Efforts were done to obtain patches whose shapes related to the normal derivatives are as good as possible. In order to measure the quality of the patches, we examined the proportionality of the area, length of the curved edges and the norm of the normal vectors. That is important because we use the same wavelet threshold for all patches. We avoid a patch where some curved boundary edges are very long while others are very short. For each patch , the quality gauge for the area is
The average of the above value is collected on Table 1 which shows that although the areas are different, their ratios with the average area are not excessive. Similar quality measurements as (131) have been utilized for the length of curved edges separating the patches  as well as the normal vectors which are computed at some tensor product samples. It is observed from Table 1 that all involved entities in the patches are practically proportional.
5.2. Double Layer Potential
For the double layer potential, one can solve the linear system iteratively without recourse to any preconditioner because the system is well conditioned. A GM-Res method serves as a solver of the system which is non-symmetric. We
Table 1. Quality of the NURBS patches for several molecules.
Table 2. BEM accuracy for the double layer and the single layer in function of the levels.
shall examine first the error reductions in term of the multiscale level. The results are collected in Table 2 which contains both the absolute errors and the relative errors. The absolute errors are obtained by comparing with the exact solution at some fixed samples in the interior . A division by the largest value of the exact solution provides the relative error.
In the first part of Table 2, we collect the convergence of the BEM simulation in function of the level ranging from 1 to 4. Those results have been obtained from a water cluster molecule containing 386 patches. The ratio between the nonzero counts of the compressed and uncompressed matrices are also exhibited. We examine now the general linear characteristic of the relative errors. We display in Figure 2(a) the error evolution where we consider five levels applied to a propane molecule using several right hand sides. We consider five exact solutions which are respectively and , , and that all have vanishing Laplacian. The right hand side is the restriction of the function on the boundary . The propane molecule admits 75 patches. The error reduction is lightly affected by the used right hand side but in general the errors reduce satisfactorily in function of the wavelet levels. In fact, they decrease linearly in logarithmic scale in function of the BEM levels.
5.3. Single Layer Potential
The single layer yields a linear system which is badly conditioned  . We use an additive domain decomposition preconditioner to remedy that bad conditioning. We briefly summarize the procedure whose comprehensive description is detailed in our former work  . In term of geometric structure, the overlapping domain decomposition is as follows
In term of linear spaces, this leads to the decomposition
Denote the orthogonal projection onto with respect to the bilinear form related to the single layer operator by
The ASM operator is defined by
The ASM procedure consists in solving a local problem which is projected by onto the subdomain . The results of the number of iterations are summarized in Figure 3(a) where it is observed that the ASM method needs significantly less iterations than the direct method. A larger view for the iterations less than 50 is exhibited in Figure 3(b). Like the double layer, the same test for the compression is detailed in the second part of Table 2 in the case of the single layer potential. Further, the decrease of the BEM error in function of the levels for different RHS is presented in Figure 2(b) where a propane molecule is used as in the double layer case.
5.4. Apparent Surface Charge
We want now to examine the values of electrostatic reaction potential which are computed with the PCM-wavelet technique. It is beyond the scope of
Figure 3. (a) Preconditioning the single layer potential for the entire iterations. (b) Zoom for preconditioning the single layer potential for iterations 0 - 60.
this article to provide a detailed comparison with other softwares as all results here are only computed for analytical comparison where we do not use any physical units. The explicit values for the reaction potential are unknown except for some simple geometries. We consider the Born ion that consists of a single sphere of radius which is centered at the origin and to which an electric charge q is assigned while the dielectric coefficient of the solute is adjusted. The analytical expression  for the reaction potential reads
. All those parameters are unitless in the sense that no physical units are used to measure them. Our next test consists in examining the effect of the dielectric coefficient . The corresponding result is depicted in Figure 4(a) in which we consider three configurations corresponding to , and respectively. We execute the wavelet-PCM simulation on multiscale level 3 while the dielectric coefficient varies from to . We observe for all three configurations that the exact reaction potential and the computed PCM results align well in Figure 4(a) where a zoom is also provided for in which the reaction potential drops down very quickly. As a next test, we investigate the variation of
Figure 4. Unitless comparison for the Born ion: (a) dielectric coefficients, (b) electric charge, (c) radius.
the electrostatic potential in function of the electric charge q. The corresponding outcome is depicted in Figure 4(b) where we consider four configurations corresponding to , , , while the dielectric coefficient is constantly . The electric charge is allowed to vary in the range . One observes in Figure 4(b) that the computed reaction potential varies quadratically in function of the electric charge. In fact, a parabolic dependence on the charge q is exhibited. The purpose of the next test displayed in Figure 4(c) is to highlight the agreement between the computed PCM result and the exact solution when the ion radius varies from to . One considers three configurations where the electric charge is respectively , and while the dielectric coefficient is consistently . For all three configurations, the electrostatic reaction potential is inversely proportional with regard to the ion radius. For all cases, the results computed by means of the PCM-wavelet align well with the theoretical expectations. As a final test, we examine the accuracy of the electrostatic potential in function of the BEM-levels which vary from 1 to 6. We consider two configurations where the radii and electric charges are and respectively. The results of the BEM accuracy are summarized in Table 3 where we consider four dielectric values , , and . The exact and the computed reaction potentials which are exhibited there imply that the BEM-approximation becomes satisfactorily precise as the multiscale levels increase.
Table 3. Computed PCM-reaction potential in function of the multi-scale levels for the Born ion.
We consider the transmission problem occurring in the interaction between the solute and the solvent media. Our method is based on an integral equation formulation which is solved by means of the wavelet Boundary Element Method. The entire molecular surface is supposed to be structured in four-sided NURBS patches onto which wavelets constructed on the unit square are embedded. For the Born ion, the exact reaction potential is known explicitly without using physical units. Our results align well with the exact solutions when we vary various parameters including the electric charge, the radius and the dielectric coefficient. The current approach is limited to the linear case where the coefficients in the transmission problem are constant parameters. Prospective future works could comprise the nonlinear Poisson-Boltzmann admitting nonsmooth Dirac load function. Further, non-constant matrix-valued coefficients which are difficult to treat by BEM might be examined. These challenging works need to be approached and presented step by step.