Assessment of Classical Force-Fields for Graphene Mechanics

Abstract

The unique properties of graphene have attracted the interest of researchers from various fields, and the discovery of graphene has sparked a revolution in materials science, specifically in the field of two-dimensional materials. However, graphene synthesis’s costly and complex process significantly impairs researchers’ endeavors to explore its properties and structure experimentally. Molecular dynamics simulation is a well-established and useful tool for investigating graphene’s atomic structure and dynamic behavior at the nanoscale without requiring expensive and complex experiments. The accuracy of the molecular dynamics simulation depends on the potential functions. This work assesses the performance of various potential functions available for graphene in mechanical properties prediction. The following two cases are considered: pristine graphene and pre-cracked graphene. The most popular fifteen potentials have been assessed. Our results suggest that diverse potentials are suitable for various applications. REBO and Tersoff potentials are the best for simulating monolayer pristine graphene, and the MEAM and the AIREBO-m potentials are recommended for those with crack defects because of their respective utilization of the electron density and inclusion of the long-range interaction. We recommend the AIREBO-m potential for a general case of classical molecular dynamics study. This work might help to guide the selection of potentials for graphene simulations and the development of further advanced interatomic potentials.

1. Introduction

Graphene, a two-dimensional material, has numerous remarkable properties, such as high strength, electrical conductivity, and thermal conductivity [1,2,3,4]. Graphene has been extensively studied by researchers from various disciplines, resulting in numerous discoveries over the past two decades [5,6,7,8,9,10]. Due to its unique properties, graphene is utilized in broad applications in various fields, including energy, biomedicine, environment, electronics, material engineering, etc. [11]. The discovery of graphene was a scientific breakthrough in nanomaterials, which encouraged the discovery and engineering of other 2D materials [12,13,14,15]. There are two primary approaches to synthesizing graphene: the bottom-up construction, in which the material is built atom-by-atom to form the graphene sheet, and the top-down separation, where the graphene sheet is obtained by separating the graphite [16,17,18,19].

The difficulty of producing graphene in terms of the cost, size, and technique might cause researchers to have difficulty exploring and studying the material’s properties. Molecular dynamics (MD) simulation is a numerical tool that enables researchers to study the behavior of the material and the interaction between the particles at the microscopic level [20,21,22,23]. MD simulation is the bridge between theoretical and experimental research, enabling researchers to study materials that have not yet been synthesized, or to study them under conditions that cannot be achieved experimentally due to their complex implementation, the high risk that could be involved in such experiments, or the high cost. These advantages are typical of numerical computation and simulation in general and not only molecular dynamics.

MD is a well-established method for simulating various processes, properties, and applications of materials [24,25,26]. The accuracy of the MD simulation depends on multiple factors, such as the size of the system; especially when there are defects in the structure, it is required to ensure that the system is large enough. Another factor is the empirical potential used to represent the force between the atoms [27,28,29]. Different potential types are described with various formulae and functions.

For example, the Lennard–Jones potential is a commonly used pairwise interaction potential in MD simulations. It describes the interaction between two neutral atoms or molecules as a function of the distance between them. The Lennard–Jones potential [30,31,32] can be used to describe the interactions between carbon atoms in graphene, as well as interactions between graphene and other species (e.g., adsorbates and solvents). As a non-bonded potential, this pairwise potential describes the van der Waals interactions between non-bonded atoms. It accounts for both the attraction at long distances and the repulsion at short distances. The Embedded Atom Method (EAM) potential is a many-body potential that takes into account the electronic structure of the atoms. Although primarily used for metals, the EAM potential [33] or its variants can be adapted to describe the interactions within graphene, especially when considering the metallic nature of graphene’s electrons. The Tersoff potential is a bond-order potential specifically designed for simulating covalent materials such as silicon and carbon. It accounts for both the bond length and bond angle dependencies. The Tersoff potential [34,35,36] is widely used to describe the interactions between carbon atoms in graphene, as it accurately captures the covalent bonding characteristics of this material. As a three-body potential that includes bond-order effects, it is important for accurately describing the carbon–carbon interactions in graphene. It accounts for the variation in the bond strength depending on the local environment. The REBO potential, short for the Reactive Empirical Bond-Order interatomic potential of the second generation, is an advanced interatomic potential used in atomistic simulations to describe the motion and interactions of individual atoms [37]. The REBO potential uses the bond-order parameter to calculate the bond energy. As the bond-order changes due to bond stretching, bending, or other deformations, the bond energy is adjusted accordingly. The AIREBO-m potential modifies the repulsive part of the original AIREBO (Atomic Interaction Potential with Bond-Order) potential [38]. By replacing the LJ potential with a Morse potential and optimizing the functional form, the AIREBO-M potential improves the accuracy of simulations involving intermolecular interactions. Neural Network Potentials [39,40] are potentials derived from machine learning algorithms trained on DFT data. They can be very accurate and computationally efficient for large-scale MD simulations of graphene. Machine learning potentials have been widely introduced in recent years [41,42,43,44,45,46]. We focus on classical force-fields in this study, since they are a reasonably good compromise of accuracy and computation cost.

