Hoover Thermostat Research Articles

Ad hoc continuum-atomistic thermostat for modeling heat flow in molecular dynamics simulations

An ad hoc thermostating procedure that couples a molecular dynamics (MD) simulation and a numerical solution to the continuum heat flow equation is presented. The method allows experimental thermal transport properties to be modeled without explicitly including electronic degrees of freedom in a MD simulation. The method is demonstrated using two examples, heat flow from a constant temperature silver surface into a single crystal bulk, and a tip sliding along a silver surface. For the former it is shown that frictional forces based on the Hoover thermostat applied locally to grid regions of the simulation are needed for effective feedback between the atomistic and continuum equations. For fast tip sliding the thermostat results in less surface heating, and higher frictional and normal forces compared to the same simulation without the thermostat.

From molecular dynamics to hydrodynamics: A novel Galilean invariant thermostat

This paper proposes a novel thermostat applicable to any particle-based dynamic simulation. Each pair of particles is thermostated either (with probability P) with a pairwise Lowe-Andersen thermostat [C. P. Lowe, Europhys. Lett. 47, 145 (1999)] or (with probability 1-P) with a thermostat that is introduced here, which is based on a pairwise interaction similar to the Nosé-Hoover thermostat. When the pairwise Nosé-Hoover thermostat dominates (low P), the liquid has a high diffusion coefficient and low viscosity, but when the Lowe-Andersen thermostat dominates, the diffusion coefficient is low and viscosity is high. This novel Nosé-Hoover-Lowe-Andersen thermostat is Galilean invariant and preserves both total linear and angular momentum of the system, due to the fact that the thermostatic forces between each pair of the particles are pairwise additive and central. We show by simulation that this thermostat also preserves hydrodynamics. For the (noninteracting) ideal gas at P = 0, the diffusion coefficient diverges and viscosity is zero, while for P > 0 it has a finite value. By adjusting probability P, the Schmidt number can be varied by orders of magnitude. The temperature deviation from the required value is at least an order of magnitude smaller than in dissipative particle dynamics (DPD), while the equilibrium properties of the system are very well reproduced. The thermostat is easy to implement and offers a computational efficiency better than (DPD), with better temperature control and greater flexibility in terms of adjusting the diffusion coefficient and viscosity of the simulated system. Applications of this thermostat include all standard molecular dynamic simulations of dense liquids and solids with any type of force field, as well as hydrodynamic simulation of multiphase systems with largely different bulk viscosities, including surface viscosity, and of dilute gases and plasmas.

Constant temperature molecular dynamics simulations of energetic particle–solid collisions: comparison of temperature control methods

Appropriate temperature control methods must be incorporated into simulations that maintain constant temperature. A specific role played by these methods in modeling energetic particle–solid collisions is to absorb the excess energy wave generated by the collision that, if left unchecked, propagates through the material and reflects from the system boundaries. In this study, five temperature control methods are investigated for use in molecular dynamics simulations of carbon cluster deposition on diamond surfaces. These five methods are the Nosé–Hoover thermostat, the generalized Langevin equation (GLEQ) approach, the Berendsen method, a modified GLEQ approach where an extra damping mechanism is introduced, and a combination of the GLEQ and Berendsen methods. The temperature control capability and the effectiveness of these methods at reducing the amplitude of the reflected energy wave produced in these systems are compared and discussed. It is found that the realistic performance of these methods depends on the specifics of the system, including incident energy and substrate size. Among the five methods considered, the Berendsen method is found to be effective at removing excess energy in the early stages of the deposition process, but the resultant final temperatures are relatively high in most cases. Additionally, the performance of the Nosé–Hoover thermostat is similar to that of the Berendsen method. If the diamond substrate is large enough and the incident energy is not too high, the GLEQ approach is found to be sufficient for removal of excess energy. However, the modified GLEQ approach and the new, combined thermostat perform moderately better than the other approaches at high incident energies.

Grand canonical molecular dynamic simulations for polar and interfacial systems

