Argon force field revisited: a molecular dynamic study

There are many argon force fields that reproduce the liquid-vapor phase diagram, but the surface tension is not obtained with the same force field [1–5]. In figure 1 we can appreciate different force fields; the force field named Vr, developed by Vrabec et al [1], was obtained by considering an extra quadrupolar term to the LJP, this force field was improved to obtain both the liquid-vapor phase diagram and surface tension, adding to last property the effect of three-body interactions, [6], model developed by Goujon et al (Go). The force field developed by White (Wh) is based on a renormalization group theory [2], the phase diagram fits excellently with the experiment except for surface tension. Similar behaviors to the one mentioned above are the force fields developed by Barker et al [3], Rowley et al [4], and Rahman [5], labeled as Ba, Ro, and Ra, respectively.

Zoom In Zoom Out Reset image size

The LJP is employed as the force field for argon, its functional form is

$$\begin{eqnarray}&&U(r)=4\epsilon \left[{\left(\displaystyle \frac{\sigma }{r}\right)}^{12}-{\left(\displaystyle \frac{\sigma }{r}\right)}^{6}\right],\end{eqnarray} \tag{ 1 }$$

where is the depth potential well and r = σ is where the potential change of sign.

The main objective of this work is to improve some argon force fields to reproduce surface tension and liquid-vapor phase diagram, simultaneously, with molecular dynamics by applying a renormalization method to the force field parameters without adding any extra terms to the surface tension like Goujon et al [6] does. The molecular dynamics simulation was carried out using GROMACS software [8] with the Velocity-Verlet algorithm, an NVT ensemble was used to obtain the surface tension and density profile. A V-rescale thermostat was used (τ = 0.5 ps) to keep the temperature constant. Furthermore, a time step of Δt = 0.002 ps and periodic boundary conditions were established. Zero charge and molecular weight of 39.948 were used. All simulations were run 5 ns of production after 5 ns of equilibration. Additionally, the surface tension was calculated from the mechanical definition,

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{1}{2}{L}_{z}\left[{P}_{{zz}}-\displaystyle \frac{1}{2}({P}_{{xx}}+{P}_{{yy}})\right],\end{eqnarray} \tag{ 2 }$$

where P_xx , P_yy , and P_zz are the diagonal elements of the pressure tensor. To calculate liquid-vapor density we considered the plateau density profile in liquid or vapor regions.

Reparameterization is an empirical (fast and easy) method to obtain better properties of modeled compounds by scaling charges, LJP energy and distance parameters, and bond distance to reach the dielectric constant, surface tension, density [9] or micelle radius [10], and self-diffusion constant [11], respectively. A new reparametrization method is applied to argon force field by considering these methodologies. A new methodology is developed as follows:

• (1)  
  The distance LJP parameter, σ, is scaled to obtain the liquid density of the system.
• (2)  
  The energy LJP parameter, , is modified to obtain the surface tension.
• (3)  
  The (1) and (2) steps are repeated to obtain 5% or less of percentage error respect to experiment values.

The temperature to perform reparameterization was kept constant, after obtaining surface tension and density of the liquid; the temperature was varied to obtain the complete liquid-vapor phase diagram and the temperature dependence of the surface tension. To perform the reparameterization, a temperature value of 110K was chosen to ensure that the simulation data fit as best as possible over the entire temperature range. We reached that conclusion by taking the force field from Vrabec et al [1], then reparameterizing the force field at different temperatures, figure 2 shows the reparameterized data at 90 K, 110 K, and 130 K. The data obtained at a temperature of 110 K are closer to the experimental liquid density than other temperatures, so the Wh, Ba, Ro, Ra, and Go force fields were reparameterized at that temperature.

Zoom In Zoom Out Reset image size

The adequate form to calculate surface tension is to analyze the system to avoid any finite size effects , we analyze different interfacial areas and cut radius (r_cut) as illustrated in figure 3, where the surface tension values for an area of 4 nm × 4 nm are not stable, these values increase as the r_cut increases, for the area of 6 nm × 6 nm they are stable from r_cut = 2.8 nm, and for 8 nm × 8 nm stability is reached from 2.6 nm. Then, molecular simulation was performed at r_cut = 2.8 nm with an interfacial area of 6 nm × 6 nm and 30 nm in z-direction containing N = 6750 argon atoms. This analysis and the scaling process improved the value of the surface tension.

