Integrating diffusion maps with umbrella sampling: application to alanine dipeptide

Nonlinear dimensionality reduction techniques can be applied to molecular simulation trajectories to systematically extract a small number of variables with which to parametrize the important dynamical motions of the system. For molecular systems exhibiting free energy barriers exceeding a few k(B)T, inadequate sampling of the barrier regions between stable or metastable basins can lead to a poor global characterization of the free energy landscape. We present an adaptation of a nonlinear dimensionality reduction technique known as the diffusion map that extends its applicability to biased umbrella sampling simulation trajectories in which restraining potentials are employed to drive the system into high free energy regions and improve sampling of phase space. We then propose a bootstrapped approach to iteratively discover good low-dimensional parametrizations by interleaving successive rounds of umbrella sampling and diffusion mapping, and we illustrate the technique through a study of alanine dipeptide in explicit solvent.     

On the role of inherent structures in glass-forming materials: I. The vitrification process

In this work, we investigate the role of inherent structures in the vitrification process of glass-forming materials by using a two-component Lennard-Jones mixture. We start with a simplified model that describes the dynamics of the atomistic system as a Poisson process consisting of a series of transitions from one potential energy minimum (inherent structure) to another, the rate of individual transitions being described by a first-order kinetic law. We investigate the validity of this model by comparing the mean square displacement resulting from atomistic molecular dynamics (MD) trajectories with the corresponding mean square displacement based on inherent structures. Furthermore, in the case of vitrification via stepwise cooling, we identify the role of the potential energy landscape in determining the properties of the resulting glass. Interestingly, the cooling rate is not sufficient to define the resulting glass in a stepwise cooling process, because the time spent by the system at different temperatures (length of the steps) has a highly nonlinear impact on the properties of the resulting glass. In contrast to previous investigations of supercooled liquids, we focus on a range of temperatures close to and below the glass transition temperature, where the use of MD is incapable of producing equilibrated samples of the metastable supercooled state. Our aim is to develop a methodology that enables mapping the dynamics under these conditions to a coarse-grained first-order kinetic model based on the Poisson process approximation. This model can be used in order to extend our dynamical sampling ability to much broader time scales and therefore allow us to create computer glasses with cooling rates closer to those used experimentally. In a continuation to this work, we provide the mathematical formulation for lifting the coarse-grained Poisson process model and reproducing the full dynamics of the atomistic system.     

Computational study of the one and two dimensional infrared spectra of a vibrational mode strongly coupled to its environment: beyond the cumulant and condon approximations

The effect of the commonly employed Condon and second-order cumulant approximations on one- and two-dimensional infrared spectra is examined in the case of a vibrational mode which is strongly coupled to its environment. The analysis is performed within the context of the hydrogen stretch of a moderately strong hydrogen-bonded complex dissolved in a dipolar liquid. The IR spectra are calculated using an adiabatic mixed quantum-classical approach that treats the hydrogen quantum-mechanically, while the remaining degrees of freedom are treated classically. While the cumulant and Condon treatments are seen to produce extremely broad and rather structureless spectra, the non-Condon spectra are found to consist of several relatively narrow bands that can be traced back to subsets of bath configurations with large transition dipole moments. Thus, although the cumulant and Condon approximations can capture some general qualitative spectral trends and are able to reproduce some highly averaged quantities such as the photon-echo peak shift, they fail to reproduce many important features of the spectra. We show that the great sensitivity of the transition dipole moment to the bath configuration provides new means for decongesting the spectra, probing statistically unfavorable bath configurations, and obtaining unique information regarding the dynamics of individual subsets of bath configurations and of the rates of transitions between them.     

Anharmonic infrared and Raman spectra in Car-Parrinello molecular dynamics simulations

The infrared and Raman spectra of naphthalene crystal with inclusion of anharmonic effects have been calculated by adopting the generalized variational density functional perturbation theory in the framework of Car-Parrinello molecular dynamics simulations. The computational approach has been generalized for cells of arbitrary shape. The intermolecular interactions have been analyzed with and without the van der Waals corrections, showing the importance of such interactions in the naphthalene crystal to reproduce the structural, dynamical, and spectroscopic properties.     

Two-dimensional Raman and infrared vibrational spectroscopy for a harmonic oscillator system nonlinearly coupled with a colored noise bath

