Converged quantum dynamics with modified Shepard interpolation and Gaussian wave packets

Quasi-direct Quantum Molecular Dynamics: The Time-Dependent Adaptive Density-Guided Approach for Potential Energy Surface Construction

We present a new quasi-direct quantum molecular dynamics computational method which offers a compromise between quantum dynamics using a precomputed potential energy surface (PES) and fully direct quantum dynamics. This method is termed the time-dependent adaptive density-guided approach (TD-ADGA) and is a method for constructing a PES on the fly during a dynamics simulation. This is achieved by acquisition of new single-point (SP) calculations and refitting of the PES, depending on the need of the dynamics. The TD-ADGA is a further development of the adaptive density-guided approach (ADGA) for PES construction where the placement of SPs is guided by the density of the nuclear wave function. In TD-ADGA, the ADGA framework has been integrated into the time propagation of the time-dependent nuclear wave function and we use the reduced one-mode density of this wave function to guide when and where new SPs are placed. The PES is thus extended or updated if the wave function moves into new areas or if a certain area becomes more important. Here, we derive equations for the reduced one-mode density for the time-dependent Hartree (TDH) method and for multiconfiguration time-dependent Hartree (MCTDH) methods, but the TD-ADGA can be used with any time-dependent wave function method as long as a density is available. The TD-ADGA method has been investigated on molecular systems containing single- and double-minimum potentials and on single-mode and multi-mode systems. We explore different approaches to handle the fact that the TD-ADGA involves a PES that changes during the computation and show how results can be obtained that are in very good agreement with results obtained by using an accurate reference PES. Dynamics with TD-ADGA is essentially a black box procedure, where only the initialization of the system and how to compute SPs must be provided. The TD-ADGA thus makes it easier to carry out quantum molecular dynamics and the quasi-direct framework opens up the possibility to compute quantum dynamics accurately for larger molecular systems.

Benchmarking non-adiabatic quantum dynamics using the molecular Tully models

On-the-fly non-adiabatic dynamics methods are becoming more important as tools to characterise the time evolution of a system after absorbing light. These methods, which calculate quantities such as state energies, gradients and interstate couplings at every time step, circumvent the requirement for pre-computed potential energy surfaces. There are a number of different algorithms used, the most common being Tully Surface Hopping (TSH), but all are approximate solutions to the time-dependent Schrödinger equation and benchmarking is required to understand their accuracy and performance. For this, a common set of systems and observables are required to compare them. In this work, we validate the on-the-fly direct dynamics variational multi-configuration Gaussian (DD-vMCG) method using three molecular systems recently suggested by Ibele and Curchod as molecular versions of the Tully model systems used to test one-dimensional non-adiabatic behaviour [Ibele et al., Phys. Chem. Chem. Phys. 2020, 22, 15183-15196]. Parametrised linear vibronic potential energy surfaces for each of the systems were also tested and compared to on-the-fly results. The molecules, which we term the Ibele-Curchod models, are ethene, DMABN and fulvene and the authors used them to test and compare several versions of the Ab Initio Multiple Spawning (AIMS) method alongside TSH. The three systems present different deactivation pathways after excitation to their ππ* bright states. When comparing DD-vMCG to AIMS and TSH, we obtain crucial differences in some cases, for which an explanation is provided by the classical nature and the chosen initial conditions of the TSH simulations.

LVC/MM: A Hybrid Linear Vibronic Coupling/Molecular Mechanics Model with Distributed Multipole-Based Electrostatic Embedding for Highly Efficient Surface Hopping Dynamics in Solution

