Calculation of the Melting Entropy of Argon at Constant Volume Using Molecular Dynamics

Abstract

The melting entropy values of argon at constant volume (Δ_vS) as a function of density were calculated using molecular dynamics simulations. These calculations employed cells with N = 864 molecules. The resulting entropy of melting per particle, Δ_vS/N, was found to be essentially constant (with some slight effect of the crystal form) and nearly equal to 0.5 k, where k is the Boltzmann constant. The melting entropy values at constant pressure (Δ_pS) were also determined and it was observed that Δ_pS > Δ_vS. In addition, variations in the vibrational frequency in the liquid were compared with fluctuations in the solid near the melting point as a means of further investigating the melting entropy.

1 INTRODUCTION

The melting entropy of a substance, ΔS, is typically determined under constant pressure, thus giving Δ_pS [1]. However, such analyses will include the effect of volume changes other than those associated with the transformation of the structure due to melting. For this reason, the present work calculated melting entropy values for argon at constant volume using molecular dynamics (MD) simulations [2] and assuming that intermolecular interactions can be represented by the Lennard-Jones (LJ) function [1]. To the best of the author’s knowledge, there are no simple rules associated with the determination of Δ_pS, although a new rule is presented herein regarding Δ_vS. Specifically, this study found that the melting entropy at constant volume will be less than the melting entropy under constant pressure: Δ_pS > Δ_vS. In the present work, the Δ_vS of argon was also determined to depend somewhat on the crystal form (face centered cubic (FCC), hexagonal close packed (HCP) or body centered cubic (BCC)) and to have a value of approximately 0.5 k, where k is the Boltzmann constant.

Based on calculations using a standard number of molecules, N( = 864), the limiting value of Δ_vS was estimated to be 0.49 k as N approaches infinity at a typical density.

Herein, the entropy of super-cooled liquid argon is also compared with that of solid argon, and the temperature-dependence of the entropy of the solid phase is discussed based on the mean square displacement (MSD). This work also examined pre-melting near the melting point, again using the concept of MSD. Finally, fluctuations in the vibrational frequencies in the liquid state were compared with those occurring in the solid near the melting point, as a means of assessing the melting entropy.

2  MOLECULAR DYNAMICS

In this study, periodic boundary conditions were assumed, and a time increment of 1 fs was employed. The SCIGRESS-ME program was used for MD calculations, and each standard MD run consisted of 1,000,000 steps [3]. An initial FCC configuration was used during standard calculations and the Lennard-Jones (LJ) function was always assumed, expressed as:   

$u\left(r\right)=4\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]$(1)

, where ε is the depth of the potential well and σ is the separation at which u(σ) = 0. The constants ε and σ have units of energy and length, respectively, and numerical values for these constants in the case of argon are provided in Table 1 [3].

Table 1. Lennard–Jones parameters for argon [3].

(ε/k)/K	ε/10−21 J	σ/10−10 m
125	1.73	3.428

In these calculations, the cut off distance was half the length of the basic cell, and an NTV ensemble was used to determine variations in the internal energy, Δ_vU, and the melting temperature, T_m, using standard MD simulations [2]. The Δ_vS values were calculated as [1]:   

${\Delta }_{v}S=\frac{{\Delta }_{v}U}{T_m}$(2)

.

3  MELTING ENTROPY AT CONSTANTVOLUME

Figure 1 presents plots of the average potential energy per particle, <E_p/N > , as functions of the temperature at a density, N/V, of 0.97 σ⁻³. In obtaining these plots, the internal energy, U, was calculated as:   

$U=\frac{3}{2}NkT+⟨E_p⟩$(3)

.

Figure 1.

 The average potential energy per particle, <E_p>/N, as a function of the temperature, T, at a number density, N/V, of 0.97 σ⁻³.

One of these plots corresponds to solid and liquid phases starting from the FCC configuration, while the other is for argon as a super-cooled liquid. The super-cooled liquid data were obtained using a random initial configuration representing the final configuration at T = 2.4 ε/k. In the first plot, melting is observed at approximately T = 1.39 ε/k, as indicated by the abrupt increase in the average potential energy.

The corresponding plots of pressure, p, against temperature are presented in Figure 2. In addition to freezing, the cooled structure also exhibits an unexpected temperature dependence. Here, the average potential energy of the cooled structure has almost the same value above and below the melting temperature (Figure 1). In contrast to the results in Figure 1, it is also evident that the two pressure plots deviate from one another in the same temperature regions. Nevertheless, a sharp rise is again observed around the freezing point.

