Many-Body Dissipative Particle Dynamics Simulation of Liquid–Vapor

Article pubs.acs.org/jced

Many-Body Dissipative Particle Dynamics Simulation of Liquid− Vapor Coexisting Curve in Sodium Maryam Atashafrooz and Nargess Mehdipour* Department of Chemistry, College of Sciences, Persian Gulf University, Boushehr 75168, Iran ABSTRACT: It is known that pair potentials are not adequate to correctly describe the structure and dynamics of liquid metals. In this work many-body interactions in fluid sodium are taken into account through the generic form of the forces acting between fluid particles in the context of the many-body dissipative particle dynamics (DPD) simulation. It is shown that this method accurately takes into account the role of many-body interactions in the improvement of the equation of state and structural properties of sodium. A comparison of our calculated coexisting liquid and vapor densities with experiment indicates that this model well predicts the vapor pressure and the densities over a broad range of temperatures. Also the calculated radial distribution functions are in very close agreement with experimental values. However, the calculated diffusion coefficient of sodium at low temperature deviates considerably from experiment. Although the soft-core “pair” potential in the original DPD method results in an equation of state exhibiting no fluid−fluid phase transition, the present “many-body” DPD approach successfully predicts the liquid−vapor coexistence curve of sodium. This occurs because, in the many-body DPD model, the force acting on a pair of atoms depends not only on their positions but also on the positions of all their other neighboring atoms. Also calculations are done using “pair” potentials and using the potential of mean force (in which the effect of many-body interactions are implicitly taken into account). It is shown that the results obtained using the many-body DPD approach are more accurate than those obtained employing both pair potentials and potential of mean force.

1. INTRODUCTION Because liquid alkali metals are used in engineering applications and are fundamentally important, the knowledge of their properties is of particular interest. They are used for applications as liquid metal-cooled reactors and for converting thermal energy to electrical power, for space applications, in the Rankine cycles, and in thermoelectric power generation systems. The wide range of the liquid state, over which the physical properties change substantially (for example, the occurrence of metal−nonmetal transition) as well as their high reactivity is an indication of the complexity of liquid alkali metals.1 Because of such complexities, experimental measurements on liquid alkali metals, especially at high temperatures and pressures are difficult. Theoretically, in many developed models for fluid metals, attention has been paid to the concept of “simple liquid”, in which the particles interact through spherically symmetric pairwise additive potentials.2,3 Also molecular simulation studies have been performed based on such pair potentials. However, it is now a well-established fact that pair potentials fail to accurately describe the microscopic structure of liquid metals.4,5 For example, it is known that the distinct collective modes, extending to wavelengths comparable to interatomic separations, affect the dynamics in molten alkali metals. In contrast, such excitations in simple liquids are only well-defined at long wavelengths and are considered as macroscopic density fluctuations.4 The nature of such distinct behaviors in liquid metals, compared to that of simple liquids, is in the specific interparticle interactions. © XXXX American Chemical Society

The above-cited descriptions indicate that the densityindependent pair potentials, such as the Lennard-Jones (LJ) potential, can poorly describe the interaction of liquid sodium atoms. It is worth mentioning that there exist modified forms of pair potentials6 for liquid alkali metals, with a significant improvement over their original simplified pair potentials,7 in the literature. Recently we have developed a LJ (5−4) pair potential for liquid sodium which successfully predicts the vapor−liquid coexistence curve.8 In fact, it is shown that a LJ pair potential with a softer repulsive branch and a longer-range attractive branch better predicts the pressure−volume−temperature properties of liquid metals than the conventional LJ (12− 6) potential. Meanwhile, we would like to emphasize that even in “pair” potentials, the effect of many-body interactions are implicitly taken into account. This is due to the fact that such potentials are parametrized against experimental bulk liquid densities, where many body interactions play an important role at the state of parametrization. To take into account the manybody interactions, Raabe and Sadus4 have added an empirically determined (temperature-dependent) contribution to the pair potential for liquid mercury. Perhaps difficulties in rationalizing the many-body effects and the relative ease of calculations using pairwise additive forces have led them to use an almost pairwise additive potential, in which the temperature is included explicitly, in their potential. However, their calculations indicate Received: July 4, 2016 Accepted: September 2, 2016

A

