A Computational Assessment of the Efficacy of Halides as Shape-Directing Agents in Nanoparticle Growth

We report a comprehensive study of aqueous halide adsorption on nanoparticles of gold and palladium that addresses several limitations hampering the use of atomistic modeling as a tool for understanding and improving wet-chemical synthesis and related applications. A combination of thermodynamic modeling with density functional theory (DFT) calculations and experimental data is used to predict equilibrium shapes of halide-covered nanoparticles as a function of the chemical environment. To ensure realistic and experimentally relevant results, we account for solvent effects and include a large set of vicinal surfaces, several adsorbate coverages as well as decahedral particles. While the observed stabilization is not significant enough to result in thermodynamic stability of anisotropic shapes such as nanocubes, non-uniformity in the halide coverage indicates the possibility of obtaining such shapes as kinetic products. With regard to technical challenges, we show that inclusion of surface-solvent interactions lead to qualitative changes in the predicted shape. Furthermore, accounting for non-local interactions on the functional level yields a more accurate description of surface systems.


I. INTRODUCTION
Wet-chemical synthesis has emerged as one of the primary routes for obtaining metal nanoparticles (NPs) with tailored properties. One of the most important processes in the aqueous synthesis environment is the adsorption of ions on the NP surface. Detailed understanding of the influence of ion adlayers on growth, morphology, and physicochemical properties of NPs is therefore crucial for the continued development of nanotechnology.
Halides are one of the prototypical ionic species found in aqueous NP solutions as components of common metal salts, surfactants, and cleaning agents 1,2 . In aqueous solution, halides tend to form adlayers on metal surfaces 3,4 , resulting in lower surface energies and stabilization of higher-index facets. The propensity of halides to adsorb is highly sensitive to the atomic structure of the surface, resulting in facet-dependent adlayer structures. Since this can in turn stabilize different NP morphologies [5][6][7] and guide NP growth, halides are commonly considered as "shape-directing agents". For instance, during the growth of gold nanorods bromide plays an integral role in promoting the growth of, initially spherical, seed particles into rods [8][9][10][11][12] . Similarly, several studies have found the presence of bromide conducive to Pd nanocube growth 7,13,14 . Halides can also influence NP growth through other pathways, including the formation of metal complexes with the precursor and surface smoothing through oxidative etching 15,16 .
While the importance of parameters such as halide species, concentration, and counter ion is generally acknowledged, a quantitative understanding of their impact has yet to be achieved. This can partly be attributed to a gap between experimental and theoretical work. Experiments take place in a complex chemical environment where the metal precursor commingles with other reactants such as surfactants and reducing agents, which is often neglected completely in modeling work. This holds true in particular for density functional theory (DFT) calculations, which despite being the main workhorse of computational surface science face technical challenges, e.g, with respect to obtaining accurate surface energies 17,18 and accounting for solvation effects [19][20][21] . A comprehensive analysis of halide adsorption must therefore address multiple aspects. Here, as a step toward this goal, we present an extensive quantitative investigation of halide adsorption in aqueous solution on Au and Pd surfaces as well as NPs, which are of interest for a wide range of applications in nanomedicine [22][23][24] , clean energy 25,26 and catalysis [27][28][29] . We employ a thermodynamic model that uses a combination of data from DFT calculations and experimental free energies to create Wulff constructions. Our model enables quantitative predictions for equilibrium shape and surface area density of halide-covered NPs as a function of concentration in the range relevant for synthesis. In contrast, previous DFT-based investigations into the thermodynamics of halide adsorption on metal surfaces have largely approached the problem from an electrochemical perspective, describing the stable adsorbate geometry on single surfaces as a function of the electrode potential [30][31][32] . The combination DFT calculations with experimental data allows us to compensate for short comings in DFT calculations, including ionization and solvation, without increasing the computational burden. The inclusion of a large number of high-index facets renders the predicted shapes more realistic. The next two sections introduce the key features of the thermodynamic model and provide a summary of the computational methodology, which is followed by results and discussion. Additional material, including a detailed description of the modeling approach as well as complementary results, is presented in the Supplementary Information (SI).

