SDSPT2s: SDSPT2 with Selection
Yibo Lei Yang Guo Bingbing Suo Wenjian Liu


§ ABSTRACT 
As an approximation to SDSCI [static-dynamic-static (SDS) configuration interaction (CI), a minimal MRCI; Theor. Chem. Acc. 133, 1481 (2014)], SDSPT2 [Mol. Phys. 115, 2696 (2017)]
is a CI-like multireference (MR) second-order perturbation theory (PT2)
that treats single and multiple roots on an equal footing. This feature permits the use of configuration selection
over a large complete active space (CAS) P to end up with a much reduced reference space P̃,
which is connected only with a portion (Q̃_1) of the full first-order interacting space Q connected to P.
The effective interacting Q̃ space can further be truncated by an integral-based cutoff threshold.
With marginal loss of accuracy,
the selection-truncation procedure, along with an efficient evaluation and storage of internal contraction coefficients,
renders SDSPT2s (SDSPT2 with selection) applicable to systems that cannot be handled by the parent CAS-based SDSPT2,
as demonstrated by several challenging showcases.


Diffusion-limited settling of highly porous particles in density-stratified fluids
    Daniel M. Harris
    September 3, 2024
==================================================================================







§ INTRODUCTION

A system of strongly correlated electrons is characterized by a large number of energetically adjacent and singly occupied frontier orbitals.
The strong static/nondynamic correlation among such orbitals renders the many-electron wave function
a heavy mixture of a huge number of Slater determinants or
equivalently configuration state functions (CSF). As such, such systems (especially of low spins) go beyond the capability of
single reference methods. Instead, the use of a multireference (MR) method is mandatory. The available MR methods can be classified into
three families<cit.>: static-then-dynamic (SD), dynamic-then-static (DS), and static-dynamic-static (SDS).
Briefly, the SD type of methods start with a diagonalization of the bare molecular Hamiltonian projected onto
a reference/active space P={Φ_R; R=[1,N_R]}, so as to
obtain a reference state Ψ^(0)_k=∑_R=1^N_RΦ_RC̅_R k^(0).
This step captures static correlation, especially when Ψ^(0)_k is only qualitatively or semi-quantitatively correct.
The remaining dynamic correction to Ψ^(0)_k can be accounted for in a number of ways, even just to first order Ψ^(1)_k.
Classic examples of this family of methods include complete active space second-order perturbation theory<cit.>,
multiconfigurtion quasi-degenerate perturbation theory<cit.>, and n-electron valence
second-order perturbation theory (NEVPT2)<cit.>.
A common feature of such methods lies in that the coefficients C̅_R k^(0) of the reference state Ψ^(0)_k are used to construct Ψ^(1)_k
for dynamic correlation but are not relaxed in the presence of dynamic correlation. Colloquially,
dynamic correlation sees static correlation but static correlation
does not see dynamic correlation. This is true even for their multi-state (MS) variants<cit.>
because of their insufficient revision of the coefficients C̅_R k^(0) of any Ψ^(0)_k.
In contrast, the DS family of methods incorporate dynamic correction to each reference function Φ_R and then
construct and diagonalize an effective Hamiltonian in the space P, thereby producing a wave function that can be (very) different from Ψ^(0)_k
obtained by diagonalizing the bare Hamiltonian in the same space. Classic examples
are the (shifted) B_k type of methods<cit.>.
Since the zero-order coefficients C̅_R k^(0) are not involved in the dynamic correlation step at all,
it can be said that dynamic correlation does not see static correlation in such family of methods.
It should be clear that both the SD and DS families of methods cannot achieve balanced treatments of the static and dynamic components of
the overall correlation, especially when the two components are strongly entangled and even interchangeable.

This situation warrants a SDS type of treatment, where the predetermined coefficients C̅_R k^(0) are used
in dynamic correlation but are then sufficiently or even fully relaxed in the presence of dynamic correlation.
That is, the static and dynamic components of the overall correlation do see each other.
Several variants of multireference second-order perturbation theory (MRPT2)<cit.>,
internally contracted multireference configuration interaction (ic-MRCI)
<cit.>,
and internally contracted multireference coupled-cluster methods
<cit.>
belong to this family. However, such methods usually require the construction and diagonalization of a large Hamiltonian matrix
even just for one state.






At variance with this, one of the present authors proposed<cit.> a restricted SDS framework for
constructing the many-electron wave functions, where no matter how many electrons and how many orbitals are to be correlated,
only a 3N_P-by-3N_P Hamiltonian matrix is constructed and diagonalized for N_P states.
The framework leads to a series of methods, including SDSCI<cit.>, SDSPT2<cit.>, iterative configuration interaction (iCI)<cit.>,
iCI with configuration and perturbation (iCIPT2)<cit.>, and extended variants<cit.> of SDSCI and SDSPT2. Albeit a minimal MRCI, SDSCI
is very close in accuracy to ic-MRCI<cit.>, with a computational cost being only that
of one iteration of ic-MRCI. As an approximation to SDSCI, SDSPT2 is a CI-like MRPT2
and treats single and multiple roots in the same way.
In particular, SDSPT2 gives rise to MS-NEVPT2 for free. While SDSPT2 is usually very similar to MS-NEVPT2 in accuracy
<cit.>, it does outperform MS-NEVPT2
for situations with multiple nearly degenerate states<cit.> .
As an iterative version of SDSCI, iCI is an exact solver of full CI, whereas iCIPT2 is one of the most efficient
near-exact methods<cit.>. Taking iCI as the CASCI solver, we obtain iCISCF<cit.>, which can handle
active spaces as large as CAS(60,60) (i.e, 60 electrons in 60 orbitals).










One major problem associated with the above MR methods lies in that a large reference space P
leads to an exceedingly large first-order interacting space (FOIS) Q  that is intractable.
The only way to go is to reduce the reference space from P to P̃,
so as to reduce the FOIS from Q to Q̃.
The Q̃ space can further be treated approximately by separating it into an important subset
that is treated rigorously and an unimportant subset that can be treated approximately.
Several approaches have been proposed along this line
<cit.>, which differ from
each other in the choice of reduced reference space P̃, zeroth-order Hamiltonian H_0, perturbers spanning Q̃,
and effective Hamiltonian to be diagonalized for the final solutions.
In the present work, we apply such approximations to the
complete active space self-consistent field (CASSCF)-based SDSPT2<cit.>, so as to render it
applicable to systems that require very large active spaces.

The paper is organized as follows. The essential features of SDSPT2<cit.> are first recapitulated in Sec. <ref>. The
hole-particle symmetry-based graphic unitary group approach (HPS-GUGA)<cit.> is then employed in Sec. <ref> to construct
the spin-adapted, internally contracted configurations (ICC), starting with the orbital configurations (oCFG) contained
in the iCISCF/CASSCF<cit.> wave functions Ψ^(0)_k (=∑_R^N_RΦ_RC̅_R k^(0)).
Such ICCs span the full FOIS Q. However, the ICCs belonging to the Q-Q̃ portion of Q and hence having zero Hamiltonian matrix elements
with the reduced reference wave functions
Ψ̃^(0)_k=∑_R=1^Ñ_RΦ_RC̃̅̃_R k^(0) (resulting from the truncation of Ψ^(0)_k)
can be screened out automatically. An integral-based cutoff is further introduced to truncate the most expensive
subspaces of Q̃
The efficacy of SDSPT2s (SDSPT2 with selection) is illustrated in Sec. <ref> with
several challenging systems, including a simplified model of heme, chromium dimer, Cu_2O^2+_2 core, and transition metal complex [Co(TC-3,3)(NO)]. The paper is closed with a summary in Sec. <ref>.










§ SDSPT2

Unless otherwise stated, the notations documented in Table <ref> are to be used for the orbitals and states, under the Einstein summation
convention over repeated indices.

The SDS framework<cit.> starts with the following wave functions for N_P lowest states

    |Ψ_I⟩
    =∑_k^N_P |Ψ_k^(0)⟩C̃_kI + ∑_k^N_P|Ψ_k^(1)⟩C̃_(k+N_P)I + ∑_k^N_P|Ψ_k^(2)⟩C̃_(k+2N_P)I,

where Ψ_k^(0) (=∑^N_R_R=1Φ_RC̅_Rk^(0)), Ψ_k^(1), and Ψ_k^(2) represent the zeroth-order, first-order, and secondary functions, respectively. With the introduction of the primary (P_m) and secondary (P_s) parts of the P space,

    P_m   =   ∑_k=1^N_P|Ψ^(0)_k⟩⟨Ψ_k^(0)|,
    
    P_s   =   P-P_m
       =   ∑^N_R_R=1|Φ_R⟩⟨Φ_R|-P_m=∑_l=N_P+1^N_R|Ψ^(0)_l⟩⟨Ψ_l^(0)|,

the first-order and secondary functions can be defined as

    |Ψ^(1)_k⟩   =    Q1/E_k^(0)-H_0QH|Ψ^(0)_k⟩=∑_q∈ Q|Φ̅_q⟩C̅^(1)_qk,
    
    Q    =    1-P=∑_q∈ Q|Φ̅_q⟩⟨Φ̅_q|,
    
    |Ψ^(2)_k⟩   =    P_sH|Ψ^(1)_k⟩
       ≈    P_s^' H|Ψ^(1)_k⟩,   P_s^'=∑_l=N_P+1^N_P+M_P|Ψ^(0)_l⟩⟨Ψ_l^(0)|,
       =   ∑_R=1^N_R|Φ_R⟩C̅^(2)_Rk,

