Recent developments in the PySCF program package

Hirshfeld atom refinement based on projector augmented wave densities with periodic boundary conditions

Hirshfeld atom refinement (HAR) is an X-ray diffraction refinement method that, in numerous publications, has been shown to give H-atom bond lengths in close agreement with neutron diffraction derived values. Presented here is a first evaluation of an approach using densities derived from projector augmented wave (PAW) densities with three-dimensional periodic boundary conditions for HAR. The results show an improvement over refinements that neglect the crystal environment or treat it classically, while being on a par with non-periodic approximations for treating the solid-state environment quantum mechanically. A suite of functionals were evaluated for this purpose, showing that the SCAN and revSCAN functionals are most suited to these types of calculation.

DQC: A Python program package for differentiable quantum chemistry

Automatic differentiation represents a paradigm shift in scientific programming, where evaluating both functions and their derivatives is required for most applications. By removing the need to explicitly derive expressions for gradients, development times can be shortened and calculations can be simplified. For these reasons, automatic differentiation has fueled the rapid growth of a variety of sophisticated machine learning techniques over the past decade, but is now also increasingly showing its value to support ab initio simulations of quantum systems and enhance computational quantum chemistry. Here, we present an open-source differentiable quantum chemistry simulation code and explore applications facilitated by automatic differentiation: (1) calculating molecular perturbation properties, (2) reoptimizing a basis set for hydrocarbons, (3) checking the stability of self-consistent field wave functions, and (4) predicting molecular properties via alchemical perturbations.

Molecular Dynamics with Constrained Nuclear Electronic Orbital Density Functional Theory: Accurate Vibrational Spectra from Efficient Incorporation of Nuclear Quantum Effects

Nuclear quantum effects play a crucial role in many chemical and biological systems involving hydrogen atoms yet are difficult to include in practical molecular simulations. In this paper, we combine our recently developed methods of constrained nuclear-electronic orbital density functional theory (cNEO-DFT) and constrained minimized energy surface molecular dynamics (CMES-MD) to create a new method for accurately and efficiently describing nuclear quantum effects in molecular simulations. By use of this new method, dubbed cNEO-MD, the vibrational spectra of a set of small molecules are calculated and compared with those from conventional ab initio molecular dynamics (AIMD) as well as from experiments. With the same formal scaling, cNEO-MD greatly outperforms AIMD in describing the vibrational modes with significant hydrogen motion characters, demonstrating the promise of cNEO-MD for simulating chemical and biological systems with significant nuclear quantum effects.

Correlation-Consistent Gaussian Basis Sets for Solids Made Simple

The rapidly growing interest in simulating condensed-phase materials using quantum chemistry methods calls for a library of high-quality Gaussian basis sets suitable for periodic calculations. Unfortunately, most standard Gaussian basis sets commonly used in molecular simulation show significant linear dependencies when used in close-packed solids, leading to severe numerical issues that hamper the convergence to the complete basis set (CBS) limit, especially in correlated calculations. In this work, we revisit Dunning's strategy for construction of correlation-consistent basis sets and examine the relationship between accuracy and numerical stability in periodic settings. We find that limiting the number of primitive functions avoids the appearance of problematic small exponents while still providing smooth convergence to the CBS limit. As an example, we generate double-, triple-, and quadruple-ζ correlation-consistent Gaussian basis sets for periodic calculations with Goedecker-Teter-Hutter (GTH) pseudopotentials. Our basis sets cover the main-group elements from the first three rows of the periodic table. Especially for atoms on the left side of the periodic table, our basis sets are less diffuse than those used in molecular calculations. We verify the fast and reliable convergence to the CBS limit in both Hartree-Fock and post-Hartree-Fock (MP2) calculations, using a diverse test set of 19 semiconductors and insulators.

Compressing Many-Body Fermion Operators under Unitary Constraints

The most efficient known quantum circuits for preparing unitary coupled cluster states and applying Trotter steps of the arbitrary basis electronic structure Hamiltonian involve interleaved sequences of Fermionic Gaussian circuits and Ising interaction-type circuits. These circuits arise from factorizing the two-body operators generating those unitaries as a sum of squared one-body operators that are simulated using product formulas. We introduce a numerical algorithm for performing this factorization that has an iteration complexity no worse than single particle basis transformations of the two-body operators and often results in many times fewer squared one-body operators in the sum of squares, compared to the analytical decompositions. As an application of this numerical procedure, we demonstrate that our protocol can be used to approximate generic unitary coupled cluster operators and prepare the necessary high-quality initial states for techniques (like ADAPT-VQE) that iteratively construct approximations to the ground state.

Coupled- and Independent-Trajectory Approaches Based on the Exact Factorization Using the PyUNIxMD Package

We present mixed quantum-classical approaches based on the exact factorization framework. The electron-nuclear correlation term in the exact factorization enables us to deal with quantum coherences by accounting for electronic and nuclear nonadiabatic couplings effectively within classical nuclei approximation. We compare coupled- and independent-trajectory approximations with each other to understand algorithms in description of the bifurcation of nuclear wave packets and the correct spatial distribution of electronic wave functions along with nuclear trajectories. Finally, we show numerical results for comparisons of coupled- and independent-trajectory approaches for the photoisomerization of a protonated Schiff base from excited state molecular dynamics (ESMD) simulations with the recently developed Python-based ESMD code, namely, the PyUNIxMD program.

Pure State v-Representability of Density Matrix Embedding Theory

