Ionic Vapor Composition in Critical and Supercritical States of Strongly

Article pubs.acs.org/JPCB

Ionic Vapor Composition in Critical and Supercritical States of Strongly Interacting Ionic Compounds Vitaly V. Chaban† and Oleg V. Prezhdo*,‡ †

Instituto de Ciência e Tecnologia, Universidade Federal de São Paulo, 12231-280, São José dos Campos, SP, Brazil Department of Chemistry, University of Southern California, Los Angeles, California 90089, United States

‡

ABSTRACT: The critical point, CP (T, P), of the phase diagram quantifies the minimum amount of kinetic energy needed to prevent a substance from existing in a condensed phase. Therefore, the CP is closely related to the properties of the fluid far below the critical temperature. Approaches designed to predict thermophysical properties of a system necessarily aim to provide reliable estimates of the CP. Vice versa, CP estimation is impossible without knowledge of the vapor phase behavior. We report ab initio Born−Oppenheimer molecular dynamics (BOMD) simulations of sodium and potassium chlorides, NaCl and KCl, at and above their expected CPs. We advance the present knowledge regarding the existence of ionic species in the vapor phase by establishing significant percentages of atomic clusters: 29−30% in NaCl and 34−38% in KCl. A neutral pair of counterions is the most abundant cluster in the ionic vapors (ca. 35% of all vaporized ions exist in this form). Unexpectedly, an appreciable fraction of clusters is charged. The ionic vapor composition is determined by the vapor density, rather than the nature of the alkali ion. The previously suggested CPs of NaCl and KCl appear overestimated, based on the present simulations. The reported results offer essential insights into the ionic fluid properties and assist in development of thermodynamic theories. The ab initio BOMD method has been applied to investigate the vapor phase composition of an ionic fluid for the first time.

■

INTRODUCTION Ionic compounds (ICs) are omnipresent in chemistry. They are generally known for lower volatility and flammability and a larger liquid-state temperature range, as compared to molecular compounds.1−19 ICs can be inorganic, organic, and mixed. For instance, many room-temperature ionic liquids (RTILs)20−28 comprise an organic cation and inorganic anion, whereas inorganic crystals are essential for the Earth crust. Most ICs exhibit hydrophilic properties, to a larger or smaller extent. Furthermore, ICs are significantly miscible with each other, due to electrostatic forces. The latter feature permits tunability of the physicochemical properties, and, thus, opens up a wide avenue of possible applications. Purely inorganic ICs have high melting, boiling, and critical temperatures,2,3,6,14 and are used predominantly as solid materials. In turn, mixed and organic ICs can be low melting salts due to unfavorable ionic packing for steric and symmetry reasons. They serve as versatile solvents and tunable reaction media. Because of their organic nature, these ICs thermally decompose well below their critical and even normal boiling points. An accurate knowledge of the critical points (CPs), triple points, and all their parameters constitutes a starting point of many investigations aiming to predict thermophysical properties of ICs, most important of which are liquid state density and surface tension. The thermodynamic behavior of the majority of simple fluids is quite universal, provided that the corresponding state variables are expressed in the reduced form. This reduction is achieved by scaling all parameters by their respective values at the CP. Specific intermolecular and interionic interactions, such as hydrogen bonding, dipole− dipole association, and π-electron interactions between © XXXX American Chemical Society

relatively large species, complicate correlations of the thermophysical properties and give rise to deviations that are difficult to predict, although major trends still persist. Roman and co-workers29 published an encouraging theoretical advance, demonstrating a universal behavior of the thermodynamic properties of a wide class of pure fluids along the entire vapor−liquid coexistence curve. The CP and triple point parameters are used to convert surface tension (ST), enthalpy of vaporization (EV), and coexistence densities difference (CEDD) into the corresponding reduced values. Once reduced, the experimental data can be described by a relatively simple expression containing only two empirical parameters, which are specific for every thermodynamic property: the critical-point exponent and the slope at the triple point. Noteworthy, saturation pressure (SVP) does not obey the established universal behavior. SVP is conventionally fitted by the three-parameter Antoine equation, whereby the parameters are different below and above the normal boiling point. The analyzed fluids include noble gases, short-chain alkanes, molecular nitrogen, carbon dioxide, and aprotic polar fluids, which may be expected to exhibit more sophisticated trends. Hydrogen-bonded fluids, such as water, hydrofluoric acid, methanol, carboxylic acids, etc., were omitted from the analysis, together with other fluids exhibiting various types of specific site−site intermolecular interactions. Recently, Weiss30 published an analytic solution of the reformulated problem for selected molecular and ionic Received: March 7, 2016 Revised: April 19, 2016

A

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

