Mean spherical approximation for the Lennard-Jones-like two Yukawa model: Comparison against Monte Carlo data

Monte Carlo simulation studies are performed for the Lennard-Jones like two Yukawa (LJ2Y) potential to show how properties of this model fluid depend on the replacement of the soft repulsion by the hard-core repulsion. Different distances for the positioning of hard core have been explored. We have found, that for temperatures that are slightly lower and slightly higher of the critical point temperature for the Lennard-Jones fluid, placing the hard core at distances that are shorter than zero-potential energy is well justified by thermodynamic properties that are practically the same as in original LJ2Y model without hard core. However, going to extreme conditions with the high temperature one should be careful since presence of the hard core provokes changes in the properties of the system. The later is extremely important when the mean spherical approximation (MSA) theory is applied to treat the Lennard-Jones-like fluid.


Introduction
The so-called simple fluid models cannot be always applied to the study of realistic systems. However, this type of models is very important in basic science while studying fundamental problems in the liquid state theory. One of these problems refers to the role played by repulsive and attractive forces. While comprehension of the repulsive forces is due to the progress reached in the theoretical and computer modeling of hard sphere fluids, the simplest nontrivial model that makes it possible to study the liquid phase of a matter as well as the vapor-liquid coexistence seems to be the one that consists of a hard-core plus attractive Yukawa (HCAY) potential where is the depth of the potential energy well, R is a hard-core diameter, κ −1 is a measure of the range of the attractive tail. An important feature due to which the HCAY model is intensively employed in the studies of simple fluids is that approximate analytical and semianalytical solutions for this model are available owing to the pioneering paper by Waisman [1] on the analytic solution of the Ornstein-Zernike equation using the mean-spherical approximation (MSA). These were Henderson et al. [2] who by using this MSA solution suggested that having R = σ and κR = 1.8, the HCAY fluid is qualitatively similar to argon for the densities and temperatures of the liquid in equilibrium with its vapour and with the potential parameters approximately the same as the Lennard-Jones potential (i.e. /k 120 K, σ ≈ 3.4 Å). Figure 1  Recently, Kadiri et al. [3] have extended the mapping of the LJ system into the HCAY system over a wide domain of a phase diagram. Specifically, at every density and temperature the thermodynamics of the LJ fluid was reproduced from the equation of state of the HCAY fluid with suitable values of its two parameters R and κ. An advantage of such an approach is that it provides analytical equations for the thermodynamics of the LJ system with only two parameters that are of precision comparable to that of the simulation data. It can be compared to the fully empirical equations of state by Johnson et al. [4] that uses 33 parameters. However, one can see from figure 1 that the shapes of two potential profiles, LJ and 1Y, are rather different and one cannot be sure that the thermodynamics predicted by two models will be always the same. Thus, this approach should be used with great care while exploring the thermodynamic states outside of the mapping domain.
More possibilities are offered by the MSA solution that has been obtained by Blum and Hoye [5] for a linear superposition of the Yukawa tails i exp(−κ i r)/r. This was quite an important step since superposition of the attractive and repulsive Yukawa tails makes it possible to mimic practically any potential profile. It is natural that this MSA solution has been applied to Lennard-Jones fluid represented by the hard core and sum of two Yukawas (HC2Y) [6,7], where again R = σ determines the separation of zero energy, u(σ) = 0, and is, as usual, the depth of the potential minimum. The detailed course of the 2Y potential (with a specific set of the decay parameters κ i which will be discussed in the following section 2) is shown in figure 1 where it is also compared with the Lennard-Jones (LJ) potential. Indeed, two Yukawa terms allow us to significantly improve the mimicking of LJ potential at short and intermediate distances. This is a rather important achievement since in this case the differences in the properties of these two model systems can be attributed to the long-ranged asymptote of the LJ interaction only. Another important issue is that in general, the LJ2Y potential (like the parent LJ potential) does not require the presence of hard-core, e.g., for computer simulation studies; the hard-core is necessary if one intends to apply the MSA theory to the treatment of the problem. However, even in this case placing the hard-core at R = σ (as it usually is assumed) seems not to be always appropriate and justified since it may modify the properties of the original model. To clarify this issue, we will report here the Monte Carlo simulation data for the original LJ2Y potential without hard-core as well as for the LJ2Y potential with a few different choices of the hard-core diameter R in the range 0.8 < R/σ < 1. One more question that we are seeking to answer in this study concerns the performance of MSA theory for different choices of the hard-core diameters 0.8 < R/σ < 1. The MSA studies reported so far for the LJ2Y model are dealing mainly with the case R = σ. The structure of this paper is as follows. In section 2 we outline the way the parameters of the LJ2Y potential have been chosen and present a brief description of the Monte Carlo computer simulations and MSA computations that we have carried out for the LJ2Y model. The results are collected and discussed in section 3 while section 4 contains conclusions.

