The conformational sensitivity of iterative stockholder partitioning schemes

An information-theoretic approach to basis-set fitting of electron densities and other non-negative functions

The numerical ill-conditioning associated with approximating an electron density with a convex sum of Gaussian or Slater-type functions is overcome by using the (extended) Kullback-Leibler divergence to measure the deviation between the target and approximate density. The optimized densities are non-negative and normalized, and they are accurate enough to be used in applications related to molecular similarity, the topology of the electron density, and numerical molecular integration. This robust, efficient, and general approach can be used to fit any non-negative normalized functions (e.g., the kinetic energy density and molecular electron density) to a convex sum of non-negative basis functions. We present a fixed-point iteration method for optimizing the Kullback-Leibler divergence and compare it to conventional gradient-based optimization methods. These algorithms are released through the free and open-source BFit package, which also includes a L2-norm squared optimization routine applicable to any square-integrable scalar function.

Constrained iterative Hirshfeld charges: A variational approach

We develop a variational procedure for the iterative Hirshfeld (HI) partitioning scheme. The main practical advantage of having a variational framework is that it provides a formal and straightforward approach for imposing constraints (e.g., fixed charges on certain atoms or molecular fragments) when computing HI atoms and their properties. Unlike many other variants of the Hirshfeld partitioning scheme, HI charges do not arise naturally from the information-theoretic framework, but only as a reverse-engineered construction of the objective function. However, the procedure we use is quite general and could be applied to other problems as well. We also prove that there is always at least one solution to the HI equations, but we could not prove that its self-consistent equations would always converge for any given initial pro-atom charges. Our numerical assessment of the constrained iterative Hirshfeld method shows that it satisfies many desirable traits of atoms in molecules and has the potential to surpass existing approaches for adding constraints when computing atomic properties.

Multi-center decomposition of molecular densities: a mathematical perspective

The aim of this article is to analyze from a mathematical perspective some existing schemes to partition a molecular density into several atomic contributions, with a specific focus on Iterative Stockholder Atom (ISA) methods. We provide a unified mathematical framework to describe the latter family of methods and propose a new scheme, named L-ISA (for linear approximation of ISA), which generalizes the so-called Additive Variational Hirshfeld method. We prove several important mathematical properties of the ISA and L-ISA minimization problems and show that the so-called ISA algorithms can be viewed as alternating minimization schemes, which in turn enables us to obtain new convergence results for these numerical methods. Specific mathematical properties of the ISA decomposition for diatomic systems are also presented. Different schemes are numerically compared on different molecules and we discuss the advantages and drawbacks of each approach. Numerical results on diatomic systems illustrate the proven mathematical properties.

Using atomic charges to model molecular polarization

We review different models for introducing electric polarization in force fields, with special focus on methods where polarization is modelled at the atomic charge level. While electric polarization has been included in several force fields, the common approach has been to focus on atomic dipole polarizability. Several approaches allow modelling electric polarization by using charge-flow between charge sites instead, but this has been less exploited, despite that atomic charges and charge-flow is expected to be more important than atomic dipoles and dipole polarizability. A number of challenges are required to be solved for charge-flow models to be incorporated into polarizable force fields, for example how to parameterize the models and how to make them computational efficient.

Protein C-GeM: A Coarse-Grained Electron Model for Fast and Accurate Protein Electrostatics Prediction

The electrostatic potential (ESP) is a powerful property for understanding and predicting electrostatic charge distributions that drive interactions between molecules. In this study, we compare various charge partitioning schemes including fitted charges, density-based quantum mechanical (QM) partitioning schemes, charge equilibration methods, and our recently introduced coarse-grained electron model, C-GeM, to describe the ESP for protein systems. When benchmarked against high quality density functional theory calculations of the ESP for tripeptides and the crambin protein, we find that the C-GeM model is of comparable accuracy to ab initio charge partitioning methods, but with orders of magnitude improvement in computational efficiency since it does not require either the electron density or the electrostatic potential as input.

Analysis and visualization of energy densities. II. Insights from linear-response time-dependent density functional theory calculations

Inspired by the analysis of Kohn-Sham energy densities by Nakai and coworkers, we extended the energy density analysis to linear-response time-dependent density functional theory (LR-TDDFT) calculations. Using ethylene-tetrafluoroethylene and oxyluciferin-water complexes as examples, distinctive distribution patterns were demonstrated for the excitation energy densities of local excitations (within a molecular fragment) and charge-transfer excitations (between molecular fragments). It also provided a simple way to compute the effective energy of both hot carriers (particle and hole) from charge-transfer excitations via an integration of the excitation energy density over the donor and acceptor grid points.

First-Principles Many-Body Nonadditive Polarization Energies from Monomer and Dimer Calculations Only: A Case Study on Water

