Fully implicit ADI schemes for solving the nonlinear Poisson-Boltzmann equation

Abstract The Poisson-Boltzmann (PB) model is an effective approach for the electrostatics analysis of solvated biomolecules. The nonlinearity associated with the PB equation is critical when the underlying electrostatic potential is strong, but is extremely difficult to solve numerically. In this paper, we construct two operator splitting alternating direction implicit (ADI) schemes to efficiently and stably solve the nonlinear PB equation in a pseudo-transient continuation approach. The operator splitting framework enables an analytical integration of the nonlinear term that suppresses the nonlinear instability. A standard finite difference scheme weighted by piecewise dielectric constants varying across the molecular surface is employed to discretize the nonhomogeneous diffusion term of the nonlinear PB equation, and yields tridiagonal matrices in the Douglas and Douglas-Rachford type ADI schemes. The proposed time splitting ADI schemes are different from all existing pseudo-transient continuation approaches for solving the classical nonlinear PB equation in the sense that they are fully implicit. In a numerical benchmark example, the steady state solutions of the fully-implicit ADI schemes based on different initial values all converge to the time invariant analytical solution, while those of the explicit Euler and semi-implicit ADI schemes blow up when the magnitude of the initial solution is large. For the solvation analysis in applications to real biomolecules with various sizes, the time stability of the proposed ADI schemes can be maintained even using very large time increments, demonstrating the efficiency and stability of the present methods for biomolecular simulation.


