Immobilization of Au nanoparticles on graphite tunnels through nanocapillarity

Abstract

Atomistic computer simulations on the generation of nanotunnels on graphite and the subsequent immobilization of gold nanoparticles are presented in this work. A Morse potential dependent on the coordination of carbon atoms was parameterized based on density functional theory including long dispersion forces. The set up chosen is such that a direct comparison with the experiments is possible. The model is able to reproduce crucial experimental aspects such as the phenomena of capillarity and the final height of the immobilized nanoparticle. Results presented here can inspire the design of new platforms for protein immobilizations.

---

1 Introduction

Immobilization of nanostructures on selective surfaces is an important technical issue for several reasons.¹ For the design of polymer nano-composites, polymers are filled with particles in order to improve their mechanical, optical, and chemical affinity properties, and it is required for the particles to be attached strongly to the substrate in such a way that agglomeration is avoided and the nano-composite is kept dispersed on the surface.² Iron magnetic nanoparticles have been used to immobilize proteins acting as active agents in enzyme-catalyzed reactions,³ and thiolated gold nanoparticles have been kept fix in a matrix using multidentate imidazole surface ligands, with the aim of using the particles as a platform for bioconjugation with oligonucleotides and peptides in biosensing applications.⁴ In the design of nano-electronic circuits, it is a requirement to keep nanoparticles (NP) strongly attached to the surface even if the system is subjected to mechanical stress, thermal fluctuations, or external electric fields; in order to reach this mechanical stability, Au nanoparticles have been immobilized in SiO₂ surfaces by embedding them in octadecylsiloxane⁵ or by covering them with amines and thiols.⁶ Another possible way of immobilizing naked metal nanoparticles is to deposit them onto a porous polymer support or on polymer brushes. Nevertheless, in applications such as catalysis it is desirable to have the surface of the metal particles exposed to the reacting agents. In this sense, these conventional methods that cover the nanoparticle with elongated molecules are not an efficient route for immobilization.

Palmer et al.⁷ proposed a novel methodology which consists in bombarding small Au clusters on a highly oriented pyrolytic graphite (HOPG) surface. The clusters are generated by means of radio frequency magnetron sputtering and a gas condensation cluster source, combined with a lateral time-of-flight (TOF) mass filter. These clusters arrive to the graphite surface with enough energy (1.5 keV) to penetrate the first carbon layers, creating “nano-tunnels” (holes).

The number of graphite layers penetrated by the Au cluster and the amount of damage, may depend not only in the kinetic energy, but also on the shape and the relative orientation of the cluster with respect to the surface. In their experiments, Palmer⁷ used clusters of 20 atoms in size (Au₂₀). As it is very well known, such a cluster has a tetrahedral structure, which is a fragment of the face-centered cubic lattice of bulk gold with a small structural relaxation.^8,9

Once the tunnels are created, larger Au nanoparticles, around 1 to 2 nm in size, are projected onto the graphite surface with a kinetic energy corresponding to a soft-landing regime (i.e. less than 1 eV per atom). The deposition at this soft-landing regime, preserve not only the size but also the general shape of the nanoparticles.¹ A particle that lands in the vicinity of a tunnel gets immobilized by partially inserting itself into the tunnel, a phenomenon known as “nanocapillarity”. The specific details of the immobilization process in the experiment were not explicitly known, since the resolution of the measurements made by atomic force microscopy (AFM) is not enough to investigate the final shape of the particles nor the amount of volume of the particles that got trapped into the tunnels.