temperature range. Simple empirical equations were utilized without guidance from either theory or statistical mechanical calculations considering atomistic and molecular parameters. Pitzer also predicted the CP for KCl based on analogy with NaCl. Any correlation method applied to predict the phase diagram in the vicinity of the CP requires reliable information regarding the vapor phase composition, because it affects the thermodynamic quantities used to establish the needed correlation. For instance, ion or molecule pairing in the vapor phase would decrease vapor pressure. Such pairing should depend on temperature. Deviations from the theoretical models are expected in the presence of a significant percentage of unaccounted vapor structures. For instance, NaCl is known to contain a substantial component of the ion pair dimer Na2Cl2, in addition to the monomer NaCl.32,33 While the partial pressure of Na2Cl2 is correct over the range in which it was determined experimentally, the extrapolation to higher temperatures is not certain. Gillan33 computed populations of various vapor species using a simple ionic model. The model used the direct Coulomb law at distances exceeding a predefined constant and an infinitely strong repulsion for smaller distances. Neither van der Waals forces nor electronic polarizability effects were accounted for. All ion pair dimers Na2Cl2 exhibited a ring structure in the experimental temperature range. However, there was a substantial population of open-chain structures at higher temperatures. Various structures were also proposed for Na3Cl3. For the first time, we report ab initio BOMD simulations of the vapor phases of ICs. In particular, we simulated NaCl and KCl at the CPs, hypothesized by Pitzer,32 and at a few supercritical conditions. We analyze cluster compositions, structure, dynamics, and electronic properties. There is still no experimental information on the supercritical states of these simple salts. Reliable atomically precise results are required for applications and for development and refinement of thermodynamic theories.

substances, including RTILs, namely, prediction of the CP based on the experimental triple point parameters and the experimental thermodynamic data at temperatures between the triple point and the CP. Compare, CP prediction from the data measured at water’s melting point, 273 K, amounts to 891 K (based on ST), to 686 K (based on EV), to 1887 K (based on CEDD). Measurements near the normal boiling point, at 370 K, allow for more successful predictions: 736 K (based on ST), 695 K (based on EV), 804 K (based on CEDD). The experimentally determined critical temperature for water is 647.1 K. One concludes that, the closer experimental results are to the CP, the more trustworthy is the established correlation. We suppose that the observed discrepancies occur largely due to the neglected microscopic changes in the liquid and vapor phases upon heating. Measurement of the CP parameters may be challenging for a significant fraction of the known fluids, both molecular and ionic. Unlike for gases and volatile solvents without hydrogen bonding, normally consisting of tiny particles, recording the upper part of the phase diagram for bulky organic molecules and most ICs is inaccessible for experimental techniques. One reason is an insufficiently high thermal decomposition temperature. Indeed, many organic substances decompose or oxidize prior to boiling. Another reason is a very high CP, a few thousands of degrees Celsius. Although certain compounds do not decompose upon such an extreme heating, e.g., alkali metals halides, working at high temperatures is technically problematic, furthermore subjecting the results to large uncertainties. Measurements of SVPs of ICs at room and moderately elevated temperatures are not efficient, since the temperature range is too small, giving rise to large error bars. The presently available experimental data for ICs are likely not sufficient for reliable estimation of the CPs. Emerging over the last decades, computer simulation methods, representing the systems at the atomistic or coarse-grained levels, with particles interacting via ab initio or carefully calibrated Hamiltonian, are able to overcome the above difficulties. Classical potential-based molecular dynamics (MD) allows one to prevent chemical reactions; therefore, a molecule of interest can be maintained as a stable particle at any thermodynamic conditions. Alternatively, density functional theory (DFT) based Born− Oppenheimer molecular dynamics (BOMD) allows one to carry out the simulations at very high temperatures including quantum mechanical effects at the ab initio level of description. In principle, as discussed above, the critical temperature can be obtained through extrapolation of certain thermodynamic properties over the region of higher temperatures. ST and EV reach zero, whereas liquid and vapor densities become identical at the critical temperature. Fifty years ago, Kirshenbaum and co-workers31 observed a phase coexistence curve for sodium chloride (NaCl) and potassium chloride (KCl). The authors extrapolated empirical vapor density equations and assumed that the law of rectilinear diameters (LRD) held. Recall, LRD states that the dependence of the mean density on temperature is linear, while the liquid density near the CP is reduced below the straight line by the amount equal to the vapor density. Accordingly, the density must be a smooth continuous curve. The CP can be located with a reasonably moderate uncertainty. The approach employed by Kirshenbaum and co-workers31 was critically reviewed by Pitzer,32 who concluded that the reported densities were clearly overestimated above 2000 K. According to Pitzer, these discrepancies were not surprising, because the corresponding estimates arose from extrapolations over a large

■

METHODOLOGY

The electronic structure calculations and adiabatic BOMD simulations were conducted with the Vienna Ab initio Simulation Package (VASP).34 A pure DFT method with a converged plane-wave basis set works well to generate wave functions for the three-dimensional periodic systems. Periodic systems in BOMD eliminate undesirable boundary effects, allowing a virtually infinite system to be represented. The simulated systems (Table 1) include 20 ions of either NaCl or Table 1. Simulated Systems, Their Composition, and Parametersa no. 1 2 3 4 5 6 7 8

