Development of Parallel On-the-Fly Crystal Algorithm for Global Exploration of Conical Intersection Seam Space.

Conical intersection (CI) seams are configuration spaces of a molecular system where two or more (spin) adiabatic electronic states are degenerate in energy. They play essential roles in photochemistry because nonradiative decays often occur near the minima of the seam, i.e., the minimum energy CIs (MECIs). Thus, it is important to explore the CI seams and discover the MECIs. Although various approaches exist for CI seam exploration, most of them are local in nature, requiring reasonable initial guesses of geometries and nuclear gradients during the search. Global search algorithms, on the other hand, are powerful because they can fully sample the configurational space and locate important MECIs missed by local algorithms. However, global algorithms are often computationally expensive for large systems due to their poor scalability with respect to the number of degrees of freedom. To overcome this challenge, we develop the parallel on-the-fly Crystal algorithm to globally explore the CI seam space, taking advantage of its superior scaling behavior. Specifically, Crystal is coupled with on-the-fly evaluations of the excited and ground state energies using multireference electronic structure methods. Meanwhile, the algorithm is parallelized to further boost its computational efficiency. The effectiveness of this new algorithm is tested for three types of molecular photoswitches of significant importance in material and biomedical sciences: photostatin (PST), stilbene, and butadiene. A rudimentary implementation of the algorithm is applied to PST and stilbene, resulting in the discovery of all previously identified MECIs and several new ones. A refined version of the algorithm, combined with a systematic clustering technique, is applied to butadiene, resulting in the identification of an unprecedented number of energetically accessible MECIs. The results demonstrate that the parallel on-the-fly Crystal algorithm is a powerful tool for automated global CI seam exploration.

"On-the-fly" Crystal : How to reliably and automatically characterize and construct potential energy surfaces.

In this work, the Crystal  code, developed previously by the authors to find "holes" as well as legitimate transition states in existing potential energy surface (PES) functions [JPC Lett. 11, 6468 (2020)], is retooled to perform on-the-fly "direct dynamics"-type PES explorations, as well as automatic construction of new PES functions. In all of these contexts, the chief advantage of Crystal  over other methods is its ability to globally map the PES, thereby determining the most relevant regions of configuration space quickly and reliably-even when the dimensionality is rather large. Here, Crystal  is used to generate a uniformly spaced grid of density functional theory (DFT) or ab initio points, truncated over the relevant regions, which can then be used to either: (a) hone in precisely on PES features such as minima and transition states, or; (b) create a new PES function automatically, via interpolation. Proof of concept is demonstrated via application to three molecular systems: water (H 2 O), (reduced-dimensional) methane (CH 4 ), and methylene imine (CH 2 NH).

Dwell Times, Wavepacket Dynamics, and Quantum Trajectories for Particles with Spin 1/2.

The theoretical connections between quantum trajectories and quantum dwell times, previously explored in the context of 1D time-independent stationary scattering applications, are here generalized for multidimensional time-dependent wavepacket applications for particles with spin 1/2. In addition to dwell times, trajectory-based dwell time distributions are also developed, and compared with previous distributions based on the dwell time operator and the flux-flux correlation function. Dwell time distributions are of interest, in part because they may be of experimental relevance. In addition to standard unipolar quantum trajectories, bipolar quantum trajectories are also considered, and found to relate more directly to the dwell time (and other quantum time) quantities of greatest relevance for scattering applications. Detailed calculations are performed for a benchmark 3D spin-1/2 particle application, considered previously in the context of computing quantum arrival times.

Quantum Dynamical Investigation of Dihydrogen-Hydride Exchange in a Transition-Metal Polyhydride Complex.

