A charge-scaling implementation of the variational electrostatic projection method

Two new charge-scaling methods for efficient modeling of the solvated macromolecular environment in hybrid QM/MM calculations of biological reactions are presented. The methods are extensions of the variational electrostatic projection (VEP) method, and allows a subset of atomic charges in the external environment to be adjusted to mimic, in the active dynamical region, the electrostatic potential and field due to the large surrounding macromolecule and solvent. The method has the advantages that it offers improved accuracy, does not require the use of a three-dimensional grid or auxiliary set of fitting points, and requires only minor molecular simulation code modifications. The VEP-cs and VEP-RVM+cs methods are able to attain very high accuracy (relative force errors of 10(-7) or better with appropriate choice of control parameters), and take advantage of a recently introduced set of high-order discretization schemes and Gaussian exponents for boundary element solvation and VEP methods. The methods developed here serve as potentially powerful tools in the arsenal of computational techniques used in multiscale computational modeling problems.     

Solution of the linearized Poisson-Boltzmann equation

Improved methods are formulated for solution of the linearized Poisson-Boltzmann equation, to be used in conjunction with electronic structure calculation on a solute together with dielectric continuum representation of the salt-containing solvent. Volume polarization effects due to quantum mechanical penetration of solute charge density outside the cavity that excludes solvent are treated by exact and by approximate methods analogous to those previously developed for the salt-free case. With boundary element approaches, exact solutions lead to coupled equations for a pair of cavity surface distributions that mimic the polarization of the solvent dielectric and the ionic atmosphere. A novel means is found to effectively decouple these equations, yielding more efficient practical methods for their numerical solution. Detailed comparisons are given to related boundary element formulations previously reported in the literature, which neglect volume polarization, and analogous decoupling is also found for the pair of surface distributions invoked there. Illustrative results are provided for a simple spherical example.     

Mixed homopolymer brushes grafted onto a nanosphere

Microphase separation of mixed A∕B polymer brushes grafted onto a nanosphere with its radius comparable to the size of polymers is investigated by numerical implementation of the self-consistent field theory. The idea is to embed the sphere within a larger cubic computational cell and use a "masking" technique to treat the spherical boundary. The partial differential equations for the chain propagator on the sphere can thus be readily solved with an efficient and high-order accurate pseudospectral method involving fast Fourier transform on a cubic cell. This numerical technique can circumvent the "pole problem" due to the use of a spherical coordinate system in conventional finite difference or finite element grid. We systematically investigate the effect of the total grafting density, composition, chain length asymmetry between two grafted homopolymers as well as spherical radius, i.e., substrate curvature on the formation of island structure with specific arrangement in a regular lattice. A series of island structures with different island numbers representing specific structure symmetry ranging from 2 to 12 except for 11 are found, in contrast to conventional hexagonal arrangement for polymer brushes on a planar substrate. Among these parameters, the spherical radius plays a significant role in determining the type of island structures, i.e., the morphology formed on the sphere.     

Binding energies of on-axis hydrogenic and nonhydrogenic donors in GaAs/Ga1−x Al x As

In a previous work [Phys. Rev. B43(1991)9262], the binding energies of hydrogenic and nonhydrogenic on-axis donors in GaAs/Ga_1?x A1 _X As quantum-well wires of circular cross section have been calculated as functions of the radius of the quantum-well wire. In both the hydrogenic and nonhydrogenic cases, a variational trial wave function was chosen that could be written as the product of an “envelope function”, and a function containing the variational parameter. It was assumed in these calculations that the potential barrier that exists at the surface of the GaAs cylinder and the surrounding Ga_1?x A1 _X As matrix is infinite. For the envelope function, an ordinary Bessel function of the first kind and of order zero was chosen. This envelope function satisfies the boundary condition, the vanishing of the trial function at the surface of the quantum-well wire. The question arises: how sensitive are the calculated binding energies to the choice of the envelope function? In the present work, we attempt to provide a partial answer to this question by choosing another envelope function, a spherical Bessel function of order zero. This function also satisfies the boundary condition and makes the trial wave function vanish at the interface between the GaAs cylinder and the Ga_1?x A1 _X As matrix. Our calculations show that the binding energies of both the hydrogenic and the nonhydrogenic on-axis donors depend on the choice of the envelope function.     

