Implementations of Nosé-Hoover and Nosé-Poincaré thermostats in mesoscopic dynamic simulations with the united-residue model of a polypeptide chain

Molecular dynamics (MD) simulations generate a canonical ensemble only when integration of the equations of motion is coupled to a thermostat. Three extended phase space thermostats, one version of Nose-Hoover and two versions of Nose-Poincare, are compared with each other and with the Berendsen thermostat and Langevin stochastic dynamics. Implementation of extended phase space thermostats was first tested on a model Lennard-Jones fluid system; subsequently, they were implemented with our physics-based protein united-residue (UNRES) force field MD. The thermostats were also implemented and tested for the multiple-time-step reversible reference system propagator (RESPA). The velocity and temperature distributions were analyzed to confirm that the proper canonical distribution is generated by each simulation. The value of the artificial mass constant, Q, of the thermostat has a large influence on the distribution of the temperatures sampled during UNRES simulations (the velocity distributions were affected only slightly). The numerical stabilities of all three algorithms were compared with each other and with that of microcanonical MD. Both Nose-Poincare thermostats, which are symplectic, were not very stable for both the Lennard-Jones fluid and UNRES MD simulations started from nonequilibrated structures which implies major changes of the potential energy throughout a trajectory. Even though the Nose-Hoover thermostat does not have a canonical symplectic structure, it is the most stable algorithm for UNRES MD simulations. For UNRES with RESPA, the "extended system inside-reference system propagator algorithm" of the RESPA implementation of the Nose-Hoover thermostat was the only stable algorithm, and enabled us to increase the integration time step.     

Periodic boundary condition induced breakdown of the equipartition principle and other kinetic effects of finite sample size in classical hard-sphere molecular dynamics simulation

We examine consequences of the non-Boltzmann nature of probability distributions for one-particle kinetic energy, momentum, and velocity for finite systems of classical hard spheres with constant total energy and nonidentical masses. By comparing two cases, reflecting walls (NVE or microcanonical ensemble) and periodic boundaries (NVEPG or molecular dynamics ensemble), we describe three consequences of the center-of-mass constraint in periodic boundary conditions: the equipartition theorem no longer holds for unequal masses, the ratio of the average relative velocity to the average velocity is increased by a factor of [N/(N-1)]1/2, and the ratio of average collision energy to average kinetic energy is increased by a factor of N/(N-1). Simulations in one, two, and three dimensions confirm the analytic results for arbitrary dimension.     

Enhanced sampling via strong coupling to a heat bath: relationship between Tsallis and multicanonical algorithms

We show that Tsallis and multicanonical statistical mechanics are equivalent under specific conditions and that they describe a system strongly coupled to a heat bath. The concept of the strong coupling to a heat bath, in which energy fluctuation is larger than that in the canonical ensemble [J. Chem. Phys. 119, 7075 (2003)], plays a key role in relating Tsallis formalism to multicanonical formalism. The equivalence between these formalisms allows us to obtain an appropriate q parameter in the Tsallis algorithm to enhance the sampling in the phase space in a manner similar to the multicanonical algorithm. An enhanced sampling in the configurational space by use of the strong coupling formalism is demonstrated in a Lennard-Jones fluid.     

Configurational constant pressure molecular dynamics

