Predicting kinetic nanocrystal shapes through multi-scale theory and simulation: Polyvinylpyrrolidone-mediated growth of Ag nanocrystals

Colloidal metal nanocrystals with well-defined shapes exhibit properties that are attractive for a diverse array of applications. With robust shape-control techniques, morphology-dependent properties of the nanocrystals can be tailored towards specific applications in areas such as catalysis,^1–4 plasmonics,^5–9 and photovoltaics.^10–12 Although many syntheses yielding monodisperse samples of nanocrystals have been reported,^13,14 the molecular understanding of their shape-control mechanism still remains elusive.

There are many accounts of solution-phase syntheses of Ag nanocrystals via reduction of AgNO₃ salt in an ethylene glycol (EG) solution in the presence of polyvinylpyrrolidone (PVP),^15–19 making this a convenient model system for studying the solution-phase shape control mechanism. By varying synthesis parameters, different Ag nanostructures such as cubes,¹⁹ bipyramids,²⁰ and wires²¹ can be obtained in these solution-phase syntheses. While Ag nanocrystals with well-defined shapes can be successfully obtained, the atomistic details of how such shapes emerge are rather ambiguous.

In the solution-phase synthesis of Ag nanocrystals, PVP likely plays the role of a structure-directing agent (SDA) which helps to mediate the final nanocrystal shape. The structure-directing capability of PVP has been attributed to preferential binding of PVP to Ag(100) over Ag(111) facets.^22,23 In the case of fivefold-twinned Ag nanowires, it was hypothesized that PVP passivates the {100} side facets and leaves the {111} end facets essentially uncovered.²³ This allows Ag atoms to deposit selectively at the {111} end facets relative to the {100} side facets, to grow nanowires. Preferential binding of PVP to Ag(100) has been confirmed by evidence from scanning electron microscopy (SEM) that shows the replacement of PVP by dithiol molecules at the {111}-faceted nanowire ends, but not at the {100}-faceted sides.²³ Preferential binding of PVP to {100} over {111} facets is also seen in dispersion-corrected density-functional theory (DFT) calculations of PVP binding energies to Ag surfaces.^24–29 

A deeper question arises of how the preferential binding of an SDA to a particular metal facet could lead to shape-controlled synthesis in which that facet is expressed. Surface-adsorbed SDAs, such as PVP, prevent random nanocrystal aggregation and allow for well-defined shapes to emerge. They can influence the surface free energies of the nanocrystal facets, which could promote various thermodynamic shapes. SDAs can also influence the kinetics of Ag atom deposition and surface diffusion, which could lead to kinetic shapes.¹⁸ When the rate of deposition is much greater than the rate of inter-facet surface diffusion, growth is kinetically controlled. In this case, a facet’s growth rate is directly proportional to the rate of deposition onto the facet.²⁹ Under kinetic control, the facets with the slowest linear growth rates (normal to the facet) are expressed in the final nanocrystal shape. On the other hand, growth is thermodynamically controlled when the rate of inter-facet surface diffusion is significantly greater than the rate of deposition. Atoms diffuse between facets to minimize the surface free energy of the nanocrystal. Given enough time for equilibration, the final nanocrystal shape will be made up of facets that minimize the surface free energy under thermodynamic control.

One likely example of kinetic nanocrystal shapes can be found in experimental studies of the PVP-mediated growth of large Ag nanocrystals from cubic Ag seed crystals in an EG solution.¹⁹ In these studies, they found that {100}-faceted cubic seeds continue to grow as cubes if the concentration of PVP is sufficiently high. Below a critical PVP concentration, the cubic seeds evolve into cubo-octahedra and eventually into {111}-faceted octahedra. The edge lengths of the seeds used in the experimental study are ∼40 nm and the seeds grow to ∼100 nm dimensions,¹⁹ which is too large for inter-facet diffusion to be significant.

In recent work, we computed the fluxes of solution-phase Ag atoms to Ag(100) and Ag(111) surfaces with adsorbed PVP using direct molecular dynamics (MD) simulations.³⁰ These simulations showed that there is a higher flux of solution-phase Ag atoms to Ag(111) than to Ag(100), due to the stronger binding of PVP to Ag(100). The ratio of Ag(111)/Ag(100) fluxes was large enough to predict cubic Ag(100) nanostructures for sufficiently high PVP surface coverages within the framework of the kinetic Wulff construction, consistent with experiment.¹⁹ 

We note that solution-phase Ag fluxes can also be obtained from the mean-first passage time (MFPT) t_m, by solving the Smoluchowski equation for the diffusion of the Ag atom in the potential of mean force (PMF).³¹ For an atom beginning at a height of z above a surface (see Fig. 1) that is absorbed irreversibly onto the surface at a height of z_f, t_m is given by