In this work, we present an investigation based on Molecular Dynamics (MD) where we simulate both the generation of nano-tunnels in graphite, and the immobilization of Au nanoparticles. In order to do so, it is essential to count with a potential that correctly describes the gold–carbon interaction. Despite being the most used potential for gold and graphite, the Lennard-Jones fitting of Lewis et al.,¹⁰ severely overestimates the adsorption energy, and underestimates the equilibrium height of gold clusters on such a surface. In fact, this fitting was “determined rather loosely from various two-body models for Ag–C and Pt–C interactions” and it was expected to “provide a qualitatively correct description of the system”.¹⁰ Nevertheless, several simulations were performed using this potential, and quantitative descriptions were made despite the limitations of the model.^11,12 Another limitation of this potential is that it does not distinguish between carbon atoms at tunnel edges and at pristine graphene, although their reactivities are radically different. For all these reasons, a new parametrization of a Morse potential based on Density Functional Density (DFT) considering dispersion forces was developed. In the next sections, the procedure of this parametrization as well as the settings of the MD simulations are described.

2 Computational details

2.1 Interaction models

The interaction between Au atoms is modelled using the Embedded Atom Model (EAM), due to Daw and Baskes.¹³ The EAM is a many-body potential based on DFT calculations, and it considers a background electronic density and the electrostatic interaction due to the overlap of atomic nuclei. The carbon–carbon interaction is modelled using the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO), due to Stuart,¹⁴ developed for hydrocarbons. This potential function considers non-bonded interactions, and it includes terms for dihedrals and torsion angles.

Ab initio calculations were performed using the Quantum Espresso/PWSCF code.¹⁵ Vanderbilt ultrasoft pseudopotentials¹⁶ were used with the Perdew–Burke–Ernzerhof¹⁷ (PBE) and the local density approximation¹⁸ (LDA). A 36 Ry kinetic energy cutoff and a 300 Ry charge density cutoff was used. The reciprocal space of a (1 × 1 × 1) supercell was sampled with a (12 × 21 × 1) k-point grid using the Monkhorst–Pack method¹⁹ and reduced inversely with the supercell size. All these variables were carefully parameterized. A Gaussian broadening of 0.01 Ry was applied. The convergence of the forces for DFT relaxations are equal to 0.02 eV Å⁻¹. The supercells were constructed in such a way that the periodic images of the NP were at least 10 Å from each other. There is a vast spectrum of DFT methods that attempt to describe van der Waals (vdW) interactions. We performed several tests with different vdW-DF functionals. In each case a notorious increase of the adsorption energy was found with respect to the PBE functional. The calculations shown in this work correspond to the vdW-DF2-c09x functional^20,21 as it is the most adequate choice due to its precision in both long and short range interactions and its fast convergence.^22,23

2.2 C–Au potential parametrization

In order to successfully simulate the experiments reported by R. Palmer,⁷ two processes have to be correctly described by the potential chosen: the physisorption and the chemisorption of gold nanoparticles on graphite. In the first case the cluster is interacting with a pristine layer of carbon atoms, while in the second case the nanoparticle is interacting with the carbon edges of the tunnel. The criterion chosen to distinguish these two scenarios is the coordination number of carbon atoms (excluding gold). In the case of physisorption, the NP interacts with pristine graphene, where each carbon atom has a coordination of three. In the case of chemisorption, the NP interacts with the tunnel edge, where the coordination of carbon atoms with dangling bonds can be either one or two. Consequently, three Morse potentials for these three coordination numbers were built, conforming what we called the Coordination Dependent Morse Potential (CDMP).

Let us describe the physisorption first. As a first approximation the gold–carbon interaction was obtained from the interaction between graphene (not graphite) and gold clusters (not nanoparticles) of increasing size, using DFT. Such clusters consisted in two-layer structures showing the (111) face to the carbon surface (see Fig. 1a). The nomenclature for metal clusters is Au_m−n, where m is the number of atoms facing the graphene surface and n is the total number of cluster atoms. The difference between the adsorption on graphite and on graphene was found to be negligible, which allow us to choose the latter for the rest of the DFT calculations. After full relaxation of the interacting parts, single point calculations were performed for clusters at several heights from graphene. The corresponding physisorption curves are shown in Fig. 1b. Note that the equilibrium energy becomes less negative with increasing cluster size and tend to a plateau. In fact, the adsorption energy becomes invariant at clusters larger than Au₆₋₁₆. This should hold also for nanoparticles. The reason for this behavior can be attributed to the bond order conservation. As the clusters grow in size the metal atoms become more coordinated with each other, being less available to the carbon surface. Since nanoparticles from 147 atoms were employed in the experiments, we focused on the two largest clusters, where the equilibrium energy is already converged. These outcomes were used to parametrize the physisorption of gold clusters on graphite. The resulting parameters are in Table 1 (for the case of n_C = 3).

