An Extended van der Waals Equation of State Based on Molecular Dynamics Simulation

Yosuke KATAOKA and Yuri YAMADA


Return

1 Introduction

The van der Waals EOS is well known in physical chemistry because it has a very simple form and explains the characteristics of a real gas [1]. The following pVT relation also explains the vapor-liquid critical point:

here a and b are constants that characterize the molecules. The variables N, V, and T are the number of molecules, the volume, and the temperature of the system. The Boltzmann constant is denoted as k. The vapor-liquid transition point is estimated by the Maxwell construction in the van der Waals loop region in the pVT relation, where the unstable state is also found [1]. This loop appears because van der Waals theory assumes a single phase. Molecular simulation also treats the phase as a single phase in a conventional method. This loop is found in molecular simulations with 100-10,000 molecules [2 - 6]. The loop appears in molecular simulation because the molecular system has uniform structure on the macroscopic scale under the periodic boundary condition. In this sense, the van der Waals EOS is important in the understanding of the molecular simulation results.
On the other hand, the UVT relation in the van der Waals EOS is too simple [6]. The internal energy U is described as follows in the van der Waals EOS [7]:

The heat capacity under the constant volume CV of real liquid depends on temperature, as confirmed by macroscopic and molecular observations [1].
The present study proposes an extended form of the van der Waals EOS to explain the pVT and UVT relations of the gas and liquid phase of the spherical molecules by a simple treatment of the attraction part in the canonical ensemble partition function. The formulation follows the perturbation theory with the hard sphere system as a reference system [7]. Although the original van der Waals EOS includes only the averaged value of the attraction energy, the extended van der Waals EOS represents the attraction part by a cluster expansion in an intuitive manner, whereby the Lennard-Jones interaction is assumed as typical in the spherical molecules.
The extended van der Waals EOS for the Lennard-Jones system has only three parameters, two of which are the Lennard-Jones parameters. The remaining parameter is the effective volume of the hard sphere of the reference system. The pVT and UVT relations will be compared with those of the molecular dynamics [MD] simulation on the Lennard-Jones system.
The present paper follows the statistical mechanical theory on the original van der Waals EOS [7], which is then extended to explain, at least qualitatively, the UVT relation in a broad region in the temperature-density domain. Finally, the theoretical and simulation results are compared.

2 Theory of van der Waals EOS

The theory of van der Waals EOS [7] is reviewed in this section and will be extended in the following section. A molecular system with N spherical molecules of mass m at temperature T and volume V is considered. The canonical partition function Q(N,V,T) is given as follows:

Here, the potential energy of the system is written as U(r1, r2, ..., rN). The quantity ZN is called the configuration integral. In the case of a perfect gas, ZN is given as follows, because no interaction is assumed:

The thermodynamic quantities are written by the partition function as follows:

According to this formula, in the case of a perfect gas, the pressure is given as follows:

The interaction energy is approximated by the van der Waals approximation. The statistical mechanical theory on van der Waals EOS adopts the perturbation approximation. The total potential energy is divided into two parts:

where and are the potential energy of the reference system and that of the perturbed system. The quantity b is the inverse of temperature: b = 1/kT. The configuration integral is as follows:

This equation can be rewritten as follows:

The reference system is the hard sphere system, and its pair potential function uHS(r) is as follows:

The pair correlation function g(r) of the hard sphere system is approximated as follows:

When is sufficiently small, the term in the configuration integral ZN can be rewritten as follows:

The average is estimated by the above pair correlation function as follows:

Here u(1)(r) is the pair potential:

In this manner, the configuration integral can be written as follows:

In addition, is approximated as follows:

Using this configuration integral, one can write the pressure as follows:

The internal energy U is derived from the partition function as follows:

3 Extension of van der Waals EOS

In this section, van der Waals EOS is extended by a cluster expansion of the attractive term in the canonical partition function in the perturbation scheme described in the previous section. There may be monomers, dimers, trimers, etc., in the system according to the temperature and density. The largest cluster in the system is assumed to be an M-mer. The structure and energy of the J-mer cluster is simplified as follows. The J-mer cluster has (J - 1) bonds that start from the central molecule being considered. Each bond has energy -e. The energy of the J-mer cluster is then . The factor 1/2 appears because the present theory is a one-body approximation, and this factor is necessary in order to obtain a reasonable total energy for the system. Therefore, in the present theory, the J-mer cluster has the Boltzmann factor according to the temperature of the system. The geometrical probability of its being found in the system is proportional to the pre-exponential factor . Here, n is the volume per molecule, and n0 is that in the case of the closest packing. The exponent (J - 1)/3 is assumed because the finding probability of the J-mer cluster is proportional to the ratio of the length of the closest packing structure to that of the characteristic length of the volume per molecule of the system. Based on this assumption the partition function is written as follows:

