Simulating the Excited-State Dynamics of Polaritons with Ab Initio Multiple Spawning.

Over the past decade, there has been a growth of interest in polaritonic chemistry, where the formation of hybrid light-matter states (polaritons) can alter the course of photochemical reactions. These hybrid states are created by strong coupling between molecules and photons in resonant optical cavities and can even occur in the absence of light when the molecule is strongly coupled with the electromagnetic fluctuations of the vacuum field. We present a first-principles model to simulate nonadiabatic dynamics of such polaritonic states inside optical cavities by leveraging graphical processing units (GPUs). Our first implementation of this model is specialized for a single molecule coupled to a single-photon mode confined inside the optical cavity but with any number of excited states computed using complete active space configuration interaction (CASCI) and a Jaynes-Cummings-type Hamiltonian. Using this model, we have simulated the excited-state dynamics of a single salicylideneaniline (SA) molecule strongly coupled to a cavity photon with the ab initio multiple spawning (AIMS) method. We demonstrate how the branching ratios of the photodeactivation pathways for this molecule can be manipulated by coupling to the cavity. We also show how one can stop the photoreaction from happening inside of an optical cavity. Finally, we also investigate cavity-based control of the ordering of two excited states (one optically bright and the other optically dark) inside a cavity for a set of molecules, where the dark and bright states are close in energy.

Ultrafast internal conversion and photochromism in gas-phase salicylideneaniline.

Salicylideneaniline (SA) is an archetypal system for excited-state intramolecular proton transfer (ESIPT) in non-planar systems. Multiple channels for relaxation involving both the keto and enol forms have been proposed after excitation to S1 with near-UV light. Here, we present transient absorption measurements of hot gas-phase SA, jet-cooled SA, and SA in Ar clusters using cavity-enhanced transient absorption spectroscopy (CE-TAS). Assignment of the spectra is aided by simulated TAS spectra, computed by applying time-dependent complete active space configuration interaction (TD-CASCI) to structures drawn from nonadiabatic molecular dynamics simulations. We find prompt ESIPT in all conditions followed by the rapid generation of the trans keto metastable photochrome state and fluorescent keto state in parallel. Increasing the internal energy increases the photochrome yield and decreases the fluorescent yield and fluorescent state lifetime observed in TAS. In Ar clusters, internal conversion of SA is severely hindered, but the photochrome yield is unchanged. Taken together, these results are consistent with the photochrome being produced via the vibrationally excited keto population after ESIPT.

Efficient quantum analytic nuclear gradients with double factorization.

Efficient representations of the Hamiltonian, such as double factorization, drastically reduce the circuit depth or the number of repetitions in error corrected and noisy intermediate-scale quantum (NISQ) algorithms for chemistry. We report a Lagrangian-based approach for evaluating relaxed one- and two-particle reduced density matrices from double factorized Hamiltonians, unlocking efficiency improvements in computing the nuclear gradient and related derivative properties. We demonstrate the accuracy and feasibility of our Lagrangian-based approach to recover all off-diagonal density matrix elements in classically simulated examples with up to 327 quantum and 18 470 total atoms in QM/MM simulations with modest-sized quantum active spaces. We show this in the context of the variational quantum eigensolver in case studies, such as transition state optimization, ab initio molecular dynamics simulation, and energy minimization of large molecular systems.

TeraChem protocol buffers (TCPB): Accelerating QM and QM/MM simulations with a client-server model.

