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.

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.

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.

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.

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.

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.

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.

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.

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. 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.

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.

Effect of Nonplanarity on Excited-State Proton Transfer and Internal Conversion in Salicylideneaniline.

Salicylideneaniline (SA) is a prototype for excited-state intramolecular proton transfer (ESIPT) reactions in nonplanar molecules. It is generally understood that the dominant photochemical pathway in this molecule is ESIPT followed by nonradiative decay due to twisting about its phenolic bond. However, the presence of a secondary internal conversion pathway resulting from frustrated proton transfer remains a matter of contention. We perform a detailed nonadiabatic dynamics simulation of SA and definitively identify the existence of both reaction pathways, thereby showing the presence of a secondary photochemical pathway and providing insight into the nature of ESIPT dynamics in molecules with nonplanar ground-state geometries.

Atomistic non-adiabatic dynamics of the LH2 complex with a GPU-accelerated ab initio exciton model.

We recently outlined an efficient multi-tiered parallel ab initio excitonic framework that utilizes time dependent density functional theory (TDDFT) to calculate ground and excited state energies and gradients of large supramolecular complexes in atomistic detail - enabling us to undertake non-adiabatic simulations which explicitly account for the coupled anharmonic vibrational motion of all the constituent atoms in a supramolecular system. Here we apply that framework to the 27 coupled bacterio-chlorophyll-a chromophores which make up the LH2 complex, using it to compute an on-the-fly nonadiabatic surface-hopping (SH) trajectory of electronically excited LH2. Part one of this article is focussed on calibrating our ab initio exciton Hamiltonian using two key parameters: a shift δ, which corrects for the error in TDDFT vertical excitation energies; and an effective dielectric constant ε, which describes the average screening of the transition-dipole coupling between chromophores. Using snapshots obtained from equilibrium molecular dynamics simulations (MD) of LH2, we tune the values of both δ and ε through fitting to the thermally broadened experimental absorption spectrum, giving a linear absorption spectrum that agrees reasonably well with experiment. In part two of this article, we construct a time-resolved picture of the coupled vibrational and excitation energy transfer (EET) dynamics in the sub-picosecond regime following photo-excitation. Assuming Franck-Condon excitation of a narrow eigenstate band centred at 800 nm, we use surface hopping to follow a single nonadiabatic dynamics trajectory within the full eigenstate manifold. Consistent with experimental data, this trajectory gives timescales for B800→B850 population transfer (τB800→B850) between 650-1050 fs, and B800 population decay (τ800→) between 10-50 fs. The dynamical picture that emerges is one of rapidly fluctuating LH2 eigenstates that are delocalized over multiple chromophores and undergo frequent crossing on a femtosecond timescale as a result of the atomic vibrations of the constituent chromophores. The eigenstate fluctuations arise from disorder that is driven by vibrational dynamics with multiple characteristic timescales. The scalability of our ab initio excitonic computational framework across massively parallel architectures opens up the possibility of addressing a wide range of questions, including how specific dynamical motions impact both the pathways and efficiency of electronic energy-transfer within large supramolecular systems.

Excited-State Dynamics of a Benzotriazole Photostabilizer: 2-(2'-Hydroxy-5'-methylphenyl)benzotriazole.

A large number of common photostabilizers are based on the 2-(2'-hydroxyphenyl)benzotriazole structure. One common example is 2-(2'-hydroxy-5'-methylphenyl)benzotriazole, or TINUVIN-P. The excited-state dynamics of this molecule have been extensively characterized by ultrafast spectroscopies. These experiments have established that upon photoexcitation TINUVIN-P exhibits excited-state proton transfer followed by a remarkably fast internal conversion. We simulate the excited-state dynamics using ab initio multiple spawning (AIMS) and a complete active space configuration interaction (CASCI) wave function with a correction from density functional theory (DFT) to generate the potential energy surfaces. We predict ultrafast proton transfer on the order of 20 fs followed by simultaneous twisting and pyramidalization until a seam of conical intersection is reached. Near the intersection seam population transfer to the ground state is highly efficient. The process is best described as ballistic wavepacket motion from the Franck-Condon point along a barrierless coordinate leading to the seam of intersection. Internal conversion is primarily mediated by a minimum-energy conical intersection (MECI) with a high degree of pyramidalization. We posit that the presence of a nitrogen atom in the bond linking the phenyl to the benzotriazole allows for the rapid pyramidalization and the short excited-state lifetime.

Excited-State Dynamics of 2-(2'-Hydroxyphenyl)benzothiazole: Ultrafast Proton Transfer and Internal Conversion.

One of the most widely studied model systems for excited-state proton transfer (ESPT) is the 2-(2'-hydroxyphenyl)benzothiazole (HBT) molecule. This compound undergoes ultrafast ESPT followed by internal conversion to return to the ground state. In the present work, we simulate the nonadiabatic photochemistry of HBT using ab initio multiple spawning (AIMS) nuclear dynamics and a complete active space configuration interaction (CASCI) method in conjunction with wave function-in-DFT embedding to obtain ground- and excited-state potential surfaces on-the-fly. Our simulation predicts ultrafast ESPT with a time constant of 48-54 fs and an excited-state lifetime of 1.7-1.8 ps. Following proton transfer, HBT becomes trapped in a metastable keto structure on the S1 state. Eventually, the molecule begins to twist and proceeds toward a seam of intersection with the ground state where internal conversion is highly efficient.

Complete active space configuration interaction from state-averaged configuration interaction singles natural orbitals: Analytic first derivatives and derivative coupling vectors.

A new complete active space configuration interaction (CASCI) method was recently introduced that uses state-averaged natural orbitals from the configuration interaction singles method (configuration interaction singles natural orbital CASCI, CISNO-CASCI). This method has been shown to perform as well or better than state-averaged complete active space self-consistent field for a variety of systems. However, further development and testing of this method have been limited by the lack of available analytic first derivatives of the CISNO-CASCI energy as well as the derivative coupling between electronic states. In the present work, we present a Lagrangian-based formulation of these derivatives as well as a highly efficient implementation of the resulting equations accelerated with graphical processing units. We demonstrate that the CISNO-CASCI method is practical for dynamical simulations of photochemical processes in molecular systems containing hundreds of atoms.

A direct-compatible formulation of the coupled perturbed complete active space self-consistent field equations on graphical processing units.

We recently developed an algorithm to compute response properties for the state-averaged complete active space self-consistent field method (SA-CASSCF) that capitalized on sparsity in the atomic orbital basis. Our original algorithm was limited to treating small to moderate sized active spaces, but the recent development of graphical processing unit (GPU) based direct-configuration interaction algorithms provides an opportunity to extend this to large active spaces. We present here a direct-compatible version of the coupled perturbed equations, enabling us to compute response properties for systems treated with arbitrary active spaces (subject to available memory and computation time). This work demonstrates that the computationally demanding portions of the SA-CASSCF method can be formulated in terms of seven fundamental operations, including Coulomb and exchange matrix builds and their derivatives, as well as, generalized one- and two-particle density matrix and σ vector constructions. As in our previous work, this algorithm exhibits low computational scaling and is accelerated by the use of GPUs, making possible optimizations and nonadiabatic dynamics on systems with O(1000) basis functions and O(100) atoms, respectively.