The potential models and computational details
Originally, the LJ2Y potential is a model made up of two Yukawa tails without any hard-core, where r m is the position of the potential energy well, 1 > 0 and 2 > 0 are the strengths of the repulsive and attractive contributions, respectively, while κ −1 1 and κ −1 2 are the measures of the range of the corresponding tails. Superposition of two Yukawa tails in the form of equation (3) to mimic the LJ potential has been already used by other authors [8][9][10]. However, all previous studies have been concerned with the medium and long distances without paying much attention to short separations governed by repulsive forces.
To map the properties of the LJ fluid onto those of a LJ2Y fluid, the above four parameters of potential function (3) must be determined. Focussing first on the course of the LJ function itself, the following three conditions seem evident: the coincidence of the location and depth of the potential minimum of the LJ and LJ2Y functions, and the location of the potential zero, Then there remains to impose one more condition on the 2Y potential parameters to complete the set of equations. Some authors attempt to follow the (12-6) LJ curve to certain intermediate distance r * by equalizing the integrals from two potential functions in the interval from r min to r * . By contrast, in this study we will set the condition that both potentials attain the same value at the location of the inflection point r = r inf = 1.244455σ that follows immediately after the position of the potential well of the LJ potential, i.e., Proceeding in this way we obtain a 2Y potential that very accurately reproduces the LJ potential profile at distances up to inflection point, but showing some discrepancies after inflection point due to an exponential decay at large distances. The resulting parameters are as follows: 1 = 1954325.046σ, 2 = 50.26984765σ, κ 1 σ = 13.66462 and κ 2 σ = 3.10147. These values are consistent with those in equation (2) and in figure 1 (b). Figure 2 shows the vapour-liquid phase diagram for the LJ2Y potential model. To determine the vapour-liquid envelope we used the common Gibbs ensemble with the total number of particles N = 512 and applied the long-range correction to truncate the potential at r c = 3 N/ρ * . For the purpose of comparison we are presenting computer simulation data for the liquid-vapor coexistence in a parent LJ fluid. The differences that are observed in figure 2 concern the lowering of the critical point temperature and are caused by the differences between two potential models, LJ and LJ2Y, at distances r larger than the position of inflection point r inf . These differences are consistent with what has been already learned from the computer simulation studies of the HCAY model, namely, the shortening of the range of attraction brings about a decrease of the critical point temperature.
The corresponding MC data are collected in table 1. To discuss the role played by repulsive forces of the hard-core origin we employed both Monte Carlo (MC) computer simulations and MSA approaches. First, we performed MC computer simulations for the original LJ2Y model given by equation (3) for a set temperature and density conditions that are representative for LJ fluid model. Then similar MC computer simulations at the same thermodynamic states have been repeated for the LJ2Y model with an artificially embedded hard-core at different interparticle separations smaller than LJ diameter σ.
In all these studies, the standard Monte Carlo simulations for the LJ2Y model fluids were carried out in an NVT ensemble with N = 512 particles using the largest possible cutoff for each isotherm and with appropriate long range corrections applied [12]. In addition to the common thermodynamic properties, the internal energy and pressure (evaluated by the virtual volume change), we also determined, for all thermodynamic state conditions considered, the excess chemical 33005-4 potential using the standard Widom's particle insertion method [12]. Our MSA study is based on the general solution obtained by Blum and Hoye [5] for a mixture of hard spheres interaction via pair potential that is the sum of Yukawa tails with different decay parameters. The numerical solution of the resulting set of nonlinear algebraic equation was obtained using a relatively simple and effective iterative procedure proposed recently by Kalyuzhnyi and Cummings [13].
Firstly, we are expecting that the placement of a hard-core at different distances may modify the properties set up for the initial LJ2Y potential function given by equation (3). This issue can be resolved by analyzing the MC simulation data. Secondly, the MSA performance could depend on the position of the hard-core, and this can be concluded by comparing with computer simulation data.

