Journal of the Korean Physical Society, Vol. 59, No. 1, July 2011, pp. 196∼199
Brief Reports
Ab initio Calculations with van der Waals Corrections: Benzene-benzene Intermolecular Case and Graphite Jinwoo Park Department of Physics and Graphene Research Institute, Sejong University, Seoul 143-747, Korea
Byung Deok Yu Department of Physics, University of Seoul, Seoul 130-743, Korea
Suklyun Hong∗ Department of Physics, Graphene Research Institute, and Institute of Fundamental Physics, Sejong University, Seoul 143-747, Korea (Received 29 April 2011) For the benzene-benzene intermolecular case and graphite, we have performed ab initio density functional theory (DFT) calculations with van der Waals corrections. We use two DFT methods: the Vienna ab initio simulation package (VASP) and the Fritz Haber Institute ab initio molecular simulations package (FHI-aims). For the van der Waals description, we adopt Grimme’s DFT-D2 approach in VASP while the Tkatchenko-Scheffler (TS) scheme is used in FHI-aims. First, we calculate the equilibrium distance between the in-plane benzene-dimer molecules by using the two methods. The equilibrium intermolecular distance obtained using the DFT-D2 scheme of Grimme is shorter than that obtained using the TS scheme. The equilibrium structural parameters of graphite are calculated by varying the lattice parameters a and c simultaneously and are compared with previous results. PACS numbers: 71.15.Nc, 61.46.-w Keywords: Van der Waals interation, Benzene, Graphite, Ab initio calculations DOI: 10.3938/jkps.59.196
plemented in two computational packages, VASP [3–6] and FHI-aims [7,8], respectively. These simulation packages with the vdW corrections were developed recently, so scientists can now use the vdW schemes to investigate organic and carbon-based materials. For reliability of these vdW schemes, we need to check how much different the calculation results using the Grimme and the TS vdW schemes are from each other. To check the validity of van der Waals corrections for the in-plane benzene dimer and graphite, we preformed DFT calculations as given in VASP and FHI-aims within the GGA [9]. In the VASP code, plane waves up to a cutoff energy (Ecut ) of 400 eV are included to expand the wave functions, and the ions are represented by using projector-augmented-wave (PAW) potentials [10,11] and Grimme’s DFT-D2 vdW as implemented in VASP. On the other hand, in the FHI-aims code, we use DFT with an all electron full-potential scheme and a “tight” numerical atom-centered orbitals (NAO) “tier2” basis set and adopt a van der Waals description implemented in the code [2]. In both calculations, the exchange-correlation energy is described by using the functional of Perdew, Burke, and Ernzerhof (PBE) based on the GGA [9].
Most theoretical scientists investigating electronic structure and chemical properties of materials use density functional theory (DFT) methods. For DFT calculations, proper choices of exchange-correlation functionals [generalized gradient approximation (GGA) or local density approximation (LDA)], potentials, basis sets, etc., should be made to obtain “good and reliable” computational results that are compatible with experimental ones. Despite their great success in most cases, DFT methods tend to fail in treating the long-range behavior of interatomic London dispersion interactions like van der Waals (vdW) interactions. The vdW forces play a very important role in descriptions of organic and carbon-based materials. Recently, two appropriate modelings of vdW forces were made by Grimme [1] and by Tkatchenko and Scheffler (TS) [2]. One, as given by Grimme’s approach, is based on a semi-empirical GGAtype theory, while the other is based on a non-empirical theory. The Grimme and TS vdW schemes are now im∗ Corresponding
author; E-mail:
[email protected]; Fax: +82-2-
3408-4316
-196-
Ab initio Calculations with van der Waals Corrections: Benzene-benzene · · · – Jinwoo Park et al.
-197-
Table 1. Calculated equilibrium lattice constants (a, c, and c/a) and bulk modulus (B0 ) for graphite in comparison with previous calculations and experiments.
VASP
FHI-aims Exp [14]
Ecut (eV) 400 500 600 700 800 900 1000 1000 [13] -
V (˚ A3 ) 33.81 33.94 33.95 33.93 33.92 33.91 33.91 33.80 34.94 35.2
a (˚ A) 2.46 2.39 2.46 2.46 2.46 2.46 2.46 2.46 2.46 2.462
First, we consider in-plane benzene-dimer molecules to investigate the validity of van der Waals corrections in describing the intermolecular interactions. We use a supercell of dimensions of 22 ˚ A × 30 ˚ A × 15 ˚ A (see Fig. 1). To calculate their binding energy as a function of distance d between two benzene molecules, we fix the y and z coordinates of two neighboring carbon atoms of two benzene molecules. Except for these eight coordinates, all coordinates are fully relaxed until the Hellmann-Feynman forces are lower than 0.01 eV/˚ A. For the Brillouin-zone (BZ) integration, we use the gamma point. A Gaussian broadening with a width of 0.01 eV is used to accelerate the convergence. To find the equilibrium distance between the in-plane benzene-dimer molecules, we calculate their binding energy as a function of distance (see Fig. 2). As a result, the equilibrium distances are 6.597 ˚ A and 6.682 ˚ A for the VASP and FHI-aims codes, respectively, while their corresponding binding energies are -0.0426 eV and -0.0453 eV. The FHI-aims code gives a slightly longer equilibrium distance (about 1.3%) and a more stable binding energy (6.3%). Without considering the vdWforce, the line shapes of FHI-aims and VASP are almost the same and do not give an equilibrium distance. The above results imply that TS’s vdW is more effective in dealing with long-range interactions than Grimme’s vdW because the tail shape of the FHI-aims vdW (red line) is lower and longer than that of the VASP vdW (green line in Fig. 2). Next, we consider graphite, which consists of layers of carbon with hexagonal crystal structures. Using the PBE-D2 approach in VASP, we determine its structural parameters (a = 2.46 ˚ A, c = 6.45 ˚ A). In the calculations, the kinetic energy cutoff is first set to be 600 eV, and a 16 × 16 × 8 k-point mesh in the BZ is used for the k-space integration. The optimized equilibrium lattice constants a and c are obtained simultaneously by considering the total energies at about 440 points over the a-c configuration space. The results are also tested by increasing a cutoff energy from 400 eV to 1000 eV (see Table 1).
c (˚ A) 6.43 6.87 6.45 6.46 6.46 6.45 6.45 6.45 6.67 6.707
c/a 2.61 2.87 2.62 2.63 2.63 2.62 2.62 2.62 2.71 2.724
B0 (GPa) 43.09 39.74 40.22 40.18 40.44 40.60 40.45 37.8 64.73 34-42
Ecoh (eV) 9.34 9.33 9.33 9.34 8.49 9.34 9.34 8.099 9.13 7.371
Fig. 1. (Color online) Atomic structure of an in-plane benzene dimer.
The obtained values for the structural parameters are compared to the PBE-D2 results by Bu˘cko et al. [12], in which a cutoff energy of 1000 eV and a 16 × 16 × 8 kpoint mesh in the BZ were used. Our equilibrium lattice constants a and c are almost the same as those of Bu˘cko et al. The difference in B0 and Ecoh may come from the fitting method over the a - c configuration space through the simultaneous optimization. Using the TS vdW scheme in FHI-aims, we also determine the structural parameters for graphite. The equilib-
-198-
Journal of the Korean Physical Society, Vol. 59, No. 1, July 2011
Fig. 2. (Color online) Binding energy of an in-plane benzene dimer as a function of distance d between the centers of the benzene molecules.
˚ is almost the same as rium lattice parameter a of 2.46 A that of the PBE-D2 in VASP, but the equilibrium lattice parameter c of 6.67 ˚ A is larger than that of the PBE-D2 in VASP. Especially, the structural parameters obtained by using FHI-aims are much closer to the experimental values, except for the bulk modulus B0 . The long-range Coulomb screening and many-body effects play a significant role for graphite binding. The TS vdW scheme is known to overestimate the interlayer binding in graphite [14,15], which may explain the large deviation in B0 . To determine the equilibrium lattice constants for graphite, we have investigated the potential energy surface E(a, c) with respect to the lattice constants a and c. Figures 3(a) and 3(b) show the contour plots of E(a, c) with the PBE-D2 in VASP and the TS vdW scheme in FHI-aims, respectively. We find a very long trough along the c axis, which clearly reflects weak vdW interactions between graphite layers. In summary, we performed DFT calculations with van der Waals corrections for in-plane benzene-dimer molecules and graphite. We checked the validity of vdW corrections such as the Grimme and TS vdW schemes implemented in VASP and FHI-aims, respectively. Compared to the FHI-aims code, the VASP gives a shorter equilibrium distance between the in-plane benzene dimer and a shorter equilibrium lattice parameter c of graphite by a few percent. The structural parameters obtained by using FHI-aims are closer to the experimental values. As a result, we conclude that both the Grimme and TS vdW schemes implemented in computational codes are valid for describing van der Waals corrections.
ACKNOWLEDGMENTS We would like to thank A. Tkatchenko for valuable discussions. This research was supported by the Converging Research Center Program (2010K001069) and the Priority Research Centers Program (2011-0018395)
Fig. 3. Contour plots of the potential energy surface E(a, c) with respect to the lattice constants a and c for graphite: (a) the PBE-D2 in VASP and (b) the TS vdW scheme in FHI-aims.
through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (MEST). Calculations were performed by using the Partnership & Leadership for the nationwide Supercomputing Infrastructure (PLSI) resources and Korea Research Environment Open NETwork (KREONET) of the Korea Institute of Science and Technology Information.
REFERENCES [1] S. Grimme, J. Comput. Chem. 27, 1787 (2006). [2] A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). [3] G. Kresse and J. Hafner, Phys. Rev. B 48, 13115 (1993). [4] G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994). [5] G. Kresse and J. Furthm¨ uller, Comput. Mater. Sci. 6, 15 (1996). [6] G. Kresse and J. Furthm¨ uller, Phys. Rev. B 54, 11169 (1996).
Ab initio Calculations with van der Waals Corrections: Benzene-benzene · · · – Jinwoo Park et al. [7] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun. 180, 2175 (2009). [8] V. Havu, V. Blum, P. Havu and M. Scheffler, J. Comput. Phys. 228, 8367 (2009). [9] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996); 78, 1396 (E) (1997). [10] P. E. Bl¨ ochl, Phys. Rev. B 50, 17953 (1994).
-199-
[11] G. Kresse and J. Joubert, Phys. Rev. B 59, 1758 (1999). ´ [12] T. Bu˘cko, J. Hafner, S. Leb´egue and J. G. Angy´ an, J. Phys. Chem. 114, 11814 (2010). [13] Y. X. Zhao and I. L. Spain, Phys. Rev. B 40, 993 (1989). [14] O. A. von Lilienfeld and A. Tkatchenko, J. Chem. Phys. 132, 234109 (2010). [15] F. Hanke, J. Comput. Chem. 32, 1424 (2011).
