NONLINEAR POISSONBOLTZMANN: 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 solutesolvent 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,
 Aposteriori 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
Ω=Ω^{u}UΩ
^{v}:
The coefficients ε(x) and κ(x) are
spacedependent functions related to the dielectric
value and the modified DebyeHueckle 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 counterions of type i and q_{i} is its charge
while β:=1/(k_{B}T) with k_{B} and T being the
Boltzmann constant and the temperature.


For real chemical simulation, the righthand 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)
PoissonBoltzmann equation:
Such an approximation of
sinh yields
very inaccurate results for highly charged
electrostatic potentials. As a consequence, we consider the
nonlinear PoissonBoltzmann 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_{1}(m_{1},r_{1}) and B_{2}(m_{2},r_{2}), 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^{4} 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:
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
PoissonBoltzmann 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
FEbases. 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
 QuasiNewton based on BFGS,
 Efficient and accurate line search,
 FastFourier 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 quasiNewton 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 PoissonBoltzmann
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
PoissonBoltzmann 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
interlevel 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 (BROYDENFLETCHERGOLDFARBSHANNO / FAST HARTLEY TRANSFORM):

The usual quasiNewton 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 matrixvector
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 PoissonBoltzmann, 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 BFGSmatrices are not sparse. Hence, one would need
O(n^{2}) memory space for the allocation.
 Although the matrix from the BFGS
updates is symmetric positive definite, its nonsparsity
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 matrixvector multiplication for a dense BFGSmatrix 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 FFTBFGS 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 FFTBFGS update to the
nonsparse 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 QuasiNewton behaves very efficiently. The resulting FFTQUASINewton 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 CooleyTuckey algorithm already implicitly 'halfen'
the data at every iteration.
ADAPTIVE REFINEMENT AND HIERARCHY:

We consider aposteriori error estimates which can be computed very efficiently.
Although the initial problem is a nonlinear one, the aposteriori error estimator
is a linear one. Based upon the strengthned CauchySchwartz inequality, one
can deduce an aposteriori error estimator which is both
efficient
and
reliable. In other words, one can derive the next property
where the constants
c_{1} and
c_{2} 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 CauchySchwarz 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:
 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 aposteriori 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.