Introduction

Since the discovery of topological insulators1,2,3,4,5,6,7, topological band theory has been widely applied to solid state materials8,9,10,11 and leads to the further classification of both insulators12,13,14,15,16,17,18 and semimetals19,20,21,22,23,24,25,26,27,28,29 based on the quantum topological invariants. The development of the topological semimetals stimulated extensive research in understanding a variety of topological properties. Very recently, some comprehensive databases for the symmetry protected nonmagnetic topological states were established by three independent works30,31,32. These useful databases contain not only the topological insulators, but also the semimetals with degeneracy points located at high-symmetric lines or points. However, in these databases, an important category of topological materials, e.g., the Weyl semimetal (WSM) whose degeneracy point can locate at any position in the entire Brillouin zone, is still absent.

Generally speaking, topological semimetals are the systems whose conduction and valence bands cross with each other in the Brillouin zone and cannot be removed by perturbations of Hamiltonian without breaking specific symmetries. Topological WSM is such a material with double degeneracy Weyl points. Compared with other topological states, the WSM is extraordinary since the topology is not protected by any crystalline symmetry. It can only appear in the system whose spin degeneracy is removed by breaking time reversal or spacial inversion symmetry. In this case, the bands that create the Weyl points take the 2 × 2 Weyl equation form, \(H=\pm {v}_{F}\overrightarrow{k}\cdot \overrightarrow{\sigma }\), which can be viewed as half of the massless Dirac equation with defined chirality in three-dimensional space33. Because a two-bands model Hamiltonian can be described by a linear combination of three Pauli matrices generally, the linear band crossing can be obtained by tuning the parameters without considering any additional conditions. Therefore, besides time reversal or spacial inversion symmetry, the lattice periodic boundary condition is the unique requirement for Weyl points. Accordingly, the identification of topology based on the characteristics of high-symmetry points does not work anymore for WSMs, and that is why the WSMs are still absent in the comprehensive databases30,31,32.

In this work, we only focus on the experimental structures to complete the databases by including nonmagnetic WSMs. The total information about their topological properties is listed in the Supplementary Information. At the same time, two effective tools used for generating Wannier function automatically and searching Weyl points are developed in Fig. 1. The accuracy of them is essential for finding materials with observable properties of Weyl points, Some of them have already been theoretically predicted and experimentally confirmed, such as the family of TaAs23,24,25,26. Besides type-I and type-II WSMs, nine chiral WSMs with inversion and mirror symmetries breaking are also listed in our database. For chiral WSMs, the truly quantized circular photogalvanic effect (CPGE)34,35,36,37 could be induced by the energy difference of Weyl points with opposite chirality. Moreover, we explored the bulk photovoltaic effect (BPVE)38 of all these noncentrosymmetric WSMs. A giant BPVE which reach up to 700 μA/V2 is predicted, and some of them can even extend to the regime of visible light.

Fig. 1: Workflow of the high-throughput calculations for the WSMs identification.
figure 1

a The experimental noncentrosymmetric crystal structures are loaded from the Inorganic Crystal Structure Database (ICSD) after removing alloys and repeated structures. All the compounds containing strongly correlated f-electrons are also removed from the list. With the final input list, we perform a nonmagnetic calculation to obtain the band gap and the charge carrier density. Here we label the compounds with magnetic elements V, Mn, Fe, Co, and Ni or magnetic moments above 0.05 μB per unit cell as magnetic. Furthermore, for the remaining 232 semimetal phases, the Bloch wave functions are automatically projected into high-symmetry atomic-orbital-like Wannier functions and generate the corresponding tight-binding (TB) model Hamiltonians. Using these TB model Hamiltonians, we search the Weyl points around the Fermi energy. b Sub-branch of automatic Wannier function generation. c Sub-branch of automatic Weyl points search with the Berry phase approach method.

Results

