Full dimensional (15-dimensional) quantum-dynamical simulation of the protonated water dimer. II. Infrared spectrum and vibrational dynamics

Resource Optimization for Quantum Dynamics with Tensor Networks: Quantum and Classical Algorithms

The exponential scaling of the quantum degrees of freedom with the size of the system is one of the biggest challenges in computational chemistry and particularly in quantum dynamics. We present a tensor network approach for the time-evolution of the nuclear degrees of freedom of multiconfigurational chemical systems at a reduced storage and computational complexity. We also present quantum algorithms for the resultant dynamics. To preserve the compression advantage achieved via tensor network decompositions, we present an adaptive algorithm for the regularization of nonphysical bond dimensions, preventing the potentially exponential growth of these with time. While applicable to any quantum dynamical problem, our method is particularly valuable for dynamical simulations of nuclear chemical systems. Our algorithm is demonstrated using ab initio potentials obtained for a symmetric hydrogen-bonded system, namely, the protonated 2,2'-bipyridine, and compared to exact diagonalization numerical results.

A time averaged semiclassical approach to IR spectroscopy

We propose a new semiclassical approach to the calculation of molecular IR spectra. The method employs the time averaging technique of Kaledin and Miller upon symmetrization of the quantum dipole-dipole autocorrelation function. Spectra at high and low temperatures are investigated. In the first case, we are able to point out the possible presence of hot bands in the molecular absorption line shape. In the second case, we are able to reproduce accurate IR spectra as demonstrated by a calculation of the IR spectrum of the water molecule, which is within 4% of the exact intensity. Our time averaged IR spectra can be directly compared to time averaged semiclassical power spectra as shown in an application to the CO2 molecule, which points out the differences between IR and power spectra and demonstrates that our new approach can identify active IR transitions correctly. Overall, the method features excellent accuracy in calculating absorption intensities and provides estimates for the frequencies of vibrations in agreement with the corresponding power spectra. In perspective, this work opens up the possibility to interface the new method with the semiclassical techniques developed for power spectra, such as the divide-and-conquer one, to get accurate IR spectra of complex and high-dimensional molecular systems.

Near-infrared spectroscopy of H3O+⋯Xn (X = Ar, N2, and CO, n = 1-3)

Near-infrared (NIR) spectra of H3O+⋯Xn (X = Ar, N2, and CO, n = 1-3) in the first overtone region of OH-stretching vibrations (4800-7000 cm-1) were measured. Not only OH-stretching overtones but also several combination bands are major features in this region, and assignments of these observed bands are not obvious at a glance. High-precision anharmonic vibrational simulations based on the discrete variable representation approach were performed. The simulated spectra show good agreement with the observed ones and provide firm assignments of the observed bands, except in the case of X = CO, in which higher order vibrational mode couplings seem significant. This agreement demonstrates that the present system can be a benchmark for high precision anharmonic vibrational computations of NIR spectra. Band broadening in the observed spectra becomes remarkable with an increase of the interaction with the solvent molecule (X). The origin of the band broadening is explored by rare gas tagging experiments and anharmonic vibrational simulations of hot bands.

QuTree: A tree tensor network package

We present QuTree, a C++ library for tree tensor network approaches. QuTree provides class structures for tensors, tensor trees, and related linear algebra functions that facilitate the fast development of tree tensor network approaches such as the multilayer multiconfigurational time-dependent Hartree approach or the density matrix renormalization group approach and its various extensions. We investigate the efficiency of relevant tensor and tensor network operations and show that the overhead for managing the network structure is negligible, even in cases with a million leaves and small tensors. QuTree focuses on providing simple, high-level routines while retaining easy access to the backend to facilitate novel developments. We demonstrate the capabilities of the package by computing the eigenstates of coupled harmonic oscillator Hamiltonians and performing random circuit simulations on a virtual quantum computer.

Eigenstate calculation in the state-averaged (multi-layer) multi-configurational time-dependent Hartree approach

A new approach for the calculation of eigenstates with the state-averaged (multi-layer) multi-configurational time-dependent Hartree (MCTDH) approach is presented. The approach is inspired by the recent work of Larsson [J. Chem. Phys. 151, 204102 (2019)]. It employs local optimization of the basis sets at each node of the multi-layer MCTDH tree and successive downward and upward sweeps to obtain a globally converged result. At the top node, the Hamiltonian represented in the basis of the single-particle functions (SPFs) of the first layer is diagonalized. Here p wavefunctions corresponding to the p lowest eigenvalues are computed by a block Lanczos approach. At all other nodes, a non-linear operator consisting of the respective mean-field Hamiltonian matrix and a projector onto the space spanned by the respective SPFs is considered. Here, the eigenstate corresponding to the lowest eigenvalue is computed using a short iterative Lanczos scheme. Two different examples are studied to illustrate the new approach: the calculation of the vibrational states of methyl and acetonitrile. The calculations for methyl employ the single-layer MCTDH approach, a general potential energy surface, and the correlation discrete variable representation. A five-layer MCTDH representation and a sum of product-type Hamiltonian are used in the acetonitrile calculations. Very fast convergence and order of magnitude reductions in the numerical effort compared to the previously used block relaxation scheme are found. Furthermore, a detailed comparison with the results of Avila and Carrington [J. Chem. Phys. 134, 054126 (2011)] for acetonitrile highlights the potential problems of convergence tests for high-dimensional systems.

