Density Functional Theory Study of Pd Aggregation on a Pyridine-Terminated Self-Assembled Monolayer.

By using density functional theory calculations, the initial steps towards Pd metal cluster formation on a pyridine-terminated self-assembled monolayer (SAM) consisting of 3-(4-(pyridine-4-yl)phenyl)propane-1-thiol on an Au(1 1 1) surface are investigated. Theoretical modelling allows the investigation of structural details of the SAM surface and the metal/SAM interface at the atomic level, which is essential for elucidating the nature of Pd-SAM and Pd-Pd interactions at the liquid/solid interface and gaining insight into the mechanism of metal nucleation in the initial stage of electrodeposition. The structural flexibility of SAM molecules was studied first and the most stable conformation was identified, planar molecules in a herringbone packing, as the model for Pd adsorption. Two binding sites are found for Pd atoms on the pyridine end group of the SAM. The strong interaction between Pd atoms and pyridines illustrates the importance of SAM functionalisation in the metal nucleation process. Consistent with an energetic driving force of approximately -0.3 eV per Pd atom towards Pd aggregation suggested by static calculations, a spontaneous Pd dimerisation is observed in ab initio molecular dynamic studies of the system. Nudged elastic band calculations suggest a potential route with a low energy barrier of 0.10 eV for the Pd atom diffusion and then dimerisation on top of the SAM layer.

Introduction

Self‐assembled monolayers (SAMs), comprised of highly ordered and densely packed chemisorbed molecules,1 provide a versatile approach to modify surfaces and tune interfacial properties. Among diverse technological applications of SAMs2, 3, 4, 5 their potential as active and passive materials in electronic and spintronic devices6, 7, 8 is of great interest. However, one of the challenges in electronic applications of SAMs is the reproducible fabrication of top contacts in electrode/SAM/electrode sandwich structures. Studies on the metallisation of SAMs have shown that, in general, metal easily penetrates the monolayer and diffuses laterally at the SAM/substrate interface.9, 10, 11 In an electrochemically controlled process10, 11, 12, 13, 14 deposition starts at defects in the SAM where ions reach the substrate, that is, nucleation occurs at the SAM/substrate interface and growth proceeds from there in a mushroom‐like fashion, thus resulting in a short‐circuiting of bottom and top contacts.9, 10, 12 Therefore, a different approach is required if the advantages of electrodeposition, namely simplicity, nanoscale precision and scalability, are to be harnessed.13, 15

A strategy to restrict electrochemical metal deposition to the top of a SAM has been successfully demonstrated by Kolb and co‐workers employing a two‐step process as illustrated in Scheme 1.14 The key point is the reduction of metal ions coordinated to the SAM by electron transfer through the monolayer (step II). Upon reduction of the 2D layer of ions, the neutralised metal species diffuse at the SAM/electrolyte interface, nucleate to clusters and subsequently grow into nanoparticles. Originally demonstrated by Pd metallisation of a mercaptopyridine SAM on Au(1 1 1), the scheme has been extended to different metals (e.g., Pd, Pt, Au, Rh), different coordinating moieties (‐SH, thiazole or amino), and SAMs of different molecular architectures.14, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28 Determined by a molecular property instead of a defect in the SAM, the coordination mediated deposition scheme not only provides routine access to metal/SAM/metal sandwich structures but also offers new opportunities for metal deposition on SAMs including bimetallic systems through the combined deposition of coordinated metal and metal from the bulk electrolyte.25, 29

Scheme 1

Illustration of coordination controlled metal deposition on a SAM and molecular structures of various SAM forming pyridine thiols. The process comprises complexation of metal ions to the SAM in a non‐electrochemical environment (I) followed by their electrochemical reduction (II). Light and dark grey circles represent metal ions and metal atoms, respectively.

