Density functional theory for carbon dioxide crystal Yiwen Chang, Jianguo Mi, and Chongli Zhong Citation: The Journal of Chemical Physics 140, 204706 (2014); doi: 10.1063/1.4878413 View online: http://dx.doi.org/10.1063/1.4878413 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/20?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Pressure dependent stability and structure of carbon dioxide—A density functional study including long-range corrections J. Chem. Phys. 139, 174501 (2013); 10.1063/1.4826929 High pressure ionic and molecular crystals of ammonia monohydrate within density functional theory J. Chem. Phys. 137, 064506 (2012); 10.1063/1.4737887 Density-functional theory for polymer-carbon dioxide mixtures: A perturbed-chain SAFT approach J. Chem. Phys. 137, 054902 (2012); 10.1063/1.4742346 Density functional theory, molecular dynamics, and differential scanning calorimetry study of the RbF–CsF phase diagram J. Chem. Phys. 130, 134716 (2009); 10.1063/1.3097550 Density functional approach for modeling CO 2 pressurized polymer thin films in equilibrium J. Chem. Phys. 130, 084902 (2009); 10.1063/1.3077861
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
THE JOURNAL OF CHEMICAL PHYSICS 140, 204706 (2014)
Density functional theory for carbon dioxide crystal Yiwen Chang, Jianguo Mi,a) and Chongli Zhong State Key Laboratory of Organic-Inorganic Composites, Beijing University of Chemical Technology, Beijing 100029, China
(Received 27 February 2014; accepted 5 May 2014; published online 27 May 2014) We present a density functional approach to describe the solid−liquid phase transition, interfacial and crystal structure, and properties of polyatomic CO2 . Unlike previous phase field crystal model or density functional theory, which are derived from the second order direct correlation function, the present density functional approach is based on the fundamental measure theory for hard-sphere repulsion in solid. More importantly, the contributions of enthalpic interactions due to the dispersive attractions and of entropic interactions arising from the molecular architecture are integrated in the density functional model. Using the theoretical model, the predicted liquid and solid densities of CO2 at equilibrium triple point are in good agreement with the experimental values. Based on the structure of crystal-liquid interfaces in different planes, the corresponding interfacial tensions are predicted. Their respective accuracies need to be tested. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4878413] I. INTRODUCTION
Crystallization phenomena are ubiquitous in both natural and artificially synthesized systems. During the nucleation process, the microstructure of new crystal is strongly influenced by the property of bulk fluid and the complex spatial patterns of solid. In order to understand crystallization from microscopic scale, some essential issues, such as the solid−liquid phase transition, the structure and energy of crystal and interface, need to be clarified. However, it has so far been difficult to obtain these quantitative descriptions theoretically. The difficulty becomes more serious when dealing with polyatomic fluid−solid systems. For simple monatomic fluids, only enthalpic interactions affect these properties. Whereas for polyatomic fluids, although enthalpic interactions are usually dominant, entropic interactions arising from the molecular architecture can play a critical role in dictating the structure and properties. A successful theoretical model should consider their cooperative contributions. There are two different routes to develop crystallization theory. One is derived from liquid-state theory, and another is developed from field theory. Due to the success of liquid-state theory, there has been much effort in describing the solid state with liquid properties, leading to the development of density functional theory (DFT). In recent years, investigations in condensed matter have given unambiguous evidence that DFT is a very versatile and powerful tool to study the structural and thermodynamic properties of crystallization.1, 2 Meanwhile, the phase-field crystal (PFC) model developed from phasefield theory has proven to be prospective in a variety of problems, such as interfaces,3 polycrystalline pattern formation,4, 5 or crystal nucleation and growth.6–8 Of later years, PFC combined with classical DFT has been developed and applied to solve this problem in a highly efficient way.9, 10 It has shown a) Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2014/140(20)/204706/7/$30.00
that the combined theory is an ideal tool to study nucleation process. In the above theoretical models, the energy expressions are unexceptionally derived from the direct correlation function with a second-order expansion form. It is well-known that such approximation gives poor results even for describing monatomic fluids in equilibrium state. The deficiency is mainly due to the inaccurate direct correlation function inside the hard-sphere. Without accurate free energy expression, both the PFC model and the DFT are a kind of coarse-grained method, which focus mainly on phenomenological description. On the other hand, the development of the free energy functional for hard-sphere is presented through its milestones in the weighted density approximation and the fundamental measure theory.11 The extensions of these approaches to hard-sphere mixtures include the treatment of polydisperse systems. Very recently, we have presented a three-dimensional (3D) DFT approach to describe the crystal−liquid interfaces of Lennard-Jones fluid.12 In the theoretical model, the fundamental measure theory is applied to describe the free energy functional of hard-sphere repulsion, and a weighted density method13 is used to describe the free energy contribution arising from the van der Waals attractions. The liquid-solid equilibrium, density profiles within crystal cell and at liquid−solid interface, and interfacial tensions were calculated for facecentered-cubic (fcc) and body-centered-cubic (bcc) nucleus. The results are in good agreement with simulation data, showing that the model is quantitatively accurate in describing nucleation thermodynamics of the Lennard-Jones fluid. In the recent years, DFT was further developed and applied to research polyatomic fluids or polymer systems. For instances, a Wertheim’s first-order thermodynamic perturbation theory (TPT1) based DFT was proposed for chain conformation and bond orientation correlation function of polyatomic fluid;14 Frischknecht et al.15 employed a DFT based on interfacial statistical associating fluid theory (SAFT)
140, 204706-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-2
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
equation to probe the first-order phase transition in athermal blend of nanoparticle and homopolymer near a wall. The developments of DFT are conducive to obtain the structure and properties for polyatomic fluids more accurately. The main purpose of this work is to extend the above 3DDFT approach to polyatomic fluids. CO2 is taken as an example. Compared with monatomic molecules, a CO2 molecule contains three sites, and its structure has a marginal deviation from the linear geometry. The contribution of enthalpy to the free energy functional is represented with the van der Waals interactions. Meanwhile, the entropic contribution is correlated to the architecture of CO2 molecule including bond length, bond angle, and connectivity. Similar to our previous work, the van der Waals interactions are split into hard-sphere repulsions and dispersive attractions, and the free energy functional for hard-sphere repulsions is still represented by the fundamental measure theory for solid phase. In particular, by solving the reference interaction site model integral equation, the nonbonded dispersive attractions, bond length, and bond angle are integrated into the direct correlation functions of different sites. Once the direct correlation functions are derived, the corresponding free energy functional expressions can be easily constructed. Finally, the connectivity of linear molecule can be viewed as associating spheres, and the free energy functional follows naturally the associating theory.16 In the 3D-DFT calculations, fast Fourier transform to evaluate the convolutions of the free energy functional is employed. This algorithm can extensively save computational time and memory usage.17 In order to understand crystal nucleation comprehensively, the corresponding liquid−solid phase coexistence and the melting enthalpy are calculated with the new free energy functional, and the results are com-
pared with experimental values to test its reliability. Accordingly, density profiles in crystal cell and interface, and interfacial tension are predicted for fcc crystal orientation. II. THEORY AND EQUATIONS
For inhomogeneous system, DFT is applied to obtain the density distribution in the equilibrium state by minimizing the grand potential [ρ α (r)], or equivalently the intrinsic Helmholtz free energy F[ρ α (r)], as a functional of the density distribution ρ α (r). In general, the grand potential [ρ α (r)] can be expressed as the following form: [ρα (r)] = F [ρα (r)] + dr [ρα (r) (Vα (r) − μα )], α
(1) where α is a site in one molecule, μα is the chemical potential in the ensemble, Vα (r) is the external potential, and in bulk system, Vα (r) = 0. For polyatomic fluid, we begin by separating the total free energy into four components F [ρα (r)] = kB T drρα (r) ln ρα (r)3α − 1 α
+F [ρα (r)] + F att [ρα (r)] +F chain [ρα (r)] (2) hs
representing the contributions of ideal reference, hard-sphere repulsion, nonbonded van der Waals attractions and architecture, and molecular connectivity to the total free energy. Here kB is the Boltzmann constant, T is the absolute temperature, is the thermal de Broglie wavelength. The minimization of the grand potential functional leads to the Euler-Lagrange equation
−3
ρα (r) =
δ F hs [ρα (r)] + F att [ρα (r)] + F chain [ρα (r)] − βVα (r) exp βμα − β δρα (r)
with β = 1/kB T. The equilibrium 3D density profile of the system can be obtained by solving this equation. For hard-sphere contribution, the fundamental measure theory18 is postulated to have the form hs (4) βF [ρα (r)] = hs [nm (r)] dr, where the free energy for the hard sphere system hs [nm (r)] is a functional of the weighted densities of the system nm (r) with m = 0, 1, 2, 3, V 1, V 2, T . To adequately describe the solid, the tensor version of White Bear II for hard spheres6 is applied to modify the free energy hs [nm (r)] n1 n2 −nV 1 · nV 2 = −n0 ln(1−n3 )+ ϕ1 (n3 )+ϕ2 (n3 ) 1−n3 3 nV 2 ·nT ·nV 2 −n2 nV 2 ·nV 2 −tr n3T +n2 tr n2T × 16π (1−n3 )2
(3)
with ϕ1 (n3 ) = 1 +
2n3 − n23 + 2 (1 − n3 ) ln (1 − n3 ) , 3n3
(6)
2n3 − 3n23 + 2n33 + 2 (1 − n3 )2 ln (1 − n3 ) , 3n23 (7) where n0 , n1 , n2 , n3 are scalar weighted densities, nV 1 , nV 2 are vector weighted densities, and nT is tensor weighted density. These weighed densities are based on the deconvolution of the Mayer f-function19 for hard spheres, which is a function of the radii of both spheres, into a sum of weighed functions based on the individual particle radii.18 The weighted densities are written as (8) nm (r) = ρα r wα(m) r − r dr , ϕ2 (n3 ) = 1 −
α
(5)
where the weight functions wα(m) are given by wα(2) = δ(Rα − |r|), wα(3) = θ (Rα − |r|), wα(0) = wα(2) /4π Rα2 ,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-3
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
(2) 1) w (1) w(V 2) = rδ(Rα − |r|)/|r|, w(V α = wα /4π Rα , α (V 2) (T ) 2 = w α /4π Rα , and wα = rrδ(Rα − |r|)/r , respectively. Here δ(r)and θ (r) denote the Dirac delta function and Heaviside step function. Accordingly, the local chemical potential of hard-sphere repulsion can be calculated by δβF hs [ρα (r)] ∂hs (m) = w (r − r )dr . (9) δρα (r) ∂n m m
To construct the free energy functional due to attraction, several methods can be used, such as mean-field theory,20 perturbation theory,21, 22 and weighted density approximation.13 Furthermore, the more accurate weight functions can be obtained by a Heaviside step function with a temperaturedependent diameter.23 In the work, Fatt [ρ α (r)] is expressed with the weighted density approximation method13
N att βF [ρα (r)] = ρα (r) a att [ρ¯αα (r)] α
+
N 1 att a ρ¯αχ (r) 2 χ=1,χ=α
att ρ¯αχ (r ) att δa ωαχ (r−r )dr dr + ρχ (r ) δ ρ¯αχ (r )
The free energy contributions from the molecular connectivity are based on the association16, 28 N 1 βF chain [ρα (r)] = ρα r dr 2 α=1
{α} δ(|r − r | − σαα ) × 1−ln dr 2 4π σαα α ×yαα ( ξγ )ρα r , (15) where yαα ({ξ γ }) is the cavity correlation function between the sites α and α in one molecule,29 and σ αα = (σ α + σ α )/2 is the interaction diameter of two site. The notation {α} on the innermost summation of this functional indicates this sum is to be performed over all sites α , which are connected to the site α. The cavity correlation function is expressed as a set of weighed nonlocal densities ξ γ , which resemble the Rosenfeld hard sphere densities. ξ γ is written as ξγ (r) = dr ωα(γ ) (r − r )ρα (r ), γ = 0, 1, 2, 3 (16) α
with γ
ωα(γ ) (r) =
(10) att (r), in which ρ¯αχ (r) is defined via a weight function ωαχ att (r) r − r dr (11) ρ¯αχ (r) = ρα r ωαχ
and att ωαχ
(r) =
att (r) cαχ att (r) dr cαχ
σα θ (σα − |r|). 8σα3
(17)
Accordingly, the corresponding local chemical potential for molecular connectivity can be calculated by δβF chain [ρα (r)] δρα (r) {α}
,
(12)
att (r) is the direct correlation function of attractive where cαχ part between site α in one molecule and site χ in another att hs hs (r) = c(r) − cαχ (r), and cαχ (r) is the molecule; hence cαχ 24 hard-sphere reference term. c(r) can be obtained by solving the reference interaction site model integral equation25 h(r) = dr dr ω(|r − r |)c(|r − r |)[ω(r ) + ρh(r )],
(13) where h(r) is the total correlation function, ω(r) is the intramolecular correlation function, the details are given elsewhere.26 To solve the integral equation, the hypernetted chain approximation is applied, and the force field parameters27 are: εc-c /kB = 29.0 K, εo-o /kB = 83.0 K, εc-o /kB = 49.06 K, σ c = 2.785 Å, σ o = 3.064 Å, lc-o = 1.161 Å, and θ = 174.2◦ . The local chemical potential for the attractive interaction is then written as δβF att [ρα (r)] δρα (r) N δa att [ρ¯αχ (r )] att ω a att [ρ¯αχ (r)]+ ρα (r ) . = (r−r )dr αχ δ ρ¯αχ (r ) χ=1 (14)
1 = [1 − zαα ({ξγ } − tαα (r)] 2 α {κ} N ∂zκκ 1 (γ ) − ρκ (r) ω (r − r ) dr 2 ∂ξγ α γ κ κ=1 −
1 2
dr
{α}
ρα (r)
α
∂tαα (r) δ(|r − r | − σαα ) , 2 ∂nαα 4π σαα
(18) δ(|r−r |−σ ) where nαα (r) = dr 4πσ 2 αα ρα r , zαα ({ξ γ }) = αα ln yαα ({ξ γ }) and tαα (r) = ln nαα (r). On the other hand, the crystal phase should be regarded as a self-structured fluid, for which, even in the absence of any external potential, the minimum of [ρ(r)] is not achieved by a homogeneous density ρ(r) = ρ 0 , but rather by a strongly modulated density ρ(r), with the symmetry of the crystal lattice. The bulk number density of the solid phase is an inhomogeneous function, which is characterized by the lattice vectors {R}:6 2 ρsol (r) = ρsol + (θ/π )3/2 e−θ|R−r| , (19)
R=0
where ρ sol denotes the averaged density of the solid under equilibrium state, and the parameter θ is calculated by a restricted minimization of free energy functional of crystal. The
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-4
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
lattice vectors are R = h1 b1 + h2 b2 + h3 b3 for integer h1 , h2 , h3 , and b1 , b2 , b3 are primitive translation vectors of the lattice. Using Eq. (3), one can determine the interfacial density profile by a full minimization of the free energy in threedimensional real space. For such a complicated problem, the initial input is important to obtain a stable solution. In the solid–liquid interfacial region, the initial density ρ(r) can be assumed as the following form:3 ρ(r) = ρliq + (ρsol − ρliq )0 (z) 2 G (z)(θ/π )3/2 e−θ|R−r| +
(20)
R=0
in which ρ liq denotes the density of liquid under equilibrium state, and G (z) characterizes the interface varying from the crystal phase to the liquid phase ⎧ 1, |z| < z0 ⎪ ⎨ z−z0 1 G (z)= 2 1 + cos π zG , z0 ≤ |z| ≤ zG , (21) ⎪ ⎩ 0, |z| > zG
zG = |zG − z0 | = (R1 /R)υ z,
0 ≤ υ ≤ 1, R ≥ R1 , (22) where R1 is the magnitude of the smallest nonzero lattice vectors, the parameter z is the width of the interface, and the parameter υ controls the rate of broadening of the solid density peaks. z and υ can be obtained by a restricted minimization of free energy.30, 31 In the numerical calculations, the system is discretized into 256Lx × 256Ly × 256Lz grids within a rectangular cuboid box with side lengths Lx , Ly , and Lz . The widths Lx and Ly are chosen to be a few lattices spacing and Lz must be large enough to describe interfaces. As a result, the length of Lx and Ly is about 6 Å, and the length of Lz is about 25 Å. To calculate the bulk thermodynamic properties of solid CO2 from Eq. (2), we employ a simple Picard iteration procedure. One unit cell with periodic boundary conditions is used in these calculations. The initial configuration is constructed using Eq. (19) to represent a sum of discrete delta functions at the lattice sites in the unit cell. Based on the density profile, the free energy is calculated. At each average density, the calculation is repeated to minimize the free energy with respect to the lattice spacing, and an optimized θ is achieved. We determine the final density profile through the full minimization of the free energy in three-dimensional real space. To solve Eq. (3), we employ δFex (r)/δρ(r) but not ρ(r) during the iterations, which is convenient to the convergence. The iteration procedure is repeated until the average fractional difference over any 3D grid point between the old and the new δFex (r)/δρ(r) is less than 1.0 × 10−4 , and the vacancy concentration is nvac < 10−2 .
FIG. 1. FCC crystal structure of CO2 under moderate pressure.
node stands for a molecule of CO2 . As shown in Fig. 1, CO2 molecules are located at eight corners and six face centers of the cubic. In the following discussion, the structure and properties are calculated at vapor-liquid−solid triple point (T = 216.96K). Fig. 2 presents the free energy profiles for liquid and crystal phases. As can be seen, the free energy densities of liquid and solid increase with the increasing density. Most importantly, the free energy of solid is lower than that of liquid in the region of high density, leading to formalism of crystal. Coexistence densities are obtained from the free energy curves by the standard common tangent construction. Compared with the available experimental data in parentheses,35 the equilibrium liquid and solid densities are 1197 (1178) kg/m3 and 1492 (1512) kg/m3 , respectively. Both of the results show that the theoretical predictions are quite close to the actual data. This example strongly suggests that the present density functional approach accurately captures the enthalpic and entropic
III. RESULTS AND DISCUSSION
There are several kinds of crystal structures of CO2 at different pressures.32–34 At moderate pressure, the crystal structure is fcc lattice, where each face has five nodes and each
FIG. 2. Free energy as a function of density. Two points denote the equilibrium densities of liquid and fcc crystal, respectively.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-5
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
contributions of polyatomic CO2 , and its applications to crystal and crystal−fluid interface will be comfortable. Fig. 3 depicts the calculated density profiles of carbon and oxygen in (100) crystal plane. According to the lattice vectors, the peaks correspond to the particles located about their lattice positions. In one unit cell, there are five peaks in fcc. Comparing the two density fields, one sees that the density peak of carbon is much lower than that of oxygen, but slightly higher than its half value. Since the carbon atom is located in the center CO2 molecule, the contribution of molecular connectivity to the free energy functional of carbon is slightly higher than that of one oxygen, leading to larger entropic interactions. According to the density distribution of carbon and oxygen, we calculate the distribution of free energy, entropy, and enthalpy in one unit cell. The results are plotted depicted in Fig. 4. It is shown that the energy distribution is coherently related to the density distribution. The locations of the highest peaks in the energy distribution curves are the locations of the highest local densities. At high density, the hard-sphere repulsion plays the dominate role, whereas at low density, the attractive contribution is more important to make free energy lower or even negative. The distributions of entropy and enthalpy are applied to calculate the molar enthalpy of fusion. In order to investigate the thermodynamic properties of the crystal-fluid interface, it is important to calculate the
FIG. 4. Energy distribution of CO2 in (100) crystal plane. (a) Free energy, (b) entropy, and (c) enthalpy.
FIG. 3. Local density of carbon (a) and oxygen (b) in (100) crystal plane of fcc coexisting with liquid at the triple point.
density distribution of interface. The density profiles of different sites at coexistence are constructed under the constraint of solid and liquid densities. The density distributions of interface are computed by minimizing the grand potential. The results in equilibrium (100) interface are given in Fig. 5. The considerable variations in three-dimensional structure can be seen clearly. While the corresponding
equilibrium planaraveraged density profiles ρ(z) = A−1 ρ(r)dxdy are presented in Fig. 6. After careful inspection of the profiles in Fig. 5, one sees that the peaks of density profile are slowly reduced from
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-6
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
FIG. 5. Two-dimensional cuts of density profiles for carbon (a) and oxygen (b) at the (100) fcc–liquid interface.
the crystal side to the liquid side. Near the liquid side, the density variations are still obvious with narrow interfaces. In Fig. 6, the peak of carbon is slightly higher than half that of oxygen until converged to the liquid density. The result corresponds to the structure of crystal.
FIG. 6. The averaged interfacial density of carbon and oxygen at z direction.
FIG. 7. The platform of density profile of carbon. (a) (100) – liquid interface, (b) (110) – liquid interface, (c) (111) – liquid interface.
Fig. 7 reflects the density distribution of carbon in the equilibrium (100), (110), and (111) interfaces, respectively. It is shown that (100) crystal-liquid interface has the highest peak with 8 interfacial layers and presents a square grid distribution, (111) crystal−liquid interface has the same 8 distribution layers but presents a hexagonal grid distribution, (110) crystal−liquid interface has the highest peak with 4 interfacial layers and presents a rectangular grid distribution. According to the density distribution of interfaces, the interfacial tensions for different orientations are predicted. The results of γ 100 , γ 110 , and γ 111 are 12.27 mN/m, 11.85 mN/m, and 11.09 mN/m, respectively. Unfortunately, the corresponding experimental or computational values are now unavailable to test their reliability. Whereas we can calculate the melting enthalpy using the entropies of the equilibrium liquid and
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-7
Chang, Mi, and Zhong
crystal with H = TS, and compare it with the corresponding experimental value (8.347 kJ/mol).36 The theoretical result is 7.138 kJ/mol, showing a reasonable prediction of the theory. IV. CONCLUSION
We have presented a three-dimensional density functional approach to quantitatively describe the thermodynamic properties of crystal−liquid interface of CO2 . In the theoretical framework, the fundamental measure theory for hardsphere repulsion of solid is applied as the basis, and the contributions of enthalpic interactions due to dispersive attractions and of entropic interactions arising from the molecular architecture are integrated in the free energy functional. As a result, the total free energy functional contains the contributions of ideal reference, hard-sphere repulsion, nonbonded van der Waals attractions and architecture, and molecular connectivity. In particular, the third term is derived from the reference interaction site model integral equation, which is particularly suitable for structure description of real fluids. The liquid−solid phase equilibrium and the melting enthalpy at triple point has been correctly predicted, showing that the present theoretical model is quantitatively accurate to describe the inhomogeneous properties of polyatomic fluids. The density morphology of fcc crystal and of liquid−solid interface have been presented, but the accuracy of the predicted corresponding interfacial tensions is unknown, since the corresponding experimental or computational values are still unavailable. More importantly, we expect the present theoretical approach can be extended to deal with the structure and properties of crystallization in more complicated fluids, such as polymer melts. ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China (Grant Nos. 21076006 and 21276010). 1 H.
Löwen, Lecture Notes, 3rd Warsaw School of Statistical Physics (Warsaw University Press, 2010), p. 87. 2 P. Harrowell, J. Phys.: Condens. Matter 22(36), 364106 (2010). 3 W. Curtin, Phys. Rev. B 39(10), 6775 (1989).
J. Chem. Phys. 140, 204706 (2014) 4 Y.
Rosenfeld, M. Schmidt, H. Löwen, and P. Tarazona, Phys. Rev. E 55(4), 4245 (1997). 5 P. Tarazona, Phys. Rev. Lett. 84(4), 694 (2000). 6 M. Oettel, S. Görig, A. Härtel, H. Löwen, M. Radu, and T. Schilling, Phys. Rev. E 82(5), 051404 (2010). 7 G. Kahl and H. Löwen, J. Phys.: Condens. Matter 21(46), 464101 (2009). 8 K. Elder, M. Katakowski, M. Haataja, and M. Grant, Phys. Rev. Lett. 88(24), 245701 (2002). 9 K. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, Phys. Rev. B 75(6), 064107 (2007). 10 A. Jaatinen, C. Achim, K. Elder, and T. Ala-Nissila, Phys. Rev. E 80(3), 031602 (2009). 11 Y. Rosenfeld, J. Phys.: Condens. Matter 8(50), L795 (1996). 12 X. Wang, J. Mi, and C. Zhong, J. Chem. Phys. 138(16), 164704 (2013). 13 S.-C. Kim and S. H. Lee, J. Phys.: Condens. Matter 16(36), 6365 (2004). 14 Y.-X. Yu, G.-H. Gao, and X.-L. Wang, J. Phys. Chem. B 110(29), 14418 (2006). 15 A. L. Frischknecht, V. Padmanabhan, and M. E. Mackay, J. Chem. Phys. 136(16), 164904 (2012). 16 S. Tripathi and W. G. Chapman, Phys. Rev. Lett. 94(8), 087801 (2005). 17 M. G. Knepley, D. A. Karpeev, S. Davidovits, R. S. Eisenberg, and D. Gillespie, J. Chem. Phys. 132(12), 124101 (2010). 18 Y. Rosenfeld, Phys. Rev. Lett. 63(9), 980 (1989). 19 H. T. Davis, Statistical Mechanics of Phases, Interfaces, and Thin Films (VCH, New York, 1996). 20 F.-Q. You, Y.-X. Yu, and G.-H. Gao, J. Chem. Phys. 123(11), 114705 (2005). 21 Y. Rosenfeld, J. Chem. Phys. 98(10), 8126 (1993). 22 B. Peng and Y.-X. Yu, J. Chem. Phys. 131(13), 134703 (2009). 23 Y.-X. Yu, J. Chem. Phys. 131(2), 024704 (2009). 24 J. K. Percus and G. J. Yevick, Phys. Rev. 110(1), 1 (1958). 25 D. Chandler, H. C. Anderson, and J. D. Weeks, “Optimized Cluster Expansions for Classical Fluids. 3. Applications to Ionic Solutions and Simple Liquids” (Department of Chemistry, California University San Diego, La Jolla, 1972). 26 K. G. Honnell, J. G. Curro, and K. S. Schweizer, Macromolecules 23(14), 3496 (1990). 27 Y. Zhang, J. Yang, and Y.-X. Yu, J. Phys. Chem. B 109(27), 13375 (2005). 28 S. Tripathi and W. G. Chapman, J. Chem. Phys. 122, 094506 (2005). 29 E. McGarrity, A. Frischknecht, and M. Mackay, J. Chem. Phys. 128(15), 154904 (2008). 30 A. Härtel, M. Oettel, R. E. Rozas, S. U. Egelhaaf, J. Horbach, and H. Löwen, Phys. Rev. Lett. 108(22), 226101 (2012). 31 R. Ohnesorge, H. Löwen, and H. Wagner, Phys. Rev. E 50(6), 4801 (1994). 32 B. Olinger, J. Chem. Phys. 77(12), 6255 (1982). 33 K. Aoki, H. Yamawaki, M. Sakashita, Y. Gotoh, and K. Takemura, Science 263(5145), 356–358 (1994). 34 C. Yoo, H. Cynn, F. Gygi, G. Galli, V. Iota, M. Nicol, S. Carlson, D. Häusermann, and C. Mailhiot, Phys. Rev. Lett. 83(26), 5527 (1999). 35 A. Yokozeki, Fluid Phase Equilibria 222–223, 55 (2004). 36 O. Maass and W. Barnes, Proc. R. Soc. London Ser. A 111(757), 224 (1926).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
THE JOURNAL OF CHEMICAL PHYSICS 140, 204706 (2014)
Density functional theory for carbon dioxide crystal Yiwen Chang, Jianguo Mi,a) and Chongli Zhong State Key Laboratory of Organic-Inorganic Composites, Beijing University of Chemical Technology, Beijing 100029, China
(Received 27 February 2014; accepted 5 May 2014; published online 27 May 2014) We present a density functional approach to describe the solid−liquid phase transition, interfacial and crystal structure, and properties of polyatomic CO2 . Unlike previous phase field crystal model or density functional theory, which are derived from the second order direct correlation function, the present density functional approach is based on the fundamental measure theory for hard-sphere repulsion in solid. More importantly, the contributions of enthalpic interactions due to the dispersive attractions and of entropic interactions arising from the molecular architecture are integrated in the density functional model. Using the theoretical model, the predicted liquid and solid densities of CO2 at equilibrium triple point are in good agreement with the experimental values. Based on the structure of crystal-liquid interfaces in different planes, the corresponding interfacial tensions are predicted. Their respective accuracies need to be tested. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4878413] I. INTRODUCTION
Crystallization phenomena are ubiquitous in both natural and artificially synthesized systems. During the nucleation process, the microstructure of new crystal is strongly influenced by the property of bulk fluid and the complex spatial patterns of solid. In order to understand crystallization from microscopic scale, some essential issues, such as the solid−liquid phase transition, the structure and energy of crystal and interface, need to be clarified. However, it has so far been difficult to obtain these quantitative descriptions theoretically. The difficulty becomes more serious when dealing with polyatomic fluid−solid systems. For simple monatomic fluids, only enthalpic interactions affect these properties. Whereas for polyatomic fluids, although enthalpic interactions are usually dominant, entropic interactions arising from the molecular architecture can play a critical role in dictating the structure and properties. A successful theoretical model should consider their cooperative contributions. There are two different routes to develop crystallization theory. One is derived from liquid-state theory, and another is developed from field theory. Due to the success of liquid-state theory, there has been much effort in describing the solid state with liquid properties, leading to the development of density functional theory (DFT). In recent years, investigations in condensed matter have given unambiguous evidence that DFT is a very versatile and powerful tool to study the structural and thermodynamic properties of crystallization.1, 2 Meanwhile, the phase-field crystal (PFC) model developed from phasefield theory has proven to be prospective in a variety of problems, such as interfaces,3 polycrystalline pattern formation,4, 5 or crystal nucleation and growth.6–8 Of later years, PFC combined with classical DFT has been developed and applied to solve this problem in a highly efficient way.9, 10 It has shown a) Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2014/140(20)/204706/7/$30.00
that the combined theory is an ideal tool to study nucleation process. In the above theoretical models, the energy expressions are unexceptionally derived from the direct correlation function with a second-order expansion form. It is well-known that such approximation gives poor results even for describing monatomic fluids in equilibrium state. The deficiency is mainly due to the inaccurate direct correlation function inside the hard-sphere. Without accurate free energy expression, both the PFC model and the DFT are a kind of coarse-grained method, which focus mainly on phenomenological description. On the other hand, the development of the free energy functional for hard-sphere is presented through its milestones in the weighted density approximation and the fundamental measure theory.11 The extensions of these approaches to hard-sphere mixtures include the treatment of polydisperse systems. Very recently, we have presented a three-dimensional (3D) DFT approach to describe the crystal−liquid interfaces of Lennard-Jones fluid.12 In the theoretical model, the fundamental measure theory is applied to describe the free energy functional of hard-sphere repulsion, and a weighted density method13 is used to describe the free energy contribution arising from the van der Waals attractions. The liquid-solid equilibrium, density profiles within crystal cell and at liquid−solid interface, and interfacial tensions were calculated for facecentered-cubic (fcc) and body-centered-cubic (bcc) nucleus. The results are in good agreement with simulation data, showing that the model is quantitatively accurate in describing nucleation thermodynamics of the Lennard-Jones fluid. In the recent years, DFT was further developed and applied to research polyatomic fluids or polymer systems. For instances, a Wertheim’s first-order thermodynamic perturbation theory (TPT1) based DFT was proposed for chain conformation and bond orientation correlation function of polyatomic fluid;14 Frischknecht et al.15 employed a DFT based on interfacial statistical associating fluid theory (SAFT)
140, 204706-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-2
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
equation to probe the first-order phase transition in athermal blend of nanoparticle and homopolymer near a wall. The developments of DFT are conducive to obtain the structure and properties for polyatomic fluids more accurately. The main purpose of this work is to extend the above 3DDFT approach to polyatomic fluids. CO2 is taken as an example. Compared with monatomic molecules, a CO2 molecule contains three sites, and its structure has a marginal deviation from the linear geometry. The contribution of enthalpy to the free energy functional is represented with the van der Waals interactions. Meanwhile, the entropic contribution is correlated to the architecture of CO2 molecule including bond length, bond angle, and connectivity. Similar to our previous work, the van der Waals interactions are split into hard-sphere repulsions and dispersive attractions, and the free energy functional for hard-sphere repulsions is still represented by the fundamental measure theory for solid phase. In particular, by solving the reference interaction site model integral equation, the nonbonded dispersive attractions, bond length, and bond angle are integrated into the direct correlation functions of different sites. Once the direct correlation functions are derived, the corresponding free energy functional expressions can be easily constructed. Finally, the connectivity of linear molecule can be viewed as associating spheres, and the free energy functional follows naturally the associating theory.16 In the 3D-DFT calculations, fast Fourier transform to evaluate the convolutions of the free energy functional is employed. This algorithm can extensively save computational time and memory usage.17 In order to understand crystal nucleation comprehensively, the corresponding liquid−solid phase coexistence and the melting enthalpy are calculated with the new free energy functional, and the results are com-
pared with experimental values to test its reliability. Accordingly, density profiles in crystal cell and interface, and interfacial tension are predicted for fcc crystal orientation. II. THEORY AND EQUATIONS
For inhomogeneous system, DFT is applied to obtain the density distribution in the equilibrium state by minimizing the grand potential [ρ α (r)], or equivalently the intrinsic Helmholtz free energy F[ρ α (r)], as a functional of the density distribution ρ α (r). In general, the grand potential [ρ α (r)] can be expressed as the following form: [ρα (r)] = F [ρα (r)] + dr [ρα (r) (Vα (r) − μα )], α
(1) where α is a site in one molecule, μα is the chemical potential in the ensemble, Vα (r) is the external potential, and in bulk system, Vα (r) = 0. For polyatomic fluid, we begin by separating the total free energy into four components F [ρα (r)] = kB T drρα (r) ln ρα (r)3α − 1 α
+F [ρα (r)] + F att [ρα (r)] +F chain [ρα (r)] (2) hs
representing the contributions of ideal reference, hard-sphere repulsion, nonbonded van der Waals attractions and architecture, and molecular connectivity to the total free energy. Here kB is the Boltzmann constant, T is the absolute temperature, is the thermal de Broglie wavelength. The minimization of the grand potential functional leads to the Euler-Lagrange equation
−3
ρα (r) =
δ F hs [ρα (r)] + F att [ρα (r)] + F chain [ρα (r)] − βVα (r) exp βμα − β δρα (r)
with β = 1/kB T. The equilibrium 3D density profile of the system can be obtained by solving this equation. For hard-sphere contribution, the fundamental measure theory18 is postulated to have the form hs (4) βF [ρα (r)] = hs [nm (r)] dr, where the free energy for the hard sphere system hs [nm (r)] is a functional of the weighted densities of the system nm (r) with m = 0, 1, 2, 3, V 1, V 2, T . To adequately describe the solid, the tensor version of White Bear II for hard spheres6 is applied to modify the free energy hs [nm (r)] n1 n2 −nV 1 · nV 2 = −n0 ln(1−n3 )+ ϕ1 (n3 )+ϕ2 (n3 ) 1−n3 3 nV 2 ·nT ·nV 2 −n2 nV 2 ·nV 2 −tr n3T +n2 tr n2T × 16π (1−n3 )2
(3)
with ϕ1 (n3 ) = 1 +
2n3 − n23 + 2 (1 − n3 ) ln (1 − n3 ) , 3n3
(6)
2n3 − 3n23 + 2n33 + 2 (1 − n3 )2 ln (1 − n3 ) , 3n23 (7) where n0 , n1 , n2 , n3 are scalar weighted densities, nV 1 , nV 2 are vector weighted densities, and nT is tensor weighted density. These weighed densities are based on the deconvolution of the Mayer f-function19 for hard spheres, which is a function of the radii of both spheres, into a sum of weighed functions based on the individual particle radii.18 The weighted densities are written as (8) nm (r) = ρα r wα(m) r − r dr , ϕ2 (n3 ) = 1 −
α
(5)
where the weight functions wα(m) are given by wα(2) = δ(Rα − |r|), wα(3) = θ (Rα − |r|), wα(0) = wα(2) /4π Rα2 ,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-3
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
(2) 1) w (1) w(V 2) = rδ(Rα − |r|)/|r|, w(V α = wα /4π Rα , α (V 2) (T ) 2 = w α /4π Rα , and wα = rrδ(Rα − |r|)/r , respectively. Here δ(r)and θ (r) denote the Dirac delta function and Heaviside step function. Accordingly, the local chemical potential of hard-sphere repulsion can be calculated by δβF hs [ρα (r)] ∂hs (m) = w (r − r )dr . (9) δρα (r) ∂n m m
To construct the free energy functional due to attraction, several methods can be used, such as mean-field theory,20 perturbation theory,21, 22 and weighted density approximation.13 Furthermore, the more accurate weight functions can be obtained by a Heaviside step function with a temperaturedependent diameter.23 In the work, Fatt [ρ α (r)] is expressed with the weighted density approximation method13
N att βF [ρα (r)] = ρα (r) a att [ρ¯αα (r)] α
+
N 1 att a ρ¯αχ (r) 2 χ=1,χ=α
att ρ¯αχ (r ) att δa ωαχ (r−r )dr dr + ρχ (r ) δ ρ¯αχ (r )
The free energy contributions from the molecular connectivity are based on the association16, 28 N 1 βF chain [ρα (r)] = ρα r dr 2 α=1
{α} δ(|r − r | − σαα ) × 1−ln dr 2 4π σαα α ×yαα ( ξγ )ρα r , (15) where yαα ({ξ γ }) is the cavity correlation function between the sites α and α in one molecule,29 and σ αα = (σ α + σ α )/2 is the interaction diameter of two site. The notation {α} on the innermost summation of this functional indicates this sum is to be performed over all sites α , which are connected to the site α. The cavity correlation function is expressed as a set of weighed nonlocal densities ξ γ , which resemble the Rosenfeld hard sphere densities. ξ γ is written as ξγ (r) = dr ωα(γ ) (r − r )ρα (r ), γ = 0, 1, 2, 3 (16) α
with γ
ωα(γ ) (r) =
(10) att (r), in which ρ¯αχ (r) is defined via a weight function ωαχ att (r) r − r dr (11) ρ¯αχ (r) = ρα r ωαχ
and att ωαχ
(r) =
att (r) cαχ att (r) dr cαχ
σα θ (σα − |r|). 8σα3
(17)
Accordingly, the corresponding local chemical potential for molecular connectivity can be calculated by δβF chain [ρα (r)] δρα (r) {α}
,
(12)
att (r) is the direct correlation function of attractive where cαχ part between site α in one molecule and site χ in another att hs hs (r) = c(r) − cαχ (r), and cαχ (r) is the molecule; hence cαχ 24 hard-sphere reference term. c(r) can be obtained by solving the reference interaction site model integral equation25 h(r) = dr dr ω(|r − r |)c(|r − r |)[ω(r ) + ρh(r )],
(13) where h(r) is the total correlation function, ω(r) is the intramolecular correlation function, the details are given elsewhere.26 To solve the integral equation, the hypernetted chain approximation is applied, and the force field parameters27 are: εc-c /kB = 29.0 K, εo-o /kB = 83.0 K, εc-o /kB = 49.06 K, σ c = 2.785 Å, σ o = 3.064 Å, lc-o = 1.161 Å, and θ = 174.2◦ . The local chemical potential for the attractive interaction is then written as δβF att [ρα (r)] δρα (r) N δa att [ρ¯αχ (r )] att ω a att [ρ¯αχ (r)]+ ρα (r ) . = (r−r )dr αχ δ ρ¯αχ (r ) χ=1 (14)
1 = [1 − zαα ({ξγ } − tαα (r)] 2 α {κ} N ∂zκκ 1 (γ ) − ρκ (r) ω (r − r ) dr 2 ∂ξγ α γ κ κ=1 −
1 2
dr
{α}
ρα (r)
α
∂tαα (r) δ(|r − r | − σαα ) , 2 ∂nαα 4π σαα
(18) δ(|r−r |−σ ) where nαα (r) = dr 4πσ 2 αα ρα r , zαα ({ξ γ }) = αα ln yαα ({ξ γ }) and tαα (r) = ln nαα (r). On the other hand, the crystal phase should be regarded as a self-structured fluid, for which, even in the absence of any external potential, the minimum of [ρ(r)] is not achieved by a homogeneous density ρ(r) = ρ 0 , but rather by a strongly modulated density ρ(r), with the symmetry of the crystal lattice. The bulk number density of the solid phase is an inhomogeneous function, which is characterized by the lattice vectors {R}:6 2 ρsol (r) = ρsol + (θ/π )3/2 e−θ|R−r| , (19)
R=0
where ρ sol denotes the averaged density of the solid under equilibrium state, and the parameter θ is calculated by a restricted minimization of free energy functional of crystal. The
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-4
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
lattice vectors are R = h1 b1 + h2 b2 + h3 b3 for integer h1 , h2 , h3 , and b1 , b2 , b3 are primitive translation vectors of the lattice. Using Eq. (3), one can determine the interfacial density profile by a full minimization of the free energy in threedimensional real space. For such a complicated problem, the initial input is important to obtain a stable solution. In the solid–liquid interfacial region, the initial density ρ(r) can be assumed as the following form:3 ρ(r) = ρliq + (ρsol − ρliq )0 (z) 2 G (z)(θ/π )3/2 e−θ|R−r| +
(20)
R=0
in which ρ liq denotes the density of liquid under equilibrium state, and G (z) characterizes the interface varying from the crystal phase to the liquid phase ⎧ 1, |z| < z0 ⎪ ⎨ z−z0 1 G (z)= 2 1 + cos π zG , z0 ≤ |z| ≤ zG , (21) ⎪ ⎩ 0, |z| > zG
zG = |zG − z0 | = (R1 /R)υ z,
0 ≤ υ ≤ 1, R ≥ R1 , (22) where R1 is the magnitude of the smallest nonzero lattice vectors, the parameter z is the width of the interface, and the parameter υ controls the rate of broadening of the solid density peaks. z and υ can be obtained by a restricted minimization of free energy.30, 31 In the numerical calculations, the system is discretized into 256Lx × 256Ly × 256Lz grids within a rectangular cuboid box with side lengths Lx , Ly , and Lz . The widths Lx and Ly are chosen to be a few lattices spacing and Lz must be large enough to describe interfaces. As a result, the length of Lx and Ly is about 6 Å, and the length of Lz is about 25 Å. To calculate the bulk thermodynamic properties of solid CO2 from Eq. (2), we employ a simple Picard iteration procedure. One unit cell with periodic boundary conditions is used in these calculations. The initial configuration is constructed using Eq. (19) to represent a sum of discrete delta functions at the lattice sites in the unit cell. Based on the density profile, the free energy is calculated. At each average density, the calculation is repeated to minimize the free energy with respect to the lattice spacing, and an optimized θ is achieved. We determine the final density profile through the full minimization of the free energy in three-dimensional real space. To solve Eq. (3), we employ δFex (r)/δρ(r) but not ρ(r) during the iterations, which is convenient to the convergence. The iteration procedure is repeated until the average fractional difference over any 3D grid point between the old and the new δFex (r)/δρ(r) is less than 1.0 × 10−4 , and the vacancy concentration is nvac < 10−2 .
FIG. 1. FCC crystal structure of CO2 under moderate pressure.
node stands for a molecule of CO2 . As shown in Fig. 1, CO2 molecules are located at eight corners and six face centers of the cubic. In the following discussion, the structure and properties are calculated at vapor-liquid−solid triple point (T = 216.96K). Fig. 2 presents the free energy profiles for liquid and crystal phases. As can be seen, the free energy densities of liquid and solid increase with the increasing density. Most importantly, the free energy of solid is lower than that of liquid in the region of high density, leading to formalism of crystal. Coexistence densities are obtained from the free energy curves by the standard common tangent construction. Compared with the available experimental data in parentheses,35 the equilibrium liquid and solid densities are 1197 (1178) kg/m3 and 1492 (1512) kg/m3 , respectively. Both of the results show that the theoretical predictions are quite close to the actual data. This example strongly suggests that the present density functional approach accurately captures the enthalpic and entropic
III. RESULTS AND DISCUSSION
There are several kinds of crystal structures of CO2 at different pressures.32–34 At moderate pressure, the crystal structure is fcc lattice, where each face has five nodes and each
FIG. 2. Free energy as a function of density. Two points denote the equilibrium densities of liquid and fcc crystal, respectively.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-5
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
contributions of polyatomic CO2 , and its applications to crystal and crystal−fluid interface will be comfortable. Fig. 3 depicts the calculated density profiles of carbon and oxygen in (100) crystal plane. According to the lattice vectors, the peaks correspond to the particles located about their lattice positions. In one unit cell, there are five peaks in fcc. Comparing the two density fields, one sees that the density peak of carbon is much lower than that of oxygen, but slightly higher than its half value. Since the carbon atom is located in the center CO2 molecule, the contribution of molecular connectivity to the free energy functional of carbon is slightly higher than that of one oxygen, leading to larger entropic interactions. According to the density distribution of carbon and oxygen, we calculate the distribution of free energy, entropy, and enthalpy in one unit cell. The results are plotted depicted in Fig. 4. It is shown that the energy distribution is coherently related to the density distribution. The locations of the highest peaks in the energy distribution curves are the locations of the highest local densities. At high density, the hard-sphere repulsion plays the dominate role, whereas at low density, the attractive contribution is more important to make free energy lower or even negative. The distributions of entropy and enthalpy are applied to calculate the molar enthalpy of fusion. In order to investigate the thermodynamic properties of the crystal-fluid interface, it is important to calculate the
FIG. 4. Energy distribution of CO2 in (100) crystal plane. (a) Free energy, (b) entropy, and (c) enthalpy.
FIG. 3. Local density of carbon (a) and oxygen (b) in (100) crystal plane of fcc coexisting with liquid at the triple point.
density distribution of interface. The density profiles of different sites at coexistence are constructed under the constraint of solid and liquid densities. The density distributions of interface are computed by minimizing the grand potential. The results in equilibrium (100) interface are given in Fig. 5. The considerable variations in three-dimensional structure can be seen clearly. While the corresponding
equilibrium planaraveraged density profiles ρ(z) = A−1 ρ(r)dxdy are presented in Fig. 6. After careful inspection of the profiles in Fig. 5, one sees that the peaks of density profile are slowly reduced from
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-6
Chang, Mi, and Zhong
J. Chem. Phys. 140, 204706 (2014)
FIG. 5. Two-dimensional cuts of density profiles for carbon (a) and oxygen (b) at the (100) fcc–liquid interface.
the crystal side to the liquid side. Near the liquid side, the density variations are still obvious with narrow interfaces. In Fig. 6, the peak of carbon is slightly higher than half that of oxygen until converged to the liquid density. The result corresponds to the structure of crystal.
FIG. 6. The averaged interfacial density of carbon and oxygen at z direction.
FIG. 7. The platform of density profile of carbon. (a) (100) – liquid interface, (b) (110) – liquid interface, (c) (111) – liquid interface.
Fig. 7 reflects the density distribution of carbon in the equilibrium (100), (110), and (111) interfaces, respectively. It is shown that (100) crystal-liquid interface has the highest peak with 8 interfacial layers and presents a square grid distribution, (111) crystal−liquid interface has the same 8 distribution layers but presents a hexagonal grid distribution, (110) crystal−liquid interface has the highest peak with 4 interfacial layers and presents a rectangular grid distribution. According to the density distribution of interfaces, the interfacial tensions for different orientations are predicted. The results of γ 100 , γ 110 , and γ 111 are 12.27 mN/m, 11.85 mN/m, and 11.09 mN/m, respectively. Unfortunately, the corresponding experimental or computational values are now unavailable to test their reliability. Whereas we can calculate the melting enthalpy using the entropies of the equilibrium liquid and
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
204706-7
Chang, Mi, and Zhong
crystal with H = TS, and compare it with the corresponding experimental value (8.347 kJ/mol).36 The theoretical result is 7.138 kJ/mol, showing a reasonable prediction of the theory. IV. CONCLUSION
We have presented a three-dimensional density functional approach to quantitatively describe the thermodynamic properties of crystal−liquid interface of CO2 . In the theoretical framework, the fundamental measure theory for hardsphere repulsion of solid is applied as the basis, and the contributions of enthalpic interactions due to dispersive attractions and of entropic interactions arising from the molecular architecture are integrated in the free energy functional. As a result, the total free energy functional contains the contributions of ideal reference, hard-sphere repulsion, nonbonded van der Waals attractions and architecture, and molecular connectivity. In particular, the third term is derived from the reference interaction site model integral equation, which is particularly suitable for structure description of real fluids. The liquid−solid phase equilibrium and the melting enthalpy at triple point has been correctly predicted, showing that the present theoretical model is quantitatively accurate to describe the inhomogeneous properties of polyatomic fluids. The density morphology of fcc crystal and of liquid−solid interface have been presented, but the accuracy of the predicted corresponding interfacial tensions is unknown, since the corresponding experimental or computational values are still unavailable. More importantly, we expect the present theoretical approach can be extended to deal with the structure and properties of crystallization in more complicated fluids, such as polymer melts. ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China (Grant Nos. 21076006 and 21276010). 1 H.
Löwen, Lecture Notes, 3rd Warsaw School of Statistical Physics (Warsaw University Press, 2010), p. 87. 2 P. Harrowell, J. Phys.: Condens. Matter 22(36), 364106 (2010). 3 W. Curtin, Phys. Rev. B 39(10), 6775 (1989).
J. Chem. Phys. 140, 204706 (2014) 4 Y.
Rosenfeld, M. Schmidt, H. Löwen, and P. Tarazona, Phys. Rev. E 55(4), 4245 (1997). 5 P. Tarazona, Phys. Rev. Lett. 84(4), 694 (2000). 6 M. Oettel, S. Görig, A. Härtel, H. Löwen, M. Radu, and T. Schilling, Phys. Rev. E 82(5), 051404 (2010). 7 G. Kahl and H. Löwen, J. Phys.: Condens. Matter 21(46), 464101 (2009). 8 K. Elder, M. Katakowski, M. Haataja, and M. Grant, Phys. Rev. Lett. 88(24), 245701 (2002). 9 K. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, Phys. Rev. B 75(6), 064107 (2007). 10 A. Jaatinen, C. Achim, K. Elder, and T. Ala-Nissila, Phys. Rev. E 80(3), 031602 (2009). 11 Y. Rosenfeld, J. Phys.: Condens. Matter 8(50), L795 (1996). 12 X. Wang, J. Mi, and C. Zhong, J. Chem. Phys. 138(16), 164704 (2013). 13 S.-C. Kim and S. H. Lee, J. Phys.: Condens. Matter 16(36), 6365 (2004). 14 Y.-X. Yu, G.-H. Gao, and X.-L. Wang, J. Phys. Chem. B 110(29), 14418 (2006). 15 A. L. Frischknecht, V. Padmanabhan, and M. E. Mackay, J. Chem. Phys. 136(16), 164904 (2012). 16 S. Tripathi and W. G. Chapman, Phys. Rev. Lett. 94(8), 087801 (2005). 17 M. G. Knepley, D. A. Karpeev, S. Davidovits, R. S. Eisenberg, and D. Gillespie, J. Chem. Phys. 132(12), 124101 (2010). 18 Y. Rosenfeld, Phys. Rev. Lett. 63(9), 980 (1989). 19 H. T. Davis, Statistical Mechanics of Phases, Interfaces, and Thin Films (VCH, New York, 1996). 20 F.-Q. You, Y.-X. Yu, and G.-H. Gao, J. Chem. Phys. 123(11), 114705 (2005). 21 Y. Rosenfeld, J. Chem. Phys. 98(10), 8126 (1993). 22 B. Peng and Y.-X. Yu, J. Chem. Phys. 131(13), 134703 (2009). 23 Y.-X. Yu, J. Chem. Phys. 131(2), 024704 (2009). 24 J. K. Percus and G. J. Yevick, Phys. Rev. 110(1), 1 (1958). 25 D. Chandler, H. C. Anderson, and J. D. Weeks, “Optimized Cluster Expansions for Classical Fluids. 3. Applications to Ionic Solutions and Simple Liquids” (Department of Chemistry, California University San Diego, La Jolla, 1972). 26 K. G. Honnell, J. G. Curro, and K. S. Schweizer, Macromolecules 23(14), 3496 (1990). 27 Y. Zhang, J. Yang, and Y.-X. Yu, J. Phys. Chem. B 109(27), 13375 (2005). 28 S. Tripathi and W. G. Chapman, J. Chem. Phys. 122, 094506 (2005). 29 E. McGarrity, A. Frischknecht, and M. Mackay, J. Chem. Phys. 128(15), 154904 (2008). 30 A. Härtel, M. Oettel, R. E. Rozas, S. U. Egelhaaf, J. Horbach, and H. Löwen, Phys. Rev. Lett. 108(22), 226101 (2012). 31 R. Ohnesorge, H. Löwen, and H. Wagner, Phys. Rev. E 50(6), 4801 (1994). 32 B. Olinger, J. Chem. Phys. 77(12), 6255 (1982). 33 K. Aoki, H. Yamawaki, M. Sakashita, Y. Gotoh, and K. Takemura, Science 263(5145), 356–358 (1994). 34 C. Yoo, H. Cynn, F. Gygi, G. Galli, V. Iota, M. Nicol, S. Carlson, D. Häusermann, and C. Mailhiot, Phys. Rev. Lett. 83(26), 5527 (1999). 35 A. Yokozeki, Fluid Phase Equilibria 222–223, 55 (2004). 36 O. Maass and W. Barnes, Proc. R. Soc. London Ser. A 111(757), 224 (1926).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 64.106.42.43 On: Wed, 01 Oct 2014 10:35:39
