The required one-and two-electron integrals (Equation ( ))-as well as the AO-overlap integrals needed in Equation (49)-were calculated using the PySCF module. In order to calculate the V , J, K, and O integrals using the above equations, the density matrix coefficients and the one-and two-electron integrals have to be given in the same basis set. However, the ordering and normalization conventions of the basis functions differs between Gaussian and PySCF. Therefore, after extracting and computing all the density-matrix coefficients from Gaussian's output, we transform them to match the PySCF convention. The AO integrals computed by PySCF were then used without any modification. The fields are gray-white shaded to ease the distinguishing between different site states. The matrices of the density-matrix coefficients that are read directly from Gaussian's fchk file are coloured blue, while those that are calculated from the CI and MO coefficients of the site states (read from Gaussian's rwf and fchk file, respectively) are coloured red. The remaining non-zero coefficients are coloured black and are calculated from either blue or red ones. The density-matrix coefficients for the pairs of states in the lower triangle of the figure are obtained by transposing the matrices of the density-matrix coefficients from the upper triangle. While transposing, the spin labels also swap positions, e.g. D βα (T 1 2 , S 1 ) = D αβ (S 1 , T 1 2 ) T . Note that the figure can trivially be generalized from S 1 , S 2 , T 1 , and T 2 to any pair of excited singlet/triplet site states.