# The mold integration method for the calculation of the crystal-fluid interfacial free energy from simulations.

## Abstract

The interfacial free energy between a crystal and a fluid, , is a highly relevant parameter in phenomena such as wetting or crystal nucleation and growth. Due to the difficulty of measuring experimentally, computer simulations are often used to study the crystal-fluid interface. Here, we present a novel simulation methodology for the calculation of . The methodology consists in using a mold composed of potential energy wells to induce the formation of a crystal slab in the fluid at coexistence conditions. This induction is done along a reversible pathway along which the free energy difference between the initial and the final states is obtained by means of thermodynamic integration. The structure of the mold is given by that of the crystal lattice planes, which allows to easily obtain the free energy for different crystal orientations. The method is validated by calculating for previously studied systems, namely the hard spheres and the Lennard-Jones systems. Our results for the latter show that the method is accurate enough to deal the anisotropy of with respect to the crystal orientation. We also calculate for a recently proposed continuous version of the hard sphere potential and obtain the same as for the pure hard sphere system. The method can be implemented both in Monte Carlo and Molecular Dynamics. In fact, we show that it can be easily used in combination with the popular Molecular Dynamics package GROMACS.

## I Introduction

A fluid and a crystal at coexistence are divided by a flat interface. The work needed to create such interface per unit area is known as the interfacial free energy. The crystal-fluid interfacial free energy, , is a highly relevant quantity due to its central role in important phenomena such as wetting or crystal nucleation and growth (1); (2); (3). Despite its importance, is still unknown for many substances given that there is not an easy and reliable way of measuring experimentally (4). This difficulty contrasts with the determination of the fluid-fluid interfacial free energy, a task for which there are well established experimental and computational techniques (5); (6). Unfortunately, these techniques are not easy to implement when one of the phases involved in the coexistence has an infinitely large viscosity, as the crystal phase has. Moreover, is anisotropic and depends on the orientation of the crystal with respect to the fluid. The situation for water, arguably the most important substance on earth, is a good example of the difficulty of measuring : while it is well known that the interfacial tension of liquid water at ambient conditions is 72 mN/m, the reported values for the ice-water interfacial free energy at ambient pressure range from to 25 to 35 mN/m (7).

Computer simulations can be used to assess experimental measurements of and to improve our understanding on the crystal-fluid interface at a molecular scale. An important effort has been devoted to develop simulation methodologies to calculate the crystal-fluid interfacial free energy. To the best of our knowledge, these are the existing computational methods for the calculation of the crystal-fluid interfacial free energy: the cleaving method, the capillary fluctuation method, the methadynamics method, the tethered Monte Carlo method, and the Classical Nucleation Theory method. The cleaving method, proposed by Broughton and Gilmer in 1986 (8), was the first method devised to directly compute in a simulation. In this scheme the reversible work needed to cleave and re-combine the crystal and the fluid is calculated by thermodynamic integration. This method is still in use and an improved version of it has been recently employed to calculate for several water models (9); (10), hard particles (11); (12), Lennard Jones (13) and dipolar fluids (14). The cleaving method has recently been further improved by sorting out some hysteresis issues (15). In the tethered Monte Carlo scheme (16) a complex order parameter is used to allow for a continuous transition between the fluid and the solid. This method has been applied to the hard sphere system (16). The Metadynamics method (17) uses the rare event simulation technique Metadynamics (18) to obtain the work of formation of the interface from a fluid at coexistence. This methodology was originally applied to a Lennard-Jones system (17) and has also been used to assess experimental measurements of for Pb (19). The crystal-fluid interfacial free energy can also be indirectly estimated by combining simulation measurements of the size of critical nuclei with classical nucleation theory (20). This approach has been used, e. g., to estimate for chlatrates (21) or water (22).

All the aforementioned methods have proved successful in the calculation of for a number of systems. However, not all methods are equally good in terms of accuracy, simplicity and computational cost. The anisotropy of for different crystal orientations is easy to study with the capillary fluctuation and with the cleaving methods, whereas dealing with the orientation of the crystal with respect to the fluid is not so trivial for the other methods. On the other hand, while the Metadynamics and the tethered Monte Carlo methods converge well for relatively small system sizes, the capillary fluctuation and the classical nucleation methods require large system sizes. From a practical point of view there are methods simple to implement, like that based on classical nucleation theory, or more cumbersome ones like the cleaving method that requires following a multi-step thermodynamic route. Moreover, all methods but the cleaving require a local-bond order parameter in order to either detect the interface or induce its formation. Such order parameter, used to distinguish liquid-like from solid-like molecules, may be difficult to conceive if the structure of the crystal lattice is complex. In this work we present a simple method for the direct calculation of that gives accurate results even for relatively small system sizes. The calculation of for different crystal orientations is trivial with this methodology. The method can be easily implemented in a bespoke Monte Carlo (MC) code or even in open access Molecular Dynamics (MD) packages as GROMACS (23). In brief, we use a mold of potential energy wells placed at the positions of the atoms in a lattice plane to induce the formation of a crystal slab in the fluid at coexistence conditions. The work of formation of the crystal slab, obtained via thermodynamic integration, is directly related to the interfacial free energy.

We test the method by calculating the interfacial free energy of hard spheres (HS) and Lennard-Jones (LJ) (for several orientations of the crystal) and by comparing the results with values published in the literature (24); (25); (11); (13); (17). Moreover, we compute for the pseudo hard-sphere potential recently proposed in Ref. (26).

## Ii The mold integration method

### ii.1 Description of the method

In this section we describe a new methodology to compute the interfacial free energy between a crystal and a fluid, , by means of computer simulations. The basic idea is to reversibly induce the formation of a thin crystalline slab in the fluid (see Fig. 1 for snapshots of a fluid and a fluid with a crystal slab). The work needed to form such crystalline slab, , is related to . Because the formation of the crystal slab is performed at coexistence conditions the fluid and the fluid plus the crystal slab have the same chemical potential. Then, is just the interfacial free energy times the area of the interface times 2. The factor of 2 is due to the fact that when the crystal slab is formed two crystal-fluid interfaces are created (see Fig. 1, bottom). Thus, can be simply obtained as:

(1) |