II. THERMODYNAMIC MODEL
In equilibrium, a particle with fixed volume minimizes its surface free energy by assuming a shape described by a set of points where γ [n] is the surface free energy 33 . This is known as Wulff's theorem and W is referred to either as the Wulff shape (WS) or Wulff construction. For a given crystalline NP of material M, it is typically sufficient to consider only a small set of orientations described by Miller indices {hkl}.
If the twin boundary energy of M is known, Wulff's theorem can be generalized to twinned particles, the most prominent examples are decahedra and icosahedra. The latter comprise five and twenty tetrahedral grains, respectively, that are patched together to fill space 34,35 .
We are interested in the case where a NP is in equilibrium with a source of adsorbates X. The free energy can then be written as the sum of the free energy of the clean surface and the change induced by adsorption of N units of X, i.e.
For an aqueous halide X − aq , we can approximate the surface free energy of adsorption on the right-hand side of Eq. (2) using total energies from DFT calculations by To avoid a charged slab calculation, a thermodynamic cycle 36 is used to replace E DFT M(hkl) : N X − aq that would otherwise appear in Eq. (3) with the energy of a neutral system E DFT [M(hkl) : N X] from which the appropriate multiple of the work function Φ [M(hkl)] is subtracted.
In Eq. (3), the chemical environment, i.e., temperature, pressure, and concentration, is represented by the chemical potential µ. The surface free energy of adsorption is a linear function of µ where the slope is equal to the adsorbate surface area density n s = N/A. For fixed (hkl), the slope is furthermore proportional to the coverage Θ = N A prim /A, where A prim is the area of the primitive surface cell. To obtain the dependence of ∆γ on the chemical environment we assume ideal solution conditions, under which Here, the •-superscript indicates that the quantity is given at standard state, defined by P • = 1 bar and c • = 1 M. The first term on the right-hand side of Eq. (4) is thus the chemical potential at the standard state conditions. Here, we calculate its value using a scheme in the spirit of Ref. 21, where DFT calculations are combined with experimental reaction energies, in the present case taken from the curated and consistency-checked data provided in Ref. 37. (Details are provided in the first section of the SI.) Furthermore, the two first terms on the righthand side of Eq. (3) contain surface-solvent interactions.
Since explicit inclusion of water molecules at the DFT level is computationally expensive, they are accounted for here using an implicit solvation model 38,39 .

III. COMPUTATIONAL DETAILS
All DFT calculations were carried out with the Vienna ab initio simulation package 40 (VASP, version 5.4.4), which uses plane-wave basis sets 41 and the projector augmented wave method 42,43 . To accurately model surfaces well as bulk systems, the van-der-Waals density functional with consistent-exchange (CX) 44 was employed. For comparison, a subset of the calculations was also carried out using the local density approximation (LDA), Perdew-Burke-Ernzerhof (PBE), 45 and Perdew-Burke-Ernzerhof for solids (PBEsol) functionals 46 . The planewave basis set was truncated at a cut-off energy of 450 eV, and the electronic levels were smeared using a Gaussian scheme with smearing parameter σ = 0.1 eV. For Brillouin zone integration, Monkhorst-Pack and Γ-centered grids were used depending on the symmetry of the unit cell, with a k-point density of 0.2Å −1 . Geometry optimization was performed for all systems, during which the atomic positions were allowed to relax until the forces were less than 10 meVÅ −1 .
For the calculations of the bare surface energies, we employed the extrapolation approach described in Refs. 47,48 using 10 slabs of different thickness for each (hkl). For the adsorbate calculations, surfaces were modeled using slabs with an approximate thickness of 5a 0 / √ 3 where a 0 is the lattice constant of the parent crystal. Spurious interactions with periodic images were minimized by introducing a vacuum region of 20Å in the direction normal to the slabs. To account for interactions between surface and solvent, the VASPsol implicit solvation model 38,39 was applied to both clean and adsorbatecovered slabs. In addition to the three low-index surfaces (111), (100), (110), 49 adsorption was also modeled on slabs of the high-index surfaces (210), (211), (221), (310), (311), (321), (322), (331), and (332). For each slab halide coverages corresponding to Θ ∈ {1, 1/2, 1/3, 1/4} ML were considered. On the (111) surface, the Θ = 1/3 coverage was realized in a √ 3 × √ 3 unit cell. For (111), (100), and (110) we included all commonly considered high-symmetry sites. For all other surfaces we employed Bayesian optimization to determine the most stable configurations 50 .
The obtained surface energies were subsequently used to construct WSs of single crystalline as well as deca-hedral and icosahedral NPs and average (facet-weighted) surface energies using the wulffpack Python package 51 .

