Adaptive multilevel finite element solution of the Poisson–Boltzmann equation I. Algorithms and examples

This article is the first of two articles on the adaptive multilevel finite element treatment of the nonlinear Poisson–Boltzmann equation (PBE), a nonlinear eliptic equation arising in biomolecular modeling. Fast and accurate numerical solution of the PBE is usually difficult to accomplish, due to the presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domain, and rapid (exponential) nonlinearity. In this first article, we explain how adaptive multilevel finite element methods can be used to obtain extremely accurate solutions to the PBE with very modest computational resources, and we present some illustrative examples using two well‐known test problems. The PBE is first discretized with piece‐wise linear finite elements over a very coarse simplex triangulation of the domain. The resulting nonlinear algebraic equations are solved with global inexact Newton methods, which we have described in an article appearing previously in this journal. A posteriori error estimates are then computed from this discrete solution, which then drives a simplex subdivision algorithm for performing adaptive mesh refinement. The discretize–solve–estimate–refine procedure is then repeated, until a nearly uniform solution quality is obtained. The sequence of unstructured meshes is used to apply multilevel methods in conjunction with global inexact Newton methods, so that the cost of solving the nonlinear algebraic equations at each step approaches optimal O(N) linear complexity. All of the numerical procedures are implemented in MANIFOLD CODE (MC), a computer program designed and built by the first author over several years at Caltech and UC San Diego. MC is designed to solve a very general class of nonlinear elliptic equations on complicated domains in two and three dimensions. We describe some of the key features of MC, and give a detailed analysis of its performance for two model PBE problems, with comparisons to the alternative methods. It is shown that the best available uniform mesh‐based finite difference or box‐method algorithms, including multilevel methods, require substantially more time to reach a target PBE solution accuracy than the adaptive multilevel methods in MC. In the second article, we develop an error estimator based on geometric solvent accessibility, and present a series of detailed numerical experiments for several complex biomolecules. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1319–1342, 2000     

Multigrid solution of the Poisson—Boltzmann equation

A multigrid method is presented for the numerical solution of the linearized Poisson–Boltzmann equation arising in molecular biophysics. The equation is discretized with the finite volume method, and the numerical solution of the discrete equations is accomplished with multiple grid techniques originally developed for twodimensional interface problems occurring in reactor physics. A detailed analysis of the resulting method is presented for several computer architectures, including comparisons to diagonally scaled CG, ICCG, vectorized ICCG and MICCG, and to SOR provided with an optimal relaxation parameter. Our results indicate that the multigrid method is superior to the preconditioned CG methods and SOR and that the advantage of multigrid grows with the problem size. © 1993 John Wiley & Sons, Inc.     

An explicit linear six-step method with vanished phase-lag and its first derivative

An explicit linear sixth algebraic order six-step method with vanished phase-lag and its first derivative is constructed in this paper. We will study the method theoretically and computationally. Theoretical investigation contains the building of the method, the calculation of the local truncation error, the comparative error analysis of the new method with the method with constant coefficients and the stability analysis of the new method using scalar test equation with different frequency than the frequency of the scalar test equation used for the development of the method. Computational investigation contains the application of the new obtained linear six-step method to the resonance problem of the radial time independent Schrödinger equation. The theoretical and computational study lead us to the summary that the new proposed linear scheme is computationally and theoretically more efficient than other well known methods for the numerical solution of the Schrödinger equation and related periodic initial or boundary value problems.     

Hybrid boundary element and finite difference method for solving the nonlinear Poisson-Boltzmann equation

A hybrid approach for solving the nonlinear Poisson-Boltzmann equation (PBE) is presented. Under this approach, the electrostatic potential is separated into (1) a linear component satisfying the linear PBE and solved using a fast boundary element method and (2) a correction term accounting for nonlinear effects and optionally, the presence of an ion-exclusion layer. Because the correction potential contains no singularities (in particular, it is smooth at charge sites) it can be accurately and efficiently solved using a finite difference method. The motivation for and formulation of such a decomposition are presented together with the numerical method for calculating the linear and correction potentials. For comparison, we also develop an integral equation representation of the solution to the nonlinear PBE. When implemented upon regular lattice grids, the hybrid scheme is found to outperform the integral equation method when treating nonlinear PBE problems. Results are presented for a spherical cavity containing a central charge, where the objective is to compare computed 1D nonlinear PBE solutions against ones obtained with alternate numerical solution methods. This is followed by examination of the electrostatic properties of nucleic acid structures.     