Figure 2.

 The pressure, p, as a function of the temperature, T, at a number density, N/V, of 0.97 σ⁻³.

A more detailed plot focusing on the region around the melting point is shown in Figure 3. The melting temperature, T_m, was calculated as the average of the temperature values at which the abrupt changes in potential energy occurred. The melting temperature values obtained using the three different crystal structures (HCP, FCC and BCC) are plotted as functions of the number density, N/V, in Figure 4. At the greatest density, the FCC plot exhibits the highest T_m, while the BCC plot has the lowest value. This Figure demonstrates that T_m varies significantly as the density is increased but that the differences resulting from changing the crystal structure are minimal. Figure 5 plots the pressure values of solid and liquid argon at the melting point as functions of density, while Figure 6 shows the average potential energy values per particle, <Ep>/N, at the melting point as functions of density.

Figure 3.

 The average potential energy per particle, <E_p>/N, as a function of temperature, T, in the vicinity of the melting point, at a number density, N/V, of 0.97 σ⁻³. The red arrow indicates the averaged T_m.

Figure 4.

 The melting temperature, T_m, as a function of the number density, N/V, for hexagonal closed packing (HCP), face centered cubic (FCC) and body centered cubic (BCC) structures.

Figure 5.

 The pressure, p, in the solid phase at the melting point as a function of the number density, N/V, for HCP, FCC and BCC structures.

Figure 6.

 The average potential energies per particle, <E_p>/N, in the solid and liquid phases at the melting point as a function of the number density, N/V, for HCP, FCC and BCC structures.

The melting entropy values at constant volume with respect to density are plotted against density in Figure 7. These values were calculated using equation (2) in conjunction with the relationship:   

${\Delta }_{v}U=⟨E_p\left(L\right)⟩-⟨E_p\left(S\right)⟩$(4)

, where L and S indicate the liquid and solid-states, respectively, as shown in Figure 6. In each case, the melting entropy at constant volume was found to be approximately 0.5 k, and the averages were as follows.   

$\begin{array}{c}{⟨{\Delta }_{v}S/N⟩}_{HCP}=0.53\text{ k}\\ {⟨{\Delta }_{v}S/N⟩}_{FCC}=0.45\text{ k}\\ {⟨{\Delta }_{v}S/N⟩}_{BCC}=0.38\text{ k}\end{array}$(5)

Figure 7.

 The melting entropy per particle at constant volume, Δ_vS/N, in the solid phase as a function of the number density, N/V, for HCP, FCC and BCC structures.

Because the HCC and FCC are both close packing structures, it was expected that the melting entropies for these configurations would be almost equal. Based on this consideration and the data in Figure 7, the melting entropy can be summarized as:   

${\left({\Delta }_{v}S/N\right)}_{FCC\&HCP}\approx \left(0.49\pm 0.04\right)\text{ k}$(6)

.

The melting entropy at constant volume of the FCC state is plotted as a function of the inverse of N in Figure 8 at a density d = 0.97 σ⁻³.The limiting value of Δ_vS/N at an infinitely large value of N is:   

$\underset{N\to \infty }{\lim }\frac{{\Delta }_{v}S}{N}=0.49$(7)

.

Figure 8.

 The melting entropy per particle at constant volume, Δ_vS/N, as a function of the inverse of the number of molecules in the basic cell, 1/N, at a number density, d, of 0.97 σ⁻³.

4 MELTING ENTROPY UNDER CONSTANTPRESSURE

The melting entropy values under constant pressure with respect to density, Δ_pS/N, are plotted against density in Figure 9. These were calculated using both the equation of state [4] and NPT-MD methods [2]. The equation of state provides a more reliable phase boundary result than the present NPT-MD, because the equation of state is compiled from many simulation results. Although the NPT-MD results deviate from the equation of state results, both plots show the same qualitative trend. It is also evident that the Δ_pS/N values are larger than the Δ_vS/N values, such that the relationship below may be written.   

${\Delta }_{p}S>{\Delta }_{V}S$(8)

Figure 9.

 The melting entropy per particle, ΔS/N, as a function of the number density, N/V, of the solid phase for an FCC structure. The equation of state (EOS), NTP and NTV results are compared. The values of EOS and NTP correspond to the phase transition under constant pressure. The results of NTV are calculated MD where N, T and V are constant.

It is also evident that Δ_pS/N decreases as the density increases and that this rate of decrease is greatest in the low density region.

