Classical density functional theory, unconstrained crystallization, and polymorphic behavior

While in principle, finite temperature density functional theory (ftDFT) should be a powerful tool for the study of crystallization, in practice this has not so far been the case. Progress has been hampered by technical problems which have plagued the study of the crystalline systems using the most sophisticated Fundamental Measure Theory models. In this paper, the reasons for the difficulties are examined and it is proposed that the tensor functionals currently favored are in fact numerically unstable. By reverting to an older, more heuristic model it is shown that all of the technical difficulties are eliminated. Application to a Lennard-Jones fluid results in a demonstration of power of ftDFT to describe crystallization in a highly inhomogeneous system. First, we show that droplets attached to a slightly hydrophobic wall crystallize spontaneously upon being quenched. The resulting crystallites are clearly faceted structures and are predominantly HCP structures. In contrast, droplets in a fully periodic calculational cell remain stable to lower temperatures and eventually show the same spontaneous localization of the density into 'atoms' but in an amorphous structure having many of the structural charactersitics of a glass. A small change of the protocol leads, at the same temperature, to the formation of crystals, this time with the FCC structure typical of bulk Lennard-Jones solids. The FCC crystals have lower free energy than the amorphous structures which in turn are more stable than the liquid droplets. It is demonstrated that as the temperature is raised, the free energy differences between the structures decreases until the solid clusters become less stable than the liquid droplets and spontaneously melt. The presence of energy barriers separating the various structures is therefore clearly demonstrated.


I. INTRODUCTION
Work on classical Density Functional Theory (cDFT) for classical systems [1,2] has long been driven by the promise of having a tool that allows for the determination of free energies and inhomogeneous density distributions for arbitary systems such as confined liquids and solids and glasses (confined or not) and the energy barriers separating these states that would then open the way to studying all aspects of nucleation: lifetimes of metastable states, transition pathways etc without problem-specific modelling bias. Indeed, these goals are so important that they have driven much recent work in other fields as well. A good example is the field of computer simulation where the search for unbiased methods of studying free energy surfaces and transitions between metastable states that do not rely on the definition of collective variables to characterize such states has been prominant (see, e.g., the contributions of Piaggi et al [3] and of Swinburn and Marinica [4] that present two recent steps in this direction). In principle, cDFT evades the need for collective variables since one works directly with the density field. A well-known and closely related alternative to cDFT are Phase-Field theories which are a more mesoscopic theoretical approach that extends the Landau theory of phase transitions by adding additional variables and terms to the gradient-expanded free energy as suggested from DFT (for an overview, see e.g. [5] where it is emphasized that phase-field theory can perhaps best be viewed as a mesoscopic, coarse-grained version of cDFT). But the truely microscopic approach based on cDFT that is capable of accurately describing atomic-scale correlations in liquids as well as solids and amorphous systems has so far remained elusive.
All forms density functional theory (classical and quantum, zero and finite temperature) are formally based on exact theorems stating that for a fixed temperature and chemical potential in the grand canonical ensemble, there exists a functional Λ [n] of the local number density n (r) (or partial densities for a multicomponent system) having the property that its value is uniquely minimized by the equilibrium density, i.e. if the equilibrium density is n eq (r) then Λ [n] ≥ Λ [n eq ] with equality if and only if n (r) = n eq (r) . Furthermore, when evaluated at this minimizing density field, one has that Λ [n eq ] = Ω, the grand canonical free energy for the system. The functional Λ depends also on any external fields present, but this dependence is trivial, being given by Λ [n; φ] = F [n] + (φ (r) − µ) n (r) dr where φ (r) is the external potential, µ is the chemical potential and F [n], called the Helmholtz free energy functional, is the part of Λ independent of the field. In general, F [n] is not known and so applications of cDFT depend on some sort of approximation for it. These are often constructed based on the rare systems for which exact results are possible such as the ideal gas, hard particles in highly confined geometries and hard-rods in one dimension. For this reason, the most highly refined theories in three dimensions are for hard-spheres and so it has long been the case that a crucial test of model cDFT's has been their ability to describe inhomogeneous hard-sphere systems, including the solid phase and the liquid-solid transition, see e.g. [1,2].
The most accurate and widely accepted models today are of a class collectively known as Fundamental Measure Theory (FMT) [6]. These arose out of the exact results of Percus and coauthors for one-dimensional systems [7][8][9] as well as an approximate description of hard-sphere statistical mechanics called Scaled Particle Theory [10]. Rosenfeld's original FMT gave a good description of fluids near a wall but failed to stabilize (i.e. to have a minimum corresponding to) the hard-sphere solid phase. After much work by several authors, Tarazona suggested a modification of Rosenfeld's theory that succeeded in stabilizing the hard-sphere solid and that was motivated by demanding that the functional reproduce certain exact results [11]. A heuristic modification of Tarazona's theory was proposed by Roth et al [12] that aimed to improve the description of the dense liquid and this version of FMT, known as the White Bear (WBI) functional, was for some time considered to be the best available model. For example, it gives a very good quantitative prediction of the hard-sphere freezing transition. A relatively recent modification, the so-called White Bear II (WBII) model, gives an even more accurate representation of the hard-sphere solid including a good description of the vacancy concentration [13]. These models are not perfect -there are difficulties in extending them to mixtures [14] and to metastable lattice structures [15] -but their value as the best current description of the hard-sphere solid cannot be questioned.
Originally, because of computational constraints, cDFT calculations for the solid phase were performed using a resticted model of the density field whereby it was represented as a sum of Gaussians centered at the fixed lattice sites of a given Bravais lattice. Only the widths of the Gaussians and sometimes their amplitudes, were varied so as to minimize Λ. This clearly accords with an intuitive idea of what a solid should look like and also has recovers the homogeneous liquid in the limit that the Gaussian widths become infinetly large. Non-Gaussian corrections were also sometimes investigated but generally seemed of minor importance. Only in the last 20 years have full-scale, uncontrained cDFT calculations been performed using either finite element or pseudo-spectral methods. In particular, Oettel et al have pioneered the investigation of the hard-sphere solid phase with finite element methods and using the WB models [16]. However, one surprising point noted in the course of these applications has been a numerical delicacy of the WB models. In order to obtain useful results and to avoid numerical divergences, the authors report having to work at constant number rather than constant chemical potential, as is more natural in the cDFT framework and in some cases being forced to carefully maintain strict cubic symmetry in the course of the numerical optimizations as well as having to use relatively fine numerical discretizations.
Because of these limitations, the results reported so far on the use of cDFT to describe the solid phase has generally been limited to homogeneous solids, to planar liquid-solid interfaces and to two-dimensional systems. For example, the homogeneous hard-sphere crystal has long been used as a test of cDFT models and, less often, the phase diagram of homogeneous solids for systems with pair potentials have been studied(see, e.g., the discussion in [17]). But in all cases, these calculations are based on pre-defined lattice structures and the creation "by hand" of localized density fields. One of the few applications beyond homogeneous crystals that have been studied in the literature by many authors is the planar liquid-solid hard-sphere interface, e.g. by Curtin [18], Lutsko [19] and Oettel and coworkers [20,21], but always starting with the introduction of pre-selected crystalline structures by hand. Similarly, cDFT has long been applied to study the free energy of glass-like systems [22,23] but, once again, based on pre-determined structures (typically the pair-distribution function of Bennett [24]). In contrast, in our work, we demonstrate the spontaneous formation of solid clusters having no relation at all to the symmetries of the calculational cell or of the applied boundary conditions and we derive a pair-distribution function from our stuctures that is the equivalent of the Bennett result. Only a few earlier works have demonstrated any spontaneous localization of a density field (e.g. the work on hard-sphere glasses of Dasgupta and Walls [25] and the extension of this work by Chaudhuri et al [26] which used very simple models limited to homogeneous systems, and that of Archer and coworkers [27,28] on two-dimensional freezing of soft-matter) and we are not aware of any such work applied to highly inhomogeneous systems such as clusters (as is necessary for applications like clustering in proteins, nucleation and crystallization) in three dimensions.
In this paper, we propose that the numerical difficulties of the WB functionals are due to a pathology of the tensor models: namely, that their approximation for Λ [n] is not bounded from below and that in fact it has no minimum. In this interpretation, the constraints mentioned above -cubic symmetry and constant number -suppress this pathology to allow results to be obtained from (artificial or actual) metastable stationary points of the functional. But while the constant number constraint is benign, the constraint of cubic symmetry clearly limits applicability of the model to high symmetry systems and thus restricts the overall utility of the model. In this paper, it is shown that using a more heuristic, but demonstrably stable model, eliminates all of these problems and leads to a much more robust description of the solid. The resulting model is not a definitive replacement for the tensor models -the exact limits satisfied by the latter are undoubtedly crucial to incorporate -but serves as a proof of concept of the feasibility of a fully robust cDFT solid.
We go on to use this robust hard-sphere functional to construct a functional for arbitrary pair potentials by treating the attractive tail of the potential using the standard mean-field approach. This model is used to perform calculations for a Lennard-Jones potential. We first determine the entire vapor-liquid-solid phase diagram. Next, we study the behaviour of small droplets in contact with a wall. As the temperature is quenched, typical layering structure develops within the droplets until at a certain point, the density field spontaneously localizes into atoms arranged in a faceted, HCP crystallite. The same protocol applied to the case of spherical droplets with no walls present shows a similar spontaneous localization but into an amorphous structure with many of the structural characteristics found in a glass. An intermediate protocol where a "wall" is present at the start of minimization but then removed results in the formation of homogeneous crystals which are again faceted but which have the FCC structure typical of bulk Lennard-Jones. We show that at intermediate temperatures, all three structures -droplet, crystal and amorphous cluster -are metastable minima of the free energy functional thus implying the presence of energy barriers separating them. These results clearly show the utility and far-reaching protential of this model for the study of nucleation, crystallization, polymorphism and, potentially, the glass transition.
Section II of this paper describes the technical problems with current hard-sphere functionals and our proposal for their origin and resolution. The result is a functional that is more numerically robust than any of the previous models. We implemented this functional in a new finite-element algorithm (described in detail in an appendix) designed to rigorouly maintain the term-by-term positivity of the contributions to the Λ-functional and validate it by detailed comparison to previous works on hard-spheres. In Section III, the new hard-sphere functional is used to construct a mean-field description of a system interacting via arbitrary pair-potentials and applied specifically to a Lennard-Jones potential. Calculations are presented for the phase diagram of the homogeneous system and of critical clusters for the vapor-liquid and vapor-solid transitions. We find that the solid clusters exhibit polymorphism depending on the environment in which they form. Section IV summarizes the results and discusses their implications for futher improvements of FMT hard-sphere models and for applications to the study of clustering, nucleation, crystallization and amorphous systems.

