Identification of strongly interacting organic semimetals

Dirac- and Weyl point- and line-node semimetals are characterized by a zero band gap with simultaneously vanishing density of states. Given a sufficient interaction strength, such materials can undergo an interaction instability, e.g., into an excitonic insulator phase. Due to generically flat bands, organic crystals represent a promising materials class in this regard. We combine machine learning, density functional theory, and effective models to identify specific example materials. Without taking into account the effect of many-body interactions, we found the organic charge transfer salts (EDT-TTF-I$_2$)$_2$(DDQ)$\cdot($CH$_3$CN) and TSeF-TCNQ and a bis-1,2,3-dithiazolyl radical conductor to exhibit a semimetallic phase in our ab initio calculations. Adding the effect of strong particle-hole interactions for (EDT-TTF-I$_2$)$_2$(DDQ)$\cdot($CH$_3$CN) and TSeF-TCNQ opens an excitonic gap in the order of 60 meV and 100 meV, which is in good agreement with previous experiments on these materials.


INTRODUCTION
Semimetals have attracted huge attention due to their striking transport properties, analogies to high energy physics phenomena, and potential for functionalization [1][2][3]. Their realization relies on a delicate combination of symmetry, electron-filling, and band ordering enforcing the existence of the nodes in the band structure at the chemical potential while having a vanishing density of states (DOS) at the crossing point [4][5][6]. It has been shown extensively for the case of Dirac semimetals, that under a sufficiently high interaction strength a dynamical mass term can be generated leading to a quantum phase transition into a gapped phase [7][8][9][10]. This quantum phase transition strongly depends on: i) the effective fine structure constant α eff , describing the ratio of the coupling of the Fermion field to its gauge field versus the kinetic energy; ii) the dimension of the system; iii) the number of fermionic flavors. Similar phenomena where also discussed in the case of Weyl semimetals [11,12] and line-node semimetals [5,13,14].
With the goal of identifying experimentally feasible materials to investigate interaction effects in nodal semimetals we focus on organic crystals. Organic crystals typically exhibit strong intramolecular forces and weak intermolecular forces, leading to tiny hopping amplitudes for electrons between molecules and resulting flat electronic bands. The flatness corresponds to a tiny quasiparticle kinetic energy and dominant interaction effects. We focus on the excitonic insulator state occurring when a weakly screened Coulomb interaction between a hole and an electron leads to an electron-hole bound state [15,16].
Even though organics seem to be promising materials for strong interaction effects, we face several major difficulties: i) the search space is massive, e.g., the crystallographic open database stores ≈ 200, 000 crystal structures containing carbon and hydrogen [17]; ii) organics are typically large band gap insulators [18]; iii) complex unit cells and strong correlation effects are challenging for an ab initio description. Hence, we apply the following procedure: first, we apply machine learning to narrow down the search space to a computationally feasible set of materials which are predicted to show a tiny band gap; second, using density functional theory we compute the band structures for these materials; third, we select the occurring semimetals, construct effective electronic models and numerically solve a self-consistent Schwinger-Dyson equation to estimate the size of an excitonic gap.

RESULTS
The general workflow of our study can be summarized as follows.
1. Machine learning. We trained a neural network (the continuous-filter convolutional neural network scheme -SchNet [19]) on 24,134 band gaps of nonmagnetic materials taken from the Organic Materials Database -OMDB [18]. We applied the model to 202,117 materials containing carbon and hydrogen stored in the crystallographic open database -COD [20].
2. Band structure calculations. We select 414 materials where the band gap is predicted small, but nonzero (0.01 eV ≤ ∆ ≤ 0.4 eV) and perform medium accuracy ab initio calculations using VASP incorporating the effect of spin-orbit interaction (SOI). Note that all calculations stored within the OMDB were peformed without SOI. Out of the 414 materials we found promising features in the band structures for 9 materials. For these 9 materials we performed high-accuracy VASP calculations taking into account structural optimization and SOI. We found the organic charge transfer salts (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN) and TSeF-TCNQ and a bis-1,2,3-dithiazolyl radical conductor which exhibit a semimetallic phase. Based on symmetry and chemistry we determine the relevant mechanisms to protect the nodal features.
We will discuss the outcome of the three steps in more detail in the following.