A SOLUTION TO THE UPWARD MIGRATION OF AN INTRUSION

The evolution of the upward migration of the magma is a nonlinear and unstable problem in mathematics. It is difficult to solve it. And using the numerical method, the solution is relatively tedious and time-consuming. This paper introduces a method of the instantaneous point source to solve the linear and unstable heat conduction equation during the infinite period of time instead of the solution of the nonlinear and unstable heat conduction equation. The results obtained by this method coincide with those by the numerical method, meaning that this method offers a simple way to solve the nonlinear and unstable heat conduction equation.     

Assessment of linear finite‐difference Poisson–Boltzmann solvers

CPU time and memory usage are two vital issues that any numerical solvers for the Poisson–Boltzmann equation have to face in biomolecular applications. In this study, we systematically analyzed the CPU time and memory usage of five commonly used finite‐difference solvers with a large and diversified set of biomolecular structures. Our comparative analysis shows that modified incomplete Cholesky conjugate gradient and geometric multigrid are the most efficient in the diversified test set. For the two efficient solvers, our test shows that their CPU times increase approximately linearly with the numbers of grids. Their CPU times also increase almost linearly with the negative logarithm of the convergence criterion at very similar rate. Our comparison further shows that geometric multigrid performs better in the large set of tested biomolecules. However, modified incomplete Cholesky conjugate gradient is superior to geometric multigrid in molecular dynamics simulations of tested molecules. We also investigated other significant components in numerical solutions of the Poisson–Boltzmann equation. It turns out that the time‐limiting step is the free boundary condition setup for the linear systems for the selected proteins if the electrostatic focusing is not used. Thus, development of future numerical solvers for the Poisson–Boltzmann equation should balance all aspects of the numerical procedures in realistic biomolecular applications. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Linear augmented Slater-type orbital method for free standing clusters

We have developed a Scalable Linear Augmented Slater-Type Orbital (LASTO) method for electronic-structure calculations on free-standing atomic clusters. As with other linear methods we solve the Schr?dinger equation using a mixed basis set consisting of numerical functions inside atom-centered spheres and matched onto tail functions outside. The tail functions are Slater-type orbitals, which are localized, exponentially decaying functions. To solve the Poisson equation between spheres, we use a finite difference method replacing the rapidly varying charge density inside the spheres with a smoothed density with the same multipole moments. We use multigrid techniques on the mesh, which yields the Coulomb potential on the spheres and in turn defines the potential inside via a Dirichlet problem. To solve the linear eigen-problem, we use ScaLAPACK, a well-developed package to solve large eigensystems with dense matrices. We have tested the method on small clusters of palladium.     

Flexible polyelectrolyte simulations at the Poisson-Boltzmann level: a comparison of the kink-jump and multigrid configurational-bias Monte Carlo methods

We present a new approach for simulating the motions of flexible polyelectrolyte chains based on the continuous kink-jump Monte Carlo technique coupled to a lattice field theory based calculation of the Poisson-Boltzmann (PB) electrostatic free energy "on the fly." This approach is compared to the configurational-bias Monte Carlo technique, in which the chains are grown on a lattice and the PB equation is solved for each configuration with a linear scaling multigrid method to obtain the many-body free energy. The two approaches are used to calculate end-to-end distances of charged polymer chains in solutions with varying ionic strengths and give similar numerical results. The configurational-bias Monte Carlo/multigrid PB method is found to be more efficient, while the kink-jump Monte Carlo method shows potential utility for simulating nonequilibrium polyelectrolyte dynamics.     

Numerical simulations of phase separation dynamics in a water-oil-surfactant system

We have studied numerically the dynamics of the microphase separation of a water-oil-surfactant system. We developed an efficient and accurate numerical method for solving the two-dimensional time-dependent Ginzburg-Landau model with two order parameters. The numerical method is based on a conservative, second-order accurate, and implicit finite-difference scheme. The nonlinear discrete equations were solved by using a nonlinear multigrid method. There is, at most, a first-order time step constraint for stability. We demonstrated numerically the convergence of our scheme and presented simulations of phase separation to show the efficiency and accuracy of the new algorithm.     

A study of one dimensional nonlinear diffusion equations by Bernstein polynomial based differential quadrature method

