Coupled-cluster methods

Coupled-cluster methods

The following coupled-cluster methods are implemented in ElemCo.jl:

• ccsd - closed-shell implementation, for open-shell systems defaults to uccsd,
• uccsd - unrestricted implementation,
• rccsd - restricted implementation (for high-spin RHF reference only),
• ccsd(t) - closed-shell implementation,
• dcsd - closed-shell implementation, for open-shell systems defaults to udcsd,
• udcsd - unrestricted implementation,
• rdcsd - restricted implementation (for high-spin RHF reference only),
• λccsd - calculation of Lagrange multipliers, closed-shell implementation,
• λccsd(t) - closed-shell implementation,
• λdcsd - calculation of Lagrange multipliers, closed-shell implementation.

The most efficient version of closed-shell CCSD/DCSD in ElemCo.jl combines the dressed factorization from ^[Kats2013] with the cckext type of factorization from ^[Hampel1992] and is given by

\[\begin{align*} \mathcal{L} &= v_{kl}^{cd} \tilde T^{kl}_{cd} + \left(\hat f_k^c + f_k^c\right) T^k_c + Λ_{ij}^{ab} \left(\hat v_{kl}^{ij} \red{+ v_{kl}^{cd} T^{ij}_{cd}}\right) T^{kl}_{ab} + Λ_{ij}^{ab} R^{ij}_{pq} δ_a^p δ_b^q \red{+Λ_{ij}^{ab} v_{kl}^{cd}T^{kj}_{ad}T^{il}_{cb}}\\ &+ Λ_{ij}^{ab} \mathcal{P}(ai;bj)\left\{\left(\hat f_a^c - \red{2\times}\frac{1}{2}v_{kl}^{cd} \tilde T^{kl}_{ad}\right)T^{ij}_{cb} - \left(\hat f_k^i + \red{2\times}\frac{1}{2}v_{kl}^{cd}\tilde T^{il}_{cd}\right)T^{kj}_{ab} \right.\\ &+ \left(\hat v_{al}^{id} + \frac{1}{2} v_{kl}^{cd}\tilde T^{ik}_{ac}\right)\tilde T^{lj}_{db} - \hat v_{ka}^{ic} T^{kj}_{cb} -\hat v_{kb}^{ic} T^{kj}_{ac} \red{-v_{kl}^{cd}T^{ki}_{da}\left(T^{lj}_{cb}-T^{lj}_{bc}\right)}\\ &\left.- R^{ij}_{pq} \left(δ_k^p δ_b^q - \frac{1}{2} δ_k^p δ_l^q T^l_b\right) T^k_a \right\} +Λ_i^a R^{ij}_{pq}\left( 2δ_a^p δ_j^q - δ_j^p δ_a^q \right) -Λ_i^a T^k_a R^{ij}_{pq}\left( 2δ_k^p δ_j^q - δ_j^p δ_k^q \right)\\ &+Λ_i^a \hat h_a^i + Λ_i^a \hat f_j^b \tilde T^{ij}_{ab} - Λ_i^a \hat v_{jk}^{ic} \tilde T^{kj}_{ca}, \end{align*}\]

where

\[R^{ij}_{pq} = v_{pq}^{rs} \left(\left(T^{ij}_{ab}+T^i_a T^j_b\right)δ_r^a δ_s^b +δ_r^i T^j_b δ_s^b + T^i_a δ_r^a δ_s^j + δ_r^i δ_s^j \right). \]

The DCSD Lagrangian is obtained by removing terms in red. Integrals with hats are dressed integrals, i.e. they are obtained by dressing the integrals with the singles amplitudes, e.g., $\hat v_{kl}^{id} = v_{kl}^{id} + v_{kl}^{cd} T^i_c$.

ao_cc_setup!(EC::ECInfo; closed_shell::Bool, orbitals=nothing) -> EHF

Set up an AO-direct CC run. Freezes core / deleted / frozen-virtual orbitals (reducing EC.space to the active space, exactly like the DF/MO path via freeze_orbitals!) and folds the frozen-core mean field into an effective one-electron AO Hamiltonian ("h1eff_AA" closed-shell; per-spin "h1eff_mm_AA"/"h1eff_MM_AA" open-shell) — the AO-direct analogue of freeze_orbs_in_dump; the frozen-core energy is added to the reference. Then builds the bare (undressed) active-space quantities the run needs — the MO Fock (f_mm[/f_MM]/e_m[/e_M]), 1-e Hamiltonian, and ⟨ij|ab⟩ — from the AO integrals and the stored MO coefficients (via ao_dressed_ints / ao_dressed_ints_unrestricted with T1=∅). Returns the reference HF energy.

closed_shell is the spin treatment of the residual the caller will run, not a property of the stored orbitals: an unrestricted residual on restricted (RHF/ROHF) orbitals simply duplicates α into β (unrestrict!, as uhf does with a restricted guess) and builds the per-spin reference from two identical orbital sets. The opposite combination (closed-shell residual, UHF orbitals) is a different calculation and is rejected by the caller (ccdriver derives an MO dump instead). orbitals lets the caller hand over the (cMO, classes) pair it already loaded.

ao_direct_orbitals_spin(EC::ECInfo) -> SpinMatrix

The correlation reference orbitals of an AO-direct run — the SINGLE orbital source every AO-direct consumer must use (the dressing, the kext rotation, the Λ generalized-Fock terms, the setup), so that no two of them can end up in different bases.

Once ao_cc_setup! has settled the reference it is read back from C_Am/C_AM (save_ao_correlation_orbitals!), so every consumer — the setup itself, the Rot the residual hoists, the Λ generalized-Fock terms, the OQV rotation — sees the SAME orbitals even when the reference is not a pure function of the orbital file (wf.dump4core_only).

build_ht_mo_blocks!(EC, names; Rv=nothing)

Build the closed-shell bare MO blocks names (e.g. ("vvvo","ovoo") for (T), plus "vovv"/"ooov" for Λ(T)) that the AO-direct triples / λ-triples read, from the occupied-bra half-transformed store "ht_oAAA" (still on file from ao_cc_setup!) and the active-space MO coefficients. The consumers pick the blocks up with load4idx.

Rv is an optional virtual-space rotation (e.g., for the pseudo-canonicalization)

build_ht_mo_blocks_unrestricted!(EC, names; Rv=nothing)

Unrestricted analogue of build_ht_mo_blocks!: build the bare MO blocks names (same-spin vvvo/vooo/VVVO/VOOO and the opposite-spin vVvO/vVoV/vOoO/oVoO the unrestricted (T) reads) from the per-spin half-transformed stores "ht_oAAA_a"/"ht_oAAA_b" (built in ao_cc_setup!) and the per-spin MO coefficients. Each block's ht_mo_block_spec says which store supplies its occupied index and which spin's coefficients go on the free slots.

Rv is an optional per-spin virtual-space rotation (e.g., for the pseudo-canonicalization)

calc_1RDM(EC::ECInfo, U1, U2, T1, T2; jacobian=false)

Calculate the 1RDM for the closed-shell CCSD or DCSD equations.

Return D1[p,q]=$D_p^q$, the 1RDM without T1 singles terms, and dD1[p,q]=$\hat D_p^q$, the 1RDM with all T1 terms included.

if jacobian=true, the T1 contributions to dD1 are calculated without the +2 T1 term, to be used in the calculation of the CCSD Jacobian.

calc_1RDM(EC::ECInfo, U1, U1os, U2, U2ab, T1, T2, T2ab, spin; jacobian=false)

Calculate the spin-1RDM for the unrestricted CCSD or DCSD equations.

U1, U2, T1, T2 are the Lagrange multipliers and amplitudes for spin∈{:α,:β}, U1os are the singles Lagrange multipliers for opposite spin, and U2ab, T2ab are the αβ Lagrange multipliers and amplitudes.

Return D1[p,q]=$D_p^q$, the 1RDM without T1 singles terms, and dD1[p,q]=$\hat D_p^q$, the 1RDM with all T1 terms included.

if jacobian=true, the T1 contributions to dD1 are calculated without the + T1 term, to be used in the calculation of the CCSD Jacobian.

calc_MP2(EC::ECInfo, addsingles=true)

Calculate closed-shell MP2 energy and amplitudes. The amplitudes are stored in T_vvoo file. If addsingles: singles are also calculated and stored in T_vo file. Return EMp2 OutDict with keys (E, ESS, EOS, EO).

calc_UMP2(EC::ECInfo, addsingles=true)