We propose two new algorithms for generating isothermal-isobaric molecular dynamics. The algorithms are based on an extended phase space dynamics where two extra degrees of freedom, representing the thermostat and the barostat, are included. These new methods adopt a totally different approach towards molecular dynamics simulation in the isothermal-isobaric ensemble. They are fully configurational in the sense that only the particle positions are required in the control of the system temperature and pressure. Following on from the works of Delhommelle and Evans [Mol. Phys., 99, 1825 (2001)] and of Braga and Travis [J. Chem. Phys., 123, 134101 (2005)] concerning configurational canonical dynamics, these new algorithms can be seen as a natural extension to the isothermal-isobaric ensemble. We have validated both of our new configurational isothermal-isobaric schemes by conducting molecular dynamics simulations of a Lennard-Jones fluid and comparing the static and dynamic properties for a single state point. We find that both schemes generate similar results compared with schemes which use kinetic temperature and pressure control. We have also monitored the response of the system to a series of isothermal compressions and isobaric quenches. We find that the configurational schemes performed at least as well as the kinetic based scheme in bringing the system temperature and pressure into line with the set point values of these variables. These new methods will potentially play a significant role in simulations where the calculation of the kinetic temperature and pressure can be problematic. A well known example resides in the field of nonequilibrium simulations where the kinetic temperature and pressure require a knowledge of the streaming velocity of the fluid in order to calculate the true peculiar velocities (or momenta) that enter into their definitions. These are completely avoided by using our configurational thermostats and barostats, since these are independent of momenta. By extending the analysis of Kusnezov et al. [Ann. Phys., 204, 155 (1990)] in order to derive a set of generalized Nose-Hoover equations of motion which can generate isothermal-isobaric dynamics in a number of different ways, we are able to show that both of our new configurational barostats and Hoover's kinetic isothermal-isobaric scheme are special cases of this more general set of equations. This generalization can be very powerful in generating constant pressure dynamics for a variety of systems.     

Role of angular momentum conservation in unimolecular translational energy release: validity of the orbiting transition state theory

The translational kinetic energy release distribution (KERD) for the halogen loss reaction of the bromobenzene and iodobenzene cations has been reinvestigated on the microsecond time scale. Two necessary conditions of validity of the orbiting transition state theory (OTST) for the calculation of kinetic energy release distributions (KERDs) have been formulated. One of them examines the central ion-induced dipole potential approximation. As a second criterion, an adiabatic parameter is derived. The lower the released translational energy and the total angular momentum, the larger the reduced mass, the rotational constant of the molecular fragment, and the polarizability of the released atom, the more valid is the OTST. Only the low-energy dissociation of the iodobenzene ion (E approximately 0.45 eV, where E is the internal energy above the reaction threshold) is found to fulfill the criteria of validity of the OTST. The constraints that act on the dissociation dynamics have been studied by the maximum entropy method. Calculations of entropy deficiencies (which measure the deviation from a microcanonical distribution) show that the pair of fragments does not sample the whole of the phase space that is compatible with the mere specification of the internal energy. The major constraint that results from conservation of angular momentum is related to a reduction of the dimensionality of the dynamics of the translational motion to a two-dimensional space. A second and minor constraint that affects the KERD leads to a suppression of small translational releases, i.e., accounts for threshold behavior. At high internal energies, the effects of curvature of the reaction path and of angular momentum conservation are intricately intermeddled and it is not possible to specify the share of each effect.     

Multiple time scale molecular dynamics for fluids with orientational degrees of freedom. II. Canonical and isokinetic ensembles

We have developed several multiple time stepping techniques to overcome the limitations on efficiency of molecular dynamics simulations of complex fluids. They include the modified canonical and isokinetic schemes, as well as the extended isokinetic Nosé-Hoover chain approach. The latter generalizes the method of Minary, Tuckerman, and Martyna for translational motion [Phys. Rev. Lett. 93, 150201 (2004)] to systems with both translational and orientational degrees of freedom. Although the microcanonical integrators are restricted to relatively small outer time steps of order of 16 fs, we show on the basis of molecular dynamics simulations of ambient water that in the canonical and isokinetic thermostats the size of these steps can be increased to 50 and 75 fs, respectively (at the same inner time step of 4 fs). Within the generalized isokinetic Nosé-Hoover chain algorithm we have derived, huge outer time steps of order of 500 fs can be used without losing numerical stability and affecting equilibrium properties.     

The flying ice cube: Velocity rescaling in molecular dynamics leads to violation of energy equipartition

