The seniority quantum number in Tensor Network States

We employ tensor network methods for the study of the seniority quantum number - defined as the number of unpaired electrons in a many-body wave function - in molecular systems. Seniority-zero methods recently emerged as promising candidates to treat strong static correlations in molecular systems, but are prone to deficiencies related to dynamical correlation and dispersion. We systematically resolve these deficiencies by increasing the allowed seniority number using tensor network methods. In particular, we investigate the number of unpaired electrons needed to correctly describe the binding of the neon and nitrogen dimer and the $D_{6h}$ symmetry of benzene.


Introduction
The quantum mechanical characterization of molecular systems is highly nontrivial due to its manybody character. The need for approximate methods arises for all but the smallest problems. When choosing a suitable approximate method, a consideration has to be made between the accuracy and the complexity of the method. The well-known Hartree-Fock (HF) method provides a mean-field solution for molecular systems rather cheaply. The deficit between the exact ground state energy and the approximate Hartree-Fock energy is an important quantity in quantum chemistry and is called the correlation energy. Often, the distinction between strong (static and nondynamical) and weak (dynamical) correlation is made [1][2][3][4]. Although both strong and weak correlations are electronic by nature, they have a different origin; the latter results mainly from the dynamical shortrange correlations of electrons, whereas static correlation originates from near-degeneracies of several rivaling electron configurations. Many approximate methods often only excel in capturing one type of correlation. For example, the complete active space self consistent field method (CASSCF) [5,6] and density matrix renormalization group (DMRG) [7][8][9] are capable of capturing strong correlations, while coupled cluster (CC) [2,10] and perturbative methods [11] are more suitable for dynamical correlations. In an effort to capture both types of correlation, combinations of these methods have also been developed such as CASPT2 (CAS with perturbation theory up to second order) [12,13], DMRG-CASPT2 [14,15], DMRG-NEVPT2 (DMRG with second-order N-electron valence state perturbation theory) [16], p-DMRG [17], MRCC (multireference coupled cluster) [18][19][20] and DMRG-TCC (DMRG-tailored coupled cluster) [21][22][23].
The majority of contemporary electronic structure methods start from a reference state, typically the single-reference HF ground state, and systematically build in correlations by considering elementary excitations from this reference. The conventional approach is to consider particle-hole (ph) excitations from the HF ground state, as is common in CC [2] or truncated configuration interaction (CI) methods [1]. This way, it is possible to construct a hierarchy of multiple n-ph excitations which are assumed to be decreasing in importance with increasing n. Although tailormade for dynamical correlations, e.g. in CC theory, it is impractical for static (or non-dynamical) correlation. It was recently observed [24] that the seniority scheme is much better suited to capture static correlations associated with the entanglement structure of single-bond breaking processes. Defined as the number of unpaired electrons in a Slater determinant, the seniority quantum number organizes the Hilbert space by the amount of broken closed-shell singlet pairs with respect to a set of (doubly degenerate) spin orbitals. For molecular systems dominated by singlet-pairs bond structures, it was shown that most of the strong static correlation in a system can already be captured in the subspace spanned by all determinants with zero seniority (no unpaired electrons) [24][25][26][27][28][29][30]. Although this tremendously reduces the dimension of the Hilbert space at hand, finding the exact doubly occupied configuration interaction (DOCI) wave function is still an exponentially scaling problem. At first glance, the seniority scheme seems only marginally more manageable than the full problem. Interestingly, the antisymmetric product of one-reference orbital geminals (AP1roG) [28,[31][32][33], also known as pair-coupled cluster doubles (pCCD) [34][35][36], appears to provide a reliable approximation to the DOCI ground state solution for a wide range of molecular systems [36] while staying computationally tractable at a mean-field scaling computational cost [28,37].
Notwithstanding its salient features, there remain several challenges that need to be overcome in order to make the AP1roG wave function quantitatively accurate. The outstanding challenges, which are shared by all methods expressed in the seniority scheme, are (i) the incorporation of dynamical correlation and (ii) the choice of a preferential orbital set, also referred to as the orbital optimization (OO) problem. Another challenge (iii) is the apparent lack of London dispersion correlations in the seniority-zero methods which are crucial to model large-size molecular systems.
The lack of dynamical correlation in the zero seniority wave functions is well illustrated by the poor description of the correlation energy of the Ne atom, as well as the near-constant parallelity error in the bond dissociation curve of the nitrogen dimer [28,38,39]. Dynamical correlation is generally encoded in a set of Slater Determinants with a few ph excitations from the HF reference state. Consequently, methods targeting these excitations, such as (multi-reference) perturbation theory [38,40], linearized CC [35,41], extended random phase approximiation (ERPA) [42] or selected configuration interaction (CI) [43], are very well equipped to capture those correlations. However, systematic generalizations of these methods in order to include dynamical correlation from higher-order ph excitations prove either technically and computationally demanding, or break the size consistency of the reference AP1roG wave function.
Size consistency of the AP1roG wave function, or by extension any seniority-zero method, is guaranteed when the spin orbitals are optimized such that the energy is in a variational minimum [44]. Unfortunately, current optimization methods result into a single unique set of spin orbitals, which can lead to nonphysical symmetry breaking effects in resonating bond structures, such as the aromatic structures in benzene [33], or incorrect characterizations of covalent tripletbond couplings, such as in the nitrogen dimer [44].
Regardless the size consistency and correct description of the static correlations associated with bond-dissociation processes, seniority-zero methods have recently been identified as essentially free from London dispersion energy [45], which is remarkable given that 2-electron systems are exactly described by (orbital optimized) seniority-zero methods, capturing the non-covalent Lennard-Jones 1/R 6 behavior of the dispersion energy in the large R → ∞ separation limit of the hydrogen dimer.
In order to obtain a global understanding of the deficiencies of the seniority-zero methods, it is quintessential to include all possible broken-pair excitations from higher seniority sectors in a systematic way. Higher seniority subspaces have been studied in the past years using CI approaches, [24,29,34,46,47], or energy renormalization group (ERG) approaches [48]. The limiting factor of these methods is the pernicious computational scaling whenever no truncation in the Slater determinants is considered. While dynamical correlation is typically included with just a few ph excitations from the HF reference state, corresponding to low-seniority quantum states, it is not clear at present how many broken pairs are needed to restore the correct symmetries or include London dispersion. As a result, there is a need for an analytic method that can assess seniority non-zero contributions in a systematic way at a favorable computational scaling.
In this paper, we use the concept of seniority in junction with tensor network states. In contrast to many other quantum many-body methods, tensor network states consider the whole collection of Slater determinants, and approximate the exact quantum states by restricting the amount of entanglement between local degrees of freedom. Tensor network states are capable of encoding local symmetries of quantum states [49]; therefore they provide a good framework to investigate broken pair excitations, as seniority can be related to the irrep label of the su (2) quasi-spin algebra [50]. In practice, the idea is to perform DMRG in a subspace of the Hilbert space up to a fixed global seniority quantum number, and increase the seniority quantum number until full convergence of the correlation energy is obtained. This procedure will be explained in detail in Section 2. In the proceeding sections, we will present results for the nitrogen dimer (Section 3.1), benzene (Section 3.2) and the neon dimer (Section 3.3), to discuss higher-seniority properties of dynamical correlation, symmetry breaking/restoration and dispersion respectively.

Tensor networks
Pioneered by Steve White in 1992 [7,8], tensor networks have proven to be a natural language for the entanglement in strongly correlated many body systems. In the tensor network state, each tensor represents a 'local' physical degree of freedom. By connecting them in a network, correlations between the different physical degrees of freedom can be encoded through their virtual degrees of freedom. The exact layout of the network influences the entanglement structure that can be represented; it is easier to correlate physical degrees of freedom that are close in the network.
First and foremost, these tensor network methods have established themselves in the field of condensed matter physics as a wide range of successful tensor networks have been developed for numerous problems. Some notable examples are matrix product states (MPS) [7,8,[51][52][53], projected entangled pair states (PEPS) [54] and the multiscale entanglement renormalization ansatz (MERA) [55]. They all provide, in their own way, an efficient representation of certain entanglement structures.
In quantum chemistry, tensor networks have also proven their worth in the study of molecules with strong correlations [9,[56][57][58][59][60][61][62][63][64][65]. Quantum chemists don't traditionally study molecules in a Hilbert space built from completely local basis functions (e.g. a grid in three dimensions), but atomic orbital basis sets such as Gaussian-type or Slater-type orbitals are used. These sets give electrons the right flexibility needed for chemistry while the basis size is kept small. On the flip side, the loss of locality in the basis functions makes a suitable network for the entanglement between the physical degrees of freedom less straightforward than for most condensed matter problems. Furthermore, in an atomic orbital basis set, the long range two-body coulomb interactions in the Hamiltonian become four-point interactions. The loss of locality and the need for an efficient evaluation of the Hamiltonian has ensured that the most simple networks are still the most preferred ones. The density matrix renormalization group (DMRG) is, by far, the most popular tensor network method in quantum chemistry and corresponds with the optimization of the linear MPS. Another option for a simple tensor network is the three-legged tree tensor network state (T3NS) [66,67]. It is a subclass of the more general tree tensor networks (TTNS) [63][64][65] and was recently introduced by some of us. In this paper, we use these two networks for the study of several chemical systems in restricted seniority subspaces. In the next sections we explain the implementation of restricted seniority for the case of DMRG. However, the ideas are readily adaptable to T3NS and were implemented for both cases in our in-house T3NS-code [68].

Seniority and tensor networks
The non relativistic quantum chemical Hamiltonian to study is given by where i, j, k and l are the indices of the orbitals and σ and τ index the spin of the electrons. This Hamiltonian showcases several symmetries, e.g. the particle conservation and total spin symmetry of the electrons. These symmetries can be easily exploited in tensor networks by writing the different tensors in an invariant form under group action of the symmetry [60,67,[69][70][71][72][73][74][75][76]. Although the seniority is not a symmetry of the quantum chemical Hamiltonian, it is still possible to apply the same idea. In this case, we write each tensor in the network in an invariant form for the seniority. For example, the tensors of rank three present in the MPS can be made invariant by imposing the following restriction for the tensor elements: or graphically where a, b and c denote the (physical or virtual) degrees of freedom of T and ν a , ν b and ν c are their respective seniority numbers which is well-defined; each state in the degrees of freedom a, b, and c are eigenstates of the seniority operator. In this example, the seniority number of the first two degrees of freedom of the tensor sum up to the seniority number of the third one. This reflects the fact that seniority is an additive feature, as unpaired electrons from different orbitals all contribute to the total seniority of the state. It is clear that this restriction implies a kind of flow for the seniority number in the network which is indicated in Eq. (3) by the directed edges. An example of an MPS wave function built from three of these invariant tensors with the flow indicated is given by The physical degrees of freedom (the occupancies of the spatial orbitals) are denoted by a, b and c in this example and have a seniority ν ∈ {0, 1}. α and β are virtual degrees of freedom. The vacuum state enters the MPS at the leftmost degree of freedom (vac) and has a seniority ν = 0, while the last degree of freedom represents the targeted state |Ψ〉. The graphical depiction implies for each connected edge a summation over its corresponding indices, as is common in the tensor network language [49,54,67,[77][78][79]. The only difference with implementing a U(1)-symmetry of the system, e.g. particle conservation or conservation of the spin projection, is the needed summation over the renormalized states of the target edge in Eq. (4). This is necessary as the seniority is not a conserved quantum number. Eigenstates of the Hamiltonian are not necessarily eigenstates of the seniority operator and the target state can be a linear combination of Slater determinants with different seniority numbers. To target such a state, the final states at the edge called target should be a set of eigenstates of the seniority operator which combine to the targeted state when summed.
The set of possible seniority numbers for the wave function is with k the number of spatial orbitals, N ↑ (N ↓ ) the number of electrons with spin up (down) and N tot the total number of electrons. For every renormalized state in the last bond, we have ν target ∈ Ω. By restricting ν target to a subset S, i.e. ν target ∈ S ⊆ Ω, ground states in seniority-restricted subspaces can be targeted.
In a similar fashion, one could use also use other non-conserved quantum numbers than the seniority. For example we could use the excitation number with respect to the Hartree Fock wave function. By only allowing Slater determinants with a certain amount of excitations, tensor networks can be used as an approximate configuration interaction (CI) solver with arbitrary allowed excitation levels.

Suboptimal decomposition
When using a wave function ansatz as shown in Eq. (4), we impose a restriction on the left renormalized states at each splitting of the network. Due to the fact that a vacuum state enters in the left most bond (vac) and all tensors used are invariant under the seniority operator, the left renormalized states need to have a well defined seniority number. This restriction does not hold for the right renormalized states since multiple states with different seniority exit at the right most bond (target).
This restriction results in the need of a possibly larger bond dimension than when discarding seniority. We illustrate this using a wave function with three electrons in three spatial orbitals: In Eq. (7), the Schmidt decomposition for a partitioning between the first two and the last orbital is given. At this partitioning only a virtual bond dimension of one is needed to represent the state. However, when we impose that the left states, i.e. the states in the first two orbitals, of the decomposition should also be eigenstates of the seniority operator, the needed bond dimension at this partitioning increases to two, confer Eq. (6).

DOCI and tensor networks
Restricting the calculation to configurations with ν = 0, i.e. all electrons are paired, is easily done with the aforementioned method. However, it is more efficient to directly implement the quantum chemical Hamiltonian projected on the DOCI-subspace where only paired electrons are allowed. The DOCI-Hamiltonian is given by where b † i and b i are the bosonic pair creation and annihilation operators and n i is the pair number operator at orbital i. They are given by and This Hamiltonian only scales quadratically with the number of orbitals in contrast with the quartic scaling of the full Hamiltonian. TNS calculations in the DOCI subspace can be performed with a lower polynomial scaling, as stated in Table 1 for the DMRG and T3NS. We find that DOCI ground state wave functions have in general lower entanglement than their corresponding FCI ground state wave function; accurate results for DOCI can be obtained with a much lower bond dimension. The synergy between the lower polynomial scaling and the lower bond dimension needed, makes DOCI-TNS very fast and a good option for initializing tensor network calculations in the FCI space. For example, DOCI calculations with 162 electron pairs and 261 spatial orbitals have been executed in a few minutes on a common laptop.

Applications
We discuss some calculations with the seniority-restricted tensor network code. As these calculations are orbital dependent, several types of orbitals are considered. The effect of allowing progressively more broken pairs is also studied within each orbital set. In Section 3.1 and Section 3.3, the dissociation of the nitrogen and neon dimer are considered, respectively. Section 3.2 discusses the benzene molecule, a system demonstrating artificial D 6h symmetry breaking in the seniority-zero subspace [33]. Coupled cluster natural orbitals and Löwdin orthogonalized atomic orbitals are obtained with PySCF [80][81][82]. DOCI-optimized orbitals are generated through an in-house DOSCF code and were carefully checked to correspond to the lowest possible DOCI energy, i.e. the global minimum [44]. The seniority-restricted tensor network calculations were executed with our own T3NS-code [68]. All seniority-restricted tensor network calculations are MPS calculations. We exploit the spin symmetry and the reported bond dimensions for the tensor networks are reduced bond dimensions; renormalized states belonging to the same multiplet are represented by one reduced renormalized state, thus reducing the needed bond dimension. Seniority-restricted tensor network calculations are, just as regular tensor network calculations, not exact; the accuracy can be controlled by the bond dimension. The following calculations use a large enough bond dimension to ensure quantitatively accurate potential energy surfaces.

Nitrogen dimer
Characterized by a triple bond breaking, the nitrogen dimer is a much visited test case for new quantum chemical methods, and has already been investigated as such in the seniority framework by Bytautas et al. [24] using an active space in the cc-pVDZ basis with D 2h -symmetry adapted MOs. Here, we study the nitrogen dimer in a cc-pVDZ basis set with all electrons correlated, however the DMRG results are qualitatively similar to the results in [24]. Seniority-restricted spin-adapted DMRG with a reduced bond dimension up to a 1000 is used to optimize the ground state in the different subspaces. The allowed seniority increases from 0 (DOCI) up to 10 for the largest calculations, allowing 5 electron pairs to be broken. In Fig. 1, the dissociation curves are given for calculations within the different seniority subspaces. Calculations were performed for canonical orbitals (Fig. 1a), DOCI-optimized orbitals (Fig. 1b) and CCSD natural orbitals (Fig. 1c). Although the DOCI-optimized orbitals are optimized for the seniority-zero subspace specifically they also perform better in higher seniority subspaces, albeit marginally. Eventually for ν ≤ 8 and onward, all orbital sets give quasi-FCI energies. (c) In Ref. [32], it is shown that the seniority-two sector decouples from the seniority-two-plus-zero sector up to first order for DOCI-optimized orbitals; only a small correction should occur due to the introduction of single broken pairs in this orbital set. Putting this first order decoupling to the test, we notice indeed a small energy correction for the DOCI-optimized orbitals, smaller than for canonical orbitals. In Fig. 2, the weights of the different seniority subspaces are plotted for the ground state in both canonical and DOCI-optimized orbitals. It is yet another illustration that for DOCI-optimized orbitals (Fig. 2b) the seniority-two subspace is less important than for canonical orbitals (Fig. 2a). However, a first order decoupling is not an exact one; there are other orbital sets possible which give even smaller energy corrections. This is illustrated by the natural orbitals (Fig. 1c) which give even smaller energy corrections when allowing single broken pairs in this system. Löwdin orthogonalized atomic orbitals (c) are given.
As a last observation we note that the largest change in energy occurs when including the seniority-four subspace, and this for all orbital sets in Fig. 1. This trend was also noticed in Ref. [24] for the nitrogen dimer in nonlocal orbitals. When including up to seniority four the energies are close to FCI around the binding distance; however, the binding energy itself is still overestimated due to missing dynamical correlation at the dissociation (values are given in Table 2 for both canonical and DOCI-optimized orbitals).
Intuitively, we would expect a much larger error when excluding the seniority-six subspace as Hund's rule dictates dissociation to two nitrogen atoms with each three unpaired electrons. However, seniority and pairing is an orbital-dependent concept [44]; we need to keep in mind 10  canonical 424  382  377  338  322  322  DOCI-optimized 383  405  367  333  321  321   Table 2: Binding energies in mE h for seniority-restricted calculations in both canonical and DOCI-optimized orbitals.
that Hund's rule applies to a nitrogen atom with orbitals localized around that atom. To study the interpretation of Hund's rule in non-local orbitals, we consider a toy model of two sets of three orbitals (p x , p y , p z ) and (p x , p y , p z ). Each set of orbitals mimics the local p-orbitals of each nitrogen atom which are singly occupied and couple together to a S = 3 /2 state, as dictated by Hund's rule. Our tensor network calculations target over the whole dissociation curve a singlet state for the dimer, so the two toy-nitrogen atoms should couple as [ 3 /2, 3 /2] 0 . Mimicking non-local orbitals, we rotate the orbitals pairwise as follows: π 1 = p x cos θ + p x sin θ , π * 1 = −p x sin θ + p x cos θ π 2 = p y cos θ + p y sin θ , π * 2 = −p y sin θ + p y cos θ σ = p z cos θ + p z sin θ , σ * = −p z sin θ + p z cos θ . In Fig. 3, the weights of the different seniority sectors is given for the [ 3 /2, 3 /2] 0 coupled toy wave function in function of the rotation angle θ . As can be seen in this model seniority-six is actually of no importance when working with delocalized orbitals (θ = π /4). Instead, the correct dissociation can be described with only seniority-zero-plus-four and both seniorities equally important.
Both canonical orbitals and the DOCI-optimized orbitals are delocalized for the 2p-orbitals in this system. This dominating importance of the seniority zero and four for the wave function at dissociation is very clear in Fig. 2a and Fig. 2b. The other seniority sectors have very small contributions in comparison. As an illustration, we also included calculations with Löwdin orthogonalized atomic orbitals in Fig. 2c. As these orbitals are localized, it corresponds with θ = 0 in Fig. 3. These orbitals do give rise to a very important seniority-six subspace at dissociation, as predicted by Hund's rule. Evenmore, all seniority sectors smaller than six express a superexponential decay. Figure 4: Benzene in a STO-6G basis set for different distortion angles. The minimal energies and the corresponding distortion angles for increasing seniority subspaces are given in the inset. A graphical depiction of the in-plane distortion of the aromatic ring in benzene is also shown.

Benzene
In this section, the in-plane distortion of benzene is studied. The exact nature of the distortion is given in the inset in Fig. 4 and is characterized by the angle θ . At θ = 60 • the D 6h symmetric equilibrium structure geometry of benzene is obtained. At other angles, the distortion introduces alternating shorter and longer carbon-carbon bonds. For this system Boguslawski et al. [33] showed that benzene (θ = 60 • ) is not the equilibrium structure within the seniority-zero subspace; an artificial symmetry breaking occurs when allowing orbital optimization.
We use DOCI-optimized orbitals in the STO-6G basis set to study this artificial symmetry breaking with all electrons correlated. We chose STO-6G as the distortion angle of the minimal energy DOCI structure is particularly large for this basis set. The tensor network calculations are executed with a reduced bond dimension of 1000.
In Fig. 4 the results for the ground state in the different seniority subspaces are given. In accordance with Boguslawski et al. [33], we notice that, indeed, the ground state is not found at θ = 60 • in the seniority zero subspace. When the breaking of one pair is allowed in this orbital set, the correction is rather small and the correct symmetry is not restored; as expected due to the aforementioned first order decoupling of the seniority-zero and seniority-two subspace in these orbitals. 57 When including progressively higher seniorities, the stable configuration moves closer to the expected D 6h symmetric benzene. The potential energy surface enjoys a large qualitative correction when including the seniority-four subspace in the calculations. However the predicted most stable configuration is still off by 0.59 • . The inclusion of seniority-six further improves the quality of the potential energy surface, but only at seniority-eight the correct symmetry seems to be recovered, at least up to the resolution of our performed calculations. At this point, the results become very close to the full seniority results. In Fig. 5, the weights of the different seniority sectors in the ground state during distortion are also given. These weights do not express the large changes as were seen during dissociation of the nitrogen dimer in Fig. 2. This is quite expected as the bond breaking of the nitrogen dimer is a far more outspoken change than the small benzene distortions in this section.

Neon dimer
The neon dimer, constituted by just two noble gas atoms, is very weakly bonding. Although the electrons do not form covalent bonds between the two atoms, it expresses some bonding character due to weak dispersion forces. In Ref. [83], an empirically fitted potential curve results in a binding energy of −134 µE h and a binding distance of 3.091 Å.
As the binding of the neon dimer is rather weak and due to dynamical correlations, it will be very sensitive to the chosen basis set size. For an accurate description of the potential energy curve, a large basis set should be chosen and basis set superposition errors (BSSE) should be taken into account appropriately [84]. A clear example of the importance of BSSE-corrections is the dissociation curve on the Hartree-Fock level. At this level of theory no binding is expected as the Hartree-Fock solution is dispersion-free. However, when using small basis sets, one would find a binding neon dimer at the Hartree-Fock level if one neglects to correct the BSSE [84].
We study the neon dimer in the aug-cc-pVDZ basis; it was found that this basis set has a favorable tradeoff between mitigating BSSE and numerical stability issues of larger basis sets. Calculations with different seniority sectors are executed while using DOCI-optimized orbitals with a frozen 1s core. Reduced bond dimensions up to 800 are used for the DMRG calculations. As the aug-cc-pVDZ basis is a rather small basis for capturing dispersion forces, appropriately removing BSSE is important. This is done by using the Boys and Bernardi counterpoise correction [85].  In Fig. 6 the raw uncorrected results are given for the different calculations. For all seniority calculations the neon dimer seems to be bonding. However, for seniority-zero and senioritytwo-plus-zero the bonding seems to be very weak; only for ν ≤ 4 calculations and higher the bonding character is qualitatively corresponding with the full seniority case. For the counterpoise correction, equivalent calculations as for the dimer are executed but where one neon atom is replaced by a chargeless, electronless ghost atom. This way, we can approximately correct for the extra stabilization each neon monomer experiences by the extra added basis functions of the other monomer. The BSSE-corrected dissociation energy for the dimer at distance r is then given by The same level of theory should be used for these ghost calculations as for the original calculation. This poses a difficulty since the seniority restricted calculations are not size consistent; E dissoc (r → ∞) does not tend to zero as is desired. Assume we have executed a dimer calculation with ν ≤ 4, using ghost calculations with the same ν ≤ 4 will over-correct, while ghost calculations at the lower ν ≤ 2 will under-correct. We try to solve this problem by both over-and undercorrecting and shift both curves to 0 in the dissociation limit. Results for these BSSE-corrections are given for different seniority sectors in Fig. 7 with the grayed area indicating where the exact BSSE-correction is expected to be. From Fig. 7a and Fig. 7b, it seems that the weak bonding character present in Fig. 6 for ν = 0 and ν ≤ 2 practically or completely disappears when taking BSSE-corrections into account. When correcting ν ≤ 4 calculations, the over-corrected dissociation underestimates the dissociation energy a bit with respect to the FCI BSSE-corrected calculations while the under-corrected dissociation overestimates the dissociation energy, as can be expected. It seems thus that calculations with seniority-zero and seniority-two do not model the needed dispersion and at least seniority-four is needed. This suggests the breaking of at least one electron pair at each Ne atom is needed, inducing polarization effects in each atom which give rise to the dispersion energy.  : Dissociation curves for the neon dimer with and without BSSE correction for different seniority subspaces. When BSSE-correction is performed, the used seniority subspace for the ghost atom calculation is given in brackets in the legend. The dissociation curve for a full-seniority calculation with full BSSE-correction is also shown in each subfigure. Results are shown for ν = 0 (a), ν ≤ 2 (b) and ν ≤ 4 (c) subspace calculations.
Finally, we notice that the FCI dissociation energy with BSSE-correction is a factor of three smaller than empirical measurements and the bond length is overestimated. This is quite normal when studying dynamical correlations in small basis sets. The limited basis set does not allow all the needed flexibility for the stabilization of the dimer.

Conclusion
In this paper, the concept of seniority is joined with tensor network methods. By using seniorityinvariant tensors in a tensor network, we can force all the renormalized states in the virtual bonds to be eigenstates of the seniority operator. This allows for arbitrary seniority-restricted calculations. For DOCI (doubly occupied configuration interaction) calculations, we can immediately implement the DOCI-projected quantum chemical Hamiltonian in Eq. (8). This results in a very fast tensor network calculation, partly because of the simpler Hamiltonian, partly because the correlations in the seniority-zero subspace for molecular systems are easily captured by tensor networks; even for very large systems, a bond dimension of less than 100 suffices for energies within chemical accuracy of the exact DOCI energy.
Several systems are studied within different seniority subspaces. For the dissociation of the nitrogen dimer, only a quantitative dissociation curve can be obtained when at least two pairs are allowed to be broken. This can be theoretically explained due to Hund's rule. The in-plane distortion of benzene and its artificial D 6h symmetry breaking in the seniority-zero subspace [33] is also studied. A large correction of the artificial symmetry breaking occurs when including seniority-four, however up to eight unpaired electrons are needed for a complete restoration of the correct benzene point group symmetry in the used basis set. Finally, also the dissociation of the neon dimer is considered. At the seniority-zero level of theory the neon dimer is non-binding; DOCI does not capture the dispersion forces needed for the weak binding characteristic of neon. Only at seniority-four and onward, the dispersion forces are adequately picked up.
For all systems, the seniority-two subspace has only a small contribution to the total wave function when using DOCI-optimized orbitals; as expected by the theoretical first order decoupling between seniority-zero and seniority-two subspaces in these types of orbitals [32]. However, a first order decoupling is not an exact decoupling and other orbital sets can be found which attribute even less importance to the seniority-two subspace. An example of this is given by the natural orbitals of the nitrogen dimer in Fig. 1c.