Embedded equation-of-motion coupled-cluster theory for electronic excitation, ionisation, electron attachment, and electronic resonances

Overview of Developments in the MRCC Program System

mrcc is a versatile suite of quantum chemistry programs designed for accurate ab initio and density functional theory (DFT) calculations. This contribution outlines the general features and recent developments of the package. The most popular features include the open-ended coupled-cluster (CC) code, state-of-the-art CC singles and doubles with perturbative triples [CCSD(T)], second-order algebraic-diagrammatic construction, and combined wave function theory-DFT approaches. Cost-reduction techniques are implemented, such as natural orbital (NO), local NO (LNO), and natural auxiliary function approximations, which significantly decrease the computational demands of these methods. This paper also details the method developments made over the past five years, including efficient schemes to approach the complete basis set limit for CCSD(T) and the extension of our LNO-CCSD(T) method to open-shell systems. Additionally, we discuss the new approximations introduced to accelerate the self-consistent field procedure and the cost-reduction techniques elaborated for analytic gradient calculations at various levels. Furthermore, embedding techniques and novel range-separated double-hybrid functionals are presented for excited-state calculations, while the extension of the theories established to describe core excitations and ionized states is also discussed. For academic purposes, the program and its source code are available free of charge, and its commercial use is also facilitated.

Coupled-cluster treatment of complex open-shell systems: the case of single-molecule magnets

We investigate the reliability of two cost-effective coupled-cluster methods for computing spin-state energetics and spin-related properties of a set of open-shell transition-metal complexes. Specifically, we employ the second-order approximate coupled-cluster singles and doubles (CC2) method and projection-based embedding that combines equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) with density functional theory (DFT). The performance of CC2 and EOM-CCSD-in-DFT is assessed against EOM-CCSD. The chosen test set includes two hexaaqua transition-metal complexes containing Fe(II) and Fe(III), and a large Co(II)-based single-molecule magnet with a non-aufbau ground state. We find that CC2 describes the excited states more accurately, reproducing EOM-CCSD excitation energies within 0.05 eV. However, EOM-CCSD-in-DFT excels in describing transition orbital angular momenta and spin-orbit couplings. Moreover, for the Co(II) molecular magnet, using EOM-CCSD-in-DFT eigenstates and spin-orbit couplings, we compute spin-reversal energy barriers, as well as temperature-dependent and field-dependent magnetizations and magnetic susceptibilities that closely match experimental values within spectroscopic accuracy. These results underscore the efficiency of CC2 in computing state energies of multi-configurational, open-shell systems and highlight the utility of the more cost-efficient EOM-CCSD-in-DFT for computing spin-orbit couplings and magnetic properties of complex and large molecular magnets.

Core Binding Energy Calculations: A Scalable Approach with the Quantum Embedding-Based Equation-of-Motion Coupled-Cluster Method

We investigated the use of density matrix embedding theory to facilitate the computation of core ionization energies (IPs) of large molecules at the equation-of-motion coupled-cluster singles doubles with perturbative triples (EOM-CCSD*) level in combination with the core-valence separation (CVS) approximation. The unembedded IP-CVS-EOM-CCSD* method with a triple-ζ basis set produced ionization energies within 1 eV of experiment with a standard deviation of ∼0.2 eV for the core65 data set. The embedded variant contributed very little systematic error relative to the unembedded method, with a mean unsigned error of 0.07 eV and a standard deviation of ∼0.1 eV, in exchange for accelerating the calculations by many orders of magnitude. By employing embedded EOM-CC methods, we computed the core ionization energies of the uracil hexamer, doped fullerene, and chlorophyll molecule, utilizing up to ∼4000 basis functions within 1 eV from experimental values. Such calculations are not currently possible with the unembedded EOM-CC method.

Projected Atomic Orbitals As Optimal Virtual Space for Excited State Projection-Based Embedding Calculations

The projected atomic orbital (PAO) technique is presented for the construction of virtual orbital spaces in projection-based embedding (PbE) applications. The proposed straightforward procedure produces a set of virtual orbitals that are used in the final, high-level calculation of the embedded active subsystem. The PAO scheme is demonstrated on intermolecular potentials of bimolecular complexes in ground and excited states, including Rydberg excitations. The results show the outstanding performance of the PbE method when used with PAO virtual orbitals compared with those produced using common orbital localization techniques. The good agreement of the resulting PbE potential curves with those from high-level ab initio dimer calculations, also in diffuse basis sets, confirms that the PAO technique can be suggested for future applications using top-down embedding methods.

A Perspective on Sustainable Computational Chemistry Software Development and Integration

