1. Introduction

Artificial lighting accounts for one-fifth of global electricity consumption, with a half of that amount consumed by inefficient incandescent and fluorescent emission sources^{[1, 2]}. Recently, metal halide perovskites have gained much attention thanks to their outstanding optoelectronic characters including high photo-absorption efficiency, tunable emission across the entire visible spectrum, exceptional defect tolerance and low-cost synthesis processing^[3-5]. Nevertheless, so far, the best-performing halides contain hazardous and bioaccumulative lead and have unsatisfactory stability against moisture and temperature, particularly for organic–inorganic hybrid lead halide perovskites^[6-9]. Nowadays replacing the toxic Pb in the perovskite structures with alternative non-toxic, environment-friendly, earth-abundant and cost-effective elements such as transition metals is of critical importance for improving light-emitting capability^[10-14].

Recently, Cu-based halide perovskites and their derivatives (CHPs) have been emerging and attracted increasing attention because of the replacement of Pb with abundant, economic and environment-friendly Cu element^[10-14]. The rich chemistry of Cu with multiple valence states and low coordination number lead to their complex crystal structures, including 2D layer, 1D chain, 0D isolated units^[15]. These low-dimensional crystal structures, resulting in low electronic dimensionality, may give rise to large exciton binding energy (E_b) and high photoluminescent quantum yield (PLQY) due to the quantum confinement effect (QCE)^{[14, 16]}. For example, E_b of 0D Cs₄PbBr₆ with ~240 meV is much larger than 3D CsPbBr₃ with ~18 meV^{[8, 9]}. Recently, Hosono et al. reported in experiment that Cs₃Cu₂I₅ (named as 325-type) films and single crystals owned high PLQY with ~60% and ~90%, respectively^[14]. In crystal structure of Cs₃Cu₂I₅ compounds, spatially isolated [Cu₂I₅]^3– anion is surrounded by Cs⁺ cations and two Cu⁺ ions possess lower coordination with a tetrahedral and a trigonal types, respectively. Spatially isolated [Cu–I] polyhedron induces enormous E_b with ~490 meV, much higher than Pb-based halide perovskites^[14]. Zhao et al. also discovered that another type of CHP Cs₂CuX₄ (X = Cl, Br, and Br/I, named as 214-type) quantum dots possessed high PLQY of ~50% for blue-green light emitting and excellent air and photo-stability, where spatially isolated [Cu–X] tetrahedron is also surrounded by Cs⁺ ions and Cu²⁺ ion possessing tetrahedral coordinations^[17]. In the Cs₂CuBr₄ synthesis processing, the phenomenon about Cu²⁺ ion partially reducing to Cu⁺ was observed through X-ray diffraction and the X-ray photoelectron spectroscopy characteristics^[17].

Apart from 325- and 214-type CHPs which have been reported in lighting application, other phases such as 123-, 113-, 213-, 327-, and 459-type, have also been reported in experiments^{[15, 18-21]}. Diverse phases of CHPs with isolated [Cu–X] building blocks provide a treasure trove for potential applications in lighting. Nevertheless, the phase stability of A_lCu_mX_n (A = Rb and Cs; X = Cl, Br, and I) have not been clearly investigated as shown in Table 1. This has hampered further development of those materials for practical applications. Meanwhile, the ample crystal structures together with variable valence states of Cu, make it difficult to control the synthesis of CHPs. Therefore, it is necessary to provide a landscape of phase stability for CHPs and suggestions of chemical environments to synthesize particular compounds per request.

In this paper, first-principles calculations are performed to comprehensively investigate the thermodynamic stability of CHPs (A_lCu_mX_n, A = Rb and Cs; X = Cl, Br, and I) based on thermodynamic equilibrium growth condition. In total, there are 42 CHPs compounds with seven kinds of stoichiometric ratio including A₃Cu₂X₅, ACu₂X₃, A₂CuX₃, A₄Cu₅X₉, A₂CuX₄, ACuX₃, and A₃Cu₂X₇, as shown in Table 1. For chemical environments and stability calculations, all possible competing phases including compositional elements Rb, Cs, Cu, Cl, Br, I, binary compounds and ternary compounds from inorganic crystal structure database (ICSD) have been considered. We have found that most of existing A_lCu_mX_n phases in experiments have a stable growth region and positive decomposition energy in our calculations, which confirms the reliability of our calculations. Moreover, we also acquire ten new and stable phases encompassing Rb₃Cu₂I₅, Rb₃Cu₂Br₅, RbCu₂Cl₃, Rb₂CuI₃, Rb₂CuBr₄, RbCuBr₃, Rb₃Cu₂Br₇, Cs₃Cu₂Br₇, Cs₃Cu₂Cl₇, and Cs₄Cu₅Cl₉, which are not yet reported experimentally. Therefore, our work provides prospective guidance for experiment to synthesize the above phases, amplifying the scope of phases of CHPs.

