Obtaining the solid-liquid interfacial free energy via multi-scheme thermodynamic integration: Ag-ethylene glycol interfaces

The solid-liquid interfacial free energy $γsl$ is an important quantity in wetting, nucleation, and crystal growth.^1–4 For example, in solution-phase, shape-controlled metal nanocrystal syntheses, thermodynamic control over the product morphology essentially depends on the values of $γsl$ for various crystal facets under the framework of the Wulff construction. Accurate and facet-specific values of $γsl$ are needed to correctly predict the Wulff shapes of crystals. Interfacial free energies are also used to calculate the nucleation free-energy barrier and predict the nucleation rate in classical nucleation theory.^3–5 

Although the interfacial free energy is a critical control parameter in crystal growth and many other phenomena, it is still difficult to directly measure this quantity in experiments. Experimentally, the interfacial free energy or interfacial tension can be indirectly estimated from the nucleation rate,^4,5 the equilibrium crystal growth rate,⁶ and the contact angle between a fluid droplet and a solid surface.^7,8 However, these methods are always associated with large errors and the incapability to resolve differences among various crystal planes with reasonable accuracy.^2,4–6 Many theoretical and simulation methods have consequently been developed to calculate the interfacial free energy $γ$ and interfacial tension $σ$⁠, in lieu of experimental means.^9–23 In this paper, we focus on estimating these quantities via atomic-scale simulations.

Kirkwood and Buff⁹ first showed the dependence of surface tension on pressure anisotropy via a mechanical route, but this relation is limited to fluids, planar surfaces, and continuous potentials. Some studies extended their method to solid-liquid interfaces.^10–13 However, unlike fluid-fluid interfaces, where $γ$ and $σ$ are identical in value, the interfacial tension $σ$ is not the same as the interfacial free energy $γ$ for solid-liquid interfaces, so that $γsl$ cannot be indirectly estimated using $σsl$⁠.^24,25 The interfacial tension can be obtained from the interfacial free energy via the relation

where A is the interfacial area. The value of $σ$ is only identical to $γ$ when changing A (i.e., creating a strain) does not change the interfacial free energy per area. For fluids, $dγ/dA$ remains invariant to changes in A because the interface is preserved when stress or strain is applied.²⁴ Solid interfaces are often sensitive to strain upon compressing or stretching, resulting in non-zero $dγ/dA$ such that $σ≠γ$⁠. Therefore, methods for calculating $σsl$ are not suitable for estimating $γsl$⁠.

Efforts from various groups have been dedicated to developing systematic methods to calculate $γsl$ from molecular dynamics^13–22 (MD) or Monte Carlo^22,23 (MC) simulations. Many of these methods^{13–19,22,23} are designed for crystal melts (i.e., a crystal coexistent with its liquid phase at the melting temperature), such as the interfacial fluctuation method¹⁴ and a series of “cleaving-wall” methods.^15–17 These methods have been successfully applied to idealized systems, such as Lennard-Jones particles^{14,16,17,20,21} or hard spheres.^22,26 Some can potentially be extended to more practical systems with proper modification.^21,27–29 For example, the “phantom-wall” method by Leroy and Müller-Plathe^20,21 can, in principle, compute the Gibbs solid-liquid interfacial free energy for any solid-liquid system. However, this method requires meticulous tuning of its parameters to avoid numerical instability. In this paper, we present a multi-scheme thermodynamic integration (TI) method that is designed to obtain $γsl$ from MD simulations based on complex force fields for realistic systems. MD simulations based on our method can be easily carried out using open-source simulation packages, such as LAMMPS.³⁰ 

Our method is inspired by a series of TI-based methods in which $γsl$ is obtained from canonical MD simulations. The idea behind these methods originated from the work of Broughton and Gilmer,¹⁵ who noted that $γsl$ can be calculated from the excess free energy, or the reversible work, required to form an interface of unit area between a solid and liquid from the individual bulk phases. They use a “cleaving potential” to gradually split the liquid and solid bulk phases into slabs. They subsequently join the crystal and liquid slabs together and remove the “cleaving potential” to create the interface. The value of $γsl$ is the total free-energy change of the entire process divided by the interfacial area. This method involves designing precise cleaving-potential functions for each facet. The uncertainty in the final results, caused primarily by hysteresis, can be too large to distinguish facet-specific $γsl$⁠, so that this method would fail to predict Wulff shapes.

There have been improvements to Broughton and Gilmer’s method in recent years. Davidchack and Laird¹⁶ modified their approach by replacing the cleaving potential with a planar wall, in the “cleaving-wall” method. The use of a physical wall simplifies the interaction between the system of interest and the cleaving potential and better resolves facet-specific $γsl$⁠, ameliorating another issue with the Broughton-Gilmer approach. Recently, Benjamin and Horbach¹⁷ pointed out that hysteresis still exists in some steps in the cleaving-wall method, which makes the process irreversible and, thus, leads to an inaccurate estimate of the free energy. They demonstrated that their improved TI scheme can successfully remove hysteresis compared to previous TI approaches. In addition, the system-wall interaction is further simplified to a repulsive Gaussian potential in their method.

The multi-scheme TI method that we will describe below is based on the Benjamin-Horbach approach. It addresses two additional features for generalized implementation. All the cleaving-wall and similar methods described above involve manipulating the periodic boundary conditions (PBC) of the simulation cell in the direction orthogonal to the interface. In previous work, the methods for manipulating the PBC were implemented in specialized simulation codes and it would require code modifications to implement these methods in open-source simulation packages. To generalize the TI procedure, we incorporate separate schemes for the solid and liquid phases to circumvent these issues. In addition, we introduce a numerical approximation to obtain the thermodynamic integrand when it is difficult to obtain analytical differentiation of the potential energy. These features make our method more user-friendly and applicable to a broad range of systems.