Our automatic workflow for searching nonmagnetic WSMs is shown in Fig. 1. It includes two crucial algorithms, the Wannier function generation and Weyl points search. Constrained by the accuracy of the density functional theory (DFT), we don’t consider the systems with strongly correlated f-elections, namely Lanthanide and Actinide series except element La, while beginning our work with the experimental noncentrosymmetric crystal structures in the Inorganic Crystal Structure Database (ICSD)39. Here we remain some systems with 3d-elections, which may also be strongly correlated systems, to make sure that we don’t omit some potential candidates. After removing alloys and repeated structures, totally 8896 compounds left. Among them, all the systems with the total magnetic moments exceeding 0.05 μB per unit cell and magnetic elements V, Mn, Fe, Co, and Ni are regarded as magnetic compounds. Furthermore, the gap and the charge carrier density (N) of the nonmagnetic phases are checked. Because we are only interested in WSMs with observable topological Weyl-points-related properties, the insulators and good metals (with N > 2 × 1022/cm3, which is obtained by the comparison with Co2VGa40) will not be considered in the further calculations.

After all the screening processes above, we obtain the list of the nonmagnetic noncentrosymmetric experimental crystal structures who may possess Weyl points near Fermi level. For these compounds, the atomic-orbital-like Wannier functions are constructed by an effective algorithm automatically as shown in Fig. 1b. The crucial problem to sole for this step is the setting of the fitting windows with appropriate projection orbital basis sets. In the nonmagnetic first-principles calculations, we select s orbitals of alkaline-earth metals, d orbitals of transition metals and s + p orbitals of p-block elements as the minimum basis sets. Thus, we can label the maximum weight value among the selected band as r. The weight represent the ratio of selected basis projection. If the weight r is smaller than 0.8, the minimum basis sets will be insufficient to produce the DFT band structure; hence, we switch the bases to the maximum basis sets which include s + p + d orbitals of all the elements. With the selected basis sets, the inner projection energy window can be easily fixed as 90% of the maximum weight value r. Subsequently, we loop over the outer projection energy window and check the mean absolute error between the DFT and Wannier band structures around the Fermi level until the error is smaller than 0.02 eV. Meanwhile, it is noteworthy that a good projection can be always guaranteed because of the large band gap between valence and core states.

Using the high-symmetry tight-binding (TB) model Hamiltonians obtained in the construction of Wannier functions, the Weyl points can be scanned. Since the existence of Weyl points in noncentrosymmetric nonmagnetic system only requires the translational symmetry, the fact that the Weyl points would locate at generic k-points without any special symmetry renders the search difficult. The typical way to search for the Weyl points use the band gap as a criterion. However, this method is much easy to reach the local minimum value, thereby losing some Weyl points. In this work, we use the Berry phase of the parallelepiped corresponding to the k-mesh of the reciprocal space as the criterion. Initially, based on the TB model Hamiltonian, we genera;te a relatively sparse 51 × 51 × 51 k-mesh along the three reciprocal vectors to calculate the Berry phase of each parallelepiped. If the value is larger than 0.01 × 2π, Weyl points might be present inside the parallelepiped. For this type of parallelepiped, we increase the k-mesh finer and repeat calculating Berry phase. Eventually, the Berry phase will increase in some fine parallelepiped and finally reach a value n × 2π, where n is an integer, corresponding to a single Weyl point. In our calculations, the parallelepiped is refined to around 1/16003 of the reciprocal lattice to obtain the Berry phase with n ≥ 0.99.

After searching the Weyl points, we have a raw list of the WSMs with all the information of Weyl points. It may also have some magnetic materials such as antiferromagnets. Thus, we process a double check here to modify the raw list, and the list of nonmagnetic WSMs are obtained.

Despite having as many as 8896 noncentrosymmetric candidates, we only find 46 nonmagnetic WSMs with Weyl points near the Fermi level within 300 meV, including nine chiral crystal structures. This result express that the good nonmagnetic WSMs are difficult to find. The details of every WSMs that we find by our algorithm are given in Supplementary Table 1. Until now, the Weyl points have been classified into type-I and type-II according to the shape of the linear band crossings. The Fermi surface would shrink to a singularity point for the type-I Weyl point, while the type-II Weyl point presents a Fermi surface with a linear crossing which is composed of electron and hole41. On the other hand, another main difference is the precondition of chiral anomaly. It can exist along any direction in the type-I WSM as long as electric and magnetic fields are not perpendicular to each other, but only exist when the electric and magnetic fields point to some selected direction in the type-II WSM.

