Transport properties of the square-well fluid from molecular dynamics simulation

In this work, we have calculated self-diffusion and shear viscosity, two of the most important transport properties, of the spherical square-well (SW) fluid interacting with potential range $\lambda = 1.5 \, \sigma$. To this end, we have used a combination of molecular dynamics simulation and the continuous version of the square-well (CSW) intermolecular potential recently proposed by Zer\'on et al. [Mol. Phys. 116, 3355 (2018)]. In addition to that, we have also determined a number of equilibrium properties, including internal energy, compressibility factor, radial distribution function, and coordination number. All properties are evaluated in a wide range of temperatures and densities, including subcritical and supercritical thermodynamic conditions. Results obtained in this work show an excellent agreement with available data reported in the literature and demonstrate that the CSW intermolecular potential can be used in molecular dynamics simulations to emulate SW transport properties with confidence.


I. INTRODUCTION
The spherical square-well (SW) force field is one of the simplest intermolecular potentials that incorporates both repulsive and attractive interactions.2][3][4][5][6][7][8] Due to its simple mathematical form, SW fluids have been extensively investigated using perturbation theory [9][10][11][12][13][14][15][16][17][18][19][20][21][22] and molecular simulation using the well-known Monte Carlo (MC) technique. 19, Struural and equilibrium properties have been routinely evaluated using the aforementioned theoretical and numerical procedures, including radial distribution function, coordination number, internal energy, compressibility factor, coexistence densities, vapor pressure, interfacial tension, high-temperature expansion coefficients of free energies, and chemical potential, among many others.Unfortunately, there is a limited number of publications in the literature focusing on the determination of transport properties in systems that interact through discontinuous potentials from computer simulation, and the SW fluid is not an exception.][47][48][49][50] The limited number of studies addressing transport properties of fluids interacting through discontinuous intermolecular potentials may be attributed to two primary factors.Firstly, the calculation of transport properties is unfeasible through Monte Carlo (MC) simulation, as it requires the continuous tracking of all particles' trajectories over time.Secondly, investigating discontinuous potentials in molecular dynamics (MD) simulations using conventional methods is not feasible, as the continuous application of forces is indispensable for solving Newton's equations.
To address this issue, authors have previously suggested two distinct approaches.The initial option involves modifying the molecular dynamics algorithm itself.A viable choice is to adopt the collision scheme, replacing the solution of Newton's equations.In this scheme, energy and momentum are conserved quantities both before and after the collision. 51,52Through this algorithm, direct particle contact is circumvented, enabling the computation of transport properties.This approach is exemplified in the works of Michels and Trappernier.Unfortunately, the adoption of this technique is not widespread.
4][55][56][57][58] Examples of this approach are the soft version of the hard sphere 55 and the continuous version of the SW intermolecular potential (CSW) recently proposed by Zerón et al. 56 In this investigation, we adopt the latter approach.Consequently, the principal aim of this study is to determine selfdiffusion, shear viscosity, and diverse structural and thermodynamic properties of the SW fluid through the integration of molecular dynamics simulations with the CSW intermolecular potential.
The organization of the paper is as follow: In Sec.II, we describe the methodology used in this work to estimate transport properties, particularly, self-diffusion coefficients and shearviscosity using standard equilibrium MD simulations.Specific details related to the MD implementation can be found in Sec.III.The results obtained as well as their discussion, are described in Sec.IV.Finally, conclusions are presented in Sec.V. Numerical data of this work has been included as Supplementary Material.

II. TRANSPORT PROPERTIES BY COMPUTER SIMULATIONS
Transport phenomena involves the study of a number of properties in systems in which mass, energy, charge, momentum or angular momentum is exchanged.MD simulations employ two frameworks to calculate transport properties: equilibrium and non-equilibrium.In both scenarios, established formalisms and equations exist for computation.For instance, equilibrium MD simulations commonly utilize the Green-Kubo expression and the Einstein formula.
When considering the equilibrium state of the system, it becomes possible to express a given transport property in terms of an appropriate correlation function of a dynamic magnitude.In the cases of self-diffusion and shear viscosity, the dynamic properties are linked to the positions of the particles constituting the system and the pressure tensor, respectively.Both are defined in the subsequent sections.We recommend the reader the manuscript of Maggin et al. 59 for a detailed account on how to calculate self-diffusivity and viscosity from equilibrium MD.

