Theoretical study of large-scale graphene on the Cu(111) surface using machine learning potential

Abstract

The properties of graphene supported on a metal surface are highly dependent on its size. However, for graphene with dimensions of several nanometers or larger, high-precision studies remain limited. This study investigates the interaction of large-scale graphene with the Cu(111) surface, focusing on the average adsorption and average formation energies of graphene configurations with varying carbon atom numbers. The force field parameters for the graphene–Cu(111) system were derived using the high-dimensional neural network potential, and its applicability and accuracy were validated. The study evaluates the average adsorption and average formation energies for graphene nanosheets, as well as zigzag and armchair graphene nanoribbon configurations and establishes relationships between these energies and the number of carbon atoms. The findings provide valuable insights into the structural evolution, stability, and adsorption behavior of graphene on Cu(111) from nanoscale to mesoscale.

---

1. Introduction

Carbon nanostructures, such as graphene and carbon nanotubes, have been widely utilized as reinforcement phases in metal matrix composites due to their excellent physical and chemical properties. The incorporation of graphene or carbon nanotubes into metal matrices can significantly enhance the mechanical strength,¹ thermal conductivity,² electrical conductivity,³ and tribological performance of the materials.⁴ Among these composites, copper-based composites exhibit excellent electrical and thermal conductivity, good processability, and superior corrosion resistance, making them highly suitable for applications in electronic packaging and power transmission devices.⁵

Currently, the main methods for fabricating graphene-reinforced copper matrix composites include electrostatic self-assembly,⁶ molecular-level mixing,⁷ chemical vapor deposition,^8–11 electrochemical deposition,¹² powder metallurgy,¹³ among others. Among these, powder metallurgy has gained the most widespread application due to its simplicity, high production efficiency, and cost-effectiveness.¹³ As a key technique within powder metallurgy, ball milling offers several advantages in the preparation of graphene-reinforced copper matrix composites, including uniform graphene dispersion in the copper matrix, ease of processing, suitability for large-scale production, and strong controllability.¹⁴ Therefore, ball milling holds great potential in this field. However, in practical applications, critical aspects such as the structural evolution of graphene, the stability of different graphene configurations, and its adsorption behavior on copper surfaces remain insufficiently studied. Additionally, during the ball milling process, factors such as graphene size variation, structural damage, and agglomeration can significantly impact the quality and performance of the composite material. Hence, an in-depth investigation into the underlying mechanisms of these factors is essential for optimizing the fabrication process.

Molecular simulations have been demonstrated as an effective research tool for comprehensively understanding these processes, providing insights into the stability of different graphene configurations on copper surfaces and analyzing their structural evolution, formation energy, and binding energy.¹⁵ Computational simulations can not only provide valuable insights into the adsorption behavior and energy variations of graphene with different sizes and configurations on copper surfaces but also serve as crucial theoretical guidance for fabricating graphene-reinforced copper matrix composites via ball milling. By exploring the size, morphology, and interactions of graphene with the copper matrix, we can optimize the fabrication process, laying a solid foundation for the controllable synthesis and efficient application of graphene-reinforced copper matrix composites.

Density functional theory (DFT) is widely regarded for its exceptional computational accuracy, particularly in comparison to classical molecular dynamics.

The calculation of formation energy and adsorption energy of materials on metal surfaces using DFT is a crucial approach for investigating surface chemistry and material properties.¹⁶ However, the computational demands of DFT are significantly higher, often surpassing those of molecular dynamics by several orders of magnitude, which presents considerable challenges in terms of computational cost. In traditional molecular dynamics simulations, force field parameters are typically derived from empirical formulas or quantum chemistry calculations. However, achieving an optimal balance between computational efficiency and accuracy in complex systems remains a significant challenge. To address this issue, recent advances in artificial intelligence have spurred the development of machine learning-based interatomic potential methods, particularly for complex material systems.

Notably, the high-dimensional neural network potential (HDNNP) leverages their ability to process nonlinear and high-dimensional data to train machine learning models using DFT results.¹⁷ These models accurately describe interatomic interactions in molecular dynamics simulations and precisely fit potential energy surfaces associated with intricate molecular interactions. By training machine learning models on high-precision quantum chemistry or DFT data, HDNNP-based force fields achieve accuracy comparable to that of quantum chemistry methods, while alleviating the computational burden and enabling large-scale molecular dynamics simulations. These advanced force fields have demonstrated exceptional performance in capturing chemical reaction dynamics,¹⁸ heterogeneous interactions,¹⁹ and the behavior of complex material systems,²⁰ making them indispensable tools for addressing the computational challenges inherent in traditional methods. Consequently, these approaches have become powerful and efficient means of advancing research in complex materials. Xu et al.²¹ trained a Cu–C–H neural network potential using a custom deep potential learning platform and applied it to identify the lowest-energy structures of carbon clusters on the Cu(111) surface, spanning cluster sizes from C₂ to C₂₄, while also calculating the corresponding formation energies for each structure. Previous studies by Xu et al.²¹ primarily focused on the configurations and formation energies of small carbon clusters on the Cu(111) surface. However, in large-scale experimental processes, such as the mixing of graphene and copper in a planetary ball mill, it becomes crucial to examine the stable configurations and adsorption energies of graphene structures, particularly those with a large number of carbon atoms (hundreds or even thousands) on the Cu(111) surface. Such analyses are essential for optimizing mixing processes and enhancing the quality of the complex material. Accordingly, this study aims to investigate the stable configurations, structural evolution, formation energies, and adsorption energies of graphene nanosheets and nanoribbons with higher carbon atom numbers on the Cu(111) surface.