Calculate unrestricted MP2 energy and amplitudes. The amplitudes are stored in T_vvoo, T_VVOO, and T_vVoO files. If addsingles: singles are also calculated and stored in T_vo and T_VO files. Return EMp2 OutDict with keys (E, ESS, EOS, EO).

calc_UMP2_energy(EC::ECInfo, addsingles=true)

Calculate open-shell MP2 energy from precalculated amplitudes. If addsingles: singles energy is also calculated. Return EMp2 OutDict with keys (E, ESS, EOS, EO).

calc_cc(EC::ECInfo, method::ECMethod)

Calculate coupled cluster amplitudes.

Exact specification of the method is given by method. Returns energies ::OutDict with the following keys:

• "E" - correlation energy
• "ESS" - same-spin component
• "EOS" - opposite-spin component
• "EO" - open-shell component (defined as $E_{αα} - E_{ββ}$)
• "EIAS" - internal-active singles (for 2D methods)
• "EW" - singlet/triplet energy contribution (for 2D methods)

calc_ccsd_vector_times_Jacobian(EC::ECInfo, U1, U2; dc=false, with_rhs=true)

Calculate the vector times the Jacobian for the closed-shell CCSD or DCSD equations.

if with_rhs is true, the right-hand side of the Lambda equations is also added. Return R1, R2

calc_ccsd_vector_times_Jacobian(EC::ECInfo, U1a, U1b, U2a, U2b, U2ab; dc=false, with_rhs=true)

Calculate the vector times the Jacobian for the unrestricted CCSD or DCSD equations.

if with_rhs is true, the right-hand side of the Lambda equations is also added. Return R1a, R1b, R2a, R2b, R2ab

calc_intermediates4Jacobian(EC::ECInfo, method::ECMethod)

Calculate intermediates required in calc_ccsd_vector_times_Jacobian for the closed-shell CCSD or DCSD equations.

calc_lm_cc(EC::ECInfo, method::ECMethod)

Solve coupled cluster Lagrange multipliers and persist the converged U_* tensors for downstream post-processing.

Exact specification of the method is given by method.

calc_pertT(EC::ECInfo, method::ECMethod; save_t3=false)

Calculate (T) correction for [Λ][U]CCSD(T)

Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict. If save_t3 is true, the T3 amplitudes are saved in T_vvvooo file (only for closed-shell).

calc_posMP2(EC::ECInfo, addsingles=true)

Calculate restricted MP2 energy and amplitudes with positrons. The amplitudes are stored in T_vvoo and T_veop files. Return EMp2 OutDict with keys (E, ESS, EOS, EO).

Fraction of deleted orbitals above which ao_cc_setup! suggests the MO-based route. AO-direct always works in the FULL AO dimension, while a derived MO basis correlates in the reduced one, so the AO-direct overhead grows roughly like (nao/(nao-ndel))^4 — ≈2.4× at 20%.

rs-block length the 4-external contractions fall back to: the default of the shared get_spaceblocks, so that kext_rs_blocksize reproduces the previous blocking exactly wherever it does not grow it.

Growth of the rs-block length per result column beyond KEXT_RS_MINBLOCK, see kext_rs_blocksize. Measured on the calc_K2 GEMM shape (MKL, 8 threads, norb = 96…232): the block GEMM saturates around k ≈ 5–8 times the number of result columns.

Result width below which the blocking is left at KEXT_RS_MINBLOCK, see kext_rs_blocksize. Deliberately conservative: block lengths in 129…511 measured slower than the old default at norb=164, ncols=144 (−6%, where the 128-block GEMM already runs at 93–100% of its shape roofline), and the sweep is not monotonic. Ramping only above ncols = 256 puts the smallest block the ramp can produce at 128 + 8·129 = 1160, so that band is unreachable by construction rather than by a tuned special case. The cost is some forgone gain for ncols ∈ (128, 256]; the benefit is that no shape can regress.

PseudoCanonicalTransform

Holds transformation matrices for pseudo-canonicalization of amplitudes and integrals.

For biorthogonal systems, the Fock matrix is non-Hermitian and requires separate left (L) and right (R) eigenvector matrices: F = R * Diagonal(ϵ) * L^†.

The transformation convention is:

• Lower indices (first half of index string) → transformed using Left eigenvectors
• Upper indices (second half of index string) → transformed using Right eigenvectors

PseudoCanonicalTransform(EC::ECInfo; restricted::Bool=true)

Construct a PseudoCanonicalTransform by checking and diagonalizing Fock matrix blocks.

If the Fock matrix is already diagonal (within threshold), returns identity transforms. For non-diagonal blocks, computes the eigenvectors for pseudo-canonicalization.

SVD_triples_to_singles_and_doubles_residuals(EC)

Calculate contributions from triples to singles and doubles residuals.

ao_core_fock(EC::ECInfo, Dcore::AbstractMatrix) -> Matrix

Closed-shell mean-field (2J−K) AO Fock contribution of the density Dcore (spatial, one particle per orbital), built directly from the exact AO integrals: $F_{μν} = Σ_{ρσ} (2⟨μρ|νσ⟩ − ⟨μρ|σν⟩) D_{ρσ}$. Used to fold the frozen core into an effective one-electron Hamiltonian in ao_cc_setup!.

ao_core_ufock(EC::ECInfo, Da, Db) -> (Fa, Fb)

Unrestricted mean-field AO Fock contributions of the α/β core densities Da/Db: the Coulomb term uses the total density, the exchange the same-spin density — $F^α = J[D_α+D_β] − K[D_α]$, $F^β = J[D_α+D_β] − K[D_β]$, with $J[D]_{pq}=⟨pr|qs⟩D_{rs}$ and $K[D]_{pq}=⟨pr|sq⟩D_{rs}$. Used to fold the (per-spin) frozen core into effective one-electron Hamiltonians in ao_cc_setup!.

ao_direct_orbitals(EC::ECInfo) -> Matrix

MO coefficients for an AO-direct run: the (α) correlation reference orbitals (ao_direct_orbitals_spin) with the linearly-dependent (deleted/redundant) columns — the last n_deleted_orbitals — dropped, so the AO-direct MO space excludes them. They carry no electrons and their (near-)zero coefficients would otherwise corrupt the Rot AO↔MO rotation in cc_kext!. The dropped columns are the highest orbital indices (appended by the canonical orthogonalization, and by the completion of a projected set).

ao_dressed_coeffs(cMO, T1, occ, virt) -> (CL, CR)

T1-dressed bra/ket MO coefficient sets: CL = [C_o | C_v − C_o·T1ᵀ] (bra/particle) and CR = [C_o + C_v·T1 | C_v] (ket/hole). Empty T1 ⇒ CL = CR = cMO.

ao_dressed_ints(EC::ECInfo, T1, cMO)

Build the T1-dressed integrals the closed-shell calc_cc_resid needs in its use_kext path — dh_mm, df_mm, d_oovo, d_oovv, d_oooo, d_voov, d_vovo — directly from the exact AO integrals (the ± supermatrix store and "h_AA") and MO coefficients cMO, without a transformed MO dump. The 4-external (vvvv) term is left to cc_kext!, which contracts the AO integrals directly with Rot=cMO.

The (non-unitary) T1 similarity and the AO→MO transform are folded into two dressed coefficient sets — bra/particle and ket/hole:

  C̃ᴸ = [ C_o | C_v − C_o·T1ᵀ ],   C̃ᴿ = [ C_o + C_v·T1 | C_v ],

giving dressed ⟨pq|rs⟩ = Σ ⟨μν|ρσ⟩ C̃ᴸ[μ,p] C̃ᴸ[ν,q] C̃ᴿ[ρ,r] C̃ᴿ[σ,s] (bra indices use C̃ᴸ, ket indices C̃ᴿ).

Every block needed has at least two occupied indices, so the transform goes to the occupied space as early as possible. Looping once over the AO integrals (one σ = ket-2 slab at a time, read straight from the triangular storage — no detri_int2), we build three intermediates that each carry exactly two occupied indices and the retained σ:

  ooAA ← v_ooAA[i,j,ρ,σ]  = Σ ⟨μν|ρσ⟩ C̃ᴸ_o[μ,i] C̃ᴸ_o[ν,j]     → d_oooo / d_oovo / d_oovv
  AooA ← v_AooA[μ,i,j,σ] = Σ ⟨μν|ρσ⟩ C̃ᴸ_o[ν,i] C̃ᴿ_o[ρ,j]     → d_voov / d_vooo
  oAoA ← v_oAoA[i,ν,j,σ] = Σ ⟨μν|ρσ⟩ C̃ᴸ_o[μ,i] C̃ᴿ_o[ρ,j]     → d_vovo

