Ab Initio Study of the Electronic and Energy Properties of Diamond Carbon

Abstract

In this chapter, we present a study on the electronic properties of diamond carbon, using band structure and density of states calculations. The calculations are based on the use of the grid-based projector-augmented wave (GPAW) and atomic simulation environment (ASE) methods. The main results of our work are the optimization of diamond energy (to −17.57 eV) and the calculation of the gap with the PBE (Perdew, Burke, and Ernzerhof) and the functional hybrid PBE0 hybrid functional, which is about 5.368 eV (the closest value to the value found in the literature). We were also able to reproduce the experimental value of the lattice constant of diamond to within 0.2% for PBE0 and 0.4% for PBE. Our results contribute to the study of the electronic properties of diamond using GPAW and ASE simulation, which is a set of Python modules, designed to facilitate the setup, execution, and analysis of atomic/electronic calculations. This tight integration of ASE and GPAW should be exploited in future research of the electronic properties of diamond, which is one of the most promising materials for the integrated electronic and photonic, radio, optoelectronic, and quantum devices industry. This chapter provides interesting information for the theoretical and experimental communities working in this field.

1. Introduction

Carbon is a unique element that appears in many different configurations of electronic states [1]. Materials such as graphite, graphene, carbon nanotubes, and diamond have basically the same chemical composition but a different crystal structure. This is called allotropy of carbon (see Figure 1 [2, 3]). Carbon has six electrons occupying the 1s², 2s², and 2p² orbitals. The 1s² electrons are two strongly bound core electrons, while the 2s² and 2p² electrons are weakly bound valence electrons. The valence electrons give rise to the 2s, 2p_x, 2p_y, and 2p_z orbitals in the crystal phase. These orbitals are important for the formation of covalent bonds in carbon materials. Since the energy difference between the upper 2p and lower 2s orbitals is relatively small compared to other group 4 elements such as silicon (Si), the electron wave functions can easily mix with each other, which is called hybridization [4]. In spⁿ hybridization, the (n + 1) bonds are responsible for the local structure. Thus, in sp. hybridization, a one-dimensional chain structure is formed, which is called carbyne. Carbyne [5] is a chain of carbon atoms held together by alternating double or single and triple atomic bonds. This makes it a true zero-dimensional material, unlike graphene atomic sheets that have a top and bottom, or hollow nanotubes that have an inside and outside. According to calculations by Liu et al., carbyne has twice the tensile stiffness of graphene and carbon nanotubes and nearly three times that of diamond.

Figure 1.

Different allotropic forms of carbon [2, 3].

An sp² hybridization forms a two-dimensional planar structure called graphene. In sp³ hybridization, a regular three-dimensional tetrahedral structure is formed, which is known as diamond.

In graphite, in addition to the three in-plane bonds, out-of-plane π orbitals are formed with highly delocalized electrons that contribute primarily to electrical and thermal conductivity. The interaction of the free π electrons with light makes graphite appear black. In addition, the weak van der Waals interaction between the graphite sheets makes the graphite flexible and cleavable, as the graphene sheets can easily move relative to each other. A carbon nanotube can be viewed as a rolled-up graphite sheet and thus a one-dimensional structure. The bonds in the nanotubes are essentially of the sp² type. However, the circular curvature causes quantum confinement and σ − π hybridization in which the three σ bonds are slightly out of plane. In this part, we are interested in the study of diamond.

Diamond has a crystal structure consisting of a face-centered cubic lattice of which one tetrahedral site out of two is occupied by a carbon atom, as represented in Figure 2 [6].

Figure 2.

Elementary crystal lattice cell of diamond [6].

The elementary lattice cell of diamond has carbon atoms, and its parameter is 3.57 Å. The shortest distance between two carbon atoms, that is, the distance corresponding to the quarter of the diagonal of the mesh is 1.54 Å, which makes the crystal lattice of diamond very dense (1.76 × 1023 atoms cm⁻³) in comparison with other materials. The compactness of its structure, associated with the high energies of the carbon-carbon bonds (E_(C-C) = 360 kJ mol⁻¹), is at the origin of diamond mechanical properties, which are well known (see Table 1) [6].

