Numerical solution of the nonlinear Poisson–Boltzmann equation: Developing more robust and efficient methods

We present a robust and efficient numerical method for solution of the nonlinear Poisson-Boltzmann equation arising in molecular biophysics. The equation is discretized with the box method, and solution of the discrete equations is accomplished with a global inexact-Newton method, combined with linear multilevel techniques we have described in an article appearing previously in this journal. A detailed analysis of the resulting method is presented, with comparisons to other methods that have been proposed in the literature, including the classical nonlinear multigrid method, the nonlinear conjugate gradient method, and nonlinear relaxation methods such as successive overrelaxation. Both theoretical and numerical evidence suggests that this method will converge in the case of molecules for which many of the existing methods will not. In addition, for problems which the other methods are able to solve, numerical experiments show that the new method is substantially more efficient, and the superiority of this method grows with the problem size. The method is easy to implement once a linear multilevel solver is available and can also easily be used in conjunction with linear methods other than multigrid. © 1995 by John Wiley & Sons, Inc.     

Multigrid methods in density functional theory

A multigrid method for real-space solution of the Kohr-Sham equations is presented. By using this multiscale approach, the problem of critical slowing down typical of iterative real-space solvers is overcome. The method scales linearly in computer time with the number of electrons if the orbitals are localized. Here, we describe details of our multigrid method, present preliminary many-electron numerical results illustrating the efficiency of the solver, and discuss its strengths and limitations. © 1997 John Wiley & Sons, Inc.     

Multigrid high-order mesh refinement techniques for composite grid electrostatics calculations

A new method for performing high-order mesh-refinement multigrid computations in real space is presented. The method allows for accurate linear scaling electrostatics calculations over composite domains with local nested fine patches. The Full Approximation Scheme (FAS) multigrid technique is utilized for a sequence of refinement patches of increasing resolution. Conservation forms are generated on coarse scales by additional defect correction terms that counter the local excess fluxes at the boundaries. Formulas are given for arbitrary order, extending the existing technique of Bai and Brandt. Test calculations are presented for a singular source in three dimensions that illustrate the multigrid convergence properties, numerical accuracy, and correct order of the approach. © 1999 John Wiley & Sons, Inc. J Comput Chem 20: 1731–1739, 1999     

Mulitstep methods with vanished phase-lag and its first and second derivatives for the numerical integration of the Schrödinger equation

A tenth algebraic order eight-step method is developed in this paper. For this method  we require the phase-lag and its first and second derivatives to be vanished. A comparative error analysis and a comparative stability analysis are also presented in this paper. The new proposed method is applied for the numerical solution of the one-dimensional Schrödinger equation. The efficiency of the new methodology is proved via the theoretical analysis and the numerical applications. General conclusions about the importance of several properties on the construction of numerical algorithms for the approximate solution of the radial Schrödinger equation are also presented.     

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     

Multi-multigrid solution of modified Poisson–Boltzmann equation for arbitrarily shaped molecules

A new multi-multigrid method is presented for solving the modified Poisson–Boltzmann equation based on the Kirkwood Hierarchy of equations, with Loeb's closure, on a three-dimensional grid. The results are compared with standard Poisson–Boltzmann calculations, which are known to underestimate the local concentration of counterions near charged parts of molecules, mainly due to neglect of fluctuations in the ionic concentrations. In the present study, the Kirkwood hierarchy of equations is discretized with the finite volume method and solved using multigrid techniques. The new possibility for solution of the three-dimensional modified Poisson–Boltzmann equation, for the first time within a model including a dielectric discontinuity, and within reasonable computational time, enables the calculation of higher valence ion distributions around arbitrarily shaped biological macromolecules. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 893–901, 1998     

A two-dimensional Haar wavelet approach for the numerical simulations of time and space fractional Fokker–Planck equations in modelling of anomalous diffusion systems

In this paper, the numerical solution for the fractional order Fokker–Planck equation has been presented using two dimensional Haar wavelet collocation method. Two dimensional Haar wavelet method is applied to compute the numerical solution of nonlinear time- and space-fractional Fokker–Planck equation. The approximate solutions of the nonlinear time- and space-fractional Fokker–Planck equation are compared with the exact solutions as well as solutions available in open literature. The present scheme is very simple, effective and convenient for obtaining numerical solution of the time and space-fractional Fokker–Planck equation.     

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.     

A fast adaptive multigrid boundary element method for macromolecular electrostatic computations in a solvent

