Ⅰ. INTRODUCTION

Two-dimensional metal-organic frameworks (MOFs) have attracted great attention for their intriguing properties. Compared with the existed 2D organic systems [1-5], especially the carbon-based sheets, many 2D MOFs not only exhibit great capacity of conductivity and light-harvesting, but also manifest applications in other regions as catalysis, gas storage, and particular magnetism attributed to the coordinating with metal atoms [6-11]. Transition metal phthalocyanine-based 2D MOFs are one of such materials that have long been studied in order to achieve tunable electronic, optical, and magnetic properties through substitution of transition metal (TM) atoms. In the structure of TMPc, the TM atoms are located at center of the host phthalocyanine (Pc) and coordinated with four nitrogen atoms, so they are prevented from clustering when the TMPc molecules form 2D structure. For the TMPc, the magnetism is owing to the unfilled d shell of the TM atom, and the optical properties are generally determined mainly by the ligands. Successfully synthesizing 2D periodic TMPc-based monolayer in experiment [12, 13] provided a novel synthetic method for exploring poly-TMPc structures with fascinating properties, and new types of 2D MOFs based on the TMPcs and their derivatives have been synthesized subsequently [14-18]. In addition, theoretical studies predicted potential applications of these 2D TMPc-based MOFs in hydrogen storage [19], CO$_2$ capture [20], light absorption [21], magnetic storage [22], and spintronics [23, 24].

In order to expand family of 2D TMPc-based MOFs and further explore their properties and applications, here we propose an intriguing TMPc-based porous monolayer with geometric frustrated kagome lattice. The spin ordering of kagome systems are of interest and are investigated intensively due to the existence of spin frustration [25-28]. The kagome lattice is composed of interlaced equilateral triangles, so that if the spins on the vertices prefer to be antiferromagnetic (AFM) coupling, the structural arrangement of spins precludes simultaneous satisfaction of the nearest-neighbor interactions. Especially, the $S$$=$1/2 AFM kagome systems may be ideal candidates for spin-liquid ground state [29]. In our proposed 2D TMPc-based kagome framework, the TM atoms embedded in the center of host Pc and dispersed uniformly on the lattice vertices can act as spin centers in the kagome system, making it an ideal platform to investigate the magnetic coupling with tunable spins $S$ relied on the metal ions and related physical properties.

In this work, by employing first-principles calculations, we have systematically investigated the electronic structure, magnetic behaviors, and optical properties of free-standing 2D kag-TMPc frameworks. Both the 2D metal-free (MF) and TM phthalocyanine frameworks are taken into account and the kag-MFPc framework could be regarded as the host structure. For the kag-TMPc frameworks, we choose 3d TM atoms varying from Cr to Zn in the periodic table. The 2D kag-MFPc framework is a nonmagnetic semiconductor with a direct band gap of 1.25 eV. When the TM atoms (Cr, Mn, Fe, Co, Ni, Cu, Zn) are introduced, the magnetic moments per unit cell of the kag-TMPc are approximate to 4, 3, 2, 1, 0, 1, and 0 μ$_\rm{B}$ respectively. Only the 2D kag-NiPc and kag-ZnPc frameworks are nonmagnetic, and the 2D kag-MnPc framework exhibits a stable ferromagnetic (FM) ordering due to the effective magnetic exchange interactions and prefers out-of-plane easy-magnetization axis. The Monte Carlo (MC) simulations based on the Heisenberg model suggest that the Curie temperature ($T_\rm{C}$) is about 125 K, close to the situation in 2D MnPc square framework [23]. The 2D kag-CrPc framework is an $S$$=$2 kagome antiferromagnet with spin arrangement of RT3 state. The optical absorption spectra indicate that the 2D kag-TMPc frameworks have even wider light absorption ranges than those of the TMPc molecules.

Ⅱ. COMPUTATIONAL PROCEDURES

