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.     

Highly accurate and simple analytical approach to nonlinear Poisson–Boltzmann equation

A novel method is suggested to analytically solve a nonlinear Poisson–Boltzmann (NLPB) equation. The method consists chiefly of reducing the NLPB equation to linear PB equation in several segments by approximating a free term of the NLPB equation by piecewise linear functions, and then, solving analytically the linear PB equation in each segment. Superiority of the method is illustrated by applying the method to solve the NLPB equation describing a colloid sphere immersed in an arbitrary valence and mixed electrolyte solution; extensive test indicates that the resulting analytical expressions for both the electrical potential distribution Ψ (r) and surface charge density/surface potential relationship (σ/Ψ ₀) are characterized with two properties that mathematical structures are much simpler than those previously reported and application scope can be arbitrarily wide by adjusting the linear interpolation range. Finally, it is noted that the method is “universal” in that its applications are not limited to the NLPB equation.     

Quadrature method of moments for aggregation-breakage processes

Investigation of particulate systems often requires the solution of a population balance, which is a continuity statement written in terms of the number density function. In turn, the number density function is defined in terms of an internal coordinate (e.g., particle length, particle volume) and it generates integral and derivative terms. Different methods exist for numerically solving the population balance equation. For many processes of industrial significance, due to the strong coupling between particle interactions and fluid dynamics, the population balance must be solved as part of a computational fluid dynamics (CFD) simulation. Such an approach requires the addition of a large number of scalars and the associated transport equations. This increases the CPU time required for the simulation, and thus it is clear that it is very important to use as few scalars as possible. In this work the quadrature method of moments (QMOM) is used. The QMOM has already been validated for crystal growth and aggregation; here the method is extended to include breakage. QMOM performance is tested for 10 different cases in which the competition between aggregation and breakage leads to asymptotic solutions.     

Adsorption integral equation via complex approximation with constraints: The Langmuir kernel

The relationship between the measured adsorption isotherm and unknown energy distribution function is described by so‐called adsorption integral equation, a linear Fredholm integral equation of the first kind. We consider the case of the Langmuir kernel when the equation can be reduced to the Stieltjes integral equation. A new method for solving the Stieltjes equation is developed. The method is based on the ideas of complex approximation with constraints. The numerical algorithms constructed on the base of this method allow reduction of the problem under consideration to linear or linear‐quadratic programming problems. The method is compared with the usual regularization methods. The obtained results can be useful for the evaluation of the experimental adsorption energy distribution from experimental data. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 191–200, 2000     

A study of the collisional fragmentation problem using the Gamma distribution approximation

The nonlinear fragmentation population balance formulation has been elevated in recent years from a prototype for studying nonlinear integro-differential equations to a vehicle for analyzing and understanding several physicochemical processes of technological interest. The so-called pure collisional fragmentation, which is the particular mode of nonlinear fragmentation induced by collisions between particles, is studied here. It is shown that the corresponding population balance equation admits large time asymptotic (self-similarity) solutions for homogeneous fragmentation and collision functions (kernels). The self-similar solutions are given in closed form for some simple kernels. Based on the shape of the self-similar solutions the method of moments with Gamma distribution approximation is employed for transient solution (from initial state to establishment of the asymptotic shape) of the collisional fragmentation equation. These solutions are presented for several sets of parameters and their behavior is discussed rather extensively. The present study is similar to the one has already been performed for the case of the much simpler linear fragmentation equation [G. Madras, B.J. McCoy, AIChE J. 44 (1998) 647].     

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 reduction method for multiple time scale stochastic reaction networks

In this paper we develop a reduction method for multiple time scale stochastic reaction networks. When the transition-rate matrix between different states of the species is available, we obtain systems of reduced equations, whose solutions can successively approximate, to any degree of accuracy, the exact probability that the reaction system be in any particular state. For the case when the transition-rate matrix is not available, one needs to rely on the chemical master equation. For this case, we obtain a corresponding reduced master equation with first-order accuracy. We illustrate the accuracy and efficiency of both approaches by simulating several motivating examples and comparing the results of our simulations with the results obtained by the exact method. Our examples include both linear and nonlinear reaction networks as well as a three time scale stochastic reaction-diffusion model arising from gene expression.     

Highly efficient and exact method for parallelization of grid‐based algorithms and its implementation in DelPhi