Density matrix embedding theory (DMET) formally requires the matching of density matrix blocks obtained from high-level and low-level theories, but this is sometimes not achievable in practical calculations. In such a case, the global band gap of the low-level theory vanishes, and this can require additional numerical considerations. We find that both the violation of the exact matching condition and the vanishing low-level gap are related to the assumption that the high-level density matrix blocks are noninteracting pure-state v-representable (NI-PS-V), which assumes that the low-level density matrix is constructed following the Aufbau principle. To relax the NI-PS-V condition, we develop an augmented Lagrangian method to match the density matrix blocks without referring to the Aufbau principle. Numerical results for the 2D Hubbard and hydrogen model systems indicate that, in some challenging scenarios, the relaxation of the Aufbau principle directly leads to exact matching of the density matrix blocks, which also yields improved accuracy.

Low-Scaling Excited State Calculation Using the Block Interaction Product State

We develop an automatic and efficient scheme for the accurate construction of the bases for excitonic models, which can enable "black-box" excited state structure calculations for large molecular systems. These new and optimized bases, which are named the block interaction product state (BIPS), can be expressed as the direct products of the local states for each chromophore. Each chromophore's local states are selected by diagonalization of its reduced density matrix, which is obtained by quantum chemical calculation of the small subsystem composed of the chromophore and its nearest neighbors. We implemented the BIPS framework with fragment-based calculations considering two- and three-body interactions. Test calculations for eight different molecular aggregates demonstrate that this framework provides an accurate description of not only the excitation energies but also the first-order wave function properties (dipole moment and transition dipole moment) of the low-lying excited states at a low-scaling computational cost.

Regularized second-order correlation methods for extended systems

Second-order Møller-Plesset perturbation theory (MP2) constitutes the simplest form of many-body wavefunction theory and often provides a good compromise between efficiency and accuracy. There are, however, well-known limitations to this approach. In particular, MP2 is known to fail or diverge for some prototypical condensed matter systems like the homogeneous electron gas (HEG) and to overestimate dispersion-driven interactions in strongly polarizable systems. In this paper, we explore how the issues of MP2 for metallic, polarizable, and strongly correlated periodic systems can be ameliorated through regularization. To this end, two regularized second-order methods (including a new, size-extensive Brillouin-Wigner approach) are applied to the HEG, the one-dimensional Hubbard model, and the graphene-water interaction. We find that regularization consistently leads to improvements over the MP2 baseline and that different regularizers are appropriate for the various systems.

Ab initio calculations in atoms, molecules, and solids, treating spin-orbit coupling and electron interaction on an equal footing

We incorporate explicit, non-perturbative treatment of spin-orbit coupling into ab initio auxiliary-field quantum Monte Carlo (AFQMC) calculations. The approach allows a general computational framework for molecular and bulk systems in which material specificity, electron correlation, and spin-orbit coupling effects can be captured accurately and on an equal footing, with favorable computational scaling vs system size. We adopt relativistic effective-core potentials that have been obtained by fitting to fully relativistic data and that have demonstrated a high degree of reliability and transferability in molecular systems. This results in a two-component spin-coupled Hamiltonian, which is then treated by generalizing the ab initio AFQMC approach. We demonstrate the method by computing the electron affinity in Pb, the bond dissociation energy in Br₂ and I₂, and solid Bi.

Equivariant representations for molecular Hamiltonians and N-center atomic-scale properties

Symmetry considerations are at the core of the major frameworks used to provide an effective mathematical representation of atomic configurations that is then used in machine-learning models to predict the properties associated with each structure. In most cases, the models rely on a description of atom-centered environments and are suitable to learn atomic properties or global observables that can be decomposed into atomic contributions. Many quantities that are relevant for quantum mechanical calculations, however-most notably the single-particle Hamiltonian matrix when written in an atomic orbital basis-are not associated with a single center, but with two (or more) atoms in the structure. We discuss a family of structural descriptors that generalize the very successful atom-centered density correlation features to the N-center case and show, in particular, how this construction can be applied to efficiently learn the matrix elements of the (effective) single-particle Hamiltonian written in an atom-centered orbital basis. These N-center features are fully equivariant-not only in terms of translations and rotations but also in terms of permutations of the indices associated with the atoms-and are suitable to construct symmetry-adapted machine-learning models of new classes of properties of molecules and materials.

Fermi.jl: A Modern Design for Quantum Chemistry

Approximating molecular wave functions involves heavy numerical effort; therefore, codes for such tasks are written completely or partially in efficient languages such as C, C++, and Fortran. While these tools are dominant throughout quantum chemistry packages, the efficient development of new methods is often hindered by the complexity associated with code development. In order to ameliorate this scenario, some software packages take a dual approach where a simpler, higher-level language, such as Python, substitutes the traditional ones wherever performance is not critical. Julia is a novel, dynamically typed, programming language that aims to solve this two-language problem. It gained attention because of its modern and intuitive design, while still being highly optimized to compete with "low-level" languages. Recently, some chemistry-related projects have emerged exploring the capabilities of Julia. Herein, we introduce the quantum chemistry package Fermi.jl, which contains the first implementations of post-Hartree-Fock methods written in Julia. Its design makes use of many Julia core features, including multiple dispatch, metaprogramming, and interactive usage. Fermi.jl is a modular package, where new methods and implementations can be easily added to the existing code. Furthermore, it is designed to maximize code reusability by relying on general functions with specialized methods for particular cases. The feasibility of the project is explored through evaluating the performance of popular ab initio methods. It is our hope that this project motivates the usage of Julia within the community and brings new contributions into Fermi.jl.

Simulating X-ray photoelectron spectra with strong electron correlation using multireference algebraic diagrammatic construction theory

A new theoretical approach for the simulations of X-ray photoelectron spectra of strongly correlated molecular systems that combines multireference algebraic diagrammatic construction theory (MR-ADC) with a core–valence separation (CVS) technique.

Multicomponent MP4 and the inclusion of triple excitations in multicomponent many-body methods

