# Finite-size corrections to the free energies of crystalline solids

## Abstract

We analyse the finite-size corrections to the free energy of crystals with a fixed center of mass. When we explicitly correct for the leading ($`\mathrm{ln}N/N`$) corrections, the remaining free energy is found to depend linearly on $`1/N`$. Extrapolating to the thermodynamic limit ($`N\mathrm{}`$), we estimate the free energy of a defect free crystal of particles interacting through an $`r^{12}`$-potential. We also estimate the free energy of perfect hard-sphere crystal near coexistence: at $`\rho \sigma ^3`$=1.0409, the excess free energy of a defect-free hard-sphere crystal is $`5.91889(4)kT`$ per particle. This, however, is not the free energy of an equilibrium hard-sphere crystal. The presence of an finite concentration of vacancies results in a reduction of the free energy that is some two orders of magnitude larger than the present error estimate.

The earliest numerical technique to compute the free-energy of crystalline solids was introduced some 30 years ago by Hoover and Ree . At present, the “single-occupancy-cell” method of Ree and Hoover is less widely used than the so-called “Einstein-crystal” method proposed by Frenkel and Ladd . The latter method employs thermodynamic integration of the Helmholtz free energy along a reversible artificial pathway between the system of interest and an Einstein crystal. The Einstein crystals serves as a reference system, as its free energy can be computed analytically. Since its introduction, Einstein-crystal method has been used extensively in studies of phase equilibria involving crystalline solids. For numerical reasons - to suppress a weak divergence of the integrand - the Einstein-crystal method calculations have to be carried out at fixed center of mass. The free energy of the reference crystal is also calculated under the center-of-mass constraint, and the final calculated free energy of the unconstrained crystal is determined by correcting for the effect of imposing the constraint in the calculations. In the original paper, the fixed center-of-mass constraint was only applied to the particle coordinates, but not to the corresponding momenta. This is irrelevant as long as one computes the free-energy difference between two structures that have either both constrained or both unconstrained centers of mass. However, when computing the absolute free energy of a crystal, one needs to transform from the constrained to the unconstrained system. In the original paper, this transformation was not performed consistently. This resulted in a small but noticeable effect on the computed absolute free energy of the crystal. Below, we describe the proper approach to calculate the free energy of arbitrary molecular crystalline solids. The derivation differs from the earlier work in two respects: first of all, we explicitly show the effect of momentum constraints. And secondly, we generalize the expression to an arbitrary crystal containing atoms or molecules with different masses.

The main point of interest involves the calculation of the partition function of a crystal with and without a constrained center of mass. The partition function for an unconstrained, $`d`$-dimensional crystalline solid of $`N_{mol}`$ molecules composed of a total of $`N`$ atoms is given by

$$Q=c_Nd^{dN}\stackrel{}{r}d^{dN}\stackrel{}{p}\mathrm{exp}[\beta (\stackrel{}{r}_i,\stackrel{}{p}_i)]$$
(1)   
where $`c_N`$=$`(h^{dN_{mol}}N_1!N_2!\mathrm{}N_m!)^1`$, where there are $`N_1`$ indistinguishable molecules of type 1, $`N_2`$ molecules of type 2, etc., where $`N_1+N_2+\mathrm{}+N_m=N_{mol}`$, and $`h`$ is Planck’s constant. It should be noted that, in all calculations of phase equilibria between systems that obey classical statistical mechanics, Planck’s constant drops out of the result. Hence, in what follows, we omit all factors $`h`$. Using the result of the appendix in an article by Ryckaert and Ciccotti , one can show that the constrained partition function $`Q^{con}`$ is given by

$`Q^{con}`$ $`=`$ $`c_N{\displaystyle d^{dN}\stackrel{}{r}d^{dN}\stackrel{}{p}\mathrm{exp}[\beta (\stackrel{}{r}_i,\stackrel{}{p}_i)]}`$ (2)   
$`\times `$ $`\delta (\stackrel{}{\sigma }(\stackrel{}{r}))\delta (𝐆^1\stackrel{}{\dot{\sigma }})`$ (3)   
where $`\stackrel{}{\sigma }(\stackrel{}{r})`$ and $`\stackrel{}{\dot{\sigma }}`$ are the constraints and time derivatives of the constraints, respectively, and