no. ion pairs 10 10 10 10 10 10 20 10

NaCl NaCl NaCl KCl KCl KCl NaCl NaCl

temperature (K)

mass (amu)

box side (Å)

3900 4100 4300 3470 3670 3870 3900 2000

584.43 584.43 584.43 745.51 745.51 745.51 1168.86 584.43

20.07 20.07 20.07 21.77 21.77 21.77 25.29 20.07

a

Equilibrated high temperature geometries obtained from the BOMD simulations are exemplified in Figure 1.

B

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B KCl. One additional system of 40 ions was simulated to identify the possible size dependence of the investigated properties. The unit cells were constructed considering mass densities at the CP predicted by Pitzer, 120 kg m−3 for both vapors. For the sake of comparison, supercritical states were prepared with the same densities but at higher temperatures. NaCl was simulated at 2000 K to observe a vapor−liquid coexistence. The exchange-correlation functional proposed by Perdew, Burke, and Ernzerhof in the generalized gradient approximation was employed.35 The projector-augmented wave method was used for pseudopotentials to simulate core (nonvalence) electrons with a higher computational efficiency.36 The initial geometries were energy-minimized to eliminate nonphysical forces before BOMD. The nuclear equations of motion were propagated with an integration time step of 1.0 ps. Every system was simulated during 100.0 ps to record trajectories for the calculation of the equilibrium microscopic properties, such as radial distribution functions (RDFs), self-diffusion coefficients (SDCs), ionic vapor composition, and selected thermodynamic quantities. The cohesive energy, Ecoh, was computed as Ecoh = Epot − 10 × (E1an + E1cat), where E1an and E1cat are potential energies (all nucleus−electron interactions) of the isolated anion and cation, respectively. An atom−atom RDF shows by what factor the local density of a given atom type at a certain distance from the reference atom exceeds the average density, normalized to unity. SCD was computed from the dependence of the meansquared displacement (MSD) on time. An ionic cluster was defined as a collection of ions, in which no ion is separated from all other ions by more than 4.0 Å (position of the first RDF minimum). Thermochemical properties were computed using the Gaussian-4 (G4) composite theory,37 which combines a few levels of theory up to coupled-cluster calculations, to predict highly accurate thermodynamic potentials. According to Curtiss and co-workers,37 an average absolute deviation of the G4 results from experiment amounts to 0.83 kcal mol−1, obtained using 454 test cases. These results represent an improvement compared to the earlier theories.38 Only global-minimum structures of the ion clusters were used in these computations. Global geometry optimization was fulfilled by the basin hopping algorithm39 with the temperature parameter of 3900 K to overcome barriers on the potential energy surface, and the maximum step width of 0.5 Å. Each ionic cluster was subjected to 50 consequent iterations (local geometry optimizations). The lowest obtained potential energy was assigned to the global minimum. Local optimizations were performed using the B3LYP hybrid DFT40,41 in conjunction with the 6-31G* basis set. The self-consistence field convergence threshold was set to 10−6 hartree. The Nosé thermostat was used to maintain constant temperature.42 Packmol43 was used to prepare initial geometries. VMD44 and Gabedit45 codes were used for preparation and analysis of system snapshots.

Figure 2. Evolution of potential energy Epot, total energy Etot, and temperature T: NaCl at 3900 K (red solid line) and KCl at 3470 K (green dashed line).

cases, confirmed by block-averaging of thermodynamic quantities over equal time intervals, such as 0−2 ps, 2−4 ps, etc. This is in concordance with high ionic self-diffusion coefficients, as discussed below. After equilibration, the energies and temperatures fluctuate significantly due to small system sizes. In NaCl at 3900 K, the standard deviation of the system’s potential and kinetic energies amounts to 1.79 eV, whereas the standard deviation of temperature amounts to 201 K. The fluctuations are smaller in KCl at 3470 K, because the temperature is lower: 1.66 eV for both energies and 146 K for temperature. Cohesive energies slightly decrease with temperature (Figure 3). Due to smaller temperatures of simulation (3470−3870 K

Figure 3. Cohesive energy in (a) NaCl vapor and (b) KCl vapor versus temperature. The shown error bars were computed from the potential energy fluctuations at the respective temperature during trajectory.

vs 3900−4300 K), the cohesive energy appears more negative in KCl. At room conditions, NaCl features a stronger bond, due to a smaller radius and a larger covalent-bond admixture. Similarity of the cohesive energies in the vicinity of the CP of each salt is expected, since the CP is determined by a relation of kinetic energy to potential energy, per degree of freedom. Note that Figure 3 depicts Ecoh per single ion, while the cohesive energy used to characterize the crystal lattice is often computed with respect to an ion pair, i.e., Ecoh = Epot − n × Eion pair, where n stands for the number of ions in the supercell. The conventional definition of Ecoh is poorly applicable to an ionic vapor, since ions are disordered at this temperature. The values in Figure 3 reflect a potential energy gain due to formation of ion pairs and larger clusters. Compare, the Gibbs free energy

