Semiclassical Molecular Dynamics Research Articles

Semiclassical Nonadiabatic Molecular Dynamics Using Linearized Pair-Density Functional Theory.

Nonadiabatic molecular dynamics is an effective method for modeling nonradiative decay in electronically excited molecules. Its accuracy depends strongly on the quality of the potential energy surfaces, and its affordability for long direct-dynamic simulations with adequate ensemble averaging depends strongly on the cost of the required electronic structure calculations. Linearized pair-density functional theory (L-PDFT) is a recently developed post-self-consistent-field multireference method that can model potential energy surfaces with an accuracy similar to expensive multireference perturbation theories but at a computational cost similar to the underlying multiconfiguration self-consistent field method. Here, we integrate the SHARC dynamics and PySCF electronic structure code to utilize L-PDFT for electronically nonadiabatic calculations and use the combined programs to study the photoisomerization reaction of cis-azomethane. We show that L-PDFT is able to successfully simulate the photoisomerization without crashes, and it yields results similar to the more expensive extended multistate complete active space second-order perturbation theory. This shows that L-PDFT can model internal conversion, and it demonstrates its promise for broader photodynamics applications.

Unveiling the electron-induced ionization cross sections and fragmentation mechanisms of 3,4-dihydro-2H-pyran.

The interactions of electrons with molecular systems under various conditions are essential to interdisciplinary research fields extending over the fundamental and applied sciences. In particular, investigating electron-induced ionization and dissociation of molecules may shed light on the radiation damage to living cells, the physicochemical processes in interstellar environments, and reaction mechanisms occurring in combustion or plasma. We have, therefore, studied electron-induced ionization and dissociation of the gas phase 3,4-dihydro-2H-pyran (DHP), a cyclic ether appearing to be a viable moiety for developing efficient clinical pharmacokinetics and revealing the mechanisms of biofuel combustion. The mass spectra in the m/z = 10-90 mass range were measured at several different energies of the ionizing electron beam using mass spectrometry. The mass spectra of DHP at the same energies were simulated using on-the-fly semi-classical molecular dynamics (MD) within the framework of the QCxMS formalism. The MD settings were suitably adjusted until a good agreement with the experimental mass spectra intensities was achieved, thus enabling a reliable assignment of cations and unraveling the plausible fragmentation channels. Based on the measurement of the absolute total ionization cross section of DHP (18.1 ± 0.9) × 10-16cm2 at 100eV energy, the absolute total and partial ionization cross sections of DHP were determined in the 5-140eV electron energy. Moreover, a machine learning algorithm that was trained with measured cross sections from 25 different molecules was used to predict the total ionization cross section for DHP. Comparison of the machine learning simulation with the measured data showed acceptable agreement, similar to that achieved in past predictions of the algorithm.

The Fate of the Formic Acid Proton on the Anatase TiO2(101) Surface.

As a prototype adsorption reaction of gas Brønsted acid on oxides, we study the adsorption of formic acid on anatase. We perform infrared spectroscopy measurements of adsorbed HCOOH and HCOOD on TiO2 nanopowders, from 13 K up to room temperature in an ultra-high vacuum chamber. We assign the IR signals via computed spectra from nuclear quantum dynamics simulations using our divide-and-conquer semiclassical ab initio molecular dynamics method. The acid proton forms an extraordinarily short and strong hydrogen bond with the surface oxygen. The strength of this hydrogen bond, that compares to H bonds in ice at high pressures, is at the root of a substantial redshift with respect to the typical free OH stretching frequency, which eludes its straightforward detection.

Computational modeling of CH4 and CO2 adsorption on monolayer graphenylene: Implications for optoelectronic properties and hydrogen production