We present a theoretical framework for a hybrid linear vibronic coupling model electrostatically embedded into a molecular mechanics environment, termed the linear vibronic coupling/molecular mechanics (LVC/MM) method, for the surface hopping including arbitrary coupling (SHARC) molecular dynamics package. Electrostatic embedding is realized through the computation of interactions between environment point charges and distributed multipole expansions (DMEs, up to quadrupoles) that represent each electronic state and transition densities in the diabatic basis. The DME parameters are obtained through a restrained electrostatic potential (RESP) fit, which we extended to yield higher-order multipoles. We also implemented in SHARC a scheme for achieving roto-translational invariance of LVC models as well as a general quantum mechanics/molecular mechanics (QM/MM) interface, an OpenMM interface, and restraining potentials for simulating liquid droplets. Using thioformaldehyde in water as a test case, we demonstrate that LVC/MM can accurately reproduce the solvation structure and energetics of rigid solutes, with errors on the order of 1-2 kcal/mol compared to a BP86/MM reference. The implementation in SHARC is shown to be very efficient, enabling the simulation of trajectories on the nanosecond time scale in a matter of days.

Micro-Solvated DMABN: Excited State Quantum Dynamics and Dual Fluorescence Spectra

In this work, we report a complete analysis by theoretical and spectroscopic methods of the short-time behaviour of 4-(dimethylamino)benzonitrile (DMABN) in the gas phase as well as in cyclohexane, tetrahydrofuran, acetonitrile, and water solution, after excitation to the La state. The spectroscopic properties of DMABN were investigated experimentally using UV absorption and fluorescence emission spectroscopy. The computational study was developed at different electronic structure levels and using the Polarisable Continuum Model (PCM) and explicit solvent molecules to reproduce the solvent environment. Additionally, excited state quantum dynamics simulations in the diabatic picture using the direct dynamics variational multiconfigurational Gaussian (DD-vMCG) method were performed, the largest quantum dynamics "on-the-fly" simulations performed with this method until now. The comparison with fully converged multilayer multiconfigurational time-dependent Hartree (ML-MCTDH) dynamics on parametrised linear vibronic coupling (LVC) potentials show very similar population decays and evolution of the nuclear wavepacket. The ring C=C stretching and three methyl tilting modes are identified as the responsible motions for the internal conversion from the La to the Lb states. No major differences are observed in the ultrafast initial decay in different solvents, but we show that this effect depends strongly on the level of electronic structure used.

Direct nonadiabatic quantum dynamics simulations of the photodissociation of phenol

Gaussian wavepacket methods are becoming popular for the investigation of nonadiabatic molecular dynamics. In the present work, a recently developed efficient algorithm for the Direct Dynamics variational Multi-Configurational Gaussian (DD-vMCG) method has been used to describe the multidimensional photodissociation dynamics of phenol including all degrees of freedom. Full-dimensional quantum dynamic calculations including for the first time six electronic states (¹ππ, 1¹ππ*, 1¹πσ*, 2¹πσ*, 2¹ππ*, 3¹ππ*), along with a comparison to an existing analytical 4-state model for the potential energy surfaces are presented. Including the fifth singlet excited state is shown to have a significant effect on the nonadiabatic photodissociation of phenol to the phenoxyl radical and hydrogen atom. State population and flux analysis from the DD-vMCG simulations of phenol provided further insights into the decay mechanism, confirming the idea of rapid relaxation to the ground state through the ¹ππ/1¹πσ* conical intersection.

Molecular excited states through a machine learning lens

Theoretical simulations of electronic excitations and associated processes in molecules are indispensable for fundamental research and technological innovations. However, such simulations are notoriously challenging to perform with quantum mechanical methods. Advances in machine learning open many new avenues for assisting molecular excited-state simulations. In this Review, we track such progress, assess the current state of the art and highlight the critical issues to solve in the future. We overview a broad range of machine learning applications in excited-state research, which include the prediction of molecular properties, improvements of quantum mechanical methods for the calculations of excited-state properties and the search for new materials. Machine learning approaches can help us understand hidden factors that influence photo-processes, leading to a better control of such processes and new rules for the design of materials for optoelectronic applications.

Improved algorithm for the direct dynamics variational multi-configurational Gaussian method