Even though various metals have been successfully deposited on top of different SAMs by employing the coordination controlled scheme, the underlying mechanism and the factors controlling metal nucleation and growth at the outer SAM interface remain largely unclear. This is highlighted by the morphology of the metal deposits for which both two‐dimensional monoatomic‐height islands14, 19, 22, 23, 24, 27, 28 and three‐dimensional particles a few nanometres in size18, 25, 26 have been reported, but the cause of the differences in morphology are not understood. As the timescale on which the initial stages of the deposition occur limits experiments to an a posteriori characterisation, tackling the problem from the computational side is an essential complementary approach to gain insight into the processes involved in the very early stages of the deposition process. However, the complexity of the problem involving SAMs with differences in molecular orientation and packing densities, as well as neutral and ionic metal species in an electrochemically controlled environment require significant simplifications. Leaving the electrochemical environment aside, computational studies performed by Keith and Jacob30, 31 examined diffusion, the prerequisite for nucleation. In density functional theory (DFT) calculations of single Pd atoms on a 4‐mercaptopyridine SAM, activation barriers were found to be too high for Pd diffusion. Adding a water molecule to the single Pd atom to mimic an aqueous environment did not result in a reduction of the barrier. However, in the case of pre‐adsorbed Pd atoms, the barriers decreased substantially, thus greatly facilitating lateral diffusion. Furthermore, the preference of the adsorption site changes from in‐between π‐systems to on top of the pyridine N atoms. Kučera and Groß32, 33, 34 focused on the geometric and electronic structure of the Pd/SAM/Au interface involving idealised, complete Pd monolayers. They explained the experimentally observed reduction of Pd density of states (DOS) at the Fermi level35 by strong N−Pd bonds and hydrogen adsorption on the Pd layer. In contrast, the presence of water does not significantly alter the Pd DOS at the Fermi level. These calculations were performed for a short and rigid linker, 4‐mercaptopyridine. We now extend these studies to a longer, more flexible one (see below). Special attention is called to the nature of the metal–SAM interaction and to the driving force for Pd–Pd aggregation on the SAM. We also go beyond simple static optimisations by performing first‐principles molecular dynamics simulations, which allow us to propose a possible diffusion mechanism for metal aggregation.

In DFT calculations, the initial stage of Pd cluster formation on a pyridine‐terminated SAM surface is studied. Instead of using the pyridine‐thiol (4‐mercaptopyridine, Py0) on Au (see Scheme 1), where several polymorphs have been found,36, 37, 38, 39 we consider thiols where a rigid aromatic unit is combined with a short aliphatic linker, a molecular architecture that has been established to yield well‐defined, densely packed SAMs.40, 41, 42 To reduce the computational cost, 3‐(4‐pyridine‐4‐yl‐phenyl)‐methane‐1‐thiol (PyP1) was studied, which, compared with experiments with SAMs of 3‐(4‐pyridine‐4‐yl‐phenyl)‐propane‐1‐thiol (PyP3),29, 43 has a reduced length of the alkane spacer.

Results and Discussion

SAM on Au(1 1 1) surface structure

Before addressing Pd adsorption on the PyP1 SAM, we discuss structural details of the native monolayer. Scanning tunnelling microscope (STM) studies proposed a (2√3×√3) R30° structure of the monolayer on Au(1 1 1) with a herringbone packing of the pyridine‐based molecules,43, 44 corresponding to the structurally analogous ω‐biphenyl‐alkanethiols with an odd number of methylene spacers on Au(1 1 1).41 Like biphenylthiolates,45, 46, 47, 48 PyP1 can exhibit coplanar or twisted conformations of the pyridine‐benzene moiety and the rotational freedom between the aromatic rings is a clear source for structural variations. Unfortunately, this twist angle is, in general, difficult to access experimentally. From commonly applied spectroscopic techniques like near‐edge X‐ray absorption fine structure (NEXAFS) spectroscopy, orientational information is derived from transition dipole moments (TDM, see Figure 1 a) through angles α, which are determined by the combination of two quantities, tilt (β) and twist (γ) angles. Therefore, an experimental value for the molecular tilt is based on assumptions about the inter‐ring twist angle (ω), that is, the deviation from coplanarity of the aromatic rings.43, 44 The situation is even more complicated if there is more than one molecule per unit cell. Experimentally, it cannot be differentiated whether a value arises from molecules of all the same orientation or represents the average of differently orientated molecules.49 This ambiguity highlights the need for structure calculations to complement experiments.

Figure 1

(a) Schematic drawing of adsorbed PyP1 (substrate omitted for clarity) with N, C, H and S represented by blue, grey, white and yellow balls, respectively. The inter‐ring twist angle is characterised by the angle ω. The orientation of the molecule is defined by the tilt angle (β) of the 4,4′‐molecular axis of the aromatic part with respect to the surface normal, and by the twist angles (γ) of the rings with respect to the plane spanned by the surface normal z and the 4,4′‐molecular axis. The orientation of a transition dipole moment (TDM) is expressed by the angle α, defined with respect to the surface normal. (b) Top view of PyP1 SAM on Au(1 1 1) with a (2√3×√3) R30° unit cell (red dashed line) and the equivalent (3×√3) cell used in the calculations (black solid line).