Multidimensional vibrational response functions of a harmonic oscillator are reconsidered by assuming nonlinear system-bath couplings. In addition to a standard linear-linear (LL) system-bath interaction, we consider a square-linear (SL) interaction. The LL interaction causes the vibrational energy relaxation, while the SL interaction is mainly responsible for the vibrational phase relaxation. The dynamics of the relevant system are investigated by the numerical integration of the Gaussian-Markovian Fokker-Planck equation under the condition of strong couplings with a colored noise bath, where the conventional perturbative approach cannot be applied. The response functions for the fifth-order nonresonant Raman and the third-order infrared (or equivalently the second-order infrared and the seventh-order nonresonant Raman) spectra are calculated under the various combinations of the LL and the SL coupling strengths. Calculated two-dimensional response functions demonstrate that those spectroscopic techniques are very sensitive to the mechanism of the system-bath couplings and the correlation time of the bath fluctuation. We discuss the primary optical transition pathways involved to elucidate the corresponding spectroscopic features and to relate them to the microscopic sources of the vibrational nonlinearity induced by the system-bath interactions. Optical pathways for the fifth-order Raman spectroscopies from an "anisotropic" medium were newly found in this study, which were not predicted by the weak system-bath coupling theory or the standard Brownian harmonic oscillator model.     

Assigning vibrational spectra of ferryl-oxo intermediates of cytochrome C oxidase by periodic orbits and molecular dynamics

Complexity is inherent in biological molecules not only because of the large number of atoms but also because of their nonlinear interactions responsible for chaotic behaviours, localized motions, and bifurcation phenomena. Thus, versatile spectroscopic techniques have been invented to achieve temporal and spatial resolution to minimize the uncertainties in assigning the spectra of complex molecules. Can we associate spectral lines to specific chemical bonds or species in a large molecule? Can energy stay localized in a bond for a substantial period of time to leave its spectroscopic signature? These longstanding problems are investigated by studying the resonance Raman spectra of ferryl-oxo intermediates of cytochrome c oxidase. The difference spectra of isotopically substituted ferryl oxygen ((16)O minus (18)O) in the cytochrome c oxidase recorded in several laboratories show one or two prominent positive peaks which have not been completely elucidated yet. By applying the hierarchical methods of nonlinear mechanics, and particularly the study of periodic orbits in the active site of the enzyme, in conjunction with molecular dynamics calculations of larger systems which include the embraced active site by the protein and selected protonated/deprotonated conformations of amino acids, we translate the spectral lines to molecular motions. It is demonstrated that for the active site stable periodic orbits exist for a substantial energy range. Families of periodic orbits which are associated with the vibrations of Fe(IV)=O bond mark the regions of phase space where nearby trajectories remain localized, as well as assign the spectral bands of the active site in the protein matrix. We demonstrate that proton movement adjacent to active site, which occurs during the P --> F transition, can lead to significant perturbations of the Fe(IV)=O isotopic difference vibrational spectra in cytochrome c oxidase, without a change in oxidation state of the metal sites. This finding links spectroscopic characteristics to protonation events occurring during enzymatic turnover.     

Marginal states in a cubic autocatalytic reaction

Marginal steady state belongs to a special class of states in nonlinear dynamics. To realize this state we consider a cubic autocatalytic reaction A + 2B → 3B in a continuous-stirred-tank-reactor, where the flow rate of the reactant A can be controlled to manipulate the dynamical behavior of the open system. We demonstrate that when the flow rate is weakly noisy the autocatalytic reaction admits of a steady state which is marginal in nature and is surrounded by infinite number of periodic trajectories. When the uncatalyzed reaction A → B is included in the reaction scheme, there exists a marginal steady state which is a critical state corresponding to the point of transition between the flow branch and the equilibrium branch, similar to gas-liquid critical point of transition. This state loses its stability in the weak noise limit.     

Molecular latent space simulators

Small integration time steps limit molecular dynamics (MD) simulations to millisecond time scales. Markov state models (MSMs) and equation-free approaches learn low-dimensional kinetic models from MD simulation data by performing configurational or dynamical coarse-graining of the state space. The learned kinetic models enable the efficient generation of dynamical trajectories over vastly longer time scales than are accessible by MD, but the discretization of configurational space and/or absence of a means to reconstruct molecular configurations precludes the generation of continuous atomistic molecular trajectories. We propose latent space simulators (LSS) to learn kinetic models for continuous atomistic simulation trajectories by training three deep learning networks to (i) learn the slow collective variables of the molecular system, (ii) propagate the system dynamics within this slow latent space, and (iii) generatively reconstruct molecular configurations. We demonstrate the approach in an application to Trp-cage miniprotein to produce novel ultra-long synthetic folding trajectories that accurately reproduce atomistic molecular structure, thermodynamics, and kinetics at six orders of magnitude lower cost than MD. The dramatically lower cost of trajectory generation enables greatly improved sampling and greatly reduced statistical uncertainties in estimated thermodynamic averages and kinetic rates.

Latent space simulators learn kinetic models for atomistic simulations and generate novel trajectories at six orders of magnitude lower cost.     