The power of quantum chemistry to predict the ground and excited state properties of complex chemical systems has driven the development of computational quantum chemistry software, integrating advances in theory, applied mathematics, and computer science. The emergence of new computational paradigms associated with exascale technologies also poses significant challenges that require a flexible forward strategy to take full advantage of existing and forthcoming computational resources. In this context, the sustainability and interoperability of computational chemistry software development are among the most pressing issues. In this perspective, we discuss software infrastructure needs and investments with an eye to fully utilize exascale resources and provide unique computational tools for next-generation science problems and scientific discoveries.

Irradiation of Plutonium Tributyl Phosphate Complexes by Ionizing Alpha Particles: A Computational Study

The PUREX solvent extraction process, widely used for recovering uranium and plutonium from spent nuclear fuel, utilizes an organic solvent composed of tributyl phosphate (TBP). The emission of ionizing particles such as alpha particles, resulting from the decay of plutonium, makes the organic solvent vulnerable to degradation. Here, we study the ultrashort time alpha irradiation of tributylphosphate (TBP) and Pu(NO3)4(TBP)2 complex formed in the PUREX process. Electron dynamics is propagated by Real-Time-Dependent Auxiliary Density Functional Theory (RT-TD-ADFT). We investigate the use of previously proposed absorption boundary conditions (ABC) in the molecular orbital space to treat secondary electron emission. Basis set and exchange correlation functional effects with ABC are reported as well as a detailed analysis of the ABC parametrization. Preliminary results on the water molecule and then on TBP show that the phenomenological nature of the ABC parameters necessitates selecting appropriate values for each system under study. Irradiation of free and complexed TBP shows an influence of the ligands on the variation of atomic charges on the femtosecond time scale. An accumulation of atomic charges in the alkyl chains of TBP is observed in the case where the nitrate groups are predominantly irradiated. In addition, we find that the Pu atom regains its electric charge very rapidly after being hit by the projectile, with the coordination sphere serving as an electron reservoir to preserve its formal redox state. This study paves the road toward a full understanding of the degradation of organic extracants employed in the nuclear industry.

Elucidating the Role of Hydrogen Bonding in the Optical Spectroscopy of the Solvated Green Fluorescent Protein Chromophore: Using Machine Learning to Establish the Importance of High-Level Electronic Structure

Hydrogen bonding interactions with chromophores in chemical and biological environments play a key role in determining their electronic absorption and relaxation processes, which are manifested in their linear and multidimensional optical spectra. For chromophores in the condensed phase, the large number of atoms needed to simulate the environment has traditionally prohibited the use of high-level excited-state electronic structure methods. By leveraging transfer learning, we show how to construct machine-learned models to accurately predict the high-level excitation energies of a chromophore in solution from only 400 high-level calculations. We show that when the electronic excitations of the green fluorescent protein chromophore in water are treated using EOM-CCSD embedded in a DFT description of the solvent the optical spectrum is correctly captured and that this improvement arises from correctly treating the coupling of the electronic transition to electric fields, which leads to a larger response upon hydrogen bonding between the chromophore and water.

Benchmarking Aspects of Ab Initio Fragment Models for Accurate Excimer Potential Energy Surfaces

While Coupled-Cluster methods have been proven to provide an accurate description of excited electronic states, the scaling of the computational costs with the system size limits the degree for which these methods can be applied. In this work different aspects of fragment-based approaches are studied on noncovalently bound molecular complexes with interacting chromophores of the fragments, such as π-stacked nucleobases. The interaction of the fragments is considered at two distinct steps. First, the states localized on the fragments are described in the presence of the other fragment(s); for this we test two approaches. One method is founded on QM/MM principles, only including the electrostatic interaction between the fragments in the electronic structure calculation with Pauli repulsion and dispersion effects added separately. The other model, a Projection-based Embedding (PbE) using the Huzinaga equation, includes both electrostatic and Pauli repulsion and only needs to be augmented by dispersion interactions. In both schemes the extended Effective Fragment Potential (EFP2) method of Gordon et al. was found to provide an adequate correction for the missing terms. In the second step, the interaction of the localized chromophores is modeled for a proper description of the excitonic coupling. Here the inclusion of purely electrostatic contributions appears to be sufficient: it is found that the Coulomb part of the coupling provides accurate splitting of the energies of interacting chromophores that are separated by more than 4 Å.

Quantum computing in pharma: A multilayer embedding approach for near future applications