d_vovo[a,i,b,j] = ⟨ai|bj⟩ = ⟨ia|jb⟩ is obtained from oAoA by electron-exchange symmetry (so it, like the others, keeps σ and avoids a separate ket-2 contraction). The remaining two AO indices of each intermediate are transformed last, only to the spaces the blocks need — the full nao⁴ tensor, the full MO integrals, and every all-virtual block are never formed. Occupied-restricted indices use i,j,k,l,m,n.

ao_dressed_ints_unrestricted(EC::ECInfo, T1a, T1b, cMOa, cMOb)

Unrestricted (UHF) analogue of ao_dressed_ints: build the T1-dressed αα/ββ/αβ integral blocks the open-shell use_kext residual and the dressed Fock need, directly from the exact AO integrals (the ± supermatrix store) and the per-spin effective 1-e Hamiltonians "h1eff_mm_AA"/"h1eff_MM_AA" (frozen core folded per spin). The 4-external (vvvv) term is left to the unrestricted cc_kext! with Rota=cMOa, Rotb=cMOb. Uses the occ-early passes ao_ss_blocks/ao_os_blocks — no dense nao⁴ tensor is formed. Empty T1a/T1b ⇒ bare blocks.

ao_lagrange_K2(EC, T1, U2, o4s, v4s, spin=:α) -> Kmmoo

AO-direct K_{mn}^{rs} = Σ_pq ⟨pq|rs⟩ dU2_{mn}^{pq} for the closed-shell / same-spin Λ kext (spin selects the orbital set for the unrestricted :β call), mirroring the amplitude kext (cc_kext!). The MO-space dressing calc_dU2 is skipped: the AO dressed-Λ2 density is folded directly, dU2_AO[μ,ν,m,n] = Σ_ab CLv[μ,a] CLv[ν,b] U2[a,b,m,n] with the dressed virtual bra CLv = C_v − C_o·T1ᵀ (same as ao_dressed_coeffs). By ⟨pq|rs⟩=⟨rs|pq⟩ this is pm_K2!'s ket-pair contraction: ½-scale the μν diagonal (scalepp), pack the pair triangularly, contract against the ± store, and rotate the kept ρσ pair back to MO with cMO. Real-valued (the bra↔ket reuse of the ± store); complex AO-direct is a follow-up.

ao_lagrange_K2ab(EC, T1a, T1b, U2ab) -> KmMoO

Opposite-spin (αβ) analogue of ao_lagrange_K2: K_{mN}^{rS} = Σ_pQ ⟨pQ|rS⟩ dU2_{mN}^{pQ}. The αβ dressed Λ2 density folds directly with BOTH dressed virtual bras, dU2_AO[μ,ν,m,N] = Σ_aB CLva[μ,a] CLvb[ν,B] U2ab[a,B,m,N] (the AO image of calc_dU2(EC,T1a,T1b,U2ab,'o','v','O','V')), so calc_dU2 is skipped. The pair is not symmetric, so — exactly as the amplitude cc_kext! αβ branch — the full square with a ½-scaled diagonal goes to pm_K2ab! (which does the explicit ± fold over both pair orders); the two kept AO externals are then folded back with the α and β orbitals respectively. Real-valued, like the same-spin variant.

ao_os_blocks(pm, La_o,La_v,Ra_o,Ra_v, Lb_o,Lb_v,Rb_o,Rb_v) -> NamedTuple

Opposite-spin (αβ) counterpart of ao_ss_blocks: one pm_os_sweep over the ± supermatrix store builds the occupied-contracted intermediates for both spins, then only the remaining AO indices are transformed into the needed spaces. No nao⁴ tensor is formed.

Transform the remaining AO indices of the five opposite-spin intermediates into the dressed αβ blocks.

ao_rotate_ints(EC::ECInfo, R::AbstractMatrix) -> Crot

AO-direct analogue of rotate_ints for the orbital-optimized QV methods: rather than re-transforming an MO integral dump, the orbital rotation is folded into the MO coefficients (Crot = cMO·R) and the (bare) integral blocks are rebuilt straight from the AO integrals in the rotated basis — the same "fold the rotation into the coefficients" route the T1 dressing and the Λ2 kext use.

The general-orbital block stays in the AO space (save_ao_AAAo!): calc_oqv_gradient only ever contracts its three general indices with the virtual rotation, so it is indifferent to whether they are MO or AO — it just receives Crot (whose virtual columns are the AO→rotated-virtual coefficients) in place of R. Returns Crot, which the caller also passes on as the AO→MO map for the residual and cc_kext!.

ao_ss_blocks(pm, Lo, Lv, Ro, Rv) -> NamedTuple

Same-spin occ-early pass (the closed-shell ao_dressed_ints kernel, reused per spin): one pm_occ_early sweep over the ± supermatrix store builds three occupied-contracted intermediates, then only the two remaining AO indices are transformed into the needed spaces. Returns the dressed oooo/oovo/oovv/voov/vovo/vooo blocks (bra columns from the dressed Lo,Lv, ket from Ro,Rv). No nao⁴ tensor and no all-virtual block is ever formed.

Transform the remaining AO indices of the occ-early intermediates into the dressed same-spin blocks.

bare_ints2!(out, EC::ECInfo, spaces::String)

Fill out with the bare MO integral block spaces, caching it on scratch so that the CC iterations extract it from the integral dump only once.

Blocks larger than a quarter of the integral dump are NOT cached: at that size the copy stops being cheap next to what it saves, and it keeps a pathological case (vvvv at nvirt^4) from being written to scratch. Such blocks only arise with cc.calc_d_vvvv and friends, which the use_kext residual does not use.

The cache file carries the plain space string as its name — bare integrals are never prefixed (only dressed blocks are, d_*), and these are the same names the AO-direct engine blocks use (save_mo_blocks!); both are temporary, so the two routes can never serve each other a stale file.

The cache is temporary (description="tmp"), so delete_temporary_files! drops it at the end of every driver run: a regenerated dump can never be served a stale block. A cached file whose dimensions no longer match (a changed space) is silently regenerated.

calc_3ext_times_T2(EC::ECInfo, T2::AbstractArray, o1::Char='o', v1::Char='v', o2::Char='o', v2::Char='v')

Calculate $\hat v_{mb}^{cd} T^{ij}_{cd}$ intermediate required in calc_ccsd_vector_times_Jacobian and store as vT_ovoo[m,b,i,j].

calc_3ext_times_T2(EC::ECInfo, T2a::AbstractArray, T2b::AbstractArray, T2ab::AbstractArray)

Calculate $\hat v_{mb}^{cd} T^{ij}_{cd}$ intermediates required in calc_ccsd_vector_times_Jacobian for αα, ββ, and αβ amplitudes and store as vT_ovoo[m,b,i,j], vT_OVOO[M,B,I,J], vT_oVoO[m,B,i,J] and vT_OvOo[M,b,I,j].

calc_4idx_T3T3_XY(EC::ECInfo, T2, UvoX, ϵX)

Calculate $D^{ij}_{ab} = T^i_{aXY} T^j_{bXY}$ using half-decomposed imaginary-shifted perturbative triple amplitudes $T^i_{aXY}$ from T2 (and UvoX)

calc_D2(EC::ECInfo, T1, T2, spin::Symbol)

Calculate $^{σσ}D^{ij}_{pq} = T^{ij}_{cd} + P_{ij}(T^i_c T^j_d +δ_{ik} T^j_d + T^i_c δ_{jl} + δ_{ik} δ_{jl})$ with $P_{ij} X_{ij} = X_{ij} - X_{ji}$. Return as D[pqij]

calc_D2(EC::ECInfo, T1, T2, scalepp=false; Rot=zeros(Float64,0,0))

Calculate $D^{ij}_{pq} = T^{ij}_{cd} + T^i_c T^j_d +δ_{ik} T^j_d + T^i_c δ_{jl} + δ_{ik} δ_{jl}$.

If scalepp: D[ppij] elements are scaled by 0.5 (for triangular summation). If Rot is provided, the pq indices of D2 are rotated (e.g., to AO basis).

Return as D[pqij]

calc_D2ab(EC::ECInfo, T1a, T1b, T2ab, scalepp=false)

Calculate $^{αβ}D^{ij}_{pq} = T^{ij}_{cd} + T^i_c T^j_d +δ_{ik} T^j_d + T^i_c δ_{jl} + δ_{ik} δ_{jl}$ Return as D[pqij]