A. FMT Models
In this and the following Sections, attention will be restricted to single component systems in three dimensions composed of hard spheres with diameter σ. The Helmholtz free energy functional is written as with the ideal gas contribution given by where k B is Boltzmann's constant and T is the temperature and, for later use, we define β = 1/(k B T ). The second term is the excess free energy and in Fundamental Measure Theory is modeled using Rosenfeld's ansatz as where n α (r; [n]) stands for a family of linear functionals of the density having the generic form and the weights w α (r) are in general short-ranged functions. The standard density functionals are the local packing fraction, which is just the density averaged over a spherical volume corresponding to the hard-sphere diameter, and the scalar, vector and tensor surface averages or, equivalently, It is sometimes interesting to separate the latter into a sum of its trace and traceless part as T (r) = 1 3 s (r) 1 + U (r). Notice that since the density itself is always greater than or equal to zero, η, s and T are also positive semi-definite. In a liquid, for which the density is by definition uniform n (r) = n, these become Different FMT models are distinguished by the form of Φ (n α ). In Rosenfeld's original proposal, this was As long as η < 1, the first term on the right is obviously positive semi-definite. The second term shares this property as can be seen from the simple inequality so s 2 (r) ≥ v 2 (r). On the other hand, the third term is not. As a somewhat realistic example that illustrates this, consider a hard wall with normal in the z-direction located at z = 0. Let the density be zero for z < 0, "outside" the volume, and a constant ρ 0 for z > 0. The weighted densities are zero for z < −σ/2, and in the domain −σ/2 < z < σ/2 they are (see Appendix A for details) Clearly, for values of z in the given domain, the first expression is indeed non-negative. However, the second is not: for example, its value at z = −σ/4 is There is no a priori physical reason that this numerator must be non-negative and it is, in fact, the reason that Rosenfeld's model fails to stabilize the solid phase: in that application, the numerator can be negative while the local value of η tends to one near the lattice sites in the solid thus making the free energy unbounded from below. Nevertheless, Rosenfeld's model does give an impressively good description of a hard-sphere fluid near a wall and has the important property that it reproduces the Percus-Yevik direct correlation function for hard-spheres via the exact relation [2] c (r 12 ; n) = − lim n(r)→n and it follows that the functional also reproduces the Percus-Yevik equation of state for the uniform fluid. Therefore, much work went into trying to improve it. The best functionals currently in use are based on Tarazona's generalization, in which the numerator of the third term is replaced by The resulting functional, Φ T (n α ), is indeed able to stablize the solid phase and like Rosenfeld's original functional, Tarazona's reproduces the Percus-Yevik direct-correlation function and equation of state for the uniform fluid. It does not, however, give very good results for liquid-solid coexistence of hard spheres. This is because the Percus-Yevik approximation is not very good at liquid densities and so the liquid equation of state is poorly approximated. Roth et al [12] suggested an empirical modification of the form with Whereas Tarazona's functional was based on demanding that FMT reproduce certain exact limits, this functional, known as the "White Bear" functional, was mostly an empirical modification aimed at reproducing the more accurate Carnahan-Starling equation of state for the liquid phase. With this modification, the theory gives a quantitatively accurate desciption of the hard-sphere solid and of liquid-solid coexistence. Since the equation of state for the uniform fluid is no longer that of the Percus-Yevik approximation, it is necessarily the case that the implied direct correlation function is also different. In fact, the WBI model gives with η = π 6 nσ 3 and a (W BI) 0 For reference, we note that the Percus-Yevik direct correlation function has the same form but with coefficients a (P Y ) 0 The WBI direct correlation function is in fact more accurate than the PY direct correlationh function at high densities [12]. Finally, a similar functional was developed using a more accurate equation of state for mixtures, the so-called "White Bear II" model [13], It is worth pausing to note that it is easy to show that φ WBII (1−η) 3 ≥ 0 so increases monotonically from 0 and is therefore always nonegative and therefore so is φ WBII 2 (η). On the other hand, φ WBII 2 (η) has the well-behaved value φ WBII 2 (1) = 2 3 so that it does not affect the postivity of the third term in the Φ and it introduces no additional singular behaviour in the high-density limit.
As noted above, it has been reported in the literature that fully three-dimensional calculations for the hard sphere solid using the tensor theories must be carefully controled to preserve cubic symmetry or numerical instabilities are encountered [16]. In our own calculations, we find that the tensorial term does indeed lead to numerical instabilities. In particular, in the solid phase with its highly localized density peaks leading to local packing fractions η(r) very close to one, the instabilities are due to the numerator becoming negative leading to an overwhelming large negative contribution to the total free energy. The result is that the total free energy tends to negative infinity and so no physical minimum exists. The question then is whether this is a numerical problem, which can be corrected, or a property of the functional itself. We therefore first ask whether the numerator of the tensorial functional is positivedefinite, like the other FMT contributions, or not. If not, this opens the possibility that the numerical instabilities are in fact of the same origin as in Rosenfeld's original model. (Note however, that in earlier calculations using Gaussians to represent the solid density, Rosenfeld's model was unstable but the tensor models were not [11].) Returning to the example of the discontinuous density at a wall, one finds (see Appendix A for details) that the tensor density in the domain −σ/2 < z < σ/2 is so that the tensor numerator (from the third term in Eqs. (15) and (20) All of the factors in this expression are positive for z in the given domain except for the last one which becomes negative for z < −0.13σ. Another example is one with zero density outside the wall and a linearlly increasing density, ρ 0 z a , for z > 0 with a > σ/2. Again, in the domain −σ/2 < z < σ/2 one finds that which obviously assumes a small but negative value at z = 0. In conclusion, these examples serve to show that the third contribution to the tensor functional is not positive definite. Having also tried, without success, to eliminate the numerical instabilities by carefully controlling the numerics (i.e. by modifying the algorithm so as to preserve the inequalities between the weighted densities, by decreasing the lattice spacing and checking for round-off errors in the calculations) we conclude that the mix of empirical and analytic evidence suggests that the functional is indeed unbounded from below and that the numerical difficulties encountered when doing free minimization of the solid are due to this fact.

B. A bounded alternative
A truely bounded (from below) alternative was in fact proposed heuristically by Rosenfeld et al prior to the introduction of the tensor densities [29]. The idea was to replace the numerator in Rosenfeld's theory, (9), is now obviously non-negative and which agrees with Rosenfeld's theory when v 2 /s 2 << 1 -e.g. in the liquid limit. The latter property is necessary if one wishes to still obain the PY direct correlation function for the liquid. There is another way to motivate this approximation starting with the tensor theory [2]. One could try to replace the traceless tensor U by a combination of the scalar and vector quantities. Given that it is traceless and must scale linearlly with the density, any such approximation must have the form is an arbitrary scalar function. Evaluating the tensor numerator with this approximation gives Interestingly, the simplest choice A (x) = 1 yields the same positive-definite numerator. In the following, we refer to this as the RSLT functional. One can further "upgrade" to the Carnahan-Starling equation of state either by making the same approximation for the tensor density in the WBI and WBII functionals or by following the logic of their original derivations but starting with the positive-definite ansatz for the numerator: both routes yield the same result which is with either φ 1 and φ 2 being the same as in the tensor case (i.e. the WBI and WBII functions). Furthermore, it turns out that the direct correlation function for the fluid is then identical to that obtained from the tensor theories. It therefore seems that the differences between the two are subtle and probably only important in highly inhomogeneous systems like a solid. The upgraded RSLT model using the WBI functions will be the main alternative discussed in the remainder of this paper and will be referred to as the modified RSLT or mRSLT model.

C. Calculations: Canonical vs Grand Canonical interpretations
The goal is to minimize Λ [n; φ] at fixed chemical potential and external field resulting in a density n * (r) satisfying An alternative is to minimize F [n] at constant total particle number, N * , by constructing a Lagrangian and solving the equations Operationally, one could guess a value of λ and solve the first equation for the density, call it n λ (r) and then vary λ until finding a value that satisfies the second equation (call it λ N * ). The resulting density will satisfy and comparison to Eq.29 shows that this is precisely the local density for a system at chemical potential µ = λ N * . There are two ways to interpret this result. Rigorously, this means that the local density so obtained is the grand canonical local density for the system at this chemical potential. The fact that, as a matter of algorithmic convenience, it was obtained by minimizations at constant particle number is irrelevant. However, it is known [30] that for sufficiently large systems, the ensembles become equivalent so that one can, heuristically, take this local density to be that of a canonical system with fixed particle number N * . Practically, speaking there are in fact methods to evaluate the corrections involved in this identification but in practice they seem to be negligable for systems of 100 or more particles [31][32][33].
Note that the fixed-N procedure can be simplified by eliminating the Lagrange multiplier. Multiplying the first of Eqs.29 by the density and integrating, one sees that this particular density also satisfies so that one can write This means that n λ N * (r) can also be characterized as the solution to which gives a more practical route that does not involve a Lagrange multiplier and it is this equation that we actually solve in numerical calculations (when working at fixed particle number) via the implementation discussed in Appendix B. So, in this case, one first fixes the lattice parameter and the number of particles and determines the minimizing density, n λ N * ,a [r], giving Λ (N * , a) = Λ [n λ N * ,a ]. This then is minimized with respect to either one of the arguments, say N * , to get the minimizing value, N * * (a), and Λ (a) = Λ (N * * (a) , a). At the same time, one gets µ (a) = λ N * * (a) which gives the mapping between lattice parameter and chemical potential.

D. Validation for the homogeneous hard-sphere crystal
In order to validate our algorithm and to make contact with previous calculations, we have performed an extensive study of the homogeneous hard-sphere FCC crystalline phase for which there is no external field and the crystal extends indefinitely in all directions. Because there are no long-ranged interactions, it is sufficient to restrict computations to a single unit cell with periodic boundary conditions. (In fact, for simplicity, we use a cubic cell rather than the minimal, non-cubic cell for this lattice.) The minimal cubic unit cell has sides of length a and contains 4 lattice points so that the density of lattice points is n FCC = 4/a 3 . If every lattice site were occupied, this would also be the physical density of the crystal but in fact, one expects a certain fraction of unoccupied lattice sites (i.e. vacancies) so that the actual average density will be n < n FCC . It is related to the local density via where the notation indicates that the integral is restricted to a single unit cell. The number of vacancies in the unit cell will then be n FCC a 3 − na 3 = 4 − na 3 and the vacancy concentration (ratio of number of vacancies to the number of lattice sites) is c vac = 1 − n/n FCC . In principle, one should allow the lattice parameter to vary during the minimization. Instead, this minimization is divided into two steps: first, a minimization at fixed values of the lattice parameter to get n * a (r) and then a second minimization with of Λ [n * a ; φ] with respect to a. Our implementation of these calculations follows in broad outline previous work: the density field is discretized on a cubic computational lattice and the weighted densities calculated using fast Fourier Transforms. The main novelty of our implementation is that the rather than using the analytic Fourier Transform of the weighting functions as has been done previously, we use a scheme designed to guarantee that the various analytic inequalities between the weighted densities discussed above are rigorously preserved. This may result, for technical reasons discussed in Appendix C in somewhat slower convergence (i.e. the need for a finer lattice spacing) but ensures that no artifacts arise due to the FFT technique. A detailed description of the algorithm is given in Appendix C.
Because the results of the calculations are rather technical, we have collected most of the results in Appendix D and only briefly discuss the main results here. First, following Ref. [16] we have minimized at fixed-N and fixed lattice parameter with subsequent minimization with respect to the lattice parameter and examined convergence of the algorithm as the lattice parameter is refined. For sufficiently fine lattice spacing (up to 256 points per hard-sphere lattice parameter) we reproduce previous results for the WB models to several significant figures thus validating the implementation. Doing this for many different densities allows us to determine the phase coexistence with the liquid phase. For the WB models, we are again in close numerical agreement with previous calculations.
Concerning the numerical stability of the algorithm, we found, as reported previously, that numerical instabilities arose when we tried to minimize the WB functionals at constant chemical potential. On the other hand, at fixed particle number, we did not find it necessary to maintain strict cubic symmetry during the minimizations and we were able to use much coarser calculational lattices than were possible with previous implementations [34]. We speculate this may be due to the careful treatment of the FFT calculations and the preservation of the analytic inequalities.
We then extend the calculations to the proposed mRSLT model. In this case, we found that -as expected -we were able to minimize at fixed chemical potential with no instabilities as in the WB models. Also as expected, due to the more heuristic nature of hte model, the results for liquid-solid coexistence were, while reasonable, somewhat poorer. The liquid(solid) packing fractions at coexistence from simulation [35] are 0.492(0.545), the WBII model gives 0.495(0.545) and the mRSLT yields 0.513(0.546). While an error of some 3% is not excessive, there is nevertheless room for improvement.
A final test was to determine the vacancy concentration as a function of the average density. Oettel et al [16] report this as being a rather sensitive test of the models and is one for which the WBII funcitonal is clearly superior to the WBI functional (indeed, WBI does not produce physical results for this quantity because minimization with respect to the lattice spacing is not possible). The mRSLT functional does produce physically reasonable results that compare well with WBII as well as simulation.

A. cDFT Functional
To further illustrate the capabilities of a model that is free of divergences we consider a simple fluid. The development below is generic and can be adapted to any pair potential but for the sake of calculations, we have used a Lennard-Jones potential, which is cutoff at a separation r c and shifted to give the interaction potential All results reported here were obtained using a cutoff of r c = 3σ. The potential is separated into a short-ranged repulsive part, v 0 (r), and a long-ranged attractive part, w (r), using the WCA prescription [36,37] according to which v (r) = v 0 (r) + w (r) with where r 0 = 2 1/6 σ is the minimum of the potential. An effective hard-sphere diameter, d, is computed using the Barker-Henderson prescription [37,38], and it should be noted that the effective hard-sphere diameter depends on temperature but not on the density. The model DFT Helmholtz functional is then constructed as For a uniform fluid, n (r) = n, all of the Helmholtz free energy assumes akind of generalized van der Waals form, with a = w (r) dr.
To implement this functional, only two changes from the hard-sphere calculations are necessary. The first is that one must include the extra term coming from the attractive part of the potential. This is efficiently evaluated using standard FFT techniques. The second is that care must be taken with the periodic boundaries. Whereas in the hard-sphere case, it sufficed to take the calculational cell to be the cubic unit cell (and in fact one can invoke cubic symmetry to reduce this to an octant of the cubic cell, although we did not do this) the attractive tail will in general have a range greater than the lattice spacing and so a larger calculational cell must be used. (Again, other, more efficient but more complicated possibilities exploiting cubic symmetry and not using the minimum image prescription are possible but not pursued here.)

B. Phase diagram
As a first application, the model was used to calculate the vapor-liquid-FCC phase diagram for a homogeneous system. Since the typical lattice spacing in an FCC Lennard-Jones solid is on the order of σ, we used a calculational cell consisting of 4×4×4 primitive cubic unit cells with a discretization lattice spacing of ∆ = 0.05σ. The vapor-liquid binodal is easily determined from the model. Indeed, for a uniform system at densityn, one has that where the Carnahan-Starling excess free energy per unit volume is which is the usual thermodynamic relation between the chemical potential and the free energy. Coexistence is therefore determined by finding the (two) densities,n 1 andn 2 satisfying Eq. (44) and, since the fundamental theorem of DFT says the equilibrium state minimizes Λ[n], the two states must also satisfy Λ(n 1 ) = Λ(n 2 ). These two relationsequivalent to the usual thermodynamic criteria of equality of chemical potential and pressure -uniquely determine the coexisting states. For the solid, the calculation is conceptially the same. First, a chemical potential µ is chosen. Then, Λ[n] is minimized for unit cells of size k∆ with typical values being 127 ≤ k ≤ 135. Interpolation of the resulting values of Λ are used to determine a minimum and, hence, the cell size (i.e. lattice density) and grand-canonical free energy for the solid phase for the given chemical potential. The same interpolation also yields the average density of the solid phase. From Eq. (44), the density of the liquid (vapor) for the chemical potential is also determined. Finally, the procedure is repeated for different chemical potentials until equality of the grand canonical free energies is found, thus giving the coexisting densities for the solid and liquid (vapor) phases. The resulting phase diagram is shown in Fig.1 where the reduced temperature is T * ≡ k B T / . The results are typical of previous calculations and serve as a baseline for the application of this model to crystallization.

C. Heterogeneous Crystallization
As an example of the capabilities of a well-behaved DFT functional we consider the rapid quench of a small Lennard-Jones droplet attached to a wall. The calculational cell is rectilinear with dimensions 20σ × 20σ × 10σ with periodic boundaries in the first two (x,y) directions. In the z-direction, the cell is bounded by impenatrable walls. The lower wall at z = 0 is hard meaning that the particles do not interact with it except that they cannot pass through it: it can be represented as an external potential which is zero for z > 0 and infinite for z < 0. The interaction with the upper wall is also taken to be infinite for z > z wall and for z < z wall it is the Steele 4-9 potential [39,40], which is a continuum potential intended to mimic the average potential generated by the (100) surface of an FCC crystal, where the position of the wall is z wall = 9σ. (The wall is not at the limit of the cell, 10σ, for technical reasons since the FMT weighted densities are non-zero for one hard-sphere radius outside the calculational domain.) The calculations were performed using wall = 0.15 , σ wall = σ and a discretization of ∆ = σ/7 which for a typical solid density of nσ 3 = 0.9 corresponds to approximately 11.5 lattice points per unit cell. The initial density was taken as n(r) = n vap + (n liq − n vap )e −βv wall (r) Θ (R − (r 0 − r)) with n vap , n liq being the coexistence densities of liquid and vapor at the given temperature, R is the radius of the (spherical) droplet and r 0 = (10σ, 10σ, z wall ) is its center. This is therefore a crude model of a spherical drop resting on the upper wall. The initial temperature was taken to be k B T = 0.8 and the cDFT functional minimized at constant particle number. Because of the rather large system sizes needed for the calculations, to ensure that there is no direct self- interaction of the droplets, we have used a rather coarse lattice with a lattice spacing of σ/7 in this and all subsequent calculations described below.(There are approximately 297 particles in the system with almost all concentrated in the droplet.) As discussed above, the resulting configurations can be viewed in two different ways. On the one hand, the droplets are expected to approximate the density distribution for a canonical system, with the approximation generally becoming more accurate as the number of particles in the system increases. On the other hand, the result of minimizing at constant particle number has the exact interpretation of being the stationary density distribution for the grand canonical system with a chemical potential determined by the calculation. In the present case, it corresponds to an unstable stationary point of the grand-canonical functional and, so physically, to the grand-canonical critical droplet.
Once this configuration is relaxed, the temperature is lowered and the system again relaxed, etc. Since there are no fluctuations in the system (except those induced by numerical noise and the discreteness of the calculational lattice), this corresponds physically to a rapid quench. Figure 2 shows the evolution of the density as the temperature is lowered: i.e. according to the two interpretations it is either the evolution of a quenched system in the canonical ensemble or a sequence of critical droplets in the grand canonical ensemble. As the temperature is lowered, one sees the development of more and more structure within the droplet corresponding to layering of the dense fluid near the wall until, at the lowest temperature, the density spontaneously localizes and the system crystallizes into a hexagonal close-packed (HCP) structure. This crystallization is completely unconstrained and untemplated: it is not induced artificially and the resulting structure is distinct from the cubic calculational lattice. It is also surprising since the minimum energy state for bulk Lennard-Jones is an FCC crystal. However, the difference in energy between HCP and FCC is small and in fact it is known from simulations that small droplets attached to walls will crystallize into HCP structures [41]. The crystal can be reheated, with minimization of the energy functional at each temperature, and the free energy difference between the droplet and the crystallite determined until at k B T = 0.55 the crystal spontaneously melts. For the range 0.45 < k B T / < 0.55 both the droplet and the crystallite are metastable thus establishing that there is a free energy barrier separating them: the free energy functional therefore gives a mean-field model of crystallization. At k B T / = 0.50 the total free energy of the system with crystallite is about 7k B T smaller than that with the droplet. Further properties of this model, including a demonstration of the spontaneous formation of the bulk FCC phase, will be discussed in a later publication.

D. Homogeneous Solidification: Polymorphic Behaviour
Encouraged by this, we turned to the same problem in the absence of a wall with the aim of studying homogeneous crystals. We began with a free-floating droplet that was equilibrated at a temperature above the triple point. We then lowered the temperature and equilibrated resulting in the series of droplets shown in Figure 3. At a sufficiently low temperature, we again saw the spontaneous formation of a localized structure. To characterize it, localized peaks obtained in DFT are identified as particle coordinates using particle tracking method that are commonly employed for confocal microscopy analysis [42] [See Fig.4]. From there, we deduced the pair distribution function which exhibit a split of the second peak [See Fig. 5.a]. Furthermore, topological cluster classification analysis as developed by Mallins et al. [43] is carried out in order to identify local structures. The amorphous cluster is characterized by an increasing emergence of icosahedral order as the temperature is decreased. Icosahedral order has been conjectured as a signature of the glass transition. In this picture, geometrical frustrations are induced because icosahedrons are not able to tessellate 3D space which then lead to the dynamical arrest observed in glass transition [44]. While both experimental [45,46] and numerical evidences [47][48][49][50] have been found for such phenomenon, our results demonstrate that the icosahedral rich structure is not only accompanying or causing the dynamical arrests but it is a genuine minimum in the free energy landscape and should be therefore considered as a thermodynamically stable state. The conclusion that the system has found a free energy minimum corresponding to a glassy structure is plausible. In order to check that this was not a pathology of the model, we performed the same calculation but at the beginning of the minimization at the lowest temperature, we introduced a simple planar field at the center of the droplet which was intended to break the spherical symmetry. (The specific field was arbitrarily chosen to be φ(r)/ = −0.1e |z|/10 where the droplet is centered at x = y = z = 0. ) This field was removed after a few thousand minimization steps and the minimization was allowed to continue as before. In this case, we obtained a nearly perfect FCC crystal with a free energy some 222k B T below that of the amorphous cluster. The crystal is shown in Fig. 6 where its faceting into (111) and (100) planes is evident. As illustrated in Fig.5, both hcp and fcc clusters are found within the crystal with far higher populations than for the amorphous cluster. From these calculations, we already know that the crystalline and amorphous clusters are both minima of the free energy functional and so there must necessarily be a free energy barrier between them. To investigate their stability, we then increased the temperature and minimized starting in each case with the previous (lower temperature) structure. The resulting free energies are shown in Fig.7. The amorphous structure spontaneously melts at k B T / = 0.45 indicating that the energy barrier separating it from the liquid state vanishes. The solid is stable up to k B T / = 0.45 at which point it also melts.
We note that it is widely known that glasses are not normally seen in simulations of one-component Lennard-Jones systems and so our amorphous clusters might seem suspicious. However, we have verified that using the particle coordinates as input to a (constant-N, canonical) molecular dynamics simulation, the clusters immediately crystallize so that the energy barrier between the amorphous and crystalline states must be very small compared to the thermal energy, thus explaining why they are not normally observed. We also note that our crystalline cluster does not have the geometry of the minimum energy cluster found by Xing et al [51] and reported in the Cambridge Cluster Database [52]. To investigate this, we have again used our structures as the starting state for the minization of the Hamiltonian in a particle-based algorithm (using facilities in the LAMMPS package [53]) and we have performed the same calculation using the database structure (but using our cut and shifted Lennard-Jones potential). We find that the energy per particle of the amorphous cluster converges to −6.20 , our FCC crystal gives −6.48 and the database structure gives −6.50 so that while the latter is, technically, a more stable state, the energy difference is very small. Finally, it is important to be clear that our calculations do not, per se, involve nucleation because there are no thermal fluctuations: we quench to sufficiently low temperatures that the barrier to crystallization either vanishes or becomes so small that numerical noise is enough to drive the transition. But such a model can obviously be used together with techniques for exploring free energy surfaces to give a complete picture of crystallization. It should be emphasized that DFT by itself is not sufficient to resolve all issues associated with the nucleation of stable solid clusters. This is because nucleation is a fundamentally non-equilibrium process and, has been recently empahsized [54,55], requires a dynamical description in addition to the free energies provided by DFT. However, it has recently become clear how to combine DFT and dynamics -e.g. fluctuating hydrodynamics -so as to describe first order phase transitions [55]. The relation between this dynamical description and the more coarse-grained approach of Classical Nucleation Theory has been established as has the development of reduced description based on coarsegrained order parameters [55,56]. With the addition of free energy functionals capable of describing crystalline phases without constraint, the path is now open to attach theoretical questions that have until now been inaccessible.

IV. CONCLUSIONS
In this paper, the limits of current FMT models for hard-sphere crystallization have been re-examined. It was demonstrated that the tensor theories are not positive definite and, based on the numerical evidence, it was suggested that they are indeed not bounded from below thus opening the possibility that the hard-sphere crystalline states they predict are in fact metastable states. The true minima, in this case, consist of some sort of highly localized densities with local packing fractions approaching unity that give an unbounded negative total free energy becoming lower as the density peaks become more localized. In other words, the functionals are unbounded from below. As such, they can only be viewed as very limited realizations of the "true" DFT functional which should give the (known) stable crystalline state as an absolute minimum of the functional. The correct behaviour has been illustrated by using an older, positive-definite, and therefore bounded, excess free-energy functional that is adapted to give the Carnahan-Starling equation of state. This functional is certainly not the final word in hard-sphere functionals as it lacks several advantages of the tenor theories: in particular, the description of hard-sphere freezing is not as quantitatively accurate as that of the tensor theories and (probably related) it lacks consistency with some known exact results that the tensor theories respect [11]. Nevertheless, the virtues of being demonstrably bounded are also very clear: namely, the ability to minimize with no special regard to maintaining cubic symmetry and even the possibility to minimize at constant chemical potential with no additional effort.
Ours is not the first work to encounter problems with the tensorial functionals. We note that extensions to nonspherical hard particles have also proven difficult as discussed, e.g. by Wittman et al [57] who observed that some of these problems are solved by using another, older, functional of Tarazona and Rosenfeld [58] which is also demonstrably bounded from below and which could therefore be an interesting alternative to the RSLT functional used here. In any case, none of this definitively shows that the (newer) tensorial functionals are unphysical. The present work only shows that (a) the numerator of the funcitonal is not positive-definite and (b) that numerical evidence suggests it could be unbounded from below. The second point is supported by the fact that simply replacing the functional by the bounded mRSLT functional with no other change to the calculations gives numerical stability. While our algorithmic implementation is nonstandard and could also be responsible for some of the improved stability, we did find the same instabilities with our algorithm as with older ones.
The mRSLT functional was then used as the basis for a mean-field treatment of a Lennard-Jones system. The full vapor-liquid-solid phase diagram was computed including full minimization of the density. Then, as a first illustation of the novel applications that become possible with a stable functional, the crystallization of a droplet on an attractive wall was illustrated. The density of the droplet was calculated (at constant particle number) for lower and lower temperatures until at sufficiently low temperature, the droplet crystallzes into an HCP structure. This localization of the density and crystallization is completely spontaneous and takes place with no externally-imposed lattice structure. Because these are static calculations, the crystallization is not the result of nucleation: rather, the temperature is so low that any energy barrier between the liquid and solid states is so small that either (a) there is in fact no barrier or (b) there is a barrier but it is so small that numerical noise in the calculation is enough to allow the system to overcome it. Nevertheless, despite being performed at constant particle number, the relevance to nucleation at constant chemical potential was noted so that the present results do in fact have a bearing on nucleation.
The study of clusters was continued by replacing the wall with periodic boundaries. In this case, an amorphous cluster was produced. We also succeeded in generating a lower energy FCC crystalline cluster by introducing a symmetry-breaking external field in the early stages of minimization and then removing it. These clusters were heated and we demonstrated a region of three-phase coexistence until at higher temperatures the amorphous and crystalline clusters eventually melted. This work continues and we are in the process of studying similar clusters of different sizes. We emphasized that the calculations can be viewed as either approximate (but for such large systems, probably very good) results for a fixed particle number canonical system or as exact results in the grand-canonical ensemble. For the canonical systems, the droplets and clusters are stable states due to the constraint of fixed N : if the clusters grow, the density of the vapor outside of them decreases and the system becomes understaturated leading to cluster sublimation while if the clusters shrink, the vapor becomes more highly supersaturated leading to cluster growth [59]. Viewed as exact grand-canonical results, the clusters are metastable critical clusters and as such give a direct indication of the nucleation barrier.
These calculations show that cDFT can be a powerful tool in the study of inhomogeneous classical systems with the capability of producing unexpected results. It is hoped that the potential of this application will serve to inspire more work on how to improve the FMT theories so that they retain their quantitative accuracy while correcting their weaknesses which seem to be tied to the issue of boundedness of the functional. Promising starting points for such developments could be the older Tarazona-Rosenfeld tensorial functional [58] discussed also by Wittman et al [57] and the recently proposed scalar functional of Hansen-Goos et al [60] built upon a new consistency condition discovered by Santos [61]. The potential applications to nucleation, pre-melting of surfaces, cluster formation in proteins, wetting, crystallization and perhaps even the glass transition are evident. In particular, they open the door to a complete approach to the description of first order phase transitions in which the free energy surfaces produced can be combined with a stochastic-dynamical framework [54,55] and tools such as the dimer method [62] and the string method [63], as well as stochastic-processes theory [64][65][66], so as to predict nucleation rates, nucleation pathways and multi-step nucleation mechanisms. and Appendix B: Calculations at fixed particle number An alternative to the introduction of a Lagrange multiplier to fix the particle number is to introduce an auxilliary field x (r) and to write thus constructing densities with fixed particle numberber. Minimizing F [n N * ] + φ (r) n N * (r) dr with respect to x gives and using This is exactly the same as the Eq.(33) derived using a Lagrange multiplier.
Appendix C: Algorithm

Introduction
This Appendix explains the details of the algorithm used in our calculations. The cDFT functional can be written as with the ideal gas contribution To evaluate this numerically, we discretize space into a set of points R = i∆ x + j∆ y + k∆ z where ∆ is the lattice spaceing and i, j, k ∈ Z. We assume periodic boundaries so that there are N x unique points in the x-direction, 0 ≤ i < N x ,etc. The value of the density at lattice point R will be denoted n R . The thre contributions can be evaluated using, e.g., the simplest trap-rule expressions and the only remaining question is the evaluation of local weighted densities.

Weighted Densities
The weighted densities all have the form of convolutions; where the weighting functions, w (α) (r), are relatively localized: w (η) is only nonzero inside a volume of π 6 σ 3 and the others are nonzero only on the shell defined by r = σ 2 . So one possibility would be to directly evaluate them in real space. The cost would be on the order of N x N y N z σ ∆ 3 operations for η (r) and N x N y N z σ ∆ 2 for the others. Note that in the latter case, there is the complication that the domain of integration -the shells of radius σ/2 -will, in general, contain no lattice points so this is not so simple as for the case of the volume integral needed for η which can be approximated by simply summing over the lattice points in the volume. An alternative would be to rewrite the convolution in Fourier space as The analytic Fourier transform of the weight functions is easily calculated and the that of the density can be obtained by means of Fast Fourier Transform (FFT) at a cost of on the order of (N x N y N z ) log 2 (N x N y N z ) operations. The product involves (N x N y N z ) operations and the inverse FFT needed to get the real-space weighted densities gives a total cost of 2 (N x N y N z ) log 2 (N x N y N z ) + (N x N y N z ) operations. The FFT method is advantageous provided that or, if there are the same number of points in each direction, N < 2 ( σ ∆ ) 3 /3 . Clearly, the advantage of the FFT method increases as the lattice is refined but even for as few as 8 points per hard-sphere diameter the FFT method is advantageous provided the number of latticepoints in each direction is N < 2 171 .

Real space or Fourier space weighting functions?
Thus, the FFT method must clearly be the method of choice in practical calculations. However, if implemented as described here, there are potential problems. If the FT of the weighting functions is determined analytically and that of the density numerically, then there is no guarantee that the resulting weighted densities will be non-negative or that the inequalities discussed in the main text will be respected (e.g. that s 2 R > v 2 R ). Since one of our main goals was to investigate the stability of the models, a more careful implementation of this approach -one guaranteed to respect these properties -was used. Returning to the real-space version of the calculation, let us assume that there is some set of points r a and weights u a with 0 ≤ a < M that allow us to approximate an integral over a spherical shell of radius t as Then, the values of the weighted density s (r) at the lattice points could be evaluated as Since we only know the density on the lattice points, the value of the density on the shell must be evaluated with some sort of interpolation. The only requirement for the following is that the method of interpolation be linear in the density: we make the simplest choice of trilinear interpolation. For any point r in the domain, a lattice point S (r) = i (r) ∆ x + j (r) ∆ y + k (r) ∆ zcan be identified having the property that etc. and the value of the density at r is estimated by doing linear interpolation with resepct to the eight corners of this cube. Defining etc., the result is that Finally, noting that for any lattice vector R, we have that S (R + r) = R + S (r), it follows that ∆ x (R + r) = ∆ x (r) so that n (R + r) and Analogous expressions for the other shell-averaged densities are easily found. For the local packing fraction, we obtain something similar by integrating this shell expression with respect to the radial variable using a one-dimensional Gauss-Legendre scheme. These real-space, discretized expressions for the weighted densities satisfy all of the physical inequalities and the discrete convolutions can be efficiently (and, up to numerical noise, exactly) evaluated using discrete FFT.

Spherical Integrals
One question not addressed so far is what points are used to do the spherical integrals. If one naively applies one-dimensional integration schemes to spherical coordinates, its easy to see that spatial asymmetries are created that could bias the calculations. The standard method [67], used for example in quantum-chemical calculations, are to use points and weights some subset of spherical harmonics exactly. The two most common schemes are Lebedev quadrature and Spherical t-schemes [67]. The difference between the two is that the former strictly enforce octahedral and inversion symmetry whereas the latter are constructed by demanding that the weights of all points are equal. The efficency of the two schemes is similar [67]. Both have been tested and the differences in the present calculations are negligable. When using spherical t-schemes, we have enforced cubic and reflection symmetry by expanding the list of points first so that besides each point (x, y, z) we also have all permutations of the same values of x, y, z (thus going from one to 6 points). We then take this expanded list of points and include all reflections: (x, y, z), (−x, y, z)...(−x, −y, z)...(−x, −y, −z). In total, this expands the list of N points to N × 6 × 8 = 48N . In our calculations, we have used point sets from the website Ref. [68,69]. The results seem to be robust with respect to the number of points used: we have typically used sets containing 21,118 points. The volumetric spherical integrals are performed by supplementing the integration on the spherical shell with Gauss-Legendre integration in the radial coordinate using up to 512 points. Numerical calculations were implemented using the GSL scientific library [70] and the Armadillo linear algebra library [71] and Fourier transforms perfomed using the FFTW library [72].

Minimization
After discretization, the ftDFT functional becomes a function of the densities at each lattice point. The gradients of this function are easily calculated "analytically" (i.e. analytic expressions that are calculated at the same time as the functional itself). Indeed, the discretized Λ functional has the form where the sums are over all lattice points, w ik is the attractive part of the potential evaluated at the lattice sites, v ext i is the external potential evaluated on the lattice sites and n is the collection of values of the density on the lattice points. Minimization requires that or, upon rearrangement, These forces were used to implement either conjugate gradient (CG) minimization [73] or the Fast Inertial Relaxation Engine (FIRE) of Bitzek et al [74]. For the liquid, both methods acheived convergence quickly and with similar effort. For the solid phase, the CG would sometimes fail to converge or converge very slowly while FIRE always converged and generally with much less effort. Minimization was terminated using the criterion where the last approximation is valid when the forces are small (the case of interest) and avoids the expense of exponentiating over all of the lattice points. We found that convergence was generally achieved for = 10 −2 and, to be sure, the criterion = 10 −4 was used throughout the calculations. We have implemented these algorithms by discretizing the density on a cubic grid (details are given in Appendix C). Aside from the physical parameters, the only numerical parameter is the grid spacing, ∆. The convergence of the calculations as the lattice spacing is decreased (at fixed particle number and using the WBII model) is illustrated by the results in Table I which includes an extrapolation to the continuum limit. Even the coarsest grid, with only 8 points in each Cartesisan direction, corresponding to about 5 points per hard sphere diameter, gives free energies accurate to about 6%. At 256 3 grid points, the results are within 0.2% of the continuum limit. Compared to Ref. [16], which are presumeably converged to several significant figures, the present calculations converge more slowly. Most likely this is due to specific design goals in our algorithm described in Appendix C which lead to the use of linear interpolation and simple trap-rule integration. Table II shows results from calculations using the WBII functional at constant particle number at different densities. In each case, the average density,n, and the vacancy concentration, c vac , was fixed thus determining the total number of particles in a unit cell as N = 4 × (1 − c vac ) and thus the lattice parameter via 4(1 − c vac )/a 3 =n. Varying the vacancy concentration while holding the average density constant is therefore equivalent to varying the lattice parameter. Also shown are results reported in Ref. [16], where the authors state that they used 64 3 to 256 3 grid points per unit cell (increasing as the density increases) and presumeably, the numbers given are converged to a several significant digits. The present calculations using 256 3 grid points are consistent with these values with the typical difference being on the order of 0.1%. Table III shows the coexistence conditions predicted by the various theories. While the values are somewhat dependent on the grid spacing, they are largely insensitive to the vacancy concentration. The values of the postivedefinite (mRSLT) theory are also given. Even though the latter incorporates the Carnahan-Starling equation of state, and is therefore equivalent to the WBI theory for the liquid phase, it is less accurate in predicting coexistence due to higher overall predictions for the solid pressure. It therefore lacks the qualitative accuracy of the tensor theories.

Vacancy concentration
The results for freezing are relatively insensitive to the vacancy concentration. Nevertheless, DFT formally requires that the Λ-functional be fully minized with respect to the density. In the present context, this means that it must be minimized with respect to the vacancy concentration. As discussed in detail in Ref. [16], such a minimization is not even possible with functionals such as WBI and the fact that WBII can be minimized and gives physically reasonable results is one of its strongest features. In the case of the RSLT functional, there is no doubt that it can be minimized  [16], these calculations were performed at fixed particle number with a constant vacancy concentration of cvac = 1×10 −4 . The columns are the lattice density, the lattice packing fraction, the number of points in each Cartesian direction of the unit cell and the Helmoltz free energy per particle reduced by the temperature. All calculations were performed using the WBII model. with respect to the average density (at fixed lattice constant): this follows simply because as the local density increases, it must at some point cause the local packing fraction to increase beyond unity and this will necessarily lead to a (positive) divergence in the functional. Figure 8 shows the results of such a minimization, together with the WBII results and data from simulations. The mRSLT theory gives physically reasonable results at all densities.
Even more impressive in this context is the fact that there are no technical difficulties in minimizing at constant chemical potential rather than constant average density (or, equivalently, constant vacancy concentration). Even for WBII, Oettel et al state that such a minimization at constant chemical potential is not technically feasible and we have confirmed with the present code that attempts to do so typically lead to (negative) divergence of the energy functional. This is not the case with the mRSLT model and the Fig. 8 shows that the resulting vacancy concentrations are consistent with those obtained from the constant density minimizations.  8. Vacancy concentration as a function of average density determined using the mRSLT model both by minimizing the free energy at constant particle number and at constant chemical potential with subsequent minimization with respect to the lattice spacing in both cases. Simulations results due to Kwak et al [75] and Bennett and Alder [76] are shown as are the WBII results of Oettel et al [16].