DOI: 10.1021/acs.jced.6b00586 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

that even the inclusion of such an indirect estimate of manybody forces considerably improves the accuracy of their results. In fact the potential energy of interaction between alkali metal particles results from the collective motions of electrons and ions. Therefore, the effect of electron gas on the ionic core must be taken into account in the potential energy of alkali metals. This means that such a potential must be a function of the electron density, and hence, of the metal density. On the basis of this concept, a number of effective density-dependent potentials have been suggested in the literature for alkali metals.4,9,10 It is known that potentials based on embedded atom models (EAM), in which the potential energy of metal is composed of a pairwise additive contribution plus an embedding contribution, depending on the effective electron density, are successful in predicting physical properties of liquid metals.11−13 In this respect, Monte Carlo and molecular dynamics simulations on liquid sodium are done using the EAM potentials.11−13 However, a disadvantage of EAM potentials is the large number of their parameters; in the case of liquid sodium, which is the subject of the present work, the number of parameters vary from 511,15 to 20.13 The proper way of introducing manybody effects into the interactions depends on the nature of the system of interest. In this work, our key to address the manybody interactions in sodium is through comparing the simulation results with the microscopic structure of the fluid (the radial distribution function) and the experimental pressure−volume−temperature (equation of state) data. In principle, we need a method to link between the macroscopic equation of state of the fluid sodium and its effective interparticle forces. A proper solution to this problem is known to be through deriving the generic form of the forces acting between fluid particles as a function of the excess free energy that characterizes the system (as a macroscopic property). Such a density-dependent force is, in general, nonpairwise additive, and its role in improving the equation of state of a soft repulsive fluid is greeted in the context of the dissipative particle dynamics (DPD) simulation.14,15 In fact the original DPD simulation approach, with a soft repulsive potential, has been improved by introducing the local-densitydependent (many-body) interactions to the expression for interparticle forces.16 It is the purpose of this work to employ this approach to study the equation of state and the structural properties of fluid sodium, as a typical example of a metallic system, and to examine the role of many-body interactions in the improvement of its equation of state.

Fi =

∑ (FCij + FijD + FijR ) (1)

j≠i

The conservative force is a soft repulsive force, defined as

FCij = AwC(rij)eij

(2)

where rij is the distance between particles i and j, A is the maximum repulsion between them, and eij = rij /rij. The dissipative and random forces are expressed as15 FijD = −γw D(rij)(rij. vij)eij

(3)

and FijR = σw R(rij)θij eij

(4)

Here γ is the dissipation strength, σ is the noise strength, vij = vi − vj (vi and vj being the velocities of particles i and j, respectively), θij is a random variable with a Gaussian distribution and unit variance. The weight functions wD and wR are defined as follow:

wC(rij) = w R (rij) =

rij ⎧ 1− (rij ≤ rC) ⎪ rC w D(rij) = ⎨ ⎪ (rij > rC) ⎩0

(5)

where rC is the cutoff diatnce (rC = 1 in dimensionless units). The random force heats up the system, whereas the dissipative force reduces the relative velocities of particles, tending to cool down the system. Therefore, the dissipative and random forces act together as a thermostat. Although the DPD method has received particular attention for simulation of fluids of a variety of molecular shapes and complexities, its soft-core conservative pair potential results in an equation of state with fixed quadratic density dependence, exhibiting no fluid−fluid phase transition, which is quite different from the equation of state of a real fluid. To overcome this difficulty, Pagonabarraga and Frenkel16 developed the many-body DPD model to describe the thermodynamic behavior of the fluids. In this method the amplitude A in the conservative force expression depends somehow on the density. Later Trofimov et al.17 and Warren18 introduced a local density expression into the conservative force. This many-body DPD differs from the original DPD model in that the definition of conservative force depends not only on the interparticle separation, but also on the instantaneous local particle density, which in turn depends on the positions of all other neighboring particles. In other words, the conservative force in the manybody DPD is a many-body force. On the basis of the desired thermodynamic behavior of the fluid, Pagonabarraga and Frenkel16 derived the functional form of the density-dependent conservative forces. According to the Warren method,18 manybody DPD is an extended version of standard DPD in which the soft pair potential is attractive (A < 0) and a repulsive many-body contribution with a shorter range is added; that is, eq 2 for the conservative forces is modified as17,18