Vibrational Spectra of Highly Anharmonic Water Clusters: Molecular Dynamics and Harmonic Analysis Revisited with Constrained Nuclear-Electronic Orbital Methods

Vibrational spectroscopy is widely used to gain insights into structural and dynamic properties of chemical, biological, and materials systems. Thus, an efficient and accurate method to simulate vibrational spectra is desired. In this paper, we justify and employ a microcanonical molecular simulation scheme to calculate the vibrational spectra of three challenging water clusters: the neutral water dimer (H4O2), the protonated water trimer (H7O3+), and the protonated water tetramer (H9O4+). We find that with the accurate description of quantum nuclear delocalization effects through the constrained nuclear-electronic orbital framework, including vibrational mode coupling effects through molecular dynamics simulations can additionally improve the vibrational spectrum calculations. In contrast, without the quantum nuclear delocalization picture, conventional ab initio molecular dynamics may even lead to less accurate results than harmonic analysis.

Approaching the "Zundel" Limit: Tuning the Vibrational Coupling in N2H+Ng, Ng = {He, Ne, Ar, Kr, Xe, and Rn}

The diazenylium ion (N₂H⁺) is a ubiquitous ion in dense molecular clouds. This ion is often used as a dense gas tracer in outer space. Most of the previous works on diazenylium ion have focused on the shared-proton stretch band, ν_H⁺. In this work, we have performed reduced-dimensional calculations to investigate the vibrational structure of N₂H⁺Ng, Ng = {He, Ne, Ar, Kr, Xe, and Rn}. We demonstrate a few interesting things about this system. First, the vibrational coupling in N₂H⁺ can be tuned to switch on interesting anharmonic effects such as Fermi resonance or combination bands by tagging it with different noble gases. Second, a comparison of the vibrational spectrum from N₂H⁺He to N₂H⁺Rn shows that the ν_H⁺ can be swept from an "Eigen-like" to a "Zundel-like" limiting case. Anharmonic calculations were performed using a multilevel approach, which utilized the MP2 and CCSD(T) levels of theories. Binding energies for the elimination of Ng in N₂H⁺Ng are also reported.

Calculating Vibrational Excited State Absorptions with Excited State Constrained Minimized Energy Surfaces

The modeling and interpretation of vibrational spectra are crucial for studying reaction dynamics using vibrational spectroscopy. Most prior theoretical developments focused on describing fundamental vibrational transitions while fewer developments focused on vibrational excited state absorptions. In this study, we present a new method that uses excited state constrained minimized energy surfaces (CMESs) to describe vibrational excited state absorptions. The excited state CMESs are obtained similarly to the previous ground state CMES development in our group but with additional wave function orthogonality constraints. Using a series of model systems, including the harmonic oscillator, Morse potential, double-well potential, quartic potential, and two-dimensional anharmonic potential, we demonstrate that this new procedure provides good estimations of the transition frequencies for vibrational excited state absorptions. These results are significantly better than those obtained from harmonic approximations using conventional potential energy surfaces, demonstrating the promise of excited state CMES-based methods for calculating vibrational excited state absorptions in real systems.

Optimal Mode Combination in the Multiconfiguration Time-Dependent Hartree Method through Multivariate Statistics: Factor Analysis and Hierarchical Clustering

The multiconfiguration time-dependent Hartree (MCTDH) method and its multilayer extension (ML-MCTDH) are powerful algorithms for the efficient computation of nuclear quantum dynamics in high-dimensional systems. By providing time-dependent variational orbitals and an optimal choice of layered effective degrees of freedom, one is able to reduce the computational cost to an amenable number of configurations. However, choices related to selecting properly the mode grouping and tensor tree are strongly system dependent and, thus far, subjectively based on intuition and/or experience. Therefore, herein we detail a new protocol based on multivariate statistics─more specifically, factor analysis and hierarchical clustering─for a reliable and convenient guiding in the optimal design of such complex "system-of-systems" tensor-network decompositions. The advantages of employing the new algorithm and its applicability are tested on water and two floppy protonated water clusters with large amplitude motions.

Proton-transfer spectroscopy beyond the normal-mode scenario

A stochastic theory is developed to predict the spectral signature of proton-transfer processes and is applied to infrared spectra computed from ab initio molecular-dynamics simulations of a single [Formula: see text] cation. By constraining the oxygen atoms to a fixed distance, this system serves as a tunable model for general proton-transfer processes with variable barrier height. Three spectral contributions at distinct frequencies are identified and analytically predicted: the quasi-harmonic motion around the most probable configuration, amenable to normal-mode analysis, the contribution due to transfer paths when the proton moves over the barrier, and a shoulder for low frequencies stemming from the stochastic transfer-waiting-time distribution; the latter two contributions are not captured by normal-mode analysis but exclusively reported on the proton-transfer kinetics. In accordance with reaction rate theory, the transfer-waiting-contribution frequency depends inversely exponentially on the barrier height, whereas the transfer-path-contribution frequency is rather insensitive to the barrier height.

The coupling of the hydrated proton to its first solvation shell