Calculation of the Maxwell stress tensor and the Poisson-Boltzmann force on a solvated molecular surface using hypersingular boundary integrals

The electrostatic interaction among molecules solvated in ionic solution is governed by the Poisson-Boltzmann equation (PBE). Here the hypersingular integral technique is used in a boundary element method (BEM) for the three-dimensional (3D) linear PBE to calculate the Maxwell stress tensor on the solvated molecular surface, and then the PB forces and torques can be obtained from the stress tensor. Compared with the variational method (also in a BEM frame) that we proposed recently, this method provides an even more efficient way to calculate the full intermolecular electrostatic interaction force, especially for macromolecular systems. Thus, it may be more suitable for the application of Brownian dynamics methods to study the dynamics of protein/protein docking as well as the assembly of large 3D architectures involving many diffusing subunits. The method has been tested on two simple cases to demonstrate its reliability and efficiency, and also compared with our previous variational method used in BEM.     

Computation of electrostatic forces between solvated molecules determined by the Poisson-Boltzmann equation using a boundary element method

A rigorous approach is proposed to calculate the electrostatic forces among an arbitrary number of solvated molecules in ionic solution determined by the linearized Poisson-Boltzmann equation. The variational principle is used and implemented in the frame of a boundary element method (BEM). This approach does not require the calculation of the Maxwell stress tensor on the molecular surface, therefore it totally avoids the hypersingularity problem in the direct BEM whenever one needs to calculate the gradient of the surface potential or the stress tensor. This method provides an accurate and efficient way to calculate the full intermolecular electrostatic interaction energy and force, which could potentially be used in Brownian dynamics simulation of biomolecular association. The method has been tested on some simple cases to demonstrate its reliability and efficiency, and parts of the results are compared with analytical results and with those obtained by some known methods such as adaptive Poisson-Boltzmann solver.     

A boundary element method for molecular electrostatics with electrolyte effects

A boundary element method is developed to compute the electrostatic potential inside and around molecules in an electrolyte solution. A set of boundary integral equations are derived based on the integral formulations of the Poisson equation and the linearized Poisson-Boltzmann equation. The boundary integral equations are then solved numerically after discretizing the molecular surface into a number of flat triangular elements. The method is applied to a spherical molecule for which analytical solutions are available. Use is made of both constant and linearly varying unknowns over the boundary elements, and the method is tested for various values of parameters such as the dielectric constant of the molecule, ionic strength, and the location of the interior point charge. The use of the boundary integral method incorporating the nonlinear Poisson-Boltzmann equation is also briefly discussed.     

Analyzing polymer crazing as quasifracture

This paper deals with a viscoelastic boundary element method for analyzing a polymer quasifracture usually called a craze in polymers. A time-dependent boundary stiffness is considered on the quasifracture envelope surface. The viscoelastic property of the glassy polymer is represented by a generalized Kelvin model with multiple retardation times. According to the linear viscoelastic correspondence principle, the associated elasticity solution can be solved by applying the general integral boundary element method. Then the viscoelastic solution in the time domain can be obtained by applying a collocation Laplace inversion transformation. Using these methods, the quasifracture problem composed of an isolated craze opening with time-dependent stiffness traction in a stressed rectangular plate is analyzed. The displacement profile and the stress distribution around the craze envelope surface are computed.     

Properties of nearly one-electron molecules. I. An iterative Green function approach to calculating the reaction matrix