This study implements the full multicomponent third-order (MP3) and fourth-order (MP4) many-body perturbation theory methods for the first time. Previous multicomponent studies have only implemented a subset of the full contributions, and the present implementation is the first multicomponent many-body method to include any connected triples contribution to the electron-proton correlation energy. The multicomponent MP3 method is shown to be comparable in accuracy to the multicomponent coupled-cluster doubles method for the calculation of proton affinities, while the multicomponent MP4 method is of similar accuracy as the multicomponent coupled-cluster singles and doubles method. From the results in this study, it is hypothesized that the relative accuracy of multicomponent methods is more similar to their single-component counterparts than previously assumed. It is demonstrated that for multicomponent MP4, the fourth-order triple-excitation contributions can be split into electron-electron and electron-proton contributions and the electron-electron contributions ignored with very little loss of accuracy of protonic properties.

Removing the Deadwood from DFT/MRCI Wave Functions: The p-DFT/MRCI Method

The combined density functional theory and multireference configuration interaction (DFT/MRCI) method is a powerful tool for the calculation of excited electronic states of large molecules. There exists, however, a large amount of superfluous configurations in a typical DFT/MRCI wave function. We show that this deadwood may be effectively removed using a simple configuration pruning algorithm based on second-order Epstein-Nesbet perturbation theory. The resulting method, which we denote p-DFT/MRCI, is shown to result in orders of magnitude saving in computational timings, while retaining the accuracy of the original DFT/MRCI method.

Chemical bonding theories as guides for self-interaction corrected solutions: Multiple local minima and symmetry breaking

Fermi-Löwdin orbitals (FLOs) are a special set of localized orbitals, which have become commonly used in combination with the Perdew-Zunger self-interaction correction (SIC) in the FLO-SIC method. The FLOs are obtained for a set of occupied orbitals by specifying a classical position for each electron. These positions are known as Fermi-orbital descriptors (FODs), and they have a clear relation to chemical bonding. In this study, we show how FLOs and FODs can be used to initialize, interpret, and justify SIC solutions in a common chemical picture, both within FLO-SIC and in traditional variational SIC, and to locate distinct local minima in either of these approaches. We demonstrate that FLOs based on Lewis theory lead to symmetry breaking for benzene-the electron density is found to break symmetry already at the symmetric molecular structure-while ones from Linnett's double-quartet theory reproduce symmetric electron densities and molecular geometries. Introducing a benchmark set of 16 planar cyclic molecules, we show that using Lewis theory as the starting point can lead to artifactual dipole moments of up to 1 D, while Linnett SIC dipole moments are in better agreement with experimental values. We suggest using the dipole moment as a diagnostic of symmetry breaking in SIC and monitoring it in all SIC calculations. We show that Linnett structures can often be seen as superpositions of Lewis structures and propose Linnett structures as a simple way to describe aromatic systems in SIC with reduced symmetry breaking. The role of hovering FODs is also briefly discussed.

Scalable Molecular GW Calculations: Valence and Core Spectra

We present a scalable implementation of the GW approximation using Gaussian atomic orbitals to study the valence and core ionization spectroscopies of molecules. The implementation of the standard spectral decomposition approach to the screened-Coulomb interaction, as well as a contour-deformation method, is described. We have implemented both of these approaches using the robust variational fitting approximation to the four-center electron repulsion integrals. We have utilized the MINRES solver with the contour-deformation approach to reduce the computational scaling by 1 order of magnitude. A complex heuristic in the quasiparticle equation solver further allows a speed-up of the computation of core and semicore ionization energies. Benchmark tests using the GW100 and CORE65 data sets and the carbon 1s binding energy of the well-studied ethyl trifluoroacetate, or ESCA molecule, were performed to validate the accuracy of our implementation. We also demonstrate and discuss the parallel performance and computational scaling of our implementation using a range of water clusters of increasing size.

Dipole Moment Calculations Using Multiconfiguration Pair-Density Functional Theory and Hybrid Multiconfiguration Pair-Density Functional Theory

The dipole moment is the molecular property that most directly indicates molecular polarity. The accuracy of computed dipole moments depends strongly on the quality of the calculated electron density, and the breakdown of single-reference methods for strongly correlated systems can lead to poor predictions of the dipole moments in those cases. Here, we derive the analytical expression for obtaining the electric dipole moment by multiconfiguration pair-density functional theory (MC-PDFT), and we assess the accuracy of MC-PDFT for predicting dipole moments at equilibrium and nonequilibrium geometries. We show that MC-PDFT dipole moment curves have reasonable behavior even for stretched geometries, and they significantly improve upon the CASSCF results by capturing more electron correlation. The analysis of a dataset consisting of 18 first-row transition-metal diatomics and 6 main-group polyatomic molecules with a multireference character suggests that MC-PDFT and its hybrid extension (HMC-PDFT) perform comparably to CASPT2 and MRCISD+Q methods and have a mean unsigned deviation of 0.2-0.3 D with respect to the best available dipole moment reference values. We explored the dependence of the predicted dipole moments upon the choice of the on-top density functional and active space, and we recommend the tPBE and hybrid tPBE0 on-top choices for the functionals combined with the moderate correlated-participating-orbitals scheme for selecting the active space. With these choices, the mean unsigned deviations (in debyes) of the calculated equilibrium dipole moments from the best estimates are 0.77 for CASSCF, 0.29 for MC-PDFT, 0.24 for HMC-PDFT, 0.28 for CASPT2, and 0.25 for MRCISD+Q. These results are encouraging because the computational cost of MC-PDFT or HMC-PDFT is largely reduced compared to the CASPT2 and MRCISD+Q methods.

Excited States of Crystalline Point Defects with Multireference Density Matrix Embedding Theory