A. Clean surface properties
We first assess the efficacy of various exchangecorrelation (XC) functionals for accurately representing quantities important in the modeling of surface systems. To this end, the lattice constant, clean surface energy as well as the work function are reported for the late non-magnetic face-centered cubic transition metals (Table I). The set of XC functionals includes LDA, PBE, and PBEsol, all of which are commonly used in surface calculations. This set is extended by the non-local CX functional, since several recent studies have highlighted the advantages of including non-local correlations for obtaining accurate results for surfaces 18,52 as well as bulk systems 53,54 . For each combination of functional and quantity, mean absolute error (MAE) and mean absolute percentage error (MAPE) were calculated relative to available experimental data [55][56][57][58] . For the lattice constants, PBEsol and CX both offer comparable improvement over LDA and PBE, respectively. The calculated values are in quantitative agreement with previous benchmarks for bulk systems 54 , where it was found that CX outperforms PBEsol also when a wider range of bulk properties and materials is considered.
Experimentally, surface energies are commonly derived from surface tension measurements on liquid droplets, which provide a facet-weighted average surface energy 59 . Here, a significant improvement is seen in the MAPE values in going from conventional generalized gradient approximation (GGA) functionals to CX. In particular, while PBE gives reasonable results for surface energies of simple metals 59 , the MAPE of 35% found here suggests that caution should be exercised in the application of PBE to late transition metals. It is well-known that LDA, on the other hand, yields relatively accurate surface energies 60 , which can, however, be attributed to an error-cancellation effect 18 . As can be seen from the MAPE values the systematic extension of DFT to include non-local interactions represented by CX leads to a slight improvement over LDA. The systematic increase in surface energy seen when comparing CX to either PBE or PBEsol reflects the contribution of non-local correlation to the energetic cost of cleaving a crystal to create a surface.
The {111} work functions exhibit similar trends as the lattice constants, with PBE and LDA systematically under and overestimating experimental data, respectively. The MAPE is reduced by about 2% when either PBEsol or CX are used instead, with the latter being slightly closer to the experimental values.
Overall, CX outperforms the other functionals over the range of elements and properties included in this comparison. Since both surface energy and work function appear explicitly in the surface free energy of adsorption in Eq. (3), CX is thus the preferred functional for the present purpose.
While XC functionals can exhibit large differences in predictions of surface energies, the equilibrium shape of a NP is sensitive to surface energy ratios rather than absolute values by virtue of Wulff's theorem Eq. (1). For this reason, we also compare actual WSs obtained via the above set of functionals for Au and Pd. To give an accurate representation of the WS, the set of three lowindex facets is further augmented by several high-index facets.
The WSs for single crystalline particles (Fig. 1) reveal that, while both Au and Pd particles have shapes derived from the regular truncated octahedron (RTO) with additional faceting provided by high-index surfaces, the fractional area occupied by the facets varies not only in magnitude, but also in ordering when the XC functional is changed. The differences are most pronounced in the case of Au, where a large portion of the surface is occupied by high index facets. In conclusion, the choice of XC functional has ramifications that go beyond the degree of numerical agreement with experiment, and can lead to qualitative differences in predicted NP shapes.