The routine use of electronic structures in many chemical simulation applications calls for efficient and easy ways to access electronic structure programs. We describe how the graphics processing unit (GPU) accelerated electronic structure program TeraChem can be set up as an electronic structure server, to be easily accessed by third-party client programs. We exploit Google's protocol buffer framework for data serialization and communication. The client interface, called TeraChem protocol buffers (TCPB), has been designed for ease of use and compatibility with multiple programming languages, such as C++, Fortran, and Python. To demonstrate the ease of coupling third-party programs with electronic structures using TCPB, we have incorporated the TCPB client into Amber for quantum mechanics/molecular mechanics (QM/MM) simulations. The TCPB interface saves time with GPU initialization and I/O operations, achieving a speedup of more than 2× compared to a prior file-based implementation for a QM region with ∼250 basis functions. We demonstrate the practical application of TCPB by computing the free energy profile of p-hydroxybenzylidene-2,3-dimethylimidazolinone (p-HBDI-)-a model chromophore in green fluorescent proteins-on the first excited singlet state using Hamiltonian replica exchange for enhanced sampling. All calculations in this work have been performed with the non-commercial freely-available version of TeraChem, which is sufficient for many QM region sizes in common use.

InteraChem: Exploring Excited States in Virtual Reality with Ab Initio Interactive Molecular Dynamics.

InteraChem is an ab initio interactive molecular dynamics (AI-IMD) visualizer that leverages recent advances in virtual reality hardware and software, as well as the graphical processing unit (GPU)-accelerated TeraChem electronic structure package, in order to render quantum chemistry in real time. We introduce the exploration of electronically excited states via AI-IMD using the floating occupation molecular orbital-complete active space configuration interaction method. The optimization tools in InteraChem enable identification of excited state minima as well as minimum energy conical intersections for further characterization of excited state chemistry in small- to medium-sized systems. We demonstrate that finite-temperature Hartree-Fock theory is an efficient method to perform ground state AI-IMD. InteraChem allows users to track electronic properties such as molecular orbitals and bond order in real time, resulting in an interactive visualization tool that aids in the interpretation of excited state chemistry data and makes quantum chemistry more accessible for both research and educational purposes.

Rank-reduced coupled-cluster. III. Tensor hypercontraction of the doubles amplitudes.

We develop a quartic-scaling implementation of coupled-cluster singles and doubles (CCSD) based on low-rank tensor hypercontraction (THC) factorizations of both the electron repulsion integrals (ERIs) and the doubles amplitudes. This extends our rank-reduced (RR) coupled-cluster method to incorporate higher-order tensor factorizations. The THC factorization of the doubles amplitudes accounts for most of the gain in computational efficiency as it is sufficient, in conjunction with a Cholesky decomposition of the ERIs, to reduce the computational complexity of most contributions to the CCSD amplitude equations. Further THC factorization of the ERIs reduces the complexity of certain terms arising from nested commutators between the doubles excitation operator and the two-electron operator. We implement this new algorithm using graphical processing units and demonstrate that it enables CCSD calculations for molecules with 250 atoms and 2500 basis functions using a single computer node. Furthermore, we show that the new method computes correlation energies with comparable accuracy to the underlying RR-CCSD method.

GPU acceleration of rank-reduced coupled-cluster singles and doubles.

We have developed a graphical processing unit (GPU) accelerated implementation of our recently introduced rank-reduced coupled-cluster singles and doubles (RR-CCSD) method. RR-CCSD introduces a low-rank approximation of the doubles amplitudes. This is combined with a low-rank approximation of the electron repulsion integrals via Cholesky decomposition. The result of these two low-rank approximations is the replacement of the usual fourth-order CCSD tensors with products of second- and third-order tensors. In our implementation, only a single fourth-order tensor must be constructed as an intermediate during the solution of the amplitude equations. Owing in large part to the compression of the doubles amplitudes, the GPU-accelerated implementation shows excellent parallel efficiency (95% on eight GPUs). Our implementation can solve the RR-CCSD equations for up to 400 electrons and 1550 basis functions-roughly 50% larger than the largest canonical CCSD computations that have been performed on any hardware. In addition to increased scalability, the RR-CCSD computations are faster than the corresponding CCSD computations for all but the smallest molecules. We test the accuracy of RR-CCSD for a variety of chemical systems including up to 1000 basis functions and determine that accuracy to better than 0.1% error in the correlation energy can be achieved with roughly 95% compression of the ov space for the largest systems considered. We also demonstrate that conformational energies can be predicted to be within 0.1 kcal mol-1 with efficient compression applied to the wavefunction. Finally, we find that low-rank approximations of the CCSD doubles amplitudes used in the similarity transformation of the Hamiltonian prior to a conventional equation-of-motion CCSD computation will not introduce significant errors (on the order of a few hundredths of an electronvolt) into the resulting excitation energies.