$$G_{ij}=\underset{k=1}{\overset{N}{}}\frac{1}{m_k}\frac{\sigma _i}{\stackrel{}{r}_k}\frac{\sigma _j}{\stackrel{}{r}_k}$$
(4)   
The same integration limits implicit in Eqn. (1) are also used in Eqn. (3). To constrain the center of mass (CM), we take $`\stackrel{}{\sigma }(\stackrel{}{r})=_{i=1}^N\mu _i\stackrel{}{r}_i`$, and thus, $`\stackrel{}{\dot{\sigma }}=_{i=1}^N(\mu _i/m_i)\stackrel{}{p}_i`$, where $`\mu _im_i/_im_i`$. Note that in Eqns. (1) and (3) we have assumed that there are no additional internal molecular constraints, such as fixed bond lengths or bond angles.

We first consider the case of an Einstein crystal, which has a potential energy function given by $`U_{Ein}=(\alpha /2)_{i=1}^N(\stackrel{}{r}_i\stackrel{}{r}_i^{(0)})^2`$, where $`\stackrel{}{r}_i^{(0)}`$ are the equilibrium lattice positions. Note that the particles in a crystal are associated with specific lattice points and therefore behave as if they are distinguishable - thus, $`c_N=1`$ (as we omit the factor $`1/h^{d(N1)}`$). It is easy to show that

$$Q_{Ein}^{CM}=Z_{Ein}^{CM}P_{Ein}^{CM},$$
(5)   
with

$`Z_{Ein}^{CM}`$ $`=`$ $`{\displaystyle d^{dN}\stackrel{}{r}\underset{i=1}{\overset{N}{}}\mathrm{exp}[(\beta \alpha /2)r_i^2]\delta (\underset{i=1}{\overset{N}{}}\mu _i\stackrel{}{r}_i)}`$ (6)   
$`=`$ $`\left({\displaystyle \frac{\alpha \beta }{2\pi _i\mu _i^2}}\right)^{d/2}\left({\displaystyle \frac{2\pi }{\alpha \beta }}\right)^{Nd/2}`$ (7)   
$`=`$ $`\left({\displaystyle \frac{\alpha \beta }{2\pi _i\mu _i^2}}\right)^{d/2}Z_{Ein}`$ (8)   
and

$`P_{Ein}^{CM}`$ $`=`$ $`{\displaystyle d^{dN}\stackrel{}{p}\underset{i=1}{\overset{N}{}}\mathrm{exp}[(\beta /2m_i)p_i^2]\delta (\underset{i=1}{\overset{N}{}}\stackrel{}{p}_i)}`$ (9)   
$`=`$ $`\left({\displaystyle \frac{\beta h^2}{2\pi M}}\right)^{d/2}{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{2\pi m_i}{\beta h^2}}\right)^{d/2}`$ (10)   
$`=`$ $`\left({\displaystyle \frac{\beta h^2}{2\pi M}}\right)^{d/2}P_{Ein}\text{ ,}`$ (11)   
where $`M=_im_i`$, while $`Z_{Ein}`$ and $`P_{Ein}`$ are the corresponding contribution to $`Q_{Ein}`$,the partition function of the unconstrained Einstein crystal. Clearly,

$$Q_{Ein}^{CM}=(\underset{i}{}m_i/\underset{i}{}m_i^2)^{d/2}(\beta ^2\alpha h^2/4\pi ^2)^{d/2}Q_{Ein}$$
(12)   
Similarly, one can show that the partition function for an arbitrary crystalline system subject to the CM constraint is given by

$$Q^{CM}=Z^{CM}(\beta h^2/2\pi M)^{d/2}\underset{i=1}{\overset{N}{}}(2\pi m_i/\beta h^2)^{d/2},$$
(13)   
with