2. THEORY As stated in the preceding section, we are going to address the many-body interactions in fluid sodium in the framework of the development of a DPD simulation approach to modify the equation of state of a soft-repulsive fluid. Therefore, we give a historical background of the development of the DPD simulation with the main emphasize on the role and manifestation of the many-body forces. The DPD is a simulation method, introduced by Hoogerbrugge and Koelman14 and later reformulated by Españoland Warren,15 to study the hydrodynamic behavior of complex fluids. In this method the fluid is composed of a number of particles, the motions of which obey Newton’s law. The force is a contribution of pairwise conservative, dissipative, and random forces; that is,

FCij = AwC(rij)eij + B(ρi + ρj )wαC(rij)eij

(6)

where B is the strength of the repulsive forces, wC(rij) is defined in eq 5, and B

DOI: 10.1021/acs.jced.6b00586 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data rij ⎧ 1− (rij ≤ αrC) ⎪ αrC wαC(rij) = ⎨ ⎪ (rij > αrC) ⎩0

Article

constant temperature, T, and constant pressure, P (the NPT ensemble), with the same method and the same parameters as explained in section 1, above. (3) The many-body DPD simulations in which the parameters in eqs 6−8) were tuned against the experimental pressure−volume−temperature data of fluid sodium. Because of the higher accuracy of the experimental liquid and vapor densities at lower temperatures, the DPD parameters were tuned against about 10 experimental data points in the low temperature regime.25 The parameters in reduced and real units are given in Table 1. The neighbor list length was 2.1 nm, and the time step was 4.0 fs. The DPD simulations were performed with a homemade code.

(7)

Here the parameter α (α < 1) determines the range of the attractive many-body forces, and ρi indicates the average local density around particle i and is expressed as ρi =

∑ j≠i

2 rij ⎞ 15 ⎛ 1 − ⎜ ⎟ αrC ⎠ 2πα 3rC3 ⎝

(8)

Although the original DPD formulation, as explained above, has been applied extensively to study systems of varieties of complexities,14−18 for liquid metals the inclusion of many-body interactions in the potential energy of interaction is known to be necessary. Therefore, in this work we examine the effect of inclusion of many-body forces in the original DPD formulation to predict the structural and dynamical properties of fluid sodium. To this end we did our calculations with a simplified many-body DPD model, as explained above, and compared the results with the predictions from simulations employing “pair” potentials. Also the results have been compared with the predictions from the potential of mean force (obtained by inversion of the radial distribution function at 378 K), in which the “many-body” force, at the state point it is derived (378 K), is included.

Table 1. Many-Body DPD Parameters in Reduced and Real Unitsa value in Parameter

dimensionless units

rC kT m α γ A B

1 1 1 0.25 0.212 −1.534 14.13

real units

conversion factor

2.0 23.0 0.1 −0.95 300.0

unit nm kJ mol−1 g mol−1

(mkT)0.5/rC4 kT/rC2 kTrC/α

kJ mol−1 nm−4 ps kJ mol−1 nm−2 kJ nm

a

k and m are the Boltzmann constants and the atomic mass of a sodium atom.

3. SIMULATION AND METHOD In this work all simulations of the liquid phase were done for 2500 sodium atoms. In the case of the gaseous phase the number of sodium particles ranged from 100 to 2000 particles, depending on the temperature (vapor pressure). Three types of simulations were done: (1) Molecular dynamics (MD) simulations, employing the pair potentials for which the force field parameters are given in the literature. In this case a LJ (12−6) pair potential, with a set of potential energy parameters tuned by Bhansali et al.19 (ε = 4.98 kJ mol−1 and σ = 0.324 nm), a new pair potential, developed by Ramana6 based on the improvement of the Morse potential,7 and our recently proposed LJ (5−4) potential8 were used. Simulations were done with the YASP simulation package.20 The cutoff distance was 2.0 nm and neighbors were included up to a distance of 2.2 nm. The temperature and pressure were kept constant using a Berendsen thermostat and a Berendsen barostat.21 The time constants for temperature and pressure couplings were 0.2 ps, and 1.0 ps, respectively. (2) MD simulations in which a potential of mean force, constructed using the so-called iterative Boltzmann inversion method,22 was used to describe the interaction of sodium atoms. The details of the method can be found in the literature.22,23 Here the method is just briefly described. To construct the potential of mean force, one must be aware of the target radial distribution function (in this case we have used the experimental radial distribution function24 reported at T = 378 K). Inversion of such a radial distribution function gives an initial guess for the potential of mean force. The initial guess is then improved in an iterative cycle so that the target and the calculated (from simulations using the potential of mean force) radial distribution functions closely resemble each other.22,23 Simulations were done for a constant number of particles, N, at