Predictions of Pre-edge Features in Time-Resolved Near-Edge X-ray Absorption Fine Structure Spectroscopy from Hole-Hole Tamm-Dancoff-Approximated Density Functional Theory.

Time-resolved near-edge X-ray absorption fine structure (TR-NEXAFS) spectroscopy is a powerful technique for studying photochemical reaction dynamics with femtosecond time resolution. In order to avoid ambiguity in TR-NEXAFS spectra from nonadiabatic dynamics simulations, core- and valence-excited states must be evaluated on equal footing and those valence states must also define the potential energy surfaces used in the nonadiabatic dynamics simulation. In this work, we demonstrate that hole-hole Tamm-Dancoff-approximated density functional theory (hh-TDA) is capable of directly simulating TR-NEXAFS spectroscopies. We apply hh-TDA to the excited-state dynamics of acrolein. We identify two pre-edge features in the oxygen K-edge TR-NEXAFS spectrum associated with the S2 (ππ*) and S1 (nπ*) excited states. We show that these features can be used to follow the internal conversion dynamics between the lowest three electronic states of acrolein. Due to the low, O(N2) apparent computational complexity of hh-TDA and our GPU-accelerated implementation, this method is promising for the simulation of pre-edge features in TR-NEXAFS spectra of large molecules and molecules in the condensed phase.

Nonadiabatic Dynamics Simulation of the Wavelength-Dependent Photochemistry of Azobenzene Excited to the nπ* and ππ* Excited States.

Azobenzene is one of the most ubiquitous photoswitches in photochemistry and a prototypical model for photoisomerizing systems. Despite this, its wavelength-dependent photochemistry has puzzled researchers for decades. Upon excitation to the higher energy ππ* excited state instead of the dipole-forbidden nπ* state, the quantum yield of isomerization from trans- to cis-azobenzene is halved. The difficulties associated with unambiguously resolving this effect both experimentally and theoretically have contributed to lasting controversies regarding the photochemistry of azobenzene. Here, we systematically characterize the dynamic photoreaction pathways of azobenzene by performing first-principles simulations of the nonadiabatic dynamics following excitation to both the ππ* and the nπ* states. We demonstrate that ground-state recovery is mediated by two distinct S1 decay pathways: a reactive twisting pathway and an unreactive planar pathway. Increased preference for the unreactive pathway upon ππ* excitation largely accounts for the wavelength-dependent behavior observed in azobenzene.

Ab Initio Nonadiabatic Molecular Dynamics with Hole-Hole Tamm-Dancoff Approximated Density Functional Theory.

The study of photoinduced dynamics in chemical systems necessitates accurate and computationally efficient electronic structure methods, especially as the systems of interest grow larger. The linear response hole-hole Tamm-Dancoff approximated (hh-TDA) density functional theory method was recently proposed to satisfy such demands. The N-electron electronic states are obtained by means of double annihilations on a doubly anionic (N + 2)-electron reference state, allowing for the ground and excited states to be formed on the same footing and thus enabling the correct description of conical intersections. Dynamic electron correlation effects are incorporated by means of the exchange-correlation functional. The accuracy afforded by the simultaneous treatment of static and dynamic correlation in addition to the relatively low computational cost, comparable to that of time-dependent density functional theory (TDDFT), makes it a promising ab initio electronic structure method for on-the-fly generation of potential energy surfaces in nonadiabatic dynamics simulations of photochemical systems, particularly those for which the nπ* and ππ* electronic excitations are most relevant. Here, we apply the hh-TDA method to nonadiabatic dynamics simulations of prototypical photochemical processes. First, we demonstrate the ability of hh-TDA to adequately describe conical intersection geometries. We next examine its ability to describe the ultrafast excited state dynamics of photoexcited ethylene through an ab initio multiple spawning (AIMS) dynamics simulation. Finally, we present an alternative variant of the hh-TDA method, which uses orbitals from a fractional occupation number Kohn-Sham (FON-KS) calculation applied to an ensemble with N-electrons. The resulting method is termed floating occupation molecular orbital hh-TDA (FOMO-hh-TDA). This scheme allows us to combine hh-TDA with global hybrid functionals and allows us to avoid unbound valence orbitals that may result from an (N + 2)-electron self-consistent field (SCF) procedure. FOMO-hh-TDA-BHLYP faithfully reproduces the nonadiabatic dynamics of trans-azobenzene (TAB) and is used to characterize the excited state decay pathways from the first (nπ*) excited state.