Our first-principles calculations based on spin-polarized density functional theory (DFT) are carried out in the framework of the projector augmented wave (PAW) method [30, 31] with the Perdew-Burke-Ernzerhof type generalized gradient approximation [32], implemented in Vienna Ab initio Simulation Package (VASP) [33]. All self-consistent calculations are performed with an energy cutoff of 400 eV for the plane-wave expansion of electron wavefuctions. The Brillouin zone integration is generated in a grid of 5$\times$5$\times$1 according to the Monkhorst-Pack scheme [34] and 3$\times$6$\times$1 for 2$\times$1 supercell. Periodic boundary condition is applied to simulate the 2D sheets, and the vacuum space along the $z$ direction is larger than 20 Å to avoid interactions between neighbouring slabs. In order to treat the strongly correlated d orbitals of TM atoms, GGA+$U$ scheme [35] is adopted to describe the Coulomb and exchange corrections, while the delocalized s and p electrons are treated by the standard GGA method. The correlation energy $U$ and the exchange energy $J$ of the 3d electrons are set to 4 eV and 1 eV respectively. These values had been tested and used in previous studies of similar TM-organic complex systems [23, 36, 37]. Our tests also confirm that the changes of $U$ and $J$ values will not change our conclusions qualitatively. All the structures have been fully relaxed through conjugated-gradient algorithm with criteria of 10$^{-5}$ eV and 0.01 eV/Å of total energy and Hellman-Feymann force, respectively.

Ⅲ. RESULTS AND DISCUSSION

The MF phthalocyanine is always regarded as the host material of abundant TMPc derivates. Comprehension of intrinsic properties of organic ligand parts separated from the TM atoms in this type of system is of significance and can inform us the specific characters related to H, C, and N. For an MF phthalocyanine monomer, i.e., H$_2$Pc molecule, two H atoms are bonded to opposite nitrogen atoms in the Pc ring exhibiting an aromatic, 18 $\pi$-electrons conjugation. When building the 2D framework, we remove the outmost 16 C atoms in four benzene groups of H$_2$Pc molecules and link them periodically in an arrangement of kagome lattice as shown in FIG. 1(a). The fully structural relaxation result suggests that the lowest energy structure is planar without any buckling and preserves the six-fold symmetry, displaying characteristics of a perfect kagome lattice. To check thermal stability of the 2D kag-H$_2$Pc framework, we have performed $ab$ $initio$ molecular dynamics (MD) simulation at 300 K with time step of 1 fs using the Nose heat bath scheme [38] and the phonon spectrum is also calculated to examine the dynamic stability. FIG. 1(c) gives snapshot of the geometric structure of kag-H$_2$Pc framework after 4.5 ps MD simulations at 300 K that the geometry can still be maintained. The variation of total energy with time (FIG. 1(d)) shows that the total energy remains almost invariant during the MD simulation and all the vibrational modes are real in whole Brillouin zone (FIG. 1(e)), suggesting the framework is thermally and dynamically stable. The optimized lattice constant is 19.00 Å, while the band structure shows that this extended $\pi$-conjugated system is a nonmagnetic semiconductor, and has a direct band gap of 1.25 eV which is located at the $\Gamma$ point. Furthermore, our calculations under the bixial strain show that the energy gap of 2D kag-H$_2$Pc keeps almost unchanged (with variation within 0.02 eV, not shown) as the lattice constant is gradually increased by 10%.

Next we will discuss the situations of 3d TM atoms coordinated frameworks. The TM atoms are located at the center of the Pc pores, and each TM atom donates two 4s electrons to the Pc ligands so that it is bonded with the adjacent nitrogen and remains in the +2 valence state. Structural optimization shows that the TM atoms are dispersed uniformly on the kagome vertices, and the lattice constants of the 2D kag-TMPc frameworks (TM atom varying from Cr to Zn) are shown in FIG. 6 and Table Ⅰ. The bond length between the TM atom and its adjacent N atom increases along the association direction of TMPc units ($d_\rm{TM-N}$$^\rm{a}$) and decreases in the direction perpendicular to the former ($d_\rm{TM-N}$$^\rm{b}$) compared to the TMPc molecules ($d_\rm{TM-N}$$^\rm{M}$). Moreover, the bond lengths of all of the TM--N bonds in the kag-TMPc frameworks and individual TMPc molecules decrease from TM$=$Cr to TM$=$Ni and then increase until TM$=$Zn, which well agrees with the previous study [39].