The main research objectives of this study are as follows: first, the force field parameters for the graphene–Cu(111) system were fitted using the HDNNP, and its applicability and accuracy were validated. Second, the study of distance variations between different graphene configurations and the Cu(111) surface is conducted to provide quantitative and qualitative insights into the intermolecular interactions. Additionally, the average formation and the average adsorption energies of graphene nanosheets, zigzag and armchair graphene nanoribbon configurations with different carbon atom numbers on the Cu(111) surface were calculated, and relationships between these energies and the number of carbon atoms were derived. This study provides insights into the interaction of graphene with the Cu(111) surface, facilitating the design of graphene-based materials.

2. Computational details

2.1 Density functional calculations

All DFT calculations were performed using the CP2K software,²² using Gaussian and plane wave (GPW) formalism.²³ The Perdew–Burke–Ernzerhof (PBE) functional with the generalized gradient approximation (GGA) method was used for evaluating the exchange correlation functional,²⁴ combined with the molecular optimized DZVP local basis sets and Goedecker–Teter–Hutter (GTH) pseudopotentials.^25,26 Moreover, the DFT-D3 (BJ) dispersion correction was included to account for long-range van der Waals interactions. The plane wave density cutoff and energy convergence criterion were set at 400 Ry and 1 × 10⁻⁶ Ha, respectively.

The average formation energies E_f of the graphene/Cu(111) system were calculated by

E_f = (E(C_n/Cu(111)) − E(Cu(111)) − n × E(C))/n (1)

where E(C_n/Cu(111)), E(Cu(111)), and E(C) is the energy of graphene/Cu(111), Cu(111) surface, and single C atom derived from graphite, respectively. n is the atom number of C.

The average adsorption energies E_ads of graphene/Cu(111) were calculated by

E_ads = (E(C_n/Cu(111)) − E(Cu(111)) − E(C_n))/n (2)

where E(C_n/Cu(111)), E(Cu(111)), and E(C_n) is the energy of graphene/Cu(111), Cu(111) surface, and graphene, respectively. n is the atom number of C. C_n denotes a graphene nanosheet or nanoribbon consisting of n C atoms.

2.2 Training of the HDNNP

The machine learning potential was trained using the n2p2 v2.2.0 software.²⁷ The HDNNP consists of 2 hidden layers. Each hidden layer contains 10 nodes, allowing the model to capture complex, non-linear relationships in the input data. The architecture is designed to efficiently approximate the potential energy surface by learning the intricate interactions between atoms in the system. A cutoff radius of 14 bohr was applied to limit the range of atomic interactions considered by the symmetry functions, ensuring computational efficiency while retaining accuracy. A total of 280 symmetry functions were employed to encode the local atomic environments. The details about the type and parameters of symmetry functions are listed in the ESI.† The DFT training data were generated using the minimum hopping method.²⁸ This method efficiently explores different structures by initializing velocities along low-curvature directions, allowing rapid transitions over small energy barriers to neighboring basins.

The committee method reported by Schran et al.²⁹ was applied to enhance the accuracy of HDNNP. The committee disagreement helps track and control the model's accuracy compared to its ab initio method and training set. The computational cost of committee-HDNNP calculations is nearly the same as that of a single HDNNP, except for the initial requirement of training HDNNP multiple times. Eight HDNNPs were trained using the same training data, with the random number generator seed varied for each. Different graphene models are considered in the HDNNP. The dataset comprised 12 415 structures, including graphene/copper models with and without defects, and structure models of graphene/copper systems of the graphene nanosheet and nanoribbon (Fig. 1).

Fig. 1 Structural models of graphene/copper systems utilized in the HDNNP framework. (a) Graphene nanosheet, (b) defective graphene, (c) periodic graphene, and (d) graphene nanoribbon.

During the HDNNP fitting process, 90% of the data was allocated for training, while the remaining 10% was reserved for testing. In this study, a systematic evaluation of the HDNNP model was conducted by comparing its predicted energies and forces with reference DFT calculations to assess its accuracy and reliability. As shown in Fig. 2a, the energies of 12 415 data points obtained from HDNNP and DFT exhibit strong agreement. A similar trend is observed in the force calculations, as illustrated in Fig. 2b, further confirming the consistency between HDNNP and DFT results. The final HDNNP model achieved a root mean square error (RMSE) of 4.02 meV per atom for predicted energies and 155 meV Å⁻¹ for forces. By quantifying the consistency between HDNNP and DFT results and evaluating the RMSE, the findings demonstrate that the HDNNP model possesses the capability to accurately capture the potential energy surface and atomic interactions within the system.