Discussion

Representative materials

Here, we choose two typical candidates covering each class of WSM as examples. For type-I WSM, we highlight an orthorhombic sulfide TaAgS3 with space group Cmc21 (No. 36) as shown in Fig. 2a. The energy dispersion calculated by DFT is given in Fig. 2d. It presents a semimetallic structure with small electron and hole pockets located around Γ and Y points, respectively. With inversion symmetry breaking, the spin degeneracy is lifted. Finally, two pairs of Weyl points, a minimum number of Weyl points in time-reversal-symmetric systems, would be generated. Because of the C2z rotation and the time reversal symmetry, the Weyl points are constrained in the kz = 0 plane as shown in Fig. 2b. All the Weyl points are related with each other and located at the same energy, only ~30 meV below the charge neutral point, owing to the existence of the mirror planes mx and my, On the other hand, low charge carrier density can make the Weyl points observable by transport measurements. The Weyl points are further confirmed by the (001) surface state calculation. When these Weyl points are projected into (001) surface, the Weyl points with different chirality are projected into different position. If the chemical potential is set at the energy of the Weyl points, the Fermi arcs connecting the Weyl points can be observed as shown in Fig. 2e. The Fermi arc length can be ~15% of the reciprocal lattice vector, which makes it easy to detect by surface detection methods, such as angle-resolved photoemission spectroscopy (ARPES) and scanning tunneling microscopy (STM).

Fig. 2: Type-I WSM TaAgS3.
figure 2

a Orthorhombic crystal structure of TaAgS3. b Brillouin zone of TaAgS3 and Weyl points distribution which is obtain from Weyl point search algorithm. c Energy dispersion along one pair of Weyl points parallel to the kx axis. d Band structure with (red dot lines) and without (black solid lines) SOC along high-symmetry lines. The band inversion appears around Γ point. e Surface Fermi arcs with energy fixed at the Weyl points. The purple and blue dots represent the Weyl points with positive and negative chirality, respectively.

The other example Ag2S is a hybrid WSM42 containing both two kinds of Weyl points. It also has an orthorhombic crystal structure with space group P212121 (No. 19), as shown in Fig. 3a. Different from TaAgS3, without considering spin-orbital coupling (SOC), the band inversion between the p orbital of S and the d orbital of Ag forms two Dirac points on the high-symmetric lines ΓX and ΓY which is protected by the C2 rotation symmetry30,31,32. Furthermore, turning on SOC leads to the separation of the Dirac points and four pairs of Weyl nodes appear, since SOC breaks the spin rotation symmetry. The locations of Weyl points in the Brillouin zone shown in Fig. 3b connected with each other by time reversal, C2x and C2y rotation symmetries. Through the band structure along the k-path crossing one pair of Weyl points in Fig. 3d, we find these Weyl points belong to two different classes, the type-I Weyl points with chirality +1 (W1) and the type-II Weyl points with chirality −1 (W2). They are distinguished by different color in Fig. 3. At the same time, the type-II Weyl points W2 can also be located by the linear-crossing-like Fermi surface in Fig. 3e if we put the chemical potential at W2.

Fig. 3: Hybrid WSM Ag2S with chiral crystal symmetry.
figure 3

a Orthorhombic crystal structure of Ag2S. b Brillouin zone and Weyl points distribution which is obtain from Weyl point search algorithm. c Band structure with (red dot lines) and without (black solid lines) SOC along high-symmetry lines. The obvious band crossing points are located at ΓX and ΓY without SOC. d Energy dispersion along one pair of Weyl points with opposite chirality. W1 and W2 represent the type-I Weyl point and the type-II Weyl point, respectively. e The constant energy surface in the energy of the type-II Weyl point W2. the linear crossing point of electron and hole pockets is the exact location of W2. f, g (001) surface states with energy fixed at W1 and W2. The purple and blue dots represent the type-I Weyl points with positive chirality and type-II Weyl points with negative chirality, respectively.

