POLARIZABLE CONTINUUM MODEL FOR WAVELETS ON NURBS PATCHES

We are interested in the application of wavelet techniques on four-sided patches. The Polarizable ContinuumModel which is a special case of Poisson-Boltzmann equation is treated by means of the boundary integral equations. The patches are isotropically shaped and they admit similar surface areas. The integral equation is decoupled into two linear equations which involve respectively the single layer and the double layer operators. A wavelet thresholding technique is used to compress the linear systems into quasi-sparse ones. We examine the accuracy of the different approximations involved in the integral equations. Domain decomposition preconditioner severs as acceleration of the linear solver of the single layer which is badly conditioned.


Introduction
Potential applications of solute-solvent interactions include synthetic medical design, charge transfer, simulation of membranes and nanotubes, studies of biological process.We focus on the Boundary Element Method (BEM) using wavelets on NURBS (Nonuniform Rational B-Splines) patches to solve the Polarizable Continuum Model (PCM) which is a particular case of the Poisson-Boltzmann equation when the coefficient related to the modified Debye-Hückel parameter vanishes.BEM becomes a preferred tool for solving solvent problems because it is difficult to apply 3D solvers such as FEM (Finite Element Method) for several reasons.Most FEM programs and analysis assume that the right hand side (RHS) is a sufficiently smooth function whereas one has in the solvent problem a RHS which is a sum of nonsmooth Diracs defined only in the sense of distribution leading to extremely fine adaptive mesh refinements.In addition, the related transmission problem is solved by FEM in the whole 3D space while only the solution on the molecular surface is required.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 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 is subsequently used as a RHS to another equation involving the single layer operator.That process yields several numerical approximations.For one, the replacement of the RHS by the solution to the second kind integral equation.In addition, the Galerkin approximation which originates from solving the equation in a finite dimensional space which is the piecewise constant space in our case.For that matter, there exists the approximation obtained from substituting the original operators by their compressed counterparts.The linear systems from BEM [12,24,20,28,45] are fully populated by using the classical polynomial basis.Wavelet basis [10,27,12,11] serves well as a compression technique to obtain quasi-sparse matrices [5,25,26,34].By annihilating the small matrix coefficients in function of certain well chosen threshold values, only very few coefficients remain.The objective is to significantly reduce the number of nonzero matrix entries without substantially deteriorating the accuracy.The wavelet technique compresses the resulting matrix where the ratio of compression becomes better as the multiscale level grows.It is frequently observed in practice that at level four or five, only an order of 0.1 percent of the Galerkin matrices are essential to achieve the convergence.That thresholding approach would not function to anisotropic patches where some curved edges are relatively much longer in comparison to other edges.In addition, when the surface areas of the patches are considerably dissimilar, the thresholding technique would equally fail.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.We will briefly report on the results of such a decomposition.The system for the double layer operator has a bounded condition number and no preconditioner is required.In the opposite, the single layer potential is severely badly conditioned and a domain decomposition techniques is used as a preconditioner to reduce the conditioner which substantially decrease the required number of iterations to solve the linear system.Our past works: a splitting method for CAD surfaces has been proposed in [37] for BEM simulation.Additionally, methods for checking the regularity of the mappings has been proved in [39].While approximations are required to obtain global continuity in [39] for CAD objects, it can be achieved exactly for molecular surfaces in [19,38].Furthermore, a real chemical simulation by using wavelet BEM is described in [44] for the quantum computation.The surface structure which is required by the wavelet-BEM is unfortunately very complicated to construct in contrast to the standard mesh generation [41].In this work, we compare only the results with analytical functions.We display numerical results pertaining to the single and double layer

