Two-dimensional core-softened model with water like properties. Study by thermodynamic perturbation theory

Thermodynamic properties of the particles interacting through smooth version of Stell-Hemmer interaction were studied using Wertheim's thermodynamic perturbation theory. The temperature dependence of molar volume, heat capacity, isothermal compressibility and thermal expansion coefficient at constant pressure for different number of bonding sites on particle were evaluated. The model showed water-like anomalies for all evaluated quantities, but thermodynamic perturbation theory does not properly predict the dependence of these properties at a fixed number of bonding points.


Introduction
They were Stell and Hemmer who first proposed core-softened potentials in 1970 [1]. In their early work, they stressed that negative curvature in interaction potential might lead to a second critical point in addition to a standard liquid-gas critical point. In different works [2,3] it has been shown that coresoftened potentials and similar shouldered potentials can reproduce various fluid anomalies that are typical of water and other substances with angular dependent interactions, such as silica [3], silicon [4], BeF 2 [5]. Core-softened potentials were also used to study single-component liquid metal systems [6][7][8][9][10] and as solvent for studying ions [11]. Poole et al. [12] proposed liquid-liquid phase transition as an explanation for anomalous properties of water. After that there was an increased interest to the studies of these liquid-liquid phase transitions. Franseze et al. [13] suggested that the liquid-liquid phase transition and its critical point might be caused by the potential with two characteristic distances (hard core and soft core). In their work, they reported the existence of the low-density liquid phase and the high-density liquid phase obtained for 3D model using molecular dynamics (MD) simulations. On the other hand, 2D MD produced only a density anomaly but no liquid-liquid phase transition [14,15]. Scala et al. [16] carried out MD simulations of 2D discrete and smoothed version of potential to study liquid anomalies. These studies were continued by Buldyrev et al. [17] to explore liquid-liquid phase transition for 2D and 3D version of potentials and by Almudallal et al. [18]. They both produced phase diagrams for a discrete version of potential with liquid anomalies, and no liquid-liquid critical point in stable liquid region was obtained.
Our aim here is to apply Wertheim's thermodynamic perturbation theory (TPT) to capture the physics of the model of 2D molecules interacting by Stell-Hemmer potential. In recent years, a theory has been developed for fluids comprised of molecules that associate into dimers and higher clusters due to the presence of highly directional attractive forces [19][20][21]. In the present work, we apply the thermodynamic perturbation theory (TPT) [19][20][21][22] to central symmetric attractive potential.

Model
The smooth version of the core-softened potential proposed by Scala et al. [16] is used in this work.
The interaction potential U (r ) is a sum of a Gaussian well and the Lennard-Jones (LJ) part of the potential U (r ) = U LJ (r ) +U a (r ), ε is here the well-depth and σ is the distance, where LJ part of the potential is zero. The Gaussian part of the interaction is as follows: This part of potential is stronger than the LJ part and is the reason that the particles make strong association. We also refer to this part of potential as association potential. We apply the units and values of model parameters as used before by Scala et al. [16] as ε = 1.0, σ = 1.0, λ = 1.7, a = 25.0, r 0 = 1.5σ. Figure 1 shows the shape of the smooth version of the core-softened potential used here.

Monte Carlo simulation details
We performed Monte Carlo simulations in the isothermal-isobaric (NpT) ensemble to obtain thermodynamic properties of the model. At each step, the displacements in the x, y coordinates were chosen randomly. We used periodic boundary conditions and the minimum image convention to mimic an infinite system of particles. The starting configurations were selected at random. Every 10 moves of particles an attempt is made to scale the dimensions of the box and all of its component particles in order to hold the pressure constant. 5×10 4 moves per particle were needed to equilibrate the system. The statistics were gathered over the next 1×10 6 moves to obtain well converged results. All simulations were performed with N =200 or N =400 molecules. The maximum change of dimensions of the box was calibrated during equilibration simulations. The physical properties of the system such as enthalpy and volume were calculated as the statistical averages of these quantities over the course of simulations [24]. The heat capacity, C p , the isothermal compressibility, κ, and the thermal expansion coefficient, α are computed from the fluctuations [25] of enthalpy, H , and volume, V .

