Domain Decomposition Method for Implicit Solvation Models

Get Complete Project Material File(s) Now! »

Domain decomposition methods for implicit sol-vation models

Implicit solvation models are widely-used in the chemistry community, but little interaction with applied mathematics can be observed despite the fact that solutions to partial differential equations need to be approximated. The underlying physical law is provided by electrostatic interaction involving elliptic operators.
In this section, we first introduce two state-of-art Schwarz domain decomposition (dd) methods respectively for solving the COSMO and the PCM: the ddCOSMO and the ddPCM. Then, we introduce two new Schwarz domain decomposition methods respectively for solving an SES-based PCM and the PB solvation model, which are contributions of this thesis. The methods generally consist of two steps:
[1] if the electrostatic problem is defined on an unbounded domain, transform the original problem equivalently to problems that are defined on a bounded domain and might be coupled by a non-local condition;
[2] develop a classical Schwarz domain decomposition method that decomposes the bounded domain into balls and only solve local coupled sub-problems restricted to balls.
For the SES-based PCM, we will first construct the model and then propose the corresponding domain decomposition method for solving it, called the ddPCM-SES method. Further, for the PB solvation model which is already established in Section 2.2.3, we will propose a domain decomposition method for solving it, which is called the ddLPB method.

ddCOSMO

To solve the electrostatic problem of the COSMO, the finite element method or the finite difference method can be used. However, the computational cost is too expensive for a large realistic molecule. In particular, meshing the solute cavity Ω of a complicated molecule is already too costly. The state-of-art COSMO solver is the ddCOSMO [25, 72, 69], a Schwarz domain decomposition method [96, 107] developed for solving the COSMO in the past several years, which has attracted much attention due to its impressive efficiency, that is, it performs about two orders of magnitude faster than other equivalent methods [69].
The reaction potential is indeed the electrostatic potential that is additionally created by the pressure of the solvent with respect to vacuum.
Using the Schwarz decomposition of Ω, this Laplace equation can be recast as the following group of coupled sub-equations, each restricted to Ωj:
Figure 6: 2D schematic diagram of Γij (red) and Γej (blue) associated with Ωj.
then applied to solve the coupled sub-equations. For instance, the idea of the Ja-cobi algorithm is to solve each local problem based on the boundary condition of the neighboring solutions derived from the previous iteration. During this iterative procedure, the computed value of ψr|Γij is updated step by step and finally converges to the exact value.

ddPCM

Latter on, a similar Schwarz domain decomposition method for solving the clas-sical PCM with the PDE (1.5) of electrostatic potential and the definition (2.2) of ε(x) was proposed in [108] (called ddPCM, see also website [70]), which is based on the integral equation formalism (IEF) of PCM [26].
In the above two equations, the second one is equivalent to the Laplace equation of the COSMO and therefore, can be solved by the ddCOSMO solver. Solving the first equation involves decomposing the solute cavity Ω into balls and solving a group of sub-equation each restricted to a ball. To do this, a similar Schwarz decomposition method to the ddCOSMO is applied, see details in [108].
Both the ddCOSMO and the ddPCM fully take advantage of the geometrical structure of the VdW-cavity or SAS-cavity, i.e., the union of balls. These two methods give a fast resolution of respestively the COSMO and the classical PCM. Then, we would like to generalize them to an SES-based PCM and the PB solvation model.

ddPCM-SES

As mentioned, the previous two methods (ddCOSMO and ddPCM) are based on the VdW-cavity or the SAS-cavity, due to their simpler geometrical structure and easier computation. However, this might not be physically appropriate and the choice of the cavity is indeed important as pointed out in [110, Section II. C.]: The shape and size of the cavity are critical factors in the elaboration of a method. An ideal cavity should reproduce the shape of the solute M, with the inclusion of the whole charge distribution ρM and with the exclusion of empty spaces which can be filled by the solvent continuous distribution. If the cavity is too large the solvation effects are damped; if it is too small serious errors may arise in the evaluation of the interaction energy for the portions of ρM (atoms or bonds) near the solute-solvent boundary.
The SES-cavity, the region enclosed by the SES, is thought to be a more appropri-ate choice of solute cavity, since it has a stronger physical meaning in the sense that it represents the region where solvent molecules (represented by idealized spheres) can not touch. As a consequence, we would like to study the possibility of constructing and then solving an SES-based solvation model, such as an SES-based PCM.