Accurate and affordable methods to characterize the electronic structure of solids are important for targeted materials design. Embedding-based methods provide an appealing balance in the trade-off between cost and accuracy─particularly when studying localized phenomena. Here, we use the density matrix embedding theory (DMET) algorithm to study the electronic excitations in solid-state defects with a restricted open-shell Hartree-Fock (ROHF) bath and multireference impurity solvers, specifically, complete active space self-consistent field (CASSCF) and n-electron valence state second-order perturbation theory (NEVPT2). We apply the method to investigate the electronic excitations in an oxygen vacancy (OV) on a MgO(100) surface and find absolute deviations within 0.05 eV between DMET using the CASSCF/NEVPT2 solver, denoted as CAS-DMET/NEVPT2-DMET, and the nonembedded CASSCF/NEVPT2 approach. Next, we establish the practicality of DMET by extending it to larger supercells for the OV defect and a neutral silicon vacancy in diamond where the use of nonembedded CASSCF/NEVPT2 is extremely expensive.

Third-Order Many-Body Expansion of OSV-MP2 Wave Function for Low-Order Scaling Analytical Gradient Computation

We present a many-body expansion (MBE) formulation and implementation for efficient computation of analytical energy gradients from the orbital-specific-virtual second-order Møllet-Plesset perturbation theory (OSV-MP2) based on our earlier work (Zhou et al. J. Chem. Theory Comput. 2020, 16, 196-210). The third-order MBE(3) expansion of OSV-MP2 amplitudes and density matrices was developed to adopt the orbital-specific clustering and long-range termination schemes, which avoids term-by-term differentiations of the MBE energy bodies. We achieve better efficiency by exploiting the algorithmic sparsity that allows us to prune out insignificant fitting integrals and OSV relaxations. With these approximations, the present implementation is benchmarked on a range of molecules that show an economic scaling in the linear and quadratic regimes for computing MBE(3)-OSV-MP2 amplitude and gradient equations, respectively, and yields normal accuracy comparable to the original OSV-MP2 results. The MPI-3-based parallelism through shared memory one-sided communication is further developed for improving parallel scalability and memory accessibility by sorting the MBE(3) orbital clusters into independent tasks that are distributed on multiple processes across many nodes, supporting both global and local data locations in which selected MBE(3)-OSV-MP2 intermediates of different sizes are distinguished and accordingly placed. The accuracy and efficiency level of our MBE(3)-OSV-MP2 analytical gradient implementation is finally illustrated in two applications: we show that the subtle coordination structure differences of mechanically interlocked Cu-catenane complexes can be distinguished when tuning ligand lengths; and the porphycene molecular dynamics reveals the emergence of the vibrational signature arising from softened N-H stretching associated with hydrogen transfer, using an MP2 level of electron correlation and classical nuclei for the first time.

Constructing and representing exchange-correlation holes through artificial neural networks

One strategy to construct approximations to the exchange-correlation (XC) energy E_XC of Kohn-Sham density functional theory relies on physical constraints satisfied by the XC hole ρ_XC(r, u). In the XC hole, the reference charge is located at r and u is the electron-electron separation. With mathematical intuition, a given set of physical constraints can be expressed in a formula, yielding an approximation to ρ_XC(r, u) and the corresponding E_XC. Here, we adapt machine learning algorithms to partially automate the construction of X and XC holes. While machine learning usually relies on finding patterns in datasets and does not require physical insight, we focus entirely on the latter and develop a tool (ExMachina), consisting of the basic equations and their implementation, for the machine generation of approximations. To illustrate ExMachina, we apply it to calculate various model holes and show how to go beyond existing approximations.

Approaching the basis set limit in Gaussian-orbital-based periodic calculations with transferability: Performance of pure density functionals for simple semiconductors

Simulating solids with quantum chemistry methods and Gaussian-type orbitals (GTOs) has been gaining popularity. Nonetheless, there are few systematic studies that assess the basis set incompleteness error (BSIE) in these GTO-based simulations over a variety of solids. In this work, we report a GTO-based implementation for solids and apply it to address the basis set convergence issue. We employ a simple strategy to generate large uncontracted (unc) GTO basis sets that we call the unc-def2-GTH sets. These basis sets exhibit systematic improvement toward the basis set limit as well as good transferability based on application to a total of 43 simple semiconductors. Most notably, we found the BSIE of unc-def2-QZVP-GTH to be smaller than 0.7 mEh per atom in total energies and 20 meV in bandgaps for all systems considered here. Using unc-def2-QZVP-GTH, we report bandgap benchmarks of a combinatorially designed meta-generalized gradient approximation (mGGA) functional, B97M-rV, and show that B97M-rV performs similarly (a root-mean-square-deviation of 1.18 eV) to other modern mGGA functionals, M06-L (1.26 eV), MN15-L (1.29 eV), and Strongly Constrained and Appropriately Normed (SCAN) (1.20 eV). This represents a clear improvement over older pure functionals such as local density approximation (1.71 eV) and Perdew-Burke-Ernzerhof (PBE) (1.49 eV), although all these mGGAs are still far from being quantitatively accurate. We also provide several cautionary notes on the use of our uncontracted bases and on future research on GTO basis set development for solids.

Towards multistate multimode landscapes in singlet fission of pentacene: the dual role of charge-transfer states

Singlet fission duplicates triplet excitons for improving light harvesting efficiency. The presence of the interaction between electronic and nuclear degrees of freedom complicates the interpretation of correlated triplet pairs. We report a quantum chemistry study on the significance and subtleties of multistate and multimode pathways in forming triplet pair states of the pentacene dimer through a six-state vibronic-coupling Hamiltonian derived from many-electron adiabatic wavefunctions of an ab initio density matrix renormalization group. The resulting spin values of the singlet manifolds on each pentacene center are computed, and the varying spin nature can be distinguished clearly with respect to dimer stacking and vibronic progression. Our monomer spin assignments reveal the coexistence of both lower-lying weak and higher-lying strong charge transfer states which interact vibronically with the triplet pair state, providing important implications for its generation and separation occurring in vibronic regions. This work conveys the importance of the many-electron process requiring close low-lying singlet manifolds to determine the subtle fission details, and represents an important step for understanding vibronically resolved spin states and conversions underlying efficient singlet fission.