This article describes an unexpected phenomenon encountered during MD simulations: velocity rescaling using standard protocols can systematically change the proportion of total kinetic energy (KE) found in motions associated with the various degrees of freedom. Under these conditions, the simulation violates the principle of equipartition of energy, which requires a mean kinetic energy of RT/2 in each degree of freedom. A particularly pathological form of this problem occurs if one does not periodically remove the net translation of (and rotation about) the center of mass. In this case, almost all of the kinetic energy is converted into these two kinds of motion, producing a system with almost no kinetic energy associated with the internal degrees of freedom. We call this phenomenon “the flying ice cube.” We present a mathematical analysis of a simple diatomic system with two degrees of freedom, to document the origin of the problem. We then present examples from three kinds of MD simulations, one being an in vacuo simulation on a diatomic system, one involving a low resolution model of DNA in vacuo, and the third using a traditional all-atom DNA model with full solvation, periodic boundary conditions, and the particle mesh Ewald method for treating long-range electrostatics. Finally, we discuss methods for avoiding the problem. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 726–740, 1998     

A microscopic model of chemical reactions with reorganization. High-energy approximation

A dynamical model of a chemical reaction, accompanied by reorganization of the immediate environment of the isolated chemical subsystem, is proposed. The model enables studying the emergence of nonequilibrium distribution functions as a combined result of the interaction within the dynamical subsystem and the energy exchange with a subsystem of inactive degrees of freedom (thermal bath). The study is based on the quasiclassical high-energy approximation for nonadiabatic effects in the energy exchange within the dynamical subsystem, for strong and weak coupling of the oscillator mode with the thermal bath. Such an approximation allows for the important statement that nonequilibrium effects in thermal reactions are absent if the initial translational distribution along the reaction coordinate and the initial vibrational distribution in transversal degrees of freedom are Boltzmann-like with the same temperature. The results obtained in the absence of the initial equilibrium distribution have been used for interpreting the kinetics of endothermic plasmochemical reactions proceeding under nonequilibrium conditions.     

Generalized simulated tempering for exploring strong phase transitions

An extension of the simulation tempering algorithm is proposed. It is shown to be particularly suited to the exploration of first-order phase transition systems characterized by the backbending or S-loop in the statistical temperature or a microcanonical caloric curve. A guided Markov process in an auxiliary parameter space systematically combines a set of parametrized Tsallis-weight ensemble simulations, which are targeted to transform unstable or metastable energy states of canonical ensembles into stable ones and smoothly join ordered and disordered phases across phase transition regions via a succession of unimodal energy distributions. The inverse mapping between the sampling weight and the effective temperature enables an optimal selection of relevant Tsallis-weight parameters. A semianalytic expression for the biasing weight in parameter space is adaptively updated "on the fly" during the simulation to achieve rapid convergence. Accelerated tunneling transitions with a comprehensive sampling for phase-coexistent states are explicitly demonstrated in systems subject to strong hysteresis including Potts and Ising spin models and a 147 atom Lennard-Jones cluster.     

The role of long-range forces in the determination of translational kinetic energy release. Loss of C4H4+ from the benzene and pyridine cations

Kinetic energy release distributions (KERDs) for the benzene ion fragmenting into C 4H 4 (+) and C 2H 2 have been recorded by double-focusing mass spectrometry in the metastable energy window and by a retarding field experiment up to an energy of 5 eV above the fragmentation threshold. They are compared with those resulting from the HCN loss reaction from the pyridine ion. Both reactions display a similar variation of the kinetic energy release as a function of the internal energy: the average release is smaller than statistically expected, with a further restriction of the phase-space sampling for the C 5H 5N (+) dissociation. Ab initio calculations of the potential-energy profile have been carried out. They reveal a complicated reaction mechanism, the last step of which consists in the dissociation of a weakly bound ion-quadrupole or ion-dipole complex. The KERDs have been analyzed by the maximum entropy method. The fraction of phase space effectively sampled by the pair of fragments has been determined and is similar for both dissociations. Both reactions are constrained by the square root of the released translational energy, epsilon (1/2). This indicates that in the latter stage of the dissociation process, the reaction coordinate is adiabatically decoupled from the bath of the bound degrees of freedom. For the C 6H 6 (+) fragmentation, the analysis of the experimental results strongly suggests that, just as for a spherically symmetric interaction potential, the translational motion is confined to a plane. For the dissociation of the pyridine ion, the main dynamical constraint is also a restriction to a two-dimensional subspace. This dimensionality reduction of the translational phase space is due to the fact that the Hamiltonian of both weakly bound complexes contains a cyclic coordinate.     