The Zundel (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{5}{O}_{2}^{+}$$\end{document}H5O2+) and Eigen (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{9}{O}_{4}^{+}$$\end{document}H9O4+) cations play an important role as intermediate structures for proton transfer processes in liquid water. In the gas phase they exhibit radically different infrared (IR) spectra. The question arises: is there a least common denominator structure that explains the IR spectra of both, the Zundel and Eigen cations, and hence of the solvated proton? Full dimensional quantum simulations of these protonated cations demonstrate that two dynamical water molecules and an excess proton constitute this fundamental subunit. Embedded in the static environment of the parent Eigen cation, this subunit reproduces the positions and broadenings of its main excess-proton bands. In isolation, its spectrum reverts to the well-known Zundel ion. Hence, the dynamics of this subunit polarized by an environment suffice to explain the spectral signatures and anharmonic couplings of the solvated proton in its first solvation shell.

The Zundel [H(H₂O)₂]⁺ and Eigen [H(H₂O)₄]⁺ cations exhibit radicallly different infrared spectra and are the limiting dynamical structures involved in proton mobility in liquid water. Here, the authors find through quantum dynamics simulations that two polarized water molecules and a proton suffice to explain the key spectroscopic features connected to proton mobility for both species.

State-resolved infrared spectrum of the protonated water dimer: revisiting the characteristic proton transfer doublet peak

The infrared (IR) spectra of protonated water clusters encode precise information on the dynamics and structure of the hydrated proton. However, the strong anharmonic coupling and quantum effects of these elusive species remain puzzling up to the present day. Here, we report unequivocal evidence that the interplay between the proton transfer and the water wagging motions in the protonated water dimer (Zundel ion) giving rise to the characteristic doublet peak is both more complex and more sensitive to subtle energetic changes than previously thought. In particular, hitherto overlooked low-intensity satellite peaks in the experimental spectrum are now unveiled and mechanistically assigned. Our findings rely on the comparison of IR spectra obtained using two highly accurate potential energy surfaces in conjunction with highly accurate state-resolved quantum simulations. We demonstrate that these high-accuracy simulations are important for providing definite assignments of the complex IR signals of fluxional molecules.

Adaptive Fitting of Potential Energy Surfaces of Small to Medium-Sized Molecules in Sum-of-Product Form: Application to Vibrational Spectroscopy

A semi-automatic sampling and fitting procedure for generating sum-of-product (Born-Oppenheimer) potential energy surfaces based on a high-dimensional model representation is presented. The adaptive sampling procedure and subsequent fitting relies on energies only and can be used for re-fitting existing analytic potential energy surfaces in sum-of-product form or for direct fits from ab initio computa- tions. The method is tested by fitting ground electronic state potential energy surfaces for small to medium sized semi-rigid molecules, i.e., HFCO, HONO, and HCOOH, based upon ab initio computations at the CCSD(T)-F12/cc-pVTZ-F12 or MP2/aug-cc-pVTZ levels of theory. Vibrational eigenstates are computed using block improved relaxation in the Heidelberg MCTDH package and compared to available experimental and theoretical data. The new potential energy surfaces are compared to the best ones currently available for these molecules, in terms of accuracy, including of resulting vibrational states, required numbers of sampling points, and number of fitting parameters. The present procedure leads to compact expansions and scales well with the number of dimensions for simple potentials such as single or double wells.

A non-hierarchical correlation discrete variable representation

The correlation discrete variable representation (CDVR) facilitates (multi-layer) multi-configurational time-dependent Hartree (MCTDH) calculations with general potentials. It employs a layered grid representation to efficiently evaluate all potential matrix elements appearing in the MCTDH equations of motion. The original CDVR approach and its multi-layer extension show a hierarchical structure: the size of the grids employed at the different layers increases when moving from an upper layer to a lower one. In this work, a non-hierarchical CDVR approach, which uses identically structured quadratures at all layers of the MCTDH wavefunction representation, is introduced. The non-hierarchical CDVR approach crucially reduces the number of grid points required, compared to the hierarchical CDVR, shows superior scaling properties, and yields identical results for all three representations showing the same topology. Numerical tests studying the photodissociation of NOCl and the vibrational states of CH₃ demonstrate the accuracy of the non-hierarchical CDVR approach.

Exploring the Origins of Spectral Signatures of Strong Hydrogen Bonding in Protonated Water Clusters

The effects of anharmonicity on the spectral features of strong ionic hydrogen bonds are explored through reduced dimensional studies of the couplings between the hydrogen bonding OH and the donor-acceptor OO stretching vibrations in protonated water clusters with 2-4 water molecules. Specifically, this study focuses on how the anharmonicities and couplings in these ions are reflected in the vibrational spectra by exploring the intensities of the transitions to states with excitation in both the OH and the OO stretching vibrations and changes in the frequency of the OO stretching vibration when the OH stretching vibration is excited. These questions are addressed through the application of several approximate treatments that are based on an adiabatic separation of the high-frequency OH and low-frequency OO stretching vibrations as well as low-order expansions of the potential and dipole surfaces. While an adiabatic approximation captures most of the trends found in the spectra and from an analysis of the two-dimensional model, a vibrational Franck-Condon approach fails to capture the intensities of these transitions. Of the terms in the expansion of the dipole moment function, those that are proportional to Δr_OH and Δr_OH² are found to provide the largest contributions to the calculated intensities of the transitions involving excitation of both the OH and the OO stretches. This leads to the conclusion that the intensities of these transitions encode information about the frequency and anharmonicity of the OH stretching vibration and how they are affected by changes in the OO distance. The anharmonicity of the potential also leads to changes in the OO stretching frequency with excitation of the OH stretching vibration. The direction of this change in frequency encodes additional information about the strength of the ionic hydrogen bond.