Fig. 1 (a) Bottom view of two of the five clusters used to study physisorption (Au₃₋₉ and Au₇₋₁₉). Metal atoms facing the surface and further from it are colored yellow and orange, respectively. (b) Physisorption energy (per metal atom facing the carbon surface) for five clusters on graphene at several heights. At least one metal atom is in a top site. Points correspond to DFT values while lines were obtained by the CDMP. (c) Top site view of structures used to study chemisorption. Atoms are colored yellow for gold, cyan for carbon, white for hydrogen, and red for carbon atoms involved in the chemical bond with gold. (d) Chemisorption energy (per C–Au bond) of a gold cluster (Au₇₋₁₉) at several distances from the edge of graphene. Points correspond to DFT values while lines were obtained by the CDMP. To ease the comparison, values of physisorption of the same cluster on the carbon surface are also shown.

Table 1 Parameters of the coordination dependent Morse potential for the gold–carbon interaction: coordination number (n_C), well depth (D_e), equilibrium bond distance (r_e), and well width (α)

nC	De (eV)	re (Å)	α
1	1.38185	1.95452	2.25875
2	1.46368	2.07302	1.93911
3	0.00832	3.87035	1.25707

---

Chemisorption processes were simulated by the adsorption of cluster Au₇₋₁₉ at the edges of a graphene flake, shown in Fig. 1c. These edges are planar structures closed with hydrogen atoms and open at one side to form chemical bonds with the gold structure. The two most common edges observed in our simulations of graphite tunnels are zigzag (with coordination number n_C = 2) and Klein kink²⁴ (with n_C = 1) edges. We chose these two to perform the DFT study. In this case, the structures were relaxed by keeping the second gold layer (further from the edge) fixed. Then, single point calculations were performed for the cluster and the graphene edge at several distances from each other (frozen at the configuration of minimum energy). The resulting chemisorption curves are shown in Fig. 1d, where they are compared with the physisorption curve of the same cluster on graphene. The equilibrium adsorption energies for both edges are almost equal and one order of magnitude larger than that of the physisorption (an average of −1.42 eV for chemisorption and −0.17 eV for physisorption). However, at distances larger than 3.6 Å physisorption becomes dominant, manifesting the typical behavior of long range interactions. These DFT data were used to parametrize the cases of coordination one and two, resulting in the corresponding values shown in Table 1.

3 MD simulations

We performed two kinds of molecular dynamics simulations using the LAMMPS code,²⁵ an open source program designed to simulate a large variety of systems at atomistic resolution. For all the sets of simulations made in this study we employed a time step of 1 fs, and the canonical simulations were made using a Nosé–Hoover thermostat with a damping parameter of 100 ps. The graphite plates were thermalized with a short NVT run of 1000 steps at 300 K, previous to any interaction with the metal clusters. In the starting stage, we simulated the implantation of Au₂₀ clusters in graphite and studied the damage provoked at the surface. The implantation of Au₂₀ clusters was made at several kinetic energies, covering a range from 0.5 to 2.5 keV in 0.25 keV intervals (in the experiments of Palmer the energy was set at 1.5 keV). The reason to investigate this range was to find the minimum energy required to make enough damage for the creation of the tunnels in the surface, and to compare the amount of damage as function of energy and relative orientation of the cluster. We performed 10 independent runs for each case, varying only the initial distribution of velocities at each run, assigning the velocities through a random Gaussian distribution at 300 K. Since the Au₂₀ has pyramidal shape, three high-symmetry orientations were used in the simulations: exposing a face, an edge, and a vertex towards the graphene surface. These orientations are represented in Fig. 2.

