Improving the Description of Nonmagnetic and Magnetic Molecular Crystals via the van der Waals Density Functional
aa r X i v : . [ phy s i c s . c o m p - ph ] J a n Improving the Description of Nonmagnetic and Magnetic Molecular Crystalsvia the van der Waals Density Functional
Masao Obata and Makoto Nakamura
Graduate School of Natural Science and Technology,Kanazawa University, Kakuma, Kanazawa 920-1192, Japan
Ikutaro Hamada
International Center for Materials Nanoarchitectonics (WPI-MANA),National Institute for Materials Science (NIMS), Tsukuba 305-0044, Japan
Tatsuki Oda
Graduate School of Natural Science and Technology,Kanazawa University, Kakuma, Kanazawa 920-1192, Japan andInstitute of Science and Engineering, Kanazawa University, Kakuma, Kanazawa 920-1192, Japan
We have derived and implemented a stress tensor formulation for the van der Waals densityfunctional (vdW-DF) with spin-polarization-dependent gradient correction (GC) recentlyproposed by the authors [J. Phys. Soc. Jpn. , 093701 (2013)] and applied it to non-magnetic and magnetic molecular crystals under ambient condition. We found that the cellparameters of the molecular crystals obtained with vdW-DF show an overall improvementcompared with those obtained using local density and generalized gradient approximations.In particular, the original vdW-DF with GC gives the equilibrium structural parameters ofsolid oxygen in the α -phase, which are in good agreement with the experiment. I. INTRODUCTION
The oxygen molecule has a spin-triplet ground state, and the competing magnetic and molecularinteractions result in a rich variety of structural, electronic, and magnetic phases under widetemperature and pressure ranges. Much effort has been devoted to clarifying the magnetism andatomic structure of solid oxygen.
The magnetism of solid oxygen is particularly interesting, asthe magnetic moment is sensitive to the relative distance between molecules owing the short-rangefeature of the direct magnetic interaction. Thus, accurate prediction of the crystal structure isone of the key issues for understanding the magnetism in solid oxygen. Theoretical works on solidoxygen in the literature are mostly for those under a high pressure, presumably because undera relatively low pressure, the vdW interaction plays an important role, where the conventionalsemilocal density functional theory (DFT) fails to describe it properly.An accurate description of the vdW interaction is, however, challenging when using the elec-tronic structure method, particularly, DFT within the semilocal approximation, which is knownto describe the dispersion interaction poorly. In quantum chemistry, there are well-establishedmethods, such as the second-order Møller-Plesset perturbation theory (MP2), coupled-cluster cal-culations with single, double, and perturbative triple excitations [CCSD(T)], and the symmetry-adopted perturbation theory (SAPT), which can describe the vdW interaction. The quantumMonte-Carlo (QMC) method and the adiabatic connection fluctuation dissipation theorem withinthe random phase approximation (ACFDT-RPA) are also known to describe the vdW interactionaccurately. However, because these highly accurate methods are computationally very demandingand applicable only to small systems, practical methods based on DFT are desirable for the simula-tions of complex systems. In this sense, semi-empirical dispersion-corrected DFT is widely used,and recent developments allow more accurate calculations of the vdW interaction with analmost negligible overhead. A more elaborate approach based on the maximally localized Wannierfunction was also proposed to correct for the missing dispersion interaction in DFT. Among other approaches, the van der Waals density functional (vdW-DF), which was developedby Rydberg et al. and Dion al. , has attracted much attention, because it does not requireempirical parameters and depends only on charge density and its gradient, thereby allowing one toperform large-scale calculation at a modest computational cost. There have been applications ofvdW-DF to several systems, and to extend the applicability of the methods, there have been severalproposals for more accurate vdW-DFs. However, because the vdW interaction is not directlyrelated to spin polarization, a truly spin-polarized version of vdW-DF has not yet been developed,except for the one proposed by Vydrov and Van Voorhis. Thus, the application of vdW-DFs tomagnetic systems is quite limited. In our previous work, we proposed a practical approach to amagnetic system within the framework of vdW-DF, in which spin-polarization-dependent gradientcorrection is added to the local correlation energy and potential, instead of developing the spin-polarized version of the nonlocal correlation. The method has been applied to an oxygen dimer,and shown to be applicable to magnetic systems.In this work, we have studied solid oxygen at ambient pressure in the α phase ( α -O ) usingour vdW-DF for magnetic materials. We implemented the stress tensor arising from the nonlo-cal correlation to optimize the cell parameters of α -O and used several vdW-DFs with differentexchange and nonlocal correlation. We found that one of the vdW-DFs shows a significant im-provement on the structural properties over the conventional local density approximation (LDA)or generalized gradient approximation (GGA). The remaining part of this paper is organized asfollows. In Sect. II, we describe our implementation of vdW-DF and the computational detail. Totest our implementation of vdW-DF, we first applied our code to solid argon, graphite, selenium,and dry ice (solid carbon dioxide), in which vdW forces play an important role in the binding, andthe results are given in Sects. III A–III D. The results for solid oxygen are discussed in detail inSect. III E. The summary for this paper is given in Sect. IV. II. METHODA. van der Waals density functional
The exchange-correlation energy in vdW-DF is given by E xc = E x + E locc + E nlc , (1)where E x , E locc , and E nlc are the exchange energy, short-range local correlation energy, and nonlocalcorrelation energy, respectively. E nlc is expressed as a double integral of the nonlocal interactionkernel φ ( q , q , r ) over the spatial coordinates r and r as E nlc = 12 Z Z n ( r ) φ ( q , q , r ) n ( r ) d r d r , (2)where r = | r − r | , q = q ( r ), q = q ( r ), and q ( r ) is a function of the charge density n ( r ) andits gradient |∇ n ( r ) | . Rom´an-P´erez and Soler (RPS) proposed an efficient method for evaluating E nlc , as the direct computation of the double integral in Eq. (2) is extremely time-consuming. Inthe RPS scheme, the vdW kernel φ is expanded in terms of the interpolating function p α ( q ), whichsatisfies p α ( q β ) = δ αβ , so that φ is a fixed value at a given q β point, allowing the use of the fastFourier transform (FFT) in the evaluation of E nlc . However, φ has a divergence when q α , q β → φ with the softform, and an LDA-like approximation was employed for E nlc near the origin. Alternatively, Wuand Gygi (WG) proposed a simplified implementation to avoid the divergence in φ as follows.The nonlocal interaction kernel, multiplied by q and q , is expanded as q q φ ( q , q , r ) ≃ X αβ q α p α ( q ) q β p β ( q ) φ αβ ( r ) , (3)so that the divergence at q α , q β → η α ( r ) = q α n ( r ) p α ( q ( r )) /q ( r ), the nonlocal correlation is calculated in the reciprocal space as E nlc = 12 X αβ Z Z η α ( r ) η β ( r ) φ αβ ( r ) d r d r = Ω2 X G X αβ η ∗ α ( G ) η β ( G ) φ αβ ( G ) , (4)where G is the reciprocal lattice vector G = | G | , η α ( G ) and φ αβ ( G ) are the Fourier componentsof η α ( r ) and φ αβ ( r ), respectively, and Ω is the volume of a unit cell. The Fourier component φ αβ ( G ) is calculated as φ αβ ( G ) = 1 q α F αβ (cid:18) Gq α (cid:19) , (5)where F αβ ( k ) = 4 πk Z ∞ uφ (cid:18) u, q β q α u (cid:19) sin( ku ) du. (6)Note that the G = 0 term in Eq. (4) vanishes when α = β . We discuss the behaviors of the φ function in Appendix B.The nonlocal correlation potential within the White-Bird algorithm is given by v nlc ( r i ) = δE nlc δn ( r i ) ≃ N Ω dE nlc dn ( r i )= X α u α ( r i ) ∂η α ( r i ) ∂n ( r i ) + X j u α ( r j ) ∂η α ( r j ) ∂ ∇ n ( r j ) ∂ ∇ n ( r j ) ∂n ( r i ) , (7)where u α ( r i ) = Ω N X βj φ αβ ( r ij ) η β ( r j ) = X β G φ αβ ( G ) η β ( G ) e i G · r i , (8)with N being the number of FFT grids. B. Pressure tensor
To optimize the lattice parameter, it is necessary to calculate the pressure tensor, and there isan additional contribution from the nonlocal correlation in vdW-DF.
In our formulation, anadditional contribution to the pressure tensor Π kℓ , ( k, ℓ = x, y, z ) is calculated as(Π nlc ) kℓ = − X m ∂E nlc ∂h km ( h t ) mℓ (9)= − E nlc δ kℓ − N X m X j X α u α ( r j ) ∂η α ( r j ) ∂n ( r j ) ∂n ( r j ) ∂h km ( h t ) mℓ + 1 N X j X α u α ( r j ) ∂η α ( r j ) ∂ |∇ n ( r j ) | ( ∇ n ( r j )) k ( ∇ n ( r j )) ℓ |∇ n ( r j ) |− N X m X j X α u α ( r j ) ∂η α ( r j ) ∂ |∇ n ( r j ) |× ∇ n ( r j ) |∇ n ( r j ) | · X G ∂n ( G ) ∂h km ( h t ) mℓ ( i G ) e i G · r j + 12 X αβ X G η ∗ α ( G ) η β ( G ) ∂φ αβ ( G ) ∂ | G | G k G ℓ | G | , (10)where the matrix h = ( a , a , a ) is a set of primitive lattice vectors a k ( k = 1 , ,
3) and h t is theinverse of h . C. Spin-polarization-dependent gradient correction of the local correlation
In the original vdW-DF, LDA was employed for the short-range local correlation energy, to avoidpossible double counting of the contribution from |∇ n | contained in E nlc . In a spin-polarized sys-tem, spin-polarization-dependent gradient correction should be included in the correlation energyand potential, but such a contribution is missing in the nonlocal correlation, as it is mostly for-mulated for the spin-unpolarized system. In our previous work , we propose a simple scheme toinclude such spin-polarization-dependent gradient contribution as follows. We define the local partof the correlation functional as E locc = E LSDAc [ n ↑ , n ↓ ] + ∆ E c [ n, ζ ] , (11)where ∆ E c [ n, ζ ] = E GGAc [ n ↑ , n ↓ ] − E GGAc [ n/ , n/ Z n { H ( r s , ζ, t ) − H ( r s , , t ) } d r , (12) H , r s , ζ , and t are the gradient contribution, Seitz radius ( n = 3 / πr s ), relative spin polarization,and dimensionless density gradient proportional to |∇ n | , respectively. We use the functional formfor H proposed by Perdew, Burke, and Ernzerhof (PBE). The present functional is reduced tothe original one in the absence of spin polarization.
D. Technicalities
We have implemented vdW-DF in our in-house DFT code, which uses a plane-wave basis setand ultrasoft pseudopotentials.
In the construction of the pseudopotentials, we neglected thenonlocal correlation ( E nlc ) and employed the semilocal exchange and correlation functionals, whichare consistent with the solid-state calculation. For instance, the pseudopotential constructed usingthe refit Perdew Wang (PW86R) exchange and the LDA correlation was used in vdW-DF2calculation. The use of the pseudopotential generated with the semilocal exchange correlationfunctional has been justified. We also performed the calculations using the pseudopotentialgenerated with the PBE functional and found that there is an insignificant difference in bindingenergy ( ∼ k -points. In this work, we employed the Perdew-Zunger exchange correlation for LDA and the PBE functional for GGA .For the vdW-DF calculations, we employed the following functionals: the original vdW-DF withthe revPBE exchange and LDA correlation, the second version of the vdW-DF (vdW-DF2), whichuses PW86R exchange, vdW-DF paired with Cooper’s (C09) exchange (vdW-DF C09x ), andvdW-DF2 paired with C09 exchange (vdW-DF2
C09x ). We also implemented the revised Vydrov-Van Voorhis functional (rVV10) with the WG implementation (see Appendix A), which usesthe PW86R exchange, PBE correlation, and the nonlocal correlation. The functional is designed TABLE I. Optimized lattice parameter ( a fcc ), equilibrium volume ( V ), and cohesive energy (∆ E coh ) forsolid argon in the fcc structure, obtained using different density functionals. The CCSD(T), ACFDT-RPAusing PBE orbitals, and experimental values are also listed for comparison.Functional a fcc V ∆ E coh (˚A) (˚A /atom) (meV/atom)LDA 4.97 30.61 135PBE 6.04 55.03 14vdW-DF 5.53 42.34 142vdW-DF2 5.29 37.05 116vdW-DF C09x
C09x a b c d d ea Ref. 45. b Ref. 46. c Ref. 47. d Ref. 48. e Ref. 49. to give accurate C coefficients and to vanish in the uniform electron gas limit. In the case ofantiferromagnetic α -O , we employed the spin-polarization-dependent gradient correction (GC) asdescribed previously, and “-GC” is appended to the functional name when GC is considered .Note that GC is not necessary for rVV10, because the PBE-GGA correlation is used for the localcorrelation.To evaluate φ in the vicinity of q α = 0, we use a linear mesh up to q = q m with N lin grid points,and the logarithmic mesh is used from q m to q c with N log grid points, where q c is the cutoff for the q -mesh. We use q c = 8 a . u ., q m = 0 .
01 a . u ., N lin = 3, and N log = 28, except for rVV10, in which q c = 3 a . u . , and q m = 0 .
005 a . u . are used. We confirmed that the cutoff q c is sufficiently large,by monitoring the distribution of q and the nonlocal correlation energy as a function of q (seeAppendix C). TABLE II. Optimized lattice parameters ( a hcp and c hcp ), equilibrium volume ( V ), and cohesive energy(∆ E con ) for hcp argon. Functional a hcp c hcp /a hcp V ∆ E coh (˚A) (˚A /atom) (meV/atom)LDA 3.51 1.633 30.45 135PBE 4.16 1.633 50.92 13vdW-DF 3.89 1.633 41.63 142vdW-DF2 3.74 1.636 37.08 116vdW-DF C09x
C09x a a Ref. 50.
III. RESULTSA. Fcc and hcp argons
We begin with solid argon in the face-centered-cubic (fcc) and hexagonal-close-packed (hcp)structures, as a typical vdW bonded solid. It is known that at ambient pressure, argon in the fccstructure is the most stable, whereas that in the hcp structure is metastable. For both phases, weoptimized lattice parameters using several vdW-DFs as well as LDA and GGA. The energy cutoffsof 60 and 500 Ry were used for wave function and electron density, respectively. A 14 × × × ×
8) MP special k -point set was used for the fcc (hcp) phase. The optimized lattice constantand cohesive energy for fcc argon are summarized in Table I, along with CCSD(T) and experimentalvalues. The lattice constant obtained using LDA (PBE) is underestimated (overestimated),whereasthose with vdW-DFs are improved except for vdW-DF2 C09x , and most of the vdW-DFs improvecohesive energy. Our results for fcc argon are in good agreement with the recent results by Tranand Hutter. The results for hcp argon are summarized in Table II. The difference between the cohesiveenergies for the fcc and hcp phases is very small (at most 2 meV/atom) and these phases are almostdegenerated. The lattice parameters for the hcp structure satisfy the ideal ratio a fcc ∼ √ a hcp , c hcp /a hcp = 1 . a fcc , a hcp , and c hcp are the lattice constant for the fcc structure, thein-plane lattice constant for the hcp structure, and the out-of-plane-lattice constant for the hcpstructure, respectively. The deviations of a fcc from the ideal value are 0.15 ˚A for PBE, 0.36˚A for vdW-DF2 C09x , and 0.03 ˚A for other functionals. Our result indicates that these two phasesare almost identical and cannot reproduce the experimentally observed stable fcc phase at ambientpressure. This apparent discrepancy may be solved by introducing an entropic effect from phonons.Indeed, Ishikawa et al . performed lattice dynamics calculations and successfully predicted a stablefcc argon by taking into account the entropic contributions. B. Graphite
Graphite consists of the two-dimensional carbon allotrope, graphene, and graphene layers boundin the out-of-plane direction with a vdW interaction. Graphite has been studied extensively bothexperimentally and theoretically, and is often used to assess the accuracy of a new method for weakinteractions, as benchmark calculations with a highly accurate electronic structure method such asQMC and ACFDT-RPA. In our calculation, the Brillouin zone was sampled using an 8 × × k -point set, and planewave cutoffs of 80 and 480 Ry were used for wave functions and electron density, respectively. Wefully optimized the lattice parameters according to the calculated internal pressure, and bindingenergy was determined from the difference in total energy at the equilibrium and that at aninterlayer distance larger than 13 ˚A. In Table III, optimized lattice parameters and binding energiesobtained using different exchange-correlation functionals are summarized.In accordance with the literature, LDA gives lattice parameters, which are in good agreementwith the experimental results, whereas GGA yields a much larger out-of-plane lattice parameter( c ) and negligible binding energy, suggesting that the latter cannot predict the binding of graphite.The binding energy obtained with LDA is smaller than that obtained from the experiment and bya highly accurate theoretical method. In general, our vdW-DF results are in good agreement withthose of previous studies, that employed different vdW-DFs. vdW-DF and vdW-DF2predict the binding energy (∆ E b ), in good agreement with the experimental results. They also im-prove the description of structural properties, but the equilibrium volumes ( V ) are overestimated. V is improved using vdW-DF C09x and rVV10, but they overestimate ∆ E b . We note that our ∆ E b obtained using rVV10 is almost twice as large as that reported by the original authors. At present,the origin of the discrepancy is yet to be clarified, but similar larger values ( ∼
70 meV/atom) wereobtained using different implementations of rVV10 and different potentials. On the otherhand, vdW-DF2
C09x predicts that lattice parameters and ∆ E b are in good agreement with the0 TABLE III. In-plane lattice constant ( a ), out-of-plane lattice constant ( c ), equilibrium volume ( V ), and thebinding energy ∆ E b for graphite. The numbers in parentheses were obtained using the experimental value.Functional a c V ∆ E b (˚A) (˚A) (˚A /atom) (meV/atom)LDA 2.45 6.69 8.71 24PBE 2.47 8.76 11.61 2vdW-DF 2.49 7.18 9.61 55vdW-DF2 2.48 7.05 9.41 53vdW-DF C09x
C09x a b (2.46) 6.777 (8.88) 71vdW-DF-cx c d (2.46) 6.852 (8.98) 56 ± e (2.46) 6.68 (8.75) 48Expt. 2.46 f f f ± ga Ref. 26. b Ref. 54. c Ref. 27. d Ref. 52. e Ref. 53. f Ref. 55. g Ref. 56.FIG. 1. (Color online) Trigonal structure of selenium. experimental results, in line with a previous study. TABLE IV. Optimized lattice parameters ( a and c ), equilibrium volume ( V ), internal parameter ( x ), shortestSe-Se distance in a chain ( ℓ ), shortest Se-Se distance in different chains ( ℓ ), and binding energy (∆ E b ) fortrigonal selenium.Functional a c V x ℓ ℓ ∆ E b (˚A) (˚A) (˚A /atom) (˚A) (˚A) (meV/atom)LDA 3.89 5.04 21.94 0.2540 2.40 3.06 395PBE 4.51 4.97 29.16 0.2148 2.36 3.58 55vdW-DF 4.70 5.04 32.15 0.2091 2.39 3.74 189vdW-DF2 4.57 5.12 30.81 0.2177 2.42 3.62 228vdW-DF C09x
C09x a a a b a aa Ref. 64. b Ref. 65.
C. Trigonal selenium
Trigonal selenium is formed by a bundle of one-dimensional covalently bonded atomic chiralchains. As already shown in previous works, the gradient correction to LDA improves the de-scription of selenium, but is still less satisfactory, presumably because of the lack of vdW interactionbetween the chiral chains in GGA. Buˇcko et al . demonstrated using the semi-empirical dispersioncorrection that the structural parameters are significantly improved, and the agreement with theexperimental results becomes much better. Here, we use vdW-DF to address the importance ofvdW interaction.The trigonal selenium belongs either to the P
21 or P
21 space group and has three atoms perunit cell (see Fig. 1). The atomic position in a cell can be described using the internal parameter x , which scales the distance from the screw axis to the atomic position. The structure of thetrigonal selenium is also characterized by the shortest Se-Se distance in the same chain ( ℓ ) andthe shortest Se-Se distance in different chains ( ℓ ). We used an 8 × × k -point set andplane wave cutoffs of 60 and 500 Ry for wave functions and electron density, respectively. Bindingenergy was calculated from the difference between the total energies of a solid in equilibrium andthe isolated chain.Calculated structural parameters and binding energy are summarized in Table IV. In general,2 FIG. 2. (Color online) Cubic
P a ). the lattice constant along the chiral chain ( c ) is overestimated by all the functionals used in thisstudy, whereas the accuracy for the lattice constant a varies. LDA significantly underestimates a ,whereas PBE significantly improves the description of a , in good agreement with previous studies.Both vdW-DF and vdW-DF2 overestimate a and c , whereas the use of the C09 exchange (vdW-DF C09x and vdW-DF2
C09x ) leads to the underestimation of a . On the other hand, although c isslightly overestimated, rVV10 provides a balanced description of both the lattice constants. It isfound from our results that there is a correlation between calculated a and c : if a is overestimated,the error in c tends to be smaller. Thus, to obtain accurate structural parameters for the trigonalselenium, it is very important to describe both covalent and vdW interactions accurately, suggestingthat this material can be a good benchmark system for a vdW inclusive functional. D. Dry ice
The solid form of carbon dioxide is called dry ice. Carbon dioxide has no dipole moment, andthus attractive vdW forces play an important role in the condensation. We chose this material asa representative application of our vdW-DF implementation to non magnetic molecular crystals.Crystalline dry ice has a cubic symmetry with the space group of
P a × × k -point set was used for Brillouin zone sampling.Optimized structural parameters and binding energy for the dry ice are given in Table V,along with the MP2 and experimental (150K) results. Our PBE equilibrium volume is in goodagreement with the theoretical value obtained by Bonev et al. (52.9 ˚A / CO ) using the samefunctional. On the other hand, both the equilibrium volume and binding energy obtained usingvdW-DFs are in better agreement with the highly accurate MP2 and experimental (150K) results,3 TABLE V. Equilibrium lattice parameter ( a ), volume ( V ), C-O bond length ( ℓ b ), and binding energy (∆ E b )for dry ice. Functional a V ℓ b ∆ E b (˚A) (˚A / CO ) (˚A) (meV/CO )LDA 5.28 36.88 1.165 370PBE 6.04 55.11 1.175 103vdW-DF 5.77 47.94 1.180 364vdW-DF2 5.61 44.04 1.178 346vdW-DF C09x
C09x a b b b ca Ref. 69. b Ref. 70. c Ref. 71. suggesting the improvement over LDA and PBE. The maximum deviation of the equilibrium volumeobtained using vdW-DF is 9.1% compared with that obtained in the low-temperature experiment.The above results suggest that PBE overestimates the equilibrium volume because of the lack ofdispersion forces, and a more accurate description of the structure and energetics of dry ice is madepossible by considering the dispersion forces with vdW-DF. Regarding the severe underestimationof the binding energy using vdW-DF2
C09x , it has been found from a systematic assessment usingthe S22 dataset that although it predicts reasonable intermolecular separation (distance), thefunctional severely underestimates the binding energy, because the C09 exchange is too repulsiveat a relatively large density gradient relevant to the intermolecular region. This result implies thatvdW-DF2 C09x inaccurately describes the intermolecular vdW interaction.
E. Solid oxygen in the α phase α -O has an antiferromagnetic ground state with the crystal structure of the C /m space group,as shown in Fig. 3. There are four lattice parameters, a , b , c , and β , and the internal parameters of ℓ b (bond length of O molecule) and θ (tilted angle of molecular axis). The molecular axis is almostperpendicular to the ab -plane and slightly tilted within the ac -plane. The angle θ was reported tobe ∼ ◦ in the experiment. . In the calculation, we used plane wave cutoffs of 160 and 960 Ry for4 FIG. 3. (Color online) Crystal structure of solid oxygen in α phase. The arrows on molecules represent themagnetic moment.TABLE VI. Optimized lattice parameters ( a , b , c , and β ), nearest-neighbor distance ( ℓ mol ≡ √ a + b / V ), bond length ( ℓ b ), magnetic moment ( M a ) on oxygen atom, and binding energy ofmolecule (∆ E b ). Experimental values are shown for comparison.Functional a b c β ℓ mol V ℓ b M a ∆ E b (˚A) (˚A) (˚A) (deg) (˚A) (˚A ) (˚A) ( µ B ) (meV)LDA 3.29 3.28 4.05 113.9 2.32 19.93 1.202 0.30 552PBE 4.59 3.93 5.05 122.1 3.02 38.54 1.218 0.66 41vdW-DF 4.68 3.68 4.70 125.2 2.98 33.05 1.231 0.66 221vdW-DF-GC 4.94 3.57 4.91 128.4 3.05 33.85 1.231 0.66 213vdW-DF2 3.83 3.85 4.29 118.6 2.72 27.75 1.235 0.60 225vdW-DF2-GC 3.91 3.88 4.30 119.0 2.76 28.56 1.235 0.62 209vdW-DF C09x
C09x -GC 3.59 3.47 4.17 115.4 2.50 23.51 1.215 0.50 255vdW-DF2
C09x
C09x -GC 3.79 3.62 4.36 115.8 2.62 26.92 1.217 0.57 71rVV10 3.71 3.62 4.18 117.2 2.59 24.94 1.225 0.56 240Expt. a a Ref. 4. wave functions and electron density, respectively, and an 8 × × k -point set was used for theBrillouin zone sampling. During the structural optimization, the molecular axis was fixed in thedirection perpendicular to the ab -plane, but the residual forces on the atom were small, typicallyless than 0.001 eV/˚A.In Table VI, we report the structural parameters and binding energies obtained using different5 TABLE VII. Magnetic energy for α -O (∆ E mag ).Functional ∆ E mag (meV)PBE 95vdW-DF 118vdW-DF-GC 87vdW-DF2 304vdW-DF2-GC 251 functionals. As expected, LDA severely underestimates the equilibrium volume by 42.8%, whereasPBE overestimates it, but the error is marginal (10.6%). However, this unexpectedly small erroris due to the cancellation of the errors in a ( − b (14.6%). The error in c is surpris-ingly small, but the trend in the lattice parameters estimated using PBE is not systematic. Onthe other hand, all the vdW-DFs underestimate a and c , whereas b is overestimated, and as aresult, the equilibrium volumes are consistently underestimated. The angle β is also consistentlyunderestimated. Among vdW-DFs, the original vdW-DF proposed by Dion et al . predicts themost accurate structural parameters, and inclusion of GC further improves the structural parame-ters, which are in good agreement with the experimental results, suggesting the importance of thespin-polarization-dependent gradient correction of the local correlation.The nearest-neighbor distance between O molecules ( ℓ mol ≡ √ a + b /
2) is a very impor-tant quantity in α -O , as it strongly correlates with the antiferromagnetic interaction between O molecules. We found that ℓ mol is underestimated by LDA and PBE, but is increased using vdW-DF-GC to a distance of 3.05 ˚A which is close to the experimental value (3.20 ˚A). We also note thatthere is a correlation between β and c estimated using vdW-DF, i.e., the smaller the c , the smallerthe β . As discussed in previous works, the distance between magnetic molecules is determinedby a subtle balance between magnetic and vdW interactions. In the solid state, not only the bal-ance between magnetic and vdW interactions within the ab -plane, but also that between ab -planes(out-of-plane direction) plays a decisive role: a complicated interplay among the antiferromagnetic,ferromagnetic, and vdW interactions in three dimensions determines the structure of α -O . Themolecules in the ab -plane can be considered to form a triangular configuration. In the experimentalstructure, this configuration is very similar to an equilateral triangle, where the antiferromagneticinteraction between the neighboring molecules could destabilize the magnetic structure or distort6the triangle to reduce the magnetic frustration. On the basis of this consideration and the knowl-edge acquired in a previous study , the underestimation of the nearest-neighbor distance ( ℓ mol )using vdW-DF may be ascribed to the overestimation of the antiferromagnetic interaction betweenmolecules with respect to the ferromagnetic one. The structural properties depend strongly on theexchange energy functional used, as shown in Table VI. Such a behavior has also been shown inthe previous study on the pair of oxygen molecules . The revPBE exchange functional in the orig-inal vdW-DF gives a more repulsive intermolecular potential than the other exchange functionals.The other vdW-DFs consisting of different exchange functionals (C09 and PW86R) underestimate ℓ mol , indicating that these exchange functionals enhance the overstabilization of the antiferromag-netic state. Thus, to improve the description of α -O , it is necessary to develop more accurateexchange and correlation functionals that predict a balanced description of antiferromagnetic andferromagnetic interactions.The magnetic interaction J is proportional to − t / ∆ E , where t and ∆ E are the transfer integraland the energy gap between the spin-up and spin-down states immediately above and below theFermi level, respectively. Because GGA is known to underestimate the energy gap for insulatingand molecular systems, the antiferromagnetic states are overstabilized as a result of the underesti-mation of ∆ E (overestimation of | t | as well). To capture a trend in the functionals on the magneticinteraction, we calculated magnetic energy, defined as the energy difference between the ferromag-netic (F) and antiferromagnetic(AF) states at the same AF crystal structure (∆ E mag = E F − E AF ).The results obtained using PBE, vdW-DF, and vdW-DF2 (with and without GC) are summarizedin Table VII. It is found that vdW-DF, which predicts more accurate lattice parameters (larger ℓ mol ) than vdW-DF2, yields a smaller ∆ E mag , suggesting that the intermolecular distance playsan important role in the magnetic interaction. Inclusion of GC leads to the enlargement of thenearest O distance ( ℓ mol ) and reduction of ∆ E mag , because GC destabilizes the antiferromagneticstates, which works as a repulsive force between the nearest O molecules . IV. SUMMARY
In this work, we have studied solid oxygen in the α -phase with the antiferromagnetic configura-tion, using several variants of vdW-DF and the correction scheme for magnetic systems proposedin our previous study. To test the vdW-DF energy, force, and stress implemented in our programcode, we performed the calculations for solid argon in the fcc and hcp phases, graphite, trigonalselenium, and dry ice, in which the vdW interaction is considered to be important. We have found7that vdW-DF shows an overall improvement in the structure and energy of these materials overLDA and GGA, in good agreement with previous studies. In the case of antiferromagnetic solidoxygen, we have found that vdW-DFs consistently underestimate the lattice parameters, and theoriginal vdW-DF proposed by Dion et al . with our spin-polarization-dependent gradient correctionscheme (vdW-DF-GC) provides the most accurate structural parameters among other functionals.We point out the competing nature of ferromagnetic, antiferromagnetic, and vdW interactions in asolid oxygen, and for a more accurate prediction of the structural properties of solid oxygen, a bal-anced description of these interactions is decisively important. Thus, we suggest that solid oxygencan be a critical test material for an electronic structure method, which can account for magneticand vdW interactions accurately. Applications of the present vdW-DF scheme for the magneticsystems to the molecular magnet or metalorganic molecular systems are highly anticipated, whichpromote research on molecular spintronics. ACKNOWLEDGMENTS
The computation in this work was performed using the facilities of the Supercomputer Center,Institute for Solid State Physics, University of Tokyo, and the facilities of the Research Cen-ter for Computational Science, National Institutes of Natural Sciences, Okazaki, Japan. Thiswork was partly supported by Grants-in-Aid for Scientific Research from JSPS/MEXT (GrantNos. 22104012 and 26400312), the Strategic Programs for Innovative Research (SPIRE), MEXT,Japan, the Computational Materials Science Initiative (CMSI), Japan, and the World PremierInternational Research Center Initiative (WPI) for Materials Nanoarchitectonics, MEXT, Japan.
Appendix A: rVV10 with the Wu-Gygi implementation
The rVV10 functional has been proposed for implementing the VV10 functional with theefficient RPS algorithm. We note that, very recently, the original VV10 has been implementedusing the WG scheme by Corsetti et al . Here, we briefly describe our implementation of rVV10with the algorithm proposed by WG. The rVV10 nonlocal correlation functional is given by E nlc = 12 Z Z n ( r ) k / ( r ) φ rVV10 (˜ q , ˜ q , r ) n ( r ) k / ( r ) d r d r , (A1)where φ rVV10 (˜ q , ˜ q , r ) is the rVV10 nonlocal correlation kernel, ˜ q = ˜ q ( r ), ˜ q = ˜ q ( r ), and k ( n ( r )) = 3 πb ( n/ π ) / /
2. ˜ q i is defined as ˜ q ( r i ) = ω ( n ( r i ) , |∇ n ( r i ) | ) /k ( n ( r i )) with ω defined inRef. 23. The parameter b is an adjustable parameter, which is determined by fitting on a training8 f ab ( G ) ( e V Å ) G(Å −1 )( q a , q b )=(0.053, 0.053)( q a , q b )=(0.579, 0.579)( q a , q b )=(2.425, 2.425)−400 0 400 800 1200 1600 0 0.2 0.4 0.6 0.8 1 f ab ( G ) ( e V Å ) G(Å −1 )( q a , q b )=(2.425, 0.138)( q a , q b )=(2.425, 1.504)( q a , q b )=(2.425, 2.425) FIG. 4. (Color online) Fourier component of the kernel function φ αβ ( G ) for several sets of q α and q β .Diagonal (upper panel) and off-diagonal (lower panel) components are shown. set (see Refs. 23 and 26 for details). By expanding the nonlocal correlation kernel in the samemanner as Eq. (3), the nonlocal correlation energy is written as E nlc = Ω2 X G X αβ ˜ η ∗ α ( G )˜ η β ( G ) φ rVV10 αβ ( G ) , (A2)where ˜ η α ( G ) is the Fourier component of ˜ η α ( r ), which is defined by ˜ η α ( r ) = q α n ( r ) p α (˜ q ( r )) /k / ( r )˜ q ( r ),and the Fourier transform of the nonlocal correlation kernel φ rVV10 is given by φ rVV10 αβ ( G ) = 1 q α / F rVV10 αβ (cid:18) G √ q α (cid:19) , (A3) F rVV10 αβ ( k ) = 4 πk Z ∞ u sin( ku ) φ rVV10 (cid:18) u, r q β q α u (cid:19) du. (A4)The nonlocal correlation potential is calculated similarly to Eq. (7). Appendix B: Fourier component of the kernel function φ αβ ( G ) In the present implementation, the Fourier components of the kernel function φ αβ ( G ) are cal-culated from the tabulated data for φ ( d , d ) ( d φ ( d , d ) in the case of WG implementation) usingEq. (5). The typical functions of φ αβ ( G ) are presented in Fig. 4. For the small q α ’s, φ αα ( G )deviates slightly from zero at G = 0 in the present computation, which may cause a numericalerror. However, the error in the absolute value of E nlc associated with the deviation is typically ∼ a = 5.29Å a =10.58Å D ( q a ) D ( q a ) e cnl e cnl a = 5.29Å a =10.58Å D i s t r i bu t i on o f q N on l o c a l c o rr e l a t i on ene r g y ( e V / a t o m ) q a D ( q a ) D ( q a ) e cnl e cnl FIG. 5. (Color online) Distributions of q function (solid curves) and ε nlc ( α ) (dotted curves) for fcc argoncalculated using vdW-DF (upper panel) and vdW-DF2 (lower panel). The red open and blue filled symbolsare for the lattice parameters ( a ) of 5.29 and 10.58 ˚A, respectively. Appendix C: Distribution of q function To expand Eq. (3), an appropriate cutoff q c and mesh for the q function ( q α ) are needed. Todetermine q c and q α , we investigated the distribution of the q function by sampling the valueson each FFT grid point. The distribution D ( q ) was calculated as a histogram by counting thenumber of samples while satisfying q α ≤ q ( r ) < q α +1 . In Fig. 5, D ( q )’s for the vdW-DF andvdW-DF2 are plotted for different lattice parameters. Furthermore, we analyzed the contributionof the nonlocal correlation energy from each q α value, by defining E nlc = X α ε nlc ( α ) ,ε nlc ( α ) = Ω2 X G | η α ( G ) | φ αα ( G ) + 2 X β ( <α ) Re( η ∗ α ( G ) η β ( G )) φ αβ ( G ) , (C1)where ε nlc ( α ) is a positive value. ε nlc ( α )’s for fcc argon calculated using vdW-DF and vdW-DF2are shown in Fig. 5. There is almost no contribution of the nonlocal correlation energy at thesmallest q α , although the samples of D ( q α ) exist there. Note that the samples at the smallest q α come from diluted electron densities. In the case of a = 5 .
29 ˚A, ε nlc ( α )’s have a peak at q = 2 . q α increases. In the case of a = 10 .
58 ˚A, the shrunken electron cloudaround atomic cores causes the distributions of D ( q α ) at larger q α ’s( > E nlc term. This is presumably because the vdW interaction is both0canceled out in the atom and fully damped at the atomic distances. From the above result, it isconfirmed that the q α mesh defined in the present work covers the energy range in the system. Y. A. Freiman and H. J. Jodl, Physics Reports , 1 (2004). G. C. DeFotis, Phys. Rev. B , 4714 (1981). S. Klotz, Th. Str¨assle, A. L. Cornelius, J. Philippe, and Th. Hansen, Phys. Rev. Lett. , 115501 (2010). R. J. Meier and R. B. Helmholdt, Phys. Rev. B , 1387 (1984). R. D. Etters, K. Kobashi, and J. Belak, Phys. Rev. B , 4097 (1985). M. Otani, K. Yamaguchi, H. Miyagi, and N. Suzuki, J. Phys.: Condens. Matter , 11603 (1998). K. Kususe, Y. Hori, S. Suzuki, and K. Nakao, J. Phys. Soc. Jpn. , 2692 (1999). K. Nozawa, N. Shima, and K. Makoshi, J. Phys. Soc. Jpn. , 377 (2002). K. Nozawa, N. Shima, and K. Makoshi, J. Phys.: Condens. Matter , 335219 (2008). M. Santoro, F. A. Gorelli, L. Ulivi, R. Bini, and H. J. Jodl, Phys. Rev. B , 064428 (2001). S. Grimme, J. Comput. Chem. , 1787 (2006). S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys. , 154104 (2010). A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. , 073005 (2009). A. Tkatchenko, R. A. DiStasio, Jr., R. Car, and M. Scheffler, Phys. Rev. Lett. , 236402 (2012). P. L. Silvestrelli, Phys. Rev. Lett. , 053002 (2008). P. L. Silvestrelli, J. Chem. Phys. , 054106 (2013). H. Rydberg, M. Dion, N. Jacobson, E. Schr¨oder, P. Hyldgaard, S. I. Simak, D. C. Langreth, andB. I. Lundqvist, Phys. Rev. Lett. , 126402 (2003). M. Dion, H. Rydberg, E. Schr¨oder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. , 246401(2004) [Erratum , 109902 (2005)]. O. A. Vydrov and T. Van Voorhis, Phys. Rev. Lett. , 063004 (2009). J. Klimeˇs, D. R. Bowler, and A. Michaelides, J. Phys.: Condens. Matter , 022201 (2010). V. R. Cooper, Phys. Rev. B , 161104(R) (2010). K. Lee, ´E. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, Phys. Rev. B , 081101(R)(2010). O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. , 244103 (2010). J. Klimeˇs, D. R. Bowler, and A. Michaelides, Phys. Rev. B , 195131 (2011). J. Wellendorff, K. T. Lundgaard, A. Møgelhøj, V. Petzold, D. D. Landis, J. Nørskov, T. Bligaard, andK. W. Jacobsen, Phys. Rev. B , 235149 (2012). R. Sabatini, T. Gorni, and S. de Gironcoli, Phys. Rev. B , 041108(R) (2013). K. Berland and P. Hyldgaard, Phys. Rev. B , 035412 (2014). I. Hamada, Phys. Rev. B , 121103(R) (2014). M. Obata, M. Nakamura, I. Hamda, and T. Oda, J. Phys. Soc. Jpn. , 093701 (2013). G. Rom´an-P´erez and J. M. Soler, Phys. Rev. Lett. , 096102 (2009). J. Wu and F. Gygi, J. Chem. Phys. , 224107 (2012). J. A. White and D. M. Bird, Phys. Rev. B , 4954 (1994). R. Sabatini, E. K¨u¸c¨ukbenli, B. Kolb, T. Thonhauser, and S. de Gironcoli, J. Phys.: Condens. Matter , 424209 (2012). J. P. Perdew, K. Burke, and M. Ernzehof, Phys. Rev. Lett. , 3865 (1996). R. Car and M. Parrinello, Phys. Rev. Lett. , 2471 (1985). T. Oda, A. Pasquarello, and R. Car, Phys. Rev. Lett. , 3622 (1998). T. Oda, J. Phys. Soc. Jpn. , 519 (2002). D. Vanderbilt, Phys. Rev. B , 7892 (1990). K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vanderbilt, Phys. Rev. B , 10142 (1993). J. P. Perdew and Y. Wang, Phys. Rev. B , 8800(R) (1986). ´E. D. Murray, K. Lee, and D. C. Langreth, J. Chem. Theory Comput. , 2754 (2009). H. J. Monkhorst and J. D. Pack, Phys. Rev. B , 5188 (1976). J. P. Perdew and A. Zunger, Phys. Rev. B , 5048 (1981). Y. Zhang and W. Yang, Phys. Rev. Lett. , 890 (1998). F. Tran and J. Hutter, J. Chem. Phys. , 204103 (2013). K. Ro´sciszewski, B. Paulus, P. Fulde, and H. Stoll, Phys. Rev. B , 5482 (2000). J. Harl and G. Kresse, Phys. Rev. B , 045136 (2008). O. G. Peterson, D. N. Batchelder, and R. O. Simmons, Phys. Rev. , 703 (1966). L. A. Schwalbe, R. K. Crawford, H. H. Chen, and R. A. Aziz, J. Chem. Phys. , 4493 (1977). J. Wittlinger, R. Fischer, S. Werner, J. Schneider, and H. Schulz, Acta Cryst. B , 745 (1997). T. Ishikawa, M. Asano, M. Obata, N. Suzuki, T. Oda, I. Hamada, and K. Shimizu, private communica-tion. L. Spanu, S. Sorella, and G. Galli, Phys. Rev. Lett. , 196401 (2009). S. Leb`egue, J. Harl, T. Gould, J. G. ´Angy´an, G. Kresse, and J. F. Dobson, Phys. Rev. Lett. , 196401(2010). T. Bj¨orkman, A. Gulans, A. V. Krasheninnikov, and R. M. Nieminen, Phys. Rev. Lett. , 235502(2012). Y. Baskin and L. Meyer, Phys. Rev. , 544 (1955). R. Zacharia, H. Ulbricht, and T. Hertel, Phys. Rev. B , 155406 (2004). I. Hamada and M. Otani, Phys. Rev. B , 153412 (2010). R. E. Mapasha, A. M. Ukpong, and N. Chetty, Phys. Rev. B , 205402 (2012). G. Graziano, J. Klimeˇs, F. Fernandez-Alonso, and A. Michaelides, J. Phys.: Condens. Matter , 424216(2012). rVV10 was also implemented in the JuNoLo code using the Rom´an-P´erez-Soler (RPS) algorithm. Theimplementation of the RPS algorithm is described by Callsen et al . P. Lazi´c, N. Atodiresei, M. Alaei, V. Caciuc, S. Bl¨ugel, and R. Brako, Comput. Phys. Commun. ,371 (2010). M. Callsen, N. Atodiresei, V. Caciuc, and S. Bl¨ugel, Phys. Rev. B , 085439 (2012). M. Callsen and I. Hamada, private communication. H. E. Swanson, N. T. Gilfrich, and G. M. Ugrinic , National Bureau of Standerds Circular 539, Vol. 5(U. S. Government Printing Office, Washington, D. C., 1955) p. 54. P. Cherin and P. Unger, Inorg. Chem. , 1589 (1967). A. D. Corso and R. Resta, Phys. Rev. B , 4327 (1994). G. Kresse, J. Furthm¨uller, and J. Hafner, Phys. Rev. B , 13181 (1994). T. Buˇcko, J. Hafner, S. Leb`egue, and J. G. ´Angy´an, J. Phys. Chem. A , 11814 (2010). O. Sode, M. Ke¸c, K. Yagi, and S. Hirata, J. Chem. Phys. , 074501 (2013). A. Simon and K. Peters, Acta Cryst. B , 2750 (1980). A. Otero-de-la-Roza and E. R. Johnson, J. Chem. Phys. , 054103 (2012). S. A. Bonev, F. Gygi, T. Ogitsu, and G. Galli, Phys. Rev. Lett. , 065501 (2003). I. Hamada, private communication. F. Corsetti, E. Artacho, J. M. Soler, S. S. Alexandre, and M.-V. Fernandez-Serra, J. Chem. Phys.139