High-throughput first-principle prediction of collinear magnetic topological materials

The success of topological band theory and symmetry-based topological classification significantly advances our understanding of the Berry phase. Based on the critical concept of topological obstruction, efficient theoretical frameworks, including topological quantum chemistry and symmetry indicator theory, were developed, making a massive characterization of real materials possible. However, the classification of magnetic materials often involves the complexity of their unknown magnetic structures, which are often hard to know from experiments, thus, hindering the topological classification. In this paper, we design a high-throughput workflow to classify magnetic topological materials by automating the search for collinear magnetic structures and the characterization of their topological natures. We computed 1049 chosen transition-metal compounds (TMCs) without oxygen and identified 64 topological insulators and 53 semimetals, which become 73 and 26 when U correction is further considered. Due to the lack of magnetic structure information from experiments, our high-throughput predictions provide insightful reference results and make the step toward a complete diagnosis of magnetic topological materials.

Utilizing group theory, band degeneracy and connectivity at/ between different maximal K-points can be precisely determined, which lays the theoretical basis for automating the topological classification of real material. Before the comprehensive understanding of the relationship between topological invariant and group symmetry, there already exists a simplified algorithm that can quickly diagnose the three-dimensional (3D) topological insulators with inversion symmetry 75 . At each time-reversal invariant momentum (TRIM), Bloch wave functions carry definite eigenvalues of the inversion operator. The topological invariant of 3D TIs can be calculated as the product of parity eigenvalues of all occupied bands at all eight TRIM. The connection of topological invariant to the discrete crystalline symmetry eigenvalues at a finite number of TRIM greatly simplifies the classification of topological materials, as tedious integration of wave functions in the Brillouin zone (BZ) is no more required.
The Fu-Kane scheme for topological invariant for inversion invariant systems greatly inspires the search for quick diagnosis schemes for systems with other symmetries. There are two closely related schemes proposed simultaneously in 2017, i.e., topological quantum chemistry (TQC) 76 and symmetry indicator (SI) theory 77,78 , which are highly efficient in diagnosing nearly all known topological phases. Consequently, high-throughput characterization of inorganic materials was quickly carried out with astonishing accuracy and efficiency 79 . Most nonmagnetic stable materials tabulated in ICSD (Inorganic Crystal Structure Database) and Materials Project have been examined and topologically classified, representing a great success of the TQC and SI theories [76][77][78] .
The same idea has also been extended to the magnetic space group (MSG) 80,81 . With the carefully designed computer tools, 549 magnetic materials from MAGNDATA database were first topologically classified by Benevig's group 82 , and several attractive magnetic topological material candidates were discovered. However, unlike the classification of nonmagnetic materials, large-scale diagnosis of magnetic materials is further subject to the uncertainty of the magnetic configurations. Experimental techniques for magnetic structures are available only for high-quality single crystals, and the determination of a magnetic structure usually requires much more effort than just characterizing a crystal structure. For this reason, MAGNDATA only contains the magnetic configurations of a tiny portion of magnetic materials. Magnetic structures of many more materials are not yet available experimentally. Thus, although the theoretical foundation and the right computer tools have been developed, massive magnetic-material diagnosis is apparently hindered by the lack of knowledge of magnetic structures.
Without experimental insights, to make progress on the highthroughput topological diagnosis of magnetic materials, theoretically, we have to first predict a candidate magnetic structure before topological characterization. Determining magnetic configuration is a highly nontrivial task, and one can only locate the magnetic ground state within the space of a limited number of initial configurations. The reliability of such theoretical prediction highly depends on the available initial configurations. In this work, we restrict ourselves to the collinear magnetic structures to feature a high-throughput calculation. We design a workflow to examine the possible collinear magnetic structures theoretically and, subsequently, determine their topological natures. Because of the theoretical nature of the proposed magnetic structures, our diagnosis no longer provides a definite answer to the topological nature of the examined magnetic materials. It only offers a possible candidate topological characterization, which awaits further experimental verification. Nevertheless, our work screens the search area and provides valuable insights into future experimental exploration. Several interesting magnetic topological materials have been proposed, making a first step toward the complete diagnosis of magnetic topological materials.

