Coupled-cluster tools

CCTools

A collection of tools for working with coupled cluster theory.

calc_HF_energy(EC::ECInfo, closed_shell)

Calculate HF energy from FCIDump and EC info.

calc_ab_doubles_dot(T2, T2_)

Calculate dot of unrestricted αβ doubles amplitudes.

calc_contra_cs_doubles_dot(T2, T2_)

Calculate dot of contravariant closed-shell doubles amplitudes.

calc_contra_cs_singles_dot(T1, T1_, state=0)

Calculate dot product of contravariant closed-shell singles amplitudes.

calc_contra_doubles_norm(T2a, T2b, T2ab)

Calculate squared norm of unrestricted doubles amplitudes (the same as calc_doubles_norm)

calc_contra_doubles_norm(T2)

Calculate squared norm of closed-shell contravariant doubles amplitudes.

calc_contra_singles_norm(T1a, T1b)

Calculate squared norm of unrestricted singles amplitudes (same as calc_singles_norm(T1a, T1b)).

calc_contra_singles_norm(T1)

Calculate squared norm of closed-shell contravariant singles amplitudes.

calc_cs_doubles_dot(T2, T2_)

Calculate dot of closed-shell doubles amplitudes.

calc_cs_singles_dot(T1, T1_, state=0)

Calculate dot product of closed-shell singles amplitudes.

calc_cs_triples_dot(T3, T3_)

Calculate dot of closed-shell triples amplitudes.

calc_deco_doubles_norm(T2, tT2=Float64[])

Calculate squared norm of doubles (for decomposed doubles: without contravariant!) T2 are decomposed doubles amplitudes T2[X,Y]=$T_{XY}$ or full doubles amplitudes T2[a,b,i,j]=$T^{ij}_{ab}$.

If the contravariant amplitude tT2 is provided, the norm will be calculated as $T_{XY} T̃_{XY}$.

calc_deco_triples_norm(T3)

Calculate a simple norm of triples (without contravariant!)

calc_doubles_norm(T2a, T2b, T2ab)

Calculate squared norm of unrestricted doubles amplitudes.

calc_doubles_norm(T2)

Calculate squared norm of closed-shell doubles amplitudes.

calc_fock_matrix(EC::ECInfo, closed_shell, print_out=true)

Calculate fock matrix from FCIDump

calc_mixedspin_triples_dot(T3, T3_)

Calculate dot of unrestricted mixed-spin triples amplitudes.

calc_rotated_HF_energy(EC::ECInfo)

Calculate the Hartree-Fock energy in the rotated orbital basis.

calc_samespin_doubles_dot(T2, T2_)

Calculate dot of unrestricted same-spin doubles amplitudes.

calc_samespin_triples_dot(T3, T3_)

Calculate dot of unrestricted same-spin triples amplitudes.

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

Calculate coupled-cluster closed-shell singles energy using dressed fock matrix.

if fock_only is true, the energy will be calculated using only non-dressed fock matrix. Returns total energy, SS, OS, and Openshell (0.0) contributions as OutDict with keys (E, ESS, EOS, EO).

calc_singles_norm(T1a, T1b)

Calculate squared norm of unrestricted singles amplitudes.

calc_singles_norm(T1)

Calculate squared norm of closed-shell singles amplitudes.

calc_triples_norm(T3aaa, T3bbb, T3abb, T3aab)

Calculate squared norm of unrestricted triples amplitudes.

calc_triples_norm(T3)

Calculate squared norm of triples amplitudes.

calc_u_singles_dot(T1, T1_, state=0)

Calculate dot of unrestricted singles amplitudes.

check_projection_rank(P, expected_rank::Int, space_name::String)

Check if the projection matrix P has full rank for the expected dimension.

Returns (is_full_rank::Bool, actual_rank::Int). Prints a warning if the rank is less than expected.

clean_cs_triples!(T3)

Clean closed-shell triples amplitudes by setting $T^{iii}_{abc} = T^{ijk}_{aaa} = 0$.

