Introduction

The fundamental character and technological relevance of group-IVa elements (tetrels) have motivated repeated investigations and systematics on their polymorphism1,2,3,4,5. Carbon polymorphs are promising as hard and transparent materials6,7. Semiconducting silicon (Si) is versatile for micro- and nanoelectronic devices, while germanium (Ge) displays comparatively higher carrier mobility, a finer band gap tunability and good compatibility with high-dielectric constant materials8. Metallization occurs in silicon and germanium upon compression1. In Ge lowering of phonon frequencies promotes electron-phonon coupling towards superconductivity2,3. The possibility of metallic germanium under room conditions is very intriguing and intensively debated2,9, while superconductivity in elemental Ge appears under pressure3. In the lower pressure range, improved optical properties via band-gap tuning can be achieved in a different polymorph.

Engineering viable new compounds with superior properties entails a detailed understanding of structural changes10. Under pressure germanium bears similarities with silicon1 by comparatively higher transition pressures with respect to Si, due to Ge core d-electrons11. Upon compression semiconducting Ge (cubic diamond) transforms into β-tin type (space group I41/amd) at about 10 GPa12 and then to Imma phase13, simple hexagonal (P6/mmm)14, followed by orthorhombic Cmca phase15 and finally upon further compression above 180 GPa, by the hexagonal close-packed arrangement (P63/mmc)15.

The phase diagram of germanium is further complicated by a family of tetrahedral structures16,17. Type-II clathrate Ge(cF136) exists at ambient conditions17. Other germanium modifications are reported, Ge(tP12), Ge(cI16) (γ-silicon type, BC8) and Ge(hR8)16. BC8 Ge18 is accessible through decompression from β-tin Ge. In a nutshell, during the last few years, new dense and open phases of germanium have been experimentally observed or synthesized. Nevertheless, a systematic approach to including known and finding novel germanium forms is still outstanding and of top priority in order to explain recent experiments2.

In this work we explore the energy landscape of germanium, both at ambient conditions and upon moderate compression. By means of ab initio metadynamics, we efficiently sampled structural transformations along appropriate collective reaction coordinates. Besides the already known dense phases of germanium, we found two novel allotropes (green in Fig. 1). The first one is a monoclinic modification of germanium (mC16) slightly less dense than diamond, is an indirect band gap semiconductor and is unprecedented for tetrel elements. The second one is a five-coordinated (square pyramidal) metallic intermediate structure (tI4, sp. gr. I4/mmm), which incurs in the diamond (cF8) → β-tin phase (tI4) transition and which has been postulated to exist in homologue silicon (bct-5)19. We report on transformation paths, energetic, mechanical and electronic properties. For the metallic bct-5 phase, we calculated superconducting temperatures down to ambient pressure compared to β-tin phase, based on the electron-phonon coupling mechanism20.

Figure 1
figure 1

Lower pressure region of the phase diagram of Ge, augmented by two novel phases mC16 and bct-5, found by ab initio metadynamics runs.

bct-5 shows characteristic square pyramidal 5-fold coordination of Ge atoms. In monoclinic mC16 four-rings are a characteristic feature. The arrows indicate the direction of metadynamics evolution. The pressures were evaluated based on the common tangent construction (see below, Fig. 2).

Full size image

Results

The mC16 structure (Fig. 1, C2/m, a = 7.6094 Å, b = 7.9746 Å, c = 6.5668 Å, β = 104.10°) arised from a metadynamics run started from diamond (8 atoms box, p = 1 bar, T = 300 K). Ge atoms occupy three Wyckoff positions: (4i) 0.70984 0.50000 0.67434, (4i) 0.60981 0.50000 0.29080, (8j) 0.65012 0.76388 0.11596. Strikingly, mC16 is less dense than diamond (see Fig. 2), although topologically as dense as diamond or lonsdaleite. Its bulk modulus amounts to 51.2 GPa that is slightly lower than that of the diamond type-structure (60.7 GPa), estimated from the fit to the third order Birch-Murnaghan equation of state. Applying pressure to the mC16 allotrope in an additional metadynamics run resulted into a direct transition to the β-tin phase. Decompressing the latter is known to generate metastable phases typically denser than diamond Ge. Therefore, a viable route to mC16, like for other recent germanium allotropes, could rather be the oxidation of suitable germanium Zintl salt precursors, i.e. via chemical synthesis.

Figure 2
figure 2

Equation of states of Ge mC16 and bct-5, compared to cubic diamond and β-tin type.

bct-5 features a reduced volume per atom compared to diamond type (cF8), while the total energy minimum lies lower than β-tin. mC16 on the contrary is less dense and energetically close to the diamond type. The tangents used for evaluating equilibrium pressures of Fig. 1 are highlighted.

Full size image

Upon compression diamond transforms into β-tin and it subsequently follows the same transition sequence of silicon phases. Along the diamond → β-tin transition metadynamics (64 atoms box, p = 10 GPa, T = 300 K, Gaussian width δs = 4.5 (GPa Å3)(1/2), Gaussian height W = δs2) visited an intermediate of bct-5 topology (I4/mmm, a = 3.5491, c = 6.4478, Ge(4e) 0.0 0.0 0.19273). The bct-5 bulk modulus is 58.7 GPa, slightly lower than that of the β-tin phase (68.2 GPa). This five-connected structure has been proposed for silicon19, but has never been observed so far.