To calculate the equilibrium phase transition point, a method developed by Eslami and Müller−Plathe was used.26,27 The details of the method can be found in the literature;26,27 however a summary of the method is give here. A simulation of the liquid phase, in the NPT ensemble, is done and the excess chemical potential is calculated using the Widom’s test particle insertion method.28 Then a simulation of the vapor phase, in the grand canonical ensemble, is done in which the chemical potential is fixed according to the chemical potential of the liquid phase and its variation with pressure. The vapor phase then adjusts its state to the liquid phase (so that the two phases have the same chemical potential, temperature, and pressure). All samples were equilibrated for 2 ns. After this initial equilibration step, production runs were done for 5 ns.

4. RESULTS AND DISCUSSION 4.1. Liquid Vapor Densities. We have done simulations on liquid and vapor sodium at different temperatures varying from about 378 to 2200 K and at saturation pressures. The liquid and vapor densities calculated employing the many-body DPD model, as explained here, are shown in Figure 1 and are compared with experiment25 over a wide range of temperatures and pressures. The calculated liquid−vapor coexisting densities, obtained employing this method, are compared with experimental data25 in Figure 1. It is seen that the calculated liquid and vapor densities well agree with the experiment up to temperatures close to the critical temperature (≈2500 K). As it is expected, due to large fluctuations close to the critical point, the simulation results at very high temperatures (close to the critical temperature) show considerable deviations from experiment. The calculated vapor pressure curve is also compared with experiment25 in Figure 2. Again, except at very high temperatures, a very good agreement with C

DOI: 10.1021/acs.jced.6b00586 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

not shown here. However, the reported coexistence curve for sodium in ref 6 shows a substantial deviation (about 30 %) from the experimental curve. The results of LJ (5−4) potential are, however, comparable with the present results (at high temperatures the present many-body DPD results are closer to experiment that those predicted using the LJ(5−4) potential). This is expected, as in both methods potential energy functions with a soft repulsive branch and a long-range attractive branch are used. To take the role of many-body forces (implicitly) into account, a standard method in the literature is to construct the potential of mean force. The liquid densities calculated employing the potential of mean force (constructed at 378 K) are given in Figure 1 and are compared with experimental data.25 It is seen that the results are pretty accurate in the vicinity of the state point at which the potential of mean force is constructed. This is a known feature of this method.26,27 In other words, the potential of mean force implicitly includes the many-body interactions. However, as the potential of mean force is parametrized at 378 K and the many-body interactions are temperature/density-dependent, it does not reflect a correct picture of the role of many-body interactions at temperatures/ densities deviating considerably from the initial state point at which it is constructed. The results presented in this section indicate that respecting the role of many-body interactions in fluid sodium considerably improves its equation of state. The accuracy of predicted liquid and vapor densities and the vapor pressure using the manybody DPD model, as compared with experiment, is quite comparable with that using the EAM potential13 up to temperatures as high as 2000 K. However, the EAM potential better represents the phase coexistence region at high temperatures (it predicts the critical temperature at about 2300−2400 K). 4.2. Radial Distribution Function. Although the temperature-dependence of calculated saturated liquid and vapor densities, in Figure 1, indicates a good test for the ability of the present many-body DPD model to describe macroscopic properties of fluid sodium over a wide temperature range, an examination of the microscopic structure of the fluid, showing the packing of atoms around each other, is much more involved. This provides us with a direct comparison between simulation and experiment by comparing the calculated microscopic structure of liquid with that measured experimentally from X-ray or neutron diffraction studies. Here we have extracted the radial distribution functions of liquid sodium at 378 and 723 K and compared the results with experiment24 in Figure 3. It is seen that the calculated g(r) curves, using the present many-body DPD potential is in a very good agreement with experiment. Also the positions (intensities) of our calculated first g(r) peaks at 378 and 723 K, 0.37 nm (2.5) and 0.374 nm (1.99), respectively, are in close agreement with the corresponding values obtained using the EAM potential, 0.369 nm (2.5) and 0.373 nm (2.0), respectively. We have done similar calculations with the pair potential, by Ramana,6 and reported the calculated g(r) curves in Figure 3. Compared to experimental data, the g(r) curves, obtained from Ramana’s pair potential,6 are higher in intensity and are slightly shifted toward shorter distances. 4.3. Diffusion Coefficient. In this work we have calculated the diffusion coefficients of sodium over a wide range of temperatures to examine the validity of the present potential energy function in prediction of dynamics of liquid sodium.

