Crystal and electronic structures of superhard B2CN : An ab initio study

Solid State Communications 152 (2012) 71–75

Contents lists available at SciVerse ScienceDirect

Solid State Communications journal homepage: www.elsevier.com/locate/ssc

Crystal and electronic structures of superhard B2 CN: An ab initio study Quan Li a,b , Dan Zhou c , Hui Wang c,∗,1 , Wanjin Chen d,∗,2 , Baojia Wu c,e,∗∗ , Zhijian Wu b , Weitao Zheng a a

College of Materials Science and Engineering, Jilin University, Changchun 130012, China

b

State Key Laboratory of Rare Earth Resource Utilization, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun 130022, China

c

State Key Lab of Superhard Materials, Jilin University, 2699 Qianjin Str., Changchun 130012, China

d

College of Physics, Jilin Normal University, 1301 Haifeng Str., Siping 136000, China

e

Department of Physics, College of Science, Yanbian University, Yanji, Jilin 133002, China

article

info

Article history: Received 7 September 2011 Received in revised form 26 October 2011 Accepted 29 October 2011 by J. R. Chelikowsky Available online 6 November 2011 Keywords: A. Superconductors C. Structure prediction D. Superhard

abstract With the increasing demand for specific applications in high pressure and electronic devices, the search for superhard superconducting materials remains a topic of great interest. Using particle swarm optimization algorithm, we report five competitive structures with clear tetrahedrally sp3 hybridization, among which two metallic orthorhombic structures (Pmm2 and Pmma) with the maximal stable bonds (C–C + B–N) are energetically more favorable than earlier proposed structures. Further first principles calculations suggest that the predicted five structures possess simultaneously superhard and superconducting properties. The five structures are with the similar calculated hardness (56–58 GPa), but show distinct difference in superconducting critical temperature, ranging from 2 to 53 K (with the Coulomb parameter µ∗ of 0.13). © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Materials combining superhard (Vickers hardness exceeding 40 GPa) and metallic (or superconducting) properties are of considerable demand for the creation of multifunctional devices, e.g., high-pressure devices for investigating the electric properties of various materials under extreme conditions [1]. Over the past several decades, scientists mainly focused on the covalent B–C–N compounds with short and strong three dimensional covalent bonds to search for new superhard materials. Furthermore, since the discovery of superconductivity in boron-doped diamonds (4 K) [2], the superconducting behaviors in ‘‘covalent metals’’ (dominated by strong covalent bonds, rather than the metallic bonding) have also been an exotic subject [3]. It is suggested that the strongly directional covalent character in metallic compounds leads to large phonon frequencies and electron–phonon coupling potential, which contribute to increase in Tc [4]. Thus, metallic B–C–N compounds with exclusively covalent

∗ ∗∗

Corresponding authors. Correspondence to: Department of Physics, College of Science, Yanbian University, Yanji, Jilin 133002, China. Tel.: +86 433 2732221. E-mail addresses: [email protected] (Q. Li), [email protected] (H. Wang), [email protected] (W. Chen), [email protected] (B. Wu). 1 Tel.: +86 431 85168276. 2 Tel.: +86 434 3291986. 0038-1098/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ssc.2011.10.042

bond can be expected to simultaneously possess superhard and superconducting properties. Metallicity of ternary B–C–N compounds can be evaluated by the following simple rule: pZV (B) + mZV (C) + lZV (N) ̸= 4n [5]. The values p, m, l, and n are integers, and ZV (B), ZV (C), and ZV (N) are the atomic valence states (2s and 2p) for B, C, and N, respectively. Recently, diamond-like BCx N (x = 1 − 4) [6–10], BC5 [11], and B2 CN [12] crystals have been synthesized by highpressure and high-temperature (1500–2500 K) conditions, in which BC5 and B2 CN are suggested to be conductors. Diamondlike B2 CN was synthesized by boric acid and melamine as starting materials under a pressure of 5.5 GPa [12] which is relatively lower than the synthetic pressures for BCx N (7.7–50 GPa) [6–10] and BC5 (24 GPa) [11]. The superconducting properties of BC5 have been widely explored [13–15]; however, no such study has been performed for diamond-like B2 CN, since the crystal structure of the synthesized diamond-like B2 CN is not known, which is the key to reflect the physical properties of materials. Based on different atomic configurations in the eight-atom zinc-blende cubic unit cell, He et al. [16] proposed seven possible structures with the lowest energy phase adopting a tetragonal structure with space group P-4m2 (t-B2 CN). However, the proposed ground state structure is based on the knowledge of known crystal information. There is a possibility that hitherto unknown structures are stable instead. In this paper, we have adopted our newly developed particle swarm optimization (PSO) algorithm [17–19] to explore the structures of B2 CN at large pressure ranges of 0–100 GPa. Two