where W(z) is the PMF, or the free-energy profile of the system as the atom transits from the bulk solution to the surface, and D(z) is the diffusion coefficient of the atom. The diffusion coefficient is in general a function of the atom’s height above the surface z and this is evident in Fig. 1 (produced using visual molecular dynamics software³²), as the Ag atom must diffuse from the bulk solution through an adsorbed PVP layer to reach the Ag substrate below. The reflecting boundary at z₀ represents the boundary between the near-surface region and the bulk-solution region, as shown in Figure 1. We note that the PMF and D(z) can be obtained via free-energy sampling techniques. For systems in which there is a large free-energy barrier for deposition, this may be the only feasible approach, since the time scale for deposition could exceed the MD time scale.

In this work, we present a multi-scale theoretical framework for predicting kinetic nanocrystal shapes. Unlike previous approaches,^33–35 we account for the solution-phase and surface-adsorbed films explicitly, as these factors can be highly influential in determining facet growth rates. We consider the case that the crystal growth rate is limited by the flux of solution-phase atoms to the crystal surface. Using umbrella-sampling MD simulations, we calculate the PMF and the diffusion coefficient along the deposition path of an Ag atom approaching an Ag surface. We incorporate these quantities into Eq. (1) to obtain the MFPT and we compare these values to those from direct MD simulations and from mesoscopic Brownian dynamics (BD) simulations. We employ this protocol to investigate PVP’s flux-regulating capabilities in the growth of Ag nanocrystals in an EG solution. By incorporating the calculated MFPTs into the kinetic Wulff construction, we use our results to predict kinetic Ag nanocrystal shapes. Our results are consistent with experiment.¹⁹ This protocol can be applied to predict nanocrystal shapes in similar systems.

The simulation systems are constructed to mimic the analogous Ag/EG/PVP experimental system of the seed-mediated growth of Ag nanocubes.¹⁹ We include EG as an explicit solvent and we use atactic PVP oligomers of various lengths to represent the SDA.²⁹ We represent the {100} and {111} facets of the Ag nanocrystal seeds as periodic infinite slabs. The depositing Ag atom is modeled as charge neutral.

We performed all MD simulations using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code,³⁶ compiled with our extension for calculating the van der Waals interactions between Ag and organic species.³⁷ The equations of motion were integrated using the velocity-Verlet algorithm³⁸ with a 1 fs time step for equilibration. Temperature and pressure were controlled using the Nosé-Hoover thermostat and barostat.^39,40

We use the Metal-Organic Many-Body (MOMB) force field description for the Ag-EG-PVP system developed by Zhou et al.^29,37 to match results from dispersion-corrected DFT calculations. In this force field, the embedded atom method (EAM)⁴¹ is used to describe Ag-Ag interactions and the CHARMM force field^42,43 is used for organic-organic (EG and PVP) interactions. The Ag-organic interactions are described using pair-wise potentials, consisting of a Morse potential and Grimme’s form for van der Waals dispersion interactions,⁴⁴ augmented with a many-body EAM term to describe the chemical bonding of the oxygen in PVP and EG to silver.³⁷ 

All simulation boxes are constructed with periodic boundary conditions in all directions. We study PVP pentamers (PVP5), decamers (PVP10), and icosamers (PVP20). In simulations of PVP5 and PVP10, we place a (14 × 14 × 12) Ag(100) slab or a (14 × 16 × 12) Ag(111) slab in a simulation box filled with EG solvent molecules and an amount of PVP in excess of that to produce an adsorbed monolayer. For PVP20, we use Ag slabs with larger in-plane dimensions: (24 × 24 × 12) for Ag(100) and (24 × 28 × 12) for Ag(111) to accommodate the longer length of this molecule.

All simulation systems are equilibrated before production to ensure full relaxation of PVP chains on the Ag surfaces. We first relaxed systems with PVP5 and PVP10 in the NPT ensemble for 10 ns at a temperature of 433 K and a pressure of 1 bar to adjust the box size, and we further equilibrated them in the NVT ensemble at 433 K for an additional 45 ns, which includes 20 ns of parallel tempering using 8 replicas ranging from 433 K to 503 K with an interval of 10 K per replica. Systems with PVP20 are relaxed in the NPT ensemble for 20 ns followed by 80 ns of NVT equilibration, with no parallel tempering. At the end of equilibration, both Ag surfaces are covered with PVP chains and excess PVP chains reside in the solution phase. We removed the excess PVP and all PVP chains adsorbed on one of the two surfaces so that we have a monolayer of adsorbed PVP oligomers on one Ag surface in all systems. After removing the excess PVP chains, all systems are again relaxed in the NPT (5 ns) followed by the NVT ensemble (15 ns) to ensure equilibration. The number of each type of molecules is given in the supplementary material. A typical system configuration after equilibration is shown in Figure 1.