Quantum computers are special purpose machines that are expected to be particularly useful in simulating strongly correlated chemical systems. The quantum computer excels at treating a moderate number of orbitals within an active space in a fully quantum mechanical manner. We present a quantum phase estimation calculation on F₂ in a (2,2) active space on Rigetti's Aspen-11 QPU. While this is a promising start, it also underlines the need for carefully selecting the orbital spaces treated by the quantum computer. In this work, a scheme for selecting such an active space automatically is described and simulated results obtained using both the quantum phase estimation (QPE) and variational quantum eigensolver (VQE) algorithms are presented and combined with a subtractive method to enable accurate description of the environment. The active occupied space is selected from orbitals localized on the chemically relevant fragment of the molecule, while the corresponding virtual space is chosen based on the magnitude of interactions with the occupied space calculated from perturbation theory. This protocol is then applied to two chemical systems of pharmaceutical relevance: the enzyme [Fe] hydrogenase and the photosenzitizer temoporfin. While the sizes of the active spaces currently amenable to a quantum computational treatment are not enough to demonstrate quantum advantage, the procedure outlined here is applicable to any active space size, including those that are outside the reach of classical computation.

Performance of Multilevel Methods for Excited States

The performance of multilevel quantum chemical approaches, which utilize an atom-based system partitioning scheme to model various electronic excited states, is studied. The considered techniques include the mechanical-embedding (ME) of "our own N-layered integrated molecular orbital and molecular mechanics" (ONIOM) method, the point charge embedding (PCE), the electronic-embedding (EE) of ONIOM, the frozen density-embedding (FDE), the projector-based embedding (PbE), and our local domain-based correlation method. For the investigated multilevel approaches, the second-order algebraic-diagrammatic construction [ADC(2)] approach was utilized as the high-level method, which was embedded in either Hartree-Fock or a density functional environment. The XH-27 test set of Zech et al. [ J. Chem. Theory Comput., 2018, 14, 4028] was used for the assessment, where organic dyes interact with several solvent molecules. With the selection of the chromophores as active subsystems, we conclude that the most reliable approach is local domain-based ADC(2) [L-ADC(2)], and the least robust schemes are ONIOM-ME and ONIOM-EE. The PbE, FDE, and PCE techniques often approach the accuracy of the L-ADC(2) scheme, but their precision is far behind. The results suggest that a more conservative subsystem selection algorithm or the inclusion of subsystem charge-transfers is required for the atom-based cost-efficient methods to produce high-accuracy excitation energies.

Cherry-Picking Resolvents: Recovering the Valence Contribution in X-ray Two-Photon Absorption within the Core-Valence-Separated Equation-of-Motion Coupled-Cluster Response Theory

Calculations of first-order response wave functions in the X-ray regime often diverge within correlated frameworks such as equation-of-motion coupled-cluster singles and doubles (EOM-CCSD), a consequence of the coupling with the valence ionization continuum. Here, we extend our strategy of introducing a hierarchy of approximations to the EOM-EE-CCSD resolvent (or, inversely, the model Hamiltonian) involved in the response equations for the calculation of X-ray two-photon absorption (X2PA) cross sections. We exploit the frozen-core core-valence separation (fc-CVS) scheme to first decouple the core and valence Fock spaces, followed by a separate approximate treatment of the valence resolvent. We demonstrate the robust convergence of X-ray response calculations within this framework and compare X2PA spectra of small benchmark molecules with the previously reported density functional theory results.

Wavefunction embedding for molecular polaritons

Polaritonic chemistry relies on the strong light-matter interaction phenomena for altering the chemical reaction rates inside optical cavities. To explain and to understand these processes, the development of reliable theoretical models is essential. While computationally efficient quantum electrodynamics self-consistent field (QED-SCF) methods, such as quantum electrodynamics density functional theory (QEDFT) needs accurate functionals, quantum electrodynamics coupled cluster (QED-CC) methods provide a systematic increase in accuracy but at much greater cost. To overcome this computational bottleneck, herein we introduce and develop the QED-CC-in-QED-SCF projection-based embedding method that inherits all the favorable properties from the two worlds, computational efficiency and accuracy. The performance of the embedding method is assessed by studying some prototypical but relevant reactions, such as methyl transfer reaction, proton transfer reaction, as well as protonation reaction in a complex environment. The results obtained with the new embedding method are in excellent agreement with more expensive QED-CC results. The analysis performed on these reactions indicate that the electron-photon correlation effects are local in nature and that only a small region should be treated at the QED-CC level for capturing important effects due to cavity. This work sets the stage for future developments of polaritonic quantum chemistry methods and it will serve as a guideline for development of other polaritonic embedding models.

Simulating core electron binding energies of halogenated species adsorbed on ice surfaces and in solution via relativistic quantum embedding calculations