72

Q. Li et al. / Solid State Communications 152 (2012) 71–75

Fig. 1. (Color online) Structures of B2 CN predicted in this work. (a) t-B2 CN (P-4m2). Atomic positions: B at 2g (0.5, 0, 0.2701), C at 1b (0.5, 0.5, 0), and N at 1d (0, 0, 0.5). (b) P3m1. Atomic positions: B at 1a (0, 0, 0.1705) and 1b (0.3333, 0.6667, 0.6418), C at 1b (0.3333, 0.6667, 0.004), and N at 1a (0, 0, 0.5127). (c) R3m, Atomic positions: two inequivalent B at 1a (0.7322, 0.7322, 0.7322) and 1a (0.2234, 0.2234, 0.2234), respectively, C at 1a (0.3451, 0.3451, 0.3451), and N at 1a (0.8478, 0.8478, 0.8478). (d) Pmm2. Atomic positions: four inequivalent B at 1b (0, 0.5, 0.9091), 1c (0.5, 0, 0.1374), 1c (0.5, 0, 0.6811), and 1d (0.5, 0.5, 0.5092), two inequivalent C at 1a (0, 0, 0.2747) and 1b (0, 0.5, 0.3797), and two inequivalent N at 1a (0, 0, 0.7910) and 1d (0.5, 0.5, 0.0228). (e) Pmma. Atomic positions: two B at 2e (0.25, 0, 0.6889) and 2f (0.25, 0.5, 0.0845), C at 2f (0.25, 0.5, 0.5517), and N at 2e (0.25, 0, 0.1948). Table 1 The calculated equilibrium structural parameters, volumes (V ) and formation energy (1E ) per 4 atoms, bulk modulus (B0 ), and the radio of stable bonds (C–C + B–N) and the other bonds (R) for B2 CN. B2 CN

t-B2 CN P3m1 R3m Pmm2 Pmma

a (Å)

2.526 2.574 4.640 2.560 2.510

b (Å)

2.526 2.574 4.640 2.507 2.510

c (Å)

3.973 4.442 4.640 7.910 7.898

V (Å3 )

25.350 25.478 25.424 25.370 25.294

R

1:1 1:1 1:1 5:3 5:3

1E (eV) 0.85 0.97 0.95 0.55 0.53

metallic orthorhombic structures (Pmm2 and Pmma) are predicted here and found to be energetically much superior to the earlier structures [16]. Strikingly, further first principles calculations suggest that the predicted B2 CN structures possess simultaneously superhard and superconducting properties. 2. Computational method The search of high-pressure structures (0–100 GPa) was performed with variable-cell PSO simulations [17–20] with the CALYPSO code [21] in conjunction with ab initio structure relaxations within the framework of density-functional theory. The calculations of the structure relaxations, electronic density of states (DOS) were carried out within the Perdew–Burke–Ernzerhof (PBE) parameterization of the generalized gradient approximation (GGA) [22] as implemented in the Vienna ab initio simulation package (VASP) [23]. The all-electron projector augment wave (PAW) method [24,25] was adopted with 1s2 treated as core for B, C and N atoms. We used the plane-wave kinetic energy cutoff energies of 1000 eV, which was shown to give excellent convergence of the total energies (∼1 meV) and structural parameters. The lattice dynamics and electron–phonon coupling were studied using the plane-wave pseudopotential method [26,27] within the PBE–GGA, through the Quantum-ESPRESSO package [28]. Convergence tests gave the choice of kinetic energy cutoffs of 90 Ry.

B0 (GPa)

333 319 318 318 328

λ 0.56 1.25 1.25 0.52 0.42

ωlog (K) 726 526 631 1057 929