$$Z^{CM}=d^{dN}\stackrel{}{r}\mathrm{exp}[\beta U(\stackrel{}{r}_i)]\delta (\underset{i=1}{\overset{N}{}}\mu _i\stackrel{}{r}_i)$$
(14)   
while the partition function of the unconstrained crystal is given by

$$Q=Z\underset{i=1}{\overset{N}{}}(2\pi m_i/\beta h^2)^{d/2},$$
(15)   
with

$$Z=d^{dN}\stackrel{}{r}\mathrm{exp}[\beta U(\stackrel{}{r}_i)]$$
(16)   
Note that, as far as the kinetic part of the partition function is concerned, the effect of the fixed center of mass constraint is the same for an Einstein crystal as for an arbitrary “realistic” crystal.  Using Eqns. (13) and (15), the Helmholtz free energy difference between the constrained and unconstrained crystal is given by

$`\beta (FF^{CM})`$ $`=`$ $`\mathrm{ln}(Q/Q^{CM})`$ (17)   
$`=`$ $`\mathrm{ln}(Z^{CM}/Z){\displaystyle \frac{d}{2}}\mathrm{ln}(2\pi M/\beta h^2)`$ (18)   
We note that

$`{\displaystyle \frac{Z^{CM}}{Z}}`$ $`=`$ $`{\displaystyle \frac{d^{dN}\stackrel{}{r}exp[\beta U(\stackrel{}{r}_i)]\delta \left(_i\mu _i\stackrel{}{r}_i\right)}{d^{dN}\stackrel{}{r}\mathrm{exp}[\beta U(\stackrel{}{r}_i)]}}`$ (19)   
$`=`$ $`\delta ({\displaystyle \underset{i}{}}\mu _i\stackrel{}{r}_i)=𝒫(\stackrel{}{r}_{CM}=\stackrel{}{0})`$ (20)   
where $`\stackrel{}{r}_{CM}_i\mu _i\stackrel{}{r}_i`$, and $`𝒫(\stackrel{}{r}_{CM})`$ is the probability distribution function of the center of mass, $`\stackrel{}{r}_{CM}`$.

To calculate $`𝒫(\stackrel{}{r}_{CM})`$ we exploit the fact that the equilibrium crystal lattice is invariant to translations over displacements through linear combinations of integer multiples of the lattice vectors. This is true if the crystal lattice is subject to periodic boundary conditions. Consequently, the probability distribution of the center of mass of the lattice is evenly distributed over a volume equal to that of the Wigner-Seitz cell of the lattice positioned at the center of the volume over which we carry out the integration in the partition function. Since the average center of mass of the crystal is equal to the center of mass of the lattice, it follows that $`𝒫(\stackrel{}{r}_{CM})=1/V_{ws}=N_{ws}/V`$, where $`V_{ws}`$ is the volume of a Wigner-Seitz cell, and $`N_{ws}`$ is the number of such cells in the system. Thus, $`Z^{CM}/Z=𝒫(\stackrel{}{r}_{CM}=\stackrel{}{0})=N_{ws}/V`$ In the case of one molecule per cell, this implies $`Z^{CM}/Z=N_{mol}/V`$, where $`N_{mol}`$ is the number of molecules in the system.

In the Frenkel-Ladd free energy calculation, the free energy difference between the constrained crystal and the reference systems is given by

$$\beta F^{CM}=\beta F_{Ein}^{CM}\beta _0^1𝑑\lambda \mathrm{\Delta }U_\lambda ^{CM}$$
(21)   
where the statistical average of $`\mathrm{\Delta }UU_{Ein}U`$ is calculated by simulation for fixed CM as a function of $`\lambda `$ under an effective potential given by $`\stackrel{~}{U}(\lambda )=(1\lambda )U+\lambda U_{Ein}`$. Note that the center of mass must be calculated in the same manner as described in the paragraph above. Further, note that this expression is only rigorously valid for systems interacting with continuous potentials. In the case of particles with discontinuous potentials, e.g. hard particles, the internal potential energy cannot be turned off continuously. The calculation for this case differs slightly, and is discussed in detail in the original article and in Ref. .

