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 Cyclic Cluster Model (≡ SCM-Γ; 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 imposes translational symmetry not by Bloch summation but by Wigner–Seitz-supercell (WSSC) integral weighting. 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 — adds the exchange-correlation Pulay term Σ_μν (∂E_xc/∂P_μν)|... ∂(folded grid quantities)/∂R and the weight-derivative of the XC grid; the HF-exchange fraction scales the Γ four-center term. Outlook.

• Post-HF (MP2/CCSD) — relaxed-density forces: replace P,W,Γ by the method’s relaxed one- and two-particle densities (and the orbital-response Z-vector), folded and scattered identically. Outlook; cf. the CASSCF general-Γ precedent (HANDOVER_CAS_GRADIENT.md).

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.

---

2. Implementation map

Target: python/vibeqc/periodic/ccm/gradient_analytic.py, run_ccm_rhf_gradient / run_ccm_uhf_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 and UHF drivers differ only in how they build those (§1.7). 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.

---

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.

• 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).