In order to induce the formation of the crystal slab we use a mold composed of potential energy wells. The location of the wells is given by that of the particles in the lattice plane whose is calculated. In Fig. 1 we show a snapshot of the mold used for the calculation of for the 100 plane of hard spheres (red spheres). Each potential well must be small enough so that it can only accommodate one particle. When the mold is switched off, particles freely diffuse in the fluid (Fig. 1, top). On the contrary, when the mold is switched on, every well contains a particle and, if the wells are sufficiently narrow, a crystal slab is formed at coexistence conditions (Fig. 1, bottom). Typically, the mold consists of 1 or 2 crystalline planes. In Fig. 1 we show a mold composed of a single plane. When filled with particles, the mold induces the formation of crystalline planes in either side (Fig. 1, bottom) thus giving rise to two crystal-fluid interfaces. By gradually switching the interaction between the mold and the particles the work of formation of the crystal slab at coexistence conditions can be obtained by means of thermodynamic integration.

To perform thermodynamic integration we define the following potential energy:

(2) |

where is the number of particles and the number of wells; denotes the positions of all particles and the position of the wells (which are kept fixed during the simulation); is the potential energy given by the interaction between all particles and is the potential energy given by the interaction between the mold and the particles. is a parameter that varies from 0 to 1 connecting the initial state (mold switched off, Fig. 1, top) with the final state (mold switched on, Fig. 1, bottom). The interaction between the mold and the particles, , is pair additive:

(3) |

where is a square-well interaction between the particle and the well that depends on the distance between their centers, :

(4) |

Where and are the radius and the depth of the wells respectively. These are the only adjustable parameters of the method. Below we explain how to deal with the tuning of these parameters. In any case, can not be larger than the particle radius to avoid multiple filling of a single well.

By performing thermodynamic integration in (27) one can obtain the free energy difference between the fluid and the fluid plus the filled mold, :

(5) | |||||

where and are the coexistence pressure and temperature respectively. The mold coordinates and the edges of the simulation box parallel to the mold are kept fixed throughout the simulation. For that purpose, the pressure is exerted only along the axis perpendicular to the mold, the axis in our case. Thus, the edge is allowed to fluctuate to keep the pressure constant, while the and edges do not change. The mold coordinates are not rescaled when a volume move is performed. The integrand, , is evaluated in simulations for various values of and then integrated numerically to get . is the free energy change due to the appearance of the crystal slab plus that due to the interaction between the particles and the mold. The latter is simply given by (recall that the well-particle interaction is just a square well of depth ). To calculate the interfacial free energy we are just interested in the free energy change due to the generation of the crystal slab:

(6) |

This equation, combined with Eq. 1, allows in principle for the calculation of in a straightforward manner.

There is one open issue, though: the value of , and hence that of , depends on . Therefore, one has to find a priori which value of gives the right value for . We shall refer to this radius as (optimal well radius). In order to chose it is important to understand the way in which the mold affects the free energy landscape. In Fig. 2 we show a sketch of the free energy profile that separates the fluid from the crystal as a function of the crystallinity degree (XD). The latter can be measured, for instance, with the aid of a local bond order parameter that quantifies the number of crystal-like particles in the system (28); (29). The black curve in Fig. 2 corresponds to the free energy profile in the absence of any mold. The liquid and the crystal have the same free energy given that the simulations are carried out at coexistence conditions. In between both phases there is a free energy plateau corresponding to the presence of a crystal slab in the fluid at coexistence conditions. Given that when a crystal slab is present there are two interfaces of the same area (see Fig. 1, bottom), the free energy difference between the plateau and the minima is . For the free energy profile at low crystallinity degrees changes to that sketched by the blue curve. In this case, the free energy gap between the plateau and the fluid’s minimum is reduced by the mold, but there is still a free energy cost to form a crystal slab. Therefore, a fluid where a mold with is switched on can remain stable for a long time before any crystal slab arises. If the minimum given by the blue curve is shallow (values of larger than but close to ) a fluid slab can form after some induction period due to thermal fluctuations. For the free energy profile changes to that schematically shown by the red curve in Fig. 2. Accordingly, as soon as the mold is switched on a crystal slab will quickly develop in order to minimize the free energy. Therefore, the evolution of XD depends on whether is larger or smaller than . Exploiting this difference can be enclosed within a certain range by running simulations for different values of and monitoring the behaviour of XD.

Once is identified, one could in principle perform thermodynamic integration for to obtain . However, in a flat landscape as that given by (green curve in Fig. 2), XD can freely grow and may eventually fall into the crystal’s basin. Therefore, it is not advisable to perform thermodynamic integration for . Instead, it is safe to integrate to states with (blue profile in Fig. 2) for which there is a well defined minimum in the free energy profile. Although these integrations yield underestimates of , they provide a function that can be extrapolated to in order to obtain the right value of .

In summary, the method, which we call mold integration method, consists of the following steps: i) Preparation: The mold coordinates are obtained and an initial configuration of the fluid at coexistence conditions is prepared. The simulation box dimensions must be compatible with those of the mold. ii) Choice of : The mold is switched on and several simulations are run starting from the fluid configuration previously prepared. Simulations are repeated for different values of . By monitoring XD a range within which is enclosed is identified. iii) Calculation of : Thermodynamic integration is performed by gradually switching the mold on. By repeating this for several values of a function is obtained. iv) The extrapolation of to provides the definite value of .

Being the method above described completely novel, our approach shares some features with existing methodologies. For example, in the Metadynamics method (17) the work needed to create a crystal slab in a fluid at coexistence is also calculated, although in a different way (via Metadynamics as opposed to thermodynamic integration) and with the aid of local-bond order parameters that may be difficult to find for complex crystal structures. This difficulty is bypassed in our method with the use of an ad hoc mold In some recent implementations of the cleaving method a cleaving potential based on the location of the particles in the crystal plane is used (9), in resemblance to our mold of potential energy wells. However, whereas in our method the mold is used as a platform for the growth of a crystal slab in the fluid, in the cleaving method the cleaving potential is used to cleave both phases in order to subsequently recombine them. Therefore our route to a system at coexistence is more direct than that proposed in the cleaving framework. The use of potential wells is not exclusive of methods for the calculation of . Potential wells have also been used in the calculation of the free energy of amorphous and crystalline solids via thermodynamic integration (30).