B. Particle morphology and the halogen group
Before presenting any results related to halide adsorption, we note that typical wet-chemical synthesis concentrations are on the order of 0.01-1 M. Hence, we use c X − aq = 0.1 M as a representative value for aqueous states throughout the remainder of the text, unless the concentration dependence is explicitly considered. Similarly, while variations in temperature and pressure can in principle be accounted for in Eq. (3), we restrict ourselves to ambient conditions. The basic features of halogen adsorption on metal surfaces can be appreciated by observing the changes in surface adsorption energy as the adsorbate state is changed from the dimerized state to the aqueous ion. This is illustrated in Fig. 2 for the (100) surface of Au and Pd. The results for the aqueous ions are given both with and without corrections for surface-solvent interactions provided by the VASPsol implicit solvation model. 61 Note that we use the thermodynamic sign convention where more negative surface free energies indicate greater stability. Furthermore, Fig. 2 shows the surface adsorption energies as given by Eq. (3) where the contribution from the clean surface is not included. Negative (positive) energies thus correspond to exothermic (endothermic) reactions. In the dimerized state, halogen adsorption strength correlates closely with electronegativity. This is clearly observed, for example, on Au(100) (Fig. 2a, left) where the most electronegative halogens are also more strongly adsorbed. In going from dimer to solvated ion, however, the trend is reversed (Fig. 2a,b middle) and the adsorption strength decreases in the order I − > Br − > Cl − > F − . This reversal reflects differences in ion-solvent interaction, where light halide ions interact more strongly with the solvent resulting in weaker adsorption 62 . Similar results have previously been obtained from PBE calculations complemented by thermodynamic data 12 . Including the surface-solvent interaction shifts the adsorption energy upward by more than 0.2 eV. For Cland Bron {100}, this shift is as large as the change in adsorption energy in going from dimerized state to aqueous ion, and for Fit can even result in a sign change of the adsorption energy. Therefore, it is clear these interactions cannot, in general, be neglected.
We now consider the influence of halide adsorption on NP shape. Here, the use of the generalized Wulff constructions enables us to also include decahedral NPs in On Au NPs, adsorption of F - (Fig. 3b,f) is largely endothermic, and a distribution of facets primarily comprising {111} and vicinal facets (Fig. 3b, f) is obtained, identical to a clean NP under aqueous solvation conditions (Fig. 3a,e). Adsorption of Clis, on the other hand, exothermic and leads to changes in faceting. While the area fraction of {111} remains similar to the clean case, with around 20% for both single-crystalline and decahedral NPs, most high-index facets are eliminated. Indeed, the facet distribution is dominated by {311}, which accounts for 70% of the available area (Fig. 3c,g). Upon replacement of Clby Br -, the area fractions again undergo significant changes resulting in a WS with 50% {111}, 20% {100} as well as smaller amounts of {210} and {310} (Fig. 3d,h). In the case of Pd, Fadsorption can be exothermic, in particular {110} facets are stabilized relative to clean, aqueous NPs (Fig. 4a,b). For singlecrystalline NPs this yields a rhombic dodecahedron, while one obtains a truncated pentagonal bipyramid in the case of decahedral symmetry (Fig. 4e,f). In contrast to Au, the shapes obtained in the presence of Clor Brare similar to each other (single-crystalline NPs: Fig. 4c,d; decahedral NPs: Fig. 4g,h) with about 80% of their surfaces being {210} facets alongside smaller fractions of {111} and {100}. We note that the prevalence of {210} facets in this case is largely due to geometric constraints since the (100) surface free energy is of similar magnitude for a wide range of concentrations (Fig. S2).
In combination, the significant differences in faceting with halide adsorption and NP composition show that morphology prediction for halide-covered NPs is nontrivial. It is also noteworthy that without the use of an implicit-solvation model, high-index facets do not appear in any of the above WSs (Figs. S3 and S4). Hence, even if the adsorbate state is properly accounted for in the thermodynamic analysis, failure to include surfacesolvent interactions can result in qualitatively incorrect predictions. In addition to NP faceting, we can study the halide coverage, or surface area density. The analysis leads to the important observation that for both Au and Pd NPs, the area density of adsorbates at fixed bulk concentration is facet dependent. Considering for example Bron Au (Fig. 5a) at a concentration of 0.1 M, the surface area density on {100} and {111} facets is 6.0 nm −2 and 4.6 nm −2 , respectively, corresponding to coverages of 1/2 ML and 1/3 ML. This can yield anisotropic shapes since growth can be hampered in the case of high coverage. In a kinetically controlled synthesis, in which the surface diffusion is slow compared to the arrival rate of monomers at the surface, this would promote growth on {111} and turn octahedra into cubes 65 while increasing the aspect ratio of decahedra 66,67 .