2. Computational details

All of the calculations were executed with spin-unrestricted density functional theory method, as implemented in the Vienna ab-initio simulation package (VASP)^[22] by using the projector augmented wave (PAW) pseudopotential^[23]. We employ the generalized gradient approximation (GGA) parametrized by Perdew, Burker, and Ernzerhof (PBE)^[24] as electronic exchange-correlation functional. The kinetic energy cutoff with plane wave basis set is 400 eV and the k-point meshes with grid spacing of 2π × 0.025 Å^–1. All of the structures were fully relaxed until the total energy and force per atom were less than 10^–4 eV and –0.01 eV/Å, respectively.

To evaluate the thermodynamic stability of CHPs A_lCu_mX_n (A = Rb and Cs; X = Cl, Br, and I) with different types, we first calculate the chemical potential range for equilibrium growth of compound to identify the proper chemical potentials for synthesizing particular compounds in experiments^{[25, 26]}. Second, we quantitatively calculated the thermodynamic stability via considering the optimal decomposition pathways to their competing phases through linear programming. For calculation of chemical potential range, the chemical potential μ_α (α is the element that constituted to the CHPs) is constrained by the values that keep a stable host compound, and avoid the formation of other competing phases, including elemental solids. The thermodynamic equilibrium growth conditions need to satisfy the following three relations.

$$l{\mu _A} + m{\mu _{{\rm{Cu}}}} + n{\mu _{\rm{X}}} = \Delta H\left( {{{\rm{A}}_l}{\rm{C}}{{\rm{u}}_m}{{\rm{X}}_n}} \right),\;\;l,m,\;n = 1,2, \ldots ,N,$$

(1)

$${\mu _\alpha } \leqslant 0,\quad {\alpha = {\rm{A}},{\rm{Cu}},{\rm{X}}} ,$$

(2)

$${h_i}{\mu _{\rm{A}}} + {k_i}{\mu _{{\rm{Cu}}}} + {l_i}{\mu _{\rm{X}}} \leqslant \Delta H\left( {{{\rm{A}}_{{h_i}}}{\rm{C}}{{\rm{u}}_{{k_i}}}{{\rm{X}}_{{l_i}}}} \right),\;\;i = 1,2, \ldots ,N,$$

(3)

where μ_α is the chemical potential of constituent element α referring to the stable solid/gas in the growth conditions. ΔH is the formation enthalpy, A_lCu_mX_n and ${ {{\rm{A}}_{{h_i}}}{\rm{C}}{{\rm{u}}_{{k_i}}}{{\rm{X}}_{{l_i}}}}$ represent the thermodynamic equilibrium phase and all the existing competing phases, respectively. Eq. (1) is for the thermodynamic equilibrium growth condition, Eq. (2) is to avoid the atomic species that depositing to elemental phases, Eq. (3) is to hamper all the existing competing phase.

Then, the thermodynamic stability of CHPs A_lCu_mX_n was furtherly confirmed through decomposition energy calculation based on optimal decomposition pathway (ODP) using linear programming method. Specific details are as follows:

$${{\rm{A}}_l}{\rm{C}}{{\rm{u}}_m}{{\rm{X}}_n} \to \mathop \sum \limits_{i = 1}^i {x_i}{}{{\rm{A}}_{{h_i}}}{\rm{C}}{{\rm{u}}_{{k_i}}}{{\rm{X}}_{{l_i}}},$$

(4)

$$\Delta {H_d} = \mathop \sum \limits_{i = 1}^i {x_i}E\left( {{{\rm{A}}_{{h_i}}}{\rm{C}}{{\rm{u}}_{{k_i}}}{{\rm{X}}_{{l_i}}}} \right) - E\left( {{{\rm{A}}_l}{\rm{C}}{{\rm{u}}_m}{{\rm{X}}_n}} \right),$$

(5)