Fig. 2 (a) Comparison of per-atom energy obtained from HDNNP and DFT. (b) Comparison of forces obtained from HDNNP and DFT.

There are many local minimum structures of graphene ribbons/flakes adsorption on the Cu surface. To find the most stable adsorption structures, we considered several different adsorption models for geometric optimization at first. Then molecular dynamics simulations were performed from 1000 K to 0 K over 10 ps with 1 fs time step in the canonical (NVT) ensemble for the lowest adsorption structures to further confirm the most stable adsorption structures using HDNNP. The structures and adsorption properties at different temperatures were analysed. For each structure, MD simulation was performed for 30 ps with 1 fs time step in the NVT ensemble at each temperature using HNNP. For the property analysis, we utilized structural configurations sampled from the final 20 ps of the simulations.

3. Results and discussion

3.1 Comparison of results between HDNNP and DFT

3.1.1 Periodic graphene. The adsorption of graphene on the Cu(111) surface is highly dependent on the specific positions of the carbon atoms relative to the copper lattice. The carbon atoms in graphene align directly above copper atoms in the second and third layers, occupying HCP and FCC sites. The carbon atoms in graphene are aligned directly above the copper atoms in the first and third layers, corresponding to the TOP and FCC sites. The carbon atoms in graphene are situated directly above the copper atoms in the first and second layers, corresponding to the TOP and HCP sites. Fig. 3 illustrates the optimized periodic graphene structures on the Cu(111) surface with different adsorption sites. Fig. 3c depicts the most stable configuration, while the relative energies of the other two configurations, compared to the most stable structure, are also provided.

Fig. 3 Optimized periodic graphene structures on the Cu(111) surface with different adsorption sites: (a) the carbon atoms in graphene occupying the HCP and FCC sites, (b) the carbon atoms in periodic graphene occupying the TOP and FCC sites, (c) the carbon atoms in graphene occupying the TOP and HCP sites. (d) Side view of the Cu(111) structure. Relative energies of the lowest-energy structures calculated using HDNNP and DFT, with DFT results shown after the slash. All energies are given in units of eV.

As shown in Fig. 3c, the most stable configuration occurs when the carbon atoms occupy the TOP and HCP sites, highlighting the significance of these positions in stabilizing the graphene structure. The relative energy obtained by DFT is 0.004 eV for the configuration in Fig. 3b, where the carbon atoms occupy the TOP and FCC sites. It indicates that the configuration in Fig. 3b is still stable, and this is consistent with the previous publication.³⁰ The configuration in Fig. 3a, where the carbon atoms occupy the HCP and FCC, have higher energies, and the relative energy obtained by DFT is 0.452 eV, suggesting less favorable adsorption configurations.

As shown in Fig. 3, the relative energy trends derived from the two methods exhibited a high degree of consistency. Additionally, the average formation energies and average adsorption energies calculated using the two methods were also evaluated for the periodic graphene/Cu(111) in this study. The average formation energy of periodic graphene/Cu(111) is 0.034 eV per atom (HDNNP) and 0.031 eV per atom (DFT), respectively. The average adsorption energy of periodic graphene/Cu(111) is −0.086 eV per atom (HDNNP) and −0.085 eV per atom (DFT), respectively. The similarity between the computational results obtained from the HDNNP method and those from the DFT method demonstrates the high accuracy of the HDNNP approach, confirming the reliability of its outcomes.

The cell lengths of the Cu unit cell are 3.599 Å (PBE + D3BJ) and 3.601 Å (HDNNP), and they are 2.469 Å (PBE + D3BJ) and 2.466 Å (HDNNP) for the graphene unit cell. The upper two layer by layer average distances in the z direction of the Cu(111) surface from top to bottom are 2.077 Å and 2.078 Å at the PBE + D3BJ computational level, and they are 2.085 Å and 2.074 Å using HDNNP. The average distances in the z direction from graphene to the Cu(111) surface for graphene occupying the TOP and HCP sites, graphene occupying the TOP and FCC sites, and graphene occupying the HCP and FCC sites are 3.110 Å, 3.150 Å and 3.270 Å at PBE + D3BJ computational level, and they are 3.108 Å, 3.110 Å and 3.179 Å using HDNNP. These results further demonstrate that the HDNNP reproduces the geometric structures very well. For the periodic graphene/Cu(111) model, we used the parameters of the Cu cell, which results in lattice mismatch 3.16% for graphene. We calculated the distortion degree of graphene ribbon on the Cu surface compared to the isolated graphene ribbon, which is 0.89% (zigzag graphene nanoribbon) and 3.14% (armchair graphene nanoribbon), respectively, which does not increase with increasing size.

3.1.2 Graphene nanosheet C₂₄. C₂₄ is a relatively small fragment of a graphene nanosheet, containing sufficient atoms to capture the essential characteristics of a hexagonal graphene nanosheet, while avoiding the significant computational cost associated with a full graphene sheet. The adsorption behavior and structural stability of C₂₄ on the Cu(111) surface at different adsorption sites were investigated. The preferred adsorption configurations and insights into the stability of C₂₄ on the metallic surface were obtained by examining the interactions between the carbon atoms and specific copper lattice positions. Furthermore, the energies of the optimized structures of the C₂₄ graphene nanosheet on the Cu(111) surface were calculated using both HDNNP and DFT methods to evaluate the accuracy of the HDNNP in this study.