■

RESULTS AND DISCUSSION Critical points of the inorganic ionic melts, ca. 4000 K, appear well above their respective melting points, ca. 1000 K. Therefore, the corresponding liquid-state temperature range is very wide. Its applications are naturally limited due to technical complications with working at such high temperatures. A thermodynamic equilibration at these temperatures is extremely quick (Figure 2) taking less than 2 ps in all simulated C

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B for pairing Na+ and Cl− in the gas phase at 3900 K and 258 bar is −431 kJ mol−1. Since this value is over 3 times smaller than Ecoh, one should expect existence of more sophisticated ionic clusters in the vapor of both NaCl and KCl. Figure 4 depicts RDFs for all possible ion−ion combinations in NaCl and KCl. Table 2 summarizes the principal parameters

SDCs for all ions and average ionic SDCs in each system were computed from the slope of the MSD dependence on time (Figure 5). Because of a high temperature, the standard

Figure 4. Radial distribution functions computed between all possible pairs of atoms in (a) NaCl and (b) KCl: cation−anion (red solid line), cation−cation (green dashed line), anion−anion (blue dash-dotted line).

Figure 5. Translational self-diffusion coefficients, D, of the cations and anions in (a) NaCl and (b) KCl versus temperature; (c) mean-squared displacements (MSDs) versus simulated time and the fitting line, Dt=0 = 0, for KCl at 3470 K. The error bars were computed from three equal trajectory parts of the 100 ps long trajectory. No drift of D was observed, indicating proper equilibration of all systems.

of RDFs, including the position of the first peak and its height, at various temperatures. Due to a larger diameter of the potassium atom, its RDF peaks are located at somewhat larger distances, as compared to those in NaCl. Compare 2.5 Å (Na− Cl) to 2.9 Å (K−Cl) and 4.4 Å (Na−Na at 3900 K) to 5.1 Å (K−K at 3470 K). Chlorine−chlorine distances depend on the cation, since the cation separates anions (and vice versa) in the ionic clusters. An existence of the Cl−Cl peak suggests formation of ionic clusters that are larger than an ion pair. Otherwise, no peaks would be present in the gas phase. RDFs suggest that 3900 K (for NaCl) and 3470 K (for KCl) fall above the CPs for these salts, since no second peak is observed. The second peak would indicate a presence of the nascent liquid phase. Strictly speaking, it is impossible for the liquid phase to form in simulation cells that are much smaller than the thermodynamic limit. As we discuss below, it is still possible to hypothesize whether the simulated system is fully gaseous or not. Therefore, the original predictions of the CPs are somewhat overestimated.

deviations of SDCs are quite large, in a few cases exceeding the difference of SDCs at the neighboring temperatures. The SDCs are of the order of 10−6 m2 s−1, which is ca. 1000 larger than self-diffusion in most liquids at room temperature, ca. 10−9 m2 s−1.46 Transport of the gas particles is inversely proportional to their mass. This is confirmed by our results: SDC (Na+) > SDC (Cl−) in concordance with M (Na+) < M (Cl−) in NaCl. Ionic clusters existing in the ionic vapor determine vapor pressure, since vapor pressure is proportional to the concentration of particles in the vapor phase. Previous studies suggested predominant formation of ion pairs and, to a lesser extent, of neutral ion quartets with a modest admixture of neutral ion sextets. Our simulations, which account for electronic and quantum effects in the course of molecular dynamics, support this knowledge (Figure 6). In addition, they reveal a significant fraction of the charged ionic clusters comprising three, five, and seven ions. In the case of NaCl, the fraction of ions belonging to the charged cluster decays smoothly as the cluster size increases. In turn, the average number of KCl ions belonging to triplets and quintets is very

Table 2. Positions of Peaks and Their Heights in the Ion−Ion RDFs at Different Simulated Temperaturesa Na−Cl in NaCl

Na−Na in NaCl

T (K)

rmax (Å)

RDFmax

rmax (Å)

RDFmax

rmax (Å)

RDFmax

3900 4100 4300

2.5 ± 0.0 2.6 ± 0.0 2.5 ± 0.0

15.3 ± 0.3 15.0 ± 0.3 13.5 ± 0.2

4.4 ± 0.2 4.9 ± 0.2 4.5 ± 0.2

2.0 ± 0.1 1.7 ± 0.1 1.7 ± 0.1

4.5 ± 0.2 4.8 ± 0.2 5.0 ± 0.2

2.1 ± 0.1 1.9 ± 0.1 1.8 ± 0.1

K−Cl in KCl

a

Cl−Cl in NaCl

K−K in KCl

Cl−Cl in KCl

T (K)

rmax (Å)

RDFmax

rmax (Å)

RDFmax