Using Eqns. (12), (18) and (21), we find that the free energy per molecule of the unconstrained crystal is given by

$`{\displaystyle \frac{\beta F}{N_{mol}}}`$ $`=`$ $`\left({\displaystyle \frac{dN}{2N_{mol}}}\right)\mathrm{ln}(2\pi /\beta \alpha ){\displaystyle \frac{1}{N_{mol}}}\mathrm{ln}{\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle \frac{2\pi m_i}{\beta h^2}}\right]^{d/2}`$ (24)   
$`{\displaystyle \frac{\beta }{N_{mol}}}{\displaystyle _0^1}𝑑\lambda \mathrm{\Delta }U_\lambda ^{CM}`$   
$`{\displaystyle \frac{d}{2N_{mol}}}\mathrm{ln}\left({\displaystyle \frac{\alpha \beta }{2\pi _i\mu _i^2}}\right){\displaystyle \frac{\mathrm{ln}(V/N_{mol})}{N_{mol}}}`$   
If we consider the special case of a system of single-atom, identical particles ($`m_i=m`$ and $`N=N_{mol}`$), we obtain the following:

$`{\displaystyle \frac{\beta F}{N}}`$ $`=`$ $`{\displaystyle \frac{d}{2}}\mathrm{ln}\left[{\displaystyle \frac{4\pi ^2m}{\alpha \beta ^2h^2}}\right]{\displaystyle \frac{\beta }{N}}{\displaystyle _0^1}𝑑\lambda \mathrm{\Delta }U_\lambda ^{CM}`$ (26)   
$`{\displaystyle \frac{d}{2N}}\mathrm{ln}(\alpha \beta /2\pi ){\displaystyle \frac{d}{2}}{\displaystyle \frac{\mathrm{ln}N}{N}}+{\displaystyle \frac{\mathrm{ln}\rho }{N}}\text{ ,}`$   
where $`\rho N/V`$. The difference between the present result and the one obtained in ref. is in the fourth term on the right-hand side: $`d\mathrm{ln}N/2N`$. The original article implicitly gave the value +$`\mathrm{ln}N/2N`$ for a 3D crystal. While the difference between the two expressions tends to zero in the limit of large $`N`$, it is non-negligible for system sizes typically employed in the numerical calculations. However, the calculated free energy differences between two solids, such as that between the FCC and HCP hard-sphere crystals, to which the method was applied both in the original article and, more recently, in Ref. , are unaffected by this correction.

In practice, we usually need not calculate the absolute free energy of a crystal, but excess free energy, $`F_{ex}FF_{id}`$, where $`F_{id}`$ is the ideal gas free energy. Let us therefore compute the finite-size corrections to the latter quantity: Given that $`\beta F_{id}/N=\mathrm{ln}[V^N(2\pi m/\beta h^2)^{dN/2}/N!]/N`$, we find

$`{\displaystyle \frac{\beta F_{ex}}{N}}`$ $`=`$ $`{\displaystyle \frac{d}{2}}\mathrm{ln}\left[{\displaystyle \frac{2\pi }{\alpha \beta }}\right]{\displaystyle \frac{\beta }{N}}{\displaystyle _0^1}𝑑\lambda \mathrm{\Delta }U_\lambda ^{CM}`$ (29)   
$`{\displaystyle \frac{d}{2N}}\mathrm{ln}(\alpha \beta /2\pi )+{\displaystyle \frac{\mathrm{ln}\rho }{N}}`$   
$`{\displaystyle \frac{d+1}{2}}{\displaystyle \frac{\mathrm{ln}N}{N}}\mathrm{ln}\rho +1{\displaystyle \frac{\mathrm{ln}2\pi }{2N}},`$   
where we have used $`\mathrm{ln}N!N\mathrm{ln}NN+(\mathrm{ln}2\pi N)/2`$.