A. Self-diffusion coefficient
The most common way for determining the self-diffusion coefficient, D, involves the Einstein approach.In this approach, the dynamic magnitude is directly related with the positions of the particles, as functions of time, and D is proportional to the slope of the mean-square displacement (MSD) at long times: where r i (t) is the position of the i-particle at time t.Note that the ensemble average in Eq. ( 1) is taken for all N particles in the simulation to improve statistical accuracy of this property. 59,60

B. Shear viscosity
The shear-viscosity, η, can be calculated with the wellknown Green-Kubo formula: where V is the volume of the system, k B the Boltzmann constant, T the temperature, P αβ (t 0 )P αβ (t 0 + t) t 0 is the autocorrelation function of the pressure tensor, which measure the rate at which the momentum in direction β is transferred along the α direction.The components of the pressure tensor, P αβ , can be expressed as: The first term is the kinetic contribution, being m the molecular mass, v α i the α component of the velocity of the particle i, v i .The second term is the intermolecular potential contribution, where r α i j is the α component of the intermolecular vector r i j between particles i and j, and f β i j is the β component of the vector force of j over i particle, f i j .In homogeneous isotropic systems P αβ is symmetric and η can be enhance averaging multiple terms from the stress tensor.

III. SIMULATION DETAILS
All simulations have been performed in the NVT or canonical ensemble using the well-known GROMACS package (double precision). 61,62To calculate all properties, we are using argon fluids interacting with the continuous version of the SW potential proposed by Zerón et al. 56 The reduced intermolecular potential, v CSW (x) = V CSW (x)/ε, as a function of the reduced distance, x = r/σ , is written as: Here, the hard sphere diameter value is σ = 0.3405 nm, the well energy ε = 0.996078 kJ/mol, and the interaction range λ = 1.5 σ .The softness parameters, n = 2500 and m = 20000, are taken from the original work. 56The CSW potential and its derivative are input into GROMACS through a file, where numerical values are tabulated at intervals of 0.0001 nm.
We have performed simulations using two system sizes, 108 and 1000 molecules.The first election allows to compare our results obtained using the CSW intermolecular potential with simulation data taken from the literature for systems formed from the same number of molecules but interacting with the original SW potential.The second election comprises a significantly larger system.Comparison between results obtained using both system sizes will allow to analyze the possible existence of system-size effects on the transport properties calculated in this work.In both cases, the initial simulation boxes are prepared placing spherical molecules interactions via the CSW at random positions in a cubic simulation box with periodic boundary conditions applied in all directions.The studied densities range from 160 to 1450 kg/m 3 .
In all cases, we firstly run all simulations 2 ns to stabilize the system, followed by 20 ns runs to obtain averaged properties.The leap-frog integrator 63 is employed to solve the dynamic equations with a time step of 0.1 fs.Temperature is controlled using the Nosé-Hoover thermostat 64,65 with a coupling constant of τ = 2 ps in the range of temperatures studied, 80 − 600 K.The positions of all particles are saved each 50 steps and pressure components each 5 steps.
Since we are using the CSW version of the potential, GRO-MACS treats it as a standard intermolecular continuous potential.Consequently, a cut-off distance must be provided.In our calculations, U CSW (r c ) ≈ 10 −15 kJ/mol when λ = 1.5 σ and the cut-off distance is set equal to r c = 0.512 nm.Since the potential value is practically equal to zero, we do not apply long-range corrections during simulations.