{Φ̅_q} in Eq. (<ref>) are orthonormalized ICCs obtained by single and double excitations from Ψ_k^(0).
Eq. (<ref>) represents the projection of the Lanczos vector H|Ψ^(1)_k⟩ onto the secondary space,
which is further approximated to Eq. (<ref>) to simply the evaluation of the matrix elements.
It should be clear from Eq. (<ref>) that the secondary functions are linear combinations
of the reference CSFs, with the coefficients related to the first-order functions.
As such, they can facilitate the relaxation of the reference coefficients in the presence of dynamic correlation,
particularly for avoided crossings<cit.> or multiple quasi-degenerate states<cit.>.
The fact that |Ψ_k^(0)⟩, |Ψ_k^(1)⟩, and |Ψ_k^(2)⟩ have decreasing weights
in the wave function |Ψ_I⟩ justify the characterization of them as primary, external, and secondary states.
Since both |Ψ_k^(1)⟩ and |Ψ_k^(2)⟩ are specific to |Ψ_k^(0)⟩, the generalized eigenvalue problem

    𝐇̃𝐂̃=𝐒̃𝐂̃𝐄̃

for determining the expansion coefficients of |Ψ_I⟩ is only of dimension 3N_P. The Hamiltonian and metric matrices have the following structures

    H̃   =   [ P_mHP_m   P_mHQ P_mHP_s;   QHP_m     QHQ   QHP_s; P_sHP_m   P_sHQ P_sHP_s ],
       =   [         E_k^(0)δ_kl ⟨Ψ_k^(0)|H|Ψ_l^(1)⟩                   0; ⟨Ψ_l^(1)|H|Ψ_k^(0)⟩ ⟨Ψ_k^(1)|H|Ψ_l^(1)⟩ ⟨Ψ_k^(1)|H|Ψ_l^(2)⟩;                   0 ⟨Ψ_l^(2)|H|Ψ_k^(1)⟩ ⟨Ψ_k^(2)|H|Ψ_l^(2)⟩ ], k,l=1,⋯, N_P,
    S̃   =   [              δ_kl                 0                 0;                 0 ⟨Ψ_k^(1)|Ψ_l^(1)⟩                 0;                 0                 0 ⟨Ψ_k^(2)|Ψ_l^(2)⟩ ].


The above is nothing but a minimal MRCI (dubbed as SDSCI<cit.>). Taking SDSCI as the seed,
various methods can be derived<cit.>, among which the simplest variant
is SDSPT2, which amounts to replacing the QHQ block in Eq. (<ref>) with QH_0Q. Like NEVPT2<cit.>,
the Dyall CAS/A Hamiltonian<cit.> is adopted here for H_0:







To satisfy the zeroth-order Schödinger equation,

    H_0|Ψ_J^(0)⟩=E_J^(0)|Ψ_J^(0)⟩,⟨Ψ_I^(0)|Ψ_J^(0)⟩=δ_IJ, ,



    H_0=H_I^D + H_A^D, 
       H_I^D=∑_iϵ_iÊ_ii + ∑_aε_aÊ_aa +C_I^D,
       H_A^D=∑_tu f^c_tuÊ_tu+1/2∑_tuvw(tu|vw)ê_tu,vw,
       C_I^D= E^c-2∑_iε_i, E^c=∑_i(h_ii+f^c_ii), 
       f^c_pq=h_pq+∑_i[2(pq|ii)-(pi|iq)].

where the quasi-canonical orbital energies, ϵ_i and ϵ_a,
are obtained by diagonalizing the generalized Fock matrix

    F_pq=f^c_pq+∑_tu[(pq|tu)-1/2 (pu|tq)]D_tu, 
       D_tu=∑_k w_k⟨Ψ_k^(0)|Ê_tu|Ψ_k^(0)⟩

for the doubly occupied and virtual subspaces separately.
The particular choice of the constant C^D_I (<ref>) is to
make ⟨Φ_μ|H^D|Φ_ν⟩ equal to ⟨Φ_μ|H|Φ_ν⟩  for all CSFs
in P, such that

    ⟨Ψ_k^(0)|H^D|Ψ_l^(0)⟩=⟨Ψ_k^(0)|H|Ψ_l^(0)⟩=E_k^(0)δ_kl,

where

    E_k^(0)=E^c+∑_tuf^c_tu⟨Ψ^(0)_k|E^t_u|Ψ^(0)_k⟩+1/2∑_tuvw(tu|vw)⟨Ψ^(0)_k|E^tv_uw|Ψ^(0)_k⟩.












Since only the P_mHP_m and P_mHQ blocks of Eq. (<ref>) are required by NEVPT2, the free production of NEVPT2 by SDSPT2 is obvious.
However, it should be noted that SDSPT2 is not size consistent. Nevertheless, the size consistency errors
can readily be cured by the Pople correction<cit.>, as demonstrated before<cit.>.



  Notations for the orbitals and states used in this work.
  =20pt
  
   Space         Indices

      Orbital space     

   Arbitrary orbital        p,q,r,s 

   Hole ( or closed) orbital        i,j,k,l 

   Active orbital        u,v,t,w 

   External ( or virtual) orbital          a,b,c,d 

               Configuration space     

   Step vector in GUGA        |(d)_μ⟩=|d_1,d_2,⋯,d_n⟩ 

   Configuration state function (CSF)         Φ_μ,Φ_ν 

   Electronic state        Ψ_I,Ψ_J 

   Reference CSF        Φ_R 

   Reference orbital configuration        Φ_R̅^oCFG 

   Reference wavefunction        Ψ^(0)_k,Ψ^(0)_l 

   A sub-DRT in the active space formed by all vertices and arcs in DRT from vertex X̅ to Y
          X̅Y 

   




§ IMPLEMENTATION OF SDSPT2S

The SDSPT2 method can be divided into two categories based on the choice of reference space. As illustrated in Fig. <ref>, SDSPT2 with a CAS reference has been implemented using HPS-GUGA directly in our previous work<cit.>, where the molecular orbitals are optimized through CASSCF. For systems with a large active space as the reference, SDSPT2 involves a vast number of Φ_R and Φ_q for the P and Q spaces, respectively, raising several questions:

(1) How should sufficient reference configurations be selected?

(2) How can an affordable Q space be constructed?

(3) How can redundant Φ_q be eliminated?

(4) How can negligible Φ_q and Φ̅_q be pruned?

(5) How can the corresponding H̃ and S̃ be calculated?

For the second category, iCISCF calculations with the C_min parameter are performed to generate sufficiently accurate Ψ_k^(0) and molecular orbitals as preliminary steps for SDSPT2. However, in practice, the number of Φ_R for Ψ_k^(0) is too large to be used as reference CSFs for SDSPT2 with large active spaces. Therefore, the reference CSFs that make up the P space must be chosen based on the coefficients of the iCISCF expansion, subject to the condition

    min(|C̅_Rk^(0)|) > P_min,  k ∈ N_P,

where P_min is the threshold for selecting the dominant Φ_R among all N_P Ψ_k^(0). It should be noted that this selection process might result in a scenario where all the chosen reference CSFs have no electron occupation on certain active orbitals, rendering these active orbitals irrelevant within the selected active space. In such cases, C_min should be adjusted to a smaller value to ensure the selection of reference CSFs that encompass all active orbitals.

In previous works<cit.>, the construction of the Q space for large active spaces is typically accomplished by exciting the reference orbital configurations (Φ_R̅^oCFG), which are derived from the CSFs of the P space. Given that one Φ_R̅^oCFG corresponds to several Φ_R with the same occupation pattern, the Q space constructed from the single and double excitations of selected Φ_R̅^oCFG includes many redundant Φ_q that do not interact with the chosen Φ_R. To remove the Q_r space containing these redundant Φ_q, the effective First Order Interacting Space (eFOIS) method was employed in our implementation. This approach is beneficial because it allows us to discard Φ_q when both ⟨Φ_q|Ê_pq|Ψ_k^(0)⟩=0 and ⟨Φ_q|ê_pq,rs|Ψ_k^(0)⟩=0 for all N_P target states, where Ê_pq and ê_pq,rs represent one- and two-electron excitation operators, respectively.





  §.§.§ CI space construction with selection

As is well known, dealing with the Q space for large active spaces as references is complex. Therefore, efficiently constructing and compactly representing the Q space is one of the main objectives of the selection-based SDSPT2 in this study. One method to address these challenges is the graphical unitary group approach (GUGA)<cit.>, which offers a compact structure of CSFs through distinct row tableaux (DRT) as proposed by Shavitt<cit.>.

The DRT is comprised of nodes (or vertices) and arcs (or sloped line segments). It represents all the CSFs of the full-CI method by specifying only the electron number (N) and orbital number (n), along with a defined total spin value S for the target state. Each node (a_r,b_r) signifies a unique row, while an arc labeled d_r shows the step from (a_r,b_r) downward to the adjacent node (a_r-1,b_r-1). The value of d_r is calculated as 3a_r+b_r, where x_r=x_r-x_r-1 for x=a,b. In this context, a_r=N_r/2-S_r, b_r=2S_r, with N_r representing the electron count in the first r occupied orbitals, and S_r being the intermediate spin value derived from the coupling of the first r orbitals. The possible values for d_r are 0, 1, 2, or 3. A d_r value of 0, depicted by a vertical line (arc), indicates that the r-th orbital is unoccupied. Values of d_r equal to 1, 2, or 3, represented by increasing slopes (arcs), correspond to positive spin coupling (spin up), negative spin coupling (spin down), and doubly occupied r respectively. The step from (a_r,b_r) to (a_r-1,b_r-1) through different d_r follows specific rules:

