First-Principles Study of the Structural, Electronic, Dynamic, and

Article pubs.acs.org/JPCC

First-Principles Study of the Structural, Electronic, Dynamic, and Mechanical Properties of HOPG and Diamond: Influence of Exchange−Correlation Functionals and Dispersion Interactions Christoph Lechner,† Baptiste Pannier,† Philippe Baranek,*,† Nancy C. Forero-Martinez,‡,§ and Holger Vach‡ †

EDF R&D, Department Materials and Mechanics of Components (MMC), EDF Lab Les Renardières, Avenue des Renardières, F-77818 Moret-sur-Loing Cedex, France ‡ LPICM, CNRS, Ecole Polytechnique, Université Paris-Saclay, F-91128 Palaiseau, France ABSTRACT: Various properties of two polymorphs of carbon, highly oriented pyrolytic graphite (HOPG) and diamond, were investigated at the ab initio level using different Hamiltonians with all-electron Gaussian-type functions (GTF) and projector augmented wave (PAW) basis sets. Their equilibrium lattice parameters, cohesive and interlayer interaction energies, band structures, vibrational frequencies, and elastic constants were evaluated. The calculations were performed at the Hartree−Fock, density functional theory (DFT), and hybrid (B3LYP and PBE0) levels. As regards DFT, the local density and generalized gradient (PBE and PBEsol) approximations were used. For GTF, the influence of the basis set superposition error (BSSE) was assessed. Since these approaches do not take dispersion interactions correctly into account, two different versions of Grimme’s dispersion correction, D2 and D3, were evaluated. The D2 and D3 corrections were reparameterized in order to reproduce the experimental structure of HOPG. For the properties depending on the description of the covalent bonds, such as the cohesive energy, C11 and C12, the dispersion corrections have a negligible influence. The best agreement is found for B3LYP-GTF-D3, PBE providing close results. For the properties depending on the description of the interaction between different layers, such as the interlayer interaction energies, C33, and C44, using the dispersion corrections improves the results. In general, D3 performs better than D2. PBE overall provides better results than the other functionals. As regards the basis sets, PAW describes the vibrational frequencies better than GTF due to the BSSE, whereas GTF provides a better description of the bandwidths and interband separations, elastic constants, and thermodynamical data. On the basis of the overall results, PBE-GTF-D3 appears to be a good compromise for an accurate description of the properties of graphite.

■

INTRODUCTION High interest in the properties of graphite has again arisen due to the upcoming dismantlement of several nuclear power plants of the type UNGG (Uranium Naturel Graphite Gaz), which used nuclear graphite as a neutron moderator and deflector. For instance, in France, the dismantlement will lead to 23 000 tons of irradiated graphite. Therefore, a detailed understanding of the behavior of the irradiated graphite is needed in order to find an optimal solution for permanent disposal. Clearly, this type of graphite is far more complex and heterogeneous in its structure than an ideal monocrystal of graphite.1,2 Nevertheless, it is necessary to investigate this most simple system in order to understand its properties and later on study the changes caused by the perturbation of this ideal state. Graphite is also a widely used material in other fields of industry, for example, in batteries3 and fuel cells,4,5 as lubricant,6,7 or for steel making. Graphite is the thermodynamically most stable allotrope of carbon under standard conditions.8 Its semimetallic character © XXXX American Chemical Society

and highly anisotropic structure lead to very interesting properties. It is an electrical and thermal conductor similar to a metal; however, it is rather chemically inert. In addition, it has a melting point of 3800 K and lubricity that resembles more closely that of a nonmetallic compound. Graphite naturally occurs in several forms, including the highly oriented pyrolytic graphite (HOPG), which will be the system of choice for this study. The interactions in graphite are highly anisotropic: The properties within the graphene planes are mainly driven by covalent forces; however, the interactions and distance between different planes are mainly driven by dispersive forces. The anisotropy of the interactions within graphite makes it difficult to be described accurately on a theoretical level. Received: October 23, 2015 Revised: January 1, 2016

A

DOI: 10.1021/acs.jpcc.5b10396 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C However, structural,9−17 electronic,18−23 dynamical,24−27 mechanical, 1 3 , 1 7 , 2 3 , 2 8 , 2 9 and thermodynamical properties13,14,16,23,27,30−33 have been studied extensively. Nowadays, it is generally accepted that standard density functionals are not able to describe dispersion interactions correctly. While the local density approximation (LDA) usually overbinds dispersively bound systems, the generalized gradient approximation (GGA) and hybrid schemes often are not even able to describe the dispersion interactions in a qualitatively correct manner.34,35 Furthermore, the Hartree−Fock method, which completely neglects dynamical correlation, also fails at describing this type of interaction. However, a good description of these interactions is needed in order to describe the graphite system accurately as a whole. In the last 10−15 years several different schemes have been proposed to make up for this systematic lack of dispersion interaction. Among them, there are the nonlocal van der Waals density functionals (vdw-DF), 36−39 where a long-range contribution of the correlation energy, which is determined locally from the electron−gas energy evaluated with gradient corrected LDA, is included in the exchange−correlation functional. There exist also several semiempirical schemes, such as those proposed by Grimme40−42 (the two most recent being commonly named “D2” and “D3” corrections) or by Tkatchenko and Scheffler,43,44 which have a more or less simple functional dependence of the dispersion interaction on the nuclei coordinates. While Grimme’s dispersion correction D2 is applied very often to study solid state systems,16,27,45−48 his most recent version D3 is not available in all periodic codes and thus far less used.49−53 However, several modifications were introduced in D3, such as polarizability constants, which depend on the coordination number of an atom, as well as three-body interactions. These new features may result in a more accurate description of the dispersion interactions acting in a solid system. The vast majority of modeling in the field of solid state physics is carried out by using plane wave basis sets. A much less used possibility is to express the Bloch functions as a linear combination of atom-centered Gaussian-type functions (GTF). This approach originates from the field of quantum chemistry. While the orthogonality of plane waves guarantees that the properties are not biased by the basis set superposition error (BSSE), the locality of GTF assures an easy accessibility of a population analysis or projected density of states (PDOS). Additionally, the use of hybrid functionals is much easier with a GTF basis, since HF calculations are well established in quantum chemistry. However, due to the extensive use of plane wave basis sets, there is little knowledge about the impact of the basis set on the different properties. A very general overview and comparison of these two methods is given in ref 54. The literature of theoretically studied graphite is rich but scattered. However, to our knowledge a detailed investigation of the performance of different Hamiltonians, density functional theory (DFT) and HF, does not exist. This work aims at exploring the performance of different methods with respect to the description of the above-mentioned properties in HOPG and diamond. Therefore, it highlights the differences between results obtained by using different Hamiltonians and basis sets or by applying the D2 and D3 correction schemes, which account for the missing dispersion interactions. Moreover, it presents a comprehensive discussion to elucidate which method works best to investigate a given property.

The article is organized as follows: In the next section, the computational details are described. Structural data and thermodynamical, electronic, dynamic, and mechanical properties are then reported. The obtained results are discussed with respect to other theoretical work and experimental data. Finally, we conclude which method is able to describe the properties studied to a satisfying extent and which version of Grimme’s dispersion correction is more suitable for the studied system.

■

METHODOLOGICAL ASPECTS The calculations were performed with the periodic CRYSTAL0955,56 and QUANTUM-ESPRESSO57 codes. Within CRYSTAL, both the Hartree−Fock and the Kohn−Sham equation as well as hybrid schemes, such as B3LYP58−60 or PBE0,61 where the exchange operator is a linear combination of the HF and DFT ones, can be solved self-consistently. Pseudopotential and all-electron Gaussian-type basis sets can be used. The QUANTUM-ESPRESSO code, on the other hand, uses the density functional techniques as well as the hybrid schemes58−61 with different plane wave pseudopotentials.57 In this work, several exchange-correlation functionals based on the local density approximation and generalized gradient approximation were used. They will be indicated as follows: LDA for Dirac−Slater62 exchange plus Perdew−Zunger63 correlation potential; PBE for Perdew−Burke−Ernzerhof64 exchange-correlation functional, and PBEsol65 corresponding to the revised PBE improving the description of the equilibrium properties of solids. Two hybrid HF/KS functionals were also considered, namely, B3LYP58−60 and PBE0.61 All of the studied properties were determined with the primitive cells of HOPG and diamond. For the Gaussian basis sets, carbon was described by an allelectron basis set.66,67 The basis set is a contraction of 6s211sp-1d* GTFs (a contraction of 6, 2, 1, 1, and 1 Gaussiantype functions for the 1s, 2sp, 3sp, 4sp, and d shells, respectively). The outer exponents of C were optimized using an energy criterium and are reported in Table 1. The basis set optimization was carried out using the LoptCG68 script, which performs numerical gradient optimizations based on the conjugate gradient method.69 Table 1. Exponents and Coefficients of the Contracted Gaussian Basis Set Adopted in the Present Study for Ca coefficient shell

exponent

s(d)

p

sp

11.705961 2.877171 0.792605 0.234684 0.927994

−0.160393 −0.645825 1 1 1

0.296329 1.383579 1 1

d a

Only the most diffuse 211sp-1d* GTF are given (see refs 66 and 67 for a complete set of data).

Within QUANTUM-ESPRESSO, projector augmented waves (PAW) with a scalar relativistic pseudopotential were used as the basis set. In this pseudopotential, the carbon 2s and 2p electrons were treated as valence electrons. In order to take into account the contribution of dispersion interactions to the interlayer interaction of HOPG more reliably, two types of Grimme’s correction were investigated: B


Article

The Journal of Physical Chemistry C the semiempirical London dispersion-type correction40,41 (the D2 correction implemented in CRYSTAL and QUANTUMESPRESSO) and the most recent dispersion correction D3,42 which we implemented in CRYSTAL09 for this study. As regards the D2 correction, the original dispersion coefficient Cij6 and atomic van der Waals radius RvdW as proposed by Grimme40,41 were used, since the net atomic charge of C in graphite and diamond is close to that of neutral carbon according to a Mulliken population analysis. The values are 1.75 J nm6 mol−1 and 1.452 Å, respectively. The damping function steepness d was fixed at 20 Å. The cutoff distance to truncate the direct lattice summation was 30 Å. Only the scaling factor s6 was optimized in order to obtain a c parameter with an error of less than 0.5% with respect to experiment. For each level of approximation, the used values for the s6 parameter are given in Table 2. In the following, the obtained corrected results will be referred to as “-D2”.

these points leads to severe convergence issues. The meshes are listed in Table 3. Table 3. Monkhorst−Pack k-Point Meshes Used for CRYSTAL Calculations calculation optimization frequencies elastic constants cohesive energy interlayer interaction energy

basis

HF

LDA

PBE

PBEsol

B3LYP

PBE0

0.585 1.340a

−0.470 −2.740a −0.025

0.250 0.410 0.5775

−0.087 −0.460 0.2325

0.350 0.740

0.167 0.255

18 16 14 18 16

× × × × ×

18 16 14 18 16

× × × × ×

18 16 14 18 16