Ultrafast Structural Dynamics of Biomolecules Examined by Multiple‐Mode 2D IR Spectroscopy: Anharmonically Coupled Motions are in Harmony

Quantum chemical computations, molecular dynamics simulations, and linear and nonlinear infrared spectral simulations are carried out for four representative biomolecules: cellobiose, alanine tripeptide, L ‐α‐glycerylphosphorylethanolamine, and the DNA base monomer guanine. Anharmonic transition frequencies and anharmonicities for the molecules in vacuum are evaluated. Instantaneous normal‐mode analysis is performed and the vibrational frequency distribution correlations are examined for the molecules solvated in TIP3P water. Many local and regional motions of the biomolecules are predicted to be anharmonically coupled and their vibrational frequencies are predicted to be largely correlated. These coupled and correlated vibrational motions can be easily visualized by pairwise cross peaks in the femtosecond broadband two‐dimensional infrared (2D IR) spectra, which are simulated using time‐domain third‐order nonlinear response functions. A network of distinctive spectral profiles of the 2D IR cross peaks, including peak orientations and positive and negative signal patterns, are shown to be intimately connected with the couplings and correlations. The results show that the vibrational couplings and correlations, driven by solvent interactions and also by intrinsic vibrational interactions, are vibrational mode dependent and thus chemical group dependent, and form the structural and dynamical basis of the anharmonic vibrators that are ubiquitous in biomolecules.     

Interference and quantization in semiclassical response functions

Application of the Herman-Kluk semiclassical propagator to the calculation of spectroscopic response functions for anharmonic oscillators has demonstrated the quantitative accuracy of these approximate dynamics. In this approach, spectroscopic response functions are expressed as multiple phase-space integrals over pairs of classical trajectories and their associated stability matrices. Here we analyze the Herman-Kluk semiclassical approximation to a linear response function and determine the origin of the capacity of this method to reproduce quantum effects in a response function from classical dynamical information. Our analysis identifies those classical trajectories that contribute most significantly to the response function on different time scales. This finding motivates a procedure for computing the linear response function in which the interference between pairs of classical trajectories is treated approximately, resulting in an integral over a single average trajectory, as in a purely classical calculation.     

Nonequilibrium molecular dynamics simulations with a backward-forward trajectories sampling for multidimensional infrared spectroscopy of molecular vibrational modes

A full molecular dynamics (MD) simulation approach to calculate multidimensional third-order infrared (IR) signals of molecular vibrational modes is proposed. Third-order IR spectroscopy involves three-time intervals between three excitation and one probe pulses. The nonequilibrium MD (NEMD) simulation allows us to calculate molecular dipoles from nonequilibrium MD trajectories for different pulse configurations and sequences. While the conventional NEMD approach utilizes MD trajectories started from the initial equilibrium state, our approach does from the intermediate state of the third-order optical process, which leads to the doorway-window decomposition of nonlinear response functions. The decomposition is made before the second pump excitation for a two-dimensional case of IR photon echo measurement, while it is made after the second pump excitation for a three-dimensional case of three-pulse IR photon echo measurement. We show that the three-dimensional IR signals are efficiently calculated by using the MD trajectories backward and forward in time for the doorway and window functions, respectively. We examined the capability of the present approach by evaluating the signals of two- and three-dimensional IR vibrational spectroscopies for liquid hydrogen fluoride. The calculated signals might be explained by anharmonic Brownian model with the linear-linear and square-linear system-bath couplings which was used to discuss the inhomogeneous broadening and dephasing mechanism of vibrational motions. The predicted intermolecular librational spectra clearly reveal the unusually narrow inhomogeneous linewidth due to the one-dimensional character of HF molecule and the strong hydrogen bond network.     

Transition state theory for laser-driven reactions

Recent developments in transition state theory brought about by dynamical systems theory are extended to time-dependent systems such as laser-driven reactions. Using time-dependent normal form theory, the authors construct a reaction coordinate with regular dynamics inside the transition region. The conservation of the associated action enables one to extract time-dependent invariant manifolds that act as separatrices between reactive and nonreactive trajectories and thus make it possible to predict the ultimate fate of a trajectory. They illustrate the power of our approach on a driven Henon-Heiles system, which serves as a simple example of a reactive system with several open channels. The present generalization of transition state theory to driven systems will allow one to study processes such as the control of chemical reactions through laser pulses.     

2D regional correlation analysis of single-molecule time trajectories

We report a new approach of 2D regional correlation analysis capable of analyzing fluctuation dynamics of complex multiple correlated and anticorrelated fluctuations under a noncorrelated noise background. Using this new method, by changing and scanning the start time and end time along a pair of fluctuation trajectories, we are able to map out any defined segments along the fluctuation trajectories and determine whether they are correlated, anticorrelated, or noncorrelated; after which, a cross-correlation analysis can be applied for each specific segment to obtain a detailed fluctuation dynamics analysis. We specifically discuss an application of this approach to analyze single-molecule fluorescence resonance energy transfer (FRET) fluctuation dynamics where the fluctuations are often complex, although this approach can be useful for analyzing other types of fluctuation dynamics of various physical variables as well.     

