Exploration of stable compounds, crystal structures, and superconductivity in the Be-H system
Shuyin Yu, Qingfeng Zeng, Artem R. Oganov, Chaohao Hu, Gilles Frapper, Litong Zhang
vversion 0
Exploration of stable compounds, crystal structures, andsuperconductivity in the Be-H system
Shuyin Yu, Qingfeng Zeng, ∗ Artem R. Oganov,
2, 3, 4
Chaohao Hu, Gilles Frapper, and Litong Zhang Science and Technology on Thermostructural Composite Materials Laboratory,School of Materials Science and Engineering,Northwestern Polytechnical University, Xi’an, Shaanxi 710072, PR China Department of Geosciences, Center for Materials by Design,and Institute for Advanced Computational Science,State University of New York, Stony Brook, NY 11794-2100, USA Moscow Institute of Physics and Technology,Dolgoprudny, Moscow Region 141700, Russia School of Materials Science and Engineering,Northwestern Polytechnical University, Xi’an, Shaanxi 710072, PR China School of Materials Science and Engineering,Guilin University of Electronic Technology,Guilin 541004, People’s Republic of China IC2MP UMR 7285, Universit´e de Poitiers - CNRS, Poitiers 86022, France a r X i v : . [ c ond - m a t . m t r l - s c i ] J u l bstract Using first-principles variable-composition evolutionary methodology, we explored the high-pressure structures of beryllium hydrides between 0 and 400 GPa. We found that BeH remainsthe only stable compound in this pressure range. The pressure-induced transformations are pre-dicted as Ibam → P ¯3 m → R ¯3 m → Cmcm → P /nmm , which occur at 24, 139, 204 and 349GPa, respectively. P ¯3 m R ¯3 m structures are layered polytypes based on close packings of Hatoms with Be atoms filling octahedral voids in alternating layers. Cmcm and P /nmm containtwo-dimensional triangular networks with each layer forming a kinked slab in the ab -plane. P ¯3 m R ¯3 m are semiconductors while Cmcm and P /nmm are metallic. We have explored super-conductivity of both metal phases, and found large electron-phonon coupling parameters of λ =0.63for Cmcm with a T c of 32.1-44.1 K at 250 GPa and λ =0.65 for P /nmm with a T c of 46.1-62.4 Kat 400 GPa. The dependence of T c on pressure indicates that T c initially increases to a maximumof 45.1 K for Cmcm at 275 GPa and 97.0 K for P /nmm at 365 GPa, and then decreases withincreasing pressure for both phases. . INTRODUCTION The search for new high-temperature superconductors has attracted great enthusiasmin both fundamental and applied research. Owing to its low mass and high electron den-sity, “metallic hydrogen” has been predicted to possess a high superconducting transitiontemperature ( T c >
200 K) . However, hydrogen remains insulating at extremely highpressure ( >
320 GPa ), which are too high for any applications. Another feasible methodof obtaining the properties of metallic hydrogen is to form hydrogen-rich alloys with otherelements . Due to “chemical precompression”, the pressure of metallization may be reducedsignificantly.Inspired by the elusive state of matter, theoretical and experimental research has madeconsiderable progress towards exploring superconductivity in hydrogen-rich compounds, e.g.for group IVa hydrides, calculations predicted that SiH , GeH , SnH and PbH may become superconductors at high (yet lower than pure H) pressure. The origin of high-pressure superconductivity can be derived from the particular “H ” units, which are a featurecommon to hydrides of alkali metals , alkaline earth metals and group IVa elements .Experiments suggested metallization of SiH at ∼
60 GPa and its superconducting transi-tion temperature ( T c ) is 17 K at 96 and 120 GPa , though debates remain. In addition, thesuperconductivity of group IIIa hydrides (BH , AlH and GaH ) and alkaline earthmetal hydrides (CaH n , SrH n and BaH n ) have also been extensively explored.Beryllium hydrides can be an interesting subject of study, because low atomic mass ofBe may lead to very high T c values. The only known beryllium hydride is BeH . Theground-state structure of BeH is body-centered orthorhombic with Ibam symmetry. Atambient conditions, BeH is an insulator with a pronounced band gap of 5.5 eV . Vajeestonet al. proposed that BeH undergoes a series of phase transitions α → β → γ → δ → (cid:15) atpressures of 7, 51, 87 and 98 GPa, respectively, and reported that BeH remains insulatingup to 100 GPa. Zhang et al. systematically investigated the pressure-induced metallizationof alkaline earth hydrides, and found the metallization pressure of Pnma -BeH to be greaterthan 300 GPa. Wang et al. predicted that BeH reaches a metallic state by a R ¯3 m → Cmcm phase transition, instead of a direct band gap closure in R ¯3 m phase.3 I. COMPUTATIONAL METHODOLOGY
First-principles variable-composition evolutionary simulations were performed at 0, 50,100, 150, 200, 250, 300 and 400 GPa using the USPEX code , which has the capability ofdiscovering possible stoichiometries and the corresponding stable and metastable structuresat given pressure-temperature conditions, and has successfully predicted a large number ofstable structures . The initial generation of structures and compositions was producedrandomly with the use of space groups picked randomly from the total list of 230 groups.50% of the lowest-enthalpy structures were used as parents for the next generation. Inaddition, 20% of structures in each new generation were produced by lattice mutation, 15%by atomic transmutation and 15% were produced randomly. Each generation contained 60structures and runs proceeded for up to 50 generations.The underlying structure relaxations were carried out using the Vienna
Ab-initio
Sim-ulation Package (VASP) code , in the framework of density functional theory (DFT) within the Perdew Burke Ernzerhof generalized gradient approximation (PBE-GGA) . Thefrozen all-electron projected augmented wave approach (PAW) was adopted to describe thecore electrons and their effects on valence orbitals. A plane-wave kinetic energy cutoff of 600eV and dense Monkhorst-Pack k -point grids with a resolution higher than 2 π × A − were used for all structures. The most stable structures were studied further at increasedaccuracy using a reciprocal-space grid better than 2 π × A − .Phonon calculations were carried out using the supercell approach as implemented in thePHONOPY code . Electron-phonon coupling (EPC) calculations were explored using thepseudopotential plane-wave method within PBE-GGA, as implemented in the Quantum-Espresso package . In these calculations, we used the kinetic energy cutoff of 60 Ry andMonkhorst-Pack k -point grids of 20 × ×
12 for the
Cmcm phase and 16 × × P4/nmm phase with a Methfessel-Paxton smearing factor of 0.05 Ry. Additionally, q -meshes of 5 × × Cmcm and 4 × × P4/nmm were used to calculate the electron-phonon coupling matrix elements, respectively. We used the Allen-Dynes-modified McMillanequation to estimate T c , as follows: T c = ω log . (cid:34) − . λ ) λ − µ ∗ (1 + 0 . λ ) (cid:35) (1)where ω log is the logarithmic average frequency, λ is the electron-phonon coupling constantand µ ∗ is the Coulomb pseudopotential, which is assumed to be between 0.10-0.13 .4 II. RESULTS AND DISCUSSIONS
FIG. 1. (Color online). Convex hull phase diagrams for the Be-H system at 0, 50, 150, 250 and400 GPa.
Fig. 1 shows the convex hull phase diagrams for the Be-H system at selected pres-sures. The ground-state enthalpy of formation ∆H f is defined as ∆H f (Be x H y )=∆H(Be x H y )- x∆H(Be) - y∆H(H). A compound is thermodynamically stable if it has lower enthalpy thanany isochemical mixture of the elements or other compounds. Such stable compounds formthe convex hull. Based on our evolutionary searches, elemental Be adopts the P6 /mmc structure below 390 GPa and bcc Im ¯3 m structure above 390 GPa. Our findings are in goodagreement with previous calculations . Hydrogen undergoes a series of phase transitions: P6 /m (P <
105 GPa),
C2/c (105 < P <
207 GPa),
Cmca-12 (270 < P <
385 GPa),
Cmca (P >
385 GPa) , in addition to the experimentally known Ibam structure, we found aseries of pressure-induced structural transformations ( Ibam → P ¯3 m1 → R ¯3 m → Cmcm → P4/nmm ) with increasing pressure. Notably, we did not find other stable compounds be-sides BeH over the entire pressure range 0 - 400 GPa. The detailed structural parametersof these predicted phases are summarized in Tab. I.At ambient conditions, BeH crystallizes in the orthorhombic Ibam structure. This struc-ture consists of a three-dimensional network of distorted tetrahedra with Be atoms sittingat the center of the tetrahedra and H atoms at the corner in a bridged position between two5
IG. 2. (Color online). Extended crystal structures of solid BeH for (a) the P ¯3 m1 structureat 50 GPa; (b) the R ¯3 m structure at 150 GPa; (c) the Cmcm structure at 250 GPa; and (d) the
P4/nmm structure at 400 GPa. The large blue spheres represent Be atoms, while the small redand green spheres indicate two symmetrically inequivalent H atoms.
Be atoms. The orthorhombic phase transforms to a CdI -type structure ( P ¯3 m1 ; Fig. 2a) at24 GPa, then it transforms to a related CdCl -type structure ( R ¯3 m ; Fig. 2b) at 139 GPa(Fig. 3a). Both structures are made of layers of edge-sharing BeH -octahedra, but stackingsequences of these layers are different. The shortest interlayer H-H distances decrease from1.83 ˚A at 50 GPa to 1.51 ˚A at 150 GPa (Fig. 3b).We found a similar high-pressure structure in the B-H system : at P >
168 GPa, the
P6/mmm -BH structure is the most stable phase, and may be described as stacking BH-layerswith planar closely-packed arrays of boron atoms. On-top H atoms locate symmetricallybetween the boron layers. Note that if one assigns a formal charge of -1 to H (hydride-like),both the boron atoms in BH and beryllium atoms in BeH have a formal ns valence electronconfiguration. Notably, the structural transformation from Ibam to P ¯3 m1 is accompaniedby a large density jump of 9.81% while only 0.79% increase occurs at the transition from P ¯3 m1 to R ¯3 m (Fig. 3b). The structures in ref. are metastable with respect to the P ¯3 m1 structure.At 204 GPa, the layered R ¯3 m structure transforms into an orthorhombic Cmcm structure6
ABLE I. Optimized structural parameters for the predicted BeH structures at selected pressures.Pressure Space group No. Lattice parameters Atom Wyckoff positions(GPa) (˚A, deg) Sites x y z50 P ¯3 m1 a = 2.085 Be c = 3.104 H R ¯3 m a = 1.886 Be c = 8.316 H Cmcm a = 1.796 Be b = 5.503 H1 c = 2.840 H2 P4/nmm a = 1.906 Be c = 3.240 H1 (Fig. 2c). In this structure, Be atoms are coordinated by eight hydrogens, whereas hydro-gens are in the fourfold coordination (H1 atoms - in planar square coordination, H2 atoms- in tetrahedral coordination). Note that H1 atoms form flat pure-hydrogen rectangularlayers with shortest H-H distance of 1.42 ˚A at 250 GPa. The tetragonal P4/nmm structure(PbClF-type) becomes stable at 349 GPa. In this structure, Be atoms are coordinated bynine hydrogen atoms (forming a capped tetragonal antiprism); H1 atoms are in a fourfold(planar square) coordination and H2 atoms in a fivefold (square pyramid) coordination.The structure can be viewed as layered, with double layers formed by Be and H2 atoms,alternating with square layers formed by H1 atoms. The Be-Be distances at 400 GPa are1.901 and 1.906 ˚A within the double layers, and 2.328 ˚A between these layers, reinforcingcohesion of the highly delocalized covalent three-dimensional BeH structure. The shortestH-H distance is between H1 atoms, i.e. in the pure-hydrogen square layer - at 400 GPa thisdistance is 1.345 ˚A. At such distances overlap of atomic orbitals is strong enough to makethe material metallic. This is very similar to the Pbcn -SiH with the closest H-H distanceof 1.35 ˚A.Fig. 4 shows electronic densities of states (DOS) for the P ¯3 m1 and R ¯3 m structures, fromwhich it can be clearly seen that both structures are semiconductors with band gaps of7 IG. 3. (Color online). (a) Enthalpy per atom for various BeH structures as a function ofpressure with the P ¯3 m1 structure taken as the reference; (b) Computed equations of state of BeH (solid lines) and shortest H-H distance (dotted lines). p states and H- s states, suggesting large degree of covalency. The covalent bonds aremainly from the intralayer BeH octahedra while the interlayer interactions are mainly vander Waals forces. Even at very high pressures these layered structures remain insulating(e.g., the R ¯3 m structure has band gap 0.18 eV at 200 GPa).Fig. 5 shows the band structures, partial DOS and electron localization function (ELF)for the Cmcm phase at 250 GPa and
P4/nmm phase at 400 GPa, respectively. The bandstructures reveal that both structures are metallic with several bands crossing the Fermi level,and a pseudogap. The dispersed valence and conduction bands near the Fermi level signifya relatively large DOS at the Fermi level (0.098 and 0.107 electrons/eV/f.u., respectively),which may favor superconducting behavior. The valence band widths are greater than in thelow-pressure phases, which indicates enhanced electron delocalization, thus, more electronsparticipate in bonding interactions, which promotes structural stability.Distributions of the electron localization function (ELF) reveal electron accumulation onthe H atoms. For the
P4/nmm phase, the ELF between H atoms in the square H-layer isclose to 0.5 (Fig. 5e), equal to the value for the electron gas. The DOS and ELF analysis arealso in agreement with Bader analysis. In the Cmcm phase, Be atoms have charge +1.58,while the charges of H1 and H2 atoms are -0.82 and -0.76, respectively. In the
P4/nmm phase, the charges of Be, H1 and H2 atoms are +1.57, -0.94 and -0.63, respectively.8
IG. 4. (Color online). Total and Partial density of states (DOS) for the (a) P ¯3 m1 phase at 50GPa and (b) R ¯3 m phase at 150 GPa. The calculated phonon spectra for the
Cmcm and
P4/nmm phases establish dynamicalstability, as there are no imaginary phonon frequencies anywhere in the Brillouin zone (Fig.6). We further explored the superconductivity for both structures by performing electron-phonon coupling (EPC) calculations. For the
Cmcm structure at 250 GPa, the electron-phonon coupling parameter λ is 0.63, indicative of quite strong EPC. Using the calculated ω log of 1670.7 K and the commonly accepted values of the Coulomb pseudopotential µ ∗ (0.1-0.13) , we obtained T c in the range of 32.1-44.1 K using the modified Allen-Dynes-modifiedMcMillan equation . Our results are similar to those of Wang et al. , the differences beingdue to different pseudopotentials. From Fig. 7c, it is clear that T c of the Cmcm structurefirst increases and then decreases with increasing pressure, and reaches a maximum (45.1K) at ∼
275 GPa.Fig. 7 shows the total and partial phonon density of states together with the Eliash-berg phonon spectral function α F( ω ) and electron-phonon integral λ ( ω ) as a function offrequency for the P4/nmm structure at 400 GPa. Low-frequency ( <
35 THz) vibrations are9
IG. 5. (Color online). (a) and (b) band structures and partial DOS for
Cmcm phase at 250GPa and
P4/nmm phase at 400 GPa, respectively; (c), (d) and (e) electron localization function(ELF) through specific surfaces.FIG. 6. (Color online). Calculated phonon spectra along selected high symmetry points for (a)the
Cmcm phase at 250 GPa and (b)
P4/nmm phase at 400 GPa. IG. 7. (Color online). (a) Total and partial phonon density of states (PDOS) for
P4/nmm phase;(b) Eliashberg phonon spectral function α F( ω ) and electron-phonon integral λ ( ω ) as a function offrequency at 400 GPa; and (c) calculated superconducting transition temperature ( T c ) vs pressurefor Cmcm and
P4/nmm phases, triangles and squares represent µ ∗ =0.1 and 0.13, respectively. mostly related to Be atoms, while higher-frequency ( >
45 THz) modes mainly come fromthe vibration of H1 and H2 atoms. At 400 GPa, The calculated EPC parameter λ is 0.65,indicating rather strong EPC in the P4/nmm structure. Using the calculated ω log of 2170K and µ ∗ (0.1-0.13), we obtained T c in the range of 46.1-62.4 K. The vibrations of Be below35 THz contribute about 44.3% of total λ , while the vibrations of H1 and H2 above 45 THzcontribute about 55.7% with no obvious difference between H1 and H2 vibrations to λ . Inaddition, the pressure dependence of T c displays the same trend as observed in the Cmcm phase and reaches a maximum of 97.0 K at 365 GPa. This is one of the highest T c valuespredicted in literature. Note that the Allen-Dynes formula is expected to be reliable whenis less than 1-1.5 , which is the case here. IV. CONCLUSIONS
In summary, using variable-composition evolutionary simulations for crystal structureprediction, we investigated the high-pressure phases of solid beryllium hydrides in the pres-sure range of 0-400 GPa. BeH is found to be the only stable beryllium hydride. Thepressure-induced transformations are predicted to be Ibam → P ¯3 m1 → R ¯3 m → Cmcm → P4/nmm , which occur at 24, 139, 204 and 349 GPa, respectively. The layered P ¯3 m1 and R ¯3 m structures belong to the well-known CdI and CdCl types, respectively. The Cmcm
P4/nmm phases contain 8- and 9-coordinate Be atoms, respectively, and layers of Hatoms with short H-H distances, responsible for metallic conductivity. The entire phasetransformations are first-order with volume shrinkage values of 9.81%, 0.79%, 2.71% and0.43%, respectively. The P ¯3 m1 and R ¯3 m structures are semiconductors while the Cmcm and
P4/nmm phases are metallic. Electron-phonon coupling calculations show that the
Cmcm and
P4/nmm structures are phonon-mediated superconductors, with large electron-phonon coupling parameters of 0.63 for the
Cmcm phase with a T c of 32.1-44.1 K at 250GPa and 0.65 for the P4/nmm phase with a T c of 46.1-62.4 K at 400 GPa. Dependence of T c on pressure indicates that T c will increase initially to a maximum value of 45.1 K for the Cmcm phase at 275 GPa and 97.0 K for the
P4/nmm phase at 365 GPa, respectively, andthen decrease with increasing pressure for both structures.
ACKNOWLEDGMENTS
We thank the Natural Science Foundation of China (Grants No. 51372203, No. 11164005and No. 51332004), the National Basic Research Program of China (973 Program, GrantNo. 2014CB643703), the Basic Research Foundation of NWPU (Grant No. JCY20130114),the Foreign Talents Introduction and Academic Exchange Program (Grant No. B08040), theNational Science Foundation (Grants No. EAR-1114313 and No. DMR-1231586), DARPA(Grants No. W31P4Q1310005 and No. W31P4Q1210008), and the Government of the Rus-sian Federation (Grant No. 14.A12.31.0003) for financial support. The authors also acknowl-edge the High Performance Computing Center of NWPU, Shanghai Supercomputer Centre,and the National Supercomputing Center in Shenzhen and the GENCI-CINES (France) forthe allocation of computing time on their machines. ∗ [email protected] N. W. Ashcroft, Physical Review Letters , 1748 (1968). T. Barbee, A. Garc´ıa, and M. L. Cohen, Nature , 369 (1989). P. Cudazzo, G. Profeta, A. Sanna, A. Floris, A. Continenza, S. Massidda, and E. Gross,Physical Review Letters , 257001 (2008). P. Loubeyre, F. Occelli, and R. LeToullec, Nature , 613 (2002). N. Ashcroft, Physical Review Letters , 187002 (2004). J. Feng, W. Grochala, T. Jaro´n, R. Hoffmann, A. Bergara, and N. Ashcroft, Physical ReviewLetters , 017006 (2006). Y. Yao, J. Tse, Y. Ma, and K. Tanaka, EPL (Europhysics Letters) , 37003 (2007). M. Martinez Canales, A. R. Oganov, Y. Ma, Y. Yan, A. O. Lyakhov, and A. Bergara, PhysicalReview Letters , 087005 (2009). G. Gao, A. R. Oganov, A. Bergara, M. Martinez Canales, T. Cui, T. Iitaka, Y. Ma, and G. Zou,Physical Review Letters , 107002 (2008). Z. Li, W. Yu, and C. Jin, Solid State Communications , 353 (2007). J. Tse, Y. Yao, and K. Tanaka, Physical Review Letters , 117004 (2007). G. Gao, A. R. Oganov, P. Li, Z. Li, H. Wang, T. Cui, Y. Ma, A. Bergara, A. O. Lyakhov,T. Iitaka, et al. , Proceedings of the National Academy of Sciences , 1317 (2010). P. Zaleski Ejgierd, R. Hoffmann, and N. Ashcroft, Physical Review Letters , 037002 (2011). E. Zurek, R. Hoffmann, N. Ashcroft, A. R. Oganov, and A. O. Lyakhov, Proceedings of theNational Academy of Sciences , 17640 (2009). J. Hooper, B. Altintas, A. Shamp, and E. Zurek, The Journal of Physical Chemistry C ,2982 (2013). J. Hooper, T. Terpstra, A. Shamp, and E. Zurek, The Journal of Physical Chemistry C (2014). X. J. Chen, V. V. Struzhkin, Y. Song, A. F. Goncharov, M. Ahart, Z. Liu, H. K. Mao, andR. J. Hemley, Proceedings of the National Academy of Sciences , 20 (2008). M. Eremets, I. Trojan, S. Medvedev, J. Tse, and Y. Yao, Science , 1506 (2008). C. Hu, A. R. Oganov, Q. Zhu, G. R. Qian, G. Frapper, A. O. Lyakhov, and H. Y. Zhou,Physical Review Letters , 165504 (2013). I. Goncharenko, M. Eremets, M. Hanfland, J. Tse, M. Amboage, Y. Yao, and I. Trojan, PhysicalReview Letters , 045504 (2008). A. Islam, M. Ali, and M. Ali, Physica C: Superconductivity , 403 (2010). G. Gao, H. Wang, A. Bergara, Y. Li, G. Liu, and Y. Ma, Physical Review B , 064118 (2011). H. Wang, S. T. John, K. Tanaka, T. Iitaka, and Y. Ma, Proceedings of the National Academyof Sciences , 6463 (2012). G. S. Smith, Q. C. Johnson, D. K. Smith, D. Cox, R. L. Snyder, R.-S. Zhou, and A. Zalkin,Solid State Communications , 491 (1988). B. T. Wang, P. Zhang, H. L. Shi, B. Sun, and W. D. Li, The European Physical Journal B ,303 (2010). P. Vajeeston, P. Ravindran, A. Kjekshus, and H. Fjellvag, Applied physics Letters , 34(2004). C. Zhang, X. J. Chen, R. Q. Zhang, and H. Q. Lin, The Journal of Physical Chemistry C ,14614 (2010). Z. Wang, Y. Yao, L. Zhu, H. Liu, T. Iitaka, H. Wang, and Y. Ma, The Journal of ChemicalPhysics , 124707 (2014). A. R. Oganov and C. W. Glass, The Journal of chemical physics , 244704 (2006). A. R. Oganov, A. O. Lyakhov, and M. Valle, Accounts of chemical research , 227 (2011). A. O. Lyakhov, A. R. Oganov, H. T. Stokes, and Q. Zhu, Computer Physics Communications , 1172 (2013). A. R. Oganov, Y. Ma, A. O. Lyakhov, M. Valle, and C. Gatti, Reviews in Mineralogy andGeochemistry , 271 (2010). Q. Zeng, J. Peng, A. R. Oganov, Q. Zhu, C. Xie, X. Zhang, D. Dong, L. Zhang, and L. Cheng,Physical Review B , 214107 (2013). W. Zhang, A. R. Oganov, A. F. Goncharov, Q. Zhu, S. E. Boulfelfel, A. O. Lyakhov, E. Stavrou,M. Somayazulu, V. B. Prakapenka, and Z. Konˆopkov´a, Science , 1502 (2013). Q. Zhu, A. R. Oganov, and A. O. Lyakhov, Physical Chemistry Chemical Physics , 7696(2013). G. Kresse and J. Furthm¨uller, Physical Review B , 11169 (1996). P. Giannozzi, S. De Gironcoli, P. Pavone, and S. Baroni, Physical Review B , 7231 (1991). S. Baroni, P. Giannozzi, and A. Testa, Physical Review Letters , 1861 (1987). J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters , 3865 (1996). P. E. Bl¨ochl, Physical Review B , 17953 (1994). H. J. Monkhorst and J. D. Pack, Physical Review B , 5188 (1976). A. Togo, F. Oba, and I. Tanaka, Physical Review B , 134106 (2008). P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L.Chiarotti, M. Cococcioni, I. Dabo, et al. , Journal of Physics: Condensed Matter , 395502(2009). M. Methfessel and A. Paxton, Physical Review B , 3616 (1989). P. B. Allen and R. Dynes, Physical Review B , 905 (1975). Z. C. Guo, F. Luo, and Y. Cheng, Computational Materials Science , 139 (2014). C. J. Pickard and R. J. Needs, Nature Physics , 473 (2007). R. F. Bader, “Atoms in molecules: a quantum theory. international series of monographs onchemistry 22,” (1990). P. B. Allen and B. Mitrovi´c, Solid State Physics , 1 (1983)., 1 (1983).