### ii.2 Implementation

The implementation of the mold integration technique in MC is rather straightforward. A routine to evaluate the interaction between the particles and the mold, , via Eq. 3 has to be incorporated to a standard MC code. is evaluated every time a move is attempted and the change of associated to the move is added to the energy change according to which the trial move is accepted or rejected. To perform thermodynamic integration via Eq. 5 the average value of must be evaluated in the course of the simulation.

It is also possible to implement the mold integration technique in MD. We briefly discuss here how to do it for the popular MD package GROMACS (23). The trick is to consider the wells as a special kind of atom. The interaction between the wells and the particles, Eq. 4, can be approximated by the following equation:

(7) |

where , and have the same meaning as in Eq. 4 and controls the steepness of the well’s walls. This potential is continuous and differentiable and can therefore be used in MD. In Fig. 3 we compare given by Eq. 4 (black) with that given by Eq. 7 (red) for . It is evident that a square-well interaction is well approximated by Eq. 7. In GROMACS it is possible to define the well-particle interaction given by Eq. 7 in a tabular form, so there is no need to modify the source code to program the interaction between the wells and the particles. The interaction between different wells has to be also defined in GROMACS in a tabular form. Such interaction is simply 0. In order to fix the position of the wells we use the ‘frozen’ GROMACS option. To perform thermodynamic integration via Eq. 5 we need to be able evaluate for a given value of . To do that we run the simulation with a well-particle interaction given by . Since GROMACS provides average values of the potential energy for any kind of pair interaction, one can obtain as the average particle-well potential energy, and, in turn, . Finally, GROMACS also allows that the pressure is exerted only in one specific direction of the simulation box. Therefore, GROMACS includes all required tools for an easy implementation of the mold integration method in MD.

## Iii Results and discussion

### iii.1 A worked example: of hard spheres

#### Preparation

The first step is obtaining the mold coordinates and a fluid configuration at coexistence conditions. The dimensions of the simulation box must be consistent with those of the mold. The mold coordinates are obtained by replicating the unit cell and taking a plane of the resulting lattice. In our case we replicate 1x7x7 times the fcc unit cell and take a plane of atoms parallel to the plane. The resulting coordinates are shown by the red spheres in Fig. 1, top. Under periodic boundary conditions the mold is just a 100 plane of an fcc lattice. A configuration of the fluid at coexistence conditions (pressure = 11.54 (31)) is then prepared in a box whose and edges have the same length as the and sides of the mold. To achieve this we equilibrate the fluid in an simulation where pressure is exerted only along the axis (we refer to this as ensemble). In this way the length of the axis, , is allowed to fluctuate, while and are kept fixed to the desired value (7 times the unit cell side in this particular example). The resulting simulation box, that contains 1960 particles, is shown in Fig. 1, top, alongside the corresponding mold. We summarize the system size used for the study of the HS system in the top row of Table 1.

System | ST | (x)/() | ||||||
---|---|---|---|---|---|---|---|---|

HS | MC | 100 | 10.978x10.978 | 98 | 1 | 1960 | 0.315 | 0.586 (8) |

PHS | MD | 100 | 12.531x12.531 | 256 | 2 | 5632 | 0.375 | 0.588 (8) |

#### Choice of

Once the fluid is equilibrated we proceed to run simulations starting from a fluid configuration. The mold is switched on at the beginning of the simulations. If the interaction between the mold and the particles is sufficiently large all wells are quickly filled when the mold is switched on. We find this to be the case when in Eq. 4 is larger than . We monitor XD in the course of our simulations. As a measure of XD we use the following parameter, :

(8) |

where is the actual density of the system and and are the coexistence densities of the fluid and the solid respectively. Thus, fluctuates around 0 when the whole system is fluid, and around 1 when the whole system is crystalline. As a crystal slab grows in the fluid should take intermediate values between the typical ones for the fluid and the crystal. In the appendix A we show that this simple way of quantifying XD is totally equivalent to a more sophisticated one based in counting the the number of particles in the largest cluster of solid-like particles.

In Fig. 4 we show the evolution of for several values of . For a given value of we run 10 trajectories starting from the same initial configuration in order to have a statistical picture of the behaviour of the system upon switching the mold on. The trajectories differ in the seed for the random number generator. Each MC simulation consists of a million sweeps. A sweep, in turn, consists of a displacement attempt per particle plus a volume move. The displacement shifts for volume and displacement moves are tuned so that an average acceptance of 30-40 per cent is attained.

Three different types of behaviour can be seen when the trajectories are inspected for each : (a) Behaviour consistent with the presence of a deep minimum in the free energy-XD profile: in plots e) and f) of Fig. 4 XD stays low and fluctuates around a certain equilibrium value for all trajectories. This is consistent with the situation sketched by the blue curve in Fig. 2: is larger than and XD fluctuates around the minimum given by the blue curve. (b) Behaviour not consistent with the presence of a minimum the free energy-XD profile: in plots a) and b) of Fig. 4 XD readily grows as the mold is switched on and each trajectory evolves differently from the others. This corresponds to the situation illustrated by the red curve in Fig. 2: is smaller than and, due to the presence of the mold, there is no free energy penalty for the growth of XD. (c) Behaviour consistent with the presence of a shallow minimum in the free energy-XD profile: In plots c) and d) of Fig. 4 XD fluctuates around an equilibrium value for some trajectories although, stochastically, there are trajectories that visit high values of XD. For instance, for the trajectory given by the orange curve stays at low XD until it jumps around MC sweeps. This phenomenology suggests that there is a minimum in the free energy profile as indicated in Fig. 2 by the blue curve. However, the gap between the minimum and the horizontal plateau is not high and can be stochastically overcome by thermal activation.

Since the optimal radius, , must be in between the highest that shows no hint of a minimum () and the lowest that does show it () we take .

In summary, the recipe to find is to look for a value of that is comprised in between the largest one that shows no indication of the presence of a minimum in the free energy-XD profile and the smallest one that does show it. A given shows no indication of a minimum if XD can grow and evolves in a different way for different trajectories. By contrast, a given shows indication of a minimum if some trajectories show that XD fluctuates around a low constant value. In Appendix B we show how the free energy profile along the XD coordinate can be estimated for each well radius using the information contained in Fig. 4. This is quite helpful to identify which radii generate a free energy profile with a minimum and which ones do not.

