# Journal of Chemical Theory and Computation (JCTC)

## When is the Fukui Function Not Normalized? The Danger of Inconsistent Energy Interpolation Models in Density Functional Theory

### Abstract

When one defines the energy of a molecule with a noninteger number of electrons by interpolation of the energy values for integer-charged states, the interpolated electron density, Fukui function, and higher-order derivatives of the density are generally not normalized correctly. The necessary and sufficient condition for consistent energy interpolation models is that the corresponding interpolated electron density is correctly normalized to the number of electrons. A necessary, but not sufficient, condition for correct normalization is that the energy interpolant be a linear function of the reference energies. Consistent with this general rule, polynomial interpolation models and, in particular, the quadratic E vs N model popularized by Parr and Pearson, do give normalized densities and density derivatives. Interestingly, an interpolation model based on the square root of the electron number also satisfies the normalization constraints. We also derive consistent least-norm interpolation models. In contrast to these models, the popular rational and exponential forms for E vs N do not give normalized electron densities and density derivatives.

## Minimal Basis Iterative Stockholder: Atoms-in-Molecules for Force-Field Development

### Abstract

Atomic partial charges appear in the Coulomb term of many force-field models and can be derived from electronic structure calculations with a myriad of atoms-in-molecules (AIM) methods. More advanced models have also been proposed, using the distributed nature of the electron cloud and atomic multipoles. In this work, an electrostatic force field is defined through a concise approximation of the electron density, for which the Coulomb interaction is trivially evaluated. This approximate "pro-density" is expanded in a minimal basis of atom-centered s-type Slater density functions, whose parameters are optimized by minimizing the Kullback-Leibler divergence of the pro-density from a reference electron density, e.g. obtained from an electronic structure calculation. The proposed method, Minimal Basis Iterative Stockholder (MBIS), is a variant of the Hirshfeld AIM method but it can also be used as a density-fitting technique. An iterative algorithm to refine the pro-density is easily implemented with a linear-scaling computational cost, enabling applications to supramolecular systems. The benefits of the MBIS method are demonstrated with systematic applications to molecular databases and extended models of condensed phases. A comparison to 14 other AIM methods shows its effectiveness when modeling electrostatic interactions. MBIS is also suitable for rescaling atomic polarizabilities in the Tkatchenko-Sheffler scheme for dispersion interactions.

## eReaxFF: A Pseudo-Classical Treatment of Explicit Electrons within Reactive Force Field Simulations

### Abstract

We present a computational tool, eReaxFF, for simulating explicit electrons within the framework of the standard ReaxFF reactive force field method. We treat electrons explicitly in a pseudoclassical manner that enables simulation several orders of magnitude faster than quantum chemistry (QC) methods, while retaining the ReaxFF transferability. We delineate here the fundamental concepts of the eReaxFF method and the integration of the Atom-condensed Kohn–Sham DFT approximated to second order (ACKS2) charge calculation scheme into the eReaxFF. We trained our force field to capture electron affinities (EA) of various species. As a proof-of-principle, we performed a set of molecular dynamics (MD) simulations with an explicit electron model for representative hydrocarbon radicals. We establish a good qualitative agreement of EAs of various species with experimental data, and MD simulations with eReaxFF agree well with the corresponding Ehrenfest dynamics simulations. The standard ReaxFF parameters available in the literature are transferrable to the eReaxFF method. The computationally economic eReaxFF method will be a useful tool for studying large-scale chemical and physical systems with explicit electrons as an alternative to computationally demanding QC methods.

## A comparison of barostats for the mechanical characterization of metal-organic frameworks

### Abstract

In this paper, three barostat coupling schemes for pressure control, which are commonly used in molecular dynamics simulations, are critically compared to characterise the rigid MOF-5 and the flexible MIL-53(Al) metal-organic frameworks. We investigate the performance of the three barostats, the Berendsen, the Martyna-Tuckerman-Tobias-Klein (MTTK) and the Langevin coupling methods, in reproducing the cell parameters and the pressure versus volume behaviour in isothermal-isobaric simulations. A thermodynamic integration method is used to construct the free energy profiles as a function of volume at finite temperature. It is observed that the aforementioned static properties are well reproduced with the three barostats. However, for static properties depending nonlinearly on the pressure, the Berendsen barostat might give deviating results as it suppresses pressure fluctuations more drastically. Finally, dynamic properties, which are directly related to the fluctuations of the cell, such as the time to transition from the large-pore to the closed-pore phase, cannot be well reproduced by any of the coupling schemes.

## Variational optimization of the second order density matrix corresponding to a seniority-zero configuration interaction wave function

### Abstract