A fast multigrid boundary element (MBE) method for solving the Poisson equation for macromolecular electrostatic calculations in a solvent is developed. To convert the integral equation of the BE method into a numerical linear equation of low dimensions, the MBE method uses an adaptive tesselation of the molecular surface by BEs with nonregular size. The size of the BEs increases in three successive levels as the uniformity of the electrostatic field on the molecular surface increases. The MBE method provides a high degree of consistency, good accuracy, and stability when the sizes of the BEs are varied. The computational complexity of the unrestricted MBE method scales as O(N_at), where N_at is the number of atoms in the macromolecule. The MBE method is ideally suited for parallel computations and for an integrated algorithm for calculations of solvation free energy and free energy of ionization, which are coupled with the conformation of a solute molecule. The current version of the 3-level MBE method is used to calculate the free energy of transfer from a vacuum to an aqueous solution and the free energy of the equilibrium state of ionization of a 17-residue peptide in a given conformation at a given pH in ∼ 400 s of CPU time on one node of the IBM SP2 supercomputer. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18: 569–583, 1997     

The multiscale coarse-graining method. IX. A general method for construction of three body coarse-grained force fields

The multiscale coarse-graining (MS-CG) method is a method for constructing a coarse-grained (CG) model of a system using data obtained from molecular dynamics simulations of the corresponding atomically detailed model. The formal statistical mechanical derivation of the method shows that the potential energy function extracted from an MS-CG calculation is a variational approximation for the true potential of mean force of the CG sites, one that becomes exact in the limit that a complete basis set is used in the variational calculation if enough data are obtained from the atomistic simulations. Most applications of the MS-CG method have employed a representation for the nonbonded part of the CG potential that is a sum of all possible pair interactions. This approach, despite being quite successful for some CG models, is inadequate for some others. Here we propose a systematic method for including three body terms as well as two body terms in the nonbonded part of the CG potential energy. The current method is more general than a previous version presented in a recent paper of this series [L. Larini, L. Lu, and G. A. Voth, J. Chem. Phys. 132, 164107 (2010)], in the sense that it does not make any restrictive choices for the functional form of the three body potential. We use hierarchical multiresolution functions that are similar to wavelets to develop very flexible basis function expansions with both two and three body basis functions. The variational problem is solved by a numerical technique that is capable of automatically selecting an appropriate subset of basis functions from a large initial set. We apply the method to two very different coarse-grained models: a solvent free model of a two component solution made of identical Lennard-Jones particles and a one site model of SPC/E water where a site is placed at the center of mass of each water molecule. These calculations show that the inclusion of three body terms in the nonbonded CG potential can lead to significant improvement in the accuracy of CG potentials and hence of CG simulations.     

Fluorescent microscopic determination of proteins in human serum with the self-ordered ring of nuclear fast red formed on the solid support of glass slides

A self-ordered ring (SOR) technique based on the assembly of fluorescent molecules on the solid support of glass slides is presented for the detection of trace amount of proteins by using fluorescence microscope. At pH 6.62 and with the aid of poly(vinyl alcohol)-124 (PVA-124), a droplet of nuclear fast red (NFR) solution can form a fluorescent SOR on hydrophobic glass slide, and the presence of proteins can enhance its fluorescence intensity. When the volume of the droplet is 0.2 μl, linear response is observed in the range 0-100 pg, and the limit of determination is 2.44 pg for BSA and 0.91 pg for HSA (3σ). The results of determination for three human serum samples were identical with those obtained according to the Bradford method using Coomassie Brilliant Blue (CBB G-250).     

Multiscale coarse-graining and structural correlations: connections to liquid-state theory

A statistical mechanical framework elucidates the significance of structural correlations between coarse-grained (CG) sites in the multiscale coarse-graining (MS-CG) method (Izvekov, S.; Voth, G. A. J. Phys. Chem. B 2005, 109, 2469; J. Chem. Phys. 2005, 123, 134105). If no approximations are made, then the MS-CG method yields a many-body multidimensional potential of mean force describing the interactions between CG sites. However, numerical applications of the MS-CG method typically employ a set of pair potentials to describe nonbonded interactions. The analogy between coarse-graining and the inverse problem of liquid-state theory clarifies the general significance of three-particle correlations for the development of such CG pair potentials. It is demonstrated that the MS-CG methodology incorporates critical three-body correlation effects and that, for isotropic homogeneous systems evolving under a central pair potential, the MS-CG equations are a discretized representation of the well-known Yvon-Born-Green equation. Numerical calculations validate the theory and illustrate the role of these structural correlations in the MS-CG method.     

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 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.     

Numerical study of electroosmotic slip flow of fractional Oldroyd-B fluids at high zeta potentials

In this paper, an investigation of the electroosmotic flow of fractional Oldroyd-B fluids in a narrow circular tube with high zeta potential is presented. The Navier linear slip law at the walls is considered. The potential field is applied along the walls described by the nonlinear Poisson–Boltzmann equation. It's worth noting here that the linear Debye–Hückel approximation can't be used at the condition of high zeta potential and the exact solution of potential in cylindrical coordinates can't be obtained. Therefore, the Matlab bvp4c solver method and the finite difference method are employed to numerically solve the nonlinear Poisson–Boltzmann equation and the governing equations of the velocity distribution, respectively. To verify the validity of our numerical approach, a comparison has been made with the previous work in the case of low zeta potential and the excellent agreement between the solutions is clear. Then, in view of the obtained numerical solution for the velocity distribution, the numerical solutions of the flow rate and the shear stress are derived. Furthermore, based on numerical analysis, the influence of pertinent parameters on the potential distribution and the generation of flow is presented graphically.     

