A Self-Adjusting FEM-BEM Coupling Scheme for the Nonlinear Poisson-Boltzmann Equation

TL;DR

This work tackles the nonlinear Poisson–Boltzmann equation in biomolecular electrostatics by introducing a self-adjusting FEM–BEM coupling that confines nonlinearities to a near-solute region via a three-region model. A Newton–Raphson solver, augmented with cubic initial approximations and an automated relaxation factor ω^opt derived from a one-dimensional root, achieves robust, fast convergence without user intervention. Validation against APBS on spherical geometries and extensive tests on RNA and highly charged biomolecules demonstrate accurate energies and substantial speedups (up to 1.37×) with practical automation. The methodology is broadly applicable to other nonlinear Laplace-type problems and can be extended to further optimisations and other nonlinear PDEs in physics and chemistry.

Abstract

The Poisson-Boltzmann equation is widely used to model molecular electrostatics; however, it is usually solved in linearised form because the sinh nonlinearity is challenging, limiting its applicability in highly charged systems such as nucleic acids. This work presents a solution method for the nonlinear Poisson-Boltzmann equation based on a coupled finite/boundary element scheme that automatically finds an optimal relaxation parameter, ensuring fast and reliable convergence of the nonlinear solver without user intervention. We validated our solver against APBS for a spherical cavity, and used RNA-based structures to perform a thorough study of the different algorithmic choices, and to test our implementation. We found that the best alternative to solve the Poisson-Boltzmann equation was using a Newton-Raphson method where the nonlinearity was gradually introduced with a cubic approximation in the first iteration. Newton-Raphson was also the best method to find the optimal relaxation factor, reducing the number of iterations by 40%. Including other optimisation techniques, we were able to obtain a 1.37x speed-up with respect to the best hand-picked relaxation factor for 1HC8 (molecule with highest charge in our tests), avoiding any trial-and-error process to find the relaxation factor.

Figures (8)

• Figure 1: Two-region domain split according to Equation \ref{['eq:Variables']}.
• Figure 2: Domain split into three regions according to Equation \ref{['eq:Variables_2Superficie']}.
• Figure 3: Different surfaces around the molecule. The van der Waals (VdW) and solvent accessible (SAS) surfaces are not used in this work. The SES corresponds to $\Gamma_m$ and the dotted line is $\Gamma_s$. The Test Radius corresponds to the size of one water molecule ($r$=1.4 Å).
• Figure 4: Cross section of an example FEM mesh (in grey), with the SES identified within (fine mesh region close to the interface). The BEM mesh (equivalent to the outer boundary of the FEM mesh) is identified in blue.
• Figure 5: Solvation energy of a sphere for different values of $\kappa$ in linear and nonlinear FEM-BEM compared to APBS.