For QUANTUM-ESPRESSO, plane wave cutoffs of 900 and 250 au were used for HOPG and diamond, respectively. Methfessel−Paxton first-order spreading of 0.01 au was applied. For QUANTUM-ESPRESSO, all calculations were performed with 18 × 18 × 18 and 16 × 16 × 16 Monkhorst−Pack k-point meshes71 for HOPG and diamond, respectively. For both codes, the convergence criteria on total energies (and for the elastic constants and frequencies determination) were 10−9 (10−12) au. Atomic displacements and forces thresholds were 1.8 × 10−3 and 4.5 × 10−4 au, respectively. With these computational conditions, the obtained data can be considered as fully converged. Cohesive energies, Ecoh, were evaluated with respect to the isolated atoms. For both codes, the equilibrium energies of bulk graphite and diamond were corrected by the zero point energy. The cohesive energy follows this definition

Table 2. s6 (D2) and s8 (D3) Scaling Factors Used for Each Hamiltonian and Basis Set GTF-D2 GTF-D3 PAW-D2

mesh

a

To be consistent with the parametrization of the D3 correction at the PBE, PBEsol, B3LYP, and PBE0 levels, sr,870 has been set to 1.0 for HF and LDA.

Regarding the D3 correction, the influence of the Becke− Johnson and zero damping was evaluated for the PBE, PBEsol, PBE0, and B3LYP functionals with the proposed parameters by Grimme:70 For all Hamiltonians investigated, both dampings underestimate the interlayer distance by ∼4 ± 2% with respect to experiment. In order to correct this behavior, the scaling parameters of the dispersion correction functions were reparameterized. However, we noticed that an optimization of the Becke−Johnson damping parameters does not allow one to obtain an efficient correction of this discrepancy; for instance, at the PBE level, a decrease of s6 from 1.0 to 0.75 and of s8 from 0.7875 to 0.25 increases the c parameter only from 6.4 to 6.6 Å. Therefore, only zero damping was used. With s6 fixed at 1.0, the s8 parameter was optimized to obtain a c parameter with an error of less than 0.5% with respect to experiment. This approach offers the same flexibility as the reparametrization of the D2 correction while preserving the correct asymptotics of the dispersion interactions. The cutoff values for dispersion interaction and coordination number were 50 and 20 Å, respectively. The three-body terms were not included. For each level of approximation used within CRYSTAL, the used values for the s8 parameter are given in Table 2. In the following, the obtained corrected results will be referred to as “-D3”. For evaluation of the Coulomb and exchange series within CRYSTAL, the truncation thresholds for the bielectronic integrals, as defined in the CRYSTAL manual,56 were set to 10−8, 10−8, 10−8, 10−8, and 10−16 au. Smearing of the Fermi surface was not used. Due to CRYSTAL09’s limit with respect to the number of irreducible k points, different Monkhorst− Pack meshes71 were used in order to balance the increase in k points with the various reductions of symmetry, as needed for determination of the elastic constants and frequencies. For evaluation of the interlayer interaction energy, the mesh was reduced due to a degeneracy at high-symmetry points, sampling

Ecoh =

N0EC − E Bulk N0

(1)

where EC, EBulk, and N0 represent the energy of an isolated carbon, the energy of the primitive cell of the bulk graphite (or diamond), and the number of carbon atoms per primitive cell, respectively. The EC calculations have been performed in the spin-polarized framework and carbon in its triplet state. Contrary to plane wave basis sets, Gaussian-type basis sets are nonorthogonal and therefore suffer from the BSSE. Thus, for CRYSTAL, Ecoh was corrected for the BSSE according to the counterpoise method of Boys and Bernardi.72 To evaluate it a three-layer (001) HOPG surface, periodic along the [100] and [010] directions, was defined; then the energy of an isolated carbon was determined by placing it at the center of a 4 × 4 supercell of the slab with all of its neighbors treated as ghost atoms. For diamond, the BSSE was evaluated from the energy of an isolated carbon at the center of a cluster of 80 ghost atoms “cut” in the bulk material. For QUANTUM-ESPRESSO, the energy of the isolated carbon was obtained by using a cubic supercell with a lattice parameter of 15 Å. The interlayer interaction energy, EIL, follows this definition Bulk E IL = E I − E IL

(2)

EBulk IL

where EI and are the energies of an isolated graphene plane and of a graphene plane in graphite, respectively. EBulk IL is the energy difference between a seven- and a six-layer (001) surface. For CRYSTAL, EIL was also corrected for the BSSE following the scheme of Boys and Bernardi:72 To evaluate it, the energy of an isolated graphene plane with a single layer of ghost atoms above and below was determined. C


Article

The Journal of Physical Chemistry C

Figure 1. Crystal structure and vibrational modes of HOPG. The space group is P63/mmc with a, b, and c the lattice axes of the primitive hexagonal cell; green and red atoms are the two irreducible carbons with fractional coordinates (0, 0, 1/4) and (1/3, 2/3, 1/4), respectively. Given frequencies refer to the experimental data when available (see Table 9). The silent mode at 127 cm−1 can be measured with neutron scattering. Arrows describe the atomic displacements characterizing each vibration.

carbons, (0, 0, 1/4) and (1/3, 2/3, 1/4), respectively. Within the graphitic planes, the carbon atoms are covalently bound due to the sp2 hybridization, whereas due to the delocalization of the partially occupied 2pz orbitals perpendicular to the graphene layers, the weak interlayer binding is supposed to arise from dispersion or van der Waals interactions. It is convenient to use dC−C, the distance of the C−C bond within the graphitic planes (or in diamond), and dg, the interlayer spacing, as geometrical parameters to define the structure. Table 4 shows the structural parameters computed in this study with different Hamiltonians, basis sets, and types of dispersion correction for graphite. Figures 2 and 3 show the relative errors obtained for dC−C and dg without dispersion correction with respect to experiment.79 In order to make the evaluation of the performance of the various Hamiltonians and basis sets easier, the determined lattice parameters, elastic

For QUANTUM-ESPRESSO, which uses a 3D model of the surfaces, the various slabs were separated by a void of 10 Å. With QUANTUM-ESPRESSO, the ”Elastic” script73 was used for the determination of the elastic constants. For the visualization the XCRYSDEN,74 MOLDRAW,75−77 and POVRay78 software were used.

■

RESULTS AND DISCUSSION Structure Optimization. Graphite, in its most ordered phase called highly oriented pyrolytic graphite, is composed by a regular stacking of bidimensional hexagonal layers (graphene planes). The most stable stacking of carbon layers corresponds to the ABAB stacking shown in Figure 1. The corresponding crystal structure has the space group symmetry P63/mmc; the unit cell geometry is fully defined by the a and c lattice parameters and the fractional coordinates of two nonequivalent D


Article

The Journal of Physical Chemistry C Table 4. Equilibrium Geometry of HOPG Calculated with Different Hamiltonians and Basis Setsa Hamiltonian

basis

HF HF-D2 HF-D3 LDA LDA-D2 LDA-D3 LDA LDA-D2 PBE PBE-D2 PBE-D3 PBE PBE-D2 PBEsol PBEsol-D2 PBEsol-D3 PBEsol PBEsol-D2 PBE0 PBE0-D2 PBE0-D3 B3LYP B3LYP-D2 B3LYP-D3

GTF GTF GTF GTF GTF GTF PAW PAW GTF GTF GTF PAW PAW GTF GTF GTF PAW PAW GTF GTF GTF GTF GTF GTF

LDAb PW91b PBEb LDAc PW91c vdW- DF1(revPBE)d LDAe LDAe LDAf PBEf PBEg PBE-D2g vdW-DF2(revPBE)h vdW-DF2h HSE06/PBEsolI PBE-D2I PW91+vdWj exp.

NC-PP/PW NC-PP/PW NC-PP/PW FLAPW FLAPW PW US-PP/PW PAW US-PP/PW US-PP/PW PAW PAW PAW PAW PAW PAW PW

a 2.440 2.438 2.440 2.445 2.446 2.445 2.445 2.445 2.464 2.463 2.463 2.466 2.463 2.457 2.457 2.457 2.459 2.458 2.447 2.447 2.447 2.454 2.453 2.454 2.442 2.465 2.466 2.447 2.470 2.470 2.443 2.448 2.440 2.461 2.470 2.460 2.472 2.470 2.447 2.464 2.455 2.459k 2.456l 2.461m 2.464n

c 7.327 6.673 6.669 6.140 6.670 6.676 6.636 6.672 7.049 6.672 6.671 8.620 6.676 6.515 6.675 6.674 7.245 6.672 6.958 6.671 6.674 7.087 6.670 6.673 other works 6.622 6.742 6.745 6.622 10.251 7.520 6.575 6.582 6.686 8.490 8.840 6.450 7.195 6.891 7.225 6.436 6.690 6.672k 6.696l 6.706m 6.711n

c/a

dC−C

dg

V0

3.003 2.737 2.733 2.511 2.727 2.731 2.714 2.728 2.861 2.709 2.708 3.496 2.711 2.652 2.717 2.716 2.947 2.715 2.843 2.726 2.727 2.887 2.719 2.719

1.409 1.407 1.409 1.412 1.412 1.411 1.412 1.412 1.422 1.422 1.422 1.424 1.422 1.419 1.419 1.419 1.420 1.419 1.413 1.413 1.413 1.417 1.416 1.417

3.663 3.337 3.334 3.070 3.335 3.338 3.318 3.336 3.525 3.336 3.336 4.310 3.338 3.258 3.338 3.337 3.623 3.336 3.479 3.335 3.337 3.544 3.335 3.337

37.771 34.343 34.387 31.781 34.572 34.550 34.362 34.551 37.056 35.048 35.060 45.398 35.073 34.061 34.905 34.890 37.934 34.903 36.094 34.585 34.618 36.975 34.762 34.814

2.702 2.735 2.735 2.706 4.150 3.045 2.691 2.689 2.740 3.450 3.579 2.622 2.911 2.790 2.953 2.612 2.725 2.713 2.726 2.725 2.724

1.413 1.423 1.424 1.413 1.426 1.426 1.411 1.413 1.409 1.421 1.426 1.420 1.427 1.426 1.413 1.422 1.228 1.420 1.418 1.421 1.423

3.311 3.371 3.373 3.311 5.126 3.760 3.288 3.291 3.343 4.245 4.420 3.225 3.598 3.446 3.612 3.218 3.345 3.336 3.348 3.353 3.356

34.339 35.478 35.522 34.339 54.162 39.732 33.984 34.159 34.473 44.531 46.706 33.803 38.077 36.409 37.466 33.840 34.919 34.938 34.979 35.174 35.286

a

Lattice parameters (a and c), C−C bond length in the graphene plane (dC−C), and interlayer distance (dg) are in Angstroms; equilibrium volume of the unit cell (V0) is in Å3. Other theoretical results and experimental data are given for comparison. bReference 9. Hamann’s norm-conserving pseudopotential.10 cReference 11. dReference 12. eReference 14. fReference 15. gReference 16; in this work s6 = 0.75. hReference 17. IReference 27; in this work s6 = 0.75. jReferences 13 and 89. London-type dispersion with parametrization of Ortmann et al. kReference 79; experiment at 4 K. l Reference 98. mReferences 99−101; experiment at 300 K. nReference 102.

The computed dC−C are in good agreement with the experimental data, whatever the basis set and Hamiltonian; the average relative error is ∼−0.3% with a mean total deviation of ∼0.5%. HF and LDA underestimate a and dC−C by 1% but not for the same reasons: LDA systematically underestimates the lattice parameters of crystals whatever their nature, but at the HF level, due to the lack of correlation, the electronic charge is more concentrated around the atoms or around the

