A Molecular Dynamics Study of Carbon Dimerization on Cu(111) Surface with Optimized DFTB Parameters

YIN Di, QIU Zongyang, LI Pai, LI Zhenyu

Hefei National Laboratory for Physical Sciences at the Microscale, University of Science and Technology of China, Hefei 230026,P. R. China.

Abstract: Cu has been widely used as a substrate material for graphene growth.To understand the atomistic mechanism of growth, an efficient and accurate method for describing Cu-C interactions is necessary, which is the prerequisite of any possible large-scale molecular simulation studies. The semi-empirical density-functional tight-binding (DFTB) method has a solid basis from the density functional theory (DFT) and is believed to be a good tool for achieving a balance between efficiency and accuracy. However, existing DFTB parameters cannot provide a reasonable description of the Cu surface structure. At the same time,DFTB parameters for Cu-C interactions are not available. Therefore, it is highly desirable to develop a set of DFTB parameters that can describe the Cu-C system,especially for surface reactions. In this study, a parametrization for Cu-C systems within the self-consistent-charge DFTB(SCC-DFTB) framework is performed. One-center parameters, including on-site energy, Hubbard, and spin parameters,are obtained from DFT calculations on free atoms. Two-center parameters can be calculated based on atomic wavefunctions. The remaining repulsive potential is obtained as the best compromise to describe different kinds of systems. Test calculations on Cu surfaces and Cu- or C atom-adsorbed Cu surfaces indicate that the obtained parameters can generate reasonable geometric structures and energetics. Based on this parameter set, carbon dimerization on the Cu(111) surface has been investigated via molecular dynamics simulations. Since they are the feeding species for graphene growth, it is important to understand how carbon dimers are formed on the Cu surface. It is difficult to observe carbon dimerization in brute-force MD simulations even at high temperatures, because of the surface structure distortion.To study the dimerization mechanism, metadynamics simulations are performed. Our simulations suggest that carbon atoms will rotate around the bridging Cu atom after a bridging metal structure is formed, which eventually leads to the dimer formation. The free energy barrier for dimerization at 1300 K is about 0.9 eV. The results presented here provide useful insights for understanding graphene growth.

Key Words: Copper surface; Carbon dimer; DFTB; Molecular simulation

1 Introduction

Reaction dynamics of carbon species on Cu surfaces have attracted broad research interest recently due to its relevance to graphene growth. Epitaxial growth is promising for producing high quality graphene in large scale. To better control the growth process, it is desirable to understand the atomistic growth mechanisms. An intense research effort has been devoted to graphene growth mechanism studies in the past several years1,2. Recently, we revealed that carbon dimer is the dominant feeding species during graphene growth under low H2partial pressures3,4. In the minimum energy path of carbon dimerization, an important intermediate state is the bridging metal (BM) structure, where two carbon atoms are bridged by an elevated Cu atom5. BM structure can be formed spontaneously upon approaching of two carbon atoms. Then,by conquering an about 0.3 eV energy barrier, a carbon dimer can finally form.

Although the importance of carbon dimers in graphene growth and the dimerization process have already been studied from first principles, these studies are based on ideal surface models and do not take the substrate structure distortion effect into account, which can, however, be very significant due to the graphene growth temperature closing to the melting point of the substrate material Cu (about 1300 K)6. Therefore, instead of an exploration of the static potential energy surface based on the nudged elastic band method7, an explicit molecular dynamics(MD) simulation of the dimerization process on Cu surfaces at the experimental temperature is very desirable.

Considering the extremely low concentration of carbon species on the surface during graphene growth, an explicit MD simulation of the dimerization process is expected to be very time consuming. To speed up the simulation, a semi-empirical method such as density-functional tight-binding (DFTB)8,9can be used, which has a rigorous basis from density functional theory (DFT) and can thus provide a good balance between accuracy and efficiency. In DFTB, most parameters are directly calculated from electronic structure of atoms10–12, which ensures a reasonably good transferability. Recently, DFTB has been widely used to study the growth mechanisms of graphene and carbon nanotubes13–15.