Hoover has analyzed the system-size dependence of the entropy of a classical harmonic crystal with periodic boundaries. In this study, it was established that the leading finite-size correction to the free energy per particle of a harmonic crystal is equal to $`\beta ^1\mathrm{ln}N/N`$. If the harmonic approximation is valid, then this implies that the integral in Eqn. (21) should vary as $`+\mathrm{ln}N/N`$ plus higher order correction terms of the order of $`N^1`$, $`N^2`$, etc. Consequently, an inspection of Eqn. (29) suggests that $`\beta F_{ex}/N+(d1)\mathrm{ln}N/(2N)`$ will scale as $`N^1`$, if we neglect terms of order $`𝒪(1/N^2)`$.

To test this prediction we have used the Einstein-crystal method to calculate the absolute Helmholtz free energy of a (three-dimensional) FCC crystal of soft spheres interacting with pair potential of $`u(r)=ϵ(\sigma /r)^{12}`$ for systems of size $`N`$=216, 810, 1728, 5832 and 12096. The sizes were chosen such that the simulation box shape is cubic for each system. Further, the simulations were carried out at $`k_BT/ϵ`$=1.0, $`\rho \sigma ^3`$=1.1964, and employed a coupling constant of $`\alpha \sigma ^2/ϵ`$=66.0. The results are shown in Figure 1 and are clearly consistent with the predictions. The solid line is a linear fit which extrapolates to $`\beta F_{ex}/N`$=1.97298(7) at $`N=\mathrm{}`$ (where the figure between brackets is an estimate for the error in the last digit). Incidentally, we note that, at this density and temperature, the FCC phase of soft-spheres, is more stable than the HCP phase (by an amount $`\mathrm{\Delta }F_{FCCHCP}/(Nk_BT)=0.0028(8)`$). The present results suggest that we are able to correctly account for the leading ($`\mathrm{ln}N/N`$) dependence of the free energy of an arbitrary crystal. In the analytical calculation of free-energy of a harmonic crystal, it is always assumed that the center of mass of the crystal is fixed. Hence, the numerical results presented above do not provide an independent test of the validity of our expression for the contribution to the free energy due to the center-of-mass motion of the crystal.

We can perform a similar analysis for a system of hard spheres ($`\rho \sigma ^3`$=1.0409, $`\lambda _{max}=\alpha /2`$=3000). The results are shown in Figure 2. For hard spheres, $`\beta F_{ex}/N`$ extrapolates to a value of $`5.91889(4)`$ at $`N=\mathrm{}`$, well within the error margin of the original results of Hoover and Ree ($`5.924(15)`$). Note that the slopes of the fits (which are proportional to the $`1/N`$ behavior of the finite size effect) are similar, although not exactly equal. It should be stressed that none of these calculations take into account the existence of defects in the crystal, which, at these levels of precision, is significant. In fact, using the early numerical results by Bennett and Alder  we can estimate the equilibrium vacancy concentration in a hard-sphere crystal at coexistence to be 2.6 $`10^4`$. Such a vacancy concentration has an noticeable effect on the location of the melting point. For instance, the Gibbs free energy per particle at coexistence is lowered by an amount $`\mathrm{\Delta }\mu 3`$ $`10^3kT`$ . This correction is far from negligible, as it is some two orders of magnitude larger than the present numerical accuracy in the absolute free energy. It is likely that vacancies also lower the equilibrium free energy of the soft-sphere crystal. However, for that model, the equilibrium vacancy concentration has, to our knowledge, not been computed.

We would like to thank Stella Consta, Benito Groh, Jonathan Doye and Bela Mulder for stimulating and very useful discussions relating to the material discussed in this work. This work is part of the research program of the “Stichting Fundamenteel Onderzoek der Materie” (FOM) and is supported by NWO (‘Nederlandse Organisatie voor Wetenschappelijk Onderzoek’). JP acknowledges the financial support provided by the Computational Materials Science program of NWO, and by the Natural Sciences and Engineering Research Council of Canada.