Figure 1. Comparison of the calculated liquid−vapor coexisting curve for sodium with experimental data25 (open circles). The experimental vapor densities at high temperatures (open squares) are extrapolated using the experimental liquid densities via the law of rectilinear diameters. The markers indicate calculations using the LJ (12−6) potential (downward triangles), potential of mean force (upward triangles), the LJ (5−4) potential (diamonds), and many-body DPD potential (filled circles).

Figure 2. Calculated radial distribution functions of liquid sodium, at 378 and 723 K, compared with experiment (symbols).24 The full and dashed curves show many-body DPD simulation results and molecular dynamics simulation results, using the pair-potential by Ramana,6 respectively.

experimental vapor pressure results is seen. As stated above, the original DPD (without inclusion of many-body interactions) exhibits no liquid−vapor phase transition. This already reveals the role of many-body interactions in the improvement of the equation of state of fluid sodium. However, to have an estimate of the ability of “pair” potentials to predict the liquid−vapor coexisting curve of sodium, a comparison is made between the results generated using the many-body DPD approach and the those generated using three “pair” potentials, namely the LJ (12−6) potential,19 an improved pair potential by Ramana,6 and our recently proposed LJ (5−4) potential.8 In Figure 1 the temperature-dependence of the density of liquid sodium and its equilibrating vapor, calculated using the LJ (12−6) potential, is shown and is compared with experimental data.25 Unless it is at low temperature regime, over which the LJ (12−6) pair potential parameters are tuned, the agreement of the calculated densities with experimental data is poor. The predicted liquid− vapor coexisting densities, using Ramana’s pair potential,6 are D

DOI: 10.1021/acs.jced.6b00586 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

our calculated diffusion coefficients at high temperatures occurs because the present many-body DPD potential gives a density lower than the experimental values in the high temperature regime. Therefore, the trend observed for the deviation of diffusion coefficients agrees with the trend for deviations in liquid densities, compared with experiment.25

5. CONCLUSIONS The role of many-body forces in the interaction of liquid metal atoms is a long established fact. However, the inclusion of many-body interactions in the potential energy functions is not relatively easy. For example, the recently proposed accurate EAM potential for sodium consists of 20 parameters.13 This means that a large body of parametrization is required for this purpose. Moreover, the implementation of such a potential energy function in molecular simulation packages is not straightforward. In this work, a relatively simple approach is employed to take into account the many body interactions in sodium. This is done by introduction of the local-density-dependent (manybody) interactions to the original pair potential in the context of the DPD approach. It is shown that this method accurately takes into account the nature of many-body interactions in fluid sodium. An analysis of our calculated coexisting liquid and vapor densities and their comparison with experiment25 indicates that this model well predicts the vapor pressure and the densities from low (close to the triple point) to high (close to the critical point) temperatures. Also the calculated radial distribution function curves, from the results of the many-body DPD simulations well agree with experiment.24 Our results show that while the present many-body DPD approach predicts the structure of sodium in close agreement with experiment, its dynamics (examined in terms of the self-diffusion coefficient) deviates considerably from experiment both at high (T > 1200 K) and low temperatures (T < 550 K). The calculated diffusion coefficients of liquid sodium at low and high temperatures are lower and higher, respectively, than the experimental values.29,30 At the intermediate temperature regime, however, the agreement with the predictions from EAM potential13 is good. It is worth considering that the trend of observed deviations in the calculated diffusion coefficients with experiment follows the same trend in the deviation between calculated and experimental densities. Besides examining the validity of pair potentials for prediction of structural properties of sodium, calculations are done employing a simple LJ (12−6) potential,19 a pair potential parametrized by Ramana,6 and our recently proposed LJ (5−4) potential.8 The results show that only the LJ (5−4) potential, consisting a softer repulsive branch and a longer-range attractive branch than the conventional pair potentials, is adequate for description of the interactions in fluid sodium. Also we have checked the validity of a known method for implicit inclusion of many-body interactions in the potential energy function, that is, the construction of potential of mean force by tuning the calculated and experimental radial distribution functions. It is seen that in this case the liquid densities are in a good agreement with experiment just over a limited temperature range (≈200 K) around the temperature at which this potential is constructed. This means that the magnitude of many-body interactions largely depends on the temperature/density. While such a large temperature/densitydependency of the many-body interactions is taken into account implicitly in the potential of mean force at the state

