A coupled BEM-FEM for molecular

Preparing to load PDF file. please wait...

0 of 0
A coupled BEM-FEM for molecular

Transcript Of A coupled BEM-FEM for molecular

Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X
A coupled BEM-FEM for molecular electrostatics with electrolyte effects M.D. Marcozzi
?, CTmW 5W&y MYz^ry
Abstract A coupled boundary element - finite element method is developed to compute the electrostatic potential inside and around molecules in an electrolyte solution. A nonlocal boundary value problem is developed by linearizing the far field pattern of the electrostatic potential, obtaining an integral representation of the solution exterior to an auxiliary boundary. We demonstrate existence and uniqueness of the weak solution as well as the strong convergence of the Galerkin approximation. Convergence estimates are provided for a coupled BEM-FEM approach. 1 Introduction Electrostatic effects play an important role in determining the intramolecular structure and stability of biological macromolecules, while electrostatic interactions affect intermolecular energetics and kinetics. To this end, we consider a macroscopic dielectric model in which the molecular electrostatics are governed by the Poisson equation with different charge distributions for the interior and exterior regions of the molecule. For the interior region, a volume distribution of point charges, originating from charged groups in amino acid residues (for proteins) is commonly used, while for the exterior region, the number density of ionic species, considered as point charges, in the surrounding electrolyte solution is assumed to follow a Boltzmann distribution as in the Debye-Huckel theory of electolytes (cf. [14] and references therein).
The governing field equations, the Poisson equation coupled to the nonlinear Poisson-Boltzmann equation, constitute a strongly nonlinear system in an unbounded domain. In order to expedite a numerical approximation of the solution we consider a linearization of the farfieldbased upon the decay of the electrostatic potential away from the molecule. By introducing an auxiliary boundary

Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X

128 Boundary Element Technology

