Metal–Nonmetal Transition in Fluid Cesium: A Many-Body

Employing this many-body potential, simulations are done, in the grand ...... Edwards , P. P.; Lodge , M. T. J.; Hensel , F.; Redmer , R. A Metal Cond...
5 downloads 0 Views 1MB Size
Article Cite This: J. Chem. Eng. Data XXXX, XXX, XXX-XXX

pubs.acs.org/jced

Metal−Nonmetal Transition in Fluid Cesium: A Many-Body Dissipative Particle Dynamics Simulation Approach Golandam Faramarzi and Nargess Mehdipour* Department of Chemistry, Persian Gulf University, Boushehr 75168, Iran S Supporting Information *

ABSTRACT: A many-body dissipative dynamics simulation approach is presented for inclusion of many-body interactions in fluid cesium. Employing this many-body potential, simulations are done, in the grand canonical ensemble, to calculate the vapor−liquid equilibria over a wide range of temperatures (from the melting temperature to the critical temperature). The calculated coexisting liquid and vapor densities and the vapor pressures are in close agreement with experiment. The metal-nonmetal transition is examined in terms of cluster formation in low density cesium with increasing pressure. Our analysis of cluster formation along the critical isotherm shows that at pressures noticeably higher than the critical pressure very large spherical clusters are formed in the simulation box. The cluster sizes decrease with decreasing pressure to the critical pressure. At the critical pressure, only small clusters are seen in the simulation box. The cluster structures also change noticeably as the metal−nonmetal transition approaches. In the metallic state, high-density, spherical clusters are formed. Decreasing the pressure to the critical pressure causes rearrangement in the structure of spherical clusters to ellipsoidal (lower density) clusters. In the nonmetallic region, the clusters are very diffuse (low density) and do not have a well-defined structure. Adopting the fraction of cesium atoms involved in clusters larger than a threshold size as the fraction contributing to electrical conductivity, our simulation results well predict a sudden decrease in the electrical conductivity in close agreement with experimental observations.



INTRODUCTION Because of the environmental problems associated with the use of fossil fuels, attention has been focused on the renewable energy sources, among them concentrated solar power is the subject of recent interest and development.1 This technology is highly efficient at high operational temperatures, normally higher than 1000 °C. At such extreme conditions, conventional fluids cannot be employed as heat transfer agents. Due to high boiling points and high thermal conductivities, liquid metals are used as heat transfer fluids in such applications. Obviously, for such applications a molecular-level understanding of electronic and thermodynamic properties of liquid metals is essential in determining the right material for the desired application. From a scientific point of view, an important issue in the context of relationship between thermodynamics, structure, and interionic interactions in liquid metals is the metal−nonmetal (MNM) transition. Because of the low critical point of cesium (the critical temperature, pressure, and density of cesium are 1924 K, 9.2 MPa, and 380 kg m−3, respectively)3 it has been the subject of most experimental investigation of the liquid−vapor and MNM transitions.3−6 The experimental measurements3−6 and theoretical findings7−9 indicate strong interplay between the liquid−vapor coexistence and the MNM transition. Away from the critical point, the density difference between the coexisting liquid and vapor phases is large, and the MNM transition occurs at the same density as the liquid−vapor © XXXX American Chemical Society

transition occurs. Close to the critical point, however, it is difficult to experimentally determine the MNM transition. In this region, a strong effect of the phase transition on the electronic structure, and hence, on the conductivity is observed. However, a sharp (first-order) electronic transition is only observed across the liquid−vapor phase boundary. Although investigation of the phase behavior of fluid alkali metals has been the subject of several experimental and theoretical studies, the connection of phase behavior to MNM transition has remained elusive. The correlation between the density and the conductivity shows that the density decrease is an important factor governing the MNM transition. In 1940s Landau and Zeldovitch10 suggested that the first-order electronic transition and the liquid−vapor transition in liquid metals are separated. Subsequent theoretical efforts confirmed this idea.11,12 However, another idea, proposed by Krumhansl, is that the liquid−vapor and the MNM transitions coincide and there is only one critical point.13 Due to the strong interplay between fluctuations in the densities of liquid and vapor phases (close to the critical point) and the rapid change in the electronic structure at the MNM transition, it is not easy to distinguish between the two mechanisms. Experimental Received: July 14, 2017 Accepted: September 28, 2017

A

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