Results and discussions
The studies that we are reporting here have been performed for three reduced temperatures, T * = 1.25, 1.45 and 4.85. These three temperatures reflect three important temperature conditions of the parent LJ fluid, namely, subcritical, supercritical and extreme. This can be illustrated by the phase diagram shown in figure 2 with a reminder that critical point temperature of the LJ fluid is around T * c,LJ ≈ 1.35 [11]. The distances that have been explored for the position of hard-core diameter are: R = σ, 0.95σ, 0.9σ, 0.85σ and 0.8σ. The values of the internal energy U * , pressure P * and compressibility factor Z = P/ρkT for each hard-core position and at three temperature conditions are collected in tables from 3 to 7. Additionally, the density dependencies of the pressure P * and internal energy U * are illustrated in figures 3 and 4, respectively. Typical radial distribution functions are presented in figures 5 and 6. Tables and figures show both computer simulation data and the results obtained from the MSA theory. For the sake of comparison between different parts, figures 3 and 4 are using the same scale which is not always convenient for distinguishing the data within the same part of the figure. For these purposes we recommend to look for the data in tables 2-7.
For the lowest temperature, T * = 1.25, which is slightly below the critical temperature for the LJ fluid one can see (upper parts in figures 3 and 4) that MC simulation data indicate very weak dependence of the thermodynamic data on the replacement of soft repulsion by stiff hard-core repulsion in the entire range of distances r σ. This is more evident for pressure, when only the limiting case R = σ (filled triangles in figures from 3 to 6) indicates the tendency to be separated from the rest of data; in the case of internal energy, similar tendency is seen for the case R = 0.95σ. Quite similar dependencies on the hard-core position are shown by MC computer simulation data for the higher temperature, T * = 1.45 (see the middle parts in figures 3 and 4). Indeed, we can see that the case R = σ is rather special for these temperature conditions. Interestingly, this is also true for the MSA theory (thick solid lines in figures from 3 to 6). We also note that being compared with MC data for the case of a hard-core placed at zero-potential energy separation distance R = σ, the MSA theory performs quite well for the thermodynamics of the system. The MSA curves calculated with the hard-core placed at shorter distances, R = 0.95σ, for both pressure and energy (dashed lines in figures 3 and 4) are clearly separated from the MSA results obtained with the hard-core at R = σ and are very close to the MC computer simulation data. At the highest temperature, T * = 4.85, (see the bottom parts in figures 3 and 4) that we are referring to as the extremely high for the LJ fluid (more than three times higher than its critical point temperature), simulation data show that fixing a hard-core at distances r 0.9σ produces, for all the calculated thermodynamic properties, the values that are very close to those for the initial LJ2Y model without a hard-core. At the same time, hard-core R = 0.9σ seems to be the smallest one when the MSA theory performs well. Although this is not the case of the pressure  (bottom part in figure 3), but this is already seen from the behavior of the internal energy (bottom part in figure 4) and, especially, from the comparison with MC data for the radial distribution function (figures 5 and 6) when extra oscillations around the first minima start to appear.