To demonstrate how this method is applied, we use canonical MD simulations to model two possible interfaces between solid Ag and liquid-phase ethylene glycol (EG). These interfaces are relevant in the solution-phase growth of Ag nanocrystals.^31,32 We consider the Ag(100)-EG and Ag(111)-EG interfaces and we obtain $γAg(100)-EG$ and $γAg(111)-EG$ from our simulations. This paper is arranged as follows. In Sec. II A, we describe the general procedure in the multi-scheme TI method. A complete description of force fields used in the Ag-EG system is included in Sec. II B, followed by stepwise mathematical TI expressions for obtaining $γsl$ in Sec. II C. The MD simulation details can be found in Sec. III. Results are presented and discussed in Sec. IV and conclusions are given in Sec. V.

The solid-liquid interfacial free energy is the free-energy change when creating a unit area interface between a solid and a liquid from the individual bulk phases. One path to create such an interface is to cleave simulation cells holding bulk liquid and solid phases into slabs separately and bring them into contact. This approach has been implemented in prior cleaving-potential/cleaving-wall studies.^15–17 However, implementing such a pathway in MD simulations of non-ideal systems can be difficult. For example, it can be difficult to introduce a wall that avoids cutting through chemical bonds (which would cause numerical instability) in molecular systems. Alternatively, we can obtain $−γsl$ by following the reverse path — going from a cell containing a solid-liquid interface to the individual bulk phases. Doing so, we can successfully circumvent possible simulation difficulties, such as the large repulsion that results when the cleaving-wall cuts through liquid molecules.

To create individual bulk solid and liquid phases, we use the physical cleaving-wall approach, similar to that in the work of Davidchack and Laird¹⁶ and by Benjamin and Horbach.¹⁷ We follow the idea of Benjamin and Horbach to include a short-range repulsive Gaussian wall, because it is straightforward to implement and has a minimal influence on the system of interest. By using a physical wall, we can also precisely control the volumes of the solid and liquid phases, so that they remain at the desired bulk densities throughout the simulations.

A step-wise description of our method is given in Fig. 1. The multi-scheme TI procedure proceeds as follows:

• Prep step.

  To initiate the multi-scheme TI process, we equilibrate a simulation cell containing a solid-liquid interface. As suggested in previous studies,^16,17 the solid and liquid regions are designed to have similar volumes. The solid-liquid interface is in the xy-plane with dimensions $Lx≈Ly$ and $Lz≈2Lx≈2Ly$⁠. For crystals, the values of L_x and L_y can be slightly different, depending on the in-plane atom arrangement. All dimensions need to be large enough to exceed twice the force-field cutoff distance. The initial simulation cell has PBC in all directions. In setting up the initial solid slab, it should be noted that an additional solid slab is required to assist in the recovery of the PBC to achieve a bulk solid in step 5 in Fig. 1. The original solid slab should be designed so that it forms a continuous crystal phase when the additional slab is inserted.

• Step 1.

  In the equilibrated solid-liquid cell, we gradually switch on two repulsive Gaussian walls at each of the two solid-liquid interfaces — a total of four walls. Each wall only interacts with either the solid or liquid — the walls interacting with the solid are shown in black in Fig. 1 and the walls interacting with the liquid are shown in green. The repulsion can be different for the solid and liquid, as long as the walls effectively prevent the translation of molecules across them. The positions of the walls should be carefully chosen following two considerations: First, the wall should disturb the system as little as possible, and second, because the wall positions are fixed throughout the entire process, the solid volume confined by the wall should be the same as the bulk solid volume to aid in the recovery of PBC later in the process. The walls can be either perfectly flat or made up of densely packed fictitious atoms.

• Step 2.

  After the repulsive Gaussian walls are fully turned on, we gradually switch off the solid-liquid interaction. This step can be subject to hysteresis if the interaction is switched off too quickly. This is because solid surfaces can induce weak ordering of liquid molecules in the near-surface region and this ordering disappears as the solid-liquid interaction diminishes. At the end of this process, the interacting solid-liquid interfaces vanish and we have non-interacting solid and liquid slabs whose positions and volumes are fixed by the walls.

• Step 2b.

  We separate the two non-interacting slabs into two individual simulation cells. Note that while both simulation boxes are still fully periodic in all dimensions, there are vacuum gaps in the z dimension that render the systems effectively non-periodic. The z-dimension periodicity is fully recovered in later steps, to create bulk solid and liquid.

• Step 3.
  After removing the solid, the liquid slab is in contact with a vacuum gap that has the volume of the original solid. To transform this to bulk liquid, the adjacent liquid interfaces are brought into contact by decreasing the length L_z of the vacuum gap. The liquid phase is then homogeneous and in full PBC in all directions. The free-energy difference is obtained by integrating the long-range force f_p from periodic images along the pathway L using

  $ΔF=∫fpdL.$

  (2)

  We use Fig. 2 to demonstrate how Eq. (2) is implemented. Since forces are vectors, we divide the liquid slab equally into two halves (subsections 1 and 2 in Fig. 2) and their periodic images (subsections 1^′ and 2^′ in Fig. 2) to obtain a positive long-range force in the positive direction and a negative long-range force along the negative direction. Initially, liquid subsection 1 and periodic-image subsection 2^′ are L₀ apart. Along the pathway in which subsections 1 and 2^′ have each moved L₀/2 toward each other, the free energy in Eq. (2) is given by

  $ΔF3=∫0L02fpdL+∫0−L02(−fp)dL.$

  (3)

  Here f_p represents the long-range force experienced by subsection 1 from 2^′ and −f_p represents the long-range force experienced by subsection 2^′ from 1. At the end of step 3, the liquid is bulk-like except it does not have translational freedom across the two overlapping walls.
