Immobilization of Au nanoparticles on graphite tunnels through nanocapillarity

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.


Introduction
Immobilization of nanostructures on selective surfaces is an important technical issue for several reasons 1 .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. 2 Iron magnetic nanoparticles have been used to immobilize proteins acting as active agents in enzymecatalyzed reactions, 3 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. 4In 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 2 surfaces by embedding them in octadecylsiloxane 5 or by covering them with amines and thiols. 6nother 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.
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 7 used clusters of 20 atoms in size (Au 20 ).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,9nce 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/atom).The deposition at this softlanding regime, preserve not only the size but also the general shape of the nanoparticles. 1 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. 10 , 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" 10 .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

Interaction Models
The interaction between Au atoms is modelled using the Embedded Atom Model (EAM), due to Daw and Baskes 13 .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 14 , 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 15 .Vanderbilt ultrasoft pseudopotentials 16 were used with the Perdew-Burke-Ernzerhof 17 (PBE) and the local density approximation 18 (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 19 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 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

C-Au potential parametrization
In order to successfully simulate the experiments reported by R. Palmer 7 , 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 6−16 .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).
Chemisorption processes were simulated by the adsorption of cluster Au 7−19 at the edges of a graphene flake, shown in   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.

MD simulations
We performed two kinds of molecular dynamics simulations using the LAMMPS code 25 , 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 20 clusters in graphite and studied the damage provoked at the surface.The implantation of Au 20 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 20 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.
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 20 cluster.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 20 cluster loses its pyramidal shape, although according to ab-initio MD it should melt at a temperature of around 800 K. 26 It is not surprising that EAM can not describe the stability, angularity, and relativistic effects of Au 20 , 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 20 pyramid, and compared the corresponding surface damage with that of the soft cluster.The MD runs for these rigid clusters were 20500 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 20 .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 20 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 20 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.MD simulations also revealed other observables.At kinetic energies lower than 0.5 eV, in most cases soft Au 20 does not penetrate into the graphite surface, it bounces back instead.At kinetic energies higher than 0.5 eV the soft Au 20 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.

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 7 the same gold nanoparticles were considered for the simulations, namely Au 147 , Au 309 , Au 561 , and Au 923 .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 20 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. 7As expected, the average height increases with the NP size.For smaller nanopar- ticles 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).Another important outcome of the simulations is that most of the nanoparticles conserves their icosahedral structure.Indeed, only the smallest one (Au 147 ) adopts a spherical shape after immobilization.In the rest of the cases, most of the (111)

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 20 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 147 ) 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.

Fig. 1
Fig. 1 (a) Bottom view of two of the five clusters used to study physisorption (Au 3−9 and Au 7−19 ).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 7−19 ) 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.

Fig. 1c .
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 24 (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

Fig. 2
Fig. 2 Relative orientations of Au 20 with respect to the graphite surface (below the particles, not shown).a) Face; b) Edge; c) Vertex.

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

Fig. 5
Fig.5Comparison 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.

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

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 (α).