For the simulations of the PyP1 SAM, a (3×√3) unit cell containing two molecules in a herringbone arrangement (black rectangle in Figure 1 b) was chosen, as this is equivalent to the (2√3×√3) R30° unit cell, in accordance with experiments as well as calculations on related biphenylthiols.45 Both planar and twisted geometries of the adsorbed molecules were considered. Changes in tilt and inter‐ring twist angles give rise to several possible local minima, which are categorised according to the combination of conformations. The obtained geometries for the lowest energy structure in each category are shown in Figure 2, with structural details listed in Table 1.

Figure 2

Energies and top view of obtained minima structures in three different motifs: (a) planar–planar, (b) twisted–twisted, (c) twisted–planar. Only the two molecules in the unit cell are shown. The energies are given relative to the most stable planar–planar structure.

Table 1

Tilt, twist and TDM tilt angles of the lowest energy structures of the three motifs shown in Figure 2. The energies are given relative to the minimum planar–planar structure.

Geometry	Tilt angle, β [°]	Inter‐ring twist angle, ω [°]		Alpha angle of phenyl, α ph [°]	Alpha angle of pyridine, α py [°]	Relative energy [meV]
planar–planar	6.9	2.4	0.5	86.0	86.2	0
twisted–twisted	5.5	75.6	72.0	86.8	87.1	125
twisted–planar	13.0	10.3	72.1	81.9	81.4	679

As a general observation, the tilt angle of the adsorbed molecule is fairly small. The molecules are either almost planar or strongly twisted with the inter‐ring twist angle ω=72–76°, which is strongly deviated from the most stable geometry in the gas phase (ω=36°) calculated at the same level. With the tilt and the twist of aromatic rings within the molecule (for definitions of angles, see Figure 1 a), the orientation of the transition dipole moment (orthogonal to the ring) can be calculated by using the relationship50 [Eq. (1)]:

The average α values for phenyl and pyridine moieties differ only slightly even for the structures with large ω values, which is a consequence of the rather upright orientation. Both values are noticeably higher than the experimental values obtained from NEXAFS,43, 44 that is, approximately 70° for phenyl and 60–80° for pyridine. The discrepancy is explained by the fact that calculated values come from static calculations at 0 K whereas experimental values are dynamic averages at room temperature and also include molecules at grain boundaries and defects.

Polymorphs with both molecules in the unit cell of the same conformation (i.e., both planar or both twisted) are energetically more favoured as they have both aromatic rings in the preferred herringbone arrangement, where the intermolecular Coulomb repulsion is minimised. The planar–planar conformation is found to be the most stable conformation, which is in qualitative agreement with the result for biphenyl thiolate on Au,45 although the twisted–twisted one is only about 100 meV higher in energy in this case. For biphenyl thiolate, twisted–twisted forms were computed to be up to about 800 meV less stable.45 One reason for this quantitative difference may be that in the biphenyl thiolate, the S atom is directly bonded to an aryl ring, whereas in our PyP1 system it is not. The possibility of conjugation with the S atom in the biphenyl thiolate may well bestow additional stability onto the planar forms. The small energy difference for the PyP1 SAM suggests that the polymorph with strongly twisted molecules may be present under ambient conditions. Although our results suggest the possible coexistence of different polymorphs, we take the planar conformation, which has the lowest energy, as the model system for Pd adsorption.

Pd adsorption on SAM

In experiments of Pd deposition on top of SAMs,14, 18, 25, 26 PdII ions coordinated to nitrogen are electrochemically reduced. Thus, it is reasonable to assume that Pd atoms bonded to N (Figure 3 a) represent the structure immediately after metal discharge. Taking the energetically favoured structure of the PyP1 SAM with a herringbone packing of planar molecules, the structure involves two Pd atoms per unit cell, that is, a 1:1 molecule/Pd ratio. Noting that the presence of solvent molecules has been demonstrated to prevent the metal penetrating the SAM31 and also reduce the Pd diffusion barrier in Py0,30, 31 one H2O molecule per Pd was added to mimic the aqueous environment of the experiments. From model calculations on isolated, neutral [Pd(H2O)] and Pd(py)(H2O) clusters (py=pyridine, see the Supporting Information for details), we confirmed that the neutral Pd(py) moiety preferentially binds one additional water molecule (i.e., n=1) in a linear geometry, similar to the preferred coordination mode of the isoelectronic Ag+.51 Apart from the terminal N atoms, another binding site turns out to be favoured for the Pd atoms, namely in between the pyridine rings as shown in Figure 3 b. The structure with Pd atoms binding through the π‐system is energetically more stable by 0.55 eV than the one with Pd binding through N. It is noted that both binding sites are different from the most stable sites found for Py0 SAM, which is Pd binding in between two parallel rings through N. In that case, Pd and N atoms have almost the same height above the Au surface.