contra2covariant(T2)

Transform contravariant doubles amplitudes to covariant.

dump_wavefunction_with_amplitudes!(EC::ECInfo, T1, T2; orbopt=false)

Dump orbitals and CC amplitudes to the TREXIO file specified in wf.store.

Does nothing if EC.options.wf.store is empty. For closed-shell case with singles T1 and doubles T2.

If orbopt=true, the orbitals are rotated using the T1 amplitudes via rotation_matrix and stored instead of the original orbitals. T1 amplitudes are not stored in that case.

dump_wavefunction_with_amplitudes!(EC::ECInfo, T1a, T1b, T2a, T2b, T2ab; orbopt=false)

Dump orbitals and unrestricted CC amplitudes to the TREXIO file.

Does nothing if EC.options.wf.store is empty.

If orbopt=true, the orbitals are rotated using the T1 amplitudes via rotation_matrix and stored instead of the original orbitals. T1 amplitudes are not stored in that case.

dump_wavefunction_with_determinants!(EC::ECInfo, dets, coeffs; nstates=0)

Dump orbitals and CIPHI determinants with CI coefficients to TREXIO file(s).

For single-state (nstates=0 or nstates=1), writes to wf.store. For multi-state, writes each state to a separate file per TREXIO standard:

• State 1: wf.store
• State n>1: wf.store with _stateN suffix

Does nothing if EC.options.wf.store is empty.

Arguments

• dets::Vector{<:AbstractDeterminant}: Determinants
• coeffs::AbstractVecOrMat{<:Number}: CI coefficients (vector for 1 state, matrix for multi-state)
• nstates::Int=0: Number of states (0 = infer from coeffs)

print_main_singles(U1, nelem; info="", thr=1e-4)

Utility function to print the main nelem singles.

project_amplitudes(EC::ECInfo, T1_old, T2_old, cMO_old::SpinMatrix, cMO_new::SpinMatrix, 
                   basis_old::BasisSet, basis_new::BasisSet; 
                   classes_old::Tuple{Vector{String},Vector{String}}=(String[], String[]),
                   classes_new::Tuple{Vector{String},Vector{String}}=(String[], String[]))

Project CC amplitudes from an old orbital basis to a new one.

The projection is performed as:

• Singles: $T_{a}^{i,\text{new}} = P_v^{a'} T_{a'}^{i',\text{old}} P_o^{i'}$
• Doubles: $T_{ab}^{ij,\text{new}} = P_v^{a'} P_v^{b'} T_{a'b'}^{i'j',\text{old}} P_o^{i'} P_o^{j'}$

where $P$ are projection matrices between old and new orbital spaces.

If classes_old is provided, it is used to determine which orbitals in the old basis were occupied ("Inactive"/"Active") vs virtual ("Virtual"). This is essential when core orbitals are frozen, as the MO coefficient matrix includes all orbitals but amplitudes only involve active orbitals.

If classes_new is provided, it is used similarly for the new basis. Otherwise, EC.space['o'] and EC.space['v'] are used (which may have frozen core indices renumbered from FCIDUMP).

Returns (T1_new, T2_new) for closed-shell case.

project_amplitudes(EC::ECInfo, T1a_old, T1b_old, T2a_old, T2b_old, T2ab_old,
                   cMO_old::SpinMatrix, cMO_new::SpinMatrix,
                   basis_old::BasisSet, basis_new::BasisSet;
                   classes_old=(String[], String[]), occupations_old=(Float64[], Float64[]))

Project unrestricted CC amplitudes from an old orbital basis to a new one.

Arguments

• classes_old: Tuple of (alphaclasses, betaclasses) orbital class vectors from old calculation. Used to identify Core/Deleted orbitals to exclude.
• occupations_old: Tuple of (alphaoccupations, betaoccupations) from old calculation. If provided, used to determine occupied/virtual orbitals (more reliable for UHF). If empty, falls back to classes_old or assumes contiguous occupied orbitals.

Returns (T1a_new, T1b_new, T2a_new, T2b_new, T2ab_new).