Fig. 4 presents the optimized structures and corresponding relative energies of various configurations of the hexagonal graphene nanosheet composed of C₂₄. Fig. 4d shows the most stable configuration of C₂₄ on the Cu(111) surface, where the central six carbon atoms are located at the HCP and FCC sites. The results suggest that the HCP and FCC sites play a critical role in stabilizing the C₂₄ structure, particularly in the most stable configuration where the central six carbon atoms occupy these sites. This indicates that the registry between carbon atoms and the copper surface layers significantly affects the overall energy and stability of the system.

Fig. 4 Optimized graphene C₂₄ structures on the Cu(111) surface with different adsorption sites: (a) the carbon atom occupying the TOP and FCC sites, (b) the carbon atom occupying the TOP and HCP sites, (c) the carbon atom occupying the HCP and FCC sites, (d) the most stable configuration of C₂₄ on the Cu(111) surface. Relative energies of the lowest-energy structures calculated using HDNNP and DFT, with DFT results shown after the slash. All energies are given in units of eV.

In comparison to the structure shown in Fig. 4d, which corresponds to the most stable configuration, the structure seen in Fig. 4c exhibits slightly higher energy. The relative energy obtained by DFT is 0.472 eV. As depicted, during the evolution toward the most stable configuration, some carbon atoms shift from bridge sites (Fig. 4c) to top sites (Fig. 4d). The relative energy (shown in Fig. 4) trends obtained by the two methods were consistent. The formation energy of C₂₄/Cu(111) is 0.671 eV per atom (HDNNP) and 0.670 eV per atom (DFT), respectively. The average adsorption energy of C₂₄/Cu(111) is −0.584 eV per atom (HDNNP) and −0.558 eV per atom (DFT), respectively. The results indicated that the average formation and average adsorption energies calculated using the HDNNP were in close agreement with the DFT results. Meanwhile, according to the radial pair distribution functions of C–C and C–Cu shown in Fig. S1 in the ESI,† the HDNNP reproduces the geometric structures obtained in DFT calculations very well. It is worth mentioning that the C₂₄ graphene nanosheet was not included in the training dataset, which further demonstrates the reliability and accuracy of the HDNNP calculations.

3.2 Different configurations of graphene on Cu (111)

The structural configurations and energetic properties of three typical graphene structures—nanosheet, zigzag nanoribbon, and armchair nanoribbon—on the Cu(111) surface were systematically investigated. Simulations were performed on graphene samples with varying numbers of carbon atoms, revealing their trends in the average formation energy and the average adsorption energy as functions of the carbon atom number.

3.2.1 Structure properties of graphene nanosheets on the Cu(111) surface. As shown in Fig. 5, the graphene nanosheet samples were composed of structures with carbon atom counts: C₂₄, C₅₄, C₉₆, C₁₅₀, C₂₁₆, C₂₉₄, C₃₈₄, C₄₈₆, and C₆₀₀. Both top-view and side-view structural representations were obtained on the Cu(111) surface. Subsequently, the average formation energy (E_f) and average adsorption energy (E_ads) were calculated for each graphene structure.

Fig. 5 Graphene nanosheets with various carbon atom numbers on the Cu(111) surface, along with their optimized configurations, average formation energies, and average adsorption energies (a)–(i). Values after the slash represent average adsorption energies. All energies are given in units of eV.

To investigate the structural evolution of graphene nanosheets on the copper surface with the increasing carbon atom numbers, this study analyzes the variation in the distance between C atoms in nanosheets and the copper surface. The perpendicular distance d of each carbon atom of graphene to the Cu surface is computed, and the resulting distribution of these distances is denoted as f(d). The data are Gaussian broadened using a kernel density estimator with a bandwidth of 0.1. Each data point is convolved with a Gaussian kernel of this width. The resulting density estimate is automatically normalized, such that the total area under the curve equals 1. Fig. 6 illustrates the distribution of distances between the Cu(111) surface and graphene, and the distance between the periodic graphene and the copper surface was also labeled.

Fig. 6 Distance distribution of graphene nanosheets from the Cu(111) surface.

As the number of carbon atoms increases, the proportion of edge carbon atoms in the graphene nanosheet gradually decreases. As shown in Fig. 6, when the distances (d) between the carbon atoms in the graphene nanosheet and the Cu(111) surface are 1.65, 2.17, and 2.70 Å, the distribution of carbon atoms in the nanosheet shows a decreasing trend as the number of carbon atoms increases (Table S1, ESI†). The peaks at 1.65 Å and 2.1 Å reflect the strong chemical adsorption between the graphene edge and the metal surface, which is consistent with the phenomenon of chemical bonding between graphene and metal surfaces revealed by previous studies.^31,32 The peaks at distances greater than 2.70 Å indicate an equilibrium distance dominated by van der Waals interactions.