Zoom In Zoom Out Reset image size

The reparameterization procedure was applied to obtain the surface tension and the liquid-vapor phase diagram at a fixed temperature. Results shown that original force fields have low percentage error (0.19–0.74) for liquid density but Ba model, see table 1, in the case of surface tension the percentage errors are of 15–56. After reparameterizing all argon force fields was obtained low percentage errors; the maximum error was detected in the Ro force field with 0.19 % for liquid density and for Wh force field with 6.74 % for surface tension, which are low respect to the original ones. The Ba force field is the best reparameterized model due to the combined errors of both properties is lower than any other force field. These results are independent of how the original parameters were obtained and depend on the ability to adjust them to obtain the experimental properties.

The temperature dependence of these properties is shown in figure 4, it is observed that all force fields are close to the experimental surface tension as a function of temperature, the Wh force field is closer than the other force fields [2] for experimental gas density, see inset in figure 4. The phase diagram has a similar behavior for almost all force fields; the liquid density branches are well described, the surface tension and phase diagram are well reproduced by all models after reparameterization.

Zoom In Zoom Out Reset image size

The Lennard-Jones parameters from literature and the reparameterized ones (obtained in this work) are summarized in table 2. The new or reparametrized parameters are close to each other, but far from the original ones that only reproduce the liquid-vapor phase diagram. In figure 5, original and new force fields parameters represent separate regions in the -σ phase space, there is an original force field (Ba) [3] which is not part of any region, this force field did not reproduce both properties, this phase space indicates that, if we want to reproduce the liquid-vapor phase diagram, we have to choose a value of from 0.971 07 kJmol⁻¹ to 0.997 74 kJ mol⁻¹ and an σ value between 0.336 05 nm to 0.340 50 nm, if we also want to reproduce the surface tension then we have to choose a value of from 0.941 91 kJ mol⁻¹ to 0.946 37 kJ mol⁻¹ and a value of σ from 0.336 05 nm to 0.337 13 nm.

Zoom In Zoom Out Reset image size

Table 2. Left table literature and right table reparameterized and σ LJP parameters for argon compound in kJ mol⁻¹ and nm, respectively.

Force field		σ	Force field		σ
	kJ mol−1	nm		kJ mol−1	nm
Vr [6]	0.97107	0.33952	Vr rep	0.94191	0.33646
Wh [2]	0.98275	0.34050	Wh rep	0.94639	0.33713
Ba [3]	1.18144	0.33605	Ba rep	0.93570	0.33605
Ro [4]	0.99607	0.34050	Ro rep	0.93631	0.33641
Ra [5]	0.99774	0.34000	Ra rep	0.94071	0.33677

The empirical process used in this work could be improved by obtaining the LJP parameters by a renormalization group theory, as in the Wh model, which was obtained [2, 12] by applying it not only to liquid-vapor phase diagram prediction but also to the surface tension. The renormalization procedure considers the contribution of the repulsive potential for high temperatures and for low temperatures the part of the attractive potential, considering a renormalization procedure to the Helmholtz free energy, this expression and a correct choice of parameters (, σ) provide a complete prediction of the liquid-vapor phase diagram.

In conclusion, we study some argon force fields by obtaining surface tension and liquid-vapor phase diagram, a new methodology to reparameterize force fields was applied to improve them at a fixed temperature, obtaining that all force fields converge in a region in the -σ phase space, which is different from the original ones, which reproduces only the liquid-vapor phase diagram. All reparameterized force fields reproduce well the experimental surface tension temperature dependence and liquid-vapor phase diagram. Furthermore, the Ba model has low percentage errors for both reparameterized properties at 110K and the Wh force field has the better behavior for gas density. This empirical procedure could be justified by applying a renormalization method in both LJP parameters.

Acknowledgments

Data availability statement