Property	Unit	Silicon	Diamond
Density	g cm^−3	2.33	3.52
Melting point	K	1687	3773
Hardness	GPa	8.5	57–101
Coefficient of thermal expansion	K^−1	2.6 × 10^−6	1.0 × 10^−6
Young’s modulus	GPa	130–180	1050–1200
Thermal conductivity	W m^−1 K^−1	148	2500
Prohibited bandwidth GAP	eV	1.12	5.45
Electron mobility	cm^2 V^−1 s^−1	1350	4500

Table 1.

Some physical properties of diamond compared to silicon.

Its high hardness, estimated at 7000 kg mm⁻² (Knoop hardness) [6], makes diamond the most resistant material. This characteristic allows it to be used in many mechanical applications such as cutting, drilling, or polishing. Its very high Young’s modulus makes it advantageous for applications involving microelectromechanical system oscillators (MEMS)-type resonators [6]. In addition to its high thermal conductivity [7], which is five times greater than that of copper, diamond is used as a heat sink for electronic or optoelectronic components such as power transistors and laser diodes. More recently, it has been used for the realization of optical windows for high-power lasers or X-ray tubes; it is the combination of its optical, mechanical, and thermal properties, which is decisive [6].

Diamond is a three-dimensional crystalline material. Because of its structure, it is considered one of the most promising materials for the integrated electronic and photonic, radio, optoelectronic, and quantum device industry.

Following our publications on the electronic and magnetic structures of graphite using the Korringa-Kohn Rostoker (KKR) method [8] and on the electronic properties of the adsorption of CO and N₂ molecules on graphene and on a graphite plate using GPAW and ASE [7], we decided to set up a simulation system based on GPAW and ASE to study the electrical and energetic properties of diamond. The main results of our work are the optimization of diamond energy, as well as the calculation of the density of state and the band structure.

2. Calculation methods

In this chapter, we present a real-space implementation of the PAW method in the open-source software package GPAW [9]. Note that there are several software packages that currently implement the PAW method using a plane-wave basis [10, 11, 12]. The PAW method [13, 14] is formally an all-electron method that provides an exact transformation between smooth pseudo wave functions and all-electron wave functions. While in practice the PAW method resembles pseudo potential methods, it remedies several shortcomings of conservative norm or ultra-smooth pseudo-potentials.

The GPAW code is built on top of the ASE [3], which is a set of Python modules designed to facilitate the setup, execution, and analysis of atomistic/electronic calculations. This tight integration of ASE and GPAW should be maintained in the future. Interest in ASE has grown considerably over the past few years, to the point where ASE now supports about 12 different force and energy calculators.

There are of course several other open-source projects focused on DFT and time-dependent (TDDFT), such as abinit, quantum espresso and KKR [15], and Octopus. How does GPAW fit into this code market? The main feature that distinguishes GPAW from its competitors is the combination of a real-space description with the PAW method. The PAW method allows an accurate description of essentially all electrons in the frozen core, leading to smooth pseudo wave functions even for transition metals. The real-space description allows for easy and highly scalable parallelization due to the real-space decomposition, which allows for accurate calculations even for large systems.

Among the features currently implemented in GPAW, we can mention the calculation of static response functions using density-functional perturbation theory and the more general calculations of dynamic response functions in TDDFT. Force calculations as well as adiabatic and Ehrenfest dynamics are also implemented in TDDFT. In addition, work is underway to expand the number of atomic configurations to include all elements up to atomic number 86 [16].