#### Calculation of

Once we get a value for we proceed to calculate for by means of thermodynamic integration (Eq. 5). Thermodynamic integration is performed by gradually switching the mold on in such way that all wells are filled with particles when the upper integration limit is reached (i.e. when in Eq. 5). In Fig. 5 (a) we plot the average number of filled wells versus the parameter that controls the strength of the interaction between the particles and the mold via Eq. 2. Each point in Fig. 5 is obtained in an MC simulation consisting of equilibration sweeps and production sweeps. The plot in Fig. 5 corresponds to a well-particle interaction parameter (Eq. 4) and to an . The value of must guarantee that every well contains a particle for . Provided that this condition is fulfilled, can take any value. However, it is convenient that is not too large so that the integrand varies smoothly as increases. In the particular case study we present here the mold is conformed by 98 wells (see Fig. 1, top) As shown in Fig. 5 (a) all 98 wells are filled when approaches 1. On average, about wells are occupied when the mold is switched off. The curve that is actually integrated in Eq. 5 is shown in Fig. 5 (b). The integrand, , is simply given by the product between the average number of filled wells and . The integral of the curve shown in Fig. 5 (b) is , which gives the free energy difference between the system with the mold on and the system with the mold off. To simply get the free energy difference between the fluid and the fluid having the structure induced by the mold, , we need to subtract to the interaction between the mold and the fluid: . divided by two times the area gives an (under)estimate of the interfacial free energy (Eq. 1): . In this step is evaluated for some other values of in order to extrapolate to in the following step.

As previously discussed, in order for thermodynamic integration to be reversible we must avoid integrating at values of that entail any risk that the system crystallizes. Clearly, Fig. 4 shows that such risk is negligible for and , since XD fluctuates around a low, equilibrium value for all trajectories. The situation is not so clear for , where the trajectory given by the orange curve appears to have jumped to the free energy plateau from where the system could evolve towards the crystalline state. Therefore, by performing thermodynamic integration at there is a small chance that the system crystallizes in the typical simulation time required to perform thermodynamic integration. Hence, according to the study shown in Fig. 4, it is safe to perform thermodynamic integration only for . However, one can also try doing thermodynamic integration for ’s closer to and validate the integration a posteriori by checking that the system did not crystallize for any integration point. One of these checks is shown in Fig. 6, where we plot XD for the runs used to compute each integration point in Fig. 5. XD stays low for all integration points, which guarantees that the integration is reversible. It is important to do this check after performing thermodynamic integration, specially for ’s close to .

#### Extrapolation of

Above we have discussed in detail the calculation of for . The same calculation has to be repeated for other values of in order to get a function that can be extrapolated to the radius previously identified as the optimal one: . In Fig. 7 we show for the HS system. The dependency of on looks rather linear, which allows to easily extrapolate to . The extrapolation is given by the open symbol with the error bar in Fig. 7. Thus, our estimated value for for the 100 plane of HS is . The main error source in our calculation comes from the uncertainty in determining . The uncertainty in the thermodynamic integration also contributes, although to a lesser extent, to our final error bar. Our value is in very good agreement with the most recent estimate of (11) (horizontal dashed lines in Fig. 7), obtained via the cleaving method (8). Our value is also in agreement with the latest estimate of for the 100 plane of HS from capillary wave fluctuations: (24).

This excellent result proves the ability of our mold integration method to evaluate the crystal-fluid interfacial free energy. The method is simple conceptually and easy to implement. With a 10 processors machine and our bespoke MC algorithm the calculation of the crystal-fluid interfacial free energy for the 100 plane of HS took us about two days. To further validate the methodology in the following section we show results for the LJ system, for which we compute the interfacial free energy not only for the 100 plane but also for the 111 plane.

### iii.2 for the LJ system

In this section we report the calculation of for the LJ model as modified by Broughton and Gilmer (32). We perform the calculation at the triple point of the model, which was determined in Ref. (33) to be at , corresponding to a pressure ( and are the interaction parameters of the LJ system (32)). At these thermodynamic conditions the density of the fluid is 0.828 and that of the crystal 0.945 (13).

To illustrate the suitability of our methodology to deal with the anisotropy of the crystal-fluid interfacial free energy we calculate for two different crystal orientations. The orientation of the crystal with respect to the fluid is indicated by the Miller indices of the plane parallel to the interface. In this work we calculate for the 100 and the 111 orientations. Obtaining for a given crystal orientation with the mold integration method just requires using a mold coming from the lattice plane that defines such orientation. In Fig. 8 we show the molds used for the calculation of the 100 (left) and the 111 (right) interfacial free energies.

In the previous section, where we describe the calculation of for the HS system, we use a single layer mold. However, the mold can be composed of more than a single layer. We prove in this section that using molds composed of two layers (bilayer mold) one obtains results which are consistent with those obtained by using a single layer (monolayer mold). Moreover, we also compare in this section MC with MD. In the implementation section above we describe the way the mold integration method can be easily implemented in the popular MD simulation package GROMACS, which we use to calculate for the LJ system.

The first step is to prepare the mold and a fluid configuration at equilibrium in a simulation box compatible with the mold. We give details on how to do this for the HS system in the previous section. In Table 2 we list the simulations used for the calculation of for the 100 and 111 crystal orientations for the LJ system. We indicate the area of the simulation box side parallel to the mold (x), the number of wells that conform the mold, the number of layers of the mold and the number of particles of the system.

Once the initial set up is ready we run a number of trajectories (about 10) for different values of in order to find . As previously described for the HS system, can be found by looking at the behaviour of XD as a function of the simulation time for different ’s. The values of thus obtained are shown in Table 2. It is interesting to realize that changes from one crystal orientation to another. is always smaller for the 111 plane. This shows that it is necessary to adjust the value of for every crystal orientation separately. We also show in Table 2 that also depends on the number of layers that conform the mold. Monolayer molds need smaller values of to induce the formation of a crystal slab than bilayer ones. The work required to fill a well increases as its radius decreases since the smaller the well’s volume the more unfavourable it becomes confining a particle inside. Therefore, one has to supply more energy per well to a monolayer mold in order to get the same energy per unit area as in a bilayer mold. This seems reasonable since a bilayer mold has twice as many wells per unit area. Finally, from our analysis it also turns out that MD s are slightly larger than MC ones (when compared for the same crystal orientation and the same number of layers in the mold). This may be due to the fact that, although the continuous potential given by Eq. 7 closely follows a square-well interaction (see Fig. 3), the equivalence is not perfect. Nevertheless, as we show below, the small difference in between MC and MD is not reflected in the estimated value of .

