Classical Density Functional Theory applied to the solid state
CClassical Density Functional Theory applied to the solid state
James F. Lutsko ∗ and C´edric Schoonen Center for Nonlinear Phenomena and Complex Systems CP 231,Universit´e Libre de Bruxelles, Blvd. du Triomphe, 1050 Brussels, Belgium (Dated: October 1, 2020)
Abstract
The standard model of classical Density Functional Theory for pair potentials consists of a hard-sphere functional plus a mean-field term accounting for long ranged attraction. However, mostimplementations using sophisticated Fundamental Measure hard-sphere functionals suffer frompotential numerical instabilities either due to possible instabilities in the functionals themselvesor due to implementations that mix real- and Fourier-space components inconsistently. Here,we present a new implementation based on a demonstrably stable hard-sphere functional that isimplemented in a completely consistent manner. The present work does not depend on approximatespherical integration schemes and so is much more robust than previous algorithms. The methodsare illustrated by calculating phase diagrams for the solid state using the standard Lennard-Jonespotential as well as a new class of potentials recently proposed by Wang et al (Phys. Chem. Chem.Phys. 22, 10624 (2020)). The latter span the range from potentials for small molecules to thoseappropriate to colloidal systems simply by varying a parameter. We verify that cDFT is able tosemi-quantitatively reproduce the phase diagram in all cases. We also show that for these problemscomputationally cheap Gaussian approximations are nearly as good as full minimization based onfinite differences. ∗ a r X i v : . [ phy s i c s . c o m p - ph ] S e p . INTRODUCTION Classical Density Functional Theory (cDFT) is an exact theory that has become a ver-satile tool for the studying the properties of inhomogeneous systems such as fluid inter-faces and solids at length scales going down to the molecular level[1, 2]. Conceptually,cDFT calculations involve the minimization of a functional of the local density resultingin both the equilibrium density distribution and the free energy of the system. Notablerecent applications include the description of wetting phenomena[3], the calculation of hy-dration free-energies and microscopic structure of molecular solutes[4] and the descriptionof crystallization pathways[5]. An important part of this utility lies in the highly-developeddescription of correlations due to excluded-volume effects that are arise whenever moleculesinteract via potentials having divergent short-ranged repulsion as is the case, e.g., in simplefluids and is captured in such commonly used models as the Lennard-Jones, Stillinger-Weberand hard-core Yukawa potentials. This capability is most highly developed in models basedon Fundamental Measure Theory which was first introduced by Rosenfeld[6], inspired byexact results due to Percus[7–9], and that has been steadily developed over the last 30 years(see, e.g., Ref. 10).Despite this progress, several issues have slowed the application of cDFT to the most chal-lenging of inhomogeneous systems, namely the solid phase where the local density varies bymany orders of magnitude over distances of the molecular diameter. Indeed, the applica-tion to solids is sufficiently challenging that until recently such applications as existed werebased on simplifications such as the modeling of the local density as a sum of Gaussianprofiles. The application of the full machinery of cDFT to the solid phase only began withthe work of Oettel and coworkers in 2010[11] and there remain relatively few applications.Perhaps the main reason for this is that these early calculations seemed to be very delicateand plagued by numerical instabilities. This has led, e.g., to calculations only being possi-ble at fixed particle number rather than at fixed chemical potential as is more natural incDFT(see, e.g. Ref. [12]). Recently, it has been proposed that these difficulties are traceableto instabilities inherent in some of the most popular models (based on the so-called WhiteBear FMT functionals) as well as to certain inconsistencies in the implementation of thecalculations (as explained in more detail below)[13, 14]. The solution to the first problemis the careful selection of models that are provably stable which, while more heuristic than2he most advanced models, are sufficiently accurate for many purposes. One of the mainpurposes of the present work is to present a solution to the second problem: namely, a robustimplementation that is free of instabilities.The second goal of this work is to examine the transferability of the cDFT model - i.e.its robustness when applied to different interaction potentials. The recent work on solidshas mostly focused on the Lennard-Jones potential and it relatively good results have beenreported, compared to simulation. Here, we also examine a new class of potentials introducedby Wang et al[15] who thoroughly characterized their vapor-liquid-solid phase diagrams andso provide an excellent test case. We focus on the two examples studied in detail in thatwork: namely, an analog of the Lennard-Jones potential and a second that models colloids.Colloids and simple fluids have qualitatively different phase diagrams and it is therefore ofinterest to verify to what extent cDFT is able to describe such variations.In the following, we first review the basic elements of the standard cDFT model con-sisting of the sum of an ideal gas contribution, an FMT hard-sphere contribution and amean-field term that captures the details of the potential. We also present our numericalimplementation of this model. In the next Section, we present our results first for the theLennard-Jones and WHDF potential parameterized as described above for both simple flu-ids and colloids. We have calculated the vapor-liquid-solid phase diagrams, the liquid-vaporsurface tensions and we also present some details of the solid phase such as the asymmetryof the density distributions,the vacancy concentration and the difference between the HCPand FCC phases. We end with some conclusions.
II. THEORYA. The standard cDFT model
The fundamental quantity in cDFT is the local density ρ ( r ). It is important to em-phasize that for an equilibrium system this is identical to the as the one-body probabilitydistribution and as such is a microscopic quantity that involves no coarse-graining or otherapproximations and, as such, it gives a description of the equilibrium density distributionvalid down to the smallest length scales. In cDFT, the local density of an equilibrium systemis determined by minimizing a functional denoted here as Λ [ ρ ]. The formal theory underly-3ng cDFT assures us that such a functional exists, that it is unique and that when evaluatedat the equilibrium density, it is equal to the grand-canonical free energy of the system, Ω.The ”standard” cDFT model is written asΛ [ ρ ] = F (id) [ ρ ] + F (HS) [ ρ ] + 12 (cid:90) ρ ( r ) ρ ( r ) w att ( r , r ) d r d r + (cid:90) ρ ( r ) ( φ ( r ) − (cid:101) µ ) d r (1)where the terms are, in sequence, the ideal gas contribution, the hard-sphere contribution,the mean-field contribution and the external field contribution. The ideal-gas part is F (id) [ ρ ] = k B T (cid:90) (cid:0) ρ ( r ) ln l ρ ( r ) − ρ ( r ) (cid:1) d r (2)where k B is Boltzmann’s constant, T is the temperature and l is any convenient lengthscale. The next two terms depend on the intermolecular potential which we take to be apair potential v ( r ). This is separated into a repulsive part, v ( r ) and an attractive part w att ( r ) = v ( r ) − v ( r ). We use the WCA prescription whereby v ( r ) = v ( r ) − v ( r ) for r < r and zero for r ≥ r where r is the minimum of the potential. Then, an effectivehard-sphere diameter d is constructed using any convenient prescription: here, we use theBarker-Henderson recipe d = (cid:90) r (cid:0) − e − βv ( r ) (cid:1) dr (3)where β ≡ /k B T . For the hard-sphere contribution, we use the FMT ansatz F (HS) [ ρ ] = k B T (cid:90) Φ ( η ( r ) , s ( r ) , v ( r )) d r (4)where the fundamental measures are n ( α ) ( r ) = (cid:90) w ( α ) ( r − r (cid:48) ) ρ ( r (cid:48) ) d r (cid:48) (5)with the hard-sphere radius R = d/ w ( η ) ( r ) = Θ ( R − r ), w ( s ) ( r ) = δ ( R − r ) and w ( v ) ( r ) = r r δ ( R − r ) where Θ ( x ) is the step function equal to 1 if x > η, s, v ) and here we use the modified RSLT functionΦ ( η, s, v ) = − πσ s ln(1 − η ) + 12 πσ s − v − η + 124 π s (1 − ( v /s )) (1 − η ) φ ( η ) (6)with φ ( η ) = 1 − − η + 3 η − − η ) ln(1 − η )3 η (7)4hich is easily shown to be free of instabilities[13]. The mean-field term is already givenexplicitly and we just note that w att ( r , r ) = w att ( | r − r | ), the attractive part of thepotential. Finally, φ ( r ) represents any external one-body field (this plays no role in thepresent work) and (cid:101) µ = µ − ln (cid:0) Λ l (cid:1) where µ is the chemical potential and Λ is the thermalwavelength. The model functional has been presented for a single component system butgeneralization to multiple components is straightforward. B. Implementation
We discretize the density on a cubic lattice of N x × N y × N z points with spacing ∆ andour key approximation is that the density field is approximated by trilinear interpolation ofthe values of the density at the lattice positions. So, working in units of ∆ = 1, and defining ρ S = ρ ( S x (cid:98) x + S y (cid:98) y + S z (cid:98) z ) where S x , S y , S z are integers with 1 ≤ S x ≤ N x , etc. we write ρ ( r ) = (cid:88) I =0 , A I ( r − S ( r )) ρ ( S ( r ) + I ) (8)where S ( r ) is the lattice point nearest the origin which is a corner of the computational cellcontaining the point r , i.e. that for which | S x | ≤ | r x | < | S x | + 1, etc and (cid:80) I =0 , is shorthandfor (cid:80) I x =0 , (cid:80) I y =0 , (cid:80) I z =0 , . The coefficients are A I ( r ) = ( δ I x (1 − x ) + x ) (cid:0) δ I y (1 − y ) + y (cid:1) ( δ I z (1 − z ) + z ) (9)which just means that we interpolate the density linearly between the lattice positionsforming the cell containing the point r .
1. Evaluation of the ideal part of the free energy
In principle, one would like to write the exact expression for the ideal gas contribution as F (id) [ ρ ] = (cid:88) S (cid:90) C ( S ) (cid:0) ρ ( r ) ln l ρ ( r ) − ρ ( r ) (cid:1) d r (10)where C ( S ) is the cell for which the corner at lattice site S is the closest to the origin.One could then, in each of the integrals insert the trilinear interpolation and perform theintegral. However, given the non-linearity of the expression, we have found this prohibitivelydifficult to do analytically, although it seems in principle possible. We have therefore used5he simpler approximation whereby the density is treated as being constant in each cellgiving the straightforward discretization F (id) [ ρ ] (cid:39) k B T ∆ (cid:88) S (cid:0) ρ S ln l ρ S − ρ S (cid:1) . (11)Our benchmarking suggests that this remains surprisingly accurate even for relatively highlylocalized density distributions (see Supplementary material). Note that the trilinear inter-polation, in any case, gives the exact result for the average number of particles (cid:104) N (cid:105) = (cid:90) ρ ( r ) d r = (cid:88) S (cid:90) C ( S ) ρ ( r ) d r = trilinear ∆ (cid:88) S ρ S . (12)
2. Evaluation of hard-sphere contribution
The FMT contribution to the free energy is discretized in the simplest way as F (HS) [ ρ ] (cid:39) k B T ∆ (cid:88) S Φ ( η S , s S , v S ) . (13)In principle, one could use a more accurate scheme but since the fundamental measures aremore slowly varying than the density itself, this seems sufficiently accurate on a computa-tional lattice that is fine enough to resolve the density. The main effort is therefore theevaluation of the fundamental measures at the lattice points.The standard method of computing the fundamental measures is to note that they areconvolutions so that they can be efficiently evaluated by transforming the weights and densityto Fourier space, multiplying and then performing and inverse Fourier transform. Theobvious strategy is to Fourier transform the weights analytically since they are simple todo while the density requires a discrete Fourier transform. This in principle leads to aninconsistency where it is impossible to prove that the resulting approximation preserves thestability of the free energy functional. We therefore follow the example of Ref.13 and doeverything consistently on the lattice. This leads to the practical problem that the FMTweights are not well adapted to evaluation on a lattice: the step function can be handledin obvious ways but the delta-functions are less obvious. In Ref.13 this was addressed byusing pre-compiled integration points on a sphere but this is a less than optimal solution asit leads to small-scale variations due to the pseudo-random nature of the points. Here, wepresent a straightforward alternative that avoids these issues by performing the necessaryintegrals analytically giving a fast and easily coded implementation.6he fundamental measures evaluated at the lattice sites are n ( α ) ( S ) = (cid:90) w ( α ) ( S − r ) ρ ( r ) d r (14)so that substituting the trilinear interpolation for the density and making a few simpletransformations leads to n ( α ) ( S ) = (cid:88) S (cid:48) (cid:101) w ( α ) ( S − S (cid:48) ) ρ ( S (cid:48) ) (15)with (cid:101) w ( α ) ( S ) = (cid:88) K =0 , (cid:90) C (0) w ( α ) ( S + K − r ) A K ( r ) d r (16)where C (0) is the cubic cell defined by the diagonal points (0 , ,
0) and (1 , ,
3. The mean-field contribution
Evaluation of the mean-field term using the trilinear interpolation scheme for the densityresults in F (MF) [ ρ ] = 12 ∆ (cid:88) S , S (cid:48) ρ S ρ S (cid:48) (cid:101) w (MF) S − S (cid:48) (17)where the weights have the form (cid:101) w (MF) S = ∆ − (cid:88) L =0 , (cid:88) K =0 , (cid:90) C ( ) (cid:90) C ( ) A L ( r ) A K ( r ) w att ( | r − r + K − L + S | ) d r d r (18)Because the attractive part of the potential typically varies slowly over the lengthscale ofa hard sphere, it turns out to be quite accurate to use the simple approximation (cid:101) w (MF) S = w att ( S ) (see Supplementary text for illustration). In principle, one could approximate thepotential by some sort of interpolation - e.g. polynomial - within each computational celland evaluate the resulting expressions for the weights analytically but we find no advantageto this complication. Once again, once the weights (cid:101) w S are known, the evaluation of theforces and energy is efficiently coded via discrete Fourier transforms.7 . Thermodynamics of the homogeneous fluid phases The thermodynamics of the bulk fluid phases of this model are simple to describe. Whenthe local density is constant, ρ ( r ) = ρ , the fundamental measures are as well and have thevalues η ( r ) = π ρd , s ( r ) = πρd and v ( r ) = 0. The free energy functional becomes asimple function 1 V β
Λ ( ρ ) = ρ ln ρl − ρ + ρ η (4 − η )(1 − η ) + 12 βa vdw ρ − (cid:101) µρ (19)where the Van der Waals constant is a vdw = 4 π (cid:90) ∞ r w ( r ) dr → ∆ (cid:88) S (cid:101) w S . (20)The equilibrium state must minimize this so that the equilibrium density ρ satisfiesln ρl + η (3 η − η + 8)(1 − η ) + βa vdw ρ = (cid:101) µ (21)and the grand-canonical free energy is1 V β
Ω = 1
V β
Λ ( ρ ) = − ρ (1 + η + η − η )(1 − η ) − βa vdw ρ = − βp (22)where the last equality reminds us that this is just the negative of the pressure. This canbe written as a virial expansion, βpρ = 1 + (cid:18) πd βa vdw (cid:19) ρ + 10 (cid:18) πd (cid:19) ρ + 18 (cid:18) πd (cid:19) ρ + ... (23)so that in particular, the second virial coefficient is B = πd + βa vdw .Two phase coexistence is possible when there are two solutions, ρ and ρ , which are bothglobal minima of the free energy function and so give equal free energies,ln ρ l + η (3 η − η + 8)(1 − η ) + βa vdw ρ = ln ρ l + η (3 η − η + 8)(1 − η ) + βa vdw ρ (24) − ρ (1 + η + η − η )(1 − η ) − βa vdw ρ = − ρ (1 + η + η − η )(1 − η ) − βa vdw ρ The spinodals are the densities at which the derivative of the pressure vanishes, dβpdρ (cid:12)(cid:12)(cid:12)(cid:12) s = 1 + 4 η s + 4 η s − η s + η s (1 − η s ) + β s a vdw ρ s = 0 (25)8hich, incidentally, just requires solving a quintic polynomial equation to determine thespinodal densities. The critical point occurs at the inflection point of the pressure and so isdetermined by the system dβpdρ (cid:12)(cid:12)(cid:12)(cid:12) c = 1 + 4 η c + 4 η c − η c + η c (1 − η c ) + β c a vdw ρ c = 0 (26) d βpdρ (cid:12)(cid:12)(cid:12)(cid:12) c = πd η c − η c )(1 − η c ) + β c a vdw = 0Eliminating the mean-field term between these leaves an equation involving only the densityand this has only a single physical solution that is easily determined numerically resultingin ρ c = 0 . d (27) k B T c = − . a vdw d p c ρ c k B T c = 0 . a vdw . III. CALCULATIONAL PROCEDURESA. Gaussian profiles
We have calculated results both by minimizing the discretized density field, which werefer to as “full” minimization. We have also, for comparison, performed calculations bymodeling the solid as a sum of Gaussians at the lattice sites. For the latter, the density fieldis ρ ( r ) = (1 − c ) (cid:16) απ (cid:17) / (cid:88) i exp (cid:0) − α ( r − R i ) (cid:1) (28)where the sum is over the four positions of the molecules in the unit cell ((0 , , a/ , a/ , , a/ , a/ a/ , , a/
2) where a is the lattice constant). The parameter α controls the width and c is the vacancy concentration (number of vacancies per latticesite). Although this model can be implemented more efficiently, we have simply used itto determine the discretized density from which the density functional is evaluated as in9he case of full minimization. The constrained, or Gaussian, minimization thus consists ofminimizing with respect to the two parameters c and α . We carry out this minimizationusing the Nelder-Mead algorithm[16] as provided in the GSL library[17]. Note that in anyminimization of the FMT functional, it is possible that density fields are generated whichcause the local packing fraction, η ( r ), to exceed one, which is outside the domain of thefunction. This simply indicates that the minimization routine has made too large of anadjustment to the density and in this case we return a large value for the free energy thuspushing the search back into the physical region. B. Full minimization
For full minimization of the discretized density field we start with an initial guess (basede.g. on minimized Gaussian profiles or the density field minimized at some other thermody-namic parameters) and use the Fast Inertial Relaxation Engine (FIRE) algorithm[18]. Thisis a type of gradient descent with inertia and we find it reliably converges in typically a fewthousand iterations for the case of the solid (and an order of magnitude faster for the fluid).Details concerning the parameters used are given in the Supplementary text and here weonly note that as in the case of the Gaussian profiles, care must be taken to backtrack if anattempted adjustment of the density takes it outside the physical domain.
C. Determining properties of the homogeneous solid phase
A difficulty of the solid phase is that one must take explicit account of the periodicity ofthe lattice, even when working directly with the density field. This is because a homogeneoussystem is necessarily modeled using a finite computational cell with periodic boundaries sothat its size must be commensurate with the lattice. An FCC solid can be described usinga non-primitive cubic simulation cell with lattice sites at each of the 8 corners as well as onthe centers of each face. If the computational lattice spacing is ∆ and if the length of the cellis N computational nodes, then its physical length is a = N ∆ which will then be the latticespacing of the solid phase. Our procedure is to fix these quantities and then minimize thedensity profile at constant chemical potentials µ as explained in the previous subsection.The result is the free energy functional evaluated at these parameters, Ω( N, µ ; T, ∆) and the10orresponding density field. We then change the chemical potential to µ + ∆ µ and use thepreviously determined density field as the initial guess for full minimization at this chemicalpotential. This is repeated over a range of chemical potentials so that what we end up withis Ω( N, µ + m ∆ µ ; T, ∆) for different values of N and m . Then, for each value of m , i.e. foreach chemical potential in the grid, we locate the three values of N which contain a minimumof the free energy and finally, estimate the the optimal (non-integer) value of N by quadraticinterpolation on these three values. This results in a list of free energies, Ω( µ + m ∆ µ ; T, ∆),and these are used to find coexistence with the vapor and/or liquid phases (i.e. the valuesof chemical potentials where the two phases have equal free energies). This is again refinedusing quadratic interpolation. All quantities are then determined by the same quadraticinterpolation except the vacancy concentration which is sometimes determined via linearinterpolation (because the quadratic interpolation fails for this quantity in some cases).We note that an alternative would be to hold N fixed and to vary the lattice spacing ∆.We choose not to to do this for several reasons. First, the value of the VdW parameter, a VdW that determines the bulk thermodynamics of the liquid phase changes as we change∆ but not when we change N so that this introduces some unphysical variation into themodel. Furthermore, in this case, there are two ways to lower the density of the homogeneousphase: by keeping the local density fixed and changing ∆ or by holding the spacing fixedand varying the local density. This makes the limit of the homogeneous liquid ambiguous.One might also wonder about the relevance of reporting a minimum in the lattice spacingof fractional values of N determined by interpolation. In fact, the limitation to the discretevalues of lattice parameter is an artifact of trying to limit the cost of the calculations byusing a minimal cubic cell. If, e.g., the computational cell were two lattice spacings in length,so that 2 a = N ∆, then one could place one FCC lattice position at the origin, (0 , , N/ , ,
0) and one at ( N ∆ , ,
0) so that the FCC lattice spacing would be ( N/ D. Procedure to compute solid properties with the HCP lattice
In addition to the solid phase computations using the FCC lattice, we performed cal-culations for the HCP (hexagonal close-packed) lattice, using the Gaussian profiles. The11eometry of the HCP structure makes it more difficult to study than the FCC one. This isbecause our implementation of classical DFT computations uses rectangular cells with thesame grid spacing ∆ in all directions. It is possible to use a rectangular cell to constructthe HCP lattice that is compatible with periodic boundary conditions (see Fig. 1) but theratio between the side’s lengths are irrational numbers. That means we cannot constructthe HCP lattice directly as our implementation requires all rectangular side’s length to bemultiples of the same grid spacing ∆. Therefore, instead of directly constructing the HCPlattice, we compute the solid properties for regular hexagonal lattices near close-packing andinterpolate the results for the ideal close-packing (HCP) proportions. More precisely, we firstcompute the rectangle side’s lengths that are multiples of the grid spacing ∆ and that arethe closest to the ideal close-packing proportions. They define our reference cell. Then, wegenerate other rectangular cells by adding -1, 0 or 1 times the grid spacing to the side lengthof the reference cell, for two of the three axis x, y, z. We perform computations for all ofthese cells and then interpolate the results successively along each of the two selected axis,using quadratic interpolations.
FIG. 1. The HCP structure. The blue lines highlight the borders of the rectangular cell used forcomputations. The blue dots are the lattice sites which constitute the basis of the rectangular cell. V. RESULTS
We have performed calculations for three potentials: the Lennard-Jones potential andtwo potentials recently proposed by Wang et al. The Lennard-Jones potential is v LJ ( r ) = 4 (cid:15) (cid:18)(cid:16) σr (cid:17) − (cid:16) σr (cid:17) (cid:19) (29)which is then cutoff at a distance r c and shifted to give v LJ ( r ; r c ) = ( v LJ ( r ) − v LJ ( r c )) Θ ( r c − r ) . (30)The LJ potential is widely used to model simple (e.g. atomic or small-molecule) fluids andmetals. The Wang-Ramirez-Hinestrosa, Dobnikar and Frenkel (WHDF) potentials v W HDF ( r ) = α ( r c ) (cid:15) (cid:18)(cid:16) σr (cid:17) − (cid:19) (cid:18)(cid:16) r c r (cid:17) − (cid:19) Θ ( r c − r ) (31)where α ( r c ) = 2 (cid:16) r c σ (cid:17) (cid:16)(cid:0) r c σ (cid:1) − (cid:17) (32)This potential was introduced as a simplification relative to the LJ potential which is de-signed to go smoothly to zero at its cutoff and to be deformable between a LJ-like potential(when the cutoff is large, e.g. r c = 2 σ ) and a colloid-like potential when the cutoff is small(e.g. r c = 1 . σ ). The difference in the two cases is attributable to the difference in thewidth of the attractive well compared to the repulsive part of the potential and results inqualitatively different phase diagrams as illustrated below. A. Surface tension
We have calculated the liquid-vapor surface tension for our model systems using a cellconsisting of a lattice of 1 × × N points with N = 20 ,
000 and a spacing of ∆ = 0 . ρ vap and ρ liq from theequation of state and the corresponding chemical potential. We then create an initial densityin which the two phases each occupy half the computational cell with a sharp boundarybetween them. We minimize to get the equilibrium density distribution ρ ∗ ( r ) and computethe surface tension as the excess surface free energy γ = (Ω[ ρ ∗ ] − Ω( ρ vap )) / (2∆ ). The results13re shown in Fig. 2 along with simulation results for the LJ system with the same cutoff[19]and the results of Wang et al[15]. The DFT model, which has no adjustable parameters,compares reasonably well with simulation, particularly for intermediate temperatures. Beinga mean-field model, it is not expected to capture the behavior near the critical point whilethe fact that the mean-field contribution is motivated by the high-temperature limit in liquidstate theory[20] may account for the systematic deviation at lower temperatures. FIG. 2. Surface tension for the Lennard-Jones potential with cutoff R c = 3 σ and 4 σ and theWHDF potential for a simple fluid, cutoff of R c = 2 σ , and for a colloidal system, R c = 1 . σ .The calculations were performed using the DFT model with no adjustable parameters. Simulationresults for the LJ system[19] and the WHDF systems[15] are also shown with error bars of onestandard error. B. Phase diagrams
The vapor-liquid phase diagrams for our systems were already calculated in the course ofevaluating the surface tension. We have also calculated the FCC-solid phase diagrams usinga cubic cell with a fixed lattice spacing, ∆. The cell represents a non-primitive unit cell ofthe solid so that its length is the lattice spacing of the FCC solid and in a perfect solid, the14otal number of molecules in such a cell would be 4. However, when minimizing the freeenergy functional, the total number of molecules varies and in general is less than this inthe equilibrium state: the difference is a measure of the equilibrium vacancy concentrationin the solid phase which we calculate as c = − N [ ρ ∗ ]4 where N [ ρ ] = (cid:82) C ρ ( r ) d r is the totalnumber of molecules in the cell C .The resulting phase diagram for a Lennard-Jones potential is shown in Fig.3. The mean-field model reproduces the correct qualitative behavior with vapor, liquid and solid binodalsand a liquid-vapor spinodal. The Gaussian and full mininizations are in close agreementthus showing the close correspondence of the two. As seen in the Figure, the vapor-liquidcritical point becomes lower as the cutoff of the potential is reduced but the quualitativebehavior does not change. In contrast, the liquid-solid binodals show little sensitivity to thecutoff.Figure 4 shows the computed phase diagram for the WHDF potential with a cutoff of 2 σ which is quite similar to the Lennard-Jones phase diagram. Comparison to the simulationdata of Wang et al[15] shows that, while qualitatively quite realistic, the mean-field modeldoes not agree quantitatively with simulation. In particular, the fluid-solid binodals aredisplaced towards lower densities and higher temperatures than in the simulations. Thisis to be expected since no attempt is made in the mean-field model to reproduce even thefluid-phase thermodynamics quantitatively.The results for the WHDF potential with a smaller cutoff, producing a colloid-like phasediagram, are shown in Fig.5 where the typical suppression of the liquid-vapor transitioninto the metastable region of the fluid-solid transition is evident. The cDFT again faithfullyreproduces this qualitative behavior and, indeed, is in reasonable quantitative agreementwith the simulations. In summary, while the cDFT is not reliably quantitatively accurate,it does a good job of tracking the qualitative behavior resulting from variations in theinteraction potential. C. Vacancies
The vacancy concentration determined by the DFT calculations is shown as a functionof temperature in Fig. 6 for full minimization and in Fig. 7 in the Gaussian approximation.The results are qualitatively the same in the two cases and the numerical differences, while15
IG. 3. Vapor-liquid-FCC phase diagram for the LJ potential for cutoffs of 3 σ and 4 σ as obtainedby full minimization of the DFT functional. The broken lines show the spinodals. The resultsobtained using Gaussian profiles are also shown and are very close to the full minimization. real, are modest overall. In Fig. 8 we compare both calculations to some simulation data forthe LJ potential, but results are only available for higher temperatures. For the LJ system,the vacancy concentration is relatively insensitive to the cutoff and and is nearly constantfor liquid-solid coexistence, at least in the range of temperatures reported here. This isqualitatively consistent with the simulation results and is much better than older calculationsthat used more primitive models of the hard-sphere free energy functional[21, 22]. At lowertemperatures, on the solid-vapor coexistence curve, it drops sharply with temperature as onewould expect. The WHDF potential for simple fluids gives similar results near the triplepoint but drops as the temperature increases and actually becomes negative (indicatinginterstitials rather than vacancies). This may well be an unphysical artifact of the cDFTcalculations - or due to inaccuracies in the interpolations - but there are no simulationresults to compare to. We note that the Gaussian minimizations do not show this behavior.16 IG. 4. Vapor-liquid-FCC phase diagram for the WHDF potential with range R c = 2 σ whichbehaves like the Lennard-Jones potential as computed via full minimization of the DFT functionaland using Gaussian profiles. The coexistence curves based on the data of Wang et al.[15] are shownas shaded regions. Finally, the results for the WHDF colloidal potential seem reasonable but there are againno independent results for comparison.
D. Density profiles
It has been seen that the use of Gaussian profiles produces results which are quanti-tatively very similar to the results of unconstrained minimization of the cDFT functional.Nevertheless, differences are expected due to the fact that the neighborhood of a latticeposition in the solid is not spherically symmetric so that the Gaussian profiles, which forcespherical symmetry of the contribution of each lattice site, can only be an approximation.In Fig. 9 we show the density profile for a LJ solid near the triple point ( k B T = 0 . IG. 5. Vapor-liquid-FCC phase diagram for the WHDF potential with range R c = 1 . σ whichmodels a colloidal phase diagram as computed via full minimization of the DFT functional andusing Gaussian profiles. The coexistence curves based on the data of Wang et al.[15] are shown asshaded regions. a = 66 × . σ = 1 . σ ). The central region of the profile is well-fit by a Gaussian butaway from the center, the distribution is much broader than a Gaussian. The differenceis not due to contributions from the neighbors: the nearest neighbor distance is a/ √ e − r and a best fit normalized Gaussian has α = 86. Eventaking the slightly broader latter function, the contribution at the half the nearest neighbordistance is (cid:0) π (cid:1) . e − . × a/ √ ≈ × − which is much smaller than the excess in the tailof the density. Figure 10 shows the difference in density along lines running from a latticeposition in the direction (110) and (111) relative to that along the (100) direction and onesees that at intermediate distances, there is some asymmetry. This may be attributable tothe fact that moving along the (110) direction means towards a nearest neighbor and atthis density, the nearest neighbor distance, 1 . σ is slightly more than the position of the18 IG. 6. The vacancy concentration for the colloidal potential WHDF with cutoff 1 . σ on the solid-fluid coexistence curves and the simple fluid potentials on the liquid-solid (higher temperature)and vapor-solid (lower temperature) coexistence curves, obtained with the full minimization. Inthe case of the WHDF potential with cutoff 2 σ , the open symbols are values for which c is negativeand so the absolute value is shown on the logarithmic scale used here. minimum of the potential well, 1 . σ , giving a small preference for adding density in thatdirection at the expense of the others. E. Relative stability of FCC and HCP configurations
The free energy difference between the FCC and HCP configurations is reported on Fig.11as a function of the temperature. These DFT computations use the Gaussian profiles to pa-rameterize the density field and have been performed for the WHDF potential with r c = 1 . r c = 2 . IG. 7. The same as Fig. 6 but showing the results in the Gaussian approximationFIG. 8. Comparison of vacancy concentrations on the solid-liquid coexistence curve for theLennard-Jones potential with cutoff 4 σ . This graph shows our results (line and open circles)next to simulations[23] (line with error bars) and the older calculations of Singh et al[22] andMcrae[21]. IG. 9. Density distribution near a lattice site for a LJ solid with cutoff 3 σ , k B T = 0 . N = 66lattice points in each direction and βµ = − α ≈
95 and 86, respectively. Significant variation from the Gaussians in the tailregion cannot be accounted for by nearest-neighbor contributions. for the two configurations, about the order of magnitude of the expected interpolation errorwhich is estimated to be ∆ βω ∼ − for the larger cutoff and an order of magnitude largerfor the smaller cutoff. Such small differences are expected because these two structures onlydiffer at the second neighbor. Wang et al. also state that simulations for the WHDF po-tential lead to very similar free energies for the FCC and HCP structures [15]. Our resultsindicate that the presented DFT model correctly reproduce this behavior and that using theFCC lattice gives a reasonable description of the solid free energies whether or not it is themost stable configuration. 21 IG. 10. Illustration of asymmetries in the density distribution. The plot shows the differencein density near a lattice site moving along the (110) (black, full line and squares) and (111) (red,broken line and circles) directions compared to that along the (100) direction. The conditions arethe same as for Fig. 7.
V. CONCLUSIONS
We have presented a fully robust finite-elements implementation of the standard cDFTmodel. The novel elements of our implementation are (a) it is entirely formulated in real-space and only uses discrete FFT’s to evaluate discrete convolutions efficiently rather thanmixing discrete real-space quantities and analytically Fourier-transformed continuous quanti-ties in an uncontrolled manner and (b) the real-space weights needed in FMT are evaluatedanalytically thus avoiding the need for using tables of pre-compiled spherical integrationpoints as were previously used[13]. The latter fact significantly reduces numerical noisein the calculations and, together with the intrinsic stability of the methods, we routinelydo minimizations of the solid phase at constant chemical potential which, as noted in the22
IG. 11. Free energy difference per unit volume, βω = βω HCP − βω FCC as a function of temperaturecomputed using Gaussian profiles and the WHDF potential. For each temperature, the computa-tions of both FCC and HCP phases were performed at the same chemical potential, close to thecorresponding fluid-solid coexistence. We estimate that the numerical errors in the differences areon the order of 1 × − for the larger cutoff and 1 × − for the smaller one. Introduction, has previously proven infeasible.We have used this algorithm to first reproduce standard results for the Lennard-Jonespotential and then to test the robustness of cDFT in describing more general phase diagramsas result from the WHDF potentials. We confirm that the model performs well across alltests. We believe that this implementation is well suited for more challenging applicationssuch as dynamic DFT and studies of nucleation. Finally, we have verified the accuracy ofthe computationally cheap Gaussian approximation which opens the possibility to the useof pseudo-spectral methods with radial basis functions that could provide a more efficient23lternative to the finite-difference algorithms commonly used. [1] R. Evans, “The nature of the liquid-vapour interface and other topics in the statistical me-chanics of non-uniform, classical fluids,” Adv. Phys. , 143 (1979).[2] James F. Lutsko, “Recent developments in classical density functional theory,” Adv. Chem.Phys. , 1 (2010).[3] Robert Evans, Maria C. Stewart, and Nigel B. Wilding, “A unified description of hy-drophilic and superhydrophobic surfaces in terms of the wetting and drying transitionsof liquids,” Proceedings of the National Academy of Sciences , 064110 (2020),https://doi.org/10.1063/1.5142651.[5] James F. Lutsko, “How crystals form: A theory of nucleation pathways,” Sci. Adv. , eaav7399(2019).[6] Y. Rosenfeld, “Free-energy model for the inhomogeneous hard-sphere fluid mixture anddensity-functional theory of freezing,” Phys. Rev. Lett. , 980 (1989).[7] J. K. Percus, “Equilibrium state of a classical fluid of hard rods in an external field,” J. Stat.Phys. , 505–511 (1976).[8] J. K. Percus, “One-dimensional classical fluid with nearest-neighbor interaction in arbitraryexternal field,” J. Stat. Phys. , 67 (1981).[9] T. K. Vanderlick, H. T. Davis, and J. K. Percus, “The statistical mechanics of inhomogeneoushard rod mixtures,” J. Chem. Phys. , 7136 (1989).[10] Roland Roth, “Fundamental measure theory for hard-sphere mixtures: a review,” Journal ofPhysics: Condensed Matter , 063102 (2010).[11] M. Oettel, S. G¨orig, A. H¨artel, H. L¨owen, M. Radu, and T. Schilling, “Free energies, vacancyconcentrations, and density distribution anisotropies in hard-sphere crystals: A combineddensity functional and simulation study,” Phys. Rev. E , 051404 (2010).[12] Mostafa Mortazavifar and Martin Oettel, “A fundamental measure density functional for fluid nd crystal phases of the asakuraoosawa model,” J. Phys.: Condens. Matter , 44018 (2016).[13] James F. Lutsko and Julien Lam, “Classical density functional theory, unconstrained crystal-lization, and polymorphic behavior,” Phys. Rev. E , 012604 (2018).[14] James F. Lutsko, “Explicitly stable fundamental measure theory models for classical densityfunctional theory.” Http://arxiv.org/abs/2009.09390.[15] Xipeng Wang, Simn Ramrez-Hinestrosa, Jure Dobnikar, and Daan Frenkel, “The lennard-jones potential: when (not) to use it,” Phys. Chem. Chem. Phys. , 10624–10633 (2020).[16] J. A. Nelder and R. Mead, “A simplex method for function minimization,” The ComputerJournal , 308–313 (1965).[17] M. Galassi et al, “Gnu scientific library,” , accessed:2018-04-01.[18] Julien Gunol, Wolfram G. Nhring, Aviral Vaid, Frdric Houll, Zhuocheng Xie, Aruna Prakash,and Erik Bitzek, “Assessment and optimization of the fast inertial relaxation engine (fire) forenergy minimization in atomistic simulations and its implementation in lammps,” Computa-tional Materials Science , 109584 (2020).[19] Patrick Grosfils and James F. Lutsko, “Dependence of the liquid-vapor surface tension on therange of interaction: A test of the law of corresponding states,” J. Chem. Phys. , 054703(2009).[20] J.-P. Hansen and I.R. McDonald, Theory of Simple Liquids (Academic Press, San Diego, Ca,1986).[21] Robin McRae, John D. McCoy, and A. D. J. Haymet, “Density functional theory of vacancies,”Journal Chemical Physics , 4281 (1990).[22] Sunil P. Singh and Shankar P. Das, “Perturbation theory for classical solids with vacancydefects,” Phys. Rev. B , 144113 (2007).[23] Apoorva Purohit, Andrew J. Schultz, Sabry G. Moustafa, Jeffrey R. Errington, and David A.Kofke, “Free energy and concentration of crystalline vacancies by molecular simulation,”Molecular Physics , 3027–3041 (2018). ppendix A: The FMT weights In this Appendix, we give a somewhat simplified form for the analytic determination ofthe discrete FMT weights which is discussed in detail in the Supplementay material. Forthe local packing fraction, one needs to evaluate w ( η ) ( S ) = (cid:88) I =0 , (cid:90) C ( ) Θ ( R − | S + I | ) A I ( r ) d r . (A1)Before tackling this, we note that the other fundamental measures follow directly from thisone via w ( s ) ( S ) = ∂∂R w ( η ) ( S ) (A2) w ( v ) ( S ) = − ∂∂ S w ( η ) ( S )so we need only concentrate on determining the first one. The details are given in the Sup-plementary material and here we just report the result which has been somewhat simplifed.In general, for all of the weights one finds that (cid:101) w ( α ) ( S ) = (cid:88) I ∈{− , } (cid:88) K = , I Θ (cid:0) R − ( S x + K x ) − ( S x + K y ) − ( S x + K z ) (cid:1) ( − K x + K y + K z (A3) × G ( α ) ( S + I , S + K ) . For the local packing fraction, and for the positive octant S j ≥
0, the required function is G ( η ) ( T , V ) = − (cid:0) R − V (cid:1) − (cid:0) R − V (cid:1) ( T · V ) (A4) − T x T y T z V x V y V z + π T x T y T z R + P xyz J ( T , V )with J ( T , V ) = − T x T y T z R arcsin V x V y (cid:113) ( R − V x ) (cid:0) R − V y (cid:1) (A5)+ 124 T y T z (cid:16) (cid:0) R − V x (cid:1) + 4 T x V x (cid:0) R − V x (cid:1)(cid:17) (cid:32) arcsin V y (cid:112) R − V x − π (cid:33) + 1120 T z (cid:16) (cid:0) R − V x − V y (cid:1) + 5 T x V x (cid:0) R − V x − V y (cid:1) + 20 T x T y V x V y (cid:17) (cid:113) R − V x − V y − (cid:0) T x V x + 3 T x T y V x V y (cid:1) (cid:0) R − V (cid:1) − T x V x − T x T y V x V y P xyz is an operator indicating a sum over all 6 permutations of the x,y and zcomponents of the vectors (performed simultaneously: that is, one element of the sum is J ( T x , T y , T z , V x , V y , V z ), another is J ( T y , T x , T z , V y , V x , V z ), etc.). For the packing fraction,the spherical symmetry implies that the weights for other octants (e.g. with S x <
0) are thesame as for | S x | . The weights for the other fundamental measures follow from this result viadifferentiation (they are given explicitly in the Supplementary text). Note that it is easy tosee that J ( T , V ) = 0 if V = R so terms involving a derivative of the step function gives nocontribution.These results have been checked by (a) independently deriving the result for w ( s ))