A successful DFTB study is possible only if a good parameter set exists. Recently, Heine et al. has reported a parametrization of the electronic part of DFTB model for most elements in the periodic table16, but the repulsive part is just developed for some elements without copper17. In an older database of DFTB parameters for various chemical elements available online18, a complete DFTB parameter set for both copper and carbon is available in the matsci-0-3 section19,20.However, hetero-nuclear parameters to describe interactions between Cu and C are not available.

In this article, DFTB parameters for Cu and C are optimized aiming at a good description of surface reactions. Compared to the matsci-0-3 parameters, a significant improvement has been achieved in Cu parameters to describe surface systems. The obtained Cu-C parameters can predict correct adsorption geometries of Cu/C atoms on Cu surfaces. Adsorption energies also qualitatively agree with more accurate DFT calculations.These DFTB parameters are then used for MD simulations of the carbon dimerization on Cu surface. BM structures widely exist during the simulation and the process from such configurations to carbon dimers are also observed via a metadynamics simulation.

2 Parametrization details

Within the framework of self-consistent-charge (SCC) DFTB,total energy is expressed as a sum of three terms8

The first term requires a zero-order Hamiltonian matrix H0to be determined from the electronic part of DFTB parameters,which can be obtained via calculating electronic structure of neutral atoms8. As DFTB on-site parameters, diagonal elements of H0are the corresponding atomic orbital energies of a pseudo-atom confined in a wall potential.

where

The confinement potential is added to avoid the atomic orbitals becoming too diffusive21. Off-diagonal elements of H0are the corresponding two center integrals between two atomic orbitals on two different atoms10.

where Iμ∈, Jυ∈. All two center integrals can be expressed as linear combinations of Slater-Koster integrals where two atoms approach each other along a high symmetry axis22.

The second term in the right side of Eq. (1) is the energy contribution from charge transfer, where Δqaand Δqbare atomic charges and γabdepends on the DFTB Hubbard parameters8. Using the Janack’s theorem, the Hubbard parameters are calculated as the derivative of energy with respect to occupation number of the highest occupied shell23.The third term in total energy comes from remaining repulsive interactions between two atoms after the zero-order band structure energy and SCC energy are considered. It is assumed to be short-ranged with a cutoff distance Rc10.

Therefore, the whole DFTB parameters include on-site and Hubbard parameters, Slater-Koster integrals, and repulsive energy. The first three parts are also called the electronic part,which can be calculated from DFT within some approximations.All remaining parts which can’t be described by the electronic part are included in the repulsive energy and can be obtained by subtracting DFTB electronic energy from total energy of a reference system24, or via fitting multiple systems25,26.

Electronic part of the parameters was obtained using the parametrization module implemented in the Hotbit package21.For each element, the outmost valence shells (3d, 4p, 4s for copper and 2p, 2s for carbon) are considered. Atomic orbital energies of free atom calculated with the ATOM program implemented in the SIESTA package27are used as on-site parameters (Table 1) to ensure the correct limit for free atoms.For spin-polarized DFTB calculations, which is critical to get correct energies for individual atoms, additional spin parameters (Table 2) are required and they can be calculated from28:

DFTB calculations are performed using the DFTB + package29.DFT calculations are carried out with the Vienna Ab initio Simulation Package (VASP)30,31with PBE exchange-correlationfunctional32.

Table 1 On-site energies and Hubbard parameters for Cu and C.

Table 2 Spin constants for Cu and C.

2.1 Slater-Koster integrals

After parameters listed in Tables 1 and 2 are obtained from individual atom calculations, the obtained wave functions can be used as basis functions to calculate Slater-Koster integrals.The results can be affected by the confinement potential via parameter r0. Since it only depends on the electronic part of the DFTB parameters, band structure can be used to evaluate Slater-Koster integrals. For Cu, we calculate band structures for fcc, bcc, and sc phases. Then, average of energy deviation are calculated as16

The k points run over selected paths and the lowest 1-5 bands are considered. As shown in Fig. 1, the best r0for band structure occurs at about 0.222 nm. Similar to Cu, by calculating band structure for graphene, diamond and sc phases,the optimized r0for C is found to be about 0.148 nm.

