Calculation of electrostatic interaction energies in molecular dimers from atomic multipole moments obtained by different methods of electron density partitioning

Accurate and fast evaluation of electrostatic interactions in molecular systems is still one of the most challenging tasks in the rapidly advancing field of macromolecular chemistry, including molecular recognition, protein modeling and drug design. One of the most convenient and accurate approaches is based on a Buckingham-type approximation that uses the multipole moment expansion of molecular/atomic charge distributions. In the mid-1980s it was shown that the pseudoatom model commonly used in experimental X-ray charge density studies can be easily combined with the Buckingham-type approach for calculation of electrostatic interactions, plus atom-atom potentials for evaluation of the total interaction energies in molecular systems. While many such studies have been reported, little attention has been paid to the accuracy of evaluation of the purely electrostatic interactions as errors may be absorbed in the semiempirical atom-atom potentials that have to be used to account for exchange repulsion and dispersion forces. This study is aimed at the evaluation of the accuracy of the calculation of electrostatic interaction energies with the Buckingham approach. To eliminate experimental uncertainties, the atomic moments are based on theoretical single-molecule electron densities calculated at various levels of theory. The electrostatic interaction energies for a total of 11 dimers of alpha-glycine, N-acetylglycine and L-(+)-lactic acid structures calculated according to Buckingham with pseudoatom, stockholder and atoms-in-molecules moments are compared with those evaluated with the Morokuma-Ziegler energy decomposition scheme. For alpha-glycine a comparison with direct "pixel-by-pixel" integration method, recently developed Gavezzotti, is also made. It is found that the theoretical pseudoatom moments combined with the Buckingham model do predict the correct relative electrostatic interactions energies, although the absolute interaction energies are underestimated in some cases. The good agreement between electrostatic interaction energies computed with Morokuma-Ziegler partitioning, Gavezzotti's method, and the Buckingham approach with atoms-in-molecules moments demonstrates that reliable and accurate evaluation of electrostatic interactions in molecular systems of considerable complexity is now feasible.     

Toward accurate transferable electrostatic models for polypeptides: A distributed multipole study of blocked amino acid residue charge distributions

We have developed a method for building up accurate electrostatic models for polypeptides, based on a distributed multipole representation of the SCF charge densities for the dipeptides (CH₃.CO.NH.CHR.CO.NH.CH₃) of the naturally occurring amino acids. It is based on the observation that each peptide residue has almost the correct formal charge (0, ±1). We find that the variations in the backbone charge distributions (excluding proline) with sidechain have a negligible effect on the predicted electrostatic potential around the residue. However changes in the atomic multipoles with the torsion angles (?, ψ, χ) are more significant, and may need to be taken into account if electrostatic potential close to the residue is required to high accuracy. This type of DMA peptide library provides more accurate, more theoretically based, estimates of the electrostatic potential around polypeptides than current models.     

The Atomic Partial Charges Arboretum: Trying to See the Forest for the Trees

Atomic partial charges are among the most commonly used interpretive tools in quantum chemistry. Dozens of different ‘population analyses’ are in use, which are best seen as proxies (indirect gauges) rather than measurements of a ‘general ionicity’. For the GMTKN55 benchmark of nearly 2,500 main-group molecules, which span a broad swathe of chemical space, some two dozen different charge distributions were evaluated at the PBE0 level near the 1-particle basis set limit. The correlation matrix between the different charge distributions exhibits a block structure; blocking is, broadly speaking, by charge distribution class. A principal component analysis on the entire dataset suggests that nearly all variation can be accounted for by just two ‘principal components of ionicity’: one has all the distributions going in sync, while the second corresponds mainly to Bader QTAIM vs. all others. A weaker third component corresponds to electrostatic charge models in opposition to the orbital-based ones. The single charge distributions that have the greatest statistical similarity to the first principal component are iterated Hirshfeld (Hirshfeld-I) and a minimal-basis projected modification of Bickelhaupt charges. If three individual variables, rather than three principal components, are to be identified that contain most of the information in the whole dataset, one representative for each of the three classes of Corminboeuf et al. is needed: one based on partitioning of the density (such as QTAIM), a second based on orbital partitioning (such as NPA), and a third based on the molecular electrostatic potential (such as HLY or CHELPG).     