Hole-hole Tamm-Dancoff-approximated density functional theory: A highly efficient electronic structure method incorporating dynamic and static correlation.

The study of photochemical reaction dynamics requires accurate as well as computationally efficient electronic structure methods for the ground and excited states. While time-dependent density functional theory (TDDFT) is not able to capture static correlation, complete active space self-consistent field methods neglect much of the dynamic correlation. Hence, inexpensive methods that encompass both static and dynamic electron correlation effects are of high interest. Here, we revisit hole-hole Tamm-Dancoff approximated (hh-TDA) density functional theory for this purpose. The hh-TDA method is the hole-hole counterpart to the more established particle-particle TDA (pp-TDA) method, both of which are derived from the particle-particle random phase approximation (pp-RPA). In hh-TDA, the N-electron electronic states are obtained through double annihilations starting from a doubly anionic (N+2 electron) reference state. In this way, hh-TDA treats ground and excited states on equal footing, thus allowing for conical intersections to be correctly described. The treatment of dynamic correlation is introduced through the use of commonly employed density functional approximations to the exchange-correlation potential. We show that hh-TDA is a promising candidate to efficiently treat the photochemistry of organic and biochemical systems that involve several low-lying excited states-particularly those with both low-lying ππ* and nπ* states where inclusion of dynamic correlation is essential to describe the relative energetics. In contrast to the existing literature on pp-TDA and pp-RPA, we employ a functional-dependent choice for the response kernel in pp- and hh-TDA, which closely resembles the response kernels occurring in linear response and collinear spin-flip TDDFT.

TeraChem: Accelerating electronic structure and ab initio molecular dynamics with graphical processing units.

Developed over the past decade, TeraChem is an electronic structure and ab initio molecular dynamics software package designed from the ground up to leverage graphics processing units (GPUs) to perform large-scale ground and excited state quantum chemistry calculations in the gas and the condensed phase. TeraChem's speed stems from the reformulation of conventional electronic structure theories in terms of a set of individually optimized high-performance electronic structure operations (e.g., Coulomb and exchange matrix builds, one- and two-particle density matrix builds) and rank-reduction techniques (e.g., tensor hypercontraction). Recent efforts have encapsulated these core operations and provided language-agnostic interfaces. This greatly increases the accessibility and flexibility of TeraChem as a platform to develop new electronic structure methods on GPUs and provides clear optimization targets for emerging parallel computing architectures.

Performance of Coupled-Cluster Singles and Doubles on Modern Stream Processing Architectures.

We develop a new implementation of coupled-cluster singles and doubles (CCSD) optimized for the most recent graphical processing unit (GPU) hardware. We find that a single node with 8 NVIDIA V100 GPUs is capable of performing CCSD computations on roughly 100 atoms and 1300 basis functions in less than 1 day. Comparisons against massively parallel implementations of CCSD suggest that more than 64 CPU-based nodes (each with 16 cores) are required to match this performance.

Strong, Nonresonant Radiation Enhances Cis-Trans Photoisomerization of Stilbene in Solution.