Absorption Spectra of Solids from Periodic Equation-of-Motion Coupled-Cluster Theory

We present ab initio absorption spectra of six three-dimensional semiconductors and insulators calculated using Gaussian-based periodic equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD). The spectra are calculated efficiently by solving a system of linear equations at each frequency, giving access to an energy range of tens of electronvolts without explicit enumeration of excited states. We assess the impact of cost-saving approximations associated with Brillouin zone sampling, frozen orbitals, and the partitioned EOM-CCSD approximation. Although our most converged spectra exhibit line shapes that are in good agreement with experimental spectra, they are uniformly shifted to higher energies by about 1 eV, which is not explained by the remaining finite-size errors. We tentatively attribute this discrepancy to a combination of vibrational effects and the remaining electron correlation, i.e., triple excitations and above.

Fixed- and Partial-Node Approximations in Slater Determinant Space for Molecules

We present a study of fixed- and partial-node approximations in Slater determinant basis sets, using full configuration interaction quantum Monte Carlo (FCIQMC) to perform sampling. Walker annihilation in the FCIQMC method allows partial-node simulations to be performed, relaxing the nodal constraint to converge to the FCI solution. This is applied to ab initio molecular systems, using symmetry-projected Jastrow mean-field wave functions for complete active space (CAS) problems. Convergence and the sign problem within the partial-node approximation are studied, which is shown to eventually be limited in its use due to the large walker populations required. However, the fixed-node approximation results in an accurate and practical method. We apply these approaches to various molecular systems and active spaces, including ferrocene and acenes. This also provides a test of symmetry-projected Jastrow mean-field wave functions in variational Monte Carlo for a new set of problems. For trans-polyacetylene molecules and acenes, we find that the time to perform a constant number of fixed-node FCIQMC iterations scales as O(N1.44) and O(N1.75), respectively, resulting in an efficient method for CAS-based problems that can be applied accurately to large active spaces.

Multireference Algebraic Diagrammatic Construction Theory for Excited States: Extended Second-Order Implementation and Benchmark

We present an implementation and benchmark of new approximations in multireference algebraic diagrammatic construction theory for simulations of neutral electronic excitations and UV/vis spectra of strongly correlated molecular systems (MR-ADC). Following our work on the first-order MR-ADC approximation [J. Chem. Phys. 2018, 149, 204113], we report the strict and extended second-order MR-ADC methods (MR-ADC(2) and MR-ADC(2)-X) that combine the description of static and dynamic electron correlation in the ground and excited electronic states without relying on state-averaged reference wave functions. We present an extensive benchmark of the new MR-ADC methods for excited states in several small molecules, including the carbon dimer, ethylene, and butadiene. Our results demonstrate that, for weakly correlated electronic states, the MR-ADC(2) and MR-ADC(2)-X methods outperform the third-order single-reference ADC approximation and are competitive with the results from equation-of-motion coupled cluster theory. For states with multireference character, the performance of the MR-ADC methods is similar to that of an N-electron valence perturbation theory. In contrast to conventional multireference perturbation theories, the MR-ADC methods have many attractive features, such as a straightforward and efficient calculation of excited-state properties and a direct access to excitations outside of the frontier (active) orbitals.

Plasmon Character Index: An Accurate and Efficient Metric for Identifying and Quantifying Plasmons in Molecules

Plasmons, which are collective and coherent oscillations of charge carriers driven by an external field, play an important role in applications such as solar energy harvesting, sensing, and catalysis. Conventionally, plasmons are found in bulk and nanomaterials and can be described with classical electrodynamics. In recent years, plasmons have also been identified in molecules, and these molecules have been utilized to build plasmonic devices. As molecular plasmons can no longer be described by classical electrodynamics, a description using quantum mechanics is necessary. In this Letter, we develop a quantum metric to accurately and efficiently identify and quantify plasmons in molecules. A number, which we call the plasmon character index (PCI), can be calculated for each electronic excited state and describes the plasmonicity of the excitation. PCI is developed from the collective and coherent excitation picture in orbitals and shows excellent agreement with the predictions from scaled time-dependent density functional theory but is vastly more computationally efficient. Therefore, PCI can be a useful tool in identifying and quantifying plasmons and will inform the rational design of plasmonic molecules and nanoclusters.

Chebyshev Matrix Product States with Canonical Orthogonalization for Spectral Functions of Many-Body Systems

We propose a method to calculate the spectral functions of many-body systems by Chebyshev expansion in the framework of matrix product states coupled with canonical orthogonalization (coCheMPS). The canonical orthogonalization can improve the accuracy and efficiency significantly because the orthogonalized Chebyshev vectors can provide an ideal basis for constructing the effective Hamiltonian in which the exact recurrence relation can be retained. In addition, not only the spectral function but also the excited states and eigenenergies can be directly calculated, which is usually impossible for other MPS-based methods such as time-dependent formalism or correction vector. The remarkable accuracy and efficiency of coCheMPS over other methods are demonstrated by calculating the spectral functions of spin chain and ab initio hydrogen chain. For the first time we demonstrate that Chebyshev MPS can be used to deal with ab initio electronic Hamiltonian effectively. We emphasize the strength of coCheMPS to calculate the low excited states of systems with sparse discrete spectrum. We also caution the application for electron-phonon systems with dense density of states.

Computing molecular excited states on a D-Wave quantum annealer

