PHYSICAL REVIEW A VOLUME 38, NUMBER 8 Molecular-dynamics pressure thermostat OCTOBER 15, 1988 Luis F. Rull and Juan J. Morales' Departamento de F&'sica Teorica, Facultad de FI'sica, Universidad de Sevilla, 41080 Sevilla, Spain Soren Toxvaerd Chemistry Laboratory III, H. C. grsted Institute, University of Copenhagen, DE 2100 -Copenhagen g, Denmark (Received 21 April 1988) The influence of the pressure thermostat on the equilibrium properties in the Nose-Hoover (NH) method is evaluated. Three tests were made to choose the isobaric relaxation time ~~. For the present system (N =108 particles, po'=0. 75, and kT/a=0. 9) the lowest value for which ~~ is Gaussian was -0. 1, which is twice the mean collision time ~. For ~~ &0. 1 the NH method obeys the general Langevin equation giving the same Mori coefficients, and for v~ & 0. 1 the method breaks down. I. INTRODUCTION In traditional molecular-dynamics (MD) computer simulations, the trajectories of N particles are obtained for constant energy and momentum, the calculations representing a sample from a microcanonical ensemble MD(E, V, N). Recently some computer techniques have been developed which reproduce systems under con- straints by introducing reservoirs which couple with all N particles at the same time. For example, if the simulation is done at constant temperature T, one can simply keep the kinetic energy constant by scaling the velocities of the particles and letting the total energy fluctuate. If the simulation is done at constant pressure p, the volume fluctuates and the average volume is determined by the balance between the internal pressure and the externally imposed pressure p, „. One of those methods, due to Nose' and Hoover (NH), reproduces both the canonical MD(T, V, N) and isothermal-isobaric MD(T, p, N) proba- bility density in the phase space of an N-body time- reversible classical mechanical system in equilibrium. This means that the NH equations obey a generalized Liouville equation. In the NH method the coupling of the particles with the corresponding reservoirs is done by introducing iso- thermal g and isobaric e thermodynamic friction coefficients. These coefficients are related to the corre- sponding relaxation times of the system ~T and ~ . Evans and Holian derived the equilibrium fluctuation expres- sion for a linear response by an NH thermostat and found that in the thermodynamic limit this response is the same as that of the corresponding Gaussian isothermal system. In other words, the Gaussian thermostat makes use of a Gaussian constraint: g has a Gaussian distribution in the NH theory. In a previous paper we solved the NH equations at constant temperature and pressure for a two-dimensional Lennard-Jones system near the melting and freezing points. The MD(T, V, N) and MD(T, p, N) results of the solid-fluid transition showed that the T-p thermostat works well. In this article we present three tests to show how the pressure thermostat works using MD simulation with the NH procedure. The first test is the calculation of the Mori coefficients, because if the NH equations obey a Liouville equation, they also obey the general Langevin equation exactly, i. e. , the constraints have no influence on the time behavior and thus should give the same Mori coefficients. The influence on the equilibrium time behav- ior is most accurately determined in this way (instead of by self-diffusion) since the first three coefficients can be obtained directly during the MD calculations using a fifth-order predictor-corrector algorithm. Likewise, the steady equilibrium probability function derived by NH is — (Di r T)I2 NpT e where D is the dimension of the system. The isobaric friction coefficients(r ) is defined as e= VIDV, where Vis the volume of the system. As pressure is an intensive property, the relevant relax- ation time with p — pex is the thermostat relaxation time ~' defined as (4) Thus e=b, Clap where bC=(p — p, „)IpT is the excess compressibility factor. The e(rp) will have a Gaussian distribution, although some limit value of ~ must exist for which e(rp ) is no longer Gaussian. The second test will show that the thermostat breaks down when the 38 4309 1988 The American Physical Society