Ontology of temperature in nonequilibrium systems

The laws of thermodynamics provide a clear concept of the temperature for an equilibrium system in the continuum limit. Meanwhile, the equipartition theorem allows one to make a connection between the ensemble average of the kinetic energy and the uniform temperature. When a system or its environment is far from equilibrium, however, such an association does not necessarily apply. In small systems, the regression hypothesis may not even apply. Herein, it is shown that in small nonequilibrium systems, the regression hypothesis still holds, though with a generalized definition of the temperature. The latter must now be defined for each such manifestation.     

Atomistic simulations of the solid-liquid transition of 1-ethyl-3-methyl imidazolium bromide ionic liquid

Achieving melting point around room temperature is important for applications of ionic liquids. In this work, molecular dynamics simulations are carried out to investigate the solid-liquid transition of ionic liquid 1-ethyl-3-methyl imidazolium bromide ([emim]Br) by direct heating, hysteresis, void-nucleation, sandwich, and microcanonical ensemble approaches. Variations of the non-bonded energy, density, diffusion coefficient, and translational order parameter of [emim]Br are analyzed as a function of temperature, and a coexisting solid-liquid system is achieved in the microcanonical ensemble method. The melting points obtained from the first three methods are 547 ± 8 K, 429 ± 8 K, and 370 ± 6 K; while for the sandwich method, the melting points are 403 ± 4 K when merging along the x-axis by anisotropic isothermal-isobaric (NPT) ensemble, 393 ± 4 K when along the y-axis by anisotropic NPT ensemble, and 375 ± 4 K when along the y-axis by isotropic NPT ensemble. For microcanonical ensemble method, when the slabs are merging along different directions (x-axis, y-axis, and z-axis), the melting points are 364 ± 3 K, 365 ± 3 K, and 367 ± 3 K, respectively, the melting points we get by different methods are approximately 55.4%, 21.9%, 5.1%, 14.5%, 11.6%, 6.5%, 3.4%, 3.7%, and 4.3% higher than the experimental value of 352 K. The advantages and disadvantages of each method are discussed. The void-nucleation and microcanonical ensemble methods are most favorable for predicting the solid-liquid transition.     

Molecular dynamics simulation in the grand canonical ensemble

An extended system Hamiltonian is proposed to perform molecular dynamics (MD) simulation in the grand canonical ensemble. The Hamiltonian is similar to the one proposed by Lynch and Pettitt (Lynch and Pettitt, J Chem Phys 1997, 107, 8594), which consists of the kinetic and potential energies for real and fractional particles as well as the kinetic and potential energy terms for material and heat reservoirs interacting with the system. We perform a nonlinear scaling of the potential energy parameters of the fractional particle, as well as its mass to vary the number of particles dynamically. On the basis of the equations of motion derived from this Hamiltonian, an algorithm has been proposed for MD simulation at constant chemical potential. The algorithm has been tested for the ideal gas, for the Lennard-Jones fluid over a wide range of temperatures and densities, and for water. The results for the low-density Lennard-Jones fluid are compared with the predictions from a truncated virial equation of state. In the case of the dense Lennard-Jones fluid and water our predicted results are compared with the results reported using other available methods for the calculation of the chemical potential. The method is also applied to the case of vapor-liquid coexistence point predictions.     

A simulation method for the calculation of chemical potentials in small, inhomogeneous, and dense systems