The possibility of using quantum computers for electronic structure calculations has opened up a promising avenue for computational chemistry. Towards this direction, numerous algorithmic advances have been made in the last five years. The potential of quantum annealers, which are the prototypes of adiabatic quantum computers, is yet to be fully explored. In this work, we demonstrate the use of a D-Wave quantum annealer for the calculation of excited electronic states of molecular systems. These simulations play an important role in a number of areas, such as photovoltaics, semiconductor technology and nanoscience. The excited states are treated using two methods, time-dependent Hartree–Fock (TDHF) and time-dependent density-functional theory (TDDFT), both within a commonly used Tamm–Dancoff approximation (TDA). The resulting TDA eigenvalue equations are solved on a D-Wave quantum annealer using the Quantum Annealer Eigensolver (QAE), developed previously. The method is shown to reproduce a typical basis set convergence on the example \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2$$\end{document}H2 molecule and is also applied to several other molecular species. Characteristic properties such as transition dipole moments and oscillator strengths are computed as well. Three potential energy profiles for excited states are computed for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3$$\end{document}NH3 as a function of the molecular geometry. Similar to previous studies, the accuracy of the method is dependent on the accuracy of the intermediate meta-heuristic software called qbsolv.

Nucleus-electron correlation revising molecular bonding fingerprints from the exact wavefunction factorization

We present a novel theory and implementation for computing coupled electronic and quantal nuclear subsystems on a single potential energy surface, moving beyond the standard Born-Oppenheimer (BO) separation of nuclei and electrons. We formulate an exact self-consistent nucleus-electron embedding potential from the single product molecular wavefunction and demonstrate that the fundamental behavior of the correlated nucleus-electron can be computed for mean-field electrons that are responsive to a quantal anharmonic vibration of selected nuclei in a discrete variable representation. Geometric gauge choices are discussed and necessary for formulating energy invariant biorthogonal electronic equations. Our method is further applied to characterize vibrationally averaged molecular bonding properties of molecular energetics, bond lengths, and protonic and electron densities. Moreover, post-Hartree-Fock electron correlation can be conveniently computed on the basis of nucleus-electron coupled molecular orbitals, as demonstrated for correlated models of second-order Møllet-Plesset perturbation and full configuration interaction theories. Our approach not only accurately quantifies non-classical nucleus-electron couplings for revising molecular bonding properties but also provides an alternative time-independent approach for deploying non-BO molecular quantum chemistry.

GQCP: The Ghent Quantum Chemistry Package

The Ghent Quantum Chemistry Package (GQCP) is an open-source electronic structure software package that aims to provide an intuitive and expressive software framework for electronic structure software development. Its high-level interfaces (accessible through C++ and Python) have been specifically designed to correspond to theoretical concepts, while retaining access to lower-level intermediates and allowing structural run-time modifications of quantum chemical solvers. GQCP focuses on providing quantum chemical method developers with the computational "building blocks" that allow them to flexibly develop proof of principle implementations for new methods and applications up to the level of two-component spinor bases.

Improving MP2 bandgaps with low-scaling approximations to EOM-CCSD

Despite its reasonable accuracy for ground-state properties of semiconductors and insulators, second-order Møller-Plesset perturbation theory (MP2) significantly underestimates bandgaps. In this work, we evaluate the bandgap predictions of partitioned equation-of-motion MP2 (P-EOM-MP2), which is a second-order approximation to EOM coupled-cluster theory with single and double excitations. On a test set of elemental and binary semiconductors and insulators, we find that P-EOM-MP2 overestimates bandgaps by 0.3 eV on average, which can be compared to the underestimation by 0.6 eV on average exhibited by the G₀W₀ approximation with a Perdew-Burke-Ernzerhof reference. We show that P-EOM-MP2, when interpreted as a Green's function-based theory, has a self-energy that includes all first- and second-order diagrams and a few third-order diagrams. We find that the GW approximation performs better for materials with small gaps and P-EOM-MP2 performs better for materials with large gaps, which we attribute to their superior treatment of screening and exchange, respectively.

Scalable and Predictive Spectra of Correlated Molecules with Moment Truncated Iterated Perturbation Theory

A reliable and efficient computation of the entire single-particle spectrum of correlated molecules is an outstanding challenge in the field of quantum chemistry, with standard density functional theory approaches often giving an inadequate description of excitation energies and gaps. In this work, we expand upon a recently introduced approach that relies on a fully self-consistent many-body perturbation theory coupled to a nonperturbative truncation of the effective dynamics at each step. We show that this yields a low-scaling and accurate method across a diverse benchmark test set that is capable of treating moderate levels of strong correlation effects, and we detail an efficient implementation for applications involving up to ∼1000 orbitals on parallel resources. We then use this method to characterize the spectral properties of the antimalarial drug molecule artemisinin, resolving discrepancies in previous works concerning the active sites of the lowest-energy fundamental excitations of the system.

Cluster many-body expansion: A many-body expansion of the electron correlation energy about a cluster mean field reference

The many-body expansion (MBE) is an efficient tool that has a long history of use for calculating interaction energies, binding energies, lattice energies, and so on. In the past, applications of MBE to correlation energy have been unfeasible for large systems, but recent improvements to computing resources have sparked renewed interest in capturing the correlation energy using the generalized nth order Bethe-Goldstone equation. In this work, we extend this approach, originally proposed for a Slater determinant, to a tensor product state (TPS) based wavefunction. By partitioning the active space into smaller orbital clusters, our approach starts from a cluster mean field reference TPS configuration and includes the correlation contribution of the excited TPSs using the MBE. This method, named cluster MBE (cMBE), improves the convergence of MBE at lower orders compared to directly doing a block-based MBE from a RHF reference. We present numerical results for strongly correlated systems, such as the one- and two-dimensional Hubbard models and the chromium dimer. The performance of the cMBE method is also tested by partitioning the extended π space of several large π-conjugated systems, including a graphene nano-sheet with a very large active space of 114 electrons in 114 orbitals, which would require 10⁶⁶ determinants for the exact FCI solution.

Conservation laws in coupled cluster dynamics at finite temperature