An ab initio R-matrix method for determining the molecular reaction matrix of scattering theory is introduced. The method makes use of a principal-value Green function to compute the collision channel wave functions for the scattered electron, in combination with the Kohn variational scheme for the evaluation of R-matrix eigenvalues on a spherical boundary surface at short range. This technique permits the size of the bounded volume in the variational calculation to be reduced, making the computations fast and efficient. The reaction matrix is determined in a form that minimizes its energy dependence. Thus the procedure does not require modification or an increase in the computational effort to study the electronic structure and dynamics in Rydberg molecules with extremely polar ion cores. The analysis is specialized to examine the bound-state and free-electron scattering properties of nearly one-electron molecular systems, which are characterized by a Rydberg/scattering electron incident on a closed-shell ion core. However, it is shown that the treatment is compatible with all-electron/ab initio representations of open-shell and nonlinear polyatomic ion cores, emphasizing its generality. The introduced approach is used to calculate the electronic spectrum of the calcium monofluoride molecule, which has the extremely polar (Ca+2F-)+e- closed-shell ion-core configuration. The calculation utilizes an effective single-electron potential determined by M. Arif, C. Jungen, and A. L. Roche [J. Chem. Phys. 106, 4102 (1997)] previously. Close agreement with experimental data is obtained. The results demonstrate the practical utility of this method as a viable alternative to the standard variational approaches.     

Variational electrostatic projection (VEP) methods for efficient modeling of the macromolecular electrostatic and solvation environment in activated dynamics simulations

New methods for the calculation of electrostatic interactions between the active dynamical region and surrounding external solvated macromolecular environment in hybrid quantum mechanical/molecular mechanical (QM/MM) simulations are presented. The variational electrostatic projection (VEP) method, and related variational reverse-mapping procedure (VEP-RVM) utilize an expansion in Gaussian surface elements for representation of electrostatic interactions. The use of a discretized surface that surrounds the active dynamical region greatly reduces the number of interactions with the particles of the external environment. The methods are tested on two catalytic RNA systems: the hammerhead and the hairpin ribozymes. It is shown that with surface elements numbering from 302 to 1202 points the direct VEP and VEP-RVM methods are able to obtain relative force errors in the range of 0.5-0.05% and 0.09-0.0001%, respectively, using a 4.0 A projection buffer. These results are encouraging and provide an essential step in the development of new variational macromolecular solvent-boundary methods for QM/MM calculations of enzyme reactions.     

A note on atomic density

In atomic systems, electron density has a simple finite expansion in spherical harmonics times radial factors. The difficulties in the calculation of some radial factors are illustrated in the low‐lying states of the carbon atom. Single‐particle methods such as Hartree–Fock and approximate density functional theory cannot ensure the correct expansion of the density in spherical harmonics. Wave‐function methods are appropriate but, as some expansion terms are entirely due to correlation, these methods only will give correct results for high‐quality variational functions. Using full‐configuration integration (CI), all the terms predicted by the theory appear and are not negligible but the convergence of the term due to correlation toward its correct value is uncertain even for very large CI spaces. © 2012 Wiley Periodicals, Inc.     

Solution of the nonlinear Poisson-Boltzmann equation using pseudo-transient continuation and the finite element method

The nonlinear Poisson-Boltzmann (PB) equation is solved using Newton-Krylov iterations coupled with pseudo-transient continuation. The PB potential is used to compute the electrostatic energy and evaluate the force on a user-specified contour. The PB solver is embedded in a existing, 3D, massively parallel, unstructured-grid, finite element code. Either Dirichlet or mixed boundary conditions are allowed. The latter specifies surface charges, approximates far-field conditions, or linearizes conditions "regulating" the surface charge. Stability and robustness are proved using results for backward Euler differencing of diffusion equations. Potentials and energies of charged spheres and plates are computed and results compared to analysis. An approximation to the potential of the nonlinear, spherical charge is derived by combining two analytic formulae. The potential and force due to a conical probe interacting with a flat plate are computed for two types of boundary conditions: constant potential and constant charge. The second case is compared with direct force measurements by chemical force microscopy. The problem is highly nonlinear-surface potentials of the linear and nonlinear PB equations differ by over an order of magnitude. Comparison of the simulated and experimentally measured forces shows that approximately half of the surface carboxylic acid groups, of density 1/(0.2 nm2), ionize in the electrolyte implying surface charges of 0.4 C/m2, surface potentials of 0.27 V, and a force of 0.6 nN when the probe and plate are 8.7 nm apart.     