(1) d_r=0, (a_r-1,b_r-1)=(a_r,b_r);

(2) d_r=1, (a_r-1,b_r-1)=(a_r,b_r-1);

(3) d_r=2, (a_r-1,b_r-1)=(a_r-1,b_r+1);

(4) d_r=3, (a_r-1,b_r-1)=(a_r-1,b_r);

A DRT is constructed by connecting all possible nodes from the top (head) vertex (a_n,b_n) to the bottom (tail) vertex (0,0), ensuring that these connections adhere to the specified rules. The inverse sequence of d_r, ranging from the first to the n-th orbital, forms a step vector |(d)_μ⟩, which is used to sequentially document each CSF (|Φ_μ⟩).

The vast number of CSFs impedes Full-CI's capability to compute systems with numerous electrons and orbitals, despite the use of DRT to streamline the representation of the Full-CI space comprising these CSFs. The conventional technique to compress the CI space involves selecting reference CSFs and subsequently generating only their singly and doubly excited CSFs, thereby eliminating all potential higher-level excitations.
This truncation approach facilitates the construction of CSFs for the multi-reference configuration interaction with single and double excitations (MRCISD).<cit.>



In the MRCISD method, orbital space is typically partitioned into three categories: hole (closed or inactive) space, active space, and external (virtual) space. Within this framework, the reference CSFs are characterized by double occupancy in the hole space, variable occupancy in the active space, and no occupancy in the external space.
Accordingly, a DRT can be segmented into hole, active and external components as illustrated in Fig. <ref>. In this context, n_e, n_a, and n_h represent the counts of the external, active, and hole orbitals, respectively. Notably, the HPS-GUGA employs an inverse orbital order, as highlighted previously<cit.>.
In MRCISD, excitations from the hole space relative to the reference CSFs are limited to moving only one or two electrons, resulting in one or two holes in the hole space. Consequently, the external space is capable of accommodating only one or two additional electrons.

Fig. <ref> illustrates that the hole space graph is symmetrically matched with the external space graph, demonstrating hole-particle symmetry (HPS). This DRT contains two boundaries: one between the hole and active spaces, and the other between the active and external spaces. For S=0, the boundary vertices between the hole and active spaces are defined by S̅ for the CI subspace with two singlet holes, T̅ with two triplet holes, D̅ with one hole, and V̅ with no hole. In the case of S=1/2, D̅ is replaced by D̅_1/2 for spin up and D̅_-1/2 for spin down electrons on the hole space. For S ≥ 1, an additional vertex T̅_-1 appears for S_z=-1 on the hole space, while S̅ and T̅ represent S_z=0 and S_z=1, respectively. This results in six vertices: T̅_-1, D̅_-1/2, S̅, T̅, D̅ (or D̅_1/2), V̅ on the boundary between the hole and active spaces. The corresponding (a_r,b_r) values for these vertices are calculated as follows:
T̅_-1: (a-n_h+2,b-2),
D̅_-1/2: (a-n_h+1,b-1),
S̅: (a-n_h+1,b),
T̅: (a-n_h,b+2),
D̅: (a-n_h,b+1),
V̅: (a-n_h,b),
where a = N/2 - S and b = 2S, with N being the total number of electrons.
Conversely, the boundary between the active and external spaces comprises four vertices: S, T, D, and V, associated with (1,0), (0,2), (0,1), and (0,0), respectively. These correspond to the aforementioned S̅, T̅, D̅, and V̅ in the HPS-GUGA representation, as depicted in Fig. <ref>.

The whole CI space of uncontracted MRCISD (uc-MRCISD), as represented by the DRT in Fig. <ref>, can be segmented into several subspaces, each represented by sub-DRTs within the HPS-GUGA framework. Each sub-DRT corresponds to a collection of arcs on the active space, bounded above by X̅ and below by Y. Consequently, a sequence of |(d)_μ⟩ originating from the head vertex, passing through the X̅Y pair, and terminating at the tail vertex constitutes a CI subspace denoted as X̅Y.
The number of steps from X̅ to the head vertex is termed the number of up steps (NUS(X̅)), while the steps from Y to the tail vertex represent the number of down steps (NDS(Y)). The quantity D(X̅Y) represents the sum of all possible steps between X̅ and Y. Thus, the entire CI space can be categorized by a sequence of X̅Y pairs, with the total dimension expressed as:

    Dim(uc-MRCISD)=∑_X̅,YNUS(X̅) × D(X̅Y) × NDS(Y).

In this formulation, all CSFs are naturally orthonormal. Each X̅Y pair signifies a specific excitation type, as illustrated in Table <ref>. This table also highlights the correspondence between conventional CI subspace notations<cit.> and those used in HPS-GUGA under the column related to S_l^(k).


  CI subspace correspondence between S_l^(k) and sub-DRTa.
  =20pt
  t]@ccccccccccccc@
   Excitation operator (Ê^X̅Y_M)     n_hb     n_ec     S_l^(k)     sub-DRT(X̅Y)d   

   Ê_uv,ê_uv,tw     0     0     S^(0)     V̅V 

Ê_ui, ê_ui,vw      1     0     S_i^(1)       D̅V

Ê_au, ê_au,vw      0     1     S_a^(-1)      V̅D

Ê_ai, ê_ai,uv, ê_ui,av      1     1     S_a,i^(0)     D̅D

ê_ui,vj      2     0     S_ij^(2)         P̅V

ê_au,bv      0     2     S_ab^(-2)        V̅P

ê_ai,uj      2     1     S_a,ij^(1)       P̅D

ê_ai,bu      1     2     S_ab,i^(-1)      D̅P

ê_ai,bj      2     2     S_ab,ij^(0)      P̅P

   
  
  
  a S_l^(k) denotes CI subspace on ref. <cit.>, where k is the difference between the numbers of holes and particles.
  
  b Number of holes on the hole space.
  
  c Number of electrons on the virtual space.
  
  d P = S or T, S: two-electron singlet; T: two-electron triplet; D: one-electron doublet; V: void.
  


The example of water depicted in Fig. <ref> demonstrates that a mere 32 nodes and their connecting arcs suffice to characterize the 259 CSFs encompassing the total CI space, showcasing the compactness of the representation. Henceforth, the efficient identification of these nodes becomes paramount for the construction of the CI space (P and Q) in this study.
Given the structured and straightforward DRT configurations of the hole and external subspaces, the initial eight nodes pertaining to the hole space are predetermined. Consequently, the terminal eight nodes associated with the external space can also be determined effortlessly without the need for a search.
The search for nodes is specifically required only within each sub-DRT (X̅Y), focusing on the region between X̅ and Y, guided by the reference in the active space.

When the CAS comprises 6 orbitals, either Φ_R or Φ_R̅^oCFG, searching for nodes between each pair of X̅ and Y becomes straightforward. This ease arises because all possible excitations are confined to the active orbitals relative to each Φ_R within the CAS. The variations among CI subspaces (or X̅Y pairs) primarily stem from differences in hole and particle numbers, as elucidated in Table <ref>. Each sub-DRT, extending from X̅ to Y, can be conceived as a miniature DRT equivalent to a Full-CI calculation, per the discussion above. Notably, all sub-DRTs possess identical nodes on the active space; for example, nodes 9∼24 are shared in Fig. <ref>. This indicates that all Φ_R are utilized collectively in the construction of the CI space. It is crucial to highlight that nodes situated on the (n_a-1)-th row are discarded if they do not connect downwards to any Y on the n_a-th row. This removal criterion stems from the fact that the step vectors passing through these nodes represent excitations exceeding double excitations from all Φ_R in the CAS.

In the context of a large active space with specific reference configurations, not all nodes and arcs within each sub-DRT, bounded by X̅ and Y, are permissible. This restriction arises because certain Φ_R within the CAS are omitted as reference CSFs. Consequently, nodes and arcs associated with excited CSFs originating from these unselected Φ_R̅^oCFG (or Φ_R) must be eliminated. This process requires an individual examination of single and double excitations from each Φ_R̅^oCFG. For instance, Fig. <ref>(a) shows that for the sub-DRT of D̅D, 8 nodes generate 5 excited CSFs when one Φ_R̅^oCFG is used. When two Φ_R̅^oCFG are chosen, there are 9 nodes producing 8 excited CSFs, as depicted in Fig. <ref>(b). Moreover, selecting three Φ_R̅^oCFG results in 10 nodes and 9 excited CSFs, illustrated in Fig. <ref>(c). Notably, the graph connecting X̅ and Y in Fig. <ref>(c) mirrors that of the entire sub-DRT between X̅ and Y in Fig. <ref>. This similarity indicates that selecting three Φ_R̅^oCFG from the CAS can replicate the same CI subspace generated by employing all Φ_R̅^oCFG of the CAS. Concerning the selected Φ_R, certain excited CSFs become redundant and can be effectively pruned using the eFOIS method mentioned earlier.



The process of constructing the CI space using selected Φ_R̅^oCFG is evidently more complex than when the entire CAS serves as the reference. This complexity necessitates the application of two specific restrictions to eliminate nodes and arcs within the sub-DRTs pertaining to the selected configurations. Table <ref> introduces two auxiliary indices related to these restrictions. The first index, T^ex_X̅Y, corresponds to T_R̅^r,μ, which is designed to track the cumulative excitation number (a_r,b_r)_μ in relation to Φ_R̅^oCFG. Here, r signifies the r-th active orbital, and μ represents the index of node (a_r,b_r)_μ, as depicted in Fig. <ref>. Various methods can define T_R̅^r,μ, including a comparison of the occupation numbers of CSFs that pass through (a_r,b_r)_μ with those of Φ_R̅^oCFG in the r-th orbital. The definition of T_R̅^r,μ can be mathematically expressed as:

    T_R̅^r,μ=∑_v=1^r(N^occ_v,μ-N^occ_v,R̅), when N^occ_v,μ≥ N^occ_v,R̅,