If scalepp: D[ppij] elements are scaled by 0.5 (for triangular summation)

calc_E_Coe(eigenvalue_vector, q)

Calculate the coefficient matrix in QV-CCD/DCD residual calculation.

calc_K2(int2, D2, tripp; symmetrize=true)

Calculate the kext K2 contribution to the CCSD residuals by directly contracting the integrals with D2.

$K^{ij}_{pq} = v_{pq}^{rs} D^ij_rs$

Return K2pq::Array{4}.

calc_K2ab(int2, D2)

Calculate the kext K2ab contribution to the CCSD residuals by directly contracting the integrals with D2.

$K^{iJ}_{pQ} = v_{pQ}^{rS} D^iJ_rS$

Return K2pq::Array{4}.

calc_R_from_U_F(e::Vector, X::Matrix, F::Matrix, q::Float64)

Calculate intermediate R with F and eigenvalue e, eigenvectors X of corresponding U in QV-CCD/DCD.

calc_cc_resid(EC::ECInfo, T1, T2, T3; dc=false, tworef=false, fixref=false)

Calculate CCSDT or DC-CCSDT residual.

calc_cc_resid(EC::ECInfo, T1, T2; dc=false, tworef=false, fixref=false, linearized=false, Rot=zeros(Float64,0,0), singles_resid=false)

Calculate CCSD or DCSD closed-shell residual.

calc_cc_resid(EC::ECInfo, T1a, T1b, T2a, T2b, T2ab, T3aaa, T3bbb, T3abb, T3aab; dc=false)

Calculate UCCSDT or UDC-CCSDT residual.

calc_cc_resid(EC::ECInfo, T1a, T1b, T2a, T2b, T2ab; dc=false, tworef=false, fixref=false)

Calculate UCCSD or UDCSD residual.

calc_ccsd_vector_times_Jacobian4ab(EC::ECInfo, U1a, U1b, U2a, U2b, U2ab, D1a, D1b; dc=false, with_rhs=true)

Calculate the left vector times the CCSD/DCSD Jacobian for αβ component. Additionally, remaining contributions to the singles residual are calculated.

if with_rhs is true, the right-hand side of the Lambda equations is also added. Return ΔR1a, ΔR1b, R2ab

calc_ccsd_vector_times_Jacobian4spin(EC::ECInfo, U1a, U1b, U2a, U2b, U2ab, D1, dD1, dD1os, spin; dc=false, with_rhs=true)

Calculate the vector times the CCSD/DCSD Jacobian for the given spin (same-spin residual for doubles). The singles residual is missing some terms which are added in calc_ccsd_vector_times_Jacobian4ab.

if with_rhs is true, the right-hand side of the Lambda equations is also added. Return R1 and R2

calc_ccsdt(EC::ECInfo, useT3=false, cc3=false)

Calculate decomposed closed-shell DC-CCSDT amplitudes.

If useT3: (T) amplitudes from a preceding calculations will be used as starting guess. If cc3: calculate CC3 amplitudes.

calc_correlated_1rdm(EC::ECInfo, method::ECMethod,
                     T1a, T1b, T2a, T2b, T2ab,
                     U1a, U1b, U2a, U2b, U2ab)

Build the full-space unrestricted correlated 1-RDM from amplitudes and converged Lagrange multipliers.

calc_correlated_1rdm(EC::ECInfo, method::ECMethod, T1, T2, U1, U2)

Build the full-space restricted correlated 1-RDM from amplitudes and converged Lagrange multipliers.

calc_correlated_1rdm(EC::ECInfo, method::ECMethod)

Build the correlated 1-RDM from the saved CC amplitudes and the converged Lagrange multipliers stored in the temporary U_* files.

calc_correlation_norm(EC::ECInfo, U1, U2)

Calculate the norm of the correlation part of the CCSD or DCSD equations using Lagrange multipliers and amplitudes.

Return $⟨Λ|Ψ⟩ = ⟨Λ_1|T_1⟩ + ⟨Λ_2|T_2+\frac{1}{2}T_1 T_1⟩$.

calc_correlation_norm(EC::ECInfo, U1a, U1b, U2a, U2b, U2ab)

Calculate the norm of the correlation part of the UCCSD or UDCSD equations using Lagrange multipliers and amplitudes.

Return $⟨Λ|Ψ⟩ = ⟨Λ_1|T_1⟩ + ⟨Λ_2|T_2+\frac{1}{2}T_1 T_1⟩$.

calc_dU2(EC::ECInfo, T1, T12, U2, o1='o', v1='v', o2='o', v2='v')

Calculate the "dressed" $Λ_2$ for CCSD/DCSD.

T12 is the T1 amplitude for the second electron of U2 (=T1 for closed-shell and same-spin U2). Return dU2[p,q,m,n]=$Λ_{mn}^{ab}δ_a^p δ_b^q - Λ_{mn}^{ab}T^i_a δ_i^p δ_b^q - Λ_{mn}^{ab}δ_a^p T^j_b δ_j^q + Λ_{mn}^{ab}T^i_a T^j_b δ_i^p δ_j^q$.

calc_doubles_energy(EC::ECInfo, T2a, T2b, T2ab)

Calculate energy for αα (T2a), ββ (T2b) and αβ (T2ab) doubles amplitudes. Returns total energy, SS, OS and Openshell contributions as OutDict with keys (E,ESS,EOS,EO).

calc_doubles_energy(EC::ECInfo, T2, T2ep)

Calculate energy for αα (T2), and αp (T2ep) doubles amplitudes. Returns total energy, SS, OS and Openshell contributions as OutDict with keys (E,ESS,EOS,EO).

calc_doubles_energy(EC::ECInfo, T2)

Calculate coupled-cluster closed-shell doubles energy. Returns total energy, SS, OS and Openshell (0.0) contributions as OutDict with keys (E,ESS,EOS,EO).

calc_dressed_ints(EC::ECInfo, T1;
          calc_d_vvvv=EC.options.cc.calc_d_vvvv, calc_d_vvvo=EC.options.cc.calc_d_vvvo,
          calc_d_vovv=EC.options.cc.calc_d_vovv, calc_d_vvoo=EC.options.cc.calc_d_vvoo)

Dress integrals with singles.

$\hat v_{ab}^{cd}$, $\hat v_{ab}^{ci}$, $\hat v_{ak}^{cd}$ and $\hat v_{ab}^{ij}$ are only calculated if requested in EC.options.cc or using keyword-arguments.

calc_dressed_ints(EC::ECInfo, T1, T12, o1::Char, v1::Char, o2::Char, v2::Char;
          calc_d_vvvv=EC.options.cc.calc_d_vvvv, calc_d_vvvo=EC.options.cc.calc_d_vvvo,
          calc_d_vovv=EC.options.cc.calc_d_vovv, calc_d_vvoo=EC.options.cc.calc_d_vvoo)

Dress integrals with singles amplitudes.

The singles and orbspaces for first and second electron are T1, o1, v1 and T12, o2, v2, respectively. The integrals from EC.fd are used and dressed integrals are stored as d_????. $\hat v_{ab}^{cd}$, $\hat v_{ab}^{ci}$, $\hat v_{ak}^{cd}$ and $\hat v_{ab}^{ij}$ are only calculated if requested in EC.options.cc or using keyword-arguments.

calc_focklike_vT2(EC::ECInfo, T2a::AbstractArray, T2b::AbstractArray, T2ab::AbstractArray)

Calculate the fock-like intermediates $v_{kl}^{cd} T^{il}_{cd}$ and v_{kl}^{cd} T^{kl}_{ad} required in calc_ccsd_vector_times_Jacobian for αα/ββ, and αβ amplitudes. Store as vT_oo[ki], vT_vv[ac], vT_OO[ki] and vT_VV[ac].

calc_hylleraas(EC::ECInfo, T1a, T1b, T2a, T2b, T2ab, R1a, R1b, R2a, R2b, R2ab)

Calculate singles and doubles Hylleraas energy. Returns total energy, SS, OS and Openshell contributions as OutDict with keys (pE,pESS,pEOS,pEO,E,ESS,EOS,EO) where pE are projected energies.

calc_hylleraas(EC::ECInfo, T1, T2, R1, R2)

Calculate closed-shell singles and doubles Hylleraas energy. Returns total energy, SS, OS and Openshell (0.0) contributions as OutDict with keys (pE,pESS,pEOS,pEO,E,ESS,EOS,EO) where pE are projected energies.

calc_hylleraas4spincase(EC::ECInfo, o1, v1, o2, v2, T1, T1OS, T2, R1, R2, fov)