C. Sensitivity to the chemical environment
Our thermodynamic model allows for the prediction of NP energetics and shape changes as a function of thermodynamic control parameters. The most relevant parameter in this case is the halide concentration, which in a typical NP synthesis is on the order of 10-1000 mM. In applications where protecting the NPs against agglomeration is the primary concern, the reduction in the average (facet-weighted) surface energyγ with respect to the clean NP provides a measure of the overall stabilization that can be achieved by halide capping. The stabilization follows the expected energetic ordering Br − < Cl − < F − ( Fig. 6; also see Fig. 2). Increasing the halide concentration yields an exponential decrease inγ with the exception of Fon Au. In the latter case, since the adsorption is endothermic on many gold facets, no stabilization is observed. At a concentration of 0.1 M, Bradsorption leads to a decrease ofγ by 45% and 60% for Au and Pd, respectively. Thus, halide adsorption can significantly lower the average surface energy of the NP providing a quantitative measure for the efficacy of halides as capping agents.
To complement the WSs presented in the previous section, we also map the facet area fractions as a function of the halide concentration (Fig. 7). For brevity, the results presented here are restricted to Br -, since this is arguably the most important halide in shape-controlled NP synthesis. The corresponding maps for Fand Clcan be found in the SI (Figs. S5 and S6.).
We observe that the sensitivity to Br concentration is largest for Au NPs, and crucially, that qualitative changes to the faceting can occur within the range of  inated by increasing the concentration to 3.3 M. For comparison, the area fraction at the representative concentration 0.1 M is 20% (Fig. 7a, For Pd NPs, the facet area fractions (Fig. 7c,d) change monotonically with concentration and at a lower rate than for Au. Regardless of the concentration, {210} facets represent about 80% of the surface area and as the concentration increases the fraction of {100} facets grows at the expense of {111}.
The above analysis shows, in conjunction with the results presented for the halide substitutions ( Fig. 3 and  4), that the shape of Au and Pd NPs can be tuned to some extent by changing either halide species or concentration. The effect is, however, not significant enough to render anisotropic shapes such as cubes thermodynamically stable.

V. CONCLUSIONS AND OUTLOOK
In this study, we have demonstrated how DFT calculations combined with an implicit-solvation model and experimental data can be used to formulate a thermodynamic model with the ability to predict the equilibrium shape and surface area density of halides on aqueous NPs under realistic conditions, including finite temperature, pressure and halide concentration. Our contribution thus helps bridge the gap between experimental and theoretical studies of aqueous NPs.
For Au and Pd particles, halide species and concentration play an important role in determining NP faceting. In particular for Au NPs, qualitative differences in faceting can be obtained by varying the concentration of Cland Br -, the two halides most commonly encountered in wet-chemical synthesis. Here, the nonmonotonic response of the distribution of facet area fractions to changes in concentration underscores the subtlety of NP shape prediction in complex environments and underscores the need for accurate atomistic simulations. We emphasize, however, that adsorption of halides alone cannot stabilize anisotropic shapes such as cubes or rod-like particles.
Furthermore, our calculations reveal that the halide surface area density is facet dependent and for both Au and Pd NPs the highest density is found on the {100} facets. Our results are thus consistent with the picture of halides as capping-agents that can selectively passivate one or more facets and, given the right geometry and kinetic conditions, promote their growth 66 .
In light of the above conclusions regarding the relation between halides and shape anisotropy, it is evident that further investigations into the synergistic effects between halides and their cations, most notably cetyltrimethylammonium bromide (CTAB), are required. Some steps in this direction have already been taken using classical molecular dynamics (MD) 11,68 as well as DFT simulations 12,69 . In the case of classical MD, the reliability of any force-field description of the halide-metal interface must be carefully evaluated since traditional force-fields do not account for effects such as charge transfer and polarization at the interface 70,71 . In the case DFT calculations, the thermodynamic model utilized in the present work is in principle straightforward to apply to systems where both anions and cations are present, but the accuracy is limited by the quality of the available thermodynamic data.
Another aspect of NP shape prediction that was deemed to be outside the scope of the present study concerns the treatment of finite size effects. It is well-known that the Wulff construction does not account for such effects and the predicted WSs can thus only be considered realistic for large particles. Quantifying the limit for when a NP can be considered large in this sense is important given the high area fractions predicted here for high-index facets such as {210} and {310}. More precisely, from a geometric point of view a NP needs Weighted average surface energy of a NP gives a simple measure of the NPs stability against agglomeration or coalescence. The grey line indicates the concentration in Fig. 2 to 4. The average surface energy decreases with increasing concentration and with increasing atomic number of the halide. For Brconcentrations on the order of 1 M, the average surface energy can decrease by as much as 50% relative to that of clean, aqueous particles, thus suggesting the use of Bras an effective capping agent. to be large enough to accommodate stepped surfaces in the first place, and the energetic penalty incurred by the facet edges and corners needs to be considered. As a result, a small NP may end up exposing smaller amounts of vicinal facets then what is predicted by a particular WS.