IV. RESULTS
The ability of the CSW intermolecular potential in reproducing equilibrium properties of spherical SW fluids in bulk vapor or liquid phases, as well as at vapor-liquid coexistence conditions, has been addressed in the original work by Zerón et al. 56 That study has demonstrated excellent agreement between predictions obtained from MD simulations using the CSW potential and those from MC using the original SW force field.The comparison includes a number of properties, encompassing internal energies, compressibility factors, radial distribution functions, coexistence densities, vapor pressures, and vapor-liquid interfacial tensions.However, transport properties have not been previously explored using this continuous version of the intermolecular potential.
To demonstrate the suitability of the CSW intermolecular potential for standard MD simulations in accurately predicting the transport properties of the SW fluid, we perform simulations to determine the self-diffusion coefficient and shear viscosity of particles interacting via the CSW potential, with λ = 1.5 σ , under subcritical and supercritical conditions.The thermodynamic states considered in this work are illustrated in Fig. 1.Note that the majority of the chosen states align with those previously studied by Michels and Trappeniers. 46,47,50irstly, we compute the self-diffusion coefficients and shear viscosity for a system formed by108 particles in order to compare the MD results obtained using the CSW potential with values found by Michels and Trappeniers (using the standard SW potential).It is noteworthy that we employ the same system size as Michels and Trappeniers to avoid differences due to finite-size effects.Subsequently, we conduct simulations under identical thermodynamic conditions, but this time employing systems with 1000 molecules.As we have already mentioned, this allowed us to analyze the impact of the number of molecules on both transport properties and equilibrium and structural properties.
To quantify the ability of the CSW intermolecular potential in predicting transport, equilibrium, and structural properties of the standard SW potential, we also calculate the average absolute relative deviation (AAD%): and the discontinuous red line separates the subcritical and supercritical regions of the phase diagram at the critical temperature.
Here i runs for all the simulation values used in the calculation and n t is the total number of values.

A. Self-diffusion coefficient
We have performed MD-NVT simulations for systems formed by 108 molecules interacting through the CSW intermolecular potential with range λ = 1.5 σ .Particularly, we have simulated states at densities from ρ * = 0.025 to 0.86 and in the range of temperatures T * = 0.667 − 5.For each point studied, 5 independent simulations of 20 ns are running to determine the average value and to estimate the uncertainty associated to the self-diffusivity.The numerical value has been calculated using the Einstein relation of Eq. (1), i.e., we take one sixth of the slope of the linear fit applied to the 10 − 90% curve of the MSD in the diffusivity zone.
The final values of the self-diffusion coefficient of spherical SW molecules obtained from MD-NVT simulations using the CSW intermolecular potential are shown in Fig. 2. We have also included the values obtained by Michels and Trappeniers 50 using the standard SW potential model.As can be seen, there is an excellent agreement between both results at all temperatures and densities (AAD% ≈ 3%).It is important to remark that estimations obtained by Michels and Trappeniers were calculated using discontinuous Molecular Dynamics and the Green-Kubo formula for a system with the same size but interacting via the standard SW intermolecular potential. 50t this point, we have shown that it is possible to accurately predict the self-diffusion coefficient of SW molecules using the CSW intermolecular potential.However, since the self-diffusivity is system-size dependent, it is also useful to estimate D using different number of molecules.This will allow to extrapolate D to the limit of thermodynamic limit (N → ∞). 59Unfortunately, this requires to perform additional simulations.Instead of this, we choose an alternative approach that involves simulating a system formed by 1000 molecules and the use of the Yeh and Hummer 66 correction given by: where D ∞ is the self-diffusivity in the infinite system, D sim is the self-diffusion calculated simulating a system with N particles in a cubic box of side L, k B is the Boltzmann constant, T the temperature, ξ = 2.837298, and η the shear viscosity.Note that we use in Eq. ( 6) the estimated values of η obtained using the CSW intermolecular potential as it is explained in the following Section (all the self-diffusion coefficient values obtained in this work for systems formed by 108 and 1000 molecules, as well as the value obtained using Eq. ( 6), are included in the Supporting information).Fig. 3 shows the self-diffusion coefficient, as a function of density, at supercritical and subcritical temperatures as obtained from MD-NVT simulations of systems formed by 1000 molecules that interact via the CSW potential.As can be seen, the values obtained simulating 1000 molecules are slightly higher than those obtained using 108 molecules, as expected.For this reason we do not compare both sets of data in a single plot.As in Fig. 2, we also use a logarithmic D * T * − ρ * representation in order to better show the behavior of the selfdiffusion coefficient as a function of density.As can be seen, D behaves as a monotonically decreasing function of density and shows a change in the curvature for ρ * < 0.2.It is re- markable that the same change in curvature is also predicted in the work of Michels and Trappeniers. 50In all cases, λ = 1.5 σ .