The ongoing shift toward clean, sustainable energy is a primary driving force behind hydrogen fuel research. Safe and effective storage of hydrogen is a major challenge (particularly for mobile applications) and requires a detailed understanding of the atomic level interactions of hydrogen with its host materials. The light mass of hydrogen, however, implies that quantum effects are important, so a quantum dynamical treatment is required to properly account for these effects in computational simulations. As one such example, we describe herein the hydrogen exchange dynamics between a hydride and a dihydrogen ligand in the [FeH(H2)(PH3)4]+ model complex. A global three-dimensional (3D) potential energy surface (PES) was constructed by fitting to and interpolating from a discrete set of grid points computed using density functional theory; exact quantum dynamical calculations were then carried out on the 3D PES using discrete variable representation basis sets. Energy levels and their quantum tunneling splittings were computed up to 3000 cm-1 above the ground state. Within that energy range, all three fundamentals have been identified using wave function plots, as well as the first three overtones of the exchange (reaction coordinate) motion and several of its combination bands. From the tunneling splittings, the Boltzmann-averaged tunneling rates were computed. The Arrhenius plot of the total exchange rate shows a clear transition around 150 K, below which the activation energy is essentially zero and above which it is less than half of the electronic structure barrier. This indicates that exchange rates are governed by quantum tunneling throughout the relevant temperature range with the low-temperature regime dominated by a single quantum (ground) state. This work is the first-ever fully quantum dynamical study to investigate the hydrogen exchange dynamics between hydride and dihydrogen ligands coordinated to a transition-metal complex.

Adiabatic and Nonadiabatic Dynamics with Interacting Quantum Trajectories.

We present a quantum dynamics method based on the propagation of interacting quantum trajectories to describe both adiabatic and nonadiabatic processes within the same formalism. The idea originates from the work of Poirier [Chem. Phys.2010,370, 4-14] and Schiff and Poirier [J. Chem. Phys.2012,136, 031102] on quantum dynamics without wavefunctions. It consists of determining the quantum force arising in the Bohmian hydrodynamic formulation of quantum dynamics using only information about quantum trajectories. The particular time-dependent propagation scheme proposed here results in very stable dynamics. Its performance is discussed by applying the method to analytical potentials in the adiabatic regime, and by combining it with the exact factorization method in the nonadiabatic regime.

Full-dimensional Schrödinger wavefunction calculations using tensors and quantum computers: the Cartesian component-separated approach.

Traditional methods in quantum chemistry rely on Hartree-Fock-based Slater-determinant (SD) representations, whose underlying zeroth-order picture assumes separability by particle. Here, we explore a radically different approach, based on separability by Cartesian component, rather than by particle [J. Jerke and B. Poirier, J. Chem. Phys., 2018, 148, 104101]. The approach appears to be very well suited for 3D grid-based methods in quantum chemistry, and thereby also for so-called "first-quantized" quantum computing. We first present an overview of the approach as implemented on classical computers, including numerical results that justify performance claims. In particular, we perform numerical calculations with four explicit electrons that are equivalent to full-CI matrix diagonalization with nearly 1015 SDs. We then present an implementation for quantum computers for which the number of quantum gates (and to a lesser extent, the number of qubits) can be dramatically reduced, in comparison with other quantum circuitry that has been envisioned for implementing first-quantized "quantum computational chemistry" (QCC).

Hitting the Trifecta: How to Simultaneously Push the Limits of Schrödinger Solution with Respect to System Size, Convergence Accuracy, and Number of Computed States.

Methods for solving the Schrödinger equation without approximation are in high demand but are notoriously computationally expensive. In practical terms, there are just three primary factors that currently limit what can be achieved: 1) system size/dimensionality; 2) energy level excitation; and 3) numerical convergence accuracy. Broadly speaking, current methods can deliver on any two of these three goals, but achieving all three at once remains an enormous challenge. In this paper, we shall demonstrate how to "hit the trifecta" in the context of molecular vibrational spectroscopy calculations. In particular, we compute the lowest 1000 vibrational states for the six-atom acetonitrile molecule (CH3CN), to a numerical convergence of accuracy 10-2 cm-1 or better. These calculations encompass all vibrational states throughout most of the dynamically relevant range (i.e., up to ∼4250 cm-1 above the ground state), computed in full quantum dimensionality (12 dimensions), to near spectroscopic accuracy. To our knowledge, no such vibrational spectroscopy calculation has ever previously been performed.