$$\mathop \sum \limits_{i = 1}^i {x_i}{h_i} = l,\quad\mathop \sum \limits_{i = 1}^i {x_i}{k_i} = m,\quad\mathop \sum \limits_{i = 1}^i {x_i}{l_i} = n,$$

(6)

$$0 \leqslant {x_i} \leqslant {\rm{min}}\left( {\frac{l}{{{h_i}}},\frac{m}{{{k_i}}},\frac{n}{{{l_i}}}} \right),$$

(7)

where the x_i is the molar fraction of possible competing phases, unknown variables. Eq. (4) is the decomposition pathway of host compounds, Eq. (5) is the decomposition energy between existing competing phases and host phase, Eq. (6) is the mass-conservation constraints for all atomic species, Eq. (7) is the minimum value of x_i. The linear programming approach ensures that the calculated decomposition energy is based on the optimal decomposition pathway. If the value of ΔH_d is a positive number, then this decomposition reaction is an endothermic reaction, indicating that this compound is thermodynamically stable.

3. Results and discussions

3.1 A₃Cu₂X₅ phase

A₃Cu₂X₅ phase exhibits orthorhombic crystal structure with space group of Pnma (No. 62) as shown in Fig. 1(a). It has two types of Cu⁺ ion sites, a tetrahedral site and a trigonal site, each of which constitutes spatially isolated [Cu₂X₅]^3– (X = Cl, Br, I), encircled by A⁺ (A = Rb, Cs) ions. Thermodynamic stability of six kinds of 325-type CHPs are evaluated based on thermodynamic equilibrium growth conditions and decomposition energy with ODP as shown in Fig. 1(b) and Table 2, respectively. Except for Rb₃Cu₂Cl₅, all five compounds possess stability region in cyan polygon. In consistent with experiments, Cs₃Cu₂X₅ (X = Cl, Br, I) have been synthesized in experiments and Cs₃Cu₂I₅ has been applied to the luminescent equipment with high PLQY^[14]. Our results show that Rb₃Cu₂I₅ and Rb₃Cu₂Br₅ can also be stable but Rb₃Cu₂Cl₅ is unstable. The minimum decomposition energies together with corresponding decomposition pathways calculated by linear programing method have been shown in Table 2. Interestingly, the ODP of Cs₃Cu₂Cl₅ perovskite are different for Cs₃Cu₂Br₅ and Cs₃Cu₂I₅, i.e. Cs₃Cu₂X₅ → 2CsX + CsCu₂X₃, X = Br and I, their ΔH_d are 25 and 29 meV/atom for the X site of Br and I, respectively. These decomposition reactions are non-redox processes. In contrast to these two compounds, the ODP of Cs₃Cu₂Cl₅ perovskite is a disproportionation reaction with ΔH_d of 33 meV/atom, namely Cs₃Cu₂Cl₅ → Cu + CsCl + Cs₂CuCl₄, the competing phase has CsCl secondary phase, it is also consistent with experiment that CsCl additional phase observed as minor impurity^[15]. So it should be carefully controlled in the synthesis process. If the decomposition pathway of Cs₃Cu₂Cl₅ perovskite is same as the ODP of Cs₃Cu₂X₅ (X = Br and I), then its decomposition energy is 39 meV/atom, which is a little larger than the ODP by 6 meV/atom. For Rb₃Cu₂X₅ (X = Cl, Br, and I), their ODPs are all non-redox reactions, Rb₃Cu₂I₅ phase owns the similar ODP with the Cs₃Cu₂I₅ phase by ΔH_d of 12 meV/atom, but different from the ODP of Rb₃Cu₂Br₅ phase with ΔH_d of 5 meV/atom, i.e., Rb₃Cu₂Br₅ → 4/3Rb₂CuBr₃ + 1/3RbCu₂Br₃. For Rb₃Cu₂Cl₅, because of more complex competing phases, the ODP with ΔH_d of –2 meV/atom is Rb₃Cu₂Cl₅ → 7/5RbCl + 2/5Rb₄Cu₅Cl₉, suggesting its poorer stability than its bromide and iodide counterparts.

3.2 ACu₂X₃ phase