Next, we introduce the molecular partition function of the attractive term qa as follows:

This molecular partition function is compared with that of the original EOS, eβar. For this molecular partition function, the attractive term of the internal energy Ua is obtained as follows:

The contribution of the attractive interaction to pressure is as follows:

The total internal energy and pressure are as follows:

The Helmholtz free energy and entropy are calculated by the standard method from the partition function.
If we express the molecular interaction in the Lennard-Jones form, the parameters in the attractive part are given by the Lennard-Jones parameters. The Lennard-Jones potential is written as follows:

The constant e in Eq. (22) is equal to eLJ in Eq. (27). The closest packing structure in the Lennard-Jones system has the following volume:

This value is used in Eq. (22).
The last parameter in the present model is the effective volume of the hard sphere b in Eq. (22). We adjust this parameter to fit the calculated pressure and internal energy with the observed pressure and internal energy (macroscopic or molecular observations).

4 Molecular Dynamics Simulation

NVT MD simulations [8] on fluid argon were performed at several temperatures to obtain the average pressure p and internal energy U as a function of volume V. The Lennard-Jones potential is adopted as shown in Table 1. The time step is 1 fs, and the standard MD length is 1,000,000 steps. The velocity scaling method was used to control temperature. The number of molecules in the cell is 1,000. A periodic boundary condition is assumed on the cubic basic cell [8]. The assumed initial configuration is the simple cubic lattice, in which only the density varies according to the volume. The Materials Explorer (v4 and v5) program was used for the present simulation [9].
The obtained pVT and UVT relations are shown in Figure 1 and Figure 2. The abscissa shows the density r = N/V. The obtained pressure qualitatively resembles the van der Waals EOS (see Eq. (19)). The pressure is given in units of e/s3 as well as in units of atm [10]. On the other hand, the calculated internal energy depends on temperature and density in a different manner than the van der Waals EOS shown in Eq. (21). This difference is discussed in the next section.

Table 1. Lennard-Jones potential parameters for argon [9].
e/(kcal/mol)(e/k)/Ks/m
0.24831253.43E-10


Figure 1. The pVT relation calculated by MD simulation on fluid argon using the Lennard-Jones model. The ordinate shows pressure p given in Lennard-Jones units e/s3 (Table 1) and in atm [10]. The abscissa shows the number density N/V and the mass density r in units of g/cm3. The constant temperature curves are shown at T = 10, 50, 100, 175, 250, and 500 K.


Figure 2. The UVT relation calculated by MD simulation on fluid argon using the Lennard-Jones model as a function of density. The ordinate shows the internal energy per molecule U/N in Lennard-Jones units e (Table 1) and in units of kJ/mol. The abscissa shows the number density N/V and the mass density r in units of g/cm3. The constant temperature curves are shown at T = 10, 50, 100, 175, 250, and 500 K.

5 pVT and UVT relations in the extended van der Waals EOS

The adjustable parameter b is chosen by trial and error as s3 through a comparison of the extended EOS with the results of MD simulation. The fitted pVT and UVT relations are described in Figure 3 and Figure 4 on the extended van der Waals EOS. The pVT relation in Figure 3 is globally very similar to that in Figure 1. On the whole, the internal energy in Figure 4 also resembles that of the Lennard-Jones system shown in Figure 2. A more direct comparison is shown in Figure 5 and Figure 6. The pVT and UVT relations in Figure 5 and Figure 6 mean that the extended EOS can reproduce, at least qualitatively, the MD simulation results.


Figure 3. The pVT relation calculated by the extended van der Waals EOS. The ordinate is pressure p in Lennard-Jones units e/s3 (Table 1). The abscissa is the number density N/V. The constant temperature curves are shown at T = 10, 50, 100, 175, 225, 275, and 500 K.


Figure 4. The UVT relation calculated by the extended van der Waals EOS as a function of density. The ordinate shows the internal energy per molecule U/N in Lennard-Jones units e (Table 1). The abscissa shows the number density N/V. The constant temperature curves are shown at T = 10, 50, 100, 175, 225, 275, and 500 K.

6 Extended EOS compared with MD simulation

The calculated results are compared with the MD simulations in Figure 5 and Figure 6. In these figures, the present EOS provides reasonable results for pressure and internal energy, although the agreement is only semi-quantitative. Since we have only three parameters in the present EOS, the agreement is satisfactory. Although the function is not so simple as the original van der Waals EOS, the physical meanings of the three parameters are clear in the present case. The construction rule is also unambiguous in the canonical partition function.
The simulation results are usually fitted by EOS with several fitting parameters [11, 12]. The critical point is determined by such EOS [11, 12]. The critical point determined by the present EOS is compared with that shown in Table 2. The obtained result does not differ greatly from the MD simulation results even with only three parameters in the present extended EOS.