where N^occ_v,μ ranges from 0 to 2, N^occ_v,R̅ also ranges from 0 to 2, v spans from 1 to n_act, and R̅ ranges from 1 to n_ref^oCFG. Here, n_ref^oCFG represents the total count of Φ_R̅^oCFG, N^occ_v,R̅ denotes the occupation number of Φ_R̅^oCFG on the v-th row, and N^occ_v,μ signifies the occupation number of the v-th orbital following electron excitation from this orbital, considering possible excitation types d_v of 0, 1, 2, or 3. Consequently, T_R̅^r,μ serves to count the number of excited electrons from X̅ to (a_r,b_r)_μ with the Φ_R̅^oCFG. Should T_R̅^r,μ for all Φ_R̅^oCFG surpass T^ex_X̅Y for a given sub-DRT, such as V̅D, and the condition outlined below is satisfied:

    min(T_R̅^r,μ)>T^ex_X̅Y,

then the node (a_r,b_r)_μ and its downward linking arcs to Y are removed.

The second auxiliary index is the occupation index N^ex_X̅Y, which pertains to the cumulative occupation number of Φ_R̅^oCFG on the r-th row. This relationship is defined by:

    N^sum_r,R̅=N^occ_r,R̅+N^sum_r-1,R̅, N^sum_0,R̅=0,

where N^sum_r,R̅ represents the summation of occupation numbers up to the r-th row for Φ_R̅^oCFG. Consequently, the remaining active electron number (N^re_r,R̅) up to the r-th orbital of Φ_R̅^oCFG is given by:

    N^re_r,R̅=N_a-N^sum_r,R̅.

In this study, V̅V, composed of Φ_R, serves as the P space in selection-based SDSPT2. This choice stems from the observation that the CI space derived from iCISCF(2) can effectively serve as the P space for subsequent SDSPT2 calculations when the static correlation computed by iCISCF(2) closely matches that obtained through CASSCF<cit.>. Each Q subspace, generated by single and double excitations from Φ_R̅^oCFG, exhibits a characteristic where N̅^r_μ=(2a_r+b_r)_μ does not exceed the value of V̅V without excitation from Φ_R̅^oCFG by more than two. This discrepancy is captured by N^ex_X̅Y as outlined in Table <ref>. Thus, the adjusted remainder of active electrons in Φ_R̅^oCFG is expressed as:

    N_R̅^r=N^re_r,R̅+N^ex_X̅Y.

Should the candidate node (a_r,b_r)_μ meet the condition:

    N̅^r_μ=(2a_r+b_r)_μ>max(N_R̅^r),

it must be eliminated, as no electrons remain available for excitation from any Φ_R̅^oCFG. With both restrictions, as defined by Eqs.(<ref>) and (<ref>), the viable nodes and their linking relationships from X̅ to Y can be efficiently determined.


  CI subspace constructions by auxiliary excitation and occupation indices relative to the selected orbital configurations,
  which are labeled by the number of parenthesis of (T^ex_X̅Y,N^ex_X̅Y), respectively.a
  =40pt
  t]@ccccccccccccc@
   (T^ex_X̅Y,N^ex_X̅Y)     V     D     P  

   V̅     (0,0)           (1,2)     (0,2) 

   D̅     (2,2)           (1,2)     (0,2) 

   P̅     (2,2)           (1,2)     (0,2) 

   
  
  
  a P = S or T, S: two-electron singlet; T: two-electron triplet; D: one-electron doublet; V: void. P̅ = S̅ or T̅, S̅: two-hole singlet; T̅: two-hole triplet; D̅: one-hole doublet; V̅: void.
  


For the illustrated example shown in Fig. <ref>(a), when a single Φ_1^oCFG is selected, the node (2,1) for r=1 is eliminated because N̅^1_μ=5, which exceeds N_R̅^r=N_1^1=4 (calculated as 2+2), thus satisfying the condition in Eq.(<ref>). Similarly, the node (1,1) for r=2 is removed since N̅^2_μ=3 surpasses N_R̅^r=N_1^2=2. In the case of Fig. <ref>(b), where two Φ_R̅^oCFG are selected, the node (1,1) for r=2 remains because N̅^2_μ = N_R̅^r = N_2^2 = 3. Furthermore, with three selected Φ_R̅^oCFG as in Fig. <ref>(c), the node (2,1) for r=1 is retained because N̅^1_μ = N_R̅^r = N_3^1 = 5. It is worth noting that the condition specified in Eq.(<ref>) is not met in this example due to the limited active space.

It is important to note that no duplicate nodes exist across all sub-DRTs, as each newly introduced node (a_r,b_r)_μ must be verified as unique by comparison with the existing nodes in the r-th orbital. This verification process is a time-consuming aspect of the selection-based construction of sub-DRTs. For every node, a quadruple (r,a_r,b_r,λ_r)_μ and T_R̅^r,μ corresponding to n_ref^oCFG must be stored. Here, λ_r denotes the irreducible representation associated with orbital symmetry. Each (r,a_r,b_r,λ_r)_μ can be stored using a single 64-bit integer for molecular systems where n_a ≤ 255 and N_a/2-S, since four integers ranging from 0 to 255 can be encoded within a 64-bit integer using bitwise operations. Furthermore, each T_R̅^r,μ can be represented using 2 bits, as its values are constrained to 0, 1, 2, or 3. The value of 3 is significant because in MRCISD, all high excitations lead to T_R̅^r,μ≥ 3, allowing one 64-bit integer to store 32 different T_R̅^r,μ values. Indeed, two nodes that share the same (r,a_r,b_r,λ_r)_μ but have different T_R̅^r,μ values are considered distinct, not to mention those with differing (r,a_r,b_r,λ_r)_μ.

The conditions and bit representation discussed above facilitate a straightforward and efficient algorithm to address the aforementioned question (2). We propose the following steps for this algorithm to search and verify new nodes in each sub-DRT:

(1) Establish the upper boundary vertex X̅.

(2) Search for (a_r,b_r)_μ connecting from the directly preceding active orbital using one of d_r=0, 1, 2, 3, adhering to the aforementioned rules of DRT, and check if the resulting node will lead to higher excitation using the conditions in Eqs.(<ref>) and (<ref>).

(3) It is only necessary to check for duplication on the same orbital, as nodes with different orbital indices are inherently distinct. Therefore, when there is more than one node on the r-th orbital, compare the (r,a_r,b_r,λ_r)_μ of the new node with the existing (r,a_r,b_r,λ_r)_ν. If (r,a_r,b_r,λ_r)_μ≠ (r,a_r,b_r,λ_r)_ν, the new node remains.

(4) If (r,a_r,b_r,λ_r)_μ = (r,a_r,b_r,λ_r)_ν but ∑_R̅T_R̅^r,μ≠∑_R̅T_R̅^r,ν, it represents a new node.

(5) If ∑_R̅T_R̅^r,μ = ∑_R̅T_R̅^r,ν, further examination of each pair of integers with 32 (T_R̅^r,μ, T_R̅^r,ν) is required. If all integer pairs are identical, the node is a duplicate; however, any differing integer pair results in the retention of the candidate node.

(6) When the new node is located on the last active orbital n_a, it should be removed if its (a_n_a,b_n_a)_μ differs from that of the lower boundary vertex Y.

In this approach, each node between X̅ and Y is distinct, facilitating the creation of orthonormally excited CSFs from Φ_R̅^oCFG within the Q subspace of X̅Y. To store a node with the necessary information (r,a_r,b_r,λ_r)_μ, only one 64-bit integer is required. Furthermore, the indices of the next four connected nodes corresponding to d_r = 0, 1, 2, 3 are saved using four additional 64-bit integers. This efficient storage mechanism ensures that each node can be uniquely identified and linked to its adjacent nodes, facilitating the navigation and processing within the defined Q subspace.
Given the significantly lower number of nodes compared to the generated CSFs, the memory required to save these CSFs in this study is minimal when contrasted with that used for bit representation. Moreover, the generation of new nodes and the checking for duplicates are performed concurrently, thereby eliminating the substantial computational load associated with duplicate checking of generated CSFs when using other techniques like iCI<cit.> and ICE<cit.>.






















  §.§.§ Evaluating matrix elements of SDSPT2

For the SDSPT2, it is necessary to perform an internal contraction of the excited CSFs in order to compact the Q space, which consequently reduces the number of combination coefficients of Ψ_k^(1) on Eq.(<ref>).
The internally contracted functions, denoted as, Φ^X̅Y_MI, are characterized as excited functions relative to the reference state |Ψ^(0)_I⟩, and are defined as follows:<cit.>

    Φ^X̅Y_MI   =Ê^X̅Y_M|Ψ^(0)_I⟩=∑_μ∈X̅Y|Φ_μ⟩⟨Φ_μ|Ê^X̅Y_M|Ψ^(0)_I⟩
       = ∑_μ∈X̅Y|Φ_μ⟩∑_R∈V̅V⟨Φ_μ|Ê^X̅Y_M|Φ_R⟩ C^I_R
       =∑_μ∈X̅Y|Φ_μ⟩C̃^I_μ M.