For carbon atom counts of C₄₈₆ and C₆₀₀, prominent peaks emerge at 3.70 Å, mainly attributed to the curvature effect at the edges of the graphene nanosheet. As the number of carbon atoms increases to C₆₀₀, a noticeable peak appears at 3.11 Å, which coincides with the distance between periodic graphene and the Cu(111) surface. This further suggests that as the number of carbon atoms in the graphene nanosheet increases, its structure gradually approaches that of periodic graphene.

Linear regression analyses were performed to evaluate the correlations between average formation and average adsorption energies and the number of carbon atoms. Fig. 7 shows the correlation between the average formation energy and the number of carbon atoms, as well as the correlation between the average adsorption energy and the number of carbon atoms for the graphene nanosheet configuration on the Cu(111) surface. The simulation results revealed clear trends in the energies. For graphene nanosheet, the average formation energy decreased with an increasing number of carbon atoms. This trend suggests that larger graphene structures are more energetically favorable to form on the Cu(111) surface. However, the absolute value of average adsorption energy decreased as the number of carbon atoms increased, and the absolute value of total adsorption energy increased.

Fig. 7 Graphene nanosheets on the Cu(111) surface: relationship between average formation energy and the number of carbon atoms (a) and (b), as well as the relationship between average adsorption energy and the number of carbon atoms (c) and (d). The fitting method excludes data points with a small number of carbon atoms. The dashed lines represent the results of periodic graphene.

When fitting the relationship between energy and the number of carbon atoms, two methods are employed: one including data points with a small number of carbon atoms and the other excluding them. Fig. 7b and d exclude data points with a small number of carbon atoms, while Fig. S1a and b (ESI†) include these data points. For graphene nanosheets, the average formation energy exhibits a linear relationship with the square root of the inverse carbon atom number n^−1/2, and the linear correlation is strong for both fitting methods. However, the average adsorption energy fitted using data points excluding small carbon atom numbers shows a linear relationship with n^−1/2 with a linear correlation R² value closer to 1. This indicates that, compared to including small carbon atom numbers in the fitting, this method provides a more accurate prediction of the average adsorption energy of larger nanosheets on the Cu(111) surface. Furthermore, the average adsorption energy of the graphene nanosheet at the limit of n^−1/2 is −0.080 eV per atom, which is very close to that (−0.086 eV per atom) of infinite graphene. To determinate the interaction strength between graphene edge atoms and the Cu(111) surface, we calculated the adsorption energy of C₆H₅ on the Cu(111) surface, where the interaction is dominated by the single C–Cu bond, and it is −2.49 eV. Due to the smaller graphene nanosheet having a higher edge atom ratio, its average adsorption energy is higher. The linear functions with respect to the inverse square root of n may be related to the dimension of the graphene nanosheet. The graphene nanosheets can be regarded as 1-dimensional materials. For the metal nanoparticles, which can be regarded as 0-dimensional materials, several theoretical studies demonstrated that the average cohesive energy is linear with respect to n^−1/3.^33,34 The increased dimension of the graphene nanosheet changes the n^−1/3 relationship into n^−1/2.

3.2.2 Structure properties of zigzag graphene nanoribbons on the Cu(111) surface. As shown in Fig. 8, the zigzag graphene nanoribbon samples included structures with the following carbon atom counts: C₃₀₀, C₄₀₀, C₅₀₀, C₆₀₀, C₇₀₀, C₈₀₀, C₉₀₀, and C₁₀₀₀. Both top-view and side-view structural representations were obtained on the Cu(111) surface, and the average formation energy (E_f) and the average adsorption energy (E_ads) were determined for each graphene structure. Additionally, from Fig. 8(i), it can be seen that the graphene edge is of the zigzag type.

Fig. 8 Zigzag graphene nanoribbons with various carbon atom numbers on the Cu(111) surface, along with their optimized configurations, average formation energies, and average adsorption energies (a)–(h). Values after the slash represent average adsorption energies. All energies are given in units of eV. (i) Side view of the configurations.

Fig. 9 shows the distribution of distances between zigzag graphene nanoribbons and the Cu(111) surface, with the distance between periodic graphene and the copper surface also indicated. In Fig. 9, when the distance (d) between the carbon atoms in the zigzag graphene nanoribbons and the Cu(111) surface is 1.64, 2.07, 2.70 and 2.92 Å, the distribution of carbon atoms exhibits a decreasing trend as the number of carbon atoms in the zigzag graphene nanoribbons increases. When d > 3.11 Å, as the number of carbon atoms in the zigzag graphene nanoribbons increases, the peak area gradually enlarges (indicating an increased distribution of carbon atoms), and the peak gradually shifts toward 3.11 Å. This further suggests that as the number of carbon atoms in the zigzag graphene nanoribbons increases, its structure progressively transitions towards that of periodic graphene.

Fig. 9 Distance distribution of zigzag graphene nanoribbons from the Cu(111) surface.

Fig. 10 illustrates the relationship between average formation energy and the number of carbon atoms, as well as the relationship between average adsorption energy and the number of carbon atoms for the zigzag graphene nanoribbon configuration on the Cu(111) surface. For the zigzag graphene nanoribbons, the average formation energy decreased as the number of carbon atoms increased, and the absolute value of the total adsorption energy increased with the number of carbon atoms. Fig. 10b and d include data points with a small number of carbon atoms, while Fig. S2a and b (ESI†) exclude these data points.