• Step 4.

  We gradually remove the two overlapping walls inside the liquid cell by decreasing their interaction with the liquid to zero. As we proceed, liquid molecules gradually gain translational freedom in the z direction across the wall. After the wall is removed in this step, the liquid is completely in its bulk phase.

  The procedure we use for the liquid phase is not applicable to solids. This is because, when two solid periodic images are close, the strong attraction between the periodic cells can be large enough to cause distortion of the solid atoms and it can inhibit the vibrational motion of solid atoms at elevated temperatures, which is subject to undesirable hysteresis. In solids and liquids without long-range interactions, we can recover the PBC using the procedure in step 5.

• Step 5.

  In this step, the periodicity of the solid is recovered along the z dimension. To execute this step, we put the solid slab and its bound walls adjacent to a frozen bulk solid with the same L_x and L_y (shown in pink in Fig. 1). We freeze the additional bulk solid to maintain its bulk atomistic structure throughout the process. It is important to match the interfacial structure between the solid slab and the frozen bulk solid if the solid is a crystal. We gradually turn on the interaction between the frozen solid and the solid slab of interest until the interaction reaches full strength. In this way, the frozen bulk acts like an external field that helps to transform a solid surface, which was originally in contact with vacuum, to a part of the bulk. At the end of this process, we recover PBC in the z dimension and the solid slab is bulk-like.

• Step 6.

  We immediately remove the additional frozen bulk solid and one of the two walls, and we gradually remove the remaining wall inside the solid of interest. The immediate removal of the frozen region does not disturb the solid of interest because all of these atoms are bulk-like after PBC are recovered at the end of step 5. Thus, the solid at one side of the cell boundary in the z direction readily interacts with the solid at the other side of the cell boundary after the frozen part is removed.

  We retain one wall in this process instead of two to avoid double counting the solid-wall interaction. This is because in steps 1, 2, and 5, each surface of the solid slab interacts with only one of the two walls. If both walls are retained after removing the frozen bulk solid, the two walls at the cell boundaries will overlap and the solid near the walls would experience twice the repulsion as in previous steps. As we gradually remove the remaining solid wall, the center-of-mass motion of the solid needs to be removed to avoid the solid crossing the wall when the repulsion is weak, otherwise hysteresis may occur. A bulk solid is obtained at the end of step 6.

In all steps except for step 3, a parameter $λ$⁠, which varies between 0 and 1, is coupled to the relevant intermolecular interactions to achieve a smooth transformation from the initial state to the final state. The free-energy difference $ΔF$ of each step can be obtained via TI using

where $U(λ)$ is the coupled potential energy at $λ$⁠.

The sum of all free-energy changes in steps 1-6 is the negative of the free energy to create a solid-liquid interface with an area of A = L_xL_y. The interfacial free energy per unit area is given by

in which the factor of 2 in the denominator corresponds to the two interfaces in a simulation cell with PBC.

To better demonstrate how the multi-scheme TI method is implemented to obtain the Ag-EG interfacial free energy, we describe the force fields employed in MD simulations before going into mathematical details in Sec. II C.

Our MD simulations are based on a set of force fields that have high fidelity to both experiment and first-principles density functional theory (DFT) calculations.^33–39 It has been shown in other studies that this set of force fields can correctly capture the interatomic interactions in the Ag-EG system.^37,40,41 We use the embedded atom method (EAM) potential to describe Ag–Ag interactions. The EAM potential for Ag used in this work is in tabulated form provided by Williams et al.³⁴  We use the CHARMM force field³³ for liquid-phase EG interactions. Long-range Coulomb interactions are calculated via the particle-particle/particle-mesh Ewald summation in k-space.

The interactions between Ag and EG are characterized by the Metal-Organic Many-Body (MOMB) force field,^37,39 which consists of three potential models fitted by Zhou et al.,³⁷ namely, the Morse potential for short-range Pauli repulsion and short-range direct bonding, the van der Waals (vdW) potential implemented in Grimme’s DFT-D2 method,^42,43 and an EAM-like many-body potential which supplies one-way electronic density from the O atoms on EG to Ag. The validity and accuracy of these force fields are discussed in detail elsewhere.^37–39 The mathematical expression of the Morse potential is given by

where D₀ and $α$ are parameters that control the depth and width of the potential minimum, r₀ is the distance corresponding to the energy minimum, and r_ij is the distance between atoms i and j. Grimme’s expression for vdW interactions has the form

where values of the global scaling s₆ factor, d, the species-dependent dispersion coefficient C₆, and vdW radii R₀ can be found in Zhou et al.³⁷  The EAM-like one-way electron density supply from O to Ag is given by

where $ρO→Ag$ is the electron density from O to Ag, $ρAg−Ag$ is the electron density among Ag atoms and f_O (=0.55 here) is a scaling factor. The EAM potential for Ag–Ag interactions and the EAM-like potential for O–Ag interactions are integrated into a single Finnis and Sinclair EAM (EAM-fs) potential form for utilization.⁴⁴ 

