Analytic nuclear forces for the Γ-CCM (aiccm2026dev-a) – Methods

Owned by the Γ-CCM analytic-forces line (paper #3). Companion to aiccm2026dev_a_theory.md (the -a method paper) and aiccm_a_position.md (the re-derivation/adjudication memo). The implementation spec this Methods section is the rigorous backing for is ../aiccm2026dev_a_analytic_gradient.md; the finite-difference gate is ccm_numerical_gradient (properties.py).

Discipline (as the position memo): derive every link and prove it; every equation below becomes an inline comment at the implementation site (CLAUDE.md §8); the finite-difference gate is the arbiter of correctness, to ~1e-6 Ha/bohr on the H₂ chain.

This document derives the analytic nuclear gradient dE_CCM/dR of the Γ-CCM Hartree-Fock energy from first principles. The energy being differentiated is the closed-shell RHF-CCM energy per reference cell; the derivation gives the variational Hellmann-Feynman/Pulay split, proves that the WSSC fold commutes with the nuclear derivative, states the cyclic image-sum, derives the adjoint (Lᵀ) contraction that turns each folded-integral derivative into a molecular gradient call on the padded cluster, gives the nuclear-repulsion gradient in closed form, and proves the translational sum rule. RHF is treated in full; UHF, KS, and post-HF are the outlook (§1.7).

The numbering of the subsections mirrors §1 of handovers/HANDOVER_AICCM_GRADIENT.md.

---

0. Setup, notation, and the three index sets

0.1 The Γ-CCM energy

The Γ-CCM construction (Peintinger & Bredow 2014) places the crystal on a finite Born-von-Kármán torus – the reference supercell of size nrep = (N₁,N₂,N₃) – and defines its finite Hamiltonian by Wigner-Seitz-supercell (WSSC) integral weighting. A real-supercell evaluation, where used, is a representation of this already specified weighted Hamiltonian; it does not identify Γ-CCM with the separate finite-character χ-CCM construction. Every molecular HF/KS/post-HF kernel then applies unchanged to the WSSC-weighted (“effective”) integrals. The closed-shell RHF-CCM energy per reference cell is

E_CCM = Σ_μν P_μν h^CCM_μν
      + ½ Σ_μνλσ P_μν P_λσ [ (μν|λσ)^CCM − ½ (μλ|νσ)^CCM ]
      + V_nn^CCM ,                                                         (0.1)

equivalently E_CCM = ½ Σ_μν P_μν (h^CCM_μν + F^CCM_μν) + V_nn^CCM with F^CCM = h^CCM + J^CCM[P] − ½ K^CCM[P]. Here:

• P_μν = 2 Σ_i^occ C_μi C_νi is the converged closed-shell density on the home (reference-cell, nbf-dimensional) AO basis;

• h^CCM = T^CCM + V^CCM, the WSSC-folded core Hamiltonian;

• (μν|λσ)^CCM is the WSSC-folded four-center electron-repulsion tensor (Mulliken/chemist order), the symmetric bridge of ccm_eri_symmetric for method="aiccm2026dev-a";

• V_nn^CCM is the WSSC-weighted nuclear repulsion per reference cell.

It is convenient to write the two-electron energy with the closed-shell HF AO 2-particle density (2-RDM)

Γ_μνλσ = P_μν P_λσ − ½ P_μλ P_νσ ,   E_2e = ½ Σ_μνλσ Γ_μνλσ (μν|λσ)^CCM .     (0.2)

(One checks the exchange relabelling ν↔λ reproduces −¼ Σ P_μν P_λσ (μλ|νσ)^CCM.) Γ is exactly the arbitrary AO 2-RDM that two_electron_gradient_casscf(basis, mol, gamma_ao_flat) (gradient.cpp:1473, bound at bindings.cpp:2840) contracts against ∂(μν|λσ)/∂R (Helgaker Eq. 12.5.7); this is what makes the 2-e term a Python-only call (§2).

0.2 The WSSC fold L

Every effective matrix is a linear, geometry-locally-constant fold L of a padded-cluster molecular integral. The padded cluster (PaddedCluster, padded.py) is the reference supercell plus its cluster-lattice image cells g, built as one ordinary molecule on which standard molecular integrals are evaluated; the CCM matrix over the reference cell is assembled by selecting home↔image sub-blocks and applying the WSSC weights w^g_{AB} = cell_weight_matrices()[g][A,B]. Concretely:

• S^CCM = L₂(S_pad), T^CCM = L₂(T_pad) – _fold_two_center (padded.py:217);

• V^CCM = L_V(V_pad) – ccm_nuclear (padded.py:247), per-nucleus-image weighted;

• (μν|λσ)^CCM = L₄(eri_pad) – ccm_eri_symmetric (padded.py:361);

• V_nn^CCM – ccm_nuclear_repulsion (padded.py:437).

All folds draw the same weights from CCMSystem.cell_weight_matrices() (system.py:397).

0.3 Three index sets – the crux of every scatter

The derivation lives on three distinct atom-index sets; conflating them is the one way to get the adjoint wrong, so we fix the notation once.

symbol	set	size	code
β	unit-cell atom	n_basis_atoms	ccm.n_basis_atoms, ccm.unit_atom_pos
A	supercell atom (the “real” cluster atoms u_j)	n_atoms = n_cells·n_basis_atoms	ccm.n_atoms, ccm.atom_positions
(A,g)	padded-cluster atom: supercell atom A in image cell g	n_pad_atom	PaddedCluster.pad_atom_of

The supercell is built for i: for j: for k: for β, so a supercell atom factorises as A = ((i·N₂+j)·N₃+k)·n_basis_atoms + β; the inverse map β = atom_beta[A], (i,j,k) = atom_cell[A] is CCMSystem._decompose_atoms (system.py:282).

Resolved open detail – the WSSC weights are SUPERCELL-indexed. cell_weight_matrices()[g] is an (n_atoms, n_atoms) matrix W[g][A,B] indexed by supercell atoms (system.py:465: Wg = Tg[beta[:,None], beta[None,:], pidx] with beta of length n_atoms and pidx of shape (n_atoms, n_atoms)); the docstring states the same, and _fold_two_center consumes it as w_atom[ao[:,None], ao[None,:]] with ao = ccm.ao_atom mapping each AO to its supercell atom. Every adjoint scatter below is therefore indexed by A, never by β. The fold-back to the gradient output happens only at the end, by summing the per-supercell-atom (or per-pad-atom) gradient over the supercell copies of each unit-cell atom (§1.3).

The gate ccm_numerical_gradient (properties.py:585) displaces a unit-cell atom β (_displaced_unit_system shifts unit_system.unit_cell[β], then rebuilds the whole CCMSystem, so all n_cells supercell copies and all their padded images move together), and returns (n_basis_atoms, 3). The analytic gradient must match that shape and that displacement semantics.

---

1. The theory (paper §Methods)

1.1 Variational gradient and the Hellmann-Feynman/Pulay split

Let R_β be a Cartesian coordinate of unit-cell atom β. The MOs C(R) solve the CCM-RHF equations subject to orthonormality in the folded overlap, Cᵀ S^CCM(R) C = 1. Introduce the constrained Lagrangian

𝓛(C, ε; R) = E_CCM(C; R) − Σ_ij ε_ij ( Σ_μν C_μi S^CCM_μν(R) C_νj − δ_ij ).      (1.1a)

At the SCF solution the orbital-variation stationarity ∂𝓛/∂C = 0 holds and the multipliers are the canonical orbital energies, ε_ij = ε_i δ_ij. Because 𝓛 = E on the constraint manifold, the total energy derivative along the solution path C(R) equals the partial (explicit, integral-only) derivative of 𝓛:

dE_CCM/dR_β = ∂𝓛/∂R_β |_{C,ε fixed}
            = ∂E_CCM/∂R_β |_C − Σ_ij ε_ij Σ_μν C_μi (∂S^CCM_μν/∂R_β) C_νj .       (1.1b)

The orbital-response terms Σ_μi (∂E/∂C_μi)(dC_μi/dR) cancel against the constraint-response exactly at a stationary point – this is the standard Hellmann- Feynman argument (Pulay 1969; Helgaker, Jørgensen & Olsen, Molecular Electronic-Structure Theory, Ch. 10 & Eq. 12.5.7), unchanged here because the only modification is that the metric is the folded overlap S^CCM. Writing the explicit derivative out and collecting,

dE_CCM/dR_β =  Σ_μν P_μν ∂h^CCM_μν/∂R_β
            +  ½ Σ_μνλσ P_μν P_λσ ∂[(μν|λσ)^CCM − ½(μλ|νσ)^CCM]/∂R_β
            −  Σ_μν W_μν ∂S^CCM_μν/∂R_β
            +  ∂V_nn^CCM/∂R_β ,                                                    (1.1)

with the energy-weighted density (Pulay/force-relaxation matrix)

W_μν = Σ_i n_i ε_i C_μi C_νi   (n_i = 2, closed shell)  =  (C·diag(n ε)·Cᵀ)_μν .   (1.1c)

The four terms of (1.1) are the nuclear-attraction/kinetic Hellmann-Feynman term, the two-electron Hellmann-Feynman term (written compactly as ½ Σ Γ ∂(μν|λσ)^CCM with Γ from (0.2)), the Pulay overlap (basis-relaxation) term, and the nuclear-repulsion term. Equation (1.1) is formally identical to the molecular RHF gradient; the entire CCM content sits inside the folded integrals M^CCM and their derivatives, treated next.

Remark (stationarity sensitivity). The cancellation of the orbital-response terms in (1.1b) holds only at a stationary point – ∂E_CCM/∂κ_{ai} = F^CCM_{ai} = 0 for every occupied-virtual rotation. A loosely converged SCF leaves a residual F_{ai}, which leaks into the analytic gradient at first order through the Pulay term, even though the energy itself is only second-order sensitive to it. The analytic force therefore demands a tightly converged SCF (here the DIIS commutator F^CCM D S^CCM − S^CCM D F^CCM driven well below 10⁻⁶), whereas a finite- difference force is insensitive to the same non-stationarity. This asymmetry is a property of variational analytic gradients in general (Helgaker-Jørgensen-Olsen, Ch. 10), restated here because it is the one practical prerequisite for the formula to reproduce the finite-difference reference.

1.2 The fold commutes with the nuclear derivative

Each fold L (for M ∈ {S,T,V,(··|··)}) is a fixed linear map: a sum over cells g of WSSC weights w^g_{AB} times a selection of padded sub-blocks (cell_to_cols[g]), optionally bra↔ket-symmetrised. Both the weights and the selection maps are constants under the nuclear derivative (Lemma 1.2.1), so

∂M^CCM/∂R_β = L( ∂M_pad/∂R_β )   for  M ∈ {S, T, V, (··|··)} .                    (1.2)

Lemma 1.2.1 (WSSC weights are locally position-independent). The weight w^g_{AB} is determined entirely by which minimum-image Wigner-Seitz cell of the cluster lattice the class displacement d = (r_{βA} − r_{βB}) + Δcell·a falls into (minimum_image, system.py:436). This is a piecewise-constant function of the atomic positions: it is locally constant on the open interior of every WS cell and changes only across the measure-zero WS-boundary set (where minimum_image splits the weight equally among tied cells within weight_tol_bohr). Hence ∂w^g_{AB}/∂R = 0 for any geometry not on a boundary tie, and (1.2) holds there.

Caveat (boundary ties). At a boundary tie the cluster energy E_CCM(R) has a kink: the analytic gradient (1.1)-(1.2) is one element of the subgradient, and a central finite difference straddling the tie would disagree by O(1) in the tie-crossing component. The gate must therefore be run at a generic geometry (e.g. an H₂ chain at a non-commensurate bond length); there the weight is smooth and analytic == FD to ~1e-6. This is the only place where the analytic and finite-difference gradients are allowed to differ, and it is physical, not a bug. (The handover §0 statement “the WSC weights depend only on the lattice” is the interior-of-cell version of this lemma; Lemma 1.2.1 is the precise form.)

1.3 The cyclic image derivative (the image-sum)

Unit-cell atom β is realised in the padded cluster as the entire set of padded atoms

pad(β) = { (A, g) : atom_beta[A] = β,  g ∈ cells } ,                              (1.3a)

i.e. all n_cells supercell copies of β (atom_beta[A]=β) and all their image cells g. The cyclic (Born-von-Kármán) constraint – the defining content of the gate’s _displaced_unit_system + full CCMSystem rebuild – moves every member of pad(β) rigidly together. Therefore

∂/∂R_β = Σ_{A : atom_beta[A]=β}  Σ_{g}  ∂/∂R_{(A,g)} .                            (1.3)

Operationally: compute the ordinary per-pad-atom molecular gradient grad_pad[p] = ∂(·)/∂R_{pad_atom_of[p]} on the padded cluster, then image-sum

grad_cell[β] = Σ_{p : atom_beta[ pad_atom_of[p].A ] = β}  grad_pad[p] .           (1.3b)

This two-stage fold-back (over padded cells g and over supercell copies A) is the translational fold-back that ccm_numerical_gradient performs implicitly by displacing the unit-cell atom and propagating to every image.

1.4 The adjoint contraction Lᵀ (turns each fold into a molecular gradient call)

Combining (1.1), (1.2) and linearity, each energy-gradient term has the form Σ_μν C^CCM_μν ∂M^CCM_μν/∂R_β for some symmetric contraction coefficient C^CCM ∈ {P, W} (or the 2-RDM Γ^CCM for the four-center). Using M^CCM = L(M_pad) and the adjoint Lᵀ of the linear fold,

Σ C^CCM · ∂M^CCM/∂R_β = Σ C^CCM · L(∂M_pad/∂R_β) = Σ Lᵀ(C^CCM) · ∂M_pad/∂R_β .   (1.4)

So we scatter the CCM contraction coefficient onto the padded layout with the same WSSC weights (Lᵀ), then hand the scattered density to the ordinary molecular gradient routine on the padded cluster, and finally image-sum (1.3b). This is the device that makes the whole gradient a sequence of existing molecular gradient calls – no new C++. The adjoints, term by term:

(L₂ᵀ) two-center (P → P_pad, W → W_pad). The fold (_fold_two_center) is

M^CCM_μν = ½ ( out_μν + out_νμ ),   out_μν = Σ_g w^g_{a(μ)a(ν)} M_pad[μ, ν@g] ,    (1.4a)

with a(μ)=ao_atom[μ] (supercell atom of AO μ), ν@g = cell_to_cols[g][ν], home rows μ ∈ {0..nbf−1}. For symmetric C^CCM and symmetric M_pad, the energy term collapses (the ½(out+outᵀ) and the μ↔ν symmetry of C^CCM cancel the ½):

Σ_μν C^CCM_μν M^CCM_μν = Σ_g Σ_μν C^CCM_μν w^g_{a(μ)a(ν)} M_pad[μ, ν@g]
                       = Σ_{p,q} X_pad[p,q] M_pad[p,q] ,                          (1.4b)
   X_pad[μ, ν@g] += w^g_{a(μ)a(ν)} C^CCM_μν   (scatter; home rows μ),

and since only the symmetric part of X_pad survives contraction with the symmetric ∂M_pad/∂R, the padded density fed to the molecular routine is

P_pad = ½ ( X_pad + X_padᵀ ) .                                                    (1.4c)

L₂ᵀ is then: P_pad = L₂ᵀ(P) for the 1-e Hellmann-Feynman term and W_pad = L₂ᵀ(W) for the overlap Pulay term; drive one_electron_gradient_contribution(basis_pad, mol_pad, P_pad) (kinetic part) and overlap_gradient_contribution(basis_pad, mol_pad, W_pad) and image-sum.

(L_Vᵀ) nuclear attraction – the subtlest. ccm_nuclear (padded.py:247) weights each nucleus image c = (C₀, g_c) by ω_{A,c} = w^{g_c}_{A,C₀} and folds the bra/ket pair independently:

V^CCM_μν = ½( out + outᵀ ),
out_μν = Σ_c Σ_{g_ν} w^{g_ν}_{a(μ)a(ν)} · ½(ω_{a(μ),c}+ω_{a(ν),c}) · I^c_pad[μ, ν@g_ν] , (1.4d)
   I^c_pad[p,q] = ⟨ p | −Z_{C₀}/|r−R_c| | q ⟩   (signed/attractive, padded basis).

The gradient of Σ P V^CCM has two physically distinct pieces:

1. Basis-center (bra/ket) derivatives of I^c_pad: fold exactly like L₂ᵀ, but with the combined weight w^{g_ν}·½(ω_{a(μ)}+ω_{a(ν)}) and the operator evaluated for nucleus image c. The derivative lands on the unit-cell atoms owning μ and ν.

2. Hellmann-Feynman nucleus-position derivative ∫ φ_μ ∂(−Z_{C₀}/|r−R_c|)/∂R_{C₀} φ_ν, weighted by the same combined weight, summed onto the unit-cell atom that owns nucleus image c (i.e. atom_beta[C₀]).

The natural primitive is compute_external_charge_density_gradient(basis, mol, D, charges, positions) (bindings.cpp:2868), whose result carries atom_grad (the bra+ket basis-center derivatives, piece 1) and point_grad (the per-point-charge / nucleus-position derivatives, piece 2) separately – exactly the two pieces above. Per nucleus image c, build the combined-weighted density D^c_pad = L₂ᵀ-style scatter of P with weight w^{g_ν}·½(ω_{a(μ),c}+ω_{a(ν),c}), run the external-charge gradient with the single charge Z_{C₀} at R_c, route atom_grad to the AO-owner unit-cell atoms and point_grad[c] → grad_cell[atom_beta[C₀]], then sum over c. (Equivalently, mirror ccm_nuclear’s per-c compute_nuclear construction and differentiate.)

(L₄ᵀ) four-center (Γ^CCM → Γ_pad). ccm_eri_symmetric (padded.py:361) builds

(μν|λσ)^CCM = ½( eff + effᵀ_braket ),
eff[μν,λσ] = Σ_{g_c,g_e} w^0_{a(μ)a(ν)} · bridge_{μνλσ}^{g_c g_e} · w^{g_e−g_c}_{a(λ)a(σ)}
                       · eri_pad[μ, ν@0, λ@g_c, σ@g_e] ,                          (1.4e)
bridge = ¼( w^{g_c}_{a(μ)a(λ)} + w^{g_c}_{a(ν)a(λ)} + w^{g_e}_{a(μ)a(σ)} + w^{g_e}_{a(ν)a(σ)} ),
effᵀ_braket = eff.transpose(λσ ↔ μν) .

By the same adjoint identity, the contraction Σ Γ^CCM (μν|λσ)^CCM equals Σ Γ_pad · eri_pad with the scatter

Γ̄ = ½( Γ^CCM + Γ^CCM.transpose(2,3,0,1) )                 (adjoint of the closing braket transpose),
Γ_pad[μ, ν@0, λ@g_c, σ@g_e] += w^0_{a(μ)a(ν)} · bridge_{μνλσ}^{g_c g_e}
                                 · w^{g_e−g_c}_{a(λ)a(σ)} · Γ̄_μνλσ .              (1.4f)

Flatten Γ_pad row-major (((μ·nb_pad+ν)·nb_pad+λ)·nb_pad+σ) and feed two_electron_gradient_casscf(basis_pad, mol_pad, Γ_pad_flat); image-sum. The dense n_pad_ao⁴ tensor confines this to small / 1-D / 2-D clusters (the _check_padded_eri_size guard), the same regime as the energy build.

1.5 The nuclear-repulsion gradient (build-order step 1)

V_nn^CCM (padded.py:437) is the WSSC-weighted lattice sum over supercell atom pairs and cells, with the on-site (A=B, g=0) term removed:

V_nn^CCM = ½ Σ_g Σ_{A,B}' w^g_{AB} Z_A Z_B / |d_{ABg}| ,   d_{ABg} = r_A − (r_B + R_g) ,  (1.5a)

R_g = g·A_cluster, the prime excluding (A=B, g=0), and the ½ removing the ordered double-count. Differentiating (1.5a) directly (weights constant, Lemma 1.2.1), the per-supercell-atom gradient is

∂V_nn^CCM/∂r_A = − Σ_g Σ_{B}'  w^g_{AB} Z_A Z_B  d_{ABg} / |d_{ABg}|³ ,            (1.5b)

(the A-as-bra and A-as-ket contributions are equal under g↔−g, A↔B, cancelling the ½), and the gradient output for unit-cell atom β is the supercell image-sum (1.3b restricted to the supercell, since V_nn involves no padded cells beyond the R_g offsets):

∂V_nn^CCM/∂R_β = Σ_{A : atom_beta[A]=β}  ∂V_nn^CCM/∂r_A .                          (1.5c)

This is isolated and the simplest term to gate against a finite difference of ccm_nuclear_repulsion; it is build-order step 1. In code we differentiate the exact pair list the energy uses (every (g,A,B) term with w^g_{AB}≠0, (A,g)≠(B,0)), accumulating ±w Z_A Z_B d/|d|³ onto the two endpoints, so the analytic and finite-difference values agree by construction (away from boundary ties).

1.6 The translational sum rule

Theorem 1.6.1. Σ_β ∂E_CCM/∂R_β = 0 (vector), for every Cartesian component.

Proof. E_CCM depends on the atomic positions only through differences: the molecular integrals over the padded cluster depend on AO-center separations; the nuclear repulsion (1.5a) depends on d_{ABg} = r_A − r_B − R_g; and the WSSC weights depend on the class displacements (r_{βA}−r_{βB}) + Δcell·a (and, being piecewise-constant, are exactly invariant under a rigid shift – the displacement is unchanged). A rigid translation of the whole crystal by an arbitrary constant t moves every unit-cell atom by t (R_β → R_β + t for all β) and leaves all of these differences – hence E_CCM – invariant. Therefore

0 = d E_CCM(R + t·1) / dt |_{t=0} = Σ_β (∂E_CCM/∂R_β) · 1 = Σ_β ∂E_CCM/∂R_β .       (1.6)

∎ Because the WSSC weights are exactly (not just locally) invariant under rigid translation, the sum rule holds even at boundary ties – unlike the per-atom gradient (Lemma 1.2.1). It is both a paper result and the implementation’s built-in check: ccm_numerical_gradient notes the per-cell forces sum to ≈ 0, and the analytic gradient must satisfy Σ_β grad_cell[β] = 0 to working precision.

1.7 Scope of the equations

• Closed-shell RHF – fully derived above; the implementation target.

• UHF – implemented and validated (run_ccm_uhf_gradient). Equation (1.1) holds with the total density P = P^α + P^β in the 1-e term, the total energy-weighted density W = W^α + W^β (each W^σ = Σ_i ε^σ_i C^σ_i C^σ_iᵀ, occupation 1 per spin orbital) in the Pulay term, and the spin-resolved AO 2-RDM

  Γ_μνλσ = P_μν P_λσ − P^α_μλ P^α_νσ − P^β_μλ P^β_νσ                              (1.7)
  

  (total-density Coulomb, same-spin exchange; reduces to (0.2) when P^α=P^β=P/2) in the four-center term. The folds L, the adjoints Lᵀ, and the image-sum are spin-independent, so the same term cores serve both references. Validated to machine precision against the molecular UHF analytic gradient (LiH⁺/STO-3G 6.9e-11, HeH/6-31G 2.7e-9) and to the FD floor with O(h²) Richardson on the open-shell HeH chain. Edge case: a fully-occupied spin manifold (n_occ_σ = nbf) makes P^σ = S⁻¹ and the DIIS commutator vanish identically for that spin, so the SCF check is blind to it and ε^σ (hence W^σ) may be under-converged – a minimal-basis SCF pathology (HeH/STO-3G), not a gradient defect; the generic n_occ_σ < nbf case is exact.

• KS-CCM (RKS) – implemented and validated (run_ccm_rks_gradient). One structural fact decides the whole derivation: in the KS-CCM Fock F = h^CCM + J^CCM + V_xc[ρ] (− α K^CCM) (dft.py) only the Coulomb / exact exchange (the four-center) carries the WSSC fold; the exchange-correlation term is the ordinary molecular semi-local functional of the cluster density on a standard molecular Becke grid over the supercell – it does NOT fold. Hence

  dE_KS/dR =  Σ P ∂h^CCM/∂R                       (reuse: kinetic L₂ᵀ + V_ne L_Vᵀ)
           +  ½ Σ Γ_KS ∂(μν|λσ)^CCM/∂R            (reuse: L₄ᵀ, Γ_KS = P⊗P − (α/2) P⊗P_swap)
           −  Σ W ∂S^CCM/∂R                       (reuse: overlap Pulay L₂ᵀ)
           +  ∂E_xc/∂R |_supercell                (molecular XC-Pulay, no fold)
           +  ∂V_nn^CCM/∂R ,                                                       (1.7b)
  

  with α the HF-exchange fraction (α=0 pure ⇒ Coulomb-only Γ_KS=P⊗P; α=1 recovers RHF). The first three terms and V_nn reuse the HF cores verbatim; the XC term is the molecular XC-Pulay evaluated on the supercell (n_atoms) density/grid – the same supercell + Becke grid the KS-CCM energy uses – folded back by the supercell image-sum (fold_supercell_gradient_to_cell; no padded images, since the XC lives on the supercell, not the padded cluster). A small binding xc_pulay_gradient_rks (bindings.cpp) surfaces the molecular XC-Pulay kernel (gradient.cpp:805/896) – the one piece the HF gradient did not need. Validated to machine precision against the molecular compute_gradient_rks in the isolated limit for PBE (2.2e-16), PBE0 (1.4e-16), and B3LYP (5.6e-17); periodic FD on the H₂ chain agrees at the grid-fixed floor (~4e-7, the Becke weight-derivative term both the analytic and the molecular DFT gradient neglect – tighten the grid to close it). UKS is implemented too (run_ccm_uks_gradient): the same decomposition with total P=P^α+P^β, total W=W^α+W^β, the open-shell 2-RDM Γ=P⊗P − α(P^α⊗P^α + P^β⊗P^β), and the _uks XC-Pulay binding (xc_pulay_gradient_uks) on the supercell spin densities – the HF cores are already spin-general. Machine precision vs the molecular compute_gradient_uks in the isolated limit (PBE/PBE0/B3LYP); periodic FD at the grid-fixed floor on a gapped (dilute) open-shell chain. Two DFT-specific caveats, both inherited from vibe-qc’s molecular DFT gradient, not defects: (i) the DFT force satisfies the translational sum rule only to the grid-fixed floor (~1e-6) – the neglected Becke weight derivative is precisely the term that would make it exact – so Σ_β grad = 0 holds to that tolerance for KS/UKS (machine-exact for HF); (ii) a metallic (dense) open-shell chain is spin-unstable under a GGA, so a finite-difference check there is ill-posed (the displaced SCFs jump spin-broken branches) while the analytic force stays well-defined – validate periodic UKS on a gapped cell.

• MP2 + UMP2 + CCSD – implemented and validated (run_ccm_mp2_gradient / run_ccm_ump2_gradient / run_ccm_ccsd_gradient): the relaxed-density gradients, derived in §1.8 (MP2/UMP2) and §1.9 (CCSD) below. The structural statement is the same as for HF/KS: replace P, W, Γ by the method’s relaxed one- and two-particle densities (built from the CCM amplitudes, the Λ multipliers, and the CCM-CPHF/Z-vector orbital response), folded and scattered through the identical adjoints L₂ᵀ, L_Vᵀ, L₄ᵀ. Open-shell CCSD (UCCSD) and the (T) correction to the forces remain outlook; cf. the CASSCF general-Γ implementation precedent (_casscf.py).

The dense-n_pad⁴ regime confines the analytic four-center to small / 1-D / 2-D clusters; genuine 3-D analytic forces need a scalable weighted ERI-gradient C++ kernel (a build_jk_ccm_weighted sibling) – a later milestone, not part of this Methods section.

1.8 The MP2 relaxed-density gradient (post-HF forces)

The CCM-MP2 correlation energy (mp2.py) is the ordinary closed-shell RMP2 expression evaluated on the canonical CCM-RHF orbitals and the WSSC-folded four-center:

E₂ = Σ_ijab t̃_ijab g_iajb ,      g_iajb = (ia|jb)^CCM  (chemist order),
t_ijab = g_iajb / D_ijab ,       D_ijab = ε_i + ε_j − ε_a − ε_b ,
t̃_ijab = 2 t_ijab − t_ijba ,                                                    (1.8a)

i,j,k occupied, a,b,c virtual, p,q,r,s general CCM MOs. Every ingredient of (1.8a) is a folded quantity (g from L₄, ε/C from the folded Roothaan-Hall problem), so by §1.2 the AO-derivative pieces of dE₂/dR fold through the same adjoints as HF; the new work is the response of C and ε, handled by the standard MP2 Lagrangian (Pople-Krishnan-Schlegel-Binkley 1979; Handy-Schaefer 1984; Helgaker-Jørgensen-Olsen Ch. 12), assembled here with every matrix CCM-folded.

(1.8.1) Amplitude stationarity (no amplitude response). E₂ is the stationary value of the (generalized) Hylleraas functional

J[t] = 2 Σ t̃_ijab g_iajb + Σ t̃_ijab (Q t)_ijab ,
(Q t)_ijab = Σ_c [ F_ac t_ijcb + F_bc t_ijac ] − Σ_k [ F_ik t_kjab + F_jk t_ikab ] ,  (1.8b)

with F the MO Fock matrix: δJ/δt = 0 ⇔ Q t̃ = −g̃ ⇔ t = g/D in the canonical basis (F = diag(ε)), and J[t*] = E₂. Stationarity in t removes any dt/dR term from the gradient; only the explicit integral derivatives and the orbital response survive. Because (1.8b) is written with full Fock blocks F_oo, F_vv (not the diagonal ε), the functional – hence the derivation – is invariant under rotations within the occupied and within the virtual spaces, which is what licenses evaluating all final expressions in the canonical basis with the symmetric orthonormality gauge for the oo/vv response blocks.

(1.8.2) Orbital response and the Z-vector. Parametrize dC/dR = C U with U + Uᵀ = −S̃, S̃_pq ≡ (Cᵀ (∂S^CCM/∂R) C)_pq (orthonormality), and let superscript [R] denote AO-derivative-at-fixed-C MO transforms (e.g. F^[R] = Cᵀ(∂h^CCM/∂R + ∂G^CCM[P]/∂R|_P)C). Differentiating (1.8a) and collecting the U-independent, occupied-rotation, virtual-rotation, and ε-response pieces gives the intermediates (all g blocks folded):

d_ij = −2 Σ_kab t̃_ikab t_jkab ,        d_ab = +2 Σ_ijc t̃_ijac t_ijbc ,
X_pi = Σ_jab t̃_ijab g_pajb ,           Y_pa = Σ_ijb t̃_ijab g_ipjb ,             (1.8c)

(d is the unrelaxed MP2 correction density; one checks ∂E₂/∂ε_p = d_pp). The occupied-virtual rotations U_ai cannot be eliminated by the orthonormality alone; their coefficient defines the MP2 orbital-rotation gradient (the Lagrangian right-hand side)

L_ai = 4 X_ai − 4 Y_ia + 4 (Cᵀ G^CCM[C d C ᵀ] C)_ai ,                            (1.8d)

with G^CCM[·] = J^CCM − ½K^CCM the folded Fock response and d = d_oo ⊕ d_vv. Instead of solving the coupled-perturbed equations once per nuclear coordinate, one Z-vector solve (Handy-Schaefer 1984) handles all of them:

M z = L ,     M_{ai,bj} = (ε_a − ε_i) δ_ab δ_ij + A_{ai,bj} ,
A_{ai,bj} = 4(ai|bj)^CCM − (ab|ij)^CCM − (aj|ib)^CCM ,                           (1.8e)

– the CCM-CPHF equation: the same electronic Hessian A as molecular CPHF, built from the folded four-center; M is symmetric, and for the dense validation regime it is solved directly.

(1.8.3) Assembled relaxed densities. Substituting z back and moving every S̃-proportional term into the Pulay slot yields the gradient in exactly the HF form (1.1), with three density replacements. The relaxed one-particle correction (MO basis, symmetric)

P^Δ = C [ d_oo  −½z ᵀ ]
        [ −½z   d_vv ] C ᵀ ,      P_tot = P_SCF + P^Δ ;                          (1.8f)

the two-particle density Γ_tot = Γ_HF + Γ_sep + Γ_ns with the separable (Fock-response) and non-separable (amplitude) parts

Γ_sep,μνλσ = 2 P^Δ_μν P_λσ − P^Δ_μλ P_νσ ,
Γ_ns,μνλσ  = 4 Σ_ijab t̃_ijab C_μi C_νa C_λj C_σb ;                               (1.8g)

and the energy-weighted density W_tot = W_HF + W^Δ with (MO blocks, symmetric; G̃ ≡ Cᵀ G^CCM[P^Δ] C with the full P^Δ of (1.8f))

W^Δ_ij = X_ij + X_ji + ½ d_ij (ε_i + ε_j) + 2 G̃_ij ,
W^Δ_ab = Y_ab + Y_ba + ½ d_ab (ε_a + ε_b) ,
W^Δ_ia = W^Δ_ai = 2 Y_ia − ½ z_ai ε_i .                                          (1.8h)

The total CCM-MP2 gradient is then the HF contraction pipeline verbatim – (1.8f)-(1.8h) fed to the same L₂ᵀ/L_Vᵀ/L₄ᵀ adjoint cores of §1.4 plus the V_nn term of §1.5:

dE_MP2/dR_β =  Σ P_tot ∂h^CCM/∂R_β  +  ½ Σ Γ_tot ∂(μν|λσ)^CCM/∂R_β
            −  Σ W_tot ∂S^CCM/∂R_β  +  ∂V_nn^CCM/∂R_β .                          (1.8i)

This is the paper’s structural claim for post-HF forces: the WSSC fold commutes with the entire relaxed-density machinery – amplitudes, CPHF/Z-vector, and density assembly are the molecular equations on folded integrals, and the final AO contractions reuse the HF adjoints unchanged. The translational sum rule (Theorem 1.6.1) applies verbatim (its proof used only translation invariance of E_CCM, which holds for E₂ too).

Validation (all numbers from real runs, H₂ chain STO-3G unless noted).

1. Response identity (isolates d, z, and the Z-vector solve from the geometric machinery): for a random symmetric one-electron perturbation h^CCM → h^CCM + λV, the analytic Tr[P_tot V] matches the finite-difference dE_MP2/dλ to 2.8e-11 at (2,1,1) (the relaxation correction itself is 7.7e-6 – captured to 5 significant digits past it).

2. Geometric gate at (3,1,1): run_ccm_mp2_gradient vs the FD MP2 gradient (ccm_numerical_gradient with a tight-SCF MP2 runner) shows textbook O(h²) truncation – 4.3e-6 → 1.1e-6 → 2.7e-7 for h = 2e-3, 1e-3, 5e-4 – i.e. the analytic value is exact and the residual is pure FD truncation.

3. Isolated molecular limit, (1,1,1) 80-bohr box: matches a finite difference of vibe-qc’s molecular run_rhf + run_mp2 total energy (an independent implementation – C++ MP2 on canonical AO ERIs, no CCM code shared) to 3.6e-9 at h = 1e-4, the FD floor.

4. Sum rule: Σ_β dE_MP2/dR_β = 0 to 6.7e-16.

Caveat – near-degenerate SCF branch pairs break the FD reference, not the analytic force. At (2,1,1) the H₂-chain CCM-RHF ground state is translation-symmetry-broken (‖T P Tᵀ − P‖ ≈ 6e-3 for the one-cell translation T), and the two symmetry-partner solutions are split by only ~8e-12 Ha (the folded integrals’ known residual cyclic-invariance imperfection lifts the exact degeneracy). Displaced-geometry SCFs then land on either branch unpredictably; E_HF is branch-insensitive (~1e-11) but E₂ differs by ~8e-6 between branches, so the FD MP2 gradient at (2,1,1) is erratic while the FD HF gradient is clean. The analytic gradient is well-defined on the branch the reference SCF converged to. Same class as the UKS spin-branch FD caveat of §1.7: gate FD comparisons on a system with a unique SCF solution (the (3,1,1) chain here).

(1.8.4) Open-shell UMP2. The unrestricted gradient follows the same template with spin-resolved blocks (occupation 1 per spin orbital throughout; σ ∈ {α,β}, σ̄ the opposite spin). Amplitudes per the UMP2-CCM energy (ump2.py):

t^σσ_ijab = ḡ^σσ_iajb / D^σσ ,   ḡ^σσ_iajb = (ia|jb)^σσ − (ib|ja)^σσ ,
t^αβ_ijab = g^αβ_iajb / D^αβ     (i,a ∈ α; j,b ∈ β),
E₂ = ¼ Σ t^αα ḡ^αα + ¼ Σ t^ββ ḡ^ββ + Σ t^αβ g^αβ .                              (1.8j)

The intermediates absorb the spin factors (same-spin coefficient 1, cross-spin 2):

d^σ_oo = −½ t^σσ t^σσ − t^σσ̄ t^σσ̄ ,     d^σ_vv mirror (+) ,
X^σ_pi = Σ t^σσ ḡ^σσ_pajb + 2 Σ t^σσ̄ g^σσ̄_pajb ,   Y^σ_pa likewise on (ip|jb),
L^σ_ai = X^σ_ai − Y^σ_ia + 2 (C^σᵀ G^σ[D^α_ur, D^β_ur] C^σ)_ai ,                 (1.8k)

with the spin-resolved Fock response G^σ[X^α,X^β] = J[X^α+X^β] − K[X^σ]. The Z-vector is the spin-coupled CCM-CPHF solve on the stacked ov space,

M [z^α; z^β] = [L^α; L^β] ,
M^σσ_{ai,bj} = (ε^σ_a−ε^σ_i)δδ + 2(ai|bj)^σσ − (ab|ij)^σσ − (aj|ib)^σσ ,
M^σσ̄_{ai,bj} = 2 (ai|bj)^σσ̄ ,                                                   (1.8l)

and the assembled densities mirror (1.8f)–(1.8h) per spin: P^Δσ with ov block −½z^σ; W^Δσ_oo = ¼(X^σ+X^σᵀ) + ½d^σ_oo(ε+ε) + (C^σᵀ G^σ[P^Δα,P^Δβ] C^σ)_oo, W^Δσ_vv = ¼(Y^σ+Y^σᵀ) + ½d^σ_vv(ε+ε), W^Δσ_ov = ½Y^σ_ia − ½z^σ_ai ε^σ_i; Γ_ns = 2t^αα⊗(C^α)⁴ + 2t^ββ⊗(C^β)⁴ + 4t^αβ⊗(C^α)²(C^β)²; Γ_sep = 2P^Δ_tot⊗P_tot − 2P^Δα⊗P^α − 2P^Δβ⊗P^β (exchange slot). The final contraction is (1.8i) with the UHF Γ_HF of (1.7) and W = W^α+W^β+W^Δα+W^Δβ. Every block’s closed-shell reduction reproduces the RMP2 set exactly (t^αβ = t, t^σσ = t − t_swap, X^σ = 2X_cs, z^σ = ½z_cs, …), which is also pinned numerically below.

Validation (UMP2). (1) Closed-shell reduction: the UMP2 pipeline on the closed-shell H₂ chain (3,1,1) equals the RMP2 gradient to 6.2e-15 (machine ε; every spin factor pinned). (2) Response identity on the open-shell HeH chain (3,1,1) (doublet, 9 e⁻): 1.7e-10. (3) Geometric gate (HeH (3,1,1)): O(h²): 4.9e-6 → 1.2e-6 → 3.1e-7 for h = 2e-3, 1e-3, 5e-4. (4) Isolated limit (HeH/6-31G, 80-bohr box): matches the FD of vibe-qc’s molecular run_uhf + run_ump2 to 1.1e-9 at h = 1e-4 (the FD floor). (5) Sum rule: 1.6e-14.

Open-shell UMP2 is implemented and validated. CCSD relaxed-density forces are next (§1.9).

1.9 The CCSD relaxed-density gradient (post-HF forces)

CCSD closes the post-HF forces line. The energy (ccsd.py) is a spin-orbital CCSD (Stanton-Gauss-Watts-Bartlett 1991) evaluated on the canonical CCM-RHF orbitals and the WSSC-folded four-center; every MO integral is a folded quantity, so by §1.2 the fold commutes with the entire relaxed-density machinery (Λ, the relaxed 1-/2-particle densities, and the CCM-CPHF orbital response). The gradient is again the HF contraction pipeline (1.1) verbatim, with the CCSD relaxed densities:

dE_CCSD/dR =  Σ (P_SCF+P^Δ) ∂h^CCM/∂R  +  ½ Σ (Γ_HF+Γ_sep+Γ_ns) ∂(μν|λσ)^CCM/∂R
           −  Σ (W_HF+W^Δ) ∂S^CCM/∂R   +  ∂V_nn^CCM/∂R .                          (1.9a)

(1.9.1) The CCSD Lagrangian and its densities. With the converged amplitudes T = {t₁,t₂} and the de-excitation multipliers Λ = {λ₁,λ₂} (Handy-Schaefer 1984), the CCSD gradient is the relaxed-density gradient of Scheiner, Scuseria, Lee, Rice & Schaefer, J. Chem. Phys. 87, 5361 (1987) and Gauss & Cremer, Chem. Phys. Lett. 138, 131 (1987). Every ingredient of (1.9a) is a partial derivative of the CCSD correlation Lagrangian

L(T,Λ; f, W) = E_corr(T; f, W) + ⟨Λ, Ω(T; f, W)⟩ ,                               (1.9b)

(Ω the projected T-amplitude residual, f/W the folded 1-/2-electron MO integrals), taken at the stationary point (∂L/∂T = ∂L/∂Λ = 0):

• the response 1-PDM γ = ∂L/∂f (blocks d_oo, d_vv, d_ov = λ₁, d_vo);

• the non-separable 2-PDM (cumulant) Γ_ns = ∂L/∂W;

• the generalized Fock GF = ∂L/∂U (orbital variation C → C(1+U) with the reference Fock rebuilt from h_core, so the G[P_ref] reference-density response under the rotation is included).

(1.9.2) The Z-vector (CCM-CPHF). The orbital-rotation gradient is the antisymmetric ov block of the generalized Fock, X_ai = GF_ai − GF_ia, and the orbital relaxation solves the same CCM-CPHF equation as MP2/HF (Eq. 1.8e): M z = X, M_{ai,bj} = (ε_a−ε_i)δδ + 4(ai|bj)^CCM − (ab|ij)^CCM − (aj|ib)^CCM. (The X_GF computed from a rigid Fock transform misses the reference response; rebuilding f from h_core supplies it, equivalently the CCSD analog of MP2’s +4G^CCM[γ] term, verified numerically to 1.1e-9.) The relaxed 1-PDM is then P^Δ = C[γ + (−½z in ov/vo)]Cᵀ, P_tot = P_SCF + P^Δ.

(1.9.3) Assembled 2-PDM and energy-weighted density. As for MP2, the total 2-PDM splits Γ_tot = Γ_HF + Γ_sep + Γ_ns with Γ_sep = 2P^Δ⊗P − P^Δ⊗P_swap (separable coupling of P^Δ with the reference) and Γ_ns the CCSD cumulant. The energy-weighted (Pulay) density is W_tot = W_HF + W^Δ with

W^Δ = ¼ C(GF + GFᵀ)Cᵀ ,                                                          (1.9c)

which is −dE_corr/dS^CCM via the metric re-orthonormalization U^S = −½ Cᵀ(∂S)C (the generalized Fock already carries the orbital-relaxation weighting, so no separate z·ε block appears, unlike MP2’s explicit W^Δ_ov).

(1.9.4) Implementation: the robust complex-step route. The CCSD energy is a dense spin-orbital n_so⁶ path (small / 1-D / 2-D clusters only); the gradient matches that regime and takes the most provably correct route to the relaxed densities. Because L (1.9b) is a polynomial in the integrals, the complex-step derivative Im L(x + iε)/ε is exact to machine precision (Squire & Trapp, SIAM Rev. 40, 110 (1998)), so γ = ∂L/∂f, Γ_ns = ∂L/∂W, and GF = ∂L/∂U are obtained without any hand-derived CC density algebra, and Λ from the exact complex-step transpose-Jacobian of Ω (no hand-coded Λ equation). The final AO contractions reuse the identical L₂ᵀ/L_Vᵀ/L₄ᵀ cores and the same CCM-CPHF A. This is O(n_so²)/O(n_so⁴) Lagrangian evaluations, affordable for the validation clusters and philosophically consistent with the dense CCSD energy; an iterative Λ solver with closed-form densities is the scalable-kernel follow-on.

Validation (H₂ chain STO-3G, CCSD without (T)). (1) Relaxed-1-PDM response identity (h^CCM → h^CCM + λV, re-solving reference + CCSD): Tr[P_tot V] vs FD dE_CCSD/dλ to 1.9e-9; this isolates γ + the Z-vector orbital relaxation. (2) 2-PDM response ∂E_corr/∂W: Σ Γ_ns V₄ to 4.3e-12 (frozen-orbital). (3) W response ∂E/∂S: −Tr[W_tot V_S] to 9.2e-9. (4) Geometric gate (3,1,1): run_ccm_ccsd_gradient vs FD, textbook O(h²): 4.4e-6 → 1.1e-6 → 3.4e-7 for h = 2e-3, 1e-3, 5e-4. (5) Isolated (1,1,1) limit: matches the FD of vibe-qc’s molecular run_ccsd (C++ spin-adapted CCSD, no shared code) to 3.6e-9 at h = 1e-4. (6) Sum rule: 2.4e-15. The same (2,1,1) branch-pair FD caveat as MP2 (§1.8) applies; gate on (3,1,1).

Open-shell (UCCSD) and the (T) correction to the forces follow the same relaxed-density template; they remain follow-on work.

---

2. Implementation map

Target: python/vibeqc/periodic/ccm/gradient_analytic.py, run_ccm_rhf_gradient / run_ccm_uhf_gradient / run_ccm_rks_gradient / run_ccm_uks_gradient / run_ccm_mp2_gradient / run_ccm_ump2_gradient / run_ccm_ccsd_gradient (ccm, *, method="aiccm2026dev-a", ...) -> (n_basis_atoms, 3), exported from vibeqc.periodic.ccm. The four term cores (kinetic, nuclear attraction, overlap Pulay, two-electron) take explicit P/W/Γ arrays, so the RHF/UHF/KS/MP2/CCSD drivers differ only in how they build those (§1.7, §1.8, §1.9). Build order (each step finite-difference-checkable against ccm_numerical_gradient; the gradient is all-or-nothing, so each term is gated in isolation before the total):

1. ∂V_nn^CCM/∂R (§1.5) – direct pair-list differentiation; gate vs FD of ccm_nuclear_repulsion.

2. 1-e Σ P ∂h^CCM/∂R (§1.4 L₂ᵀ for T; L_Vᵀ for V incl. the HF-nucleus term).

3. Overlap Pulay −Σ W ∂S^CCM/∂R (§1.4 L₂ᵀ).

4. 2-e Σ Γ ∂(μν|λσ)^CCM/∂R (§1.4 L₄ᵀ + two_electron_gradient_casscf).

5. Total vs ccm_numerical_gradient on the H₂ chain (2,1,1)/(4,1,1) STO-3G to ~1e-6 Ha/bohr; per-cell forces sum ≈ 0 (Theorem 1.6.1).

Molecular primitives (no new C++): nuclear_repulsion_gradient, one_electron_gradient_contribution, overlap_gradient_contribution, two_electron_gradient_casscf, compute_external_charge_density_gradient (bindings.cpp 2753-2870), all driven on the padded Molecule+BasisSet with the adjoint-scattered P_pad, W_pad, Γ_pad, then image-summed via pad_atom_of + atom_beta (§1.3).

Two implementation subtleties worth recording, both load-bearing for exactness:

1. Kinetic isolation. T^CCM and V^CCM carry different folds, so they cannot be driven by the combined T+V one-electron primitive together. ∂T is obtained from one_electron_gradient_contribution(..., effective_charges=0), which zeroes the V_ne derivative kernel and leaves the pure kinetic gradient.

2. 8-fold symmetrization of Γ_pad. two_electron_gradient_casscf sums over unique shell quartets ×degeneracy and assumes the supplied AO 2-RDM carries the full 8-fold ERI permutation symmetry. The L₄ᵀ scatter places Γ̄ into distinct padded slots, breaking μ↔ν/λ↔σ symmetry on the padded layout; since ∂eri_pad/∂R is itself 8-fold symmetric, Σ Γ_pad ∂eri = Σ sym₈(Γ_pad) ∂eri, so Γ_pad is 8-fold symmetrized (then scaled by ½) before the call.

Validation (RHF, implemented)

python/vibeqc/periodic/ccm/gradient_analytic.py · tests/test_ccm_gradient_analytic.py. Three independent checks on STO-3G:

• Per-term, isolated. Each of the five terms (V_nn, T, V_ne, overlap Pulay, two-electron) matches a finite difference of its own fixed-density contraction to ≈ 1e-9 Ha/bohr on the H₂ chain (2,1,1)-(4,1,1).

• Total vs the full-energy FD gate. run_ccm_rhf_gradient matches ccm_numerical_gradient with textbook O(h²) convergence – the residual halves by 4× per halving of h (4.25e-6 → 1.06e-6 → 2.66e-7 → 6.65e-8 → 1.06e-8 for h = 2e-3 … 1e-4 on (3,1,1)), i.e. the analytic value is exact and the gap is pure finite-difference truncation. (Requires the tight SCF of the remark above.)

• Isolated molecular limit (analytic vs analytic). In a (1,1,1) 80-bohr box the Γ-CCM analytic gradient reproduces vibe-qc’s molecular analytic RHF gradient (compute_gradient) to machine precision (5.6e-16) – no finite difference involved. The Γ-CCM force reduces exactly to the molecular force in the molecular limit.

MP2 validation (response identity, geometric O(h²) gate, independent molecular cross-check, sum rule, and the (2,1,1) branch-pair FD caveat): §1.8.

---

References

• T. Helgaker, P. Jørgensen, J. Olsen, Molecular Electronic-Structure Theory (Wiley, 2000) – analytic gradient theory, Ch. 10 and Eq. 12.5.7 (the general AO 2-RDM ERI-gradient contraction).

• P. Pulay, “Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules. I.”, Mol. Phys. 17, 197 (1969) – the Hellmann-Feynman/Pulay split and the energy-weighted density.

• J. A. Pople, R. Krishnan, H. B. Schlegel, J. S. Binkley, “Derivative studies in Hartree-Fock and Møller-Plesset theories”, Int. J. Quantum Chem. Symp. 13, 225 (1979) – the MP2 gradient (§1.8).

• N. C. Handy, H. F. Schaefer III, “On the evaluation of analytic energy derivatives for correlated wave functions”, J. Chem. Phys. 81, 5031 (1984), doi:10.1063/1.447489 – the Z-vector elimination of the orbital response (§1.8, Eq. 1.8e; §1.9).

• A. C. Scheiner, G. E. Scuseria, J. E. Rice, T. J. Lee, H. F. Schaefer III, “Analytic evaluation of energy gradients for the single and double excitation coupled cluster (CCSD) wave function”, J. Chem. Phys. 87, 5361 (1987) – the CCSD relaxed-density gradient (§1.9).

• J. Gauss, D. Cremer, “Analytical evaluation of energy gradients in quadratic configuration interaction theory”, Chem. Phys. Lett. 138, 131 (1987) – the CC/QCI relaxed-density gradient (§1.9).

• J. F. Stanton, J. Gauss, J. D. Watts, R. J. Bartlett, “A direct product decomposition approach for symmetry exploitation in many-body methods. I. Energy calculations”, J. Chem. Phys. 94, 4334 (1991) – the spin-orbital CCSD amplitude equations (§1.9, the folded CCSD kernel).

• W. Squire, G. Trapp, “Using complex variables to estimate derivatives of real functions”, SIAM Rev. 40, 110 (1998) – the exact complex-step differentiation used for the CCSD relaxed densities (§1.9.4).

• M. F. Peintinger, T. Bredow, “The Cyclic Cluster Model at Hartree-Fock Level”, J. Comput. Chem. 35, 839 (2014), doi:10.1002/jcc.23550 – the CCM/WSSC weighting (eqs 5-6, 12-13, 18).