When choosing a potential for graphene simulations, researchers must balance the trade-off between computational efficiency and the accuracy required for the specific application. For example, while harmonic potentials are computationally efficient, they may not accurately capture the behavior of graphene under extreme conditions. In contrast, DFT-based methods provide high accuracy but are computationally more demanding and thus limited in the time and length scales that they can access. The choice of potential will depend on the specific properties of graphene that are being investigated and the computational resources available. Thus, the potential should correctly resemble the experimental behavior of the forces among the atoms, which is what we are trying to accomplish in this paper.

In this paper, tensile tests of pristine graphene are simulated, followed by comparing the results produced by MD simulation using various empirical potentials with the experimental results. Additionally, we evaluated the agreement of the potentials with experiments for pre-cracked graphene by examining its fracture toughness.

2. Model and Methods

2.1. Illustration of the Structure

In this work, all MD simulations are performed using LAMMPS [47]. The graphene unit cell is constructed as shown in Figure 1, where the positions of the atoms are determined relative to the lattice vectors (note that there are multiple choices for the unit cell and the lattice vector). The unit cell is replicated in the x- and y-directions.

The final structure, shown in Figure 1c, has dimensions of 5.11 nm × 5.16 nm × 0.334 nm, and the total number of atoms is 1008. The model’s boundary conditions are periodicity in the x- and y-directions and a fixed boundary along the z-direction. This represents the intermolecular force of the bonds between the atoms. With the initial model, we need to equilibrate the system to ensure that there is no overlap among the atoms and no deformation in the structure. The equilibration is performed under the temperature of 300 K and pressure of 0 GPa. In the equilibration process, the system’s energy is minimized: we need to find the configuration that results in the minimum energy, since the less energy the configuration has, the more stable it is.

The uniaxial tensile test was conducted to study the mechanical properties. Loading was applied in the horizontal direction (the zigzag direction) with an engineering strain rate of 0.001 per picoseconds. The load ceased when the material failed. Then, the stress–strain relationship was obtained.

2.2. Description of the Used Potentials [47]

The first potential to be discussed is the Tersoff potential (after the researcher who formulated it, Jerry Tersoff); then, the REBO (Reactive Empirical Bond-Order) and the AIREBO (Adaptive Intermolecular REBO) potentials are discussed. The Tersoff potential only considers the nearest-neighbor interaction of two atoms, and the system’s total energy is equal to the sum of all these interactions. It is given by the following expression [48,49]:

$$E={\displaystyle \frac{1}{2}}{\sum }_i^N{\sum }_{j\ne i}^N{f}_C\left(r_{ij}\right)[V_R\left(r_{ij}\right)+b_{ij}{V}_A\left(r_{ij}\right)],$$

(1)

where N is the number of atoms; f_c is the smooth cutoff function to switch off the effect of the potential after a certain distance; r_ij is the distance between atom i and atom j; V_R is the repulsive part of the potential; V_A is the attractive part of the potential; b_ij is the bond-order term, which depends on the coordination and the angle between the atoms.