One point we need to emphasized is that Ag2S does not present any mirror symmetry. Although the mirror symmetry is not necessary for the existence of Weyl points, the chirality of Weyl points must be odd with mirror operations. Therefore, the mirror symmetry forces one pair of Weyl points with opposite chirality locate at the same energy. In another words, WSM without mirror symmetry exhibits a net chirality if the Fermi energy is located between one pair of Weyl points with opposite chirality. One of the most important properties for mirror-symmetry-broken WSMs is the quantized CPGE, which is the only possible quantized signal induced by the monopole Weyl point so far. This can also be applied to Ag2S with W1 ~23 meV above the charge neutral point and W2 ~28 meV below the charge neutral point, as illustrated in Fig. 3d. Because the charge carrier density is only around 5 × 1020/cm3 and the contribution of the trivial Fermi surfaces can be neglected, Ag2S may be an ideal candidate to observe the quantized CPGE.

The chiral WSM Ag2S is further confirmed by the chiral surface Fermi arcs. Fixing the energy at W1 and W2, we can obtain the extremely long chiral Fermi arcs in (001) surface in Fig. 3f, g, respectively. Owing to the large spin splitting, the spin texture should be easily observed by both ARPES and STM. Since two W2 Weyl points overlap with each other when they are projected into (001) surface, two Fermi arcs start from same projected Weyl point. These extremely clear topological properties offer a good platform for both surface detection and bulk nonlinear optical transports.

Shift current analysis

Very recently, a giant shift current was observed experimentally in WSM TaAs43, which suggests a topological effect in BPVE and the practical applications of WSMs in optical detectors and solar energy conversion. Inspired by this work, we screened the nonlinear optical response for all the nonmagnetic WSMs in our database. Half of them can host a large shift current above 100 μA/V2 with photon energy larger than 0.1 eV (Supplementary Table 2 and Supplementary Figs. 13). Our database shows that huge shift current can be obtained in both infrared and visible spectrum in noncentrosymmetryic WSMs, and it provides exotic material candidates for the study of interplay of different functions based on WSMs.

Through perturbative analysis, the expression of shift current44 is written as

$${J}_{i}(0)=2{\sigma }_{ijk}(0;\omega ,-\omega ){Re}[{E}_{j}(\omega ){E}_{k}(-\omega )],$$
(1)

where i, j, k = x, y, z. The second order optical conductivity tensor can be computed as

$${\sigma }_{ijk}(0,\omega ,-\omega )=\frac{i{e}^{3}\pi }{\hslash V}\sum _{\overrightarrow{k}}\sum _{m,n}{f}_{nm}\left({r}_{mn}^{j}{r}_{nm}^{k;i}+{r}_{mn}^{k}{r}_{nm}^{j;i}\right)\delta \left(\hslash \omega -{E}_{mn}\right),$$
(2)

where n and m are band indexes, fnm = fnfm is the Fermi-Dirac distribution difference, \({r}_{nm}^{i}={A}_{nm}^{i}(1-{\delta }_{n,m})\) is the inter-band Berry connection, \({r}_{nm}^{i;j}={\partial }_{{k}_{j}}{r}_{nm}^{i}-i({A}_{nn}^{j}-{A}_{mm}^{j}){r}_{nm}^{i}\) depends on intermediate virtual states, and Enm = EnEm is the energy difference between two bands. The delta function \(\delta \left(\hslash \omega -{E}_{mn}\right)\) is disposed by a Gauss distribution with smearing factor 0.01 eV. To avoid the numerical divergence with zero denominator, we replace the transition-energy-related term \(\frac{1}{{E}_{mn}}\) by \(\frac{{E}_{mn}}{{E}_{mn}^{2}+{\alpha }^{2}}\) with α = 0.04 eV, which is stable for the shift current with photon energy above 0.1 eV. After a k-grid convergence test, we find that the change is less than 5% with k-grid increasing from 100 × 100 × 100 to 200 × 200 × 200. Here, we use a grid 200 × 200 × 200 in our calculations. The calculated tensors of shift current are double-checked by the symmetry. For convenient, we set the tensor at each k-point as