Purely molecular dynamical formulation of grand canonical ensemble due to Pettitt and co-workers [C.G. Lynch, B.M. Pettitt, J. Chem. Phys. 107 (1997) 8594–8610, and references therein] was extended to implement Nosé–Hoover thermostat and introduce multiple fractional particles. The algorithm was applied to simulation of μVT ensembles of pure TIP4P water and methanol molecules at 298 K, as well as the process of water penetration of a N-acetyl-morpholine (NAM) membrane, where our scheme successfully generated water molecules to compensate for those lost in the NAM mebrane. The procedure reproduced experimental density of water for input chemical potential of −24.0 kJ/mol, as well that of methanol for −19.0 kJ/mol. It was shown that by putting the fractional dynamics on hold and assigning values corresponding to desired integration points to the fractional particles, one can apply the polynomial-path method of thermodynamic integration to configurations sampled during the simulation run. The resulting free energies agreed very well with model-specific values of excess free energy of water, while showing considerable deviation in case of methanol.

Calculation of Free Energy Profiles for Elementary Bimolecular Reactions by ab Initio Molecular Dynamics: Sampling Methods and Thermostat Considerations

A study designed to refine the procedure for performing ab initio molecular dynamics calculations (AIMD) on chemical reactions is presented. Of key interest is the calculation of changes in free energy along the entire reaction path. Several simple elementary reactions are studied with the Car−Parrinello projector augmented-wave (CP−PAW) density-functional theory (DFT) methodology. The illustrative gas-phase bimolecular addition reactions are (i) a σ complexation of BH3 + H2O → H2O·BH3, (ii) the Diels−Alder reactions of butadiene with ethene, C4H6 + C2H4 → cyclohexene, 1,3-cyclopentadiene (CP) and ethene, CP + C2H4→ norbornene, and the stereoselective reaction of 5-amino-CP with ethene, amino-CP + C2H4 → amino-norbornene, (iii) the carbene cyclopropanation Cl2C + C2H4→ Cl2C3H4, and (iv) the dimerization of ketene. These reactions were used to test both the slow-growth and point-wise thermodynamic integration (STI and PTI) methods of phase-space sampling as well as the Nose−Hoover and Andersen thermostats....

Time-reversible deterministic thermostats

Seven different time-reversible deterministic thermostats are considered here and applied to a simple particle-based nonequilibrium heat-flow problem. This approach is robust. Results for all these different thermostats agree rather well for system widths of ten particle diameters or more. The simplest of the thermostats is the Gauss–Nosé–Hoover thermostat, based on kinetic-energy control. Higher moments of the particle momenta can be controlled by extensions of this idea involving as many as three additional thermostat variables. Generalizations of the deterministic thermostats suited to simulating “stochastic” and “Brownian” dynamics are discussed here too.

Large-system phase-space dimensionality loss in stationary heat flows

Thermostated tethered harmonic lattices provide good illustrations of the phase-space dimensionality loss Δ D which occurs in the strange attractor distributions characterizing stationary nonequilibrium flows. We use time-reversible nonequilibrium molecular dynamics, with two Nosé–Hoover thermostats, one hot and one cold, to study a family of square heat-conducting systems. We find a phase-space dimensionality loss which can exceed the dimensionality associated with the two driving Nosé–Hoover thermostats by as much as a factor of four. We also estimate the dimensionality loss ΔD H in the nonthermostated (Hamiltonian) part of phase-space. By measuring the projection of the total dimensionality loss there we show that nearly all of the loss occurs in the nonthermostated part. Thus this loss, which characterizes the extreme rarity of nonequilibrium states, persists in the large-system thermodynamic limit.

Melting of palladium clusters—Canonical and microcanonical Monte Carlo simulation

We present Monte Carlo simulations of single palladium clusters of 13, 34, 54, 55, 147 and 309 atoms. The clusters are modeled by a many-body potential and they have been simulated at constant temperature or constant total energy. The caloric curves of the clusters, with the exception of Pd34, exhibit an S-bend at melting which is typical for a finite system. We have also observed the typical coexistence region of solid and molten clusters both in the canonical and the microcanonical ensembles. Pd34, in contrast, melts without an accompanying peak in heat capacity and at melting the atoms become mobile without any significant change in geometric structure. For the larger clusters a free energy barrier inhibits phase switching. In some cases of phase change from molten to solid structure the barrier is of purely entropic character. By a conversion of the results in the microcanonical simulations into temperature-dependent data, the simulations at fixed temperature and fixed total energy have been compared. The agreement is in most cases good. The results are furthermore compared to earlier molecular dynamics simulations with the Nose–Hoover thermostat. These results are in good agreement with the Monte Carlo simulations as well.