Figure 3

Side and top views of Pd atom adsorption on PyP1‐SAM on Au(1 1 1) through (a) N termini and (b) aromatic rings. For clarity, only the molecules of one unit cell are shown in the top view. (c) Pd adsorption energy split into several contributions: H2O–H2O interaction energy, Pd–H2O interaction energy, Pd⋅⋅⋅Pd interaction energy, SAM–Au deformation energy and Pd(H2O) adsorption energy.

To investigate the factors that drives a Pd atom to bind at a specific site, we split the energy into several contributions: H2O–H2O, Pd–H2O and Pd⋅⋅⋅Pd interactions, SAM–Au deformation (including the top Au layer) and adsorption of the Pd(H2O) complex onto the SAM. These were evaluated through single‐point energy calculations of the optimised [PdN‐SAM]2‐Au and [PdAr‐SAM]2‐Au structures, deleting atoms as appropriate. As seen from Figure 3 c, Pd adsorption is overall exothermic but involves an endothermic term originating from changes in the structure of the SAM (and the top Au layer). When the Pd atom resides on top of the SAM, a substantial exothermic contribution comes from the complexation of the N‐bound Pd atom with a water molecule (Pd–H2O). H2O–H2O and Pd⋅⋅⋅Pd interactions make minor contributions of only about 5 % each. The situation is rather different when Pd is in between the π‐rings. Besides overall reduced interactions, Pd–H2O interactions are now minor (≈10 % of interactions). Pd⋅⋅⋅Pd interactions, even though still minor, have increased (≈15 %), owing to a smaller Pd⋅⋅⋅Pd distance (4.31 Å) compared with [PdN‐SAM]2‐Au (4.68 Å, see Table S1 in the Supporting Information). H2O–H2O interactions now contribute most (≈75 %) as a result of a more embedded position of Pd, which affords a more extended H2O network on top of the SAM with favourable H2O–H2O interactions. Pd(H2O) adsorption has a destabilising influence on the SAM layer, noticeably in the [PdAr‐SAM]2‐Au structure (see SAM–Au deformation in Figure 3 b), which is not surprising because this layer is more disrupted when Pd atoms bind through the π system and sit 1.3 Å below the N plane compared with the on‐top configuration where Pd is located 1.6 Å above (see Table S1 in the Supporting Information). However, this destabilisation is overcompensated by the strong interaction between Pd(H2O) and the SAM layer, which also accounts for the most pronounced energy difference between the two binding modes (≈2.1 eV). Arising from relatively strong Pd–π interactions, there is a clear preference for Pd atoms to be located in between the pyridine rings instead of sitting on top of the SAM and interacting via N. In fact, the difference in the Pd(H2O) adsorption energy is so large that it is the main factor determining the relative stability of the two binding modes. This indicates that the Pd–SAM interaction is the most important contribution determining the overall stability. Its substantial variation in energy as function of position of the Pd atoms should be a critical factor in the formation of clusters and, thus, the overall metallisation process as it will both affect the kinetics of Pd diffusion and define nucleation sites. According to the semiempirical dispersion corrections, dispersion interactions contribute noticeably, but not decisively to this Pd–SAM interaction (between ca. 0.4 eV to 0.8 eV, or ca. 10 %, see Table S3 in the Supporting Information). It should be noted, however, that the combined Pd–H2O and H2O–H2O interactions are substantial contributions. We note that direct comparison of these results with those of Keith and Jacob30, 31 is difficult because of the different nature of the substrates (Py0 vs. PyP1) and extent of solvation (most of Keith and Jacob's calculations were performed without water). The calculations agree that Pd0 should preferentially bind in between the pyridine rings rather than on top of the N atoms, but the precise binding site and interaction strength may depend on the system under scrutiny.

The formation of larger metal clusters arguably requires diffusion of individual Pd atoms (or smaller clusters) on top of the SAM surface. It is quite likely that N‐bound Pd atoms (or smaller clusters, see below) will be involved in this process. The energy difference between Pd atoms adsorbed via either N atoms or π‐systems of the pyridine could be critically important for the kinetics of diffusion of metal atoms towards cluster formation after discharge. Based on the above calculations, which show a clear preference for Pd intercalation (ΔE=−0.55 eV), diffusion of the metal atoms on the PyP1 SAM seems a rather hindered process, (although even higher barriers between 0.7 eV and 1.6 eV have been reported for single Pd atoms, free or hydrated, on Py0 SAM30, 31). In these structures, Pd atoms are well separated and, therefore, have little interaction (Figure 3 c). However, towards the formation of dimers and, subsequently, cluster configurations have to be considered where Pd⋅⋅⋅Pd interactions become substantial. We optimised a number of structures with different locations of the Pd atoms involving Pd⋅⋅⋅Pd contacts (i.e., smaller than sum of the van der Waal radii, which is 3.26 Å). The focus was on arrangements consisting of a mixture of atoms intercalated (located between the pyridine rings below the N‐plane) and on top of the SAM (well above the N‐plane), that is, a combination of sites, which, in the case of non‐interacting Pd, are energetically unfavourable and preferred, respectively. Figure 4 shows a compilation of the structures. The breakdown of Pd adsorption energy for these local minima can be found in the Supporting Information (Figure S2 and Table S3). For comparison, the initial structure with Pd on top of N, representing the earliest stage of metal deposition, is depicted in Figure 4 a (structure I, identical to the one in Figure 3 a). If one Pd atom is placed in between the pyridine rings, the other one remaining at the N‐terminus (Figure 4 a, II), the optimised structure yields a Pd⋅⋅⋅Pd distance of 3.69 Å and an energy lower by 0.57 eV compared with the on‐top configuration. Notably, this structure is very similar in energy to that of Figure 3 b where both Pd sit between the pyridine rings. We located two additional minima with shorter Pd⋅⋅⋅Pd distances. One of these (Figure 4 a, III), where this shortening of the Pd distance is associated with a lateral displacement of the on‐top Pd and upward movement of the intercalated Pd is more stable than structure II by about 70 meV (ΔE=−0.64 eV). The corresponding Pd–Pd bond length of 2.89 Å is close to the one of the bulk metal (2.75 Å), but significantly longer than in an isolated Pd2 cluster (2.53 Å) calculated at the same level. The structures suggest a substantial energetic driving force for dimer formation of approximately −0.3 eV per Pd atom with rather different heights of the atoms. A substantially different structure is the one depicted in Figure 4 a, IV, which is characterised by both atoms intercalated. With a Pd–Pd distance of 2.56 Å, they reside at the bottom and top of the rings, respectively. Notably, despite the quite different positions of the metal atoms the energy of the structure is the same as in structure II. The conclusions drawn from these structures are that there is a substantial energetic driving force for dimer formation and the energy landscape involves a multitude of configurations, which are very similar in energy but differ substantially in position, bond length and orientation.

Figure 4

(a) Side and top views of optimised [Pd(H2O)‐SAM]2‐Au local minima. For clarity, only the molecules of one unit cell are shown in the top view. The zoom‐in images showing the coordination of water molecules can be found in Figure S3 in the Supporting Information. (b) Local density of states (LDOS) of the two Pd atoms in [Pd(H2O)‐SAM]2‐Au. PdN refers to Pd atoms bound to nitrogen atom of pyridine in the SAM molecule. PdAr refers to Pd atoms bound to the aromatic system of the pyridine ring.

To elucidate the nature of Pd–SAM and Pd–Pd interactions, we determined the local density of states (LDOS) of Pd. As shown in Figure 4, for conformations with isolated Pd atoms, the LDOS of PdN (the Pd atom binding via N) are quite similar. The PdAr (the Pd atoms binding to the π system of the pyridine rings) exhibit a significant reduction of the LDOS near the Fermi energy (Figure 4 b). The maximum is shifted to about −1.8 eV below the Fermi level. The Pd LDOS is further altered when a Pd–Pd contact is formed. Such a large change of the LDOS when Pd changes its binding site indicates a substantial charge transfer between the involved constituents. Therefore, Bader analysis was performed to get a more quantitative view of the charge distribution within the Pd/SAM/Au interface. To estimate the influence of Pd(H2O) adsorption, we compare the results for systems with and without Pd(H2O) adsorption. The results are summarised in Table 2. In the bare SAM–Au system without Pd(H2O) adsorption, the SAM molecule is negatively charged. The electrons are transferred from Au to the SAM molecule through the S−Au bond. The SAM molecule gets more negatively charged upon Pd(H2O) adsorption on top of the nitrogen (Figure 4 a, I), which suggests an additional electron transfer from Pd(H2O) to the SAM molecule. When one Pd atom diffuses into between the pyridine rings (II), the accumulation of electrons at the SAM molecule increases. It further increases when a Pd dimer forms (III, IV). Interestingly, the electron depletion at Au remains almost unaltered upon Pd(H2O) adsorption and the change of Pd binding site. Also, a closer inspection of the charge distribution within the SAM layer reveals that the negative charge accumulation upon Pd(H2O) adsorption is mainly on the pyridine ring (see the charge of the pyridine ring listed in parentheses in Table 2). The charge distribution in corresponding Pd(py)(H2O) model systems (see Table S2 in the Supporting Information) is remarkably similar to those in the SAM minima in Table 2. This means the S–Au interaction, the aliphatic linker and even the lower phenyl ring are hardly affected by the change in Pd(H2O) adsorption. This result highlights the dominating role of the SAM terminating pyridine group and suggests studying the Pd nucleation process without taking the Au substrate and whole SAM molecule into account. We will exploit this finding in the following to study suitably truncated models.

Table 2

Atomic charges from Bader analysis of the total electron density for SAM–Au and [Pd(H2O)‐SAM]2‐Au local minima in Figure 4. The charges, given in units of the elementary charge (e), are the average values for the two Pd atoms, the two H2O molecules and the two SAM molecules, as well as the sum of all Au atoms within one unit cell.

	SAM–Au	[Pd‐SAM]2‐Au
		I	II	III	IV
Pd		0.060	0.194	0.220	0.241
H2O		0.023	0.019	0.020	0.017
SAM (pyridine)	−0.095 (−0.006)	−0.181 (−0.063)	−0.303 (−0.190)	−0.334 (−0.217)	−0.357 (−0.229)
Au	0.189	0.196	0.181	0.188	0.198

Molecular dynamics simulation of the initial Pd nucleation process

The above calculations on static structures at 0 K suggests an intrinsic preference for Pd dimers over structures with isolated atoms on top of the SAM. However, with experiments typically performed at room temperature, it is necessary to take temperature effects into account. Therefore, we employed ab initio molecular dynamics (MD) simulations to model Pd dimerisation on the pyridine‐terminated SAM. Based on the charge distributions discussed in the preceding section, the SAM was simplified to pyridine molecules. Another simplification is the omission of the substrate. To mimic attachment to a surface, the z coordinates of the anchoring C atoms (in the para position from the N atom) were fixed. With this setup, the pyridine molecules are free to move in a plane, but not to move away from it, that is, detach from a surface. As documented in section 1 of the Supporting Information, static geometry optimisations of this truncated model afford qualitatively similar structures and relative energies as the full SAM models discussed above.

Two MD trajectories were generated: one starting from a structure with two N‐binding Pd atoms (representing the scenario directly after discharge, structure I in Figure 4), and another one starting from one Pd atom intercalated between two pyridine rings and one N‐binding Pd atom (corresponding to structure II in Figure 4). The MD trajectory starting from two Pd atoms intercalated between the pyridine rings (corresponding to the structure in Figure 3 b) was also generated but is not discussed here because the structure remained unreactive in terms of Pd⋅⋅⋅Pd contact formation (see the Supporting Information for details). The time evolution of the Pd⋅⋅⋅Pd distance is plotted for both simulations in Figure 5, along with selected snapshots from the trajectories. In both cases, spontaneous Pd dimerisation is observed after approximately 10 ps. For the starting structure with two N‐binding Pd atoms, the two Pd atoms are well separated initially with a distance of approximately 4.4 Å. The Pd–N interaction is fairly strong, so that only small oscillations in the Pd⋅⋅⋅Pd distance are seen in the initial 8 ps. At approximately 10 ps, a Pd2 dimer is formed, indicated by a sharp decrease in the Pd⋅⋅⋅Pd distance. Upon Pd–Pd contact formation, the electronic structure of Pd changes; one Pd atom loses its coordinative complexed water and starts to interact with the π system of the pyridine ring (see the left of the two Pd atoms depicted in Figure 5 a). The mean Pd−Pd bond length for the remainder of the simulation is 2.65(9) Å, which is very close to the bond length of an isolated Pd2 cluster (i.e., 2.56(9) Å) calculated for pristine Pd2 by using the same box dimensions. Although results from the static model suggest that the structure with a Pd2 dimer intercalated between the pyridine rings (structure IV in Figure 4) is close in energy to the structure with a Pd2 formed on top of the rings (structure III in Figure 4), such a metal penetration is not observed on the timescale of our simulation.

Figure 5

Pd⋅⋅⋅Pd distances (in Å) and selected snapshots (side views) along two MD trajectories of truncated pyridine‐terminated SAM models starting from (a) a configuration with two Pd atoms bound to N and (b) a configuration with one Pd atom intercalated between the aromatic rings and the other one bound to N.

The product dimer structure formed in this trajectory (shown at 14.4 ps in Figure 5 a) was subject to geometry optimisation. In this minimum, as in the MD snapshot, one Pd(H2O) is bound to a N atom and the other Pd loses its H2O ligand and starts to interact with the π system of the pyridine ring. At the PBE‐D2 level it is −0.12 eV lower than the starting minimum (corresponding to structure I), see Figure S4 in the Supporting Information (where this product minimum is labelled V). This structure is thus confirmed as a viable intermediate en route to structures with more deeply intercalated Pd dimers (such as II and IV in Figure 4 a), and, possibly, to much larger clusters that cannot penetrate into the SAM layers anymore.

For the other trajectory, the simulation started with one Pd atom adsorbed in between pyridine rings and another one binding to N. With one vacant N site, the overall mobility of the N‐bound Pd atom seems to increase. After only 3 ps, this Pd atom starts to move from its end‐on position along the nitrogen lone pair into a more side‐on, perpendicular position approaching a neighbouring pyridine ring. This diffusion causes a rearrangement of the model SAM, in this case the pyridine layer, which pushes the intercalated Pd slightly “downwards” (i.e., away from the N‐terminus), but it still interacts strongly with the pyridine rings. This rearrangement gives rise to an initial increase in the Pd⋅⋅⋅Pd distance (to ca. 5 Å). After a total simulation time of approximately 10 ps, the initially N‐bound Pd atom further diffuses in between the pyridine rings and forms a Pd–Pd contact. The Pd−Pd bond is now 2.99 ű0.19 Å, somewhat longer than the bond length of the Pd2 dimer in the previous case (cf. Figure 5 b). Also, the oscillation of bond length is slightly larger, indicating that there is a competition between Pd–Pd and Pd–pyridine interactions.

The spontaneous formation of Pd–Pd contacts in unconstrained MD simulations is remarkable, because it implies very little kinetic hindrance for this process. To explore the energetics of this Pd diffusion and aggregation, we analysed the two dimerisation pathways seen in MD simulations by using the nudged elastic band (NEB) method. The calculation results reveal that the barrier for the aggregation of two N‐binding Pd is 0.78 eV (Figure 6 a). A lower barrier is computed for the dimerisation when there is one Pd atom intercalated into between the pyridine rings (Figure 6 b). The overall process was divided into two segments, i) the Pd diffusion from the on‐top‐N site into between the pyridine rings, and ii) the Pd2 dimer formation. The energetic barrier for the Pd intercalation is approximately 70 meV, which agrees with the previous report of very low Pd diffusion barrier for the Py0‐Au system with a pre‐adsorbed Pd atom in between the rings,30 although in our case, the initial intercalated Pd atom is below the N plane. In the second segment, where the dimerisation occurs, a higher barrier of 0.25 eV is computed. The overall barrier for this Pd diffusion and aggregation process is 0.10 eV. This result suggests that the Pd intercalation may play an important role in the Pd aggregation at the earliest stage of the deposition. We note that the NEB‐derived barriers for this second path (Figure 6 b) are indeed very low, consistent with very fast processes, whereas that for the first path (Figure 6 a) is too high to be compatible with the spontaneous process seen in the MD. It may well be that the simple NEB approach, which can only give an upper limit to the barrier on the actual minimum energy reaction path (MERP), failed to locate the MERP in this case. The conformational flexibility of an MD simulation may be required to find this. For a more reliable and quantitative appraisal of the actual activation barriers involved, more sophisticated simulation methods (such as constrained free‐energy MD simulations) would be needed. In any event, the energetic driving force for Pd−Pd bond formation apparent in the static optimisations agrees well with the spontaneous dimer formation in the MD simulations of the truncated model. There is thus a notable and noticeable driving force for Pd aggregation at the earliest stage of the deposition, immediately after PdII reduction.