Time evolution of ML-MCTDH wavefunctions. II. Application of the projector splitting integrator

The multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) approach can suffer from numerical instabilities whenever the wavefunction is weakly entangled. These instabilities arise from singularities in the equations of motion (EOMs) and necessitate the use of regularization of the EOMs. The Projector Splitting Integrator (PSI) has previously been presented as an approach for evolving ML-MCTDH wavefunctions that is free of singularities. Here, we will discuss the implementation of the multi-layer PSI with a particular focus on how the steps required relate to those required to implement standard ML-MCTDH. We demonstrate the efficiency and stability of the PSI for large ML-MCTDH wavefunctions containing up to hundreds of thousands of nodes by considering a series of spin-boson models with up to 10⁶ bath modes and find that for these problems, the PSI requires roughly 3-4 orders of magnitude fewer Hamiltonian evaluations and 2-3 orders of magnitude fewer Hamiltonian applications than standard ML-MCTDH and 2-3/1-2 orders of magnitude fewer evaluations/applications than approaches that use improved regularization schemes. Finally, we consider a series of significantly more challenging multi-spin-boson models that require much larger numbers of single-particle functions with wavefunctions containing up to ∼1.3×10⁹ parameters to obtain accurate dynamics.

Time evolution of ML-MCTDH wavefunctions. I. Gauge conditions, basis functions, and singularities

We derive a family of equations-of-motion (EOMs) for evolving multi-layer multiconfiguration time-dependent Hartree (ML-MCTDH) wavefunctions that, unlike the standard ML-MCTDH EOMs, never require the evaluation of the inverse of singular matrices. All members of this family of EOMs make use of alternative static gauge conditions than those used for standard ML-MCTDH. These alternative conditions result in an expansion of the wavefunction in terms of a set of potentially arbitrary orthonormal functions, rather than in terms of a set of non-orthonormal and potentially linearly dependent functions, as is the case for standard ML-MCTDH. We show that the EOMs used in the projector splitting integrator (PSI) and the invariant EOM approaches are two special cases of this family obtained from different choices for the dynamic gauge condition, with the invariant EOMs making use of a choice that introduces potentially unbounded operators into the EOMs. As a consequence, all arguments for the existence of parallelizable integration schemes for the invariant EOMs can also be applied to the PSI EOMs.

Coupled Local-Mode Approach for the Calculation of Vibrational Spectra: Application to Protonated Water Clusters

Spectroscopic studies of protonated water clusters (PWCs) have yielded enormous insights into the fundamental nature of the hydrated proton. Here, we introduce a new coupled local-mode (CLM) approach to calculate PWC OH stretch vibrational spectra. The CLM method combines a sampling of representative configurations from density functional theory (DFT)-based ab initio molecular dynamics (AIMD) simulations with DFT calculations of local-mode vibrational frequencies and couplings. Calculations of inhomogeneous OH stretch vibrational spectra for H⁺(H₂O)₄ and H⁺(H₂O)₂₁ agree well with experiment and higher-level calculations, and decompositions of the calculated spectra in terms of the coupled modes aids in the interpretation of the spectra. This observation is consistent with the idea that capturing anharmonicity and coupling is as important to accuracy as the underlying level of electronic structure theory. The CLM calculations can easily discern the configuration that dominates the experimental measurement for H⁺(H₂O)₅, which can adopt several low-energy conformations.

Using Diffusion Monte Carlo Wave Functions to Analyze the Vibrational Spectra of H7O3+ and H9O4. J Phys Chem A 2021;125:7185-7197. [PMID: 34433268 DOI: 10.1021/acs.jpca.1c05025]

An approach for evaluating spectra from ground state probability amplitudes (GSPA) obtained from diffusion Monte Carlo (DMC) simulations is extended to improve the description of excited state energies and allow for coupling among vibrational excited states. This approach is applied to studies of the protonated water trimer and tetramer, and their deuterated analogs. These ions provide models for solvated hydronium, and analysis of these spectra provides insights into spectral signatures of proton transfer in aqueous environments. In this approach, we obtain a separable set of internal coordinates from the DMC ground state probability amplitude. A basis is then developed from products of the DMC ground state wave function and low-order polynomials in these internal coordinates. This approach provides a compact basis in which the Hamiltonian and dipole moment matrix are evaluated and used to obtain the spectrum. The resulting spectra are in good agreement with experiment and in many cases provide comparable agreement to the results obtained using much larger basis sets. In addition, the compact basis allows for interpretation of the spectral features and how they evolve with cluster size and deuteration.

Electron binding energies and Dyson orbitals of OnH2n+1 +,0,- clusters: Double Rydberg anions, Rydberg radicals, and micro-solvated hydronium cations