ACu₂X₃ phase exhibits orthorhombic space group of Cmcm (No. 63), the 1D chain with two [Cu–X] tetrahedra through edge-sharing in a row, isolated by A⁺ ion, as shown in Fig. 2(a). Fig. 2(b) shows the stability region of ACu₂X₃ phase in chemical potential map. The results show that all the ACu₂X₃ (A = Rb and Cs; X = Cl, Br, and I) compounds have thermodynamic stability region with cyan polyhedron, surrounded by the boundary composed of competing phases, which is in agreement with existing experiments that five of them have been synthesized successfully in experiments^[15], except for RbCu₂Cl₃. The minimum decomposition energies together with corresponding decomposition pathways have been shown in Table 2, well matching the thermodynamic stability region calculations. CsCu₂X₃ possess the same ODP, i.e., CsCu₂X₃ → 4/3CuX + 1/3Cs₃Cu₂X₅ (X = Cl, Br, and I). Their corresponding ΔH_d are 0, 22, and 24 meV/atom for CsCu₂I₃, CsCu₂Br₃, CsCu₂Cl₃, respectively. The ODPs of RbCu₂I₃, RbCu₂Br₃, and RbCu₂Cl₃ are all different, i.e., RbCu₂I₃ → RbI + 2CuI, RbCu₂Br₃ → 3/2CuBr + 1/2Rb₂CuBr₃, RbCu₂Cl₃ → 3/4CuCl + 1/4Rb₄Cu₅Cl₉. Their decomposition energies are 23, 27, and 6 meV/atom for RbCu₂I₃, RbCu₂Br₃, and RbCu₂Cl₃, respectively. Considering that CsCu₂I₃ has a small stability region and tiny ΔH_d, experimental synthesis of its Rb counterpart may be easier since RbCu₂Cl₃, has much larger stable region and corresponding larger decomposition energies, which may be more useful in photoluminescence fields than CsCu₂X₃.

3.3 A₂Cu₁X₃ phase

In contrast from the ACu₂X₃ phase with wider 1D chain by two Cu⁺ ions in one row, the A₂CuX₃ phase with the symmetry of Pnma (No. 62) also owns a 1D chain with only one Cu⁺ ion via vertex-sharing one line in Fig. 3(a). Their phase stability regions versus the chemical potential of Cu and X element are shown in Fig. 3(b). Only Rb₂CuI₃ and Rb₂CuBr₃ possess slim stability region with cyan polyhedron. Meanwhile, Rb₂CuBr₃ phase has been synthesized successfully in experiments^[15]. Even though Rb₂CuCl₃ perovskite has no stability region from simulation, the decomposition energy calculations confirmed that its ΔH_d with ODP is a positive number with 8 meV/atom, namely Rb₂CuCl₃ → 1/2Cu + RbCl + 1/2Rb₂CuCl₄, as shown in Table 2. Indeed, Rb₂CuCl₃ exists in experiments and the additional competing phase RbCl was also observed, consistent with our predictions in the ODP^[15]. In addition, Rb₂CuI(Br)₃ phases own the same ODP with positive ΔH_d, i.e. Rb₂CuX₃ → 3/2RbX + 1/2RbCu₂X₃, X = Br, I. While the 213-type CHP with the A site of Cs element also possess the same ODP, namely Cs₂CuX₃ → 1/2CsX + 1/2Cs₃Cu₂X₅, X = Cl, Br, and I, their ΔH_d are all negative indicating their instability.

3.4 A₄Cu₅X₉ phase

A₄Cu₅X₉ compound is also a class of complex Cu-based compounds in Fig. 4(a). The crystal structure with the space group of Pc (No. 7) owns three types of Cu⁺ site, a tetragonal site, a trigonal site, and a 2-fold coordination site, which forms 0D isolated [Cu₅X₉]^4– anion, isolated by A⁺ ions. Their phase stability regions versus the chemical potential of Cu and X element are shown in Fig. 4(b). Only two compounds, Rb₄Cu₅Cl₉ and Cs₄Cu₅Cl₉ own slim stability region with cyan polygon, surrounded by competing phases. The decomposition energies are 21 and 8 meV/atom for Rb₄Cu₅Cl₉ and Cs₄Cu₅Cl₉, respectively, in Table 2. Other 459-type CHP have negative ΔH_d with the ODP. It is observed that the ODP for Rb₄Cu₅Cl₉ phase is a disproportionation reaction, i.e., Rb₄Cu₅Cl₉ → 12/5Cu + 1/5Rb₂CuCl₃ + 6/5Rb₃Cu₂Cl₇, while for the Cs₄Cu₅Cl₉, it is a non-redox reaction, i.e., Cs₄Cu₅Cl₉ → 3/4Cs₃Cu₂Cl₅ + 7/4CsCu₂Cl₃. Experimentally, Rb₄Cu₅Cl₉ perovskite had been able to synthesize successfully, matches well with our predictions^[15].