We compute the PMF using umbrella sampling⁴⁵ and umbrella integration.⁴⁶ Umbrella sampling was performed using the collective variables module in LAMMPS.⁴⁷ An exhaustive description of the procedure is given by Kästner and Thiel.^46,48 To initialize the process, an Ag atom is added into the bulk solvent and it is held in the bulk solution using a biasing harmonic potential to generate initial states for umbrella-sampling simulations. From each initial state, the Ag atom is “pushed” through the solution onto the Ag surface using a moving harmonic potential to generate initial states for windows along the path of Ag atom deposition. The orthogonal distance of the Ag atom from the Ag surface is the reaction coordinate, which is denoted as z. When a “pushing” simulation is performed from different initial states, different deposition paths are observed. These dissimilar sampling paths reflect the inhomogeneity of the PVP film.

We note that it is also possible to generate initial trajectories using “pulling” simulations, in which the Ag atom begins at the Ag surface. In fact, we performed both types of simulations (we reported the results of a “pulling” simulation for an Ag atom from a PVP10 film previously³⁰) and we observed that the PMF profiles in the “pulling” simulation have a longer tail in the bulk solvent phase—especially with PVP20. This occurs when the Ag atom is pulled because PVP chains tend to attach to the atom and follow it into solution. This long-range Ag-PVP attachment does not occur in most trajectories for direct MD simulations,³⁰ where solution-phase Ag atoms frequently bypass the PVP tails. By doing “pushing” simulations, this artifact is minimized. In addition, initialization of the “pushing” simulations can be de-phased to generate different deposition paths more efficiently than in “pulling” simulations. This is because the Ag atom diffuses faster in the solution phase than on the Ag surface. The PMF profiles computed from “pulling” simulations are not suitable for the current work because the artifacts introduced will significantly influence the quantitative predictions made here, although previous qualitative interpretations³⁰ of the PMF profiles are not as sensitive to these artifacts.

After the initial “pushing” simulation, we begin umbrella sampling by applying a biasing harmonic potential centered on the Ag atom for every window along the deposition path. At each window, we run canonical MD simulations initiating at the “pushing” coordinates and we sample the reaction coordinate z every 100 time steps. From each window of the umbrella sampling MD simulation, we obtain a biased histogram of the reaction coordinate. Using the umbrella integration technique,⁴⁶ the PMF profile can be extracted from the biased probability-density histogram combined from all windows.

During umbrella sampling, we use a 1.5 fs time step with C–H and O–H bonds constrained using the SHAKE algorithm.⁴⁹ The umbrella-sampling parameters are chosen based on guidelines given by Kästner and Thiel.⁵⁰ These parameters consist of the force constant of the biasing harmonic potential, the window spacing, and the number of windows. These values are reported in Table I. By varying the force constant for the PVP5-Ag(100) system, we observe that the PMF profiles obtained from different force constants converge, as shown in Fig. S1 in the supplementary material. We sampled 3 ns per window, which leads to 120–270 ns per sampling path, combining all windows sampled. Using 3 different sampling paths to ensure the PMF profiles are converged, we have 360–810 ns of simulation time for each PMF profile computed.

The PMF profiles for Ag atom deposition onto Ag(100) and Ag(111) are shown in Figure 2 (plotted using Python’s matplotlib package⁵¹). We note, as discussed above, that each of the PMF profiles in Fig. 2 is an average over three individual profiles. Due to the heterogeneity of the PVP films, which is evident in Fig. 1, there is variation among the individual profiles, as can be seen in Fig. 3. Nevertheless, the individual pathways exhibit similar general features and, as discussed above, the averages in Fig. 2 do represent considerable sampling of Ag atom deposition pathways. We ensured that the Ag atom has sampled the entire x–y space within the PVP film by mapping a 2D histogram to show regions that the Ag atom had sampled, as shown in Fig. S2 in the supplementary material.

From Fig. 2, we see that in all cases, the PMF profile is flat in the bulk solvent phase far away from the Ag surface. The bulk solution value is set to be zero in all systems, as the reference value. As the Ag atom begins to contact the edge of the PVP film, the PMF decreases from its solution-phase value. For PVP5 and PVP10 films, the location where the PMF begins to decrease is further away from Ag(111) than Ag(100). This suggests enhanced trapping of the Ag atom by the extended PVP film on Ag(111), which we previously proposed as one of the flux-regulating capabilities of PVP.³⁰ The decrease of the PMF from its bulk value is more gradual and there are less discernable differences between Ag(100) and Ag(111) for PVP20 films. We will return to this point below.

