NONLINEAR POISSON-BOLTZMANN: modeling, adaptivity, multilevel.

OVERVIEW ABOUT IONIC PROBLEM:

We consider the interaction of solute and solvent media which are respectively denoted by Ω^u and Ω^v. The surface Γ represents the solute-solvent boundary which is in our case the Connolly surface. The solute is located inside the cavity composed of immobile charges while the external domain is constituted of mobile ions. Our ongoing and future objectives:

Modeling of data acquired from PDB (Protein Data Bank),

Coarse curved tetrahedral patches from the Connolly surface,

Nonlinear correction with fast convergence,

FFT and FHT(Fast Hartley Transform) based nonlinear solver,

A-posteriori error estimator and adaptive refinement,

Higher order elements using Dubiner bases.

We consider the following nonlinear PDE admitting the unknown u, related to the electrostatic potential, on the ionic solution Ω=Ω^uUΩ ^v:

The coefficients ε(x) and κ(x) are space-dependent functions related to the dielectric value and the modified Debye-Hueckle parameter. Those coefficients might be highly discontinuous between Ω^u and Ω^v but the jump of solution u is required to be zero across the interface Γ. The value of n_i is the number density of counter-ions of type i and q_i is its charge while β:=1/(k_BT) with k_B and T being the Boltzmann constant and the temperature.

For real chemical simulation, the right-hand expression f(x) is related to the internal electric charges while the boundary function g is zero. But for our numerical models, we let them general for comparison reason. In the monovalent case, the above nonlinear expression is replaced by sinh. If the solvent ions are assumed to occupy sites on a lattice, the nonlinear expression is

In the monovalent case, by considering only the first term of the Taylor expansion, one obtains the linear (or linearized) Poisson-Boltzmann equation:

Such an approximation of sinh yields very inaccurate results for highly charged electrostatic potentials. As a consequence, we consider the nonlinear Poisson-Boltzmann equation (PBE).

MODELING AND DATA ACQUISITION FROM PDB:

For our simulation, we use two kinds of geometric discretization:

Unstructured meshes composed of fine tetrahedra. In this case, the simulation is performed at the physical domains.
Hierarchical representation: we decompose the domains into curved parametric tetrahedral patches and we simulate at the parametric domains.

Both methods require some geometric preparation. The geometric information are acquired from the Van der Waals radii and the PDB files containing the coordinates of the nuclei. We need a weighted Voronoi decomposition. For two spheres B₁(m₁,r₁) and B₂(m₂,r₂), we use the power distance

With respect to the power distance, one can tessellate the whole space into weighted Voronoi cells (Laguerre). This decomposition is obtained by using the convex hull in R⁴ of the initial points added with weights. The weighted Delaunay is used to derive the decomposition where you use the orthocenter and the orthoradius of the projection of the weighted convex hull. After the initial decomposition process, the molecular surfaces are represented as composite trimmed surfaces. Each one of them is composed of one base NURBS surface and a few NURBS curves which trim the base NURBS surface. Thus, the raw trimmed surfaces typically have several sides and internal holes. For the unstructured meshes, we need to apply usual techniques for the generation of the tetrahedra. As for the hierarchical case, more care must be considered.

In most interesting cases where the probe radius is about 1.4 Angstrom and the number of atoms are more than 50, there are typically many trimmed surfaces (toroidal and spherical) which are very thin. As a consequence, the usual approach about decomposing the trimmed surfaces into fine triangles which are then tetrahedralized does not function anymore. More accurately, that is still valid but it produces tetrahedra which are undesirable for subsequent numerical simulation. What we need are block tetrahedra which are large enough and well shaped. For some related works, follow the next links:

Unstructured models: unstrucured 1 , unstrucured 2 , unstrucured 3 , unstrucured 4 .
Hierarchical models: hierarchical 1 , hierarchical 2 , hierarchical 3 .
Clementine for Connolly modeling
Molecular modeling from PDB
Publication

For the unstructured case, we treat a single mesh for the whole simulation domain. The simulation is thus performed directly on the Ω. As for the multiblock version, our treatment has the next properties: We have a nested sequence of finite sequence of finite dimensional spaces:

We do not use a single mesh for the whole model but rather many parametric meshes (one for each curved tetrahedra),
The results in a coarser structure is used to enhance those in a fine one,
In the course of the refinement, we keep track of the the parent/children hierarchy: which tetrahedra of V_k correspond to elements of V_k+1,
We organize our data structure and code structure efficiently so as to treat the hierarchical refinements patch by patch.

SEQUENCE OF NONLINEAR PROBLEMS:

A related work is the method of Xie et al:
Xie and Zhou: A new minimization protocol for solving nonlinear Poisson-Boltzmann mortar finite element equation, (BIT Numerical Mathematics (2007) 47: 853�871)
This is a very good source for the mortar case but it has serious drawbacks. For Xie et al, in order to solve the nonlinear optimization, they evaluate on the fly the integrals

by using some numerical quadrature where ψ_j is some FE-bases. As a consequence, the evaluations of the objective functional, gradient and the Hessian are either very expensive or inaccurate. That method is still feasible in the 2D case if the model is very small as they considered. To remedy those drawbacks, our PBE nonlinear solver is featured by

Quasi-Newton based on BFGS,

Efficient and accurate line search,

Fast-Fourier Transform Hessian,

Nonlinear correction with fast convergence.

We use a method avoiding numerical quadrature on the fly for functional, gradient or Hessian evaluations. We replace the initial problem by a sequence of nonlinear problems which are simpler to solve. If we assume that we have an approximation u_h^(k), the next approximation is obtained by solving for the correction function