The issue with the Tersoff potential is that it does not distinguish between the different types of bonds and only considers the bond’s length. As a result, it overestimates the system’s energy (this might be why most of the Tersoff potentials in

Table 1. Comparison between the experimental results and simulation results for different interatomic potentials in terms of the mechanical properties.

Potential Name	Young’s Modulus (GPa)	Failure Strength (GPa)	Fracture Strain	Toughness (GPa)
Experiment [50,51,52]	1000 ± 100	130 ± 10	0.25	21.1
CH.airebo-m	1025.11	93.93	0.169	10.5
CH.rebo	1043.13	131.53	0.282	18.08
C.meam	1022.26	83.81	0.400	7.42
CCu	798.76	109.33	0.400	14.45
C.lcbop	950.93	175.37	0.400	39.83
BNC	1009.93	125.65	0.235	18.36
SiCGe	1164.12	148.07	0.270	18.03
SiC-a	1151.20	150.19	0.275	18.64
SiC-b	1046.34	142.19	0.284	18.48
SiC-c	1164.12	148.07	0.270	18.03
SiC-d	583.48	106.43	0.400	12.77
SiC-e	300.08	42.24	0.400	5.59
SiC-gw	310.06	41.62	0.400	5.59
FeC	828.16	131.78	0.400	22.13
SiC-f	574.29	98.10	0.400	10.99

have a higher failure strength than the experiment potentials). The REBO potential is introduced to take into account the conjugated and the radical character of the system or the configuration (the conjugated configuration is defined as the ability of the system to resonate, which can be due to the existence of double or triple bonds, or the radical characteristic; a radical configuration contains unpaired electrons) [49].

The radical and the conjugation effects are included in the REBO potential as the correction to the bond-order term; consequently, instead of only having b_ij, it will be $b_{ij}+{\sum }_{ij}^{RC}+b_{ij}^{DH}$, where the second term describes the radical and conjugated characteristic and the third term describes the dihedral angle for the double bonds between carbon atoms. However, the REBO potential is restricted only to a short-range interaction, which is why the AIREBO potential was introduced, which describes the long-range interaction as the Lennard–Jones potential and the torsional potential due to the rotation of the bonds. Note that we have to use the potential CH.airebo-m, where m stands for the Morse potential [53]. The Morse potential is expected to produce more accurate results, since it converges more quickly as you leave the atom and diverges more slowly than the LJ potential as you go toward the atom [49].

$$E^{AIREBO}={\displaystyle \frac{1}{2}}{\sum }_i^N\sum _{j\ne i}^N[E_{ij}^{REBO}+E_{ij}^{LJ}+{\sum }_{k\ne i,j}{\sum }_{l\ne i,j,k}{E}_{ijkl}^{TORSION}]$$

(2)

The next potential is the BOP (bond-order potential) [54], which can be described as an improved version of the Tersoff potential that uses the theory of quantum mechanics to distinguish between the single and double bonds in the structure or the configuration. An advantage of the BOP potential is that it can be used for different phases, since it is formulated using quantum mechanics. However, the difficulty lies in its parameterization (that might be why it yields inaccurate values for the pristine graphene, since it was parameterized for the CCu). The energy of the bond in the BOP potential considers the bond-order and the bond integral, where the first is defined as follows [55]:

Bond order = (number of bonding electrons − number of antibonding electrons)/2

When two atoms move close to one another, their orbital overlap and the electrons within the overlap region contribute to the formation of bonds between the atoms. These electrons are called bonding electrons; the electrons outside of the overlap region are known as antibonding electrons. Additionally, the bond integral is defined in terms of the probability for the electron to jump from one orbital in the atom to another atom’s orbital. The expression below describes the energy of the system that uses the BOP potential [55]:

$$E={\displaystyle \frac{1}{2}}{\sum }_{i=1}^N\sum _{j=i_1}^{i_N}{\Phi }_{ij}\left(r_{ij}\right)-{\sum }_{i=1}^N\sum _{j=i_1}^{i_N}{\beta }_{\sigma ,ij}\left(r_{ij}\right)·{\Theta }_{\sigma ,ij}-{\sum }_{i=1}^N\sum _{j=i_1}^{i_N}{\beta }_{\pi ,ij}\left(r_{ij}\right)·{\Theta }_{\pi ,ij}$$

(3)