The spin-free excitation operator, Ê^X̅Y_M, represents either Ê_pq or ê_pq,rs, where M is a collective index that specifies the active orbital indices participating in these excitation operators. The target state index is denoted by I, and the specific definitions of Ê^X̅Y_M are provided in Table <ref> for distinct sub-DRTs (or CI subspaces). Furthermore, the contraction function can be described as a linear combination of the complete set of Φ_μ functions within X̅Y, as evidenced by the identity ∑_μ∈X̅Y|Φ_μ⟩⟨Φ_μ|=1.

In the HPS-GUGA representation, as described by Wang et al.<cit.>, |Φ_μ⟩ is defined as a step vector |(d)_μ⟩, which can be further detailed as |((d)_h(d)_a(d)_e)_μ⟩. Here, (d)_h, (d)_a, and (d)_e represent the step values for the hole, active, and external spaces, respectively.
The coupling coefficient ⟨(d)_μ|Ê^X̅Y_M|(d)_ν⟩ can be rewritten as:

    ⟨(d)_μ|Ê^X̅Y_M|(d)_ν⟩   =⟨ ((d)_e(d)_a(d)_h)_μ|Ê^X̅Y_M| ((d)_h(d)_a(d)_e)_ν⟩
       =∑_J=0,1ω_J HLS(X̅X̅^')_μν·  ALS(X̅Y,X̅^'Y^')_μν· ELS(YY^')_μν,

Here, |(d)_μ⟩ and |(d)_ν⟩ belong to sub-DRTs X̅Y and X̅^'Y^', respectively.
The one-electron coupling coefficient is simplified compared to the two-electron case. For the one-electron coupling coefficient, the summation and index J are not needed, and the values are fixed at J=0 and ω_J=1. This means that for one-electron coupling, there is no spin coupling to consider, and the factor ω_J is always 1, indicating no exchange or direct type coupling complexity.
For the two-electron coupling coefficient, the index J accounts for the spin coupling between two generators, E_pq and E_rs, each having a spin of 1/2. The factor ω_J is determined by the intersection of the two generator lines: ω_J=1 for exchange type coupling and ω_J=-1 for direct type coupling. If E_pq and E_rs have no orbital overlap, then ω_J=0.

The coupling coefficients are reexpressed as a product of three segmental factors: HLS, ALS, and ELS, which are defined as follows:

   HLS(X̅X̅^')_μν=∏^n_h_r=1W(Q_r(X̅X̅^');(d_r)_μ(d_r)_ν,▵b_r,(b_r)_ν,J) 

   ALS(X̅Y,X̅^'Y^')_μν=∏^n_h+n_a_r=n_h+1W(Q_r(X̅Y,X̅^'Y^');(d_r)_μ(d_r)_ν,
▵b_r,(b_r)_ν,J) 

   ELS(YY^')_μν=∏^n_r=n_h+n_a+1W(Q_r(YY^');(d_r)_μ(d_r)_ν,▵b_r,(b_r)_ν,J)

Here, W(Q_r;(d_r)_μ(d_r)_ν,▵ b_r,b_r,J) represents the segment factors that depend on the segment type Q_r, with ▵b_r=(b_r)_ν-(b_r)_μ. The detailed descriptions of these types and their corresponding segment factors are provided by Paldus, Boyle, and Payne <cit.>.
HLS(X̅X̅^')_μν, termed as the hole loop shape, represents a graphical illustration commencing from the loop head to X̅. This may manifest as either a complete or partial loop, contingent upon whether the loop is closed in the hole space. Similarly, ELS(YY^')_μν signifies the external (partial) loop shape corresponding to the graph extending from Y to the loop tail.
These two types of shapes are associated with specific formulae, enabling their pre-calculation as described by Wang et al. <cit.>. In contrast, the active partial and complete loops, denoted as ALS(X̅Y,X̅^'Y^')_μν, necessitate determination through loop searching originating from X̅ to Y. It is noteworthy that confining the search to the active space results in fewer partial loops compared to implementations not utilizing the HPS-GUGA approach. This reduction not only conserves computational time but also mitigates computing bottlenecks. Moreover, multiple ALS(X̅Y,X̅^'Y^')_μν possessing identical active orbital indexes can be simultaneously searched. Subsequently, these can be multiplied in parallel with block-stored molecular integrals sharing the same orbital indexes, thereby expediting the evaluation of the Hamiltonian matrix.


The contraction coefficient, C̃^I_μ M on Eq.(<ref>), can be derived from the interaction between Φ_μ and the reference wavefunction. This coefficient can alternatively be expressed as

    C̃^I_μ M   =∑_R∈V̅V⟨Φ_μ|Ê^X̅Y_M|Φ_R⟩ C^I_R
       =∑_R∈V̅V⟨ (d)_μ|Ê^X̅Y_M| (d)_R⟩ C^I_R
       =∑_R∈V̅V⟨ ((d)_e(d)_a(d)_h)_μ|Ê^X̅Y_M| ((d)_h(d)_a(d)_e)_R⟩ C^I_R
       =ELS(YV) · A^I_μ M· HLS(X̅V̅)

where

    A^I_μ M=∑_R∈V̅VALS(X̅Y,V̅V)_μ R C^I_R.

HLS(X̅V̅) and ELS(YV) are precomputed and do not influence the C̃^I_μ M values. These coefficients pertain exclusively to A^I_μ M, which relates to the step vectors connecting X̅ to Y within the active space, involving Φ_μ and Φ_R. The linear dependence of the internally contracted functions Φ^X̅Y_MI is also determined by A^I_μ M.
For example, contraction functions Φ^D̅D_MI and Φ^D̅D_NI, which involve different orbitals u and v but share the same indices i and a, may be nonorthogonal. To eliminate this linear dependence, one can compute the overlap matrix S^I_MN = ∑_μν⟨ A^I_μ M | A^I_ν N⟩. Subsequently, Löwdin's symmetric orthogonalization<cit.> can be applied to this matrix to ensure orthogonality among the functions.









The generated Φ^X̅Y_MI serve as basis functions for the CI Hamiltonian matrix, e.g. Dyall Hamiltonian in Eq.(<ref>), allowing the coupling coefficients to be adapted as follows:

    ⟨Φ^X̅Y_MI|Ê^X̅Y_M|Φ^X̅^'Y^'_NJ⟩   =∑_J=0,1ω_J ∑_μ∈X̅Y,ν∈X̅^'Y^'C̃^I_μ M⟨Φ_μ|Ê^X̅Y_M|Φ_ν⟩C̃^J_ν N
       =∑_J=0,1ω_J ELS(YY^') · HLS(X̅X̅^') ·∑_μ∈X̅Y,
    ν∈X̅^'Y^'A^I_μ M· ALS(X̅Y,X̅^'Y^')_μν· A^J_ν N
       =∑_J=0,1ω_J ELS(YY^') · HLS(X̅X̅^') · [𝐀^I^†·ALS·𝐀^J]_MN,
       =∑_J=0,1ω_J ELS(YY^') · HLS(X̅X̅^') ·ALS^'_MN.

In this equation, ALS^' represents the assembled matrix of internally contracted partial loops, while 𝐀^I and 𝐀^J are matrices of contraction coefficients. The expression in Eq.(<ref>) differs from that in Eq.(<ref>) primarily by the linear transformation of ALS from uncontracted to internally contracted forms, mediated by A^I_μ M.

In the context of large active spaces with numerous active orbitals, internal contraction may become infeasible. Specifically, the two CI subspaces D̅V and V̅D contain a substantial number of Φ^X̅Y_MI (or Φ̅_q as discussed in subsection <ref>). This abundance leads to a computational bottleneck, necessitating the storage of A^I_μ M, where the largest data block encompasses approximately n_a^3 × D(V̅D) or n_a^3 × D(D̅V) 64-bit real numbers.
To mitigate memory consumption, we employ the FOIS method to prune Φ^D̅V_MI and Φ^V̅D_MI. For example, Φ_μ^D̅V does not contribute to Φ^D̅V_MI and can be pruned if it satisfies the condition:

    max|⟨Φ_μ^D̅V|(iu|vw)ê_iu,vw|Φ_R⟩ C^I_R|^R∈V̅V_i=1,n_h<Q_min,

where Q_min is the pruning threshold. This criterion is analogous to that used for pruning generated CSFs in methods such as SBCI<cit.> and iCI<cit.>.
Through this pruning procedure, many of Φ^D̅V_MI are filtered due to their minor contraction coefficients. Similarly, V̅D can be compressed using the specified Q_min. The number of Φ_μ^D̅V or Φ_μ^V̅D can be reduced by two orders of magnitude or more, further reducing an order of magnitude of Φ^D̅V_MI or Φ^V̅D_MI. This approach effectively addresses the issue mentioned in question (4).
Furthermore, the sparse nature of 𝐀^J significantly reduces both computational cost and memory requirements, resolving the issue presented in question (5).









































§ RESULTS AND DISCUSSION

The SDSPT2 algorithm described in the previous section has been implemented in a developing version of the Xi’an-CI module <cit.> within the BDF program package <cit.>.  Scalar relativistic effects were considered in all calculations using the spin-free exact two-component (sf-X2C) relativistic Hamiltonian <cit.> and the core orbitals were also frozen.