ST | (x)/() | ||||||

MC | 100 | 11.323x11.323 | 98 | 1 | 1960 | 0.305 | 0.372(8) |

MD | 100 | 11.323x11.323 | 98 | 1 | 1960 | 0.315 | 0.372(8) |

MD | 100 | 14.543x14.543 | 324 | 2 | 6480 | 0.385 | 0.373(8) |

MC | 111 | 13.726x9.906 | 120 | 1 | 2160 | 0.285 | 0.350(8) |

MD | 111 | 13.726x9.906 | 120 | 1 | 2160 | 0.295 | 0.354(8) |

MD | 111 | 13.726x9.906 | 240 | 2 | 2160 | 0.385 | 0.348(8) |

100 | 12158 | 0.371(3) (13) | |||||

100 | 38740 | 0.369(8) (34) | |||||

100 | 2352 | 0.370(2) (17) | |||||

100 | 1790 | 0.34(2) (35) | |||||

111 | 8916 | 0.347(3) (13) | |||||

111 | 38740 | 0.355(8) (34) | |||||

111 | 1674 | 0.35(2) (35) |

Once is identified for each system the next step is to perform thermodynamic integration for at least a couple of values of in order to obtain a function that can be extrapolated to . In Fig. 9 we show for all systems investigated. Red symbols correspond to the results for the 100 orientation and black ones to the 111 orientation. Let us start by discussing the results for the 100 orientation. Filled symbols correspond to the calculation of via thermodynamic integration and empty ones to the extrapolation of to . The black squares correspond to MC and the black circles to MD simulations, both with a monolayer mold. It is clear from Fig. 9 that both MC and MD yield consistent results for the calculation of via thermodynamic integration. A linear extrapolation of the MC and MD data to their corresponding values of (see Table 2) provides an estimate for , indicated by the open symbols in Fig. 9 and reported in Table 2. Within the error of the method both MC and MD give the same for the 100 orientation. By using a bilayer mold (black diamonds) we get also, within error, the same value for . Red squares and circles correspond to the results for a 111 monolayer mold as obtained from MC and MD respectively. Again, a good agreement between both simulation techniques is obtained. Moreover, the results for the bilayer (red diamonds) give the same as the monolayer mold. In summary, in Fig. 9 we show that the mold integration technique gives consistent results regardless the simulation technique (MC or MD), or the number of layers in the mold (1 or 2). The accuracy of the technique is sufficient to distinguish between the interfacial free energy of two different crystal orientations (100 and 111).

As discussed above, in order to compute we recommend to obtain first two or three under-estimates of for , where thermodynamic integration is reversible, and then extrapolate the results to . A close inspection of Fig. 9 shows that the MC estimate of for the 100 orientation was directly performed at . This allows to directly estimate without the need of any extrapolation, but apparently contradicts the advice of performing thermodynamic integration for . In fact, we had to resort to a tailored type of MC move in order to perform thermodynamic integration for , where the reversibility of thermodynamic integration is compromised by the possibility that the system fully crystallizes. Such move consisted in performing blocks of thousands of MC sweeps that are accepted or rejected according to the final value of XD. If XD increases beyond the point at which the system is committed to fully crystallize, the whole block of MC sweeps is rejected and re-started with a different seed for the random number generator. Otherwise, the MC simulation continues normally. In order to set both the length of the simulation blocks and the XD threshold it is necessary to get some experience first by examining several unbiased runs at . In this way we have two MC estimates of for the 100 interface: one coming from the ‘direct’ calculation of at (as described in this paragraph), and another coming from the extrapolation of estimates for . As shown in Fig. 9 both estimates coincide pretty well, which gives us confidence in the extrapolation procedure described in previous sections. Although both ways of estimating are equally valid, we recommend the use of the extrapolation method because it is more general as it does not require the implementation of the MC-block moves described in this paragraph.

In Table 2 we show that one can obtain consistent results for molds with one or two layers. In principle any number of layers can be used. However, one must take into account that increases as the number of layers increases (see table 2) and that can not be larger than 0.5 in order to avoid multiple filling of the wells. For this reason, in practise, we could not use a mold with more than two layers to compute . In any case, there is no practical advantage in using bi-layer over mono-layer molds. In fact, with a mono-layer mold the number of well-particle interactions is half as many and the code runs faster.

The interfacial free energy for the LJ model at the triple point has been directly determined by Broughton and Gilmer in 1986 using the cleaving method (35), by Laird and Davidchack using a more accurate variant of the same methodology (13), by Morris and Song using a capillary fluctuation approach (34) and by Angioletti-Uberti et al. using a Metadynamics-based approach (17). In the bottom part of table 2 we report the values of obtained in these works. The agreement between our data and those obtained in Refs. (13); (34) is very good. The comparison with Ref. (35) is not bad either, particularly taking into account that in 1986 the computational resources did not allow for accurate enough calculations to distinguish between different crystal orientations. In Ref. (20) was indirectly estimated via classical nucleation theory (). The discrepancy with our data may be partly due to the fact that in Ref. (20) a value of averaged over several crystal orientations is provided. In summary, in Table 2 we show that the values obtained from our method are in good agreement with the general consensus reached for the of LJ for two different crystal orientations.

The extensive work done on the LJ system allows for a comparison between different simulation methods. In table 2 we report the number of molecules used for the calculation of for each simulation method. The number we report for Ref. (34), where the capillary fluctuation method was used, is an average over all systems employed. In terms of computational cost, the capillary fluctuation method is expensive. Moreover it requires the evaluation of the spectrum of capillary waves for at least three different crystal orientations in order to provide a value of . The method based on Classical Nucleation Theory also requires a large number of particles because such theory works best for large cluster sizes (20). The cleaving method used in Ref. (35) used relatively small systems, but the results were not accurate enough to distinguish between different crystal orientations. Many more particles were used in Ref. (13) in order to gain accuracy. Both the Metadynamics method and our mold integration method are capable of producing accurate results for systems of less than 2000 particles. We have checked that there are no significant system size effects present in our simulations. In Table 2 we show that a calculation with 6480 particles gives the same result as one with 1960. Therefore, the possibility of using small systems is a positive aspect of our method.