3.5 A₂CuX₄ compounds

Since all the above CHPs own monovalent Cu⁺ ion on B site, their maximum coordination number is four, induced by higher 3d¹⁰ energy level and smaller ion radius of Cu, in consistent with the report by Xiao et al.^[16]. The CHPs with divalent Cu²⁺ ion on B site also possess abundant phases, such as 214-, 113-, and 327-type, most of which are constituted with elongated [Cu–X] octahedra, induced by Jahn-Teller distortion, except for 214-type with space group of Pnma (No. 62), which is composed of 0D [Cu–X] tetrahedra, isolated by A⁺ ions, as shown in Fig. 5(a). In addition to Pnma crystal structure, Rb₂CuCl₄ compounds can exist in Cmca (No. 64) structure with 2D layers of corner-sharing [CuCl₄]^2– octahedra^[17-19]. The comparisons of total energy for non-existing compounds with two crystal structures are calculated in Table 3. The equilibrium growth region of A₂CuX₄ compounds with stable crystal phases are assessed in Fig. 5(b). Cs₂CuCl₄ undergoes non-redox ODP with ΔH_d of 34 meV/atom: Cs₂CuCl₄ → CsCl + CsCuCl₃, while Cs₂CuBr₄ goes through a redox ODP with ΔH_d of 14 meV/atom, i.e., Cs₂CuBr₄ → 1/2CsBr + 1/4CsBr₃ + 1/2CsCuBr₃ + 1/4Cs₃Cu₂Br₅, due to more competing phases. Experimentally, the Cu²⁺ ions of Cs₂CuBr₄ perovskite can be partially reduced to Cu⁺ ions. Our calculations reveal that the Cu⁺ may exist in the occurrence of Cs₃Cu₂Br₅ compound^[17]. Although Rb₂CuCl₄ has been reported in experiments, we haven’t found its stability region and its ΔH_d is also negative (–9 meV/atom) with ODP, Rb₂CuCl₄ → 1/2RbCl + 1/2Rb₃Cu₂Cl₇. The discrepancies may be ascribed to the computational errors. Meanwhile, we suggested to double check the experimental results for Rb₂CuCl₄, in particular to possible existence of secondary phases including RbCl and Rb₃Cu₂Cl₇. For A₂CuX₄ compounds, we discover a novel and stable 214-type CHPs Rb₂CuBr₄ via phase stability region in cyan polygon and decomposition energy with ODP by ΔH_d of 10 meV/atom, i.e. Rb₂CuBr₄ → 2RbBr + CuBr₂.

3.6 ACuX₃ compounds

Due to the divalent Cu²⁺ ions, the crystal structure of ACuX₃ compounds also own elongated [Cu–X] octahedra via corner- and face-sharing connection, forming spatially isolated 1D coplanar chain and 3D network through vertex-sharing, surrounded by A⁺ ion in Fig. 6(a). Experimentally, CsCuBr₃, CsCuCl₃, and RbCuCl₃ compounds have been synthesized successfully, where CsCuBr₃ phase with space group of C2221 symmetry (No. 20) owns 3D network via vertex-sharing connection by coplanar double [CuX₆] octahedron unit, while CsCuCl₃ and RbCuCl₃ phases possess 1D chain through face-sharing [CuX₆] octahedron, they have different space group P6₁22 (No. 178) and Pbcn (No. 60), respectively^{[20, 21]}. The stable structures of other non-existing 113-type compounds in experiment are calculated according to the above existing crystal structures, as summed in Table 4. These stable structures with the lowest energy all own C2221 symmetry, same as CsCuBr₃ phase. Then the phase stability regions of the above six 113-type CHP are evaluated via thermodynamic equilibrium growth conditions, as shown in Fig. 6(b). In consistent with the experiment, CsCuBr(Cl)₃ compounds all possess slim stability region. Theoretical calculations did not find stability region for experiment existing RbCuCl₃, which may due to the computational errors or experimental ignorance of secondary phases, similar to the case of Rb₂CuCl₄. Meanwhile, we discover a new stable RbCuBr₃, with long and narrow stability region and positive decomposition energies of 20 meV/atom, i.e., RbCuBr₃ → RbBr + CuBr₂, as shown in Table 2. This means that RbCuBr₃ compound is not prone to disintegrating their competing phases RbBr and CuBr₂. Based on the above research, we also find a new 113-type CHPs with C2221 symmetry RbCuBr₃ to be stable.