Tc (K)

µ∗ = 0.10

µ∗ = 0.13

13.4 49.5 59.5 14.9 5.1

9.0 44.2 53.1 9.4 2.4

the ratio of stable bonds (C–C + B–N) over less stable bonds (B–C + C–N + B–B) of the obtained B2 CN structures have been presented in Table 1. It is noted that the tetragonal P-4m2 structure is identical to the earlier t-B2 CN structure [16]. To the best of our knowledge, the other four structures are firstly reported here. T -B2 CN, Pmm2 and Pmma can be viewed as a periodic sandwichlike structure formed by stacking CBNB, BNBBCCBN and CCBNBBNB along the ⟨100⟩ direction of diamond structure, respectively. Instead, the stacking sequences of P3m1 and R3m follows the ⟨111⟩ direction. With the fixed number of 1 and 2 formula units in the primitive cell, it is found that the energetically best structures are the diamond-[100] structures: t-B2 CN and Pmma, respectively. Pmma structure with ‘‘boron chains’’ and ‘‘carbon chains’’ possesses the lowest energy containing 1 and 2 f.u. per primitive cell, and 320 meV per unit lower than the earlier t-B2 CN structure. The C–C and B–N chemical bonds are more favorable than B–B, B–C and C–N to form diamond-like structure with sp3 hybridized, thus we define C–C and B–N bonds as stable bonds and the other relatively weaker bonds as less stable bonds. In the two most stable structures (Pmma and Pmm2), it is interesting to note that the ratio of the stable chemical bonds over less stable bonds reaches 5:3, which is the largest among the five structures. Further increasing the formula units (3 or 4 f.u./cell) in the simulation cell, we find the produced structures tend to separate into bulk diamond + c-BN + boron layers. The formation energy (1E ) is defined as

3. Results and discussion

1E = E (B2 CN) − E (C) − E (BN) − E (B).

We performed variable-cell structure prediction simulations at a large pressure range 0–100 GPa with systems containing one to two formula units (f.u.) in the simulation cell. As shown in Fig. 1, five competitive structures with clear tetrahedrally sp3 hybridization were predicted: tetragonal P-4m2 (1 f.u./cell), hexagonal P3m1 (1 f.u./cell), rhombohedral R3m (1 f.u./cell), and orthorhombic Pmm2 (2 f.u./cell) and Pmma (2 f.u./cell). The theoretical lattice parameters, volumes, formation energies, and

The diamond, c-BN, and α -B12 were chosen as the reference phases. We find that all the B2 CN structures have positive formation energy (Table 1), indicating that these structures are tend to phase decomposition. The calculated projected electronic DOS for the predicted structures are presented in Fig. 2. The valence electrons in B2 CN populate sp3 hybridized orbitals and give rise to σ -bonds between nearest neighbors. However, the valence bands are partially

Q. Li et al. / Solid State Communications 152 (2012) 71–75

73

Table 2 µ µ The calculated bond lengths (dµ ), bond numbers in the primitive cell (N µ ), valence electron density (Ne ), Mulliken overlap population (P µ ), Phillips ionicity (fi ), and Vickers hardness of µ-type bonds and crystal (Hvµ and Hv ) for B2 CN.

t-B2 CN

P3m1

R3m

Pmm2

Pmma

µ

µ

Bond type

dµ (Å)

Nµ

Ne

Pµ

fi

Hvµ (GPa)

B–N B–C B–N B–N B–C B–C B–N B–N B–C B–C C–C B–N B–N B–N B–N B–C B–C B–B C–C B–N B–N B–C B–B

1.558 1.657 1.520 1.592 1.607 1.660 1.524 1.593 1.604 1.652 1.503 1.547 1.547 1.563 1.564 1.639 1.678 1.850 1.497 1.545 1.555 1.674 1.832

4 4 1 3 1 3 1 3 1 3 2 2 2 2 2 2 2 2 2 4 4 4 2

0.785 0.635 0.824 0752 0.675 0.632 0.818 0.749 0.676 0.637 0.853 0.805 0.805 0.788 0.787 0.656 0.626 0.465 0.859 0.807 0.796 0.629 0.474