Machine Learning
According to the OMDB, organic crystals show a mean ab initio band gap of ≈ 3 eV with a standard deviation of 1 eV [18]. Fig. 1(a) shows a comparison of the band gap distribution of the training set with the band gap distribution obtained using machine learning (ML). While the amount of data is much bigger for the materials taken from the COD, the general shape of the histogram of calculated and predicted gaps agrees well, i.e., the ML model successfully reproduced the band gap statistics. Due to the highly complex structures of organic crystals and the relatively small data set, our trained ML model has a large mean absolute error (MAE) of 0.406 eV. Compared to the average band gap of ≈ 3 eV, this value represents a sufficient accuracy. However, as we are interested in the tiny gap regime, the order of the MAE is similar to our acceptance range 0.01 eV ≤ ∆ ≤ 0.4 eV. Nevertheless, the general guidance of our model was sufficient to identify a few final example materials. We note that a more sophisticated prediction scheme reaching a high accuracy on small and complex data sets would significantly advance the outcome of our approach in the future.
The ABABAB π-stacked bis-1,2,3-dithiazolyl radical was synthesized by Yu et al. [29]. It crystallizes in the non-centrosymmetryc space group P2 1 2 1 2 1 (19) and orders in a canted antiferromagnetic structure with T N ≈ 4.5 K. It was reported to undergo a spin-flop transition to a field-induced ferromagnetic state at 2 K and a magnetic field strength of H ≈ 20 kOe.
For all materials the calculated crossing points occur at the Fermi level with simultaneously vanishing DOS. The selected k-path follows the convention of pymatgen [30]. (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN) shows a crossing of 4 bands at the Brillouin zone center (Γ; Fig. 2(d)) exhibiting a tiny splitting of two two-fold band degeneracies away from Γ due to a small magnetization of two central carbon atoms with ≈ 0.65 µ B . For (TSeF-TCNQ) we tization for all sites within the (TSeF-TCNQ) unit cell. For the π-stacked bis-1,2,3-dithiazolyl radical we find a 2-fold degenerate Weyl node within the ferromagnetic phase along the path ΓX (X = (0.5, 0.0, 0.0), Fig. 2 (f)). A comparison with the metallic band structure of the nonmagnetic and the fully-gapped band structure of an antiferromagnetic phase can be found in the supplementary material.