The system-wall interaction is described by a general pairwise-repulsive Gaussian potential that is symmetrical at the wall, given by

where r_ik is the distance between atom i and wall atom k, and parameters H and $σh$ describe the height and width of the Gaussian peak, respectively. The wall is composed of fictitious atoms arranged as a two-dimensional square lattice with a lattice constant of 0.486 $Å$⁠. The value of maximum repulsion H is set to be $25kBT$ for atoms in EG and $30kBT$ for Ag. This combination is sufficient to prohibit molecules from moving across the wall.

Ideally, free energies are calculated by integrating analytical derivatives $⟨∂U(λ)∂λ⟩$ over $λ$⁠. However, analytical derivatives for complex $λ$-coupling, such as coupling $λ$ to $ρ$ in the EAM-fs potential, can be involved. When analytical derivatives are not available, we can obtain $⟨∂U(λ)∂λ⟩$ numerically by perturbing configurations sampled at $λ$ with interactions scaled by $λ±δ$⁠, and retrieve $⟨∂U(λ)∂λ⟩$ using the central difference, or forward and backward difference for $λ=0$ and $λ=1$⁠, respectively. For example, the central-difference numerical approximation of $⟨∂U(λ)∂λ⟩$ is given by

where $𝐑λ$ is the configuration sampled at $λ$⁠. The perturbing parameter $δ$ should be small compared to $λ$⁠. In our Ag-EG simulations, we chose to implement a linear coupling of $λ$ to the relevant interactions. As discussed elsewhere, the shortcoming of using linear coupling is insufficient precision at the beginning or the end of the TI process.⁴⁵ We resolve this problem by including more $λ$ points at necessary places. Detailed mathematical expressions for all steps are given below.

In step 1, two sets of overlapping walls for the solid (W_S) and liquid (W_L) are switched on and $λ$ is linearly coupled to the wall potential as

where U₁ is the $λ$-coupled potential in step 1, $ϕGauss(r)$ is the unscaled Gaussian potential in Eq. (9), $ϕGauss(λ,r)$ is the scaled Gaussian potential, i represents system atoms, and j represents wall atoms. Equation (11) implies that $λ$ is coupled to H in Eq. (9). The thermodynamic integrand at each $λ$ can be calculated as

It is important to note that the configurations ${𝐑λ}$ in Eq. (12) are sampled at each $λ$ value, not at the unscaled $λ=1$⁠. The value of $U1(𝐑λ)$ can be obtained by mapping the unscaled Gaussian potential onto configurations sampled at $λ$-scaled points. The free-energy difference in step 1 is obtained via

This integration goes from 0 to 1 as the walls are switched on.

In step 2, the Ag-EG interactions are switched off by coupling $λ$ to the Morse potential, the vdW potential, and the electron density $ρO→Ag$⁠. In this step, $λ$ is coupled to D₀ in the Morse potential in Eq. (6), C₆ in the vdW potential in Eq. (7), and $ρO→Ag$ in the EAM-like potential in Eq. (8). The scaled potential energy is given by

in which $φAg(∑ρO→Ag(r))$ is the embedding function that describes the one-way electron density supply from O to Ag, and the sums run over i Ag atoms, j atoms in EG, and k O atoms in EG. Analytical derivatives of potential energy for the first two terms in Eq. (14) can be obtained as described in step 1 because of the explicit coupling of $λ$⁠. Numerical treatment is needed for $φAg(∑ρO→Ag(r))$ due to the complexity of $φAg$⁠. The thermodynamic integrand is integrated over $λ$ from 1 to 0 as we turn off the Ag-EG interactions. The free-energy difference in step 2 can therefore be obtained via

In step 3, we transform the repeating liquid slabs and vacuum gaps to a bulk-like liquid, and the free-energy change is obtained by integrating the long-range force along the path using Eq. (3).

Step 4 is similar to step 1 except that $λ$ goes from 1 to 0 as we gradually remove the wall. The free-energy difference $ΔF4$ is calculated as

where $U4(λ,𝐑λ)$ is the total Gaussian potential measured at each $λ$ point.

In step 5, we recover the PBC of the Ag slab assisted by an additional frozen bulk Ag crystal (Ag_FZ). At the beginning of step 5, when the interaction between Ag and Ag_FZ is 0, the frozen bulk Ag resembles vacuum to the Ag slab, serving as a smooth transient state from the end of step 2. We design the unscaled Ag_FZ →Ag interaction to be the same as the Ag–Ag interaction, so that at the end of this process, the frozen bulk Ag resembles Ag next to the Ag slab of interest.

The EAM description of Ag–Ag interactions contains two parts: the embedding function and pairwise interaction. As step 5 proceeds, we gradually increase both the electron density supply from Ag_FZ to Ag and the pairwise interaction between Ag and Ag_FZ. It is not necessary to supply electron density from Ag to Ag_FZ because the additional frozen bulk Ag is treated as an external field. The scaled potential energy in step 5 is given by

where $ρAgFZ→Ag$ is the electron density from Ag_FZ to Ag, $VAg−AgFZ$ is the pairwise interaction, i represents the Ag atoms in the slab and j represents the Ag_FZ atoms. The pairwise interaction in Eq. (17) can be differentiated analytically in a similar way to the Gaussian, Morse, and vdW potentials, but the derivative of energy resulting from the embedding function is obtained numerically for the same reason as in step 2. The free-energy difference of step 5 has the form

in which U_5,em is the energy from the embedding part and U_5,pair is from the pairwise part in Eq. (17).