0.64 0.66 0.48 0.69 0.51 0.70 0.46 0.70 0.49 0.72 0.81 0.66 0.63 0.63 0.66 0.68 0.65 0.54 0.83 0.67 0.61 0.65 0.56

0.258 0.220 0.537 0.140 0.487 0.140 0.572 0.140 0.521 0.095 0.144 0.220 0.276 0.276 0.220 0.181 0.239 0.435 0.173 0.201 0.312 0.239 0.400

67 48 53 70 40 53 51 69 39 57 91 73 68 65 69 54 45 21 89 75 64 46 23

Hv (GPa) 57

57

58

56

56

It is known that the bulk modulus measures the resistance to volume change, which is not always a good indicator of hardness of a substance. It is known that hardness determination involves both elastic and plastic deformation. Therefore, hardness is too complex to be described by first principles. Fortunately, several hardness models are proposed to evaluate the intrinsic hardness of ideal crystals and are quite successful for the applications to covalent compounds, and even for metals [29–36]. For electronic conductors, the electrons delocalized to contribute to the conduction should be excluded from the hardness calculation. However, the major carriers are holes and the valence electrons are mainly localized to form covalent bonds for hole conductors. It is thus unnecessary to include the metallic correction in the hardness calculation for the hole conductors. We explicitly estimated the intrinsic hardness using Gao’s model as following [30–32]: µ

Hvµ = 350(Neµ )2/3 (dµ )−2.5 e−1.191fi , Hv =

Fig. 2. (Color online) The calculated projected electron density of states for (a) P-4m2, (b) P3m1, (c) R3m, (d) Pmm2, and (e) Pmma at zero pressure, respectively. The horizontal dashed lines represent the Fermi level.

filled since the charges of this system are not balanced due to the electron deficiency (existing holes) of B2 CN in comparison with diamond. Thus, t-B2 CN and Pmma structures show metallic behavior with hole nature. The bulk modulus (B0 ) is obtained by fitting the energy–volume data points to the third-order Birch–Murnaghan equation of state as listed in Table 1. All the predicted structures of B2 CN have large bulk modulus (318–333 GPa), indicating that those structures are very incompressible as expected from their sp3 chemical bondings.

 µ ∏

µ Nµ

(Hv )

1/ ∑ N µ

.

The calculated bond lengths dµ , bond numbers N µ , valence elecµ µ tron density Ne , Phillips ionicity fi [37], and Vickers hardness Hv for B2 CN are listed in Table 2. Here, superscript µ denotes µ-type bond. The previous hardness calculations have shown that the accuracy of the calculated hardness is within 5% [30–32]. Significantly, it is found that the Vickers hardness for the current predicted structures can reach 56–58 GPa, indicative of the superhard property. Since the current theory slightly overestimates the volume as expected by GGA, larger bulk modulus and hardness might be expected in real applications. Phonon calculations give a criterion for the crystal stability and indicate, through soft modes, structural instability. Therefore, the phonon dispersion curves for the five structures were calculated. As shown in Fig. 3, no imaginary phonon frequency is observed for the five structures, indicating the structural stability. It is suggested that strong covalent bonding could lead to large phonon frequencies and a large electron–phonon coupling potential [4]. To explore the superconductivity, the theoretical spectral function α 2 F (ω)/ω [38] and the integrated EPC parameter λ(ω) as a function of frequency have been investigated at ambient pressure. Superconducting Tc can be estimated from the Allen–Dynes

74

Q. Li et al. / Solid State Communications 152 (2012) 71–75

Fig. 3. Phonon dispersion for (a) P-4m2, (b) P3m1, (c) R3m, (d) Pmm2, and (e) Pmma at zero pressure, respectively.

Fig. 4. (Color online) Eliashberg phonon spectral function α 2 F (ω)/ω, and the electron–phonon integral λ(ω) for (a) P-4m2, (b) P3m1, (c) R3m, (d) Pmm2, and (e) Pmma at zero pressure, respectively.

modified McMillan equation Tc =

ωlog 1.2



1.04(1+λ)



exp − λ−µ∗ (1−0.62λ) [38]. ∗ To estimate Tc , the effective Coulomb repulsion   parameter µ must ω