$${\sigma }_{ijk}^{\prime}(\overrightarrow{k};0,\omega ,-\omega )=\frac{i{e}^{3}\pi }{\hslash V}\sum _{m,n}{f}_{nm}\left({r}_{mn}^{j}{r}_{nm}^{k;i}+{r}_{mn}^{k}{r}_{nm}^{j;i}\right)\delta \left(\hslash \omega -{E}_{mn}\right).$$
(3)

In the experimental measurement of WSM TaAs (ref. 43), an apparent measured value σxxz = 34 ± 3.7 μAV−2 with photon energy ћω = 0.117 eV is obtained. Our calculation gives a well-fitting theoretical value σxxz = 27 μAV−2 with the same photon energy. Furthermore, ref. 43 added the reflectance of TaAs at 10.6 μm to the apparent measured value manually and the giant shift current σxxz = 154 ± 17 μAV−2 is obtained in the end. At the same time, σzzz has the biggest value for TaAs in both our calculation and ref. 43.

The largest shift current in WSMs was found in BiTeI with space group P3, in which the component of σxxz can reach up to 750 μA/V2 at a photon energy 0.66 eV. It is around five times larger than that in TaAs. Moreover, four independent tensor components in BiTeI are above 200 μA/V2 in a large photon energy window of 0.3–0.7 eV (see Fig. 4a), which allows large shift current in different experimental setup. Hence the new WSM BiTeI offers a more ideal candidate for the generation of strong BPVE. Considering the layered structure of BiTeI stacking along z, this component provides a natural advantage to integrate the material into devices.

Fig. 4: Shift current tensor σijk.
figure 4

It is a function of photon energy for a BiTeI, b LaIrP, c LaRhP and d SrHgSn. They are calculated based on TB model Hamiltonians generated by step b of our algorithm. All of them have large shift current with tensor σijk > 100 μA/V2 and huge range of photon energy which can even reach visible light (>1.65 eV).

Beside the largest shift current candidate BiTeI, we also find another two kind of candidates. One is multifunctional candidates combining topology, strong shift current, and superconductivity. As indicated in Fig. 4b, c, the shift current tensor can reach up to σzzz = 564 μA/V2 and σzzz = 250 μA/V2 for the new found WSMs LaIrP and LaRhP, respectively. In addition, since both LaIrP and LaRhP are superconductors45 with transition temperatures Tc = 5.3 and 2.5 K, these two materials will provide good platform for the study of the interplay among topology, superconductivity and nonlinear optical responses. The other one is the candidates with large photon energy. For all these WSMs with large shift current in Supplementary Table 2 and Supplementary Figs. 13, we find that the peak value of σijk locates at an energy below 1 eV, except SrHgSn, CaBiAg and CaCdPb. Fig. 4d shows the shift current tensor of SrHgSn with space group P63mc (No. 186) as an example. One can find a peak value σzzz = 165 μA/V2 at photon energy ~1.8 eV, which is already located at the range of red light.