B. Shear-viscosity
We have also considered the behavior of the shear-viscosity, as a function of density, of spherical SW molecules in the same range of temperatures previously studied.To this end, we use the Green-Kubo formula given by Eq. ( 2).This requires the calculation the stress tensor throughout MD-NVT simulations saving information for short periods of time to correctly capture the behavior of the auto-correlation function in Eq. (2).According to the literature, 59 it is recommendable to save the components of the pressure tensor each 5 to 10 fs for systems interacting through continuous intermolecular potentials.However, we have checked that systems of spherical SW molecules that interaction via the CSW intermolecular potential need to save the components of the pressure tensor every 0.5 fs.
It is possible to determine the diagonal, P xx , P yy , P zz , and off-diagonal, P xy , P xz , P yz , components of the pressure tensor from the same simulations performed to determine the self-diffusion coefficients.In particular, we average the autocorrelation function of several components of the pressure tensor, the three off-diagonal components and the equivalent terms (P xx − P yy )/2 and (P yy − P zz )/2. 67e use a home-made program to analyze and calculate the shear viscosity curves, as functions of density, at several subcritical and supercritical temperatures.As for the case of diffusivity, the reported data (averages and uncertainties) correspond to the average of 5 independent runs for each thermodynamic state studied.In each simulation, the value of the viscosity corresponds to the time-average value of η(t) from times at which the plateau is reached until the upper time limit in Eq. ( 2).In all cases, this time limit is fixed at 100 ps.The time at which the plateau occurs depends on the thermodynamic conditions: ∼ 35 ps, for low temperatures and densities, ∼ 20 ps at low temperatures and high densities or vice versa, and ∼ 8 ps for high temperatures and densities.It is interest- ing to mention that it is possible to estimate each value of η in a different manner.Particularly, one could fit the data using some analytic function, as suggested by Maginn et al., 59 and to obtain results close to those obtained using in this work.Fig. 4 shows the shear viscosity, as a function of density, at different temperatures as obtained from MD-NVT simulation of systems formed by 108 and 1000 molecules interacting through the CSW intermolecular potential with range λ = 1.5 σ .We have also included the simulation data obtained by Michels and Trappeniers. 46,47As can be seen, agreement between results obtained using the standard SW potential and the CSW implementation of Zerón 56 is excellent in the whole range conditions of temperature and density.Particularly, AAD% ∼ 3% for simulations using 108 molecules.The simulation results obtained using 108 and 1000 molecules are very similar, indicating that shear-viscosity is a collective property that it does not depend on the system size.A detailed account of the results obtained in this section are tabulated and included in the Supplementary Material.