Calculate singles and doubles Hylleraas energy for one spin case.

Returns OutDict with keys (pE2,pE1,pE1_2,E2,E1,E1_2) where pE are projected energies and E are Hylleraas energies, and 2, 1, 1_2 are doubles, singles and quadratic singles contributions.

calc_intermediates4triples(EC::ECInfo)

Calculate intermediates for decomposed triples independent of the amplitudes.

calc_lagrange_K2(int2, D2; symmetrize=true)

Calculate the kext K2 contribution to the CCSD Lagrange residuals by directly contracting the integrals with D2.

$K_{ij}^{pq} = v^{pq}_{rs} D_ij^rs$

Return Kmmoo::Array{4}.

calc_lagrange_K2ab(int2, D2)

Calculate the kext K2 contribution to the UCCSD Lagrange residuals for the αβ part by directly contracting the integrals with D2.

$K_{mN}^{rS} = v^{rS}_{pQ} D_{mN}^{pQ}$

Return KmMoO::Array{4}.

calc_lagrange_Koox(int2, D2)

Calculate 3-index intermediate in the kext in CCSD Lagrange residuals.

Return Koox::Array{3}, where $Koox_{mn}^x = v^{x}_{pq} D_{mn}^{pq}$.

calc_oqv_gradient(EC::ECInfo, q1T, q2T, R)

Calculate the gradient for orbital optimization in OQV-CCD/DCD. q1T and q2T are the transformed amplitudes. R is the rotation matrix. Return the gradient as a matrix.

calc_pertT_closed_shell(EC::ECInfo; save_t3=false)

Calculate (T) correction for closed-shell CCSD. If the Fock matrix is not diagonal in the occ-occ and virt-virt blocks, performs pseudo-canonicalization before the (T) calculation.

Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict.

calc_pertT_mixedspin(EC::ECInfo, T1, T2, T1os, T2mix, pct, spin::Symbol)

Calculate mixed-spin (T) correction for UCCSD(T) (i.e., ααβ or ββα).

spin ∈ (:α,:β) T1 and T2 are same-spin amplitudes, T1os are opposite-spin amplitudes, and T2mix are mixed-spin amplitudes with the second electron being spin, i.e., Tβα for spin == :α and Tαβ for spin == :β. Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict.

calc_pertT_samespin(EC::ECInfo, T1, T2, pct, spin::Symbol)

Calculate same-spin (T) correction for UCCSD(T) (i.e., ααα or βββ). spin ∈ (:α,:β)

Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict.

calc_pertT_unrestricted(EC::ECInfo)

Calculate (T) correction for UCCSD.

Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict.

calc_qG_cc(EC::ECInfo, qV, T2, T2t, Ae, AX, Be, BX, Ce, CX, Ye, YX, We, WX, AU1, BU1, CU1, Y1, W1, q)

Calculate QV-DCD closed-shell residuals qG for CC. qV is the specially defined integrals in QV-CCD. T2 and T2t are the transformed doubles amplitudes. Ae, AX, Be, BX, Ce, CX, Ye, YX, We, WX are the eigenvalues and eigenvectors of the corresponding U in QV-CCD. AU1, BU1, CU1, Y1, W1 are the AU^(-q/2), BU^(-q/2), CU^(-q/2), Y^(-q/2) and W^(-q/2) matrices. Return qG as a matrix.

calc_qG_dc(EC::ECInfo, qV, T2, T2t, qVD, Ae, AX, Be, BX, Ye, YX, AU1, BU1, Y1, q)

Calculate QV-DCD closed-shell residuals qG. qV and qVD are the specially defined integrals in QV-DCD. T2 and T2t are the transformed doubles amplitudes. Ae, AX, Be, BX, Ye, YX are the eigenvalues and eigenvectors of the corresponding U in QV-DCD. AU1, BU1, Y1 are the AU^(-q/2), BU^(-q/2) and Y^(-q/2)matrices. Return qG as a matrix.

calc_qT_cc(AU1::Matrix, BU1::Matrix, CU1::Array, Y1::Array, W1::Array, T2, T2t)

Calculate QV-DCD closed-shell transformed amplitudes qT.

calc_qT_dc(AU1::Matrix, BU1::Matrix, Y1::Array, T2, T2t)

Calculate QV-DCD closed-shell transformed amplitudes qT.

calc_qvcc_resid(EC::ECInfo, T1, T2; dc=false, orbopt=false, brueckner=false)

Calculate QV-CCD or QV-DCD closed-shell residual.

calc_rings_vT2(EC::ECInfo, T2a::AbstractArray, T2b::AbstractArray, T2ab::AbstractArray; dc=false)

Calculate the ring intermediates required in calc_ccsd_vector_times_Jacobian $\hat y_{am}^{ie} = \hat v_{am}^{ie} - \hat v_{am}^{ei} + 2x_{am}^{ie} + v_{mL}^{eD} T^{iL}_{aD}$ and $\hat y_{Bn}^{Jf} = \hat v_{nB}^{fJ} + 2x_{Bn}^{Jf} + v_{nL}^{fD} T^{LJ}_{DB}$ and the spin-flip version of them, with $2x_{am}^{ie} = T^{il}_{ad} (v_{lm}^{de} \red{- v_{ml}^{de}})$ and $2x_{Am}^{Ie} = T^{Il}_{Ad} (v_{lm}^{de} \red{- v_{ml}^{de}})$

The intermediates are stored as vT_voov[amie], vT_VOOV[AMIE], vT_VoOv[BnJf], vT_vOoV[bNjF],

calc_rotated_fock(EC::ECInfo, into, R::Matrix)

Calculate the Fock matrix in the rotated orbital basis with the original into is the mmmo integral. R is the rotation matrix. Return the Fock matrix.

calc_singles_energy(EC::ECInfo, T1a, T1b; fock_only=false)

Calculate energy for α (T1a) and β (T1b) singles amplitudes. Returns total energy, SS, OS and Openshell contributions as OutDict with keys (E,ESS,EOS,EO).

calc_singles_energy(EC::ECInfo, T1; fock_only=false)

Calculate coupled-cluster closed-shell singles energy. Returns total energy, SS, OS and Openshell (0.0) contributions as OutDict with keys (E,ESS,EOS,EO).

calc_space4project_voXL(EC::ECInfo, T2)

Calculate space for project_voXL=true approximation.

It is a combination of spaces for triples and contravariant doubles.

calc_triples_decomposition(EC::ECInfo)

Decompose $T^{ijk}_{abc}$ as $U^{iX}_a U^{jY}_b U^{kZ}_c T_{XYZ}$.

calc_triples_decomposition_without_triples(EC::ECInfo, T2)

Decompose $T^{ijk}_{abc}$ as $U^{iX}_a U^{jY}_b U^{kZ}_c T_{XYZ}$ without explicit calculation of $T^{ijk}_{abc}$.

Compute perturbative $T^i_{aXY}$ and decompose $D^{ij}_{ab} = (T^i_{aXY} T^j_{bXY})$ to get $U^{iX}_a$.

calc_triples_residuals!(EC::ECInfo, R1, R2, T2)

Calculate decomposed triples DC-CCSDT residuals.

calc_vT2_intermediates(EC::ECInfo, T2a, T2b, T2ab; dc=false)

Calculate intermediates required in unrestricted calc_ccsd_vector_times_Jacobian

calc_vT2_intermediates(EC::ECInfo, T2; dc=false)

Calculate intermediates required in closed-shell calc_ccsd_vector_times_Jacobian

calc_ΛpertT_closed_shell(EC::ECInfo)

Calculate (T) correction for closed-shell ΛCCSD(T).

The amplitudes are stored in T_vvoo file, and the Lagrangian multipliers are stored in U_vvoo file. If the Fock matrix is not diagonal in the occ-occ and virt-virt blocks, performs pseudo-canonicalization before the (T) calculation. Return ( "ET3"=(T) energy, "ET3b"=[T] energy) OutDict.

calc_ΛpertT_mixedspin(EC::ECInfo, T2, T2mix, U1, U2, U1os, U2mix, pct, spin::Symbol)

Calculate mixed-spin (T) correction for ΛUCCSD(T) (i.e., ααβ or ββα).

spin ∈ (:α,:β) U1 and U2/T2 are same-spin Lagrange multipliers/amplitudes, U1os are opposite-spin Lagrange multipliers, and U2mix/T2mix are mixed-spin Lagrange multipliers/amplitudes with the second electron being spin, i.e., Tβα for spin == :α and Tαβ for spin == :β. Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict.

