11814
J. Phys. Chem. A 2010, 114, 11814–11824
Improved Description of the Structure of Molecular and Layered Crystals: Ab Initio DFT Calculations with van der Waals Corrections ´ ngya´n¶ Toma´sˇ Bucˇko,*,†,§ Ju¨rgen Hafner,‡ Se´bastien Lebe`gue,*,¶ and Ja´nos G. A Department of Physical and Theoretical Chemistry, Faculty of Natural Sciences, Comenius UniVersity, Mlynska´ Dolina, SK-84215 BratislaVa, SloVakia, Fakulta¨t fu¨r Physik and Center for Computational Materials Science, UniVersita¨t Wien, Sensengasse 8/12, Wien 1090 Austria, and E´quipe de Mode´lisation Quantique et Cristallographique, CRM2, Institut Jean Barriol, UMR 7036, CNRS - Nancy-UniVersite´, B.P. 239, F-54506 VandœuVre-le`s-Nancy, France ReceiVed: July 13, 2010; ReVised Manuscript ReceiVed: September 3, 2010
The implementation of technique for full structural optimizations of complex periodic systems in the DFTPAW package VASP, including the volume and shape of the unit cell and the internal coordinates of the atoms, together with a correction that allows an appropriate modeling of London dispersion forces, as given by the DFT-D2 approach of Grimme [Grimme, S. J. Comp. Chem. 2006, 27, 1787], is reported. Dispersion corrections are calculated not only for the forces acting on the atoms, but also for the stresses on the unit cell. This permits a simultaneous optimization of all degrees of freedom. Benchmark results on a series of prototype systems are presented and compared to results obtained by other methods and experimental data. It is demonstrated that the computationally inexpensive DFT-D2 scheme yields reasonable predictions for the structure, bulk moduli, and cohesive energies of weakly bonded materials.
Density functional theory (DFT) has become an important tool in investigating the structure of various types of materials.1–5 Whereas in small molecular systems it is the chemical bonding between nearest neighbors that determines the structure, with the increasing size of the system, and in particular in molecular assemblies, weak intermolecular forces start to play a moreand-more determining role. This is particularly true in crystalline molecular materials held together either by electrostatic and polarization interactions or, in their absence or in case of their weaker intensity, by London dispersion forces. To cope with this challenge, theory should have at its disposal not only appropriate total energy functionals, which take into account correctly these weak intermolecular forces, but also adapted tools to perform geometry optimizations on shallow potential energy surfaces with many local minima, which allow us to optimize not only atomic positions but also the lattice/ cell parameters in crystalline materials. Popular density functionals are unable to describe correctly van der Waals interactions resulting from dynamical correlations between fluctuating charge distributions,6 therefore they are intrinsically unadapted to optimize lattice parameters in molecular materials. In spite of the considerable progress achieved in this domain during the recent few years,7 a fully ab initio treatment of the dispersion forces is still too expensive to become a routine tool for larger systems. A pragmatic method to work around this problem has been given by the DFT-D approach,8 which consists in adding a semiempirical dispersion potential to the conventional Kohn-Sham DFT energy. Among
numerous attempts to provide a consistent parametrization compatible with the exchange-correlation functional, the approach of Grimme9–11 seems to be the most successful one. In this method the van der Waals interactions are described via a simple pairwise force field optimized for popular DFT functionals, including PBE12 and B3LYP.13,14 The performance of the method to optimize extended systems has been tested recently by other researchers, see for example, refs 15–19. After completion of our work, Grimme and his co-workers published a considerably more sophisticated parametrization, called DFTD3.20 In the present article we focus on the previous version of the method,10 referred to as DFT-D2. In this study we describe our implementation of Grimme’s DFT-D2 approach in the popular solid-state PAW code VASP,21–24 including not only forces but also stresses on the unit cell, making possible full geometry relaxations of atomic positions, cell volume, and shape. Full relaxation of molecular crystal structures has been attempted by Ferenczy et al.25,26 by semiempirical model Hamiltonians combined with an empirical force field. The general framework for full structural optimizations of periodic systems has been discussed in ref 27. Furthermore we test the performance of Grimme’s DFT-D2 method for the optimization of atomic and lattice degrees of freedom on a set of selected extended systems. In the following section the reader finds a brief description of the methodology, including working expressions for the calculation of stresses. Section 3 describes results for various classes of materials: rare gas solids, molecular crystals, layered materials and finally some covalently bound and partially ionic systems. The main conclusions are presented in Section 4.
* To whom correspondence should be addressed. E-mail: tomas.bucko@ univie.ac.at (T.B.);
[email protected] (S.L.). † Comenius University. ‡ Vienna University. ¶ Nancy-University. § Previous address: Fakulta¨t fu¨r Physik and Center for Computational Materials Science, Universita¨t Wien, Sensengasse 8/12, Wien 1090 Austria.
2. Methodologies 2.1. Electronic Structure Calculations. Periodic DFT calculations have been performed using the VASP code.21–24 The Kohn-Sham equations have been solved variationally in a plane wave basis set using the projector-augmented-wave (PAW)
1. Introduction
10.1021/jp106469x 2010 American Chemical Society Published on Web 10/05/2010
Structure of Molecular and Layered Crystals
J. Phys. Chem. A, Vol. 114, No. 43, 2010 11815
TABLE 1: Parameters for Grimme’s PBE-D2 Force Field Used in This Study (after ref 10) element
Ci6
Ri
element
Ci6
Ri
H B C N O Ne Na Si S Cl
0.14 3.13 1.75 1.23 0.70 0.63 5.71 9.23 5.57 5.07
1.001 1.485 1.452 1.397 1.452 1.243 1.144 1.716 1.683 1.639
Ar K V Se Kr Nb Mo Te Xe
4.61 10.80 10.80 12.64 12.01 24.67 24.67 31.74 29.99
1.595 1.485 1.562 1.771 1.727 1.639 1.639 1.892 1.881
method of Blo¨chl,28 as adapted by Kresse and Joubert.29 The exchange-correlation energy was described by the functional of Perdew, Burke, and Ernzerhof (PBE)12 based on the generalized gradient approximation. Computational details such as plane wave cutoff and Brillouin zone sampling are discussed for each system separately in Section 3. For all the computations presented in this paper, standard PAW data sets have been used. 2.2. Dispersion Energy Correction. van der Waals interactions have been computed using the semiempirical correction of Grimme.10 In this method, the total energy of the system is defined as a sum of the self-consistent Kohn-Sham energy (EKS-DFT) and a semiempirical correction (Edisp):
EDFT-D ) EKS-DFT + Edisp
Nat
Nat
Cij
∑ ∑ ∑′ |ri,0 -6 rj,L|6 f(|ri,0 - rj,L|) i)1 j)1
saR ) mod(
a is the R-component of the Cartesian position vector where xa,L R of atom a in the L ) (l1, l2, l3) unit cell of the solid. It follows from eq 2 that the position derivatives of the potential energy in the DFT-D approach are:
∂EDFT-D ∂sµi
(2)
∂Edisp ∂sµi
where the summations are over all atoms Nat and all translations of the unit cell L ) (l1,l2,l3), the prime indicates that i * j for L ) 0, s6 is a global scaling factor, C6ij denotes the dispersion coefficient for the atom pair ij, and ri,L is a position vector of atom i after performing L translations of the unit cell along lattice vectors. The term f(rij) is a damping function
f(r ) )
1 + e-d(r /R -1) ij
ij
(3)
whose role is to scale the force field such as to minimize contributions from interactions within typical bonding distances. Combination rules for dispersion coefficients Cij6 and vdW radii Rij are
Cij6 ) √Ci6Cj6
(4)
Nat
)
∑ ∑ ∂ |r
∂EKS-DFT ∂sµi
+
∂Edisp
(7)
∂sµi
j)1
) -s6
∑ ∑ ( |r j)1
i,0
L
′
i,0
L
∂ |ri,0 - rj,L | -r | ∂sµi
∂Edisp
′
j,L
)
i,0 j,L ij 6 d - ij e-d(|r -r |/R -1)f(|ri,0 - rj,L|) × - rj,L | R
Cij6 |ri,0 - rj,L |6
f(|ri,0 - rj,L |)
-r | ∂sµi
i,0
∂ |r
j,L
(8) The lattice-vector derivatives of the potential energy are computed using:
∂EDFT-D ) -Ω ∂hRβ
∑ σRγDFT-D(ht)γβ-1
(9)
γ
Here Ω is the volume of the unit cell, and σDFT-D is the stress tensor that can be expressed as a sum of σDFT computed at ab initio level plus a correction due to the force field σD. The latter term is computed as: N
and
D σµ,ν )-
Rij ) Ri + Rj
)
where the term (∂EKS-DFT)/(∂sµi ) is computed using the Hellmann-Feynman theorem and (∂Edisp)/(∂siµ) is the contribution from Grimme’s potential:
L
1
(6)
a
β
Nat
ij
∑ hRβ-1xβa,L , laR)
(1)
The dispersion energy for periodic systems is defined as
Edisp ) -s6
of the unit cell. In practice, however, terms corresponding to interactions over distances longer than a certain suitably chosen cutoff radius contribute only negligibly to Edisp and can be ignored. In this study, the cutoff radius was set to 40 Å. 2.3. Dispersion Energy Gradients. Most of the modern optimization algorithms use potential energy gradients to minimize the energy with respect to the geometry of the system. In this section we briefly discuss the calculation of all the quantities needed for a full structural optimization of periodic systems: total energy derivatives with respect to atomic positions and with respect to lattice vector components. Let the structure of a periodic system be characterized by three lattice vectors, arranged in the matrix h ) [a1, a2, a3] and by the 3N fractional coordinates, s ) {saR} defined as:
N
∑ ∑ ∑′ ∂ |ri,0 -disprj,L| xµi,0 ∂ |r ∂x-i,0 r
1 2Ω i)1
j)1
∂E
L
i,0
j,L
|
ν
(5)
(10)
respectively. All force field parameters were set as proposed in ref 10 for the PBE functional: s6 ) 0.75, d ) 20. Values of C6i and Ri are compiled in Table 1. In principle, the summation ΣL′ in eq 2 should run over all possible (i.e., infinitely many) translations
2.4. Structural Optimizations and Bulk Modulus Calculations. Atomic positions, lattice volume, and lattice parameters have been optimized simultaneously using the conjugate gradient method.30 Bulk moduli have been determined by fitting the energies from a series of constant volume optimizations (i.e.,
11816
J. Phys. Chem. A, Vol. 114, No. 43, 2010
Bucˇko et al. TABLE 2: Computed Equilibrium Lattice Constants, Bulk Moduli, and Cohesive Energies for Ne, Ar, Kr, and Xe fcc Crystals, Compared with Results from Other Theoretical Studies and with Experimental Data compound
method
a (Å)
B0 (GPa)
∆Ecoh (meV/atom)
expt PBE PBE-D2 GGA+vdW40 RPA41
4.46432 4.56 4.23 4.562 4.5
1.132 1 4 2.1
2038 20 58
expt PBE PBE-D2 GGA+vdW40 RPA41
5.30033 5.92 5.38 6.002 5.3
2.733