Spin polarization calculations are carried out to investigate the magnetism in the 2D kag-TMPc frameworks. The calculation results show that except for the kag-NiPc and kag-ZnPc frameworks which are nonmagnetic, all the other kag-TMPc frameworks have spin-polarized ground states. The magnetic moments per unit cell of the 2D kag-TMPc frameworks are 12, 9, 6, 3, and 3 μ$_\rm{B}$ in their FM states for TM$=$Cr, Mn, Fe, Co, and Cu, respectively. In order to understand origination of the magnetism in these frameworks, we have analyzed the spatial distribution of their spin-polarized electron densities ($\rho$$=$$\rho$$^\uparrow$$-$$\rho$$^\downarrow$), and found that the magnetic moment mainly originates from the TM atoms. Moreover, the adjacent N atoms are slightly spin-polarized, and there exists AFM coupling between the TM atoms and their adjacent N atoms, except for the kag-CuPc framework in which the FM coupling between two kinds of atoms is found in the ground state. As an example, the isosurface and 2D slice of the spin-polarized electron density of the 2D FM kag-MnPc framework are illustrated in FIG. 2 (a) and (b). Considering that there are three TM atoms in one unit cell, we can deduce that every TM atom has magnetic moment of about 4, 3, 2, 1, and 1 μ$_\rm{B}$ in the 2D kag-TMPc framework with TM$=$Cr, Mn, Fe, Co, and Cu, respectively.

Since the 4s electrons of the central TM atoms have been involved in the chemical bonding between the TM atoms and the nearest N atoms, the local magnetism in the 2D kag-TMPc framework is considered to originate from spin-splitting of the 3d electrons of the central TM atoms. Due to the kagome crystalline field and D$_\rm{4h}$ symmetry of the TMPc molecule, degeneracy of five d orbitals is removed and they are split into three groups : two approximately degenerate orbitals d$_{xz}$/d$_{yz}$, two approximately degenerate orbitals d$_{xy}$/d$_{{{x}^2-{y}^2}}$, and a d$_{{z}^2}$ orbital. This can be clearly confirmed by projected density of states (PDOS) of the five 3d orbitals in the kag-TMPc frameworks, as shown in FIG. 2(c) and FIG. 3. The d$_{{z}^2}$ orbital is extremely localized in the energy scale compared to the other d orbitals, because some orbital hybridizations between the latter d orbitals and the Pc ligands result in that the latter d orbitals exhibit wider energy broadening. After analyzing the 3d electron configurations in detail, we can see that in the 2D kag-CrPc framework, the four d electrons of Cr occupy the d$_{xz}$/d$_{yz}$, d$_{{z}^2}$ and partial d$_{xy}$/d$_{{{x}^2-{y}^2}}$ orbitals in spin up channel, giving rise to a magnetic moment of about 4 μ$_\rm{B}$. After that, as the TM atom changes from Mn to Ni, the increased d electrons are sequentially filled into the d$_{xz}$/d$_{yz}$, d$_{{z}^2}$ and partial d$_{xy}$/d$_{{{x}^2-{y}^2}}$ orbitals in the spin down channel, resulting in the magnetic moments of about 3 μ$_\rm{B}$ (Mn), 2 μ$_\rm{B}$ (Fe), 1 μ$_\rm{B}$ (Co), and 0 μ$_\rm{B}$ (Ni), respectively. In the 2D kag-CuPc framework, only the partial d$_{xy}$/d$_{{{x}^2-{y}^2}}$ orbitals in the spin down channel are yet unoccupied so that 1 μ$_\rm{B}$ for each Cu ion is observed. All the 3d orbitals of Zn are fully occupied in the 2D kag-ZnPc framework, so it is nonmagnetic.