Figure 6

Energy diagram for Pd2 dimer formation starting from (a) a configuration with two Pd atoms bound to N and (b) a configuration with one Pd atom intercalated between the aromatic rings and the other one bound to N. The initial state, transition state (TS), intermediate and final states are shown.

Conclusion

We have applied DFT calculations to study the initial stage of Pd metallisation of a (2√3×√3) R30° pyridine‐terminated SAM on Au(1 1 1). Considering the structural flexibility of the PyP1 molecule, the most stable configuration in the respective SAM was identified to consist of planar molecules arranged in a herringbone packing. Using this packing as a starting structure for the Pd adsorption studies, various binding sites were found for Pd at the pyridine end group above and in between the SAMs. The Pd atom prefers to reside in between the aromatic rings rather than on top of the ring. The strong interaction between Pd and pyridine highlights the importance of SAM functionalisation in metal nucleation processes. The dimerisation of Pd is favourable on the SAM surface with an energetic driving force of approximately −0.3 eV per Pd atom. Consistent with this driving force, and indicative of low kinetic barriers, spontaneous Pd dimerisation is seen in ab initio MD studies of the system. By using ab initio MD and NEB calculations, we have identified a dimerisation path with a small energy barrier of 0.10 eV for the system with one pre‐adsorbed Pd atom.