On neglecting Coriolis and related couplings in first-principles rovibrational spectroscopy: Considerations of symmetry, accuracy, and simplicity. II. Case studies for H2O isotopologues, H3+, O3, and NH3.

For centuries, it has been known that vibrational and rotational degrees of freedom are in general not separable. Nevertheless, surprisingly little is known about the best strategies for approximately separating these degrees of freedom in practice-even in the case of semirigid molecules, where the separation is most meaningful. There is also some confusion in the literature about the proper way to quantify the magnitude of the Coriolis (i.e., rotation-vibration) coupling in rovibrational Hamiltonians or its effect on the rovibrational eigenenergies. In this study, a vibrational-coordinate-independent metric is proposed to quantify the magnitude of the Coriolis contribution to the rovibrational Hamiltonian. The impact of Coriolis coupling on the rovibrational eigenenergies is computed numerically exactly, using both full and various truncated Hamiltonians. The role played by the choice of the vibrational coordinate system-and especially by the choice of "embedding" or body-fixed frame-is examined extensively, both numerically and analytically. This investigation targets several molecular prototypes, all of which serve as important benchmarks for the high-resolution spectroscopic community. Most of these are triatomic molecules, including water (H216O), its deuterated isotopologues (D216O and HD16O), H3+, and ozone (16O3), but the tetratomic ammonia molecule (14NH3) is also investigated. These studies provide important insight into the nature of Coriolis coupling under various circumstances. The findings of this study also have significant practical ramifications, vis-à-vis the use of simplifying numerical approximation techniques in nuclear-motion computations.

Plumbing Potentials for Molecules with Up To Tens of Atoms: How to Find Saddle Points and Fix Leaky Holes.

Potential energy surfaces (PESs) play an indispensable role in molecular dynamics but are notoriously difficult to flesh out properly in large-dimensional spaces. In particular, the undetected presence of PES holes, i.e., unphysical saddle points beyond which the potential energy drops arbitrarily, can have devastating effects on both classical and quantum dynamics calculations. In this study, the Crystal algorithm is developed as a tool for efficiently and accurately finding PES holes, as well as legitimate saddle points, even in very large-dimensional configuration spaces. The approach is applied to three large-dimensional PESs for molecular systems of current interest: uracil, naphthalene, and formic acid dimer. Low-lying PES holes are discovered and located for the first two systems-including naphthalene, for which no holes were previously suspected, to the best of our knowledge. Likewise, the double-well, double-proton-transfer isomerization saddle point for formic acid dimer is also located.

An algorithm to find (and plug) "holes" in multi-dimensional surfaces.

We have developed an algorithm to detect holes in multi-dimensional real-valued surfaces-such as the potential energy surfaces (PESs) that describe the nuclear motion of molecules in the context of the Born-Oppenheimer approximation. For our purposes, a PES "hole" is defined as an unphysical saddle point, beyond which the potential energy drops (typically) without limit to negative infinity. PES holes are numerical artifacts that can arise when fitting PES functional forms to discrete ab initio data-even when the data is of high quality, and/or for comparatively few degrees of freedom (DOF). Often undetected, PES holes can have devastating effects on subsequent dynamical calculations, especially if they occur at low energies. In this paper, we present a highly efficient algorithm designed to systematically identify hole configurations and energies. The method is applied to a variety of molecular PESs ranging up to 30 DOF. A number of evidently previously undetected PES holes are reported here-surprisingly, even for PESs that have been available for decades. The code itself (Crystal) is presented together with a user manual. These tools may be of great benefit for PES developers, who can use the information they provide to fix holes, once identified. More generally, the methodology can be applied in any context involving multi-dimensional surfaces.

On neglecting Coriolis and related couplings in first-principles rovibrational spectroscopy: considerations of symmetry, accuracy, and simplicity.

