PCCP View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 3062
View Journal | View Issue
Absolute thermodynamic properties of molten salts using the two-phase thermodynamic (2PT) superpositioning method Jin Wang, Brahmananda Chakraborty and Jacob Eapen* We show that the absolute thermodynamic properties of molten salts (mixtures of KCl and LiCl) can be accurately determined from the two-phase thermodynamic (2PT) method that is based on superpositioning of solid-like and gas-like (hard-sphere) vibrational density of states (DoS). The 2PT predictions are in
Received 25th June 2013, Accepted 28th October 2013
excellent accordance with those from the thermodynamic integration method; the melting point of KCl evaluated from the free energy and the absolute entropy shows close conformity with the experimental/ NIST data. The DoS partitioning shows that the Li+ ions in the eutectic LiCl–KCl molten mixture are largely
DOI: 10.1039/c3cp52632a
solid-like, unlike the K+ and Cl ions, which have a significant gas-like contribution, for temperatures
www.rsc.org/pccp
employed for chemical and nuclear reprocessing applications.
ranging from 773 K to 1300 K. The solid-like states of the Li+ ions may have practical implications when
1. Introduction There is an engaging interest in determining the absolute thermodynamic properties such as free energy, enthalpy, entropy and specific heat in the liquid and disordered states.1–7 The term absolute is used in the sense that the thermodynamic properties can be evaluated directly through the partition function for a given state in the canonical ensemble; there is however, an arbitrary reference energy which is associated with the potential energy of the system. Knowing the absolute properties is particularly useful for optimizing the back-end of the nuclear fuel cycle, which entails the crucial step of separation of actinides from the spent nuclear fuel, followed by the development of new fuel through recycling. The advantage of having a closed fuel cycle includes a reduction of high-level nuclear waste and the optimal utilization of the fissionable isotopes. While aqueousbased (wet) processes are well-established there is a current interest for high temperature pyrometallurgical processes8–10 that utilize molten salts such as LiCl or KCl. A eutectic mixture of LiCl and KCl is particularly favored because of its lower melting point (626 K) relative to that of the constituent salts (883 K for LiCl and 1043 K for KCl). Originally developed for spent metallic fuel from the fast reactors, pyro processes can also be used for spent fuel from the current light water reactors (LWRs) using electroreduction, electrorefining, and pyropartitioning; they are also regarded to be more proliferation-resistant than the traditional wet processes. Department of Nuclear Engineering, North Carolina State University, Raleigh, NC 27695, USA. E-mail: [email protected]; Tel: +1 919-515-5952
3062 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
Reprocessing of nuclear fuel involves a large number of fissionable actinides8,11 as well as non-fissionable fission products; the separation processes typically depend on the chemical activity of the species, which in turn depends on the relative free energy.10 Thermodynamic integration (TI), which is the method of choice for a variety of applications,3,6,12–15 especially with molecular dynamics (MD) simulations, may not be the most optimal choice for nuclear reprocessing applications given the inordinate number of species that are targeted for separation. In several instances, TI can be limited by the optimal choice of the integration path that is not known a priori, and is also constrained in having a reverse path that implicitly excludes a first order phase change. Another challenge arises in proscribing appropriate reference systems, which can allow computation of absolute thermodynamic properties. While TI can successfully predict free energy differentials with actinide transmutation,10 integration paths typically require trajectories that span over 1 ns, which is expensive by current standards. A problem of a different nature may arise with actinide dissolution in molten salts – depending on the volume fraction, the interaction strength and the thermodynamic state, the dissolved species may form glassy states or even undergo a first-order liquid–liquid phase transition,16–20 which is somewhat poorly understood currently. Very strong attractive interactions have been shown to cause an amorphous transition in a model nanocolloidal system with a volume fraction as low as 5%.21 Thus it would be profitable to identify a methodology that can obviate some or all of the aforesaid possible drawbacks for systems with large number of species in complex liquid forms. Originally developed by Lin, Blanco and Goddard,22 the twophase thermodynamic (2PT) model has recently gained wide
This journal is © the Owner Societies 2014
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
Paper
PCCP
attention and success in predicting absolute thermodynamic properties for several types of fluids.1,2,4,5,7,23–26 Motivated by the approach of evaluating the partition function in the solid state using idealized harmonic oscillators, the central idea behind the 2PT methodology is to construct the partition function of a complex liquid state through a superposition of solid-like and gas-like vibrational density of states (DoS). An exact representation of the idealized states can be derived by treating the gas-like state through hard-sphere interactions with weighting functions provided by the Carnahan–Starling equation of state,27 and assuming the solid-like state to be a collection of independent harmonic oscillators. The thermodynamic properties of a liquid system thus can be shown to be a combination of the properties of a set of harmonic oscillators (solid-like) and a set of hard spheres (gas-like). Since the 2PT method thus far has not been employed to evaluate the properties of molten salts and mixtures, its efficacy or accuracy is currently unknown. In this investigation, we show that the absolute thermodynamic properties of LiCl/KCl molten salts can be accurately determined from the 2PT method. The 2PT free energy is in excellent agreement with that determined from the thermodynamic integration (TI) method, with relative errors of 1% or less, typically. The melting point of ionic KCl (1025 K) compares favorably with the experimental data (1041 K); the entropy of the melt phase also shows good conformity with the recommended data from NIST. We then determine the absolute free energy and entropy which are currently unknown for a eutectic mixture of LiCl and KCl for temperatures ranging from 773 K to 1300 K. From the partitioning of the DoS, which is uniquely possible in the 2PT method, we show that the Li+ ions in the LiCl–KCl mixture are dominantly solid-like even at temperatures as high as 1300 K. For chemical and nuclear applications the solid-like states can potentially impose practical limits, for example, to the amount of dissolved species that can be practically accommodated during reprocessing.
2. The 2PT methodology The thermodynamic properties in the canonical ensemble can be evaluated from the canonical partition function (Q); for example, the internal energy (E), entropy (S), specific heat at constant volume (Cv), and the Helmoholtz free energy (A) can be evaluated as28 @ ðln QÞ (1a) E ¼ kB T 2 @T @ S ¼ kB T ðln QÞ þ ln Q @T
Cv ¼
@ @ kB T 2 ðln QÞ @T @T A = kBT(ln Q)
This journal is © the Owner Societies 2014
(1b)
(1c)
(1d)
where kB is the Boltzmann constant and T is the absolute temperature. The thermodynamic properties are seldom evaluated through the partition function as it is known exactly only for idealized systems such as harmonic oscillators and ideal gas. Using a normal mode analysis, a solid state can be approximated as a system of non-interacting harmonic oscillators (HO); for this case, the total canonical partition function (QHO) can be expressed as29 QHO ¼
3N Y
qi
(2)
i¼1
where qi is the partition function of the ith harmonic oscillator mode, and N is the total number of oscillators. The total partition function can be shown to be related to the density of states as29 ð1 ln QHO ¼ dnGðnÞ ln½qðnÞ (3) 0
where G(n) is the density of states, which is given by22 ð 4 1 GðnÞ ¼ dtCðtÞei2pnt kB T 0
(4)
where C(t) is the mass-weighted velocity autocorrelation function,28 which can be directly computed from a molecular dynamics (MD) simulation (ab initio or classical). Since an exact quantum-mechanical expression is known for q(n), the total partition function can be evaluated using eqn (3) and (4), and the thermodynamic properties such as given by eqn (1a) to (1d) can be assessed for the solid states with reasonable accuracy.29 For a liquid state, the system cannot be considered as a collection of harmonic oscillators and a direct application of the above approximation will lead to non-physical results. The 2PT method22 extends the above idea to the liquid state by assuming a superposition of two idealized states – hard spheres (HS) that correspond to a gas-like state, and harmonic oscillators that correspond to a solid-like state. First it is hypothesized that the total density of states of the system that is now regarded to be a non-interacting mixture of hard spheres and harmonic oscillators can be decomposed as HO Gk(n) = GHS k (n) + Gk (n)
(5)
where k denotes the different species in the mixture. A key step in the 2PT method is in defining and evaluating a fluidicity parameter ( f ) given by the following expression for partitioning the system into the idealized states. Ð1 dnGHS 3NkHS NkHS k ðnÞ fk ¼ Ð01 ¼ (6) ¼ 3Nk Nk 0 dnGk ðnÞ where NHS k is the effective number of hard sphere atoms or ions of the kth component of the mixture. Note that the integral of the total density of states over all frequencies is simply 3Nk. Thus the 3Nk degrees of freedom of the kth component is assumed to be composed of 3NHS or equivalently, 3fkNk hard k sphere or diffusive degrees of freedom, and 3(Nk NHS k ) or equivalently, 3Nk(1 fk) solid-like or non-diffusive degrees of freedom.25 With this approximation, the thermodynamic
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3063
View Article Online
PCCP
Paper
h i 1 sHS ðnÞ HS HS WAHS k ðnÞ ¼ W ðnÞ W ðnÞ ¼ Ek Sk k 2 3kB
properties are formally expressed as25 Ek ¼ Ek0 þ kB T
ð 1 0
HS dnGHS k ðnÞWEk ðnÞ þ
ð1
HO dnGHO k ðnÞWE ðnÞ
0
(7a) Sk ¼ kB
ð 1
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
0
Cv;k ¼ kB
HS dnGHS k ðn ÞWSk ðn Þ þ
ð 1 0
ð1 0
HO dnGHO k ðn ÞWS ðn Þ
HS dnGHS k ðnÞWCv;k ðnÞ þ
ð1 0
(7b)
HO dnGHO k ðnÞWCv ðnÞ
(7c) Ak ¼
Ek0
þ kB T
ð 1 0
HS dnGHS k ðnÞWAk ðnÞ
ð1
þ
0
HO dnGHO k ðnÞWA ðnÞ
(7d) where W is a weight function for the appropriate property, and E0k is the reference energy which is given by ref. 25 fk 0 MD Ek ¼ Ek 3Nk kB T 1 (8) 2 is the total energy from MD simulations for the kth where EMD k component. Thus each extensive thermodynamic property, for each mixture component, comprises of two contributions – one from diffusive (hard spheres) degrees of freedom, and the other from vibrating (harmonic oscillators) degrees of freedom. The distinguishing part of the 2PT method is in deriving an exact expression for GHS k given by ref. 22. " #1 pGk ð0Þn 2 HS Gk ðnÞ ¼ Gk ð0Þ 1 þ (9) 6fk Nk Thus the density of states for the hard sphere phase can be calculated as a function of the fluidicity parameter and Gk(0) which is just the zero-frequency value of the system (total) density of states. MD simulations are employed to calculate Gk(n) from the velocity autocorrelation function. With a known GHS k , the density of states of the HO phase is then evaluated as HS GHO k (n) = [Gk(n) Gk (n)]. The weight functions are known exactly for the HO and HS idealized states, and are given by5,25 WEHO ¼ WSHO
bhn bhn þ bhn 2 e 1
bhn ln 1 ebhn ¼ bhn e 1 WCHO ¼ v
ðbhnÞ2 ebhn ðebhn
1Þ
1 ebhn WAHO ðnÞ ¼ ln ebhn=2
WSHS ðnÞ ¼ k
1 2
sHS k ðnÞ 3kB
3064 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
where b1 = kBT, h is the Planck constant and sHS is the entropy of the hard sphere phase (per atom/ion), which in the 2PT model is evaluated as ( ) yk 4Þ y^k ð3^ HS IG sk ¼ sk þ kB ln½zk ðy^k Þ þ (12) ð1 y^k Þ2 The first term of the above equation denotes the ideal gas entropy, however, weighted only for the HS phase. It is given by ref. 5: ( " #) 5 2pmk kB T 3=2 Vk þ ln (13) ¼ k sIG B k 2 h2 f k Nk where mk and Nk are the mass and the total number of ions/ atoms of the kth species, respectively. The term zk PVk/(NkkBT) in eqn (12) delineates the compressibility of the hard sphere phase, which is determined from the Carnahan–Starling equation given by zk ðy^k Þ ¼
1 þ y^k þ y^k2 y^k3
(10b)
(10d) (11a) (11b)
(14)
ð1 y^k Þ3
where ˆyk is the weighted hard sphere packing fraction (by the fluidicity parameter) for the kth component, and is given by ˆyk = fkyk, where yk denotes the packing fraction which is expressed as 3 prk sHS k yk ¼ (15) 6 is the hard sphere where rk is the number density and sHS k diameter of the kth species that in turn influences the fluidicity parameter (fk). In the 2PT methodology it is assumed that the fluidicity parameter is simply a ratio of the self-diffusivity (Dk) to the hard sphere diffusivity at zero pressure (D0HS ), which can k be written as f ¼
DðT; PÞ HS D0HS k ðT; P; s Þ
(16)
The self-diffusivity can be evaluated from MD simulations as the time integral of the velocity autocorrelation, and is expressible as Dk ðT; PÞ ¼
(10c)
2
WEHS ðnÞ ¼ WCHS ðnÞ ¼ k v;k
(10a)
(11c)
kB TGk ð0Þ 12mk Nk
(17)
where Gk(0) is evaluated from the expression ð 4 1 dtCk ðtÞei2pnt Gk ðnÞ ¼ kB T 0
(18)
We invoke the stationary principle of classical correlation functions,30 and thus only the real part of the Fourier transform needs to be used for calculating the above integral. Note that C(t) is the mass-weighted velocity autocorrelation function defined as Nk ð 1 Nk ð 1 X X Ck ðtÞ ¼ mk dthvk ðtÞ:vk ð0Þi ¼ mk dthck ðtÞi (19) i¼1
0
i¼1
0
This journal is © the Owner Societies 2014
View Article Online
Paper
PCCP
where ck(t) is the velocity autocorrelation function of the kth component. The hard sphere diffusivity in the limit of zero pressure/density is computed exactly from the Chapman– Enskog’s theory, which is given by 3 kB T 1=2 HS D0HS ¼ T; P; s (20) k k 2 pm k 8r sHS
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
k
k
The velocity autocorrelation for HS phase decays exponentially in time and hence the HS diffusivity can be expressed as22 DHS k ðT; fk rk Þ ¼
kB TGk ð0Þ 12mk fk Nk
(21)
The above expression can be related to the zero-pressure hard sphere diffusivity through the following relationship 0HS T; fk rk ; sHS DHS k ðT; fk rk Þ ¼ Dk k
4^ yk ½zðy^k Þ 1
(22)
where ˆyk = fky. After a few algebraic steps the following relationship can be established25 2Dk9/2fk15/2 6Dk3fk5 Dk3/2fk7/2 + 6Dk3/2fk5/2 + 2fk 2 = 0 (23) where 2Gk ð0Þ pkB T 1=2 Nk 1=3 6 2=3 Dk ¼ 9Nk mk p Vk
where sk and xk are the ionic radius and mole fraction, respectively, of the kth component. The Kirkwood and Buff (KB) theory31 gives a better approximation for partial volumes but it entails the evaluation of somewhat ill-convergent KB integrals. We have used a fixed ionic radius of 0.152 nm, 0.167 nm and 0.09 nm for K+, Cl and Li+ ions, respectively.32 As shown later, the ionic radius approximation works quite well for the current application. We have also observed that the calculated thermodynamic properties are not sensitive to moderate changes in the ionic radii. Finally, the molar (m) thermodynamic properties of the mixture are determined as xk Ekm
k¼1
m A ¼
n X
xk Skm R^
n X
xk lnðxk Þ
(26b)
k¼1
^ x k Am k þ RT
k¼1
n X
xk lnðxk Þ
(26c)
k¼1
ˆ is the universal gas constant, and n is the total number where R of components in the mixture.
3. Molecular dynamics simulations MD simulations are principally employed for calculating the velocity autocorrelation function. The molecular system typically consists of 1000 ions interacting through a long-ranged electrostatic potential and a Born–Huggins–Mayer (BHM) short-ranged potential. Originally developed by Fumi and Tosi,33 the rigid-ion potential is able to predict several thermodynamic, structural and transport properties of molten LiCl–KCl mixtures34,35 and other alkali halides33 with reasonable accuracy. The effect of polarization, however, is important for multicomponent systems,36–39 although it tends to be small for monovalent systems. The functional form of short-ranged part of the BHM potential is given by
(24)
k¼1
n X
n X
Fshort ¼ Bij earij ij
Eqn (24) contains the details of the physical state as well as the zero-frequency density of states that can be evaluated from MD simulations. Eqn (23) now can be iteratively solved to get the fluidicity parameter (fk) for each component of the mixture. Once f is known, the density of states and the weight functions of the HS phase can be computed, and the HO density of states can be determined as GHO = (Gk GHS k k ). Thus all the properties given by eqn (7a) to (7d) can be determined, for each component of the mixture. In this investigation, the partial volume (Vk) in eqn (24) is calculated from the relative ionic radii, which is given by sk3 V (25) Vk ¼ n P N 3 xk sk
m
m S ¼
Cij Dij rij6 rij8
(27)
where rij is the interionic separation between two ions, i and j, and a, B, C and D are constants.34,35 While the first term corresponds to the electron cloud repulsion, the second and third terms account for the dipole–dipole and the dipole–quadrupole dispersion interactions, respectively. The MD simulations are performed using the LAMMPS software40 with periodic boundary conditions at zero pressure using a timestep of 1 fs. Particle–particle–particle– mesh is used for columbic interactions, and the temperature and ´–Hoover algorithm. From pressure are controlled using the Nose the velocity data collected from MD simulations, the velocity autocorrelation function (VACF) is constructed, typically with a correlation time of 25 ps, using an overlapped data structure.41 From the Fourier transform of the VACF, the 2PT thermodynamic properties are computed using an in-house computer program.
4. Results and discussion We will first discuss (1) the numerical convergence of the fluidicity parameter with correlation time, followed by (2) a benchmark investigation of the 2PT method against published results from thermodynamic integration (TI) method for KCl, (3) computation of the molar free energy and the entropy change across melting in KCl along with a comparison to NIST data,42 and (4) an assessment of the density of states and thermodynamic properties of a eutectic mixture of LiCl–KCl for different temperatures. 4.1
Numerical convergence of the fluidicity parameter
k¼1
The fluidicity parameter for each component of the mixture is iteratively solved using eqn (23). Further the partitioning of the
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3065
E ¼
(26a)
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
PCCP
Paper
Fig. 1 Relative error in reconstructing the fluidicity parameter for the K+ ions in molten KCl.
system into solid-like and gas-like states depends on resolving the zero-frequency value of the density of states. To assess the optimal correlation time (for VACF), we have first computed fk using eqn (23) and then compared with the reconstructed ˆfk from the equation Ð1 ð HS 1 1 0 dnGk ðnÞ ^ Ð dnGHS ¼ fk ¼ 1 k ðnÞ 3Nk 0 0 dnGk ðnÞ
(28)
Insufficient sampling and numerical errors in the evaluation of VACF and the density of states will be reflected as a difference between ˆfk and fk. Fig. 1 shows the error in the reconstruction, expressed as a percentage, for the K+ ions in molten KCl at a temperature of 1100 K and at a density of 20 nm3. At short correlation times that are less than 5 ps, the numerical errors are significant; they however, reduce almost exponentially with increasing correlation time. It is interesting to note that the rate of decrease of the relative error is considerably weaker beyond B20 ps; in our simulations, we have used a fixed correlation time of 25 ps, which generally results in a reconstruction error of 3 to 4%. 4.2
Fig. 2 Comparison of absolute (molar) free energy of molten KCl from the 2PT and TI14 methods at a temperature of 1100 K and a density of 20 nm3.
The free energy computed by the 2PT method (circles) and coupling TI method14 (squares) shown in Fig. 2 confirms the accuracy of the 2PT method against a well-established method; the relative error between the predictions is typically less than 1%. Unlike the coupling TI method, the 2PT method requires only modest computing time to derive the free energy with comparable precision. Also shown in the figure is the reconstructed free energy (dotted line) from a thermodynamic integration approach using the internal energy computed from the 2PT method.5 First the entropy at a certain temperature (T2) is computed as ð T2 1 dEk dT (29) Sk ðT2 Þ ¼ Sk ðT1 Þ þ T dT 2PT T1 where Ek is the internal energy determined from the 2PT method as shown in eqn (7a) for the kth component of the mixture. The internal energy is next linearly fitted to temperature which gives a constant derivative for Ek (the specific heat at constant volume). Then eqn (29) is integrated (numerically or exactly) to evaluate the entropy at different states, followed by the reconstruction of the free energy as Ak = Ek TSk. As evident from Fig. 2, the reconstructed free energy is in close agreement with the original estimate; somewhat better agreement can be obtained from a higher order fit to the energy variation with temperature.
Benchmarking of the 2PT method
We will now turn to benchmarking the predictions from the 2PT method with published data on free energy of molten KCl determined from the thermodynamic integration (TI) method. Absolute free energy has been computed by Rodrigues and Fernandes using a coupling TI method with the BHM potential.14 The Einstein crystal method is used to compute the absolute free energy of the solid state, while a two-step procedure is employed for the liquid state. First the BHM potential is converted to a half-wing repulsive Gaussian potential with the repulsive barrier kept close to that of the original potential. Next the Gaussian potential is converted into a null potential to transition into the ideal gas reference state with known thermodynamic properties. For the liquid states, checks have been made to ensure that the TI paths are reversible and no phase change occurs along them.14
3066 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
4.3
Free energy and entropy change across melting in KCl
A first-order phase transition occurs at the thermodynamic melting point, which is exemplified by a slope change for the free energy along with a discontinuity in properties such as entropy. In Fig. 3, the free energies of the solid and liquid phases are plotted as a function of temperature – a discernible change in the slope occurs at B1025 K, which corresponds to the melting point of KCl. This 2PT estimate compares favorably with the experimental melting point of 1043 K with a relative error of 1.5%, approximately. It may be noted that the fluidicity parameter for the solid states is negligible. The melting point is also measurable from other properties from a MD simulation such as from the change in the structure function; the 2PT method however, is particularly advantageous as it can compute the properties such as entropy and specific
This journal is © the Owner Societies 2014
View Article Online
Paper
PCCP
close agreement again highlights the self-consistency of the 2PT methodology. From the top inset, it can be observed that the fluidity parameter ( f ) shows a jump upon melting; as expected, the fluidicity values are close to zero for the solid state.
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
4.4 Molar free energy and entropy of eutectic LiCl–KCl mixture
Fig. 3 Molar free energy across melting in KCl with the 2PT method. (bottom inset) Entropy jump from the 2PT method and comparison to NIST standard entropy data (open circles).42 (top inset) Change in the fluidicity parameter (f) across melting.
heat which undergo a discontinuity across the phase change. The bottom inset shows an entropy jump of 0.0273 kJ(mol K)1 across melting which agrees very well with the NIST data of 0.0269 kJ(mol K)1 (ref. 42) at standard conditions, incurring a relative error of 1.5%. However, there is a systematic overprediction by the 2PT method leading to a relative error of 2–3% between the absolute entropy values for the liquid states and somewhat higher for the solid states; the difference can partly arise from the inaccuracy of the interionic potential. The dotted lines in the bottom inset correspond to the reconstructed entropy values obtained by integrating eqn (29) – the
Lithium ions in the LiCl–KCl mixtures behave very differently from the potassium or chlorine ions; our recent investigation shows that the molar flux correlations for the Li+ group (Li+–Li+, Li+–Cl and Li+–K+) behave very differently from the K+ group.43 Of particular interest is the formation of long lived dynamical cage for the Li+–Li+ and Li+–Cl interactions, even at high temperatures, which has been ascribed to the intermediate-ranged ordering of the Li+ ions.35,36 An incomplete mixing of LiCl and KCl has been put forward as a plausible mechanism for the intermediateranged ordering, which manifests as a pre-peak in the structure function at 1 Å, approximately. In the current study, we investigate the density of states, particularly the partitioning into solid-like and gas-like components, and probe the dynamical characteristics that determine the thermodynamic properties. The total DoS for the K+ ions (see Fig. 4) is similar to that observed in typical liquids with a discernible zero-frequency component that augments with increasing temperature, and a single peak at a frequency of B2 THz. In contrast the total DoS for the Li+ ions portrays a plateau-like peak region spanning a broad range of frequencies, which is well-preserved even at high temperatures. Interestingly, the magnitude of the total DoS at zero frequency for the Li+ ions is almost zero indicating a solid-like behavior, which is further affirmed by the DoS partitions; the HO contribution practically accounts for the total DoS with negligible
Fig. 4 Density of states (DoS): total (left column), solid-like (HO, middle column), and gas-like (HS, right column) for Li+ (top row) and K+ ions (bottom row) at different temperatures (773 K to 1300 K). The DoS plots for different temperatures – 773 K, 850 K, 950 K, 1043 K, 1100 K, 1150 K, 1200 K, 1300 K – are off-set by 40 ps for visual clarity.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3067
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
PCCP
Fig. 5 Free energy and entropy for LiCl–KCl mixture at eutectic composition (xKCl = 0.42).
contribution from the hard spheres. The Li+ behavior is consistent with the persistent and prolonged backscattering observed in the molar flux correlations,43 as well as the presence of the intermediate-ranged ordering.36 The total DoS for the K+ ions, on the other hand, is dominated by the contribution from the harmonic oscillators at low temperatures; as temperature increases, the contribution from the hard spheres also increases significantly. Interestingly, the total and partial DoS characteristics for the Cl ions are very similar to those of the K+ ions; however, with a more prominent contribution from the hard sphere phase. Thus in the molten LiCl–KCl state, even at high temperatures, Li+ ions depict a predominantly solid-like characteristic, while the K+ and Cl ions portray a more conventional fluid-like behavior with a significant contribution from the hard sphere phase. Next we show the free energy and entropy of LiCl–KCl mixture at eutectic composition for different temperatures in Fig. 5 (for 1 mole of LiCl and 1 mole of KCl). With the addition of LiCl to KCl, the melting point reduces from 1043 K and attains the lowest value of 626 K at the eutectic composition (xKCl = 0.42). As expected the entropy increases with increasing temperature and shows a somewhat linear behavior from B1050 K, which perhaps, is related to the occurrence of melting of pure KCl around this temperature (1043 K); our recent work also shows an instability in the K+–Li+ Maxwell–Stefan (MS) diffusivity at B1100 K.43 Further investigations are needed to ascertain the physical significance, if any, for this behavior.
5. Concluding remarks The two-phase thermodynamic (2PT) method developed by Lin, Blanco and Goddard22 is emerging as a robust and accurate methodology for determining the absolute thermodynamic properties of fluids and fluid mixtures. Although the method has only been recently developed, the idea that the liquid state can be considered as an intimate mixture of solid-like and gas-like components dates back to the early work of Eyring and Ree44 who proposed that a liquid molecule briefly exhibits solid-like properties before showing a gas-like behavior as it jumps into neighboring vacancies. In the 2PT method, the partition
3068 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
Paper
function of a complex liquid state is constructed through a superposition of solid-like and gas-like vibrational density of states (DoS). We have applied this method to calculate the absolute free energy and entropy of molten KCl, and LiCl–KCl mixture at eutectic composition. The 2PT predictions show excellent agreement with those from a coupling thermodynamic integration (TI) method; for comparable states, the relative error is typically less than 1%. Further the 2PT method accurately predicts the melting point of molten KCl and the entropy change across melting with errors not exceeding few percents. The 2PT method is less reliant on prior judgment – as pointed out,13 the TI method generally benefits from the judicious choice of potential parameters/optimal paths for determining the absolute properties, which are not always known, a priori. As shown in this investigation, the 2PT method is robust and can handle complex liquids and mixtures with relative ease. These desirable features make the 2PT method very appealing for chemical and nuclear reprocessing applications where the processes entail a large number of dissolved species. Further there is a satisfying physical basis for the 2PT method, which allows the extraction of dynamical characteristics that determine the thermodynamic properties; as elucidated before, the 2PT method is unique in this aspect. We have shown that the Li+ ions in molten LiCl–KCl have very dominating solidlike and practically negligible fluid-like characteristics, unlike the K+ and Cl ions. This may have practical implications – for example, it is easier for the dissolved species in the presence of solid-like ions to change the solvent characteristics, say to highly viscous, glassy states.21 Thus there may be a trade-off in having a eutectic mixture of LiCl–KCl – on one hand, a significantly lower melting point can be achieved that is of immense value from a practical point of view, while on the other, the presence of solidlike ions, particularly close to the melting point, may not be optimal for reprocessing and waste recovery applications.
Acknowledgements Funding from DOE-NEUP program is gratefully acknowledged. The authors are thankful for the useful discussions with Dane Morgan at the University of Wisconsin.
References 1 T. A. Pascal, D. Scharf, Y. Jung and T. D. Kuhne, J. Chem. Phys., 2012, 137, 244507. 2 A. Debnath, B. Mukherjee, K. G. Ayappa, P. K. Maiti and S.-T. Lin, J. Chem. Phys., 2010, 133, 174704–174714. 3 T. Schilling and F. Schmid, J. Chem. Phys., 2009, 131, 231102–231104. 4 S.-T. Lin, P. K. Maiti and W. A. Goddard, J. Phys. Chem. B, 2010, 114, 8191–8198. 5 S.-N. Huang, T. A. Pascal, W. A. Goddard, P. K. Maiti and S.-T. Lin, J. Chem. Theory Comput., 2011, 7, 1893–1901. 6 V. Ovchinnikov, M. Cecchini and M. Karplus, J. Phys. Chem. B, 2013, 117, 750–762.
This journal is © the Owner Societies 2014
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
Paper
7 A. M. Teweldeberhan and S. A. Bonev, Phys. Rev. B, 2011, 83, 134120. 8 T. Kato, T. Inoue, T. Iwai and Y. Arai, J. Nucl. Mater., 2006, 357, 105–114. 9 J. J. Katz, G. T. Seaborg and L. R. Morss, The Chemistry of the Actinide Elements (II ed.), Chapman and Hall, New York, 1986. 10 M. Salanne, C. Simon, P. Turq and P. A. Madden, J. Phys. Chem. B, 2008, 112, 1177–1183. 11 P. Masset, R. J. M. Konings, R. Malmbeck, J. Serp and J.-P. Glatz, J. Nucl. Mater., 2005, 344, 173–179. 12 M. Nayhouse, A. M. Amlani, V. R. Heng and G. Orkoulas, J. Phys.: Condens. Matter, 2012, 24, 155101. 13 J. Anwar, D. Frenkel and M. G. Noro, J. Chem. Phys., 2003, 118, 728–735. 14 P. C. R. Rodrigues and F. M. S. S. Fernandes, J. Chem. Phys., 2007, 126, 024503–024510. 15 S. Habershon and D. E. Manolopoulos, J. Chem. Phys., 2011, 135, 224111–224116. 16 Y. Katayama, T. Mizutani, W. Utsumi, O. Shimomura, M. Yamakata and K.-i. Funakoshi, Nature, 2000, 403, 170–173. 17 F. Mallamace, M. Broccio, C. Corsaro, A. Faraone, D. Majolino, V. Venuti, L. Liu, C.-Y. Mou and S.-H. Chen, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 424. 18 M. Beye, F. Sorgenfrei, W. F. Schlotter, W. Wurth and ¨hlisch, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 16772–16776. A. Fo 19 G. N. Greaves, M. C. Wilding, S. Fearn, D. Langstaff, F. Kargl, ´rus, C. J. Benmore, R. Weber, C. M. S. Cox, Q. V. Van, O. Maje Martin and L. Hennet, Science, 2008, 322, 566–570. 20 X. Mei and J. Eapen, Phys. Rev. B, 2013, 87, 134206. 21 J. Eapen, J. Li and S. Yip, Phys. Rev. E, 2007, 76, 062501. 22 S.-T. Lin, M. Blanco and W. A. Goddard III, J. Chem. Phys., 2003, 119, 11792. 23 T. A. Pascal, S.-T. Lin and W. A. Goddard III, Phys. Chem. Chem. Phys., 2011, 13, 169–181.
This journal is © the Owner Societies 2014
PCCP
24 C. Zhang, L. Spanu and G. Galli, J. Phys. Chem. B, 2011, 115, 14190–14195. 25 P.-K. Lai, C.-M. Hsieh and S.-T. Lin, Phys. Chem. Chem. Phys., 2012, 14, 15206–15213. 26 T. A. Pascal, W. A. Goddard, P. K. Maiti and N. Vaidehi, J. Phys. Chem. B, 2012, 116, 12159–12167. 27 N. F. Carnahan and K. E. Starling, J. Chem. Phys., 1970, 53, 600–603. 28 D. A. McQuarrie, Statistical Mechanics, University Science Books, California, 2000. 29 P. H. Berens, D. H. J. Mackay, G. M. White and K. R. Wilson, J. Chem. Phys., 1983, 79, 2375–2389. 30 J. P. Boon and S. Yip, Molecular Hydrodynamics, Dover, New York, 1991. 31 J. G. Kirkwood and F. P. Buff, J. Chem. Phys., 1951, 19, 774–777. 32 R. D. Shannon, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Cryst., 1976, 32, 751. 33 M. P. Tosi and F. G. Fumi, J. Phys. Chem. Solids, 1964, 25, 45–52. 34 F. Lantelme and P. Turq, J. Chem. Phys., 1982, 77, 3177–3187. 35 M. C. C. Ribeiro, J. Phys. Chem. B, 2003, 107, 4392–4402. 36 M. Salanne, C. Simon, P. Turq and P. A. Madden, J. Phys.: Condens. Matter, 2008, 20, 332101. 37 N. Ohtori, M. Salanne and P. A. Madden, J. Chem. Phys., 2009, 130, 104507–104513. 38 M. Salanne and P. A. Madden, Mol. Phys., 2011, 109, 2299–2315. 39 N. Galamba, C. A. N. d. Castro and J. F. Ely, J. Chem. Phys., 2004, 120, 8676. 40 S. J. Plimpton, J. Comput. Phys., 1995, 117, 1–19. 41 D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge, II edn, 2004. 42 NIST, http://webbook.nist.gov/chemistry/, 2013. 43 B. Chakraborty, J. Wang and J. Eapen, Phys. Rev. E, 2013, 87, 052312. 44 H. Eyring and T. Ree, Proc. Natl. Acad. Sci. U. S. A., 1961, 47, 526–537.
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3069
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 3062
View Journal | View Issue
Absolute thermodynamic properties of molten salts using the two-phase thermodynamic (2PT) superpositioning method Jin Wang, Brahmananda Chakraborty and Jacob Eapen* We show that the absolute thermodynamic properties of molten salts (mixtures of KCl and LiCl) can be accurately determined from the two-phase thermodynamic (2PT) method that is based on superpositioning of solid-like and gas-like (hard-sphere) vibrational density of states (DoS). The 2PT predictions are in
Received 25th June 2013, Accepted 28th October 2013
excellent accordance with those from the thermodynamic integration method; the melting point of KCl evaluated from the free energy and the absolute entropy shows close conformity with the experimental/ NIST data. The DoS partitioning shows that the Li+ ions in the eutectic LiCl–KCl molten mixture are largely
DOI: 10.1039/c3cp52632a
solid-like, unlike the K+ and Cl ions, which have a significant gas-like contribution, for temperatures
www.rsc.org/pccp
employed for chemical and nuclear reprocessing applications.
ranging from 773 K to 1300 K. The solid-like states of the Li+ ions may have practical implications when
1. Introduction There is an engaging interest in determining the absolute thermodynamic properties such as free energy, enthalpy, entropy and specific heat in the liquid and disordered states.1–7 The term absolute is used in the sense that the thermodynamic properties can be evaluated directly through the partition function for a given state in the canonical ensemble; there is however, an arbitrary reference energy which is associated with the potential energy of the system. Knowing the absolute properties is particularly useful for optimizing the back-end of the nuclear fuel cycle, which entails the crucial step of separation of actinides from the spent nuclear fuel, followed by the development of new fuel through recycling. The advantage of having a closed fuel cycle includes a reduction of high-level nuclear waste and the optimal utilization of the fissionable isotopes. While aqueousbased (wet) processes are well-established there is a current interest for high temperature pyrometallurgical processes8–10 that utilize molten salts such as LiCl or KCl. A eutectic mixture of LiCl and KCl is particularly favored because of its lower melting point (626 K) relative to that of the constituent salts (883 K for LiCl and 1043 K for KCl). Originally developed for spent metallic fuel from the fast reactors, pyro processes can also be used for spent fuel from the current light water reactors (LWRs) using electroreduction, electrorefining, and pyropartitioning; they are also regarded to be more proliferation-resistant than the traditional wet processes. Department of Nuclear Engineering, North Carolina State University, Raleigh, NC 27695, USA. E-mail: [email protected]; Tel: +1 919-515-5952
3062 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
Reprocessing of nuclear fuel involves a large number of fissionable actinides8,11 as well as non-fissionable fission products; the separation processes typically depend on the chemical activity of the species, which in turn depends on the relative free energy.10 Thermodynamic integration (TI), which is the method of choice for a variety of applications,3,6,12–15 especially with molecular dynamics (MD) simulations, may not be the most optimal choice for nuclear reprocessing applications given the inordinate number of species that are targeted for separation. In several instances, TI can be limited by the optimal choice of the integration path that is not known a priori, and is also constrained in having a reverse path that implicitly excludes a first order phase change. Another challenge arises in proscribing appropriate reference systems, which can allow computation of absolute thermodynamic properties. While TI can successfully predict free energy differentials with actinide transmutation,10 integration paths typically require trajectories that span over 1 ns, which is expensive by current standards. A problem of a different nature may arise with actinide dissolution in molten salts – depending on the volume fraction, the interaction strength and the thermodynamic state, the dissolved species may form glassy states or even undergo a first-order liquid–liquid phase transition,16–20 which is somewhat poorly understood currently. Very strong attractive interactions have been shown to cause an amorphous transition in a model nanocolloidal system with a volume fraction as low as 5%.21 Thus it would be profitable to identify a methodology that can obviate some or all of the aforesaid possible drawbacks for systems with large number of species in complex liquid forms. Originally developed by Lin, Blanco and Goddard,22 the twophase thermodynamic (2PT) model has recently gained wide
This journal is © the Owner Societies 2014
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
Paper
PCCP
attention and success in predicting absolute thermodynamic properties for several types of fluids.1,2,4,5,7,23–26 Motivated by the approach of evaluating the partition function in the solid state using idealized harmonic oscillators, the central idea behind the 2PT methodology is to construct the partition function of a complex liquid state through a superposition of solid-like and gas-like vibrational density of states (DoS). An exact representation of the idealized states can be derived by treating the gas-like state through hard-sphere interactions with weighting functions provided by the Carnahan–Starling equation of state,27 and assuming the solid-like state to be a collection of independent harmonic oscillators. The thermodynamic properties of a liquid system thus can be shown to be a combination of the properties of a set of harmonic oscillators (solid-like) and a set of hard spheres (gas-like). Since the 2PT method thus far has not been employed to evaluate the properties of molten salts and mixtures, its efficacy or accuracy is currently unknown. In this investigation, we show that the absolute thermodynamic properties of LiCl/KCl molten salts can be accurately determined from the 2PT method. The 2PT free energy is in excellent agreement with that determined from the thermodynamic integration (TI) method, with relative errors of 1% or less, typically. The melting point of ionic KCl (1025 K) compares favorably with the experimental data (1041 K); the entropy of the melt phase also shows good conformity with the recommended data from NIST. We then determine the absolute free energy and entropy which are currently unknown for a eutectic mixture of LiCl and KCl for temperatures ranging from 773 K to 1300 K. From the partitioning of the DoS, which is uniquely possible in the 2PT method, we show that the Li+ ions in the LiCl–KCl mixture are dominantly solid-like even at temperatures as high as 1300 K. For chemical and nuclear applications the solid-like states can potentially impose practical limits, for example, to the amount of dissolved species that can be practically accommodated during reprocessing.
2. The 2PT methodology The thermodynamic properties in the canonical ensemble can be evaluated from the canonical partition function (Q); for example, the internal energy (E), entropy (S), specific heat at constant volume (Cv), and the Helmoholtz free energy (A) can be evaluated as28 @ ðln QÞ (1a) E ¼ kB T 2 @T @ S ¼ kB T ðln QÞ þ ln Q @T
Cv ¼
@ @ kB T 2 ðln QÞ @T @T A = kBT(ln Q)
This journal is © the Owner Societies 2014
(1b)
(1c)
(1d)
where kB is the Boltzmann constant and T is the absolute temperature. The thermodynamic properties are seldom evaluated through the partition function as it is known exactly only for idealized systems such as harmonic oscillators and ideal gas. Using a normal mode analysis, a solid state can be approximated as a system of non-interacting harmonic oscillators (HO); for this case, the total canonical partition function (QHO) can be expressed as29 QHO ¼
3N Y
qi
(2)
i¼1
where qi is the partition function of the ith harmonic oscillator mode, and N is the total number of oscillators. The total partition function can be shown to be related to the density of states as29 ð1 ln QHO ¼ dnGðnÞ ln½qðnÞ (3) 0
where G(n) is the density of states, which is given by22 ð 4 1 GðnÞ ¼ dtCðtÞei2pnt kB T 0
(4)
where C(t) is the mass-weighted velocity autocorrelation function,28 which can be directly computed from a molecular dynamics (MD) simulation (ab initio or classical). Since an exact quantum-mechanical expression is known for q(n), the total partition function can be evaluated using eqn (3) and (4), and the thermodynamic properties such as given by eqn (1a) to (1d) can be assessed for the solid states with reasonable accuracy.29 For a liquid state, the system cannot be considered as a collection of harmonic oscillators and a direct application of the above approximation will lead to non-physical results. The 2PT method22 extends the above idea to the liquid state by assuming a superposition of two idealized states – hard spheres (HS) that correspond to a gas-like state, and harmonic oscillators that correspond to a solid-like state. First it is hypothesized that the total density of states of the system that is now regarded to be a non-interacting mixture of hard spheres and harmonic oscillators can be decomposed as HO Gk(n) = GHS k (n) + Gk (n)
(5)
where k denotes the different species in the mixture. A key step in the 2PT method is in defining and evaluating a fluidicity parameter ( f ) given by the following expression for partitioning the system into the idealized states. Ð1 dnGHS 3NkHS NkHS k ðnÞ fk ¼ Ð01 ¼ (6) ¼ 3Nk Nk 0 dnGk ðnÞ where NHS k is the effective number of hard sphere atoms or ions of the kth component of the mixture. Note that the integral of the total density of states over all frequencies is simply 3Nk. Thus the 3Nk degrees of freedom of the kth component is assumed to be composed of 3NHS or equivalently, 3fkNk hard k sphere or diffusive degrees of freedom, and 3(Nk NHS k ) or equivalently, 3Nk(1 fk) solid-like or non-diffusive degrees of freedom.25 With this approximation, the thermodynamic
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3063
View Article Online
PCCP
Paper
h i 1 sHS ðnÞ HS HS WAHS k ðnÞ ¼ W ðnÞ W ðnÞ ¼ Ek Sk k 2 3kB
properties are formally expressed as25 Ek ¼ Ek0 þ kB T
ð 1 0
HS dnGHS k ðnÞWEk ðnÞ þ
ð1
HO dnGHO k ðnÞWE ðnÞ
0
(7a) Sk ¼ kB
ð 1
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
0
Cv;k ¼ kB
HS dnGHS k ðn ÞWSk ðn Þ þ
ð 1 0
ð1 0
HO dnGHO k ðn ÞWS ðn Þ
HS dnGHS k ðnÞWCv;k ðnÞ þ
ð1 0
(7b)
HO dnGHO k ðnÞWCv ðnÞ
(7c) Ak ¼
Ek0
þ kB T
ð 1 0
HS dnGHS k ðnÞWAk ðnÞ
ð1
þ
0
HO dnGHO k ðnÞWA ðnÞ
(7d) where W is a weight function for the appropriate property, and E0k is the reference energy which is given by ref. 25 fk 0 MD Ek ¼ Ek 3Nk kB T 1 (8) 2 is the total energy from MD simulations for the kth where EMD k component. Thus each extensive thermodynamic property, for each mixture component, comprises of two contributions – one from diffusive (hard spheres) degrees of freedom, and the other from vibrating (harmonic oscillators) degrees of freedom. The distinguishing part of the 2PT method is in deriving an exact expression for GHS k given by ref. 22. " #1 pGk ð0Þn 2 HS Gk ðnÞ ¼ Gk ð0Þ 1 þ (9) 6fk Nk Thus the density of states for the hard sphere phase can be calculated as a function of the fluidicity parameter and Gk(0) which is just the zero-frequency value of the system (total) density of states. MD simulations are employed to calculate Gk(n) from the velocity autocorrelation function. With a known GHS k , the density of states of the HO phase is then evaluated as HS GHO k (n) = [Gk(n) Gk (n)]. The weight functions are known exactly for the HO and HS idealized states, and are given by5,25 WEHO ¼ WSHO
bhn bhn þ bhn 2 e 1
bhn ln 1 ebhn ¼ bhn e 1 WCHO ¼ v
ðbhnÞ2 ebhn ðebhn
1Þ
1 ebhn WAHO ðnÞ ¼ ln ebhn=2
WSHS ðnÞ ¼ k
1 2
sHS k ðnÞ 3kB
3064 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
where b1 = kBT, h is the Planck constant and sHS is the entropy of the hard sphere phase (per atom/ion), which in the 2PT model is evaluated as ( ) yk 4Þ y^k ð3^ HS IG sk ¼ sk þ kB ln½zk ðy^k Þ þ (12) ð1 y^k Þ2 The first term of the above equation denotes the ideal gas entropy, however, weighted only for the HS phase. It is given by ref. 5: ( " #) 5 2pmk kB T 3=2 Vk þ ln (13) ¼ k sIG B k 2 h2 f k Nk where mk and Nk are the mass and the total number of ions/ atoms of the kth species, respectively. The term zk PVk/(NkkBT) in eqn (12) delineates the compressibility of the hard sphere phase, which is determined from the Carnahan–Starling equation given by zk ðy^k Þ ¼
1 þ y^k þ y^k2 y^k3
(10b)
(10d) (11a) (11b)
(14)
ð1 y^k Þ3
where ˆyk is the weighted hard sphere packing fraction (by the fluidicity parameter) for the kth component, and is given by ˆyk = fkyk, where yk denotes the packing fraction which is expressed as 3 prk sHS k yk ¼ (15) 6 is the hard sphere where rk is the number density and sHS k diameter of the kth species that in turn influences the fluidicity parameter (fk). In the 2PT methodology it is assumed that the fluidicity parameter is simply a ratio of the self-diffusivity (Dk) to the hard sphere diffusivity at zero pressure (D0HS ), which can k be written as f ¼
DðT; PÞ HS D0HS k ðT; P; s Þ
(16)
The self-diffusivity can be evaluated from MD simulations as the time integral of the velocity autocorrelation, and is expressible as Dk ðT; PÞ ¼
(10c)
2
WEHS ðnÞ ¼ WCHS ðnÞ ¼ k v;k
(10a)
(11c)
kB TGk ð0Þ 12mk Nk
(17)
where Gk(0) is evaluated from the expression ð 4 1 dtCk ðtÞei2pnt Gk ðnÞ ¼ kB T 0
(18)
We invoke the stationary principle of classical correlation functions,30 and thus only the real part of the Fourier transform needs to be used for calculating the above integral. Note that C(t) is the mass-weighted velocity autocorrelation function defined as Nk ð 1 Nk ð 1 X X Ck ðtÞ ¼ mk dthvk ðtÞ:vk ð0Þi ¼ mk dthck ðtÞi (19) i¼1
0
i¼1
0
This journal is © the Owner Societies 2014
View Article Online
Paper
PCCP
where ck(t) is the velocity autocorrelation function of the kth component. The hard sphere diffusivity in the limit of zero pressure/density is computed exactly from the Chapman– Enskog’s theory, which is given by 3 kB T 1=2 HS D0HS ¼ T; P; s (20) k k 2 pm k 8r sHS
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
k
k
The velocity autocorrelation for HS phase decays exponentially in time and hence the HS diffusivity can be expressed as22 DHS k ðT; fk rk Þ ¼
kB TGk ð0Þ 12mk fk Nk
(21)
The above expression can be related to the zero-pressure hard sphere diffusivity through the following relationship 0HS T; fk rk ; sHS DHS k ðT; fk rk Þ ¼ Dk k
4^ yk ½zðy^k Þ 1
(22)
where ˆyk = fky. After a few algebraic steps the following relationship can be established25 2Dk9/2fk15/2 6Dk3fk5 Dk3/2fk7/2 + 6Dk3/2fk5/2 + 2fk 2 = 0 (23) where 2Gk ð0Þ pkB T 1=2 Nk 1=3 6 2=3 Dk ¼ 9Nk mk p Vk
where sk and xk are the ionic radius and mole fraction, respectively, of the kth component. The Kirkwood and Buff (KB) theory31 gives a better approximation for partial volumes but it entails the evaluation of somewhat ill-convergent KB integrals. We have used a fixed ionic radius of 0.152 nm, 0.167 nm and 0.09 nm for K+, Cl and Li+ ions, respectively.32 As shown later, the ionic radius approximation works quite well for the current application. We have also observed that the calculated thermodynamic properties are not sensitive to moderate changes in the ionic radii. Finally, the molar (m) thermodynamic properties of the mixture are determined as xk Ekm
k¼1
m A ¼
n X
xk Skm R^
n X
xk lnðxk Þ
(26b)
k¼1
^ x k Am k þ RT
k¼1
n X
xk lnðxk Þ
(26c)
k¼1
ˆ is the universal gas constant, and n is the total number where R of components in the mixture.
3. Molecular dynamics simulations MD simulations are principally employed for calculating the velocity autocorrelation function. The molecular system typically consists of 1000 ions interacting through a long-ranged electrostatic potential and a Born–Huggins–Mayer (BHM) short-ranged potential. Originally developed by Fumi and Tosi,33 the rigid-ion potential is able to predict several thermodynamic, structural and transport properties of molten LiCl–KCl mixtures34,35 and other alkali halides33 with reasonable accuracy. The effect of polarization, however, is important for multicomponent systems,36–39 although it tends to be small for monovalent systems. The functional form of short-ranged part of the BHM potential is given by
(24)
k¼1
n X
n X
Fshort ¼ Bij earij ij
Eqn (24) contains the details of the physical state as well as the zero-frequency density of states that can be evaluated from MD simulations. Eqn (23) now can be iteratively solved to get the fluidicity parameter (fk) for each component of the mixture. Once f is known, the density of states and the weight functions of the HS phase can be computed, and the HO density of states can be determined as GHO = (Gk GHS k k ). Thus all the properties given by eqn (7a) to (7d) can be determined, for each component of the mixture. In this investigation, the partial volume (Vk) in eqn (24) is calculated from the relative ionic radii, which is given by sk3 V (25) Vk ¼ n P N 3 xk sk
m
m S ¼
Cij Dij rij6 rij8
(27)
where rij is the interionic separation between two ions, i and j, and a, B, C and D are constants.34,35 While the first term corresponds to the electron cloud repulsion, the second and third terms account for the dipole–dipole and the dipole–quadrupole dispersion interactions, respectively. The MD simulations are performed using the LAMMPS software40 with periodic boundary conditions at zero pressure using a timestep of 1 fs. Particle–particle–particle– mesh is used for columbic interactions, and the temperature and ´–Hoover algorithm. From pressure are controlled using the Nose the velocity data collected from MD simulations, the velocity autocorrelation function (VACF) is constructed, typically with a correlation time of 25 ps, using an overlapped data structure.41 From the Fourier transform of the VACF, the 2PT thermodynamic properties are computed using an in-house computer program.
4. Results and discussion We will first discuss (1) the numerical convergence of the fluidicity parameter with correlation time, followed by (2) a benchmark investigation of the 2PT method against published results from thermodynamic integration (TI) method for KCl, (3) computation of the molar free energy and the entropy change across melting in KCl along with a comparison to NIST data,42 and (4) an assessment of the density of states and thermodynamic properties of a eutectic mixture of LiCl–KCl for different temperatures. 4.1
Numerical convergence of the fluidicity parameter
k¼1
The fluidicity parameter for each component of the mixture is iteratively solved using eqn (23). Further the partitioning of the
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3065
E ¼
(26a)
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
PCCP
Paper
Fig. 1 Relative error in reconstructing the fluidicity parameter for the K+ ions in molten KCl.
system into solid-like and gas-like states depends on resolving the zero-frequency value of the density of states. To assess the optimal correlation time (for VACF), we have first computed fk using eqn (23) and then compared with the reconstructed ˆfk from the equation Ð1 ð HS 1 1 0 dnGk ðnÞ ^ Ð dnGHS ¼ fk ¼ 1 k ðnÞ 3Nk 0 0 dnGk ðnÞ
(28)
Insufficient sampling and numerical errors in the evaluation of VACF and the density of states will be reflected as a difference between ˆfk and fk. Fig. 1 shows the error in the reconstruction, expressed as a percentage, for the K+ ions in molten KCl at a temperature of 1100 K and at a density of 20 nm3. At short correlation times that are less than 5 ps, the numerical errors are significant; they however, reduce almost exponentially with increasing correlation time. It is interesting to note that the rate of decrease of the relative error is considerably weaker beyond B20 ps; in our simulations, we have used a fixed correlation time of 25 ps, which generally results in a reconstruction error of 3 to 4%. 4.2
Fig. 2 Comparison of absolute (molar) free energy of molten KCl from the 2PT and TI14 methods at a temperature of 1100 K and a density of 20 nm3.
The free energy computed by the 2PT method (circles) and coupling TI method14 (squares) shown in Fig. 2 confirms the accuracy of the 2PT method against a well-established method; the relative error between the predictions is typically less than 1%. Unlike the coupling TI method, the 2PT method requires only modest computing time to derive the free energy with comparable precision. Also shown in the figure is the reconstructed free energy (dotted line) from a thermodynamic integration approach using the internal energy computed from the 2PT method.5 First the entropy at a certain temperature (T2) is computed as ð T2 1 dEk dT (29) Sk ðT2 Þ ¼ Sk ðT1 Þ þ T dT 2PT T1 where Ek is the internal energy determined from the 2PT method as shown in eqn (7a) for the kth component of the mixture. The internal energy is next linearly fitted to temperature which gives a constant derivative for Ek (the specific heat at constant volume). Then eqn (29) is integrated (numerically or exactly) to evaluate the entropy at different states, followed by the reconstruction of the free energy as Ak = Ek TSk. As evident from Fig. 2, the reconstructed free energy is in close agreement with the original estimate; somewhat better agreement can be obtained from a higher order fit to the energy variation with temperature.
Benchmarking of the 2PT method
We will now turn to benchmarking the predictions from the 2PT method with published data on free energy of molten KCl determined from the thermodynamic integration (TI) method. Absolute free energy has been computed by Rodrigues and Fernandes using a coupling TI method with the BHM potential.14 The Einstein crystal method is used to compute the absolute free energy of the solid state, while a two-step procedure is employed for the liquid state. First the BHM potential is converted to a half-wing repulsive Gaussian potential with the repulsive barrier kept close to that of the original potential. Next the Gaussian potential is converted into a null potential to transition into the ideal gas reference state with known thermodynamic properties. For the liquid states, checks have been made to ensure that the TI paths are reversible and no phase change occurs along them.14
3066 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
4.3
Free energy and entropy change across melting in KCl
A first-order phase transition occurs at the thermodynamic melting point, which is exemplified by a slope change for the free energy along with a discontinuity in properties such as entropy. In Fig. 3, the free energies of the solid and liquid phases are plotted as a function of temperature – a discernible change in the slope occurs at B1025 K, which corresponds to the melting point of KCl. This 2PT estimate compares favorably with the experimental melting point of 1043 K with a relative error of 1.5%, approximately. It may be noted that the fluidicity parameter for the solid states is negligible. The melting point is also measurable from other properties from a MD simulation such as from the change in the structure function; the 2PT method however, is particularly advantageous as it can compute the properties such as entropy and specific
This journal is © the Owner Societies 2014
View Article Online
Paper
PCCP
close agreement again highlights the self-consistency of the 2PT methodology. From the top inset, it can be observed that the fluidity parameter ( f ) shows a jump upon melting; as expected, the fluidicity values are close to zero for the solid state.
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
4.4 Molar free energy and entropy of eutectic LiCl–KCl mixture
Fig. 3 Molar free energy across melting in KCl with the 2PT method. (bottom inset) Entropy jump from the 2PT method and comparison to NIST standard entropy data (open circles).42 (top inset) Change in the fluidicity parameter (f) across melting.
heat which undergo a discontinuity across the phase change. The bottom inset shows an entropy jump of 0.0273 kJ(mol K)1 across melting which agrees very well with the NIST data of 0.0269 kJ(mol K)1 (ref. 42) at standard conditions, incurring a relative error of 1.5%. However, there is a systematic overprediction by the 2PT method leading to a relative error of 2–3% between the absolute entropy values for the liquid states and somewhat higher for the solid states; the difference can partly arise from the inaccuracy of the interionic potential. The dotted lines in the bottom inset correspond to the reconstructed entropy values obtained by integrating eqn (29) – the
Lithium ions in the LiCl–KCl mixtures behave very differently from the potassium or chlorine ions; our recent investigation shows that the molar flux correlations for the Li+ group (Li+–Li+, Li+–Cl and Li+–K+) behave very differently from the K+ group.43 Of particular interest is the formation of long lived dynamical cage for the Li+–Li+ and Li+–Cl interactions, even at high temperatures, which has been ascribed to the intermediate-ranged ordering of the Li+ ions.35,36 An incomplete mixing of LiCl and KCl has been put forward as a plausible mechanism for the intermediateranged ordering, which manifests as a pre-peak in the structure function at 1 Å, approximately. In the current study, we investigate the density of states, particularly the partitioning into solid-like and gas-like components, and probe the dynamical characteristics that determine the thermodynamic properties. The total DoS for the K+ ions (see Fig. 4) is similar to that observed in typical liquids with a discernible zero-frequency component that augments with increasing temperature, and a single peak at a frequency of B2 THz. In contrast the total DoS for the Li+ ions portrays a plateau-like peak region spanning a broad range of frequencies, which is well-preserved even at high temperatures. Interestingly, the magnitude of the total DoS at zero frequency for the Li+ ions is almost zero indicating a solid-like behavior, which is further affirmed by the DoS partitions; the HO contribution practically accounts for the total DoS with negligible
Fig. 4 Density of states (DoS): total (left column), solid-like (HO, middle column), and gas-like (HS, right column) for Li+ (top row) and K+ ions (bottom row) at different temperatures (773 K to 1300 K). The DoS plots for different temperatures – 773 K, 850 K, 950 K, 1043 K, 1100 K, 1150 K, 1200 K, 1300 K – are off-set by 40 ps for visual clarity.
This journal is © the Owner Societies 2014
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3067
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
PCCP
Fig. 5 Free energy and entropy for LiCl–KCl mixture at eutectic composition (xKCl = 0.42).
contribution from the hard spheres. The Li+ behavior is consistent with the persistent and prolonged backscattering observed in the molar flux correlations,43 as well as the presence of the intermediate-ranged ordering.36 The total DoS for the K+ ions, on the other hand, is dominated by the contribution from the harmonic oscillators at low temperatures; as temperature increases, the contribution from the hard spheres also increases significantly. Interestingly, the total and partial DoS characteristics for the Cl ions are very similar to those of the K+ ions; however, with a more prominent contribution from the hard sphere phase. Thus in the molten LiCl–KCl state, even at high temperatures, Li+ ions depict a predominantly solid-like characteristic, while the K+ and Cl ions portray a more conventional fluid-like behavior with a significant contribution from the hard sphere phase. Next we show the free energy and entropy of LiCl–KCl mixture at eutectic composition for different temperatures in Fig. 5 (for 1 mole of LiCl and 1 mole of KCl). With the addition of LiCl to KCl, the melting point reduces from 1043 K and attains the lowest value of 626 K at the eutectic composition (xKCl = 0.42). As expected the entropy increases with increasing temperature and shows a somewhat linear behavior from B1050 K, which perhaps, is related to the occurrence of melting of pure KCl around this temperature (1043 K); our recent work also shows an instability in the K+–Li+ Maxwell–Stefan (MS) diffusivity at B1100 K.43 Further investigations are needed to ascertain the physical significance, if any, for this behavior.
5. Concluding remarks The two-phase thermodynamic (2PT) method developed by Lin, Blanco and Goddard22 is emerging as a robust and accurate methodology for determining the absolute thermodynamic properties of fluids and fluid mixtures. Although the method has only been recently developed, the idea that the liquid state can be considered as an intimate mixture of solid-like and gas-like components dates back to the early work of Eyring and Ree44 who proposed that a liquid molecule briefly exhibits solid-like properties before showing a gas-like behavior as it jumps into neighboring vacancies. In the 2PT method, the partition
3068 | Phys. Chem. Chem. Phys., 2014, 16, 3062--3069
Paper
function of a complex liquid state is constructed through a superposition of solid-like and gas-like vibrational density of states (DoS). We have applied this method to calculate the absolute free energy and entropy of molten KCl, and LiCl–KCl mixture at eutectic composition. The 2PT predictions show excellent agreement with those from a coupling thermodynamic integration (TI) method; for comparable states, the relative error is typically less than 1%. Further the 2PT method accurately predicts the melting point of molten KCl and the entropy change across melting with errors not exceeding few percents. The 2PT method is less reliant on prior judgment – as pointed out,13 the TI method generally benefits from the judicious choice of potential parameters/optimal paths for determining the absolute properties, which are not always known, a priori. As shown in this investigation, the 2PT method is robust and can handle complex liquids and mixtures with relative ease. These desirable features make the 2PT method very appealing for chemical and nuclear reprocessing applications where the processes entail a large number of dissolved species. Further there is a satisfying physical basis for the 2PT method, which allows the extraction of dynamical characteristics that determine the thermodynamic properties; as elucidated before, the 2PT method is unique in this aspect. We have shown that the Li+ ions in molten LiCl–KCl have very dominating solidlike and practically negligible fluid-like characteristics, unlike the K+ and Cl ions. This may have practical implications – for example, it is easier for the dissolved species in the presence of solid-like ions to change the solvent characteristics, say to highly viscous, glassy states.21 Thus there may be a trade-off in having a eutectic mixture of LiCl–KCl – on one hand, a significantly lower melting point can be achieved that is of immense value from a practical point of view, while on the other, the presence of solidlike ions, particularly close to the melting point, may not be optimal for reprocessing and waste recovery applications.
Acknowledgements Funding from DOE-NEUP program is gratefully acknowledged. The authors are thankful for the useful discussions with Dane Morgan at the University of Wisconsin.
References 1 T. A. Pascal, D. Scharf, Y. Jung and T. D. Kuhne, J. Chem. Phys., 2012, 137, 244507. 2 A. Debnath, B. Mukherjee, K. G. Ayappa, P. K. Maiti and S.-T. Lin, J. Chem. Phys., 2010, 133, 174704–174714. 3 T. Schilling and F. Schmid, J. Chem. Phys., 2009, 131, 231102–231104. 4 S.-T. Lin, P. K. Maiti and W. A. Goddard, J. Phys. Chem. B, 2010, 114, 8191–8198. 5 S.-N. Huang, T. A. Pascal, W. A. Goddard, P. K. Maiti and S.-T. Lin, J. Chem. Theory Comput., 2011, 7, 1893–1901. 6 V. Ovchinnikov, M. Cecchini and M. Karplus, J. Phys. Chem. B, 2013, 117, 750–762.
This journal is © the Owner Societies 2014
View Article Online
Published on 07 January 2014. Downloaded by University of Waterloo on 31/10/2014 14:42:20.
Paper
7 A. M. Teweldeberhan and S. A. Bonev, Phys. Rev. B, 2011, 83, 134120. 8 T. Kato, T. Inoue, T. Iwai and Y. Arai, J. Nucl. Mater., 2006, 357, 105–114. 9 J. J. Katz, G. T. Seaborg and L. R. Morss, The Chemistry of the Actinide Elements (II ed.), Chapman and Hall, New York, 1986. 10 M. Salanne, C. Simon, P. Turq and P. A. Madden, J. Phys. Chem. B, 2008, 112, 1177–1183. 11 P. Masset, R. J. M. Konings, R. Malmbeck, J. Serp and J.-P. Glatz, J. Nucl. Mater., 2005, 344, 173–179. 12 M. Nayhouse, A. M. Amlani, V. R. Heng and G. Orkoulas, J. Phys.: Condens. Matter, 2012, 24, 155101. 13 J. Anwar, D. Frenkel and M. G. Noro, J. Chem. Phys., 2003, 118, 728–735. 14 P. C. R. Rodrigues and F. M. S. S. Fernandes, J. Chem. Phys., 2007, 126, 024503–024510. 15 S. Habershon and D. E. Manolopoulos, J. Chem. Phys., 2011, 135, 224111–224116. 16 Y. Katayama, T. Mizutani, W. Utsumi, O. Shimomura, M. Yamakata and K.-i. Funakoshi, Nature, 2000, 403, 170–173. 17 F. Mallamace, M. Broccio, C. Corsaro, A. Faraone, D. Majolino, V. Venuti, L. Liu, C.-Y. Mou and S.-H. Chen, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 424. 18 M. Beye, F. Sorgenfrei, W. F. Schlotter, W. Wurth and ¨hlisch, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 16772–16776. A. Fo 19 G. N. Greaves, M. C. Wilding, S. Fearn, D. Langstaff, F. Kargl, ´rus, C. J. Benmore, R. Weber, C. M. S. Cox, Q. V. Van, O. Maje Martin and L. Hennet, Science, 2008, 322, 566–570. 20 X. Mei and J. Eapen, Phys. Rev. B, 2013, 87, 134206. 21 J. Eapen, J. Li and S. Yip, Phys. Rev. E, 2007, 76, 062501. 22 S.-T. Lin, M. Blanco and W. A. Goddard III, J. Chem. Phys., 2003, 119, 11792. 23 T. A. Pascal, S.-T. Lin and W. A. Goddard III, Phys. Chem. Chem. Phys., 2011, 13, 169–181.
This journal is © the Owner Societies 2014
PCCP
24 C. Zhang, L. Spanu and G. Galli, J. Phys. Chem. B, 2011, 115, 14190–14195. 25 P.-K. Lai, C.-M. Hsieh and S.-T. Lin, Phys. Chem. Chem. Phys., 2012, 14, 15206–15213. 26 T. A. Pascal, W. A. Goddard, P. K. Maiti and N. Vaidehi, J. Phys. Chem. B, 2012, 116, 12159–12167. 27 N. F. Carnahan and K. E. Starling, J. Chem. Phys., 1970, 53, 600–603. 28 D. A. McQuarrie, Statistical Mechanics, University Science Books, California, 2000. 29 P. H. Berens, D. H. J. Mackay, G. M. White and K. R. Wilson, J. Chem. Phys., 1983, 79, 2375–2389. 30 J. P. Boon and S. Yip, Molecular Hydrodynamics, Dover, New York, 1991. 31 J. G. Kirkwood and F. P. Buff, J. Chem. Phys., 1951, 19, 774–777. 32 R. D. Shannon, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Cryst., 1976, 32, 751. 33 M. P. Tosi and F. G. Fumi, J. Phys. Chem. Solids, 1964, 25, 45–52. 34 F. Lantelme and P. Turq, J. Chem. Phys., 1982, 77, 3177–3187. 35 M. C. C. Ribeiro, J. Phys. Chem. B, 2003, 107, 4392–4402. 36 M. Salanne, C. Simon, P. Turq and P. A. Madden, J. Phys.: Condens. Matter, 2008, 20, 332101. 37 N. Ohtori, M. Salanne and P. A. Madden, J. Chem. Phys., 2009, 130, 104507–104513. 38 M. Salanne and P. A. Madden, Mol. Phys., 2011, 109, 2299–2315. 39 N. Galamba, C. A. N. d. Castro and J. F. Ely, J. Chem. Phys., 2004, 120, 8676. 40 S. J. Plimpton, J. Comput. Phys., 1995, 117, 1–19. 41 D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge, II edn, 2004. 42 NIST, http://webbook.nist.gov/chemistry/, 2013. 43 B. Chakraborty, J. Wang and J. Eapen, Phys. Rev. E, 2013, 87, 052312. 44 H. Eyring and T. Ree, Proc. Natl. Acad. Sci. U. S. A., 1961, 47, 526–537.
Phys. Chem. Chem. Phys., 2014, 16, 3062--3069 | 3069