Step 6 is carried out in a similar manner as step 4, and the free-energy change $ΔF6$ is given by

where $U6(λ,𝐑λ)$ is the total Gaussian potential measured at each $λ$ point.

We perform all our simulations using a locally modified LAMMPS package that contains the function for Grimme’s vdW potential. The modified parts of LAMMPS for this work are to include the vdW interaction in the force field, which is tailored to solve a specific problem in metal nanocrystal synthesis. The choice of force fields does not affect the implementation of this method in any version of LAMMPS. All systems are simulated in the canonical ensemble, so we obtain Helmholtz free energies. Usually, the strain, stress, and energy in crystals, especially in metals, are sensitive to the lattice constant — simulating at lattice constants other than the correct value can result in inaccurate energies in NVT simulations. Since the lattice constant varies with temperature, we first perform simulations to determine the bulk Ag lattice constant at 413 K to estimate the surface free energy of Ag with high accuracy. The lattice constant corresponds to the value that gives the desired bulk pressure, which is 4.1246 $Å$ for Ag at 413 K (Fig. 3). This value is used to determine L_x and L_y in all configurations.

We generate two initial configurations, Ag(100)-EG and Ag(111)-EG, and equilibrate them at a temperature of 413 K and a pressure of 1 bar to adjust the density of the liquid. The EG-Ag(100) simulation box contains an $18×18×24$ (atoms) Ag(100) slab and the EG-Ag(111) simulation box contains an $18×20×21$ (atoms) Ag(111) slab. These supercell dimensions yield the same L_x and L_y for Ag(100), a similar L_x and L_y for Ag(111), and similar thicknesses for the two slabs of the different facets. The rest of the box is filled with EG molecules. We design the two systems to have the same solution volume and number of EG molecules, so that the calculations to acquire $γsl$ can share the same $ΔF3$ and $ΔF4$ shown in Fig. 1.

We equilibrate both systems in the NPT ensemble using the Nosé-Hoover^46,47 thermostat and barostat for 5 ns. We have justified the validity of this thermostat on small and stiff materials⁴⁸ like Ag in the supplementary material. To ensure the correct lattice constant, we fix L_x and L_y and only apply pressure coupling to the z-dimension. At the end of the NPT equilibration, the liquids in both simulation boxes are at the correct density. If the volumes of the liquid regions have subtle discrepancies, we slightly adjust the box volume so that two liquid regions have exactly the same volume for the reason stated above. Systems are then further equilibrated in the NVT ensemble using the Nosé-Hoover thermostat for an additional 5 ns. To avoid hysteresis, the system should be equilibrated extensively at each $λ$ point and the production run needs to be long enough to provide reliable ensemble average. To do so, at each $λ$⁠, we first equilibrate the system in the NVT ensemble using the Berendsen⁴⁹ thermostat for 3 ns, followed by 5 ns equilibration using the Nosé-Hoover thermostat. The production run is 5 ns in the NVT ensemble using the Nosé-Hoover thermostat. The number of $λ$ points varies with steps; details are given in Sec. IV. A time step of 1 fs is used in all simulations. Throughout the TI process, we fix the positions of the walls by excluding the wall atoms from position and velocity update, as well as temperature control. The center-of-mass motion of the system is removed at all times, otherwise hysteresis may result if it shifts in the z-direction when the wall repulsion is not sufficient.

To proceed with step 1, we carefully insert two walls composed of fictitious atoms, W_S and W_L, at each of the two solid-liquid interfaces, and two walls superimpose with each other. Since wall positions are fixed throughout the entire procedure, the choice of the wall position should consider the volume of bulk Ag, that is, the volume between two W_S should have the bulk Ag volume. The insertion of walls should also disturb the system as little as possible. At interfaces, there is a small gap between the solid and liquid molecules and this gap is about 2 $Å$ for Ag and EG on average (Fig. 4). With our chosen parameters, the Gaussian potentials of both W_S and W_L damp to 0 near r = 0.6 $Å$⁠, which is less than half of the interface gap. The results in Sec. IV also indicate that our choice of wall position based on the consideration of bulk Ag volume only marginally disturbs the system.

When we turn off the solid-liquid interaction in step 2, we set up simulations of energy perturbation around each $λ$ to numerically differentiate the embedding function of the Ag–O interaction. At each $λ$⁠, we generate a series of configurations, and rerun the same configurations using potentials scaled with $λ±δ$⁠, in which $δ$ is at least two orders of magnitude less than $λ$⁠. We apply separate thermostats to the solid and liquid to ensure correct temperature control as the solid-liquid interaction decreases.

To evaluate Eq. (3) in step 3, we need to extract the long-range force f_p (cf. Fig. 2) from all other forces within the same simulation box. To get f_p, we ran simulations at each distance L and output the total force f_p,L on each subsection, as well as the configuration of EG and the walls, at each output time step. We then re-ran the same set of configurations of EG and the walls but with a distance of L₀ between the two periodic images to obtain the non-periodic force $fp,LL0$ for the two subsections. The long-range force f_p is calculated by subtracting the non-periodic force from the total force on a subsection, given by

In step 4, we gradually switch off the liquid-wall interaction. Since in step 3, the two walls that bound the liquid are brought to overlap, we keep both walls in step 4 for a smooth transition. Simulations of step 4 are carried out in a reverse manner as in step 1.

