An approach for generating trajectory-based dynamics which conserves the canonical distribution in the phase space formulation of quantum mechanics. I. Theories

We have reformulated and generalized our recent work [J. Liu and W. H. Miller, J. Chem. Phys. 126, 234110 (2007)] into an approach for generating a family of trajectory-based dynamics methods in the phase space formulation of quantum mechanics. The approach (equilibrium Liouville dynamics) is in the spirit of Liouville's theorem in classical mechanics. The trajectory-based dynamics is able to conserve the quantum canonical distribution for the thermal equilibrium system and approaches classical dynamics in the classical (? → 0), high temperature (β → 0), and harmonic limits. Equilibrium Liouville dynamics provides the framework for the development of novel theoretical∕computational tools for studying quantum dynamical effects in large∕complex molecular systems.     

NVU dynamics. II. Comparing to four other dynamics

In the companion paper [T. S. Ingebrigtsen, S. Toxvaerd, O. J. Heilmann, T. B. Schr?der, and J. C. Dyre, "NVU dynamics. I. Geodesic motion on the constant-potential-energy hypersurface," J. Chem. Phys. (in press)] an algorithm was developed for tracing out a geodesic curve on the constant-potential-energy hypersurface. Here, simulations of NVU dynamics are compared to results for four other dynamics, both deterministic and stochastic. First, NVU dynamics is compared to the standard energy-conserving Newtonian NVE dynamics by simulations of the Kob-Andersen binary Lennard-Jones liquid, its WCA version (i.e., with cut-off's at the pair potential minima), and the Lennard-Jones Gaussian liquid. We find identical results for all quantities probed: radial distribution functions, incoherent intermediate scattering functions, and mean-square displacement as function of time. Arguments are presented for the equivalence of NVU and NVE dynamics in the thermodynamic limit; in particular, to leading order in 1∕N these two dynamics give identical time-autocorrelation functions. In the final part of the paper, NVU dynamics is compared to Monte Carlo dynamics, to a diffusive dynamics of small-step random walks on the constant-potential-energy hypersurface, and to Nose?-Hoover NVT dynamics. If time is scaled for the two stochastic dynamics to make single-particle diffusion constants identical to that of NVE dynamics, the simulations show that all five dynamics are equivalent at low temperatures except at short times.     

Langevin dynamics in constant pressure extended systems

The advantages of performing Langevin dynamics in extended systems are discussed. A simple Langevin dynamics scheme for producing the canonical ensemble is reviewed, and is then extended to the Hoover ensemble. We show that the resulting equations of motion generate the isobaric-isothermal ensemble. The Parrinello-Rahman ensemble is then discussed and we show that despite the presence of intrinsic probability gradients in this system, a Langevin dynamics approach samples the extended phase space in the correct fashion. The implementation of these methods in the ab initio plane wave density functional theory code CASTEP [M. D. Segall, P. L. D. Lindan, M. J. Probert, C. J. Pickard, P. J. Hasnip, S. J. Clarke, and M. C. Payne, J. Phys.: Condens. Matter 14, 2717 (2003)] is demonstrated.     

Rotationally resolved reactive scattering: imaging detailed Cl+C2H6 reaction dynamics

The hydrogen atom abstraction reaction of Cl (2P3/2) with ethane has been studied using the crossed molecular beam technique with dc slice imaging at collision energies from 3.2 to 10.4 kcal/mol. The products HCl (v,J) (v = 0, J = 0-5) were state-selectively detected using 2+1 resonance enhanced multiphoton ionization. The images were used to obtain the center-of-mass frame product angular distributions and translational energy release distributions. Two general features were found in all probed HCl quantum states at 6.7 kcal/mol collision energy, and these features have distinct translational energy release and angular distributions, as described for HCl (v = 0, J = 2) in a recent preliminary report [Li et al., J. Chem. Phys. 124, 011102 (2006)]. The results for HCl (v = 0, J = 2) at four collision energies were also compared to investigate the energy-dependent dynamics. We discuss the reaction in terms of a variety of models of polyatomic reaction dynamics. The dynamics of this well studied system are more complicated than can be accounted for by a single mechanism, and the results call for further theoretical and experimental investigations.     

Coriolis-coupled wave packet dynamics of H + HLi reaction

We investigated the effect of Coriolis coupling (CC) on the initial state-selected dynamics of H+HLi reaction by a time-dependent wave packet (WP) approach. Exact quantum scattering calculations were obtained by a WP propagation method based on the Chebyshev polynomial scheme and ab initio potential energy surface of the reacting system. Partial wave contributions up to the total angular momentum J=30 were found to be necessary for the scattering of HLi in its vibrational and rotational ground state up to a collision energy approximately 0.75 eV. For each J value, the projection quantum number K was varied from 0 to min (J, K(max)), with K(max)=8 until J=20 and K(max)=4 for further higher J values. This is because further higher values of K do not have much effect on the dynamics and also because one wishes to maintain the large computational overhead for each calculation within the affordable limit. The initial state-selected integral reaction cross sections and thermal rate constants were calculated by summing up the contributions from all partial waves. These were compared with our previous results on the title system, obtained within the centrifugal sudden and J-shifting approximations, to demonstrate the impact of CC on the dynamics of this system.     

Quantum dynamics study of H+NH3-->H2+NH2 reaction

We report in this paper a quantum dynamics study for the reaction H+NH3-->NH2+H2 on the potential energy surface of Corchado and Espinosa-Garcia [J. Chem. Phys. 106, 4013 (1997)]. The quantum dynamics calculation employs the semirigid vibrating rotor target model [J. Z. H. Zhang, J. Chem. Phys. 111, 3929 (1999)] and time-dependent wave packet method to propagate the wave function. Initial state-specific reaction probabilities are obtained, and an energy correction scheme is employed to account for zero point energy changes for the neglected degrees of freedom in the dynamics treatment. Tunneling effect is observed in the energy dependency of reaction probability, similar to those found in H+CH4 reaction. The influence of rovibrational excitation on reaction probability and stereodynamical effect are investigated. Reaction rate constants from the initial ground state are calculated and are compared to those from the transition state theory and experimental measurement.     

Low-temperature nucleation in a kinetic Ising model under different stochastic dynamics with local energy barriers

Using both analytical and simulational methods, we study low-temperature nucleation rates in kinetic Ising lattice-gas models that evolve under two different Arrhenius dynamics that interpose between the Ising states a transition state representing a local energy barrier. The two dynamics are the transition-state approximation [T. Ala-Nissila, J. Kjoll, and S. C. Ying, Phys. Rev. B 46, 846 (1992)] and the one-step dynamic [H. C. Kang and W. H. Weinberg, J. Chem. Phys. 90, 2824 (1989)]. Even though they both obey detailed balance and are here applied to a situation that does not conserve the order parameter, we find significant differences between the nucleation rates observed with the two dynamics, and between them and the standard Glauber dynamic [R. J. Glauber, J. Math. Phys. 4, 294 (1963)], which does not contain transition states. Our results show that great care must be exercised when devising kinetic Monte Carlo transition rates for specific physical or chemical systems.     

Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems

Pade? spectrum decomposition is an optimal sum-over-poles expansion scheme of Fermi function and Bose function [J. Hu, R. X. Xu, and Y. J. Yan, J. Chem. Phys. 133, 101106 (2010)]. In this work, we report two additional members to this family, from which the best among all sum-over-poles methods could be chosen for different cases of application. Methods are developed for determining these three Pade? spectrum decomposition expansions at machine precision via simple algorithms. We exemplify the applications of present development with optimal construction of hierarchical equations-of-motion formulations for nonperturbative quantum dissipation and quantum transport dynamics. Numerical demonstrations are given for two systems. One is the transient transport current to an interacting quantum-dots system, together with the involved high-order co-tunneling dynamics. Another is the non-Markovian dynamics of a spin-boson system.     

IBIsCO: a molecular dynamics simulation package for coarse-grained simulation

IBIsCO is a parallel molecular dynamics simulation package developed specially for coarse-grained simulations with numerical potentials derived by the iterative Boltzmann inversion (IBI) method (Reith et al., J Comput Chem 2003, 24, 1624). In addition to common features of molecular dynamics programs, the techniques of dissipative particle dynamics (Groot and Warren, J Chem Phys 1997, 107, 4423) and Lowe-Andersen dynamics (Lowe, Europhys Lett 1999, 47, 145) are implemented, which can be used both as thermostats and as sources of friction to compensate the loss of degrees of freedom by coarse-graining. The reverse nonequilibrium molecular dynamics simulation method (Müller-Plathe, Phys Rev E 1999, 59, 4894) for the calculation of viscosities is also implemented. Details of the algorithms, functionalities, implementation, user interfaces, and file formats are described. The code is parallelized using PE_MPI on PowerPC architecture. The execution time scales satisfactorily with the number of processors.     

Communication: Universal Markovian reduction of Brownian particle dynamics

Non-Markovian processes can often be turned Markovian by enlarging the set of variables. Here we show, by an explicit construction, how this can be done for the dynamics of a Brownian particle obeying the generalized Langevin equation. Given an arbitrary bath spectral density J(0), we introduce an orthogonal transformation of the bath variables into effective modes, leading stepwise to a semi-infinite chain with nearest-neighbor interactions. The transformation is uniquely determined by J(0) and defines a sequence {J(n)}(n∈N) of residual spectral densities describing the interaction of the terminal chain mode, at each step, with the remaining bath. We derive a simple one-term recurrence relation for this sequence and show that its limit is the quasi-Ohmic expression provided by the Rubin model of dissipation. Numerical calculations show that, irrespective of the details of J(0), convergence is fast enough to be useful in practice for an effective Ohmic reduction of the dissipative dynamics.     

Understanding Molecular Dynamics with Stochastic Processes via Real or Virtual Dynamics

Molecular dynamics with the stochastic process provides a convenient way to compute structural and thermodynamic properties of chemical, biological, and materials systems. It is demonstrated that the virtual dynamics case that we proposed for the Langevin equation[J. Chem. Phys. 147 , 184104 (2017)] in principle exists in other types of stochastic thermostats as well. The recommended "middle" scheme[J. Chem. Phys. 147 , 034109 (2017)] of the Andersen thermostat is investigated as an example. As shown by both analytic and numerical results, while the real and virtual dynamics cases approach the same plateau of the characteristic correlation time in the high collision frequency limit, the accuracy and efficiency of sampling are relatively insensitive to the value of the collision frequency in a broad range. After we compare the behaviors of the Andersen thermostat to those of Langevin dynamics, a heuristic schematic representation is proposed for understanding efficient stochastic thermostatting processes with molecular dynamics.     

Reactive scattering dynamics in atom+polyatomic systems: F+C2H6-->HF(v,J)+C2H5

State-to-state scattering dynamics of F+C2H6-->HF(v,J)+C2H5 have been investigated at Ecom=3.2(6) kcalmol under single-collision conditions, via detection of nascent rovibrationally resolved HF(v,J) product states with high-resolution infrared laser absorption methods. State-resolved Doppler absorption profiles are recorded for multiple HF(v,J) transitions originating in the v=0,1,2,3 manifold, analyzed to yield absolute column-integrated densities via known HF transition moments, and converted into nascent probabilities via density-to-flux analysis. The spectral resolution of the probe laser also permits Doppler study of translational energy release into quantum-state-resolved HF fragments, which reveals a remarkable linear correlation between (i) HF(v,J) translational recoil and (ii) the remaining energy available, Eavail=Etot-E(HF(v,J)). The dynamics are interpreted in the context of a simple impulsive model based on conservation of linearangular momentum that yields predictions in good agreement with experiment. Deviations from the model indicate only minor excitation of ethyl vibrations, in contrast with a picture of extensive intramolecular vibrational energy flow but consistent with Franck-Condon excitation of the methylene CH2 bending mode. The results suggest a relatively simple dynamical picture for exothermic atom+polyatomic scattering, i.e., that of early barrier dynamics in atom+diatom systems but modified by impulsive recoil coupling at the transition state between translationalrotational degrees of freedom.     

Vibrational state-dependent predissociation dynamics of ClO (A (2)Pi(3/2)): Insight from correlated fine structure branching ratios

We have studied the v'-dependent predissociation dynamics of the ClO A (2)Pi(3/2) state using velocity-map ion-imaging. Experimental final correlated state branching ratios, i.e. Cl((2)P(J=3/2,1/2)) + O((3)P(J=2,1,0)) channels, have been measured for v' = 6-11. We find that the branching ratios are highly variable and depend strongly on v', providing a window into the v'-dependent predissociation mechanism. A comparison of the experimental results with the recent model of Lane et al. (I. C. Lane, W. H. Howie and A. J. Orr-Ewing, Phys. Chem. Chem. Phys., 1999, 1, 3087) in both the diabatic and adiabatic limits suggests that the dynamics are closer to the diabatic limit. The overall Cl((2)P(J)) branching ratios are in good agreement with the diabatic model results. There are significant differences, however, between theory and experiment at the correlated state level, demonstrating the sensitivity of correlated measurements to the role of the exit channel coupling in the predissociation dynamics. The results highlight the need for more sophisticated quantum dynamical calculations to describe the correlated fine structure branching ratios in this system.     

Describing the component dynamics in miscible polymer blends: towards a fully predictive model

We have recently proposed [D. Cangialosi et al., J. Chem. Phys. 123, 144908 (2005)] an extension of the Adam-Gibbs [J. Chem. Phys. 43, 139 (1965)] theory, combined with the concept of self-concentration, to describe the temperature dependence of the relaxation time for the component segmental dynamics in miscible polymer blends. Thus, we were able to obtain the dynamics of each component in the blend starting from the knowledge of the dynamic and thermodynamic data of the pure polymers, with a single fitting parameter (alpha) which had to be obtained from the fitting of the experimental data. In the present work we demonstrate that this model is also suitable to describe the polymer segmental dynamics in concentrated polymer solutions. From this result we have developed a new route for determining the value of the alpha parameter associated with any given polymer. Once this value is known for the two components of a possible polymer blend, our model for polymer blends dynamics becomes fully predictive.     

Stereospecific NMR assignments of prochiral methyls, rotameric states and dynamics of valine residues in malate synthase G

Near complete stereospecific assignments of the prochiral methyl carbons of Leu and Val residues in malate synthase G, a 723 residue enzyme, are reported. Assignments were obtained on the basis of a 10% fractional (13)C-labeling strategy developed by Wüthrich and co-workers [Neri, D; Szyperski, T; Otting, G; Senn, H; Wüthrich, K. Biochemistry 1989, 28, 7510-7516] and, in the case of Val residues, supplemented with results from a series of new methyl-TROSY quantitative J experiments for measuring (3)J(C)(gamma)(N) and (3)J(C)(gamma)(C)' couplings. The measured (3)J couplings were also used to probe Val side chain dynamics. A strong correlation is observed between rotamer averaging established on the basis of the couplings and side chain millisecond time scale dynamics measured using methyl-TROSY based (1)H-(13)C multiple quantum relaxation dispersion experiments.     

Polarizable and flexible model for ethanol

A polarizable, flexible model for ethanol is obtained based on an extensive series of B3LYP/6-311++G(d,p) calculations and molecular dynamics simulations. The ethanol model includes electric-field dependence in both the atomic charges and the intramolecular degrees of freedom. Field-dependent intramolecular potentials have been attempted only once previously, for OH and HH stretches in water [P. Cicu et al., J. Chem. Phys. 112, 8267 (2000)]. The torsional potential involving the hydrogen-bonding hydrogen in ethanol is found to be particularly field sensitive. The methodology for developing field-dependent potentials can be readily generalized to other molecules and is discussed in detail. Molecular dynamics simulations of bulk ethanol are performed and the results are assessed based on comparisons with the self-diffusion coefficient [N. Karger et al., J. Chem. Phys. 93, 3437 (1990)], dielectric constant [J. T. Kindt and C. A. Schmuttenmaer, J. Phys. Chem. 100, 10373 (1996)], enthalpy of vaporization [R. C. Wilhoit and B. J. Zwolinski, J. Phys. Chem. Ref. Data, Suppl. 2, 2 (1973)], and experimental interatomic distributions [C. J. Benmore and Y. L. Loh, J. Chem. Phys. 112, 5877 (2000)]. The simultaneous variation of the atomic charges and the intramolecular potentials requires modified equations of motion and a multiple time step algorithm has been implemented to solve these equations. The article concludes with a discussion of the bulk structure and properties with an emphasis on the hydrogen bonding network.     

Recovering the Crooks equation for dynamical systems in the isothermal-isobaric ensemble: a strategy based on the equations of motion

The Crooks equation [Eq. (10) in J. Stat. Phys. 90, 1481 (1998)] relates the work done on a system during a nonequilibrium transformation to the free energy difference between the final and the initial state of the transformation. Recently, the authors have derived the Crooks equation for systems in the canonical ensemble thermostatted by the Nose-Hoover or Nose-Hoover chain method [P. Procacci et al., J. Chem. Phys. 125, 164101 (2006)]. That proof is essentially based on the fluctuation theorem by Evans and Searles [Adv. Phys. 51, 1529 (2002)] and on the equations of motion. Following an analogous approach, the authors derive here the Crooks equation in the context of molecular dynamics simulations of systems in the isothermal-isobaric (NPT) ensemble, whose dynamics is regulated by the Martyna-Tobias-Klein algorithm [J. Chem. Phys. 101, 4177 (1994)]. Their present derivation of the Crooks equation correlates to the demonstration of the Jarzynski identity for NPT systems recently proposed by Cuendet [J. Chem. Phys. 125, 144109 (2006)].     

Nature of the water molecules in the palisade layer of a triton X-100 micelle in the presence of added salts: A solvation dynamics study

The effect of added electrolytes on the nature of water molecules in the palisade layer of a Triton X-100 (TX-100) micelle has been investigated using solvation dynamics studies of C153 dye in the presence of different concentrations of NaCl, KCl, and CsCl salts. In all of the cases, the solvation dynamics is found to be biexponential in nature. It is seen that in the presence of added salts the solvation dynamics becomes slower. As previously reported (Charlton et al. J. Phys. Chem. B 2000, 104, 8327; Molina-Bolivar et al. J. Phys. Chem. B 2002, 106, 870), the presence of salt increases micellar hydration (and also size) for TX-100, mainly due to enhancement in the mechanically trapped water content in the palisade layer. Under normal circumstances, increased micellar hydration was expected to cause faster solvation dynamics (Kumbhakar et al. J. Phys. Chem. B 2004, 108, 19246), though in the present work, a reverse trend is in fact observed with the added salts. In accordance with solvation dynamics results, fluorescence anisotropy studies also indicate an increase in microviscosity for the palisade layer of the TX-100 micelle with the added salts. The present results have been rationalized assuming that the ions reside in the palisade layer, and due to the hydration of the ions, especially the cations, the water molecules in the palisade layer undergo a kind of clustering, causing the microviscosity to in fact increase rather than decrease as expected due to increased micellar hydration. A partial collapse of the surfactant chains due to their dehydration as caused by the hydration of the ions in the palisade layer may also add to the increase in microviscosity and the consequent retardation in relaxation dynamics in the presence of salts.     

Brownian dynamics simulations of a flexible polymer chain which includes continuous resistance and multibody hydrodynamic interactions

Using methods adapted from the simulation of suspension dynamics, we have developed a Brownian dynamics algorithm with multibody hydrodynamic interactions for simulating the dynamics of polymer molecules. The polymer molecule is modeled as a chain composed of a series of inextensible, rigid rods with constraints at each joint to ensure continuity of the chain. The linear and rotational velocities of each segment of the polymer chain are described by the slender-body theory of Batchelor [J. Fluid Mech. 44, 419 (1970)]. To include hydrodynamic interactions between the segments of the chain, the line distribution of forces on each segment is approximated by making a Legendre polynomial expansion of the disturbance velocity on the segment, where the first two terms of the expansion are retained in the calculation. Thus, the resulting linear force distribution is specified by a center of mass force, couple, and stresslet on each segment. This method for calculating the hydrodynamic interactions has been successfully used to simulate the dynamics of noncolloidal suspensions of rigid fibers [O. G. Harlen, R. R. Sundararajakumar, and D. L. Koch, J. Fluid Mech. 388, 355 (1999); J. E. Butler and E. S. G. Shaqfeh, J. Fluid Mech. 468, 204 (2002)]. The longest relaxation time and center of mass diffusivity are among the quantities calculated with the simulation technique. Comparisons are made for different levels of approximation of the hydrodynamic interactions, including multibody interactions, two-body interactions, and the "freely draining" case with no interactions. For the short polymer chains studied in this paper, the results indicate a difference in the apparent scaling of diffusivity with polymer length for the multibody versus two-body level of approximation for the hydrodynamic interactions.     

Modeling real dynamics in the coarse-grained representation of condensed phase systems

This work presents a systematic multiscale methodology to provide a more faithful representation of real dynamics in coarse-grained molecular simulation models. The theoretical formalism is based on the recently developed multiscale coarse-graining (MS-CG) method [S. Izvekov and G. A. Voth, J. Phys. Chem. B. 109, 2469 (2005); J. Chem. Phys. 123, 134105 (2005)] and relies on the generalized Langevin equation approach and its simpler Langevin equation limit. The friction coefficients are determined in multiscale fashion from the underlying all-atom molecular dynamics simulations using force-velocity and velocity-velocity correlation functions for the coarse-grained sites. The diffusion properties in the resulting CG Brownian dynamics simulations are shown to be quite accurate. The time dependence of the velocity autocorrelation function is also well-reproduced relative to the all-atom model if sufficient resolution of the CG sites is implemented.