calc_ΛpertT_samespin(EC::ECInfo, T2, U1, U2, pct, spin::Symbol)

Calculate same-spin (T) correction for ΛUCCSD(T) (i.e., ααα or βββ). spin ∈ (:α,:β)

Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict.

calc_ΛpertT_unrestricted(EC::ECInfo)

Calculate (T) correction for ΛUCCSD(T).

Return ( "ET3"=(T)-energy, "ET3b"=[T]-energy)) OutDict.

cc_kext!(EC::ECInfo, R1a, R1b, R2a, R2b, R2ab, T1a, T1b, T2a, T2b, T2ab, Rota, Rotb)

Calculate the 4-external (and some more) contribution to the UCCSD residuals.

R1* and R2* are the singles and doubles residuals to be updated. T1* and T2* are the singles and doubles amplitudes. Rot* are the rotation matrices for orbital optimization (if any).

cc_kext!(EC::ECInfo, R1, R2, T1, T2, Rot)

Calculate the 4-external (and some more) contribution to the CCSD residuals.

R1 and R2 are the singles and doubles residuals to be updated. T1 and T2 are the singles and doubles amplitudes. Rot is the rotation matrix for orbital optimization (if any).

cc_lagrange_kext!(EC::ECInfo, R1, R2, T1, U2, spin=:closed)

Calculate the contribution of the 4-external integrals to the Λ equations for CCSD/DCSD.

spin can be :closed, :α, or :β.

cc_lagrange_kext4ab!(EC::ECInfo, R1a, R1b, R2, T1a, T1b, U2)

Calculate the contribution of the 4-external integrals to the Λ αβ equations for CCSD/DCSD.

check_fock_diagonal_blocks(EC::ECInfo, spin::Symbol=:α)

Check if the Fock matrix occ-occ and virt-virt blocks are diagonal.

For the check the Fock matrix block is symmetrized (to account for non-Hermitian cases with already diagonalized blocks). Returns (is_diagonal_oo, is_diagonal_vv, F_oo, F_vv).

compute_pseudocanonical_transform(F_block::Matrix; skip::Bool=false, hermitian::Bool=true)

Diagonalize a potentially non-Hermitian Fock block. Returns (ϵ, Cr, Cl): eigenvalues, right/left eigenvector matrices.

If skip=true, skips the diagonalization and returns identity transforms. If hermitian=true, assumes the Fock block is Hermitian and uses eigen(Hermitian(...)).

Uses rotate_eigenvectors_to_real for complex pairs.

dD1_fock_vo(EC::ECInfo, dD1::AbstractMatrix) -> Matrix

The virtual-occupied block the closed-shell Λ residual needs from the general-orbital density dD1: $R^e_m += Σ_pq dD1_p^q (2⟨qm|pe⟩ − ⟨qm|ep⟩)$ — structurally the v,o block of a generalized (2J−K) Fock built with dD1. AO-direct: rotate D_AO = cMO·dD1ᵀ·cMOᵀ and read the [e,m] block straight off the half-transformed store (dD1_ht_vo with the Coulomb density 2·D_AO and the exchange density D_AO), so no nao×nao Fock matrix is formed and no general-orbital nocc·norb³ block (the ints2(EC,"momm") read of the MO path) either. Returns an [nvirt,nocc] block.

dD1_ht_vo(EC, htkey, DJ, DK, Cv) -> Matrix

The [e,m] generalized-Fock block the Λ residual needs, read off the half-transformed store htkey (built once per orbital set in ao_cc_setup!) instead of building a full nao×nao Fock:

  R^e_m = Σ_pq ( dD1ᴶ_p^q ⟨qm|pe⟩ − dD1ᴷ_p^q ⟨qm|ep⟩ )

with the two general-orbital densities supplied already rotated to the AO basis, D[x,y] = Σ_pq cMO[x,q] dD1_p^q cMO[y,p] (i.e. cMO·dD1ᵀ·cMOᵀ, the convention of the callers).

The store's bra index IS the (active) occupied m this term carries, and its two roles hold exactly the two orderings the term needs — with the column ν in the ket-2 slot both times:

  B_ν[m,x,y] = ⟨x m|y ν⟩    (B-role, occ on bra-2)  → Coulomb, contracted with `DJ[x,y]`
  A_ν[m,x,y] = ⟨m x|y ν⟩    (A-role, occ on bra-1)  → exchange, contracted with `DK[x,y]`

(the exchange uses the conjugation-free particle symmetry ⟨mq|pe⟩ = ⟨qm|ep⟩, which is what puts its virtual on the same ket-2 slot as the Coulomb's, so it too is a plain contraction of the slab's two free AO slots with an AO density instead of a partial one). Both are then done for every ν at once by ht_jk_columns! — one sequential pass over the store, each stored element read once — and the resulting t[m,ν] only has to have its ket-2 slot transformed to the virtual e. Cost: 2·nocc·nao³ MACs and nocc·nao³ elements read, against the nao⁴/2 MACs and nao⁴/4 streamed integrals of a full 2J−K build of which only the [v,o] block was kept.

Element-type generic: the store applies plain C (the AO-direct/FCIDUMP detri convention) and the contraction uses no conjugation anywhere, so unlike the ao_JK! route this is complex-correct.

dD1_ufock_vo(EC::ECInfo, dD1, dD1os, spin::Symbol) -> Matrix

Unrestricted analogue of dD1_fock_vo: the [e,m] block the open-shell Λ residual adds for spin from its own general-orbital density dD1 and the opposite-spin dD1os. The three general-orbital reads it replaces — same-spin ⟨qm|pe⟩−⟨qm|ep⟩ (a momm/MOMM block) plus the opposite-spin Coulomb ⟨Qm|Pe⟩ (oMvM/mOmV) — are the virtual-occupied block of the UHF generalized Fock F^σ = J(D^α+D^β) − K(D^σ), i.e. dD1_ht_vo with the Coulomb density D^α_AO + D^β_AO and the exchange density D^σ_AO. AO integrals are spin-free, so the opposite-spin Coulomb rides along on the SAME (spin) store — one pass over "ht_oAAA_"*spin, no nao×nao Fock matrices and no nocc·norb³ general-orbital blocks.

dress_fock_closedshell(EC::ECInfo, T1)

Dress the fock matrix (closed-shell). The dressed fock matrix is stored as df_mm.

dress_fock_oppositespin(EC::ECInfo)

Add the dressed opposite-spin part to the dressed Fock matrix.

dress_fock_samespin(EC::ECInfo, T1, o1::Char, v1::Char)

Dress the fock matrix (same-spin part).

Unrestricted dress_lambda_ints!: AO-direct dressing (with the 3-external blocks) or the dump. Empty T1a/T1b ⇒ bare blocks (Λ-UCCD/Λ-UDCD), AO-direct only.

dress_lambda_ints!(EC, T1)

Closed-shell integral dressing for the Λ residual / EOM Jacobian (dressed d_* blocks incl. the 3-external d_vovv). AO-direct: ao_dressed_ints with calc_d_vovv=true (built from the half-transformed store) instead of the MO-fcidump calc_dressed_ints. An EMPTY T1 (methods without singles, Λ-CCD/Λ-DCD) leaves the coefficients undressed and hence writes the BARE blocks under the same d_* names — the AO-direct analogue of pseudo_dressed_ints. Only supported AO-direct (the MO dressing has no empty-T1 path).

get_pseudo_orbital_energies(pct::PseudoCanonicalTransform, spin::Symbol=:α)

Get pseudo-canonical orbital energies from the transform object. Returns (ϵo, ϵv) for the specified spin.

ht_mo_block(EC, htkey, CX, CY, CZ) -> W

General half-transformed-store → MO-block kernel: build a 4-index MO integral block that carries the store's (occupied) bra index i plus three MO indices transformed by the coefficient matrices CX/CY/CZ, in ONE sweep over the store htkey (Σ_μ⟨μν|ρσ⟩C[μ,i], read column-by-column with ht_column_A!):

  W[i,X,Y,Z] = Σ_μνρσ ⟨μν|ρσ⟩ C[μ,i] CX[ν,X] CY[ρ,Y] CZ[σ,Z]   ( = ⟨iX|YZ⟩ )

Every block a CC method needs (they all have ≥1 occupied index) is one call plus an output permutation — see ht_mo_block_spec. Per σ-column two GEMMs transform the free bra ν and ket-1 ρ into the [(i,X),Y] intermediate; stacking those over σ and one final GEMM contracts ket-2 σ into Z. Generic over the element type (the store applies plain C, matching the AO-direct/FCIDUMP detri convention). Cost: one pass over the store + GEMMs; peak RAM ≈ the output.