investigations on the magnetic susceptibility measurements14 for alkali metal atoms in the vapor phase show that a high concentration of molecular aggregates forms as the vapor density increases. The electron transfer between atoms and clusters and between clusters is facilitated due to cluster formation. The behavior of experimentally measured magnetic susceptibilities confirms formation of a large fraction of aggregates in the critical region of fluid alkali metals.14 In conjunction with these experimental measurements, there is a need for accurate molecular models for a better understanding of the molecular-level processes governing the macroscopic properties. So far, simplified theoretical approaches for the estimation of electrical conductivity in metal models have been developed.15−18 In principle, ab inito molecular dynamics simulations, with the main emphasis on the role of electrons in metal atom interactions, are proper tools for electrical conductivity calculation. Such simulations are, however, computationally expensive. Our main emphasis here is to develop a simple model for simulation of liquid− vapor phase transition and MNM transition in cesium, as a typical example. To calculate liquid−vapor phase coexistence points, molecular dynamics simulations in the grand canonical ensemble19 are done. The MNM transition is calculated based on an analysis in terms of clustering of cesium atoms.

Warren24 showed that the dissipative and random forces have to satisfy a certain relation, i.e., rij ⎧ 1− (rij ≤ rC) ⎪ rC w D(rij) = wC(rij) = [w R (rij)]2 = ⎨ ⎪ (rij > rC) ⎩0

where rC is the cutoff distance. The relationship between the amplitudes and the thermal energy is expressed as δ 2 = 2γkBT

THEORY A severe limitation of the theoretical approaches for studying fluid metals is the lack of accurate pair potential energy functions describing their interatomic interactions.20,21 In fact, no reliable theoretical method is available to derive an effective potential function that describes liquid alkali metals accurately. Even, inconsistencies between predictions of various potential models are reported.22,23 It is known, however, that a class of potentials based on embedded atom models (EAM) are successful in predicting physical properties of liquid metals.24 The EAM potential is composed of a pairwise additive contribution plus an effective electron density-dependent embedding contribution. Recently, we have proposed a manybody dissipative particle dynamics (DPD) simulation method for a proper inclusion of the many-body interactions in the liquid sodium.25 In what follows we describe the method, with the main emphasis on the role of many-body interaction. The particle interactions in DPD are formulated in terms of three types of forces; a conservative force FC, a dissipative force FD, and a random force, FR. The conservative force is a soft repulsive force, defined as26

where α (α < 1) determines the range of the attractive manybody forces,ρi indicates the average local density around particle i, and B is the strength of the repulsive forces. The local density around particle i is defined as ρi =

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

(7)

According to eq 6 the force between a pair of particles, depends on the local densities around both particles, and hence, is density/temperature-dependent. In our recent study we showed that this improved version of DPD successfully takes into account the role of many-body interactions in fluid sodium.25 It is worth mentioning that even for polar fluids many-body interactions may considerably improve the accuracy of the phase coexistence results.30 In the case of liquid metals, it is known that many-body interactions play an important role in their physical properties.25 In this work we employ the many-body DPD potential, as explained above, to perform simulations to calculate the liquid− vapor phase equilibria and the MNM transition in cesium. There are a number of well-known techniques for the calculation of free energies,31 and hence, the phase equilibria32 using molecular simulation methods. A review of the existing methods is given by Kofke.33 In this work, for the calculation of phase equilibria, simulations are done in the grand canonical ensemble.19 Here the temperature, volume, V, and chemical potential, μ, are set to predefined values and the number of particles is varied to reach equilibrium. To speed up the calculations we used a method developed by Vrabec and Hasse34 and applied for phase equilibrium calculations for sorption of gases in polymers.35,36 According to this method, to

(1)

(2)

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

∑ j≠i

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 as26 FijD = −γw D(rij)(rij. vij)eij

(5)

where kB is the Boltzmann constant and T is the temperature. The soft-core conservative pair potential of the DPD results in an equation of state with fixed quadratic density dependence, exhibiting no fluid−fluid phase transition. To solve this problem, Pagonabarraga and Frenkel27 developed the manybody DPD model to describe the thermodynamic behavior of the fluids. In the many-body DPD, the amplitude A in the conservative force expression depends on the density. Later Trofimov et al.28 and Warren29 introduced a local density expression into the conservative force. According to their method, many-body DPD is an extended version of standard DPD in which a repulsive many-body contribution is added, i.e.,



FCij = AwC(rij)eij

(4)

(3)

where γ and δ are the dissipation and noise strengths, respectively, 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. In order to have a DPD thermostat mentioning the canonical distribution, Español and B

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

find the phase coexistence point, a simulation of the liquid phase in the NPT ensemble is done. The excess chemical potential in the liquid phase is calculated using the Widom’s test-particle insertion method, i.e.,37 μex = −kBT ln

PV exp( −ΔU /kBT ) (N + 1)kBT

according to the excess chemical potential of the liquid phase (see eq 11). All simulations were carried out using the simulation package YASP.44 The temperature-independent many-body DPD parameters for cesium in both dimensionless and real units are reported in Table 1. The velocity-Verlet Table 1. Many-Body DPD Parameters in Reduced and Real Units