The purpose of this work is to study numerical solutions of nonlinear diffusion equations such as Fisher’s equation, Burgers’ equation and modified Burgers’ equation by applying Bernstein differential quadrature method (BDQM). These nonlinear diffusion equations occur in many applications of theoretical, engineering and environmental sciences. Therefore, finding numerical solutions of these equations is very important. In BDQM, Bernstein polynomials are used as base functions to find weighting coefficients of differential quadrature method. We have applied BDQM on eleven different test problems taken from the literature and the computed results confirm that BDQM is an efficient method for finding solution of nonlinear partial differential equations. It is found that BDQM produces very good results even at small number of grid points.     

Multiple grid methods for classical molecular dynamics

Presented in the context of classical molecular mechanics and dynamics are multilevel summation methods for the fast calculation of energies/forces for pairwise interactions, which are based on the hierarchical interpolation of interaction potentials on multiple grids. The concepts and details underlying multigrid interpolation are described. For integration of molecular dynamics the use of different time steps for different interactions allows longer time steps for many of the interactions, and this can be combined with multiple grids in space. Comparison is made to the fast multipole method, and evidence is presented suggesting that for molecular simulations multigrid methods may be superior to the fast multipole method and other tree methods.     

The wavelet methods to linear and nonlinear reaction–diffusion model arising in mathematical chemistry

In this paper, we have applied an accurate and efficient wavelet scheme (due to Legendre polynomial) to find the numerical solutions for a set of coupled reaction–diffusion equations. This technique provides the solutions in rapid convergence series with computable terms for the problems with high degree of non linear terms appearing in the governing differential equations. The highest derivative in the differential equation is expanded into wavelet series, this approximation is then integrated while the boundary conditions are applied by using integration constants. With the help of operational matrices, the nonlinear reaction–diffusion equations are converted into a system of algebraic equations. Finally, some numerical examples to demonstrate the validity and applicability of the method have been furnished. The use of Legendre wavelets is found to be accurate, efficient, simple, and computationally attractive. This wavelet method can be used for obtaining quick solution in many chemical Engineering problems.     

A new explicit Bessel and Neumann fitted eighth algebraic order method for the numerical solution of the Schrödinger equation

An explicit eighth algebraic order Bessel and Neumann fitted method is developed in this paper for the numerical solution of the Schrödinger equation. The new method has free parameters which are defined in order the method is fitted to spherical Bessel and Neumann functions. A variable-step procedure is obtained based on the newly developed method and the method of Simos [17]. Numerical illustrations based on the numerical solution of the radial Schrödinger equation and of coupled differential equations arising from the Schrödinger equation indicate that this new approach is more efficient than other well known methods.     

Computation of dynamic adsorption with adaptive integral,finite difference,and finite element methods

Analysis of diffusion-controlled adsorption and surface tension in one-dimensional planar coordinates with a finite diffusion length and a nonlinear isotherm, such as the Langmuir or Frumkin isotherm, requires numerical solution of the governing equations. This paper presents three numerical methods for solving this problem. First, the often-used integral (I) method with the trapezoidal rule approximation is improved by implementing a technique for error estimation and choosing time-step sizes adaptively. Next, an improved finite difference (FD) method and a new finite element (FE) method are developed. Both methods incorporate (a). an algorithm for generating spatially stretched grids and (b). a predictor-corrector method with adaptive time integration. The analytical solution of the problem for a linear dynamic isotherm (Henry isotherm) is used to validate the numerical solutions. Solutions for the Langmuir and Frumkin isotherms obtained using the I, FD, and FE methods are compared with regard to accuracy and efficiency. The results show that to attain the same accuracy, the FE method is the most efficient of the three methods used.     

A new time dependent approach for solving electrochemical interfaces Part I: theoretical considerations using lie group analysis

In this paper we introduce a nonlinear partial differential equation (nPDE) of the third order to the first time. This new model equation allows the extension of the Debye-Hückel Theory (DHT) considering time dependence explicitly. This also leads to a new formulation of the meaning of the nonlinear Poisson-Boltzmann Equation (PBE) and therefore we call it the modified Poisson-Boltzmann Equation (mPBE). In the present first part of this extensive study we derive the equation from the electromagnetics from a quasistatic perspective, or more precisely the electroquasistatic approximation (EQS). Our main focus will be the analysis via the Lie group formalism and since that up to now no symmetry calculation is available we believe that it seems indispensable to apply this method yielding a deeper insight into the behaviour of the solution manifold of this new equation following electrochemical considerations. We determine the classical Lie point symmetries including algebraic properties. Similarity solutions in a most general form and suitable nonlinear transformations are obtained. In addition, a note relating to potential and generalized symmetries is drawn. Moreover we show how the equation leads to approximate symmetries and we apply the method to the first time. The second part appearing shortly after will deal with algebraic solution methods and we shall show that closed-form solutions can be calculated without any numerical methods. Finally the third part will consider appropriate electrochemical experiments proving the model under consideration.     