There is no second "occupied on bra-2" kernel: by the particle-exchange symmetry ⟨μν|ρσ⟩=⟨νμ|σρ⟩,

  ⟨Xi|YZ⟩ = ⟨iX|ZY⟩   ⟹   B(CX,CY,CZ) = permutedims(ht_mo_block(CX,CZ,CY), (2,1,4,3))

i.e. such a block is THIS kernel with the two ket coefficients exchanged. That identity is exact for complex integrals and for a non-Hermitian (similarity-transformed, CL≠CR) transformation, because it never exchanges bra with ket — unlike the swap relation in ht_mo_block_spec.

ht_mo_block_spec(name) -> (store, X, Y, Z, perm, swap)

Resolve a block name (a 4-character space string ⟨s₁s₂|s₃s₄⟩, lowercase = α, uppercase = β) into the ht_mo_block call that produces it: which per-spin store supplies the block's occupied index, which space goes on each free slot, and how to permute the result.

The store holds its occupied index on a BRA slot, so where the block's occupied index sits decides everything. With A ≡ ⟨iX|YZ⟩ and using ⟨pq|rs⟩ = ⟨qp|sr⟩ (particle exchange — always valid) and ⟨pq|rs⟩ = ⟨rs|pq⟩ (bra↔ket — swap, see below):

| occ at | rewrite | A-form | perm | swap | |:–––:|:–––––––––––––-|:––––––|:––––––|:––:| | 1 | ⟨is₂\|s₃s₄⟩ | (s₂,s₃,s₄)| (1,2,3,4) | no | | 2 | ⟨s₁i\|s₃s₄⟩ = ⟨is₁\|s₄s₃⟩ | (s₁,s₄,s₃)| (2,1,4,3) | no | | 3 | ⟨s₁s₂\|is₄⟩ = ⟨is₄\|s₁s₂⟩ | (s₄,s₁,s₂)| (3,4,1,2) | YES | | 4 | ⟨s₁s₂\|s₃i⟩ = ⟨is₃\|s₂s₁⟩ | (s₃,s₂,s₁)| (4,3,2,1) | YES |

The first occupied slot is taken, so a bra one wins whenever the block has one — the point being that rows 1–2 need only particle exchange, which holds for complex and for non-Hermitian (CL≠CR) transformations. Rows 3–4 have the occupied index on a KET slot, which no particle exchange can move to the bra, so they need the bra↔ket relation: valid only for Hermitian integrals with CL=CR and real coefficients (guarded in save_mo_blocks!; complex would need a ket-transformed store).

Blocks sharing an A-form differ only in perm/swap and so share one sweep — vovv/vvvo, ovoo/vooo, VOVV/VVVO, vOvV/vVvO, oVvV/vVoV.

ints2_t(EC::ECInfo, spaces::AbstractString)

A bare MO integral block for the (T) kernels. AO-direct reads the block prebuilt from the half-transformed store — the scratch files are named by their space string, see build_ht_mo_blocks! / build_ht_mo_blocks_unrestricted! — otherwise the block is cut from the MO fcidump by ints2.

ints2_t_oovv(EC::ECInfo, spaces::AbstractString)

As ints2_t for the oovv-class blocks (oovv/OOVV/oOvV), which the AO-direct dressing already stores as bare blocks — read via load_bare_int2 instead of the block engine.

kext_rs_blocksize(nrs, ncols; scratch_per_rs=0, nout=0)

Block length for the rs loop of the 4-external (kext) contractions (calc_K2).

Each block is a single GEMM whose contraction (k) dimension is the block length and which read-modify-writes the whole result, so the shared get_spaceblocks default of 128 streams the result nrs/128 times (212 times at norb=232) on top of the one unavoidable pass over the integrals. The ratio of the two is 2*ncols/blocklength, i.e. it is governed by ncols, the width of the GEMM's right-hand side (nocc1*nocc2, or ntri_oo for the ± variant) — and so is the measured payoff: for narrow results the result traffic does not exceed the integral stream and enlarging the block loses (up to 17% at ncols=64), while for wide ones the gain reaches 1.35× (at ncols=400, i.e. 20 correlated occupied orbitals, norb=232).

The rule is deliberately conservative (see KEXT_RS_RAMP_FROM): the blocking is left exactly as it was unless ncols exceeds that threshold, and only then grows by KEXT_RS_PER_COL per result column. So every shape either keeps its previous, bit-for-bit identical blocking or gets a block in a range where enlarging was measured to win — the non-monotonic band in between is never selected. The payoff is not governed by ncols alone (the same ncols=144 loses 6% at norb=164, where the 128-block GEMM is already at roofline, and gains 1.08× at norb=232), which is the reason for the threshold rather than a finer fit.

scratch_per_rs (elements of per-rs scratch, if any) and nout (elements of the result the scratch is measured against) cap the block so the scratch stays below a quarter of a result array that is allocated anyway — no new memory regime, at most a few percent on the peak. The result never exceeds nrs.

lm_cc_iterations!(LMs1, LMs2, EC::ECInfo, method::ECMethod)

Perform the CCSD or UCCSD iterations using Lagrange multipliers.

LMs1 are the Lagrange multipliers for the singles equations, LMs2 are the Lagrange multipliers for the doubles equations.

load_bare_int2(EC::ECInfo, name::AbstractString) -> Array

Return the bare (undressed) MO 2-e integral block name ("oovv", "OOVV", "oOvV", …), cached in the scratch file "d_"*name: extracted from the integral source (ints2) on first use and loaded thereafter, so it is materialized at most once per calculation instead of on every use. The AO-direct path has "d_"*name prebuilt by ao_cc_setup!/the AO dressing (EC.fd is empty there), so it is simply loaded. The cache is a temporary file (rebuilt each calculation), so it never goes stale. This unifies the FCIDUMP/DF and AO-direct paths — no EC.ao_direct branch.

oss_active_orbitals(EC::ECInfo)

Return the four active orbitals of an (2e,2o) open-shell singlet problem based on a single determinant reference.

pm_K2!(pm, D2, tripp)

kext K2 from the persisted ± supermatrix store — the amortized replacement for the per iteration. Reuses the same ij/rs ±-fold of the density and 4-quadrant output scatter, but obtains the ± integral action as zero-copy panel GEMMs $s\!K2 = V_s·D_s$ / $a\!K2 = V_a·D_a$ (pm_matmul!) over the stored lower block-triangle — halved flops and streaming, no per-iteration ± build. D2[tri(pq),i,j] must already carry the ½ rs-diagonal (calc_D2(...; scalepp=true), as cc_kext! passes).

The ±-fold of the density (pm_fold_ij!) and the 4-quadrant unpacking of the products (pm_scatter_K2!) are pure memory traffic around the two GEMMs; both run as one fused multi-threaded pass over the i ≤ j pairs instead of the CartesianIndex-cut gather/scatter broadcasts.

$K^{ij}_{pq} = v_{pq}^{rs} D^{ij}_{rs}$

Return K2pq::Array{4}.

pm_K2ab!(pm, D2ab_full, tripp)

αβ kext from the persisted ± store: the full contraction $K2ab_{pq}^{iJ} = Σ_{rs} ⟨pq|rs⟩ D^{iJ}_{rs}$ for the (non-rs-symmetric) αβ density. The rs-± fold is explicit — Ds/Da = ½(D[pq] ± D[qp]), with the ½ rs-diagonal already carried by calc_D2ab(...; scalepp=true) — then two pm_matmul! panel GEMMs and the ± unscatter to both pq orders (pm_scatter_K2ab!). Halved flops and streaming vs the two joint-store calc_K2 passes.

Because this fold transposes the AO pair (unlike the occupied-pair fold of pm_K2!), it is exactly calc_tri_sym_antisym! with fac = ½ — one fused, threaded, row-buffered pass instead of two CartesianIndex gathers plus two broadcasts over a pair of npp × nanb temporaries.

pm_fold_ij!(Ds, Da, D2, trioo)

±-fold of the kext density over the occupied pair: Ds/Da[rs,ij] = ½(D2[rs,i,j] ± D2[rs,j,i]) for the i ≤ j list trioo, in a single fused pass (Ds, Da are npp × length(trioo)).

Both source columns and both destination columns are contiguous in rs, so this is four streaming vectors; the previous D2[:,trioo] ± D2[:,trioo_swap] form instead paid two CartesianIndex gathers plus two broadcasts over a pair of npp × ntri_oo temporaries.