(8)

where P is the pressure, ΔU is the potential energy of interaction of the test particle with the rest of particles. Here μex is the excess chemical potential; the difference between chemical potential of the system and that of ideal gas, μid, defined as ⎡⎛ 2πmk T ⎞ B ⎟ μid = kBT ln⎢⎜ ⎢⎣⎝ h2 ⎠

3/2

kBT ⎤ ⎥ P ⎥⎦

value in

1 1 1 0.21 0.444 −1.804 5.61

(10)

kBT P

conversion factor

3.0

nm kJ mol−1 g mol−1

133.0 0.1 −0.5 200.0

unit

(mkBT)0.5/rC4 kBT/rC2 kBTrC/α

kJ mol−1 nm−5ps kJ mol−1 nm−2 kJ nm

N

(11)

d 2λ W 2 =− dt

with

v lex = v l −

real units

algorithm was used to integrate the equations of motion. The time step was 3 fs. All samples were equilibrated for 3 ns. After this initial equilibration step, a production run was done for 10 ns. For calculation of statistical uncertainties a total of 10000 blocks, with block size 100, were used for averaging. The block averages were independent. For performing simulation in the grand canonical ensemble, one of the particles of the system was chosen as a “fractional” particle, whose potential energy with the rest of the system was scaled by a fractional number λ (0 ≤ λ ≤ 1). During the course of the simulation, the fractional particle may grow to a full particle (to be added to the system), or may be removed from the system. The governing equation on the motion of the fractional number is19,39

where v is the molar volume and subscript l refers to the liquid phase. At the coexistence point (P*), the chemical potential of gas phase, μg(P*), is equal to that of the liquid phase. Therefore, performing a simulation of the gas phase in the pseudogrand canonical ensemble, the target excess chemical potential is set, as a function of pressure, according to the following equation: μgex (P*) = μ lex (P) + v lex (P* − P)

rC kBT m α γ A B

(9)

where N is the number of particles, m is the mass of particles, and h is the Planck’s constant. The excess chemical potential of the liquid phase is calculated, using the Widom’s test particle insertion method,37 from the NPT ensemble at a pressure P preferably close to the equilibrium vapor pressure (P*). Since the chemical potential of the liquid phase only slightly depends on the pressure, its value at the coexistence pressure, P* (not known), can be obtained by expanding it in terms of pressure as μ l (P*) = μ l (P) + v l(P* − P)

parameter

dimensionless units

∑ i=1

∂U (rif , λ) ∂λ

+ μex (13)

where W is the mass associated with the additional coordinate, λ, and U(rif, λ) is the potential energy of interaction between the fractional and ordinary particles. Upon growing a fractional particle to a full particle, or deletion of fractional particle from the system, decision was made either to insert a new fractional particle in the system or to choose a preexisting particle as the new fractional particle. Repeating such a procedure, equilibrium was achieved, in which the density was fluctuating around an average value, corresponding to the target chemical potential, temperature, and volume. The coexistence point was calculated from the results of this simulation.

(12)

It is shown that this later simulation forces the chemical potential of the gaseous phase to that of the liquid phase, and hence, the phase equilibrium point is found.34−36 This method is different from the Gibbs ensemble Monte Carlo simulation,38 in that it does not need simultaneous simulations of both phases; rather two separate simulations are needed to find the phase coexistence point.34 The method has been employed successfully for the calculation of phase equilibria in systems of varieties of complexities, ranging from the simple LennardJones fluid34,39 to polar fluids40 and gas solubility in polymers.41−43





RESULTS AND DISCUSSION Radial Distribution Function. Many-body DPD simulations on the liquid and vapor phases of cesium at temperatures ranging from 350 to 2000 K (critical point of cesium) were done. To adjust the many-body DPD parameters against experimental data, the calculated radial distribution functions, g(r), were compared with experiment,45 and the potential energy parameters were tuned to obtain a close match of our calculated radial distribution function curve with the experimental curve. In Figure 1 a comparison of the calculated and experimental45 radial distribution function curves at T =

SIMULATIONS In this work DPD simulations were performed for 2500 cesium atoms. The temperature was kept constant by adjusting the dissipation strength. Performing a simulation in the liquid phase in the NPT ensemble, the (excess) chemical potential is calculated, using the Widom’s test particle insertion method,37 by inserting 10000 test particles in 10000 independent configurations and averaging the Boltzmann factor. Then another simulation of the gaseous phase in the μVT ensemble is done, in which the target excess chemical potential is set C

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

Figure 1. A comparison of the our calculated (curves) and experimental (open markers)45 radial distribution functions for liquid cesium. The filled squares represent the ab initio molecular dynamics simulation results of Yuryev and Gelchinski.46 The temperatures are shown in the legend.