rmax (Å)

RDFmax

3470 3670 3870

2.9 ± 0.0 2.9 ± 0.0 2.9 ± 0.0

14.7 ± 0.2 13.0 ± 0.2 13.0 ± 0.2

5.1 ± 0.2 5.3 ± 0.2 5.3 ± 0.2

1.9 ± 0.1 1.9 ± 0.1 1.7 ± 0.1

5.4 ± 0.2 5.4 ± 0.2 5.6 ± 0.2

1.8 ± 0.1 1.8 ± 0.1 1.7 ± 0.1

Standard uncertainties were estimated by considering three equal parts of each trajectory. D

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

To complete the characterization of the critical and supercritical vapor phases, we provide global-minimum ionic configurations (Figure 8) of the most abundant clusters in

Figure 8. Ideal (optimized) geometries of certain ionic clusters observed in the NaCl vapor: (a) Na2Cl; (b) NaCl2; (c) Na2Cl2; (d) Na2Cl3; (e) Na3Cl2. In all clusters, the bonds are predominantly ionic. The lines connecting atoms are drawn to guide the eye. Figure 6. Distribution of the vaporized ions among clusters of different size (one to eight atoms). The errors bars were computed from three equal trajectory parts of the 100 ps long trajectory.

NaCl. These configurations are quite different from those in the immediate snapshot of the equilibrated state (Figure 1). This is

similar, while heptets are observed rarely. Note that Figure 6 depicts percentages of ions existing as clusters of given sizes, rather than the percentages of these clusters themselves. Thus, the number of ions existing as triplets and quintets in KCl is similar, whereas triplets are more abundant, by ca. 5/3. Interestingly, the percentages of ions, existing as pairs (ca. 35%) and quartets (ca. 16%), are equal in both NaCl and KCl. KCl was simulated at a significantly lower temperature than NaCl. An evident hypothesis is that ion clusters are predominantly determined by vapor density, rather than the nature of the alkali cation and corresponding kinetic energies per degree of freedom. Five percent of ions are unpaired, and the same fraction of ions exist as large clusters (seven and eight ions). Larger clusters were also sporadically detected in the simulations, but the corresponding probabilities do not exceed 1%. Note that observation of large ionic clusters is prone to larger uncertainties. Longer simulations and large simulation cells are required to investigate them with a trustworthy accuracy. In concordance with the distribution of the ion cluster sizes, the average size, shown versus time in Figure 7, is ca. 2.5 ions, irrespective of the alkali cation and temperature. Interestingly, fluctuations of this property in the case of NaCl are higher, being in line with higher kinetic energies and higher translational speeds of the ions, because of the lower atomic mass, M (Na, 23 amu) < M (K, 39 amu).

Figure 1. Ionic configurations of the simulated critical states after equilibration: (a) 10 KCl at 3470 K and (b) 20 NaCl at 3900 K. Sodium is red, potassium is orange, and chlorine is blue.

expected, since high kinetic energies greatly perturb geometries, and the bond lengths increase as a result of thermal expansion. The Na−Cl distance in the ion pair is equal to 2.36 Å. The distance increases to 2.48 Å in the Cl−Na−Cl triplet and to 2.47 Å in the Na−Cl−Na triplet. The Na−Cl bond length continues to increase in the ion quartet, 2.54 Å, and in the ion quintet (Na2Cl3), 2.81 and 2.49 Å. The global-minimum ionic configuration of Na3Cl2 is qualitatively different from the configuration of the other clusers. Na3Cl2 is linear with the cation−anion distances of 2.44 and 2.53 Å. Our geometries agree with the results of Gillan,33 including linear structures and the ring structure of the quartet. The bond lengths are in excellent agreement with the RDFs (Figure 4 and Table 2). Further insight into stabilities of the ionic clusters is provided by the free energies of the cluster formation reactions taking place in the vapor phase of NaCl (Table 3). Expectedly, formation of the ion pair from the individual cation and anion is most thermodynamically favorable, −431 kJ mol−1. Formation of Na2Cl2 from NaCl2 and Na2Cl is somewhat less favorable. In turn, formation of the charged clusters is systematically less likely, whereas formation of the positively charged clusters, Nax+1Clx, is even prohibited. Nevertheless, a significant fraction of these clusters exist in the vapor, thanks to large kinetic energies involved at the critical state (16.2 kJ mol−1 per degree of freedom). The fractions of ions belonging to all possible charged clusters is systematically higher in the case of KCl: 34% (3470 K), 38% (3670 K), and 36% (3870 K), whereas the analogous fractions in NaCl are 30% (3900 K), 29% (4100 K), and 30%

Figure 7. Average size of ionic clusters in vapor versus simulated time at different temperatures: (a) NaCl; (b) KCl. E

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