Graphenylene (GP) is a two–dimensional carbon allotrope with a hexagonal lattice structure containing periodic pores. The unique arrangement of GP offers potential applications in electronics, optoelectronics, energy storage, and gas separation. Specifically, its advantageous electronic and optical properties, make it a promising candidate for hydrogen production and advanced electronic devices. In this study, we employ a computational chemistry–based modeling approach to investigate the adsorption mechanisms of CH4 and CO2 on monolayer GP, with a specific focus on their effects on optical adsorption and electrical transport properties at room temperature. To simulate the adsorption dynamics as closely as possible to experimental conditions, we utilize the self–consistent charge tight–binding density functional theory (SCC–DFTB). Through semi–classical molecular dynamics (MD) simulations, we observe the formation of H2 molecules from the dissociation of CH4 and the formation of CO + O species from carbon dioxide molecules. This provides insights into the adsorption and dispersion mechanisms of CH4 and CO2 on GP. Furthermore, we explore the impact of molecular adsorption on optical absorption properties. Our results demonstrate that CH4 and CH2 affects drastically the optical adsorption of GP, while CO2 does not significantly affect the optical properties of the two–dimensional material. To analyze electron transport, we employ the open–boundary non–equilibrium Green's function method. By studying the conductivity of GP and graphene under voltage bias up to 300 mV, we gain valuable insights into the electrical transport properties of GP under optical absorption conditions. The findings from our computational modeling approach might contribute to a deeper understanding of the potential applications of GP in hydrogen production and advanced electronic devices.

Semiclassical approach to spin dynamics of a ferromagnetic S = 1 chain

Abstract Motivated by recent experimental progress on the quasi-one-dimensional quantum magnet NiNb2O6, we study the spin dynamics of an S = 1 ferromagnetic Heisenberg chain with single-ion anisotropy by using a semiclassical molecular dynamics approach. This system undergoes a quantum phase transition from a ferromagnetic to a paramagnetic state under a transverse magnetic field, and the magnetic response reflecting this transition is well described by our semiclassical method. We show that at low temperature the transverse component of the dynamical structure factor depicts clearly the magnon dispersion, and the longitudinal component exhibits two continua associated with single- and two-magnon excitations, respectively. These spin excitation spectra show interesting temperature dependence as effects of magnon interactions. Our findings shed light on the experimental detection of spin excitations in a large class of quasi-one-dimensional magnets.

Non-Hermitian molecular dynamics simulations of exciton-polaritons in lossy cavities.

The observation that materials can change their properties when placed inside or near an optical resonator has sparked a fervid interest in understanding the effects of strong light-matter coupling on molecular dynamics, and several approaches have been proposed to extend the methods of computational chemistry into this regime. Whereas the majority of these approaches have focused on modeling a single molecule coupled to a single cavity mode, changes to chemistry have so far only been observed experimentally when very many molecules are coupled collectively to multiple modes with short lifetimes. While atomistic simulations of many molecules coupled to multiple cavity modes have been performed with semi-classical molecular dynamics, an explicit description of cavity losses has so far been restricted to simulations in which only a very few molecular degrees of freedom were considered. Here, we have implemented an effective non-Hermitian Hamiltonian to explicitly treat cavity losses in large-scale semi-classical molecular dynamics simulations of organic polaritons and used it to perform both mean-field and surface hopping simulations of polariton relaxation, propagation, and energy transfer.

Modeling pH-Dependent Biomolecular Photochemistry.

The tuning mechanism of pH can be extremely challenging to model computationally in complex biological systems, especially with respect to the photochemical properties. This article reports a protocol aimed at modeling pH-dependent photodynamics using a combination of constant-pH molecular dynamics and semiclassical nonadiabatic molecular dynamics simulations. With retinal photoisomerization in Anabaena sensory rhodopsin (ASR) as a testbed, we show that our protocol produces pH-dependent photochemical properties, such as the isomerization quantum yield or decay rates. We decompose our results into single-titrated residue contributions, identifying some key tuning amino acids. Additionally, we assess the validity of the single protonation state picture to represent the system at a given pH and propose the most populated protein charge state as a compromise between cost and accuracy.

ANT 2023: A program for adiabatic and nonadiabatic trajectories