The rotation-vibration (Coriolis) coupling contribution to variationally computed rovibrational energy levels is investigated, employing triatomic AB[Formula: see text] molecules as models. In particular, calculations are performed for H[Formula: see text][Formula: see text]O, across a range of vibrational and rotational excitations, both with and without the Coriolis contribution. A variety of different embedding choices are considered, together with a hierarchy of increasingly severe approximations culminating in a generalized version of the so-called "centrifugal sudden" method. Several surprising and remarkable conclusions are found, including that the Eckart embedding is not the best embedding choice.

Exact bound rovibrational spectra of the neon tetramer.

Exact quantum dynamics calculations are performed for the bound rovibrational states of the neon tetramer (Ne4) in its ground electronic state, using pair-wise Lennard-Jones potentials and the ScalIT suite of parallel codes. The vibrational states separate into a low-lying group mostly localized to a single potential well and a higher-energy delocalized group lying above the isomerization threshold-with such a structure serving as a testament to the "intermediate" quantum nature of the Ne4 system. To accurately and efficiently represent both groups of states, the phase-space optimized discrete variable representation (PSO-DVR) approach was used, as implemented in the ScalIT code. The resultant 1D PSO effective potentials also shed significant light on the quantum dynamics of the system. All vibrational states were computed well up into the isomerization band and labeled up to the classical isomerization threshold-defined as the addition of the classical energy of a single bond, ϵ = 24.7 cm-1, to the quantum ground state energy. Rovibrational energy levels for all total angular momentum values in the range J = 1-5 were also computed, treating all Coriolis coupling exactly.

Using phase-space Gaussians to compute the vibrational states of OCHCO.

In this article, we use momentum-symmetrized phase-space Gaussians to calculate the vibrational energy eigenstates of the OCHCO+ cation. A potential energy surface provided by Bowman's group was used, albeit refit to a sixth order anharmonic force field. We have developed the "Crystal" algorithm to implement various basis set truncation strategies for our calculations. These calculations were performed using the SwitchBLADE code, designed by our group, which constructs and diagonalizes the vibrational Hamiltonian matrix in order to compute energy eigenstates. A theoretical overview of the Crystal algorithm is presented, as is a discussion of how best to truncate the basis for a floppy double-well system. Specific vibrational transition frequency results for OCHCO+ are also presented and analyzed.

Accurate characterization of the lowest triplet potential energy surface of SO2 with a coupled cluster method.

The near-equilibrium potential energy surface (PES) of the ã 3B1 state of SO2 is developed from explicitly correlated spin-unrestricted coupled cluster calculations with single, double, and perturbative triple excitations with an augmented triple-zeta correlation-consistent basis set. The lowest-lying ro-vibrational energy levels of several sulfur isotopologues have been determined using this PES. It is shown that the new ab initio PES provides a much better description of the low-lying vibrational states than a previous PES determined at the multi-reference configuration interaction level. In particular, the theory-experiment agreement for the three lowest-lying vibrational transitions is within 1-3 cm-1.

Communication: Adiabatic quantum trajectory capture for cold and ultra-cold chemical reactions.

The Langevin capture model is often used to describe barrierless reactive collisions. At very low temperatures, quantum effects may alter this simple capture image and dramatically affect the reaction probability. In this paper, we use the trajectory-ensemble reformulation of quantum mechanics, as recently proposed by one of the authors (Poirier) to compute adiabatic-channel capture probabilities and cross-sections for the highly exothermic reaction Li + CaH(v = 0, j = 0) → LiH + Ca, at low and ultra-low temperatures. Each captured quantum trajectory takes full account of tunneling and quantum reflection along the radial collision coordinate. Our approach is found to be very fast and accurate, down to extremely low temperatures. Moreover, it provides an intuitive and practical procedure for determining the capture distance (i.e., where the capture probability is evaluated), which would otherwise be arbitrary.

Accurate rovibrational energies of ozone isotopologues up to J = 10 utilizing artificial neural networks.