SES-based PCM

In some chemical calculations, it has been confirmed that taking the SES-cavity Ωses into account can yield more accurate results, such as in [98, 91]. The shape of the cavity is represented by the dielectric permittivity function which in this case equals to one (the dielectric permittivity of vacuum) within the SES-cavity.
In the classical PCM, the permittivity function ε(x) in the form of (2.2) is dis-continuous and equals to the bulk solvent dielectric permittivity outside the solute cavity. A solver for such an SES-based PCM has been proposed in [50, 17], using the integral equation formulation of PCM and an efficient mesh generator of the SES. However, it has been argued that treating the solvent dielectric permittivity as con-stant is not sufficient and as a remedy, continuous permittivity functions εs(x) have been proposed [47, 13], based on the VdW-cavity or the SAS-cavity. But a method containing a continuous permittivity function based on the SES-cavity does not ex-ist, to our knowledge. We are convinced that introducing the SES-based continuous permittivity function is the next logical step to further refine the PCM.
In our model, we assume that the dielectric permittivity in the region Ω∞ far away from the solute molecule is equal to the bulk solvent dielectric permittivity εs, which is a reasonable assumption since at the position far from the solute molecule, the solvent molecules are influenced little by the solute molecule. Between the solute cavity and Ω∞, an intermediate dielectric boundary layer (the switching region) L := Ωc∞ ∩ Ωcses is constucted, to obtain a continuous dielectric permittivity function ε(x) of the following form (see the left of Figure 7)
See Figure 8 for an example of ε(x), where ε(x) is a distance-dependent function. The “distance” here represents the signed distance to the SAS, denoted by fsas, which also characterizes the SES-cavity as follows where the bounded domain Ω0 and the unbounded domain Ω∞ are complementary in R3, [u] and [∂nu] denote the jump of u respectively its normal derivative on Γ0 := ∂Ω0.
The bounded domain Ω0 is taken as a union of balls, inspired by the geometrical structure of the solute molecule. Since Ω0 consists of a union of balls, we propose to further apply a classical domain decomposition algorithm in order to solve the two problems (3.14) by only solving local sub-problems restricted to balls.
As a consequence, a Laplace solver and a Generalized Poisson (GP) solver are developed respectively for solving the Laplace equation and the GP equation in each ball, which allows us to use the spherical harmonics as basis functions in the angular direction to propose an efficient spectral method within each ball. It is important to notice that this algorithm does not require any meshing and only involves problems in balls that are coupled to each other within the domain decomposition paradigm.
In Chapter 3 of Part II, we will present the details of this Schwarz domain decom-position method for the SES-based PCM. We will first remind different solute-solvent boundaries including the VdW surface, the SAS and the SES. We will also provide more details on the construction of the continuous dielectric permittivity function ε(x) of PCM, ensuring that the SES-cavity always has the dielectric constant of vac-uum as explained above. Then, we will present the problem formulation of the PCM as well as its equivalent transformation, and therefore, a global strategy for solving the problem. Later, we will introduce the scheme of the domain decomposition method for solving the associated partial differential equations iteratively. This requires to develop the Laplace solver and the GP-solver in a ball, which are presented. After that, we will give a series of numerical results of the proposed method.

ddLPB