Figure 5. The pVT relation calculated by the extended van der Waals EOS compared with that obtained by MD simulation. The ordinate shows pressure p in Lennard-Jones units e/s3 (Table 1). The abscissa shows the number density N/V. The constant temperature curves are shown at T = 50, 200, and 500 K.


Figure 6. The UVT relation calculated by the extended van der Waals EOS compared with that obtained by MD simulation. The ordinate shows the internal energy per molecule in Lennard-Jones units e (Table 1). The abscissa shows the number density N/V. The constant temperature curves are shown at T = 50, 200, and 500 K.

Table 2. The critical point determined by the extended van der Waals EOS is compared with that in the fitted EOS of MD simulations on Lenard-Jones system.
Tc/(e/k)pc/(e/s3)rc/s-3
Nicolas et al [11]1.350.1420.35
Kolafa & Nezbeda [12]1.33960.14050.3108
Present Study1.80.2130.33
Tc/Kpc/atmVc/(cm3/mol)
Nicolas et al [11]16960.069
Kolafa & Nezbeda [12]16759.478
Present Study22590.074

Table 3. The van der Waals coefficients for fluid argon [1].
a/(atm.L2/mol2)b/(L/mol)
1.3453.22E-02

As shown in Figure 7, the extended EOS pressure is different from the MD results at low temperatures. The effect of attractive interaction is not well incorporated in this case. This is due to the simplified cluster structure, which is determined by the short-range term in the molecular interaction. The original van der Waals EOS also yields poor results at low temperatures, as shown in Figure 7. This curve shows that the attractive force is too strong in this EOS. The van der Waals coefficients used in Figure 7 and Figure 8 are shown in Table 3 [1]. A comparison of internal energy values is shown in Figure 8 at low temperatures.


Figure 7. Pressures are compared for the extended van der Waals EOS, the original van der Waals EOS, and the MD results on the Lennard-Jones system as a model of argon at T = 10 K = 0.08 e/k. The abscissa shows the number density N/V.


Figure 8. Average potential energies per molecule Ue/N = U/N - 1.5kT are compared for the extended van der Waals EOS, the original van der Waals EOS, and the MD results on the Lennard-Jones system as a model of argon at T = 10 K = 0.08 e/k. The abscissa shows the number density N/V.

The value -6e shown in Figure 8 of the present study is obtained from the 12 nearest neighbors in the FCC lattice. With respect to the internal energy, the extended van der Waals EOS is better than the original van der Waals EOS, as shown in Figure 8. The shortcomings in the present study (as shown in Figure 7 and Figure 8) may be overcome by the inclusion of the long-range part of the attractive interaction.

7 Conclusion

An extended van der Waals EOS is obtained by a statistical mechanical treatment on the Lennard-Jones system. The new EOS reproduces the pVT and UVT relations calculated by molecular dynamics simulations, at least qualitatively. As the effective Lennard-Jones parameters are known for many molecular systems, the new EOS is useful to estimate the pVT and UVT relations on such systems.

The authors would like to thank the Research Center for Computing and Multimedia Studies, at Hosei University for use of computers.

References

[ 1] P. W. Atkins, Physical Chemistry, 6th edition, Oxford Univ. Press (1998).
[ 2] R. Yamamoto, H. Tanaka, K. Nakanishi, and X. C. Zeng, Chem. Phys. Lett., 231, 401 (1994).
[ 3] R. Yamamoto and K. Nakanishi, Phys. Rev. B, 49, 14958 (1994).
[ 4] R. Yamamoto, O. Kitao, and K. Nakanishi, Mol. Phys., 84, 757 (1995).
[ 5] R. Yamamoto and K. Nakanishi, Phys. Rev. B, 51, 2715 (1995).
[ 6] R. Yamamoto and K. Nakanishi, Mol. Simulation, 16, 119 (1996).
[ 7] D. A. McQuarrie, Statistical Mechanics, Harper Collins (1976).
[ 8] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford Univ. Press, Oxford (1987).
[ 9] http://software.fujitsu.com/jp/materials-explorer/
[10] The relation between the units of pressure is given as follows: 1 atm = 101,325 Pa.
[11] J. J. Nicolas, K. E. Gubbins, W. B. Streett, and D. J. Tildesley, Mol. Phys., 37, 1429 (1979).
[12] J. Kolafa and I. Nezbeda, Fluid Phase Equil., 100, 1 (1994).


Return