Solid-liquid interface free energy through metadynamics simulations

The solid-liquid interface free energy ${\gamma }_{sl}$ is a key parameter controlling nucleation and growth during solidification and other phenomena. There are intrinsic difficulties in obtaining accurate experimental values, and the previous approaches to compute ${\gamma }_{sl}$ with atomistic simulations are computationally demanding. We present an approach which is to obtain ${\gamma }_{sl}$ from a free-energy map of the phase transition reconstructed by metadynamics. We apply this to the benchmark case of a Lennard-Jones potential, and the results confirm the most reliable data obtained previously. We demonstrate several advantages of our approach: it is simple to implement, robust and free of hysteresis problems, it allows a rigorous and unbiased estimate of the statistical uncertainty, and it returns a good estimate of the thermodynamic limit with system sizes of a just a few hundred atoms. It is therefore attractive for applications which require more realistic and specific models of interatomic forces.

Figure 1

(Color online) Schematic representation of the flattening of the effective FES by the repulsive bias of Eq. (4), during a metadynamics simulation in a one-dimensional collective-variable space. We show the underlying FES $G\left(s\right)$ and the bias accumulated at chosen times (arbitrary units). For a sufficiently long simulation, $G\left(s\right)+V\left(s,t\right)\to \text{constant}$ so that one can obtain an accurate estimate of the free-energy surface simply by taking the negative of the bias.

Figure 2

(Color online) Angular function $c_{\alpha }\left(\hat{\mathbf{x}}\right)$ as defined in Eq. (8). The function is shown as a polar plot, centered on an fcc lattice. $c_{\alpha }$ has well-defined peaks in the directions of the nearest neighbors.

Figure 3

(Color online) Cutoff function used to define the regions $A$ and $B$ for the calculation of the two collective variables. The function varies smoothly from 0 to 1 so as to avoid discontinuities when atoms transit between the two regions.

Figure 4

(Color online) Time-averaged value of the order parameter $\overline{\varphi }$ for an individual atom in the bulk solid and liquid phases, as evaluated for the LJ system at different temperatures and across the solid-liquid transition. The bounding lines delimit the range one standard deviation above and below the mean value. Dashed lines correspond to the values of the order parameter for metastable solid and liquid. Note that the parameters $r_1$ and $r_0$ for radial cutoff (9) are scaled from the values used at $T_m$ according to the changes of the equilibrium density.

Figure 5

(Color online) Two-dimensional (2D) FES reconstructed by well-tempered metadynamics, together with selected snapshots of configurations obtained during the simulation. Atoms participating in $s_A$ are colored in orange, those in $s_B$ in blue, and the region of CV space corresponding to each snapshot is marked. The negative peaks in the FES clearly correspond to the two single-phase regions. They are separated by a very wide plateau, corresponding to the presence of well-defined interfaces between solid and liquid phases at various different positions relative to the $A/B$ partition [insets (a), (b), and (e)]. In inset (f) we report the projection of the FES along the single CV $s=\left(s_A+s_B\right)/2$.

Figure 6

(Color online) Radial pair-correlation function $g\left(r\right)$ for the liquid in the presence of an interface during our metadynamics simulations (red, dashed) and for a normal molecular-dynamics simulation of a bulk liquid (blue). It is clear that the two curves are very similar, thus ruling out quantitatively the presence of metastable structures during our simulation. In the inset, the kinetic-energy distribution during our simulation is plotted in comparison to that expected for a canonical ensemble at $T=T_m$. Again the two are very close demonstrating that the metadynamic bias does not induce any systematic deviation from the correct ensemble, and that quasiequilibrium conditions hold.

Figure 7

(Color online) A snapshot of the solid-liquid interface taken from the final part of a well-tempered metadynamics simulation. The scaled order parameter $\overline{\varphi }$ has been used to color atoms. The atoms with a liquidlike configuration, with $\overline{\varphi }<0.45$ have been hidden, the atoms with $\overline{\varphi }>0.65$ have been colored in blue. Finally, atoms in intermediate configurations, with $0.45<\overline{\varphi }<0.65$ have been made translucent. It is clear that—whatever threshold is used to ascertain the solid from the liquid state—the interface is not flat on the atomic scale.

Figure 8

(Color online) Convergence of the FES [inset (a)] and its average error [inset (b)] with respect to time for our one-dimensional (1D) simulations. Ten simulations have been performed on a $7\times 7\times 12$ supercell, with the single-CV setup described in the text. The FESs in inset (a) are constructed by averaging equal-times biases of the independent runs, and the error bars correspond to the standard deviation. In inset (b) we plot the relative error, averaged between ${\overline{\varphi }}_l$ and ${\overline{\varphi }}_s$, as a function of simulation time. The error is plotted on a log-log scale and the least-squares linear fit shows that the angular coefficient is close to the theoretical value of $-1/2$ predicted for a simple Langevin model.

Figure 9

(Color online) (a) The plateau region corresponding to the presence of the solid/liquid interface, with different relative amounts of the two phases, is drawn for different supercell sizes in the $z$ direction. The curves are shifted for display purpose, and they are only meant to demonstrate how the flat portion of $G\left(s\right)$ becomes more extended as larger supercells are considered. Finite-size effects are present also in the region for $s<0.1$. In inset (b) we compare results for simulations with different in-plane sizes, showing a further finite-size effect which depends on the change in $T_m$ and causes a residual slope of $G\left(s_A\right)/A$ even after the complete formation of an interface.