ANT 2023 is a program for quasiclassical and semiclassical trajectories, both single-surface trajectories for which the Born-Oppenheimer approximation is valid and multi-surface calculations with electronic state changes, i.e., for electronically adiabatic and electronically nonadiabatic trajectories. There are several methods available for multisurface problems: surface hopping with or without time uncertainty and with or without decoherence, semiclassical Ehrenfest, self-consistent decay of mixing, and coherent switching with decay of mixing (CSDM). The potential surface for single-surface problems may be an analytic potential function supplied by the user, or one may use direct dynamics. To use the adiabatic representation (i.e., electronically adiabatic basis functions) for electronically nonadiabatic dynamics, the user may either provide the adiabatic surfaces and nonadiabatic couplings by direct dynamics, or the program may calculate them from diabatic surfaces and diabatic couplings, which are usually analytic. One can also use analytic fits to the surfaces and couplings to carry out calculations entirely in the diabatic representation. The curvature-driven approximation is available as an option for use with CSDM and trajectory surface hopping. The ANT 2023 program is especially recommended for calculations with analytic potential energy surfaces and couplings because of its high efficiency for such calculations. Program SummaryProgram Title: ANT 2023, revision BCPC Library link to program files:https://doi.org/10.17632/gr74wwckcp.1Developer's repository link: https://doi.org/10.5281/zenodo.10011563Licensing Provisions: Apache-2.0Programming language: Fortran 90Interfaces to other programs: For direct dynamics, the code has well tested interfaces to the Gaussian 09 [1], Gaussian 16 [2], MOPACmn [3], and Molpro [4] electronic structure codes. The interfaces are included, but these electronic structure programs must be obtained separately. Users may also provide their own interfaces to other electronic structure programs.[1] M. J. Frisch et al. Gaussian 09, Gaussian, Inc., Wallingford CT, 2016.[2] M. J. Frisch et al., Gaussian 16, Gaussian, Inc., Wallingford CT, 2016.[3] J. J. P. Stewart et al., MOPAC 5.022mn, University of Minnesota, Minneapolis, MN, 2015; Zenodo. https://doi.org/10.5281/zenodo.8153079[4] H.-J. Werner et al., J. Chem. Phys. 152 (2020) 144107.Nature of problem: classical, quasiclassical, or semiclassical molecular dynamics for electronically adiabatic or electronically nonadiabatic collisions and unimolecular processes of atoms, molecules, and clustersSolution method: The program integrates the classical or semiclassical equations of motion for the motion of the nuclei. ANT 2023 can use anharmonic initial-state selection for initial states of diatomics or harmonic initial-state selection for general polyatomics. The equations of motion for electronically adiabatic processes are Hamilton's equations with potential energy functions corresponding to the Born-Oppenheimer approximation. The equations of motion for electronically nonadiabatic processes (non-Born-Oppenheimer processes) include trajectory surface hopping (fewest switches with or without time uncertainty and with or without decoherence) or several possible mean-field methods, including the semiclassical Ehrenfest method (also called the time-dependent Hartree method), self-consistent decay of mixing, and coherent switching with decay of mixing. Electronically nonadiabatic processes may be treated using either the adiabatic representation or a diabatic representation, and the potential surfaces and couplings may be obtained from analytic fits in a user-supplied subroutine or by direct dynamics, which involves requesting energies, gradients, and couplings as needed from a quantum chemistry electronic structure program. In the adiabatic representation, the coupling may be nonadiabatic coupling vectors (based on gradients of electronic wave functions), time-derivative couplings (based on overlaps of electronic wave functions), or curvature-driven couplings (based on the curvatures of the adiabatic potential energy surfaces).Additional comments: The ANT 2023 program authors are Y. Shu, L. Zhang, J. Zheng, Z. H. Li, A. W. Jasper, D. A. Bonhommeau, R. Valero, R. Meana-Pañeda, S. L. Mielke, Z. Varga, and D. G. Truhlar [5]; the authors of the article [6] describing ANT 2023 are Y. Shu, L. Zhang, and D. G. Truhlar.