Transmission problem on NURBS patches
We present in this section the required geometric structure needed by the wavelet formulation of the PCM.We will see also the problem formulation of the transmission problem which is related to the PCM.The pertaining integral equation formulation is introduced because we need only the values on the interface Γ.We consider M atoms {x k } M k=1 which are located in the solute region Ω solute .We suppose the geometry Γ satisfies the following conditions.
• We have a covering of the surface by four-sided patches Γ = N p=1 Γ p , • The intersection of two different patches Γ p and Γ q is supposed to be either empty, a common curvilinear edge or a common vertex, • Each patch Γ p where p = 1, 2, . . ., N is the image by γ p : := [0, 1] 2 → Γ p 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 Γ p , Γ q sharing a curvilinear edge, the parametric representation is subject to a matching condition.That is, a bijective affine mapping Ξ : → exists such that for all x = γ p (s) on the common curvilinear edge, one has γ p (s) = (γ q •Ξ)(s).In other words, the images of the functions γ p and γ q agree pointwise at common edges after some reorientation, • The manifold Γ is orientable and the normal vector n(x) is consistently pointing outward for any x ∈ Γ, • The patches have similar surface areas.
An illustration of the above surface structure is depicted in Fig. 1.The CAD representation of the former mappings γ p uses the concept of B-spline and NURBS [35,37,21].Consider two integers n, k such that n ≥ k ≥ 1.The interval [0,1] is subdivided by a knot sequence τ = (τ i ) n+k i=0 such that τ i < τ i+1 for i = k − 1, ..., n − 1 and such that the initial and the final entries of the knot sequence are clamped One defines the B-splines [13,6,35] basis functions as where we employ the divided difference [τ i , τ i+1 , ..., τ p ]f in which we use the truncated power functions (• − t) k + given by (x − t) k + := (x − t) k if x ≥ t, while it is zero otherwise.The integer k controls the polynomial degree k − 1 of the B-spline which admits an overall smoothness of C k−2 while the integer n controls the number of B-spline functions for which each B-spline basis N τ ,k i is supported by [τ i , τ i+k ].The NURBS patch γ p admitting the control points d i,j ∈ R 3 and weights w i,j ∈ R + is expressed as The similarity of the surface areas are needed in practice for the wavelet threshold to functions in a unified manner.We will consider only geometries which are globally smooth and which admit moderate curvature.For each patch Γ p , the Gram determinant is denoted by The Gram determinant G p and its partial derivatives are assumed to be bounded for α = (α 1 , α 2 ) where |α| = α 1 + α 2 ≤ η for η sufficiently large.We consider the transmission problem for a parameter ε ≥ 1 related to the dielectric coefficient: where q k denotes the electric charge of the k-th atom while δ(• − x k ) is the Dirac distribution centered at the coordinates of the nucleus x k .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 V, the double layer K and the adjoint double layer K * operators for x ∈ Γ Consider the components (u i , u e ) of the solution u inside the solute Ω solute and the solvent Ω solvent respectively.We recall very briefly the transformation of the transmission problem (2.5) -(2.9) into integral equations which follows the procedure of [7,8] where one replaces (u i , u e ) by (w i , w e ) such that (2.10) The apparent surface charge is defined as σ := V −1 w i so that one obtains (2.11) which shows that σ behaves as a charge distribution over Γ while the second term is a sum over the actual charges {q k }.The representation formula [29] yields of which the full detail is described in [7,8].The fact that KV = VK * , w i = Vσ and the Newton potential Apply the representation formula to the Newton potential which is harmonic to obtain In practice, one is not interested in the entire solution u of the transmission problem but in the electrostatic reaction potential (2.18) 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 (2.5)-(2.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 (2.19) is decoupled as two integral equations: That is, the solution to (2.22) is used as a RHS in (2.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 V by using a finite dimensional space V L (Γ), • Similar Galerkin error but for the double layer potential K, • Replacing the RHS for Vσ = f by the approximate solution from the double layer equation, • Replacing the operator V by its compressed version V, • Similar compression but for substituting K by K.
We will consider the space of piecewise constants V L (Γ) whose construction will be specified later on.Concerning the discrete Galerkin variational formulation in V L (Γ) for discretizing the integral equations (2.21) and (2.22), one searches for σ where K SL (x, y) and K DL (x, y) are respectively the kernel functions for the single layer V and double layer K.

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 ℓ = 0, 1, ..., L 1D , (3.25) The internal knots on the next level (ℓ + 1) 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 [0, 1] on level ℓ: ℓ where χ D designates the characteristic function with respect to D. By using the two scale relation and the inclusion ζ ℓ ⊂ ζ ℓ+1 , the spaces V ℓ [0, 1] form a nested sequence of subspaces: On account of the nestedness (3.28), the space V ℓ [0, 1] is expressed as an orthogonal sum with respect to the L 2 -scalar product where W ℓ [0, 1] is the complementary wavelet space For the explicit expression of the wavelet functions ψ ℓ i , we use the Haar wavelet defined on [0, 1] by whose relation with the single scale basis is . By using dilation and shift, one obtains for ℓ = 1, ..., L 1D and i = 1, ..., The wavelet functions ψ ℓ i constitute an orthonormal basis (3.33) where the first Dirac δ ℓ 1 ,ℓ 2 comes from the inter-level orthogonality while the second Dirac δ i 1 ,i 2 is justified by the non-overlapping of Support(ψ ℓ i 1 ) and Support(ψ ℓ i 2 ) on the same level ℓ.By applying the decomposition (3.29) recursively, one obtains on the maximal level so that we have the dimensionalities Thus, a function u ∈ V L 1D [0, 1] has two representations: in the single-scale basis and in the multiscale basis, we have respectively From here onward, we use the standard notation X Y if there is a constant c such that X ≤ cY in which c is independent on the level ℓ and the maximal level.Besides, the notation X ≃ Y amounts to Y X Y .Moreover, we denote 2 X by exp 2 {X} when the expression of X is complicated.The next norm equivalences related to the coefficients are valid [12,34] (3.39) Due to the property (3.33) and The 2D-wavelet space on the unit square is defined for any level ℓ = 0, 1, ..., L := 2L 1D in term of tensor product as follows On each patch Γ p (p = 1, ..., N) and on the whole surface Γ = ∪ N p=1 Γ p , we define for the level ℓ = 0, 1, ..., L The space V L (Γ) 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 ψ (ℓ,i) ∈ W ℓ (Γ) are as follows: They are easily derived from the 1D definitions where ℓ = ℓ 1 + ℓ 2 and the boundedness of the mappings γ p .
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 k is where the differentiation ∂ α f is interpreted in the sense of distribution [43] such that for all compactly supported smooth functions g.The Sobolev space H k (Γ) is endowed with the norm Concerning a real positive order p = k + θ such that k is an integer and θ ∈]0, 1[, the Sobolev space H p (Γ) 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 We will denote the orthogonal projection with respect to the L 2 -scalar product onto THEOREM.On the 2D level L, consider the uncompressed discrete Galerkin PCM problems for seeking Then, the accuracy between the apparent surface charge σ from (2.21) and σ L from (3.60) for r ≥ 1/2 and t ≥ −1/2 satisfies In addition, the accuracy of the reaction potential P REAC from (2.18) and PROOF.The single layer operator V can be decomposed [29] as and L is a compact operator.Therefore, it is a Fredholm operator of zero index which is in addition injective [31].As a consequence, equations (2.21) and (3.60) are uniquely solvable.On account of the Strang lemma [9] with perturbed right hand side, the orthogonal projection results in On the 2D-level L for piecewise constant setup, one has In addition, apply the inverse estimate A combination of the above relations yields therefore On the other hand, for the estimation of f − f L L 2 (Γ) , one uses the same reasoning as above where K 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 {x k } and the molecular surface Γ, By using the orthogonal projection Q L and an inverse inequality, one obtains As a consequence, one concludes the following accuracy by applying (3.71),

Truncated multiscale PCM
In this section, we consider the compression scheme for the PCM.The matrix entries where the distance between the supports is beyond a level dependent threshold are annihilated.For a matrix M = M ij j=1,...,Q i=1,...,P , define For two levels ℓ, ℓ ′ = 0, 1, ..., L, a threshold δ ℓ,ℓ ′ whose expression will be specified subsequently is used.We consider the matrix coefficients They correspond to the matrix entries in (2.23) and (2.24) for the single layer kernel K = K SL and the double layer kernel K = K DL .The following level depended truncation procedure for each pair of levels (ℓ, ℓ ′ ) is applied: The block matrices A ℓ,ℓ ′ and A ℓ,ℓ ′ have respectively the entries [A (ℓ,i),(ℓ ′ ,i ′ ) ] i,i ′ and [ A (ℓ,i),(ℓ ′ ,i ′ ) ] i,i ′ .The uncompressed and compressed matrices are blockwise: For the single layer and the double layer, define respectively (4.83) For the involved kernels, one has on account of the Calderon-Zygmund inequality : LEMMA.Consider two levels ℓ, ℓ ′ = 0, ..., L. Consider any positive real parameters (s, p, q).Suppose where λ > 1 and ∆ ℓ denotes the diameter of the supports of any wavelet ψ (ℓ,i) on level ℓ.Then, the following error estimates in matrix norms are fulfilled: By using the thresholds (δ ℓ,ℓ ′ ) for all ℓ, ℓ ′ = 0, ..., L as in (4.85) such that the parameters (s, p, q) satisfy (4.89) s + p + q = θ 0 + 2 one has for u ∈ H γ (Γ) and v L ∈ V L (Γ) in the case the single layer compression while the double layer compression yields PROOF.The value of K p,q (u, v) from (4.80) will be examined on a pair of patches Γ p × Γ q .The indices (p, q) are dropped to simplify the notation.Fix ( u, v) = (ũ 1 , ũ2 ), (ṽ 1 , ṽ2 ) as the midpoints of S (ℓ,i) and S (ℓ ′ ,i ′ ) .By considering any (u, v) = (u 1 , u 2 ), (v 1 , v 2 ) ∈ S (ℓ,i) × S (ℓ ′ ,i ′ ) , the Taylor expansion leads to for some (u * , v * ) ∈ S (ℓ,i) × S (ℓ ′ ,i ′ ) .For the first term, by multiplication with a tensor product wavelet basis and by taking the integration over × , one obtains (4.93) which is zero due to (3.33).As for the summation over |α| = |β| = 1, supposing that the Jacobians of the mappings γ p and γ q are bounded: Concerning the estimation of the far-field where the inter-support distance is larger than the diameters of the supports: Hence, one has the next inequalities for x = γ p (u) and y = γ q (v) On account of the properties (3.48)-(3.51),one deduces (4.99) Use the expression of δ ℓ,ℓ ′ and note that θ 0 + 2 + 1/L > 0 to obtain and hence, (4.100) For any x ∈ Γ, we have

By using the measure property (3.48), obtain
The 1-norm is processed in an analogous manner and therefore Consider the orthogonal projection Q ℓ onto V ℓ (Γ) and introduce Blockwise according to the levels (4.104) By utilizing the next estimates one obtains based on the Cauchy-Schwarz inequality and (4.88), The use of the sum of geometric sequence leads to and hence For piecewise constant setup, obtain for the constant parameters s, p, q For the single layer V, set γ ′ = −1/2 and use s Similarly for the double layer, setting γ ′ = 0 results in THEOREM.Consider the truncated PCM equations on level L: where the operators V and K are obtained from V and K by using the threshold δ SL ℓ,ℓ ′ and δ DL ℓ,ℓ ′ respectively.Suppose the constant parameters (s, p, q) are chosen such that (4.115) Then, the numbers of nonzero matrix coefficients are reduced from In addition, the accuracy between σ in (2.21) and σ L in (4.113) for r ≥ 1/2 and t ≥ −1/2 verifies Further, the error between the reaction potential P REAC from (2.18) and We consider first the number of nonzero entries of the compressed matrices V and K.For computing N ( A ℓ,ℓ ′ ) for each fixed level pair (ℓ, ℓ ′ ), consider a basis function ψ (ℓ,i) on level ℓ.According to [30], the number of basis functions For that, Lemma 3.2. of [30] uses an argument using a sphere of radius ∆ ℓ + ∆ ℓ ′ + δ ℓ,ℓ ′ centered at a point located upon the support S(ℓ, i).There are 2 ℓ basis functions ψ (ℓ,i) on level ℓ for each patch Γ p where p = 1, ..., N. The number of the interlevel nonzero entries on all patches is thus (4.119) By summing over all levels The first sum on all patches verifies (4.120) For the second sum N 2 , the expression of the threshold δ ℓ,ℓ ′ leads to By using sum of geometric sequence, deduce In a similar manner, one obtains Hence, obtain for the second sum by using s Concerning the first equation (4.113), the next Strang relation [9] holds .
By inserting (4.111) in the last estimate, the next inequality is fulfilled In the same manner as (3.67), one has On the other hand, for the second equation (4.114) by applying the Strang lemma [9] for the second time By using (4.112), deduce Based on the last estimate and (4.123), deduce Concerning the accuracy with respect to the electrostatic reaction potential P REAC , one proceeds as in (3.72) and (3.73) in order to obtain

Practical implementation and numerical results
We want now to present in this section some results of the former method.We present separate results pertaining to the double layer potential and the single layer potential because if they are solved individually, analytical results for arbitrary bounday Γ exist for comparisons.The only well-known analytical result for the reaction potential P REAC is a simple geometry composed of a single atom.1. Quality of the resulting patches.

L-value
5.1.Molecular models.We employ different sorts of quantum models including molecules which are acquired from PDB 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 hydrogen and oxygen atoms contained in that large sphere constitute the components of the water clusters.The formerly proposed method assumes that the sizes of the patches are proportional.We made efforts to have geometric processing which results in patches of similar surface areas.In addition, efforts were done to obtain patches whose shapes related to the normal derivatives are as good as possible.The ideal quadrilateral which is a perfect square corresponds to L ideal = 0.25π ≃ 0.785398.In Tab.For the double layer potential, one can solve the system iteratively without recourse to any preconditioner because the system is well conditioned.A GMRes method serves as a solver of the system which is non-symmetric.We examine first the error reductions in term of the multiscale level for large molecules.
The results are collected in Tab. 2 which contains both the absolute errors and the relative errors.The models consist of a water cluster and a DNA section which admit respectively 1109 patches and 2119 patches.We consider two exact solutions which are respectively U 1 (x) = .2x 2 − .15y 2 − .05z 2 and U 2 (x) = exp(.5x) cos(.5 y) that have vanishing Laplacian.The right hand side g(x) is the restriction of the function U on the boundary Γ.The error reduction is affected by the exact solutions but in general the errors reduce satisfactorily in function of the wavelet levels.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. Tose rsults have been done for a water cluster molecule containing 3xx patches.The ration between the nonzero counts of the compressed and uncompressed matrixes are also exhibited.We examine now the general  In fact, they decrease linearly in logarithmic scale in function of the BEM levels.
5.3.Single layer potential.The single layer yields a system which is badly conditioned.We use an additive domain decomposition preconditioner to remedy that bad conditioning.In term of geometric structure, the overlapping domain decomposition is as follows Γ i where (5.127) Ω p ∩ Ω q is not necessarily empty for p = q. (5.128) In term of linear spaces, this leads to the decomposition (5.129) Denote the orthogonal projection onto V(Ω p ) = V•, • with respect to the bilinear form V(•, •) related to the single layer operator by (5.130) The ASM operator is defined by (5.131)The ASM solver consists in solving a local problem which is projected by P Ω p in each subdomain 'Ω p .A comprehensive description is detailed in our former work [42].The results of the number of iterations are summarized in Figure 2(a) where it is observed that the ASM method needs a lot less iterations than the direct method.A larger view for the iterations less than 50 is exhibited in Fig. 2(b).The same test fo the compression is detailed in the secind 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 Fig. 3(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 P REAC which are computed with the PCM-wavelet technique.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 R > 0 which is centered at the origin and which is assigned an electric charge q while the dielectric coefficient ε ≥ 1 of the solute can be adjusted.The analytical expression for the reaction potential [1] of the Born ion is P REAC = −(q 2 /8πR)(1 − 1/ε).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 Fig. 4(a) in which we consider three configurations corresponding to (R, q) = (1.2,2.0), (R, q) = (2.2,2.0) and (R, q) = (1.5, 1.0) respectively.We execute the wavelet-PCM simulation on multiscale level 3 while the dielectric coefficient varies from ε = 1 to ε = 200.We observe for all three configurations that the exact reaction potential and the computed PCM results align well in Fig. 4(a) where a zoom  for ε ∈ [1,20] in which the reaction potential drops down very quickly is also provided.As a next test, we investigate the variation of the electrostatic potential in function of the electric charge q.The corresponding outcome is depicted in Fig. 4(b) where we consider four configurations corresponding to R = 4, R = 3, R = 2, R = 1 while the dielectric coefficient is constantly ε = 200.The electric charge is allowed to vary in the range q ∈ [1. 5, 4.25].One observes in Fig. 4(b) that the computed reaction potential varies quadratically in function of the electric charge.The electrostatic reaction potential varies in a parabolic manner in function of the charge q.The next test for the Born ion is displayed in Fig 4(c) where one compares the computed PCM result with the exact solution when the ion radius varies from R = 0.5 to R = 3.5.One (R, q)=(1.2,2.0)(R, q)=(2.Table 3. Computed reaction potential in term of the multi-scale levels for the Borm ion. considers three configurations where the electric charge is respectively q = 1.5, q = 2.0 and q = 2.5 while the dielectric coefficient is consistently ε = 12.5.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 agree well with the expected theoretical expectations.As a last test for the Born ion, we examine the accuracy of the electrostatic potential in function of the BEM-levels.The levels are set to vary from 1 to 6.We consider two configurations where the radii and electric charges are (R, q) = (1.2,2.0) and (R, q) = (2.3,2.5).The results of the BEM accuracy are summarized in Table 3 where we consider four dielectric values ε = 15.0,ε = 5.0, ε = 100, ε = 45.0.The exact and the computed reaction potentials are shown there.It is observed that the BEM-approximation become satisfactorily accurate as the multiscale levels increase.
Most results in this article are computed for analytical comparison where we do not use any physical units.A detailed comparison with well established quantum softwares will be reported in a separate forthcoming article.We mention only that we have made a few comparisons with real chemical simulations.For that, we use the APBS (Adaptive Poisson-Boltzmann Software) which is a reliable American software [3,4,23] Physical constants values  widely used in Universities, Laboratories and Institutes across the USA and beyond the USA.The result of the comparison is depicted in Fig. 4(d) where the median exhibits the agreement between the method using the PCM-wavelet and the FDM (Finite Difference Method) for the reaction potential measured in kJ/mol.For that we use a propane molecule and water cluster assigned with random electric charges.We depict in Fig. 5 an instance of the apparent surface charge on a real molecular model acquired from water cluster.Along with the APBS packet, we use also the format PQR which can be deduced from PDB by using PDB2PQR [15,16].The electrostatic charges are measured in elementary charge while the other physical units are summarized in Table 4.

Conclusion
We consider the transmission problem occuring 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 structure in four-sided NURBS patches for which wavelets on the unit square are constructed.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.A comparison to a reliable software exhibits also a good agreement when evaluated in physical units.

Figure 1 .
Figure 1.Patch representation of a Water Cluster with 1089 NURBS.

2 w i
where ∂ represents the conormal derivative.Combine them with the transmission jump conditions [[w]] := w e − w i = f e − f i and [[∂w]] := ∂w e − ∂w i = ∂f e − ∂f i to obtain(2.13)

Figure 2 .
Figure 2. Practical implementation of preconditioning the single layer potential V. Number of iterations vs. residual error.

Figure 3 .
Figure 3. Relative errors for several RHS in function of increasing levels using a propane molecule: (a)Double layer, (b)Single layer.

Figure 4 .
Figure 4. Unitless comparison with Born ion: (a) dielectric coefficients, (b) electric charge, (c) radius.Comparison of real molecular surfaces with APBS in (d) using unit kJ/mol.

Figure 5 .
Figure 5. Apparent surface charge on a molecular surface.

Table 2 .
1, we gather the results of our test which consists in computing the average values of L over the whole patches.We obtain a satisfactory quality because the average values of L do not tend to zero.BEM accuracy for the double layer and the single layer in function of the levels.