Fig. 1 Average energy deviation in band structure with respect to r0 in the confinement potential.

Fig. 2 Band structures for fcc Cu (upper plate) and graphene(below plate).DFT results are in black, matsci-0-3 in red, our parameterization in blue.

Notice that the purpose of this parametrization is for molecular simulations, therefore it is not necessary for us to insist on the best band structure. Tests on adsorption systems suggest to use a smaller r0about 0.127 nm for carbon. With such r0adjustment, we can still obtain reasonable band structures as exemplified in fcc Cu and graphene using the final Slater-Koster integral parameters (Fig. 2). The difference between our results and that from matsci-0-3 DFTB parameters are small, and our results are slightly better than matsci-0-3 within the region near the Fermi level.

2.2 Repulsive energy

Repulsive energy is typically positive and monotonically decreases with atom distance33. A cutoff radius Rcis used,beyond which the energy difference between DFT and DFTB is assumed to be zero or a rigid shift. We set the cutoff radius of Cu―Cu repulsive potential to be 0.265 nm, larger than the energy minimum location for both DFT and DFTB band structure energy (Fig. 3).

Fig. 3 Energy with respect to lattice parameter of fcc Cu.Black curve is DFT total energy while blue curve is DFTB band structure energy.

After the cutoff radius is determined, the simplest way to obtain repulsive energy is subtracting the DFTB band energy from the total DFT energy of a dimer with different inter-atom distances, as shown in the red curve of Fig. 4. Beside dimer, we have also checked other high-symmetry systems which can be used to directly calculate the repulsive energy, including fcc Cu and Cu slab with an adsorbed Cu atom. Four Cu atom layers are used to construct a Cu(111) slab. A vacuum layer thicker than 1.5 nm was inserted between neighboring slabs in the c-axis direction to avoid their interaction.

Cu dimer gives the highest repulsive energy and fcc Cu gives the lowest repulsive energy, where even negative values are found in the region close to the cutoff radius. The two adsorption geometries produce similar repulsive curves between the dimer and bulk curves. A genetic algorithm based fitting of the repulsive energy usually leads to local oscillations in the curve, which will generate fluctuating forces if the obtained repulsive energy is used in the final DFTB parameters.Our test calculations suggest that simply using the repulsive energy generated by the surface system with a Cu atom adsorption at top site leads to reasonable DFTB parameters for our interested systems.

For C-C and Cu-C pairs, the cutoff radius is chosen to be 0.212 and 0.270 nm, respectively. They are smaller than the next nearest neighbor distances25and larger than the nearest neighbor distance. The C-C repulsive energy curve is calculated from graphene with different lattice parameters. The Cu-C curve is obtained from the slab models with one C atom adsorbed at the top site of Cu(111) surface.

2.3 Benchmark

Based on our DFTB parameters, the optimized lattice constant of fcc Cu is 0.375 nm, which is 3.74% and 3.02%larger than the experimental value (0.362 nm)34and the DFT result (0.364 nm). The matsci-0-3 parameters of Cu give a similar lattice parameter (0.371 nm) as what we obtained.Cohesive energy is obtained based on the optimized geometry,defined as

Fig. 4 Different Cu-Cu repulsive curves obtained from different configuration sets.Black curve for the slab model with one Cu atom adsorbed at top site, green curve for Cu atom adsorbed at fcc hollow site, blue curve for fcc bulks, and red curve for dimers.

where Eatomand Etotare the energy of an individual atom and bulk system with N atoms in the computational cell. Cohesive energy for bulk Cu obtained from our parameters is 4.53 eV. It is 3.49 and 4.18 eV from DFT5and matsci-0-3 parameters,respectively. Test calculations indicate that the overestimation of lattice parameter and cohesive energy can be improved by increasing the r0value in the confinement potential when generating electronic part of the parameters. However, such an adjustment leads to overall worse performance of the parameters. For example, if large r0value is used, a large hole will be formed in C atom adsorbed Cu surface, which is not observed in DFT calculations.

