Molecular Simulation

ISSN: 0892-7022 (Print) 1029-0435 (Online) Journal homepage: http://www.tandfonline.com/loi/gmos20

Efficient embedded atom method interatomic potential for graphite and carbon nanostructures V. E. Zalizniak & O. A. Zolotov To cite this article: V. E. Zalizniak & O. A. Zolotov (2017): Efficient embedded atom method interatomic potential for graphite and carbon nanostructures, Molecular Simulation, DOI: 10.1080/08927022.2017.1324957 To link to this article: http://dx.doi.org/10.1080/08927022.2017.1324957

Published online: 30 May 2017.

Submit your article to this journal

View related articles

View Crossmark data

Full Terms & Conditions of access and use can be found at http://www.tandfonline.com/action/journalInformation?journalCode=gmos20 Download by: [The UC San Diego Library]

Date: 31 May 2017, At: 03:47

Molecular Simulation, 2017 https://doi.org/10.1080/08927022.2017.1324957

Efficient embedded atom method interatomic potential for graphite and carbon nanostructures V. E. Zalizniak and O. A. Zolotov Institute of Mathematics and Computer Science, Siberian Federal University, Krasnoyarsk, Russia

ABSTRACT

A new interatomic potential for graphite and graphene based on embedded atom method is proposed in this paper. Potential parameters were determined by fitting to the equilibrium lattice constants, the binding energy, the vacancy formation energy and elastic constants. The agreement between the calculated properties of graphite and experimental data is very good. In addition, the proposed potential quite accurately reproduces the surface energy of graphite and the binding energies of carbon atom in fullerene C60 and in SWNTs. The proposed potential is computationally more efficient than the existing widely used carbon potentials. It is intended for use in large-scale molecular dynamics simulations of carbon structures.

1. Introduction Carbon nanostructures (graphite, graphene, nanotubes and fullerenes) attract a wide interest because they possess special properties. These properties open a wide field of uses ranging from structural materials, and quantum dots, to implanted electronics. For instance, nanotube and graphene structures can be used to construct miniaturised super fast computers and transistors by reducing conductive regions, power consumption and heat accumulation. In spite of considerably increased computer performance, the application of ab initio methods for simulation of atomistic structures is still limited to relatively small systems of atoms and relatively short simulation times. By contrast, the method of molecular dynamics based on empirical interatomic potentials makes it possible to simulate much larger systems for much longer times. For this reason there is a demand for realistic and efficient interatomic potentials. A variety of carbon interatomic potentials have been developed. Examples of the carbon potentials are reactive empirical bond-order potential (REBO) developed by Brenner [1], environment dependent interatomic potential (EDIP) developed by Marks [2], reactive force field (ReaxFF) developed by van Duin et al. [3] and long-range carbon bond order potential (LCBOP) developed by Los et al. [4,5]. Zhou et al. [6] have developed an analytical bond order potential (BOP) for carbon, and performed studies to compare the proposed BOP, REBO, EDIP and ReaxFF carbon potentials with the available experimental data. All these potentials include explicit threebody interactions.

CONTACT  V. E. Zalizniak 

[email protected]

© 2017 Informa UK Limited, trading as Taylor & Francis Group

ARTICLE HISTORY

Received 23 December 2016 Accepted 25 April 2017 KEYWORDS

Interatomic potential; embedded atom method; graphite; graphene; carbon nanotube

The modified embedded atom method (MEAM) potential, originally proposed by Baskes et al. [7], was parameterised for carbon atoms by Uddin et al. [8]. The purpose of this paper is to present simple and efficient potential for carbon atoms in graphite and graphene based on embedded atom method (EAM) [9,10]. EAM potentials provide an adequate description of metallic systems and computational speed essential for computer simulations of complex materials science problems. They are much more computationally efficient as compared to three-body potentials and MEAM potentials [11]. In Section 2 of the paper, we briefly consider embedded atom method. Basic structural and physical parameters of graphite are presented in Section 3. Sections 4 and 5 describe the interatomic potential for graphite and a strategy of parameterisation. In Section 6, we present the results of fitting and testing the potential. In Section 7, we summarise our results and make some conclusions.