The radiative surface hopping (RSH) algorithm: Capturing fluorescence events in molecular systems within a semi-classical non-adiabatic molecular dynamics framework.

In this article, we present the radiative surface hopping algorithm, which enables modeling fluorescence within a semi-classical non-adiabatic molecular dynamics framework. The algorithm has been tested for the photodeactivation dynamics of trans-4-dimethylamino-4'-cyanostilbene (DCS). By treating on equal footing the radiative and non-radiative processes, our method allows us to attain a complete molecular movie of the excited-state deactivation. Our dynamics rely on a semi-empirical quantum mechanical/molecular mechanical Hamiltonian and have been run for hundreds of picoseconds, both in the gas phase and in isopropyl ether. The proposed approach successfully captures the first fluorescence processes occurring in DCS, and it succeeds in reproducing the experimental fluorescence lifetime and quantum yield, especially in the polar solvent. The analysis of the geometrical features of the emissive species during the dynamics discards the hypothesis of a twisted intramolecular charge transfer state to be responsible for the dual emission observed experimentally in some polar solvents. In a nutshell, our method opens the way for theoretical studies on early fluorescence events occurring up to hundreds of picoseconds in molecular systems.

Commensurate and incommensurate double moiré interference in twisted trilayer graphene

Twisted graphene multi-layers have been recently demonstrated to share several correlation-driven behaviours with twisted bilayer graphene. In general, the van Hove singularities (VHSs) can be used as a proxy of the tendency for correlated behaviours. In this paper, we adopt an atomistic method by combining tight-binding method with the semi-classical molecular dynamics to investigate the electronic structures of twisted trilayer graphene (TTG) with two independent twist angles. The two independent twist angles can lead to the interference of the moir\'e patterns forming a variety of commensurate/incommensurate complex supermoir\'e patterns. In particular, the lattice relaxation, twist angle and angle disorder effects on the VHS are discussed. We find that the lattice relaxation significantly influence the position and magnitude of the VHSs. In the supermoir\'e TTG, the moir\'e interference provides constructive or destructive effects depending on the relative twist angle. By modulating the two independent twist angles, novel superstructures, for instance, the Kagome-like lattice, could constructed via the moir\'e pattern. Moreover, we demonstrate that a slight change in twist angles (angle disorder) provides a significant suppression of the peak of the VHSs. Apart from the moir\'e length, the evolution of the VHSs and the LDOS mapping in real space could be used to identify the twist angles in the complicated TTG. In practice, our work could provide a guide for exploring the flat band behaviours in the supermoir\'e TTG experimentally.

Modeling of second sound in carbon nanostructures

The study of thermal transport in low-dimensional materials has attracted a lot of attention recently after discovery of high thermal conductivity of graphene. Here we study numerically phonon transport in low-dimensional carbon structures being interested in the hydrodynamic regime revealed through the observation of second sound. We demonstrate that correct numerical modeling of such two-dimensional systems requires semi-classical molecular dynamics simulations of temperature waves that take into account quantum statistics of thermalized phonons. We reveal that second sound can be attributed to the maximum group velocity of bending optical oscillations of carbon structures, and the hydrodynamic effects disappear for $T>200$K, being replaced by diffusive dynamics of thermal waves. Our numerical results suggest that the velocity of second sound in such low-dimensional structures is about 6 km/s, and the hydrodynamic effects are manifested stronger in carbon nanotubes rather than in carbon nanoribbons.

An atomistic approach for the structural and electronic properties of twisted bilayer graphene-boron nitride heterostructures