The pristine Cu(111) surface optimized with our DFTB parameters has a similar structure compared to DFT results,except for a systematical overestimation of the Cu―Cu bond length originating from the overestimation of Cu bulk lattice parameter. Notice that, when the matsci-0-3 parameters are used to optimize Cu(111) surface, the slab will break into two parts with an extraordinary large distance between the two top layers and the two bottom layers. Therefore, it cannot be used to study surface systems.

Adsorption structures and adsorption energies defined as

are tested with our DFTB parameters, where Esub, Eatom, and Etotare energies of optimized substrate, individual atoms, and the adsorption system containing N adatoms, respectively. As shown in Table 3, for different adsorption sites, the agreement between our DFTB parameters and DFT is reasonably well.

Since the Cu(111) surface is what we used to generate the repulsive energy of our parameters, we also check the performance of our parameters for the Cu(100) surface. Good consistency with DFT results for adsorption geometries isachieved as shown in Table 4, except for more overbinding compared with the Cu(111).

Table 3 Properties of Cu adatom on Cu(111) surface.

Table 4 Properties of Cu adatom on Cu(100) surface.

Before going to Cu-C interaction, we use graphene to test the performance of the parameters of the carbon part. The obtained lattice parameter is 0.247 nm, which is equal to that given by matsci-0-3 parameters and very close to the DFT value (0.248 nm). The cohesive energy is calculated to be 7.90 eV, which is very close to the DFT value of 7.91 eV35and smaller than 8.08 eV given by the matsci-0-3 parameters.

Carbon atom adsorption on Cu(111) is tested. As shown in Table 5, the agreement with DFT results of adsorption structures is reasonably well. However, we also notice that the Cu-C distance differences between top site and other sites are much larger than those in DFT results. The agreement of adsorption energy between DFTB and DFT is even better than the Cu atom adsorption case, where a systematic overbinding exists. As shown in Table 6, our DFTB parameters also give reasonable results for C atom adsorption on Cu(100) surface.

3 C dimerization on Cu(111)

Based on the DFTB parameterization described in the above section, we study carbon adsorption on Cu(111). C2occupies a fcc site and its neighboring hcp site. C3forms a chain and occupies three hollow sites. The stable structure for C6adsorption is a hexagon with each C atom adsorbed on a hollow site. In these structures, a larger simulation cell is adopted toavoid artificial interactions with neighboring cells. The relaxed geometries and adsorption energies (Table 7) are in good agreement with DFT results36.

Table 5 Properties of C adatom on Cu(111) surface.

Table 6 Properties of C adatom on Cu(100) surface.

Table 7 Adsorption energy per atom for C clusters on Cu(111) surface.

Among these surface carbon species, carbon dimer plays an dominant role in graphene growth3. The dimerization process has been theoretically studied based on ideal Cu(111) surface3,5.A BM structure is found to be spontaneously formed upon approaching of two carbon adatoms5. However, graphene is typically grown at a temperature as high as 1300 K6, where significant surface structure distortion, which is not contained in an ideal surface model, is expected. To take this effect into account, molecular dynamics simulations of carbon dimerization on Cu(111) surface is performed with a time step of 1 fs. In a 6 × 6 simulation cell with 144 Cu atoms, two carbon atoms initially 0.81 nm apart are put on the surface. As shown in Fig. 5, significant surface distortion is indeed observed after thermostating at 1000 K for 10 ps in the MD trajectory. Notice that the bottom layer is fixed at its bulk structure during MD simulations.

Test MD simulations at room temperature give a well ordered Cu surface and C atoms oscillate around surface fcc site. As the temperature increases to 500 K, C atoms can diffuse between fcc and hcp sites and approach each other. But dimerization is still difficult and no dimer forms in a 30-ps trajectory. This is caused by the insertion of a Cu atom between two C atoms. At even higher temperature (1000 K), significant surface distortion is observed, which leads to big holes sometimes and C atoms can thus diffuse into subsurface (Fig.5). C-Cu-C bridging structures similar to the BM structure reported previously5can easily form. Therefore, the dimerization process can’t be speeded up simply by increasing temperature and it is not practical to observe carbon dimerization using brute-force MD simulations.