Ab initio electron propagator methods are employed to predict the vertical electron attachment energies (VEAEs) of OH₃ ⁺(H₂O)_n clusters. The VEAEs decrease with increasing n, and the corresponding Dyson orbitals are diffused over exterior, non-hydrogen bonded protons. Clusters formed from OH₃ ^- double Rydberg anions (DRAs) and stabilized by hydrogen bonding or electrostatic interactions between ions and polar molecules are studied through calculations on OH₃ ^-(H₂O)_n complexes and are compared with more stable H^-(H₂O)_n+1 isomers. Remarkable changes in the geometry of the anionic hydronium-water clusters with respect to their cationic counterparts occur. Rydberg electrons in the uncharged and anionic clusters are held near the exterior protons of the water network. For all values of n, the anion-water complex H^-(H₂O)_n+1 is always the most stable, with large vertical electron detachment energies (VEDEs). OH₃ ^-(H₂O)_n DRA isomers have well separated VEDEs and may be visible in anion photoelectron spectra. Corresponding Dyson orbitals occupy regions beyond the peripheral O-H bonds and differ significantly from those obtained for the VEAEs of the cations.

Symmetries in the multi-configurational time-dependent Hartree wavefunction representation and propagation

In multi-configurational time-dependent Hartree (MCTDH) approaches, different multi-layered wavefunction representations can be used to represent the same physical wavefunction. Transformations between different equivalent representations of a physical wavefunction that alter the tree structure used in the multi-layer MCTDH wavefunction representation interchange the role of single-particle functions (SPFs) and single-hole functions (SHFs) in the MCTDH formalism. While the physical wavefunction is invariant under these transformations, this invariance does not hold for the standard multi-layer MCTDH equations of motion. Introducing transformed SPFs, which obey normalization conditions typically associated with SHFs, revised equations of motion are derived. These equations do not show the singularities resulting from the inverse single-particle density matrix and are invariant under tree transformations. Based on the revised equations of motion, a new integration scheme is introduced. The scheme combines the advantages of the constant mean-field approach of Beck and Meyer [Z. Phys. D 42, 113 (1997)] and the singularity-free integrator suggested by Lubich [Appl. Math. Res. Express 2015, 311]. Numerical calculations studying the spin boson model in high dimensionality confirm the favorable properties of the new integration scheme.

Evolution of the linker in microhydrated hydrogen dinitrate anions: From H+ to H5 O2. J Comput Chem 2021;42:1514-1525. [PMID: 33990989 DOI: 10.1002/jcc.26560] [Citation(s) in RCA: 0] [Impact Index Per Article: 0]

Hydrogen dinitrate anion, HNO₃ (NO₃ ^- ), is a proton-bound dimer with a very strong hydrogen bond. By employing ab initio molecular dynamics (AIMD) method, we studied the effects of the proton transfer and the rotation of the nitrates on the vibrational profiles of HNO₃ (NO₃ ^- )(H₂ O)_n (n = 0-2). The AIMD results indicate that the structure of the n = 0 cluster is very flexible, even though its hydrogen bond is quite strong. Significant rotations around the hydrogen bond and frequent transfers of proton from HNO₃ to NO₃ ^- are observed in AIMD simulations. Dynamic changes are therefore an important factor in understanding the broadening of vibrational features. For n = 1, the extent of structural fluctuation increases further, as H₂ O could move around the anion while the HNO₃ (NO₃ ^- ) core also goes through structural changes. Its vibrational spectrum can be understood as a mixture of many isomers visited during AIMD simulations. By n = 2, the structure is stabilized around one isomer, with the linker between the two nitrates being H₅ O₂ ⁺ , rather than H⁺ . Due to strong hydrogen bonds between nitrates and water molecules, this H₅ O₂ ⁺ takes the extraordinary structure with the H⁺ localized on one H₂ O, rather than being shared. While this novel structure is stable during AIMD simulations, the dynamic fluctuations in hydrogen bond distances still produce significant broadening in its vibrational profile.

Comparison of the multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) method and the density matrix renormalization group (DMRG) for ground state properties of linear rotor chains

We demonstrate the applicability of the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) method to the problem of computing ground states of one-dimensional chains of linear rotors with dipolar interactions. Specifically, we successfully obtain energies, entanglement entropies, and orientational correlations that are in agreement with the Density Matrix Renormalization Group (DMRG), which has been previously used for this system. We find that the entropies calculated by ML-MCTDH for larger system sizes contain nonmonotonicity, as expected in the vicinity of a second-order quantum phase transition between ordered and disordered rotor states. We observe that this effect remains when all couplings besides nearest-neighbor are omitted from the Hamiltonian, which suggests that it is not sensitive to the rate of decay of the interactions. In contrast to DMRG, which is tailored to the one-dimensional case, ML-MCTDH (as implemented in the Heidelberg MCTDH package) requires more computational time and memory, although the requirements are still within reach of commodity hardware. The numerical convergence and computational demand of two practical implementations of ML-MCTDH and DMRG are presented in detail for various combinations of system parameters.

Weighted-Graph-Theoretic Methods for Many-Body Corrections within ONIOM: Smooth AIMD and the Role of High-Order Many-Body Terms