For the sake of simplicity, we focus on solving the linearized Poisson-Boltzmann equation (2.8) defined in R3. This problem can be divided into a Poisson equation defined in the bounded solute cavity Ω and a homogeneous screened Poisson (HSP) equation defined in the unbounded solvent region Ωc as follows Here, κ is the Debye-Hückel screening constant. We remind that the solute cavity Ω in the context of the PB equation consists of a union of balls (i.e., we take the VdW-cavity or the SAS-cavity due to the simple geometry).
Further, we use the integral equation formulation to represent the electrostatic po-tential ψ|Ωc in the solvent region, which simultaneously gives the potential ψe of an extended screened Poisson equation in the solute cavity (this is an interior Dirichlet problem) where ψe satisfies the same HSP equation as ψ|Ωc with the same Dirichlet boundary condition on Γ, but defined in Ω.
At this moment, the initial problem has been transformed equivalently into two cou-pled problems (3.18)–(3.19) both defined on the bounded domain Ω with a coupling condition (3.20). Taking advantage of the fact that Ω is a union of overlapping balls, the Schwarz domain decomposition can be applied to solve these two equations by respectively solving a group of coupled sub-problems each defined on a ball.
Ultimately, a Laplace solver and a HSP solver are developed respectively for solv-ing the Laplace equation and the HSP equation in each ball, which allows us to use the spherical harmonics as basis functions in the spherical direction to propose an efficient spectral method within each ball. Similar to the ddCOSMO, the ddPCM and the ddPCM-SES, the ddLPB does not require any meshing and only involves problems in balls that are coupled to each other.
Remark 3.1. The ddPCM-SES method and the ddLPB method are inspired by the previous ddCOSMO method and ddPCM method which run impressively fast. The ddLPB has the same computational complexity as the ddPCM.
In Chapter 4 of Part II, we will present the details about this Schwarz domain decomposition method for the PB solvation model. We will first introduce different definitions of the solute-solvent boundary as in Chapter 3 and also the Poisson-Boltzmann equation defined in the whole space in the case where the solvent is an ionic solution. Then, we will transform equivalently the original Poisson-Boltzmann defined in R3 to two coupled equations both defined on the bounded solute cavity. Next, we will present the global strategy for solving these two coupled equations, using the domain decomposition method. The domain decomposition scheme requires to develop two single-domain solvers respectively for the Laplace equation and the HSP equation defined in a ball. After that, we will give a reformulation of the coupling conditions that should be discretized and consequently derive a global linear system to be solved. We will also present a series of numerical tests about the calculation of the electrostatic contribution to the solvation energy using the ddLPB method.
This chapter has been published as a journal paper [93]. As mentioned in the introduction of this thesis, we present in this chapter a complete characterization of the Solvent Excluded Surface (SES) for molecular systems, including a complete characterization of singularities of the surface. The theory is based on an implicit representation of the SES, which, in turn, is based on the signed distance function to the Solvent Accessible Surface (SAS). All proofs are constructive so that the theory allows for efficient algorithms in order to compute the area of the SES and the volume of the SES-cavity, or to visualize the surface. Further, we propose to refine the notion of SAS and SES in order to take inner holes in a solute molecule into account or not.

Introduction

As mentioned in Section 1 of the introduction of this thesis, the majority of chemically relevant reactions take place in the liquid phase and the effect of the environment (solvent) is important and should be considered in various chemical computations. The implicit solvation model (or continuum solvation model) is a model in which the effect of the solvent molecules on the solute are described by a continuous model [109]. In continuum solvation continuum models, the notion of molecular cavity and molecular surface is a fundamental part of the model. The molecular cavity occupies the space of the solute molecule where a solvent molecule cannot touch and the molecular surface, the boundary of the corresponding cavity, builds the interface between the solute and the solvent.
A precise understanding of the nature of the surface is essential for the implicit sol-vation model and as a consequence for running numerical computations. The van der Waals (VdW) surface, the Solvent Accessible Surface (SAS) and the Solvent Excluded Surface (SES) are well-established concepts. The VdW surface is more generally used in chemical calculations, such as in the recent developments [25, 72] for example, of numerical approximations to the COnductor-like Screening MOdel (COSMO) due to the simplicity of the cavity. Since the VdW surface is the topological boundary of the union of spheres, the geometric features are therefore easier to understand. However, the SES, which is considered to be a more precise description of the cavity, is more complicated and its analytical characterization remains unsatisfying despite a large number of contributions in literature.

Previous Work

In quantum chemistry, atoms of a molecule can be represented by VdW-balls with VdW-radii obtained from experiments [97]. The VdW surface of a solute molecule is consequently defined as the topological boundary of the union of all VdW-balls. For a given solute molecule, its SAS and the corresponding SES were first introduced by Lee & Richards in the 1970s [67, 99], where the solvent molecules surrounding a solute molecule are reduced to spherical probes [109]. The SES is also called “the smooth molecular surface” or “the Connolly surface”, due to Connolly’s fundamental work [30]. He has divided the SES into three types of patches: convex spherical patches, saddle-shaped toroidal patches and concave spherical triangles. But the self-intersection among different patches in this division often causes singularities de-spite that the whole SES is smooth almost everywhere. This singularity problem has led to difficulty in many associated researches on the SES, for example, failure of SES meshing algorithms and imprecise calculation of molecular areas or volumes, or has been circumvented by approximate techniques [28]. In 1996, Sanner treated some special singularity cases in his MSMS (Michel Sanner’s Molecular Surface) soft-ware for meshing molecular surfaces [102]. However, to our knowledge, the complete characterization of the singularities of the SES remains unsolved.