2.  Embedded-atom method The potential developed in this work is based on the formalism of the embedded-atom method (EAM) [9, 10]. In the framework of EAM, the total energy of a system can be written as

Etot =

N N ∑ ( ) ( ) ( ) 1∑ 𝜌 rnm (1) 𝜑 rnm , 𝜌n = En , En = F 𝜌n + 2 m=1 m=1 n=1

N ∑

m≠n

m≠n

here Etot is the total energy of the system of N atoms, En is the potential energy associated with atom n, ρn – the electron density

2 

 V. E. ZALIZNIAK AND O. A. ZOLOTOV

at atom n due to all other atoms, ρ (rnm) – the contribution to the electron density at atom n due to atom m at the distance rnm from atom n, F(ρn) – the embedding energy of the atom into the electron density ρn, φ(rnm) – the two body central potential between atoms n and m separated by rnm. Interpretation and functional form of φ(r), ρ(r) and F(ρ) depend on a particular method. The popularity of the EAM model results from its quantum mechanical justification, as well as its mathematic simplicity, which makes this model conducive to large-scale computer modelling. In practical applications of EAM potentials, it is also desirable to employ a switching function in order to terminate the potential and forces smoothly at the cut-off distance because the energy conservation is sensitive to the truncation of the force field. For this purpose, a switching function fc(r) can be applied to the electron density distribution and to the pair potential in a region just below the cut off distance rc:

𝜑(r) → 𝜑(r)fc (r), 𝜌(r) → 𝜌(r)fc (r)

3.  Properties of graphite The physical properties of graphite are a consequence of its layer structure. Carbon atoms within a plane form a two-dimensional hexagonal lattice (graphene). Three-dimensional graphite is obtained by a stacking of the graphene planes in a …ABAB… sequence. The ideal graphite structure is shown in Figure 1. This stable hexagonal lattice has the lattice constant a, and the interlayer distance c. Strong chemical bonding exists within the layer planes. The weaker bonds between layers are often explained to be the result of van der Waals forces. The layer structure of graphite is responsible for the anisotropies in its physical properties. Structural and some physical parameters of graphite are shown in Table 1. We also use the experimental value of the interlayer binding energy of highly oriented pyrolytic graphite εb = 0.031 eV/atom [12].

Figure 1. Schematic of graphite structure.

Table 1.  Structural and physical parameters of graphite; c11, c12, c13, c33, c44 experimental elastic constants; Egraphite – experimental cohesive energy per atom; Evf - unrelaxed vacancy formation energy. Parameter a, Å c, Å c11, eV/Å3 c12, eV/Å3 c13, eV/Å3 c33, eV/Å3 c44, eV/Å3 Egraphite, eV Evf, eV

Experimental or DFT value 2.463a [18] 3.356 [18] 6.62 [18], 6.91 [19] 1.11 [18], 0.87 [19] 0.09 [18], 0 [19] 0.22 [18], 0.24 [19] 0.025 [18], 0.03 [19] 7.37 [20] 7.6b [21]

a

The interatomic distance within a layer plane is 1.42 Å. Result of ab initio calculations.

b

4.  Potential for graphite Using (1), the potential energy of carbon atom n is written as

En = Eg,n + Eintl,n , � � 1� � � 𝜑 r + F 𝜌n , = 2 m=1 0 nm

(2)

m≠n

� � � � Eintl,n = 𝜑1 ⟨zn, 1 ⟩ + 𝜑1 ⟨zn, −1 ⟩ , where Eg,n is the energy of the in-plane bonding in graphene, Eintl,n is the inter-layer binding energy, 〈zn, ±1〉 is the average local separation between graphite layers at the location of atom n defined as

∑N ⟨zn,p ⟩ = � � fg r, z, p =

�

1, 0,

� � − zm �fg rnm , zn − zm , p , � ∑N � m=1 fg rnm , zn − zm , p