Long-range electrostatic interactions in hybrid quantum and molecular mechanical dynamics using a lattice summation approach

A method based on a lattice summation technique for treating long-range electrostatic interactions in hybrid quantum mechanics/molecular mechanics simulations is presented in this article. The quantum subsystem is studied at the semiempirical level, whereas the solvent is described by a two-body potential of molecular mechanics. Molecular dynamics simulations of a (quantum) chloride ion in (classical) water have been performed to test this technique. It is observed that the application of the lattice summations to solvent-solvent interactions as well as on solute-solvent ones has a significant effect on solvation energy and diffusion coefficient. Moreover, two schemes for the computation of the long-range contribution to the electrostatic interaction energy are investigated. The first one replaces the exact charge distribution of the quantum solute by a Mulliken charge distribution. The long-range electrostatic interactions are then calculated for this charge distribution that interacts with the solvent molecule charges. The second one is more accurate and involves a modified Fock operator containing long-range electron-charge interactions. It is shown here that both schemes lead to similar results, the method using Mulliken charges for the evaluation of long-range interactions being, however, much more computationally efficient.     

Gaussian induced dipole polarization model

A new induced dipole polarization model based on interacting Gaussian charge densities is presented. In contrast to the original induced point dipole model, the Gaussian polarization model is capable of finite interactions at short distances. Aspects of convergence related to the Gaussian model will be explored. The Gaussian polarization model is compared with the damped Thole-induced dipole model and the point dipole model. It will be shown that the Gaussian polarization model performs slightly better than the Thole model in terms of fitting to molecular polarizability tensors. An advantage of the model based on Gaussian charge distribution is that it can be easily generalized to other multipole moments and provide effective damping for both permanent electrostatic and polarization models. Finally, a method of parameterizing polarizabilities is presented. This method is based on probing a molecule with point charges and fitting polarizabilities to electrostatic potential. In contrast to the generic atom type polarizabilities fit to molecular polarizability tensors, probed polarizabilities are significantly more accurate in terms of reproducing molecular polarizability tensors and electrostatic potential, while retaining conformational transferability.     

Electrostatic Energies in the 1,4‐Dichlorobenzene Polymorph Crystals: the Role of Charge Density Overlap Effects in Crystal Packing Analysis

Electrostatic and polarization energies for the three known polymorphic crystal structures of 1,4‐dichlorobenzene, as well as for one particularly stable virtual crystal structure generated by computer search, were calculated by a new accurate numerical integration method over static molecular charge densities obtained from high level ab initio molecular‐orbital calculations. Results are compared with those from standard empirical atom‐atom force fields. The new electrostatic energies, which include charge density overlap (penetration) effects, are seen to be much larger than and sometimes of opposite sign to those derived from point‐charge models. None of the four polymorphs is substantially more stable than the others on electrostatic‐energy grounds. Molecule‐molecule electrostatic energies have been calculated for the more important molecular pairs in each of the four structures; trends are found to be very different from those indicated by point‐charge energies or by total energies estimated with a parametric atom‐atom force field. Conclusions based exclusively on analysis of intermolecular atom contacts and point‐charge electrostatics may need to be modified in the light of the new kind of calculation.     

Spatial partitioning of the molecular wave function: Reexamination of the bond‐charge model of covalent binding

Owing to the complexity of many‐body quantum interactions, theoretical problems in chemistry almost always necessitate numerical methods of solution. In turn, and somewhat paradoxically, this reliance on computation often precipitates a need for approximate models that can provide a conceptual framework for the interpretation of the accurate numerical solutions. In this report, we reexamine and improve upon the simple bond‐charge (SBC) model of Parr and Borkman with the express purpose of examining its role as a suitable intuitive model for bonding energetics. We first discuss the role of spatial partitioning in the SBC model and the molecular‐structure approximations that are implicit within the model. We then develop an improved version of the model by including electrostatic repulsion between shared electron pairs; we refer to this enhanced version as a semiempirical charge‐cloud (SCC) model. Applying the SCC model to a broad range of diatomics, we find that it is able to predict dissociation energies on a par with simple variational calculations. As an illustrative application of the SCC model, we use it to estimate unknown dissociation energies and rotational constants for a number of alkali dimers. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