Prior to step 5, an additional frozen bulk Ag with L_x and L_y dimensions exactly the same as the Ag slab is prepared. We use an eight-layer bulk oriented in the {100} direction for the Ag(100) slab and nine-layer bulk oriented in the {111} direction for the Ag(111) slab. The number of bulk layers should take the periodicity of the facet into consideration so that the surface of the slab aligns perfectly with the adjacent additional frozen bulk. To numerically obtain the derivative of the embedding energy $⟨∂U5,em∂λ⟩$⁠, we perform similar simulations to those in step 2.

Simulations in step 6 are carried out in a similar manner as in step 4, except that only one wall is kept in the simulation cell. The center-of-mass motion of the Ag is removed at all times during this process.

The facet-specific interfacial free energies $γAg(100)−EG$ and $γAg(111)−EG$ are obtained using Eq. (5). We use the same values, of $ΔF3$ and $ΔF4$ for both systems because their liquid components are of the same material and volume. Figures 5(a)–5(f) show the corresponding thermodynamic integrands per unit area in steps 1–6. Table I summarizes the free-energy, potential-energy, and entropy change per unit area in each step. The potential-energy change in each step is obtained using

where $⟨Ui,final⟩$ and $⟨Ui,initial⟩$ are the average potential energies of the system at the final and initial $λ$ points in step i. The change in entropy in each step is calculated using the following thermodynamic relation:

with a T = 413 K for all steps. Note that these steps proceed from solid-liquid interfaces to individual bulk phases, so that the total free-energy change is the negative of the interfacial free energy. Thermodynamic integrands are integrated numerically using the trapezoidal rule. The error at each $λ$ point reflects the standard error of the mean, and the errors in $ΔF$⁠, $ΔU$⁠, and $ΔS$ are calculated based on error propagation.

In general, the free-energy contributions from switching on/off the wall potential (steps 1, 4, and 6) to the total free-energy change are less significant than those from other steps, and are negligible in steps 1 and 6. The largest free-energy contribution is from step 5, which corresponds to making bulk Ag from two Ag interfaces, followed by steps 2 and 3.

In step 1, $ΔF1$ for Ag(100) is marginally higher than Ag(111) because the gap between the solid and liquid is slightly smaller in the EG-Ag(100) system (cf. Fig. 4). The thermodynamic integrands for both facets from $λ=0$ to 1 in Fig. 5(a) and $ΔF1$ in Table I are small compared to other steps. We therefore confirm that employing a short-range Gaussian repulsive potential for the wall-system interaction has a negligible effect on the system. If the repulsion is too large or the Gaussian is too wide, the wall can inhibit the normal vibrations of the Ag atoms, which leads to incorrect estimation of $ΔF$⁠.

The magnitude of $ΔF2$ for Ag(111) is higher than for Ag(100) — a result that is not obvious because EG binds slightly stronger to Ag(100) according to DFT calculations⁵⁰ and MD simulations.³⁷ To understand this result, we analyzed the conformations of the EG molecules on the Ag surfaces. We quantified the alignment of EG molecules along the axis normal to the surface plane. Since EG is a relatively small molecule and its preferred adsorption on Ag(100) is primarily due to stronger Ag–O interactions, we chose the vector connecting the two O atoms on an EG molecule to characterize the conformations of EG on both surfaces. For each surface, we construct a second-order tensor S for EG molecules within a 6 $Å$ distance normal to the Ag surface (i.e., EG molecules in the first major and minor peaks in Fig. 4) and average over all N molecules included using⁵¹ 

in which the sum runs over i EG molecules, $uα(β)i$ is one of the (x,y,z) components of the unit O–O vector uⁱ, and $α$ and $β$ are indices selected from (x,y,z). We calculate the eigenvalues and eigenvectors of this second-order tensor, and find that the O–O vectors are isotropic in the xy-plane but have a distinct alignment perpendicular to the z-axis, with an eigenvector of (−0.005, −0.003, 1.0) for Ag(100) and (0.007, 0.003, 1.0) for Ag(111). The eigenvalues corresponding to the eigenvectors that indicate significant z-axis perpendicular alignment are −0.204 and −0.212 for Ag(100) and Ag(111), respectively. This suggests that EG molecules align more orthogonally to the z-axis (i.e., flatter on the surface) on Ag(111) than on Ag(100), which may result in a tighter network parallel to the surface that causes a larger free-energy change when we turn off the Ag-EG interaction. In a previous study,^40,41 we observed that the deposition flux of Ag atoms is smaller for EG-covered Ag(111) than for Ag(100) surfaces, which is also an evidence for “flatter” alignment of EG on Ag(111).

The entropy changes in step 2 for both systems are also counterintuitive. At an interacting Ag-EG interface, EG forms ordered layers near the surface (cf. Fig. 4) and imposes vibrational constraints on the Ag surface atoms. When the Ag-EG interactions are turned off, the near-surface EG molecules gain translational freedom and the surface Ag atoms gain vibrational freedom, both of which intuitively suggest an increase of entropy, whereas $ΔS2$ in Table I are negative. To understand the entropy loss, we again look at the distribution of EG near the Ag surfaces before and after the Ag-EG interactions are switched off. We plot the distribution of O atoms on EG in the direction orthogonal to the wall on Ag(100) and Ag(111) in both cases (Fig. 6) and compare. The O-atom distribution is more representative than the center-of-mass distribution because of the relatively strong Ag–O bonding. We note that the EG density within 3 $Å$ of the wall is much lower when the Ag-EG interactions are off than when the interactions are on. The lower probability of EG to be sampled at these regions corresponds to a volume decrease, which is associated with an entropy loss,⁵² and when this entropy loss is greater than the translational and vibrational entropy gain, the overall entropy change is negative.