be known. It is defined as µ1∗ = µ1 + ln ω el , where µ is the direct ph Coulomb repulsion between paired electrons, ωel corresponds to a plasma frequency, while ωph the high-frequency cutoff in electron–phonon coupling spectral function α 2 F (ω)/ω [39]. Due to the difficulty in calculating µ and ωph , it remains challenging to directly derive an accurate µ∗ by theory. Fortunately, it is found that the µ∗ value of 0.13–0.10 gives reasonable estimated results for estimating the superconducting Tc for diamond-like B–C–N compounds (e.g., boron-doped diamond [40,41] and diamond-like c-BC5 [13,14]). As shown in Fig. 4 and Table 1, the resulting λs for diamond-[111] are much stronger than those of diamond[100]. It is found that the contributions of the relatively higher

frequency (>20 THz) to the EPC λ are not with large difference for different structures. However, the relatively lower frequency (<20 THz) in diamond-[111] structures contribute 1.12 and 1.06 to the EPC λ for P3m1 and R3m, respectively, which are much larger than those of diamond-[100] structures (0.22–0.35). with the Coulomb parameter µ∗ of 0.13, the estimated Tc for diamond[111] structures reaches a very high value, 44–53 K, which is much higher than those of diamond-[100] structures (2–9 K). We therefore conclude that the large difference in superconductivity of B2 CN is mainly related to the relatively lower frequency (<20 THz). 4. Conclusion In summary, crystal structures of B2 CN have been extensively investigated using ab initio particle swarm optimization algorithm

Q. Li et al. / Solid State Communications 152 (2012) 71–75

at a large pressure range of 0–100 GPa. With the fixed number of 1 and 2 formula units in the primitive cell, it is found that the energetically best structures are the diamond-[100] structures: t-B2 CN and Pmma, respectively. Due to the electron deficiency of sp3 hybridized B2 CN in comparison with diamond, the five diamond-like B2 CN structures show metallic with hole nature. The calculated large bulk modulus (318–333 GPa) and high hardness (56–58 GPa) reveal that they are ultra-incompressible and superhard materials. Further electron–phonon coupling calculations predict that all the five phases are superconductors. The estimated Tc for diamond-[111] structures (R3m and P-3m1) reaches a very high value, 44–53 K, in contrast to the smaller value (2–9 K) in diamond-[100] structures. In view of the excellent properties of diamond-like B2 CN, it is necessary to pay much attention on the experimental synthesis and conductivity measurements. Acknowledgments The authors gratefully acknowledge the National Natural Science Foundation of China (NSFC), Grant Nos. 10874054, 91022029, and 11025418, the China 973 Program, Grant No. 2011CB808200. Parts of the calculations were performed in the High Performance Computing Center (HPCC) of the Jilin University. References [1] G.A. Dubitskiy, V.D. Blank, S.G. Buga, E.E. Semenova, V.A. Kul’bachinskii, A.V. Krechetov, V.G. Kytin, JETP Lett. 81 (2005) 260. [2] E.A. Ekimov, V.A. Sidorov, E.D. Bauer, N.N. Mel’nik, N.J. Curro, J.D. Thompson, S.M. Stishov, Nature 428 (2004) 542. [3] V.H. Crespi, Nature Mater. 2 (2003) 650. [4] X. Blase, E. Bustarret, C. Chapelier, T. Klein, C. Marcenat, Nature Mater. 8 (2009) 375. [5] S.F. Matar, M. Mattesini, C. R. Math. Acad. Sci. Paris 4 (2001) 255. [6] V.L. Solozhenko, D. Andrault, G. Fiquet, M. Mezouar, D.C. Rubie, Appl. Phys. Lett. 78 (2001) 1385. [7] Y. Zhao, D.W. He, L.L. Daemen, T.D. Shen, R.B. Schwarz, Y. Zhu, D.L. Bish, J. Huang, J. Zhang, G. Shen, J. Qian, T.W. Zerda, J. Mater. Res. 17 (2002) 3139.

75

