First-principles structural optimization and electronic structure of the superconductor picene for various potassium doping levels
Taichi Kosugi, Takashi Miyake, Shoji Ishibashi, Ryotaro Arita, Hideo Aoki
aa r X i v : . [ c ond - m a t . s up r- c on ] N ov First-principles structural optimization and electronic structure ofthe superconductor picene for various potassium doping levels
Taichi Kosugi , Takashi Miyake , , Shoji Ishibashi , Ryotaro Arita , , , and Hideo Aoki Nanosystem Research Institute “RICS”, AIST, Umezono, Tsukuba 305-8568, Japan Japan Science and Technology Agency (JST), CREST, Honcho, Kawaguchi, Saitama 332-0012, Japan Department of Applied Physics, University of Tokyo, Hongo, Tokyo 113-8656, Japan Japan Science and Technology Agency (JST), PRESTO, Kawaguchi, Saitama 332-0012, Japan and Department of Physics, University of Tokyo, Hongo, Tokyo 113-0033, Japan
We theoretically explore the crystal structures of K x picene, for which a new aromatic supercon-ductivity has recently been discovered for x = 3, by systematically performing first-principles fullstructural optimization covering the concentration range x = 1-4. The crystal symmetry (spacegroup) of the pristine picene is shown to be preserved in all the optimized structures despite sig-nificant deformations of each picene molecule and vast rearrangements of herringbone array ofmolecules. For K x picene ( x = 1-4) optimization indicates that (i) multiple structures exist in somecases, and (ii) dopants can enter in the intralayer region as well as in the interlayer region betweenthe stack of herringbone structures. In the electronic structure obtained with the local densityapproximation for the optimized structures, the rigid-band approximation is invalidated for multi-fold reasons: the dopants affect the electronic properties not only through the rearrangement anddistortion of molecules, but also through hybridizations between the molecules and metal atoms.As a consequence the resultant Fermi surface exhibits a variety of multiband structures which takediverse topology for K picene and K picene. PACS numbers: 74.20.Pq, 74.70.Kn, 74.70.Wz
I. INTRODUCTION
The discovery of superconductivity in potassium-doped solid picene, K x picene, with T c = 7-18 K for adoping level x ≃
3, recently reported by Mitsuhashi etal. [1], is seminal as the first aromatic superconductor .This is expected to open a new avenue in quest for or-ganic superconductors. The occurrence of superconduc-tivity in the doped picene should not be an exceptionalcase, since the discovery has been followed by a more re-cent report of superconductivity in another of aromaticcompounds, coronene, by Kubozono et al. [2] with T c up to 15 K for K x coronene with x ≃
3, and also indoped phenanthrene with T c up to 5 K for x ≃ et al. [3]. These aromatic compounds are nowforming a new class in the family of carbon-based su-perconductors that include graphite-intercalation com-pounds (GICs) and doped fullerenes. For GIC, KC , thefirst carbon-based superconductor with T c = 0 . with the highest T c = 11 . C , a superconductor with T c = 20 K [6], is followed by Cs RbC with T c = 33K[7], and Cs C with T c = 40 K under 15 kbar[8] or T c = 38 K under 7 kbar[9].The electronic structures of pristine and doped picenewere first obtained by the present authors[10] as an ini-tial step to elucidate the mechanism of superconduc-tivity. It has been revealed that the conduction bandgroup is derived from the lowest unoccupied molecularorbital (LUMO) and LUMO+1 of a picene molecule.The conduction-band structure of solid picene, where themolecules are arrayed in a herringbone structure withtwo molecules per unit cell, is indeed shown to be accu- rately reproduced by a four-orbital tight-binding modeldownfolded in terms of the maximally-localized Wannierfunctions [11]. We have also determined the dopant po-sitions for K picene, where we have employed the exper-imental lattice parameters. Andres et al. [12], on theother hand, provided two, fully-optimized geometries ofK picene, one with a herringbone structure and the otherwith a laminar structure. Giovannetti and Capone[13]performed an electronic structure calculations with dif-ferent exchange-correlation functionals, and found that ahybrid functional leads to an antiferromagnetic state inK picene. Kim et al. [14] demonstrated that an antiferro-magnetic state is stabilized in K picene when the volumeis increased. They, as well as Giovannetti and Capone,suggested that a strong electron correlation exists in thesystem. This is reasonable in that the electronic correla-tion energy in solid picene is estimated to be U ≃ . et al. [16], who have estimated the phonon frequenciesand the electron-phonon coupling for various hydrocar-bon molecules. For picene they found that the electron-phonon coupling is large enough to attain T c ≃
10 K.Subedi and Boeri[17] calculated the electron-phonon in-teractions in doped picene within the rigid-band approxi-mation, and obtained coupling constants large enough toexplain the experimental T c ≃
18 K in K picene. Casula et al. [18] performed a more elaborate calculation, andobtained the electron-phonon coupling large enough toexplain the experimental T c ≃ picene as opposed to other values of x ? Thequestion is important, since it should be intimately re-lated to another question: to what extent the aromaticsuperconductivity is general. We can for instance re-member that the problem of why A fullerene ( A : alkalimetal elements) is special played an important role inelucidating its electronic properties. One experimentallyunresolved situation is that, in the doped picene, both itsatomic structure and even the number of dopant atomsactually accommodated in the samples have yet to beexperimentally determined, although there are some es-timates as discussed in Ref.[2]. Given the circumstance, atheoretical prediction for the geometry of doped picene isthus desirable, especially for systematically varied dop-ing levels. This will also give a starting point to dis-cuss whether the superconductivity is driven by a phononmechanism, or by an electronic mechanism as in someclasses of organic superconductors[19].This is exactly the purpose of the present work. Fororganic materials, not only the atomic positions but alsothe unit cell should be relaxed for a reliable optimizationof the geometry. Hence we have performed here the fullstructural optimization of K x picene, covering the concen-tration range x = 1-4. II. METHOD
For the electronic structure calculation and the struc-tural optimization in the present work, we adopt thedensity-functional theory with the projector augmented-wave (PAW) method [20] with the Quantum MAterialsSimulator (QMAS) package [21] in the local density ap-proximation (LDA)[22]. The Bloch pseudofunctions areexpanded in plane waves, with an energy cutoff of 40Ry throughout the present work. All the self-consistentfield calculations are performed with 4 × × k points.The structural optimization allows the atomic coordi-nates and the cell parameters to be relaxed with no con-straint imposed on them. In each iteration of the self-consistent field calculation, the forces and stress[23] feltby the system are calculated. The optimization is re-garded as converged when the largest force and stressare smaller than 5 × − Ht/ a B and 5 × − Ht/ a , re-spectively. We employ the maximally-localized Wannierfunction technique[11], which has been recognized to bea standard method for describing electronic orbitals insolid. There, Wannier functions are constructed fromthe Bloch wave functions belonging to the specified elec-tronic bands. In the present work we extract the transferintegrals between the localized Wannier functions con- structed for the solid picene. III. RESULTS
Pristine solid picene comprises molecular layersstacked in the c direction, as depicted in Fig. 1. Ineach layer, picene molecules are in a herringbone arrange-ment, which means a unit cell contains two molecules.The crystal lattice is monoclinic with space group P .Its lattice parameters are reported to be a = 8 . , b =6 . , c = 13 .
515 in ˚A, and the angle β = 90 . ◦ be-tween a and c axes from a single-crystal X-ray diffractionexperiment.[25] The samples used by Mitsuhashi et al. [1]had the same lattice parameters.Let us first classify the candidate positions into whichdopants can be accommodated. A conventional organicinsulator, solid pentacene, where each molecule has thesame chemical formula C H as picene but five ben-zene rings are connected on a straight line with zigzagedges, is known to accommodate dopants in its inter-layer region when doped with alkali elements,[26] with its c axis significantly expanded. By contrast, the c axis ofpicene is reported to shrink when potassium-doped with x = 2 . ab plane of the pristine picene as initial geometries,and found that the intralayer insertion is in fact energet-ically favorable. In the present work, we hence exploreoptimized geometries of the dopants with intralayer in-sertion.The optimized lattice parameters and the internal co-ordinates of potassium atoms are summarized in Table I.The crystal symmetry of pristine picene was found to bepreserved in all the optimized doped structures, despitesignificant deformations of each molecule from the planarshape along with significant rearrangements of moleculesupon doping. The preserved symmetry contrasts withthe situation in K coronene, for which we have shownthat doping lowers the symmetry.[27] The geometries,the electronic band structures, the densities of states,and the Fermi surfaces for the optimized structures areshown for K picene in Fig. 2, for K picene in Fig. 3, forK picene in Fig. 4, and for K picene in Fig. 5. In eachcase, the conduction band group comprises four bands,which have main characters of LUMO and LUMO+1 ofan isolated picene molecule. Since the herringbone struc-ture has two molecules per unit cell, we have constructedfour maximally-localized wave functions (WFs) from thefour bands around the Fermi level of the individual struc-tures. For each of them, two WFs, w h and w l , localizedat each molecule with higher and lower orbital energieswere obtained as shown in Fig. 6. ac molecular layer(intralayer region)interlayer regionΓ X AYZ CD E a * b * c * (a)(b) FIG. 1: (Color online) (a) Crystal structure of pristine solidpicene (drawn with VESTA[24]), which comprises molecularlayers stacked in the c direction. Solid lines delineate unitcells. (b) Typical first Brillouin zone of the doped picene op-timized in the present work. Due to the crystal symmetry, wehave a ∗ ⊥ b ∗ and b ∗ ⊥ c ∗ , while c ∗ and a ∗ are not orthogonalto each other. A. K picene Let us start with K picene. Two structures, K (A) andK (B), have been obtained. K (A) appears when we putone potassium atom per molecule in the intralayer regionas the initial position of the optimization, while K (B)appears when we put the dopant in interlayer. In bothstructures, the dopants end up with intralayer positionsafter the optimization. The insertion of the dopants intothe intralayer region supports the result of our previouswork.[10] In K (B), with a total energy higher than inK (A) by 0 .
316 eV per molecule, the picene moleculesare more strongly twisted than in K (A). Both systemshave metallic electronic band structures. The width ofthe band group consisting of four around their Fermi levelis 0 .
36 for K (A), and 0 .
52 eV for K (B). B. K picene K picene structure was obtained by putting one potas-sium atom in the intralayer region and another in the TABLE I: Optimized lattice constants a, b, c (in ˚A), angle β between a and c axes, volume V of a unit cell, and internalcoordinates of dopants in parentheses for K x picene with x =1-4. a b c β V K (A) 7 .
267 7 .
387 12 .
698 95 .
79 688 . . , . , . (B) 8 .
234 6 .
554 12 .
807 94 .
65 678 . . , . , . .
237 7 .
506 12 .
624 94 .
20 683 . . , . , . . , . , . K .
776 6 .
394 13 .
346 94 .
03 747 . . , . , . . , . , . . , . , . (A) 7 .
421 7 .
213 14 .
028 104 .
53 726 . . , . , . . , . , . . , . , . (B) 7 .
408 7 .
223 14 .
116 105 .
93 726 . . , . , . . , . , . . , . , . .
511 7 .
058 14 .
230 102 .
96 735 . . , . , . . , . , . . , . , . . , . , . interlayer region as the initial positions of the optimiza-tion, where we confirmed that an additional atom dopedinto K picene again prefers the intralayer region by mov-ing to an intralayer position after the optimization. Theobtained K picene is an insulator, with a band gap of0 .
038 eV at X. The widths of the highest occupied andthe lowest unoccupied band groups are 0 .
21 and 0 .
29 eV,respectively. By analyzing the projected density of statesfor the WFs, we find that each of the two band groupshas comparable weights for the two types of the Wannierfunctions w h and w l . C. K picene For K picene we first obtaine K picene(A) structure byputting three potassium atoms in the intralayer regionas the initial positions for the optimization. Andres etal. [12] obtained the monoclinic K picene structure with a = 7 . , b = 7 . , c = 14 .
018 in ˚A and β = 105 . ◦ ,in which the potassium atoms are all in the intralayer re-gion. We have also tried a structural optimization start- E n e r gy ( e V ) DOS (eV ) -1
10 200 E n e r gy ( e V ) DOS (eV ) -1
10 200Γ X A Y ZΓ D E C Z Γ X A Y ZΓ D E C Z a b acc ac bc a a * b * c * b * a * b * a * c * a * b * K picene (A) K picene (B) -0.200.20.40.60.8 -0.200.20.40.60.8 FIG. 2: (Color online) Geometries (large spheres: K atoms), electronic band structures, densities of states, and Fermi surfacesof the K structures. The origins of energy are set to the respective Fermi level. E n e r gy ( e V ) DOS (eV ) -1
10 200Γ X A Y ZΓ D E C Z abcac K picene -0.200.20.40.6 FIG. 3: (Color online) Geometry (large spheres: K atoms),electronic band structure, and density of states of theK picene. The origin of energy is set to the Fermi level. ing from their lattice parameters and atomic coordinates.Then the geometry resulted in another structure, whichwe call K picene(B) here. Although the lattice param-eters of K (A) and K (B) are close to each other withthe difference in the total energy being only 1 meV permolecule, they are distinct in the positions of dopants, asseen in Table I and in Fig. 4. K picene(B) has the lat-tice parameters and the dopant positions close to thoseof the K picene structure obtained by Andres et al. [12],and we therefore consider the two structures to be thesame. K (A) and K (B) have metallic electronic bandstructures with the same width 0 . x ≃ K picene, by putting two potassium atomsin the intralayer region and one in the interlayer re-gion as the initial positions for the optimization. Thisstructure is higher in the total energy than the K (A)structure by 0 .
465 eV per molecule. However, the lat-tice parameters are closer to the experimental results[1]of a = 8 . , b = 5 . , c = 12 .
97 in ˚A, and β = 92 . ◦ for K . picene. The electronic band structure is againmetallic for the K K structure. The four bands aroundthe Fermi level, whose width is 0 .
51 eV, are close to thehigher bands as compared with those in the K (A) andK (B) structures.Each of K (A) and K (B) has multiple Fermi sur-faces comprising a pocket and complicated, connectedsurfaces. The existence of a Fermi pocket implies three-dimensional electron transfer integrals despite the layeredstructure. On the other hand, K K has one-dimensional(a pair of planar) Fermi surfaces along with connectedsurfaces. We can also note that, while the Fermi sur-face of K coronene[27] also consists of multiple sur-faces, the topology of the Fermi surfaces for K coronene(consisting of a pair of quasi-one-dimensional, planarsurfaces and another pocket) differs from those foundhere for K picene(A) and K picene(B) (consisting ofmultiple surfaces of two- or higher dimensional charac-ters), and for K K picene (consisting of a pair of quasi-one-dimensional, planar surfaces and another, multiply-connected one).The largest interlayer transfer integral of the WFs inK (A) (K K ) is 13 . .
6) meV, while that in undopedpicene is 14 . K than in K (A) originate in the presenceof dopants in the interlayer region. The widths of thebands around the Fermi level in these doped structuresare larger than in undoped picene (0.39 eV[10]), thoughK (A) has the smaller interlayer transfer integrals. Itimplies that their larger band widths come from the en-hanced intralayer transfers of electrons due to the hy-bridization of the wave functions of potassium atoms andthe molecular orbitals. In fact, the largest intralayertransfer integral in K (A) (K K ) is 82 . .
6) meV,which is much larger than that in undoped picene, 57 . (A)picene and K K picene are shown in Fig.7. It is seen that K (A) allows for the transfers betweenmore distant WFs in the a direction, which leads to itsmore two-dimensional and wider band dispersion thanK K . D. K picene We have obtained the K picene structure by puttingthree potassium atoms in the intralayer region and onein the interlayer region as the initial positions of theoptimization, because the three dopants sitting on ev-ery other benzene rings out of five rings in each picenemolecule may tend to push the fourth dopant to the in-terlayer region. In the optimized structure, all the fourdopants end up with sitting in the intralayer region, buttwo of them are rather close to the interlayer region dueto the crowded space. The system is an insulator with aband gap of 0 .
36 eV. The highest occupied band groupconsists of two bands, separated from the second highestgroup of two bands by 0 .
01 eV. The widths of the higherand lower band groups are 0 .
17 and 0 .
20 eV, respectively. By analyzing the projected density of states for the WFs,we found that the higher and lower occupied band groupsconsist mainly of w h and w l , respectively. This contrastswith the situation in another insulator, K picene. E. Discussions
For each concentration x , the lowest-energy structureshown above contains the dopants in the intralayer re-gion. This behavior is consistent with the observedtendency[2] of the c axis of K x picene, in which the c axis shrinks monotonically with increasing x , implyingthat the dopants move into the intralayer region regard-less of the concentration. This is to be contrasted withRb x picene, whose c axis was observed to expand up to x = 2 and shrink afterwards.In Fig. 6 we have displayed the maximally-localizedWannier functions ( w h and w l with higher and lowerorbital energies, respectively) in the optimized struc-tures for all the cases of K picene to K picene. Wecan see that the WFs have basically the same featuresfor all the structures in that all the WFs have LUMOand LUMO+1 orbital characters, although small but sig-nificant amplitudes found around the potassium atomswhose magnitudes depend on the doping level x . Moreprecisely, we can show that w h and w l have charac-ters of the linear combinations of LUMO and LUMO+1: w h , l = φ LUMO ± φ LUMO+1 , as in the undoped picene[10].However, a key message in the present band struc-ture results is that, although the orbital characters ofthe WFs are common to all the optimized geometries,their band structures very much differ. This is not sur-prising, since the molecular shape as well as the herring-bone arrangement of molecules distort significantly withdoping. More importantly, the lobes of molecular WFsextend outward the molecule via the potassium atoms,with significant hybridizations of WFs around the potas-sium atoms. Specifically, we find that the small contribu-tions from the potassium wave functions have a p orbitalcharacter, as seen in Fig. 6.In other words, the rigid-band picture is invalidatedfor multifold reasons, i.e., through the rearrangementand distortion of molecules, but also via hybridizationof the potassium p orbitals with the molecular orbitalsof picene. Indeed, for all the optimized geometries wehave also calculated their band structures (not shown)with potassium atoms removed but the doping level isfixed by shifting the chemical potential to find that theband structures are drastically changed from those of thedoped systems. As a consequence the resultant Fermisurface exhibits a variety of multiband structures whichtake diverse topology for K picene and K picene. Thisobviously provides an intriguing starting point for explor-ing mechanisms of superconductivity as a future work.Final point of interest, about the molecular species de-pendence, is the following: we can note that the cal-culated band structures for all x (including x = 0[10]) E n e r gy ( e V ) DOS (eV ) -1
10 200
DOS (eV ) -1
10 200 E n e r gy ( e V ) -0.4-0.200.20.40.6 E n e r gy ( e V ) DOS (eV ) -1
10 200Γ X A
Y ZΓ D E C Z Γ X A Y ZΓ D E C ZΓ X A Y ZΓ D E C Z abc abcabcac cc a a a * c * b *K picene (A) K picene (B)K K picene b * b * b * a * a * a * c * c * b * b * a * a * -0.4-0.200.20.40.6 -0.4-0.200.20.40.6 FIG. 4: (Color online) Geometries (large spheres: K atoms), electronic band structures, densities of states, and Fermi surfacesof the K picene. The origins of energy are set to the respective Fermi level. have significant dispersion along the c direction as wellas along the a and b directions, which indicates isotropiccomponents in the electron transfers. If we compare thiswith the situation in K x coronene, we find that the band dispersion is more anisotropic in both undoped coroneneand K coronene[27], which reflects the differences be-tween the intermolecular distances in each direction.
10 200
DOS (eV ) -1 Γ X A Y ZΓ D E C Z E n e r gy ( e V ) c a abc K picene -0.4-0.200.20.4 FIG. 5: (Color online) Geometry (large spheres: K atoms),electronic band structure, density of states, and Fermi sur-faces of the K picene. The origin of energy is set to the Fermilevel. IV. SUMMARY
We have systematically performed a full structuraloptimization for K x picene, covering the doping level x = 1-4. The crystal symmetry of pristine picene of P was found to be preserved in all the optimizedstructures despite the deformation and rearrangementsof the molecules upon doping. The optimized systemsfor K picene and K picene have metallic band struc-tures, while band gaps open around their Fermi levelsfor K picene and K picene. In each of the optimized sys-tems, the electronic bands in the vicinity of the Fermilevel are derived from LUMO and LUMO+1 of the aro-matic molecule as in the undoped picene. We found that the presence of dopants affects the electronic propertiesof all the doped structures. The electronic structuresobtained with the local density approximation for theoptimized structures reveal that the rigid-band approxi-mation is invalidated due to the rearrangement and de-formation of molecules and the accompanied charge re-distribution between the molecules and dopants. Theresultant Fermi surface exhibits a variety of multibandstructures which take diverse topology for K picene andK picene.In all the optimized systems, the widths of the LUMO-and LUMO+1-derived bands do not exceed 0 . U ≃ .
85 eV,[15] indicating that the mate-rial may be a strongly correlated electronic system. Fur-ther investigations for clarifying the material propertiesand the mechanism of superconductivity in doped piceneshould thus take into account the magnetic instabilityand the strong Coulomb interaction, as pointed out byGiovanetti and Capone[13] and Kim et al. [14]Rb . picene and Ca . picene have also been reportedto become superconducting with T c ≃ x picene is ex-pected to accommodate dopants in the interlayer regiondue to a larger atomic radius of Rb than in K x picene. Acknowledgments
HA thanks Yoshihiro Kubozono for discussions. Thepresent work is partially supported by the Next Genera-tion Supercomputer Project, Nanoscience Program fromMEXT, Japan, and by Grants-in-aid No. 19051016 and22104010 from MEXT, Japan and the JST PRESTOprogram. The calculations were performed with com-putational facilities at TACC, AIST as well as at Super-computer Center of ISSP and at Information TechnologyCenter, both at University of Tokyo. [1] R. Mitsuhashi, Y. Suzuki, Y, Yamanari, H. Mitamura, T.Kambe, N. Ikeda, H. Okamoto, A. Fujiwara, M. Yamaji,N. Kawasaki, Y. Maniwa, and Y. Kubozono, Nature ,76 (2010).[2] Y. Kubozono, M. Mitamura, X. Lee, X. He, Y. Yamanari,Y. Takahashi, Y. Suzuki, Y. Kaji, R. Eguchi, K. Akaike,T. Kambe, H. Okamoto, A. Fujiwara, T. Kato, T. Kosugi,and H. Aoki, Phys. Chem. Chem. Phys. , 16476 (2011).[3] X. F. Wang, R. H. Liu, Z. Gui, Y. L. Xie, Y. J . Yan, J.J. Ying, X. G. Luo, and X. H. Chen, Nature Communi-cations , 507 (2011).[4] N. B. Hannay, T. H. Geballe, B. T. Matthias, K. Andres,P. Schmidt, and D. MacNair, Phys. Rev. Lett. , 225(1965). [5] T. E. Weller, M. Ellerby, S. S. Saxena, R. P. Smith andN. T. Skipper, Nature Phys. , 39 (2005); N. Emery,C. H´erold, M. d’Astuto, V. Garcia, Ch. Bellin, J. F.Marˆech´e, P. Lagrange and G. Loupias, Phys. Rev. Lett. , 087003 (2005).[6] A. F. Hebard, M. J. Rosseinsky, R. C. Haddon, D. W.Murphy, S. H. Glarum, T. T. M. Palstra, A. P. Ramirezand A. R. Kortan, Nature , 600 (1991).[7] K. Tanigaki, T. W. Ebbesen, S. Saito, J. Mizuki, J.S. Tsai, Y. Kubo and S. Kuroshima, Nature , 222(1991).[8] T. T. M. Palstra, O. Zhou, Y. Iwasa, P. E. Sulewski, R.M. Fleming and B. R. Zegarski, Solid State Commun. , 327 (1995). w h w l w h w l w h w l w h w l w h w l (a) K picene(A) (b) K piceneK (c) picene(A) K K (d) picene K (e) picene ac FIG. 6: (Color online) Maximally-localized Wannier functions in the optimized structures for K picene (a), K picene (b),K picene(A) (c), K K picene (d), and K picene (e). w h and w l stand for the higher- and lower-energy states, respectively. ba K (a) picene(A)K (b) piceneK baw l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h w l w h FIG. 7: (Color online) Significant ( | t | >
30 meV) transfer in-tegrals t between Wannier functions, two at each molecule, aredepicted with arrows for K picene(A) (a), and K K picene(b). [9] A. Y. Ganin, Y. Takabayashi, Y. Z. Khimyak, S. Mar-gadonna, A. Tamai, M. J. Rosseinsky, and K. Prassides,Nature Materials , 367 (2008).[10] T. Kosugi, T. Miyake, S. Ishibashi, R. Arita, and H. Aoki,J. Phys. Soc. Jpn. (2009) 113704[11] N. Marzari and D. Vanderbilt, Phys. Rev. B , 12847(1997); I. Souza, N. Marzari, and D. Vanderbilt, Phys.Rev. B , 035109 (2001).[12] P. L. de Andres, A. Guijarro, and J. A. Verg´es, Phys.Rev. B , 245113 (2011).[13] G. Giovannetti and M. Capone, Phys. Rev. B , 134508(2011).[14] M. Kim, B. I. Min, G. Lee, H. J. Kwon, Y. M. Rhee, andJ. H. Shim, Phys. Rev. B , 214510 (2011).[15] F. Roth, B. Mahns, B. B¨uchner, and M. Knupfer, Phys.Rev. B , 165436 (2011).[16] T. Kato and T. Yamabe, J. Chem. Phys. , 8592(2001); T. Kato, K. Yoshizawa and K. Hirao, J. Chem.Phys. , 3420 (2002); T. Kato and T. Yamabe, Chem.Phys. , 437 (2006).[17] A. Subedi and L. Boeri, Phys. Rev. B , 020508(R)(2011).[18] M. Casula, M. Calandra, G. Profeta, and F. Mauri, Phys.Rev. Lett. , 137006 (2011).[19] K. Kuroki, J. Phys. Soc. Jpn. , 051013 (2006).[20] P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994); G. Kresseand D. Joubert, ibid. , 566(1980); J. P. Perdew and A. Zunger, Phys. Rev. B ,5048 (1981). [23] O. H. Nielsen and Richard M. Martin, Phys. Rev. B ,3780 (1985).[24] K. Momma and F. Izumi, J. Appl. Crystallogr. , 653(2008).[25] A. De, R. Ghosh, S. Roychowdhury and P. Roychowd-hury, Acta Cryst. C , 907 (1985).[26] T. Minakata, I. Nagoya and M. Ozaki, J. Appl. Phys. , 7354 (1991); T. Minakata, H. Imai and M. Ozaki,J. Appl. Phys. , 4178 (1992); T. Ito, T. Mitani, T. Takenobu and Y. Iwasa, J. Phys. Chem. Solids , 609(2004); Y. Matsuo, T. Suzuki, Y. Yokoi and S. Ikehata, J.Phys. Chem. Solids , 619 (2004); Y. Matsuo, S. Sakaiand S. Ikehata, Phys. Lett. A , 62 (2004); B. Fang,H. Zhou and I. Homma, Appl. Phys. Lett. , 261909(2005).[27] T. Kosugi, T. Miyake, S. Ishibashi, R. Arita, and H. Aoki,Phys. Rev. B84