For all PVP-covered surfaces, the PMF decreases until it reaches a local minimum, where it is held in the PVP film prior to reaching the Ag substrate, which are global minima in Fig. 2. To reach the Ag substrate, the atom must overcome a free-energy barrier. We denote this barrier as the “hole-opening barrier” because it characterizes the Ag atom’s search for an opening within the adsorbed layer that allows it to reach the underlying substrate. In direct MD simulations,³⁰ we observed that the Ag atom finds a site on the Ag substrate either by diffusing within the adsorbed layer to find an opening in the film or by rearrangement of the adsorbed layer to open a hole underneath the Ag atom. We note that the hole-opening barrier for the Ag atom to reach the PVP-covered Ag surfaces is much smaller than the free-energy difference between the local minimum and the bulk solution. Thus, the Ag atom is held essentially irreversibly in the PVP film before it continues to the Ag substrate.

We can characterize the hole-opening barrier by the difference between the PMF value at the local minimum and that at the adjacent local maximum as the Ag atom approaches the surface. Table II shows the values of these barriers for the different PVP systems. For PVP10 and PVP20, the hole-opening barriers for Ag(111) are smaller than for Ag(100). This implies that the Ag atom can more easily access Ag(111) than Ag(100). It is consistent with a larger linear growth rate for Ag(111) facets and {100}-faceted nanostructures, and consistent with experiment.¹⁹ For PVP5, the hole-opening barriers for both Ag(100) and Ag(111) are similar in value. We can deduce that the critical chain length which allows PVP to regulate Ag fluxes is above five units.

The PMFs for the PVP-covered Ag surfaces are distinctly different from those for Ag surfaces in EG solvent only [Fig. 2(a)]. In the absence of PVP, the most prominent feature in the PMF is the free-energy barrier that the solution-phase Ag atom encounters in its approach to the Ag surface. Unlike the PVP layer, which attracts a solution-phase Ag atom and holds it before it reaches the surface, the pure EG solvent layer repels the Ag atom with a barrier to adsorption. As we will elaborate later, a solution-phase Ag atom approaching an Ag substrate with the pure EG layer frequently bounces off the free-energy barrier and returns into the solution phase instead of successfully reaching the Ag substrate.

Major features in the PMF profiles in Fig. 2 coincide with minima and maxima in the solvent and PVP density profiles. For the EG-only systems [Fig. 2(a)], minima and maxima in the PMF are due to the layering of EG induced by the Ag surfaces. Evidence for EG layers can be seen in the EG density profile shown in Fig. 4(a). Here, we see that the location of the highest maximum in the PMF in Fig. 2(a) coincides with the largest maximum in the EG density profile in Fig. 4(a). This maximum coincides with an adsorbed layer of EG immediately adjacent to the Ag surfaces. There is also a correspondence between the locations of the three, smaller maxima in the free-energy profile and those in the EG density profile and this is associated with the progressively weaker EG layering that occurs in the transition from the near-surface to the bulk-solution environment.

In the EG density profiles for the PVP layers [Figs. 4(b)-4(d)], we see that PVP displaces EG from the surface, so its near-surface density becomes less than its bulk-solution density. To understand the PMFs for PVP-covered surfaces, we examine the density profiles for the oxygen atoms in PVP, which are shown in Fig. 5. Since each segment of PVP contains one O atom, the O-atom density profiles are akin to PVP segment-density profiles. In Fig. 5, we see that adsorbed PVP chains in a monolayer have a distinct structure, with two major peaks near the surface. In previous work, we saw similar density profiles for PVP chains adsorbed in a monolayer on Au surfaces.⁵² In that study, as well as in this study, the two peaks could be attributed to two general orientations of the 2-pyrrolidone side rings with respect to the alkyl backbone of atactic PVP: a fraction of the rings is “downward facing,” to form the peak closest to the surface, and a fraction is “upward facing” and oriented away from the surface, to form the second peak. We note that the most prominent local minimum in the PMF profiles for the PVP-covered surfaces falls around the location of the second peak. Since the O atom in the repeat unit of PVP has a particularly strong interaction with Ag,^24,25,37 we conclude that O atoms in these upward-facing 2-pyrrolidone rings are primarily responsible for trapping incoming Ag atoms into the local minimum in the PMF profile.

We note that the location of the hole-opening barrier in the PMF profiles lies around 5 Å from the surface and coincides with the minimum in the density profiles in Fig. 5. Thus, the hole-opening barriers in the PMF profiles can be linked to the transit of Ag atoms from the “upward-facing” 2-pyrrolidone rings through the alkyl backbones, interstitial EG (cf. Fig. 4), and “downward-facing” pyrrolidone rings to the surface.

Examining the O-atom density peaks closest to the surface in Fig. 5, we see that these are higher for Ag(100) than for Ag(111). Considering the hole-opening barriers in Table II, we see that these are higher on Ag(100) than on Ag(111), so a higher PVP segment density near the surface correlates with a higher hole-opening barrier. We also see that the O-atom density profiles are more extended for Ag(111) than for Ag(100). Overall, the O-atom density profiles indicate denser and more compact PVP films on Ag(100) than on Ag(111).