In this work, we investigate the effects of the environment on the X-ray photoelectron spectra of hydrogen chloride and chloride ions adsorbed on ice surfaces, as well as of chloride ions in water droplets. In our approach, we combine a density functional theory (DFT) description of the ice surface with that of halogen species using the recently developed relativistic core-valence separation equation of motion coupled cluster (CVS-EOM-IP-CCSD) via the frozen density embedding formalism (FDE), to determine the K and L_1,2,3 edges of chlorine. Our calculations, which incorporate temperature effects through snapshots from classical molecular dynamics simulations, are shown to reproduce experimental trends in the change of the core binding energies of Cl^- upon moving from a liquid (water droplets) to an interfacial (ice quasi-liquid layer) environment. Our simulations yield water valence band binding energies in good agreement with experiment, which vary little between the droplets and the ice surface. For halide core binding energies there is an overall trend for overestimating experimental values, though good agreement between theory and experiment is found for Cl^- in water droplets and on ice. For HCl on the other hand there are significant discrepancies between experimental and calculated core binding energies when we consider structural models that maintain the H-Cl bond more or less intact. An analysis of models that allow for pre-dissociated and dissociated structures suggests that experimentally observed chemical shifts in binding energies between Cl^- and HCl would require that H⁺ (in the form of H₃O⁺) and Cl^- are separated by roughly 4-6 Å.

Theory of electronic resonances: fundamental aspects and recent advances

Electronic resonances are states that are unstable towards loss of electrons. They play critical roles in high-energy environments across chemistry, physics, and biology but are also relevant to processes under ambient conditions that involve unbound electrons. This feature article focuses on complex-variable techniques such as complex scaling and complex absorbing potentials that afford a treatment of electronic resonances in terms of discrete square-integrable eigenstates of non-Hermitian Hamiltonians with complex energy. Fundamental aspects of these techniques as well as their integration into molecular electronic-structure theory are discussed and an overview of some recent developments is given: analytic gradient theory for electronic resonances, the application of rank-reduction techniques and quantum embedding to them, as well as approaches for evaluating partial decay widths.

Molecular Auger Decay Rates from Complex-Variable Coupled-Cluster Theory

The emission of an Auger electron is the predominant relaxation mechanism of core-vacant states in molecules composed of light nuclei. In this non-radiative decay process, one valence electron fills the core vacancy while a second valence electron is emitted into the ionization continuum. Because of this coupling to the continuum, core-vacant states represent electronic resonances that can be tackled with standard quantum-chemical methods only if they are approximated as bound states, meaning that Auger decay is neglected. Here, we present an approach to compute Auger decay rates of core-vacant states from coupled-cluster and equation-of-motion coupled-cluster wave functions combined with complex scaling of the Hamiltonian or, alternatively, complex-scaled basis functions. Through energy decomposition analysis, we illustrate how complex-scaled methods are capable of describing the coupling to the ionization continuum without the need to model the wave function of the Auger electron explicitly. In addition, we introduce in this work several approaches for the determination of partial decay widths and Auger branching ratios from complex-scaled coupled-cluster wave functions. We demonstrate the capabilities of our new approach by computations on core-ionized states of neon, water, dinitrogen, and benzene. Coupled-cluster and equation-of-motion coupled-cluster theory in the singles and doubles approximation both deliver excellent results for total decay widths, whereas we find partial widths more straightforward to evaluate with the former method. We also observe that the requirements towards the basis set are less arduous for Auger decay than for other types of resonances so that extensions to larger molecules are readily possible.

Spin-Orbit Natural Transition Orbitals and Spin-Forbidden Transitions

Natural transition orbitals (NTOs) are in widespread use for visualizing and analyzing electronic transitions. The present work introduces the analysis of formally spin-forbidden transitions with the help of complex-valued spin-orbit (SO) NTOs. The analysis specifically focuses on the components in such transitions that cause their intensity to be nonzero because of SO coupling. Transition properties such as transition dipole moments are partitioned into SO-NTO hole-particle pairs, such that contributions to the intensity from specific occupied and unoccupied orbitals are obtained. The method has been implemented within the restricted active space (RAS) self-consistent field wave function theory framework, with SO coupling treated by RAS state interaction. SO-NTOs have a broad range of potential applications, which is illustrated by the T₂-S₁ state mixing in pyrazine, spin-forbidden versus spin-allowed 4f-5d transitions in the Tb³⁺ ion, and the phosphorescence of tris(2-phenylpyridine) iridium [Ir(ppy)₃].

Variational Solutions for Resonances by a Finite-Difference Grid Method

We demonstrate that the finite difference grid method (FDM) can be simply modified to satisfy the variational principle and enable calculations of both real and complex poles of the scattering matrix. These complex poles are known as resonances and provide the energies and inverse lifetimes of the system under study (e.g., molecules) in metastable states. This approach allows incorporating finite grid methods in the study of resonance phenomena in chemistry. Possible applications include the calculation of electronic autoionization resonances which occur when ionization takes place as the bond lengths of the molecule are varied. Alternatively, the method can be applied to calculate nuclear predissociation resonances which are associated with activated complexes with finite lifetimes.