[8] E. Knittle, R.B. Kaner, R. Jeanloz, M.L. Cohen, Phys. Rev. B 51 (1995) 12149. [9] T. Komatsu, M. Nomura, Y. Kakudate, S. Fujiwara, J. Mater. Chem. 6 (1996) 1799. [10] S. Nakano, M. Akaishi, T. Sasaki, S. Yamaoka, Chem. Mater. 6 (1994) 2246. [11] V.L. Solozhenko, O.O. Kurakevych, D. Andrault, Y. Le Godec, M. Mezouar, Phys. Rev. Lett. 102 (2009) 015506. [12] J.L. He, Y.J. Tian, D.L. Yu, T.S. Wang, S.M. Liu, L.C. Guo, D.C. Li, X.P. Jia, L.X. Chen, G.T. Zou, Chem. Phys. Lett. 340 (2001) 431. [13] M. Calandra, F. Mauri, Phys. Rev. Lett. 101 (2008) 016401. [14] Y. Yao, J.S. Tse, D.D. Klug, Phys. Rev. B 80 (2009) 094106. [15] Q. Li, H. Wang, Y. Tian, Y. Xia, T. Cui, J. He, Y. Ma, G. Zou, J. Appl. Phys. 108 (2010) 023507. [16] J.L. He, L.C. Guo, E. Wu, X.G. Luo, Y.J. Tian, J. Phys.: Condens. Matter 16 (2004) 8131. [17] Y. Wang, J. Lv, L. Zhu, Y. Ma, Phys. Rev. B 82 (2010) 094116. [18] P. Li, G. Gao, Y. Wang, Y. Ma, J. Phys. Chem. C 114 (2010) 21745. [19] J. Lv, Y. Wang, L. Zhu, Y. Ma, Phys. Rev. Lett. 106 (2011) 015503. [20] L. Zhu, H. Wang, Y. Wang, J. Lv, Y. Ma, Q. Cui, Y. Ma, G. Zou, Phys. Rev. Lett. 106 (2011) 145501. [21] Y. Ma, Y. Wang, J. Lv, L. Zhu, http://nlshm-lab.jlu.edu.cn/~calypso.html. [22] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [23] G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169. [24] P.E. Blöchl, Phys. Rev. B 50 (1994) 17953. [25] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. [26] S. Baroni, P. Giannozzi, A. Testa, Phys. Rev. Lett. 58 (1987) 1861. [27] P. Giannozzi, S. de Gironcoli, P. Pavone, S. Baroni, Phys. Rev. B 43 (1991) 7231. [28] S. Baroni, A.D. Corso, S.d. Gironcoli, P. Giannozzi, C. Cavazzoni, G. Ballabio, S. Scandolo, G. Chiarotti, P. Focher, A. Pasquarello, K. Laasonen, A. Trave, R. Car, N. Marzari, A. Kokalj, http://www.pwscf.org. [29] A. Šimůnek, J. Vackář, Phys. Rev. Lett. 96 (2006) 085501. [30] X. Guo, L. Li, Z. Liu, D. Yu, J. He, R. Liu, B. Xu, Y. Tian, H.T. Wang, J. Appl. Phys. 104 (2008) 023503. [31] F. Gao, J. He, E. Wu, S. Liu, D. Yu, D. Li, S. Zhang, Y. Tian, Phys. Rev. Lett. 91 (2003) 015502. [32] J. He, E. Wu, H. Wang, R. Liu, Y. Tian, Phys. Rev. Lett. 94 (2005) 015504. [33] K. Li, X. Wang, F. Zhang, D. Xue, Phys. Rev. Lett. 100 (2008) 235504. [34] Y. Zhang, H. Sun, C.F. Chen, Phys. Rev. Lett. 93 (2004) 195504. [35] F.M. Gao, Phys. Rev. B 73 (2006) 132104. [36] Q. Li, H. Wang, Y. Ma, J. Superhard Mater. 32 (2010) 192. [37] J.C. Phillips, Rev. Modern Phys. 42 (1970) 317. [38] P.B. Allen, Phys. Rev. B 6 (1972) 2577. [39] P.B. Allen, R.C. Dynes, Phys. Rev. B 12 (1975) 905. [40] M. Hoesch, T. Fukuda, J. Mizuki, T. Takenouchi, H. Kawarada, J.P. Sutter, S. Tsutsui, A.Q.R. Baron, M. Nagao, Y. Takano, Phys. Rev. B 75 (2007) 140508. [41] J.E. Moussa, M.L. Cohen, Phys. Rev. B 77 (2008) 64518.