To learn about the magnetic coupling interactions and magnetic ground state of the 2D kag-TMPc systems, we have considered and calculated their electronic states with both the FM and AFM coupling between the TM atoms. In the FM state, all the spins on the TM atoms are parallel. But in the AFM state, it is noteworthy that there possibly exists spin frustration, since each unit cell contains three TM atoms located on the equilateral triangle and the AFM coupling between arbitrarily two adjacent TM atoms can't be satisfied simultaneously. So we reverse one spin of each unit cell in the FM state and obtain a spin frustrated (SF) magnetic state which is considered to represent the AFM state. The calculated magnetic exchange energy $E_\rm{ex}$ ($E_\rm{ex}$$=$$E_\rm{SF}$$-$$E_\rm{FM}$) per unit cell of all the 2D kag-TMPc frameworks are also listed in the FIG. 6 and Table Ⅰ. We can observe that only the 2D kag-CrPc framework exhibits a spin frustrated ground state among seven 2D kag-TMPc frameworks. It can be explained as follows: as shown in PDOS of the 2D kag-CrPc framework (FIG. 3(a)), the highest occupied (d$_{xz}$/d$_{yz}$) and the lowest dominant unoccupied (d$_{{z}^2}$) d orbitals are respectively located in the spin-up and spin-down channels, thus the virtual hopping between them results in the AFM coupling between the Cr atoms. On the other hand, among all the concerned systems, the 2D kag-MnPc framework displays the most stable FM ground state with a magnetic exchange energy of 61.1 meV. Detailed analyses suggest that it is mainly attributed to the virtual hopping between the d$_{xz}$/d$_{yz}$ orbitals in the same spin down channel near the Fermi level, and this FM coupling is quite enhanced through the strong p-d hybridization as shown in FIG. 2(c). The d orbitals (mainly d$_{xz}$/d$_{yz}$) overlap with the conjugated p orbitals of the linked pyrroles in proximity of the Fermi level, which provides effective p-d exchange and enhances the exchange interactions so as to leading to stable long-range FM ordering. And the smaller band gap of 2D kag-MnPc also contributes to the reinforced virtual hopping and then the enhanced FM ordering. In contrast, the other 2D kag-TMPc frameworks possess larger band gaps so that their magnetic couplings are weaker relatively. In addition, the electronic band structure (FIG. 2(c)) indicates that the 2D kag-MnPc framework displays a characteristic that both the conduction band minimum (CBM) and valence band maximum (VBM) are dominated by spin down orbitals, so the robust half-metal nature with $-$100% spin polarization can be realized by adjusting the Fermi level via carrier doping.

Since the 2D kag-MnPc framework exhibits stable FM ordering, we have explored the magnetic behaviors under finite temperature by performing the MC simulations. Without an external field, the Hamiltonian of the Heisenberg model can be written as:

$ \begin{eqnarray} H = - \sum\limits_{\langle {i, j}\rangle } {{J_\alpha} {\textbf{S}_i}\cdot{\textbf{S}_j}} \end{eqnarray} $

(1)

where exchange integral $J_{\alpha}$ is defined in FIG. 4(a) with $\alpha$ referring to different types, ${\textbf{S}_i}$ and ${\textbf{S}_j}$ are the magnetic moments at site $i$ and $j$ respectively. If the next-nearest magnetic interactions are ignored, the nearest-neighbor exchange integral $J_1$ can be obtained through the formula:

$ \begin{eqnarray*} {J_1}={E_\rm{ex}}/(8|S|^2) \end{eqnarray*} $

Therefore, we estimate that the exchange integral $J_1$ in the 2D kag-MnPc framework is 0.85 meV. In the MC simulations, we employ a 100$\times$100 supercell and the simulation has lasted for 10$^6$ loops. In each loop, all the spins on the Mn atoms change their spatial direction under the law of Boltzmann distribution. The results of the MC simulations are shown in FIG. 4(b). In order to display the magnetic phase transition more clear, we have also calculated the heat capacity ($C_\rm{V}$) of the system using the expression:

$\begin{eqnarray*} {C_\rm{V}}={\mathop {\lim }\limits_{\Delta T \to 0} \frac{\Delta E_T}{\Delta T}} \end{eqnarray*} $

Here $\Delta$$E_T$ is the change of total energy of the system when the temperature raises from $T$ to $T$$+$$\Delta$$T$. Obviously, the transition from the FM state to paramagnetic state of the 2D kag-MnPc framework is a second order phase transition and the corresponding Curie temperature $T_\rm{C}$ is about 125 K, close to the analogous system of 2D MnPc square framework [23]. Moreover, the magnetic anisotropy of the 2D FM framework materials should be of importance for their applications. Our non-collinear calculations with the spin-orbit coupling suggested that the magnetic anisotropy energy (MAE) of the 2D kag-MnPc framework is 1.18 meV per unit cell (defined as MAE$=$$E_{\|}$$-$$E_{\bot}$), indicating the easy-magnetization axis is perpendicular to the kag-MnPc framework. Although this MAE value is not large enough, the preferred out-of-plane easy-magnetization axis is still interesting and valuable for the applications of the 2D kag-MnPc framework in the future spintronics and magnetic storage.