The Gauss–Seidel (GS) method is a standard iterative numerical method widely used to solve a system of equations and, in general, is more efficient comparing to other iterative methods, such as the Jacobi method. However, standard implementation of the GS method restricts its utilization in parallel computing due to its requirement of using updated neighboring values (i.e., in current iteration) as soon as they are available. Here, we report an efficient and exact (not requiring assumptions) method to parallelize iterations and to reduce the computational time as a linear/nearly linear function of the number of processes or computing units. In contrast to other existing solutions, our method does not require any assumptions and is equally applicable for solving linear and nonlinear equations. This approach is implemented in the DelPhi program, which is a finite difference Poisson–Boltzmann equation solver to model electrostatics in molecular biology. This development makes the iterative procedure on obtaining the electrostatic potential distribution in the parallelized DelPhi several folds faster than that in the serial code. Further, we demonstrate the advantages of the new parallelized DelPhi by computing the electrostatic potential and the corresponding energies of large supramolecular structures. © 2012 Wiley Periodicals, Inc.     

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.     

An Exact Propagator for Solving the Triatomic Reactive Schrödinger Equation

The exact short time propagator, in a form similar to the Crank-Nicholson method but in the spirit of spectrally transformed Hamiltonian, was proposed to solve the triatomic reactive time-dependent schrödinger equation. This new propagator is exact and unconditionally convergent for calculating reactive scattering processes with large time step sizes. In order to improve the computational efficiency, the spectral difference method was applied. This resulted the Hamiltonian with elements confined in a narrow diagonal band. In contrast to our previous theoretical work, the discrete variable representation was applied and resulted in full Hamiltonian matrix. As examples, the collision energy-dependent probability of the triatomic H+H₂ and O+O₂ reaction are calculated. The numerical results demonstrate that this new propagator is numerically accurate and capable of propagating the wave packet with large time steps. However, the efficiency and accuracy of this new propagator strongly depend on the mathematical method for solving the involved linear equations and the choice of preconditioner.     

Non-Markovian stochastic Schrödinger equations in different temperature regimes: a study of the spin-boson model

Stochastic Schrodinger equations are used to describe the dynamics of a quantum open system in contact with a large environment, as an alternative to the commonly used master equations. We present a study of the two main types of non-Markovian stochastic Schrodinger equations, linear and nonlinear ones. We compare them both analytically and numerically, the latter for the case of a spin-boson model. We show in this paper that two linear stochastic Schrodinger equations, derived from different perspectives by Diosi, Gisin, and Strunz [Phys. Rev. A 58, 1699 (1998)], and Gaspard and Nagaoka [J. Chem. Phys. 13, 5676 (1999)], respectively, are equivalent in the relevant order of perturbation theory. Nonlinear stochastic Schrodinger equations are in principle more efficient than linear ones, as they determine solutions with a higher weight in the ensemble average which recovers the reduced density matrix of the quantum open system. However, it will be shown in this paper that for the case of a spin-boson system and weak coupling, this improvement does only occur in the case of a bath at high temperature. For low temperatures, the sampling of realizations of the nonlinear equation is practically equivalent to the sampling of the linear ones. We study further this result by analyzing, for both temperature regimes, the driving noise of the linear equations in comparison to that of the nonlinear equations.     

On the equivalence between the microscopic frictional theory and transition-state theory

The generalized Langevin equations are presented by considering such microscopic motions of molecules described by the microscopic Hamiltonian whose potential function is quadratic and internal degree of freedom is multidimensional. Considering the long time behavior of the reactive mode, the Grote-Hynes equation has been derived from the generalized Langevin equations. Furthermore, we have proved that solving the Grote-Hynes equation is equivalent to solving the eigenvalue problem for the whole system, and then the Grote-Hynes treatment coincides with the transition-state theory for the whole system. © 1996 John Wiley & Sons, Inc.     

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     

Rheological properties and particle behaviors of a non-dilute colloidal dispersion composed of ferromagnetic spherocylinder particles subjected to a simple shear flow (analysis by means of mean-field approximation)

We have investigated the rheological properties and the orientational distributions of particles of a non-dilute colloidal dispersion, which is composed of ferromagnetic spherocylinder particles, subject to a simple shear flow. The mean-field approximation is applied to take into account the interactions between spherocylinder particles. The basic equation of the orientational distribution function has been derived from the balance of the torques (including the term due to the mean field approximation) acting on the particle in an applied magnetic field; this is an integrodifferential equation. Then, the governing equation has been solved by means of the method of successive approximation and Galerkin's method. The results obtained here are summarized as follows. For the case of strong magnetic interactions between particles, the particle exhibits a sharp peak of the orientational distribution even for a weak applied magnetic field. In this case, the mean magnetic moment of the particle becomes large, which leads to strong interactions between the applied magnetic field and the particle. Thus, the particle tends to point to the magnetic filed direction under these situations. Also, in this case, a large increase in viscosity is obtained due to such a restriction concerning the particle orientation.     

Adsorption integral equation via complex approximation with constraints: kernel of general form