3.7 A₃Cu₂X₇ phase

Last but not least, A₃Cu₂X₇ phase with Cu²⁺ ion on B site owns spatially isolated 2D [Cu₂Cl₇]^3– anion layer with double [Cu–X] octahedron layer by A⁺ ions, disparity with Rb₂CuCl₄ by single [Cu–Cl] octahedron layer, as illustrated in Fig. 7(a). Their stability growth regions in cyan polygon are also assessed against μ_Cu and μ_X in Fig. 7(b). Except for Cs(Rb)₃Cu₂I₇ compounds, others possess stability region with cyan polygon, consistent with the predictions of decomposition energy in Table 2. Experimentally, Rb₃Cu₂Cl₇ compound has been synthesized successfully, our prediction signifies the experimental discovery^[15]. In the 327-type CHP, we additionally discover three stable compounds, named as Cs₃Cu₂Br(Cl)₇ and Rb₃Cu₂Br₇. Their ODPs are all non-redox reaction, where Cs₃Cu₂Br(Cl)₇ compounds own similar ODP with ΔH_d of 9 and 9 meV/Å for the X site of Br and Cl element, respectively, i.e. Cs₃Cu₂X₇ → Cs₂CuX₄ + CsCuX₃, X = Br, Cl. For Rb₃Cu₂Br₇ phase, the ODP with ΔH_d of 23 meV/Å is Rb₃Cu₂Br₇ → 3RbBr + 2CuBr ₂. Therefore, three new compounds in the 327-type CHPs, namely Cs₃Cu₂Br(Cl)₇ and Rb₃Cu₂Br₇ are predicted to be stable.

Table 5 summarizes the phase stability cases encompassing phase stability region and decomposition energy of ODP, whether or not they exist at room temperature. Most of existing phases are stable according to our predictions, except for Rb₂CuCl₄ and RbCuCl₃ perovskites because they have more binary and ternary secondary phases and computational errors. Surprisingly, we also discover 10 novel CHPs with specific stability region and positive decomposition energy with ODP (i.e., Rb₃Cu₂I(Br)₅, RbCu₂Cl₃, Rb₂CuI₃, Rb₂CuBr₄, RbCuBr₃, Rb₃Cu₂Br₇, Cs₃Cu₂Br(Cl)₇ and Cs₄Cu₅Cl₉) which have not yet been reported in experiment. Our results offer importance guidance to synthesize these phases, consequently broadening the range of existing CHPs.

4. Conclusion

In summary, we have systematically studied the stability of all ternary CHPs considering thermodynamic equilibrium growth conditions and decomposition energies. They all own lower electronic dimensionality including 2D layered, 1D chained and 0D isolated unit, surrounded by A⁺ ions. The coordination number of monovalent Cu (Cu⁺) in the CHP is less than 4, namely, 2-fold, trigonal, tetragonal site. The vast majority of CHPs with Cu²⁺ ion possess elongated octahedron induced by Jahn-Teller distortion. Most of existing CHPs are predicted to be stable, which is in consistent with the experiment. Furthermore, we discovered ten novel phases with specific stability region and positive decomposition energy with ODP via element exchange method, i.e. Rb₃Cu₂I(Br)₅, RbCu₂Cl₃, Rb₂CuI₃, Rb₂CuBr₄, RbCuBr₃, Rb₃Cu₂Br₇, Cs₃Cu₂Br(Cl)₇ and Cs₄Cu₅Cl₉, which are not yet reported in experiment. Our predictions may provide insights for experimentalists to synthesize more novel inorganic CHPs, and will therefore tremendously expand the scope of existing CHP with promising applications.

Acknowledgements

The authors acknowledge funding support from National Natural Science Foundation of China (grant No. 11674237 and 51602211); National Key Research and Development Program of China (grant No. 2016YFB0700700); Natural Science Foundation of Jiangsu Province of China (grant No. BK20160299); the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD); and China Post-doctoral Foundation (grant No. 7131705619). The theoretical work was carried out at National Supercomputer Center in Tianjin and the calculations were performed on TianHe-1(A).