The spin ordering of kagome antiferromagnets is very attractive because of various unique features in their phase transitions. Messio et al. [25] provided the regular magnetic orders in frustrated kagome lattice and discussed relative stability of various possible magnetic orders, with the second- and third-neighbor exchange interactions being also taken into account. The geometries of various exchange paths in our 2D kag-CrPc framework are defined in FIG. 4(a), and we try to analyze its magnetic ordering in the ground state by adopting the $J_1$-$J_2$-$J_3$-$J_d$ model [26]. The total energy calculations have been performed for different spin states with varied spin directions of each Cr atoms on a 2$\times$1 supercell, and $J$$_{\alpha}$ ($\alpha$$=$1, 2, 3, 4, $d$) are deduced from the energy difference between them. The results indicate that $J_3$ is much smaller so as to be neglected and the exchange integrals $J_1$, $J_2$, and $J_d$ are $-$0.091, $-$0.009, and 0.003 meV, respectively. The nearest-neighbor exchange interactions are expected to dominate the magnetic ground state of the 2D kag-CrPc framework because $J_1$ is larger than the others by an order of magnitude.

So we employ the $XXZ$ Heisenberg model of the magnetic orders of kagome antiferromagnet, in which only the nearest-neighbor exchange interactions are considered [40-42]. In this model, the two possible stable magnetic orders are so-called $q$$=$0 ($Q_0$) state and $\sqrt{3}$$\times$$\sqrt{3}$ (RT3) state, corresponding to the spin states of AFM1 and AFM3 shown in FIG. 5, respectively. In the $XXZ$ Heisenberg Hamiltonian:

$ \begin{eqnarray} H = - \sum\limits_{\langle {i, j}\rangle } {J({\textbf{S}_i^x}{\textbf{S}_j^x} + {\textbf{S}_i^y}{\textbf{S}_j^y} + \Delta {\textbf{S}_i^z}{\textbf{S}_j^z})} \end{eqnarray} $

(2)

where $\Delta$ is the anisotropy of the easy-plane type, 0$\leq$$\Delta$$\leq$1. When $\Delta$$=$0, it becomes the $XY$ model [43], and while $\Delta$$=$1, it remains the same as isotropic Heisenberg model [44]. We have carried out the non-collinear calculations with spin-orbit coupling (SOC) for the 2D kag-CrPc framework in one unit cell. The calculation for the RT3 state needs a $\sqrt{3}$$\times$$\sqrt{3}$ supercell which goes beyond our computing power, so we do not consider it for the time being. Different spin states of FM, AFM1, and AFM2 (FIG. 5) with various included angles $\theta$ between the spin direction and horizon ($xy$) plane are taken into account. It is indeed that in the 2D kag-CrPc framework with one unit cell, the energy is lowered as the $\theta$ goes down and then the AFM1 state ($Q_0$) is the most stable magnetic order in which the spin directions are in the $xy$ plane of kagome vertices and the included angle between each two is 120$^\rm{o}$. However, the above calculations cannot refer to the AFM3 (RT3) state, so we adopt another method to evaluate its possibility in 2D kag-CrPc as follows. We have calculated magnetic anisotropy energy $E_\rm{MA}$ of the 2D kag-CrPc framework, which is 0.67 meV. Based on the Heisenberg Hamiltonian (Eq.(2)), the magnetic anisotropy energy $E_\rm{MA}$ of 2D kag-CrPc can be expressed as $E_\rm{MA}$$=$6$J$$|$$\Delta$$-$1$|$$|$$S$$|$$^2$, so $|$$\Delta$$-$1$|$$=$4$E_\rm{MA}$/(3$E_\rm{ex}$)=0.053. This shows that $\Delta$ is close to 1 and the 2D kag-CrPc framework can be regarded as approximately isotropic Heisenberg model. According to the ground-state phase diagram deduced in a previous theoretical study [41], the 2D kag-CrPc framework with $S$$=$2 and $\Delta$ close to 1 should exhibit a ground state of RT3 magnetic order. We also note that Müller et al. ever proposed that the kagome-lattice Heisenberg antiferromagnet with arbitrary spin $S$ favors RT3 state over $Q_0$ state for all temperatures $T$$\geq$0 K [44].