We present a weighted-graph-theoretic approach to adaptively compute contributions from many-body approximations for smooth and accurate post-Hartree-Fock (pHF) ab initio molecular dynamics (AIMD) of highly fluxional chemical systems. This approach is ONIOM-like, where the full system is treated at a computationally feasible quality of treatment (density functional theory (DFT) for the size of systems considered in this publication), which is then improved through a perturbative correction that captures local many-body interactions up to a certain order within a higher level of theory (post-Hartree-Fock in this publication) described through graph-theoretic techniques. Due to the fluxional and dynamical nature of the systems studied here, these graphical representations evolve during dynamics. As a result, energetic "hops" appear as the graphical representation deforms with the evolution of the chemical and physical properties of the system. In this paper, we introduce dynamically weighted, linear combinations of graphs, where the transition between graphical representations is smoothly achieved by considering a range of neighboring graphical representations at a given instant during dynamics. We compare these trajectories with those obtained from a set of trajectories where the range of local many-body interactions considered is increased, sometimes to the maximum available limit, which yields conservative trajectories as the order of interactions is increased. The weighted-graph approach presents improved dynamics trajectories while only using lower-order many-body interaction terms. The methods are compared by computing dynamical properties through time-correlation functions and structural distribution functions. In all cases, the weighted-graph approach provides accurate results at a lower cost.

Demystifying the Diffuse Vibrational Spectrum of Aqueous Protons Through Cold Cluster Spectroscopy

The ease with which the pH is routinely determined for aqueous solutions masks the fact that the cationic product of Arrhenius acid dissolution, the hydrated proton, or H⁺(aq), is a remarkably complex species. Here, we review how results obtained over the past 30 years in the study of H⁺⋅(H₂O)_n cluster ions isolated in the gas phase shed light on the chemical nature of H⁺(aq). This effort has also revealed molecular-level aspects of the Grotthuss relay mechanism for positive-charge translocation in water. Recently developed methods involving cryogenic cooling in radiofrequency ion traps and the application of two-color, infrared-infrared (IR-IR) double-resonance spectroscopy have established a clear picture of how local hydrogen-bond topology drives the diverse spectral signatures of the excess proton. This information now enables a new generation of cluster studies designed to unravel the microscopic mechanics underlying the ultrafast relaxation dynamics displayed by H⁺(aq).

Fermi resonance switching in KrH+Rg and XeH+Rg (Rg = Ne, Ar, Kr, and Xe)

Matrix isolation experiments have been successfully employed to extensively study the infrared spectrum of several proton-bound rare gas complexes. Most of these studies have focused on the spectral signature for the H⁺ stretch (ν₃) and its combination bands with the intermolecular stretch coordinate (ν₁). However, little attention has been paid to the Fermi resonance interaction between the H⁺ stretch (ν₃) and H⁺ bend overtone (2ν₂) in the asymmetric proton-bound rare gas dimers, RgH⁺Rg'. In this work, we have investigated this interaction on KrH⁺Rg and XeH⁺Rg with Rg = (Ne, Ar, Kr, and Xe). A multilevel potential energy surface (PES) was used to simulate the vibrational structure of these complexes. This PES is a dual-level comprising of second-order Møller-Plesset perturbation theory and coupled-cluster singles doubles with perturbative triples [CCSD(T)] levels of ab initio theories. We found that when both the combination bands (nν₁ + ν₃) and bend overtone 2ν₂ compete to borrow intensity from the ν₃ band, the latter wins over the former, which then results in the suppression of the nν₁ + ν₃ bands. The current simulations offer new assignments for the ArH⁺Xe and KrH⁺Xe spectra. Complete basis set (CBS) binding energies for these complexes were also calculated at the CCSD(T)/CBS level.

Direct product-type grid representations for angular coordinates in extended space and their application in the MCTDH approach

Multi-configurational time-dependent Hartree (MCTDH) calculations using time-dependent grid representations can be used to accurately simulate high-dimensional quantum dynamics on general ab initio potential energy surfaces. Employing the correlation discrete variable representation, sets of direct product type grids are employed in the calculation of the required potential energy matrix elements. This direct product structure can be a problem if the coordinate system includes polar and azimuthal angles that result in singularities in the kinetic energy operator. In the present work, a new direct product-type discrete variable representation (DVR) for arbitrary sets of polar and azimuthal angles is introduced. It employs an extended coordinate space where the range of the polar angles is taken to be [-π, π]. The resulting extended space DVR resolves problems caused by the singularities in the kinetic energy operator without generating a very large spectral width. MCTDH calculations studying the F·CH₄ complex are used to investigate important properties of the new scheme. The scheme is found to allow for more efficient integration of the equations of motion compared to the previously employed cot-DVR approach [G. Schiffel and U. Manthe, Chem. Phys. 374, 118 (2010)] and decreases the required central processing unit times by about an order of magnitude.

Molecular Characterization of the Surface Excess Charge Layer in Droplets

The surface excess charge layer (SECL) in droplets has often been associated with distinct chemistry. We examine the effect of the nature of ions in the composition and structure of SECL by using molecular dynamics. We find that in the presence of simple ions the thickness of SECL is invariant not only with respect to droplet size but also with respect to the nature of the ions. In the presence of simple ions, this layer has a thickness of ∼1.5-1.7 nm but in the presence of macroions it may extend to ∼2.0 nm. The proportion of ions contained in SECL depends on the nature of the ions and the droplet size. For the same droplet size, I^- and model H₃O⁺ ions show considerably higher concentration than Na⁺ and Cl^- ions. We identify the maximum ion concentration region, which, in nanodrops, may partially overlap with SECL. As the relative shape fluctuations decrease when microdrop size is approached, the overlap between SECL and maximum ion concentration region increases. We suggest the extension of the bilayer droplet structure assumed in the equilibrium partitioning model of Enke to include the maximum ion concentration region that may not coincide with SECL in nanodrops. We compute the ion concentrations in SECL, which are those that should enter the kinetic equation in the ion-evaporation mechanism, instead of the overall drop ion concentration that has been used.