We present a modification of the gauge cell Monte Carlo simulation method [A. V. Neimark and A. Vishnyakov, Phys. Rev. E 62, 4611 (2000)] designed for chemical potential calculations in small confined inhomogeneous systems. To measure the chemical potential, the system under study is set in chemical equilibrium with the gauge cell, which represents a finite volume reservoir of ideal particles. The system and the gauge cell are immersed into the thermal bath of a given temperature. The size of the gauge cell controls the level of density fluctuations in the system. The chemical potential is rigorously calculated from the equilibrium distribution of particles between the system cell and the gauge cell and does not depend on the gauge cell size. This scheme, which we call a mesoscopic canonical ensemble, bridges the gap between the canonical and the grand canonical ensembles, which are known to be inconsistent for small systems. The ideal gas gauge cell method is illustrated with Monte Carlo simulations of Lennard-Jones fluid confined to spherical pores of different sizes. Special attention is paid to the case of extreme confinement of several molecular diameters in cross section where the inconsistency between the canonical ensemble and the grand canonical ensemble is most pronounced. For sufficiently large systems, the chemical potential can be reliably determined from the mean density in the gauge cell as it was implied in the original gauge cell method. The method is applied to study the transition from supercritical adsorption to subcritical capillary condensation, which is observed in nanoporous materials as the pore size increases.     

Dynamics of quantum dissipation systems interacting with fermion and boson grand canonical bath ensembles: hierarchical equations of motion approach

A hierarchical equations of motion formalism for a quantum dissipation system in a grand canonical bath ensemble surrounding is constructed on the basis of the calculus-on-path-integral algorithm, together with the parametrization of arbitrary non-Markovian bath that satisfies fluctuation-dissipation theorem. The influence functionals for both the fermion or boson bath interaction are found to be of the same path integral expression as the canonical bath, assuming they all satisfy the Gaussian statistics. However, the equation of motion formalism is different due to the fluctuation-dissipation theories that are distinct and used explicitly. The implications of the present work to quantum transport through molecular wires and electron transfer in complex molecular systems are discussed.     

Microcanonical analysis of the kinetic method. the meaning of the "apparent entropy"

A rigorous analysis of the kinetic method is carried out using Rice-Ramsperger-Kassel-Marcus (RRKM) theory of microcanonical statistical unimolecular dissociation rates. The model employs a kinetics treatment appropriate for metastable ion dissociation. Proton-bound alkoxide dimer anions are used as model systems, with realistic vibrational and rotational parameters calculated by ab initio methods for the cluster ion and transition states leading to the competitive dissociation channels. The numerical simulations show that the kinetic method plots of ln(I2/I1) versus AAH are nearly linear but can exhibit significant curvature. The apparent entropy obtained in the extended kinetic method is not approximately equal to the thermodynamic entropy difference for dissociation, AAS(T), or for activation, deltadeltaS++(T), either at the effective temperature or at any fixed equilibrium temperature. Instead, the apparent entropy term can be related to the ratio of the microcanonical sum of states of the dissociation transition states for the kinetically selected internal energy of the dissociating ions.     

Comparisons of classical chemical dynamics simulations of the unimolecular decomposition of classical and quantum microcanonical ensembles

Previous studies have shown that classical trajectory simulations often give accurate results for short-time intramolecular and unimolecular dynamics, particularly for initial non-random energy distributions. To obtain such agreement between experiment and simulation, the appropriate distributions must be sampled to choose initial coordinates and momenta for the ensemble of trajectories. If a molecule's classical phase space is sampled randomly, its initial decomposition will give the classical anharmonic microcanonical (RRKM) unimolecular rate constant for its decomposition. For the work presented here, classical trajectory simulations of the unimolecular decomposition of quantum and classical microcanonical ensembles, at the same fixed total energy, are compared. In contrast to the classical microcanonical ensemble, the quantum microcanonical ensemble does not sample the phase space randomly. The simulations were performed for CH(4), C(2)H(5), and Cl(-)---CH(3)Br using both analytic potential energy surfaces and direct dynamics methods. Previous studies identified intrinsic RRKM dynamics for CH(4) and C(2)H(5), but intrinsic non-RRKM dynamics for Cl(-)---CH(3)Br. Rate constants calculated from trajectories obtained by the time propagation of the classical and quantum microcanonical ensembles are compared with the corresponding harmonic RRKM estimates to obtain anharmonic corrections to the RRKM rate constants. The relevance and accuracy of the classical trajectory simulation of the quantum microcanonical ensemble, for obtaining the quantum anharmonic RRKM rate constant, is discussed.