Fig. 2 Relative orientations of Au₂₀ with respect to the graphite surface (below the particles, not shown). (a) Face; (b) edge; (c) vertex.

The second stage consisted in the deposition of a Au icosahedral particle in the treated surface, in order to promote the entrapment of the particle by the tunnel created by the Au₂₀ cluster.

3.1 Au₂₀ implantation

The implantation of Au₂₀ clusters was simulated with MD. The isolated Au₂₀ clusters governed by the EAM interactions were thermalized with a short microcanonical run at 300 K of 200 steps previous to their exposure to the graphite surface. Once the system was thermalized, each implantation run was of 10 000 time steps in the microcanonical ensemble with a starting temperature of 300 K, for a total simulation time of 10 ps, time enough for the implantation, since the Au₂₀ cluster was originally positioned approximately 10 Å above the graphite surface. The simulation box consisted in a 12-layer graphite surface with dimensions 102.38 Å × 100.98 Å. At each of the 11 layers closer to the cluster, carbon atoms as far as 49 Å from the impact zone were kept frozen in space in order to avoid translation of the whole structure as effect of momentum transfer due to the impact of the gold cluster. In the last layer, carbon atoms as far as 10 Å from the impact zone were kept frozen. We found that already at 300 K, the Au₂₀ cluster loses its pyramidal shape, although according to ab initio MD it should melt at a temperature of around 800 K.²⁶ It is not surprising that EAM can not describe the stability, angularity, and relativistic effects of Au₂₀, since in such a small cluster these features are strongly influenced by quantum effects. To deal with this issue, we repeated the set of simulations using a completely rigid Au₂₀ pyramid, and compared the corresponding surface damage with that of the soft cluster. The MD runs for these rigid clusters were 20 500 steps long, for a total of 20.5 ps of simulated dynamics, since we found that energy stabilization took longer in these clusters compared against the EAM clusters. We found that the EAM cluster makes a tunnel that reaches the fifth or sixth carbon layer, while with the rigid cluster the tunnel reaches the twelfth carbon layer. It is reasonable to assume that the damage done in the implantation experiments should be between these two depths. The surface damage was characterized by measuring the area of the hole made in the first graphite layer. Fig. 3 shows the surface damage as a function of the kinetic energy of Au₂₀. These values correspond to the average of ten simulation runs for each cluster speed. As it can be observed in Fig. 3 the surface damage shows an hyperbolic growth with respect to the kinetic energy. It is also observed that Au₂₀ clusters with the face orientation are consistently the ones that produce the smaller surface damage. Apparently, the larger contact area distributes the strain on the surface more evenly than the other two orientations, so that the momentum is transferred to the surface in a less localized manner. In contrast, the other two orientations concentrate the momentum transfer in a very localized area directly below them, which explains the larger amount of damage. When Au₂₀ was considered as a rigid body, the damage was significantly higher than in the previous case, as shown in Fig. 3. Again, we expect the experimental value to be between these two situations. In this case, it is the edge orientation the one associated with smaller surface damage.

Fig. 3 Surface damage for the implantation of soft Au₂₀ on graphite, as a function of the kinetic energy at three different orientations. Results using rigid Au₂₀ at 1.5 keV are also shown for comparison purposes.

MD simulations also revealed other observables. At kinetic energies lower than 0.5 eV, in most cases soft Au₂₀ does not penetrate into the graphite surface, it bounces back instead. At kinetic energies higher than 0.5 eV the soft Au₂₀ gets partially atomized after the collision due to the associated high momenta. Consequently, some carbon and gold atoms escape from the surface, while others remain in the tunnel. It is important to remark that these simulations can help in the design of the experimental setup, in the sense that a broad range of impact energies are explored.