43605-2
Two-dimensional core-softened model T is temperature of the system and N number of particles.

Thermodynamic perturbation theory
The Helmholtz free energy of the system is the key quantity of the thermodynamic perturbation theory [19,20]. In case of the model studied in this work, this quantity is the sum of two terms N is the number of molecules, T is temperature and k B is Boltzmann's constant. The Helmholtz free energy of Lennard-Jones system, A LJ , is calculated using the Barker-Henderson perturbation theory [26] A LJ A HD is the hard-disk contribution to the Helmholtz free energy, g HD (r, η) is pair correlation function for hard disks at packing fraction η = 1 4 πd 2 ρ and ρ is the number density of molecules. d is the hard-disk diameter calculated using Barker-Henderson approximation as We used the procedure by Scalise et al. [27] to calculate the HD term of the Helmholtz free energy For g HD (r ), the expression of Gonzalez et al. [28] was used.
The association contribution to Helmholtz free energy, A a , was calculated by [19,20,29] where x is the fraction of molecules not bonded at particular interaction site and is obtained from the mass-action law [19,20] in the form ρ is the total number density. Finally, ∆ is defined by [19,20,29] ∆ = 2π g LJ (r ) f a (r )r dr .  Wertheim's theory is used. We made approximation that each particle can have N a association points in the center of a particle interacting with associating potential. We used a different number of interaction sites, from 1 to 6, the last being a coordination number of particle in a perfect hexagonal crystal. Once the Helmholtz free energy is known, other thermodynamic quantities may be calculated from standard thermodynamic relations [26] (4.12)

Results
All the results are given in reduced units; the excess internal energy and temperature are normalized to the LJ interaction parameter ε (E * = E /ε, T * = k B T /ε) and all the distances are scaled to the characteristic length σ (r * = r /σ). In figure 2, we compare the molar volume or volume per particle, V * /N , obtained from the Monte Carlo simulations, with the results of the thermodynamic perturbation theory for a different number of bonding sites on particle (N a = 1−6). The calculations were performed at a reduced pressure of P * = 0.75.
We find out that the TPT does not properly capture the results for simulations as well as it does not predict the maxima in density or minima in molar volume. In order to obtain an agreement, N a should be dynamically varied. N a should be varied with temperature and density in order to get an agreement between theoretical and simulation results. We can see from figures that we have good agreement of TPT with simulation for N a = 2 at high temperature T * = 2.0, then 3 at T * = 1, etc.

43605-4
The remaining figures show the temperature dependencies of the other thermodynamic quantities of interest: the isothermal compressibility, κ * T (figure 3), the thermal expansion coefficient, α * (figure 4), the heat capacity, C * p (figure 5), and the excess chemical potential µ * ex ( figure 6). The TPT for a fixed number of associating points is not in agreement with the Monte Carlo simulation data for all quantities. From the results we can also see that the model behaves in a way that the number of associating points is not fixed, but it changes at a constant pressure with temperature. We calculated free energy as a function of the number of boding sites, N a , and the free energy decreases within the whole range.

Conclusion
The thermodynamic perturbation theory was used to study the thermodynamics of the particles interacting through a smooth version of Stell-Hemmer interaction. The results for the molar volume, the isothermal compressibility, the thermal expansion coefficient, the heat capacity and the excess chemical potential obtained by the TPT theory for a fixed number of bonding sites on the particles are not in agreement with the computer simulation results for all the parameters studied. This is caused by the fact that the interaction is spherical symmetric, and for TPT we have to make approximations. It is crucial to 43605-5 change the spherical symmetric potential into a directional one and to have a different number of interaction points with directional forces. We cannot obtain correct thermodynamic properties with a fixed number of association points. Proper thermodynamics could be obtained if the number of association points varied with temperature and pressure.