Actually, we can estimate the equilibrium time for a Cu surface system with carbon monomers and dimers via a standard rejection-free kinetic Monte Carlo (KMC) simulation4.Here, we build a 200 × 200 grid to represent the Cu(111)surface, on which we deposit 10 to 100 C monomers. Four events, monomer diffusion, monomer combination, dimer diffusion, and dimmer decomposition, are taken into account.The corresponding energy barriers are set to 0.5, 0.9, 0.49 and 2.8 eV, respectively. KMC simulations give an equilibrium time at the order of 10?6s, which is not sensitive to the number of monomers we deposit.

Fig. 5 A snapshot of a 30-ps MD at 1000 K with two C atoms on Cu(111) surface.

In order to observe dimerization in our MD simulations at the experimental temperature (1300 K), the metadynamics technique37is applied to enhance sampling for rare events. The carbon-carbon distance is used as the collective variable. The height and width of the biased Gaussian potentials are set to 0.02 eV and 0.01 nm respectively, and they are deposited every 50 fs. To reduce simulation time, only interested regions in the phase space are sampled and half harmonic potentials are used to restrain two degrees of freedom. First, to avoid C atom diffuse into the Cu bulk, a half harmonic potential is enforced to carbon atoms and keeps C atoms stay on surface (constrain vanishes in the region 0.43 nm above the bottom Cu atom layer). At the same time, another wall potential is used to limit C―C bond being less than 0.40 nm which is larger enough to study the mechanisms of C―C bond formation and breaking.

In the metadynamics simulation, we start from a dimer on the surface. Even when the surface is melting, dimer is not decomposed until 874 Gaussian potentials are deposited. Then,one Cu atom is inserted into the two C atoms forming a linear C-Cu-C BM structure (Fig. 6a). Notice that, due to distorted surface structure at high temperature, the BM structure observed here is different from the BM structure reported previously5with copper atom slightly pull up and one carbon atom fall into the subsurface. Actually there is no distinct Cu layer except the fixed bottom layer. The two C atoms can rotate around the middle Cu atom (blue atom in Fig. 6) and the middle Cu atom prevents them from bonding again. With the simulation continued, we can see a fluctuation of the C―Cu―C angle which is coupled with the C-C distance change between 0.30 to 0.40 nm. In the next 40-ps trajectory after the decomposition of the carbon dimer, the C―Cu―C angle gradually decreases more and more and finally the two carbon atoms form a dimer again (Fig. 6b–d). Notice that such a process is similar to the minimum energy path of carbon dimerization we found previously5.

Fig. 6 Snapshots of the metadynamics simulation trajectory in which a dimer forms.C atoms are in silver and Cu atoms are in blue and yellow.

Fig. 7 Potential of mean force for carbon dimerization on copper substrate.

With more and more Gaussian potentials deposited, the dimer repeats breaking and bonding. After 4500 Gaussian potentials deposited, the C-C distance can easily change within the region from 0.12 to 0.40 nm, which means that the free energy is convergent. From the obtained free energy surface(Fig. 7), we can know that the free energy barrier of carbon dimerization is about 0.9 eV, which is significantly larger than the value (0.3 eV) from ideal surface model mainly due to the surface distortion effect. Notice that, in the breaking region(C-C distance from 0.25 to 0.40 nm), the free energy profile is less smooth. This is because most sampling configurations in this region are BM structures and the melting Cu atoms will more strongly affect the motion of C atoms.

4 Conclusions

New SCC-DFTB parameters for Cu and C are developed for surface reactions. Benchmark results indicate that this parameter set is reliable. MD simulations suggest that carbon dimerization on Cu(111) surface is not as easy as we may expect at high temperatures, which is a result of the strong surface structure distortion. Metadynamics simulations suggest that the BM structure plays an important role in carbon dimerization. A 0.9 eV free energy barrier is predicted at 1300 K.

Acknowledgment: We thank USTC-SCC, SCCAS, Tianjin,Guangzhou and Shanghai Supercomputer Centers for CPU time.