read_starting_guess4amplitudes(EC::ECInfo, ::Val{level}, spins...)

Read starting guess for excitation level.

The guess will be read from T_vo, T_VO, T_vvoo etc files. If the file does not exist, the guess will be a zeroed-vector.

rotate_U2pseudocanonical(EC::ECInfo, UaiX)

Diagonalize ϵv - ϵo transformed with UaiX (for update). Return eigenvalues and rotated UaiX

rotation_matrix(EC::ECInfo, T1; full=false, beta=false)

Make the integrals rotation matrix with T1.

If full is true, the rotation matrix will be made with core and deleted orbitals included.

save_current_doubles(EC::ECInfo, T2a, T2b, T2ab; prefix="T")

Save current doubles amplitudes T2a, T2b, and T2ab to files prefix*"_vvoo", prefix*"_VVOO", and prefix*"_vVoO"

save_current_doubles(EC::ECInfo, T2; prefix="T")

Save current doubles amplitudes T2 to file prefix*"_vvoo"

save_current_singles(EC::ECInfo, T1a, T1b; prefix="T")

Save current singles amplitudes T1a and T1b to files prefix*"_vo" and prefix*"_VO"

save_current_singles(EC::ECInfo, T1; prefix="T")

Save current singles amplitudes T1 to file prefix*"_vo"

save_or_start_file(EC::ECInfo, type, excitation_level, save=true)

Return filename and description for saving or starting amplitudes/lagrange multipliers.

type is either "T" for amplitudes or "LM" for Lagrange multipliers. excitation_level is the excitation level of the amplitudes (1, 2 etc.) If save is true, the filename for saving is returned, otherwise the filename for starting.

spin_project!(EC::ECInfo, T1a, T1b)

Spin-project singles amplitudes/residuals.

Only possible for high-spin states.

spin_project!(EC::ECInfo, T1a, T1b, T2a, T2b, T2ab)

Spin-project singles and doubles amplitudes/residuals.

Only possible for high-spin states.

spin_project_amplitudes(EC::ECInfo, with_singles=true)

Spin-project singles (if withsingles) and doubles amplitudes from files `"Tvo","TVO","Tvvoo","TVVOO"and"TvVoO"`.

transform_amplitudes2lagrange_multipliers!(Amps1, Amps2)

Transform amplitudes to first guess for Lagrange multipliers.

The amplitudes are transformed in-place.

triples_4ext!(EC::ECInfo, R3, T3)

Calculate 4-external contraction with triples amplitudes and store the result in R3.

triples_4ext_aa!(EC::ECInfo, R3a, R3aab, T3a, T3aab)

Calculate αα 4-external contractions with triples amplitudes for the unrestricted case. Uses antisymmetrized integrals $\langle cd||ab\rangle$ compressed to $a<b$.

$R^{ijk}_{cde} += \sum_{a<b} \langle cd||ab\rangle T^{ijk}_{eab}$ $R^{ijK}_{cdA} += \sum_{a<b} \langle cd||ab\rangle T^{ijK}_{abA}$

triples_4ext_aab_ab!(EC::ECInfo, R3aab, T3aab)

Calculate αβ 4-external contraction of $d_{vVvV}$ with $T^{ijK}_{abC}$ for the unrestricted case. No integral compression (different-spin contracted indices), but loops over one occupied index to reduce intermediate size.

$R^{ijK}_{bcB} -= X_{bcB}^{ijK} - X_{cbB}^{ijK}$ where $X_{bcB}^{ijK} = \sum_{a,A} \langle bB|aA\rangle T^{ijK}_{caA}$

triples_4ext_abb_ab!(EC::ECInfo, R3abb, T3abb)

Calculate αβ 4-external contraction of $d_{vVvV}$ with $T^{iIJ}_{aCB}$ for the unrestricted case. No integral compression (different-spin contracted indices), but loops over one occupied index to reduce intermediate size.