properties, indirect band gap, Raman frequency, and cohesive energy of diamond (space group Fd3¯m) are given in Table 5 for comparison. Both phases have been deeply investigated in the literature. For diamond, our results are generally in line with previous calculations. For graphite, it is also the case for the a parameter; however, the results with respect to the c parameter show major differences depending on the type of basis sets and the used computational conditions, such as the energy cutoff. E


Article


(i) The inf luence of the basis sets. As shown by the results in Table 4, whatever the Hamiltonian, the obtained data for dg are always lower for GTF with respect to those obtained with PAW. This is due to the BSSE which acts as an “additional artificial attractive force”. For instance, at the PBE level dg increases from 3.5 to 3.9 Å, when the variation of the potential energy with respect to the cell parameter c is corrected for the BSSE. When the D2 dispersion correction with the PAW s6 parameter is added to this corrected potential energy curve, the obtained dg is 3.33 Å, a good agreement with the one determined with QUANTUM-ESPRESSO. (ii) The inf luence of the Hamiltonians. Except for LDA, the interlayer distance is largely overestimated (∼+10%). The LDA data show an underestimation of dg by about −10%; this is not surprising as it is well known that LDA overbinds systems held together by dispersive forces.34,35 It is interesting to note that LDA with PAW reproduces the experimental structure rather well, which is a striking example of error cancellation. (iii) The inf luence of the dispersion correction. Grimme’s D2 and D3 dispersion corrections represent an attractive contribution to the interlayer interactions. However, the used parametrization has to be adapted to the material of interest.16,45,46 The original parametrization by Grimme, optimized for molecular systems, is known to work rather unsatisfyingly “out of the box” when using periodic boundary conditions. This is shown by the results of Bučko et al.,16 where using Grimme’s parameter leads to an underestimation of dg with respect to experiment. It is also useful to point out that for the functionals which underestimate dg without dispersion correction, such as LDA and PBEsol-GTF, the used parametrizations have to be considered from a purely empirical point of view, neglecting their physical meaning: for instance, in the case of LDA it allows one to correct the classical overbinding of the layers at this level; in the case of PBEsol-GTF it allows one to correct the BSSE (see Table 4). Finally, Figure 4 summarizes all the points mentioned above: it shows the influence of the choices of basis set, Hamiltonian, and dispersion corrections upon the obtained corrections of the lattice parameter c. It clearly illustrates the need to adapt the parametrization of Grimme’s dispersion corrections to the studied material and to the used code: The slope of GTFs is higher than that of PAW due to the BSSE; for higher initial errors, a larger s6/s8 parameter is needed to obtain the experimental structure. It is interesting to note that overall the obtained corrections vary linearly with the scaling parameters, which is probably due to the pair interaction approximation of the dispersion correction used. If this point could be verified for other molecular or lamellar crystals, it might be a simple way to parametrize these types of corrections. Cohesive and Interlayer Interaction Energies. The cohesive energies, evaluated with respect to the material in its standard state, and the interlayer interaction energies are given in Table 6 (diamond’s Ecoh is also given in Table 5). The experimental values of Ecoh have been estimated from the standard enthalpies of solids and gas, extrapolated at T = 0 K.8 Several experimental values exist for the interlayer interaction energy; however, all of those energies are not measured directly but are rather a combination of experimental data and theoretical calculations. Here, the two values of Benedict et

Figure 2. Relative errors of the distance dC−C of HOPG with respect to experimental data at 4 K79 for different Hamiltonians and basis sets without dispersion correction.

Figure 3. Relative errors of the interlayer spacing dg of HOPG with respect to experimental data at 4 K79 for different Hamiltonians and basis sets without dispersion correction.