Impact of finite-temperature and condensed-phase effects on theoretical X-ray absorption spectra of transition metal complexes

The impact of condensed-phase and finite-temperature effects on the theoretical X-ray absorption spectra of transition metal complexes is assessed. The former are included in terms of the all-electron Gaussian and augmented plane-wave approach, whereas the latter are taken into account by extensive ensemble averaging along second-generation Car–Parrinello ab initio molecular dynamics trajectories. We find that employing the periodic boundary conditions and including finite-temperature effects systematically improves the agreement between our simulated X-ray absorption spectra and experimental measurements. © 2018 Wiley Periodicals, Inc.     

Pronounced non-Condon effects in the ultrafast infrared spectroscopy of water

In the context of vibrational spectroscopy in liquids, non-Condon effects refer to the dependence of the vibrational transition dipole moment of a particular molecule on the rotational and translational coordinates of all the molecules in the liquid. For strongly hydrogen-bonded systems, such as liquid water, non-Condon effects are large. That is, the bond dipole derivative of an OH stretch depends strongly on its hydrogen-bonding environment. Previous calculations of nonlinear vibrational spectroscopy in liquids have not included these non-Condon effects. We find that for water, inclusion of these effects is important for an accurate calculation of, for example, homodyned and heterodyned three-pulse echoes. Such echo experiments have been "inverted" to obtain the OH stretch frequency time-correlation function, but by necessity the Condon and other approximations are made in this inversion procedure. Our conclusion is that for water, primarily because of strong non-Condon effects, this inversion may not lead to the correct frequency time-correlation function. Nevertheless, one can still make comparison between theory and experiment by calculating the experimental echo observables themselves.     

Massively parallel spin–orbit configuration interaction

We describe a new approach to incorporating quantum effects into chemical reaction rate theory using quantum trajectories. Our development is based on the entangled trajectory molecular dynamics method for simulating quantum processes using trajectory integration and ensemble averaging. By making dynamical approximations similar to those underlying classical transition state theory, quantum corrections are incorporated analytically into the quantum rate expression. We focus on a simple model of quantum decay in a metastable system and consider the deep tunneling limit where the classical rate vanishes and the process is entirely quantum mechanical. We compare our approximate estimate with the well-known WKB tunneling rate and find qualitative agreement.     

Non-Gaussian statistics of amide I mode frequency fluctuation of N-methylacetamide in methanol solution: linear and nonlinear vibrational spectra

By carrying out molecular dynamics simulations of an N-methylacetamide (NMA) in methanol solution, the amide I mode frequency fluctuation and hydrogen bonding dynamics were theoretically investigated. Combining an extrapolation formula developed from systematic ab initio calculation studies of NMA-(CH3OH)n clusters with a classical molecular dynamics simulation method, we were able to quantitatively describe the solvatochromic vibrational frequency shift induced by the hydrogen-bonding interaction between NMA and solvent methanol. It was found that the fluctuating amide I mode frequency distribution is notably non-Gaussian and it can be decomposed into two Gaussian peaks that are associated with two distinctively different solvation structures. The ensemble-average-calculated linear response function associated with the IR absorption is found to be oscillating, which is in turn related to the doublet amide I band shape. Numerically calculated infrared absorption spectra are directly compared with experiment and the agreement was found to be excellent. By using the Onsager's regression hypothesis, the rate constants of the interconversion process between the two solvation structures were obtained. Then, the nonlinear response functions associated with two-dimensional infrared pump-probe spectroscopy were simulated. The physics behind the two-dimensional line shape and origin of the cross peaks in the time-resolved pump-probe spectra is explained and the result is compared with 2D spectra experimentally measured recently by Woutersen et al.     

Lithium hydroxide phase transition under high pressure: an ab initio molecular dynamics study.

The high-pressure phase transition in the deuterated lithium hydroxide crystalline state has been studied by Car-Parrinello molecular dynamics simulations, in the constant-pressure, constant-temperature ensemble. The recently developed metadynamics approach has been applied to encourage the system to transform into different phases in an affordable simulation time. A previously not completely characterized high-pressure phase has been obtained. The structural and spectroscopic properties have been studied and compared with the neutron scattering, infrared and Raman measurements. It has been found that the calculated structure differs slightly from the experimental hypothesis, and that the presence of strong hydrogen bonds is the source of the red shift and of the characteristic features of the OD-stretching bands in both IR and Raman spectra.