Design of workflow
To facilitate a high-throughput prediction of magnetic topological materials, we design a workflow by combining two efficient highthroughput algorithms to examine possible collinear magnetic structures of transition-metal (TM) compounds and characterize their topological nature. The workflow is designed for FireWork by utilizing atomate 83 . In this paper, we concentrate on the TM binary and ternary compounds with strong spin-orbit coupling (SOC). Our choice of TM binary compounds consists of magnetic elements including V, Cr, Mn, Co, Fe, Ni, Ru, Mo, and Cu, and elements with strong SOC including Tl, Pb, Bi, In, Sn, Sb, Te, Ga, Ge, As, Se, Al, Si, P, and S. Alkali metals including Li, Be, Na, Mg, K, Ca, Sr, Cs, and Ba are further included in ternary compounds. We restrict the total atom numbers in a primitive cell to be less than 20 to facilitate a fast diagnosis. In total, we filtered 1049 TM compounds (TMCs) from Materials Project as our initial database to demonstrate the power of our scheme. The above choice of magnetic materials is only for demonstration purposes. Our workflow is ready to apply to any magnetic material and predict its collinear magnetic state with topological classification.
We show the workflow in Fig. 1, which consists of three major steps. In the first and the third steps, we adopt open-source implementations 82,84 , which we will also briefly discuss in the Methods section. We design the second step to synergy the two established workflows simultaneously.
In step 1, we use the high-throughput framework 84 developed by M. K. Horton et al. in 2019 to search the candidate collinear magnetic ground state of a given crystal. The proposed searching algorithm has been tested over 64 collinear antiferromagnetic materials in MAGNDATA and obtained 60% correct predictions as compared to the experimental magnetic orderings 84 . Based on this workflow, collinear magnetic orderings of more than 3000 transition-metal oxides were examined 85 , displaying a high efficiency in the theoretical screening of magnetic topological materials. We note that our workflow distinguishes from the work of Frey et al. 85 in two aspects. First, the magnetic materials we examined do not contain oxygen but elements with strong SOC to favor nontrivial topology and, at the same time, to favor possibly an easy experimental preparation. Second, we did not restrict the materials to be inversion symmetric as we further automate the topological classification with TQC algorithm 76,82 , which is a universal algorithm applicable to materials both with and without inversion symmetry. In step 2, we identify the MSG of the obtained magnetic structures in step 1. The knowledge of the magnetic space group is crucial to the third step, which works with the maximal K-points and the character table of a given MSG. In the last step, we employed the magnetic TQC routine MagVasp2trace and Check Topological Magnetic Mat. module from Bilbao Crystalline Server (BCS) to examine the irreducible representations (irreps) and, finally, the topological classification. We refer to the details of the collinear magnetic structure search and the magnetic TQC to their original works 76,82,84 . We also provide a short explanation in Supplementary Methods.
The essential step to connecting the two algorithms in steps 1 and 3 is to efficiently identify the magnetic space group of any given magnetic material. To this end, we design a function based on spglib 86 which quickly finds all symmetry operations based on the paramagnetic group. We define the paramagnetic group as the symmetry group for the same magnetic unit cell but without local moments. After removing all magnetic moments, the unit cell can be a conventional one. We keep the unit cell size and all atomic positions unchanged and determine the space group of it, which we denote as G. If the paramagnetic unit cell is a conventional one, G contains fractional lattice translations, with the translational vector being the unit vectors of its primitive cell. We further construct of a gray group of G denoted as G g ¼ G þ T G. G g contains both unitary and anti-unitary operations, and it includes all the operations of the desired magnetic space group. However, G g also contains many redundant operations with respect to the magnetic unit cell, which we will remove by applying all operations in this grep group to the magnetic cell and testing their applicability. All operations passing the symmetry check for the magnetic unit cell constitute the MSG. Comparing these symmetries with those listed in MGENPOS of BCS, we obtain the correct MSG in BNS notation 87 . At the same time, we transform the magnetic unit cell to the same convention (standard magnetic primitive) of BCS as required by the topological characterization module.  Diagnosis workflow for collinear magnetic topological insulators  Table 1. Promisingly, a large portion of the magnetic materials predicted in the first step, i.e., magnetic structure search, is found to carry certain nontrivial topology. In calculations with U correction, the portion of nontrivial materials becomes less. For the various magnetic topological phases with nonzero invariants, we categorize them mainly into two classes. One contains systems with stable topological indices. According to refs. 76,82 , it includes stable magnetic topological insulators and Smith-index semimetal (SISM). The latter are semimetals with unavoidable level crossings at generic points. The second class contains semimetals with band crossings at high-symmetry points, lines, or planes in the BZ. The band crossings of these semimetals can be either symmetry enforced, which we further term as symmetry-enforced semimetal (ESM), or the semimetals with accidental band crossings (ASM). The relative ratio of two topological class materials, i.e., TI/SISM and semimetals (ESM, ASM), is shown in (b, d) of Fig. 2. There are also materials with a magnetic ground state that cannot be topologically classified. This is due to either the poor convergence of self-consistent DFT calculation with SOC or the low symmetry of the wave functions. The irreps decomposition cannot be accomplished for the given magnetic space group.
A rough comparison of the two runs indicates that the inclusion of U corrections tends to stabilize magnetic ground states but breaks the nontrivial topology. The number of materials with certain magnetic ordering increases, but the number of TI(SISM) and ESM/ASM becomes less. Without U correction, we discovered 68 magnetic TI/SISM and 53 SMs. In the second run with U correction, these numbers become 73 and 26, respectively. We note that the calculations depend on the choices of U parameters and the trend may not lead to a general conclusion.
In the following, we take a closer look at three representative materials, i.e., SrMnSn, Rb(CrS 2 ) 2 , and Cs(MnP) 2 . Single crystals of SrMnSn and Cs(MnP) 2 are available in experiments 89,90 . SrMnSn is a G-type AFM with the magnetic unit cell doubled along the c-axis forming a tetragonal lattice. Cs(MnP) 2 is isotopic with BaZn 2 P 2 . It is paramagnetic following the Curie-Weiss law above 110 K, below which it becomes antiferromagnetic 90 . According to Table SI-2 to  Table SI antiferromagnetic with either G-type or A-type depending on the values of U.
Representative ESM material: SrMnSn. In our diagnosis, SrMnSn is identified as a trivial semiconductor in DFT calculations but a symmetry-enforced semimetal in DFT + U calculations. In the following, we focus on the DFT + U results and examine the symmetry protection in more detail. SrMnSn belongs to MSG 138.528 (P c 4 2 /ncm), which contains both inversion symmetry and S ¼ T τ z=2 . Similar to MnBi 2 Te 4 , S 2 ¼ Àe Àikz takes two values ± 1 depending on k z = 0/π. Thus, S behaves as an effective timereversal symmetry at k z = 0 plane, where a Z 2 topological invariant can be defined. Moreover, in SrMnSn, due to the simultaneous presence of inversion symmetry P, one can further define a joint operation S 0 ¼ PT τ z=2 , which is anti-unitary and squares to −1. Under the protection of S 0 , despite the loss of time-reversal symmetry, SrMnSn still has Kramer's degeneracy as in a nonmagnetic system. Every Bloch band in Fig. 3b is at least twofold degenerate.
In SrMnSn, the highest valence and the lowest conduction band cross at the Fermi level, forming a magnetic Dirac point. The antiferromagnetic Dirac semimetal has been reported in orthorhombic CuMnAs 91 , in which inversion symmetry is broken but the joint operation of PT is preserved, protecting its Kramer's degeneracy. In contrast, SrMnSn is inversion symmetric but PT is broken. The Kramer's degeneracy in SrMnSn is protected by the product of PT and the half-lattice translation τ z/2 , i.e., S 0 defined above. In SrMnSn, the magnetic Dirac point along Γ − Z is protected by S 0 and a four-fold skew rotation. Such protection mechanism is a magnetic analog of the rotational-symmetryprotected Dirac point in nonmagnetic systems 92 .
Furthermore, in SrMnSn, there is also symmetry-enforced 4-fold band degeneracies at X, M, Z. More interestingly, bands along high-symmetry lines R-A-Z form magnetic nodal loops, common in nonmagnetic systems but rarely reported in magnetic materials. The symmetry protections of X, M, Z, and R-A-Z are similar; here we take X as an illustration.
At X, the magnetic little co-group is mmm1 0 containing two orthogonal skew rotations. All Bloch bands can be further labeled by their eigenvalues. Let's take the unitary glide mirror reflectionĝ x : ðÀx þ 1=2; y; z þ 1=2Þ ðÀiσ x Þ as an example. The two separated operations ofĝ x act on the coordinate and spin spaces, respectively. We neglect the irrelevant operation on local moments for simplicity. There is alsoĝ y ¼ ðx; Ày þ 1=2; z þ 1=2Þ ðÀiσ y Þ, which anticommutes withĝ x , i.e., fĝ x ;ĝ y g ¼ 0. At X, Bloch state is an eigenstate ofĝ y withĝ y ψ j i ¼ ± i ψ j i. It is easy to prove that S 0 ψ j i has the sameĝ y eigenvalue as ψ j i. While,ĝ x ψ j i and S 0ĝ x ψ j i have opposite eigenvalues of ψ j i underĝ y . Thus, four states ð ψ j i; S 0 ψ j iÞ and ðĝ x ψ j i; S 0ĝ x ψ j iÞ stay degenerate at X, forming a stable four-fold band crossing point. A similar analysis shows that along R-A,ĝ y and g z :¼ ðx þ 1=2; y þ 1=2; ÀzÞ ðÀiσ z Þ anticommute and enforces a four-fold nodal line structure. Along A − Z, the symmetries protecting the four-fold magnetic nodal line areĉ 2xy : ðy þ 1=2; x þ 1=2; ÀzÞ Ài ffiffi ffi 2 p =2ðσ x þ σ y Þ andĝ xy : ðy; x; zÞ i ffiffi ffi 2 p =2ðσ x À σ y Þ. Our high-throughput calculations correctly predict SrMnSn to be magnetic ESM. Further analysis indicates the coexistence of multiple nodal structures in this system. The presence of Kramer's degeneracy, the magnetic Dirac point, and four-fold magnetic nodal lines make SrMnSn an appealing material platform for experiments to study magnetic ESM with multiple nodal structures.
MSG of Rb(CrS 2 ) 2 is 15.89, containing inversion symmetry but no pure fractional lattice translation. Although in Rb(CrS 2 ) 2 one cannot globally define an effective time-reversal symmetry as in SrMnSn, a local definition is still possible. MSG 15.89 contains two anti-unitary operations. Given the inversion symmetry, all k-points that respect these anti-unitary symmetries will have a Kramer's pair to degenerate in energy. We observe a double band degeneracy at M 1 and a two-fold nodal line along S − R 2 , as shown in Fig. 4b.
The valence and conduction bands of Rb(CrS 2 ) 2 separate along the k-path in Fig. 4b. The Fermi surface consists mainly of the Cr-d orbitals. Due to the small SOC, the energy separation between the valence and the conduction bands is small. When the SOC is turned off, Rb(CrS 2 ) 2 becomes an SM with a nodal ring formed by the crossing of the conduction and the valence bands, which are gapped into Weyl SM with discrete Weyl nodes when SOC is included. Due to the inversion symmetry, each pair of Weyl nodes appear at the same binding energy and locate precisely at the opposite momenta. All Weyl nodes around the Fermi level have a chiral charge of ±1. The helicoid surface states connecting two nodes with opposite chirality form the Fermi arc shown in Fig. 4c, d. We observed Fermi arcs states on both (001) and (100) surfaces. On the right plot in (c), we also project four pairs of Weyl nodes to the (001) surface. It is clear that these Weyl nodes are the sources and sinks of the Fermi arcs.  Tables 2 and  3 consistent with the experimental report 90 . In both calculations, it is characterized as an axion insulator. The overall antiferromagnetic nature of the ground state agrees well with the experimental observation 90,93 . We note that we obtained two different AFM ground-state configurations depending on the values of U. Similar controversy also exists in experiments concerning the magnetic transition temperature and the ground-state order. Here we focus on the A-type AFM structure proposed in spin-polarized DFT + U calculations, in which the magnetic moments point to [-110] and antiferromagnetically align [010]. Due to the simultaneous presence of time-reversal and inversion symmetry in Cs(MnP) 2 , every band is doubly degenerate as in nonmagnetic systems. The electronic structure shows a metallic behavior with both the lowest conduction band and the top valence band crossing the Fermi level. They overlap in energy but have no level crossings. Thus, a momentum-dependent chemical potential can separate them adiabatically, leading to a semiconductor.
We calculated the surface states at (001) surface, where timereversal symmetry is broken. The net surface magnetism is finite. The surface states show a typical gapped structure. While on the (010) surface, where half-lattice translation along [001] and timereversal define an effective surface time-reversal symmetry.
Consequently, the surface states at (010) are gapless, as shown in Fig. 5e. This is a typical feature of magnetic topological insulators in which topological surface states display gaped/ gapless structures at surface breaking/respecting time-reversal symmetry. The same behavior of topological surface states in MnBi 2 T 4 62,74,94 and the calculated topological invariant η 4I = 2 confirm the A-type AFM Cs(MnP) 2 as an axion insulator.

