Molecular Sciences Thermodynamic, Structural and Transport Properties of Lennard-jones Liquid Systems. a Molecular Dynamics Simulations of Liquid Helium, Neon, Methane and Nitrogen

Molecular dynamics calculations are carried out in order to find the properties of Lennard Jones liquids in different state points of their phase diagram. The spherical shape and the stability of the helium, neon, methane and nitrogen make the liquids easily accessible to numerical simulation. Thermodynamic, structural, and transport properties are studied and compared with both experimental data and recent theoretical investigations. In the present work, up to 22 state points are covered, some of which are near or at the triple point. It will be shown that the classical approach leads to data that are in very good agreement with experiments and other types of calculations. At high temperatures and low densities, we observe a decrease in the uncertainties in the stress autocorrelation function by increasing the number of iterations.


Introduction
The prediction of thermodynamic properties of liquids and mixtures is still an issue of interest in simulation.The literature shows that these systems are widely investigated in both experimental and theoretical ways [1][2][3][4][5].X-rays scattering of liquid methane near the triple point (90.7 K°) [6][7][8][9][10], and liquid nitrogen [11] were studied.Theoretically new simulations with path-integral formalism have been conducted to obtain thermodynamic properties, and the quantum radial functions involving Feynmann-Hibbs potentials were used for liquid helium [12].The equilibrium and non-equilibrium properties of methane were studied with Monte-Carlo (MC) [2,3,13] and Molecular Dynamics (MD) [14,15] simulations.Moreover, various potential models have been proposed for reproducing data for these liquids [16-20, 3, 21, 22].For liquid neon, several simulations involving different approaches were performed: Path-Integral Monte-Carlo (PIMC) [23,24], Path-Integral Brownian Dynamics (PIBD) [24], Monte-Carlo simulations using two effective pair potentials arising from Path-Integral formalism, the quadratic Feynmann-Hibbs(QFH), the Gaussian Feynmann-Hibbs (GFH) and the perturbation theory [25].These works focused mainly on the radial distribution function for a complete structural investigation of these liquids, and recent calculations [3,21,22], taking into account quantum effects, cover only some state points of the phase diagram and offer no way to approach the transport properties accurately.In this work, we propose a complete study over 22 state points and most structural, thermodynamic and transport properties are studied.Statistical accuracy of calculated data are carried out and compared to available data, either experimental or theoretical.

Numerical Model
Molecular Dynamic calculations were performed on a sample of size N = 864 molecules of masse m.The molecules interact by Lennard-Jones pair potential specified by the parameters ε and σ.The units of energy, length, and mass were chosen to be, respectively ε, σ, and m.The corresponding microscopic time scale is The state of our system is specified through the reduced number density and temperature ρ * =ρ σ 3 and T * = k B T/ε.The thermodynamic state points sp(ρ,T) of the systems of interests methane, helium, nitrogen and neon are summarized in Table 1.The two sp of liquid helium (Hesp1, Hesp7) were chosen for comparison with Sesé's work [12].For liquid methane, the sp (Msp1, Msp4, Msp8) were taken from Sesé's work [3].
The other state points were taken to complete our data production for methane in the liquid state.The first sp Nesp1 ( the triple point ) and the corresponding (ρ,T) data of liquid neon were taken from Singer and al studies [24].The last three sp ( Nesp2, Nesp3, Nesp4) belonging to the isotherm 35.05 K° and their data were taken close to those of Raveché and al [26].These four sp were also studied by Sesé [27].Finally, the sp for liquid nitrogen (N2sp1, N2sp2) were taken from Narten and al [11].
The equations of motion were solved using the leap-Frog integration scheme with a constant time step algorithm (∆t = 0.5×10 -14 s) and wherein the temperature T was kept constant by the constraint method [28][29].Periodic boundary conditions around the central cubic box and the minimum image truncation were included in the calculations; long-range corrections were also applied.The starting He(l) 4  5.1 0.5 0.0947 0.2380 Hesp6 He(l) 4  5.1 0.5 0.0483 0.1216 Hesp7 [10] He(l) 4  5 where ν i (t) are the velocitys of atoms at time t.The self-diffusion coefficient D of a particule is evaluated by the Einstein and Green-Kubo formulas [30]: where r i (t) are the positions of the particles at time t.D is proportional to the slope of the mean-square displacement (MSD) at long times.The shear viscosity coefficient η can be evaluated from both Green-Kubo and generalized formulas [31].For soft-body potentials, most simulators have used the Green-Kubo equation [30], which is an integral over the stress autocorrelation function, The component of the microscopic stress tensor is given by The stress autocorrelation function necessarily involves the entire system; consequently, we cannot improve the statistical precision of results for viscosity by averaging over the N particles in the system.
However, for stationary, homogeneous, uniform fluids, the statistical precision can be improved somewhat by averaging over all three terms that result from the stress tensor, ( )