The Direct Dynamics variational Multi-Configurational Gaussian (DD-vMCG) method provides a fully quantum mechanical solution to the time-dependent Schrödinger equation for the time evolution of nuclei with potential surfaces calculated on-the-fly using a quantum chemistry program. Initial studies have shown its potential for flexible and accurate simulations of non-adiabatic excited-state molecular dynamics. In this paper, we present developments to the DD-vMCG algorithm that improve both its accuracy and efficiency. First, a new, efficient parallel algorithm to control the DD-vMCG database of quantum chemistry points is presented along with improvements to the Shepard interpolation scheme. Second, the use of symmetry in describing the potential surfaces is introduced along with a new phase convention in the propagation diabatization. Benchmark calculations on the allene radical cation including all degrees of freedom then show that the new scheme is able to produce a consistent non-adiabatic coupling vector field. This new DD-vMCG version thus opens the route for effectively and accurately treating complex chemical systems using quantum dynamics simulations.

PySurf: A Framework for Database Accelerated Direct Dynamics

The greatest restriction to the theoretical study of the dynamics of photoinduced processes is computationally expensive electronic structure calculations. Machine learning algorithms have the potential to reduce the number of these computations significantly. Here, PySurf is introduced as an innovative code framework, which is specifically designed for rapid prototyping and development tasks for data science applications in computational chemistry. It comes with powerful Plugin and Workflow engines, which allows intuitive customization for individual tasks. Data is automatically stored through the database framework, which enables additional interpolation of properties in previously evaluated regions of the conformational space. To illustrate the potential of the framework, a code for nonadiabatic surface hopping simulations based on the Landau-Zener algorithm is presented here. Deriving gradients from the interpolated potential energy surfaces allows for full-dimensional nonadiabatic surface hopping simulations using only adiabatic energies (energy only). Simulations of a pyrazine model and ab initio-based calculations of the SO2 molecule show that energy-only calculations with PySurf are able to correctly predict the nonadiabatic dynamics of these prototype systems. The results reveal the degree of sophistication, which can be achieved by the database accelerated energy-only surface hopping simulations being competitive to commonly used semiclassical approaches.

A combined on-the-fly/interpolation procedure for evaluating energy values needed in molecular simulations

We propose an algorithm for molecular dynamics or Monte Carlo simulations that uses an interpolation procedure to estimate potential energy values from energies and gradients evaluated previously at points of a simplicial mesh. We chose an interpolation procedure that is exact for harmonic systems and considered two possible mesh types: Delaunay triangulation and an alternative anisotropic triangulation designed to improve performance in anharmonic systems. The mesh is generated and updated on the fly during the simulation. The procedure is tested on two-dimensional quartic oscillators and on the path integral Monte Carlo evaluation of the HCN/DCN equilibrium isotope effect.

Computation of the S1 ← S0 Vibronic Absorption Spectrum of Formaldehyde by Variational Gaussian Wavepacket and Semiclassical IVR Methods

The vibronic absorption spectrum of the electric dipole forbidden and vibronically allowed S₁(¹ A₂) ← S₀(¹ A₁) transition of formaldehyde is calculated by Gaussian wavepacket and semiclassical methods, along with numerically exact reference calculations, using the potential energy surface of Fu, Shepler, and Bowman ( J. Am. Chem. Soc. 2011, 133, 7957). Specifically, the variational multiconfigurational Gaussian (vMCG) approach and the Herman-Kluk semiclassical initial value representation (HK-SCIVR) are compared to assess the accuracy and convergence of these methods, benchmarked against numerically exact time-dependent wavepacket propagation (TDWP) on the reference potential energy surface. The vMCG calculation is shown to converge quite well with about 100 variationally evolving Gaussian functions and using a local cubic expansion instead of the conventional local harmonic approximation. By contrast, the HK-SCIVR approach with ∼10⁵ trajectories reproduces the vibrationally structured spectral envelope correctly but yields a strongly broadened spectrum. The comparison of the computed absorption spectrum with experiment shows that the relevant vibronic progressions are reasonably reproduced by all computations, but deviations of the order of 10-100 cm^-1 occur, underscoring that both electronic structure calculations and dynamical approaches remain challenging in the calculation of typical small-molecule excited-state spectra by trajectory-based methods.

Curve crossing in a manifold of coupled electronic states: direct quantum dynamics simulations of formamide