The TMPc compounds and relevant derivatives are well known for their extensive applications as optical materials. Especially, the 2D ZnPc-based metal organic framework has also been proven to have an appropriate absorber gap and interface band alignment as the donor and acceptor so that it is promising for high efficiency solar cell [21], which motivates us to study the optical properties of the 2D kag-TMPc frameworks in depth. The electronic density of states (DOS) of the H$_2$Pc/TMPc molecules which have been referred in a previous result [45] and the corresponding 2D kag-H$_2$Pc/TMPc frameworks are plotted in FIG. 6. We can see that when the H$_2$/PcTMPc molecules are connected into periodic kagome frameworks, their electronic structures are indeed restructured due to the nature of orbital interactions and the energy gaps are also reduced slightly in the 2D frameworks. For the 2D kag-TMPc framework with TM$=$Cr, Mn, Fe, and Co, the energy gap has been obviously reduced compared to that in a single molecule (also see FIG. 6 and Table Ⅰ) due to the association of TMPc. Meanwhile, in the 2D kag-NiPc, kag-CuPc, and kag-ZnPc frameworks, the energy gaps are close to those of the relevant molecules. We also find that for the 2D frameworks, there are some new states in the energy regions where the corresponding single molecules have no states, and the association induced band dispersion should be responsible for them.

To explore the specific performance in light absorption, we carry out the calculations of the absorption coefficients of the TMPc molecules and the corresponding 2D kag-TMPc frameworks. The frequency dependent dielectric function $\varepsilon$($\omega$)$=$$\varepsilon$$_1$($\omega$)$+$$i$$\varepsilon$$_2$($\omega$) includes the real and imaginary parts, in which the imaginary part $\varepsilon$$_2$($\omega$) is derived from the interband optical transition while the real part $\varepsilon$$_1$($\omega$) can be derived from $\varepsilon$$_2$($\omega$) following the Kramer-Kronig relationship [46]. Thus, the absorption coefficient $\alpha$($\omega$) can be calculated by the expression:

$ \begin{eqnarray} \alpha (\omega ) = \sqrt 2 \omega {\left[ {{\varepsilon _1} + {{\left( {{\varepsilon _1}^2(\omega ) + {\varepsilon _2}^2(\omega )} \right)}^{1/2}}} \right]^{1/2}} \end{eqnarray} $

(3)

In order to compare the intensities of optical absorptions between the TMPc molecules and TMPc based frameworks, we divide the absorption coefficient $\alpha$($\omega$) by the atoms number in one unit cell. The estimated normalized absorption densities of TMPc molecules and corresponding 2D kag-TMPc frameworks are illustrated in FIG. 7. It is noted that all the 2D kag-TMPc frameworks present better performance of light absorption in visible region (1.6--2.8 eV) than the corresponding single molecule, which may be ascribed to appearances of new states aroused by the band dispersion as shown before. On the other hand, in the 2D kag-MnPc and kag-FePc frameworks, an extra absorption peak occurs at $\sim$0.5 eV providing absorption of infrared waves, corresponding to the smallest energy gap of these two 2D frameworks (see FIG. 6 and Table Ⅰ). These results demonstrate that the 2D kag-TMPc frameworks possess the broader light absorption wavebands than the single TMPc molecule, especially those in the visible region. Considering that the 2D chemically bonded frameworks also possess higher stabilities, specific surface areas, and efficiencies than the TMPc molecule thin films, we think that the former should have the promising potentials in optical devices.

Ⅳ. SUMMARY

By employing the first-principles calculations, we have systematically studied the electronic, magnetic, and optical properties of the 2D TMPc-based kagome metal-organic frameworks. Our calculations indicate that the 2D kag-MnPc framework is an FM semiconductor with narrow band gap of 0.09 eV and the Curie transition temperature is about 125 K, and the easy-magnetization axis is perpendicular to the framework plane. The 2D kag-CrPc framework is an $S$$=$2 kagome antiferromagnet with possible magnetic order of RT3 state. The optical absorption calculations show that when the TMPc molecules form the 2D kag-TMPc frameworks, the absorption regions are broadened in the visible and even expanded to the infrared wave bands. All the conclusions above provide the evidences for bright prospect of the 2D kag-TMPc frameworks used in the future nano-electronics, spintronics, and optics.

Ⅴ. ACKNOWLEDGMENTS

This work was supported by the National Key Research & Development Program of China (No.2016YFA0200604 and No.2017YFA0204904), the National Natural Science Foundation of China (No.21473174), the Fundamental Research Funds for the Central Universities (No.WK2340000074 and No.WK2060190084). The computational resources of Super-computing Center of University of Science and Technology of China, Supercomputing Center of Chinese Academy of Sciences, Tianjing, and Shanghai Supercomputer Centers are also acknowledged.