Figure 2. Calculated saturated liquid and vapor densities from our many-body DPD simulations in the grand canonical ensemble (filled circles) compared with experimental data (open circles2 and open squares50). The open (experimental)2,50 and filled (calculated) triangles show the average densities of coexisting liquid and vapor phases. The error bars, calculated within 68% confidence limit, are smaller than the sizes of markers. The critical exponent is ≈0.9.

373, 573, and 773 K are given. In all cases the calculated g(r) curves closely match with the experimental curves.45 We have also compared our calculated g(r) curve at 573 K with that calculated (at 600 K) in a recent ab initio molecular dynamics simulation by Yuryev and Gelchinski.46 As the results show, our calculated g(r) curve is in a better agreement with experiment than that calculated by Yuryev and Gelchinski.46 The positions and intensities of our calculated g(r) curves, however, well agree with the ab initio simulation results of Matsuda et al.47 The relatively large oscillations, extending to distances as long as 1.5 nm, showing strong correlations between motions of atoms, are typical for metallic systems. Another feature of the g(r) curves in Figure 1 is that the position of maximum in the first g(r) peak (≈ 0.55 nm) changes very slightly with temperature. Both observations are confirmed by experiment45 and by previous ab initio simulation results.47−49 Liquid−Vapor Phase Coexistence Curve. The calculated liquid−vapor coexistence curve of cesium is shown in Figure 2 and is compared with experimental data.2,50 As the results show, our calculated densities perfectly recover the liquid− vapor asymmetry in the coexistence curve of cesium and are in a very good agreement with experiment.2,50 The diameter, (ρl + ρv)/2, is also shown in Figure 2 and shows a nonlinear behavior described by ⎛ρ + ρ ⎞ ⎛ T − T ⎞ξ v ⎜⎜ l ⎟⎟ − 1 = ⎜ C ⎟ ⎝ TC ⎠ ⎝ 2ρC ⎠

We have further shown the calculated vapor pressures in Figure 3 and have compared the results with experimental

Figure 3. Clausius−Clapeyron plot of vapor pressure vs reciprocal temperature. The filled and open circles represent our calculations and experimental data,2 respectively. The error bars are smaller than the size of markers.

data.2 Except in the critical region, the calculated vapor pressures agree very well with experimental data. In this article all the uncertainties in the calculated values are obtained by calculating the standard deviation (within 68% confidence limit) for n independent data points, stored during the production run. 4.3. Diffusion Coefficients. The dynamics of cesium atoms has been examined, using the present many-body DPD method, in terms of the self-diffusion coefficients. The diffusion coefficients have been calculated using the Einstein formula. As the diffusion coefficient is affected by the DPD dissipation strength,54 in this work the DPD dissipation strength is adjusted against the diffusion coefficients at moderate temperatures. In Figure 4 we have shown the Arrhenius plot of ln(D) vs 1/T and compared the results with the calculated values

(14)

where ρl and ρv are the densities of coexisting liquid and vapor phases, ρC is the critical density, TC is the critical temperature, and ξ is the critical exponent. Our calculated critical exponent for cesium is 0.93, which is very close to the experimental value, 0.9.2,50 Deviation of critical exponent from unity shows a nonclassical critical behavior for fluid cesium, which is in agreement with experiment, with previous theoretical predictions.51,52 The coexistence curve of cesium is more asymmetric than that of normal fluids. This asymmetry is due to the stronger long-range interactions in cesium, which plays a dominant role on its structure (see Figure 1). This finding is in agreement with the theoretical findings of Kim and Fisher.53 D

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

To reduce the effect of sample size on the cluster sizes and their distributions, we have simulated 100000 cesium atoms. Our analysis for the clustering of cesium atoms along the critical isotherm shows that at high pressures (P ≈ 1.5 PC, where PC is the critical pressure) almost a single crystal exists in the simulation box. With decreasing the pressure, the cluster sizes become smaller. At low pressures, only small clusters are found in the simulation box. We have shown the distribution of cluster sizes in the simulation box at the critical point in Figure 5. The results show

Figure 4. Temperature-dependence of the diffusion coefficient of saturated liquid cesium. The filled and open circles represent our calculations and calculations using the EAM potential,24 respectively. In the inset the power-law temperature-dependence of the diffusion coefficients are shown.