In this work, we present density of states and band structure calculations of diamond using the PAW method as implemented in the GPAW package (https://wiki.fysik.dtu.dk/gpaw).

Traditionally [17], software packages for this type of simulation have been implemented in compiled languages, where Fortran in its various iterations has been the most popular choice. Although interpreted dynamic languages, such as Python, can increase the efficiency of the programmer, they cannot directly compete with the raw performance of compiled languages.

However, by using an interpreted language along with a compiled language, it is possible to benefit from most of the productivity-enhancing features while still achieving good numerical performance. This approach has been implemented in the GPAW electronic structure simulation software, using a combination of Python and C programming languages.

3. Results and discussions

3.1 Optimization of the diamond structure

The diamond unit cell contains two carbon atoms located at (0.0, 0.0, 0.0) and (0.25 a, 0.25 a, 0.25 a). After optimization, the CC bond length was found to be 1.55 Å. Relaxation of all model systems was performed using an 8 × 8 × 8 Monkhorst-Pack grid sampling, used during Brillouin zone integrations (512 k-points in the irreducible part of the Brillouin zone). Tables 2 and 3 shows the results of the diamond structure optimization.

Diamond energy	−17.57 eV
Diamond atomization energy	−17.57 eV
Experimental bond length (GPAW and ASE):
Diamond energy	−17.57 eV
Bond length	1.55 Å
PBE energy minimum
Diamond energy	−17.57 eV
Bond length	1.55 Å

Table 2.

Optimization of the diamond structure.

Valence states	Energy (eV)	Radius (Å)
2s (2.00)	−13.751	0.635
2p (2.00)	−5.284	0.635
*s	13.461	0.635
*p	21.927	0.635
*d	0.000	0.635

Table 3.

Valence states, energy, and radius.

Other results provided by GPAW:

C-setup:

Name: Carbon.

Z: 6.

Valence: 4.

Core: 2.

Charge: 0.0.

Cutoffs: 0.64(comp), 1.14(filt), 1.14(core), lmax = 2.

3.2 Electronic band structure of diamond

Figure 3 shows the electronic band structure of diamond, whose unit cell is composed of two carbon atoms, located at the positions shown in Section 3.1. The band structure is plotted along the high-symmetry line ΓXWKL of the Brillouin zone. The x-axis represents the kpts = {‘path’ = ‘ΓXWKL’, ‘npoints’ = 60}, with 16 bands, 512 k-points: 8 × 8 × 8 Monkhorst-Pack grid, 60 k-dots in the irreducible part of the Brillouin zone and Coarse grid: 12 × 12 × 12 grid as well as the fine grid: 24 × 24 × 24 grid. The ordinate axis represents the energy in (eV). It is clear from Figure 3 that diamond is an insulator with a non-zero band gap since the valence and conduction bands do not meet at the k-point. This corresponds to a zero density of state when plotting the density of state as a function of energy. The band structure is then characterized by a nonlinear dispersion of the bands in the vicinity of the Fermi energy, whereas in graphene, for example, there is a linear dispersion of the bands in the vicinity of the Fermi energy [7].

Figure 3.

Band structure of diamond.

Next, we calculated the electronic band structure of diamond along the high-symmetry directions in the Brillouin zone [18]. Appendix 1 lists the calculation steps.

First, we performed a standard ground state calculation and saved the results in a .gpw file. As we are dealing with a small bulk system, the plane-wave mode is the most appropriate here.

Next, we calculated the eigenvalues along a high-symmetry path in the Brillouin zone

.

See

for the definition of the special points for a face-centered cubic (FCC) lattice.

For the band structure calculation, we fixed the density to the previously calculated ground state density, and as we wanted to calculate all k-points, we did not use symmetry (

).

Finally, we platted the band structure using ASE’s band-structure tool (

):

3.3 Density of state of diamond

As mentioned, the lattice constant used for the calculations is a = 3.57 Å, and the void used in the three-dimensional (3D) structures along the axes is equal to 0 Å. The convergence criteria used imply a total energy change of 0.0005 eV/electron, a maximum integral of absolute density change of 0.0001 electrons, and a maximum integral of absolute eigenstate change of 4.10⁻⁸ (eV)². We also used an 8 × 8 × 8 Monkhorst-Pack grid and 512 k-points in the irreducible part of the Brillouin zone. In addition, the following parameters were considered:

• Densities: coarse grid: 12 × 12 × 12 grid

• Fine grid: 24 × 24 × 24 grid;

• Total loading: 0.000000;

• Density mixing: method (separate);

• Backend (pulay);

• Linear mixing parameter (0.05);

• Mixing with 5 old densities and damping of long-wave oscillations (50)

• Hamiltonian: XC and Coulomb potentials evaluated on a 24 × 24 × 24 grid using the PBE EC function; number of atoms (2);

• Number of atomic orbitals (8);

• Number of bands in the calculation (8);

• Number of valence electrons (8).

After initialization and convergence of the graphene structure, we obtained the following results (see Table 4):

Kinetic	+16.220681 eV
Potential	−15.352109 eV
External	+0.000000 eV
XC	−19.079831 eV
Entropy (-ST)	−0.000000 eV
Local	+0.041249 eV
=
Free energy	−18.170010 eV
Extrapolated	−18.170010 eV

Table 4.

Results obtained after initialization and convergence of the graphene structure.

Figure 4 plots the density of state obtained from this calculation. The calculated electronic density of state shows that diamond is a semi-metal with a non-zero gap E_g equal to 5.5 eV, which is in agreement with other calculations [19]. The electronic spectrum near the Fermi energy has a linear form. This behavior allows diamond to offer many applications in biosensors [20] and in devices for power electronics [21], because electrons have a high mobility in diamond at Dirac point of 4500 cm² V⁻¹ s⁻¹. Diamond is a wide-band-gap isolator. It is intrinsically insulating but can become semiconductor, then metallic and superconductor if it is suitably doped. This doping, which can be of type n or p, is obtained by the substitution of a carbon atom by a donor (n) or acceptor (p) atom, respectively.

Figure 4.

Density of state (DOS) of diamond.

3.4 PBE0 and PBE calculations for diamond

PBE is a GGA function introduced by Perdew, Burke, and Ernzerhof, in which all parameters other than those of its local spin density component are fundamental constants [22, 23]. Adamo and Barone [24] realized a new hybrid Hartree–Fock /DFT model derived from the PBE GGA, called PBE0, in which the HF exchange contribution is also fixed a priori [24]. The PBE0 approach is obtained by casting the PBE function in a hybrid scheme, in which the HF ratio is fixed a priori. The integer function can be expressed as follows:

$${\mathrm{E}}_{\mathrm{xc}}^{\mathrm{PBE}0}={\mathrm{E}}_{\mathrm{xc}}^{\mathrm{PBE}}+\frac{1}{4}\left({\mathrm{E}}_{\mathrm{x}}^{\mathrm{HF}}-{\mathrm{E}}_{\mathrm{x}}^{\mathrm{PBE}}\right)$$E1

Where${\mathrm{E}}_{\mathrm{xc}}^{\mathrm{PBE}}$ and ${\mathrm{E}}_{\mathrm{x}}^{\mathrm{PBE}}$ are the correlation and GGA exchange contributions, respectively, and ${\mathrm{E}}_{\mathrm{x}}^{\mathrm{HF}}$is the HF exchange [25]. It can be interesting to examine the two approaches of PBE and PBE0 to determine the lattice constant of diamond as well as its gap value.

3.4.1 Diamond lattice constant as a function of the number of k-points

In this section, we perform a non-self-consistent PBE0 calculation based on a self-consistent PBE to represent the variation of the diamond lattice constant as a function of k-point number.

We use a Monkhorst-Pack k-point grid [26], which is essentially a uniformly spaced grid in the Brillouin zone. Another less commonly used scheme is the Chadi-Cohen k-point grid [27, 28]. Monkhorst-Pack grids are specified as n1 × n2 × n3 grids, and the total number of k-points is n1-n2-n3. The computational cost is linear with respect to the total number of k-points, so a computation on a 4 × 4 × 4 grid will be about 8 times more expensive than that on a 2 × 2 × 2 grid. Therefore, we again seek to strike a balance between convergence and ease of computation [28]. Below we examine the k-point convergence of diamond carbon. Figure 5 presents the results of the PBE0 calculation, Figure 6 presents the results of the PBE calculation, and Figure 7 presents the results of the PBE0 and PBE calculation.

Figure 5.

Calculation of the diamond lattice constant in the case of PBE0.

Figure 6.

Calculation of the diamond lattice constant in the case of PBE.

Figure 7.

Calculation of the diamond lattice constant in the case of PBE0 and PBE.

To plot the three figures, we need a grid of at least 8 × 8 × 8 k-points to achieve a convergence level of at least 0.002 Å for PBE0 and 0.0005 Å for PBE. We find a lattice constant between 3.560 and 3.562 Å in the case of PBE0 and between 3.5815 and 3.5820 Å in the case of PBE. PBE and PBE0 do not converge to the same lattice constant, so they do not give the same results. PBE0 was found to perform better than the PBE function in reproducing the experimentally found geometric features (a = 3.567 Å) [29].

The convergence of the k-points is not always monotonic as in this example, and sometimes very dense grids (e.g., up to 20 × 20 × 20) are required for highly convergent properties such as the density of states in smaller unit cells. Oscillations in the lattice constant are typical, and it can be difficult to achieve high levels of convergence. Best practice is to use the same k-point sampling grid in the calculations when possible and dense (high number of k-points) otherwise. It is important to verify convergence in these cases. The larger the unit cells, the smaller the number of k-points required. For example, if a 1 × 1 × 1 fcc unit cell shows convergent results in a 12 × 12 × 12 k-point grid, then a 2 × 2 × 2 fcc unit cell would show the same level of convergence with a 6 × 6 × 6 k-point grid. In other words, doubling the unit cell vectors results in a halving of the number of k-points [28].

Sometimes k-points are described as k-points per reciprocal atom. For example, a 12 × 12 × 12 k-point grid for a primitive fcc unit cell corresponds to 1728 k-points per reciprocal atom. A 2 × 2 × 2 fcc unit cell has eight atoms (i.e., 0.125 reciprocal atoms), and thus, a 6 × 6 × 6 k-point grid has 216 k-points, (216/0.125 = 1728 k-points per reciprocal atom) [28].

In this k-point convergence study, we use an 8 × 8 × 8 k-point grid on a unit cell containing two atoms, resulting in 1024 k-points per reciprocal atom.

3.4.2 PBE and PBE0 diamond bands

In this section, we present the calculation of band gap width of diamond, using the PBE and PBE0 calculation. We use an 8 × 8 × 8 Monkhorst-Pack grid (resolution of 512 k-points and 29 k-points in the irreducible part of the Brillouin zone). Table 5 presents the resulting Γ- Γ and Γ-X gaps obtained for PBE and PBE0 in eV:

k -points	PBE(Γ- Γ)	PBE(Γ-X)	PBE0(Γ- Γ)	PBE0(Γ-X)
2	5.003	2.996	7.653	5.368
4	5.099	3.133	7.184	5.005
6	5.109	3.151	7.102	4.937
8	5.108	3.147	7.080	4.913

Table 5.

Resulting Γ- Γ and Γ-X gaps obtained for PBE and PBE0 in eV.

4. Conclusions

We have presented here a study of the electronic properties of diamond carbon, namely the band structure and the density of states. These calculations were based on the GPAW and ASE methods. The GPAW code is built on top of the ASE, which is a set of Python modules, and designed to facilitate the configuration, execution, and analysis of atomic/electronic calculations. We performed a comparative study of PBE and PBE0 approximations on lattice constant and diamond gap calculations. Regarding the lattice constant of diamond, we were able to reproduce the experimental value to within 0.2% for PBE0 and within 0.4% for PBE. We also calculated the diamond gap and found that our results are comparable with other experimental and theoretical results. Results show that the hybrid PBE0 was performed better than PBE in reproducing the geometrical and electronic characteristics of diamond. We also optimized diamond structure and found that the calculated band structure and density of states were in agreement with other methods. This chapter provided information for the theoretical and experimental communities working on electronic properties of diamond, which is considered one of the most promising materials for the integrated electronic and photonic, radio, optoelectronic, and quantum device industry.

from ase.build import bulk

from gpaw import GPAW, PW, FermiDirac

# Perform standard ground state calculation (with plane wave basis)

c = bulk('C', 'diamond', 3.57)

calc = GPAW(mode=PW(500),

  xc='PBE',

kpts=(8, 8, 8),

  random=True, # random guess (needed if many empty bands required)

occupations=FermiDirac(0.01),

txt='Cd_gs.txt')

c.calc = calc

c.get_potential_energy()

calc.write('Cd_gs.gpw')

# Restart from ground state and fix potential:

calc = GPAW('Cd_gs.gpw',

 nbands=16,

 fixdensity=True,

  symmetry='off',

 kpts={'path': 'GXWKL', 'npoints': 60},

 convergence={'bands': 8})

calc.get_potential_energy()

# Plots the dispersion relationships

bs = calc.band_structure()

bs.plot(filename='bandstructure-diamond.png', show=True, emin=-25, emax=25.0)

plt.savefig('bandstructure-diamond.png')

plt.show ()