Continuous development of schemes for parallel computing of the electrostatics in biological systems: Implementation in DelPhi

Due to the enormous importance of electrostatics in molecular biology, calculating the electrostatic potential and corresponding energies has become a standard computational approach for the study of biomolecules and nano‐objects immersed in water and salt phase or other media. However, the electrostatics of large macromolecules and macromolecular complexes, including nano‐objects, may not be obtainable via explicit methods and even the standard continuum electrostatics methods may not be applicable due to high computational time and memory requirements. Here, we report further development of the parallelization scheme reported in our previous work (Li, et al., J. Comput. Chem. 2012, 33, 1960) to include parallelization of the molecular surface and energy calculations components of the algorithm. The parallelization scheme utilizes different approaches such as space domain parallelization, algorithmic parallelization, multithreading, and task scheduling, depending on the quantity being calculated. This allows for efficient use of the computing resources of the corresponding computer cluster. The parallelization scheme is implemented in the popular software DelPhi and results in speedup of several folds. As a demonstration of the efficiency and capability of this methodology, the electrostatic potential, and electric field distributions are calculated for the bovine mitochondrial supercomplex illustrating their complex topology, which cannot be obtained by modeling the supercomplex components alone. © 2013 Wiley Periodicals, Inc.     

Statistical investigation of surface bound ions and further development of BION server to include pH and salt dependence