using the EAM potential.24 A nearly good agreement between our calculated values and calculations based of the EAM potential is observed. In the critical region the calculated diffusion coefficients are higher than those calculated using the EAM potential. This trend is in line with the trend of deviation between calculated and experimental densities (see Figure 2). It is interesting to notice the nonlinear behavior of ln(D) vs 1/T, in Figure 4, revealing large temperature-dependencies of the activation energies for diffusion. This behavior is typical for liquid metals and can be attributed to factors, such as asymmetry in the coexistence curve and the wide temperature range over which cesium remains liquid (from the triple point ≈300 K to the critical point ≈2000 K). The error bars for calculated diffusion coefficients are obtained using the method of Pranami and Lamm.55 We have examined several functional forms, including the Vogel−Fulcher−Tammann expression56 to fit the diffusion coefficients with temperature. The diffusion coefficients of cesium over a wide range of temperatures exhibit a power-law dependence. Shown in the inset of Figure 4 is the temperaturedependency of the diffusion as the log−log plot of D vs T. The slope is 2, which is in agreement with previous molecular dynamics simulation results of Wax et al.57 4.4. Metal−Nonmetal Transition. A typical feature of the MNM transition is the structural change along the transition path. In this work, these structural changes have been investigated in terms of cluster formation, as the density of gaseous cesium increases. Here, clusters are defined as the group of atoms separated from their neighboring atoms by a threshold distance. In fact, for each cesium atom, a list of neighboring atoms is formed. Two neighboring atoms are regarded as connected atoms, if their distance is shorter than a threshold value. The search for connections further goes through the second- and higher-shell neighbors of a central atom. Here the threshold distance is chosen as 0.55 nm, corresponding to the position of maximum in the g(r) curves (see Figure 1). Our NPT simulation results show than this threshold distance is almost constant at both high and low pressures/densities. This is in agreement with the ab initio molecular dynamics simulation results in the literature.47−49

Figure 5. Distribution of cluster sizes in the simulation box at the critical point.

that a substantial amount of cesium atoms form clusters larger than 10 cesium atoms. Clearly, in our classical simulation scheme it not possible to calculate the electrical conductivity directly. However, based on experimental reports on the interplay between MNM transition and formation of molecular aggregates,14 we regard the fraction of cesium atoms involved in clusters larger than a threshold size (10 atoms) as conductive (metallic) fraction. Based on this choice, we have shown the fraction of system in the metallic/nonmetallic phase as a function of reduced pressure along the critical isotherm and compared the results with experimental4 reduced conductivities in Figure 6. Our results clearly follow the correct trend of experimental conductivities, with a sudden decrease in conductivity close to the critical pressure. In previous simulations in the literature, the change in coordination number is chosen as a criterion for MNM transition.48,58 While the coordination number of cesium atoms in the metallic region is higher than that of the nonmetallic region, no sharp change in the coordination number (characterizing the MNM boundary) occurs. Our cluster analysis method, however, reveals a sudden change in the cluster sizes and the fraction of metallic particles in the simulation box at the MNM boundary. Therefore, it provides correct features of structural changes along the MNM phase boundary. Another noticeable feature of clustering in cesium atoms is the change in clusters’ structures as the MNM transition occurs. In the metallic branch, high-density, spherical clusters are formed. In order to analyze the structure of clusters, a recent method applied by Eslami et al.59,60 for examination of local order in clusters, is applied. In this method the local order is measured in terms of the dot product of Steinhardt order parameter61 for neighboring particles. Based on our analysis, E

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

Figure 8. Pressure-dependence of the density along the critical isotherm for cesium. The filled and open markers represent the experimental data50 and our calculations, respectively. The error bars for calculated densities are smaller than the sizes of points. Figure 6. Ratio of experimental electrical conductivity, σ, to the conductivity at P = 1.5 PC, σ0, along the critical isotherm of cesium (open markers). The filled symbols indicate the fraction of cesium atoms involved in clusters of size larger than 10 atoms as a function of reduced pressure along the critical isotherm.



CONCLUSIONS A many-body DPD approach is presented for inclusion of many-body interactions in fluid cesium. Employing this manybody potential, simulations are done to calculate the vapor− liquid phase equilibria over a wide range of temperatures (from the melting temperature to the critical temperature). Our grand canonical ensemble simulations predict coexisting liquid and vapor densities and the vapor pressures in close agreement with experiment.2,50 The asymmetry in the liquid−vapor coexistence curve and the critical exponent perfectly match with experimental2,50 as well as theoretical53 findings. Also the microscopic structure of cesium atoms, expressed in terms of the radial distribution function, agrees nicely with experiment and with previous ab inito simulation results in the literature.47−49 We have further examined the dynamics of liquid cesium, along the liquid−vapor coexistence curve, by calculation the self-diffusion coefficients. The diffusion coefficients exhibit a power-law temperature dependence. The metal−nonmetal transition is examined in terms of cluster formation in low density cesium with increasing pressure/ density. Our analysis of cluster formation along the critical isotherm shows that at high pressures (P ≈ 1.5 PC) very large spherical clusters are found in the simulation box. Decreasing pressure toward the critical pressure, decreases the size of clusters. Very close to the critical pressure, only small, low density (diffuse) clusters are seen in the simulation box. In the nonmetallic side of the transition point, only very small clusters (aggregates less than 10 cesium atoms) are formed. The cluster structures also change noticeably as the MNM transition is approached. In the metallic state, high-density, spherical clusters are formed. Decreasing the pressure toward the critical pressure, rearranges the spherical clusters to ellipsoidal (lower density) clusters. In the nonmetallic region, the clusters are rather diffuse (low density) and do not have a well-defined structure. Due to the fact that the electronic charge density is localized at the nuclei in small low-density clusters, such aggregates have a very low contribution to the electrical conductivity. In contrary, the electronic charge density is delocalized, contributing to conductivity, over the larger (spherical) clusters. In this respect, our calculated fraction of conductive atoms represents a sudden change in conductivity, in close agreement with experiment.4

