key: cord-0775610-kgcbkxtc authors: Herrmann, Nils; Hanrath, Michael title: A correctly scaling rigorously spin-adapted and spin-complete open-shell CCSD implementation for arbitrary high-spin states date: 2021-11-07 journal: J Chem Phys DOI: 10.1063/5.0078020 sha: ae2abd4f2c04adcc3d9d91c6514177848b1ddcbd doc_id: 775610 cord_uid: kgcbkxtc In this paper, we report on a correctly scaling novel coupled cluster singles and doubles (CCSD) implementation for arbitrary high-spin open-shell states. The chosen cluster operator is completely spin-free, i.e. employs spatial substitutions only. It is composed of our recently developed L"owdin-type operators (doi: 10.1063/5.0026762), which ensure (1) spin completeness and (2) spin adaption i.e. spin purity of the CC wave function. In contrast to the proof-of-concept matrix-representation-based implementation presented there, the present implementation relies on second quantization and factorized tensor contractions. The generated singles and doubles operators are embedded in an equation generation engine. In the latter, Wick's theorem is used to derive prefactors arising from spin integration directly from the spin-free full contraction patterns. The obtained Wick terms composed of products of Kronecker deltas are represented by special non-antisymmetrized Goldstone diagrams. Identical (redundant) diagrams are identified by solving the underlying graph isomorphism problem. All non-redundant graphs are then automatically translated to locally -- one term at a time -- factorized tensor contractions. Finally, the spin-adapted and spin-complete (SASC) CCS and CCSD variants are applied to a set of small molecular test systems. Both correlation energies and amplitude norms hint towards a reasonable convergence of the SASC-CCSD method for a BCH series truncation of order four. In comparison to spin orbital CCSD, SASC-CCSD leads to slightly improved correlation energies with differences of up to 0.886 mE_H (0.98% w.r.t. spin orbital CCSD(T)). The CC method, pioneered by Coester and Kümmel [1, 2] , emerged as one of the most successful ab-initio methods of modern quantum chemistry. Within its framework, the cluster operatorT is applied to a reference |Ψ 0 to generate the CC wave function |Ψ CC via |Ψ CC = eT |Ψ 0 . In the most commonly used similarity-transformed CC approach, |Ψ 0 is confined to a single Hartree-Fock determinant whileT is expressed by spin orbital substitutions moving particles from the occupied to the virtual orbital space. Despite tremendous progress in the last decades , the formulation of a general coupled cluster (CC) approach to treat (high-spin) open-shell cases is still an ongoing research field. In this work, we focus specifically on single restricted open-shell Hartree-Fock (ROHF) type references resembling S = S z high-spin states. Through the application of a usual spin orbital cluster operatorT (incorporating e.g. spin orbital substitutions) to open-shell references |Ψ 0 (see e.g. [3] ), spin contamination is introduced into the CC wave function even if |Ψ 0 is represented by a genuine spin eigenfunction (see e.g. [40] [41] [42] ). This contamination arises due toT being non-commutative with theŜ 2 operator, in general. * N.Herrmann@uni-koeln.de † Michael.Hanrath@uni-koeln.de Several approaches to adapt the cluster operatorT (as originally motivated by the symmetry-adapted cluster operators of Nakatsuji and Hirao [4, 5] ) such that it commutes withŜ 2 -effectively removing all spin contamination of the CC wave function -were proposed in the literature. Naturally, the special case of closed-shell spin adaption -being fundamentally simpler than the general open-shell case -received much attention especially in the early years of spin-adapted CC theory (see e.g. [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] ). The present work, however, is concerned with the more general open-shell cases. In the literature, projective methods [6, 7] were developed, where spin-contaminated parts of the spin orbital CC expansion are removed applying an appropriate projector. Analogously, restrictive methods [8] [9] [10] [11] [12] [13] [14] [15] were proposed explicitly restricting the CC amplitudes such that spin contamination is removed. In the latter approaches, Neogrády et al. developed spin-adapted linear CCSD [8] , full CCSD [9] as well as CCSD(T) [10] for the doublet spin state. Szalay and Gauss [11] developed a restricted CCSD formalism, which was later expanded to treat excited states [12] , to the inclusion of full triples [13] as well as to an explicitly correlated R12 approach [15] . In this formalism, the correctŜ 2 expectation value of the CC wave function is enforced by the inclusion of explicit spin equations. To make their approach feasible, the proposed spin equations are followed in a truncated manifold only. The resulting CC wave functions are therefore not spin-adapted in general. Later, a comparison of spin-adapted and spin-restricted CCSDs was given by Heckert et al. [14] . Furthermore, a partially spin-adapted CCSD was presented by Knowles and Werner [16, 17] , where the cluster operator is com-posed of a mix of spin-free (spatial) and spin orbital substitutions such that the linear CCSD terms are spinadapted while the non-linear terms, hence the CC wave function eT |Ψ 0 , is not. A different approach to achieve spin adaption is followed by the so-called unitary group approach (UGA) [18] [19] [20] [21] [22] [23] [24] [25] [26] . Here, unitary group theory is applied to generate spatial, i.e. spin-free, orbital substitutions such that orthonormal configuration state functions (CSFs) are built. A cluster operator built by these substitutions is guaranteed to produce a properly spin-adapted CC wave function. Through careful choice of the selected substitutions, the completeness of the spanned spin space may also be reached. Linear CCSD [18] and full CCSD [19] theories with special emphasis on doublets, triplets and openshell singlets were developed and analyzed w.r.t. ROHF doublet [23] and open-shell singlet and triplet [24] instabilities. Further applications include the 1 A 1 , 3 B 1 separation of the methylene molecule [20, 21] , calculations of the ozone molecule [22] and an analysis of van der Waals interactions of high-spin systems [26] . As discussed by Matthews et al. [55, 56] for the closedshell case, theory and implementation of orthogonal UGA CC are quite complex compared to non-orthogonal approaches. In the latter, the cluster operator may be composed of solitary spatial substitution operators instead of linear combinations of such. This not only simplifies equation derivation, but also enables efficient implementations. Already in 1991, Janssen and Schaefer III presented an automated equation derivation engine [27] , which was also applied to a non-orthogonal spin-adapted CCSD for high-spin open shell states. Although their approach allows for spectating substitutions, terms required for spin completeness in triplet and higher spin states are missing from their cluster operators [38] . Following the early suggestion of Lindgren [57] to approximate the wave operator eT by its normal-ordered form {eT }, several CCSD [28, 29] and MRCCSD [30] implementations emerged. While the assumed normal ordering of the wave operator removes a lot of complexity from the theory (mainly by avoiding internalT contractions), its effects on the quality of the CC wave function are unknown. One of the latest developments in spin-adapted CC theory is the combinatoric open-shell CC (COSCC) approach [31] [32] [33] [34] [35] [36] [37] . Here, the wave operator is assumed to be normal-ordered while contractions among the cluster operators are reintroduced following a combinatoric cluster expansion. Theories were developed for singlereference CC [31] as well as state-universal [32] and stateselective [33] MRCC. An automated implementation engine was developed [34] and used to generate COSCCSD equations for doublet spin states. The approach was further expanded to the calculation of analytic first derivatives [35] , spin densities [36] and hyperfine coupling tensors [37] . Overall, the published spin-adapted open-shell CC methods seem to be limited by the doubles truncation and the triplet spin state. In this work, we present a newly developed CCSD implementation capable of treating arbitrary high-spin states. The presented approach utilizes our previously developed Löwdin-type non-orthogonal spatial substitution operators [38] to guarantee (i) spin adaption and (ii) spin completeness of the CC wave function. While (i) is automatically fulfilled by the usage of spin-free substitution operators, (ii) requires special care as to the selection of (usually linear-dependent) substitution operators. In our previous work [38] , this is ensured by Löwdin's projection operator method. The spin-adapted and spin-complete (SASC) CCSD implementation, presented here, utilizes a single set of CC equations. These are applicable to arbitrary S = S z states without any further equation generation or compiling. In principle, the presented methodology is directly applicable to higher CC truncation levels as well. For the sake of this publication, however, we limit ourselves to the spatial SD manifold (involving spin orbital substitutions of up to rank four). In section II, a short introduction to the SASC-CC theory is given. Section III then focuses on the automated equation generation and implementation of the SASC-CCSD method, while section IV presents proof-of-concept applications on small test molecules of spin states up to S = S z = 4. In this section, theoretical intricacies for the presented CCSD are given. Subsection II A features a short recap of the non-orthogonal, spin-adapted and spin-complete CC framework -introduced in our previous work [38] -while subsection II B introduces the spin-free CCSD equations considered in this work. A. Spin adaption and spin completeness in the CC framework As already introduced in [ where p 1 . . . p ν ∈ O ∪ A and q 1 . . . q ν ∈ A ∪ V. TheX operators denote spin orbital substitution operators given by their usual second quantized form aŝ [38] . Special care needs to be taken to construct cluster operators, which (i) are spin-adapted, (ii) produce complete spin spaces for arbitrary spin states and (iii) do not involve linear-dependent substitutions. Such operators were generated in [38] by the use of Löwdin's projection operator method. These Löwdin-type operators for the single and double substitutions are given in equations (3) and (4). In this work, we aim at solving the similarity transformed CCSD equations with forT =T 1 +T 2 as defined in equations (3) and (4) (7) Here, the indices p, q, r and s denote arbitrary spatial orbitals from orbital spaces O, A or V. In order to apply Wick's theorem toĤ, the spin-freeÊ operators may be expanded to their spin orbital analogue. Within the usual particle-hole formalism for the high-spin whereb † andb denote quasi particle/hole creation and annihilation operators, respectively. Clearly, the usual contraction rules (including the spin restriction δ σω ) hold in the occupied and virtual orbital spaces O and V. For the singly occupied active space A, the surviving contractions depend on the spin information of the contracted orbitals. For α orbitals (occupied in the spin orbital picture of the high-spin reference), only creator -annihilator contractions survive, while for β orbitals, the opposite is true. This is of particular importance when evaluating spin summations. While contractions in O or V leave at least one spin summation intact, contractions in A break down to either the α or the β part of the summation. The operatorsÊ p q andÊ pq rs may now be normal-ordered with respect to the given high-spin reference determinant |Ψ 0 (indicated by curly brackets) to yield and Interestingly, the single contractions of σ 1 and σ 2 spins occurring inÊ pq rs lead to a single spin summation for occupied indices (colored in red) and to a completely broken spin summation for active indices (colored in blue). In the latter case, this leads to pure spin orbital substitu-tionsX of α electrons. Therefore,Ĥ may be decomposed into its normal-ordered partĤ N such that two separate one-particle interrelating operatorsF (1) N andF (2) N emerge: For the fully-contracted terms of equations (10) and (11), the scalar quantity E ROHF with emerges, which corresponds to the restricted open-shell Hartree-Fock (ROHF) energy for a single high-spin determinant (c.f. [58] ). In this section, the methodology of the presented CCSD equation derivation and implementation engine is presented. In particular, differences to spin orbital CCSD implementations are highlighted. Figure 1 features a flowchart diagram showing the consecutive steps from BCH expansion to the final factorized code for one particular example. The specifics steps are explained in greater detail throughout subsections III A to III D: In subsection III A, the application of Wick's theorem [59] to the spin-free BCH terms is demonstrated. Connecting lines below the second quantized operator strings are used to abbreviate spin summations. The prefactors arising from spin summations are calculated directly from the full contraction patterns. After obtaining Wick terms composed of products of Kronecker deltas, summarizable (redundant) terms need to be collected. In this work, we apply special nonantisymmetrized Goldstone diagrams (defined throughout subsection III B) to achieve this task. These diagrams are defined such that summarizable terms are represented by isomorphic, i.e. topologically equivalent, graphs. Finally, diagram-wise index relations need to be evaluated for each Goldstone topology (c.f. subsection III C). In this work, we introduce graphs -denoted index relation graphs -to evaluate, keep track and analyze these index relations. In subsection III D, the generation of the final factorized C++ code is outlined. Factorization routes are optimized locally, i.e. diagram-wise, by a complexity estimation using index relation graphs. Following the Baker-Campbell-Hausdorff (BCH) expansion of the CC equations (5), sums of products of second quantized creators and annihilators are obtained. When applying Wick's theorem [59] to the operator strings, only fully contracted terms survive. All surviving and vanishing contraction types for the spin orbital CC as well as the spin-adapted and spin-complete CC framework (followed in this work) are shown in Table I . In the spin orbital CC formulation, the cluster operatorT is automatically normal-ordered w.r.t. the Fermi vacuum |Ψ 0 such that no internal contractions within the sameT are possible (T = 0). Furthermore, the usual spin orbital T includes hole and particle creators -â iσ andâ † aσ -only. For spin orbital CC, this results in no survivingTT and TĤ N contractions, which leads to a drastically simplified BCH expansion. Unfortunately, none of these simplification hold true in the general spin-adapted and spin-complete case. Due to the presence of active indices v, w, . . . corresponding to either occupied or virtual spin orbitals in the high-spin BCH expansion Translation to (factorized) C++ code for (i = 0; i < n occ; ++i) for (j = i + 1; j < n occ; ++j) for (a = 0; a < n vir; ++a) I : Surviving and vanishing contraction types in the spin orbital framework (implying O ∩ V = ∅) and in the spin-adapted and spin-complete framework followed in this work. Contraction Spin orbital Spin-adapted reference, some operators inT may contract from the left. Additionally, spectating indices, as e.g. inÊ av vi , lead to non-vanishing internal contractions. Despite these general considerations, Wick's theorem may be applied as in the spin orbital case. However, special care needs to be taken when evaluating spin summations. Consider, e.g., the following term contributing to E corr : From now on, we abbreviate second quantized operators by neglecting theâ such that a † pσ = (pσ) † andâ qσ = (qσ) and introduce a special notation for spin summations: Each spatial index pair (p, q) summed over the same spin variable (σ) will be represented by a connecting line below the operator expression, which we will call spin path. In contrast to spin paths, Wick contractions will be displayed above the operator strings. Using this new notation, the rhs of equation (12) now reads where the notation (. . .) F C shall denote the sum of all full contraction patterns. The expression now contains four unique spin paths containing the pairs pr, qs, ai and bj. Each operator (Wick) contraction -denoted as contraction lines above the operator strings -connecting two distinct spin paths adds a Kronecker delta for the participating spin variables such that the two spin summations are broken down to one. In other words, the spin paths of contracted operators merge. By this spin path tracing, the individual prefactors of arbitrary contraction patterns with regard to spin summation are easily derived. Every connected spin path (representing an isolated spin summation over α and β) in a full contraction pattern, contributes a prefactor of 2. For the given example, it is where the merged spin paths were highlighted in red and blue. The final terms contributing to E corr are therefore Clearly, due to equations (8) and (9), there are exceptions for active indices. Every contracted active index either involves only α (for v † w) or β (for vw † ) spins. Therefore, every path containing at least one active index contributes a prefactor of 1 instead. The expansion of Ê pq rs Ê av ij C , e.g., yields where merged spin paths of active indices are highlighted in red. Lastly, consider the example: Whenever the two active indices v and w are contracted as spin orbital particles (e.g. rv † ) and holes (e.g. p † w) on the same spin path, the whole expression vanishes. This is due to incompatible active spins (c.f. equation 9). While the creator -annihilator contraction requires α spins, the annihilator -creator contraction requires β spins. Both can never be true at the same time. An analogue argument holds for the second one-particle interrelating operatorF (2) N containing spin orbital α −→ α substitutions only. Here, any active creator -annihilator contractions on the same spin path lead to a vanishing expression. In this work, the BCH expansion of the CCSD equations (5) up to quadruply nested commutators was generated. Please note that due to the non-vanishingTT contractions, the BCH series does not naturally truncate at fourth order. We restrict our CCSD approach to BCH order four, nevertheless, since fifth, sixth, etc. order terms are expected to have negligible contributions, i.e. < 10 −7 a.u. in terms of the correlation energy. We will further discuss and motivate this issue in section IV. Application of Wick's theorem to the BCH expansion of the CCSD equations (5) leads to largely redundant products of Kronecker deltas. These bind general indices of H N to specific indices fromT or the external projections. Due to theT andĤ N tensor symmetries, many Wick terms lead to redundant equations. A much more efficient approach to CC equation derivation -established by Čížek and Paldus [43, 60, 61] -is based on the diagrammatic representations developed by Goldstone [62] and Hugenholtz [63] , who generalized Feynman's particle interaction diagrams [64] to many-body perturbation theory. In this work, we apply special spin-free Goldstone diagrams to identify summarizable terms. The general idea is to abstract algebraically redundant terms into isomorphic graphs such that the redundancy check may be performed by the optimized (sub)graph isomorphism testing algorithm VF2 [65, 66] contained in the boost graph library (BGL) [67] . In the following, all essential details of the used diagrams are given. Firstly, the new lines and features to properly represent active and/or identical lines are introduced (c.f. subsection III B 1). Then, in subsection III B 2, the embedding of tensor symmetries without implicit antisymmetrization is motivated and finally, in subsection III B 3, the incorporation of index relations into the diagram's topology is discussed. To properly represent active indices in Goldstone diagrams, a new type of line needs to be introduced. Since active index operators may contract as either holes or particles (c.f. equation 9), the vertical orientation of their diagrammatic representation not only needs to be arbitrary but also needs to be changeable. To accommodate this, we represent active indices by particle/hole lines with hollow arrows. The hollowness shall emphasize the lines not being "pinned" to the surrounding background such that they may change their vertical direction freely. For the operatorÊ va iw , e.g., the Goldstone diagram representation is given bŷ where indices v and w are represented by hollow-arrowed lines with arbitrary/fluid vertical direction. For spectating (e.g.Ê av vi ) or identical (e.g.Ê aa ii ) indices, the pairwise line affiliation needs to be part of the diagram. These pairwise identical indices will be represented by background color-coded lines, where same color indicates same indices: It is obvious that the contraction rules for active indices (c.f. equation 9) are automatically fulfilled by the connection of heads and tails of hollow-arrowed lines analogous to the connection of particle/hole lines. However, some non-trivial findings occur when building CC diagrams in the presented SASC-CC framework: (i) due to the operator generation scheme [38] , which handles spectatingÊ operators explicitly, all non-vanishingT contractions can occur among identical active indices, i.e. spectators, leading to active snakes, only. All other non-spectating operators possess either < or = relations between their active indices such that no accidental spectators are formable. Consider, e.g., the following term contributing to Ψ 0 |F Here, the two active lines corresponding to the spectating index v form an internally connected snake, which connects the two vertices of the cluster line. The final contribution to E corr therefore yields At the same time, no contribution to because indices v and w may never be iden- : (ii) when comparing topologies of spin-free Goldstone diagrams, the direction of the participating active lines must not be taken into account. This is easily motivated by a small example. Consider the terms A and B with where indices v and w connect as spin orbital holes (A: creator -annihilator) and as spin orbital particles (B: annihilator -creator), respectively. Clearly, in the full CCSD energy equation it is E corr ←− A + B = 0 due to where [. . .] + shall denote the anticommutator. In order to find these kinds of redundancies graphically, i.e. identify the Goldstone diagrams of A and B as topologically equivalent, all information about the vertical direction of the active line has to be neglected: While this certainly leads to isomorphic Goldstone diagrams A and B, it also removes information about the sign of the diagrams. In our case, however, this is of no concern because the sign (and also the prefactor) is already apparent from Wick's theorem. As stated before, we merely apply Goldstone diagram representations to identify redundant terms. (iii) Whenever identical indices are contracted to different externals, the whole expression vanishes. This again results directly from operator generation scheme [38] applied in this work, where index relations were explicitly derived (also for external indices). As an example, consider the doubles pro- Due to the contractions I † i = δ Ii and J † i = δ Ji never being compatible with I < J, the term vanishes. Please note that the same argument also holds for the particle contractions Aa † = δ Aa and Ba † = δ Ba in this case. Further remarks on the incorporation of index relations are given in subsection III B 3. (iv) Visually, the second one-particle interrelating operatorF N is therefore not required. In contrast to spin orbital substitution operatorsX, the spatial substitutionsÊ possess a pairwise permutational symmetry only. That iŝ whereP shall denote an arbitrary permutation from the symmetric group S ν . A lot of redundancy in the generated Wick terms, originates from this symmetry. Unfortunately, Goldstone diagrams can not completely incorporate the latter. The Goldstone representation ofT operators is a set of horizontally arranged vertices connected by a bold line. Each vertex is associated with an outgoing and an incoming line representing the annihilator/creator pairs of the operator. In Table II , tensor and Goldstone representations of all permutational equivalent rank [68] one to three cluster operators are shown. Goldstone diagramŝ The only symmetric feature of the shown nonantisymmetrized Goldstone diagrams is the spatial inversion with respect to the center of the cluster line. As clearly visible in Table II , for tensor ranks one and two this exactly preserves the right permutational symmetry by accident. While rank one operators do not possess any permutational symmetry, rank two operators possess just one resembling the transposition of indices 1 and 2 (highlighted by blue and gray vertices). This transposition is identical to an inversion such that the two possible diagrams are topologically equivalent, i.e. isomorphic. Rank three operators, on the other hand, possess 3! = 6 equivalent permutations. By graphical inversion, only three groups of two diagrams are found isomorphic. To resolve this problem in spin orbital theory, one may neglect all particle information by transitioning from Goldstone to Hugenholtz diagrams, as e.g. in where all cluster vertices (of the same cluster line) are merged into a single node. While this transition clearly allows for the required permutational tensor symmetry, it also destroys any information about the pairwise creator/annihilator affiliation. In the spin orbital case, this is of no concern since all participating tensors may be antisymmetrized. In the SASC approach followed in this work, however, it is particularly important to use nonantisymmetrized Goldstone diagrams (as e.g. also discussed in [34] ). Therefore, the transition to Hugenholtz diagrams is not a viable option. A different approach to achieve the correct permutational symmetry, is to replace cluster lines by a totally connected, i.e. complete, subgraph. In general, a complete graph K n is defined by a unique set of n vertices, where each unique pair of vertices is connected by a unique undirected edge. Such graphs are therefore topologically invariant to arbitrary vertex permutations. In this work, we represent all occuring substitution operators by such totally connected cluster lines. When representing the latter graphically by means of Goldstone diagrams, however, a lot of complexity is introduced through the additional cluster lines. To remove this unnecessary complexity, we abbreviate fully connected cluster lines by ring-shaped cluster lines we denote cluster super lines. In Table III , fully connected Goldstone representations as well as their cluster super line abbreviation of rank one to four substitution operators are shown. Spin-complete and linear-independent substitution operators were generated by the use of Löwdin's projection operator method [38] . One of the major consequences of this generative approach is the presence of index relations (orderings) in all substitution operators. By definition [38] , it is for the index set of each individual substitution operator. While this is beneficial in guaranteeing the exactly required amount of amplitudes, it enforces us to distinguish between e.g. t ab ij and t ab ji in Clearly, however, the unlabeled Goldstone diagrams resulting from these two operators are identical: To resolve this problem, we employed a simple solution. For each particle/active/hole line connected to a cluster super line, we place auxiliary vertices in between the cluster vertices and the contracted lines. The labels of these vertices shall then reflect the index relations within the particle/active/hole domain of the associated cluster super line (via 1 < 2 < 3 etc. contributing to the CCSD singles equation. One of the possible diagrams for this term is given by Table IV : Progression of the index relation graphs of the contraction pattern given in equation (15) . The merged indices of each step are colored. where all auxiliary vertices are explicitly shown and the external spectator V is contracted to itself. To evaluate the final expression, all index relations need to be merged according to all contracted indices. This evaluation is a non-trivial task, where special care needs to be taken. We will discuss this issue in greater detail throughout the next subsection (III C). For this particular example (collecting all redundant diagrams), the final contribution to the singles residual R AV V I reads I † J † BA a † b † ji I † J † B A a † b † j i I † J † B A a † b † j i O I J i j I = i J j I = i J = j V A B a b a A = b B B = a A = b A - - -R AV V I ←− +1 I relations at the same time -constitute an impossible situation. Therefore, valid IRGs are directed acyclic graphs in general. The relation i < j < k, l e.g. corresponds to the IRG where index l is not bound to i, j, k at all and therefore represented as an isolated (disconnected) node. In the three separated index spaces O, A and V employed in this work, only contractions within the same space survive. Therefore, three separate IRGs can be used to represent the index relations of all occupied, active and virtual indices, respectively. IRGs are particularly useful to evaluate merged index relations for a given contraction pattern. Consider, e.g., the contraction pattern (15) contributing to the CCSD doubles residual R AB IJ . In Table IV, the IRGs of zero, one and two (full) contractions in the occupied and virtual index spaces are shown. For no contractions at all, the IRGs simply reflect the index symmetry of the participating cluster and external super lines. For the given example it is i < j and I < J in the occupied space and a < b and A < B in the virtual space originating from the operatorsÊ IJ AB and E ab ij , respectively. There are no active IRGs since there are no active indices involved. The first shown contractions involve the indices I and i (displayed in red) as well as A and b (displayed in blue). Each contraction merges two distinct nodes of the IRGs such that all edges connected to either of the nodes now connect to a single merged node. For the occupied indices, nodes I and i are merged to a single node I = i, which now connects to two outgoing edges to the nodes J and j representing the index relation (I = i) < (J, j). A special case originating from the presence of = relations, emerges in the active orbital space A. Whenever such a relation is encountered, multiple IRGs are required to represent the index relations for A alone. Consider, e.g., the full contraction excerpt (assuming the two active contractions are on different spin paths) where indices v and w as well as x and y are coupled via v = w and x = y, respectively. In Table V , the uncontracted, singly and doubly contracted IRGs of the given excerpt are shown. For the uncontracted case, the unaltered index relations of the two cluster super lines are obtained. In this particular case, however, two graphs for each super line are possible. This is because the relation v = w may be expanded to v < w ∨ w < v representing two IRGs. Exactly the same holds for the indices x and y connected by x < y ∨ y < x. To graphically represent these choices, all valid IRGs are gathered in brackets such that the full index relation requires all possible choices. The first contraction (highlighted in red) merges indices v and y. Since the indices are from two separate cluster super lines with two IRG choices each, a total of four merged IRGs is possible. These include the two linear chains x < (v = y) < w and w < (v = y) < x as well as the two v-shape relations (v = y) < (w, x) and (w, x) < (v = y). The second contraction (highlighted in blue) combines indices x and w, which are now on the same IRGs. Through these unions, the two linear chain relations are rendered invalid due to occurring loops. The two v-shape relations, however, remain intact. From the original four IRG possibilities, only two survive as the final index relation (v = y) < (w = x) ∨ (w = x) < (v = y). To arrive at a final CC implementation of the generated Goldstone diagrams, the latter need to be translated to code of a given programming language. To achieve an implementation of somewhat manageable computational cost, it is mandatory to factorize the generated equations. Here, we apply a local, i.e. term/diagram-wise, factorization of each Goldstone diagram. While this is far from the global optimum, it still yields optimal scaling results in O, A and V. One major complication when evaluating the optimal factorization order of a given tensor product is the estimation of computational cost for different contraction routes. This is particularly difficult for routes including non-trivial index relations. As it turns out, the IRGs defined in the previous subsection III C can easily be used to estimate this complexity, i.e. the number of iterations, of the represented index-iterating loops. Each IRG may be decomposed into (i) its disconnected fragments and (ii) its different hierarchy levels. While decomposition (i) merely splits unconnected index relations, decomposition (ii) sorts the related indices according to their dependency from each other. Disconnected fragments scale independently of each other such that their complexities may be multiplied. Consider, e.g., the exemplary IRG in Figure 2 relating the indices i 1 to i 10 . Clearly, there are two disconnected fragments. The overall complexity of the IRG is given by the product of the complexities of its disconnected fragments 1 and 2. Within each fragment, the indices group to different hierarchy levels 1 -4 for fragment 1 and 1 -2 for fragment 2, respectively. A simple approach to sort indices into hierarchy levels is to iteratively extract nodes without incoming lines until no nodes are left. The extracted nodes during each iteration then represent the different hierarchy levels. This procedure, of course, requires the absence of any loops in the IRG. Trivially, any linear chain IRG composed of singly filled hierarchy levels only, possesses a complexity of N m for N elements of the given index space and m IRG nodes. Each additional node, attaches a singly connected, i.e. < related, index to the chain thereby increasing its complexity from N m to N m+1 . Unfortunately, valid IRGs may appear in completely different shapes and sizes as seen by the previous example. Whenever hierarchies containing multiple index nodes are encountered, a complexity evaluation must take all possible relations of these indices to each other into account. Consider, e.g., fragment 2 of the IRG shown in Figure 2 . Here, the two indices i 8 and i 9 are categorized to hierarchy level 1, i.e. have no index relation to each other. Therefore, the complexity of fragment 2 must consider all possible relations of i 8 and i 9 given by (i 8 = i 9 ) < i 10 ∨ i 8 < i 9 < i 10 ∨ i 9 < i 8 < i 10 . These in turn represent new (linearized) IRGs with trivially assignable complexities of N 2 , N 3 and N 3 , respectively. The overall complexity of fragment 2 is therefore given by In this work, arbitrary IRGs are analyzed with respect to their complexity using a recursive algorithm. This algorithm (in pseudocode) is listed as Algorithm 1. Input: disconnected IRG fragment f, size of orbital space N Data: complexity c 1 decompose f into hierarchy levels h 2 if all levels h are composed of single elements then In the presented algorithm, an arbitrary fragment f, the size of the corresponding orbital space N and the bookkeeping variable c for the calculated complexity are used. In line 1, the given fragment f is decomposed into its hierarchy levels h (as shown in figure 2 ). If all hierarchy levels are singly populated, the given IRG fragment corresponds to a linear chain with a trivially assignable complexity. In terms of the algorithm, this constitutes the recursion end. The binomial complexity is added to c (line 3) and the call is returned (line 4). In case of at least one multiple-populated, i.e. degenerate, hierarchy level d, the main recursion loop (lines 7-12) is entered. Here, all possible index relations of the degenerate index nodes are enforced. In this work, we evaluate all possible index relations as all unique combinations of = and < connections between the different nodes. This is because, both = and < possess straight forward IRG representations, which may be applied to copies of f. For each combination of index relations among the degenerate nodes, a separate copy of f is created (line 8). In these copies, all identical nodes, i.e. nodes with a = connection, are merged (line 9) while edges are added to all nodes connected by a < relation (line 10). This leaves the altered copy f_copy, where the degeneracy d was removed. To demonstrate the functionality of the algorithm, the full recursion tree of the exemplary IRG fragment given in Figure 2 as fragment 1 is shown in Figure 3 . Given the initial IRG fragment, degenerate nodes are encountered in hierarchy levels 1 and 2 highlighted in red and blue, respectively. In the first recursion level, the degeneracy of the red nodes is removed. In total, there are 13 different = and/or < relations possible. These include • one possibility for three identical nodes (i = j = k), Total complexity: (1 + 2) · 6 · N 7 + (1 + 2 + 1 + 1) · 6 · N 6 + ((1 + 2) · 1 + (1 + 1) · 6) N 5 Figure 3 : Recursion tree of the exemplary index relation graph given in Figure 2 containing two degenerate hierarchy levels (highlighted by red and blue nodes). The node labels as well as irrelevant edges are hidden to increase readability. • six possibilities for two identical nodes ((i = j) < k, k < (i = j), (i = k) < j, j < (i = k), (j = k) < i and i < (j = k)) and • six possibilities for no identical nodes (i < j < k, i < k < j, j < i < k, j < k < i, k < i < j and k < j < i). Since all the red nodes are equally connected to the next hierarchy level (the blue nodes), the particular ordering of the nodes (in < relations) is not important -it leads to identical complexities. Different graphs only emerge from the different amounts of identical nodes. This is because (1) two nodes connected by two identical lines (no loop) represent the same IRG as two nodes connected by a single line: (2) lines that skip certain hierarchy levels are irrelevant due to transitivity: Both criteria (1) and (2) are applied to the graphs of Figure 3 such that doubled as well as irrelevant lines are neglected. For each of the 13 generated graphs with removed degeneracies in the red hierarchy level, the algorithm proceeds to the degeneracies of the blue nodes representing the next hierarchy level. There are three possibilities for the two blue nodes, which include • one possibility for two identical nodes (i = j) and • two possibilities for no identical nodes (i < j and j < i). In contrast to the red nodes of hierarchy level 1, the ordering of the < relation of the blue nodes matters. If the upper blue node is placed before the lower blue node, a linear chain IRG is obtained. If the upper blue node is placed after the lower blue node, however, a new degenerate hierarchy is formed. This new degenerate hierarchy needs to be eliminated in a third step of the recursion tree leading to three linear chain IRGs. For these, the ordering of the < relations again is of no concern. Finally, the total complexity of the initial IRG fragment may be calculated by the addition of all complexities of the linear chain IRGs of the lowest recursion tree branches, i.e. the recursion endings as shown at the bottom of Figure 3 . Please note that IRGs of identical complexities originating from arbitrary < orderings are collected into groups with an assigned magnitude in Figure 3 . These need to be included multiplicatively when evaluating the total complexity. In this work, the presented complexity evaluation (c.f. Algorithm 1) is used to optimize contraction routes of tensor products. First, each generated Goldstone diagram is translated to its algebraic form. Then all possible tensor products in all possible combinations are evaluated with respect to their complexity. This is done by extracting the participating indices of the tensor product from the full IRG of the diagram. This extracted IRG may then be analyzed w.r.t. its complexity using Algorithm 1. Once all total complexities of all factorization routes per diagram are evaluated, the optimal route is determined. This is done by comparing the calculated complexities for exemplary orbital spaces of given numbers N O , N A and N V and choosing the cheapest route. The optimal routes are then automatically translated to code either adding to E corr or to a residual tensor entry. The generated code is called by a separate CC equation solver, which employs the recently published Newton-Krylov method [69] to accelerate the residual convergence. In this section, the presented spin-adapted and spincomplete (SASC) CCSD implementation is applied to several small molecular test systems. The factorization routes of all CC equations was optimized w.r.t. the examplary orbital space sizes N O = 6, N A = 3 and N V = 18. The effects of spin adaption and/or spin completeness in the non-orthogonal and orthogonal CC framework were already investigated in greater detail in our previous work [39] . Here, we merely present some proof-of-concept applications to show the feasibility of the presented implementation. In subsection IV A, the convergence of the correlation energy w.r.t. the BCH series truncation is investigated. Through several test calculations, a reasonable truncation at quadruply nested commutators is motivated. Subsection IV B features a comparison of SASC-CCSD correlation energies to (spin-contaminated) spin orbital CCSD and CCSD(T). Due to active indices in the substitution operators of the cluster operator (c.f. equations 3 and 4), non-vanishinĝ TT andTĤ N contractions are possible. These lead to a non-terminating BCH series expansion in general. In our previous work [39] , we investigated atomic systems H-N in the 6-31G basis set [70] with different orthogonal and nonorthogonal spin-adapted CC methods. We found that in Table VI : Correlation energies and Euclidean amplitude norms for SASC-CCS and SASC-CCSD for different atomic and molecular systems in the cc-pVDZ basis set [71] for a BCH series truncation of orders 1-4. The geometries of BeH and OH were taken from experimantal data [72] . ROHF calculations were conducted using the PySCF [73, 74] program package with convergence thresholds of 10 −14 a.u. in both energy and density. CC residuals were converged to 10 −7 a.u.. VII: Collection of CC correlation energies for molecular systems in different high-spin states using the cc-pVDZ basis set [71] . The given correlation energies include SASC-CCS and SASC-CCSD as presented in this work and spin orbital CCSD and CCSD(T) from the PySCF [73, 74] program package. For all CC calculations, the percentage of the recovered correlation energy w.r.t. the spin orbital CCSD(T) calculation are given. All geometries were assembled from experimental data. Geometries of BeH, CH, NH, OH, N 2 were assembled from [72] , BH 2 from [75] and CH 2 from [76] . The ground state geometry was used for all calculations. ROHF calculations were also conducted using PySCF employing convergence criteria of 10 −10 a.u. in both energy and density. The CC residuals were converged to 10 −7 a.u. for all calculations. most of the conducted calculations the correlation energy was already properly converged to 13 digits after the decimal point at a BCH truncation of quadruply nested commutators. For all cases, the maximal deviation was 10 −12 a.u. in the correlation energy w.r.t. a quintuply nested commutator truncation. Therefore, in this work, we generated CCSD equations for up to quadruply nested commutators. To support this level of truncation, the convergence of E corr as well as the Euclidean norm of theT 1 andT 2 amplitudes with increasing BCH series truncation level was investigated. The results are collected in Table VI . Already for a BCH series truncation of order 3, all presented correlation energies for both CCS and CCSD are reasonably well converged. In particular, no noticeable change from a truncation of order 3 to order 4 was detected within the first six digits after the decimal point. In comparison to the correlation energy, the amplitude norms present a much more sensitive measure for the BCH convergence. Still, only negligible deviations in bothT 1 andT 2 norms from BCH truncations of order 3 to order 4 were found. It is reasonable to assume that the same rapid convergence holds for comparable systems. Several test calculations of the presented SASC-CCSD were conducted to analyze its functionality. Initially, all calculations reported in our previous work [39] (high-spin states of atomic hydrogen to nitrogen in the 6-31G basis set [70] ), where a reference implementation of the SASC operators represented in the FCI basis was used, were reproduced. The correctness of the presented implementation was confirmed by comparing correlation energies and residuals for the converged amplitudes. Both were found in perfect agreement for up to 13 digits after the decimal point (residual mean squares were converged to 10 −14 a.u.). Therefore, it is reasonable to assume that the generated CC equations and factorizations are correct. Additionally, several test calculations for small molecular systems in the cc-pVDZ basis set [71] were performed. Table VII contains a comparison of SASC-CC to spin orbital CCSD and CCSD(T) correlation energies. In particular, the percentage of recovered correlation energy w.r.t. spin orbital CCSD(T) energies are shown. While SASC-CCS recovered between 0.00% and 25.00%, SASC-CCSD recovered between 95.70% and 99.88% of the CCSD(T) correlation energies. In direct comparison, spin orbital CCSD results of all singlet (S = 0) and all maximal high-spin results (S = N/2) are identical to the SASC-CCSD results. This is because in these cases, spin orbital CCSD is already properly spin-adapted and spin-complete. Hence, the equality of the results acts as another proof of concept for the presented implementation. In all other SASC-CCSD results, a small but noticeable improvement of the correlation energies compared to spin orbital CCSD is detected. These improvements reach from roughly 0.001 mE H to 0.886 mE H (0.01% to 0.98% w.r.t. spin orbital CCSD(T)) for octet OH and quartet BH 2 , respectively. While correlation energy differences of up to 1 mE H might be significant when aiming for chemical accuracy, the effects of the SASC framework on the correlation energy seems to be negligible in comparison to the added complexity through the incorporation of spin-adapted operators into the CC framework for most calculations. Despite these findings, we still expect SASC-CC to produce superior wave functions when compared to spin orbital CC. These may be of particular importance when calculating spin-dependent molecular properties. A similar train of thought was also followed in our previous work [38, 39] as well as e.g. by Datta and Gauss [37] , who employ the spin-adapted CC framework to predict hyperfine coupling tensors. A rigorously spin-adapted and spin-complete coupled cluster singles and doubles (CCSD) implementation capable of treating arbitrary high-spin open-shell states was reported. Following our previous work [38] , we employed the generated, Löwdin-type substitution operators (c.f. equations 3 and 4) to ensure proper spin adaption as well as spin completeness of the CC wave function. While we focus on the CCSD truncation of the full CC operator for this work, the presented scheme is directly applicable to higher substitution ranks. In section III, the generation of factorized CC equations was outlined. From the initial BCH expansion of the CCSD equations, Wick terms were generated using Wick's theorem. In subsection III A, we introduced a scheme to generate Wick terms on the basis of spatial rather than spin orbital indices. Through connecting lines below the second quantized operator string, we were able to evaluate prefactors arising from spin summations. To remove redundancies of the generated Wick terms, we represented the latter by specially crafted spin-free Goldstone diagrams such that summarizable terms corresponded to isomorphic graphs (c.f. subsection III B). In order to do so, we added new diagrammatic features for active and identical lines to our Goldstone diagram definition. We then introduced cluster super lines to incorporate the special tensor symmetry required for the spin-adapted substitution operators. Finally, we added auxiliary vertices to our topologies to represent relations between the participating cluster operator indices. In subsection III C, we introduced index relation graphs (IRGs) to evaluate, keep track and analyze index relations occurring in the tensor products of the CC equations. These IRGs may then be used to estimate complexities of different tensor products (c.f. subsection III D). Each single Goldstone diagram was translated to factorized C++ code employing an optimal factorization route determined through complexity calculations of all possible factorization routes. The final codes were then compiled in batch-wise dynamically linkable library files to be called by a separate spin-adapted and spin-complete (SASC) CC driver program. The developed SASC-CCSD method was successfully applied to several small test molecules throughout section IV. In subsection IV A, the convergence of SASC-CCS and SASC-CCSD methods w.r.t. the BCH series truncation was investigated for three examples. Due to active indices in the cluster operator (representing both occupied and virtual orbitals in the spinorbital CC picture), the BCH series does not terminate after quadruply nested commutators in general. In all presented examples, however, proper convergence of the correlation energy as well as the Euclidean amplitude norms to an accuracy of at least 10 −7 a.u. was reported. These findings agree with our previous analysis for atomic test calculations [39] . In subsection IV B, SASC-CCS and SASC-CCSD correlation energies were compared to spin orbital CCSD and CCSD(T) correlation energies. For all singlet (S = 0) and maximum high-spin (S = N/2) calculations, SASC-CCSD reproduced the spin orbital CCSD results. Since spin orbital CCSD is properly spin-adapted and spincomplete in those cases, these findings were expected and act as a proof of concept for the presented implementation. In all other calculations, small improvements of the spin orbital CCSD correlation energies were found for SASC-CCSD. While the effects of spin adaption and spin completeness seem to be small for the correlation energy (as also discussed in e.g. [37] [38] [39] ), they are expected to be significant for molecular properties, spin-dependent properties especially. 3rd IAPR-TC15 Workshop on Graph-based Representations in Pattern Recognition The Boost Graph Library: User Guide and Reference Manual Beware that rank implies tensor rank and not substitution order. A rank two operator may very well resemble a spatial substitution of order one Molecular Spectra and Molecular Structure: IV. Constants of Diatomic Molecules Electronic spectra and electronic structure of polyatomic molecules Structure of Free Polyatomic Molecules, 1