where ϕ_ij(r_ij) is the repulsive potential in the short-range interaction between two atoms; β_σ, and β_π, are the bond integrals for single and double bonds, respectively, which are a function of the coordination of the atoms in the structure, and which include the cutoff distance function as well; Θ_σ, and Θ_π, are the bond-orders for single bonds and double bonds, respectively.

Note that all these formulae for the previous potential are only general expressions, and their functions, which depend on the coordination or the local environment, are more complicated. Nonetheless, these general formulae help us to imagine the behavior of the potential.

Like the Tersoff potential, the BOP potential only considers the short-range interaction between the nearest neighbor. It ignores the long-range interaction, yielding an inaccurate result, especially in a multilayer system, where it ignores the interaction between the interplanar atoms. Instead, this is included in the LCBOP (long-range Carbon Bond-Order Potential). It is noted that this potential focuses on the interaction between carbon atoms, since the BOP potentials are difficult to parameterize, as mentioned earlier. The general formula for the LCBOP is given below [56] as follows:

$$E={\displaystyle \frac{1}{2}}{\sum }_{i,j}^N(f_{c,ij}{V}_{ij}^{SR}+S_{ij}{V}_{ij}^{LR})$$

(4)

where $V_{ij}^{SR}$ and $V_{ij}^{LR}$ are the short-range potential and long-range potential, respectively; f_c,ij is the cutoff distance function for the short-range interaction; S_ij is the switching function that is used to exclude the nearest neighbors of the long-range interaction, and it is given as S_ij = 1 − f_c,ij.

The MEAM (Modified Embedded Atom Method) [57] is another important potential, which is derived from the EAM potential [58]. The structure of the material depends on the distribution of nuclei and electrons. Therefore, the wave function of electrons can be used to describe the behavior of the atoms in the material and the energy associated with the structure of the material. However, the problem is that even a tiny sample contains many electrons. Finding the wave function for all these electrons would be nearly impossible, since the wave function of the electrons depends on the coordination of all the nuclei and the electrons in the material. Instead of using the wave function to describe a microscopic system, the EAM potential uses the electron density to describe the status of the electrons in a given system. The energy of one atom in a structure that uses the EAM potential is the following [59,60]:

$$E=F_i\left({\sum }_{j\ne i}{\rho }_i\right(r_{ij}\left)\right)+{\displaystyle \frac{1}{2}}{\sum }_{j\ne i}{\Phi }_{ij}\left(r_{ij}\right),$$

(5)

where F_i is the embedding energy function, which is a function of the ρ_i electron density; ϕ_ij is the pairwise interaction within the cutoff distance range, and it is the function of the coordination or distance between the atoms r_ij. The MEAM potential is an extension of the EAM potential that considers the existence of double and triple bonds.

3. Results and Discussion

3.1. Pristine Graphene

The analysis will focus on the following four parameters that can be extracted from the stress–strain curve: Young’s modulus [61,62], the material strength or the critical stress, the fracture strain, and the material’s toughness. Young’s modulus is obtained by calculating the slope of the linear portion of the stress–strain curve within the elastic region. It represents the material’s ability to endure or withstand the deformation or the change in length when uniaxial stress is applied to it. Young’s modulus E [63,64,65,66] can be expressed as E = FL₀/(A|L − L₀|), where F is the applied force on the surface of the material, A is the area, L₀ is the initial length, and L is the length after loading.

The material strength is the maximum stress applied to the material before it becomes fractured. The fracture strain is the strain that corresponds to the maximum stress or the strain at which the material is fractured. The toughness is the energy per unit volume that the material absorbs before it becomes cracked and is measured in (Pa). Still, usually, we are interested in the fracture toughness, since it is more realistic to have cracked graphene than pristine graphene, where the material is pre-cracked when the stress is applied, as can be seen in the discussion for the pre-cracked graphene. The fracture toughness is measured in (Pa$\sqrt{m}$). After obtaining these parameters out of the stress–strain curve, we will compare them with the experimental ones to evaluate the agreement of the results produced for each potential [24,64,67].