much less sophisticated simulation approach, appear to be overestimations for both salts. It is probably due to neglected electronic polarization, which is able to profoundly modify the shape of the vapor pressure vs temperature curve, and consequently the phase diagram. Another likely reason is that Pitzer relied only on the neutral species in the ionic vapor, whereas the contribution of the charged clusters, such as NaCl2, Na2Cl3, and Na3Cl4, and single ions, should be recognized as important. Our methodology can be used to obtain more accurate estimations of the CPs by locating the temperature as which no signs of the liquid phase exist. It is impossible to represent a complete liquid phase in the simulation box because of the size restrictions. Figure 10 exemplifies cluster

Table 3. Formal Chemical Reactions for the Formation of Ionic Clusters and the Corresponding Gibbs Free Energies, ΔG, for Sodium Chloride at 3900 K and 258 bara reactants

product

ΔG (kJ mol−1)

product charge (|e|)

Na + Cl NaCl + Na NaCl + Cl Na2Cl + Cl NaCl2 + Na Na2Cl2 + Na Na2Cl2 + Cl

NaCl Na2Cl NaCl2 Na2Cl2 Na2Cl2 Na3Cl2 Na2Cl3

−431 +5 −44 −417 −369 +22 −46

0 +1 −1 0 0 +1 −1

a

The calculations were conducted in accordance with the Gaussian-4 composite method for accurate thermochemistry. Note that the average thermal energy carried by each degree of freedom of the ionic clusters at 3900 K is 16.2 kJ mol−1.

(4300 K). No dependence on temperature was found, although one would expect that higher thermal energy fosters a larger concentration of the less thermodynamically stable formations, especially those of lone ions. Indeed, Figure 6 suggests a slight increase of the single ions at higher temperatures in both salts. We anticipate small systems (20 ions) to be sufficient for the investigation of the supercritical states. First, the estimated mass densities of the systems at TC are relatively small, 120 kg m−3, resulting in simulation cells of significant size, over 20 Å per cubic side (Table 1). Second, the number of ions in the simulation cell has to exceed the largest cluster size. As our analysis indicates, the largest observed cluster has over twice fewer atoms than the total number of ions in the system. However, a poorly predictable uncertainty in the results can be caused by large temperature and pressure fluctuations, which are inevitable for relatively small systems. To identify this and other size dependent effects, we simulated a twice larger BOMD system, involving 20 NaCl ion pairs (40 ions) (Table 1). The comparison of the ion cluster distribution in the larger and smaller systems is provided in Figure 9. The observed

Figure 10. Cluster size distribution in NaCl at 2000 K. The provided NaCl density is an average quantity computed as the mass of all ions divided by the volume of the simulation box, irrespective of the phase.

composition of NaCl at 2000 K, suggesting formation of the liquid phase (the largest ionic cluster, 20 ions) and sporadic emergence of the vapor species. One concludes that 2000 K falls below the CP of NaCl, as expected. Although the illustrated methodology may exhibit imprecisions in the vicinity of the CP due to small system size and inability to differentiate the nascent liquid phase from the large ionic clusters, it is useful to assess the CP for the strongly interacting ionic substances. Pitzer provides the minimum and maximum vapor densities separately for the NaCl, Na2Cl2, and Na3Cl3 clusters, recomputed from different data.32 These densities correspond to the vapor−liquid equilibrium, i.e., below the CP. Figure 11

Figure 9. Distribution of cluster sizes in smaller (10 ion pairs) and larger (20 ion pairs) NaCl systems at 3900 K.

differences are small and do not reveal any systematic bias. No clusters larger than eight ions were detected in either simulation. Therefore, the originally selected system size is sufficient for the investigated properties, of which the cluster size distribution is most informative. The CP has a clear thermodynamic meaning. At its temperature, the kinetic energy of the system exceeds the maximum possible interatomic potential energy, per degree of freedom. As a result, the system cannot be returned to the liquid state, despite an external pressure applied. On the basis of the cluster analysis, the simulated systems are in their supercritical states at 3900−4300 K (NaCl) and 3470−3870 K (KCl). Thus, the predictions of Pitzer,32 performed using a

Figure 11. Partial vapor densities of the neutral ion clusters in the gaseous NaCl. The data of Pitzer32 are given in color, whereas our data are designated by black symbols.

summarizes these vapor densities and compares them to the supecritical state of NaCl at 3900 K. The vapor densities in the present work were computed from the vapor composition, as discussed above. Our data are in general agreement with Pitzer’s results, with the density of ion pairs being somewhat underestimated and the density of ion sextets being overestimated. Note that Pitzer denotes the density of hextets as maximum, while we observe a smaller actual density. Our result for Na2Cl2 is between Pitzer’s maximum and minimum values. F

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

■