Contribution

In this chapter, we will characterize the above molecular surfaces with implicit functions, as well as provide explicit formulas to compute analytically the area of molecular surfaces and the volume of molecular cavities. We first propose a method to compute the signed distance function to the SAS, based on three equivalence statements which also induce a new partition of R3. As a consequence, a computable implicit function of the corresponding SES is given from the relatively simple relation-ship between the SES and the SAS. Furthermore, we will redefine different types of SES patches mathematically so that the singularities will be characterized explicitly. Besides, by applying the Gauss-Bonnet theorem [35] and the Gauss-Green theorem [38], we succeed to calculate analytically all the molecular areas and volumes, in par-ticular for the SES. These quantities are thought to be useful in protein modeling, such as describing the hydration effects [100, 19].
In addition, we will refine the notion of SAS and SES by considering the possible inner holes in the solute molecule yielding the notions of the complete SAS (cSAS) and the corresponding complete SES (cSES). To distinguish them, we call respectively the previous SAS and the previous SES as the exterior SAS (eSAS) and the exterior SES (eSES). A method with binary tree to construct all these new molecular surfaces will also be proposed in this chapter in order to provide a computationally efficient method.

Outline

We first introduce the concepts of implicit surfaces in the second section and the implicit functions of molecular surfaces are given in the third section. In the fourth section, we present two more precise definitions about the SAS, either by taking the inner holes of the solute molecule considered into account or not. Then, based on three equivalence statements that are developed, we propose a computable method to calculate the signed distance function from any point to the SAS analytically. In this process, a new Voronoi-type diagram for the SAS-cavity is given which allows us to calculate analytically the area of the SAS and the volume inside the SAS. In the fifth section, a computable implicit function of the SES is deduced directly from the signed distance function to the SAS and according to the new Voronoi-type diagram, all SES-singularities are characterized. Still within this section, the formulas of calculating the area of the SES and the volume inside the SES will be provided. In the sixth section, we explain how to construct the SAS (cSAS and eSAS) and the SES (cSES and eSES) for a given solute molecule considering the possible inner holes. Numerical results are illustrated in the seventh section and finally, we provide some conclusions of this chapter in the last section.

Introduction to implicit surfaces

The above definition of an implicit object is broad enough to include a large family of subsets of the space. In this chapter, we consider the simple case where U = R3, V = {0} and f : R3 → R is a real-valued function. As a consequence, an implicit object is represented as a zero-level set O = f−1(0), which is also called an implicit surf ace in R3, and the function f is called an implicit f unction of the implicit surface. Notice that there are various implicit functions to represent one surface in the form of a zero-level set.

Table of contents :

Introduction 
1 The larger context
2 Implicit solvation models
3 Domain decomposition methods for implicit solvation models
I Molecular Surfaces 
1 Mathematical analysis and calculation of molecular surfaces 
1.1 Introduction
1.2 Introduction to implicit surfaces
1.3 Implicit molecular surfaces
1.4 Solvent accessible surface
1.5 Solvent excluded surface
1.6 Construction of molecular surfaces
1.7 Numerical results
1.8 Conclusion
2 Meshing molecular surfaces based on analytical implicit representa-tion 
2.1 Introduction
2.2 Molecular surfaces
2.3 Construction of molecular surfaces
2.4 Molecular inner holes
2.5 Meshing
2.6 Conclusion
II Domain Decomposition Method for Implicit Solvation Models 95
3 Domain Decomposition Method for the Polarizable Continuum Model based on the Solvent Excluded Surface 
3.1 Introduction
3.2 Solute-solvent boundary
3.3 Dielectric permittivity function
3.4 Problem formulation and global strategy
3.5 Domain decomposition strategy
3.6 Single-domain solvers
3.7 Numerical results
3.8 Conclusion
4 Domain Decomposition Method for the Poisson-Boltzmann Solva-tion Model
4.1 Introduction
4.2 PB solvation model
4.3 Problem transformation
4.4 Strategy
4.5 Single-domain solvers
4.6 Global linear system
4.7 Numerical results
4.8 Conclusion

GET THE COMPLETE PROJECT