Dynamical and Statistical Mechanical Characterization of Temperature Coupling Algorithms

In this article, we investigate molecular dynamics (MD) trajectories of a butane molecule, as obtained using different types of thermostats. Results show that at low temperature, where the harmonic approximation holds, the Nose'−Hoover (NH) thermostat fails to reproduce the statistical mechanical behavior, even using simulation lengths of millions of time steps, whereas the Gaussian isokinetic (IG) thermostat reproduces quite well the expected statistical mechanical values. The Berendsen's coupling (BC) provides good results for basic properties such as the average potential and kinetic energies but fails in reproducing the canonical fluctuations. Moreover, using the speed of divergence of initially nearby trajectories in phase space as a measure of the dynamical chaoticity, we found that the NH thermostat provides very slow divergence for the physical phase space degrees of freedom, concentrating most of its chaoticity in the dynamics of the thermostat virtual degree of freedom. On the contrary, the IG thermostat provides always highly diverging trajectories in phase space, characterized by a high chaoticity of each degree of freedom. Finally, the BC thermostat provides a moderate chaotic behavior for all of the degrees of freedom. Such results suggest that even assuming for both the “rigorous” algorithms (NH and IG) a full ergodic behavior the NH thermostat could require an extremely long time to achieve convergence of the time averaged properties.

Multicanonical ensemble with Nosé–Hoover molecular dynamics simulation

We demonstrate that molecular dynamics simulations using the force scaling method with a Nosé–Hoover-chain thermostat are capable of generating multicanonical ensembles. The frequency distribution of the Nosé–Hoover-chain is broad enough to handle the energy dependent force scaling factor over a wide potential energy range, when three independent Nosé–Hoover thermostats corresponding to the three orthogonal directions are attached to each particle. The performance of this method has been tested by reproducing various equilibrium properties of one-dimensional model potential, an Ar13 cluster, and a flexible water model.

M.Dyna Mix – a scalable portable parallel MD simulation package for arbitrary molecular mixtures

A general purpose, scalable parallel molecular dynamics package for simulations of arbitrary mixtures of flexible or rigid molecules is presented. It allows use of most types of conventional molecular-mechanical force fields and contains a variety of auxiliary terms for inter- and intramolecular interactions, including an harmonic bond-stretchings. It can handle both isotropic or ordered systems. Besides an NVE MD ensemble, the simulations can also be carried out in either NVT or NPT ensembles, by employing the Nosé–Hoover thermostats and barostats, respectively. If required, the NPT ensemble can be generated by maintaining anisotropic pressures. The simulation cell can be either cubic, rectangular, hexagonal or a truncated octahedron, with corresponding periodic boundary conditions and minimum images. In all cases, the optimized Ewald method can be used to treat the Coulombic interactions. Double time-step or constrained dynamics schemes are included. An external electric field can be applied across the simulation cell. The whole program is highly modular and is written in standard Fortran 77. It can be compiled to run efficiently both on parallel and sequential computers. The inherent complexity of the studied system does not affect the scalability of the program. The scaling is good with the size of the system and with the number of processors. The portability of the program is good, it runs regularly on several common single- and multiprocessor platforms, both scalar and vector architectures included.

Grand canonical molecular dynamics for TIP4P water systems

An algorithm was developed enabling implementation of a Nosé—Hoover thermostat within the framework of grand canonical molecular dynamics [Lynch, C. G. and Pettitt, B. M., 1997, J. chem. Phys., 107, 8594]. The proposed algorithm could readily be extended to mixtures of molecular species with different chemical potentials as shown in the paper. This algorithm was first applied to simulate a μVT ensemble of TIP4P water molecules at 298 K by means of a system comprising a number of full particles and a single scaled (fractional) particle, with the scaling factor considered as a dynamic variable in its own right and chemical potential a pre-set parameter. Our finding showed that the scheme with a single fractional particle tended to freeze in metastable states as well as failed to reproduce either the real-life (−24.05 kJmol−1) or the model-specific chemical potential of water (−23.0kJ mol−1). In order to overcome the inadequacy of a single fractional particle as a chemical potential ‘probe’ the treatment of Pettitt and co-workers was extended to introduce multiple fractional particles. The extended scheme (with 4 fractional particles) was able to reproduce the actual density of water for the driving chemical potential of -24.0k mo−1. The actual behaviour of the density as a function of the chemical potential also agreed quite well with both the results of thermo-dynamic integration and the findings of Pettitt and co-workers.