Most of the curves in Figure 2 reflect the brittle nature of graphene. Tersoff potentials, including BNC, SiCGe, SiC-a, SiC-b, SiC-c, and SiC-d, exhibit similar stress–strain curves as that of CH.airebo, CH.rebo, C.meam, CCu, and C.lcbop. The potentials of SiC-gw and SiC-e do not agree with the experimental results, since they show a stress–strain curve for plastic materials with much lower stress. Table 1 shows the experimental and obtained values using the specified potential [50,51,52].

As shown in Figure 3, most of the Tersoff potentials exhibit an excellent agreement in terms of Young’s modulus, except for SiC-d, FeC, and SiC-f, where the first and the last have a Young’s modulus which is half of the experimental one, and the second has a value of 828.16 GPa, which is less than anticipated and in the range of 900–1100 GPa. However, it can be noticed that the other Tersoff potentials overestimate the failure strength by a small amount, except for BNC and SiC-b. Therefore, the best two potentials among the considered ones (in the case of pristine graphene) are the BNC potential, since the percentage difference in the fracture strain is 6%, and SiC-b potential, which has a percentage difference of 13.6% in terms of the fracture strain.

The CCu potential shows low Young’s modulus and failure strength values, and a high fracture strain. Accordingly, it is not suited for the pristine graphene simulation. The C.lcbop potential has a Young’s modulus that falls within the experimental range, while it overestimates the failure strength and the fracture strain. The C.meam and CH.airebo-m potentials have an excellent agreement with the experiments regarding Young’s modulus. However, they underestimated the failure strength; the first overestimated the fracture strain and the second underestimated it with a 32.4% percentage difference. Finally, we have the CH.rebo potential, with a Young’s modulus of 1043.13 GPa, which falls within the experimental range of 900–1100 GPa, a close failure strength value to the experimental one, with a percentage difference of 0.12%, and a fracture strain of 0.282, with percentage difference of 1.3%.

Therefore, the ideal potential for simulating pristine graphene is the CH.rebo and the BNC-Tersoff.

3.2. Pre-Cracked Graphene

In this section, the agreement of the potentials with the experiments for the pre-cracked graphene has been assessed using the stress intensity factor (SIF), the fracture toughness (KIC), crack convergence, and Griffith’s criterion as our criteria. The computational details have remained unchanged, except for the size of the system, which increased by 16 times that of the pristine graphene system. Additionally, we have included a gap at the center of the sheet to simulate the pre-cracked configuration, as illustrated in Figure 4a.

The fracture toughness is defined as the product of the critical stress and the square root of the crack’s size. Critical stress is defined as the strength of failure for a brittle material like graphene [68], as follows:

$$\mathit{fracture}\mathit{toughness}={\sigma }_c\sqrt{a_0}$$

(6)

where σ_c is the critical stress in GPa and a₀ is the half initial length of the crack. The stress intensity factor K_c (SIF) is defined as follows:

$${K_c=\sigma }_c\sqrt{{\pi a}_0}$$

(7)

The terms (SIF) and fracture toughness can be used interchangeably, since they differ by a factor of $\sqrt{\pi }$.

The following expression gives the Griffith’s criterion, as follows:

$${\sigma }_c\sqrt{a_0}=\sqrt{{\displaystyle \frac{2\gamma E}{\pi }}}\to K_c=\sqrt{2\gamma E}$$

(8)

where E is the Young’s modulus of the material and γ is the surface energy of the material.

Both Young’s modulus and the surface energy are constant, from which we expect that the product of $\sqrt{a_0}$ and σ_c should also remain a constant. To test Griffith’s criterion, we gradually increased the initial crack size; the critical stress is expected to decrease because the product should remain constant. Of course, the value should remain constant throughout the test. However, the stress intensity factor in the simulation is expected to converge to a specific value. If this value matches the experimental one, then the potential is suited for simulating graphene with crack defects [68].

In Figure 4b, inner carbon atoms in the bulk material are connected to three other atoms, while the atoms at the surface are only connected to two atoms. Since the atoms at the surface can share electrons and form bonds, they are less stable and have more energy than the inner atoms. Upon cutting the material, as in Figure 4c, two new surfaces are created with more unstable and energetically active atoms. It is anticipated that the system presented in Figure 4c will possess more energy compared to system (as shown in Figure 4a), with the energy difference being equivalent to the energy expended in dividing the system into two slabs (assuming zero friction). At this point, we have our values and are ready to compare them with the experimental ones, so we can estimate which potential is best suited for a system with cracks [69].