A final comment can be made concerning the long-range regions in the PMF profiles near the bulk-solution interface in Figure 2. For the EG+PVP20 systems, we observe a gradual decrease in the PMF that begins when the Ag atom is more than 30 Å away from the surface. This contrasts with the relatively sharp drop at 18-20 Å in the PMF for EG+PVP5 and EG+PVP10 systems. From Fig. 5, the furthest oxygen-density peak for all PVP systems is at 11-13 Å and the decline in the oxygen density extends to ∼16 Å. The distance between the furthest oxygen-density peak and the beginning of the decline in the PMF is the distance of interaction between the Ag atom and the PVP film.

In the case of PVP5, which is a rigid-rod-like molecule with a length equal to the Kuhn length of PVP in EG,³⁰ the decline in the PMF roughly correlates with the long-range region of the O-atom density profile in Fig. 5(a). However, for PVP10 and PVP20, the PMF begins to decrease over distances for which the O-atom density is negligible. This long-range decrease in the PMF is due to the interaction of the Ag atom with a few tails of PVP chains that extend into solution. The PVP20 films have the longest tails and they can be influenced by the umbrella-sampling process.

During the umbrella-sampling process, the Ag atom is held in place by a biasing harmonic potential over a 3-ns time period. At sampling windows where the Ag atom is held in bulk solution, we observed instances where the PVP oligomers shift their configurations towards an energetically favorable interaction with the depositing Ag atom. This occurs when the Ag atom is sampling the x–y space and it comes close to a PVP tail. PVP20 molecules are most susceptible to this type of configurational shift due to their long tails and flexibility.

The position-dependent diffusion coefficients of the Ag atom along the reaction coordinate [i.e., D(z) in Eq. (1)] are estimated from the simplified Woolf-Roux equation^53,54 from trajectories obtained in the umbrella-sampling MD simulations that we used to obtain the PMF profiles. The Woolf-Roux equation emerges from solving the generalized Langevin equation for a harmonic oscillator and is given by

where μ_i is the average of the reaction coordinate sampled in window i, $σz,i2=〈zi2〉−μi2$ is the variance of the reaction coordinate sampled in window i, and τ_z,i is the characteristic correlation time. The characteristic correlation time τ_z,i is obtained from the sum of autocorrelation from lag k = 0 up to the cutoff lag k_c, according to

where z(t) is the reaction coordinate along the sampled trajectories in window i. In Hummer’s reformulation of this method,⁵⁴ the cutoff lag k_c is implemented to mediate the oscillatory behavior of the autocorrelation function at long-time correlations. For all calculations, we used k_c = 18.0 ps, which gives more weight to the short-time correlations that are important for the calculation of diffusion coefficients. As discussed in the supplementary material, we verified that our diffusion-coefficient profiles are converged with respect to the umbrella-sampling force constant.

Noise in our sampling is attenuated with two serial processes. First, we employ a binning process, in which groups of five adjacent points are replaced by their average, yielding one average point every 2 Å. These values are shown as points in Fig. 6. The standard deviation between the five adjacent points is also calculated and used as the error bars in Fig. 6. A three-point rectangular smoothing is applied to the binned diffusion coefficients, in a second smoothing process, and we perform linear interpolation between these smoothed points to obtain the curves shown in Fig. 6. We conducted a sensitivity analysis to test our smoothing procedure using three-point rectangular smoothed diffusion profiles, binned diffusion coefficients, and raw diffusion coefficients from each window of data. As shown in the supplementary material, we find that our smoothing procedure provides a more regular and easily understandable diffusion-coefficient profile without altering the results.

Figure 6 shows the smoothed diffusion-coefficient profiles of the Ag atom along with binned points for the different systems investigated. In the bulk solvent phase, our diffusion coefficients converge within 21% tolerance to the bulk value of 3.3 ± 0.2 × 10⁻⁵ cm²/s that we found previously.³⁰ For our analysis, we assume the bulk diffusion coefficient to be uniformly 3.3 ± 0.2 × 10⁻⁵ cm²/s for all systems, since the bulk phase of all systems is the same, to allow for a fair comparison of the MFPTs between different systems. This assumption eliminates the fluctuation of the bulk diffusion coefficients that is introduced by the sampling noise and the oscillatory behavior of the autocorrelation function near the cutoff lag time.

For the EG-only systems [Fig. 6(a)], there is no discernable influence of the solvent environment on D(z) except perhaps at the very closest distances to the Ag surface, where there is an increase in the diffusion coefficient after the Ag atom surmounts the free-energy barrier to reach the surface [cf. Fig. 2(a)]. For systems with adsorbed PVP, the strong interaction between the Ag atom and the PVP films reduces the diffusion coefficients in the near-surface region. Similar to the PMF profiles in Fig. 2, we see that the diffusion coefficient begins its decline from the solution-phase value to its surface value further away from the surface for Ag(111) than for Ag(100) in the PVP5 and PVP10 systems. For PVP20, as for the PMF, the decline in the diffusion coefficients is similar on both Ag surfaces. The origins of this behavior are the same as those for the PMF profiles, discussed above.