3.2 Immobilization of Au nanoparticles on graphite tunnels

After simulating the impacts, we extracted the first six graphite layers and minimize their energy using the usual conjugate gradient minimization process. Following the experimental setup of Palmer et al.⁷ the same gold nanoparticles were considered for the simulations, namely Au₁₄₇, Au₃₀₉, Au₅₆₁, and Au₉₂₃. As it is very well known, these magic numbers correspond to completed onion-like geometric shells of icosahedral symmetry. The selection of these nanoparticles and the model employed allow us to make a direct comparison between theoretical simulations and experiments. The whole system was thermalized with an NVT run for 1000 time steps at 300 K. While in the experiments by Palmer et al. the soft landing of nanoparticles occurred randomly on the carbon surface, we initiated the NVT MD simulations with the particle already on graphite at the vicinity of the tunnel and moved the particle towards it with an initial velocity of 0.5 Å ps⁻¹ (this trick was used in order to accelerate the overall process). The temperature was kept at 300 K with a Nosé–Hoover thermostat. The gold nanoparticles are oriented with one (111) face resting on the carbon surface. Once the particle is immobilized by the tunnel, at the end of each run we used the last configuration to measure the degree of insertion of the particle and its maximum height with respect to the carbon surface. 40 MD simulations (10 for each particle) were made in the canonical ensemble (NVT), with the temperature fixed at 300 K and a time step of 1 fs, and the tunnels prepared as described in the previous section. Only the tunnels produced by the rigid Au₂₀ were used.

The simulations reveal that the immobilization of a nanoparticle is a very fast process, as the visual inspection of the trajectories allowed us to realize that this process take a time of the order of 0.03 ps. As soon as it reaches the tunnel, the high affinity of low-coordinated carbon atoms for gold triggers a diffusion cascade from the metal to the tunnel. As a result, an important volume of the nanoparticle gets inside the graphite hole, as shown in Fig. 4. This directly affects the relative height of the particles respect to the graphite surface. In Fig. 5, the height average for each particle is compared against the experimental results.⁷ As expected, the average height increases with the NP size. For smaller nanoparticles an excellent agreement between simulated and experimental values is reached. However, as the nanoparticle size increases, a slight deviation becomes evident. It is possible that this deviation is due to the underestimation of the affinity between gold and low coordinated carbon atoms. A greater affinity would draw more gold atoms from the NP, leading to the same outcome. Finally, it is important to stress that the EAM description of gold NPs improves with the particle size, underestimating the stability of smaller clusters. The larger the stability, the larger the resistance to capillarity. Therefore, it is possible that our simulations overestimate the height of all the immobilized NPs, but at small clusters this deviation is counteracted by a larger capillarity effect (due to the EAM underestimation of particle stability).

Fig. 4 Top and side views of four gold nanoparticles immobilized by a tunnel made by a Au₂₀ cluster. Note that with the exception of Au₁₄₇, all the nanoparticles conserve most of their icosahedral structure after immobilization.

Fig. 5 Comparison between simulations and experimental measurements of immobilized particle heights, relative to the position of the first graphite layer, as function of the size of the particle.

Another important outcome of the simulations is that most of the nanoparticles conserves their icosahedral structure. Indeed, only the smallest one (Au₁₄₇) adopts a spherical shape after immobilization. In the rest of the cases, most of the (111) faces are conserved; only the region close to the tunnel is compromised.

4 Conclusions

In the present work, molecular dynamics simulations were performed in order to understand from an atomistic perspective, the creation of nano-tunnels in graphite by Au₂₀ clusters and the subsequent immobilization of larger Au nanoparticles in such tunnels. A new C–Au potential was fitted to reproduce physical and chemical adsorption of Au nanoparticle on graphite. To do this, DFT calculations were performed and bond-order Morse-like potentials functions were fitted. This force field improves considerably the description of Au–graphite and Au–graphene dynamics compared with all previous pairwise potentials. This model is able to reproduce the phenomenon of capillarity by no other restriction than the differentiation of the Morse parameters according to the coordination of carbon atoms.