Ions are engaged in multiple biological processes in cells. By binding to the macromolecules or being mobile in the solvent, they maintain the integrity of the structure of macromolecules; participate in their enzymatic activity; or screen electrostatic interactions. While experimental methods are not always able to assign the exact location of ions, computational methods are in demand. Although the majority of computational methods are successful in predicting the position of ions buried inside macromolecules, they are less effective in deciphering positions of surface bound ions. Here, we propose the new BION algorithm ( http://compbio.clemson.edu/bion_server_ph/ ) that predicts the location of the surface bound ions. It is more efficient and accurate compared to the previous version since it uses more advanced clustering algorithm in combination with pairing rules. In addition, the BION webserver allows specifying the pH and the salt concentration in predicting ions positions. © 2015 Wiley Periodicals, Inc.     

Adaptive multilevel finite element solution of the Poisson–Boltzmann equation II. Refinement at solvent‐accessible surfaces in biomolecular systems

We apply the adaptive multilevel finite element techniques (Holst, Baker, and Wang 21 ) to the nonlinear Poisson–Boltzmann equation (PBE) in the context of biomolecules. Fast and accurate numerical solution of the PBE in this setting is usually difficult to accomplish due to presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domains, and rapid (exponential) nonlinearity. However, these adaptive techniques have shown substantial improvement in solution time over conventional uniform‐mesh finite difference methods. One important aspect of the adaptive multilevel finite element method is the robust a posteriori error estimators necessary to drive the adaptive refinement routines. This article discusses the choice of solvent accessibility for a posteriori error estimation of PBE solutions and the implementation of such routines in the “Adaptive Poisson–Boltzmann Solver” (APBS) software package based on the “Manifold Code” (MC) libraries. Results are shown for the application of this method to several biomolecular systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1343–1352, 2000     

Accuracy assessment of the linear Poisson–Boltzmann equation and reparametrization of the OBC generalized Born model for nucleic acids and nucleic acid–protein complexes

The generalized Born model in the Onufriev, Bashford, and Case (Onufriev et al., Proteins: Struct Funct Genet 2004, 55, 383) implementation has emerged as one of the best compromises between accuracy and speed of computation. For simulations of nucleic acids, however, a number of issues should be addressed: (1) the generalized Born model is based on a linear model and the linearization of the reference Poisson–Boltmann equation may be questioned for highly charged systems as nucleic acids; (2) although much attention has been given to potentials, solvation forces could be much less sensitive to linearization than the potentials; and (3) the accuracy of the Onufriev–Bashford–Case (OBC) model for nucleic acids depends on fine tuning of parameters. Here, we show that the linearization of the Poisson Boltzmann equation has mild effects on computed forces, and that with optimal choice of the OBC model parameters, solvation forces, essential for molecular dynamics simulations, agree well with those computed using the reference Poisson–Boltzmann model. © 2015 Wiley Periodicals, Inc.     

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     

Adaptive multilevel finite element solution of the Poisson–Boltzmann equation I. Algorithms and examples

This article is the first of two articles on the adaptive multilevel finite element treatment of the nonlinear Poisson–Boltzmann equation (PBE), a nonlinear eliptic equation arising in biomolecular modeling. Fast and accurate numerical solution of the PBE is usually difficult to accomplish, due to the presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domain, and rapid (exponential) nonlinearity. In this first article, we explain how adaptive multilevel finite element methods can be used to obtain extremely accurate solutions to the PBE with very modest computational resources, and we present some illustrative examples using two well‐known test problems. The PBE is first discretized with piece‐wise linear finite elements over a very coarse simplex triangulation of the domain. The resulting nonlinear algebraic equations are solved with global inexact Newton methods, which we have described in an article appearing previously in this journal. A posteriori error estimates are then computed from this discrete solution, which then drives a simplex subdivision algorithm for performing adaptive mesh refinement. The discretize–solve–estimate–refine procedure is then repeated, until a nearly uniform solution quality is obtained. The sequence of unstructured meshes is used to apply multilevel methods in conjunction with global inexact Newton methods, so that the cost of solving the nonlinear algebraic equations at each step approaches optimal O(N) linear complexity. All of the numerical procedures are implemented in MANIFOLD CODE (MC), a computer program designed and built by the first author over several years at Caltech and UC San Diego. MC is designed to solve a very general class of nonlinear elliptic equations on complicated domains in two and three dimensions. We describe some of the key features of MC, and give a detailed analysis of its performance for two model PBE problems, with comparisons to the alternative methods. It is shown that the best available uniform mesh‐based finite difference or box‐method algorithms, including multilevel methods, require substantially more time to reach a target PBE solution accuracy than the adaptive multilevel methods in MC. In the second article, we develop an error estimator based on geometric solvent accessibility, and present a series of detailed numerical experiments for several complex biomolecules. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1319–1342, 2000     

Theoretical determination of absolute free energy of reduction of plastoquinone-9 in photosystem II and of plastoquinone-n in DMF

We employed static continuum electrostatics and multi-conformation continuum electrostatics (MCCE) methods to determine the reduction potential () of PQ-9 in a section of Photosystem II (PSII). Both methods relied on the finite difference Poisson–Boltzmann (FDPB) solution. The static method brings out a value (0.01 V) that is close to the experimental one (0.05 V), thereby demonstrating that the surrounding environment critically decides the net free energy change. The value obtained from MCCE (0.04 V) is even closer to the observed value, thereby indicating the importance of protein side-chain and proton motions in the electron transfer process. Furthermore, density functional theory-dielectric polarisable continuum model (DFT-DPCM) was employed to calculate the absolute free energy of reduction of plastoquinone-n (PQ-n, where n is the number of isoprenoid units) in N,N dimethyl formamide (DMF) solvent. The DFT-DPCM method produced reduction potential values of −0.59 and −0.65 V for PQ-1 and PQ-9, respectively. These are more or less in agreement with the experimentally reported values of −0.64 and −0.62 V, respectively.     

Modeling the electrostatic potential of asymmetric lipopolysaccharide membranes: The MEMPOT algorithm implemented in DelPhi

Four chemotypes of the rough lipopolysaccharides (LPS) membrane from Pseudomonas aeruginosa were investigated by a combined approach of explicit water molecular dynamics (MD) simulations and Poisson–Boltzmann continuum electrostatics with the goal to deliver the distribution of the electrostatic potential across the membrane. For the purpose of this investigation, a new tool for modeling the electrostatic potential profile along the axis normal to the membrane, MEMbrane POTential (MEMPOT), was developed and implemented in DelPhi. Applying MEMPOT on the snapshots obtained by MD simulations, two observations were made: (a) the average electrostatic potential has a complex profile but is mostly positive inside the membrane due to the presence of Ca²⁺ ions, which overcompensate for the negative potential created by lipid phosphate groups; and (b) correct modeling of the electrostatic potential profile across the membrane requires taking into account the water phase, while neglecting it (vacuum calculations) results in dramatic changes including a reversal of the sign of the potential inside the membrane. Furthermore, using DelPhi to assign different dielectric constants for different regions of the LPS membranes, it was investigated whether a single frame structure before MD simulations with appropriate dielectric constants for the lipid tails, inner, and the external leaflet regions, can deliver the same average electrostatic potential distribution as obtained from the MD‐generated ensemble of structures. Indeed, this can be attained by using smaller dielectric constant for the tail and inner leaflet regions (mostly hydrophobic) than for the external leaflet region (hydrophilic) and the optimal dielectric constant values are chemotype‐specific. © 2014 Wiley Periodicals, Inc.     

Features of CPB: A Poisson–Boltzmann solver that uses an adaptive cartesian grid

The capabilities of an adaptive Cartesian grid (ACG)‐based Poisson–Boltzmann (PB) solver (CPB) are demonstrated. CPB solves various PB equations with an ACG, built from a hierarchical octree decomposition of the computational domain. This procedure decreases the number of points required, thereby reducing computational demands. Inside the molecule, CPB solves for the reaction‐field component (?_rf) of the electrostatic potential (?), eliminating the charge‐induced singularities in ?. CPB can also use a least‐squares reconstruction method to improve estimates of ? at the molecular surface. All surfaces, which include solvent excluded, Gaussians, and others, are created analytically, eliminating errors associated with triangulated surfaces. These features allow CPB to produce detailed surface maps of ? and compute polar solvation and binding free energies for large biomolecular assemblies, such as ribosomes and viruses, with reduced computational demands compared to other Poisson–Boltzmann equation solvers. The reader is referred to http://www.continuum‐dynamics.com/solution‐mm.html for how to obtain the CPB software. © 2014 Wiley Periodicals, 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     

Conformational sampling with Poisson–Boltzmann forces and a stochastic dynamics/Monte Carlo method: Application to alanine dipeptide

We apply a combination of stochastic dynamics and Monte Carlo methods (MC/SD) to alanine dipeptide, with solvation forces derived from a Poisson–Boltzmann model supplemented with apolar terms. Our purpose is to study the effects of the model parameters, such as the friction constant and the size of the electrostatic finite difference grid, on the rate of conformational sampling and on the accuracy of the resulting free energy map. For dialanine, a converged Ramachandran map is produced in significantly less time than what is required by stochastic dynamics or Monte Carlo alone. MC/SD is also shown to be faster, per timestep, than explicit methods. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1750–1759, 1997     

A new set of atomic radii for accurate estimation of solvation free energy by Poisson–Boltzmann solvent model

The Poisson–Boltzmann implicit solvent (PB) is widely used to estimate the solvation free energies of biomolecules in molecular simulations. An optimized set of atomic radii (PB radii) is an important parameter for PB calculations, which determines the distribution of dielectric constants around the solute. We here present new PB radii for the AMBER protein force field to accurately reproduce the solvation free energies obtained from explicit solvent simulations. The presented PB radii were optimized using results from explicit solvent simulations of the large systems. In addition, we discriminated PB radii for N‐ and C‐terminal residues from those for nonterminal residues. The performances using our PB radii showed high accuracy for the estimation of solvation free energies at the level of the molecular fragment. The obtained PB radii are effective for the detailed analysis of the solvation effects of biomolecules. © 2014 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

pKA in proteins solving the Poisson–Boltzmann equation with finite elements

Knowledge on pK_A values is an eminent factor to understand the function of proteins in living systems. We present a novel approach demonstrating that the finite element (FE) method of solving the linearized Poisson–Boltzmann equation (lPBE) can successfully be used to compute pK_A values in proteins with high accuracy as a possible replacement to finite difference (FD) method. For this purpose, we implemented the software molecular Finite Element Solver (mFES) in the framework of the Karlsberg+ program to compute pK_A values. This work focuses on a comparison between pK_A computations obtained with the well‐established FD method and with the new developed FE method mFES, solving the lPBE using protein crystal structures without conformational changes. Accurate and coarse model systems are set up with mFES using a similar number of unknowns compared with the FD method. Our FE method delivers results for computations of pK_A values and interaction energies of titratable groups, which are comparable in accuracy. We introduce different thermodynamic cycles to evaluate pK_A values and we show for the FE method how different parameters influence the accuracy of computed pK_A values. © 2015 Wiley Periodicals, Inc.     

Highly efficient implementation of pseudospectral time‐dependent density‐functional theory for the calculation of excitation energies of large molecules

We have developed and implemented pseudospectral time‐dependent density‐functional theory (TDDFT) in the quantum mechanics package Jaguar to calculate restricted singlet and restricted triplet, as well as unrestricted excitation energies with either full linear response (FLR) or the Tamm–Dancoff approximation (TDA) with the pseudospectral length scales, pseudospectral atomic corrections, and pseudospectral multigrid strategy included in the implementations to improve the chemical accuracy and to speed the pseudospectral calculations. The calculations based on pseudospectral time‐dependent density‐functional theory with full linear response (PS‐FLR‐TDDFT) and within the Tamm–Dancoff approximation (PS‐TDA‐TDDFT) for G2 set molecules using B3LYP/6‐31G*^* show mean and maximum absolute deviations of 0.0015 eV and 0.0081 eV, 0.0007 eV and 0.0064 eV, 0.0004 eV and 0.0022 eV for restricted singlet excitation energies, restricted triplet excitation energies, and unrestricted excitation energies, respectively; compared with the results calculated from the conventional spectral method. The application of PS‐FLR‐TDDFT to OLED molecules and organic dyes, as well as the comparisons for results calculated from PS‐FLR‐TDDFT and best estimations demonstrate that the accuracy of both PS‐FLR‐TDDFT and PS‐TDA‐TDDFT. Calculations for a set of medium‐sized molecules, including C_n fullerenes and nanotubes, using the B3LYP functional and 6‐31G^** basis set show PS‐TDA‐TDDFT provides 19‐ to 34‐fold speedups for C_n fullerenes with 450–1470 basis functions, 11‐ to 32‐fold speedups for nanotubes with 660–3180 basis functions, and 9‐ to 16‐fold speedups for organic molecules with 540–1340 basis functions compared to fully analytic calculations without sacrificing chemical accuracy. The calculations on a set of larger molecules, including the antibiotic drug Ramoplanin, the 46‐residue crambin protein, fullerenes up to C₅₄₀ and nanotubes up to 14×(6,6), using the B3LYP functional and 6‐31G^** basis set with up to 8100 basis functions show that PS‐FLR‐TDDFT CPU time scales as N^2.05 with the number of basis functions. © 2016 Wiley Periodicals, Inc.     

Comparative assessment of computational methods for the determination of solvation free energies in alcohol‐based molecules

The determination of differences in solvation free energies between related drug molecules remains an important challenge in computational drug optimization, when fast and accurate calculation of differences in binding free energy are required. In this study, we have evaluated the performance of five commonly used polarized continuum model (PCM) methodologies in the determination of solvation free energies for 53 typical alcohol and alkane small molecules. In addition, the performance of these PCM methods, of a thermodynamic integration (TI) protocol and of the Poisson–Boltzmann (PB) and generalized Born (GB) methods, were tested in the determination of solvation free energies changes for 28 common alkane‐alcohol transformations, by the substitution of an hydrogen atom for a hydroxyl substituent. The results show that the solvation model D (SMD) performs better among the PCM‐based approaches in estimating solvation free energies for alcohol molecules, and solvation free energy changes for alkane‐alcohol transformations, with an average error below 1 kcal/mol for both quantities. However, for the determination of solvation free energy changes on alkane‐alcohol transformation, PB and TI yielded better results. TI was particularly accurate in the treatment of hydroxyl groups additions to aromatic rings (0.53 kcal/mol), a common transformation when optimizing drug‐binding in computer‐aided drug design. © 2013 Wiley Periodicals, Inc.     

Accurate and efficient determination of higher roots in diagonalization of large matrices based in function restricted optimization algorithms

The problem of large‐scale matrix diagonalization is analyzed in the context of normal function optimization techniques with particular emphasis on the problem of obtaining high roots. New methods based on function restricted optimization algorithms are presented. The efficiency of these methods is illustrated for lowest and higher and degenerate roots of selected matrices. The diagonalization process is commonly carried out in a subspace, and involves a sort of optimization process, and the dimension of this subspace increases at each iteration. In addition, the success of a diagonalization method in obtaining a desired root strongly depends on the particular optimization procedure chosen. In this work, a rational function optimization procedure is presented that permits obtaining the lowest and higher eigenpairs in an efficient way. Update Hessian matrices formulae, routinely used in normal function optimization problems, are explored in the framework of diagonalization techniques. Finally, a diagonalization method with a fixed subspace dimension during the iterative process is presented. Some examples focused in lowest, higher and degenerate eigenpairs are discussed. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1375–1386, 2000     

Quantum supercharger library: Hyper‐parallelism of the Hartree–Fock method

We present here a set of algorithms that completely rewrites the Hartree–Fock (HF) computations common to many legacy electronic structure packages (such as GAMESS‐US, GAMESS‐UK, and NWChem) into a massively parallel compute scheme that takes advantage of hardware accelerators such as Graphical Processing Units (GPUs). The HF compute algorithm is core to a library of routines that we name the Quantum Supercharger Library (QSL). We briefly evaluate the QSL's performance and report that it accelerates a HF 6‐31G Self‐Consistent Field (SCF) computation by up to 20 times for medium sized molecules (such as a buckyball) when compared with mature Central Processing Unit algorithms available in the legacy codes in regular use by researchers. It achieves this acceleration by massive parallelization of the one‐ and two‐electron integrals and optimization of the SCF and Direct Inversion in the Iterative Subspace routines through the use of GPU linear algebra libraries. © 2015 Wiley Periodicals, Inc.     

On the achievement of high fidelity and scalability for large‐scale diagonalizations in grid‐based DFT simulations

Recent advance in high performance computing (HPC) resources has opened the possibility to expand the scope of density functional theory (DFT) simulations toward large and complex molecular systems. This work proposes a numerically robust method that enables scalable diagonalizations of large DFT Hamiltonian matrices, particularly with thousands of computing CPUs (cores) that are usual these days in terms of sizes of HPC resources. The well‐known Lanczos method is extensively refactorized to overcome its weakness for evaluation of multiple degenerate eigenpairs that is the substance of DFT simulations, where a multilevel parallelization is adopted for scalable simulations in as many cores as possible. With solid benchmark tests for realistic molecular systems, the fidelity of our method are validated against the locally optimal block preconditioned conjugated gradient (LOBPCG) method that is widely used to simulate electronic structures. Our method may waste computing resources for simulations of molecules whose degeneracy cannot be reasonably estimated. But, compared to LOBPCG method, it is fairly excellent in perspectives of both speed and scalability, and particularly has remarkably less (< 10%) sensitivity of performance to the random nature of initial basis vectors. As a promising candidate for solving electronic structures of highly degenerate systems, the proposed method can make a meaningful contribution to migrating DFT simulations toward extremely large computing environments that normally have more than several tens of thousands of computing cores.     

Highly efficient polymer‐stabilized palladium heterogeneous catalyst: Synthesis,characterization and application for Suzuki–Miyaura and Mizoroki–Heck coupling reactions

A suitable approach to stabilize palladium nanoparticles (Pd NPs), with an average diameter of 3–4 nm, on magnetic polymer is described. A new magnetic polymer containing 4′‐(4‐hydroxyphenyl)‐2,2′:6′,2″‐terpyridine (HPTPy) ligand was prepared by the polymerization of itaconic acid (ITC) as a monomer and trimethylolpropane triacrylate (TMPTA) as a cross‐linker and fully characterized. Pd NPs embedded on the magnetic polymer were successfully applied in Suzuki–Miyaura and Mizoroki–Heck coupling reactions under low palladium loading conditions, and provided the corresponding products with excellent yields (up to 98%) and high catalytic activities (TOF up to 257 hr^?1). Also, the catalyst can be easily separated and reused for at least consecutive five times with a small drop in catalytic activity.     

Electrochemical Fabrication and Properties of Highly Ordered Fe‐Doped TiO2 Nanotubes

Highly‐ordered Fe‐doped TiO₂ nanotubes (TiO₂nts) were fabricated by anodization of co‐sputtered Ti–Fe thin films in a glycerol electrolyte containing NH₄F. The as‐sputtered Ti–Fe thin films correspond to a solid solution of Ti and Fe according to X‐ray diffraction. The Fe‐doped TiO₂nts were studied in terms of composition, morphology and structure. The characterization included scanning electron microscopy, energy‐dispersive X‐ray spectroscopy, X‐ray diffraction, UV/Vis spectroscopy, X‐ray photoelectron spectroscopy and Mott–Schottky analysis. As a result of the Fe doping, an indirect bandgap of 3.0 eV was estimated using Tauc’s plot, and this substantial red‐shift extends its photoresponse to visible light. From the Mott–Schottky analysis, the flat‐band potential (E_fb) and the charge carrier concentration (N_D) were determined to be ?0.95 V vs Ag/AgCl and 5.0 ×10¹⁹ cm^?3 respectively for the Fe‐doped TiO₂nts, whilst for the undoped TiO₂nts, E_fb of ?0.85 V vs Ag/AgCl and N_D of 6.5×10¹⁹ cm^?3 were obtained.