Threaded over ij: iteration ij writes only column ij of Ds/Da, and each unordered pair {i,j} occurs once, so the writes are disjoint (i == j merely reads one column twice and yields Da[:,ij] = 0, as the ± fold requires).

pm_scatter_K2!(K2pq, sK2, aK2, trioo)

Unpack the ± kext products into K2pq[p,q,i,j], both pq orders and both ij orders: K2pq[p,q,i,j] = K2pq[q,p,j,i] = sK2[pq,ij] + aK2[pq,ij] and K2pq[p,q,j,i] = K2pq[q,p,i,j] = sK2[pq,ij] - aK2[pq,ij] (p ≤ q, ij running over trioo).

One fused pass over the output replaces the four K2pq[cut,cut] .= sK2 .± aK2 broadcasts, which re-read both products (and re-evaluated sK2 .± aK2) four times over.

The pq- and ij-diagonals are where the four quadrants alias, and they are written once, explicitly: on p == q the two pq orders coincide, on i == j the two ij orders coincide. Both carry aK2 = 0 exactly (V_a has zero pp rows; D_a has zero ii columns), so the single value s - a written there equals s + a, matching the previous last-write-wins. With no two writes of an iteration hitting one address, @simd ivdep is sound and the Threads.@threads over ij is race-free: iteration ij owns exactly the (i,j) and (j,i) planes of K2pq, and distinct ij are distinct unordered pairs, hence disjoint planes.

pm_scatter_K2ab!(K2, sK, aK)

Unpack the αβ ± kext products into both pq orders of the (pair-nonsymmetric) αβ K2: K2[p,q,iJ] = sK[pq,iJ] + aK[pq,iJ], K2[q,p,iJ] = sK[pq,iJ] - aK[pq,iJ] for p ≤ q. K2 is the output reshaped to [p, q, iJ] with the αβ pair flattened.

One fused pass over the output replaces the two K2[cut,:,:] .= sK .± aK broadcasts. The pq-diagonal, where the two orders alias, is written once (there aK = 0 exactly — V_a has zero pp rows — so s - a equals s + a, matching the previous last-write-wins). Threaded over iJ, which owns one full norb × norb plane, hence disjoint writes.

pseudo_dressed_ints(EC::ECInfo, unrestricted=false;
          calc_d_vvvv=EC.options.cc.calc_d_vvvv, calc_d_vvvo=EC.options.cc.calc_d_vvvo,
          calc_d_vovv=EC.options.cc.calc_d_vovv, calc_d_vvoo=EC.options.cc.calc_d_vvoo)

Save non-dressed integrals in files instead of dressed integrals.

pseudocan_transform!(pct::PseudoCanonicalTransform, arr::AbstractArray, indices::String; 
                  conjugate::Bool=false)

Transform an array to the pseudo-canonical basis in-place using the transformation matrices.

The indices string specifies the index type for each dimension:

• 'o'/'O': occupied α/β
• 'v'/'V': virtual α/β

The transformation convention (canonical order):

• First half of indices = lower indices (subscript/bra) → Left eigenvectors
• Second half of indices = upper indices (superscript/ket) → Right eigenvectors

For conjugate=true, the L/R roles are swapped (used for Lagrange multipliers).

Supported array dimensions: 2 and 4.

Example

# Transform T2 amplitudes stored as T_{ab}^{ij} in [a,b,i,j] order
pseudocan_transform!(pct, T2, "vvoo")

# Transform U2 Lagrange multipliers (use conjugate)
pseudocan_transform!(pct, U2, "vvoo"; conjugate=true)

# Transform mixed-spin integrals v_{aB}^{iJ}
pseudocan_transform!(pct, int_aB_iJ, "vVoO")

rotate_ints(EC::ECInfo, R::Matrix)

Update the fock matrix with rotated integrals. Rotate the orginal integrals ith the rotation matrix R. This function calculates various integrals and saves them in the ECInfo object (disk).

rotate_ints_o(EC::ECInfo, R::Matrix)

Calculate into (integral – 1 occupied) mmmo integral from the original molecular orbitals. Return as a tensor into[p',q',r',i] where p', q', r' are original molecular orbitals and i is the occupied orbital.

save_ao_AAAo!(EC::ECInfo)

Build the general-orbital block the OQV orbital gradient needs, kept in the AO space: d_AAAo[μ,ν,ρ,i] = ⟨μν|ρi⟩. Read off the half-transformed store "ht_oAAA", whose bra IS the occupied space, so the caller (ao_rotate_ints) must have rebuilt it for the current rotation.

The key is d_AAAo, not the d_mmmo of the MO route (rotate_ints_o): the two hold the same integrals over DIFFERENT one-particle bases (three AO indices here, three MO indices there), and calc_oqv_gradient is indifferent only because it contracts all three with whatever basis its rotation maps from. Separate names keep a block built in one basis from being read against the other.

The occupied index sits on a ket, which no particle exchange can move to the bra, so this needs the bra↔ket relation ⟨μν|ρi⟩ = conj⟨ρi|μν⟩ and is real-only like the other such blocks. It reads the B-role slab (ht_column_B!) — not because the A slab could not produce the block, but because the fixed [μ,ν,ρ,i] output makes [:,ν,:,:] the contiguous write.

save_ao_correlation_orbitals!(EC::ECInfo, cMO::SpinMatrix)

Persist the AO-direct correlation reference computed by ao_cc_setup! as C_Am (and C_AM for the β spin) — the single artifact every later consumer of THIS run reads, the AO-direct counterpart of the MO route's FCIDUMP. It has to be stored rather than re-derived because it is not a pure function of the orbital file: with wf.dump4core_only the frozen core is spliced in from a second file and the correlating orbitals are re-orthonormalized against it (replace_core_from_dump!). Restricted orbitals write C_Am only.

The files are temporary ("tmp"): they belong to the run, so delete_temporary_files! reclaims them at its end — like the per-run DF fd, the reference does not outlive the calculation that settled it.

Closed-shell save_mo_blocks!: the single store htkey with occupied/virtual coefficients.

save_mo_blocks!(EC, names, htkeys, coefs)

Build the bare MO blocks names and save them (mmapped, "tmp" so they are reclaimed at end-of-run) under their plain space names, in the index order the consumers read.

Blocks that resolve to the same ht_mo_block_spec A-form share ONE store sweep and differ only by the output permutation — e.g. Λ(T)'s 16 unrestricted blocks come from 12 sweeps, including a single sweep for each of the four expensive 3-external pairs.

test_UaiX(EC, UaiX)

Test UaiX

test_add_to_singles_and_doubles_residuals(R1, R2, T1, T2)

Test R1(T3) and R2(T3)

test_calc_pertT_from_T3(EC, T3)

Test [T]

test_dressed_ints(EC, T1)

Compare 3-idx dressed integrals to 4-idx dressed integrals.

transform_2idx!(arr, U1, U2)

Transform a 2-index array in place; a nothing leaves that index alone (see xform_idx!).

transform_4idx!(arr::AbstractArray{T,4}, U1, U2, U3, U4) where T

Transform a 4-index array in place: arr[p',q',r',s'] = U1[p,p'] * U2[q,q'] * U3[r,r'] * U4[s,s'] * arr[p,q,r,s].

A nothing in place of a matrix leaves that index alone — used when an index is already in the target basis (e.g. the AO-direct blocks, whose virtual indices are built from rotated coefficients), so only the remaining indices are contracted. One scratch buffer is allocated and the two are alternated, so at most one copy back into arr is needed.

warn_ao_direct_deleted(EC::ECInfo, norb)

Warn when a large fraction of the norb orbitals is deleted — whether linearly dependent or removed on purpose (e.g. by @region); both are equally costly here. AO-direct still streams the full AO dimension, so beyond AO_DIRECT_DELETED_WARN the derived-MO-dump route (which correlates in the reduced basis) is typically faster. Points at the option that switches routes.

2-index xform_idx!.

xform_idx!(dst, src, U, ::Val{k}) -> (result, spare)

Contract index k of src with U, writing into dst, and return (result, spare) so the caller can keep alternating between two buffers. Dispatch on U:

• U::Matrix — result is dst, src becomes the spare buffer;
• U::Nothing — that index is already in the target basis (e.g. the AO-direct blocks, whose virtuals are built from rotated coefficients): nothing is contracted, the result stays in src and dst stays spare.

Splitting the no-op out as its own method keeps the caller free of nothing branches and lets each contraction specialize.