The relationship between Δ_pS/N and Δ_vS/N is summarized in Table 2 and Figures. 10 and 11. Figure 10 describes the density dependence of the average potential energy, <E_p/N > , at T = 0.692 ε/k calculated by NTV MD. The pressure, p, is plotted against the density at the same temperature in Figure 11. These plots demonstrate melting near the triple point under constant pressure. In the intermediate part of each plot, melting at a constant volume is also observed around a density, N/V, of 0.86 σ⁻³. Table 2 provides the Δ_pS/N values at T_m = 0.692 ε/k and the associated numerical components. This table also provides results for another high density scenario, for which T_m = 4.5 ε/k. These data demonstrate that Δ_vS is the core part of Δ_pS.

Table 2. Numerical parameters associated with Δ_pS near the triple point and under high pressure.

Tm/(ε/k)	0.692	4.50
ΔpS/Nk	1.534	1.021
ΔvS/Nk	0.439	0.488
ΔotherU/TmNk	1.130	−0.160
ΔpV/TmNk	−0.036	0.692

Figure 10.

 The average potential energy per particle, <E_p>/N, as a function of the number density, N/V, near the triple point, T = 0.692 ε/k.

Figure 11.

 The pressure, p, as a function of the number density, N/V, near the triple point, T = 0.692 ε/k.

5 ENTROPY BELOW THE MELTING POINT

In this work, the entropy was calculated as a function of temperature using the following equation.   

$\Delta S={\int }_{Ti}^{Tf}{\left(\frac{dU}{dT}\right)}_{V}dT+{\left(\frac{\Delta U}{T}\right)}_{T_m}$(9)

The result is shown in Figure 12. Herein, the entropy values of the super-cooled liquid are plotted only in the region over which the super-cooled liquid persists in the MD simulations. At lower temperatures, incomplete solid structures were identified. It is also evident that the entropy of the super-cooled liquid at T = 0.8 ε/k is greater than that of the solid, while the entropy difference between solid and liquid is close to Δ_vS/N(see Table 3, which also includes the value calculated by thermodynamic integration [5]. Thismethod gives the entropy measured from the no-interaction system.). It should be noted that some of this difference in the entropy values results in Table 3 from errors in the numerical integration.

Figure 12.

 The entropy per particle, S/N, as a function of the temperature, T, at a number density, d, of 0.97 σ⁻³.

Table 3. The entropy values of super-cooled liquid argon, where ΔS = S,random − S,FCC, as determined using two approaches.

Method	T/(ε/k)	S/Nk, random	S/Nk, FCC	ΔS/Nk
Eq. (9)	0.80	4.50	4.10	0.40
Thermodynamic Integration	0.80	−4.65	−4.19	0.46

The temperature dependence of entropy of the solid argon can be determined based on the mean-square-displacement (MSD)$⟨{\left(\Delta r\right)}^2⟩$, which is associated with the average radius, r, of the available vibrational space as in the equation below.   

$r={\left(⟨{\left(\Delta r\right)}^2⟩\right)}^{\frac{1}{2}}$(10)

In addition, the following one-body approximation was used to estimate the entropy.   

$\Delta S=k\log \left(\frac{V}{V_0}\right)$(11)

Here the relationship between V and r is as below.   

$V\propto r^3$(12)

The results obtained from the above process are compared with the entropy values determined using equation (9) in Figure 13. It can be seen that the approximations included in equations (10) to (12) closely reproduce the S values outside of the region near the melting temperature.

Figure 13.

 The natural log of the available vibrational space, r³, as a function of temperature, T, at a number density, d, of 0.97 σ⁻³. The entropy values, S, were determined using Eq. (9) using data from Figure 12.

At a temperature value close to the melting point (T = 1.30 ε/k), the argon molecules occasionally move to other sites, as shown in Figure 14. In contrast, below this value, at T = 1.2 ε/k (which is 14% lower than the melting point), each molecule oscillates around an equilibrium position. In Figure 15, MSD values are plotted as functions of time, t, in units of τ. This unit for a Lennard-Jones particle with mass m is defined as follows.   

$\tau =\sigma {\left(\frac{m}{\epsilon }\right)}^{\frac{1}{2}}$(13)

Figure 14.

 Trajectories of molecules at T = 1.20 and 1.30 ε/k, at number density, d, of 0.97 σ⁻³.

Figure 15.

 The mean-square-displacement (MSD),$⟨{\left(\Delta r\right)}^2⟩$, as a function of time, t, in units of τ, at T = 1.44, 1.30 and 1.20 ε/k.