In step 3, we use the equivalence of the free-energy change to reversible work. The value of f_p in Fig. 5(c) reflects the total force on both subsections in Fig. 2, and the y-axis indicates the force along its path. The CHARMM force field used to describe EG includes a Lennard-Jones potential with a cutoff at 12 $Å$ and a long-range Coulomb interaction. Therefore, while the attractive contribution from the Lennard-Jones potential between two subsections is only effective when they are within 6 $Å$ from their relative center (i.e., the black dashed line in Fig. 2), the attractive contribution from the Coulomb interaction exists at all separations. As the two subsections move close to each other, f_p continuously increases, reflecting a growing attraction. However, this attraction peaks around ˜24.5 $Å$⁠, followed by a slight decrease. This is because when the vacuum gap reduces to a certain distance, EG molecules on one side of the wall begin to enter the repulsive region of the other wall. This repulsion keeps increasing as the two subsections move towards each other. At the end, when the walls overlap, EG molecules on each side of the wall feel twice the repulsion from the wall as they experienced initially. According to the design of the Gaussian potential in Eq. (9), EG from one side of the wall can experience the repulsion from the other wall when the two walls are 0.6 $Å$ apart, leading to an attractive force maximum around 24.5 $Å$⁠.

In step 4, we remove the wall from the bulk-like liquid obtained in step 3 and EG molecules gradually gain translational motion across the wall. Simulation trajectories show that the fraction of EG molecules crossing the wall remains very small until $λ$ decreases to 0.2. This is reflected in Fig. 5(d) as nearly constant small $(1A)⟨∂U∂λ⟩$ values in the region between $λ=0.2$ and 1. When the wall repulsion becomes weak, EG molecules can move across the wall or even overlap with wall atoms when $λ=0$⁠, causing large energy derivatives between $λ=0$ and 0.2.

The PBC of the Ag slab in the direction normal to the surface are recovered in step 5. In Fig. 5(e), the values of thermodynamic integrands for Ag(100) are significantly more negative than for Ag(111) at each $λ$⁠, from which we integrate a more negative $ΔF5$ for Ag(100) than Ag(111) (cf. Table I). The free-energy change in this step contributes the most to Ag-EG interfaces and represents the surface free energy of the bare Ag surfaces. We will address the significance of this quantity below.

For step 6, where we remove the wall from the bulk-like Ag, the small thermodynamic integrand values in Fig. 5(f) and the insignificant free-energy change in Table I sufficiently show that the wall does not disturb the system and has a negligible contribution towards the total free energy. The free-energy change in this step is smaller for Ag(111) than Ag(100), which is expected as the inter-layer distance for Ag(111) is larger than Ag(100).

The causes of hysteresis, as well as solutions to prevent it, have been addressed in the work of Benjamin and Horbach.¹⁷ Their six-step approach successfully resolves most hysteresis found in the four-step Davidchack-Laird approach,¹⁶ except in the final step, where the walls are gradually removed from the solid-liquid interface. The cause of hysteresis in this step is the movement of the interface orthogonal to the wall when the wall repulsion is scaled with a small $λ$⁠. Since we follow the suggestions in the Benjamin-Horbach approach, with an additional care on the center-of-mass motion relative to the wall in all steps, the paths presented in this work are free of hysteresis. However, when dealing with molecular adsorption on surfaces, the conformational change of the liquid molecule may lead to another type of hysteresis (i.e., irreversible pathways in step 2), and such hysteresis should be more pronounced for large molecules with high degrees of freedom in arranging themselves near a surface. We address this problem in the supplementary material and show that no hysteresis exists in our calculation. We will discuss more on this type of hysteresis and its solutions for large molecules in future work.

Table II summarizes the interfacial energies and entropies obtained in this work. In addition to the Ag-EG interfacial free energy, we can approximate the surface energies and entropies of bare Ag in contact with vacuum from the sum of the changes in steps 5 and 6. Our results successfully resolve the $γ$ anisotropy of Ag, and show that both the bare-surface free energy and the solid-liquid interfacial free energy of Ag(100), are higher than those of Ag(111). The interaction with liquid EG lowers the surface free energy of both facets to a similar extent, with Ag(111) being lowered slightly more than Ag(100). Moreover, we find in Table II that the excess potential energy $ΔU$ is indicative of the excess free energy $γ$ in which the facet with a higher $ΔU$ also has a higher $γ$⁠. This correspondence can be informative. Considering that it is computationally much less expensive to obtain potential energies than free energies, calculating the potential energy of a surface can provide an indication of the behavior of $γ$ for possible solid-liquid interfaces before going into extensive computation.

The signs and trends in the surface entropies obtained in this work are consistent with theories of the entropy of surface formation.⁵³ These theories indicate that surface formation for both solids and liquids is always associated with an increase in entropy. For solids, more specifically, crystals in this work, a surface atom gains entropy through having more configurations associated with defects and more vibrational freedom associated with the loss of geometrical constraints.^53,54 Although the solid surfaces in our study do not have structures such as steps, kinks, or vacancies, it is possible to extend our method to such surfaces by replacing the flat wall used in this study with a wall that accommodates these features. In the absence of surface structures, a higher entropy increase is expected for Ag(100) than for Ag(111), since the fcc(100) facet loses more nearest neighbors when forming a surface. This is seen in Table II.