Advancing beyond charge analysis using the electronic localization function: Chemically intuitive distribution of electrostatic moments

We propose here an evaluation of chemically intuitive distributed electrostatic moments using the topological analysis of the electron localization function (ELF). As this partition of the total charge density provides an accurate representation of the molecular dipole, the distributed electrostatic moments based on the ELF partition (DEMEP) allows computing of local moments located at non atomic centers such as lone pairs, sigma bonds and pi systems. As the local dipole contribution can be decomposed in polarization and charge transfer components, our results indicate that local dipolar polarization of the lone pairs and chemical reactivity are closely related whereas the charge transfer contribution is the key factor driving the local bond dipole. Results on relevant molecules show that local dipole contributions can be used to rationalize inductive polarization effects in alcohols derivatives and typical hydrogen bond interactions. Moreover, bond quadrupole polarization moments being related to a pi character enable to discuss bond multiplicities, and to sort families of molecules according to their bond order. That way, the nature of the C-O bond has been revisited for several typical systems by means of the DEMEP analysis which appears also helpful to discuss aromaticity. Special attention has been given to the carbon monoxide molecule, to the CuCO complex and to a weak intramolecular N|-CO interaction involved in several biological systems. In this latter case, it is confirmed that the bond formation is mainly linked to the CO bond polarization. Transferability tests show that the approach is suitable for the design of advanced force fields.     

A divide and conquer real-space approach for all-electron molecular electrostatic potentials and interaction energies

A computational scheme to perform accurate numerical calculations of electrostatic potentials and interaction energies for molecular systems has been developed and implemented. Molecular electron and energy densities are divided into overlapping atom-centered atomic contributions and a three-dimensional molecular remainder. The steep nuclear cusps are included in the atom-centered functions making the three-dimensional remainder smooth enough to be accurately represented with a tractable amount of grid points. The one-dimensional radial functions of the atom-centered contributions as well as the three-dimensional remainder are expanded using finite element functions. The electrostatic potential is calculated by integrating the Coulomb potential for each separate density contribution, using our tensorial finite element method for the three-dimensional remainder. We also provide algorithms to compute accurate electron-electron and electron-nuclear interactions numerically using the proposed partitioning. The methods have been tested on all-electron densities of 18 reasonable large molecules containing elements up to Zn. The accuracy of the calculated Coulomb interaction energies is in the range of 10(-3) to 10(-6) E(h) when using an equidistant grid with a step length of 0.05 a(0).     

Accurate modeling of the intramolecular electrostatic energy of proteins

The (?, ψ) energy surface of blocked alanine (N-acetyl–N′-methyl alanineamide) was calculated at the Hartree-Fock (HF)/6-31G* level using ab initio molecular orbital theory. A collection of six electrostatic models was constructed, and the term electrostatic model was used to refer to (1) a set of atomic charge densities, each unable to deform with conformation; and (2) a rule for estimating the electrostatic interaction energy between a pair of atomic charge densities. In addition to two partial charge and three multipole electrostatic models, this collection includes one extremely detailed model, which we refer to as nonspherical CPK. For each of these six electrostatic models, parameters—in the form of partial charges, atomic multipoles, or generalized atomic densities—were calculated from the HF/6-31G* wave functions whose energies define the ab initio energy surface. This calculation of parameters was complicated by a problem that was found to originate from the locking in of a set of atomic charge densities, each of which contains a small polarization-induced deformation from its idealized unpolarized state. It was observed that the collective contribution of these small polarization-induced deformations to electrostatic energy differences between conformations can become large relative to ab initio energy differences between conformations. For each of the six electrostatic models, this contribution was reduced by an averaging of atomic charge densities (or electrostatic energy surfaces) over a large collection of conformations. The ab initio energy surface was used as a target with respect to which relative accuracies were determined for the six electrostatic models. A collection of 42 more complete molecular mechanics models was created by combining each of our six electrostatic models with a collection of seven models of repulsion + dispersion + intrinsic torsional energy, chosen to provide a representative sample of functional forms and parameter sets. A measure of distance was defined between model and ab initio energy surfaces; and distances were calculated for each of our 42 molecular mechanics models. For most of our 12 standard molecular mechanics models, the average error between model and ab initio energy surfaces is greater than 1.5 kcal/mol. This error is decreased by (1) careful treatment of the nonspherical nature of atomic charge densities, and (2) accurate representation of electrostatic interaction energies of types 1—2 and 1—3. This result suggests an electrostatic origin for at least part of the error between standard model and ab initio energy surfaces. Given the range of functional forms that is used by the current generation of protein potential functions, these errors cannot be corrected by compensating for errors in other energy components. © 1995 by John Wiley & Sons, Inc.     