Twisted bilayer graphene (TBG) has taken the spotlight in the condensed matter community since the discovery of correlated phases. In this work, we study heterostructures of TBG and hexagonal boron nitride (hBN) using an atomistic tight-binding model together with semi-classical molecular dynamics to consider relaxation effects. The hBN substrate has significant effects on the band structure of TBG even in the case where TBG and hBN are not aligned. Specifically, the substrate induces a large mass gap and strong pseudo-magnetic fields that break the layer degeneracy. Interestingly, such degeneracy can be recovered with a second hBN layer. Finally, we develop a continuum model that describes the tight-binding band structure. Our results show that a real-space tight-binding model in combination with semi-classical molecular dynamics is a powerful tool to study the electronic properties of moiré heterostructures, and to explain experimental results in which the effect of the substrate plays an important role.

Photoinduced bond oscillations in ironpentacarbonyl give delayed synchronous bursts of carbonmonoxide release

Early excited state dynamics in the photodissociation of transition metal carbonyls determines the chemical nature of short-lived catalytically active reaction intermediates. However, time-resolved experiments have not yet revealed mechanistic details in the sub-picosecond regime. Hence, in this study the photoexcitation of ironpentacarbonyl Fe(CO)5 is simulated with semi-classical excited state molecular dynamics. We find that the bright metal-to-ligand charge-transfer (MLCT) transition induces synchronous Fe-C oscillations in the trigonal bipyramidal complex leading to periodically reoccurring release of predominantly axial CO. Metaphorically the photoactivated Fe(CO)5 acts as a CO geyser, as a result of dynamics in the potential energy landscape of the axial Fe-C distances and non-adiabatic transitions between manifolds of bound MLCT and dissociative metal-centered (MC) excited states. The predominant release of axial CO ligands and delayed release of equatorial CO ligands are explained in a unified mechanism based on the σ*(Fe-C) anti-bonding character of the receiving orbital in the dissociative MC states.

Quantum Vibrational Spectroscopy of Explicitly Solvated Thymidine in Semiclassical Approximation.

In this paper, we demonstrate the possibility to perform spectroscopy simulations of solvated biological species taking into consideration quantum effects and explicit solvation. We achieve this goal by interfacing our recently developed divide-and-conquer approach for semiclassical initial value representation molecular dynamics with the polarizable AMOEBABIO18 force field. The method is applied to the study of solvation of the thymidine nucleoside in two different polar solvents, water and N,N-dimethylformamide. Such systems are made of up to 2476 atoms. Experimental evidence concerning the different behavior of thymidine in the two solvents is well reproduced by our study, even though quantitative estimates are hampered by the limited accuracy of the classical force field employed. Overall, this study shows that semiclassically approximate quantum dynamical studies of explicitly solvated biological systems are both computationally affordable and insightful.

Theoretical-computational modelling of the vibrational relaxation of small inorganic species in condensed phase

In this paper we apply a recently developed theoretical-computational procedure for modelling the Vibrational Energy Relaxation (VER) of solvated chromophores represented by small inorganic species. In particular we focus our attention on four different systems, all experimentally well characterized: aqueous cyanide ion and aqueous azide ion, nitrogen dioxide and water both in chloroform. Our method essentially reconstructs the whole vibrational relaxation kinetics (with the chromophore in the ground electronic state) by: (i) determining, through semiclassical Molecular Dynamics (MD) simulation, the electric field (perturbation) produced onto the chromophore by the solvent atom motions; (ii) using the electrostatic perturbation for directly determining the chromophore quantum vibrational dynamics; (iii) calculating the rate constant for the vibrational relaxation as occurring without quantum–classical energy exchange; (iv) introducing the latter effect, a posteriori, hence obtaining the actual relaxation kinetics. Our results, in satisfactory agreement with the available experimental data, show that the VER mechanism is almost entirely determined by the fluctuating perturbation field as produced by the time-dependent motions of the environment atoms (in these cases the solvent) similarly to the well-known effects of the electromagnetic wave causing absorption/emission processes.

On-the-fly adiabatically switched semiclassical initial value representation molecular dynamics for vibrational spectroscopy of biomolecules.