The total energy/volume curves of Fig. 2 suggest bct5 as a conventional product of diamond compression (here “conventional” denotes a positive equilibrium pressure, whereby an “unconventional” product requires a “negative” pressure). However, under hydrostatic conditions β-tin is formed from diamond, while decompression leads to other germaniums, although indications of minority phases exist. Low-temperature nonequlibrium decompression further influences phase formation towards amorphous21 and novel germaniums22.

The transformation affects only one box parameter, suggesting nonhydrostatic shearing as the protocol of choice towards bct-5, also supported by the magic-stress approach which led to bct-5 in silicon19. Alternatively, low-temperature compression may be considered. The evolution of the enthalpy profile of metadynamics runs from diamond Ge (Fig. 3) shows bct-5 as a narrow plateau around metastep 160. Metadynamics runs started from bct-5 (Fig. 3, black curve) confirms it as a proper intermediate along the transition, which can be quenched down to room pressure and which is mechanically stable (see below). Mechanistically, the coordination number increases from 4 to 5 on shortening one bond, followed by flattening of the pristine tetrahedron and formation of the square pyramidal geometry of bct-5, Fig. 3b. The four bonds in the pyramid basis are 2.62 Å long, the axial one 2.48 Å.

Figure 3
figure 3

Metadynamics (DFT) runs on a 64 Ge atoms box (p1 = 10.0 GPa, T = 300 K).

bct-5 appears as a stepwise feature in the enthalpy profile of the run started from diamond (red line). Runs commenced from bct-5 (black line) evolve into β-tin. Configurations corresponding to distinct points along the runs are detailed below the graph. Energy units are eV.

Full size image

The evolution of the bonding situation from Ge diamond to β-tin over bct-5 is shown in Fig. 4. Calculation of the ELF25 shows four, one + four and two + four bond attractors, respectively. The five “bonds” in this orbital-deficient, electron-deficient metallic bct-5 result from the sp Ge valence shell. This bonding scenario is reminiscent of the recently discovered superconducting Zintl phase CaGe326, isosymmetric with bct-5.

Figure 4
figure 4

Evolution of the bonding situation from diamond to β-tin type over bct-5.

The ELF map is showing four bond attractors for diamond Ge (a, η = 0.58), one + four bond attractors for bct-5 (b, η = 0.53, transparent green isosurface η = 0.48), two + four bond attractors for β-tin (c, η = 0.51).

Full size image

The electronic band structures and phonon dispersions of mC16 and bct-5 are shown in Fig. 5. The tetrahedral mC16 phase is semiconducting with an indirect band gap of 1.43 eV (PBE-GGA27), while bct-5 is metallic and stable down to 0 GPa. Isothermic-Isobaric molecular dynamics (1 bar, 300 K, 2.5 ps) confirmed the stability of bct-5. mC16 is characteristic due to the presence of four-membered rings. However, this does not imply overall structure destabilization28. The expectation of a strained geometry is in fact not supported by total energy calculations, which place mC16 among the energetically lowest Ge allotropes. The indirect band gap and the low density (compared to the diamond type) make this germanium an attractive material. The need for a “negative” pressure makes a chemical path plausible.

Figure 5
figure 5

Band structures and phonon spectra23 (0.0 GPa) of mC16 (left) and bct-5 (right).

mC16 is a narrow-gap semiconductor (band gap = 1.43 eV), while bct-5 is metallic. Both are mechanically stable. Brillouin zone choice (bct-5) according to Ref. 24.

Full size image

Discussion

The prominent property of bct-5 is the conservation of metallic character down to ambient conditions. Calculations and experiments have shown an increase of the superconducting temperature on lowering pressure, with superconductivity still present at 6.9 GPa (Tc ≈ 6.0 K)2, outside of the existence range of β-tin. This points to the existence of another superconducting phase. The evolution of Tc as a function of pressure for bct-5 and β-tin is shown in Fig. 6. The calculated value of Tc for bct-5 at 6.9 GPa is Tc = 6.1 K (β-tin Tc = 7.7 K). The ongoing debate on the feasibility of “strange” metals other than β-tin2,9 and particularly the need for further explanations of the survival upon decompression of “metallic” metastable phases (not reliably identified as any known Ge phase up to now) open the door for a serious consideration of the role of the bct-5 in the lower pressure range, as a metallic state that qualifies for higher Tc values (via the McMillan relation).

Figure 6
figure 6

Evolution of ω, λ and Tc as a function of pressure for bct-5 and β-tin, calculated based on the electron-phonon coupling model.

The calculated equilibrium pressure between bct-5 and β-tin marks the boundary between the phases. Notice the flattening of bct-5 Tc after 5 GPa as phonons soften. The model predicts an increase of Tc in the lower pressure region.

Full size image