C. Equilibrium properties
We now concentrate in equilibrium magnitudes, including structural properties, of the SW fluid determined using the CSW version of the intermolecular potential.Particularly, we focus on internal energy, pressure, and radial distribution function.Note that simulation data in the literature is less scarce in this case since equilibrium properties can be routinely computed using standard MC simulation codes.The reason is that only the intermolecular potential, and not the intermolecular forces (that exhibit discontinuities), are needed to evaluate properties using standard MC moves.
Equilibrium properties of the SW fluid interacting via its CSW version have been previously investigated by Zerón et al. 56 However, the previous study has focused on the behavior of the system at a reduced number of temperatures and densities.Taking advantage of the simulations performed to predict transport properties of system formed by 1000 molecules, with potential range λ = 1.5 σ , we have calculated the internal energy, compressibility factor, and radial distribution function at all states analyzed in this work.This allows to compare equilibrium properties of the SW fluid obtained from MC simulations using the standard SW potential and those obtained from MD-NVT simulations using the CSW version of the intermolecular potential.The compressibility factor, Z, that measures the nonideality a fluid system (Z = 1), is defined as: where P is the pressure, V the volume, k B the Boltzmann constant, and T the temperature of the system.Fig. 5 shows the compressibility factor, as a function of the density, in a wide range of temperatures for SW fluids interacting via the CSW version of the potential with range λ = 1.5 σ .We have also included in the plot MC simulation data taken from the literature. 14,30,46,68,69As we can see, there is a good agreement between the estimations obtained using both versions of the SW intermolecular potential.Particularly, we find that AAD% < 1%.
We have also obtained the internal energy, as a function of density, in a wide range of temperatures.Isotherms for the reduced excess internal energy are represented in Fig. 6 obtained using both versions of the SW potential with a range λ = 1.5 σ .Again, agreement between results obtained using both intermolecular potentials is excellent in the whole range of temperatures and densities (AAD% ∼ 0.3%).The numerical data of u * and Z for the system formed from 1000 molecules interacting via the CSW version of the potential are tabulated in the Supporting Information.
Finally, we have determined two structural properties of the SW fluid.Particularly, we analyze the radial distribution function of the fluid.As a representative example, here we only show three g(r * ) curves at low, intermediate, and high densities and different temperatures.MD-NVT simulation predictions obtained using the CSW version of the intermolecular potential are compared with data taken from the literature 14,37 as obtained using MC simulations using the standard SW potential with range λ = 1.5 σ .As can be seen in in Fig. 7, agreement between MC and MD simulation data is excellent in all cases.
The coordination number of molecules in a fluid can be easily computed from the radial distribution function.This property, N c , accounts for the average number of neighbor particles within the well.N c can be readily calculated through the following expression: 1][72][73][74] As can be observed in Fig. 8, agreement between MC simulation data (using the standard SW potential) taken from the literature and the results obtained in this work using the CSW version of the intermolecular potential is excellent in all cases (AAD ∼ 0.4%).In all cases, λ = 1.5 σ .

V. CONCLUSIONS
We have performed MD simulations in the NVT or canonical ensemble of systems formed from 108 and 1000 molecules interacting through the continuous square-well potential with range λ = 1.5 σ .Conditions at which simulations are performed correspond to the gas and liquid phases, including subcritical and supercritical regions of the phase diagram.Transport properties (shear viscosity and self-diffusion coefficient), as well as equilibrium (internal energy and compressibility factor) and structural properties (radial distribution function and coordination number) have been determined in a wide range of temperatures and pressures.We have compared the results obtained using the combination of MD simulations and the continuous version of the square-well potential with MC simulation data obtained using the standard SW potential taken from the literature.
The absolute average relative deviation (AAD%) between our results and those taken from the literature found for transport properties is 3%, approximately, for self-diffusion coefficient and shear viscosity.Agreement in the case of equilibrium properties is even better, with AAD% is below 1%.It is important to recall here that we have not observed a trend to systematically subestimate or overestimate data, which corroborated that deviations in predictions are part of the intrinsic uncertainty associated to the simulations.
The main conclusion of this work is that the continuous version of the SW potential proposed by Zerón et al. 56 is able to reproduce accurately not only equilibrium and structural, but also transport properties of the SW fluid with potential range λ = 1.5 σ in a wide range of densities and temperatures.