Figure 3. Comparison of the calculated (open circles) and experimental25 (filled circles) vapor pressures of sodium as a function temperature.

The diffusion coefficients, calculated from the slopes of the mean squared-displacement curves, are reported in Figure 4

Figure 4. Temperature-dependence of the diffusion coefficient of liquid sodium. The filled and open circles represent our calculations using the many-body DPD potential and calculations using the LJ (5− 4) potential (diamonds) and the EAM potential13 (circles), respectively. The experimental data are the reported values by Ozelton and Swalin29 (triangles) and by Meyer and Nachtrieb30 (squares).

and are compared with the experimental values29,30 and with the results of simulations, using the EAM potential.13 A comparison of the calculated and experimental diffusion coefficients of sodium indicates that while the predictions of EAM potential well agree with experimental data, our calculated diffusion coefficients at low temperatures are considerably lower than the experimental values. As the experimental data at high temperatures are not available, we have compared our results with the calculated values from EAM potential.13 The agreement between our calculated diffusion coefficients with the calculations by Belashchenko13 is good in the intermediate temperature regime (550 K to 1200 K). At high temperatures our calculated diffusion coefficients are higher than the calculations from the EAM potential.12 Meanwhile, the calculated diffusion coefficients are compared with the results of the LJ (5−4) potential in Figure 4. The present many-body DPD results are more accurate than the predictions of LJ (5− 4) results at low temperatures. At high temperatures both methods have nearly the same accuracy. Positive deviation of E

DOI: 10.1021/acs.jced.6b00586 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

(14) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhys. Lett. 1992, 19, 155−160. (15) Español, P.; Warren, P. B. Statistical Mechanics of Dissipative Particle Dynamics. Europhys. Lett. 1995, 30, 191−196. (16) Pagonabarraga, I.; Frenkel, D. Dissipative Particle Dynamics for Interacting Systems. J. Chem. Phys. 2001, 115, 5015−5026. (17) Trofimov, S. Y.; Nies, E. L. F.; Michels, M. A. J. Thermodynamic Consistency in Dissipative Particle Dynamics Simulations of Strongly Nonideal Liquids and Liquid Mixtures. J. Chem. Phys. 2002, 117, 9383−9394. (18) Warren, P. B. Vapor-Liquid Coexistence in Many-Body Dissipative Particle Dynamics. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 066702. (19) Bhansali, A. P.; Bayazitoglu, Y.; Maruyama, S. Molecular Dynamics Simulation of an Evaporating Sodium Droplet. Int. J. Therm. Sci. 1994, 38, 66−74. (20) Müller-Plathe, F. YASP: A Molecular Simulation Package. Comput. Phys. Commun. 1993, 78, 77−94. (21) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (22) Müller-Plathe, F. Coarse Graining in Polymer Simulation, From the Atomistic to the Mesoscopic Scale and Back. ChemPhysChem 2002, 3, 754−769. (23) Eslami, H.; Karimi-Varzaneh, H. A.; Müller-Plathe, F. CoarseGrained Computer Simulation of Nanoconfined Polyamide-6, 6. Macromolecules 2011, 44, 3117−3128. (24) Waseda, Y. Structure of Non-crystalline Materials; McGraw-Hill: New York, 1980. (25) Vargaftik, N. B.; Yargin, V. S. IUPAC Handbook of Thermodynamic and Transport Properties of Alkali Metals; Ohse, R. W., Ed.; Blackwell Scientific: Oxford, 1985. (26) Eslami, H.; Müller-Plathe, F. Molecular Dynamics Simulation in the Grand Canonical Ensemble. J. Comput. Chem. 2007, 28, 1763− 1773. (27) Eslami, H.; Müller-Plathe, F. Molecular Dynamics Simulation of Sorption of Gases in Polystyrene. Macromolecules 2007, 40, 6413− 6421. (28) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808−2812. (29) Ozelton, M. W.; Swalin, R. A. Self-Diffusion in Liquid Sodium at Constant Volume and Constant Pressure. Philos. Mag. 1968, 18, 441− 451. (30) Meyer, R. E.; Nachtrieb, N. H. Self Diffusion of Liquid Sodium. J. Chem. Phys. 1955, 23, 1851−1854.