(9) Chaban, V. V.; Verspeek, B.; Khandelia, H. Novel Ultrathin Membranes Composed of Organic Ions. J. Phys. Chem. Lett. 2013, 4, 1216−1220. (10) Araque, J. C.; Hettige, J. J.; Margulis, C. J. Modern Room Temperature Ionic Liquids, a Simple Guide to Understanding Their Structure and How It May Relate to Dynamics. J. Phys. Chem. B 2015, 119, 12727−12740. (11) Hettige, J. J.; Kashyap, H. K.; Annapureddy, H. V. R.; Margulis, C. J. Anions, the Reporters of Structure in Ionic Liquids. J. Phys. Chem. Lett. 2013, 4, 105−110. (12) Chaban, V. V.; Voroshyloya, I. V.; Kalugin, O. N.; Prezhdo, O. V. Acetonitrile Boosts Conductivity of Imidazolium Ionic Liquids. J. Phys. Chem. B 2012, 116, 7719−7727. (13) Camp, P. J.; Patey, G. N. Ion association and condensation in primitive models of electrolyte solutions. J. Chem. Phys. 1999, 111, 9000−9008. (14) Panagiotopoulos, A. Z. Simulations of phase transitions in ionic systems. J. Phys.: Condens. Matter 2005, 17, S3205−S3213. (15) Pitzer, K. S. Aqueous electrolytes at near-critical and supercritical temperatures. Int. J. Thermophys. 1998, 19, 355−366. (16) Chaban, V. V.; Pereverzev, Y. V.; Prezhdo, O. V. A new model of chemical bonding in ionic melts. J. Chem. Phys. 2012, 136, 164112. (17) Sabir, A. K.; Bhuiyan, L. B.; Outhwaite, C. W. Influence of ion size and valence on classical ionic criticality. Mol. Phys. 1998, 93, 405− 409. (18) Koulinskii, V. L.; Malomuzh, N. P.; Tolpekin, V. A. Influence of charge fluctuations on the critical behavior of electrolyte solutions. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 60, 6897−6905. (19) Fuentevilla, D. A.; Sengers, J. V.; Anisimov, M. A. Critical Locus of Aqueous Solutions of Sodium Chloride Revisited. Int. J. Thermophys. 2012, 33, 943−958. (20) Yuvaraj, S. V. J.; Zhdanov, R. K.; Belosludov, R. V.; Belosludov, V. R.; Subbotin, O. S.; Kanie, K.; Funaki, K.; Muramatsu, A.; Nakamura, T.; Kawazoe, Y. Solvation Mechanism of Task-Specific Ionic Liquids in Water: A Combined Investigation Using Classical Molecular Dynamics and Density Functional Theory. J. Phys. Chem. B 2015, 119, 12894−12904. (21) Yasaka, Y.; Kimura, Y. Polarity and Nonpolarity of Ionic Liquids Viewed from the Rotational Dynamics of Carbon Monoxide. J. Phys. Chem. B 2015, 119, 15493−15501. (22) Wu, B. N.; Liang, M.; Maroncelli, M.; Castner, E. W. Photoinduced Bimolecular Electron Transfer from Cyano Anions in Ionic Liquids. J. Phys. Chem. B 2015, 119, 14790−14799. (23) So, S.; Yao, L. J.; Lodge, T. P. Permeability of Rubbery and Glassy Membranes of Ionic Liquid Filled Polymersome Nanoreactors in Water. J. Phys. Chem. B 2015, 119, 15054−15062. (24) Shimizu, Y.; Fujii, K.; Imanari, M.; Nishikawa, K. Phase Behavior of a Piperidinium-Based Room-Temperature Ionic Liquid Exhibiting Scanning Rate Dependence. J. Phys. Chem. B 2015, 119, 12552− 12560. (25) Shaikh, A. R.; Kamio, E.; Takaba, H.; Matsuyama, H. Effects of Water Concentration on the Free Volume of Amino Acid Ionic Liquids Investigated by Molecular Dynamics Simulations. J. Phys. Chem. B 2015, 119, 263−273. (26) Sahu, P. K.; Ghosh, A.; Sarkar, M. Understanding StructureProperty Correlation in Monocationic and Dicationic Ionic Liquids through Combined Fluorescence and Pulsed-Field Gradient (PFG) and Relaxation NMR Experiments. J. Phys. Chem. B 2015, 119, 14221−14235. (27) Prabhu, S. R.; Dutt, G. B. Does Addition of an Electrolyte Influence the Rotational Diffusion of Nondipolar Solutes in a Protic Ionic Liquid? J. Phys. Chem. B 2015, 119, 6311−6316. (28) del Valle, J. C.; Blanco, F. G.; Catalan, J. Empirical Parameters for Solvent Acidity, Basicity, Dipolarity, and Polarizability of the Ionic Liquids [BMIM][BF4] and [BMIM][PF6]. J. Phys. Chem. B 2015, 119, 4683−4692.