The shift current can be considered as a joint effect of the frequency-dependent density of states (DOS) and the wave function. Taking SrHgSn as an example, we can try to understand the microscopic mechanism in these two aspects. From the DOS of SrHgSn shown in Fig. 5b, one can find two peaks locating bellow and above the charge neutral point, respectively. The energy difference is around 1.8 eV, which is very close to the two peaks of σzzz in Fig. 4d. \({\sigma }_{zzz}^{\prime}(\overrightarrow{k};0,\omega ,-\omega )\) is also calculated for each k-point with the photon energy 1.8 and 1.95 eV near the two peaks. Thus we obtained that the positive contribution is mainly from the special plane kz ~0.312π/c at photon energy around 1.8 eV, while kz ~0.132π/c plane at photon energy around 1.95 eV has a negative contribution. The distribution of \({\sigma }_{zzz}^{\prime}(\overrightarrow{k};0,\omega ,-\omega )\) at photon energy 1.8 eV is dominated by twelve positive blue lines caused by the excitation from band n − 1 to bands n + 3 and n + 4, as shown in Fig. 5c, f. The distribution of \({\sigma }_{zzz}^{\prime}(\overrightarrow{k};0,\omega ,-\omega )\) changes sharply when photon energy increasing from 1.8 to 1.95 eV. At photon energy 1.95 eV, \({\sigma }_{zzz}^{\prime}(\overrightarrow{k};0,\omega ,-\omega )\) is dominated by two negative red rings as shown in Fig. 5d, e. These two red rings are derived by the light excitation from two occupied bands n and n − 1 to the non-occupied bands n + 7, n + 8, n + 9 and n + 10, as presented in Fig. 5g, h. Therefore, the two peaks of σzzz(0, ω, − ω) in Fig. 4d originate from large DOS located around EF − 1.2 eV and EF + 0.6 eV, and the sign changes because of the different band excitation at different k-points. One point need to emphasize is that the strong shift current at high photon energy is not caused by the Weyl points here.

Fig. 5: The mechanism analysis for \({\sigma }_{zzz}^{\prime}(\overrightarrow{k};0,\omega ,-\omega )\) of SrHgSn.
figure 5

a Energy dispersion along high symmetry k-path. b Density of states (DOS). ce Distribution of \({\sigma }_{zzz}^{\prime}(\overrightarrow{k};0,\omega ,-\omega )\) in reciprocal space contributed by band n − 1 in kz = 0.312π/c plane with ω = 1.8 eV, band n − 1 in kz = 0.132π/c plane with ω = 1.95 eV, and band n in kz = 0.132π/c plane with ω = 1.95 eV, respectively. fh Light excitation from band n − 1 to bands n + 3 and n + 4 in kz = 0.312π/c plane with ω = 1.8 eV, from band n − 1 to bands n + 7, n + 8, n + 9 and n + 10 in kz = 0.132π/c plane with ω = 1.95 eV, and from band n to bands n + 7, n + 8, n + 9 and n + 10 in kz = 0.132π/c plane with ω = 1.95 eV, respectively. n is the number of electrons in the primitive cell. The blue and red arrows in fh represent positive and negative contributions to σzzz.

In summary, through high-throughput calculations, we have complemented the nonmagnetic topological material databases by including the nonmagnetic WSMs. Among them, nine chiral structures, such as Ag2S, which allow the quantized CPGE and have natural long surface Fermi arcs are also recognized. Furthermore, most of them can have a large shift current, especially for BiTeI. Beside the two databases of nonmagnetic WSMs and their nonlinear optical responses, the effective algorithms for searching Weyl points by Berry phase and automatically Wannier function generation provide powerful tools for the study of magnetic topological materials.

Methods

The calculation starts from the experimental noncentrosymmetric crystal structures in the Inorganic Crystal Structure Database (ICSD, http://www2.fiz-karlsruhe.de/icsd_home.html). We take these crystal structures as the input of the full-potential local-orbital (FPLO) minimum-basis code46. Here we choose the generalized gradient approximation of the Perdew-Burke-Ernzerhof functional47 as the exchange-correlation potential. A dense 12 × 12 × 12 k-mesh is used in the static nonmagnetic calculation. SOC is also included self-consistently during the workflow. In order to get the surface states of WSMs, Green’s function method48,49 is applied with a half-infinite boundary condition. During the DFT calculations, we assume that the ground state is nonmagnetic after excluding the materials with magnetic elements and large total magnetic moment. The final database which may contains antiferromagnets is also double-checked to make sure they are indeed nonmagnetic.