DISCUSSIONS
In this paper, we presented an automated workflow to explore magnetic topological materials from the first principle. We combined two high-throughput algorithms to theoretically search collinear magnetic ground states and examined their topological nature. The first step of the workflow is to search the collinear magnetic ground state, which is accomplished by the algorithm proposed in ref. 84 . The subsequent topological characterization is done via the magnetic TQC theory 76,82 . We synchronized the two high-throughput workflows by efficiently determining the MSG, refining the standard crystal structure, and interacting with the Bilbao crystallographic server.
To demonstrate the power and feasibility of the designed workflow, we examined 1049 binary and ternary compounds containing V, Cr, Mn, Co, Fe, Ni, Ru, Mo, and Cu. We did two runs of calculations, i.e., one with U correction and the other without. Both runs predicted a considerably large number of collinear topological material candidates. Our high-throughput calculations predict highly relevant collinear topological material candidates and provide insightful guidance to further experimental investigations.
As the workflow only explores the collinear magnetic configurations, the readers need caution to understand our diagnosis results. The prediction presented in our work only offers a possible candidate topological characterization, which needs to be verified and corrected once the experimental magnetic structure is available. Our work provides a step toward a complete theoretical characterization of magnetic topological materials. Further improvement and optimization of the workflow are possible. (1) To accelerate the theoretical search, we have restricted the magnetic ground-state exploration in the non-relativistic case. The magnetic anisotropy was not taken into account. We include the SOC effect in step (3) and fix the magnetic moment to the z-direction of the standard conventional cell. The magnetic anisotropy energy (MAE) may be large in some materials such that their magnetic moments strongly prefer a specific direction. A high-throughput algorithm for determining MAE is desirable to get the moment orientation. (2) The current workflow only generates collinear magnetic configurations. Lattices of kagome, pyrochlore, etc., frustrated structures favoring noncollinear magnetic couplings, and interesting band topology cannot be correctly described by the current workflow. There is a strong demand to develop a workflow capable of efficiently exploring noncollinear magnetic structures. (3) We chose the initial ternary materials as the ones with TM elements, alkali atoms, and semiconductor elements of large SOC. However, most magnetic materials have relatively weak d − p hybridization, such that the SOC gap is too small to facilitate an excellent insulating behavior. Diagnosing materials with strong SOC and strong d − p hybridizations may lead to the discovery of excellent magnetic topological materials with both spin polarization and a large topological gap at the Fermi level. (4) In our DFT + U calculations, we chose the U parameter from first-principle estimation 88 . The resulting magnetic ordering and topological characterization, thus, depend on the U values. An unbias first-principle calculation without adjusting parameters, such as SCAN functional 95 , can eliminate the uncertainty induced by the DFT + U functional.

