In this paper, we present a new regularization scheme for the linearized Poisson–Boltzmann equation (PBE) which models the electrostatic potential of biomolecules in a solvent. This scheme is based on the splitting of the target potential into the short- and long-range components localized in the molecular region by using the range-separated (RS) tensor format [P. Benner, V. Khoromskaia, and B. N. Khoromskij, SIAM J. Comput., 2 (2018), pp. A1034–A1062] for representation of the discretized multiparticle Dirac delta [B. N. Khoromskij, J. Comput. Phys., 401 (2020), 108998] constituting the highly singular right-hand side in the PBE. From the computational point of view our regularization approach requires only the modification of the right-hand side in the PBE so that it can be implemented within any open-source grid-based software package for solving PBE that already includes some FEM/FDM disretization scheme for elliptic PDE and solver for the arising linear system of equations. The main computational benefits are twofold. First, one applies the chosen PBE solver only for the smooth long-range (regularized) part of the collective potential with the regular right-hand side represented by a low-rank RS tensor with a controllable precision. Thus, we eliminate the numerical treatment of the singularities in the right-hand side and do not change the interface and boundary conditions. And second, the elliptic PDE need not be solved for the singular part in the right-hand side at all, since the short-range part of the target potential of the biomolecule is precomputed independently on a computational grid by simple one-dimensional tensor operations. The total potential is then obtained by adding the numerical solution of the PBE for the smooth long-range part to the directly precomputed tensor representation for the short-range contribution. Numerical tests illustrate that the new regularization scheme, implemented by a simple modification of the right-hand side in the chosen PBE solver, improves the accuracy of the approximate solution on rather coarse grids. The scheme also demonstrates good convergence behavior on a sequence of refined grids.

The Poisson–Boltzmann equation (PBE) is a nonlinear elliptic parametrized partial differential equation that arises in biomolecular modeling and is a fundamental tool for structural biology. It is used to calculate electrostatic potentials around an ensemble of fixed charges immersed in an ionic solution. Efficient numerical computation of the PBE yields a high number of degrees of freedom in the resultant algebraic system of equations, ranging from several hundred thousand to millions. Coupled with the fact that in most cases the PBE requires to be solved multiple times for a large number of system configurations, for example, in Brownian dynamics simulations or in the computation of similarity indices for protein interaction analysis, this poses great computational challenges to conventional numerical techniques. To accelerate such onerous computations, we suggest to apply the reduced basis method (RBM) and the (discrete) empirical interpolation method ((D)EIM) to the PBE with a special focus on simulations of complex biomolecular systems, which greatly reduces this computational complexity by constructing a reduced order model (ROM) of typically low dimension. In this study, we employ a simple version of the PBE for proof of concept and discretize the linearized PBE (LPBE) with a centered finite difference scheme. The resultant linear system is solved by the aggregation-based algebraic multigrid (AGMG) method at different samples of ionic strength on a three-dimensional Cartesian grid. The discretized LPBE, which we call the high-fidelity full order model (FOM), yields solution as accurate as other LPBE solvers. We then apply the RBM to the FOM. DEIM is applied to the Dirichlet boundary conditions which are nonaffine in the parameter (ionic strength), to reduce the complexity of the ROM. From the numerical results, we notice that the RBM reduces the model order from $${\mathcal {N}} = 2\times 10^{6}$$to $$N = 6$$at an accuracy of $${\mathcal {O}}(10^{-9})$$and reduces the runtime by a factor of approximately 7600. DEIM, on the other hand, is also used in the offline-online phase of solving the ROM for different values of parameters which provides a speed-up of 20 for a single iteration of the greedy algorithm.