Quantum dynamics simulations are an important tool to evaluate molecular behaviour including the, often key, quantum nature of the system. In this paper we present an algorithm that is able to simulate the time evolution of a molecule after photo-excitation into a manifold of states. The direct dynamics variational multi-configurational Gaussian (DD-vMCG) method circumvents the computational bottleneck problems of traditional grid-based methods by computing the potential energy functions on-the-fly, i.e. only where required. Unlike other commonly used direct dynamics methods, DD-vMCG is fully quantum mechanical. Here, the method is combined with a novel on-the-fly diabatisation scheme to simulate the short-time dynamics of the key molecule formamide and its acid analogue formimidic acid. This is a challenging test system due to the nature and large number of excited states, and eight coupled states are included in the calculations. It is shown that the method is able to provide unbiased information on the product channels open after excitation at different energies and demonstrates the potential to be a practical scheme, limited mainly by the quality of the quantum chemistry used to describe the excited states.

Accelerating direct quantum dynamics using graphical processing units

The direct dynamics variational multi-configurational Gaussian (DD-vMCG) method is combined with electronic structure calculations accelerated by Graphical Processing Units (GPUs). This is used to identify GPU acceleration will have a significant effect for both ground and excited state simulations.

Challenges facing an understanding of the nature of low-energy excited states in photosynthesis

While the majority of the photochemical states and pathways related to the biological capture of solar energy are now well understood and provide paradigms for artificial device design, additional low-energy states have been discovered in many systems with obscure origins and significance. However, as low-energy states are naively expected to be critical to function, these observations pose important challenges. A review of known properties of low energy states covering eight photochemical systems, and options for their interpretation, are presented. A concerted experimental and theoretical research strategy is suggested and outlined, this being aimed at providing a fully comprehensive understanding.

Explicit calculation of the excited electronic states of the photosystem II reaction centre

The excited states of the photosystem II reaction centre cofactors have been calculated as a single “supermolecule”. Charge transfer states are shown to be dependent on electrostatic environment.

Basis expansion leaping: a new method to solve the time-dependent Schrödinger equation for molecular quantum dynamics

A wide variety of molecular systems that have recently come into the reach of experimental and theoretical investigation is dominated by quantum phenomena. However, even state of the art quantum propagation techniques are either unsuitable for general application to molecular systems with strong interference and tunneling characteristics or are computationally prohibitive for systems with more than a few degrees of freedom. In this Letter, we introduce a novel quantum propagation technique with wide applicability, controllable accuracy, and efficient utilization of computational resources. Its performance is validated for tunneling and dissociating systems with 1, 2, and 3 degrees of freedom, and the scaling behavior with respect to system dimensionality and requested accuracy is discussed.

Vibrations of H+(D+) in stoichiometric LiNbO3 single crystal

A first principles quantum mechanical calculation of the vibrational energy levels and transition frequencies associated with protons in stoichiometric LiNbO(3) single crystal has been carried out. The hydrogen contaminated crystal has been approximated by a model one obtains by translating a supercell, i.e., a cluster of LiNbO(3) unit cells containing a single H(+) and a Li(+) vacancy. Based on the supercell model an approximate Hamiltonian operator describing vibrations of the proton sublattice embedded in the host crystal has been derived. It is further simplified to a sum of uncoupled Hamiltonian operators corresponding to different wave vectors (ks) and each describing vibrations of a quasi-particle (quasi-proton). The three dimensional (3D) Hamiltonian operator of k=0 has been employed to calculate vibrational levels and transition frequencies. The potential energy surface (PES) entering this Hamiltonian operator has been calculated point wise on a large set of grid points by using density functional theory, and an analytical approximation to the PES has been constructed by non-parametric approximation. Then, the nuclear motion Schrödinger equation has been solved by employing the method of discrete variable representation. It has been found that the (quasi-)H(+) vibrates in a strongly anharmonic PES. Its vibrations can be described approximately as a stretching, and two orthogonal bending vibrations. The theoretically calculated transition frequencies agree within 1% with those experimentally determined, and they have allowed the assignment of one of the hitherto unassigned bands as a combination of the stretching and the bending of lower fundamental frequency.