For the monolayer adsorption on a homogeneous surface, including arbitrary range lateral interactions, the isotherm can be written as a power series of the Langmuir isotherm. If this isotherm is used as the kernel in the adsorption integral equation, this integral equation can be solved in an analytical form. Because the global isotherm is usually known as a set of experimental values, the use of a numerical method is inevitable. A new numerical method for solving the adsorption integral equation with a kernel of general form is developed. It is based on recent results concerning the structure of the local isotherm and on the ideas of complex approximation with constraints, and allows reduction of the problem under consideration to a linear‐quadratic programming problem. Results of numerical experiments are presented. The method can be useful for the evaluation of the adsorption energy distribution from experimental data. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1058–1066, 2001     

Numerical coupled Liouville approach: Application to second hyperpolarizability of molecular aggregate

In our previous study [Int. J. Quant. Chem., to appear], we have developed a novel numerical calculation scheme for a dynamics of quantum network for linear molecular aggregates under intense time‐dependent electric fields. In this approach, each molecule is assumed to be an electric dipole arranged linearly with an angle from the longitudinal axis, and the molecular interactions are taken into account by adding the radiations from these dipoles to the external electric fields. The effects of the radiations from all the dipoles involve the intermolecular distance, the speed of light, retarded polarization, and its first‐ and second‐order time derivatives at the position of each dipole. The quantum dynamics is performed by solving coupled Liouville equations composed of the Liouville equation for each dipole. In the present study, we develop a calculation approach of nonperturbative second hyperpolarizability γ in our novel approach and examine the γ of dimer models composed of two‐state molecules under the one‐photon near resonant intense laser fields. Similar phase transition‐like behavior in the field‐intensity dependence of the γ is observed. We also investigate the second hyperpolarizability spectra in the three‐photon resonant region for dimers composed of three‐state molecules, which mimic the electronic states of allyl cation. Contrary to the one‐photon resonant case, phase transition‐like behavior is not observed in the intensity dependence of γ in the three‐photon resonant region. ©1999 John Wiley & Sons, Inc. Int J Quant Chem 71: 295–306, 1999     

Nonlinear EHD Stability of the Interfacial Waves of Two Superposed Dielectric Fluids

The slow modulation of the interfacial capillary–gravity waves of two superposed dielectric fluids with uniform depths and solid horizontal boundaries, under the influence of a normal electric field and in the absence of surface charges at their interface, is investigated by using the multiple-time scales method. It is found that the complex amplitude of quasi-monochromatic traveling waves can be described by a nonlinear Schrödinger equation in a frame of reference moving with the group velocity. The stability characteristics of a uniform wave train are examined analytically and numerically on the basis of the nonlinear Schrödinger equation, and some limiting cases are recovered. Three cases appear, depending on whether the depth of the lower fluid is equal to, greater than, or less than the depth of the upper fluid. The effect of the normal electric field is determined for the three stability regions of the pure hydrodynamic case. It is found that the normal electric field has a destabilizing influence in the first stability region and a stabilizing effect in the second and third stability regions. Moreover, one new unstable region or two new stable and unstable regions appear, all of which increase when the electric field increases. On the other hand, the complex amplitude of quasi-monochromatic standing waves near the cutoff wavenumber is governed by a similar type of nonlinear Schrödinger equation in which the roles of time and space are interchanged. This equation makes it possible to estimate the nonlinear effect on the cutoff wavenumber.     

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     

Fragment-based quantum mechanical methods for periodic systems with Ewald summation and mean image charge convention for long-range electrostatic interactions

We describe an Ewald-summation method to incorporate long-range electrostatic interactions into fragment-based electronic structure methods for periodic systems. The present method is an extension of the particle-mesh Ewald technique for combined quantum mechanical and molecular mechanical (QM/MM) calculations, and it has been implemented into the explicit polarization (X-Pol) potential to illustrate the computational details. As in the QM/MM-Ewald method, the X-Pol-Ewald approach is a linear-scaling electrostatic method, in which the short-range electrostatic interactions are determined explicitly in real space and the long-range Ewald pair potential is incorporated into the Fock matrix as a correction. To avoid the time-consuming Fock matrix update during the self-consistent field procedure, a mean image charge (MIC) approximation is introduced, in which the running average with a user-chosen correlation time is used to represent the long-range electrostatic correction as an average effect. Test simulations on liquid water show that the present X-Pol-Ewald method takes about 25% more CPU time than the usual X-Pol method using spherical cutoff, whereas the use of the MIC approximation reduces the extra costs for long-range electrostatic interactions by 15%. The present X-Pol-Ewald method provides a general procedure for incorporating long-range electrostatic effects into fragment-based electronic structure methods for treating biomolecular and condensed-phase systems under periodic boundary conditions.