A new coupled wavelet-based method applied to the nonlinear reaction–diffusion equation arising in mathematical chemistry

In this paper, we have applied the wavelet-based coupled method for finding the numerical solution of Murray equation. To the best of our knowledge, until now there is no rigorous Legendre wavelets solution has been reported for the Murray equation. The highest derivative in the differential equation is expanded into Legendre series, this approximation is integrated while the boundary conditions are applied using integration constants. With the help of Legendre wavelets operational matrices, the Murray equation is converted into an algebraic system. Block pulse functions are used to investigate the Legendre wavelets coefficient vectors of nonlinear terms. The convergence of the proposed method is proved. Finally, we have given a numerical example to demonstrate the validity and applicability of the method. Moreover the use of proposed wavelet-based coupled method is found to be simple, efficient, less computation costs and computationally attractive.     

An explicit four-step method with vanished phase-lag and its first and second derivatives

In this paper we will develop an explicit fourth algebraic order four-step method with phase-lag and its first and second derivatives vanished. The comparative error and the stability analysis of the above mentioned paper is also presented. The new obtained method is applied on the resonance problem of the Schrödinger equationIn order in order to examine its efficiency. The theoretical and the computational results shown that the new obtained method is more efficient than other well known methods for the numerical solution of the Schrödinger equation and related initial-value or boundary-value problems with periodic and/or oscillating solutions.     

Adaptive multilevel finite element solution of the Poisson–Boltzmann equation II. Refinement at solvent‐accessible surfaces in biomolecular systems

We apply the adaptive multilevel finite element techniques (Holst, Baker, and Wang 21 ) to the nonlinear Poisson–Boltzmann equation (PBE) in the context of biomolecules. Fast and accurate numerical solution of the PBE in this setting is usually difficult to accomplish due to presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domains, and rapid (exponential) nonlinearity. However, these adaptive techniques have shown substantial improvement in solution time over conventional uniform‐mesh finite difference methods. One important aspect of the adaptive multilevel finite element method is the robust a posteriori error estimators necessary to drive the adaptive refinement routines. This article discusses the choice of solvent accessibility for a posteriori error estimation of PBE solutions and the implementation of such routines in the “Adaptive Poisson–Boltzmann Solver” (APBS) software package based on the “Manifold Code” (MC) libraries. Results are shown for the application of this method to several biomolecular systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1343–1352, 2000     

Exponentially - Fitted Multiderivative Methods for the Numerical Solution of the Schrödinger Equation

In this paper exponentially fitted multiderivative methods are developed for the numerical solution of the one-dimensional Schrödinger equation. The methods are called multiderivative since uses derivatives of order two and four. An application to the the resonance problem of the radial Schrödinger equation indicates that the new method is more efficient than other similar well known methods of the literature.     

A powerful truncated Newton method for potential energy minimization

With advances in computer architecture and software, Newton methods are becoming not only feasible for large-scale nonlinear optimization problems, but also reliable, fast and efficient. Truncated Newton methods, in particular, are emerging as a versatile subclass. In this article we present a truncated Newton algorithm specifically developed for potential energy minimization. The method is globally convergent with local quadratic convergence. Its key ingredients are: (1) approximation of the Newton direction far away from local minima, (2) solution of the Newton equation iteratively by the linear Conjugate Gradient method, and (3) preconditioning of the Newton equation by the analytic second-derivative components of the “local” chemical interactions: bond length, bond angle and torsional potentials. Relaxation of the required accuracy of the Newton search direction diverts the minimization search away from regions where the function is nonconvex and towards physically interesting regions. The preconditioning strategy significantly accelerates the iterative solution for the Newton search direction, and therefore reduces the computation time for each iteration. With algorithmic variations, the truncated Newton method can be formulated so that storage and computational requirements are comparable to those of the nonlinear Conjugate Gradient method. As the convergence rate of nonlinear Conjugate Gradient methods is linear and performance less predictable, the application of the truncated Newton code to potential energy functions is promising.