Another advantage is the simplicity with which it deals with different crystal orientations (one simply has to use a mold coming from the corresponding lattice plane). Also the cleaving and the capillary fluctuation methods easily deal with the anisotropy of . In both methods the fluid is brought into contact with the crystal at the desired orientation. For the capillary fluctuation method the problem is not so straightforward, though. What is obtained from the analysis of the capillary waves spectrum is the interfacial stiffness, and several orientations must be combined with a cubic harmonic expansion in order to obtain estimates of . The way to deal with the anisotropy of is more complex for the other methods. In the Metadynamics method, for instance, an order parameter has to be devised in order to induce the growth of the crystal. Finding a good order parameter may be a non-trivial task for crystal structures whose complexity goes beyond that of simple fcc or bcc lattices (36). The classical nucleation method also requires an order parameter to measure the size of the embedded clusters.

### iii.3 for the pseudo hard-sphere potential

A system composed of hard spheres (HS) is arguably the simplest non-trivial model having fluid, crystal and glass phases (37); (38). Therefore, this model is widely used by researchers on a diverse range of problems like the glass transition (39), crystal nucleation (40), or granular matter (41). There is great interest in finding a continuous potential whose kinetic and thermodynamic behaviour reproduces that of the discontinuous potential of pure HS. Finding such potential would allow to explore the physics of the HS system with simple MD simulations. This is important because MD simulation packages like GROMACS are nowadays accessible to a large scientific community. For instance, having a continuous version of the HS potential would be of great help for the investigation the crystallization of hard spheres, where large discrepancies between experimental and simulation measures of the nucleation rate have been reported (see Ref. (42) and references therein).

Quite recently, Jover et al. have proposed a continuous potential that at reduced temperature 1.5 behaves very much alike a system of pure HS in terms of the equation of state and the diffusion coefficient (26). We refer to the potential proposed by Jover et al. at reduced temperature 1.5 as the pseudo hard-sphere potential, PHS. Later on, Espinosa et al. showed that the coexistence pressure and densities for the PHS model are also very similar to those of the pure HS system (43).

The PHS potential is a good candidate for the investigation of the crystal nucleation rate of HS given that, as mentioned before, both models have a very similar thermodynamic coexistence, diffusion coefficient, and equation of state. Nothing is known, however, about the of the PHS potential, a crucial parameter in crystal nucleation (3). If was also similar to that of pure HS then the PHS model could be reliably used in MD simulations to obtain predictions about the crystallization behaviour of HS. Here, we use the mold integration method to obtain for the PHS model.

We evaluate the interfacial free energy for the 100 crystal orientation in order to compare with our results for the pure HS model. We use a two-layer mold with 128 particles in each layer. The system size is summarized in the bottom row of table 1. All simulations are performed with the GROMACS MD package.

To determine we monitor the XD as a function of time for several trajectories and for several ’s. To monitor XD we use the parameter defined in Eq. 8. The results are shown in Fig. 10. As explained above for the HS and LJ cases will be comprised in between the largest value of that shows no indication of the presence of a minimum in the free energy-XD profile and the smallest that does show it. In Fig. 10 clearly corresponds to the presence of a deep minimum (one that can not be overcome by thermal activation for any of the 10 trajectories). By contrast, 0.36 and 0.37 show no presence of a minimum because XD can grow and evolves in a different way for different trajectories. correspond to a shallow minimum because in all cases there are trajectories for which XD fluctuates for a significant period around typical values for the deep minimum case (XD). Therefore we set .

Once we get we perform thermodynamic integration for values of and obtain the solid points shown in Fig. 11. The extrapolation of these data to gives that is, within the error bar, the same value we find for the pure HS system. This result implies that the PHS model can be used with confidence for the study of the behaviour of HS. The simulation details and results for the HS and PHS models are compared in Table 1.

## Iv Summary and conclusions

We propose a novel simulation methodology for the calculation of the crystal-fluid interfacial free energy. The main idea of the method is the use of a mold of potential energy wells to induce the formation of a crystal slab in a fluid at coexistence conditions. The coordinates of the mold’s wells are given by the lattice positions of the crystal plane whose interfacial free energy is evaluated. The interaction between the wells and the fluid particles is square-well like. The free energy difference between the fluid and the fluid having the crystal slab induced by the mold is obtained by means of thermodynamic integration along a reversible path in which the wells are gradually filled. The method consists in four basic steps: (i) preparation of the mold and of the initial configuration of the fluid; (ii) estimation of the optimal radius of the wells; (iii) calculation of the interfacial free energy as a function of the well radius by thermodynamic integration; (iv) extrapolation of the function obtained in step (iii) to the optimal radius to get the final estimate of the interfacial free energy.

We validate our methodology by calculating the interfacial free energy of systems composed of hard spheres and Lennard-Jones particles. In both cases we find a very good agreement with previous estimates. Moreover, we show that our methodology is accurate enough to discriminate between different crystal orientations of the Lennard-Jones system. We also use the new method to calculate the interfacial free energy for a continuous version of the hard sphere model for which, to the best of our knowledge, the interfacial free energy had not been previously calculated. Within the statistical uncertainty of our calculations we obtain the same interfacial free energy for both the continuous and the discontinuous potentials.

One of the main advantanges of our method with respect to existing ones is that no local-bond order parameter is required to either detect the interface or induce the growth of the solid (44); (17); (20); (16). The cleaving method does not require an order parameter either (35), but it entails following a rather cumbersome thermodynamic route. Our method, gives accurate results even for relatively small systems (about 2000 particles), which can not be achieved with the capillary fluctuation (44) or the classical nucleation (20) methods. Moreover, it can potentially deal with complex crystal lattices with no extra methodological complexity. The method can be easily implemented either in Monte Carlo or in standard Molecular Dynamics packages such as GROMACS. Therefore, we hope it will be appealing to the scientific community interested in investigating the properties of the crystal-melt interface.