Designing electrostatic interactions in biological systems via charge optimization or combinatorial approaches: insights and challenges with a continuum electrostatic framework

Electrostatic interactions between biological molecules are crucially influenced by their aqueous environment, with efficient and accurate models of solvent effects required for robust molecular design strategies. Continuum electrostatic models provide a reasonable balance between computational efficiency and accurate system representation. In this article, I review two specific molecular design strategies, charge optimization and combinatorial design, paying particular attention to how the continuum framework (also briefly described herein) successfully enables both theoretical insights and molecular designs and presents a challenge in design applications due to what I call “the isostericity constraint.” Efforts to work around the isostericity constraint and other challenges are discussed. Additionally, particular emphasis is placed on using such models in the rational design of particularly tight, specific, or promiscuous interactions, in keeping with the increased sophistication of current molecular design applications.     

The hydrogen bond in ice probed by soft x-ray spectroscopy and density functional theory

We combine photoelectron and x-ray absorption spectroscopy with density functional theory to derive a molecular orbital picture of the hydrogen bond in ice. We find that the hydrogen bond involves donation and back-donation of charge between the oxygen lone pair and the O-H antibonding orbitals on neighboring molecules. Together with internal s-p rehybridization this minimizes the repulsive charge overlap of the connecting oxygen and hydrogen atoms, which is essential for a strong attractive electrostatic interaction. Our joint experimental and theoretical results demonstrate that an electrostatic model based on only charge induction from the surrounding medium fails to properly describe the internal charge redistributions upon hydrogen bonding.     

Force related atomic multipoles in planar molecules. Derivation of atomic quadrupole and octupole moments

Atomic multipoles as defined by current methods generally do not account for forces in molecules that arise from external electrostatic fields. It is pointed out that such forces and the electrostatic potential that the molecule itself generates are both determined by the molecular multipolar tensors. The latter constitute therefore the fundamental molecular constants that determine the molecular electrostatics apart from polarization. In general the multipolar tensors include contributions from the atomic multipoles and their fluxes. In planar molecules, however, the perpendicular charge flux is zero by symmetry. This gives rise to a (previously introduced) formalism that extracts analytical, force-related, atomic multipoles from the molecular multipolar tensors. This formalism is extended in this work to include force-related (FR) atomic quadrupoles and octupoles in planar molecules. The properties of the FR atomic multipoles, including their perpendicular fluxes, are discussed and some formal theoretical and computational advantages that characterize them are indicated. As an example, the electrostatics of OCS, including the molecular electrostatic potential and the forces on the nuclei due to an external point charge, is discussed.     

Electrostatic energy in the effective fragment potential method: theory and application to benzene dimer

Evaluation of the electrostatic energy within the effective fragment potential (EFP) method is presented. The performance of two variants of the distributed multipole analysis (DMA) together with two different models for estimating the charge penetration energies was studied using six homonuclear dimers. The importance of damping the higher order multipole terms, i.e. charge dipole, was also investigated. Damping corrections recover more than 70% of the charge penetration energy in all dimers, whereas higher order damping introduces only minor improvement. Electrostatic energies calculated by the numerical DMA are less accurate than those calculated by the analytic DMA. Analysis of bonding in the benzene dimer shows that EFP with inclusion of the electrostatic damping term performs very well compared to the high-level coupled cluster singles, doubles, and perturbative triples method. The largest error of 0.4 kcal/mol occurs for the sandwich dimer configuration. This error is about half the size of the corresponding error in second order perturbation theory. Thus, EFP in the current implementation is an accurate and computationally inexpensive method for calculating interaction energies in weakly bonded molecular complexes.     