I. INTRODUCTION
Numerical data of all properties analyzed in the main paper for systems of 108 and 1000 particles which interact with the continuous square well (CSW) potential of range λ = 1.5 are tabulated in this material.
FIG. 2. Self-diffusion coefficient D * , as a function of density, ρ * , of a system formed from 108 particles interacting with the CSW potential (filled circles), as obtained from MD simulations, and with the SW potential (open circles), taken from Michels and Trappeniers, 50 at several temperatures.In all cases, λ = 1.5 σ .

FIG. 3 .
FIG.3.Self-diffusion coefficient D * , as a function of density, ρ * , of a system formed from 1000 particles interacting with the CSW potential, as obtained from MD simulations at several temperatures.In all cases, λ = 1.5 σ .

5 FIG. 4 .
FIG.4.Shear-viscosity, η * , as a function of density, ρ * , of systems formed from 108 (crosses) and 1000 particles (filled circles) interacting with the CSW potential at several temperatures.We have also included the results corresponding to a system formed from 108 particles (open circles) interacting with the SW potential at the same temperature.46,47In all cases, λ = 1.5 σ .

FIG. 5 .
FIG. 5. Compressibility factor, Z, as a function of density, at several temperatures.Filled circles are the MD predictions using the CSW potential and open symbols are the MC data of Largo and Solana 68 (circles), Henderson et al. 14 (squares), Labik et al. 30 (diamond), Gil-Villegaset al. 69 (triangle up), and DMD predictions of Michels and Trappeniers 46 (triangle down) using the SW potential.In all cases, λ = 1.5 σ .The continuous curves are guides to the eye.

5 FIG. 6 .
FIG.6.Reduced internal energy of excess, u * ex = U ex /Nk B T , as a function of density, ρ * , at several temperatures.The meaning of the symbols and curves are the same as in Fig.5.

FIG. 7 . 5 FIG. 8 .
FIG.7.Radial distribution function, g(r * ), as a function of reduced distance, r * ≡ r/σ , at different densities and temperatures.Open circles represent the MC data using the SW potential taken from the the work of Largo et al.37  and the open squares from the work of Henderson et al.14The continuous line are obtained from MD simulations using the CSW potential.In all cases, λ = 1.5 σ .

TABLE I :
Self-diffusion coefficient D * for square well (λ = 1.5) fluids as function of reduced density ρ * and temperature T * from MD simulations using 108 particles interacting with the CSW potential of the same range.Number in parenthesis correspond to the uncertainty in the last digit.

TABLE II :
Self-diffusion coefficient D * for square well (λ = 1.5) fluids as function of reduced density ρ * and temperature T * from MD simulations using 1000 particles interacting with the CSW potential of the same range.Number in parenthesis correspond to the uncertainty in the last digit.

TABLE III :
Self-diffusion coefficient at infinite system D * ∞ for square well (λ = 1.5) fluids as function of reduced density ρ * and temperature T * from MD simulations using 1000 particles interacting with the CSW potential of the same range.

TABLE IV :
Shear viscosity η * for square well (λ = 1.5) fluids as function of reduced density ρ * and temperature T * from MD simulations using 108 particles interacting with the CSW potential of the same range.Number in parenthesis correspond to the uncertainty in the last digit.

TABLE V :
Shear viscosity η * for square well (λ = 1.5) fluids as function of reduced density ρ * and temperature T * from MD simulations using 1000 particles interacting with the CSW potential of the same range.Number in parenthesis correspond to the uncertainty in the last digit.

TABLE VI :
Compressibility Factor

TABLE VII :
Reduced internal energy of excess u * ex = U ex /Nǫ for square well (λ = 1.5) fluids as function of reduced density ρ * and temperature T * from MD simulations using 1000 particles interacting with the CSW potential of the same range.Number in parenthesis correspond to the uncertainty in the last digit.

TABLE VIII :
Coordination number N c for square well (λ = 1.5) fluids as function of reduced density ρ * and temperature T * from MD simulations using 1000 particles interacting with the CSW potential of the same range.Number in parenthesis correspond to the uncertainty in the last digit.