Acknowledgements

All authors thank the Spanish Ministry of Economy and Competitiviness for the financial support
through the project FIS2013-43209-P.
E. Sanz and J. R. Espinosa acknowledge financial support from the EU grant 322326-COSAAC-FP7-PEOPLE-2012-CIG
and E. Sanz from a Spanish grant Ramon y Cajal.

## Appendix A Measure of XD

Here we show that the simple order parameter used to quantify XD is totally equivalent to a more sophisticated one based on counting the the number of particles in the largest cluster of solid-like particles. Solid-like particles can be identified by means of a local-bond order parameter based on the local coordination of the particles (28). The specific parameters we use in this work are those given in Ref. (45). The largest cluster of solid-like particles corresponds to the crystal slab induced by the mold at coexistence conditions (e. g., the 6-7 crystal planes that can be seen around the mold in Fig. 1, bottom). In Fig. 12 we show the equivalent to Fig. 4 but using the number of particles in the biggest cluster of solid-like particles, , instead of as the parameter to follow the formation of the crystal. By comparing both figures it is clear that the simple parameter provides the same information as the more sophisticated .

## Appendix B Free energy profiles.

By analysing the trajectories where XD is monitored (e.g. Figs. 4 or 12) the free energy profile along the XD coordinate can be estimated as:

(9) |

where is the probability that the system takes a certain value of the order parameter, XD. Then, by simply making a histogram of XD for all trajectories performed for a given well radius it is possible to get an estimate of the corresponding free energy profile. In Fig. 13 we show the free energy profile thus calculated for the 111 plane of the LJ system using both and as measures for XD. The conclusions that can be drawn by examining either order parameter are the same: for there is a minimum whereas for there is not. Consequently, we set the optimal radius for this system to .

Eq. 9 does not give absolute free energies (there is a missing constant) but allows to determine whether there is a minimum present or not. In any case, in order to compare all curves in the same free energy scale we have shifted each minimum to the work needed to fill the wells for the corresponding well radius (calculated by thermodynamic integration via Eq. 6). For the cases where the minimum is absent () we have shifted the plateau of the curves to the work needed to fill a mold of wells with , given by the dashed horizontal line. The statistics of the free energy given by Eq. 9 are reasonably good when the system repeatedly samples configurations in the vicinity of a free energy minimum but become poor when it quickly moves along the free energy plateau. Therefore, one has to be cautious and restrict the use of Eq. 9 to small values of XD.

With the degree of accuracy we got in the present study we were able to distinguish the anisotropy between the 111 and 100 faces of the LJ system. One could in principle try to further improve the accuracy by decreasing the range within which is enclosed (by launching trajectories for more values of ). This is certainly a possibility worth exploring. However, it may require a substantial amount of trajectories to detect a mininum shallower than 1, which is the depth of the shallowest minimum we could probe (curve corresponding to in Fig. 13).

### References