These results serve to shed some light on the structural preferences and possible reaction channels of the aqueous SAM–PdII system immediately after metal ion coordination and discharge. There is a clear energetic driving force for dimerisation of two neighbouring Pd atoms and, to a lesser extent, for the intercalation of the resulting dimer in between the aromatic rings. From the absence of noticeable kinetic hindrance for this dimerisation process, it appears that the formation of larger metal clusters may not be limited by the initial nucleation process to form small clusters, but rather by the diffusion of Pd atoms (or small clusters) along the SAM surface to form larger aggregates. These diffusion barriers in turn should be related to the energetic stabilisation of small clusters intercalated between the aromatic end groups of the SAM molecules. Further calculations involving larger unit cells and larger Pd clusters are in progress to investigate this possibility.

Experimental Section

DFT calculations on SAM metallisation were performed by using the plane‐wave pseudopotential method as implemented in CASTEP,52 version 16.11. The electronic exchange and correlation were treated by using the generalised gradient approximation (GGA) formalised by Perdew, Burke and Ernzerhof53 (PBE). The long‐range dispersion forces were considered by using Tkatchenko–Scheffler's correction scheme.54 All the calculations were performed as closed shell as the test spin‐polarised calculation (on the [PdN‐SAM]2‐Au system) converged to spin 0. On‐the‐fly (OTF) pseudopotentials generated by the CASTEP software package with a kinetic cut‐off energy of 750 eV were used, which ensured a convergence of total energies to less than 1 meV per atom. Brillouin zone sampling was obtained by using a 3×7×1 Monkhorst–Pack55 k‐point grid. To model the metal nucleation on a SAM on Au substrate, we applied the repeated slab approach. A vacuum region of at least 20 Å was added to the cell to suppress interactions between repeated slabs, where a self‐consistent dipole correction56 was also applied to account for electrostatic asymmetry. In this study, the Au surface was represented by a three‐layer slab. The Au(1 1 1) nearest‐neighbour distance was set to the equilibrium value (2.907 Å) obtained from bulk metal calculation by using the same methodology. The atomic coordinates of the adsorbed molecules and the top layer of the slab were relaxed by using the Broyden, Fletcher, Goldfarb and Shannon's (BFGS) method57 (convergence criteria: maximum change in system energy=2×10−5 eV, maximum root‐mean‐square (RMS) force=0.05 eV Å−1 and maximum RMS displacement=0.002 Å). The local density of states (LDOS) was calculated with the OptaDOS code58 by using the adaptive broadening Scheme. Bader atomic charges were obtained from topological analysis of the total electron density59, 60, 61, 62 for the optimised structures. For gas‐phase calculations of isolated SAM molecules and Pd2, a 40 Å×40 Å×40 Å unit cell and a 15 Å×15 Å×15 Å unit cell were used, respectively, with a k‐point grid of 1×1×1.