Fig. 10 Zigzag graphene nanoribbons on the Cu(111) surface: relationship between average formation energy and the number of carbon atoms (a) and (b), as well as the relationship between average adsorption energy and the number of carbon atoms (c) and (d). The fitting method includes data points with a small number of carbon atoms. The dashed lines represent the results of periodic graphene.

Fig. 10b, d and Fig. S2a and b (ESI†) show that average formation energy and average adsorption energy both exhibit a strong linear relationship with the inverse carbon atom number n⁻¹ (R² ≥ 0.9997) in both fitting methods. This indicates that both approaches reliably predict the average formation and the average adsorption energies of large zigzag graphene nanoribbons on the Cu(111) surface. The different relationship of the graphene nanosheet and graphene nanoribbon can be ascribed to their different dimension. The graphene nanoribbon can be regarded as two-dimensional materials, and the increased dimension changes the n^−1/2 relationship into n⁻¹.

3.2.3 Structure properties of armchair graphene nanoribbons on the Cu(111) surface. Fig. 11 shows the optimized configurations of graphene nanoribbons with various carbon atom numbers (C₂₆₄, C₃₆₀, C₄₅₆, C₅₅₂, C₆₄₈, C₇₄₄, C₈₄₀, and C₉₃₆) in the armchair arrangement on the Cu(111) surface, along with their average formation energies and average adsorption energies. The side views in Fig. 11(a)–(h) reveal significant curvature in the edge carbon atoms of the armchair graphene nanoribbon structure. Furthermore, Fig. 11i confirms that the graphene edge adopts an armchair configuration.

Fig. 11 Armchair graphene nanoribbons with various carbon atom numbers on the Cu(111) surface, along with their optimized configurations, average formation energies, and average adsorption energies (a)–(h). Values after the slash represent average adsorption energies. All energies are given in units of eV. (i) Side view of the configurations.

Fig. 12 presents the distribution of distances between the armchair graphene nanoribbons and the Cu(111) surface, with the distance of periodic graphene from the copper surface also indicated. When the distance (d) between carbon atoms in the armchair graphene nanoribbons and the Cu(111) surface is 1.72, 2.45 and 2.96 Å, the carbon atom distribution exhibits a decreasing trend. For d > 3.11 Å, as the number of carbon atoms in the armchair graphene nanoribbons increases, the peak area expands (indicating a higher carbon atom distribution), and the peak gradually shifts toward 3.11 Å, which is the distance between periodic graphene and the Cu(111) surface.

Fig. 12 Distance distribution of armchair graphene nanoribbons from the Cu(111) surface.

Fig. 13 presents the relationships between energy and the number of carbon atoms for the armchair graphene nanoribbons. As shown in Fig. 13a and c, the average formation energy decreases with increasing carbon atom number, and the absolute value of the total adsorption energy increases. Fig. 13b and d include data points with a small number of carbon atoms, whereas Fig. S3a and b (ESI†) exclude them. Both sets of figures demonstrate that the average formation energy and the average adsorption energy exhibit a strong linear correlation with the inverse carbon atom number n⁻¹ (R² ≥ 0.9998) across both fitting methods. This indicates that the relationships between energy and carbon atom number are well captured by the regression models. The linear fitting results for the three configurations are summarized in Table 1. These findings demonstrate that the derived expressions can effectively predict the formation energy and adsorption behavior of graphene with varying carbon atom numbers on the Cu(111) surface.

Fig. 13 Armchair graphene nanoribbons on the Cu(111) surface: relationship between average formation energy and the number of carbon atoms (a) and (b), as well as the relationship between the average adsorption energy and the number of carbon atoms (c) and (d). The fitting method includes data points with a small number of carbon atoms. The dashed lines represent the results of periodic graphene.

Table 1 Fitted linear expressions and coefficients of determination (R²) for energy-carbon atom number relationships

Structure	Fitting method 1^a	R ^2	Fitting method 2^a	R ^2
a Fitting method 1 excludes data points with a small number of carbon atoms, and fitting method 2 includes them. The units of Ef and Eads are eV per atom.
Nanosheet	E f = 3.1036n^−1/2 − 0.0049	0.99857	E f = 3.3741n^−1/2 − 0.0199	0.99971
	E ads = −3.3521n^−1/2 − 0.0800	0.99832	E ads = −2.2790n^−1/2 − 0.1434	0.97910
Zigzag nanoribbons	E f = 63.6444n^−1 − 0.0073	0.99995	E f = 63.6444n^−1 − 0.0073	0.99995
	E ads = −83.4020n^−1 − 0.0631	0.99970	E ads = −82.8869n^−1 − 0.0639	0.99994
Armchair nanoribbons	E f = 63.0084n^−1 − 0.0067	0.99987	E f = 63.008n^−1 − 0.0067	0.99987
	E ads = −54.5015n^−1 − 0.0780	0.99976	E ads = −54.0864n^−1 − 0.0783	0.99982

---