Results and Discussions
The figures presented in what follows are only illustrative of the data that can be obtained through our calculations (r* = r/σ and t* = t/τ denote distances and time in reduced units).Figure 1 shows the relevant structural information for three representative points Msp1, Msp6 and Msp9 of liquid methane.The function g(r) becomes firstly zero at short distances, where repulsive forces prevent overlapping of molecules.When r is close to the collision diameter σ  g(r) increases rapidly to a maximum r = r m corresponding to the first peak.As r increases gradually, g(r) decreases showing that influence of the central molecule is disappearing and there is no order at long distances.The average number of neighbours for a methane molecule is obtained from where the integration limit R is taken as the position of the successive minima in g(r).
For the Msp1, the number of nearest neighbours is 11,05 molecules and the total number of methane molecules up to the second layer is 34,36.At a temperature of T=92K°, the g(r) from X-ray diffraction gives a number of nearest neighbours ≈12 molecules and the number of neighbours in the second layer of about ≈ 55 [6].Sesé [3] has obtained similar radial distribution function to those of the present calculations (Fig. 1) for points Msp1, Msp4, Msp8 with peak positions at different maximum R max1 = 4.05Å and R max2 = 7.75 Å and minimum R min1 = 5.75-5.85Å.The mean square displacement of atom of liquid neon at different state points is reported in Figure 3.
At low densities the curve shows two steps.In the first stage 0.0 < t* < 0.1 the system behaves like a crystal state rather then a liquid state.Arrangement of atoms at short distances shows a certain degree of dependence, which makes one think of permanent correlations existing in the solid state.In the second stage, a change in the slope is observed with a linearly rising mean square displacement; in this case, atoms possess enough kinetic energy.Points along the 35.05K° isotherm of liquid neon ( Nesp2, Nesp3, Nesp4) present the same profile.As time increases, an important gap appears between curves of the isotherm and the curve of triple point ( Nesp1).Equilibrium is reached faster at low densities and high temperatures and the overall runtime of a simulation is more important for the calculation of transport properties at high densities.
The stress autocorrelation function calculated for each points of the phase diagram and for every fluid allow us to work out the shear viscosity η.For the sp Hesp1 of liquid helium, the stress autocorrelation function is illustrated in Fig. 4a according to the three orientations.An increase of the number of iteration by a factor of 10 (from 5000 to 50000) decreases the uncertainties on η(t) as shown in Fig. 4b.The uncertainty vary about between 0.2 and -0.2 through the time (for example, for t*=0.75 it is in order of 0.2 ).At high densities and low temperatures, we observe an increase in the calculated uncertainties.On the other hand, by increasing the calculation time, a good statistical precision of the stress autocorrelation function is achieved .Table 2 shows calculated averages of different thermodynamic properties for every system of interest in different state points of their phase diagram.Different observations can be made.The fluctuations of properties such as the energy of configuration U, the total energy E and the enthalpy H for all states points do not show any changes.Fluid helium shows an unstable character because of the neglect of the quantum effects in this simplified model.We may notice very large fluctuations for liquid methane at points Msp1and Msp2, this is probably because these sp are close to the the triple point ( T=90.7K°).When moving over pressure data column of Table 2, an increase of pressure is Nesp1 Ne(l) 0.6653 0.8084 -0.587±0.079-4.925±0.018-5.923±0.017-5.651±0.1160.0395±0.00051.95±0.0001Nesp2 Ne(l) 0.9517 0.7528 -0.250±0.080-3.664±0.022-5.091±0.022-4.010±0.1310.0927±0.00131.80±0.0001Nesp3 Ne(l) 0.9517 0.7246 -0.279±0.102-3.643±0.031-5.070±0.031-4.030±0.1720.0911±0.00121.75±0.0001Nesp4 Ne(l) 0.9517 0.6877 -0.404±0.070-3.411±0.028-4.840±0.025-4.000±0.1260.1029±0.00111.52±0.0001N2sp1 N2(l) 0.6496 0.8783 0.218±0.216-5.404±0.113-6.383±0.050-5.156±0.3130.0216±0.00042.55±0.00005N2sp2 N2(l) 0.7874 0.6414 -0.630±0.053-3.484±0.021-4.665±0.021-4.465±0.0960.1040±0.00131.06±0.00035observed from point Hesp1 to Hesp6 except for state point Hesp7 where the density is smaller.For liquid methane, the pressure rises from point Msp1 to Msp2 and drops down for the point Msp3.The variation of pressure is irregular between points Msp4 and Msp9.The three points along the 35.05K° isotherm of liquid neon, the pressure drops down with decreasing densities except at the triple point where pressure makes fail at this rule.From N2sp1 to N2sp2 of liquid nitrogen, the pressure decreases remarkably with increasing temperature.The diffusion coefficients values reported in table 2 clearly show that the simulated systems are in their liquid states.At high temperatures, the molecules possess important kinetic energy to enlarge the diffusion of molecules in liquid.The diffusion coefficient fluctuations turn around a constant value.The results obtained are in close agreement with the values reported by Erpenbeck [32], Levesque and Verlet [33] and Schoen and Hoheisel [34] calculations of Lennard-Jones liquid systems.
As expected, the shear viscosity decreases with increasing temperature.It is obvious that properties calculated at high temperatures and low densities require a correlation of collective transport properties over long simulation time.One can easily work out the activation energy from the diffusion coefficient versus temperature curve, using the Arrhenius low as illustrated in Figure 5.The calculated activation energy for methane is 3.34 K.J.mol -1 and is found close to the activation energy measured by NMR (3.7 K.J.mol -1 ) [35] and Schoen and al simulations [15] (3.2 K.J.mol -1 with a spherical symmetric potential and 3.8 K.J.mol -1 for centre site-site potential).
In order to obtain a complete analysis of our results and in order to demonstrate the reliability of classic model, a comparative study of our data with experimental [12,24,27] and theoretical works [3,12,27] is performed as reported by Table 3.For this purpose, we have reported data obtained by Sesé [3,12,27] from Monte-Carlo simulation involving Feynmann-Hibbs effective potential methods; Quadratic (QFH) and Gaussian(GFH) or the path-integral method (PIMC, PIBD, PIBC, MS-PI).A good agreement is found with the quantum approach except data of reference [3] using WK(h 2 ).
The pressure value shows the chronological order: PIMC > QFH ≅ WK ≅ MCLJ > DMLJ.The FH quantum effects and PIBC method increase total and potential energies as compared with classical LJ values.On the other hand, the three methods give completely different pressures.From all these results, one can draw the following observations: with decreasing densities, (a) the FH methods lead to thermodynamic properties remarkably close to experiment, (b) QFH and GFH approaches yield similar results, QFH being somewhat much more reliable.The agreement fails as the density increases, which can eventually lead to very poor estimates of some properties (pressure at Nesp4).The large pressure fluctuations make comparisons difficult to follow, nevertheless, pressures appear to be very large for quantum models when compared to classical model in agreement with Hansen and Weis[37], Thirumalai and Hall [36] and Singer and Smith [24].From this Table 3, it appears also that the basic three methods compared (classical Molecular Dynamics, quantum Monte-Carlo and path-integral computation) deviate very strongly from each other for the case of Helium, where they show random behavior, while they lead to fairly consistent results for methane.This suggests that none of the methods succeeds in capturing the strong quantum effects that are characteristic for liquid Helium.All these approaches are semi classical finitetemperature techniques, which may provide a measure of improvement over quantum-corrected classical results (liquid neon), but which are not applicable in the strongly quantum-mechanical lowtemperature limit (liquid helium) [38,39].

Conclusion
Classical molecular dynamics simulation of Lennard Jones fluids is possible within different state point of the phase diagram.Lennard-Jones liquid energies are not affected by translational quantification because the static errors in the quantum models explain the experimental data.This is the Lennard-Jones fitting consequence [2] witch take into account one way or the other quantum effects.Quantum effects are necessary to explain the different behaviour near the triple point.Our work offers the possibility to study the transport properties with a good precision as for the diffusion coefficient or the shear viscosity.Our results compare favorably with both experiment and more sophisticated theoretical approaches that incorporate quantum effects.

Figure 4 :
Figure 4: Stress autocorrelation function of liquid helium for state point Hesp1.a) An overall runtime of 5000 steps.b) An overall runtime of 50000 steps.

Figure 5 :
Figure 5: Diffusion coefficient as a function of temperature for the pure liquid methane.

Table 1 .
Type of fluids with the value of the Lennard-Jones parameters, and thermodynamic state points investigated in this work.Densities and temperatures are in real and reduced units.

Table 3 .
A Comparative study of data obtained with different investigations for liquid helium (Hesp1, hesp7), liquid methane (Msp1, Msp4, Msp8) and liquid neon (Nesp1-Nesp4).Experimental data are reported as well.Data are reported in reduced units.