Conclusions
In this study we report the Monte Carlo (MC) simulation studies of the thermodynamics and radial distribution functions of the Lennard-Jones-like two Yukawa (LJ2Y) fluid. This fluid model has been defined to nearly exactly mimic the pair interaction at short and intermediate distances in the Lennard-Jones (12-6) fluid. In particular, the parameters of the LJ2Y have been defined from the condition that zero-potential energy distance and the slope of the potential profile at the zero-potential distance are the same. Moreover, the magnitude and position of the potential minima in both models coincide. The only differences that two potential models experience are attributed to large distances. Exactly these differences have caused quite pronounced differences in the vapour-liquid coexistence occurring in two models. However, the purpose of the present study is not to describe the LJ fluid by means of the two Yukawa model. Our goal is to find how the mean-spherical approximation (MSA) performs when it is applied to the fluid model whose soft repulsion is the same or quite similar to that of the LJ fluid. The task we set forth is twofold. First,   the application of the MSA theory requires cutting-off the soft repulsion at some interparticle separation and placing the hard-core R at this distance. However, only such a modification of the initial potential may intervene and change the properties of the system. Without this knowledge it 33005-8   will be hard to judge about the MSA performance. To fulfill this task, we performed MC computer simulations for both the initial LJ2Y model defined according to equation (3) and for hard-core modified potential (2) -the form that is usually utilized by the MSA theory. The positions of the hard-core that have been implemented in computer simulations include R = σ; 0.95σ; 0.9σ; 0.85σ and 0.8σ. Having compared this with the data obtained in computer simulations with initial 33005-9 LJ2Y potential, i.e., without hard-core, we conclude that the insertion of the hard-core into a LJlike models is extremely sensitive to the temperature conditions. Namely, at normal temperature conditions of the order and around the critical point temperature, an insertion of the hard-core at separations r 0.95σ practically does not change the thermodynamics of the initial system, while the case R = σ already shows the tendency to exhibit slightly different thermodynamic properties. The conclusions change when we explored the high temperature conditions, namely, T * = 4.85. At these temperature conditions, the presence of the hard-core at distances as short as r 0.85σ may be considered not to modify the initial system. All other cases, R = σ; 0.95σ and 0.9σ should be treated as separate model fluids that are different from the initial LJ2Y fluid model. These observations are valid for both thermodynamics and radial distribution functions.
Having taken these considerations into account we were able to fulfill the second part of our task, namely, to judge the MSA performance for the LJ2Y fluid model. First of all, we found that MSA correctly reflects the temperature conditions when the case of a hard-core R = σ starts to deviate from the other hard-core modified models. Secondly, the MSA performs rather satisfactorily for this case R = σ at all considered temperature conditions. The following conclusion seems to be rather important and reads: at normal temperature conditions (of the order and around the critical point temperature for the LJ fluid), an optimal position of the hard-core R for the MSA to be a reasonable theoretical approach is R = 0.95σ. This means that placing the hard-core at R = 0.95σ and using the MSA theory to describe the LJ-like two Yukawa model will give you the results that are of reasonable accuracy and are very close to those of the initial model. By contrast, using the hard-core for MSA description as small as R 0.9σ will result in incorrect values for initial LJ2Y fluid model and even for the LJ-like model with an embedded hard-core.
The above conclusion is valid when you are seeking the possibility to use the MSA theory at high temperatures, T * = 4.85. The only thing that changes in this case is the limiting value of the hard-core. Namely, at this high temperature you may use the hard-core position as small as R = 0.9σ and still continue to apply the MSA theory and obtain reasonably accurate description of the LJ2Y fluid model. On the contrary, using the hard-core at R = 0.85σ and 0.8σ leads to an increased inaccuracy in both thermodynamics and the radial distribution functions. We also note that at this high temperature, the hard-core modified LJ2Y models with R = σ and R = 0.95σ significantly differ from the initial LJ2Y model without hard-core and, very importantly, the MSA theory performs quite well for both these cases.