We calculated the MFPT for the various systems described above using Eq. (1). In addition to the PMF W(z) and the diffusion coefficient D(z), we need to choose values for the reflecting boundary z₀, the absorbing boundary z_f, and the initial coordinate z when integrating Eq. (1).

We take the absorbing boundary z_f to be the location of the global minimum in the PMF for all systems. We choose the reflecting boundary z₀ to match the free Ag atom concentration with the experimental bulk Ag concentration in solution.¹⁹ Here, we assume that the free Ag atom concentration is equal to the bulk AgNO₃ concentration in experiment¹⁹ and we obtain a volume of 30 174 ų per Ag atom. The distance between each Ag atom is calculated by taking the cubic root of the volume, yielding z₀ to be ∼32 Å. We note that this is likely an upper bound on the experimental concentration because it seems unlikely that every solution-phase Ag atom is charge neutral. We will return to this point below. We use the starting position of Ag atom z to also be 32 Å for all MFPT calculations.

To gain insight into the MFPT for the PVP systems, we divide it into two contributions: a diffusion time, which characterizes diffusion of the Ag atom from the bulk solution to the local minimum in the PMF, and a reaction time, which characterizes the time for an Ag atom to proceed from the local minimum in the PMF to the Ag substrate. We do not decompose the MFPT in this way for the EG system.

The MFPT for each system is shown in Figure 7. Interestingly, the MFPTs are significantly longer for the EG-only systems than for those with PVP. This trend originates from the lack of trapping (no decrease in the PMF as the Ag atom approaches the surface) in the EG-only system and the significant free-energy barrier for Ag to reach the surfaces. While the diffusion times are all similar among the PVP systems, the reaction times can vary significantly. For EG+PVP10 and EG+PVP20, the reaction time is shorter for Ag(111) than for Ag(100). This is a consequence of the smaller hole-opening barrier on Ag(111) for these systems (cf. Table II). In contrast, the reaction times are comparable on Ag(100) and Ag(111) for EG+PVP5. Comparing the diffusion and reaction times for the PVP systems, we conclude that Ag atom deposition is neither diffusion-limited nor reaction-limited.

We note that the trends in the MFPTs discussed above are in qualitative agreement with those from our earlier direct MD simulations³⁰ and the MFPTs computed here fall within a factor of two of the direct MD results. One inconsistency between the two results is that in direct MD simulations, the MFPTs increase with increasing PVP chain length and that does not occur in our calculations here. While the origins of this discrepancy are not completely clear, we note that the diffusion times tend to be longer in the direct MD simulations than ours here. This may be a dynamical effect that originates from the heterogeneity of the PVP layer.

Figure 8 shows a typical configuration of a PVP20 layer on Ag(111), where we see that there are relatively flat regions, in which PVP segments reside close to the surface, and a tail region, where PVP chains extend into solution. In many of the direct MD trajectories, the deposited Ag atom bypasses the PVP tails. While such pathways have faster diffusion of the Ag atom, they also allow the atom to turn around and head back to the bulk region, as there is no gradient in the free-energy profile pulling the atom to the surface in these regions. In the PMF profiles, on which the MFPTs calculated here are based, there is a free-energy gradient that extends to the interface between the bulk solution and the near surface region for PVP20 [cf. Fig. 2(d)]. This gradient biases diffusion toward the Ag substrates and decreases the MFPT.

Considering how heterogeneity of the PVP film could influence the MFPT, it seems that a more accurate representation of the PMF as a three-dimensional function compiled over a limited sampling time would better capture the effects of long-lived heterogeneities of the PVP films. In a direct-deposition MD trajectory, the average duration of which is on the nanosecond time scale here, the Ag atom does not sample the entire film environment, as we aim to accomplish when we compute the PMF. Thus, MFPTs based on the PMF—and especially those for PVP20—do not exactly represent trends seen in direct MD simulations. Despite this, the MFPTs predicted here are within a factor of two of those from direct MD simulations. As we will discuss below, we predict similar crystal shapes to what was predicted with direct MD simulations.³⁰ 

We can compare the calculated MFPT with the MFPT that can be deduced from the experiments by Xia et al.¹⁹ From their data for Ag nanocubes grown from an edge length of 40 nm to an edge length of 130 nm over a 20-minute period in a 1.0 mM concentration of PVP (M_W = 10 000), we can obtain the MFPT using

where ΔN is the change in the number of atoms of Ag per unit time Δt, ρ_Ag is the mass density of Ag, M_W,Ag is the molecular weight of Ag, N_A is Avogadro’s number, and ΔV is the change in the cube volume. From the experimental data, this yields t_m = 9.6 × 10⁻⁶ s per Ag atom deposition. From Fig. 7 and Table III, we see that our calculated value for the MFPT is 2.73 × 10⁻⁹ s per Ag atom deposition—a factor of ∼3000 faster.