Semiclassical (SC) vibrational spectroscopy is a technique capable of reproducing quantum effects (such as zero-pointenergies, quantum resonances, and anharmonic overtones) from classical dynamics runs even in the case of very large dimensional systems. In a previous study [Conte et al. J. Chem. Phys. 151, 214107 (2019)], a preliminary sampling based on adiabatic switching has been shown to be able to improve the precision and accuracy of semiclassical results for challenging model potentials and small molecular systems. In this paper, we investigate the possibility to extend the technique to larger (bio)molecular systems whose dynamics must be integrated by means of ab initio "on-the-fly" calculations. After some preliminary tests on small molecules, we obtain the vibrational frequencies of glycine improving on pre-existing SC calculations. Finally, the new approach is applied to 17-atom proline, an amino acid characterized by a strong intramolecular hydrogen bond.

Temperature-dependent thermal transport of single molecular junctions from semiclassical Langevin molecular dynamics

Thermal conductance of single molecular junctions at room temperature has been measured recently using picowatt-resolution scanning probes. However, fully understanding thermal transport in a much wider temperature range is needed for the exploration of energy transfer at single-molecular limit and the development of single-molecular devices. Here, employing a semiclassical Langevin molecular dynamics method, a comparative study is performed on the thermal transport of an alkane chain between Au and graphene electrodes, respectively. We illustrate the different roles of quantum statistics and anharmonic interaction in the two types of junctions. For a graphene junction, quantum statistics is essential at room temperature, while the anharmonic interaction is negligible. For a Au junction, it is the other way. Our study paves the way for theoretically understanding thermal transport of realistic single-molecular junctions in the full temperature range by including both quantum statistics and anharmonic interaction within one theoretical framework.

Unsupervised Machine Learning Neural Gas Algorithm for Accurate Evaluations of the Hessian Matrix in Molecular Dynamics.

The Hessian matrix of the potential energy of molecular systems is employed not only in geometry optimizations or high-order molecular dynamics integrators but also in many other molecular procedures, such as instantaneous normal mode analysis, force field construction, instanton calculations, and semiclassical initial value representation molecular dynamics, to name a few. Here, we present an algorithm for the calculation of the approximated Hessian in molecular dynamics. The algorithm belongs to the family of unsupervised machine learning methods, and it is based on the neural gas idea, where neurons are molecular configurations whose Hessians are adopted for groups of molecular dynamics configurations with similar geometries. The method is tested on several molecular systems of different dimensionalities both in terms of accuracy and computational time versus calculating the Hessian matrix at each time-step, that is, without any approximation, and other Hessian approximation schemes. Finally, the method is applied to the on-the-fly, full-dimensional simulation of a small synthetic peptide (the 46 atom N-acetyl-l-phenylalaninyl-l-methionine amide) at the level of DFT-B3LYP-D/6-31G* theory, from which the semiclassical vibrational power spectrum is calculated.

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.

Machine learning for vibrational spectroscopy via divide-and-conquer semiclassical initial value representation molecular dynamics with application to N-methylacetamide.

A machine learning algorithm for partitioning the nuclear vibrational space into subspaces is introduced. The subdivision criterion is based on Liouville's theorem, i.e., the best preservation of the unitary of the reduced dimensionality Jacobian determinant within each subspace along a probe full-dimensional classical trajectory. The algorithm is based on the idea of evolutionary selection, and it is implemented through a probability graph representation of the vibrational space partitioning. We interface this customized version of genetic algorithms with our divide-and-conquer semiclassical initial value representation method for the calculation of molecular power spectra. First, we benchmark the algorithm by calculating the vibrational power spectra of two model systems, for which the exact subspace division is known. Then, we apply it to the calculation of the power spectrum of methane. Exact calculations and full-dimensional semiclassical spectra of this small molecule are available and provide an additional test of the accuracy of the new approach. Finally, the algorithm is applied to the divide-and-conquer semiclassical calculation of the power spectrum of 12-atom trans-N-methylacetamide.