Protection of the nodes
To understand the nature of the crossings observed in the electronic band structures we are going to discuss the symmetry protection of the nodes. In Fig. 3 we show the molecule resolved partial DOS of (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN). While there are in total four (EDT-TTF-I 2 ) 2 molecules and two DDQ and CH 3 CN molecules in the unit cell, the inversion symmetry present in the crystal enforces pairwise degenerate contributions to the DOS. As can be verified, the main contribution to the DOS around the Fermi energy stems from (EDT-TTF-I 2 ) 2 . Due to the specific stacking structure, molecules in charge transfer salts are known to undergo a transition into a dimerized electronic state, where molecular orbitals of pairs of molecules bind significantly stronger to each other than to other molecules in the crystal [24,31,32]. In other words, we can introduce three energy scales: i) the hopping τ αβ of electrons between atoms α, β within a molecule; ii) the hopping t µν of electrons between molecules µ, ν within a molecular dimer; iii) the hopping s ij of electrons between different dimers i, j. In the limit of s ij t µν τ αβ , the problem can be separated. Assume that the Hamiltonian describing a EDT-TTF-I 2 reveals an electronic ground state with an s-like molecular orbital. Then, the Hamiltonian describing the dimerization can be written asĤ = τΨ † σ xΨ withΨ = (ĉ 1 ,ĉ 2 ), whereĉ † i (ĉ i ) creates (annihilates) an electron in molecule i. Hence, the two eigenstates are given by φ ± = 1 √ 2 (ĉ 1 ±ĉ 2 ) with the respective energies ±τ . In the crystalline unit cell these dimer states can be used to construct a basis where the φ ± are even (+) and odd (-) with respect to inversion symmetry. Furthermore, taking the center of masses of the (EDT-TTF-I 2 ) 2 dimers, the resulting lattice (without (DDQ)·(CH 3 CN)) can be approximated by a reduced lattice with one dimer site per unit cell, where the lattice constant in a direction is decreased by a factor of 1 2 . This decrease of the unit cell also introduces an effective half-filling as the initial space group symmetry is P 1, which according to its associated Bieberbach manifold is half-filled for an odd number of electrons [33]. The half-filling together with the effective bipartite lattice introduces (besides parity and time-reversal) a particle-hole symmetry. Considering a four-band model with two orbital and two spin degrees of freedom we obtain the following effective k · p Hamiltonian around the Γ point (see method section for details), The Hamiltonian reveals a symmetry protected four-fold degenerate Dirac node at the center of the Brillouin zone. By lattice symmetry this node does not carry any topological charge according to the Nielsen-Ninomiya theorem [34,35].
In general, the underlying space group symmetry for (TSeF-TCNQ) (P2 1 /c) does not allow for the high degeneracy of the nodal line observed. The nature of the crossing stems from the quasi one-dimensional nature of the charge transfer salt due to the specific molecular stacking. In Fig. 4 (a) we show the center of mass coordinates of the involved molecules which form weakly interacting one-dimensional chains. Each chain can be approximated to have a dispersion relation of E µσ = s µ 0 + 2s µ 1 cos a 2 · k, with µ being the molecule index distinguishing between TSeF and TCNQ and σ denoting the spin (note that the dispersion relation itself is independent of the spin). As there are two TSeF chains as well as two TCNQ chains present in the crystal, we observe two fourfold degenerate bands (2 chains × 2 spins per molecule) which are allowed to cross in a plane if no hybridization is taken into account (see Fig. 4 (d)). From the tilted stacking of molecules in the different chains (shown in Fig. 4 (b)+(c)) it becomes apparent that the involved hopping of electrons between different chains (interchain coupling) has to be weaker by several order of magnitude than the intrachain coupling (see also fitted parameters in methods section). In the basisΨ = (ĉ TSeF-1σ ,ĉ TCNQ-1σ ,ĉ TSeF-2σ ,ĉ TCNQ-2σ ) T and allowing for a hopping between chains we formulate an effective Hamiltonian of the form where we use ∆ 1 = t 1 cos 1 2 a 1 · k , ∆ 2,3 = t 2,3 cos 1 2 ( a 2 + a 3 ) · k . An example band structure for  Fig. 4 (e) which reproduces the ab initio band structure of Fig. 2 (e) effectively. A slightly more advanced effective Hamiltonian is described in the methods section.
The radical bis-1,2,3-dithiazolyl exhibits Weyl nodes for the low-temperature high-field ferromagnetic phase. Hence, time-reversal symmetry is broken and the crystal reflects a magnetic space group symmetry, depending on the magnetization direction. Assume we choose the magnetization in x-direction then the corresponding group is given by P 2 1 2 1 2 1 (# 19.27). To verify that the observed nodes are Weyl nodes we calculated the atom and orbital resolved weights to the band structure. We observe that the two bands forming the node belong to different atoms in the molecules within the unit cell. The main contributions stem from either the unsaturated nitrogen and the oxygen atom or the saturated nitrogen and one sulfur atom as shown in Fig. 5. Hence, we conclude that the bands belong to different orbital subspaces allowing for a crossing. In symmetry terms the two orbitals correspond to bands with different eigenvalues to the only unitary symmetry element present, (C 2x , 1/2, 1/2, 0), along the invariant line k x . Hence, a 2-band k · p Hamiltonian at a crossing point has to be invariant under the Pauli matrix τ x , H(k x , k y , k z ) = τ x H(k x , −k y , −k z )τ x . This results in an allowed low energy Hamiltonian H( K ± + k) ≈ ±t x k x τ x +t y k y τ y +t z k z τ z , at two points K ± along the x-axis. Both points reveal a monopole charge of ν ± = sign [±t x t y t z ].

Effective Models and Excitonic Gap
To explain the discrepancy between DFT semimetallic phases and experimentally observed semiconductivity in (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN) and (TSeF-TCNQ) we argue that both systems are likely to undergo an excitonic instability. Compared to known inorganic Dirac or Weyl semimetals, the band width of all three materials discussed here is smaller by at least one order of magnitude. As the decreased band width induces a small Dirac velocity (small kinetic energy) of the elementary excitations of the system, we verify that a quasiparticle interaction term becomes dominant. We estimate the size of the effect using effective band structure models. For (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN) we extended the model given in (1) to a lattice periodic version (x → sin(x)). For (TSeF-TCNQ) we construct an eight band model incorporating two spin and four orbital degrees of freedom (two belonging to molecular orbitals from TSeF and two belonging to TCNQ). The allowed dispersion relation of the model is restricted by the symmetry generators: parity, two-fold rotation about x-axis, time-reversal symmetry (details given in the methods section).
To approximate the size of the instability, we solve a BCS-like excitonic gap equation, similar to Refs [36][37][38]. Note that given the numerous degrees of freedom, multiple excitonic gap symmetries might occur, such as inter-or intranode instabilities with or without breaking of time-reversal or other lattice symmetries [38,39]. We focus on an s-wave gap. Details on the calculation are given in the methods section.
We obtained an excitonic gap of ≈ 60 meV for (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN) at T → 0 K. This value agrees with the experimentally determined gap of ≈ 105 meV. We argue that despite (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN) being a regular semiconductor, it is an excitonic insulator with an electronic structure shown in Fig. 6(a). We repeat the approach for (TSeF-TCNQ) and obtain a slightly higher value of ≈ 150 meV ( Fig. 6(b)). In contrast to the Dirac semimetal (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN), (TSeF-TCNQ) has two degenerate line nodes which gap out while the excitonic gap is formed. To verify the experimentally observed metal-insulator transition at 40 K would require a full temperature-dependent screening which is beyond our present study.

Machine Learning
The ML model is trained on a dataset of 24,134 ab initio band gaps of non-magnetic organic crystals stored in the organic materials database -OMDB [18,40]. The dataset is divided into a training, validation and test set of 15000, 3000 and 6134 materials, respectively. We used the continuous-filter convolutional neural network scheme -SchNet [19] with a batch size of 32, a cutoff radius of 5.0Å, 32 features, 3 interaction blocks, a learning rate of 10 −4 , and 50 Gaussians to expand atomic distances (see Ref. [40] for parameter choice). The initial band gap data exhibits a Wigner-Dyson like shape as shown in Fig. 1(a), having a mean of ≈ 2.9 eV and a standard deviation of ≈ 1.1 eV. Our trained ML model shows a mean absolute error (MAE) of 0.406 eV and a RMSE of 0.602 eV which is interpreted as an accuracy of ≈ 90%. Note that the underlying data is rather small (≈ 2 × 10 4 materials) and highly complex (average of 85 atoms per unit cell). The trained ML model is able to reproduce the band gap statistics on a set of 202,117 organic crystals (≈ 2 × 10 4 materials) taken from the COD [20], as shown in Fig. 1(a). These predictions also fit a Wigner-Dyson distribution ∼ x 5.62 e −0.45x 2 . The performance of our model on the test set is shown in Fig.  1(b).

ab initio Calculations
We performed ab initio calculations in the framework of the density functional theory (DFT) using a projector augmented-wave method as implemented in the Vienna ab initio simmulation package VASP [41]. Structures are taken from the COD [20] and transformed into VASP input using pymatgen [30]. During the self-consistent calculation of the electron density and the DOS we used a Γ-centered mesh with a k-mesh density of 800 points perÅ −3 for the quick materials scan of 414 candidate materials and a k-mesh density of 1500 points perÅ −3 for the refined calculations of 9 prospective semimetals. For the latter we performed an optimization of the atomic positions using a conjugate gradient algorithm. All calculations were performed with SOI and the exchange correlation functional was approximated by the strongly constrained and appropriately normed semilocal density functional SCAN [42] including van-der-Waalscorrections using the Tkatchenko-Scheffler method with iterative Hirshfeld partitioning [43,44]. The cut-off energy of the plane wave expansion was chosen to by 600 eV.

Effective Hamiltonians
We generated effective Hamiltonians for (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN) and TSeF-TCNQ. To describe the four-fold degeneracy at the Γ point observed for (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN), we construct a model with two orbital (τ i ) and two spin (σ i ) degrees of freedom, Here, each molecule in the unit cell contributes one molecular orbital. Neglecting higher and lower energy bands, a four-fold degeneracy is obtained assuming an even-and an odd-parity molecular orbital, and add respective lattice symmetries (space group P1): paritŷ P = τ 3 σ 0 ( k → − k), time-reversal −iσ 2K ( k → − k; K being the complex conjugation), and an emergent particle-hole symmetryĈ = σ 2 . Up to linear order we obtain The best fit to the ab initio calculated band structure resulted in the parameters a 1 = 0.47 eV, a 2 = −0.33 eV, a 3 = 0.11 eV, b 1 = 0.16 eV, b 2 = −0.38 eV, b 3 = 0.15 eV.
The nodes for TSeF-TCNQ occur in the interior of the Brillouin zone. Therefore, we construct a latticeperiodic model with 8 bands (4 orbital-, 2 spin degrees of freedom). The orbital degrees of freedom are obtained from respective permutations of the two TSeF and two TCNQ molecules within the primitive unit cell. We represent the generators of the factor group G/T C 2h (G is the space group P2 1 /c, T is the corresponding group of pure lattice translations) as follows,P = P(1, 2, 3, 4) × σ 0 ( k → − k),Ĉ 2x = P(2, 1, 4, 3) × iσ 1 (k x → k x , k y → −k y , k z → −k z ) and add time-reversal symmetry as before.
Here P(permutation) denotes a 4×4-dimensional permutation matrix. Note that each symmetry element comes with a corresponding double group partner. The total Hamiltonian is written as H = H orbital × H spin . We generate H orbital in the basis (|Ψ TSeF1 , |Ψ TSeF2 , |Ψ TCNQ1 , |Ψ TCNQ2 ) using

Excitonic Instability
We model quasiparticle-quasihole excitations of the material in terms of the following interaction process, where incoming (outgoing) solid lines belong to fermionic annihilation (creation) operatorsΨ q (Ψ † q ) and the dashed line to a screened scalar interaction V ( q).
The Hamiltonian is given bŷ (6) Here, ν is a corresponding quantum number and a = ± (a = −a) denotes electron and hole states.Ĥ 0 denotes the non-interacting effective Hamiltonian of the system. We derive an excitonic gap equation using a Green function approach similar to Refs [36][37][38] and define G ν,± ( p, τ − τ ) = − TΨ ν, p± (τ )Ψ † ν, p± (τ ) and F ν,± ( p, τ − τ ) = − TΨ ν, p± (τ )Ψ † ν, p∓ (τ ) . We proceed by calculating the equations of motion for G and F using the Heisenberg formalism and decomposing the fourpoint functions in terms of G and F [46]. We impose particle-hole symmetry of the electronic band structure ξ + ν ( k) = ξ − ν ( k) as well as a real excitonic gap function ∆ ν and derive the BCS-like gap equation which is in agreement with Ref. [36]. Here, E ν, p = ξ 2 ν ( p) + ∆ 2 ν ( p). We construct the exact Brillouin zone using the algorithm of Finney [47] and solve (7) using a Liouville-Neumann series, where the integration is performed using a quasi-Monte-Carlo integration scheme [48]. We calculate the zero temperature excitonic gap tanh βE ν, q 2 → 1 and approximate the interaction by a screened Coulomb interaction V ( q) = − 4π q 2 + k 2 , with k being a screening vector generally temperature temperature dependent [36]. The screening vector is approximated by the Thomas-Fermi screening, which for Dirac semimetals is given by k = 2gα π k F [38], with k F being the Fermivector and g being the Dirac cone degeneracy. Organic crystals typically exhibit purification during the crystallization process expelling dopants. We assume a clean sample with Fermi-level deviation δµ = 0.01 eV.

CONCLUSION
Typically, nodal states within the electronic structure of organic crystals are inferred by indirect measurements, e.g., via the resistivity [49], optical conductivity [50], and local spin susceptibility [51] as has been done extensively for the pressure induced semimetal phase of (BEDT-TTF) 2 I 3 . A direct observation of the electronic structure of organic materials with angle-resolved photoemission spectroscopy (ARPES) is difficult due to: usually tiny crystal sizes limiting signal; insulating behavior leading to charging of the samples; limits in available orientations and problems in preparing well-defined surface terminations. However, ARPES has been performed for organics in selected cases [52]. Here, the reported experimental maximum/median/minimum crystal sizes (mm) of the materials described by us are 0.24/0.05/0.04 for (EDT-TTF-I 2 ) 2 (DDQ)·(CH 3 CN), 0.23/0.05/0.02/ for (TSeF-TCNQ), and 0.34/0.04/0.02 for bis-1,2,3-dithiazolyl. While these crystal sizes are tiny, recent progress in the field has decreased the focus area of state-of-the-art ARPES devices tremendously [53], making a measurement of the photo electron spectrum of the three materials possible in the near future. The materials reported here provide a natural platform to investigate the rare excitonic insulator phase as a consequence of strong interaction effects within symmetry protected nodal semimetals.
While several inorganic nodal semimetals are known, the space of organics remains fairly unexplored in this regard. The realm of 3D organic semimetals is mainly composed of the single example α-(BEDT-TTF) 2 I 3 and modifications which exhibit tilted Dirac nodes under pressure (≈ 2.3 GPa) [51,54] or chemical strain [55]. Hence, the outcome of our ML and ab initio calculations are promising with respect to the identification of novel examples. We note that the choice of the s-wave excitonic instability due to strong interactions is only the simplest scenario and was chosen to estimate the size of the effect. However, other forms of gap openings and richer phase diagrams are imaginable. The ABABAB π-stacked bis-1,2,3-dithiazolyl radical was reported to undergo a transition into a canted antiferromagnet below a critical temperature of 4.5 K and to undergo a spin-flop transition to a field-induced ferromagnetic state at a temperature of 2 K and a magnetic field strength of H ≈ 20 kOe [29]. Our calculations revealed Weyl nodes along the path ΓX in the Brillouin zone for the ferromagnetic phase. To better understand the emergence of the phase and its correspondence to the underlying magnetic ordering we performed similar band structure calculations incorporating full structural optimization for an antiferromagnetic and a non-magnetic ordering (details on the ab initio calculation can be found in the methods section). Even though the nonmagnetic phase exhibits topological nodes along the path ΓY the overall behavior is metallic as there is no vanishing density of states at the Fermi level. The topological protection of the nodes for the nonmagnetic phase is a consequence of the underlying orthorhombic symmetry with space group P2 1 2 1 2 1 where the mechanism is described in Ref. [56]. The antiferromagnetic phase exhibits a clear band gap of the size of ≈ 0.15 eV. The comparison of the electronic structure for all three magnetic phases is shown in Fig. 8. The material (EDT-TTF-I)2(TCNQ) with COD-ID 4506562 exhibits topological nodes within the Brillouin zone close to the Fermi level as shown in Fig. 9. However, the material is not a semimetal as it exhibits a large density of states at the Fermi level.

Appendix D: List of predictions
We trained a machine learning model based on the continuous-filter convolutional neural network scheme SchNet [19] on 24,134 ab initio calculated band gaps stored within the organic materials database -OMDB [18] (details can be found in the method section). We applied the successfully trained model on 202,117 crystal structures stored within the crystallographic open database -COD [20]. All crystal structures considered belong to organic molecular crystals or metal organic frameworks which were synthesized before. Organic materials tend to be large band gap insulators and only 414 materials were predicted to have a band gap of 0.01 ≤ ∆ ≤ 0.4 eV, where we explicitly tried to exclude organic metals. The list of predicted band gaps for the 414 materials denoted with their COD-ID is given in Table I.  TABLE I. COD-IDs, band gap ∆, and number of atoms in the primitive unit cell Nat for the subset of 414 COD materials where the band gap was predicted to be small (0.01 ≤ ∆ ≤ 0.4 eV).
COD ID ∆ (eV) N at COD ID ∆ (eV) N at COD ID ∆ (eV) N at COD ID ∆ (eV) N at COD ID ∆ (eV) N at COD ID ∆ (eV) N at