into the exterior of the molecule, a representation of the exterior linearized problem may be obtained and coupled to the nonlinear one in the annular region bounded internally by the surface of the molecule. In a similar fashion, electrostatic effects interior to the molecule governed by the Poisson equation may be represented nonlocally on the molecule's boundary. The theory of monotone operators is utilized to treat the transmission problem on the bounded annular region and approximate solutions are sought through Galerkin procedures.
The problem of determining the electric potential and the consideration of electostatic effects is not a new one. Computations based upon the finite difference method [12] and the finite element method [11] have been made assuming far-field conditions imposed along a computational auxiliary boundary. Boundary element methods have been applied to the linearized problem as well [14]. These methods, however, are formal in their execution. We propose then a coupled boundary element - finite element scheme, demonstrating the unique existence of a generalized solution and the strong convergence of the method.
The outline of the paper is as follows. Section 2 is devoted to the formal development of the transmission problem in the bounded annular region. The development of the generalized problem and it's representation as an operator equation is considered in section 3. The following section 4 is composed of technical lemmas and the main result; namely, the strong convergence of the Galerkin approximates to the unique generalized solution. The final section 5 concerns the formulation and convergence properties of the proposed coupled BEM-FEM.
2 Problem statement
Let (7 C 3ft3 be a bounded annular region with smooth interior and exterior boundaries TQ and F, respectively. We label HO the region enclosed by TQ while the region exterior to F will be denoted by Og. Formally, we wish to determine the triple (UQ,U,U<>) in appropriately defined function spaces such that

= go6(x-xo); Vx€ %, fo


Au (x) = K* sinh u (x); V x e H,


Aw,(x) = *&, (x)i VxEH,,


4- %=) %,(x) = o (|x|"i) as |x| -+ oo,

subject to the so-called transmission conditions across the frontiers of 17; namely,

UQ = u* and eo<9u0 /<9n = tdu* /dn on FQ



if" = %+ and &r/9n = g^/^n on F.


Here, XQ denotes a point in HO, n the unit normal vector pointing away from the respective bounded enclosed region, and u^ the limit (trace) of any function u on

Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X
Boundary Element Technology 129
<90 counter to/in the direction of n. In particular, the Cauchy data are not independent, but will be seen related through (2) by a Dirichlet-Neumann mapping. Physically, the transmission conditions are statements of solution regularity.
By the Green's formula, we may obtain a representation of the solution to the interior problem (la) in terms of the Cauchy data on TO such that

- <-> -?*<*- *»> + L * <- o f£ (3)
lor all x G HO, where EQ (x-f) = l/47r|x-(| is the fundamental solution to the three-dimensional Laplacian. Letting x tend to a boundary point nontangentially, it follows as a consequence of the jump relations (cf. [3]) that
«JT (x) = g& (x - xo) + (V^ (x) + (I/ - V.) («-) (x) , (4
for all x G FQ. Here, / is the identity while \\ and V£ are the single and double layer potential operators given by
(x) = £„ (x - 0 ao (0 d (x) =
respectively, for any x G FQ. Considering similarly the normal derivative of (3) as x — » To from within OQ, we obtain likewise,

- Ss* (- -> + (' + 4 for x G TO, where

w - K (-.I w . <*)

are the adjoint double layer and hyper-singular potential operators [8], respectively. The equation (4a) (or 4b) suffices to represent the solution to the interior problem (la) and so, when coupled with (Ib) through the transmission conditions (2a), constitute a nonlocal boundary condition for w, relating the Cauchy data on FQ.
A similar analysis holds for the exterior problem (Ic) for which we find




Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X

130 Boundary Element Technology

for all x G F. With respect to the exterior problem and the representations (5), we have that E (x — £) = e~*'*~^/47r|x — £|, which is the fundamental solution for the Yokawa (cf. Helmholtz with k = IK, K > 0) equation in $?. Utilizing the expressions (4) and (5) and the transmission conditions (2), the original problem (1) may be reduced to the nonlocal transmission problem

Au (x) = t? sinh u (x) ; V x

/ 4- It %+ (x)-- % €Q-- (*)== -& (x - %o); V x E To, / - V, tT (x) + Vi - (x) = 0; V x 6 T,
where the Neumann data have the representation

(6b) (6c)

60 c/n and

= CQ (7n * - \z ' / (\^fo tc/rn y

(x) = % («-) (x) + / - v;

W = 0; Vxe r. (7b)

We note that as opposed to the formulation (1), the nonlocal problem (6) is posed on the domain fi of finite extent. This circumstance avoids the complications inherent for problems posed on 3P (cf. [15] and references therein).
3 Weak formulation
The weak formulation for the solution of problem (6) can be determined by the usual variational approach. To this end, we introduce the subspace D (0) = {v G H* (H) | sinh v (x) • v (x) € LI (H)} and define the generalized problem associated with (6): Let O-Q = du+/dn on TO and a = du~ /dn on F. We seek U = («,<7o,<7) eXD = D(tyxH~* (Fo) x H~* (F) such that

K= (sinh tz, 1;)^ + +, / + ^ ^o --%+, %^^


/ /To fO

' ^ o

= -(«+, f !;£o (x - xo)\ ,


/ 1o


--^(£o(x-xo),xo}r», (8b)




/To CO

(^^%)r - ((^ - ^2) ^,%)^ = 0,


for all test functions V = (v,Xo,x) € y = //^ (fi) x H'* (To) x #"& (F), where

a (u, v) = Vu • Vv dx ; (w, u)^ = u • v dx,


«/ 17

Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X

Boundary Element Technology 131
and (-, -} denotes the respective duality pairings on #»/' x H~W with respect to To and T. Formally, we obtain (8a) by multiplication of (6a) with v € C~ (ft) integrating by parts, and substituting in the relations (7). The construction of the set D(Sl) here ensures the existence of (sinh«,t,)^ and (sinh«,wL for all u e D (ft) and v e H * (ft). We remark that the space of test functions over ft is chosen to be more regular than that of the solution space. This is done to compensate for the restriction of the solution space from H* (ft) to D (ft). As such, we employ a Galerkin method in the space of greater regularity H* (SI) which will converge in the weaker space H* (ft).
We consider then the reduction of the variational problem (8) to an operator equation. To this end, we first note that the bilinear form a is bounded by the Holder mequality, in which case there exists a unique linear continuous operator A : H (ft) -> H (ft) such that (Au,v)^^ = a(u,u), for all u £ H* (SI) and v € H*(£l). With respect next to u € D(Sl), we consider the boundedness of (smh«,«)o for any v € H' (ft). By the Sobolev (Kondrasov) embedding
theorems, it follows that IP (ft) A C (ft). Moreover, |sinh%| = |u-i.sinhu-u| < sinhu-u, for |u| > 1. It follows then for u € D (ft) that x .— » smhu(x) € L, (SI) From this, we find

(sinhu,v)n < IHIc(n) ^ |sinhw|dx < con^. (Ml^,,, ,


for all v € H^(Sl). There exists then a uniquely determined operator B • D (ft) C

H (U) -+ ^-' (ft) such that (Bu,v)^) = (sinh u, «)„ , for all u € D(ft) and

w € H (ft). Finally, the regularity of the layer potentials and their derivatives

iollows from the properties of pseudodifferential operators and the order of their

symbols [4].Employing the mapping properties of the trace, we conclude that

there exists a unique linear continuous operator C • X = H* (SI) x #-1/2 fr 1 v

H-W (T) -,y = H- (ft) x #1/2 (r,) xffV2(T) such that


1 0 for all V G 3^.
Similar arguments may be applied to determine a unique operator 7" E AT* that

Defining A : XD C X -> y by .4 = .4_+_2 + C,


Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X

132 Boundary Element Technology

the generalized problem (8) corresponds to the operator equation

AU = F,


for all % E Ab. In particular; #% (H) C C H C D (0) C #* (H).We consider then the operator equation (11) along with the Galerkin method

(^/n - f, W&) = 0, 6 = 1, 2, . . . , n


where % G X = span {Wi, Wg, . . . ,WJ C ^.
4 Main result
In order to demonstrate our primary result, we need to first verify that the operator A, specified by (10), satisfies certain technical conditions. We provide for this by showing that A is strongly monotonic and globally coercive, satisfies a generalized monotonicity condition, and is quasi-bounded. Our main result, theorem 3 follows.
Clearly, the forms a (v, v) and (sinh v, v)^ are non-negative for all v G H* (0) . Moreover, there exist (generic) constants c such that

)r > c||x||/2^ (13a) and

- („+, TW)^ > c ||v||^(n, ; - (tT, V&~\ > c \\v\\^ ,


for all V &y. The coerciveness of Vi was demonstrated in [7], while the identities

VVo (x) = (nx x Vx) • Vi • (ny x Vy^o) (x) ; for x € TO, (14a)


x) = (nx x Vx) • li • (ny x Vy/u) (x)


follow from [8], for all ^ € H*l* (To) and p € H*'* (T). The results (13) are then a consequence of (14), integration by parts on 0H, the coerciveness of V\ and the trace theorem. We have therefore shown the global coercivity of A\ i.e. there exists a constant c > 0 such that



for all V € y. Moreover, since sinhu is strongly monotone for u G ^, it follows that there exists a number d > 0 for which



Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X

Boundary Element Technology 133

for all U - V ^ 0 € XD, i.e. A is strongly monotonic on XD* The following result is a technical generalization of the monotonicity condi-
Lemma 1 (Generalized monotonicity condition). Let f G X* and (%) be a sequence m ^ W^ % -" &/ m ^ as 7% -» oo aW



Proof. We follow here the work of Zeidler (cf. [15]). The operators A and
C from X into y* are linear and continuous and therefore weakly sequentially continuous, i.e.



as n -4 oo, for ail V G ];. Noting that ^^ ^) ^ (sinhw,u)^, along with the estimate (9),it suffices to show that

sinh (un (x)) -> sinh (u (x)) in LI (0)


as n -> oo. Moreover, note that it is sufficient to prove (19) for a subsequence of (%n).
We begin by showing that sinh (%(•))%(•) e ^i(H), i.e. % E D(H). In particular, the embedding #* (H) -. 63 (H) is compact. This implies %* -» % in LI (ft) as n -^ oo. Thus, there exists a subsequence, again denoted by (%*), such that ^ (x) -> u (x) as n -> oo, for almost all x G H. Moreover, from % — &/ in A', it follows that sup* IJ^II^ < oo, in which case (4%*,%.)^^ < c|K||^ -
^ and (C%,%)y < c||%||^ < con6(. Relation (17b) provides that


%n • W*dx < C,


where it follows from the continuity of the hyperbolic sine that sinh^(z) • ^(a:) -» sinh i/(a;) • u(z) as n -^ oo, for almost all x E ft. Therefore, by the lemma of Fatou, we have sinh u (•) • u (•) E LI (ft).
We now prove (19). Let d > 0 be Axed. For each x E ft, we have either K (x)| < d or |sinh%*(x)| < cf'i sinh^ (x) • w* (x). Note that sinhu = i/-^ . sinh w • % if %y 0. Moreover, we get |sinh%| < c((f) for all % such that |w| < &
Let now H be a measurable subset offt,in which case

Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X

134 Boundary Element Technology

By (20), the constant ci is independent of n, d, and H. Hence, J~ sinhundx < |, for all n, if d is sufficiently large and 0 is sufficiently small. By the absolute continuity of the integral, it follows that for each e > 0, there exists a 8 (e) > 0, such that

/ |sinh%n — sinhu| dx < / (|sinhun| -f |sinhw|) dx < e,



for all n and for all subsets 0 of £1 with 0 < 8 (e). By the Vitali convergence

theorem, (19) follows.


Lemma 2 (Quasi-boundedness of A.} Let (Un) be a sequence in the space y with Un —* U in X, as n —•» oo, and suppose that

for all n, then the sequence (AUn) isbounded iny .

Proof. The weak convergence of %to U in X implies that (Un) is bounded in X

and so lim^oo (>AUn,Mn}y < const. In order to obtain a contradiction, suppose

that (AUn)is unbounded in y*. Then there exists a subsequence, again denoted

by (Un) such that

||"4%n||y. -» oo as n -» oo.


By the same argument as in lemma 1, we obtain that (AUn, V)y — + (AU, V)y as

n — > oo, for all V € y, by passing to a subsequence, if necessary. The uniform

boundedness principle then implies that the sequence (AUn) is bounded. This

contradicts (21).


We now establish our main result on the unique existence of the electrostatic


Theorem 3 The transmission problem (6) has a unique, bounded generalized solution u G D (H). Moreover, for each n £ N, there exists a unique solution Un such that Un —•> u in H* (0) as n — > oo. Proof. Existence of U and Un follows from the lemmas 1 and 2 and theorem 27. B in Zeidler (cf. [15]). Uniqueness is an immediate consequence of global coercivity.
Strong convergence of the Galerkin approximation is shown as follows: Given % -+ % as n -» oo, then X%* -»./%/ in y and (>%4,%)y n— » oo. Therefore, by the uniform monotonicity of A (cf. eqn 16),

as n —> oo, in which case Un — » U.

Finally, we demonstrate the boundedness of U by noting that since A is

globally coercive (cf. eqn 15), there exists an R > 0 for which (AU — F,U)y > 0,

for all U such that \\U\\x > R. However, if U is a solution of (6), then MA = JF,

in which case ||ZY||, < R.


Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X

Boundary Element Technology 135

5 Convergence estimates
Towards quantifying the quality of the Galerkin approximation, we consider the following abstract error estimate. Lemma 4 TVtere ezz'afs o comgfomf C, 2Wepe%de%2 o/A, awcA ^Aa(


Proof. Let V* be an arbitrary element of X- It follows from (11) and (12) that - ,424, V,Jy = 0, in which case, by (16),

where the latter is a consequence of the boundedness of the operators A and C For bounded % and % (cf. theorem 3), there exists a constant C > 0 such that

by the monotonicity of the hyperbolic sine function.


We will consider the convergence of the Galerkin method in the case of con-

formal finite elements. To this end, we assume that the solution % of (11) is in

the more regular space ^ (H) x JT (Fo) x #' (F) ( and note that 6 E D') and

suppose n to be polygonal. For the analysis of nonconforming methods, we note

[6], p], [10], and [9], to name a few. Constructing then an afEne family of regular

tetrahedron of type 1 (and class C^) over W and piecewise constant splines over

9ft, we obtain the following convergence estimate.

Theorem 5 6ef% E D(ft)n#^) % #^0) x ^F 6e

Proof. As a consequence of (22), we have that

where n.% denotes the interpolant of the function %. Standard interpolation

error estimates then imply the result (cf. [2], [1], [13]).


1. Babuska, I., & Aziz, A.K. 7%e MofAemofW Fo%%
Transactions on Modelling and Simulation vol 9, © 1995 WIT Press, www.witpress.com, ISSN 1743-355X
136 Boundary Element Technology
2. Ciarlet, P.B., The Finite Element Method for Elliptic Problems, NorthHolland, Amsterdam, 1978.
3. Colton, D. & Kress, R. Integral Equation Methods in Scattering theory, Wiley Interscience, New York, 1983.
4. Costabel, M. Boundary integral equations on Lipschitz domains: Elementary results, S/AM J. MofA /IW., J&&9, 19, 613-626.
5. Gatica, G.N. & Hsiao, G.C. On the coupled BEM and FEM for a nonlinear exterior Dirichlet problem in %\ Numer. Math, 1992,61 171-214.
6. Hsiao, G.C. & Marcozzi, M.D. Variational methods for the potential flow past an airfoil, submitted to SI AM J. Math Anal.
7. Hsiao, G.C. & Wendlarid, W.L. A finite element method for some integral equations of thefirstkind, J. Math. Anal. Appli.,1977, 58, 449-481.
8. Hsiao, G.C. & Wendland, W.L. Variational Methods for Boundary Integral Equations and Mathematical Foundations of Boundary Element Methods, in preparation.
9. Johnson, C. & Nedelec, J.C. On the coupling of boundary integral and finite element methods, Math. Comp., 1980, 35, 1063-1079.
10. LeRoux, M.N. Method d'elements finis pour la resolution numerique de problemes exterieurs en dimension 2, R.AJ.R.O. Anal. Numer., 1977, 11 27-60.
11. Orttung, W.H. Direct solution of the Poisson equation for biomolecules of arbitrary shape, polarization density, and charge distribution, Ann. N.Y. /lead. Sc2'.,1977, 303, 22-37.
12. War wicker, J. & Watson, H.C. Calculation of the electic potential in the active site cleft due to a-helix dipoles, J of Mol. BioL, 1982, 157, 671-679.
13. Wendland, W.L. On asymptotic error estimates for the combined BEM and FEM,Finite Element and Boundary Element Techniques from Mathematical and Engineering point of view, ed E. Stein, W.L. Wendland, CISM Lecture Notes 301, pp 273-333, Udine, Springer-Verlag, Wein-New-York, 1988.
14. Yoon, B.J. & LenofF, A.M. A boundary element method for molecular electrostatics with electrolyte effects, J. Comp. Chem., 1990, 11, 1080-1086.
15. Zeidler, E. Nonlinear Functional Analysis and its Applications II, SpringerVerlag, New York, 1990.