(4)

𝜌(r) = 𝜌0 fc (r)(1 + 𝛽r)2 e−𝛼r

N

Eg,n

The electron density distribution is taken in the following form [13] This approximation of electron density distribution was suggested from the basic principles of quantum mechanics. The form of this distribution defines the two body potential [13] 6 ∑

𝜑0 (r) = 𝜀fc (r) exp (−𝛼r)

where parameters an depend on α and β. The embedding energy function F(ρ) is assumed to be in the following form:

{

0.9c ≤ r ≤ 1.2c , 0.7c ≤ �z� ≤ 1.3c , pz > 0 otherwise.

The interlayer potential φ1〈zn, p〉 describes the interaction between a carbon atom and the nearest graphite layer.

(5)

n=−1

m=1 �zn

(3)

an (𝛼r)n

∗

∗ 2

F(𝜌) = c0 + c1 𝜌 + c2 (𝜌 ) +

c3 (𝜌∗ )3 , ∗ 3

c4 (𝜌 ) ,

𝜌 ≤ 𝜌e 𝜌 > 𝜌e

(6)

𝜌 𝜌∗ = − 1, 𝜌e where ρe is the equilibrium electron density and cn are some coefficients. The following simple polynomial switching function is used

MOLECULAR SIMULATION 

rsw ≤ r ≤ rc

(7)

where rsw is the distance at which the switching function is applied. The values of coefficients d1, d2, d3 follow from the conditions

( ) d2 fc ( ) df ( ) rsw = 0. fc rsw = 1, c rsw = 0, dr dr 2 We set rsw = 2.5 Å and rc = 3 Å. Then parameters of the switching function are d1 = −2160, d2 = −19,440, d3 = −46,656. The following form of the inter-layer potential is adopted ) ( z −1 𝜑1 (z) = −𝜀b (1 + x) exp (−x) , x = 𝛾 (8) c

5.  Parameterisation procedure In order to define the potential of interaction between carbon atoms one need to fit two parameters of the electron density distribution α and β, coefficients of the embedding energy function (6) and parameter γ in (8). The condition of lattice equilibrium and properties of graphite can be calculated from (2) [10]. These properties include the binding energy per atom, the unrelaxed vacancy formation energy and elastic constants. Elastic constants c11, c12, c44 are defined by Eg,n but elastic constants c13, c33 are defined by Eintl,n. The approximate value of c33 is given by the following relation 2 𝜀b 𝛾 2 c 2 d 𝜑1 = (c) Va dz 2 Va

( ) 𝜕F ( ) 1 𝜕2F ( ) 𝜌e , c2 = 𝜌2e 2 𝜌e c0 = F 𝜌e , c1 = 𝜌e 𝜕𝜌 2 𝜕𝜌

(10)

Coefficient c3 is calculated from the condition F(0) = 0, and coefficient c4 is determined so that the proposed potential would adequately reproduce the binding energies of small fullerenes. The results of fitting are presented in Table 2. The value of parameter γ = 8.08 is calculated from relation (9), using experimental value of elastic constant c33.

6. Results 6.1.  Graphite properties The calculated properties of graphite are compared with the experimental values, to which they were fitted. The results are presented in Table 4. The match between experiment and the proposed EAM model is almost perfect. Functions φ(r), ρ(r) and F(ρ) of the proposed potential are shown in Figures 2–4. Pair potential function and embedding energy function have a positive curvature for the equilibrium lattice. This ensures the stability of graphene structure. In the presented model, the value of binding energy per carbon atom in graphene Egraphene = Egraphite – 2εb = 7.308 eV. As additional characteristic of the potential, the surface energy of (0 0 0 1) surface of graphite was calculated by taking

(9) 0.4

where undeformed atomic volume for carbon atom Va = 8.8 Å3. The experimental data used in fitting procedure consist of the equilibrium lattice constants, the cohesive energy, the vacancy formation energy and five elastic constants. The values of these constants were taken as the average of values given in Table 1. They are presented in the last row of Table 4. All computations are performed for zero temperature with hydrostatic pressure set to zero. General description of parameterisation procedure is given in [13]. Here, the fitting procedure suggests that the equilibrium lattice constant, the cohesive energy, the vacancy formation energy and elastic constant c11 are reproduced exactly. The

3

(a) c33 =

optimal pair (α, β) provides the minimal discrepancy between calculated and experimental values of elastic constants c12 and c44. To calculate coefficients cn of embedding energy function (6) we require five conditions. First of all there are values of ρe, F(ρe), F′(ρe) and F″(ρe) that correspond to the optimal pair (α, β). Then

Electron density , e/A

⎧ 1 , r < rsw ⎪ fc (r) = ⎨ d1 (r ∗ )3 + d2 (r ∗ )4 + d3 (r ∗ )5 , ⎪ 0, r > r ⎩ c r r ∗ = − 1, rc

 3

0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0

0.5

1

1.5

r,A

2

2.5

3

3.5

Figure 2. Electron density as a function of distance.

Table 2. Parameters of the proposed potential for graphite and graphene. Parameters of the atomic electron density distribution ρ0, e/Å3 α, 1/Å β, 1/Å

0.11232 3.0845 6.247

Parameters of pair potential ε, eV a-1 a0 a1 a2 a3 a4 a5 a6

314.6943 1 0.6502 0.3484 −1.9225E-03 −1.1794E-02 −1.3239E-03 −9.8169E-05 −3.8489E-06

Parameters of the embedding energy function ρ0, e/Å3 0.5063 c0, eV −6.6927 c1, eV −6.2968 c2, eV 1.3888 c3, eV 0.9929 c4, eV 498

4 

 V. E. ZALIZNIAK AND O. A. ZOLOTOV

Table 3. Binding energies of carbon atom in fullerenes.

4

Pair potential , eV

3

Fullerene C24 C60

2 1

Energies calculated from (2), eV/atom −6.45 −6.91

Experimental and/or DFT energies, eV/atom −6.45b [22] −(6.94 – 6.98)a [23], −7.06b [22]

a

Experimental result. DFT result.

b

0 −1 −2 1.4

1.6

1.8

2

2.2

2.4 2.6 r,A

2.8

3

3.2

3.4

Figure 3. Pair interaction energy as a function of distance between atoms.

1

Embedding energy , eV

0 −1 −2 −3 −4 −5 −6 −7 0

0.1

0.2

0.3 0.4 Electron density , e/A3

0.5

0.6

Figure 4. Embedding energy as a function of background electron density.

the energy difference between the total energy of a periodic slab and an equivalent bulk reference amount: ) ( 1 Esurf = Eslab − N ⋅ Egraphite 2S where Eslab is the total energy of an N-atom slab, Egraphite is the energy of one carbon atom in the bulk, S is the area of the slab surface, and the factor 1/2 accounts for the two surfaces in the slab. Calculations were done using the equilibrium lattice parameters. In the proposed model the force acting on a carbon atom from adjacent graphene sheets is equal to zero at the equilibrium state. Then there is no relaxation or reconstruction at the surface. Calculated graphite (0 0 0 1) surface energy is 0.188 J/m2. Experimental values of pyrolytic graphite surface energy are in the 0.1–0.2 J/m2 range [14,15].

Figure 5. The value of binding energy per carbon atom in a bare SWNT versus tube radius. Circles are results calculated using potential (2). Squares are results of ab initio calculations [17]. The solid line corresponds to Egraphene . Table 4.  Graphite properties calculated with the use of various potentials and ­experimental properties of graphite. Potential c11, eV/Å3 c12, eV/Å3 c13, eV/Å3 c33, eV/Å3 c44, eV/Å3 Egraphite, eV Evf, eV ReaxFF REBO EDIP LCBOP BOP [6] EAM (2) Exp

17.2 8.4 6.5 6.55 6.0 6.775 6.775

15.2 3.1 1.2 1.03 3.1 0.98 0.98

0 0 – – 0 0 0

0.2 0 – 0.23 0 0.23 0.23

The first term in (2) presents the potential for graphene. It was applied to calculate the binding energy of carbon atom in zigzag single wall carbon nanotubes (SWNT’s) and fullerenes. Calculated binding energies of carbon atom in fullerenes in comparison with experimental result and results of DFT calculations are presented in Table 3. Results of ab initio calculations of the binding energy of carbon atom in zigzag SWNTs were presented [16], showing curvature effects. Zigzag SWNTs are rolled graphene sheets that are characterised by integer n defining the rolling vector of

9.06 7.81 – 7.37 7.4 7.37 7.37

9.5 8.9 – 7.9 5.5 7.6 7.6

graphene [17]. The nanotube radius R = 0.5a⋅n/π. It was found that the value of binding energy per carbon atom in a bare SWNT depends on the nanotube radius as

ESWNT = Egraphene + 6.2.  Properties of carbon nanotubes and small fullerenes

0 0 – – 0 0.0279 0.0275

2.41 R2

(11)

for 4 < n < 15. Calculated binding energies of carbon atom in zigzag SWNTs in comparison with the results of ab initio calculations are presented in Figure 5. The maximum discrepancy is found to be less than 3.5%. Graphite properties obtained from different potential models are presented in Table 4. ReaxFF and REBO potentials significantly overestimate the experimental and DFT values. BOP [6] improves over ReaxFF and REBO potentials. Only LCBOP rather closely reproduces the experimental and DFT values. This potential is rather complexly built and includes long-range interactions with the cutoff at 6 Å. All mentioned potentials include explicit

MOLECULAR SIMULATION 

three-body interactions. The inclusion of more detailed effects comes at additional computational cost in comparison with the proposed potential. 6.3.  Comparison of computational costs Empirical many-body potentials for carbon may be represented in the form containing two- and three-body interactions separately. Binding energy Eb is written as a sum over nearest neighbours in the form ∑∑[ ( ) ( )] Eb = UR rnm − bnm UA rnm . n m>n

The functions UR(r) and UA(r) are pair-additive interactions that represent repulsion and attraction, respectively. The quantity rnm is the distance between pairs of nearest-neighbour atoms n and m. The term bnm determines the bond strength and it is a monotonically decreasing function of the number of competing bonds, the strength of the competing bonds and the cosines of the angles with competing bonds. To benchmark the relative computational cost of implementation of various potentials we compute energy of carbon atoms in graphene layer. The computational cost of implementation of the REBO potential is comparable to that of the proposed EAM potential because the REBO potential includes only first-neighbour interactions. The proposed EAM potential is approximately three times more efficient than bond order potential [6]. LCBOP [5] is even less computationally efficient then BOP [6] because it is rather complexly built and includes not only short-range interactions but also the long-range nonbonded interactions.

7. Conclusions This paper presents a new EAM potential model to describe graphite and graphene. The potential model has a simple function form and it is easy to use in computer simulations. The potential parameters were determined by fitting to the equilibrium lattice constants, the binding energy, the vacancy formation energy and elastic constants. The agreement between the calculated properties of graphite and experimental data is very good. In addition, the proposed potential quite accurately reproduces the surface energy of graphite and the binding energies of carbon atom in fullerene C60 and in SWNTs. The potential does not capture the properties of diamond because EAM cannot properly describe the anisotropic spatial distribution of the electron density in diamond. The proposed potential is computationally more efficient than the existing widely used carbon potentials, and it is intended for use in large-scale simulations. The proposed potential is fully compatible with EAM potential [13] parameterised for 27 metals. As a first approximation, one can assume that the form of electron density distribution of an atom depends only slightly on surrounding atoms, and the embedding energy function is the unique property of an atom because it depends only on electron density due to all other atoms. Then the proposed potential for carbon atom allows one to create a parameter database for binary carbon-metal interactions. All one need to do is to fit the parameters of pair cross-interaction functions against experimental or ab initio data for the relevant carbide system. This database will allow

 5

us to simulate various carbon-metal systems such as graphite intercalation compounds, endohedral metallofullerenes and metallacarbohedrenes. We believe that the proposed EAM potential will find ­applications in MD studies of various carbon structures.

Disclosure statement No potential conflict of interest was reported by the authors.

References   [1] Brenner DW, Shenderova OA, Harrison JA, et al. A second generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. J Phys Condens Matter. 2002;14:783–802.  [2]  Marks NA. Generalizing the environment-dependent interaction potential for carbon. Phys Rev B. 2001;63:035401.   [3] van Duin ACT, Dasgupta S, Lorant F, et al. ReaxFF: a reactive force field for hydrocarbons. J Phys Chem. 2001;105:9396–9409.   [4] Los JH, Fasolino A. Intrinsic long-range bond-order potential for carbon: performance in Monte Carlo simulations of graphitization. Phys Rev B. 2003;68:024107.  [5] Los JH, Ghiringhelli LM, Meijer EJ, et al. Improved long-range reactive bond-order potential for carbon. I. Construction. Phys Rev B. 2005;72:214102–14.   [6] Zhou XW, Ward DK, Foster ME. An analytical bond-order potential for carbon. J Comput Chem. 2015;36:1719–1735.   [7] Baskes MI. Modified embedded-atom potentials for cubic materials and impurities. Phys Rev B. 1992;46:2727–2742.   [8] Uddin J, Baskes MI, Srinivasan SG, et al. Modified embedded atom method study of the mechanical properties of carbon nanotube reinforced nickel composites. Phys Rev B. 2010;81:104103.   [9] Daw MS, Baskes MI. Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals. Phys Rev Lett. 1983;50:1285. [10]  Daw MS, Baskes MI. Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals. Phys Rev B. 1983;29:6443–6453. [11] Plimpton SJ, Thompson AP. Computational aspects of many-body potentials. MRS Bull. 2012;37(5):513–521. [12] Liu Z, Liu JZ, Cheng Y, et al. Interlayer binding energy of graphite: a mesoscopic determination from deformation. Phys Rev B. 2012;85:205418–5. [13] Zalizniak VE, Zolotov OA. Towards a universal embedded atom method interatomic potential for pure metals. J Sib Federal Univ Math Phys. 2015;8:230–240. Available from: http://elib.sfu-kras.ru/ handle/2311/16811. [14]  Eustathopoulos N, Nicholas MG, Drevet B. Pergamon materials series: wettability at high temperatures. Amsterdam: Elsevier Science; 1999. [15] Zaidi H, Robert F, Paulmier D. Influence of adsorbed gases on the surface energy of graphite: consequences on the friction behavior. Thin Solid Films. 1995;264:46–51. [16]  Gülseren O, Yildirim T, Ciraci S. Systematic ab initio study of curvature effects in carbon nanotubes. Phys Rev B. 2002;65:153405–4. [17] Dresselhaus MS, Dresselhaus G, Eklund PC. Science of fullerenes and carbon nanotubes. San Diego (CA): Academic Press; 1996. [18] Bosak A, Krisch M, Mohr M, et al. Elasticity of single-crystalline graphite: inelastic X-ray scattering study. Phys Rev B. 2007;75:153408-4. [19]  Blackslee OL, Proctor DG, Seldin EJ, et al. Elastic constants of compression‐annealed pyrolytic graphite. J Appl Phys. 1970;41:3373– 3382. [20] Brewer L. The cohesive energies of the elements; 1977. (Lawrence Berkeley National Laboratory Reports no. LBL-3720). [21] Li L, Reich S, Robertson J. Defect energies of graphite: densityfunctional calculations. Phys Rev B. 2005;72:184109–10. [22] Hua X, Çağın T, Che J, et al. QM(DFT) and MD studies on formation mechanisms of C60 fullerenes. Nanotechnology. 2000;11:85–88. [23] Chen HS, Kortan AR, Haddon RC, et al. Reactivity of C60 in pure oxygen. Appl Phys Lett. 1991;59:2956–2958.