point it is constructed, its effect is not transferable to state points deviating substantially from the original point. In other words, such a method is not efficient for such systems as liquid sodium with a broad liquid domain. A comparison of the predicted results from the present many-body DPD approach and the predictions from the pair potentials and the potential of mean force with experiment indicates that the many-body DPD method is successful in taking into account the many-body interactions (by introducing the local density expression into the conservative forces) to improve the equation of state of fluid sodium. Also compared to the accurate EAM potentials, the many-body DPD approach employed here has nearly the same accuracy in predicting the structural properties of sodium up to temperatures as high as 2000 K. At higher temperatures, however, the EAM potential gives more accurate densities than the present many-body DPD potential. Also the calculated diffusion coefficients of liquid sodium using the EAM potential are more accurate than those from our calculations.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Funding

The support of this work by the research committee of Persian Gulf University is acknowledged. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Mehdipour, N. Correlation for the Viscosity of Liquid Alkali Metals. Fluid Phase Equilib. 2013, 355, 8−11. (2) Mehdipour, N.; Boushehri, A. An Analytical Equation of State for Mercury. Int. J. Thermophys. 1997, 18, 1329−1334. (3) Pilgrim, W. C.; Morkel, Ch. State Dependent Particle Dynamics in Liquid Alkali Metals. J. Phys.: Condens. Matter 2006, 18, R585− R633. (4) Raabe, G.; Sadus, R. J. Molecular Simulation of the Vapor−Liquid Coexistence of Mercury. J. Chem. Phys. 2003, 119, 6691−6697. (5) Gupta, R. P. Lattice Relaxation at a Metal Surface. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 6265−6270. (6) Sai Venkata Ramana, A. Molecular Dynamics Simulation of Liquid−Vapor Phase Diagrams of Metals Modeled using Modified Empirical Pair Potentials. Fluid Phase Equilib. 2014, 361, 181−187. (7) Singh, J. K.; Adhikari, J.; Kwak, S. K. Vapor−Liquid Phase Coexistence Curves for Morse Fluids. Fluid Phase Equilib. 2006, 248, 1−6. (8) Mehdipour, N.; Hasheminasab, F. Molecular Dynamics Simulation of Fluid Sodium. Fluid Phase Equilib. 2016, 427, 161−165. (9) Price, D. L.; Singwi, K. S.; Tosi, M. P. Lattice Dynamics of Alkali Metals in the Self Consistent Screening Theory. Phys. Rev. B 1970, 2, 2983−2999. (10) Eslami, H.; Mozaffari, F.; Moghadasi, J. Molecular Dynamics Simulation of Potassium Along the Liquid-Vapor Coexistence Curve. J. Iran. Chem. Soc. 2010, 7, 308−317. (11) Cai, J.; Ye, Y. Y. Simple Analytical Embedded-Atom-Potential Model Including a Long Range Force for fcc Metals and Their Alloys. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 8398−8410. (12) Ramana, A. S. V. Molecular Dynamics Simulation of LiquidVapor Coexistence Curves of Metals. J. Phys. Conf. Ser. 2012, 377, 012086. (13) Belashchenko, D. K. Computer Simulation of Liquid Metals. Phys.-Usp. 2013, 56, 1176−1216. F

DOI: 10.1021/acs.jced.6b00586 J. Chem. Eng. Data XXXX, XXX, XXX−XXX