Representing molecular ground and excited vibrational eigenstates with nuclear densities obtained from semiclassical initial value representation molecular dynamics

We present in detail and validate an effective Monte Carlo approach for the calculation of the nuclear vibrational densities via integration of molecular eigenfunctions that we have preliminary employed to calculate the densities of the ground and the excited OH stretch vibrational states in the protonated glycine molecule [Aieta et al., Nat Commun 11, 4348 (2020)]. Here, we first validate and discuss in detail the features of the method on a benchmark water molecule. Then, we apply it to calculate on-the-fly the ab initio anharmonic nuclear densities in the correspondence of the fundamental transitions of NH and CH stretches in protonated glycine. We show how we can gain both qualitative and quantitative physical insight by inspection of different one-nucleus densities and assign a character to spectroscopic absorption peaks using the expansion of vibrational states in terms of harmonic basis functions. The visualization of the nuclear vibrations in a purely quantum picture allows us to observe and quantify the effects of anharmonicity on the molecular structure, also to exploit the effect of IR excitations on specific bonds or functional groups, beyond the harmonic approximation. We also calculate the quantum probability distribution of bond lengths, angles, and dihedrals of the molecule. Notably, we observe how in the case of one type of fundamental NH stretching, the typical harmonic nodal pattern is absent in the anharmonic distribution.

The "Hole" Story in Ionized Water from the Perspective of Ehrenfest Dynamics

The radiolysis of liquid water and the radiation-matter interactions that happen in aqueous environments are important to the fields of chemistry, materials, and environmental sciences, as well as the biological and physiological response to extreme conditions and medical treatments. The initial stage of radiolysis is the ultrafast response, or hole dynamics, that triggers chemical processes within complex energetic landscapes that may include reactivity. A fundamental understanding necessitates the use of theoretical methods that are capable of simulating both ultrafast coherence and non-adiabatic energy transfer pathways. In this work, we carry out an ab initio Ehrenfest dynamics study to provide a more complete description of the ultrafast dynamics and reactive events initiated by photoionization of water. After sudden ionization, a range of processes, including hole trapping and transfer, large OH oscillations, proton transfer and subsequent relay, formation of the metastable Zundel complex, and long-lived coherence, are identified and new insight into their driving forces is elucidated.

Decoding the 2D IR spectrum of the aqueous proton with high-level VSCF/VCI calculations

The aqueous proton is a common and long-studied species in chemistry, yet there is currently intense interest devoted to understanding its hydration structure and transport dynamics. Typically described in terms of two limiting structures observed in gas-phase clusters, the Zundel H₅O₂ ⁺ and Eigen H₉O₄ ⁺ ions, the aqueous structure is less clear due to the heterogeneity of hydrogen bonding environments and room-temperature structural fluctuations in water. The linear infrared (IR) spectrum, which reports on structural configurations, is challenging to interpret because it appears as a continuum of absorption, and the underlying vibrational modes are strongly anharmonically coupled to each other. Recent two-dimensional IR (2D IR) experiments presented strong evidence for asymmetric Zundel-like motifs in solution, but true structure-spectrum correlations are missing and complicated by the anharmonicity of the system. In this study, we employ high-level vibrational self-consistent field/virtual state configuration interaction calculations to demonstrate that the 2D IR spectrum reports on a broad distribution of geometric configurations of the aqueous proton. We find that the diagonal 2D IR spectrum around 1200 cm^-1 is dominated by the proton stretch vibrations of Zundel-like and intermediate geometries, broadened by the heterogeneity of aqueous configurations. There is a wide distribution of multidimensional potential shapes for the proton stretching vibration with varying degrees of potential asymmetry and confinement. Finally, we find specific cross peak patterns due to aqueous Zundel-like species. These studies provide clarity on highly debated spectral assignments and stringent spectroscopic benchmarks for future simulations.

Anharmonic quantum nuclear densities from full dimensional vibrational eigenfunctions with application to protonated glycine

The interpretation of molecular vibrational spectroscopic signals in terms of atomic motion is essential to understand molecular mechanisms and for chemical characterization. The signals are usually assigned after harmonic normal mode analysis, even if molecular vibrations are known to be anharmonic. Here we obtain the quantum anharmonic vibrational eigenfunctions of the 11-atom protonated glycine molecule and we calculate the density distribution of its nuclei and its geometry parameters, for both the ground and the O-H stretch excited states, using our semiclassical method based on ab initio molecular dynamics trajectories. Our quantum mechanical results describe a molecule elongated and more flexible with respect to what previously thought. More importantly, our method is able to assign each spectral peak in vibrational spectroscopy by showing quantitatively how normal modes involving different functional groups cooperate to originate that spectroscopic signal. The method will possibly allow for a better rationalization of experimental spectroscopy.

Inexpensive modeling of quantum dynamics using path integral generalized Langevin equation thermostats