Molecular orientations at interfaces by extended polarizable continuum model

This study presents an investigation into orientation of molecular solutes at the interface of liquid water and other media. The calculation of electrostatic free energy of molecular solute is based on an extension of the polarizable continuum model (PCM) to interfacial system. The extended PCM computational scheme is incorporated with the self‐consistent field procedure which is necessary to obtain more accurate electrostatic free energy and charge density distribution. The computation of non‐electrostatic energy for interfacial system is also realized. Applying the numerical procedure to molecular systems, N,N′‐diethyl‐p‐nitroaniline (DEPNA) at air/water interface and p‐nitrophenol (PNP) at cyclohexane/water interface, the average orientational angles are in reasonable agreement with the experimental results. Taking both the electrostatic and the non‐electrostatic energies into account, the analysis on the energy profiles shows that the electrostatic solvation energy is the dominant factor in determining the orientation angle for PNP, whereas for DEPNA, the orientation angle mainly depends on the cavitation energy. This suggests that, in addition to the electrostatic energy, taking the cavitation energy into account may provide a more complete view when we survey the molecular orientation at interface. © 2010 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

Study of peptide conformation in terms of the ABEEM/MM method

The ABEEM/MM model (atom-bond electronegativity equalization method fused into molecular mechanics) is applied to study of the polypeptide conformations. The Lennard-Jones and torsional parameters were optimized to be consistent with the ABEEM/MM fluctuating charge electrostatic potential. The hydrogen bond was specially treated with an electrostatic fitting function. Molecular dipole moments, dimerization energies, and hydrogen bond lengths of complexes are reasonably achieved by our model, compared to ab initio results. The ABEEM/MM fluctuating charge model reproduces both the peptide conformational energies and structures with satisfactory accuracy with low computer cost. The transferability is tested by applying the parameters of our model to the tetrapeptide of alanine and another four dipeptides. The overall RMS deviations in conformational energies and key dihedral angles for four di- or tetrapeptide, is 0.39 kcal/mol and 7.7 degrees . The current results agree well with those by the accurate ab initio method, and are comparable to those from the best existing force fields. The results make us believe that our fluctuating charge model can obtain more promising results in protein and macromolecular modeling with good accuracy but less computer cost.     

Electrostatic energies and forces computed without explicit interparticle interactions: a linear time complexity formulation

A rapid method for the calculation of the electrostatic energy of a system without a cutoff is described in which the computational time grows linearly with the number of particles or charges. The inverse of the distance is approximated as a polynomial, which is then transformed into a function whose terms involve individual particles, instead of particle pairs, by a partitioning of the double sum. In this way, the electrostatic energy that is determined by the interparticle interactions is obtained without explicit calculation of these interactions. For systems of positive charges positioned on a face-centered cubic lattice, the calculation of the energy by the new method is shown to be faster than the calculation of the exact energy, in many cases by an order of magnitude, and to be accurate to within 1-2%. The application of this method to increase the accuracy of conventional truncation-based calculations in condensed-phase systems is also demonstrated by combining the approximated long-range electrostatic interactions with the exact short-range interactions in a "hybrid" calculation. For a 20-A sphere of water molecules, the forces are shown to be six times as accurate using this hybrid method as those calculated with conventional truncation of the electrostatic energy function at 12 A. This is accomplished with a slight increase in speed, and with a sevenfold increase in speed relative to the exact all-pair calculation. Structures minimized with the hybrid function are shown to be closer to structures minimized with an exact all-pair electrostatic energy function than are those minimized with a conventional 13-A cutoff-based electrostatic energy function. Comparison of the energies and forces calculated with the exact method illustrate that the absolute errors obtained with standard truncation can be very large. The extension of the current method to other pairwise functions as well as to multibody functions, is described.