Introduction
With the development of theoretical methods and computational techniques in the past few decades, molecular modeling has become a more and more effective and practical approach to mimic the behavior of molecules with biological significance. However, the atomic level description of molecular systems involving both a specific molecule and its aquatic and ionic surrounding is prohibitively expensive when it comes to investigating large conformational reorganization of biomolecules, such as in protein folding [11]. Alternatively, with the implicit Poisson-Boltzmann (PB) model one may explicitly include atomic details of the solute while implicitly treat the solvent by its mean influence and the mobile ions with a statistical mechanics distribution [2]. Therefore, the PB model substantially reduces the degree of freedom of the system. In the PB model, the governing PB equation of electrostatic potential is a nonlinear elliptic equation defined on multiple domains with discontinuous coefficients across the solute-solvent interface or molecular surface [5,6]. The PB equation admits analytical solutions only for simple geometries such as spheres [22] or rods [18]. For molecules with complex geometry, the PB equation can only be solved through numerical approximations. Solving the PB equation numerically is challenging due to the consideration of non-smooth solution, discontinuous coefficients (dielectric constants and ionic constants), irregular interface, singular source terms, infinite domain, as well as nonlinearity when ionic effects are strong. There exist various numerical PB solvers and these solvers can be roughly but not completely classified into two categories: 1) The mesh-based finite difference/finite element methods [1,10,20,25,27,29]; and 2) The boundary element method [7, 15-17, 21, 24, 35]. These PB solvers are designed based on different numerical algorithms and are aimed to applications from different perspectives. For example, the PB solvers embedded in molecular modeling packages such as Delphi [29], CHARMM [20], AMBER [25], and APBS [1] all employ finite difference discretizations and smoothened molecular surfaces. Although arguably these solvers have reduced accuracy, particularly on or near the molecular surface, the efficient, robust and user-friendly features of these PB solvers brought them popularity among the bio-oriented community. Mesh-based interface methods such as the Immersed Interface Method (IIM) [27] and Matched Interface and Boundary Poisson-Boltzmann (MIBPB) solver [10] can significantly improve the accuracy by rigorously treating the numerical difficulties such as discontinuities and singularities. However, the complexity of the algorithms to a certain extent reduces the efficiency. The boundary element methods have the potential to circumvent many of the numerical difficulties, and accelerate the solver with fast algorithms such as the fast multipole method (FMM) [7,24,35] and treecode [15]. However, the nonlinearity of the PB equation is cumbersome to address under the boundary element formulation. A hybrid boundary element approach which solves the nonlinear term by the finite difference is reported in [8]. In a usual approach for solving the nonlinear PB (NPB) equation, a nonlinear algebraic system is typically formed, after a spatial discretization by means of finite difference or finite element AU + N(U) = b (1) where U is the vector of electrostatic potentials, b is the vector of source terms, and A and N are, respectively, discrete operators for the linear and nonlinear parts of the NPB equation. Various algebraic methods can then be employed to solve Eq. (1) efficiently. A straightforward way is using the nonlinear relaxation method, in which the nonlinear term will be evaluated according to the previous approximations. One then needs to solve a linear system in each relaxation step. The commonly used linear iterative methods, such as Gauss-Seidal or Jacobi, can be employed as the internal iterative solver. The widely used biomolecular simulation packages such as Delphi [29] and CHARMM [20] employ the nonlinear relaxation methods. The nonlinear conjugate gradient method solves Eq. (1) approximately by calculating a minimization problem with a nonlinear integral functional, and is built in the UHBD package [26]. In the inexact Newton method [19], the Jacobian system of Eq. (1) is solved inexactly to find a descent direction first. A line search is then conducted along the descent direction to ensure global convergence. The performance of the inexact Newton method is greatly accelerated in the APBS package, where a fast algebraic multigrid solver is combined with the inexact Newton method [18,19]. A comprehensive assessment of various algebra-based NPB solvers is carried out in [9]. The nonlinearity treatments related to the present study are pseudo-transient continuation approaches [31,32] for solving the NPB equation. In such approaches, a time dependent NPB equation is constructed by introducing a pseudo-time.
The solution of the original nonlinear boundary value system is essentially recovered by the steady state solution of the pseudo-time dependent process. The first pseudo-transient continuation approach for solving the NPB equation is introduced in [32], in which the pseudo-time convection diffusion process is solved by means of a diffusion module of an existing finite element software. To guarantee a large time step so that the long time integration can be computed quickly, the implicit Euler scheme is employed together with a linearization technique based on the first order Taylor expansion. The linearization gives rise to a coefficient matrix of the linear system at each time step by evaluating the nonlinear term at the previous time instant [32]. Since the nonlinear term of the NPB equation is treated explicitly, the overall time integration is semi-implicit. A more efficient semi-implicit approach is introduced in [ the time dependent NPB process is decomposed into diffusion, convection, and nonlinear parts. The convection and nonlinear parts are also treated explicitly, while the diffusion term is integrated implicitly by a regularized alternating direction implicit (ADI) method [13]. The regularization technique allows an even larger time increment [13]. Moreover, the Thomas algorithm [33] can be employed to solve the tridiagonal finite difference systems efficiently in each ADI step. Consequently, the semi-implicit approach [31] becomes even faster. However, a very large time increment is still prohibited in semi-implicit time integrations of pseudo-transient continuation approaches, because the nonlinear term of the NPB equation is typically a hyperbolic sine function which could be exponentially large. The objective of this paper is to introduce fully implicit ADI methods for solving the NPB equation in a pseudo-transient continuation process, so that a very large time increment can be employed. The proposed ADI methods are built based on an operator splitting or time splitting framework, in which the nonlinear subsystem of the NPB equation can be analytically integrated. Such a time splitting procedure is often found to be a powerful tool in solving various time dependent partial differential equations (PDEs) [3,4,36,39] due to its capability in enlarging the Courant-Friedrich-Lewy (CFL) number. The proposed fully implicit ADI schemes are different from the semi-implicit ADI scheme of [31] in terms of two features. First, the time-dependent NPB equation is rewritten into a convection-diffusion process in [31], while it is treated as a nonhomogeneous diffusion process in the present study. In particular, standard finite difference discretization is conducted for the entire nonhomogeneous diffusion term in the present implicit ADI integration. The resulting linear algebraic system is still tridiagonal. Second, the nonlinear term is evaluated at the previous time step in [31], while it is integrated analytically in the present operator splitting ADI schemes. Because of these two distinct features, the proposed ADI schemes become fully implicit.
In a related study, we have also developed operator splitting ADI schemes [38] for solving a generalized NPB equation [37]. Even though the schemes presented in [38] are similar to the present ones, the underlying generalized NPB equation [37] differs significantly from the classical NPB equation. The pursuit of the operator splitting ADI schemes on the classical NPB equation distinguishes the present schemes from the previous work on generalized NPB in [38]. The most distinguishing feature of the generalized NPB equation [37] is the use of a smoothly defined solute-solvent interface via some characteristic functions in the Eulerian formulation [12,34]. This smooth definition, on the one hand, makes the stability issue of the nonlinear term more severe, because the nonlinear term now takes nonzero values within the so-called ion-exclusion layer [18,37,38]. Note that in the classical NPB equation, the nonlinear term is defined outside the ion-exclusion layer through a Heaviside step function. Within the ion-exclusion layer, the electrostatic potential could be large by Coulomb's law, so that an extremely large value results after evaluating the nonlinear term. Consequently, the explicit Euler scheme is found to be unconditionally unstable in solving the time dependent generalized NPB equation [37,38], whereas the explicit Euler scheme remains to be conditionally stable in the present study for the time dependent classical nonlinear PB equation. On the other hand, the smooth definition of solute-solvent interface induces a smooth dielectric profile; thus the central finite difference can be simply employed for the spatial discretization of the generalized NPB equation with good stability. Nevertheless, the dielectric profile is discontinuous in the classical NPB equation and the potential function loses its regularity across the material interface. In the present study, a standard finite difference discretization is adopted for simplicity and the non-smooth feature of the potential function is believed to introduce additional instability in our computation. Sophisticated interface treatments [14,40] to account for the discontinuous dielectric profile might reduce the instability associated with the spatial discretization. However, such studies are beyond the scope of the present paper. The rest of this paper is organized as follows. In Section 2, we first introduce a pseudo- has the form of: − ∇ · ( (r)∇ (r)) +κ 2 (r) sinh ( (r)) = ρ (r) (2) where ρ is a singular source The dielectric constant is assumed to be piecewise values such that (r) = for r ∈ Ω and (r) = for r ∈ Ω , in the units of fundamental charge is the partial charge on the th atom of the solute located at r , and B is the Boltzmann constant. Hereκ is the modified Debye-Huckel parameter, which is defined as for r ∈ Ω where N A is Avogadro's number and I is the ionic strength in the unit of mole. Numerically, when T = 298K, we haveκ 2 = 8 486902807 Å −2 I and 2 B T = 332 06364/0 592183 Å [18]. It is noted the form of the nonlinear PB equation (2) varies on different units. The present form (2) is the dimensionless form as explained in [19]. The dimensionless electrostatic potential can be conveniently converted to the electrostatic potential of units / / by multiplying the constant 0.592183 subject to room temperature (T = 298K ) [18]. Physically, the boundary condition for the dimensionless electrostatic potential is defined at infinity: For mesh-based numerical computations, a finite domain Ω is required. The approximated analytical condition (6) can be employed to assign the boundary condition on ∂Ω Equation (6) is actually a linear superposition of Coulomb's law for a series of N partial charges at positions r and can be used to approximate the potential solved from Eq. (2) when ∂Ω is sufficiently distanced from the macromolecule subdomain. On the other hand, we have the interface jump condition across Γ, where [ ] Γ is defined to be the difference between two limiting values approaching from outside and inside, and n is the outward normal direction.

Pseudo-time governing equation and numerical difficulties
Theoretical modeling and computational simulation based on the full nonlinear PB system are highly nontrivial. This motivates the development of pseudo-transient continuation approaches for solving the nonlinear PB equation [31,32,37,38]. Essentially, a pseudo-transient variation is introduced to convert Eq. (2) from the time-independent form to a time-dependent form Typically, one needs to first specify an initial solution, which could be the electrostatic potential solved from a linearized PB equation [37] or trivially = 0. One then numerically integrates Eq. (8) for a sufficiently long time period. The solution of the original nonlinear PB equation (2) is essentially recovered by the steady state solution of the pseudo-time dependent process (8). However, there exist great difficulties in the numerical integration of the time dependent NPB equation (8). Generally speaking, because a long time integration is required for (8), explicit time stepping methods are usually not efficient for pseudo-transient continuation approaches [31,32,37,38]. In the literature, semi-implicit time stepping methods [31,32] are commonly employed to solve time dependent NPB equation (8) so that a large time step could be used for a stable simulation. Nevertheless, a fully implicit time integration method has never been constructed for solving the classical nonlinear PB equation (2). Such a development is mainly hindered by the implicit treatment of the hyperbolic sine nonlinear term.

Numerical discretization
Consider a uniform mesh partition of the computational domain Ω. Without the loss of generality, we assume the equal grid spacing in all , and directions. We denote the time increment to be ∆ . To facilitate the following discussions, we apply the notation = ( ) to denote the electrostatic potential at node ( ). To gain an in-depth understanding, we first discuss a standard numerical discretization to the time dependent NPB equation (8). In particular, the central finite difference and explicit Euler scheme are used for spatial and temporal discretization, respectively. The detailed discretization is given as ) is the fractional charge at grid point ( ), which is obtained by using the trilinear interpolation to distribute all charges in the source term ρ of Eq. (2). It is noted that the spatial discretization presented in (9) is a standard central difference scheme for the nonhomogeneous Laplacian term ∇ · ( (r)∇ ). This central scheme is of second order of accuracy when the dielectric profile (r) is smooth. Accuracy reduction is inevitable for a discontinuous dielectric profile as in the present study. To restore the second order of accuracy, sophisticated interface treatments are required to enforce interface conditions (7) into the spatial discretization [14,40]. Nevertheless, since the main objective of the present study is to develop fully implicit time integration schemes, we will continue using the simple finite difference scheme shown in (9). This ensures that our current development can be conveniently adopted in other finite difference schemes originated from different physical backgrounds. In the present study, the value of on half grid nodes, such as ( , is determined by its location ( , which is either inside/on or outside the molecular surface Γ. In particular, ( ∈ Ω , and ( A small ∆ is usually required in the discretization (9) so that the explicit Euler scheme is very inefficient. Essentially, in order to maintain the time stability, the difference between +1 and shall be small. Thus, near the solutesolvent interface Γ where the dimensionless electrostatic potential could take a large value [37,38], the nonlinear term sinh( ) is exponentially large. Consequently, a very small ∆ has to be used in the explicit Euler scheme for a stable simulation. Therefore, for a long time integration to reach the steady state, the explicit Euler scheme must be computationally very expensive. This difficulty can be partially relieved when using semi-implicit time integration schemes [31,32], but is not completely resolved because the nonlinear term is still treated explicitly in [31,32], i.e., evaluated at the present time level = . The fully implicit schemes constructed by simply evaluating the nonlinear term at = +1 will generate a nonlinear algebraic system to be solved at each time step, which is obviously time-consuming, too.

Operator splitting alternating direction implicit (ADI) schemes
To overcome the aforementioned difficulties, we will develop novel operator splitting alternating direction implicit (ADI) schemes for solving the time dependent NPB equation (8), so that the nonlinear term can be handled fully implicitly. We note that the proposed ADI schemes are similar to those developed in the pseudo-transient continuation solution of the generalized NPB equation [38]. Nevertheless, it is the first time in the literature that such schemes are constructed for solving the classical nonlinear PB equation (2), to the best of our knowledge.  (8), we will assume time-dependence only for . We note that all other functions in (8), i.e., ,κ 2 , and ρ are time independent. 2.4.1. ADI scheme 1 In the ADI scheme 1 or ADI1, at each time step from to +1 , the time dependent NPB equation (8) will be solved by a first order time splitting method in two stages [36] ∂ We then have U +1 = V +1 . For the first stage, Eq. (10) can be analytically solved, In other words, with W at , W +1 can be calculated analytically. Thus the difficulties caused by the nonlinear term sinh(·) in Eq. (2) are circumvented. The right-hand side of (12) is essentially a function of W and ∆ . We thus rewrite Eq. (12) as W +1 = F (W ; ∆ ) to facilitate the following discussions. A Douglas-Rachford type ADI scheme is then proposed to solve the nonhomogeneous diffusion equation (11). The discretization of Eq. (11) using backward-Euler integration in time and central differencing in space results in where δ 2 , δ 2 , and δ 2 are the central difference operators in the , , and directions, respectively, The Douglas-Rachford ADI scheme for Eq. (11) is formulated as [33] ( In other words, the three dimensional linear algebraic system in the implicit scheme (13) is decomposed into one dimensional linear algebraic systems in (14). Moreover, each of these linear systems has a tridiagonal structure and thus can be efficiently solved by the Thomas algorithm [5,33]. By eliminating * and * * and solving for +1 in (14), we obtain Thus, the Douglas-Rachford scheme (14) is a higher order perturbation of the backward Euler scheme (13). Since both the backward Euler scheme and the time spitting scheme (10) - (11) are first order in time, the proposed ADI1 scheme is of first order accuracy in time. In the numerical simulations, the same Dirichlet boundary values are assumed for , * and * * as for . The entire ADI1 time integration is fully implicit. In the ADI scheme 2 or ADI2, at each time step from to +1 , the time dependent NPB equation (8) will be solved by a second order time splitting method in three stages [36] ∂ We then have U +1 =W +1 . As in the ADI1 scheme, in the first and last stage of the ADI2 scheme, an analytical integration can be conducted. Symbolically, we have W +1 = F (W ; ∆ 2 ) andW +1 = F (W ; ∆ 2 ), where F is defined as in Eq. (12). A Douglas type ADI scheme is proposed to solve the nonhomogeneous diffusion equation (17) in the second stage. The discretization of (17) using Crank-Nicolson integration in time and central differencing in space results in This can be decomposed into , , and directions to give a Douglas ADI scheme for (17) [33] (1 − ∆ 2 Similarly, the tridiagonal systems in (20) can be efficiently solved by the Thomas algorithm. The elimination of * and * * in (20) leads to It can be shown that the resulting equation is a higher order perturbation of the Crank-Nicolson scheme (19) in three space variables. Thus, the temporal order of accuracy of the Douglas scheme (20) is two, which is consistent with the present time splitting (16) - (18). Therefore, the proposed fully implicit ADI2 scheme is of second order of accuracy in time. Again, the boundary values of , * and * * are the same as those for . In terms of the computational complexity of one cycle, the ADI2 scheme involves extra flops compared with the ADI1 scheme. Nevertheless, the ADI2 scheme is usually more accurate than the ADI1 scheme.

Semi-implicit ADI1 scheme
For a comparison, we also consider the counterparts of the proposed two ADI schemes with semi-implicit time integration. The semi-implicit ADI1 (semi-ADI1) scheme is constructed based on the backward-Euler integration of the time dependent NPB equation (8), but with the nonlinear term evaluated at By treating the nonlinear term in (22) as an additional source term, the same Douglas-Rachford ADI discretization used in the ADI1 scheme can be utilized to form the semi-ADI1 scheme

Semi-implicit ADI2 scheme
Similarly, the semi-implicit ADI2 (semi-ADI2) scheme is constructed based on the Crank-Nicolson integration of the time dependent NPB equation (8), By treating the nonlinear term in (24) as a source term, the same Douglas ADI discretization used in the ADI2 scheme leads to the semi-ADI2 scheme We note that the present semi-ADI1 and semi-ADI2 schemes are different from the existing semi-implicit approaches [31,32]. The primary difference is in spatial discretization. The simplest central difference discretization is used in the present semi-implicit ADI schemes, while more complex finite difference and finite element formulas are used in [31,32]. However, the explicit treatment of the nonlinear term is the same in the existing and present semi-implicit ADI schemes. We will illustrate how this explicit treatment will affect the stability in the next section.

Convergence criteria
We propose to approximate the solution to the nonlinear PB equation (2) by using the steady state solution to the time-dependent nonlinear PB equation (8), which can be achieved in theory when → ∞. Numerically we will need a stop criterion when is sufficiently large. Two convergence criteria will be adopted simultaneously in this paper to maintain the numerical efficiency.

Criterion 1
Time integration stops when the change in the solvation energy of two consecutive time steps is less than a given tolerance. This criterion is particularly useful when calibrating the proposed schemes in constructed examples with analytical solutions.

Criterion 2
Time integration stops when ≥ T , where T is a user specified time. After many numerical trials, we find T = 10 is a reasonable choice for terminating the time integration in the solution of the time-dependent NPB equation (8). This usually produces adequately converged solutions for real protein systems.

Numerical validations
In this section, we validate our proposed ADI schemes numerically. We first solve the nonlinear PB equation on a sphere and compare the numerical results with the available analytical solution for tests of stability as well as both spatial and temporal accuracy and convergence. We next apply our new algorithms to compute the solvation energy for a series of proteins with various sizes and geometric structures. All simulations are compiled with Intel visual Fortran compiler and run on a Dell 6600 with i7-2760QM 2.4GHz CPU and 16GB memory.

A spherical cavity with analytical solution
For the nonlinear PB equation, the analytical solutions are available only for simple and special geometries. We here design a case for a spherical cavity based on the reference [30] (r) = where ε = / and R is the radius of the sphere. By defining the source term ρ of Eq. (2) according to Eq. (27), it can be verified that the solution in Eq. (26) satisfies Eq.
(2) as well as the interface jump conditions in Eq. (7). Note that there are both singularity and non-smoothness involved in the present analytical solution. The singularity is from the source term as seen in Eq. (27) and the non-smoothness is due to the interface jump conditions in Eq. (7). Both features will significantly reduce the accuracy of the spatial discretization and introduces additional instability in temporal discretization. In our study of the spherical case, the chosen sphere has radius R = 1 Å with one centered charge of 1 . The dielectric constants are set as = 80 and = 1. The nonlinear constantκ is set as 1.

Convergence and stability
We first examine the convergence of the pseudo-transient continuation approach. For this purpose, we solve the time dependent nonlinear PB equation (8) with different initial solutions (r 0) until the stopping criteria are satisfied. We then compare our steady state solutions with the time independent analytical solution (26) to verify if the convergence is achieved. Besides the trivial initial solution (r 0) = 0, we specifically construct a family of initial solutions (r 0) = H cos( with different magnitude H for the study of the sensitivities of the schemes subject to the nonlinearity (the sinh( ) term). It can be seen that the initial solution (28) satisfies the boundary conditions for the computational domain We considered various H values ranging from 1 to 20 for explicit Euler, semi-ADI1, semi-ADI2, ADI1, and ADI2 simulations for a comparison. It is found that under different initial solutions, all tested schemes converge to the analytical solution when is sufficiently large, except for a few extreme cases where the convergence is ruined by the instability. This validates the present pseudo-transient continuation approach. We next quantitatively analyze the stability of the explicit Euler, semi-ADI1, semi-ADI2, ADI1, and ADI2 schemes. Because of the singularity and non-smoothness features of the solution, all time integration schemes here are conditionally stable, i.e., there is a maximum stable time increment ∆ = 2 at different spatial mesh size . We report the critical values of these schemes in Table 1. From this table, we see that with the increase of H in the initial solutions (28), and thus the increase of the nonlinearity, the critical ∆ values of the explicit Euler scheme and semi-implicit ADI schemes become exponentially smaller and approach zero eventually, while those of ADI schemes are essentially unchanged. Therefore, the fully-implicit schemes demonstrated much improved stability compared with the explicit and semi-implicit schemes.

Spatial convergence rate
We then study the spatial convergence rate of ADI schemes with the same example. The simulation results for the spherical case can be found in Table 2. Denoting as the numerical solution, we use the following measures to estimate errors in a relative manner:   In producing results in Table 2, we vary spatial increment and set ∆ = 2 /20 to fulfill the stability requirement. The convergence rate is nonuniform because the charge singularities (the delta function) and the non-smoothness of the solution across the interface Γ are treated in an approximated sense for the present standard finite difference discretization. However, we could still see in average approximately the first-order convergence in L ∞ and L 2 norms of both ADI1 and ADI2 schemes. This approximately first order convergence rate can also be found in reference [14] when the linearized PB equation on a sphere is solved using a similar spatial discretization. To further demonstrate this convergence behavior, we provide detailed error-mesh plots in Fig. 1. In these plots, we perform a least square fitting (the solid lines) of the errors against the number of mesh points N, which is inversely proportional to . The slopes of the solid straight lines are calculated and listed in the legends of both plots. The magnitude of reflects the average convergence rate. It can be seen that both ADI1 and ADI2 schemes deliver the approximately first-order convergence in both L ∞ and L 2 norms. Visually, the errors of both ADI1 and ADI2 schemes are almost the same when N is large. Nevertheless, there are some minor differences when N is small.

Temporal convergence rate
We finally focus on the temporal order of convergence, which is reported in Table 3. Since the spatial discretization error is large due to the singularity and non-smoothness as described previously, a direct comparison between numerical solutions at different and the analytical solution will not be able to reveal the temporal order. Instead, we fix = 0 125 and use the electrostatic potentials solved at the finest time increment ∆ =2.5E-05 as a reference solution. The potentials computed based on other ∆ values are benchmarked with this reference solution. This enables us to compute the temporal discretization error through essentially the cancellation of the spatial discretization error. We observe second order temporal convergence from the table for both ADI1 and ADI2 schemes. For a fixed ∆ , the ADI2 is more accurate than the ADI1. We note some better-than-expected performance for the ADI1 scheme, i.e., the numerical order   is higher than the theoretical one. In fact, the same finding has been observed for the ADI1 scheme when it is applied to solve the generalized PB equation defined on a smooth solute-solvent boundary and a conventional analysis of the temporal convergence is considered [38].

Solvation energy of proteins
In this subsection, we validate the proposed time splitting ADI schemes by solving the nonlinear PB equation for real protein systems. In particular, we will compute solvation energies on a series of proteins with different size and geometric structures. For a comparison, we also report the results produced by the explicit Euler method.  The solvation energy or free energy of solvation is the energy released when the solute in free space is dissolved in solvent. In molecular simulation, it can be computed by evaluating the difference between total free energy of the solute in the vacuum and in the solvent. In the context that only the electrostatic effects are considered, the solvation energy can be defined as where φ and φ 0 in the units of kcal/mol/ are electrostatic potentials in the presence of the solvent and the vacuum, respectively. They are obtained by scaling the dimensionless potentials with the constant 0.592183 subject to room temperature (T = 298K ) [18]. Based on the trilinear interpolation of the singular charges in ρ , the solvation energy can be calculated as We solve the nonlinear PB equation (2) on a set of 24 proteins by using explicit Euler, ADI1 and ADI2 schemes. For simulation on proteins, the dielectric constants are set as = 80 and = 1. The ionic strengths I is set as 0.15M, which makes the nonlinear constantκ = 0 1257. To quantitatively justify the accuracy, efficiency and stability, we compute the solvation energy based on electrostatic potentials solved from the nonlinear PB equation. The numerical results of the ADI1, ADI2, and Euler methods are listed in Table 4. The corresponding parameters are listed in the caption of the table and the choice of the ∆ is subject to the stability requirement. From this table, we can see that all three nonlinear methods produce consistent solvation energies. In particular, the explicit Euler scheme is expected to produce some of the most accurate results, because a much smaller ∆ is applied and the same spatial discretization is used. The second order time discretization ADI2 is supposed to yield more accurate results than the first order ADI1. Indeed, the results in the table adequately showed this, i.e., the results of the ADI2 are closer to those of the Euler scheme. On the other hand, the solution of the Linearized PB equation (LPBE) is carried out by using the same finite difference discretization and a biconjugate gradient solver. The solvation free energies of LPBE results are also listed in Table 4. It can be seen that there are some large differences between linear and nonlinear PB models, reflecting the importance of the nonlinear term in Eq. (2). We also depict the simulation results in Fig. 2. To avoid the scaling problem caused by the large variation of the solvation energy value, in Fig. 2(a) we plot the differences between nonlinear PB results generated by ADI1, ADI2, and Euler schemes and linearized PB results. This figure shows the differences between linear and nonlinear models. It also demonstrates the consistency among three nonlinear PB solvers. We then plot the CPU time used by three algorithms for solving the nonlinear PB equation in Fig. 2(b). It is seen that the explicit Euler method costs about 2-3 times more CPU time than both ADI schemes. Finally, we plot the surface potentials of a protein (1mbg) in Fig. 3. Besides the nonlinear PB results generated by the ADI1 scheme, we also compare with the linear PB results obtained by the same finite difference discretization. From

Conclusion
This paper presents the first fully implicit time integration scheme for solving the nonlinear Poisson-Boltzmann (NPB) equation in a pseudo-transient continuation framework. In such a framework, the boundary value solution of the NPB equation is essentially recovered by the steady state solution to the pseudo-time NPB process. The long time integration in this process necessitates the development of implicit schemes that allow a large time increment ∆ . However, all existing implicit methods are semi-implicit, because of the difficulty in the implicit treatment of the nonlinear hyperbolic sine term. This bottleneck is overcome in the present paper by using a novel operator splitting procedure. After splitting, the nonlinear instability can be completely avoided through analytical integration, while both Douglas and Douglas-Rachford alternating direction implicit (ADI) schemes are developed for solving the nonhomogeneous diffusion subsystem. A standard finite difference scheme is employed for spatial discretization. Other more accurate finite difference methods can be applied, as long as the tridiagonal structure is preserved in each coordinate direction. The Thomas algorithm is utilized to solve linear algebraic systems efficiently in the proposed fully implicit ADI schemes. The proposed operator splitting ADI schemes are first validated by solving the NPB equation on spherical cavities where time invariant analytical solutions are available subject to specially designed source terms. By considering different initial values, it is demonstrated that all steady state solutions of the fully-implicit ADI methods converge to the analytical solution, while those of the explicit Euler scheme and semi-implicit ADI methods blow up when the magnitude of the initial solution is large. The spatial convergence rate of the proposed ADI schemes is found to be nearly first order, because the accuracy of the standard finite difference is reduced by the singularities from the source term and non-smoothness of the solution across the molecular surface. Special cancellation is considered in order to reveal the true temporal order of accuracy. It is interesting to observe that both the proposed ADI schemes then yield a second order in time. The application of both ADI schemes to calculate solvation free energies for a series of 24 proteins is considered next. In all of our experiments, it is found that the fully implicit ADI schemes allow much larger ∆ values than that of the explicit Euler scheme, so that the present biomolecular simulation becomes faster. Moreover, since the proposed ADI schemes can withstand a very strong nonlinear effect, the full prediction power of the nonlinear solvation model can be numerically accomplished. In particular, we illustrate the theoretical difference between the linearized Poisson-Boltzmann model and the present NPB model in electrostatic analysis. The development of more robust spatial discretization associated with the proposed temporal schemes are under our consideration.