The many-body polarization energy is the major source of nonadditivity in strongly polar systems such as water. This nonadditivity is often considerable and must be included, if only in an average manner, to correctly describe the physical properties of the system. Models for the polarization energy are usually parametrized using experimental data, or theoretical estimates of the many-body effects. Here we show how many-body polarization models can be developed for water complexes using data for the monomer and dimer only using ideas recently developed in the field of intermolecular perturbation theory and state-of-the-art approaches for calculating distributed molecular properties based on the iterated stockholder atoms (ISA) algorithm. We show how these models can be calculated, and we validate their accuracy in describing the many-body nonadditive energies of a range of water clusters. We further investigate their sensitivity to the details of the polarization damping models used. We show how our very best polarization models yield many-body energies that agree with those computed with coupled-cluster methods, but at a fraction of the computational cost.

Describing Molecular Polarizability by a Bond Capacity Model

We propose a bond capacity model for describing molecular polarization in force field energy functions at the charge-only level. Atomic charges are calculated by allowing charge to flow between atom pairs according to a bond capacity and a difference in electrostatic potential. The bond capacity is closely related to the bond order and decays to zero as the bond distance is increased. The electrostatic potential is composed of an intrinsic potential, identified as the electronegativity, and a screened Coulomb potential from all other charges. The bond capacity model leads to integer fragment charges upon bond dissociation and displays linear scaling of the polarizability with system size. Bond capacity parameters can be derived from reference molecular polarizabilities, while electronegativity parameters can be derived from reference atomic charges or a reference molecular electrostatic potential. Out-of-plane polarization for planar systems is modeled by off-nuclei charge sites. The model is shown to be able to reproduce anisotropic reference molecular polarizabilities with an accuracy of ∼10% using a limited set of bond capacity parameters and can describe both inter- and intramolecular polarization.

Hydration Free Energies in the FreeSolv Database Calculated with Polarized Iterative Hirshfeld Charges

Computer simulations of biomolecular systems often use force fields, which are combinations of simple empirical atom-based functions to describe the molecular interactions. Even though polarizable force fields give a more detailed description of intermolecular interactions, nonpolarizable force fields, developed several decades ago, are often still preferred because of their reduced computation cost. Electrostatic interactions play a major role in biomolecular systems and are therein described by atomic point charges. In this work, we address the performance of different atomic charges to reproduce experimental hydration free energies in the FreeSolv database in combination with the GAFF force field. Atomic charges were calculated by two atoms-in-molecules approaches, Hirshfeld-I and Minimal Basis Iterative Stockholder (MBIS). To account for polarization effects, the charges were derived from the solute's electron density computed with an implicit solvent model, and the energy required to polarize the solute was added to the free energy cycle. The calculated hydration free energies were analyzed with an error model, revealing systematic errors associated with specific functional groups or chemical elements. The best agreement with the experimental data is observed for the AM1-BCC and the MBIS atomic charge methods. The latter includes the solvent polarization and presents a root-mean-square error of 2.0 kcal mol-1 for the 613 organic molecules studied. The largest deviation was observed for phosphorus-containing molecules and the molecules with amide, ester and amine functional groups.

Introducing DDEC6 atomic population analysis: part 3. Comprehensive method to compute bond orders

A new method to compute accurate bond orders for metallic, covalent, polar-covalent, ionic, multi-centered, aromatic, dative, dispersion, and hydrogen bonding.

Atom and Bond Fukui Functions and Matrices: A Hirshfeld-I Atoms-in-Molecule Approach

The Fukui function is often used in its atom-condensed form by isolating it from the molecular Fukui function using a chosen weight function for the atom in the molecule. Recently, Fukui functions and matrices for both atoms and bonds separately were introduced for semiempirical and ab initio levels of theory using Hückel and Mulliken atoms-in-molecule models. In this work, a double partitioning method of the Fukui matrix is proposed within the Hirshfeld-I atoms-in-molecule framework. Diagonalizing the resulting atomic and bond matrices gives eigenvalues and eigenvectors (Fukui orbitals) describing the reactivity of atoms and bonds. The Fukui function is the diagonal element of the Fukui matrix and may be resolved in atom and bond contributions. The extra information contained in the atom and bond resolution of the Fukui matrices and functions is highlighted. The effect of the choice of weight function arising from the Hirshfeld-I approach to obtain atom- and bond-condensed Fukui functions is studied. A comparison of the results with those generated by using the Mulliken atoms-in-molecule approach shows low correlation between the two partitioning schemes.

Introducing DDEC6 atomic population analysis: part 1. Charge partitioning theory and methodology

We introduce a new atomic population analysis method that performs exceptionally well across an extremely broad range of periodic and non-periodic material types.

Introducing DDEC6 atomic population analysis: part 2. Computed results for a wide range of periodic and nonperiodic materials

DDEC6 atomic population analysis gives excellent performance for small and large molecules, porous solids, dense solids, solid surfaces, organometallic complexes, nanoclusters, and magnetic materials.

Chemical Assignment of Symmetry-Adapted Perturbation Theory Interaction Energy Components: The Functional-Group SAPT Partition