Ab initio molecular dynamics (MD) simulations and nudged elastic band (NEB) calculations were performed by using the CASTEP package for a truncated‐SAM model (pyridine molecules). A cut‐off energy of 600 eV for MD simulations and 700 eV for NEB calculations was used to construct the plane‐wave basis, with OTF ultrasoft pseudopotentials and the PBE exchange‐correlation functional.53 The van der Waals forces were considered by using Grimme's dispersion correction scheme of 200663 (DFT‐D2). K‐point sampling was performed by using a Monkhorst–Pack55 grid of 3×1×1 for MD simulations and a 7×3×1 grid for NEB calculations. The reduced model size, cut‐off and k‐point grid are designed to save CPU time in the expensive MD simulations (as documented in section 1 of the Supporting Information, these truncations do not affect the qualitative findings such as relative energies of isomers). During each SCF cycle, the electronic structure was minimised to a tolerance of 2×10−6 eV. The equations of motion were integrated based on the velocity Verlet algorithm with a time step of 1.0 fs. MD simulations were performed at 298 K, and controlled by a Langevin thermostat and a canonical ensemble (NVT). The NEB method64 was used to calculate the activation barriers of Pd diffusion and dimerisation. For NEB calculations, the initial and final states were constructed based on the aggregation paths observed in the MD trajectories.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.