the clusters are liquid-like ordered. As the pressure/density decreases, the spherical clusters rearrange to ellipsoidal (lower density) clusters. In the nonmetallic region, the clusters are rather diffuse (low density) and do not have a well-defined structure. Snapshots of clusters in the metallic side, in the transition region, and in the nonmetallic side of the MNM transition are shown in Figure 7.

Figure 7. Snapshots of largest clusters, along the critical isotherm, in the simulation box. From left to right, the clusters belong to P = 1.1PC, 0.95PC, and 0.91PC, respectively.

Denser (spherical) clusters have a higher contribution to electrical conductivity than diffuser clusters. This is due to the fact that electronic charge density becomes inhomogeneous (localized) in diffuse clusters. In other words, the electrons localize at the nuclei in smaller and diffuser clusters. In contrary, the electrons are delocalized, contributing to conductivity, over the larger (spherical) clusters. Consequently, transition from a metallic to a nonmetallic phase occurs with the change in cluster sizes and structures. We have also shown in Figure 8 the variation of reduced density with the reduced pressure along the critical isotherm. The close correlation between the behavior of density and that of the conductivity (see Figure 6) with pressure along the critical isotherm shows that the density fluctuation is the main factor for observation of fluctuating clusters, and hence, the governing factor in MNM transition. F

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data



Article