The effect of an external electric field on the structure of liquid water using molecular dynamics simulations

Using molecular dynamics simulations with the rigid TIP4P water model, we have analyzed the structural change of liquid water induced by an external electric field. The temperature was controlled with a Nosé–Hoover thermostat. In this paper, we report the acquisition of liquid water with the enhanced structural regularity by applying an electric field. From the simulations under various strengths of the electric field, we can see that the threshold for the significant structural change is thought to be between 0.2 and 0.15 V/Å. When the number of six-membered rings is increased by the external electric field, so that water is forced to have structural regularity, we calculate the diffusion coefficients and discover that water we make in the simulations is not solid but still liquid under the electric field.

Computational approaches to the study of some lanthanide (III)-polyazamacrocyclic chelates for magnetic resonance imaging

A set of parameters consistent with the CHARMM force field has been determined for molecular dynamics simulations of several DOTA– and DOTP–Ln(III) chelates. Bonding and van der Waals parameters were derived from the available experimental data and analogy to similar ones in the existing force field. Net atomic charges were derived from ab initio calculations at the Hartree–Fock level to reproduce molecular electrostatic potentials (ESPs), with an effective core potential (ECP) basis set for the metal ion and the 6-31G* basis set for the ligand atoms. The charges are consistent with the TIP3P water model. Preliminary molecular dynamics simulations of the lanthanide chelates in aqueous solution were performed using the Nose–Hoover thermostat at 300 K. The new parameters correctly predicted the molecular structures and stability of the chelates major and minor isomers. ©1999 John Wiley & Sons, Inc. Int J Quant Chem 73: 237–248, 1999

Thermostats: Analysis and application.

Gaussian isokinetic and isoenergetic deterministic thermostats are reviewed in the correct historical context with their later justification using Gauss' principle of least constraint. The Nose-Hoover thermostat for simulating the canonical ensemble is also developed. For some model systems the Lyapunov exponents satisfy the conjugate pairing rule and a Hamiltonian formulation is obtained. We prove the conjugate pairing rule for nonequilibrium systems where the force is derivable from a potential. The generalized symplectic structure and Hamiltonian formulation is discussed. The application of such thermostats to the Lorentz gas is considered in some detail. The periodic orbit expansion methods are used to calculate averages and to categorize the generic transitions in the structure of the attractor. We prove that the conductivity in the nonequilibrium Lorentz gas is non-negative. (c) 1998 American Institute of Physics.

Non-equilibrium molecular dynamics integrators using Maple

The conjugate pairing rule states that for Hamiltonian systems to which a Gaussian or Nosé–Hoover thermostat is applied, the spectrum of Lyapunov exponents should exhibit a simple symmetry: that the sum of conjugate pairs of exponents should be a constant independent of the pair index. This symmetry is important for it enables one to calculate the entropy production and thereby the non-equilibrium transport coefficient, by summing just two exponents, say the maximal exponents. However, the conjugate pairing rule has never been proved for the standard computer simulation algorithm for shear flow namely the Sllod algorithm. We use the symbolic algebra program, Maple, to examine the validity of the conjugate pairing rule for the Sllod algorithm. We do this by using the multiple precision facility of Maple to carry out computer simulations for Sllod to very high accuracy. We also use the symbolic algebra capability of Maple to prove the conjugate pairing rule for ideal gases subject to thermostatted shear flow and also for unthermostatted shear flow of fluids at arbitrary density.

Melting and evaporation of argon clusters

Molecular dynamics simulation with a Nosé–Hoover thermostat was used to study melting and evaporation of free icosahedral argon clusters containing 13 to 1415 atoms. Clusters of 147 atoms or less were found to melt at temperatures clearly below the bulk melting temperature in agreement with previous results. Clusters containing 309 atoms or more were observed to desorb atoms at temperatures where the core of the cluster is solid. As a consequence of this a reliable determination of their melting temperatures using molecular dynamics was found to be complicated.

Molecular dynamics simulation for the cluster formation process of Lennard-Jones particles: Magic numbers and characteristic features