The multi-reference wave functions of all molecules were computed using iCAS <cit.> and iCISCF <cit.> with the default C_min threshold of 1.0 × 10^-4. In SDSPT2 calculations, the iCISCF reference wave function is further truncated based on thresholds P_min and Q_min to reduce computational costs for computing correlation energies. The final results are computed by adding up the SDSPT2 correlation energy and the reference energy by iCISCF(2) <cit.>. In each SDSPT2 calculation, multi-state NEVPT2 (abbreviated as NEVPT2) <cit.> result is obtained simultaneously as a side product. All calculations were performed on a computer equipped with 2 Intel(R) Xeon(R) Platinum 8375C CPUs (64 cores in total) and 1 TB of memory.



 §.§ Truncation parameter P_min and Q_min

In principle, the present SDSPT2 algorithm is flexible and could employ any MCSCF wave functions as reference. However, the computational costs of SDSPT2 for references with more than 10^6 CSFs are hardly manageable, whereas the iCISCF method can handle references with about 10^9 CSFs. Thus, the reference CSFs generated by iCISCF should be further truncated before SDSPT2 calculations using P_min. Additionally, in the SDSPT2 calculation, the CSFs in D̅V and V̅D subspaces are reduced by Q_min, as detailed in Section <ref>.

To determine the optimal values of the two parameters, a simplified model of heme (denoted Fe^IIL_2) <cit.> is used as
the model system, which is considered as a strongly correlated system<cit.>. The Cartesian coordinates of Fe^IIL_2 were taken from previous work by Radoń<cit.> and are provided in the supporting information. The iCISCF calculation with respect to three low-lying states, ^5A_g, ^3B_3g, and ^1A_g, is performed with an active space of CAS(14,17), which includes Fe 3d4s4d, two 4p orbitals within the molecular plane, as well as a frontier σ ligand orbital of nitrogens. The 1s of C, N, and Fe's 1s2s2p orbitals are frozen in SDSPT2 calculations. The def2-TZVP and the corresponding auxiliary basis sets are employed for the resolution of identity (RI) approximations in both iCISCF as well as SDSPT2 calculations.

The C_min parameter influences the accuracy of the iCISCF wave function. In previous work, we found that setting C_min to 1.0×10^-4 produced sufficiently accurate results and orbitals.<cit.> However, its impact on the accuracy of subsequent MRPT2 calculations has not been investigated. Therefore, before examining the effects of P_min and Q_min, we first test the influence of C_min on SDSPT2, using the settings P_min=C_min and Q_min=0.0. The SDSPT2 energies computed using references with different C_min values are provided in the SI. The results indicate that the energy changes for the selected three electronic states are only 0.1 kcal/mol when C_min≤ 10^-4. Hence, for the rest of the work, a C_min setting of 10^-4 is used in all calculations unless otherwise specified.

For calculations involving large active spaces, the number of CSFs generated by iCISCF can be too large to perform SDSPT2 calculations. Therefore, the size of the reference used in SDSPT2 is further reduced based on the magnitude of the CI coefficients of the iCISCF wave function. This refinement process utilizes the parameter P_min. Comparisons of both absolute and relative energies for various low-lying excited states of the Fe^IIL_2 complex, calculated employing different P_min thresholds, are provided in Table <ref>. By adopting P_min=10^-3, the resultant absolute energy deviations compared to the calculations with P_min=10^-4 are less than 10^-3 Hartree for all states. The deviations in excitation energies are less than 0.3 kcal/mol.  However, the numbers of reference CSFs in SDSPT2 calculations are reduced from 7204 to 2727 for the ^5A_g state.

It is well known that in internally contracted MRPT2 or MRCI calculations with large active space reference, the D̅V and V̅D subspaces pose computational bottlenecks. As discussed in Section <ref>, the truncation parameter Q_min could be introduced to reduce the number of functions in FOIS. The results for Fe^IIL_2 using P_min=10^-3 and different Q_min values are given in Table S4. These results indicate that by setting Q_min=10^-5, there is almost no loss in accuracy. The number of CSFs in  Φ^V̅D_q is reduced by an order of magnitude. Thus, to balance computational costs and accuracy, the settings of P_min=10^-3 and Q_min=10^-5 are used throughout the rest of the work.


  SDSPT2 energies (+1727.0 E_h) of Fe^IIL_2 with iCISCF(14,17) providing optimized orbitals and reference configurations and
  Q_min=0.0.
  =5pt
  t]@ccccccccccccc@
   P_min       number of Φ_R of ^5A_g      ^5A_g        ^3B_3g      ^1A_g      E(^5A_g → ^3B_3g) (kcal/mol)    
   E(^5A_g→ ^1A_g) (kcal/mol) 

   1.0×10^-3     2727     -0.506706         -0.500335        -0.445045      4.0 	   38.7 

   5.0×10^-4     5051     -0.505963         -0.499160  	     -0.444192      4.3 	   38.8 

   1.0×10^-4     7204     -0.505751         -0.499040        -0.444233      4.2 	   38.6 

   



  SDSPT2 energies (+1727.0 E_h) of Fe^IIL_2 with iCISCF(14,17) providing optimized orbitals and reference configurations and
  P_min=10^-3.
  =5pt
  t]@ccccccccccccc@
   Q_min      number of Φ^V̅D_q of ^5A_g      ^5A_g     ^3B_3g      ^1A_g      E(^5A_g → ^3B_3g) (kcal/mol)     E(^5A_g→ ^1A_g) (kcal/mol) 

   1.0×10^-3       58800     -0.502386        -0.492765 	     -0.438914      6.0 	   39.8 

   1.0×10^-4       71330     -0.505569        -0.498354 	     -0.443763      4.5 	   38.8 

   1.0×10^-5      214000     -0.506566        -0.500099          -0.444933      4.1 	   38.7 

   0.0                  1746717     -0.506706        -0.500335          -0.445045      4.0 	   38.7 

   


The energy gaps of Δ E(^5A_g→^3B_3g) and Δ E(^5A_g→^1A_g) of Fe^IIL_2 computed by SDSPT2, using these default thresholds, are 4.1 and 38.7 kcal/mol, respectively. These results are in good agreement with the 3.8 and 40.7 kcal/mol obtained from CCSD calculations by Radoń<cit.>.


 §.§ A test on computing resources

To analyze the computational costs associated with selection-based SDSPT2 relative to CAS-based SDSPT2, we conducted SDSPT2 calculations for the ground state of anthracene, both with and without the selection of CAS(14,14). Table <ref> presents the CPU time (in minute) for the three most time-consuming components of SDSPT2, corresponding to different values of P_min. The data reveals that constructing the Q space using CAS as a reference is highly efficient, yet the time cost escalates as P_min decreases for SDSPT2 with selection. Moreover, the time cost ratios for constructing QH_0Q and P_mHQ as defined in Eq.(<ref>) align closely with the ratios of the numbers of Φ_q and Φ_R, respectively.
It becomes evident that the primary computational task in selection-based SDSPT2 is the construction of QH_0Q, which requires approximately an order of magnitude more resources than P_mHQ. Conversely, in the non-selective SDSPT2 approach, P_mHQ emerges as the most resource-intensive component, especially when the number of Φ_R surges from P_min=1.0×10^-3 to P_min=0.0. The peak memory consumption occurs during the storage of C̃^I_μ M, as depicted in Eq.(<ref>). Notably, employing a reference selection with P_min=1.0×10^-3 reduces the memory usage to roughly one-tenth compared to using CAS as a reference. In conclusion, the selection-based SDSPT2 markedly diminishes the demand for computational resources, making it an advisable choice for multi-reference calculations involving large active spaces.


  Number of Φ_R and Φ_q, CPU time (in minute) and memory Usage (in GBtype) for the SDSPT2 calculations of the ground state of anthracene with and without selection of CAS(14,14).
  =5pt
  t]@ccccccccccccc@
   P_min     number of Φ_R     number of Φ_q     time for Q space      time for QH_0Q      time for P_mHQ     maximal used memory 

   1.0×10^-3    9638    3428202     4 	     46 	         4         1.1 

   1.0×10^-4    43826    14848192     50      273 	        28         2.3 

   0.0                 2760615    67113618     0 	    967         1429    	  11.2 

   





 §.§ Size Consistency