(17) Reinaldo-Falagán, M.; Tarazona, P.; Chacón, E.; Velasco, E.; Hernandez, J. P. Hard-sphere Fluid with Tight-Binding Electronic Interactions: A Glue Model Treatment. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 67, 024209. (18) Tarazona, P.; Chacón, E.; Vergés, J. A.; Reinaldo-Falagán, M.; Velasco, E.; Hernandez, J. P. Electrical Conductivity of a Tight-binding Hard-sphere Model for Hot Fluid Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 024203. (19) Eslami, H.; Müller-Plathe, F. Molecular Dynamics Simulation in the Grand Canonical Ensemble. J. Comput. Chem. 2007, 28, 1763− 1773. (20) 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. (21) Desgranges, C.; Delhommelle, J. Thermodynamics of Phase Coexistence and Metal−Nonmetal Transition in Mercury: Assessment of Effective Potentials via Expanded Wang−Landau Simulations. J. Phys. Chem. B 2014, 118, 3175−3182. (22) Hasheminasab, F.; Mehdipour, N. Molecular Dynamics Simulation of Fluid Sodium. Fluid Phase Equilib. 2016, 427, 161−165. (23) Mehdipour, N.; Moosavi, F. A Perturbed Hard-Sphere Equation of State for Liquid Metals. Phys. Chem. Liq. 2011, 49, 347−354. (24) Belashchenko, D. K.; Nikitin, N. Yu. The Embedded Atom Model of Liquid Cesium. Russ. J. Phys. Chem. A 2008, 82, 1283−1289. (25) Atashafrooz, M.; Mehdipour, N. Many-Body Dissipative Particle Dynamics Simulation of Liquid−Vapor Coexisting Curve in Sodium. J. Chem. Eng. Data 2016, 61, 3659−3664. (26) Español, P.; Warren, P. B. Statistical Mechanics of Dissipative Particle Dynamics. Europhys. Lett. 1995, 30, 191−196. (27) Pagonabarraga, I.; Frenkel, D. Dissipative Particle Dynamics for Interacting Systems. J. Chem. Phys. 2001, 115, 5015−5026. (28) 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. (29) 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. (30) Wierzchowski, S. J.; Fang, Z. H.; Kofke, D. A.; Tilson, J. L. Three-body Effects in Hydrogen Fluoride: Survey of Potential Energy Surfaces. Mol. Phys. 2006, 104, 503−513. (31) Lu, N.; Kofke, D. A.; Woolf, T. B. Improving the Efficiency and Reliability of Free Energy Perturbation Calculations Using Overlap Sampling Methods. J. Comput. Chem. 2004, 25, 28−39. (32) Mehta, M.; Kofke, D. A. Coexistence Diagrams of Mixtures by Molecular Simulation. Chem. Eng. Sci. 1994, 49, 2633−2645. (33) Kofke, D. A. Free Energy Methods in Molecular Simulation. Fluid Phase Equilib. 2005, 228, 41−48. (34) Vrabec, J.; Hasse, H. Grand Equilibrium: Vapour-Liquid Equilibria by a New Molecular Simulation Method. Mol. Phys. 2002, 100, 3375−3383. (35) Eslami, H.; Müller-Plathe, F. Molecular Dynamics Simulation of Sorption of Gases in Polystyrene. Macromolecules 2007, 40, 6413− 6421. (36) Eslami, H.; Mehdipour, N. Grand Canonical Ensemble Molecular Dynamics Simulation of Water Solubility in Polyamide6,6. Phys. Chem. Chem. Phys. 2011, 13, 669−673. (37) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808−2812. (38) Panagiotopoulos, A. Z. Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a New Ensemble. Mol. Phys. 1987, 61, 813−826. (39) Eslami, H.; Mojahedi, F.; Moghadasi, J. Molecular Dynamics Simulation with Weak Coupling to Heat and Material Baths. J. Chem. Phys. 2010, 133, 084105. (40) Eslami, H.; Dargahi, A.; Behnejad, H. Molecular Dynamics Simulation of Liquid−Vapor Phase Equilibria in Polar Fluids. Chem. Phys. Lett. 2009, 473, 66−71.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jced.7b00639. All properties calculated in simulations (with uncertainties) (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Nargess Mehdipour: 0000-0002-3570-0870 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) Lorenzin, N.; Abanades, A. A Review on the Application of Liquid Metals as Heat Transfer Fluid in Concentrated Solar Power Technologies. Int. J. Hydrogen Energy 2016, 41, 6990−6995. (2) Vargaftik, N. B.; Yargin, V. S. IUPAC Handbook of Thermodynamic and Transport Properties of Alkali Metals; Ohse, R. W., Ed.; Blackwell Scientific: Oxford, 1985. (3) Renkert, H.; Hensel, F.; Franck, E. U. Metal-Nonmetal Transition in Dense Cesium Vapour. Phys. Lett. A 1969, 30, 494−495. (4) Hensel, F.; Uchtmann, H. The Metal-Insulator Transition in Expanded Fluid Metals,. Annu. Rev. Phys. Chem. 1989, 40, 61−83. (5) Edwards, P. P.; Ramakrishnan, T. V.; Rao, C. N. R. The MetalNonmetal Transition: A Global Perspective,. J. Phys. Chem. 1995, 99, 5228−5239. (6) Edwards, P. P.; Lodge, M. T. J.; Hensel, F.; Redmer, R. A Metal Conducts and a Non-metal Doesn’t. Philos. Trans. R. Soc., A 2010, 368, 941−965. (7) Tarazona, P.; Chacón, E.; Hernandez, J. P. Simple Model for the Phase Coexistence and Electrical Conductivity of Alkali Fluids. Phys. Rev. Lett. 1995, 74, 142−145. (8) Likalter, A. A. Equation of State of Metallic Fluids Near the Critical Point of Phase Transition. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53, 4386−4392. (9) Keshavarzi, E.; Parsafar, G. Prediction of the Metal-Nonmetal Transition using the Linear Isotherm Regularity. J. Phys. Chem. B 1999, 103, 6584−6589. (10) Landau, L.; Zeldovich, G. On the Relation between the Liquid and the Gaseous States of Metals. Acta Physicochim USSR 1943, 18, 194. (11) Odagaki, T.; Ogita, N.; Matsuda, H. Percolation Approach to the Metal-Insulator Transition in Super-Critical Fluid Metals. J. Phys. Soc. Jpn. 1975, 39, 618−624. (12) Yonezawa, F.; Ogawa, T. Metal-Nonmetal Transitions in Expanded Fluids−Electronic and Thermodynamic Properties. Suppl. Prog. Theor. Phys. 1982, 72, 1−80. (13) Krumhansl, J. A. In Physics of Solids at High Pressures; Tomizuka, C. T., Emrick, R. M., Eds.; Academic: New York, 1965; p 425. (14) Freyland, W. Magnetic Susceptibility of Metallic and Nonmetallic Expanded Fluid Cesium. Phys. Rev. B: Condens. Matter Mater. Phys. 1979, 20, 5104−5110. (15) Nield, V. M.; Howe, M. A.; McGreevy, R. L. The Metal-nonmetal Transition in Expanded Caesium. J. Phys.: Condens. Matter 1991, 3, 7519−7525. (16) Carpena, P.; Bernaola-Galvan, P.; Ivanov, P. Ch.; Stanley, H. E. Metal−insulator Transition in Chains with Correlated Disorder. Nature 2002, 418, 955−959. G

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