We extend the finite-temperature Keldysh non-equilibrium coupled cluster theory (Keldysh-CC) [A. F. White and G. K.-L. Chan, J. Chem. Theory Comput. 15, 6137-6253 (2019)] to include a time-dependent orbital basis. When chosen to minimize the action, such a basis restores local and global conservation laws (Ehrenfest's theorem) for all one-particle properties while remaining energy conserving for time-independent Hamiltonians. We present the time-dependent Keldysh orbital-optimized coupled cluster doubles method in analogy with the formalism for zero-temperature dynamics, extended to finite temperatures through the time-dependent action on the Keldysh contour. To demonstrate the conservation property and understand the numerical performance of the method, we apply it to several problems of non-equilibrium finite-temperature dynamics: a 1D Hubbard model with a time-dependent Peierls phase, laser driving of molecular H₂, driven dynamics in warm-dense silicon, and transport in the single impurity Anderson model.

Efficient Evaluation of Two-Center Gaussian Integrals in Periodic Systems

By using Poisson's summation formula, we calculate periodic integrals over Gaussian basis functions by partitioning the lattice summations between the real and reciprocal space, where both sums converge exponentially fast with a large exponent. We demonstrate that the summation can be performed efficiently to calculate two-center Gaussian integrals over various kernels including overlap, kinetic, and Coulomb. The summation in real space is performed using an efficient flavor of the McMurchie-Davidson recurrence relation. The expressions for performing summation in the reciprocal space are also derived and implemented. The algorithm for reciprocal space summation allows us to reuse several terms and leads to a significant improvement in efficiency when highly contracted basis functions with large exponents are used. We find that the resulting algorithm is only between a factor of 5 and 15 slower than that for molecular integrals, indicating the very small number of terms needed in both the real and reciprocal space summations. An outline of the algorithm for calculating three-center Coulomb integrals is also provided.

Decomposed Mean-Field Simulations of Local Properties in Condensed Phases

The present work demonstrates a robust protocol for probing localized electronic structure in condensed-phase systems, operating in terms of a recently proposed theory for decomposing the results of Kohn-Sham density functional theory in a basis of spatially localized molecular orbitals. In an initial application to liquid, ambient water and the assessment of the solvation energy and the embedded dipole moment of H₂O in solution, we find that both properties are amplified on average-in accordance with expectation-and that correlations are indeed observed to exist between them. However, the simulated solvent-induced shift to the dipole moment of water is found to be significantly dampened with respect to typical literature values. The local nature of our methodology has further allowed us to evaluate the convergence of bulk properties with respect to the extent of the underlying one-electron basis set, ranging from single-ζ to full (augmented) quadruple-ζ quality. Albeit a pilot example, our work paves the way toward future studies of local effects and defects in more complex phases, e.g., liquid mixtures and even solid-state crystals.

Molecular vibrational frequencies from analytic Hessian of constrained nuclear-electronic orbital density functional theory

Nuclear quantum effects are important in a variety of chemical and biological processes. The constrained nuclear-electronic orbital density functional theory (cNEO-DFT) has been developed to include nuclear quantum effects in energy surfaces. Herein, we develop the analytic Hessian for cNEO-DFT energy with respect to the change in nuclear (expectation) positions, which can be used to characterize stationary points on energy surfaces and compute molecular vibrational frequencies. This is achieved by constructing and solving the multicomponent cNEO coupled-perturbed Kohn-Sham (cNEO-CPKS) equations, which describe the response of electronic and nuclear orbitals to the displacement of nuclear (expectation) positions. With the analytic Hessian, the vibrational frequencies of a series of small molecules are calculated and compared to those from conventional DFT Hessian calculations as well as those from the vibrational second-order perturbation theory (VPT2). It is found that even with a harmonic treatment, cNEO-DFT significantly outperforms DFT and is comparable to DFT-VPT2 in the description of vibrational frequencies in regular polyatomic molecules. Furthermore, cNEO-DFT can reasonably describe the proton transfer modes in systems with a shared proton, whereas DFT-VPT2 often faces great challenges. Our results suggest the importance of nuclear quantum effects in molecular vibrations, and cNEO-DFT is an accurate and inexpensive method to describe molecular vibrations.

Efficient Open-Source Implementations of Linear-Scaling Polarizable Embedding: Use Octrees to Save the Trees

We present open-source implementations of the linear-scaling fast multipole method (FMM) within the polarizable embedding (PE) model for efficient treatment of large polarizable environments in computational spectroscopy simulations. The implementations are tested for accuracy, efficiency, and usability on model systems as well as more realistic biomolecular systems. We explain how FMM parameters affect the calculation of molecular properties and show that PE calculations employing FMM can be carried out in a black-box manner. The efficiency of the linear-scaling approach is demonstrated by simulating the UV/vis spectrum of a chromophore in an environment of more than 1 million polarizable sites. Our implementations are interfaced to several open-source quantum chemistry programs, making computational spectroscopy simulations within the PE model and FMM available to a large variety of methods and a broad user base.

Metal Oxo-Fluoride Molecules OnMF2 (M = Mn and Fe; n = 1-4) and O2MnF: Matrix Infrared Spectra and Quantum Chemistry

On reacting laser-ablated manganese or iron difluorides with O₂ or O₃ during codeposition in solid neon or argon, infrared absorptions of several new metal oxo-fluoride molecules, including OMF₂, (η¹-O₂)MF₂, (η²-O₃)MF₂, (η¹-O₂)₂MF₂ (M = Mn and Fe), and O₂MnF, have been observed. Quantum chemical density functional and multiconfiguration wavefunction calculations have been applied to characterize these new products by their geometric and electronic structures, vibrations, charges, and bonding. The assignment of the main vibrational absorptions as dominant symmetric or antisymmetric M-F or M-O stretching modes is confirmed by oxygen isotopic shifts and quantum chemical calculations of frequencies and thermal stabilities. The tendency of Fe to form polyoxygen complexes in lower oxidation states than the preceding element Mn is affirmed experimentally and supported theoretically. The M-F stretching frequencies of the isolated metal oxo-fluorides may provide a scale for the local charge on the MF₂ sites in active energy conversion systems. The study of these species provides insights for understanding the trend of oxidation state changes across the transition-metal series.

Correlation-driven phenomena in periodic molecular systems from variational two-electron reduced density matrix theory

Correlation-driven phenomena in molecular periodic systems are challenging to predict computationally not only because such systems are periodically infinite but also because they are typically strongly correlated. Here, we generalize the variational two-electron reduced density matrix (2-RDM) theory to compute the energies and properties of strongly correlated periodic systems. The 2-RDM of the unit cell is directly computed subject to necessary N-representability conditions such that the unit-cell 2-RDM represents at least one N-electron density matrix. Two canonical but non-trivial systems, periodic metallic hydrogen chains and periodic acenes, are treated to demonstrate the methodology. We show that while single-reference correlation theories do not capture the strong (static) correlation effects in either of these molecular systems, the periodic variational 2-RDM theory predicts the Mott metal-to-insulator transition in the hydrogen chains and the length-dependent polyradical formation in acenes. For both hydrogen chains and acenes, the periodic calculations are compared with previous non-periodic calculations with the results showing a significant change in energies and increase in the electron correlation from the periodic boundary conditions. The 2-RDM theory, which allows for much larger active spaces than are traditionally possible, is applicable to studying correlation-driven phenomena in general periodic molecular solids and materials.

Convergence to the fixed-node limit in deep variational Monte Carlo

Variational quantum Monte Carlo (QMC) is an ab initio method for solving the electronic Schrödinger equation that is exact in principle, but limited by the flexibility of the available Ansätze in practice. The recently introduced deep QMC approach, specifically two deep-neural-network Ansätze PauliNet and FermiNet, allows variational QMC to reach the accuracy of diffusion QMC, but little is understood about the convergence behavior of such Ansätze. Here, we analyze how deep variational QMC approaches the fixed-node limit with increasing network size. First, we demonstrate that a deep neural network can overcome the limitations of a small basis set and reach the mean-field (MF) complete-basis-set limit. Moving to electron correlation, we then perform an extensive hyperparameter scan of a deep Jastrow factor for LiH and H₄ and find that variational energies at the fixed-node limit can be obtained with a sufficiently large network. Finally, we benchmark MF and many-body Ansätze on H₂O, increasing the fraction of recovered fixed-node correlation energy of single-determinant Slater-Jastrow-type Ansätze by half an order of magnitude compared to previous variational QMC results, and demonstrate that a single-determinant Slater-Jastrow-backflow version of the Ansatz overcomes the fixed-node limitations. This analysis helps understand the superb accuracy of deep variational Ansätze in comparison to the traditional trial wavefunctions at the respective level of theory and will guide future improvements of the neural-network architectures in deep QMC.

Sulfur Molecules in Space by X-rays: A Computational Study

X-ray astronomy lacks high resolution spectra of interstellar dust analogues and molecules, severely hampering interstellar medium studies based on upcoming X-ray missions. Various theoretical approaches may be used to address this problem, but they must first be shown to reproduce reliable spectra compared to the experiment. In this work, we calculate the sulfur K edge X-ray absorption spectra of H₂S, SO₂, and OCS, whose spectra are already known from X-ray experiments and predict the X-ray spectrum of CS, which as far as we are aware has not been measured, thereby hampering its detection by X-ray telescopes. We chose these four molecules as the astrochemistry of sulfur is an unsolved problem and as the four molecules are already known to exist in space. We consider three types of methods for modeling the X-ray spectra: more accurate calculations with the algebraic-diagrammatic construction (ADC) and the CC2, CCSD, and CC3 coupled cluster (CC) approaches as well as more affordable ones with transition potential density functional theory (TP-DFT). A comparison of our computational results to previously reported experimental spectra shows that the core-valence separation (CVS) approaches CVS-ADC(2)-x and CVS-CC3 generally yield a good qualitative level of agreement with the experiment, suggesting that they can be used for interpreting measured spectra, while the TP-DFT method is not reliable for these molecules. However, quantitative agreement with the experiment is still outside the reach of the computational methods studied in this work.

Benchmarking Magnetizabilities with Recent Density Functionals

We have assessed the accuracy of the magnetic properties of a set of 51 density functional approximations, including both recently published and already established functionals. The accuracy assessment considers a series of 27 small molecules and is based on comparing the predicted magnetizabilities to literature reference values calculated using coupled-cluster theory with full singles and doubles and perturbative triples [CCSD(T)] employing large basis sets. The most accurate magnetizabilities, defined as the smallest mean absolute error, are obtained with the BHandHLYP functional. Three of the six studied Berkeley functionals and the three range-separated Florida functionals also yield accurate magnetizabilities. Also, some older functionals like CAM-B3LYP, KT1, BHLYP (BHandH), B3LYP, and PBE0 perform rather well. In contrast, unsatisfactory performance is generally obtained with Minnesota functionals, which are therefore not recommended for calculations of magnetically induced current density susceptibilities and related magnetic properties such as magnetizabilities and nuclear magnetic shieldings. We also demonstrate that magnetizabilities can be calculated by numerical integration of magnetizability density; we have implemented this approach as a new feature in the gauge-including magnetically induced current (GIMIC) method. Magnetizabilities can be calculated from magnetically induced current density susceptibilities within this approach even when analytical approaches for magnetizabilities as the second derivative of the energy have not been implemented. The magnetizability density can also be visualized, providing additional information that is not otherwise easily accessible on the spatial origin of magnetizabilities.