$R^{iIJ}_{bBC} -= X_{bBC}^{iIJ} - X_{bCB}^{iIJ}$ where $X_{bBC}^{iIJ} = \sum_{a,A} \langle bB|aA\rangle T^{iIJ}_{aCA}$

triples_4ext_bb!(EC::ECInfo, R3b, R3abb, T3b, T3abb)

Calculate ββ 4-external contractions with triples amplitudes for the unrestricted case. Uses antisymmetrized integrals $\langle CD||AB\rangle$ compressed to $A<B$.

$R^{IJK}_{CDE} += \sum_{A<B} \langle CD||AB\rangle T^{IJK}_{EAB}$ $R^{iIJ}_{aCB} += \sum_{A<B} \langle CD||AB\rangle T^{iIJ}_{aAB}$

try2save_amps!(EC::ECInfo, ::Val{excitation_level}, amps...; type="T")

Save amplitudes (type="T") or Lagrange multipliers (type="LM") to file EC.options.cc.save[_lm]*"_excitation_level".

try2save_doubles!(EC::ECInfo, doubles...; type="T")

Save doubles amplitudes (type="T") or Lagrange multipliers (type="LM") to file EC.options.cc.save[_lm]*"_2".

try2save_singles!(EC::ECInfo, singles...; type="T")

Save singles amplitudes (type="T") or Lagrange multipliers (type="LM") to file EC.options.cc.save[_lm]*"_1".

try2start_amps(EC::ECInfo, ::Val{excitation_level}; type="T")

Read amplitudes (type="T") or Lagrange multipliers (type="LM") from file EC.options.cc.start[_lm]*"_excitation_level".

try2start_doubles(EC::ECInfo; type="T")

Read doubles amplitudes (type="T") or Lagrange multipliers (type="LM") from file EC.options.cc.start[_lm]*"_2".

try2start_singles(EC::ECInfo; type="T")

Read singles amplitudes (type="T") or Lagrange multipliers (type="LM") from file EC.options.cc.start[_lm]*"_1".

try_fetch_restricted_starting_amplitudes(EC::ECInfo)

Try to read and project restricted amplitudes from a TREXIO file.

The logic is:

• If wf.start is not empty: read amplitudes, MOs, and basis from wf.start, then project amplitudes to the current MO basis (obtained via fetch_orbitals from wf.dump).
• If wf.start is empty: try to read amplitudes from wf.dump (no projection needed).

The occupied space projection rank is checked and a warning is printed if not full rank.

Returns (T1, T2, success::Bool) for closed-shell case.

try_fetch_starting_determinants(EC::ECInfo; OPattern=UInt64, nstates=1)

Try to read determinants and CI coefficients from a TREXIO file for CIPHI restart.

The logic follows CC amplitude restart:

• If wf.start is not empty: read from wf.start file(s)
• If wf.start is empty: try to read from wf.dump file(s)

For multi-state, reads from separate files per TREXIO standard.

Returns (dets, coeffs, success::Bool).

try_fetch_unrestricted_starting_amplitudes(EC::ECInfo)

Try to read and project unrestricted amplitudes from a TREXIO file.

The logic is:

• If wf.start is not empty: read amplitudes, MOs, and basis from wf.start, then project amplitudes to the current MO basis (obtained via fetch_orbitals from wf.dump).
• If wf.start is empty: try to read amplitudes from wf.dump (no projection needed).

The occupied space projection rank is checked and a warning is printed if not full rank.

Returns (T1a, T1b, T2a, T2b, T2ab, success::Bool) for unrestricted case.

update_deco_doubles(EC, R2; use_shift=true)

Update decomposed doubles amplitudes.

If R2 is $R^{ij}_{ab}$, the update is calculated using update_doubles(EC, R2, use_shift=use_shift).

update_deco_triples(EC, R3, use_shift=true)

Update decomposed triples amplitudes.

Note that the sign of the residual is opposite to the usual definition of the triples residual and therefore the update is calculated using a positive denominator...

update_doubles(R2, ϵo1, ϵv1, ϵo2, ϵv2, shift)

Calculate update for doubles amplitudes.