Previously, it has been demonstrated that external electric fields may be used to exert control over chemical reactivity. In this study, the impact of a strong, nonresonant IR field (1064 nm) on the photoisomerization of cis-stilbene is investigated in cyclohexane solution. The design of a suitable reaction vessel for characterization of this effect is presented. The electric field supplied by the pulsed, near-IR radiation (εl = 4.5 × 107 V/cm) enhances the cis → trans photoisomerization yield at the red edge of the absorption spectrum (wavelengths between 337 and 340 nm). Within the microliter focal volume, up to 75% of all cis-stilbene molecules undergo isomerization to trans-stilbene in the strong electric-field environment, indicating a significant increase relative to the 35% yield of trans-stilbene under field-free conditions. This result correlates with a 1-3% enhancement in the trans-stilbene concentration throughout the bulk solution. Theoretical analysis suggests that the observed change is the result of dynamic Stark shifting of the ground and first excited states, leading to a significant redshift in cis-stilbene's absorption spectrum. The predicted increase in the absorption cross section in this range of excitation wavelengths is qualitatively consistent with the experimental increase in trans-stilbene production.

Nonadiabatic Dynamics of Photoexcited cis-Stilbene Using Ab Initio Multiple Spawning.

The photochemistry of cis-stilbene proceeds through two pathways: cis-trans isomerization and ring closure to 4a,4b-dihydrophenanthrene (DHP). Despite serving for many decades as a model system for photoisomerization, the photodynamics of cis-stilbene is still not fully understood. We use ab initio multiple spawning on a SA-2-CASSCF(2,2) potential energy surface to simulate the nonadiabatic dynamics of isolated cis-stilbene. We find the cyclization (to DHP and cis-stilbene) and isomerization (to trans- and cis-stilbene) reaction coordinates to be orthogonal; branching between the two pathways is determined on the S1 excited state within 150 fs of photoexcitation. Trajectory basis functions (TBFs) undergoing cyclization decay rapidly to the ground state in 250 fs, while TBFs moving along the isomerization coordinate remain on the excited state longer, with the majority decaying between 300 and 500 fs. We observe three avoided crossing regions in the dynamics: two along the isomerization coordinate (displaying pyramidalization and migration of an ethylenic hydrogen or phenyl group), and one DHP-like conical intersection along the cyclization coordinate. The isomeric form of the vibrationally hot photoproducts (as determined by measurement 2 ps after photoexcitation) is determined within less than 50 fs of decay to the ground state mediated by passage through a conical intersection. Excess vibrational energy of ground state cis- and trans-stilbene is channelled into phenyl torsions (with mostly opposing directionality). Our simulations are validated by direct comparison to experiment for the absorption spectrum, branching ratio of the three photoproducts (44:52:4 cis-stilbene:trans-stilbene:DHP), and excited state lifetime (520 ± 40 fs).

Psi4 1.4: Open-source software for high-throughput quantum chemistry.

PSI4 is a free and open-source ab initio electronic structure program providing implementations of Hartree-Fock, density functional theory, many-body perturbation theory, configuration interaction, density cumulant theory, symmetry-adapted perturbation theory, and coupled-cluster theory. Most of the methods are quite efficient, thanks to density fitting and multi-core parallelism. The program is a hybrid of C++ and Python, and calculations may be run with very simple text files or using the Python API, facilitating post-processing and complex workflows; method developers also have access to most of PSI4's core functionalities via Python. Job specification may be passed using The Molecular Sciences Software Institute (MolSSI) QCSCHEMA data format, facilitating interoperability. A rewrite of our top-level computation driver, and concomitant adoption of the MolSSI QCARCHIVE INFRASTRUCTURE project, makes the latest version of PSI4 well suited to distributed computation of large numbers of independent tasks. The project has fostered the development of independent software components that may be reused in other quantum chemistry programs.

Rank reduced coupled cluster theory. II. Equation-of-motion coupled-cluster singles and doubles.

Equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) is a reliable and popular approach to the determination of electronic excitation energies. Recently, we have developed a rank-reduced CCSD (RR-CCSD) method that allows the ground-state coupled-cluster energy to be determined with low-rank cluster amplitudes. Here, we extend this approach to excited-state energies through a RR-EOM-CCSD method. We start from the EOM-CCSD energy functional and insert low-rank approximations to the doubles amplitudes. The result is an approximate EOM-CCSD method with only a quadratic number (in the molecular size) of free parameters in the wavefunction. Importantly, our formulation of RR-EOM-CCSD preserves the size intensivity of the excitation energy and size extensivity of the total energy. Numerical tests of the method suggest that accuracy on the order of 0.05-0.01 eV in the excitation energy is possible with 1% or less of the original number of wavefunction coefficients; accuracy of better than 0.01 eV can be achieved with about 4% or less of the free parameters. The amount of compression at a given accuracy level is expected to increase with the size of the molecule. The RR-EOM-CCSD method is a new path toward the efficient determination of accurate electronic excitation energies.

Quantum Computation of Electronic Transitions Using a Variational Quantum Eigensolver.

We develop an extension of the variational quantum eigensolver (VQE) algorithm-multistate contracted VQE (MC-VQE)-that allows for the efficient computation of the transition energies between the ground state and several low-lying excited states of a molecule, as well as the oscillator strengths associated with these transitions. We numerically simulate MC-VQE by computing the absorption spectrum of an ab initio exciton model of an 18-chromophore light-harvesting complex from purple photosynthetic bacteria.

Rank reduced coupled cluster theory. I. Ground state energies and wavefunctions.

We propose a compression of the opposite-spin coupled cluster doubles amplitudes of the form τij ab≡Uia VTVWUjb W, where Uia V are the nV-highest magnitude eigenvectors of the MP2 or MP3 doubles amplitudes. Together with a corresponding parameterization of the opposite-spin coupled cluster Lagrange multipliers of the form λab ij≡Uia VLVWUjb W, this yields a fully self-consistent parameterization of reduced-rank coupled cluster equations in terms of the Lagrangian L0TVW,LVW. Making this Lagrangian stationary with respect to the LVW parameters yields a perfectly determined set of equations for the TVW equations and coupled cluster energy. These equations can be solved using a Lyapunov equation for the first-order amplitude updates. We test this "rank-reduced coupled cluster" method for coupled cluster singles and doubles in medium sized molecules and find that substantial compression of the T^2 amplitudes is possible with acceptable accuracy.

Stereoisomer specific reaction of hexabromocyclododecane with reduced sulfur species in aqueous solutions.

The individual degradation rates of the three dominant stereoisomers (α, ß, γ) of hexabromocyclododecane (HBCDD) with bisulfide and polysulfides were investigated at pH 9 in methanol/water solutions at two different temperatures (25 °C and 40 °C). Under all conditions investigated, α-HBCDD reacts 10 to 20 times slower with bisulfide than ß-HBCDD and γ-HBCDD. The difference in reactivity of HBCDD isomers can be explained by the different populations of stable conformers with large dihedral angle between the vicinal bromine atoms. It was also observed that the reaction of HBCDD with polysulfides is about six times faster than with bisulfide. The experiments performed in solvent mixtures with increased water content at 40 °C indicated that the reaction of HBCDD with bisulfide is faster with higher percentage of water. The much slower abiotic reaction of α-HBCDD compared to ß-HBCDD and γ-HBCDD could potentially contribute to the fact that α-HBCDD is more persistent in the environment than γ-HBCDD. Only one isomer of tetrabromocyclododecene (TBCDe-5) was identified as a degradation product of the reaction of HBCDD with reduced sulfur species. TBCDe-5 itself reacts about ten times slower with bisulfide and twenty times slower with polysulfide than HBCDD. The study demonstrates that polysulfides and bisulfides can reduce HBCDD sufficiently in natural anoxic environments and the dominant pathway for the degradation of HBCDD by reduced sulfur species is very likely to be the reductive debromination of vicinal dibromides via concerted anti-elimination.