Cluster formation of Lennard-Jones particles (65 536 atoms in a unit cell with an overall number density equal to 0.0149) was simulated by molecular dynamics. The temperature was set to decrease linearly with time by various thermostats, starting from a gas state temperature and ending at zero temperature. With the Nosé–Hoover thermostat, it was found that the translational temperature of the clusters suddenly decreased almost to zero when the cluster formation drastically increased around a reduced temperature (T*) of 0.5, while the internal temperature decreased linearly. Using the Andersen thermostat, which could simulate the aggregation of particles in an inert gas, both the internal and translational temperatures decreased almost linearly with time. When these thermostats were used, cluster–cluster and cluster–atom collisions did not give any magic number peaks in the size distribution up to 250 atoms/cluster at any temperature. Careful tracing of the cluster growth of 13-atom clusters showed no difference in reactivity between icosahedral and nonicosahedral clusters. To simulate cooling in a supersonic jet, a thermostat which controlled only the translational temperature was introduced. After the clusters were formed by cooling the system with this thermostat, their internal temperature stayed at T*≊0.5, while the translational temperature decreased linearly to zero with time as it was controlled. A long-time evaporation from these high-temperature clusters gave peaks at 13 and 19 (and less significantly at 23 and 26) which are magic number sizes corresponding to single, double, triple, and quadruple icosahedra, respectively. The internal temperatures of 13- and 19-atom clusters were higher than those of other size clusters. Higher evaporation energy was observed for the clusters of 13, 19, 23, and 26 atoms than for other size clusters after the long-time evaporation, but only the 13-atom clusters had the higher evaporation energy after cooling by the Andersen thermostat. These results suggest that magic number clusters were formed by evaporation to be trapped at the magic number sizes, and not by either cluster–atom or cluster–cluster collisions. Analyses of the radial distribution functions and the overall shapes of the generated solidlike clusters consisting of many isomers revealed the following characteristic features: The clusters around 13 and 26 atoms were close to being spherical, and the clusters around 19 atoms were oblate. Clusters around 13 atoms had an icosahedron-based structure. The clusters around 55 atoms formed by the Nosé–Hoover and the Andersen thermostats were close to spherical and had an ordered structure. Clusters from 30 to 50 atoms had a disordered structure or a mixture of the different series of structures.

Dynamics of molecules with internal degrees of freedom by multiple time-step methods

In molecular dynamics simulations of flexible molecules, the treatment of motions that take place on different time scales (e.g., bond stretching and bond angle bending vs dihedral angle librations and overall translation and rotation) or of forces that vary differently with distance (e.g., bonding forces vs electrostatic forces) is a major consideration. The most rapidly varying quantity limits the integration time step, while the most slowly varying quantity usually determines the simulation time required. The reversible multiple time-step methods introduced by Tuckerman, Berne, and Martyna [J. Chem. Phys. 97, 1990 (1992)] allow one to use different integration time steps in the same simulation so as to treat the time development of the slow and fast motions most effectively. As a first step toward extending multiple time-step methods to macromolecules, we investigate their utility for determining the dynamics of n-alkanes as neat liquids and in aqueous solution. Because the dynamics is done in Cartesian rather than internal coordinates, the different time steps are implemented in terms of hard and soft forces acting on individual atoms; e.g., hard forces are contributed by bond stretching (and bond angle bending) terms and soft forces by the remaining intramolecular and intermolecular (van der Waals and electrostatic) terms in the additive empirical potential function used to describe these systems. Integrations with and without a Nose–Hoover thermostat are performed. Comparisons are made with ordinary single time step molecular dynamics and dynamics with constraints applied to the bonds. The stability of the integration methods and the dynamic properties are considered. It is shown that stable integrations that yield agreement with standard dynamics can be achieved with a saving of up to a factor of 5 in computer time.

Control of the adiabatic electronic state in ab initio molecular dynamics

The problem of keeping the electronic state at its adiabatic value in the course of ab initio molecular dynamics runs with the Car–Parrinello method is discussed. Attention is focused on the difficulties which arise when chemical bonds are broken and formed in the course of the atomic motion. A solution to the problems is the coupling of independent Nosé–Hoover thermostats to the up and down spin orbitals. The choice of suitable parameters to control the thermostats is discussed. The success of the method is illustrated in studies of the dissociation/recombination of Na2 and the isomerization of the Na3 molecule by a process of pseudorotation.