METHOD
Our first-principle calculations are based on the generalized-gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) form [96][97][98][99][100] within the framework of the density-functional theory (DFT) using projector-augmented-wave (PAW) 101 wave functions as implemented in the Vienna Ab-Initio Simulation Package (VASP) 102,103 . The cut-off energy for the wave functions was set at 500 eV. We performed two sets of calculations for comparison, i.e., one with spin-polarized DFT and the other one with spin-polarized DFT + U. In the latter, we took the ab-initio Hubbard correction parameters for V (3.909 eV), Cr (2.982 eV), Mn (4.710 eV), Co (5.237 eV), Fe (4.545 eV), Ni (5.847 eV), Cu (7.59 eV), Ru (2.972 eV), and Mo (2.431 eV) as suggested in ref. 88 . Fixed reciprocal space density was used for the Brillouin zone (BZ) sampling with (k 1 , k 2 , k 3 ) = mul × (1/a, 1/b, 1/c), where mul = (NVabc) 1/3 . V is the volume of the reciprocal cell. N is the grid density per reciprocal cell (in the unit of Å −3 ). N takes two different values, i.e., 64 in structure relaxation calculations and 100 in the charge self-consistent calculations. The self-consistent convergence threshold for total energy was set at 10 −5 eV. During the magnetic structure search, we turned off the SOC to accelerate the calculations, which was included in the topological analysis. With SOC, the local frame of the magnetic moment becomes essential for systems with strong magnetic anisotropy. In our high-throughput calculations, the magnetic moments are aligned to the z-direction of the magnetic cell by default (see the Data availability section for the mcif files.).
The magnetic moments will be further rotated to a new direction after transforming the magnetic cell to the standard primitive cell required by MVasp2trace and BCS. Finally, the edge states were calculated by using our in-house code LTM (L ibrary for T opological M aterial calculations) with the iterative Green's function approach 104 based on the maximally localized Wannier functions 105 obtained through the VASP2WANNIER90 106 . The chirality of Weyl nodes in mp-1238796 Rb(CrS 2 ) 2 was evaluated by using the WannierTools pacakge 107 .

DATA AVAILABILITY
All data are available from the corresponding author upon reasonable request. The magnetic structures files and trace files generated by Mvasp2trace and band structures are available at https://github.com/GangLi_SHT/Collinear_MTQC.

CODE AVAILABILITY
The codes for steps 1 and 3 of the workflow are open-source and free of charge except for the VASP. The codes for automatic analysis of the magnetic space group are available at https://github.com/GangLi_SHT/Collinear_MTQC.