The fast multipole boundary element method for molecular electrostatics: An optimal approach for large systems

We propose a fast implementation of the boundary element method for solving the Poisson equation, which approximately determines the electrostatic field around solvated molecules of arbitrary shape. The method presented uses computational resources of order O(N) only, where N is the number of elements representing the dielectric boundary at the molecular surface. The method is based on the Fast Multipole Algorithm by Rokhlin and Greengard, which is used to calculate the Coulombic interaction between surface elements in linear time. We calculate the solvation energies of a sphere, a small polar molecule, and a moderately sized protein. The values obtained by the boundary element method agree well with results from finite difference calculations and show a higher degree of consistency due to the absence of grid dependencies. The boundary element method can be taken to a much higher accuracy than is possible with finite difference methods and can therefore be used to verify their validity. © 1995 by John Wiley & Sons, Inc.     

Surface-induced morphologies of ABC star triblock copolymer in spherical cavities

The morphologies and phase diagrams exhibited by symmetric ABC star triblock copolymer nanoparticles are investigated on the basis of real-space self-consistent field theory. The ABC star triblock copolymers were chosen to be tiling-forming with fixed polymer parameter and the spherical boundaries were modeled using the masking technique. We first study a number of examples where the ABC triblock copolymers confined in spherical cavities with neutral surface. Then, two types of spherical cavity distinct preferential surfaces are considered, including both A-block attractive and repulsive preferential surfaces. We aim at the effects due to various spherical cavity diameters and the degree of interactions between the polymer and the spherical surface. A variety of morphologies, such as ring-like structures, concentric sphere, and irregular cylinder, were identified in phase diagrams. The results show that both the degree of interactions and spherical diameters can influence the formation of morphologies so that ring-like structures and other novel structures could be obtained.     

Comparison of solvent reaction field representations

 Alternative ways are examined for representing a reaction field to treat the important effects of long-range electrostatic interaction with a solvent in electronic structure calculations on the properties of a solute. Several extant boundary element methods for approximate representation of the solvent reaction field in terms of surface charge distributions are considered, and analogous new methods for approximate representation in terms of surface dipole distributions are introduced. Illustrative computational results are presented on representative small neutral and ionic solutes to evaluate the relative accuracy of various methods. Received: 2 July 2001 / Accepted: 10 September 2001 / Published online: 19 December 2001     

New hybrid method for reactive systems from integrating molecular orbital or molecular mechanics methods with analytical potential energy surfaces

A computational approach to calculating potential energy surfaces for reactive systems is presented and tested. This hybrid approach is based on integrated methods where calculations for a small model system are performed by using analytical potential energy surfaces, and for the real system by using molecular orbital or molecular mechanics methods. The method is tested on a hydrogen abstraction reaction by using the variational transition-state theory with multidimensional tunneling corrections. The agreement between the calculated and experimental information depends on the quality of the method chosen for the real system. When the real system is treated by accurate quantum mechanics methods, the rate constants are in excellent agreement with the experimental measurements over a wide temperature range. When the real system is treated by molecular mechanics methods, the results are still good, which is very encouraging since molecular mechanics itself is not at all capable of describing this reactive system. Since no experimental information or additional fits are required to apply this method, it can be used to improve the accuracy of molecular orbital methods or to extend the molecular mechanics method to treat any reactive system with the single constraint of the availability of an analytical potential energy surface that describes the model system.     

Numerical comparison between Maxwell stress method and equivalent multipole approach for calculation of the dielectrophoretic force in single-cell traps

This paper presents detailed numerical calculations of the dielectrophoretic force in traps designed for single-cell trapping. A trap with eight planar electrodes is studied for spherical and ellipsoidal particles using the boundary element method (BEM). Multipolar approximations of orders one to three are compared with the full Maxwell stress tensor (MST) calculation of the electrical force on spherical particles. Ellipsoidal particles are also studied, but in their case only the dipolar approximation is available for comparison with the MST solution. The results show that a small number of multipolar terms need to be considered in order to obtain accurate results for spheres, even in the proximity of the electrodes, and that the full MST calculation is only required in the study of non-spherical particles.     