A New Fluorescein Isothiocyanate-Piperidine Self-ordered Ring Phosphorescent Probe for the Determination of Trace Protein

The ? COOH in fluorescein isothiocyanate (FITC) reacted with ? NH? in piperidine (P) to form FITC‐P on the center of indentation of polyamide membrane (PAM) when drying for 2 min at (92±1)°C. Then, the FITC‐P diffused outward from the indentation center and formed the round SOR‐P‐FITC (containing the FITC‐P self‐ordered rings). Thus, multi‐FITC accumulated on SOR‐P‐FITC, leading to the enhancement of RTP signal on bio‐target, whose I_p increased 2.0 times compared with non‐generated SOR. When bovine serum albumin (BSA) was added to the center of SOR‐P‐FITC, ? NCS of FITC in SOR‐P‐FITC reacted with ? NH₂ of BSA to form SOR‐P‐FITC‐BSA, which caused the RTP signal of FITC to enhance sharply. The ΔI_p of the system was 3.4 times higher than that without β‐CD and 4.0 times higher than that without SOR‐P‐FITC formed. Its ΔI_p was linear to the content of BSA. Therefore, a new solid substrate‐room temperature phosphorimetry (SS‐RTP) for the determination of trace protein was established using SOR‐P‐FITC as a phosphorescent probe. Under the optimum condition, the linear range of this method was 0.040–16.0 ag·spot^?1 with a detection limit (LD) of 8.5 zg·spot^?1 (0.40 µL sample solution per spot, the corresponding concentration was 2.1×10^?17 g·mL^?1), and the regression equation of working curve was ΔI_p=3.848+4.240m_BSA (ag·spot^‐1), n=6, correlation coefficient (r) was 0.9993. This method with high sensitivity had been applied to determining the content of trace protein in the water samples, and the results coincided well with those obtained with pyrocatechol violet‐Mo(VI) method (P.V.M.M.). At the same time, the mechanism of SS‐RTP using SOR‐P‐FITC as a phosphorescent probe (SOR‐P‐FITC‐SS‐RTP) was discussed.     

Potential distribution in the solution interface of a scanning tunneling microscope

The linearized Poisson—Boltzmann equation is solved in the region between a sphere and a plane, which is modelling the electrolyte solution interface between the tip and the substrate in a scanning tunneling microscope. A series expansion in modified Bessel functions and Legendre polynomials, which are solutions to the linearized Poisson—Boltzmann equation, is used to fit the boundary conditions. Another numerical method of finite difference is also used with the domain transformed into bispherical coordinates. Results for cases of different potential values on the boundary surfaces and different distances of the sphere from the plane are presented.     

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.     

The multiscale coarse-graining method. VIII. Multiresolution hierarchical basis functions and basis function selection in the construction of coarse-grained force fields

The multiscale coarse-graining (MS-CG) method is a method for determining the effective potential energy function for a coarse-grained (CG) model of a molecular system using data obtained from molecular dynamics simulation of the corresponding atomically detailed model. The coarse-grained potential obtained using the MS-CG method is a variational approximation for the exact many-body potential of mean force for the coarse-grained sites. Here we propose a new numerical algorithm with noise suppression capabilities and enhanced numerical stability for the solution of the MS-CG variational problem. The new method, which is a variant of the elastic net method [Friedman et al., Ann. Appl. Stat. 1, 302 (2007)], allows us to construct a large basis set, and for each value of a so-called "penalty parameter" the method automatically chooses a subset of the basis that is most important for representing the MS-CG potential. The size of the subset increases as the penalty parameter is decreased. The appropriate value to choose for the penalty parameter is the one that gives a basis set that is large enough to fit the data in the simulation data set without fitting the noise. This procedure provides regularization to mitigate potential numerical problems in the associated linear least squares calculation, and it provides a way to avoid fitting statistical error. We also develop new basis functions that are similar to multiresolution Haar functions and that have the differentiability properties that are appropriate for representing CG potentials. We demonstrate the feasibility of the combined use of the elastic net method and the multiresolution basis functions by performing a variational calculation of the CG potential for a relatively simple system. We develop a method to choose the appropriate value of the penalty parameter to give the optimal basis set. The combined effect of the new basis functions and the regularization provided by the elastic net method opens the possibility of using very large basis sets for complicated CG systems with many interaction potentials without encountering numerical problems in the variational calculation.