In recent years, ozone and its isotopologues have been a topic of interest in many fields of research, due to its importance in atmospheric chemistry and its anomalous isotopic enrichment-or the so-called "mass-independent fractionation." In the field of potential energy surface (PES) creation, debate over the existence of a potential barrier just under the dissociation threshold (referred to as a "potential reef") has plagued research for some years. Recently, Dawes and co-workers [Dawes, Lolur, Li, Jiang, and Guo (DLLJG) J. Chem. Phys. 139, 201103 (2013)] created a highly accurate global PES, for which the reef is found to be replaced with a (monotonic) "plateau." Subsequent dynamical calculations on this "DLLJG" PES have shown improved agreement with experiment, particularly the vibrational spectrum. However, it is well known that reaction dynamics is also highly influenced by the rovibrational states, especially in cases like ozone that assume a Lindemann-type mechanism. Accordingly, we present the first significant step toward a complete characterization of the rovibrational spectrum for various isotopologues of ozone, computed using the DLLJG PES together with the ScalIT suite of parallel codes. Additionally, artificial neural networks are used in an innovative fashion-not to construct the PES function per se but rather to greatly speed up its evaluation.

A quantum dynamical study of the rotation of the dihydrogen ligand in the Fe(H)2(H2)(PEtPh2)3 coordination complex.

Progress in the hydrogen fuel field requires a clear understanding and characterization of how materials of interest interact with hydrogen. Due to the inherently quantum mechanical nature of hydrogen nuclei, any theoretical studies of these systems must be treated quantum dynamically. One class of material that has been examined in this context are dihydrogen complexes. Since their discovery by Kubas in 1984, many such complexes have been studied both experimentally and theoretically. This particular study examines the rotational dynamics of the dihydrogen ligand in the Fe(H)2(H2)(PEtPh2)3 complex, allowing for full motion in both the rotational degrees of freedom and treating the quantum dynamics (QD) explicitly. A "gas-phase" global potential energy surface is first constructed using density functional theory with the Becke, 3-parameter, Lee-Yang-Parr functional; this is followed by an exact QD calculation of the corresponding rotation/libration states. The results provide insight into the dynamical correlation of the two rotation angles as well as a comprehensive analysis of both ground- and excited-state librational tunneling splittings. The latter was computed to be 6.914 cm-1-in excellent agreement with the experimental value of 6.4 cm-1. This work represents the first full-dimensional ab initio exact QD calculation ever performed for dihydrogen ligand rotation in a coordination complex.

Two-body Schrödinger wave functions in a plane-wave basis via separation of dimensions.

Using a combination of ideas, the ground and several excited electronic states of the helium atom and the hydrogen molecule are computed to chemical accuracy-i.e., to within 1-2 mhartree or better. The basic strategy is very different from the standard electronic structure approach in that the full two-electron six-dimensional (6D) problem is tackled directly, rather than starting from a single-electron Hartree-Fock approximation. Electron correlation is thus treated exactly, even though computational requirements remain modest. The method also allows for exact wave functions to be computed, as well as energy levels. From the full-dimensional 6D wave functions computed here, radial distribution functions and radial correlation functions are extracted-as well as a 2D probability density function exhibiting antisymmetry for a single Cartesian component. These calculations support a more recent interpretation of Hund's rule, which states that the lower energy of the higher spin-multiplicity states is actually due to reduced screening, rather than reduced electron-electron repulsion. Prospects for larger systems and/or electron dynamics applications appear promising.

Final State Resolved Quantum Predissociation Dynamics of SO2(C̃1B2) and Its Isotopomers via a Crossing with a Singlet Repulsive State.

The fragmentation dynamics of predissociative SO2(C̃1B2) is investigated on an accurate adiabatic potential energy surface (PES) determined from high level ab initio data. This singlet PES features non-C2v equilibrium geometries for SO2, which are separated from the SO(X̃3Σ-) + O(3P) dissociation limit by a barrier resulting from a conical intersection with a repulsive singlet state. The ro-vibrational state distribution of the SO fragment is determined quantum mechanically for many predissociative states of several sulfur isotopomers of SO2. Significant rotational and vibrational excitations are found in the SO fragment. It is shown that these fragment internal state distributions are strongly dependent on the predissociative vibronic states, and the excitation typically increases with the photon energy.