The properties of molecules and materials containing light nuclei are affected by their quantum mechanical nature. Accurate modeling of these quantum nuclear effects requires computationally demanding path integral techniques. Considerable success has been achieved in reducing the cost of such simulations by using generalized Langevin dynamics to induce frequency-dependent fluctuations. Path integral generalized Langevin equation methods, however, have this far been limited to the study of static, thermodynamic properties due to the large perturbation to the system's dynamics induced by the aggressive thermostatting. Here, we introduce a post-processing scheme, based on analytical estimates of the dynamical perturbation induced by the generalized Langevin dynamics, which makes it possible to recover meaningful time correlation properties from a thermostatted trajectory. We show that this approach yields spectroscopic observables for model and realistic systems that have an accuracy comparable to much more demanding approximate quantum dynamics techniques based on full path integral simulations.

Reduced rovibrational coupling Cartesian dynamics for semiclassical calculations: Application to the spectrum of the Zundel cation

We study the vibrational spectrum of the protonated water dimer, by means of a divide-and-conquer semiclassical initial value representation of the quantum propagator, as a first step in the study of larger protonated water clusters. We use the potential energy surface from the work of Huang et al. [J. Chem. Phys. 122, 044308 (2005)]. To tackle such an anharmonic and floppy molecule, we employ fully Cartesian dynamics and carefully reduce the coupling to global rotations in the definition of normal modes. We apply the time-averaging filter and obtain clean power spectra relative to suitable reference states that highlight the spectral peaks corresponding to the fundamental excitations of the system. Our trajectory-based approach allows for the physical interpretation of the very challenging proton transfer modes. We find that it is important, for such a floppy molecule, to selectively avoid initially exciting lower energy modes, in order to obtain cleaner spectra. The estimated vibrational energies display a mean absolute error (MAE) of ∼29 cm^-1 with respect to available multiconfiguration time-dependent Hartree calculations and MAE ∼ 14 cm^-1 when compared to the optically active experimental excitations of the Ne-tagged Zundel cation. The reasonable scaling in the number of trajectories for Monte Carlo convergence is promising for applications to higher dimensional protonated cluster systems.

Disentangling the Complex Vibrational Mechanics of the Protonated Water Trimer by Rational Control of Its Hydrogen Bonds

The vibrational spectrum of the protonated water trimer, H⁺(H₂O)₃, is surprisingly complex, with many strong features in the expected region of the fundamentals associated with two H-bonded OH groups on the H₃O⁺ core ion. Here we follow how the bands in this region of the spectrum evolve when the energies of the fundamentals in the H-bonded OH stretches are systematically increased by the attachment of increasingly strongly bound "tag" molecules (He, Ar, D₂, N₂, CO, and H₂O) to the free OH position on the hydronium core ion of H⁺(H₂O)₃, as well as by replacement of the hydrogen atom in the nonbonded OH group on hydronium with methyl and ethyl groups. This allows for the incremental transformation of the complex band pattern observed in H⁺(H₂O)₃ into that of the "Eigen" structure of the protonated water tetramer. Differences among the trajectories of the various bands provide an empirical way to disentangle features primarily due to the displacements of the OH stretches bound to the hydronium core from those arising from anharmonic coupling to states involving one or more quanta in lower frequency modes. The latter are found to be dramatically enhanced when the nominal frequencies of the intermolecular OH stretching modes approach those of the intramolecular bends of the H₃O⁺ and H₂O constituents in both H and D isotopologues.

Semiclassical vibrational spectroscopy with Hessian databases

We report on a new approach to ease the computational overhead of ab initio "on-the-fly" semiclassical dynamics simulations for vibrational spectroscopy. The well known bottleneck of such computations lies in the necessity to estimate the Hessian matrix for propagating the semiclassical pre-exponential factor at each step along the dynamics. The procedure proposed here is based on the creation of a dynamical database of Hessians and associated molecular geometries able to speed up calculations while preserving the accuracy of results at a satisfactory level. This new approach can be interfaced to both analytical potential energy surfaces and on-the-fly dynamics, allowing one to study even large systems previously not achievable. We present results obtained for semiclassical vibrational power spectra of methane, glycine, and N-acetyl-L-phenylalaninyl-L-methionine-amide, a molecule of biological interest made of 46 atoms.

Hydrated excess protons and their local hydrogen bond transport network as measured by translational, librational, and vibrational frequencies

A clear molecular description of excess hydrated protons and their local hydrogen bond transport network remains elusive. Here, the hydrogen bond network of excess hydrated protons in water bridges was probed by measuring their Raman spectra and comparing them to the spectra of protons in ice and water. The proton vibrational spectrum and the hydrogen bond network translational and librational spectra were recorded. The spectra of the water bridge and water exhibit clear differences, indicating the presence of a structure in water bridges when subjected to an electric field of ∼10⁶ V/m that has not been previously reported. The intermolecular Raman spectrum of the floating water bridge exhibits a hydrogen bond stretching band at 150-250 cm^-1, librational bands within the 300-1000 cm^-1 spectral range, and a large band at 1500-3000 cm^-1, which corresponds to the vibrational signature of excess hydrated protons in the water bridge structure. The excess protons are shown to move predominantly at the air/water interface, and the effect of this distribution is a measurable change in the air/water interfacial tension from ∼80 to ∼32 N/m. Therefore, hydrated protons must have a unique water arrangement that enables them to propagate without sinking into bulk water. This local polarized hydrogen bond network in the interfacial water region is characterized by a translational spectrum similar to that of ice V.