As a CI-like MRPT2 method, SDSPT2 inherently lacks size consistency. However, the size-consistency errors of SDSPT2 when using CAS as a reference can be effectively corrected through the Pople correction<cit.>, as previously validated by our research<cit.>. In this study, we also examined the errors associated with selection-based SDSPT2, by analyzing the energies of low-lying singlet and triplet states of anthracene...Rg complexes (where Rg represents rare gas atoms such as He, Ne, Ar, Kr) relative to isolated anthracene.
For these complexes, the rare gas atom was positioned 100 Å away from the center of anthracene, with the geometry of anthracene obtained from <cit.>. We selected CAS(14,14) that includes all π and π^* orbitals of anthracene and utilized SA3-iCISCF(14,14)/cc-pVDZ to optimize the orbitals and generate Φ_R.
Table <ref> demonstrates that the size-consistency errors in SDSPT2 applying the Pople correction remain within an acceptable range for both the ground and excited states of anthracene...Rg complexes, with the maximum error recorded being only 0.10 eV. It is also noteworthy that NEVPT2, when tested with selected reference CSFs from a large active space, exhibited fundamental size consistency, as evidenced in Table <ref>.
To assess the accuracy of selection-based SDSPT2 (iCISCF/SDSPT2), we additionally performed calculations using CAS-based SDSPT2 (CASSCF/SDSPT2) with C_min=P_min=Q_min=0.0, which is expected to yield more precise results. Table <ref> reveals that the size-consistency errors for iCISCF/SDSPT2 and CASSCF/SDSPT2, as well as for the analogous iCISCF/NEVPT2 and CASSCF/NEVPT2 computations, are nearly identical, suggesting consistent trends across different methods and parameter settings.


  Size-consistency errors (in meV) of SA3-iCISCF(14,14), NEVPT2 and SDSPT2 and basis set of cc-pVDZ for the low-lying singlets and triplets energies of anthracene⋯Rg (Rg = He, Ne, Ar, Kr; 100 Å away from the center of anthracene) with respect to anthracene.
  =10pt
  t]@ccccccccccccc@
   Species     iCISCF     iCISCF/NEVPT2     CASSCF/NEVPT2     iCISCF/SDSPT2a     CASSCF/SDSPT2a 

               Singlet        

   anthracene⋯He 1^1A_1     0.0    	 0.1    	 0.0     	 1.9    	 1.8 

   anthracene⋯He 2^1A_1     0.0    	 2.5    	 2.6     	 4.0    	 4.0 

   anthracene⋯He 3^1A_1     0.0    	 2.4    	 1.9     	 4.4    	 3.8 

   anthracene⋯Ne 1^1A_1     0.0    	 0.3    	 0.0      -2.0      -2.3 

   anthracene⋯Ne 2^1A_1     0.0    	 1.9    	 2.6      -1.7      -1.5 

   anthracene⋯Ne 3^1A_1     0.0    	 2.7    	 1.9       1.8       1.0 

   anthracene⋯Ar 1^1A_1     0.0    	 0.2    	 0.0     	30.2    	29.7 

   anthracene⋯Ar 2^1A_1     0.0    	 2.4    	 2.7     	31.4    	30.7 

   anthracene⋯Ar 3^1A_1     0.0    	 2.3    	 2.7     	33.8    	33.5 

   anthracene⋯Kr 1^1A_1     0.0    	 0.3    	 0.0     	95.6    	94.6 

   anthracene⋯Kr 2^1A_1     0.0    	 2.9    	 2.5     	97.2    	94.9 

   anthracene⋯Kr 3^1A_1     0.0    	 2.7    	 1.7     100.3      97.2 

               Triplet         

   anthracene⋯He 1^3A_1     0.2    	-5.3    	 1.3     	-3.1    	 3.0 

   anthracene⋯He 2^3A_1     0.1    	 0.0    	 2.1     	 1.7    	 3.8 

   anthracene⋯He 3^3A_1     0.0    	 6.2    	 2.0     	 7.7    	 3.6 

   anthracene⋯Ne 1^3A_1     0.0    	 1.3    	 1.3      -1.1      -1.4 

   anthracene⋯Ne 2^3A_1     0.0    	 6.2    	 2.1       3.8      -0.3 

   anthracene⋯Ne 3^3A_1     0.0    	 8.0    	 2.1       5.5      -0.4 

   anthracene⋯Ar 1^3A_1     0.2    	-6.4    	 1.2     	24.3    	30.6 

   anthracene⋯Ar 2^3A_1     0.1    	-0.1    	 2.1     	30.3    	31.7 

   anthracene⋯Ar 3^3A_1     0.0    	 8.8    	 2.1     	38.0    	30.8 

   anthracene⋯Kr 1^3A_1     0.3       -12.6       1.3     	84.4    	95.5 

   anthracene⋯Kr 2^3A_1     0.1    	-5.6    	 2.1     	90.7    	96.4 

   anthracene⋯Kr 3^3A_1     0.0    	 7.6    	 2.2     101.5      94.5 

   
  
  
  a SDSPT2 with Pople correction.
  
 





 §.§ The Potential Energy Curve of Cr_2


The chromium dimer, Cr_2, has been a challenging system in quantum chemistry for decades.<cit.> Early on, CAS(12,12) was one of the most commonly used references for MRPT2 or MRCI. However, various theoretical results from MRPT2 and MRCISD suggested that the minimal CAS(12,12), consisting of 3d and 4s orbitals of Cr atom, is insufficient for a quantitative description of the Cr_2 potential energy curve (PEC).
With the development of MR methods, larger active spaces such as CAS(12,22)<cit.> (3d4s4d), CAS(12,28)<cit.> (3d4s4d4p), and even CAS(12,42)<cit.> (3d4s4d4p4f) have been utilized as references to perform subsequent dynamic correlation energy calculations at the MRPT2 and MRCISD levels. Most recently, Chan et al. reported a PEC that is in quantitative agreement with experimental data.<cit.>

In the present work, the CAS(12,22) active space suggested by Chan and coworkers is adopted in our NEVPT2 and SDSPT2 calculations using cc-pwCV5Z-DK basis sets. The NEVPT2 and SDSPT2 results are depicted in Fig. <ref>, alongside the CASPT2 results with a CAS(12,28) reference<cit.> and NEVPT2 results with a CAS(12,22) reference<cit.>. Additionally, the theoretical best estimation (TBE) and experimental curves reported by Chan et al. are presented for comparison<cit.>. Our results show that both NEVPT2 and SDSPT2 methods slightly underestimate the binding energy compared to experimental data, yet offer better accuracy than the DMRG-NEVPT2 results using the same active space<cit.>. Notably, in the bond length region from 2.0 to 2.4 Å, our SDSPT2 and NEVPT2 PECs are more consistent with the latest experimental and TBE results<cit.>, whereas a plateau is observed in the DMRG-NEVPT2 PEC<cit.>.
The absolute energies, and the reference information, are detailed in Table S5. The number of reference CSFs used in the SDSPT2 calculations exhibits fluctuation. However, the total weights of reference determined by P_min=10^-3, account for no less than 98.9% of those by iCISCF(12,22) across all points on the PEC. The size-consistent errors at both NEVPT2 and SDSPT2 levels are less than 0.5 kcal/mol for Cr_2.

The PECs obtained from NEVPT2 and SDSPT2 calculations are sufficiently smooth to enable the fitting of spectroscopic constants. The spectroscopic constants of Cr_2 fitted at the iCISCF(12,22)-NEVPT2 and -SDSPT2 levels are presented in Table <ref>. The binding energies predicted by NEVPT2 and SDSPT2 are 1.486 and 1.449 eV, respectively, which are close to the earlier experimental value of 1.47(5) eV obtained from negative ion photoelectron spectroscopy measurements by Casey and Leopold<cit.>. The bond length and vibrational frequency by our methods are also in good agreement with experimental values<cit.>. The spectroscopic constants produced by relevant MR methods with large active spaces are also listed in Table <ref>. The results from DMRG-SC-NEVPT2 and DMRG-CASPT2 yield similar equilibrium bond lengths to those from SDSPT2, whereas the value from DMRG-ec-MRCISD+Q is slightly longer.







  Spectroscopic constants of Cr_2 by various MR methods.
  =20pt
  t]@ccccccccccccc@
       R_e(Å)     D_e(eV)     ω _e(cm^-1) 

   DMRG(12,28)-CASPT2/cc-pwCV5Zd     1.681     1.610              470 

   DMRG(12,22)-SC-NEVPT2/CBSe     1.655     1.435              469 

   DMRG(12,42)-ec-MRCISD+Q/ANO-RCC-VQZPf     1.71     1.62              479 

   iCISCF(12,22)-NEVPT2/cc-pwCV5Z-DKi     1.658     1.486              469  

   iCISCF(12,22)-SDSPT2/cc-pwCV5Z-DKi     1.659     1.449              462  



   TBEh     1.685     1.58±0.02     495 

   Experiment                 1.679a     1.47(5)b, 1.56±0.06c       481b  


   
  
  
   ^aExp.(1983)<cit.>, ^bExp.(1993)<cit.>,
  ^cExp.(1998)<cit.>, ^dTheo.(2011)<cit.>, ^eTheo.(2016)<cit.>, ^fTheo.(2018)<cit.>,
  ^gTheo.(2020)<cit.>, ^hTheo.(2022)<cit.>, ^iTheoretical computations in this work.
  






 §.§ [Cu_2O_2]^2+ core

The [Cu_2O_2]^2+ complex, a binuclear copper activating O_2 model system, is another popular model for validating the performance of strongly correlated methods. The isomerization energy between the bis(μ-oxo) (F = 0.0) and μ-η^2:η^2 peroxo (F = 1.0) configurations has been extensively studied by various approaches<cit.>. Without sufficiently large active spaces, such as CAS(16,14)<cit.>, CASPT2 might predict an incorrect geometric minimum. Gagliardi and colleagues<cit.> emphasized the importance of including the double-shell effect of copper for producing quantitative correct isomerization energies, which necessitates active spaces comprising 28-32 orbitals.

To compare with results from Chan's DMRG-SC-CTSD<cit.> and Sharma's SC-NEVPT2(s)<cit.>, the same active space CAS(28,32), comprising Cu 3d4d and O 2p3p orbitals, was employed with ANO-RCC-VQZP basis sets. The relative energies of different [Cu_2O_2]^2+ structures, as predicted by iCISCF(28,32) followed by NEVPT2 and SDSPT2 calculations, are presented in Table <ref>. The iCISCF(28,32) results show a shallow well at F = 0.8, indicating that calculating static correlation is insufficient. The NEVPT2 and SDSPT2 calculations yield nearly identical relative energies with respect to the μ-η^2:η^2 peroxo isomer results. The relative energies from NEVPT2 and SDSPT2 closely align with those obtained using SC-NEVPT2(s), the deviations of relative energies do not exceed 4.0 kcal/mol across all examined points. Notably, our SDSPT2 method yields larger relative energies from F = 0.4 to F=0.8 structures compared to those from DMRG-SC-CTSD and SC-NEVPT2(s).
Most recently, Ma and coworkers investigated the isomerization energy between bis(μ-oxo) and μ-η^2:η^2 peroxo isomers as well<cit.>. Utilizing renormalized-residue-based fic-MRCISD+Q with a CAS(24,24) active space, they predict the isomerization energy to be 37.44 kcal/mol. Nonetheless, with sufficiently large active spaces, all these MR methods predict similar isomerization energies between F=0.0 and F=1.0, ranging from 37.4 to 41.3 kcal/mol.


  [Cu_2O_2]^2+ energies (in kcal/mol) from various MR methods, including iCISCF, NEVPT2 and SDSPT2 with CAS(28,32) active space. Energies are relative to the μ-η^2:η^2 peroxo isomer (F = 1.0).
  =20pt
  t]@ccccccccccccc@
   Method      F=0.0     F=0.2     F=0.4     F=0.6     F=0.8     F=1.0 

   iCISCF                        18.5         11.1        5.6         1.1         -1.3       0.0 