(41) Eslami, H.; Müller-Plathe, F. Water Permeability of Poly(ethylene terephthalate): A Grand Canonical Ensemble Molecular Dynamics Simulation Study. J. Chem. Phys. 2009, 131, 234904. (42) Eslami, H.; Müller-Plathe, F. Molecular Dynamics Simulation of Water Influence on Local Structure of Nanoconfined Polyamide-6,6. J. Phys. Chem. B 2011, 115, 9720−9731. (43) Eslami, H.; Kesik, M.; Karimi-Varzaneh, H. A.; Müller-Plathe, F. Sorption and Diffusion of Carbon Dioxide and Nitrogen in Poly(methyl methacrylate). J. Chem. Phys. 2013, 139, 124902. (44) Müller-Plathe, F. YASP: A Molecular Simulation Package. Comput. Phys. Commun. 1993, 78, 77−94. (45) Winter, R.; Hensel, F.; Bodensteiner, T.; Gläser, W. The Static Structure Factor of Cesium over the Whole Liquid Range up to the Critical Point. Ber. Bunsenges. Phys. Chem. 1987, 91, 1327−1330. (46) Yuryev, A. A.; Gelchinski, B. R. Ab initio Molecular Dynamics Study of Liquid Sodium and Cesium up to Critical Point. AIP Conf. Proc. 2013, 1673, 020009. (47) Matsuda, N.; Mori, H.; Hoshino, K.; Watabe, M. Universality in the Structural Change of Expanded Liquid Alkali Metals Along the Liquid-Vapour Coexistence Curve. J. Phys.: Condens. Matter 1991, 3, 827−836. (48) Kietzmann, A.; Redmer, R.; Hensel, F.; Desjarlais, M. P.; Mattsson, T. R. Structure of Expanded Fluid Rb and Cs: A Quantum Molecular Dynamics Study. J. Phys.: Condens. Matter 2006, 18, 5597− 5605. (49) Falconi, S.; Ackland, G. J. Ab initio Simulations in Liquid Caesium at High Pressure and Temperature. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 184204. (50) Jüngst, S.; Knuth, B.; Hensel, F. Observation of Singular Diameters in the Coexistence Curves of Metals. Phys. Rev. Lett. 1985, 55, 2160−2163. (51) Goldstein, R. E.; Ashcroft, N. W. Origin of the Singular Diameter in the Coexistence Curve of a Metal. Phys. Rev. Lett. 1985, 55, 2164. (52) Goldstein, R. E.; Parola, A.; Smith, A. P. Fluctuating Pseudoatoms in Metallic Fluids. J. Chem. Phys. 1989, 91, 1843−1854. (53) Kim, Y. C.; Fisher, M. E. Singular Coexistence-curve Diameters: Experiments and Simulations. Chem. Phys. Lett. 2005, 414, 185−192. (54) Groot, R.; Warren, P. B. Dissipative Particle Dynamics: Bridging the Gap between Atomistic and Mesoscopic Simulation. J. Chem. Phys. 1997, 107, 4423−4435. (55) Pranami, G.; Lamm, M. H. Estimating Error in Diffusion Coefficients Derived from Molecular Dynamics Simulations. J. Chem. Theory Comput. 2015, 11, 4586−4592. (56) Vogel, H. The Law of the Relationship between Viscosity of Liquids and the Temperature. Phys. Z. 1921, 22, 645−646. (57) Wax, J.-F.; Albaki, R.; Bretonnet, J.-L. Temperature Dependence of the Diffusion Coefficient in Liquid Alkali Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 65, 014301. (58) Sumi, T.; Miyoshi, E.; Tanaka, K. Molecular-dynamics study of liquid mercury in the density region between metal and nonmetal. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 6153−6158. (59) Eslami, H.; Khanjari, N.; Müller-Plathe, F. A Local Order Parameter-Based Method for Simulation of Free Energy Barriers in Crystal Nucleation. J. Chem. Theory Comput. 2017, 13, 1307−1316. (60) Khanjari, N.; Eslami, H.; Müller-Plathe, F. Adaptive-NumericalBias Metadynamics. J. Comput. Chem. 2017, 13, 1307−1316. (61) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Icosahedral Bond Orientational Order in Supercooled Liquids. Phys. Rev. Lett. 1981, 47, 1297−300.

H

DOI: 10.1021/acs.jced.7b00639 J. Chem. Eng. Data XXXX, XXX, XXX−XXX