The height of immobilized NP from our simulations are in good agreement with that of the experiments. A small deviation its observed at larger cluster size, probably as a consequence of the subestimation of the tunnel size, or C–Au affinity. Among the immobilized nanoparticles tested only the smallest one (Au₁₄₇) lost its icosahedral structure. The rest of the particles presented restructuring only in the region close to the tunnel due to capillarity.

We believe this study can not only improve experimental designs to prepare new platforms for protein immobilization, but also assist in the design of new nanodevices for catalysis, nanoelectronics and nanophotonics.

Acknowledgements

This work was supported by FONCyT PICT-2015-2191, CONICET PIP 11220110100992, Secyt, Universidad Nacional de Córdoba. The authors thank the Texas Advanced Computer Center (TACC) for access to the Stampede supercomputer. Thanks also to the High Performance Computational Center (UNC, Córdoba, Argentina) for access to the Cristina supercomputer. J. A. d. l. R.-A., thanks Conacyt for the fellowship granted.

References

1. R. E. Palmer, S. Pratontep and H. G. Boyen, Nat. Mater., 2003, 2, 443–448.
2. M. Alexandre and P. Dubois, Mater. Sci. Eng., R, 2000, 28, 1–63.
3. J. Xu, J. Sun, Y. Wang, J. Sheng, F. Wang and M. Sun, Molecules, 2014, 19, 11465–11486.
4. E. Petryayeva and U. J. Krull, Langmuir, 2012, 28, 13943–13951.
5. R. Resch, S. Meltzer, T. Vallant, H. Hoffmann, B. E. Koel, A. Madhukar, A. Requicha and P. Will, Langmuir, 2001, 17, 5666–5670.
6. M. Ben Haddada, J. Blanchard, S. Casale, J.-M. Krafft, A. Vallée, C. Méthivier and S. Boujday, Gold Bull., 2013, 46, 335–341.
7. P. Rodríguez-Zamora, F. Yin and R. E. Palmer, J. Phys. Chem. A, 2014, 118, 8182.
8. J. Li, X. Li, H. J. Zhai and L. S. Wang, Science, 2003, 299, 864–867.
9. P. Gruene, D. M. Rayner, B. Redlich, A. F. G. van der Meer, J. T. Lyon, G. Meijer and A. Fielicke, Science, 2008, 321, 674–676.
10. L. J. Lewis, P. Jensen, N. Combe and J.-L. Barrat, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 16084–16090.
11. S. Kawai, A. Benassi, E. Gnecco, H. Söde, R. Pawlak, X. Feng, K. Müllen, D. Passerone, C. A. Pignedoli, P. Ruffieux, R. Fasel and E. Meyer, Science, 2016, 351, 957–961.
12. R. Guerra, U. Tartaglino, A. Vanossi and E. Tosatti, Nat. Mater., 2010, 9, 634–637.
13. M. S. Daw and M. I. Baskes, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 29, 6443.
14. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472–6486.
15. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
16. D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 41, 7892–7895.
17. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
18. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
19. J. H. Monkhorst and D. J. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192.
20. K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
21. V. R. Cooper, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 161104.
22. I. Lončarić and V. Despoja, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 075414.
23. I. Hamada and M. Otani, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 153412.
24. P. Wagner, V. V. Ivanovskaya, M. Melle-Franco, B. Humbert, J.-J. Adjizian, P. R. Briddon and C. P. Ewels, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 094106.
25. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
26. S. Krishnamurty, G. S. Shafai, D. G. Kanhere, B. S. de Bas and M. J. Ford, J. Phys. Chem. A, 2007, 111, 10769–10775.

---