- W. Boettinger, S. Coriell, A. Greer, A. Karma, W. Kurz, M. Rappaz, and R. Trivedi, “Solidification microstructures: recent developments, future directions,” Acta Materialia, vol. 48, no. 1, pp. 43 – 70, 2000.
- D. P. Woodruff, The solid-liquid interface. Cambridge: Cambridge University Press, 1973.
- K. F. Kelton, Crystal Nucleation in Liquids and Glasses. Boston: Academic, 1991.
- K. F. Kelton, Solid State Physics. Dordrecht: Academic Press, 1991.
- J. S. Rowlinson and B. Widsom, Molecular theory of capillarity. Clarendon Press, Oxford, 1982.
- G. Gloor, G. Jackson, F. Blas, and E. de Miguel, “Test-area simulation method for the direct determination of the interfacial tension of systems with continuous or discontinuous potentials,” J. Chem. Phys., vol. 123, p. 134703, 2005.
- H. R. Pruppacher, “A new look at homogeneous ice nucleation in supercooled water drops,” J. Atmosph. Sci., vol. 52, p. 1924, 1995.
- J. Q. Broughton and G. H. Gilmer, “Molecular dynamics investigation of the crystal–fluid interface. vi. e xcess surface free energies of crystal–liquid systems,” J. Chem. Phys., vol. 84, no. 10, pp. 5759–5768, 1986.
- R. Handel, R. L. Davidchack, J. Anwar, and A. Brukhno, “Direct calculation of solid-liquid interfacial free energy for molecular systems: Tip4p ice-water interface,” Phys. Rev. Lett., vol. 100, p. 036104, Jan 2008.
- R. L. Davidchack, R. Handel, J. Anwar, and A. V. Brukhno, “Ice ih-water interfacial free energy of simple water models with full electrostatic interactions,” Journal of Chemical Theory and Computation, vol. 8, no. 7, pp. 2383–2390, 2012.
- R. L. Davidchack, “Hard spheres revisited: Accurate calculation of the solid–liquid interfacial free energy,” J. Chem. Phys., vol. 133, no. 23, p. 234701, 2010.
- Y. Mu and X. Song, “Crystal-melt interfacial free energies of hard-dumbbell systems,” Phys. Rev. E, vol. 74, p. 031611, Sep 2006.
- R. L. Davidchack and B. B. Laird, “Direct calculation of the crystal–melt interfacial free energies for continuous potentials: Application to the lennard-jones system,” J. Chem. Phys., vol. 118, no. 16, pp. 7651–7657, 2003.
- J. Wang, P. A. Apte, J. R. Morris, and X. C. Zeng, “Freezing point and solid-liquid interfacial free energy of stockmayer dipolar fluids: A molecular dynamics simulation study,” The Journal of Chemical Physics, vol. 139, no. 11, p. 114705, 2013.
- R. Benjamin and J. Horbach, “Crystal-liquid interfacial free energy via thermodynamic integration,” The Journal of Chemical Physics, vol. 141, no. 4, p. 044715, 2014.
- L. A. Fernandez, V. Martin-Mayor, B. Seoane, and P. Verrocchio, “Equilibrium fluid-solid coexistence of hard spheres,” Phys. Rev. Lett., vol. 108, p. 165701, Apr 2012.
- S. Angioletti-Uberti, M. Ceriotti, P. D. Lee, and M. W. Finnis, “Solid-liquid interface free energy through metadynamics simulations,” Phys. Rev. B, vol. 81, p. 125416, Mar 2010.
- A. Laio and M. Parrinello, “Escaping free-energy minima,” Proc. Natl. Acad. Sci., vol. 99, p. 12562, 2002.
- S. Angioletti-Uberti, “The solidâliquid interface free-energy of pb: comparison of theory and experiments,” Journal of Physics: Condensed Matter, vol. 23, no. 43, p. 435008, 2011.
- X.-M. Bai and M. Li, “Calculation of solid-liquid interfacial free energy: A classical nuclea tion theory based approach,” J. Chem. Phys., vol. 124, no. 12, p. 124707, 2006.
- B. C. Knott, V. Molinero, M. F. Doherty, and B. Peters, “Homogeneous nucleation of methane hydrates: Unrealistic under realistic conditions,” J. Am. Chem. Soc., vol. 134, pp. 19544–19547, 2012.
- E. Sanz, C. Vega, J. R. Espinosa, R. Caballero-Bernal, J. L. F. Abascal, and C. Valeriani, “Homogeneous ice nucleation at moderate supercooling from molecular simulation,” Journal of the American Chemical Society, vol. 135, no. 40, pp. 15008–15017, 2013.
- B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, “Algorithms for highly efficient, load-balanced, and scalable molecular simulation,” J. Chem. Theory Comput., vol. 4, pp. 435–447, 2008.
- R. L. Davidchack, J. R. Morris, and B. B. Laird, “The anisotropic hard-sphere crystal-melt interfacial free energy from fluctuations,” J. Chem. Phys., vol. 125, p. 094710, 2006.
- R. L. Davidchack, J. R. Morris, and B. B. Laird, “The anisotropic hard-sphere crystal-melt interfacial free energy from fluctuations,” The Journal of Chemical Physics, vol. 125, no. 9, p. 094710, 2006.
- J. Jover, A. J. Haslam, A. Galindo, G. Jackson, and E. A. Muller, “Pseudo hard-sphere potential for use in continuous molecular-dynamics simulation of spherical and chain molecules,” The Journal of Chemical Physics, vol. 137, no. 14, p. 144505, 2012.
- D. Frenkel and B. Smit, Understanding Molecular Simulation, ch. 7.1. Academic Press, 2002.
- P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, “Numerical calculation of the rate of crystal nucleation in a Lennard-Jones system at moderate undercooling,” J. Chem. Phys., vol. 104, p. 9932, 1996.
- W. Lechner and C. Dellago, “Accurate determination of crystal structures based on averaged local bond order parameters,” The Journal of Chemical Physics, vol. 129, no. 11, p. 114707, 2008.
- T. Schilling and F. Schmid, “Computing absolute free energies of disordered structures by molecular simulation,” The Journal of Chemical Physics, vol. 131, no. 23, p. 231102, 2009.
- E. G. Noya, C. Vega, and E. de Miguel, “Determination of the melting point of hard spheres from direct coexistence simulation methods,” The Journal of Chemical Physics, vol. 128, no. 15, p. 154507, 2008.
- J. Q. Broughton and G. H. Gilmer, “Molecular dynamics investigation of the crystalâfluid interface. i. bulk properties,” The Journal of Chemical Physics, vol. 79, no. 10, pp. 5095–5104, 1983.
- J. Q. Broughton and G. H. Gilmer, “Molecular dynamics investigation of the crystalâfluid interface. iv. free energies of crystalâvapor systems,” The Journal of Chemical Physics, vol. 84, no. 10, pp. 5741–5748, 1986.
- J. R. Morris and X. Song, “The anisotropic free energy of the lennard-jones crystal-melt interface,” J. Chem. Phys., vol. 119, no. 7, pp. 3920–3925, 2003.
- J. Q. Broughton and G. H. Gimer, “Molecular dynamics investigation of the crystal-fluid interface. vi. excess surface free energies of crystal-liquid systems,” The Journal of Chemical Physics, vol. 84, p. 5759, 1986.
- E. E. Santiso and B. L. Trout, “A general set of order parameters for molecular crystals,” J. Chem. Phys., vol. 134, no. 6, p. 064109, 2011.
- B. J. Alder and T. E. Wainwright, “Phase transition for a hard sphere system,” J. Chem. Phys., vol. 27, no. 5, pp. 1208–1209, 1957.
- P. N. Pusey and W. van Megen, “Phase behaviour of concentrated suspensions of nearly hard colloidal spheres,” Nature, vol. 320, p. 340, 1986.
- K. N. Pham, A. M. Puertas, J. Bergenholtz, S. U. Egelhaaf, A. Moussaid, P. N. Pusey, A. B. Schofield, M. E. Cates, M. Fuchs, and W. C. K. Poon, “Multiple glassy states in a simple model system,” Science, vol. 296, p. 104, 2002.
- S. Auer and D. Frenkel, “Prediction of absolute crystal-nucleation rate in hard-sphere colloids,” Nature, vol. 409, p. 1020, 2001.
- C. Song, P. Wang, and H. A. Makse, “A phase diagram for jammed matter,” Nature, vol. 453, no. 7195, pp. 629–632, 2008.
- L. Filion, R. Ni, D. Frenkel, and M. Dijkstra, “Simulation of nucleation in almost hard-sphere colloids: The discrepancy between experiment and simulation persists,” The Journal of Chemical Physics, vol. 134, no. 13, p. 134901, 2011.
- J. R. Espinosa, E. Sanz, C. Valeriani, and C. Vega, “On fluid-solid direct coexistence simulations: The pseudo-hard sphere model,” The Journal of Chemical Physics, vol. 139, no. 14, p. 144502, 2013.
- J. J. Hoyt, M. Asta, and A. Karma, “Method for computing the anisotropy of the solid-liquid interfacial free energy,” Phys. Rev. Lett., vol. 86, pp. 5530–5533, Jun 2001.
- P. N. Pusey, E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon, and M. E. Cates, “Hard spheres: crystallization and glass formation,” Phyl. Trans. Roy. Soc. A, vol. 367, pp. 4993–5011, 2009.