We perform a direct variational determination of the second-order (two-particle) density matrix corresponding to a many-electron system, under a restricted set of the two-index $N$-representability $\mathcal{P}$-, $\mathcal{Q}$-, and $\mathcal{G}$-conditions. In addition, we impose a set of necessary constraints that the two-particle density matrix must be derivable from a doubly-occupied many-electron wave function, i.e.\ a singlet wave function for which the Slater determinant decomposition only contains determinants in which spatial orbitals are doubly occupied. We rederive the two-index $N$-representability conditions first found by Weinhold and Wilson and apply them to various benchmark systems (linear hydrogen chains, He, $\text{N}_2$ and $\text{CN}^-$). This work is motivated by the fact that a doubly-occupied many-electron wave function captures in many cases the bulk of the static correlation. Compared to the general case, the structure of doubly-occupied two-particle density matrices causes the associate semidefinite program to have a very favorable scaling as $L^3$, where $L$ is the number of spatial orbitals. Since the doubly-occupied Hilbert space depends on the choice of the orbitals, variational calculation steps of the two-particle density matrix are interspersed with orbital-optimization steps (based on Jacobi rotations in the space of the spatial orbitals). We also point to the importance of symmetry breaking of the orbitals when performing calculations in a doubly-occupied framework.

## Non-Variational Orbital Optimization Techniques for the AP1roG Wave Function

### Abstract

We introduce new nonvariational orbital optimization schemes for the antisymmetric product of one-reference orbital geminal (AP1roG) wave function (also known as pair-coupled cluster doubles) that are extensions to our recently proposed projected seniority-two (PS2-AP1roG) orbital optimization method [ J. Chem. Phys. 2014, 140, 214114)]. These approaches represent less stringent approximations to the PS2-AP1roG ansatz and prove to be more robust approximations to the variational orbital optimization scheme than PS2-AP1roG. The performance of the proposed orbital optimization techniques is illustrated for a number of well-known multireference problems: the insertion of Be into H2, the automerization process of cyclobutadiene, the stability of the monocyclic form of pyridyne, and the aromatic stability of benzene.

## Multicenter Bonding in Ditetracyanoethylene Dianion: A Simple Aromatic Picture in Terms of Three-Electron Bonds

### Abstract

The nature of the multicenter, long bond in ditetracyanoethylene dianion complex [TCNE]22– is elucidated using high level ab initio Valence Bond (VB) theory coupled with Quantum Monte Carlo (QMC) methods. This dimer is the prototype of the general family of pancake-bonded dimers with large interplanar separations. Quantitative results obtained with a compact wave function in terms of only six VB structures match the reference CCSD(T) bonding energies. Analysis of the VB wave function shows that the weights of the VB structures are not compatible with a covalent bond between the π* orbitals of the fragments. On the other hand, these weights are consistent with a simple picture in terms of two resonating bonding schemes, one displaying a pair of interfragment three-electron σ bonds and the other displaying intrafragment three-electron π bonds. This simple picture explains at once (1) the long interfragment bond length, which is independent of the countercations but typical of three-electron (3-e) CC σ bonds, (2) the interfragment orbital overlaps which are very close to the theoretical optimal overlap of 1/6 for a 3-e σ bond, and (3) the unusual importance of dynamic correlation, which is precisely the main bonding component of 3-e bonds. Moreover, it is shown that the [TCNE]22– system is topologically equivalent to the square C4H42– dianion, a well-established aromatic system. To better understand the role of the cyano substituents, the unsubstituted diethylenic Na+2[C2H4]22– complex is studied and shown to be only metastable and topologically equivalent to a rectangular C4H42– dianion, devoid of aromaticity.

## A New Mean-Field Method Suitable for Strongly Correlated Electrons: Computationally Facile Antisymmetric Products of Nonorthogonal Geminals

### Abstract

We propose an approach to the electronic structure problem based on noninteracting electron pairs that has similar computational cost to conventional methods based on noninteracting electrons. In stark contrast to other approaches, the wave function is an antisymmetric product of nonorthogonal geminals, but the geminals are structured so the projected Schrödinger equation can be solved very efficiently. We focus on an approach where, in each geminal, only one of the orbitals in a reference Slater determinant is occupied. The resulting method gives good results for atoms and small molecules. It also performs well for a prototypical example of strongly correlated electronic systems, the hydrogen atom chain.

## Hirshfeld-E partitioning: AIM charges with an improved trade-off between robustness and accurate electrostatics

### Abstract

For the development of ab-initio derived force fields, atomic charges must be computed from electronic structure computations, such that (i) they accurately describe the molecular electrostatic potential (ESP) and (ii) they are transferable to the force-field application of interest. The Iterative Hirshfeld (Hirshfeld-I or HI) scheme meets both requirements for organic molecules. For inorganic oxide clusters, however, Hirshfeld-I becomes ambiguous because electron densities of nonexistent isolated anions are needed as input. Herein, we propose a simple Extended Hirshfeld (Hirshfeld-E or HE) scheme to overcome this limitation. The performance of the new HE scheme is compared to four popular atoms-in-molecules schemes, using two tests involving a set of 248 silica clusters. These tests show that the new HE scheme provides an improved trade-off between the ESP accuracy and the transferability of the charges. The new scheme is a generalization of the Hirshfeld-I scheme and it is expected that its improvements are to a large extent applicable to molecular systems containing elements from the entire periodic table.