bond (or, in other words, the system is more ionic or covalent (see, for instance, refs 80 and 81). Among the GGA functionals, PBE overestimates a and dC−C while PBEsol corrects this tendency to reach a good agreement with the experimental data. As for other materials,82,83 the obtained error with the PBEsol functional is less than ∼0.1%. Among the hybrid functionals, PBE0 provides results with the same accuracy as HF and LDA; however, for ionic materials81,84 this functional is better adapted. The best description of dC−C is obtained with the B3LYP functional: the deviation with respect to experimental data is ∼0.2% for graphite and diamond. In the following, the obtained theoretical results on the interlayer distances will be described. Since dg is driven by the weak interlayer binding, the obtained theoretical data characterize the difficulties of taking the dispersion interactions into account. For instance, without dispersion correction, this parameter is underestimated or overestimated by 0.3−1 Å at the LDA-GTF and PBE-PAW levels, respectively (the other Hamiltonians giving intermediate results). This behavior is linked to the dispersive character of the interlayer interaction which is not correctly accounted for at the various levels of approximation (as it was shown for other materials34,35,80): as it will be discussed in the next section, the interlayer interaction energy EIL is small (less than 50 meV/ atom, see Table 6) and varies smoothly with the interlayer distances.34,35,80 The obtained data for dg deserve to be discussed in detail. The determined interlayer distance depends on three factors, namely, the choices of the basis sets, Hamiltonians, and dispersion correction. F


Article


Table 5. Equilibrium Geometry, Elastic Properties, Band Gap, and Raman Frequency of Diamond Calculated with Different Hamiltonians and Basis Setsa Hamiltonian

basis

HF LDA LDA PBE PBE PBEsol PBEsol PBE0 B3LYP

GTF GTF PAW GTF PAW GTF PAW GTF GTF

HFb LDAb HFc LDAd B3LYPe LDAf LDAg LDAh PW91h LDAI PW91I PBEI PBEj PBEk HSE06/PBEsoll exp.

LMTO LMTO GTF GTF GTF NC-PP/PW US-PP/PW FLAPW FLAPW NC-PP/PW NC-PP/PW NC-PP/PW US-PP/PW US-PP/PW PAW

a 3.550 3.535 3.535 3.570 3.572 3.554 3.555 3.545 3.568

B 509 464 467 444 422 460 452 481 451

C11 1287 1118 1115 1070 1057 1088 1071 1156 1118 other works

C12

C44

Eg

νRaman

Ecoh

EBSSE coh

120 153 161 131 127 145 141 144 118

686 604 591 572 561 587 573 623 590

12.35 4.18 4.16 4.14 4.13 4.05 4.03 6.01 5.91

1535 1333 1320 1303 1289 1319 1305 1390 1352

5.08 8.86 8.87 7.58 7.59 8.13 8.15 7.47 6.93

5.63 9.54 8.24 8.78 8.06 7.57

12.10 4.01 3.580 3.570

471

3.530 3.539 3.535 3.579 3.530 3.570 3.570 3.568 3.568 3.538 3.567m 3.556o

473

1484 4.01 5.80

465 430 455 438 439 432 444 442n 452r

1324

1060 1098

125 118

562 589

1079n

124n

578n

1289 4.75 5.48q 6.00s

1356 1332p

7.34t

Lattice parameter (a) in Angstroms; bulk modulus (B) and elastic constants (Cij) in GPa; indirect Γ−Δ1 band gap (Eg) in eV; Raman frequency (νRaman) in cm−1; cohesive energy including the zero point energy (Ecoh) in eV/atom (the energy not corrected for BSSE (EBSSE coh ) is also given). Other theoretical results and experimental data are given for comparison. bReference 103 with the experimental lattice parameter at 300 K. cReference 104. d Reference 105. eReference 106 with the experimental lattice parameter at 300 K. fReference 107. gReference 108. hReference 11. IReference 9. Hamann’s norm-conserving pseudopotential.10. jReference 15. kReference 109. lReference 27. mReference 99; experiment at 300 K. nReference 110 o Reference 8; experiment at 4 K. pReference 111. qReferences 112−114. rReference 115. sReference 116. tDeduced from the standard enthalpies of solids and gas extrapolated at T = 0 K of ref 8. a

al.85 and Zacharia et al.86 have been taken as the experimental lower and upper bound of EIL, respectively. For Ecoh, the contribution of covalent interactions is much more important. The order of magnitude of Ecoh is eV/atom, whereas for EIL it is meV/atom. The obtained results for graphite and diamond share the same trend: As expected, HF understimates Ecoh by about 30%, while LDA, on the other hand, largely overbinds the two systems by about 20%. The other functionals provide data between HF and LDA and closer to experiment: The best agreement is found for B3LYP-GTFD3. In Table 6, for GTF, Ecoh not corrected for the BSSE are also supplied. As it is shown, the average contribution of the BSSE to Ecoh is about 6%. As expected, the dispersion corrections have a negligible effect on Ecoh: their average contribution is less than 1.5%. Globally, for Ecoh, PAW and GTF results are in perfect agreement. One of the direct applications of the knowledge of Ecoh is the determination of the relative stability of diamond with respect to graphite, ΔE (EDiamond − EGraphite). Experimentally, ΔE is really small, ∼0.02 eV/atom at T = 0 K, as deduced from ref 8. However, the theoretical evaluation of ΔE is still under discussion, mostly fed by the difficulties to evaluate the influence of the dispersion interactions, zero point energy

(ZPE), and computational parameters on the properties of these two polymorphs (see, for instance, refs 11, 13, 20, 27, 87−89 and references therein). These previous works show that ΔE ranges from −0.03 to 0.31 eV/atom, depending on the Hamiltonian, basis set, dispersion, and ZPE corrections. Some authors, such as Ortmann et al.13 or Teobaldi et al.,89 studied the effect of accounting for the dispersion interaction on the diamond description at the GGA level: Adding the dispersion correction worsens the cohesive energies of graphite and diamond from 7.72 to 7.86 eV/atom and from 7.55 to 7.80 eV/ atom, respectively; the induced improvement of ΔE from 0.17 to 0.06 eV/atom is simply a cancellation of errors. Grochala27 studied the various thermodynamic properties of carbon with the hybrid HSE06 functional. As stated by Paier et al. in ref 90 this functional should give the same results as PBE0. However, Grochala obtained ΔE = −0.001 eV/atom, which is contradictory to our findings for PBE0. This seems to be due to the use of PBEsol exchange−correlation in the hybrid HSE06 functional (PBEsol underestimates ΔE with respect to PBE) and the absence of any dispersion correction for graphite. From these various works, the only clear insight is that LDA systematically determines diamond as the ground state of carbon. Our results on ΔE are summarized in Table 7. Further clear trends appear. First, except for LDA-PAW, graphite is G


Article

The Journal of Physical Chemistry C Table 6. Cohesive Energy Including the Zero Point Energy (Ecoh, in eV/atom) and Interlayer Interaction Energy (EIL, in meV/atom) of HOPG Obtained with Different Hamiltonians and Basis Setsa Hamiltonian

basis



LDAb LDAb LDAc PBEd RPAe PBEf PBE-D2f vdW-DF2g PBE-TSh PBE-TS+SCSh PW91+vdWI exp.

PAW US-PP/PW LCGTO-FF NC-PP/PW PAW PAW PAW US-PP/PW PAW PAW PAW

Ecoh

EBSSE coh

5.28 5.65 5.35 5.74 5.39 5.78 8.91 9.57 8.86 9.49 8.90 9.53 8.86 8.86 7.82 8.28 7.84 8.32 7.88 8.35 7.72 7.80 8.23 8.71 8.23 8.69 8.26 8.73 8.17 8.20 7.60 8.01 7.62 8.04 7.65 8.07 7.25 7.70 7.29 7.75 7.32 7.78 other works 8.89 8.79 8.87 7.63

7.86 7.36j

EIL −6 31 62 19 −26 4 25 23 −2 13 43 1 39 −5 −11 14 4 23 −3 7 33 −7 16 42

EBSSE IL 3 50 81 42 −10 19

Figure 4. Correction (Δc) of the determined HOPG lattice c parameter as a function of the scaling parameters of the dispersion corrections (in order to obtain c with an error less than 0.5% with respect to experiment at 4 K79). s6 and s8 are given in Table 2 and correspond to the scaling parameters of the D2 and D3 dispersion corrections, respectively. The results are those obtained with the GTFD2, GTF-D3, and PAW-D2 methods using the HF (▶), LDA (●), PBE (◆), PBEsol (■), PBE0 (▲), and B3LYP (▼) Hamiltonians.

9 30 59

13 5 30

expected, the dispersion corrections have no impact on ΔE and show the same trends as for Ecoh, as discussed above. Except for LDA-PAW and PBEsol-PAW, ΔE is systematically overestimated. Considering the interlayer interaction energy, as shown in Table 6, the influences of the BSSE and of the dispersion corrections are decisive for the calculation of EIL. For the GTF results, the obtained data clearly show that the BSSE is far from being negligible: Depending on the Hamiltonian, its contribution ranges from 50% to 300%. For instance, the results without dispersion correction show that at the HF, PBE, PBEsol, or hybrid levels there is simply no cohesion between graphitic layers without considering the BSSE. These results are coherent with the PAW ones (free of BSSE) which overestimate dg and yield small values of EIL. As expected, LDA overbinds the planes. For EIL, the choice of the dispersion correction is equally important. D2 and D3 tend to not only correct but even overcorrect this energy. The most extreme variation is obtained for HF, where it changes from −6 meV/atom, an underestimation of at least −117%, to 62 meV/atom, an overestimation of at least 19% (when considering the lower and upper bounds as references). While for Ecoh no significant advantage is gained by using D3, the interlayer interaction energies show the limits of D2 and the improvement offered by D3: They are in general underestimated by D2 due to the s6 parameter being largely decreased with respect to Grimme’s original parameters; additionally, the asymptotics of the dispersion interaction are not correct for all Hamiltonians (since s6 ≠ 1). By using D3 both of these problems are solved. Keeping s6 = 1 assures the correct asymptotics, tuning the s8 parameter alters the force, but the total dispersion correction remains negative. Finally, considering the D3 results, HF overestimates EIL by about 45%, while LDA on the other hand underestimates it by about 90%. The other functionals provide data between HF and LDA and closer to experiment. PBE-GTF-D3 and hybrid-GTFD3 have the best agreement with the experimental data. Band Structure. In this section, since graphite is a semimetallic material, a set of interband separations for which experimental data are available, has been used to evaluate the performance of the different methods. Tables 5 and 8 sum up the results for the band gap Eg of diamond and the bandwidths

8 23 49 5 34 59

48 1 55 53 82 55 84 35 ± 15k, 52 ± 5l

a For the GTF basis sets, the energies not corrected for BSSE (EBSSE coh , EBSSE IL ) are also given. Other theoretical results and experimental data are given for comparison. bReference 14. cReference 23. dReference 30. eReference 31. fReference 16. gReference 32. hReference 33. I References 13 and 89; London-type dispersion with the Ortmann et al. parametrization. jDeduced from the standard enthalpies of solids and gas extrapolated at T = 0 K of ref 8. kReference 85. lReference 86.

systematically found as the carbon’s ground state. Omitting the zero point energies has a negligible effect on ΔE, since they are roughly the same for both polymorphs. However, except for LDA, the BSSE is the most important source of error for ΔE; not correcting for the BSSE changes it on average by about −0.2 eV/atom. This is due to the fact that the BSSE is ∼0.2 eV/atom higher for diamond than for graphite. The trend of LDA is just due to a cancellation of errors: Indeed, at this level the cohesive energies not corrected for BSSE, EBSSE coh , are overestimated with respect to experiment and are 9.57 and 9.54 eV/atom for graphite and diamond, respectively. As well for PBEsol and PBE0, not taking into account the BSSE can cause a relative stability reversal between these two phases. As H


Article

The Journal of Physical Chemistry C Table 7. Calculated Energy Difference (ΔE = EDiamond − EGraphite) between Diamond and Graphite (in eV/atom) Obtained with Different Hamiltonians and Basis Sets; ΔE Includes the Zero Point Energya basis

ΔE



0.20 0.27 0.31 0.05 0.00 0.04 −0.01 −0.01 0.24 0.26 0.30 0.13 0.21 0.10 0.10 0.13 0.02 0.05 0.13 0.15 0.18 0.32 0.36 0.39 other works

LDAb LDAc PW91c PW91d PW91+vdWd LDAe PBEsole PBE-D2e HSE06/PBEsole exp.

PW LAPW LAPW PW PW PAW PAW PAW PAW

Hamiltonian

−0.00(1) 0.03 0.04 −0.00(3) 0.02f

ΔEnoZPE

ΔEBSSE noZPE

0.20 0.27 0.31 0.05 0.00 0.05 −0.01 −0.01 0.24 0.26 0.30 0.13 0.21 0.10 0.10 0.13 0.03 0.06 0.14 0.15 0.06 0.33 0.36 0.40

0.02 0.11 0.15 0.03 −0.05 0.00

best agreement is found for the LDA, PBE, and PBEsol functionals, which give similar results. Since Grimme’s dispersion corrections are applied a posteriori, their influence on the bandwidths are minor, as they propagate only through the geometry optimization. Figure 5 shows the band structure obtained with PBE-GTF and PBE-PAW as well as the density of states obtained with PBE-GTF. They are typical of a semimetallic material: no or a small gap at specific k points (here at the K and H points) and no density of states at the Fermi level. The projected DOS show that near the Fermi level the states are mainly π (pz) states in the valence band (between 0 and −5 eV) and π* (pz) in the conduction band (between 0 and 7.5 eV); as expected, the σ (spxpy) states are lower in energy (between −2.5 and −17 eV) and the σ* (spxpy) states higher up in the conduction band (over 7.5 eV). The obtained band structures show that GTF and PAW basis sets give different descriptions of certain valence and conduction bands. This is illustrated in Figure 6, which shows the relative errors with respect to the bandwidths and interband separations. • Considering the valence bands, the same results are obtained for the Γ1+−Γ5+ and Γ1+−Γ3+ interband separations for each Hamiltonian, whatever the basis set. Both dispersion corrections have no influence on the calculated data (the induced variations are less than 1%). In spite of the variance of the experimental data of the Γ1+−Γ3+ transition, the general trends given previously are verified: HF and hybrid functionals overestimate it; the best agreement with experiment is found for the LDA, PBE, and PBEsol functionals. However, the obtained errors are much more important for the splitting of the two π bands, Γ3+−Γ2− (see Figure 6). First, the error discrepancies between the GTF and the PAW results are due to numerical artifacts due to the small value of this interband transition energy. Nevertheless, a general trend can be recognized: Except for PBE-PAW, this splitting is systematically overestimated by 50−250%. Except for LDA (GTF and PAW) and PBEsol-GTF, applying the dispersion corrections leads to an increase of this overestimation of ∼25−90%. In fact, the corrected dg (see Table 4) with dispersion corrections are smaller than those obtained without correction. The interactions between graphene planes are then repulsive. Since these interactions are mainly between the π orbitals of the planes, this increases the splitting. This is coherent with the reverse behavior of the LDA (GTF and PAW) and PBEsol-GTF results, where the splittings are lowered due to the increase of dg with the dispersion corrections applied. • The different band gaps at the M and K points show exactly the same trends with respect to the variation of dg as the Γ3+−Γ2− separation (as dicussed previously), since the bands are π bands, too. The only difference is the range of variation, which here is ∼10−50%. The M1+−M5+ gap is well reproduced at the LDA-D2/D3, PBE-D2/D3, and PBEsol-D2/D3 levels with GTF (PAW underestimate it); the hybrid functionals and HF overestimate it whatever the basis set and dispersion correction. However, for the K2−K4 gap the best agreement is found at the B3LYP-D2/D3 levels whatever the basis set. The LDA and GGA functionals under-

0.04 0.08 0.11

−0.08 −0.10 −0.06

−0.05 −0.03 −0.12 0.13 0.17 0.21

−0.00 −0.02 0.15 0.17 0.31 −0.01 0.02 0.03 −0.01

a

For both GTF and PAW, the results not including the zero point energy (ΔEnoZPE) are provided; for GTF, the relative energies not corrected for BSSE and not including the zero point energy (ΔEBSSE noZPE) are also given. Other theoretical results and experimental data are given for comparison. bReference 20. cReference 11. dReferences 13 and 89. London-type dispersion with the Ortmann et al. parametrization. eReference 27. fDeduced from the standard enthalpies of solids and gas extrapolated at T = 0 K of ref 8.

and interband separations of graphite obtained with the various Hamiltonians and basis sets, respectively. From a general point of view, the determined data are in good agreement with previous calculations. The obtained band gaps for diamond show the characteristic trends of the ones calculated for ionic or semiconducting compounds: HF overestimates Eg; LDA, PBE, and PBEsol give similar results and underestimate it; the best agreement is found with the hybrid functionals B3LYP and PBE0. However, the obtained bandwidths for graphite are characteristic for metallic or semimetallic materials: HF as well as the hybrids B3LYP and PBE0 fail to describe this type of system91 and consequently overestimate the bandwidth; in that case, the I


Article


Table 8. Computed Band Structures (in eV) at the Γ, K, and M Points of HOPG Obtained with Different Hamiltonians and Basis Setsa valence bands Hamiltonian

basis

Γ1+−Γ5+



22.88 23.02 22.99 16.79 16.64 16.66 16.67 16.62 16.59 16.66 16.65 16.48 16.67 16.76 16.72 16.72 16.56 16.73 18.40 18.46 18.45 17.76 17.83 17.82

NC-PP/PW LAPW US-PP/PW FLAPW LMTO LCGTO-FF LCGTO-FF

16.7 16.6 16.7 15. 15.8 16.21 16.28 16.00h

LDAb LDAc LDAd LDAe LDAf LDAg LDAg exp.

Γ1+−Γ3+ 16.48 16.03 16.02 10.64 10.95 10.96 10.84 10.87 11.40 11.21 11.21 11.70 11.13 11.05 11.14 11.14 11.28 11.03 12.60 12.43 12.43 12.33 12.11 12.11 other works 11.2 10.9 10.9 10.9 11.4 10.99 10.95 11,j 12.5h

a

b

intervalence−conduction bands Γ3+−Γ2−

Γ5+−Γ1+

M3+−M4−

K2−K4

1.46 2.68 2.70 3.12 2.09 2.07 2.16 2.10 1.47 2.00 2.00 0.45 2.02 2.31 2.03 2.04 1.32 2.07 1.69 2.16 2.16 1.51 2.15 2.15

22.60 22.60 22.58 11.42 11.43 11.45 6.87 6.87 11.38 11.38 11.37 6.11 7.32 11.34 11.35 11.35 6.53 7.02 14.09 14.09 14.08 13.71 13.71 13.70

11.72 11.89 11.89 4.97 4.65 4.66 4.31 4.29 4.52 4.63 4.63 3.25 4.02 4.71 4.64 4.64 4.08 4.29 6.37 6.45 6.45 5.96 6.07 6.07

1.55 2.74 2.75 2.26 1.44 1.43 1.47 1.43 1.00 1.39 1.39 1.81 0.24 1.63 1.42 1.42 0.85 1.42 1.34 1.72 1.72 1.16 1.67 1.67

2.1 2.1 2.1 2. 1.4 2.03 1.98 0.9j

7.1 6.9 6.8 8.4 9.1 7.08 7.2 11.5I c

d

4.6k e

1.6l f

Other theoretical results and experimental data are given for comparison. Reference 18. Reference 19. Reference 20. Reference 21. Reference 22. gReference 23. hReference 117. IReference 118. jReference 119. kReference 120. lReference 121.

Figure 5. Band structures and density of states (DOS) of HOPG obtained at the PBE level. On the left, the calculated band structures with GTF (full lines) and PAW (dashed lines) are shown; on the right, the projected density of states on the 2s2px2py AOs (PDOS spxpy) and on the pz AOs (PDOS pz) as well as the total density of states (DOS) are shown.

J


Article


Figure 6. Relative errors of the bandwidths with respect to the experimental data92,123 for different Hamiltonians and basis sets without dispersion correction.

Figure 7. Raman and IR spectra of HOPG obtained at the PBE-GTFD3 level. Intensities are only schematic and obtained from MOLDRAW.75−77

estimate it by ∼13%, but the obtained results may depend on various numerical errors, since this gap is small (exp. 1.6 eV). The obtained results for the Γ5+−Γ1+ gap require detailed comments. The best agreement is found for GTF at the LDA level with similar results for the PBE and PBEsol functionals. However, the PAW approaches systematically underestimate this data by ∼50%. This well-known discrepancy21−23,92 is still under discussion and was never explained consistently: Some attributed it to a poor sampling of the k space;21 Boettger23 claimed a numerical instability in the px and py states near the Γ point induced by the basis set; Holzwarth et al.,92 using a mixed basis set of PW and LCAO, found the same discrepancy and attributed it to the inability of the experimental measurements to detect these states because of their low density. Considering the band structures of Figure 5, for PAW, a projection of the wave function on the localized atomic states shows that the major contribution stems from s orbitals; the same as in ref 92 but with GTF, these states are far higher in energy, in good agreement with the available experimental data. Interestingly, Tatar et al.93 qualitatively obtained the same band structure as our PAW band structure with a muffin-tin potential and a modified Korringa−Kohn− Rostoker (KKR) approach; however, when they applied the full crystal potential, the Γ1+ state moved further up in the conduction band, and the final result is in good agreement with our GTF band structure. Therefore, considering our and published results, the Γ5+−Γ1+ gap cannot be used as a reference to establish the reliability of a method. Finally, taking into account the difficulties to reproduce the π bands and the Γ5+−Γ1+ gap, the best agreement on the band structure is obtained for the LDA-D2/D3, PBE-D2/D3, and PBEsol-D2/D3 functionals, which give similar results whatever the basis set. Frequencies and Elastic Constants. Frequencies. HOPG possesses nine transverse optical vibrational modes with the following irreducible representations at the Γ point

modes are inactive; however, one was determined by neutron scattering at 127 cm−1,94 whereas the other is expected to be close to the A2u mode at 868 cm−1.95 Table 9 gives the calculated frequencies of the various vibrations with different Hamiltonians and basis sets. Our results are coherent with the established influences of the Hamiltonian on vibrational properties.97 The obtained relative errors are presented in Figure 8. It is not an easy task to estimate the influence of the BSSE on these data; in this section, the comparison between GTF and PAW frequencies will allow us to understand it better. The modes can be subdivided into three classes: the interlayer, parallel intralayer, and perpendicular intralayer modes. The parallel intralayer modes (with E1u and E2g representations) correspond to the active Raman and IR vibrations of the covalent bonds parallel to (or in) the graphitic planes with experimental frequencies around 1585 cm−1. The obtained parallel intralayer frequencies show the same characteristic trends as those calculated for semiconducting compounds, as illustrated by the results for the Raman-active mode of diamond in Table 5: as shown in Figure 8, HF overestimates it by ∼15% and the hybrid functionals B3LYP and PBE0 by ∼2% and 5%, respectively; the error of LDA is around 1%, and the best agreement is found at the PBE and PBEsol levels with less than 0.6% error with respect to experiments. For theses modes the BSSE has no important effect: the GTF frequencies are ∼3% higher than the PAW ones. As expected, the dispersion corrections have also a negligible influence on these data (∼0.5% and 2% for GTF and PAW, respectively). The interlayer modes (with B2g and E2g representations) correspond to the Raman-active and silent vibrations of the weak interlayer bonds with frequencies measured experimentally at 42 and 127 cm−1, respectively. Since these modes are directly associated with dg, the reliability of the calculated data is driven by the BSSE and the difficulties of taking the dispersion interactions into account. (i) The inf luence of BSSE. In order to evaluate the impact of the BSSE on the final results, we consider now the data without dispersion correction. As it can be noticed, the BSSE has an important effect on the determined frequencies: due to the BSSE, which overbinds the graphitic planes, GTF systematically overestimates the frequencies by 25−90% with respect to the PAW ones. (ii) The inf luence of the dispersion. The same as dg, except for LDA and PBEsol-GTF, adding the dispersion corrections (or adding an attractive potential between planes) leads

Γ = A 2u + 2B2g + E1u + 2E 2g

The different Γ-phonon’s modes are depicted in Figure 1. Figure 7 shows the IR and Raman spectra obtained with PBED3-GTF. The two E2g modes are Raman active (exp. 42 and 1582 cm−1, respectively);94,95 the A2u and E1u modes are IR active (exp. 868 and 1588 cm−1, respectively).96 The two B2g K


Article


Table 9. Frequencies (in cm−1) of the Vibrational Modes of HOPG Obtained with Different Hamiltonians and Basis Setsa Hamiltonian

basis

E2g



25.1 48.9 51.1 60.3 38.0 36.1 41.6 40.5 25.5 37.4 37.5 28.8 47.7 44.0 37.5 37.1 28.4 47.2 29.4 39.5 39.9 25.0 38.5 39.3

LDAb LDAc LDAd PBEd LDAe PBEsole HSE06e exp.

PAW NC-PP/PW NC-PP/PW NC-PP/PW PAW PAW PAW

23.1 18.9 42.1 42f

B2g 85.2 160.1 183.1 190.3 118.5 113.9 100.1 94.4 95.4 126.0 139.4 55.9 133.9 123.7 101.8 101.4 31.4 111.2 94.8 121.8 130.9 99.5 147.9 160.9 other works

108.8 61.0 97.6 127f

B2g

A2u

E2g

E1u

853.2 823.4 824.3 590.6 646.2 642.1 883.0 883.4 705.0 682.1 684.0 873.1 867.8 647.0 661.0 660.2 872.7 867.8 728.4 707.7 709.8 744.6 720.5 722.8

851.8 814.9 815.7 594.4 648.1 644.9 888.5 888.6 706.5 683.8 685.3 873.3 872.4 648.9 662.7 661.9 874.5 872.6 727.9 706.3 708.2 745.6 720.7 722.7

1760.7 1758.1 1757.4 1606.7 1608.3 1614.0 1593.6 1593.3 1575.5 1575.6 1573.8 1527.2 1557.4 1591.5 1592.0 1592.7 1563.4 1574.9 1649.1 1649.4 1647.1 1614.4 1615.2 1611.5

1763.2 1765.9 1765.0 1621.6 1614.9 1620.4 1596.9 1596.3 1579.0 1582.2 1580.3 1527.2 1560.5 1560.0 1598.6 1599.3 1564.6 1578.1 1653.3 1656.1 1653.9 1617.8 1622.2 1618.5

1593.6 1580.0 1652.7 1582h

1597 1555 1597 1569 1624.7 1608.7 1665.2 1588g

885.1 872.4 910.0

890 890 893 884 889.5 874.2 913.2 868g

a In Figure 1, it is described which modes are silent, Raman, or infrared active. Other theoretical results and experimental data are given for comparison. bReference 24. cReference 25. dReference 26. eReference 27. fReference 94. gReference 96. hReference 95.

Figure 8. Relative errors of the frequencies with respect to the experimental data94−96 for different Hamiltonians and basis sets: (a) without, (b) with D2, and (c) with D3 dispersion correction.

with less than 5% and 10% of difference for the E2g and B2g vibrations, respectively. To conclude on the interlayer modes, the best agreement with respect to experiment is found at the PBE-GTF/PAW-D2/D3 and PBE0-GTF-D3 levels. The perpendicular intralayer (or out-of-plane) modes (with B2g and A2u representations) correspond to the IR-active and silent vibrations, mainly associated with the stretching of the

to a blue shift of these frequencies. Conversely, for LDA and PBEsol-GTF, which underestimate dg, the dispersion corrections lead to a red shift of the frequencies. Except for LDA-PAW and PBEsol-GTF without dispersion correction, the error ranges from 15% to 50%, and from 20% to 75% for the E2g and B2g modes, respectively. The dispersion corrections lead to a halving of these errors. The D2 and D3 corrections give similar results L


Article


Table 10. Elastic Constants (Cij) and Reuss Average of the Bulk Modulus (KR, in GPa) of HOPG Obtained with Different Hamiltonians and Basis Setsa Hamiltonian

basis

C11



1156.0 1264.6 1314.6 1211.3 1119.0 1112.1 1103.1 1097.4 1024.4 1080.8 1104.9 814.7 1058.3 1120.3 1094.1 1102.2 985.7 1072.2 1106.6 1153.4 1174.3 1062.2 1126.8 1154.2

C12

LDAb LDAc vdW-DF2(revPBE)d vdW-DF2(PBE)d PW91+vdWe PBEf LDAf PBEDFPTf LDADFPTf exp.

PAW LCGTO-FF PAW PAW US-PP/PW US-PP/PW US-PP/PW US-PP/PW US-PP/PW

1104.8

188.2 206.6 226.2 229.0 211.2 207.5 217.2 216.4 182.1 192.1 202.1 156.5 194.1 207.0 202.2 205.0 192.4 206.0 198.8 207.4 215.7 181.1 192.4 203.8 other works 203.9 1298j 1214j 1227j 1265j 976j 1283j

1079 1118 1060 ± 20g

C13

C33

C44

KR

0.5 −5.5 0.2 0.9 5.2 1.3 −1.2 −0.9 0.5 −1.9 0.1 0.8 −6.4 0.0 0.7 0.3 0.4 −3.0 0.1 −1.7 0.0 0.5 −2.5 0.6

17.5 54.2 71.4 67.0 31.0 26.2 31.4 29.7 20.6 32.8 39.4 1.9 28.9 30.7 21.5 19.2 8.9 26.1 20.0 31.1 34.7 22.8 46.1 54.0

1.9 3.6 3.0 9.1 5.3 5.4 4.3 4.1 1.6 2.3 2.2 0.2 3.8 4.6 3.7 3.7 1.7 4.1 2.2 3.1 3.0 1.6 2.2 2.2

17.1 49.7 65.4 61.5 30.0 25.3 29.9 28.4 20.0 31.0 37.2 1.9 27.1 29.4 20.8 18.7 8.8 24.9 19.4 29.6 33.1 22.0 42.8 50.1

−2.5 −0.5

30.9 40.8 25.0 36.4 41.7 2.4 29.0 42.2 29.5 36.5 ± 0.1g

5.6

32.5 38.3

−0.5 −2.8 217 235 180 ± 20g

a

15 ± 5g b

c

d

2.0 4.1

3.9 4.5 4.5 ± 0.05g

33.8h, 35.8I

e

Other theoretical results and experimental data are given for comparison. Reference 28. Reference 23. Reference 17. References 13 and 89. London-type dispersion with the Ortmann et al. parametrization. fReference 15. gReference 122. hReference 100. IReference 101. jThese data correspond to C11 + C12.

covalent bonds perpendicular to the graphitic planes, with −1

experimentally measured frequencies near 870 cm

(ii) The inf luence of the dispersion. Due to the BSSE, the dispersion corrections do not have the same influences on the GTF and PAW results. Considering PAW, the dispersion corrections have a negligible effect of less than 1%. However, the induced corrections on the GTF frequencies are much more important, ranging from 5% to 10%. D2 and D3 dispersion give similar results. As it was discussed above, the unimportant influence of the dispersion corrections on HF frequencies is symptomatic for the lack of correlation. To conclude the study on the perpendicular intralayer modes, the best agreement is found for the PBE/PBEsol-PAW-D2/D3 approaches. Finally, considering the overall frequencies, whatever the basis set, the PBE-D2/D3 approaches give the best description of the vibrational properties. Elastic Constants. In the Voigt notation, the elastic properties of HOPG are characterized by six elastic constants: C11, C12, C13, C33, C44, and C66. In this section, the performance of the various Hamiltonians and basis sets with respect to C66

(see Figure

1). They are also strongly influenced by the BSSE and dispersion corrections. (i) The inf luence of BSSE. As previously found, in order to evaluate the impact of the BSSE on the final results, we consider now the data without dispersion correction. As can be noticed, except for HF, GTF’s frequencies are systematically underestimated with respect to the PAW ones by ∼25%. These trends are not in contradiction with the results on interlayer modes: In fact, as shown previously, the BSSE overbinds the atoms of two differents planes; the induced effect is equivalent to adding an “artificial attractive force” between these atoms which softens these stretching modes. HF does not show this trend: Because of the lack of correlation, the electronic cloud is overconcentrated around the bonds which enhance their covalent character and attenuate the BSSE effects. M


Article


Figure 9. Relative errors of the elastic constants with respect to the experimental data at 4 K79 for different Hamiltonians and basis sets: (a) without, (b) with D2, and (c) with D3 dispersion correction.

will not be discussed since C66 = (C11 − C12)/2. The elastic constants (Cij) and Reuss average of the bulk modulus (KR) are given in Table 10. The obtained errors for these data with respect to experimental results are shown in Figure 9. Our results are in line with previous calculations. Likewise as for the vibrational modes, the elastic constants can be subdivided in three families: Constants associated with the response to the deformations parallel to the graphitic planes (C11 and C12), those associated with the response to the deformations along the c axis (C33), and the elastic constants associated with the shears (C13 and C44). For C11 and C12, which are determined from the deformation of the covalent bonds parallel to the graphitic layers, the obtained results are coherent with those for diamond in Table 5: Without dispersion correction, HF, LDA, PBEsol, and PBE0 overestimate them by ∼10%; PBE underestimates them by 4%; the best agreement is found for B3LYP. It can be seen that PBE-PAW and PBEsol-PAW results compare poorly to experimental data: As it is shown by Mounet and Marzari,15 certainly due to error propagation linked to the complexity of these functionals, these data depend strongly on the used method to determine it; for instance, if they are calculated from the phonon dispersions, C11 and C12 are in better agreement with experiment (see in Table 10, data indexed with DFPT). Mounet and Marzari show that these two constants, which are deduced from the high-frequency vibrations (around 1580 cm−1), are already well described without dispersion correction: thus, taking into account the dispersion interactions yields an average overestimation of these data by ∼10% with respect to experiment, whatever the Hamiltonian and basis set. However, C13, C33, and C44 are very sensitive to the various computational parameters, and they do not show the same dependencies on the chosen Hamiltonians and basis sets. We consider now C33. First, the role of the basis set is clear: the values obtained with GTF, which overbind the layers, are higher than the PAW ones; this discrepancy can be really huge; for example, at the PBE level the difference between basis sets reaches a ratio of 20. Without dispersion correction, except for LDA-GTF due to the already discussed overbinding, all Hamiltonians underestimate it by 14−50%; LDA overestimates it by ∼85%. As discussed above, PBE-PAW and PBEsol-PAW results compare poorly to experimental data. As it can be noticed, whatever the Hamiltonian and basis set, adding the dispersion corrections not only corrects but even overcorrects the obtained data for C33 in certain cases. For instance, the HF correction reaches an overestimation of C33 by 100% with respect to experiment. The magnitude of the D2/

D3 corrections depends on the Hamiltonian but without any clear trend. For this elastic constant, the best agreement is found for PBE-GTF-D2/D3. The C44 constants have globally the same dependencies on the Hamiltonians, basis sets and dispersion corrections as C33. Since this property is a function of the second derivative of the total energy of the material, due to its low value (4.5 GPa exp.) this data is strongly influenced by the various parameters used to determine it.15 Thus, when the obtained data without dispersion correction is in good agreement with experimental data (such as PBEsol-GTF), adding the dispersion corrections worsens the results. The best results are obtained for LDAPAW and PBEsol-GTF without dispersion corrections, with LDA-PAW-D2, and PBEsol-PAW-D2. As it was pointed out in ref 15, whatever the Hamiltonian, basis set, and dispersion correction, C13 is more than poorly reproduced and in certain cases its value is negative. Obviously, this disagreement cannot be only associated with the wrong description of the interlayer distance: The D3 correction allows us to obtain positive values of C13 between 0.1 and 1; however, the D2 correction can yield positive as well as negative values of C13 without clear tendency. The work of Mounet and Marzari15 tends to show that using phonon’s dispersions cannot provide a good evaluation of this data. Actually, to the best of our knowledge, there exists no method which is able to correct it. Since KR is a linear combination of Cij, the obtained results combine all errors discussed above. Thus, at this point the different data for KR will not be commented on in detail: the PBE-GTF-D2/D3 and PBE0-GTF-D3 approaches have the best agreement with experimental data. Finally, taking into account the difficulties to reproduce the shears elastic constants, PBE-GTF-D3 represents the best compromise to obtain a reasonable description of the elastic properties of graphite.

■

CONCLUSIONS In this paper, a systematic study of the performance of six Hamiltonians (HF, LDA, PBE, PBEsol, PBE0, and B3LYP) with two types of basis sets (GTF and PAW) has been performed with respect to various properties of HOPG and diamond. The investigated physical and chemical data of these polymorphs of carbon are the equilibrium geometry, some thermodynamic data (such as the cohesive and the interlayer interaction energies), the band structure, the Raman and IR frequencies, and the elastic constants. For GTF, the influence of the basis set superposition error on the obtained results was determined. As the used methods do not take dispersion N


Article


PAW are more suitable to describe these compounds: In certain cases the BSSE allows one to improve the description of given properties; in others cases it worsens the results. However, if dispersion corrections are included, the results are of comparable quality. With respect to the dispersion correction, D3 globally performs better than D2 (as it was also mentioned for other types of materials51,52) and should therefore be the method of choice for the HOPG system. To conclude, based on the overall results, PBE-GTF-D3 is the best compromise to describe the properties of HOPG.

interactions correctly into account, two different versions of Grimme’s dispersion correction, D2 and D3, were evaluated. Their scaling parameters were optimized to reproduce the experimental interlayer spacing with an error of less than 0.5%. For the cohesive energy of HOPG, the best agreements with experimental data were obtained for B3LYP-GTF-D3 and PBE0-GTF. The relative stability of graphite in comparison to diamond was determined. Except for LDA-PAW, whatever the Hamiltonian and basis set, graphite is the ground state of carbon; the energy difference between the two phases, graphite and diamond, is systematically overestimated with respect to experiment. The dispersion corrections have a negligible effect on these properties mainly driven by the covalent bonds. However, the interlayer interaction energy, depending on the proper description of the interlayer interaction, is the property which is influenced the most by the dispersion corrections. D3 outperforms D2; the best agreements with experiment were found for PBE with either basis set and B3LYP-GTF-D3. Considering the band structures, the bandwidths and interband separations associated with the σ bonds are not influenced by the dispersion corrections. However, the obtained errors are important for the splitting of the π bands. This splitting is systematically overestimated. Applying the dispersion corrections leads to an increase of this overestimation. In fact, except for LDA and PBEsol-GTF, the dg obtained with dispersion corrections are smaller than those obtained without dispersion correction. When dispersion corrections are applied, the interactions between graphene planes obtained within the SCF cycle are repulsive, while LDA and PBEsol-GTF show the reverse trend. Taking into account the difficulties to reproduce the π bands, the best agreements on the band structure are obtained for the LDA-D2/D3, PBED2/D3, and PBEsol-D2/D3 functionals, which give similar results whatever the basis set. The vibrational frequencies of graphite can be separated into three classes: The interplane modes, the parallel intraplane modes, and the perpendicular intraplane modes. The parallel intraplane modes are reproduced rather well without any kind of dispersion correction. However, interplane modes are largely improved upon using a dispersion correction. Overall the GGA functionals perform best with respect to experimental data, PBE-GTF/PAW-D2/D3 yielding the best results. Considering the elastic properties, for C11 and C12, the best agreement is found for B3LYP, PBE providing data close to the former mentioned. The shear’s constant C44 and the two constants associated with deformations along the c axis, C33 and C13, are the most difficult to reproduce: For C44 and C33, adding the dispersion corrections not only corrects but even overcorrects the obtained data in certain cases, the magnitude of the change induced by the D2/D3 corrections show no clear trend with respect to the different Hamiltonians, and the C13 constant is systematically poorly reproduced without clear explanation. Globally, PBE-GTF-D3 allows one to obtain a reasonable description of the elastic properties of HOPG. The various obtained results illustrate some well-known characteristics of the various functionals: LDA overbinds and the other Hamiltonians underbind the graphitic planes. For GTF, the BSSE allows one to correct artificially for the lack of dispersion interaction. Therefore, the interlayer interactions are stronger than for PAW; thus, with respect to PAW and without dispersion corrections, the c parameter is underestimated, the interlayer frequencies are overestimated, and elastic constants, such as C33, are overestimated. It is difficult to state if GTF or

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +33 (0)160737608. Fax: +33 (0)160736889. Present Address §

N.C.F.-M.: Theory Group, Max Planck Institute for Polymer Research, Ackermannweg 10, G-55128 Mainz, Germany. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank the ANRT (French National Association for Research and Technology) for its financial support within the CIFRE agreement 2014/1036 (industrial convention for training through research).

■

ABBREVIATIONS BSSE basis set superposition error DFPT density functional perturbation theory GGA generalized gradient approximation GTF Gaussian-type functions HOPG highly oriented pyrolytic graphite PW plane waves PAW projector augmented waves FLAPW full-potential linearized augmented plane wave LDA local density approximation LAPW linearized augmented plane wave LCAO linear combination of atomic orbitals NC-PP norm conserving pseudo potential US-PP ultrasoft pseudo potential LMTO linear muffin-tin-orbital LCGTO-FF full-potenial linear combinations of Gaussian-type orbitals fitting function RPA random phase approximation TS Tkatchenko−Scheffler van der Waals correction method TS+SCS TS with self-consistent screening UNGG Uranium Naturel Graphite Gaz ZPE zero point energy

■

REFERENCES

(1) Kane, J.; Karthik, C.; Butt, D.; Windes, W.; Ubic, R. Microstructural Characterization and Pore Structure Analysis of Nuclear Graphite. J. Nucl. Mater. 2011, 415, 189−197. (2) Telling, R. H.; Heggie, M. I. Radiation Defects in Graphite. Philos. Mag. 2007, 87, 4797−4846. (3) Barker, J.; Gao, F. Carbonaceous Electrode and Compatible Electrolyte Solvent, 1998. http://www.google.com/patents/US5712059, U.S. Patent 5,712,059, 1998.

O


Article


Nanotubes. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 67, 035401. (25) Miyamoto, Y.; Cohen, M. L.; Louie, S. G. Ab Initio Calculation of Phonon Spectra for Graphite, BN, and BC2N Sheets. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 52, 14971−14975. (26) Wirtz, L.; Rubio, A. The Phonon Dispersion of Graphite Revisited. Solid State Commun. 2004, 131, 141−152. (27) Grochala, W. Diamond: Electronic Ground State of Carbon at Temperatures Approaching 0 K. Angew. Chem., Int. Ed. 2014, 53, 3680−3683. (28) Qi, Y.; Guo, H.; Hector, L. G.; Timmons, A. Threefold Increase in the Young’s Modulus of Graphite Negative Electrode during Lithium Intercalation. J. Electrochem. Soc. 2010, 157, A558−A566. (29) Teobaldi, G.; Tanimura, K.; Shluger, A. L. Structure and Properties of Surface and Subsurface Defects in Graphite Accounting for van der Waals and Spin-Polarization Effects. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 174104. (30) Letardi, S.; Celino, M.; Cleri, F.; Rosato, V. Atomic Hydrogen Adsorption on a Stone-Wales Defect in Graphite. Surf. Sci. 2002, 496, 33−38. (31) Lebègue, S.; Harl, J.; Gould, T.; Á ngyán, J. G.; Kresse, G.; Dobson, J. F. Cohesive Properties and Asymptotics of the Dispersion Interaction in Graphite by the Random Phase Approximation. Phys. Rev. Lett. 2010, 105, 196401. (32) Berland, K.; Borck, Ø.; Hyldgaard, P. Van der Waals Density Functional Calculations of Binding in Molecular Crystals. Comput. Phys. Commun. 2011, 182, 1800−1804. (33) Bučko, T.; Lebègue, S.; Hafner, J.; Á ngyán, J. G. TkatchenkoScheffler van der Waals Correction Method with and without SelfConsistent Screening Applied to Solids. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 064110. (34) Kristyán, S.; Pulay, P. Can (Semi)local Density Functional Theory Account for the London Dispersion Forces? Chem. Phys. Lett. 1994, 229, 175−180. (35) Pérez-Jordá, J. M.; Becke, A. A Density-Functional Study of van der Waals Forces: Rare Gas Diatomics. Chem. Phys. Lett. 1995, 233, 134−137. (36) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (37) Thonhauser, T.; Cooper, V. R.; Li, S.; Puzder, A.; Hyldgaard, P.; Langreth, D. C. Van der Waals Density Functional: Self-Consistent Potential and the Nature of the van der Waals Bond. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 125112. (38) Román-Pérez, G.; Soler, J. M. Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes. Phys. Rev. Lett. 2009, 103, 096102. (39) Lee, K.; Murray, E. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy van der Waals Density Functional. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 081101. (40) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (41) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (42) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (43) Tkatchenko, A.; Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and FreeAtom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. (44) Tkatchenko, A.; DiStasio, R. A.; Car, R.; Scheffler, M. Accurate and Efficient Method for Many-Body van der Waals Interactions. Phys. Rev. Lett. 2012, 108, 236402. (45) Civalleri, B.; Zicovich-Wilson, C. M.; Valenzano, L.; Ugliengo, P. B3LYP Augmented with an Empirical Dispersion Term (B3LYP-

(4) Bessel, C. A.; Laubernds, K.; Rodriguez, N. M.; Baker, R. T. K. Graphite Nanofibers as an Electrode for Fuel Cell Applications. J. Phys. Chem. B 2001, 105, 1115−1118. (5) Park, D. H.; Zeikus, J. G. Improved Fuel Cell and Electrode Designs for Producing Electricity from Microbial Degradation. Biotechnol. Bioeng. 2003, 81, 348−355. (6) Gupta, B.; Janting, J.; Sørensen, G. Friction and Wear Characteristics of Ion-Beam Modified Graphite Coatings. Tribol. Int. 1994, 27, 139−143. (7) Shaji, S.; Radhakrishnan, V. An Investigation on Surface Grinding Using Graphite as Lubricant. Int. J. Mach. Tool. Manu. 2002, 42, 733− 740. (8) Lide, D. Handbook of Chemistry and Physics; CRC Press: Boston, 1990−1991; Vol. 71. (9) Lee, I.-H.; Martin, R. M. Applications of the GeneralizedGradient Approximation to Atoms, Clusters, and Solids. Phys. Rev. B: Condens. Matter Mater. Phys. 1997, 56, 7197−7205. (10) Hamann, D. R. Generalized-Gradient Functionals in Adaptive Curvilinear Coordinates. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1568−1574. (11) Janotti, A.; Wei, S.-H.; Singh, D. J. First-Principles Study of the Stability of BN and C. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 64, 174107. (12) Rydberg, H.; Dion, M.; Jacobson, N.; Schröder, E.; Hyldgaard, P.; Simak, S. I.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for Layered Structures. Phys. Rev. Lett. 2003, 91, 126402. (13) Ortmann, F.; Bechstedt, F.; Schmidt, W. G. Semiempirical van der Waals Correction to the Density Functional Description of Solids and Molecular Structures. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 205101. (14) Ooi, N.; Rairkar, A.; Adams, J. B. Density Functional Study of Graphite Bulk and Surface Properties. Carbon 2006, 44, 231−242. (15) Mounet, N.; Marzari, N. First-Principles Determination of the Structural, Vibrational and Thermodynamic Properties of Diamond, Graphite, and Derivatives. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 205214. (16) Bučko, T.; Hafner, J.; Lebègue, S.; Á ngyán, J. G. Improved Description of the Structure of Molecular and Layered Crystals: Ab Initio DFT Calculations with van der Waals Corrections. J. Phys. Chem. A 2010, 114, 11814−11824. (17) Gulans, A.; Krasheninnikov, A. V.; Puska, M. J.; Nieminen, R. M. Bound and Free Self-Interstitial Defects in Graphite and Bilayer Graphene: A Computational Study. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 024114. (18) Charlier, J.-C.; Gonze, X.; Michenaud, J.-P. First-Principles Study of the Electronic Properties of Graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 4579−4589. (19) Schabel, M. C.; Martins, J. L. Energetics of Interplanar Binding in Graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46, 7185−7188. (20) Furthmüller, J.; Hafner, J.; Kresse, G. Ab Initio Calculation of the Structural and Electronic Properties of Carbon and Boron Nitride Using Ultrasoft Pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 15606−15622. (21) Jansen, H. J. F.; Freeman, A. J. Structural and Electronic Properties of Graphite via an All-Electron Total-Energy Local-Density Approach. Phys. Rev. B: Condens. Matter Mater. Phys. 1987, 35, 8207− 8214. (22) Ahuja, R.; Auluck, S.; Trygg, J.; Wills, J. M.; Eriksson, O.; Johansson, B. Electronic Structure of Graphite: Effect of Hydrostatic Pressure. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 4813− 4819. (23) Boettger, J. C. All-Electron Full-Potential Calculation of the Electronic Band Structure, Elastic Constants, and Equation of State for Graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1997, 55, 11202− 11211. (24) Dubay, O.; Kresse, G. Accurate Density Functional Calculations for the Phonon Dispersion Relations of Graphite Layer and Carbon P


Article

The Journal of Physical Chemistry C D*) as Applied to Molecular Crystals. CrystEngComm 2008, 10, 405− 410. (46) Ugliengo, P.; Zicovich-Wilson, C. M.; Tosoni, S.; Civalleri, B. Role of Dispersive Interactions in Layered Materials: A Periodic B3LYP and B3LYP-D* Study of MgOH2, CaOH2 and Kaolinite. J. Mater. Chem. 2009, 19, 2564−2572. (47) Canepa, P.; Ugliengo, P.; Alfredsson, M. Elastic and Vibrational Properties of α- and β-PbO. J. Phys. Chem. C 2012, 116, 21514− 21522. (48) Delbé, K.; Mansot, J.-L.; Thomas, P.; Baranek, P.; Boucher, F.; Vangelisti, R.; Billaud, D. Contribution to the Understanding of Tribological Properties of Graphite Intercalation Compounds with Metal Chloride. Tribol. Lett. 2012, 47, 367−379. (49) Brandenburg, J. G.; Grimme, S.; Jones, P. G.; Markopoulos, G.; Hopf, H.; Cyranski, M. K.; Kuck, D. Unidirectional Molecular Stacking of Tribenzotriquinacenes in the Solid State: A Combined X-Ray and Theoretical Study. Chem. - Eur. J. 2013, 19, 9930−9938. (50) Reckien, W.; Eggers, M.; Bredow, T. Theoretical Study of the Adsorption of Benzene on Coinage Metals. Beilstein J. Org. Chem. 2014, 10, 1775−1784. (51) Ilawe, N. V.; Zimmerman, J. A.; Wong, B. M. Breaking Badly: DFT-D2 Gives Sizeable Errors for Tensile Strengths in PalladiumHydride Solids. J. Chem. Theory Comput. 2015, 11, 5426−5435. PMID: 26574331. (52) Skelton, J. M.; Tiana, D.; Parker, S. C.; Togo, A.; Tanaka, I.; Walsh, A. Influence of the Exchange-Correlation Functional on the Quasi-Harmonic Lattice Dynamics of II-VI Semiconductors. J. Chem. Phys. 2015, 143, 064710. (53) George, J.; Reimann, C.; Deringer, V. L.; Bredow, T.; Dronskowski, R. On the DFT Ground State of Crystalline Bromine and Iodine. ChemPhysChem 2015, 16, 728−732. (54) Pisani, C.; Dovesi, R.; Erba, A.; Giannozzi, P. In Modern ChargeDensity Analysis; Gatti, C., Macchi, P., Eds.; Springer: Netherlands, 2012; pp 79−132. (55) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. R.; Zicovich-Wilson, C. M. CRYSTAL: A Computational Tool For Ab Initio Study of the Electronic Properties of Crystals. Z. Kristallogr. Cryst. Mater. 2005, 220, 571−573. (56) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J. et al. CRYSTAL09 User’s Manual; University of Torino: Torino, 2009. (57) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys.: Condens. Matter 2009, 21, 395502. (58) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (59) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula Into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (60) Moreira, I. d. P. R.; Illas, F.; Martin, R. L. Effect of Fock Exchange on the Electronic Structure and Magnetic Coupling in NiO. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 65, 155102. (61) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110, 6158−6170. (62) Dirac, P. A. M. Note on Exchange Phenomena in the Thomas Atom. Math. Proc. Cambridge Philos. Soc. 1930, 26, 376−385. (63) Perdew, J. P.; Zunger, A. Self-Interaction Correction to DensityFunctional Approximations for Many-Electron Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 5048−5079. (64) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (65) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett. 2008, 100, 136406.

(66) Dovesi, R.; Causà, M.; Orlando, R.; Roetti, C.; Saunders, V. R. Ab initio Approach to Molecular Crystals: A Periodic Hartree-Fock Study of Crystalline Urea. J. Chem. Phys. 1990, 92, 7402. (67) Catti, M.; Pavese, A.; Dovesi, R.; Saunders, V. R. Static Lattice and Electron Properties of MgCO3 (Magnesite) Calculated by Ab Initio Periodic Hartree-Fock Methods. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 9189−9198. (68) Zicovich-Wilson, C. LoptCG (Shell Procedure for Numerical Gradient Optimization); Instituto de Technologia Quimica Valencia: Valencia, 2006. (69) Hestenes, M.; Stiefel, E. Methods of Conjugate Gradients for Solving Linear Systems. J. Res. Natl. Bur. Stand. 1952, 49, 409−436. (70) Grimme, S. DFT-D3 - A Dispersion Correction for Density Functionals, Hartree-Fock and Semi-Empirical Quantum Chemical Methods, http://www.thch.uni-bonn.de/tc/index.php?section/ downloads&subsection/DFT-D3&lang/english, 2015, online; accessed July 30, 2015. (71) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188−5192. (72) Boys, S.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (73) Golesorkhtabar, R.; Pavone, P.; Spitaler, J.; Puschnig, P.; Draxl, C. ElaStic: A tool for Calculating Second-Order Elastic Constants from First Principles. Comput. Phys. Commun. 2013, 184, 1861−1873. (74) Kokalj, A. XCrySDen-A New Program for Displaying Crystalline Structures and Electron Densities. J. Mol. Graphics Modell. 1999, 17, 176−179. (75) Ugliengo, P.; Borzani, G.; Viterbo, D. MOLDRAW - Program for the Graphical Manipulation of Molecules on Personal Computers. J. Appl. Crystallogr. 1988, 21, 75. (76) Ugliengo, G.; Borzani, P.; Viterbo, D. New Developments in the MOLDRAW Graphics Program. Z. Kristallogr. 1988, 185, 712. (77) Ugliengo, P.; Viterbo, D.; Chiari, G. MOLDRAW: Molecular Graphics on a Personal Computer. Z. Kristallogr. - Cryst. Mater. 1993, 207, 9−24. (78) Persistence of Vision Raytracer, http://www.povray.org/, 2015, online; accessed Aug 16, 2015. (79) Baskin, Y.; Meyer, L. Lattice Constants of Graphite at Low Temperatures. Phys. Rev. 1955, 100, 544−544. (80) Baranek, P.; Lichanot, A.; Orlando, R.; Dovesi, R. Structural and Vibrational Properties of Solid Mg(OH)2 and Ca(OH)2 - Performances of Various Hamiltonians. Chem. Phys. Lett. 2001, 340, 362−369. (81) Labat, F.; Baranek, P.; Domain, C.; Minot, C.; Adamo, C. Density Functional Theory Analysis of the Structural and Electronic Properties of TiO2 Rutile and Anatase Polytypes: Performances of Different Exchange-Correlation Functionals. J. Chem. Phys. 2007, 126, 154703. (82) Ropo, M.; Kokko, K.; Vitos, L. Assessing the Perdew-BurkeErnzerhof Exchange-Correlation Density Functional Revised for Metallic Bulk and Surface Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 195445. (83) Csonka, G. I.; Perdew, J. P.; Ruzsinszky, A.; Philipsen, P. H. T.; Lebègue, S.; Paier, J.; Vydrov, O. A.; Á ngyán, J. G. Assessing the Performance of Recent Density Functionals for Bulk Solids. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 155107. (84) Sophia, G.; Baranek, P.; Sarrazin, C.; Rérat, M.; Dovesi, R. FirstPrinciples Study of the Mechanisms of the Pressure-Induced Dielectric Anomalies in Ferroelectric Perovskites. Phase Transitions 2013, 86, 1069−1084. (85) Benedict, L. X.; Chopra, N. G.; Cohen, M. L.; Zettl, A.; Louie, S. G.; Crespi, V. H. Microscopic Determination of the Interlayer Binding Energy in Graphite. Chem. Phys. Lett. 1998, 286, 490−496. (86) Zacharia, R.; Ulbricht, H.; Hertel, T. Interlayer Cohesive Energy of Graphite from Thermal Desorption of Polyaromatic Hydrocarbons. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 155406. (87) Ribeiro, F. J.; Tangney, P.; Louie, S. G.; Cohen, M. L. Structural and Electronic Properties of Carbon in Hybrid Diamond-Graphite Q


Article

The Journal of Physical Chemistry C Structures. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 214109. (88) Ribeiro, F. J.; Tangney, P.; Louie, S. G.; Cohen, M. L. Hypothetical Hard Structures of Carbon with Cubic Symmetry. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 172101. (89) The Effect of van der Waals Interactions on the Properties of Intrinsic Defects in Graphite. Carbon 2010, 48, 4145−4161.10.1016/ j.carbon.2010.07.029 (90) Paier, J.; Marsman, M.; Hummer, K.; Kresse, G.; Gerber, I. C.; Á ngyán, J. G. Screened Hybrid Density Functionals Applied to Solids. J. Chem. Phys. 2006, 124, 154709. (91) Paier, J.; Marsman, M.; Kresse, G. Why Does the B3LYP Hybrid Functional Fail for Metals. J. Chem. Phys. 2007, 127, 024103. (92) Holzwarth, N. A. W.; Louie, S. G.; Rabii, S. X-Ray Form Factors and the Electronic Structure of Graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1982, 26, 5382−5390. (93) Tatar, R. C.; Rabii, S. Electronic Properties of Graphite: A Unified Theoretical Study. Phys. Rev. B: Condens. Matter Mater. Phys. 1982, 25, 4126−4141. (94) Nicklow, R.; Wakabayashi, N.; Smith, H. G. Lattice Dynamics of Pyrolytic Graphite. Phys. Rev. B 1972, 5, 4951−4962. (95) Nemanich, R. J.; Solin, S. A. First- and Second-Order Raman Scattering from Finite-Size Crystals of Graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 1979, 20, 392−401. (96) Nemanich, R.; Lucovsky, G.; Solin, S. Infrared Active Optical Vibrations of Graphite. Solid State Commun. 1977, 23, 117−120. (97) Scott, A. P.; Radom, L. Harmonic Vibrational Frequencies: An Evaluation of Hartree-Fock, Møller-Plesset, Quadratic Configuration Interaction, Density Functional Theory, and Semiempirical Scale Factors. J. Phys. Chem. 1996, 100, 16502−16513. (98) Wyckoff, R. Crystal Structures; Interscience Publisher: New York, 1963; Vol. 1, p 27. (99) In The Structures of Elements; Donohue, J., Ed.; Wiley: New York, 1974. (100) Hanfland, M.; Beister, H.; Syassen, K. Graphite under Pressure: Equation of State and First-Order Raman Modes. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 39, 12598−12603. (101) Zhao, Y. X.; Spain, I. L. X-Ray Diffraction Data for Graphite to 20 GPa. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 40, 993−997. (102) Trucano, P.; Chen, R. Structure of Graphite by Neutron Diffraction. Nature 1975, 258, 136−137. (103) Svane, A. Hartree-Fock Band-Structure Calculations with the Linear Muffin-Tin-Orbital Method: Application to C, Si, Ge, and α -Sn. Phys. Rev. B: Condens. Matter Mater. Phys. 1987, 35, 5496−5502. (104) Causà, M.; Dovesi, R.; Roetti, C. Pseudopotential HartreeFock Study of Seventeen III-V and IV-IV Semiconductors. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 11937−11943. (105) Rohlfing, M.; Krüger, P.; Pollmann, J. Quasiparticle BandStructure Calculations for C, Si, Ge, GaAs, and SiC using GaussianOrbital Basis Sets. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 17791−17805. (106) Muscat, J.; Wander, A.; Harrison, N. On the Prediction of Band Gaps from Hybrid Functional Theory. Chem. Phys. Lett. 2001, 342, 397−401. (107) Pavone, P.; Karch, K.; Schütt, O.; Strauch, D.; Windl, W.; Giannozzi, P.; Baroni, S. Ab Initio Lattice Dynamics of Diamond. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 48, 3156−3163. (108) Milman, V.; Winkler, B.; White, J. A.; Pickard, C. J.; Payne, M. C.; Akhmatskaya, E. V.; Nobes, R. H. Electronic Structure, Properties, and Phase Stability of Inorganic Crystals: A Pseudopotential PlaneWave Study. Int. J. Quantum Chem. 2000, 77, 895−910. (109) Wen, B.; Zhao, J.; Bucknum, M. J.; Yao, P.; Li, T. Firstprinciples Studies of Diamond Polytypes. Diamond Relat. Mater. 2008, 17, 356−364. (110) McSkimin, H. J.; Andreatch, P. Elastic Moduli of Diamond as a Function of Pressure and Temperature. J. Appl. Phys. 1972, 43, 2944− 2948. (111) Warren, J. L.; Yarnell, J. L.; Dolling, G.; Cowley, R. A. Lattice Dynamics of Diamond. Phys. Rev. 1967, 158, 805−808.

(112) Saslow, W.; Bergstresser, T. K.; Cohen, M. L. Band Structure and Optical Properties of Diamond. Phys. Rev. Lett. 1966, 16, 354− 356. (113) Clark, C. D.; Dean, P. J.; Harris, P. V. Intrinsic Edge Absorption in Diamond. Proc. R. Soc. London, Ser. A 1964, 277, 312− 329. (114) Clark, C. The Absorption Edge Spectrum of Diamond. J. Phys. Chem. Solids 1959, 8, 481−484. (115) Geschneider, K., Jr. Solid State Physics; Academic: New York, 1964; Vol. 16, p 275. (116) Calzaferri, G.; Rytz, R. The Band Structure of Diamond. J. Phys. Chem. 1996, 100, 11122−11124. (117) Eberhardt, W.; McGovern, I. T.; Plummer, E. W.; Fisher, J. E. Charge-Transfer and Non-Rigid-Band Effects in the Graphite Compound LiC6. Phys. Rev. Lett. 1980, 44, 200−204. (118) Willis, R. F.; Fitton, B.; Painter, G. S. Secondary-Electron Emission Spectroscopy and the Observation of High-Energy Excited States in Graphite: Theory and Experiment. Phys. Rev. B 1974, 9, 1926−1937. (119) McGovern, I.; Eberhardt, W.; Plummer, E.; Fischer, J. The Bandstructures of Graphite and Graphite Intercalation Compounds as Determined by Angle Resolved Photoemission Using Synchrotron Radiation. Physica B+C 1980, 99, 415−419. (120) McClure, J. W. Band Structure of Graphite and de Haas-van Alphen Effect. Phys. Rev. 1957, 108, 612−618. (121) Willis, R. F.; Feuerbacher, B.; Fitton, B. Experimental Investigation of the Band Structure of Graphite. Phys. Rev. B 1971, 4, 2441−2452. (122) Blakslee, O. L.; Proctor, D. G.; Seldin, E. J.; Spence, G. B.; Weng, T. Elastic Constants of Compression-Annealed Pyrolytic Graphite. J. Appl. Phys. 1970, 41, 3373−3382. (123) Tatar, R.; Holzwarth, N.; Rabii, S. Energy Band Structure of Three Dimensional Graphite. Synth. Met. 1981, 3, 131−138.

R


First-Principles Study of the Structural, Electronic, Dynamic, and

Recommend Documents