Table 2. The stress intensity factor (SIF) convergence and other parameters like Young’s modulus and the surface energy.

Cracked Graphene Data for Different Potentials
Potential Name	Crack’s Length a0 (nm)	Fracture Stress σc (GPa)	$\mathbf{Stress}\mathbf{Intensity}\mathbf{Factor}(\mathit{M}$$\mathbf{Pa}\sqrt{\mathbf{m}}$$){\mathit{K}}_{\mathit{c}}={\sigma }_c\sqrt{\mathit{\pi }{\mathit{a}}_0}$	Young’s Modulus E (GPa)	Surface Energy $\mathit{\gamma }\left(\frac{\mathbf{J}}{{\mathbf{m}}^2}\right)$	$\mathbf{SIF}=\sqrt{2\mathit{E}\mathit{\gamma }}$$\left(\mathbf{MPa}\sqrt{\mathbf{m}}\right)$
CH.airebo-m	0.615	70.39	3.09	978.72	8.859	4.16
	1.107	57.87	3.41	1001.38		4.21
	1.599	50.48	3.58	952.17		4.11
	2.091	51.32	4.16	936.53		4.07
	2.583	45.42	4.09	899.32		3.99
CH.rebo	0.615	60.92	2.65	873.14	9.985	4.18
	1.107	56.33	3.23	872.51		4.17
	1.599	53.41	3.51	851.74		4.12
	2.091	49.12	3.52	820.72		4.05
	2.583	43.66	3.61	794.49		3.98
C.meam	0.615	60.92	2.68	1002.78	8.674	4.17
	1.107	56.33	3.32	983.72		4.13
	1.599	53.41	3.79	956.01		4.07
	2.091	49.12	3.98	933.93		4.03
	2.583	43.66	3.93	905.29		3.96
BNC	0.615	87.37	3.84	988.04	1.573	1.76
	1.107	73.22	4.32	1001.08		1.77
	1.599	69.16	4.90	972.47		1.75
	2.091	66.71	5.41	971.18		1.75
	2.583	62.33	5.61	910.52		1.69
SiC-b	0.615	98.73	5.53	1054.56	2.852	2.45
	1.107	81.75	5.43	1022.52		2.41
	1.599	73.75	5.23	999.83		2.39
	2.091	66.94	4.82	956.89		2.34
	2.583	61.39	4.34	947.21		2.32

shows that the surface energy and the critical stress depend significantly on the type of potential used in the simulation, even if the potential yields an excellent value for Young’s modulus. The BNC potential’s SIF value does not converge to the experimental value, which is 4 MPa$\sqrt{\mathrm{m}}$ [68,69,70,71,72,73], even though it was one of the best potentials to simulate the pristine graphene. Tersoff potentials do not produce excellent data because they do not account for the long-range interaction, which is required in the case of the pre-cracked atoms, where the atoms across the fracture could interact with each other. Figure 5 and Figure 6 exhibit that the CH.airebo-m potential agrees with the experimental data, since it accounts for the long-range interaction, including the Lennard–Jones or Morse potential. The C.meam agrees with experiments because it does not use a switching function for the cutoff distance and uses the electron density to estimate the system’s energy [74,75]. The CH.rebo potential produces better results than the Tersoff potential but is not as accurate as the CH.airebo-m potential, since it accounts only for the conjugation and radical characteristic, but not the long-range interaction.

4. Conclusions

The accuracy of the MD simulation is greatly affected by the type of interatomic interactions, namely, the force-field. Depending on the structure and configuration of the simulated material, different types of potentials may be utilized for various applications. The CH.rebo and BNC.tersoff potentials are best suited for simulating monolayer pristine graphene, where the long-range interaction is less significant in this case. The C.meam and the CH.airebo-m potentials are recommended for the graphene with crack defects, since they use the electron density and include the long-range interaction term, respectively. We recommend AIREBO-m for a general case MD study for mechanical properties. Although most of the potentials do not agree well with the experimental results of the mechanical properties for pristine or pre-cracked graphene, it is possible that these potentials might be suitable for the other properties and/or structures.