leading to

For the solving of the simplified nonlinear problem for c^(k)_h, we use an improved quasi-Newton method which is not computationally expensive since there are no integrals to be estimated. At first sight, the above sum looks expensive but it turns out that only very few terms of that sum are necessary to obtain a good accuracy.
Nonlinear solver for unstructured mesh:
On the right plots, you see some instances of the convergence for the case of unstructured meshes. We consider various values of the parameters seen in the type of Poisson-Boltzmann equation.

The convergence has the following typical feature: it is dragging down slowly until it reaches a certain closeness and then it drops very quickly. In case of hierarchy the convergence is much faster because the starting slow convergence is removed. Below are some results for the hierarchical case.

In the above figures, we show some results on a nonadaptive unstructured mesh composed only of two atoms for illustration purpose. The first figure show the exact solution. The following one is the result of the linear Poisson-Boltzmann equation which is very inaccurate. The remaining figures show some intermediate nonlinear results and the ultimate nonlinear approximation. As for the hierarchical case, the nonlinear solver follow the same line just as in the multiblock case with the exception that you have the cascading advantage and the multiple level results.

The right plots contain some outcomes for a model when enhanced with the hierarchical strategy. The different curves are obtained for various number of blocks on some levels. The increase factor indicates the increase of the number of tetrahedra on the inter-level hierarchical refinement. It is observed that on all tests that we have done the converge is much faster than what it observed above in the case of dispensing from hierarchical structure.

BFGS / FHT (BROYDEN-FLETCHER-GOLDFARB-SHANNO / FAST HARTLEY TRANSFORM):

The usual quasi-Newton use only an approximate Hessian. One of the widely used approximate Hessian updates is the BFGS which are of two kinds: direct one and the inverse one. They update respectively the Hessian and the inverse of the Hessian at every iteration. The direct BFGS still necessitates a solving of the approximate Hessian at each iteration while the inverse one needs only one matrix-vector multiplication. Still, the inverse BFGS is very much flawed by the initialization because you need to have a pretty good approximation of the initial inverse Hessian. In some problems, such an initial good inverse is easy to obtain. In our case of Poisson-Boltzmann, it is very expensive to assemble even the inverse of the initial Hessian. Therefore, a direct BFGS is more reasonable in our case. The drawbacks of using the direct BFGS update formula are:

The resulting BFGS-matrices are not sparse. Hence, one would need O(n²) memory space for the allocation.

Although the matrix from the BFGS updates is symmetric positive definite, its non-sparsity makes it difficult to solve using fast solvers like CG or GMRES.

Although the matrices are easy to update, the resulting matrices are usually badly conditioned resulting to many iteration counts

A single matrix-vector multiplication for a dense BFGS-matrix of order n which represents the number of Finite Element unknowns is very expensive.

One would like to keep the advantage of BFGS (easy update, no need to assemble full Hessian) while trying to avoid the above drawbacks.

The main idea of the FFT-BFGS is to obtain an updating formula whose results are easy to treat. First you need only O(n) memory for the storage. On the other hand, the matrix to be solved is diagonal after a few DFT or DHT transforms which can be done very fast. The usual FFT based BFGS completely ignores the initial matrices and makes an update. In our problem some part of the Hessian is available and sparse. As a consequence, we apply only the FFT-BFGS update to the non-sparse information.

Line search: the most usual ways of line searches are Armijo, Goldstein, Powell, Cauchy, etc. For our simplified nonlinear method, the line search equation is given as a univariate polynomials of low degree. Searching the line search solution amounts therefore to selecting between the few roots, which can be computed exactly, the one minimizing the objective functional in the current search direction.

On the left plot, we see some convergence history of the resulting method. Minimizing the functional necessarily amounts to detecting the position of zero gradients. As a consequence, measuring the gradient norm an in the plot is a good gauge of assessing the quality of the nonlinear solver. We observe that the Quasi-Newton behaves very efficiently. The resulting FFT-QUASI-Newton almost admits the convergence property of the exact Newton.

The above result admits the same property for both the unstructured and hierarchical settings. Recall that the implementation of the DFT/DHT by the FFT method is already hierarchical. The Cooley-Tuckey algorithm already implicitly 'halfen' the data at every iteration.

ADAPTIVE REFINEMENT AND HIERARCHY:

We consider a-posteriori error estimates which can be computed very efficiently. Although the initial problem is a nonlinear one, the a-posteriori error estimator is a linear one. Based upon the strengthned Cauchy-Schwartz inequality, one can deduce an a-posteriori error estimator which is both efficient and reliable. In other words, one can derive the next property

where the constants c₁ and c₂ are independent of the discretization properties of the underlying mesh for two local bases V(T) and Z(T). The stiffness matrix with respect to the bases of V(T) and Z(T) is given in M as follows. In function of its block matrices, the original strengthned Cauchy-Schwarz constant is expressed otherwise

Therefore, γ is given by the square root of the largest eigenvalue of the generalized eigenproblem

Here are some instance of the eigenfunctions of some well chosen local decomposition:

FORTHCOMING TASKS:

Parallel implementation: this extension is not difficult (but it takes time) since everything in the sequencial version is already set blockwise. Hence, the parallelization seems to be only a distribution of the blocks to the processors.
Currently, the choice of polynomial degree is done by heuristic. More mathematical way is needed to choose the degree as well as to decide between mesh refinement and polynomial degree elevation.
The a-posteriori error estimator always work in practice but the analysis is only partially finished for the piecewise linear case.
Improvement and optimization of the code.
Papers.