Isothermal molecular-dynamics calculations

PHYSICAL REVIEW BVOLUME 32, NUMBER 91NOVEMBER 1985 Isothermal molecular-dynamics calculations Luis F. Rull Departamento Ffsica Teorica, Facultad de Fi'sica, Universidad de Sevilla, 41012Sevilla, Spain Juan J. Morales and Francisco Cuadros Seccion de FI'sicas, Facultad de Ciencias, Universidad de Extremadura, 06071 Badajoz, Spain (Received 30 January 1985; revised manuscript received 28 May 1985) %e have performed long-time runs of molecular-dynamics computer simulations of atwo-dimensional Lennard-Jones system, without any scaling procedure. The thermodynamic properties show spontaneous fluctuations except when the system is far from the melting zone. In arecent paper, 'Toxvaerd studied critically aspecial computational technique: "isothermal-isobaric" molecular dynamics (T-p MD) (used by Abraham and Koch' to investigate melting in atwo-dimensional Lennard-Jones system). By using the calculation of the isothermal compressibility from the fluctuations of volume, Toxvaerd established that the scaling procedure was incorrect and cannot be used to obtain definitive statements about the nature of the phase transition. Toxvaerd noted that the T-p MD calculations, inspired by the N-p-T-ensemble Monte Carlo method, are obtained from Laplace transforms of the canonical partition function, and this transformation has to be performed for constant extensive state variables, excluding the volume. In an analysis of the ensembles used in moleculardynamics simulations, Lado found that, in the classical limit, for interacting systems with hard-core interactions, Monte Carlo and molecular-dynamics calculations should yield the same equation of state without corrections, i.e., BlnQ~ I (PP )NvT+-- V()V (IS/ )N- v-Ei where the configurational integral Q~ vis defined by Q~ r=Vdr" 8'(r"), where ~( ~) 0, U(r )=oo, 1, U(r )=0, and, as usual, P=(ksT) ', with the only assumptions being the basic postulates of equilibrium statistical mechanics. The left-hand side (N-V-T) corresponds to the Monte Carlo calculation, that is, with the number of particles, volume, and temperature constant, and the right-hand side (N VE) correspo-nd-s to the molecular dynamics calculation, with the number of particles, volume, and energy constant. Nose4 has recently made an exhaustive study of the ensembles in molecular dynamics to reproduce both the canonical and the isothermal-isobaric probability densities in phase space. The physical system of interest consists of N particles, to which an external macroscopic variable and its LLI 1 LLI l— ~OQ ~rl p++M+ 01~+0 oe++I ~ ~~~~~~ 0.7,. 20 40 60 '80 10 TIME STEPS FIG. 1. Reduced temperature vs time for atwo-dimensional liquid at pr~=1.00. 32 6050 1985 The American Physical Society BRIEF REPORTS 6051 1,3' eee4ee ~eeeeeeeeev v~e eeeee4eee eeee eee ee eeeeeeeeeeeeeee eel $ eeeee ee.ee ~ee ee eeee ~ T1 0.7 20 60 80 10 TINE STEPS FIG. 2. Same as Fig. 1, but for atwo-dimensional solid very close to the melting zone. The density is prm=1. 14. conjugate momentum are added, permitting the energy of the system to fluctuate in away similar to Haile and Gupta's results with athermal bath. ' In this Brief Report-, we present some results of two long runs in two dimensions close to and far from the melting zone, which were performed in order to investigate the fluctuations of the thermodynamic variables during the time of. evolution of the system. Our simulations have been carried out without any scaling procedures (for pressure and temperature). Just as it is incorrect to perform asimulation using the constant-pressure method in the melting zone, it is also incorrect to use the constant-temperature method, because there are problems related to fluctuations in space and time independent of the nature of the phase transition (first or second order). In our calculations we have chosen asystem of %=256 particles interacting with atruncated Lennard-Jones potential at r, =2.5r~ (ris the position of the potential minimum) . The initial positions of the particles formed atwodimensional triangular lattice of densities pr2=1.00 and pr~=1.14. To integrate the equations of motion we used the algorithm proposed by Toxvaerd. This algorithm is very accurate, but introduces a drift in the energy because the equations are not reversible in time; this energy-drift problem can be solved by rescaling the velocities of the particles through the heat capacity. In our simulations this was not necessary, because the maximum drift is less than 10 4, and consequently the increment in the temperature is kshT/e — =10 ', and, in addition, the usual periodic boundary conditions were used. Both systems evolved in exactly the same way. After every ten time steps (updating the table of the nearest neighbors), the velocities were renormalized to give the re- ~e ee ~e~eeeee ~.e~ee ee ~ ee ~~ee ~~ee~ee ee e 40 10 TINE STEPS 60 80 FIG. 3. Reduced pressure vs time for atwo-dimensional liquid at pr~ =1.00. 60S2 BRIEF REPORTS 32 10 I 9"~ ~Oy ~y~~ ~0 ~y~O ~~0 ~~~y ~0~~y ~4~01 ~l40 ~y$0 ~y ~g~~~O QQ ~~O+ eqi ~OO ~ee ~~QO ~I 20 1Q TI&E STEPS FIG, 4. Same as Fig. 3, but for atwo-dimensional solid very close to the melting zone. The density is prm=1. 14. duced temperature AT/a=1. 00 (e is the minimum of the I.ennard-Jones potential). The rescaling procedure was performed during the first 8000 time steps. After that the evolution of the system was free, without any scaling procedure, the thermodynamic properties being then obtained as time averages for each 800 time steps. The virial pressure was obtained for the truncated Lennard-Jones potential, not for afull Lennard-Jones potential. The state points that we have chosen correspond to the fluid state (pr =1.00) and the solid state (pr' =1.14), the latter very close to the melting zone located thermodynamically by Toxvaerd. 'Toxvaerd calculated the point where fluid and solid have equal chemical potentials corresponding to the temperature ks T/e =1.00, the difference in our case being that we have performed long runs (40000 time steps for afluid and 80000 time steps for asolid). Figures 1-4 show the results of our simulations. Each point in these figures corresponds to an average over 800 time steps. In the liquid system (Figs. 1and 3) the temperature and pressure show the "regular" fluctuation typical of the liquid system, and the system is in equilibrium once free evolution commences. But for the system close to the melting zone, equilibrium is reached after 20000 time steps (Figs. 2and 4), and the temperature spontaneously drops twice (arrows), and consequently the pressure rises. Toxvaerd' wondered if the equilibrium state of the system could be obtained in atime interval of 10000 time steps. -Our answer is no, because in the melting zone (or close to it) the system spontaneously goes back and forth between two different points of state. It seems evident that far from atransition in liquid or solid systems, computer simulations using the T-p MC method are reliable because the fluctuations are small. However, fluctuations play an important role in the melting zone; so to obtain adefinitive statement about the nature of the phase transition one should use MD simulations of free evolution, without any scaling procedure, i.e.,without temperature or pressure constraints. These constraints imply that the system is not isolated, but is in contact with an energy reservoir and that the energy is transferred by ageneralized force.' The cooperation of the University of Sevilla Computer Center is gratefully acknowledged. S.Toxvaerd, Phys. Rev. B29, 2821 (1984). F. F. Abraham and S. %.Koch, Phys. Rev. B29, 2824 (1984), and references cited therein. F. Lado, J. Chem. Phys. 75, 5461 (1981). 4S. Nose, Mol, Phys. 52, 255 (1984); J. Chem. Phys. 81, 511 (1984). 5J. M. Haile and S. Gupta, J. Chem. Phys. 79, 3067 (1983). S. Toxvaerd, J. Comput. Phys. 47, 444 (1982). 7S. Toxvaerd, Phys. Rev. A24, 2735 (1981).