The surface free energy of bare Ag in vacuum has been discussed in the literature and estimated using different approaches.^55–59 For bare Ag surfaces, Ag(111) is naturally more stable by possessing a lower surface free energy than Ag(100). A crude way to rationalize this is by counting the number of broken bonds when forming a surface from the bulk — more missing nearest neighbors means higher work to create a surface thus higher surface energy.⁵⁶ Galanakis et al. calculated $γAg$ for various facets based on the broken-bond ratio relative to the 111 plane. $γAg$ has also been obtained at 0 K using DFT,^56,57 the modified analytical EAM (MAEAM) method,⁵⁹ and the EAM potential³⁴ employed in this work. In the work of Vitos et al.,⁵⁷ Wen and Zhang,⁵⁹ and Williams et al.,³⁴ the zero-K surface energy is obtained using

where U_slab is the potential energy of a relaxed surface slab, N is the number of atoms in the slab, U_bulk is the potential energy per atom of the bulk atoms, and A is the surface area in the simulation cell.

Table III lists values of the surface energies of bare Ag(100) and Ag(111) slabs at zero K. Despite the difference in temperature, our multi-scheme TI method agrees well in the trend in literature values that $γAg(100)$ is higher than $γAg(111)$⁠. It should be noted that surface energies of Ag calculated using EAM potentials are lower than those calculated using DFT, which is a common feature in EAM potentials.³⁴ DFT overestimates the surface energy compared to experiments.⁶⁰ We note that a comparison of the potential energies for the formation of bare Ag surfaces in Table II indicates that these are very close to the zero-temperature values from our EAM potential in Table III. Thus, differences between the surface energies at 413 K and 0 K can largely be attributed to entropy. It is known from experimental⁶¹ and MC simulation^62,63 studies that the surface free energy decreases with increasing temperature. Our multi-scheme TI method is consistent with the literature by estimating the $γAg$ at 413 K about 40 mJ/m² lower than the zero K $γAg$ in Table III for each facet.

Although we find that EG does not significantly influence the trend in the interfacial free energies of the bare Ag surfaces, our methodology can be used to investigate broader topics, such as the thermodynamic role of structure directing agents (SDAs) in solution-phase shape-controlled metal nanocrystal syntheses. In solution-phase shape-controlled metal nanocrystal syntheses, SDAs are the key chemical species that direct shape evolution during all three stages: nucleation, seed formation, and growth, even at a low concentration. In the case of Ag nanocrystal formation, EG is used as the solvent, and polyvinylpyrrolidone (PVP) is a well-known effective SDA that assists to produce 100-faceted structures.^31,32 Because the bare Ag(111) surface is more energetically stable than Ag(100) and Ag(110); a hypothesis for observing 100-faceted structures in the presence of PVP is that, due to a stronger binding affinity of PVP on Ag(100) facets,^35–38 it can stabilize Ag(100) surfaces by making the solid-liquid interfacial free energy of Ag(100) lower than that of Ag(111). Our method can be applied to estimate the interfacial free energies of Ag facets in EG-PVP solution since it is not limited to pure liquid. From our current results, we can exclude the thermodynamic effect from the EG solvent because the trend and free-energy difference are similar to those of bare Ag surfaces.

Inspired by a series of interfacial free-energy studies that share the same “cleaving” idea,^15–17 we developed a six-step multi-scheme TI method to calculate the solid-liquid interfacial free energy using open-source simulation packages. The composition in the liquid phase has no restriction and can be either a pure liquid or a solution of multiple species. As in previous studies, we use TI with $λ$-coupling to calculate the free-energy changes for most steps of our method. However, we utilize the equivalence of the free-energy difference to the reversible work to obtain the free-energy change in going from a liquid-interfacial system to bulk liquid, thereby avoiding special coding to manipulate the PBC of the simulation cell. Another feature of this method is a numerical approximation to obtain the thermodynamic integrands when analytical differentiation is cumbersome to carry out. This feature allows us to deal with a wide range of force fields appropriately, such as the EAM potential used in this work.

Using our method, we obtained the interfacial free energy of Ag-EG interfaces and surface free energies of bare Ag(100) and Ag(111). We successfully resolved the anisotropy by obtaining significantly different $γAg(100)$ and $γAg(111)$⁠. The surface free energies of bare Ag surfaces are consistent with theoretical predictions at zero K that $γAg(100)$ is higher than $γAg(111)$⁠, and also exhibits the temperature dependency expected from other studies. The addition of liquid EG lowers the surface free energy of Ag but does not change the difference between Ag(100) and Ag(111) facets significantly. For this system, we find that the largest contribution to $γsl$ is the free energy to create the bare metal surface. The second-largest contribution to $γsl$ is from the liquid-solid interaction. By comparing the free energy to the potential energy, we find that the trend in free energy mirrors that of the potential energy. However, the entropy contribution is still significant and cannot be neglected. The signs and trends in the entropy of surface formation agree with prior studies that an increase of entropy is always associated with surface formation and its magnitude depends on how much geometric constraint is lost for the surface atoms.

The multi-scheme TI method can be applied to estimate the solid-liquid interfacial free energy in a wide range of systems of interest. One area of interest that motivated the development of this method is understanding the morphologies that develop in solution-phase nanocrystal growth, where crystal shapes can be significantly influenced by additive molecules, or structure-directing agents.^31,32

See supplementary material for a summary of the number of stages of each TI step, properties of the liquid phase, a test of the thermostat validity, and a test for hysteresis.

This work is funded by the Department of Energy, Office of Basic Energy Sciences, Materials Science Division, Grant No. DE-FG02-07ER46414. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) supported by Grant No. NSF/OCI-1053575.