Recently, we introduced an effective atom-pairwise partition of the many-body symmetry-adapted perturbation theory (SAPT) interaction energy decomposition, producing a method known as atomic SAPT (A-SAPT) [Parrish, R. M.; Sherrill, C. D. J. Chem. Phys. 2014, 141, 044115]. A-SAPT provides ab initio atom-pair potentials for force field development and also automatic visualizations of the spatial contributions of noncovalent interactions, but often has difficulty producing chemically useful partitions of the electrostatic energy, due to the buildup of oscillating partial charges on adjacent functional groups. In this work, we substitute chemical functional groups in place of atoms as the relevant local quasiparticles in the partition, resulting in a functional-group-pairwise partition denoted as functional-group SAPT (F-SAPT). F-SAPT assigns integral sets of local occupied electronic orbitals and protons to chemical functional groups and linking σ bonds. Link-bond contributions can be further assigned to chemical functional groups to simplify the analysis. This approach yields a SAPT partition between pairs of functional groups with integral charge (usually neutral), preventing oscillations in the electrostatic partition. F-SAPT qualitatively matches chemical intuition and the cut-and-cap fragmentation technique but additionally yields the quantitative many-body SAPT interaction energy. The conceptual simplicity, chemical utility, and computational efficiency of F-SAPT is demonstrated in the context of phenol dimer, proflavine(+)-DNA intercalation, and a cucurbituril host-guest inclusion complex.

Spatial assignment of symmetry adapted perturbation theory interaction energy components: The atomic SAPT partition

We develop a physically-motivated assignment of symmetry adapted perturbation theory for intermolecular interactions (SAPT) into atom-pairwise contributions (the A-SAPT partition). The basic precept of A-SAPT is that the many-body interaction energy components are computed normally under the formalism of SAPT, following which a spatially-localized two-body quasiparticle interaction is extracted from the many-body interaction terms. For electrostatics and induction source terms, the relevant quasiparticles are atoms, which are obtained in this work through the iterative stockholder analysis (ISA) procedure. For the exchange, induction response, and dispersion terms, the relevant quasiparticles are local occupied orbitals, which are obtained in this work through the Pipek-Mezey procedure. The local orbital atomic charges obtained from ISA additionally allow the terms involving local orbitals to be assigned in an atom-pairwise manner. Further summation over the atoms of one or the other monomer allows for a chemically intuitive visualization of the contribution of each atom and interaction component to the overall noncovalent interaction strength. Herein, we present the intuitive development and mathematical form for A-SAPT applied in the SAPT0 approximation (the A-SAPT0 partition). We also provide an efficient series of algorithms for the computation of the A-SAPT0 partition with essentially the same computational cost as the corresponding SAPT0 decomposition. We probe the sensitivity of the A-SAPT0 partition to the ISA grid and convergence parameter, orbital localization metric, and induction coupling treatment, and recommend a set of practical choices which closes the definition of the A-SAPT0 partition. We demonstrate the utility and computational tractability of the A-SAPT0 partition in the context of side-on cation-π interactions and the intercalation of DNA by proflavine. A-SAPT0 clearly shows the key processes in these complicated noncovalent interactions, in systems with up to 220 atoms and 2845 basis functions.

Numerical study on the partitioning of the molecular polarizability into fluctuating charge and induced atomic dipole contributions

In order to carry out a detailed analysis of the molecular static polarizability, which is the response of the molecule to a uniform external electric field, the molecular polarizability was computed using the finite-difference method for 21 small molecules, using density functional theory. Within nine charge population schemes (Löwdin, Mulliken, Becke, Hirshfeld, CM5, Hirshfeld-I, NPA, CHELPG, MK-ESP) in common use, the charge fluctuation contribution is found to dominate the molecular polarizability, with its ratio ranging from 59.9% with the Hirshfeld or CM5 scheme to 96.2% with the Mulliken scheme. The Hirshfeld-I scheme is also used to compute the other contribution to the molecular polarizability coming from the induced atomic dipoles, and the atomic polarizabilities in eight small molecules and water pentamer are found to be highly anisotropic for most atoms. Overall, the results suggest that (a) more emphasis probably should be placed on the charge fluctuation terms in future polarizable force field development and (b) an anisotropic polarizability might be more suitable than an isotropic one in polarizable force fields based entirely or partially on the induced atomic dipoles.

FOHI-D: an iterative Hirshfeld procedure including atomic dipoles

In this work, a new partitioning method based on the FOHI method (fractional occupation Hirshfeld-I method) will be discussed. The new FOHI-D method uses an iterative scheme in which both the atomic charge and atomic dipole are calculated self-consistently. In order to induce the dipole moment on the atom, an electric field is applied during the atomic SCF calculations. Based on two sets of molecules, the atomic charge and intrinsic atomic dipole moment of hydrogen and chlorine atoms are compared using the iterative Hirshfeld (HI) method, the iterative Stockholder atoms (ISA) method, the FOHI method, and the FOHI-D method. The results obtained are further analyzed as a function of the group electronegativity of Boyd et al. [J. Am. Chem. Soc. 110, 4182 (1988); Boyd et al., J. Am. Chem. Soc. 114, 1652 (1992)] and De Proft et al. [J. Phys. Chem. 97, 1826 (1993)]. The molecular electrostatic potential (ESP) based on the HI, ISA, FOHI, and FOHI-D charges is compared with the ab initio ESP. Finally, the effect of adding HI, ISA, FOHI, and FOHI-D atomic dipoles to the multipole expansion as a function of the precision of the ESP is analyzed.