Double layer interactions between charge-regulated colloidal surfaces: pair potentials for spherical particles bearing ionogenic surface groups

The variational approach of Reiner and Radke (1990) is employed to investigate the effect of surface charge regulation upon the double layer interaction free energy V_e of pairs of colloidal particles immersed in an electrolyte. A model for dissociating surface groups that permits the consideration of an arbitrary number ofion-complexation reactions is introduced. The variational method is then used to derive (in the Poisson-Boltzmann approximation) the configurational free energy functional Ω of an ensemble of particles bearing such groups. The Debye-Hückel (DH) linearization process is applied to this functional, and ensuing consistency issues are examined.The DH free energy is extremized for a configuration of two interacting flat plates, and Derjaguin's (1934 and 1939) method is used to obtain an approximate analytical form for V_e for two different-sized spherical particles bearing different surface groups. This second problem is next considered from the perspective ofLevine's (1934, 1939b) exact multipole expansion of the electrostatic potential surrounding two axisymmetric particles. It is shown that the linear superposition approximation (LSA) for V_e developed by Levine (1939c) and Verwey and Overbeek (1948) emerges rigorously from this formulation in the limit of large interparticle separations. The interaction free energy from Levine's expansion is calculated te a six digit accuracy for identical spheres over the range of regulated behavior from fixed surface charge density q_s to fixed surface potential ψ_s for surface-surface separation h to Debye length λ ratios from 0 to 2 and ratios of the particle radius a to λ of 0.1, 1, and 10. These results are compared to those obtained from Derjaguin's method and the linear superposition approximation. Derjaguin's method is only quantitatively accurate (in error by less than 10%) for the largest value of a/λ and becomes progressively less so as the boundary is changed from perfectly regulating (constant ψ_s) to unregulated (constant q_s). Agreement of the LSA with the exact V_e is good over a wide range of parameters, but worsens for large a/λ and small h/λ. Appendices present extensions of our approach to surfaces bearing more than one type of complexing group and to the consideration of Stern layer formation at the particle-electrolyte boundary in the context of a standard model for metal oxide-aqueous interfaces.     

Finite-element implementation for electron transport in nanostructures

We have modeled transport properties of nanostructures using Green's-function method within the framework of the density-functional theory. The scheme is computationally demanding, so numerical methods have to be chosen carefully. A typical solution to the numerical burden is to use a special basis-function set, which is tailored to the problem in question, for example, the atomic-orbital basis. In this paper we present our solution to the problem. We have used the finite-element method with a hierarchical high-order polynomial basis, the so-called p elements. This method allows the discretation error to be controlled in a systematic way. The p elements work so efficiently that they can be used to solve interesting nanosystems described by nonlocal pseudopotentials. We demonstrate the potential of the implementation with two different systems. As a test system a simple Na-atom chain between two leads is modeled and the results are compared with several previous calculations. Secondly, we consider a thin hafnium dioxide (HfO2) layer on a silicon surface as a model for a gate structure of the next generation of microelectronics.     

Encapsulation of a polyelectrolyte chain by an oppositely charged spherical surface

Using the ground state dominance approximation and a variational theory, we study the encapsulation of a polyelectrolyte chain by an oppositely charged spherical surface. The electrostatic attraction between the polyelectrolyte and the surface and the entropy loss of the encapsulated polyelectrolyte chain dictate the optimum conditions for encapsulation. Two scenarios of encapsulation are identified: entropy-dominated and adsorption-dominated encapsulation. In the entropy-dominated encapsulation regime, the polyelectrolyte chain is delocalized, and the optimum radius of the encapsulating sphere decreases with increasing the attraction. In the adsorption-dominated encapsulation regime, the polyelectrolyte chain is strongly localized near the surface, and the optimum radius increases with increasing the attraction. After identifying a universal encapsulation parameter, the dependencies of the optimum radius on the salt concentration, surface charge density, polymer charge density, and polymer length are explored.