To emphasis the critical role in the interaction between graphene and the Cu(111) surface, we fitted the relationship between the edge atom ratio (n_edge/n) and average adsorption energy, and very good relationships are found for graphene nanosheets and nanoribbons as shown in Fig. S5–S7 (ESI†). It is noted that the graphene nanosheet has a good linear relationship between the square root of edge atom ratio (n_edge/n) and average adsorption energy, but the average adsorption energy is positive at the limit of (n_edge/n)^−1/2 to 0, which can be ascribed to the more significant deformation of the graphene nanosheet.

3.3 Temperature effect of the adsorption of graphene nanosheets and nanoribbons on the Cu surface

This section analyzes the structural evolution and adsorption properties of graphene nanosheets and nanoribbons on Cu surfaces at 300 K, 600 K, and 1000 K. All the graphene nanosheets and nanoribbons are stably adsorbed on the Cu surface, and no significant changes are observed even simulated at 1000 K. As shown in Fig. S8 (ESI†), the average adsorption energies of the graphene nanosheet and nanoribbon remain nearly constant from 0 K to 600 K, and very slight decreases are observed from 600 K to 1000 K. For the largest graphene nanosheet and nanoribbon models considered, the adsorption energy decrease is within 0.01 eV per carbon atom. In contrast, smaller graphene structures exhibit more pronounced reductions of 0.01–0.02 eV per carbon atom. The average formation energies of Cu-graphene slightly increase with increasing temperature (Fig. S9, ESI†). Similar to the average adsorption energies, the smaller graphene is more sensitive to the temperature. The insignificant changes of average adsorption and formation energies can be ascribed to the very strong interaction between Cu and graphene edge carbon atoms and low proportion of edge carbon atom in graphene, which leads to the graphene maintaining intimate contact with Cu. The majority of carbon atoms interact with the Cu surface through van der Waals forces, which exhibit minimal temperature dependence in the studied range (0–1000 K). According to the radial pair distribution functions of C–C and C–Cu (Fig. S10–S18, ESI†), no considerable shift of the first C–C peak position is observed, though the intensity is reduced. However, the first peak in the C–Cu radial distribution function shifts slightly to higher distances, suggesting a weakening of the C–Cu chemical interaction with increasing temperature. As shown in Fig. 14, all the distances from C to Cu surface in graphene nanosheet C₂₄ increase with increasing temperature. But for larger graphene nanosheets C₂₁₆ and C₆₀₀, the distances below 2.5 Å shift to higher values, and those beyond 2.5 Å shift to lower values. The decreased distances indicate that the protrusions near graphene edge become smoother as temperature increases. This graphene geometric reconstruction is also observed for zigzag and armchair graphene nanoribbons (Fig. S19 and S20, ESI†).

Fig. 14 Distance distribution of graphene nanosheets from the Cu(111) surface at different temperatures: (a) C₂₄, (b) C₂₁₆, and (c) C₆₀₀.

4. Conclusions

In this study, the HDNNP is utilized to study the adsorption of graphene with different configurations on the Cu(111) surface. The consistency between the energies obtained from NNP and DFT calculations, as well as the agreement of the force fields computed by both methods, demonstrates the accuracy of the NNP approach. The analysis of distance variations between different graphene configurations and the Cu(111) surface offers both quantitative and qualitative insights into intermolecular interactions at the graphene–metal interface. The analysis of distance variations is beneficial to understanding the structural properties between graphene and the Cu surface. It can also be used to explain the different graphene heights on the Cu surface at different graphene sizes observed from experimental measurements.

Moreover, the average formation energy and the average adsorption energy as functions of the number of carbon atoms in graphene are also investigated. The average formation energy decreases and the absolute value of the total adsorption energy increases with an increasing number of carbon atoms for all configurations. This implies that larger graphene structures are more stable thermodynamically and exhibit stronger adsorption on the Cu(111) surface. Zigzag and armchair graphene nanoribbon configurations demonstrate strong linear correlations between energies and n⁻¹, while the nanosheet configuration shows a linear relationship with n^−1/2. These relationships are supported by high R² values, indicating their robustness. The zigzag and armchair nanoribbon configurations showed higher adsorption capacities and better stability compared to the nanosheet structure, which exhibited smaller energy variations with increasing carbon atom numbers.

This study offers valuable insights into the interaction strength between graphene and the Cu(111) surface, as well as the stability of the graphene/Cu(111) system from nanoscale to mesoscale. According to the linear relationships found in this work, the graphene adsorption strength and stability can be quickly evaluated.

Author contributions

Jingli Han: resources, conceptualization, data curation, visualization, writing – original draft, writing – review and editing. Rubén Cabello: conceptualization and data curation. Jordi Bonet Ruiz: writing – review and editing, and formal analysis. Alexandra Elena Plesu Popescu: conceptualization, writing – review and editing. Sergi Dosta Parras: conceptualization, writing – review and editing. Camila Barreneche: formal analysis. Yongpeng Yang: resources, conceptualization, methodology, visualization, writing – review and editing, and funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the ESI.†

Acknowledgements