CONCLUSIONS Ab initio Born−Oppenheimer simulations of NaCl and KCl at critical and supercritical conditions were reported over a significant time scale. Thermodynamic, structural, and dynamic properties, partial charge distributions, and vapor phase compositions were characterized. The vapor phase of the alkali metal halide ionic fluids consists of a variety of ionic species, ranging from the single ions to the neutral ionic octets. The previous knowledge about domination of ion pairs (primary component), quartets (to a less extent), and sextets (modest fraction) was supported. In addition, a significant percentage of the charged ionic clusters was detected. This information is important for development and verification of phenomenological theories, which predict phase diagrams at high temperatures. The ionic vapor composition at the critical and supercritical states of NaCl and KCl was found to be determined by the vapor density, rather than the nature of the alkali ions. In turn, the diffusion coefficients are in inverse proportion to the ionic mass. Small systems show signs of vapor condensation at sufficiently low temperatures, although no liquid phase could be represented in the simulation cell due to size restrictions. It may be fruitful to attempt applying such systems to predict CPs of simple substances directly.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: (213) 821-3116. Fax: (213) 740-270. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are grateful to Andrew Sifain for comments on the manuscript. O.V.P. acknowledges financial support of the US Department of Energy, Grant No. DE-SC0014429.

■

REFERENCES

(1) Panagiotopoulos, A. Z. Molecular simulation of phase-equilibria simple, ionic and polymeric fluids. Fluid Phase Equilib. 1992, 76, 97− 112. (2) Stell, G. Critical-behavior of ionic-fluid models. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 45, 7628−7631. (3) Panagiotopoulos, A. Z. Monte-Carlo simulations of phase coexistence for polymeric and ionic fluids. Fluid Phase Equilib. 1995, 104, 185−194. (4) Tkachev, N. C. Analysis of liquid-vapour phase transitions for restricted primitive model of ionic system in mean spherical approximation. Phys. Chem. Liq. 1996, 31, 71−81. (5) Chaban, V. V.; Prezhdo, O. V. Ionic and Molecular Liquids: Working Together for Robust Engineering. J. Phys. Chem. Lett. 2013, 4, 1423−1431. (6) Hynninen, A.-P.; Panagiotopoulos, A. Z. Simulations of phase transitions and free energies for ionic systems. Mol. Phys. 2008, 106, 2039−2051. (7) Araque, J. C.; Yadav, S. K.; Shadeck, M.; Maroncelli, M.; Margulis, C. J. How Is Diffusion of Neutral and Charged Tracers Related to the Structure and Dynamics of a Room-Temperature Ionic Liquid? Large Deviations from Stokes-Einstein Behavior Explained. J. Phys. Chem. B 2015, 119, 7015−7029. (8) Roy, D.; Patel, N.; Conte, S.; Maroncelli, M. Dynamics in an Idealized Ionic Liquid Model. J. Phys. Chem. B 2010, 114, 8410−8424. G

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (29) Roman, F. L.; White, J. A.; Velasco, S.; Mulero, A. On the universal behavior of some thermodynamic properties along the whole liquid-vapor coexistence curve. J. Chem. Phys. 2005, 123, 124512. (30) Weiss, V. C. Predicting critical temperatures of ionic and nonionic fluids from thermophysical data obtained near the melting point. J. Chem. Phys. 2015, 143, 144503. (31) Kirshenbaum, A. D.; Cahill, J. A.; McGonigal, P. J.; Grosse, A. V. The Density of Liquid NaCl and KCl and an Estimate of Their Critical Constants Together with Those of the Other Alkali Halides. J. Inorg. Nucl. Chem. 1962, 24, 1287−1296. (32) Pitzer, K. S. Critical-Point and Vapor-Pressure of Ionic Fluids Including NaCl and KCl. Chem. Phys. Lett. 1984, 105, 484−489. (33) Gillan, M. J. Liquid-Vapor-Equilibrium in the Restricted Primitive Model for Ionic Liquids. Mol. Phys. 1983, 49, 421−442. (34) Kresse, G.; Furthmuller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 1996, 6, 15−50. (35) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (36) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758−1775. (37) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 theory using reduced order perturbation theory. J. Chem. Phys. 2007, 127, 124105. (38) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K.; Rassolov, V.; Pople, J. A. Gaussian-3 theory using reduced M?ller-Plesset order. J. Chem. Phys. 1999, 110, 4703−4709. (39) Wales, D. J.; Doye, J. P. K. Global optimization by basinhopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms. J. Phys. Chem. A 1997, 101, 5111−5116. (40) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic-Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (41) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the ElectronDensity. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (42) Nose, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255−268. (43) Martinez, L.; Andrade, R.; Birgin, E. G.; Martinez, J. M. PACKMOL: a package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157−2164. (44) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (45) Allouche, A. R. Gabedit-A Graphical User Interface for Computational Chemistry Softwares. J. Comput. Chem. 2011, 32, 174−182. (46) CRC Handbook of Chemistry and Physics, 96th ed.; Taylor and Francis Group: Abingdon, 2015−2016; http://www.hbcpnetbase. com/.

H

DOI: 10.1021/acs.jpcb.6b02405 J. Phys. Chem. B XXXX, XXX, XXX−XXX