MSD values were compared at several temperatures (T = 1.20, 1.30 and 1.44 ε/k) and an excited state was identified at T = 1.30 ε/k. This state was found to have a larger entropy value than that determined using the approximations in equations (10) to (12), as shown in the equation below.   

${\left(\frac{\Delta S}{N}\right)}_{excitation}=0.6\text{ k}$(14)

This different entropy value is tentatively attributed to the effect of pre-melting excitation entropy. Despite the difference, it is still quite close to the value determined for the melting entropy at constant volume.

6 ANALYSIS OF VIBRATIONAL SPECTRA

Vibrational motion is frequently observed in the liquid state [6]. If we assume a harmonic oscillator with frequency ω at temperature T, the deviation, Δx, from the equilibrium point can be described as a function of mass, m, as follows.   

${\left(\Delta x\right)}^2=\frac{kT}{m{\omega }^2}$(15)

If there is a range of frequencies and the deviation, δω, from the average, <ω>, is small, Δx can be approximated by the following expression.   

${\left(\Delta x\right)}^2=\frac{k_{B}T}{m{\left(⟨\omega ⟩+\delta \omega \right)}^2}\approx \frac{k_{B}T}{m{⟨\omega ⟩}^2}{\left(1-\frac{\delta \omega }{⟨\omega ⟩}\right)}^2\approx \frac{k_{B}T}{m{⟨\omega ⟩}^2}\left(1-2\frac{\delta \omega }{⟨\omega ⟩}\right)$(16)

In addition, the change in the square of Δx can be written as a function of ω as follows.   

$\begin{array}{c}⟨{\left({\left(\Delta x\right)}^2\right)}^2⟩-{⟨{\left(\Delta x\right)}^2⟩}^2\approx 4{\left(\frac{k_{B}T}{m{⟨\omega ⟩}^2}\right)}^2\frac{⟨{\left(\delta \omega \right)}^2⟩}{{⟨\omega ⟩}^2}\approx 4{\left(\frac{k_{B}T}{m{⟨\omega ⟩}^2}\right)}^2\frac{⟨{\omega }^2⟩-{⟨\omega ⟩}^2}{{⟨\omega ⟩}^2}\\ \end{array}$(17)

If the value of <ω> in the solid and liquid states at a constant volume and temperature is essentially constant, the following expression holds true.   

$\frac{{\left\{⟨{\left(\Delta x^2\right)}^2⟩-{⟨\Delta x^2⟩}^2\right\}}_L}{{\left\{⟨{\left(\Delta x^2\right)}^2⟩-{⟨\Delta x^2⟩}^2\right\}}_S}\approx \frac{{\left\{⟨{\omega }^2⟩-{⟨\omega ⟩}^2\right\}}_L}{{\left\{⟨{\omega }^2⟩-{⟨\omega ⟩}^2\right\}}_S}$(18)

The frequency spectra at T = 1.39 ε/k and d = 0.97 σ⁻³ are shown in Figure 16. Here, the difference between <ω> in the solid and liquid states is approximately 1%. The approximation in equations (10) to (12) predicts the following difference in entropy values.   

${\left(\frac{S}{N}\right)}_{Liquid}-{\left(\frac{S}{N}\right)}_{Solid}\approx 0.6\text{ k}$(19)

Figure 16.

 The frequency spectra based on I (ω) values at T = 1.39 ε/k and d = 0.97 σ⁻³.

This value is close to the melting entropy at constant volume.

7 CONCLUSIONS

A new rule regarding Δ_vS(the melting entropy at constant volume) was identified during this work; the melting entropy per particle is almost constant as a function of density for a Lennard-Jones system and can be approximately by 0.5 k, where k is the Boltzmann constant. The author suggests that this rule should be validated by macroscopic experiments in future.

Acknowledgment

The author would like to thank the Research Center for Computing and Multimedia Studies at Hosei University for the use of computer resources.

REFERENCES

• [1]   P. W. Atkins, Physical Chemistry, Oxford Univ. Press, Oxford (1998).
• [2]   M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford (1992).
• [3]   http://www.scigress.com/ (accessed 2017)
• [4]   M. A. Barroso, A. L. Ferreira, J. Chem. Phys., 116, 7145 (2002). doi:10.1063/1.1464828
• [5]   D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, San Diego (1996)
• [6]   Y. Kataoka, Y. Yamada, J. Comput. Chem. Jpn., 13, 130 (2014). doi:10.2477/jccj.2014-0002