NEVPT2                        40.9         35.4       30.5        23.9         13.1       0.0 

SDSPT2                        40.6         35.1       30.2        23.6         12.7       0.0 


SC-NEVPT2(s)a         41.3(8)      33.5(9)     26.3(9)     19.9(8)     9.3(9)     0.0 






DMRG-SC-CTSDb         37.4         29.0        22.0        14.4         6.1       0.0 

  
  
  
  a Ref..
  
  b Ref..
  






 §.§ Cobalt Tropocoronand Complex

To further validate the performance of the newly developed SDSPT2, low-lying states of the cobalt tropocoronand complex, [Co(TC-3,3)(NO)], were studied (as depicted in Fig. <ref>). Both experimental and theoretical evidence has confirmed that the ground state of [Co(TC-3,3)(NO)] is a diamagnetic singlet<cit.>. However, the triplet state assignments remain ambiguous. Ghosh and coworkers, utilizing DFT, characterized T_1 as a low-spin S = 1/2 cobalt(II) ferromagnetically coupled to a NO radical, and T_2 as a high-spin S = 2 cobalt(III) antiferromagnetically coupled to an S=1 NO^- anion. In contrast, recent studies by Freitag et al. employing DMRG-SC-NEVPT2 with a CAS(22,22) reference<cit.> elucidated that T_1 features S = 3/2 cobalt(II) antiferromagnetically coupled to a NO radical, whereas T_2 exhibits S = 1/2 cobalt(II) ferromagnetically coupled to an S=1/2 NO radical. Upon reviewing and repeating the DFT calculations of Ghosh et al., we found that, for the T_2 state at the B3LYP level, the Mulliken spin population on cobalt is 2.636, and the spin population of the CoN_4 center is 2.96. More importantly, as Ghosh et al. reported, the NO bears a spin population of -1.092<cit.>. As a matter of fact, our broken symmetry calculation for the T_1 state started from a quintet high-spin state instead of a septet state. Thus, the T_2 state calculated by B3LYP may be better described as a S = 3/2 Co(TC-3,3) complex antiferromagnetically coupled to a NO radical, which is the T_1 state produced by DMRG-SC-NEVPT2. The results reported by Ghosh et al.<cit.> imply that the broken symmetry B3LYP method might incorrectly predict the ordering of the two triplet states, whereas the excitation energies calculated using pure functionals PW91 and OLYP align with the results of DMRG-SC-NEVPT2.

The electronic structure of [Co(TC-3,3)(NO)] was revisited using the SDSPT2 method in this work, employing the adiabatic molecular geometries<cit.> and the CAS(22,22) active space utilized by DMRG-SC-NEVPT2 computations<cit.>. The def2-TZVP basis set was utilized
alongside the corresponding auxiliary basis for the RI approximation. The initial guess for the iCISCF active space was automatically generated by iCAS<cit.>, including five Co 3d orbitals, a σ ligand MO (denoted as σ_3d_xy), eight π, π^∗ MOs of TC-3,3, four π, π^∗ MOs of NO, and four Co 4d MOs. The converged active MOs for the S_0, T_1, and T_2 states are displayed in Fig. S1-S3, respetively. The occupation numbers for the Co 3d MOs and two π^∗ MOs of NO in the active space are presented in Fig. <ref>. For comparison, the occupation numbers reported by Freitag et al. at the DMRG-SCF level are also provided<cit.>. The occupation numbers predicted by DMRG-SCF and iCISCF for both the S_0 and T_2 states are in good agreement. However, there are slight differences for the occupation numbers of 3d_xy and π^∗_NO,x in the T_1 state, leading to different excitation energies of T_1 reference.
In comparison to the iCISCF outcomes, the T_1 excitation energy is notably overestimated by DMRG-SCF, as depicted in Table <ref>. This discrepancy might stem from issues encountered during the iterative diagonalization process inherent to DMRG. To validate this hypothesis, we reassessed the iCISCF(22,22) with a focus on two roots, optimizing solely the second root. The adiabatic excitation energy for the second root, when active-active rotations were either included or not, was determined to be 38.0 kcal/mol and 38.8 kcal/mol, respectively. These values align closely with the 38.6 kcal/mol obtained from DMRG-SCF computations. Such consistency suggests that the initial root, possessing an excitation energy of 20.7 kcal/mol, was inadvertently omitted in the preceding DMRG-SCF study.

The excitation energies for the T_1 and T_2 states, as determined using iCISCF, NEVPT2, and SDSPT2 methods, are summarized in Table <ref>. Notably, the energy gap S_0 → T_2 computed via SDSPT2 (iCISCF) concurs with the findings from DMRG-SC-NEVPT2 (DMRG-SCF). However, due to disparities in the reference wave functions for T_1, the excitation energy of T_1 derived from SDSPT2 is 19.9 kcal/mol lower than that from DMRG-SC-NEVPT2, which reported degenerate excitation energies of 35.0 and 36.1 kcal/mol for T_1 and T_2, respectively. To reexamine the excitation energies from DMRG-SC-NEVPT2, additional NEVPT2 and SDSPT2 calculations were conducted based on the aforementioned second root of T_1. Both methods yielded an excitation energy of 39.6 kcal/mol, closely matching the 40.1 (40.3) kcal/mol value for T_2, thereby elucidating the degeneracy observed in prior work. For comparative analysis, DFT results by Ghosh et al. are also included<cit.>. It is acknowledged that the S-T gap in metal complexes cannot be accurately predicted by the B3LYP (PW91) functional due to a inherent bias towards triplet (singlet) states<cit.>. Conversely, S-T gaps estimated using the OLYP functional tend to be more reliable<cit.>. The excitation energies for T_1 and T_2 obtained with the OLYP functional align qualitatively with those predicted by SDSPT2 and NEVPT2.






  The excitation energies of two triplet states (in kcal/mol) of [Co(TC-3,3)(NO)] calculated by various methods.
  =10pt
  t]@ccccccccccccc@
   State    iCISCF     SDSPT2     NEVPT2     DMRG-SCFa     DMRG-SC-NEVPT2a     PW91b     OLYPb     B3LYP-D3b 

   T_1     20.7      19.9       19.7        38.6         35.0         24.4     18.0      15.0     

   T_2     29.7      40.1       40.3        29.6         36.1         25.1     23.8      10.4     

   
  
  
  a From ref.. The excitation energies were obtained from DMRG-SC-NEVPT2 (DMRG-SCF) with m=512.
  
  b re-assigned states at the DFT level are from ref..
  





§ CONCLUSIONS AND OUTLOOK


In this study, we have proposed an efficient selection based SDSPT2 method,
thereby extending its applicability to molecular systems characterized by large active space references.
This innovative SDSPT2 approach utilizes molecular orbitals and reference configurations obtained from iCISCF calculations, offering the capability to concurrently and accurately compute both static and dynamic correlations.
A pivotal challenge addressed in this work concerned the construction of the Q using HPS-GUGA.
We accomplished this efficiently by integrating advanced computer storage technologies with distinct processing techniques for different excitation subspaces.
Furthermore, we tackled the issue of dimensionality reduction in the two expansive Q subspaces of D̅V (or S_i^(1))
and V̅D (or S_a^(-1)), as well as minimizing the number of required internal contraction coefficients.
These hurdles were surmounted by amalgamating the FOIS method with sparse matrix storage technique.

Utilizing these strategies, our novel SDSPT2 method has exhibited impressive computational efficiency without compromising on calculational precision, as evidenced by tests on both model and real systems. This method has been effectively applied to molecular systems characterized by large active spaces and a high number of atoms. Notably, the complex Cr_2 system was successfully computed using SDSPT2 in conjunction with iCISCF(12,22) at the cc-pwCV5Z-DK level, yielding accurate potential energy curves and spectroscopic constants. Furthermore, SDSPT2, when paired with the substantial active space CAS(28,32), has been adeptly employed to determine the energy differences across various geometries of the Cu_2O^2+_2 core. Additionally, SDSPT2 with CAS(22,22) has accurately determined the correct energy ordering of three low-lying electronic states for the [Co(TC-3,3)(NO)] transition metal complex, which lacks spatial symmetry. Looking ahead, we aim to integrate the HPS-GUGA method based on configuration selection within other MRCI approaches, such as SDSCI and ic-MRCISD, to further enhance the computational accuracy for complex molecular systems.



§ ACKNOWLEDGEMENT

This work was supported by 
the National Natural Science Foundation of China (Grant Nos. 21833001, 21973054, 22273071 and 22273052) and
Mount Tai Climbing Program of Shandong Province.



§ DATA AVAILABILITY STATEMENT