This work is supported by Hefei Advanced Computing Center and Zhengzhou National Supercomputing Center. Sergi Dosta Parras acknowledges funding from the European Union's Horizon Europe research and innovation program under grant agreement no. 101046835 (Thermodust project). Jingli Han and Yongpeng Yang acknowledge financial support from the China Scholarship Council (no. 202308410535 and 202308410412).

References

1. W. Chen, T. Yang, L. Dong, A. Elmasry, J. Song, N. Deng, A. Elmarakbi, T. Liu, H. B. Lv and Y. Q. Fu, Nanotechnol. Precis. Eng., 2020, 3, 189–210.
2. M. Jalali, N. Hassanzadeh, R. Alizadeh and T. G. Langdon, J. Mater. Sci., 2024, 59, 20900–20928.
3. M. Dadkhah, A. Saboori and P. Fino, Materials, 2019, 12, 2823.
4. M. Yang, Y. Liu, T. Fan and D. Zhang, Prog. Mater. Sci., 2020, 110, 100652.
5. P. Hidalgo-Manrique, X. Lei, R. Xu, M. Zhou, I. A. Kinloch and R. J. Young, J. Mater. Sci., 2019, 54, 12236–12289.
6. G. Shao, P. Liu, K. Zhang, W. Li, X. Chen and F. Ma, Mater. Sci. Eng., A, 2019, 739, 329–334.
7. J. Hwang, T. Yoon, S. H. Jin, J. Lee, T. Kim, S. H. Hong and S. Jeon, Adv. Mater., 2013, 25, 6724–6729.
8. X. Li, S. Yan, X. Chen, Q. Hong and N. Wang, J. Alloys Compd., 2020, 834, 155182.
9. X. Sun, X. Liu, Z. Sun, X. Zhang, Y. Wu, Y. Zhu, Y. Song, K. Jia, J. Zhang, L. Sun, W.-J. Yin and Z. Liu, Nano Res., 2024, 17, 4259–4269.
10. B. Jiang, D. Liang, Z. Sun, H. Ci, B. Liu, Y. Gao, J. Shan, X. Yang, M. H. Rümmeli, J. Wang, T. Wei, J. Sun and Z. Liu, Adv. Funct. Mater., 2022, 32, 2200428.
11. B. Liu, Z. Sun, K. Cui, Z. Xue, Z. Li, W. Wang, W. Gu, K. Zheng, R. Liu, Y. Zhao, M. H. Rümmeli, X. Gao, J. Sun and Z. Liu, Adv. Funct. Mater., 2023, 33, 2210771.
12. X. Zhao, J. Tang, F. Yu and N. Ye, J. Alloys Compd., 2018, 766, 266–273.
13. M. Anthony Xavior and H. G. Prashantha Kumar, Mater. Today Proc., 2017, 4, 3334–3341.
14. T. Han, J. Li, N. Zhao and C. He, Carbon, 2020, 159, 311–323.
15. W. Peng, K. Sun, R. Abdullah, M. Zhang, J. Chen and J. Shi, Appl. Surf. Sci., 2019, 487, 22–31.
16. H. Jia, L. Zhang, Y. Lu, J. Wu, X. Tan and Y. Cao, Mater. Today Commun., 2024, 38, 107624.
17. C. Liang and N. R. Aluru, Tuning interfacial water friction through moiré twist, ACS Nano, 2024, 18, 16141–16150.
18. Q. Chu, X. Chang, K. Ma, X. Fu and D. Chen, Phys. Chem. Chem. Phys., 2022, 24, 25885–25894.
19. P. Hu, M. Lei, Z. Sui, X. Zhou, D. Chen and Y. Zhu, ACS Catal., 2024, 14, 8602–86187.
20. R. Otsuka, K. Shimizu, H. Wakabayashi and S. Watanabe, Appl. Phys. Express, 2024, 17, 115501.
21. S. Xu, L. Wu, Y. Fan, Y. Liu, X. Zeng and Z. Li, J. Phys. Chem. C, 2024, 128, 5697–5707.
22. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
23. B. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–487.
24. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
25. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
26. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
27. A. Singraber, T. Morawietz, J. Behler and C. Dellago, J. Chem. Theory Comput., 2019, 15, 3075–3092.
28. S. Goedecker, J. Chem. Phys., 2004, 120, 9911–9917.
29. C. Schran, K. Brezina and O. Marsalek, J. Chem. Phys., 2020, 153, 104105.
30. M. A. Stoodley, L. A. Rochford, T. Lee, B. P. Klein, D. A. Duncan and R. J. Maurer, Phys. Rev. Lett., 2024, 132, 196201.
31. J. Wintterlin and M.-L. Bocquet, Surf. Sci., 2009, 603, 1841–1852.
32. W. O. Morais, J. P. C. Felix, G. R. D. Silva, C. M. Oliveira Bastos, A. C. Dias, E. M. Flores, C. R. C. Rêgo, V. S. R. de Sousa, D. Guedes-Sobrinho and M. J. Piotrowski, Sci. Rep., 2025, 15, 15496.
33. L. Vega, F. Viñes and K. M. Neyman, Nanomaterials, 2022, 12, 17.
34. F. Viñes and A. Görling, Chem. – Eur. J., 2020, 26, 11478–11491.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp02042e

---