update_doubles!(EC::ECInfo, T2, R2)

Update doubles amplitudes in T2 with R2.

update_doubles!(EC::ECInfo, T2a, T2b, T2ab, R2a, R2b, R2ab)

Update doubles amplitudes in T2a, T2b, T2ab with R2a, R2b, R2ab.

update_doubles(EC::ECInfo, R2; spincase::Symbol=:α, antisymmetrize=false, use_shift=true)

Calculate update for doubles amplitudes for a given spincase∈{:α,:β,:αβ}.

update_singles!(EC::ECInfo, T1, R1; damp=1.0)

Update singles amplitudes in T1 with R1. The update is scaled by damp (1.0 = full step; used e.g. as the orbital-rotation damping in OQV methods).

update_singles!(EC::ECInfo, T1a, T1b, R1a, R1b; damp=1.0)

Update singles amplitudes in T1a, T1b with R1a, R1b, scaled by damp.

update_singles(R1, ϵo, ϵv, shift)

Calculate update for singles amplitudes.

update_singles(EC::ECInfo, R1; spincase::Symbol=:α, use_shift=true)

Calculate update for singles amplitudes for a given spincase∈{:α,:β}.

update_triples!(EC::ECInfo, T3, R3)

Update triples amplitudes in T3, with R3.

update_triples!(EC::ECInfo, T3a, T3b, T3aab, T3abb, R3a, R3b, R3aab, R3abb)

Update triples amplitudes in T3a, T3b, T3aab and T3abb with R3a, R3b, R3aab and R3abb.

add_singles2doubles!(T2aa, T2bb, T2ab, T1a, T1b)

Add singles to doubles amplitudes.

add_singles2doubles!(T2, T1; make_contravariant=true)

Add singles to doubles amplitudes.

If make_contravariant is true, the amplitudes will be made contravariant.

antisymmetrize_34!(d, trivv)

Antisymmetrize dims 3 and 4 of d in-place and compress to triangular. Returns d[:,:,trivv] where d[p,q,a,b] → d[p,q,a,b] - d[p,q,b,a] for a < b.

correlation_reference_orbital_data(EC::ECInfo)

Pre-fetch the orbital data to store with a restarted wavefunction, using the correlation reference orbitals — the ones the FCIDUMP was actually built from (load_orbitals_for_correlation) — rather than the raw orbitals on file.

For a plain restart (the dump already holds the reference in the current AO basis) the two coincide. But for a dump=""+start restart across a geometry (or basis) change, the reference is the projected and re-orthonormalized orbital set expressed in the current AO basis, while the raw orbitals on file are still in the old basis. Storing the raw orbitals then writes old-basis coefficients tagged with the new AO basis — inconsistent with both the current geometry and the stored (reference-basis) amplitudes — so a subsequent same-geometry restart cannot resume at the stored solution and re-optimizes over several iterations. Storing the correlation reference keeps the orbitals, amplitudes, and AO basis mutually consistent, so the restart resumes in one iteration.

Returns the raw data unchanged for an FCIDUMP-only run (no molecular system) or a rotation-only dump (no AO basis stored), where there are no orbital coefficients to re-express — the caller then stores a rotation as before. Must be called before opening the store file for writing (dump == store case).

empty_restricted_amplitudes(EC::ECInfo)

Return empty restricted amplitudes of the correct size.

Returns (T1, T2, false) where T1 and T2 are zero arrays.

empty_unrestricted_amplitudes(EC::ECInfo)

Return empty unrestricted amplitudes of the correct size.

Returns (T1a, T1b, T2a, T2b, T2ab, false) where all arrays are zero.

update_triples(R3, ϵo1, ϵv1, ϵo2, ϵv2, ϵo3, ϵv3, shift)

Calculate update for triples amplitudes.

update_triples(EC::ECInfo, R3; spincase::Symbol=:α, antisymmetrize=false, use_shift=true)

Calculate update for triples amplitudes for a given spincase∈{:α,:β,:ααβ,:αββ}.