Several factors may contribute to the theoretical-experimental discrepancy. First, we estimated the concentration of Ag from the solution-phase concentration of AgNO₃, which is an upper bound—lower values for the concentration will lead to longer t_m. There is also likely free PVP in solution in the experimental studies and we did not include this effect here. Our calculations indicate that Ag interacts strongly with PVP, such that solution-phase PVP would sequester Ag atoms and hinder their diffusion. This effect would also lead to larger t_m.

Another source of the theoretical-experimental discrepancy is that the longest PVP chain considered here has 20 segments, while the shortest chain in the experimental study has approximately 100 segments.¹⁹ Although the MFPT calculations in Fig. 7 do not show a clear trend in the MFPT with increasing chain length, the results from our direct MD simulations show an increase in t_m with increasing chain length.³⁰ We discussed likely origins of this discrepancy above. An increase in the reaction time with increasing chain length goes in a direction that would reduce the theoretical-experimental discrepancy.

In the interest of coarse-graining the dynamics in the simulation of Ag nanocrystal growth, we simulate Ag atom deposition using 1-D Brownian dynamics (BD). We use the numerical scheme derived by Grassia et al.⁵⁵ and presented clearly by Ljubetič et al.⁵⁶ The numerical scheme is evaluated on a continuous length scale along the PMF profiles. Within this scheme, the position of the Ag atom is updated using

Using Eq. (5), the magnitude and direction of each displacement step at time t and position z incorporate the contribution from the PMF W(z), the diffusion coefficient D(z), and R_t(z) which is a random noise following an uncorrelated Gaussian distribution with zero mean and 2D(z) variance given by

and

We run all BD simulations at 433 K with a 100 fs time step. The reaction coordinate domain z goes from 0 to 32 Å where the atom commences its trajectory, as discussed above. We take advantage of the BD simulation’s drastically lower computational cost by simulating the Ag atom deposition process 3000 times for each system. Thus, the variability in the deposition times can be determined and the dynamics of the Ag atom deposition process can be observed. The mean deposition time can also be obtained from BD simulations and compared to the MFPTs calculated above.

Figure 9 shows six sample trajectories from BD simulations for EG-only and the EG+PVP10 system on Ag(100). These sample trajectories show Ag atom diffusion in solution and in the PVP film before it is absorbed onto the Ag surface. When the Ag atom reaches the Ag surface, we do not observe any desorption as we continue the BD simulations up to 0.5 μs. With no adsorbed PVP, the Ag atom can reflect off the free-energy barrier in the PMF [cf. Fig. 2(a)] into the bulk solvent phase multiple times before it crosses the barrier, as observed in two of the trajectories in Fig. 9(a). Table III shows the mean deposition times obtained from BD simulations and MFPT calculations using Eq. (1). For EG-only systems, the mean deposition times are 8.40 ns and 7.16 ns for Ag(100) and Ag(111), respectively. These mean deposition times are similar in magnitude with MFPTs, which are 8.21 ns and 6.98 ns for Ag(100) and Ag(111). The mean deposition times obtained from BD simulations and MFPT are very close, which verifies the veracity of Eq. (1).

When there is a PVP film, the Ag atom tends to be trapped in the film at a distance of z = 7-8 Å, as observed in the trajectories of Figure 9(b). For all of the 18 000 BD simulations with PVP-covered Ag surfaces, we do not observe any escape of the Ag atom back into the bulk solvent once it reached the local minimum prior to the hole-opening barrier. The mean deposition times for surfaces with PVP films range from 0.95 ns to 3.46 ns using BD simulations. The MFPTs obtained from Eq. (1) for surfaces with PVP films ranging from 0.93 ns to 3.75 ns are similar in magnitude, which again verifies the correctness of Eq. (1).

We observed that the distribution of t_m from BD simulations has a good empirical fit to the log-normal distribution for systems investigated, with an example shown in Fig. 10. Log-normal first-passage-time distributions were also observed for other systems with diffusion over an energy barrier, e.g., homogeneous nucleation,⁵⁷ percolation,⁵⁸ and active transport in live cells.⁵⁹ In this distribution, there is a wide spread in the deposition times. For systems with adsorbed PVP, the deposition times ranged from 0.11–24.38 ns and the ratio of standard deviation to the mean ranges from 0.59–0.89. For the EG-only systems, the deposition times can range from 0.12–67.45 ns and the ratio of standard deviation to mean is 0.97 for both Ag(100) and Ag(111) surfaces.