In conclusion, with metadynamics runs we have found and characterized two novel Ge polymorphs. The first one, mC16, is an indirect gap semiconductor and its structure is unprecedented for the tetrel family. The second bct-5 polymorph was suggested for Si. Our simulations lean strong relevance to bct-5 in the lower pressure range, as a new metallic superconducting phase capable of stability at room conditions. We expect our predictions to stimulate further experimental work.

Methods

Metadynamics

Different approaches are successful in the prediction of novel polymorphs of the elements29,30. Important discoveries have been achieved by means of random techniques, genetic (evolutionary) algorithms, or accelerated molecular dynamics10. Metadynamics31,32,33 allows for the exploration of the energy surface along one or more collective reaction coordinates. The method is independent of the level of theory used, it does not require prior knowledge of the energy landscape and its sampling efficiency can be enhanced by parallel runs started from different configurations. The time-evolution of the system is biased by a history-dependent potential, which discourages the system from visiting already harvested regions of the potential34. Efficiency is achieved in metadynamics also through dimensionality reduction. Instead of studying the problem in the full 3N dimensional configuration space of N particles, a relatively small number of collective coordinates s = (s1..sm) is used instead, which provide a coarse-grained description of the system and are able to distinguish between different free energy minima, i.e. different phases. The inclusion in the space of collective variables of slow degrees of freedom positively impacts the performance of the method.

All metadynamics runs were performed with at least eight atoms in the simulation box which served as a collective (6-dimensional) variable. This minimal box approach was successfully used in the prediction of novel carbon polymorphs, recently7. The size of the minimal box ensured commensurability of all already known phases (except for the clathrate II phase, which requires a minimum of 34 atoms) either open or dense. Each metadynamics metastep consisted of a molecular dynamics run in the NVT ensemble for 0.5 ps (timestep 2 fs) at 300 K.

Density functional computational layers

Metadynamics was performed with different molecular dynamics layers. A Density Functional Tight Binding (DFTB)35 level of theory, as implemented in the Γ-point-only DFTB module of the CP2K code36,37, ensured rapid and accurate sampling in the low-pressure regime, characterized by four-connected Ge atoms. For higher-pressures SIESTA38 was used as the DFT molecular dynamics layer, allowing for k-point runs. Electronic states were expanded by a single-ζ basis set constituted of numerical orbitals with a norm-conserving Troullier-Martins39 pseudopotential description of the core electrons. Single-ζ basis set dramatically reduces computational times providing nonetheless, the right topology and energy differences of all the Ge allotropes under study. The charge density was represented on a real-space grid with an energy cutoff38 of 200 Ry. A Monkhorst-Pack k-point mesh of 2 × 2 × 2 ensured the convergence of the electronic part. High-pressure metadynamics was performed based on DFT. Lower pressure regions were initially explored by DFTB, followed by DFT metadynamics upon discovery of interesting novel polymorphs. In the lower pressure range the transferability between DFTB and DFT is unflawed. To validate the method a set of preliminary test runs were performed. Similarly to metadynamics applied to silicon40, the correct sequence of reported phases3 of Ge could be reproduced (not shown).

Band structures, phonons and electron-phonon coupling

Electronic structure, phonon dispersion curves and superconducting properties were calculated with the Quantum Espresso (QE)20,23 package. The superconducting critical temperature Tc was evaluated based on the Allen and Dynes modification of the McMillan formula. This required calculating the electron-phonon coupling strength λ via the Eliashberg function. The Coulomb potential value was μ = 0.1. A q-mesh of 12 × 12 × 12 was used for the evaluation of the dynamical matrix, while a Monkhorst-Pack k-point mesh of 16 × 16 × 16 ensured convergence of the electronic part.

Structure characterization

The structures visited during each run were characterized by their vertex symbols, which contain the information on all the shortest rings meeting at each atom and coordination sequences, as implemented in the TOPOS package41. In case of new structures ideal space group and asymmetric units were identified with the Gavrog Systre package42. Subsequently a variable-cell geometry optimization was performed (DFT-GGA, PBE functional27) in a plane-wave pseudopotential framework23 using Vanderbilt ultrasoft pseudopotential with non-linear core correction43. A k-point mesh of 8 × 8 × 8 ensured convergence of the electronic part, while a plane-wave basis set with an energy cut-off of 30 Ry was applied.

Electron localization function

The electron localization function (ELF)25 is a widely used tool to study chemical bonding in molecules and solids. The ELF monitors the correlation of the movement of parallel spin electrons in real space44. In this study the ELF calculations were performed by using the ELF module as implemented45 in the all-electron, full-potential local orbital (FPLO) minimal basis method46. In FPLO each atomic orbital nl with principal quantum number n and angular momentum l is represented by one basis function only. The basis functions are obtained by solving an effective Schrödinger equation which contains the spherically averaged crystal potential and an artificial confining potential47. The confining potential makes the basis functions more localized than the atomic orbitals. The FPLO method does not have any atomic (or muffintin) spheres so that the whole space is treated in a uniform manner. In the scalar relativistic calculations within the LDA scheme (Perdew and Wang48) Ge(3d, 4s, 4p, 4d) represented the basis sets. Lower-lying states were treated as core states.