Chapter 8. TVN Dynamics

8.1. Introduction

Questions exist as to whether constant-temperature, constant-volume (TVN) dynamics properly reproduces physical properties of systems. Although the equations of motion guarantee that a proper canonical distribution will be generated eventually, the primary concern is whether the dynamics reaches ergodicity, and thus an adequate sampling of the distribution, in a reasonable amount of time.

The Nosé-Hoover formulation has one free parameter, \tau_s. This variable controls the rate of transfer of kinetic energy from the system to the bath and vice versa and thus the rate of equilibration. There have been no clear guidelines as to how this parameter should be chosen [1].

We investigated how \tau_s affects dynamics in a wide variety of systems, including argon and methane clusters and periodic poly(ethylene) models.

8.2. Procedure

Five structures were generated for testing. Argon atoms were arranged in a Mackay icosahedron [2] of five shells (561 atoms). An additional argon system was built using a 256 atom face centered cubic (fcc) unit cell, 21.6204 Å on a side. A methane cluster was built by replacing each argon atom in a three-shell Mackay icosahedral structure with a methane molecule; the resulting cluster of 735 atoms was then minimized using BIOGRAF [3]. Two poly(ethylene) (PE) systems were also built. One was an infinite system composed of a minimized single chain fragment of 398 atoms in an 18 Å cubic unit cell. The ends of the chain fragment were connected through the unit cell boundary to form the infinite chain. The other PE system was a 2x6x4 supercell of 576 atoms (16 chain fragments, each of 6 monomer units).

The argon systems are a very simple, finite cluster with only heavy atoms and weak forces and an extension of that system to periodic boundary conditions. The methane system is similar to the argon cluster, but it also includes high-frequency C-H bonds. The PE cases have both periodic boundary conditions and high-frequency bonds, first using only one chain and then with inter-chain interactions.

Each system was simulated for 100 ps, using 1 fs timesteps. \tau_s values between 0.01 ps and 1.0 ps were used. The bath temperature was set to 20 K for the argon and methane cases; it was fixed at 300 K for the PE systems.

The kinetic energy for each run was plotted against time to observe how quickly the system converged to the desired temperature. In addition, the kinetic energy distribution was computed by dividing the range of kinetic energies into 100 bins and counting the number of values falling into each bin. In the ideal case, these distributions should be of Gaussian form [1].

8.3. Results

The kinetic energy distributions for various values of \tau_s are plotted in Figures 8-1 through 8-5. The argon cluster results were essentially the same as those for the fcc argon system and are not shown; the single-chain PE results similarly were essentially identical to the multi-chain PE results and are omitted.

Two particular features are notable. At longer values of \tau_s, the distributions have two distinct peaks and are often skewed so that the lower-energy peak is sharper and taller. For the argon systems, the distributions for very short values of \tau_s in Figure 8-3 also show double-peaked behavior.

The first feature is explained by looking at sample plots of the kinetic energy versus time in 8-12. The system is coupled to the bath via what is essentially a harmonic oscillator of period 2\pi\tau_s. For long values of \tau_s (Figures 8-7 and 8-9), relatively few cycles of this oscillator have occurred by the end of the 100 ps simulation. The oscillator frequency is outside the range of the modes of the system, so transfer of energy from the system to the bath occurs only through the anharmonicities of the system's oscillators. This process is slow. Thus, by the end of the simulation the system is still equilibrating and has not yet reached ergodicity. The resulting distribution is dominated by the turning points of the Nosé oscillator, where the system spends the most time.

At shorter values of \tau_s (Figures 8-6, 8-8, and 8-10), the oscillator frequency couples well with the modes of the system. Energy is transferred easily between the system and the bath, and ergodicity is reached more quickly, both in terms of simulation time and in terms of the number of oscillator periods.

At very short values of \tau_s, as seen in the argon cases (Figure 8-3), the distributions again become double-peaked. Note, however, that the range of kinetic energies in these cases is much smaller than for long \tau_s. Again, reference to the kinetic energy versus time plots (Figure 8-9) shows that the system is now highly constrained to a very limited range of KE values. In these cases, with heavy atoms, energy is transferred too rapidly from the system to the bath or vice versa, and the system has no opportunity to reach ergodicity.

The short \tau_s cases for argon were run with an unusually small timestep of 1 fs. Using a more typical longer timestep caused numerical instabilities as the Nosé oscillator failed to be properly integrated.

[IMAGE]
Figure 8-1. KE distribution for multi-chain PE crystal for varying \tau_s .

[IMAGE]
Figure 8-2. KE distribution for methane cluster system for varying \tau_s .

[IMAGE]
Figure 8-3. KE distribution for fcc argon with small \tau_s .

[IMAGE]
Figure 8-4. KE distribution for fcc argon with intermediate \tau_s .

[IMAGE]
Figure 8-5. KE distribution for fcc argon with large \tau_s .

[IMAGE]
Figure 8-6. KE vs. time for multi-chain PE crystal with small \tau_s .

[IMAGE]
Figure 8-7. KE vs. time for multi-chain PE crystal with large \tau_s .

[IMAGE]
Figure 8-8. KE vs. time for methane cluster with small \tau_s .

[IMAGE]
Figure 8-9. KE vs. time for methane cluster with large \tau_s .

[IMAGE]
Figure 8-10. KE vs. time for fcc argon with small \tau_s .

[IMAGE]
Figure 8-11. KE vs. time for fcc argon with intermediate \tau_s .

[IMAGE]
Figure 8-12. KE vs. time for fcc argon with large \tau_s

. 8.4. Conclusions

Since rapid convergence is almost always desired over slower convergence, the shortest possible \tau_s should be used. From the argon cases, it is apparent that there is a limit to how short this can be, however. Between the twin considerations of integrability of the Nosé oscillator and allowing the system to vary somewhat in KE to maintain ergodicity, a suitable choice of \tau_s appears to be 10 times the normal dynamics timestep, or 0.01 ps for typical systems involving hydrogen atoms. Using such a value, convergence to the desired temperature can be achieved in a few hundred timesteps.
Next / Previous / Table of Contents
Kian-Tat Lim