Kinetic Wulff shapes can be predicted from computed facet growth rates using methods described by Gadewar⁶⁰ and Zhang.⁶¹ In the framework of the kinetic Wulff construction, the steady-state shape is predicted as a function of the linear facet growth rates {G_i} of the i facets. With the assumptions that facet growth rates are limited by the solution-phase flux of Ag atoms to the facets, smooth layer-by-layer facet growth, and insignificant inter-facet diffusion, the linear growth rate of facet i is given by²⁹ 

where F_i is the molar flux of Ag to facet i. Kinetic Wulff shapes are defined by the condition that

where $x̄i$ is the relative distance of facet i from the nanocrystal center, $Ḡi$ is the relative growth rate of facet i—both are relative to a reference facet, and N is the number of facets in the crystal.

Since in the analogous experimental system,¹⁹ they observed nanocrystals with {100} and {111} facets, we consider the evolution of kinetic Wulff shapes ranging from {111}-faceted octahedra to {100}-faceted cubes. Figure 11 shows the kinetic Wulff shapes predicted by Eq. (9). We predict octahedra when G₁₁₁/G₁₀₀ ≤ 0.58, cubes when $G111/G100≥3$⁠, and shapes with a mixture of {111} and {100} facets between these limits.

Since $Gi∼tm,i−1$⁠, we have that G₁₁₁/G₁₀₀ = t_m,100/t_m,111 and we can use our calculated MFPTs to predict kinetic Wulff shapes. These shapes are shown in Fig. 12. All of the calculated growth-rate ratios from the MFPTs and BD simulations are comparable in value. From the BD simulations, the standard errors of the mean growth rate ratios are obtained from error propagation rules and are relatively small. We observe that both PVP10 and PVP20 films produce {100}-faceted cubes with no {111} corner facets. On the other hand, we find truncated octahedra with PVP5.

For the EG-only systems, our results suggest the formation of truncated cubes. However in experiments, we expect the formation of irregular nanostructures from aggregation because EG molecules do not adsorb strongly enough to Ag surfaces to act as a capping agent and prevent aggregation. The same trends discussed here are observed with the previously computed MFPT from direct MD simulations.³⁰ Our results agree with the experiments by Xia et al.,¹⁹ who find that Ag nanocubes are obtained with sufficient PVP concentration.

We investigated the atomistic details of how PVP actuates shape control in the solution-phase synthesis of Ag nanocubes. We quantitatively clarified the kinetic influence of PVP on the Ag atom deposition process by calculating the MFPTs and implementing BD simulations of Ag atom deposition. Our results show that PVP films with chain lengths of ten or greater can alter the flux towards Ag(111) to be greater than the flux towards Ag(100). PVP’s preferential binding to Ag(100) over Ag(111) gives PVP its flux-regulating capabilities due to two mechanisms: a reduced free-energy barrier for Ag atoms to cross the lower-density PVP film on Ag(111) and enhanced Ag trapping in the extended PVP film on Ag(111).

The predicted kinetic Wulff shapes are in agreement with the analogous experimental system,¹⁹ which supports our hypothesis that the kinetic influences of PVP dominate the mechanism of shape control for the seed-mediated growth of large nanocrystals. The BD simulations showed that there is a large variation in the deposition times. Our results provide direct evidence of how the preferential binding of PVP can lead to shape-controlled synthesis through kinetic influence of the atom deposition process.

The multi-scale framework presented here is especially useful for systems with large energy barriers, in which the deposition time-scale is not accessible in direct MD simulations and the problem of rare events can be tackled with umbrella sampling. We note that the average computational expenses per system of using this framework and direct MD simulations³⁰ are 450 ns and 350 ns, respectively. The computational expense for both methods is in the same order of magnitude for our system of interest. However, we anticipate that using this framework will be much more efficient than direct MD simulations for systems with large energy barriers. Our studies provide evidence that if certain precautions are followed in the umbrella-sampling procedure (i.e., obtaining averages from several short “pushing” simulations instead of one large “pulling” simulation), the multi-scale approach is equivalent to the direct MD approach, which is exact. Moreover, we demonstrate that MFPTs from BD simulations based on the PMFs and the computed diffusion-coefficient profiles are identical to those obtained from Eq. (1). The distribution of deposition times can be easily obtained in BD simulations, which is an added benefit of that approach.

Broadening our scope, the multi-scale theoretical framework we proposed here can be applied to other kinetically controlled crystal growth systems, including those in which growth is limited by diffusion and nucleation on the crystal surfaces. Since the chemical environment of the crystal system can be accounted for explicitly within MD simulations, the framework can be used to further the understanding of solution-phase nanocrystal synthesis. The capability to predict the growth morphologies in nanocrystal synthesis will allow for precise shape-control strategies, which are crucial in applications such as reaction selectivity in catalysis⁶² and optical sensing of biomolecules.^63,64

See the supplementary material for a table of system sizes, PMF, and diffusion coefficient profiles at different biasing force constants, 2D histograms of sampling in x–y space, and sensitivity analysis of diffusion coefficient smoothing.

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