CCResPy: Coupled-Cluster Response in Python

• CCResPy is an open‐source Python/NumPy framework implementing coupled‐cluster (CCSD) response theory for both molecular and 1D periodic systems.
• It leverages collective indices and tensor-based formulations to compute frequency-dependent electric dipole polarizabilities, addressing challenges in k-space convergence and dipole operator ambiguities.
• Designed as a method-development platform, CCResPy operates on externally generated mean-field data to facilitate benchmark-quality optical response calculations with transparent, extensible code.

Searching arXiv for CCResPy and closely related periodic CCSD response work. CCResPy, short for Coupled-Cluster Response in Python, is an open-source Python/NumPy implementation of coupled-cluster response theory centered on coupled cluster with single and double excitations (CCSD), its left-eigenvector equations, and linear-response CCSD for general one-body perturbations. In its current documented form, it supports both molecular systems and 1D periodic systems with periodic boundary conditions (PBCs), and has been extended to compute the frequency-dependent electric dipole-electric dipole polarizability tensor in the length gauge for 1D periodic chains. The code is positioned as a specialized method-development platform rather than an all-in-one quantum-chemistry package: it assumes an external mean-field calculation and operates on imported one- and two-electron quantities, dipole integrals, and PBC data (Caricato et al., 2 Aug 2025).

1. Scope and conceptual role

CCResPy is designed around a deliberately narrow but technically demanding target: high-accuracy correlated electronic-structure calculations based on CCSD and CCSD response theory. Its stated capabilities include ground-state CCSD energies and amplitudes, solution of the left-hand $\Lambda$ equations, and linear-response CCSD for one-body perturbations, with the present implementation using these ingredients to evaluate the frequency-dependent electric dipole-electric dipole polarizability. The same framework is used for molecules and for 1D periodic systems through a collective-index formalism that combines orbital and k-point labels (Caricato et al., 2 Aug 2025).

Aspect	CCResPy implementation	Role
Language	Python + NumPy	Tensor-based method development
Wavefunction level	CCSD	Ground state and response
System classes	Molecular and 1D PBC	Unified tensor equations
Response capability	LR-CCSD for one-body perturbations	Here applied to $\boldsymbol{\alpha}(\omega)$
External dependencies	Mean-field and integrals generated elsewhere	CCResPy does not perform SCF

The code is explicitly not described as a full electronic-structure environment. It does not generate integrals, perform SCF, or carry out geometry optimization. Instead, the reference Hartree-Fock calculation is done externally, in the reported work with the GAUSSIAN development version, and the necessary data are extracted and passed to CCResPy. This places it in the category of a method-development and high-accuracy prototype code, complementary to production packages such as Gaussian, CRYSTAL, and VASP, especially for benchmark-quality optical-response calculations in solids where DFT and TDDFT are usually the default (Caricato et al., 2 Aug 2025).

2. Coupled-cluster and linear-response formulation

The ground-state theory follows the standard single-reference CCSD ansatz with Hartree-Fock reference determinant $|\Phi_0\rangle$,

$$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$

with single and double excitation amplitudes $t_I^A$ and $t_{IJ}^{AB}$. For periodic systems, CCResPy introduces collective indices $P \equiv p k_p$, so that molecular and periodic equations can be written in the same algebraic form. In this notation, the periodic CCSD correlation energy per unit cell is

$$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$

with $N_k$ the number of k-points, $f_{PQ}$ Fock-matrix elements in the crystal-orbital/k-point basis, and $\boldsymbol{\alpha}(\omega)$0 antisymmetrized two-electron integrals. The dressed doubles amplitudes are

$\boldsymbol{\alpha}(\omega)$1

$\boldsymbol{\alpha}(\omega)$2

The amplitude equations are written in the usual residual-equals-zero form with orbital-energy denominators

$\boldsymbol{\alpha}(\omega)$3

Their working expressions follow Gauss-Stanton-Bartlett style CCSD gradients and intermediates, adapted to PBC in the sense described in the source work. Intermediates such as $\boldsymbol{\alpha}(\omega)$4, $\boldsymbol{\alpha}(\omega)$5, $\boldsymbol{\alpha}(\omega)$6, and $\boldsymbol{\alpha}(\omega)$7-type tensors are computed once per CCSD iteration and then reused in the tensor contractions (Caricato et al., 2 Aug 2025).

For response properties, CCResPy also solves the $\boldsymbol{\alpha}(\omega)$8-CCSD equations for the de-excitation amplitudes $\boldsymbol{\alpha}(\omega)$9 and $|\Phi_0\rangle$0. These left eigenvectors enter the symmetric linear-response formalism and the construction of $|\Phi_0\rangle$1-type intermediates and $|\Phi_0\rangle$2-type one-particle transition-density-like tensors. For a one-electron perturbation $|\Phi_0\rangle$3, the linear-response function is expressed as

$|\Phi_0\rangle$4

but CCResPy does not form the Jacobian $|\Phi_0\rangle$5 explicitly. Instead, it solves directly for perturbed amplitudes $|\Phi_0\rangle$6 from

$|\Phi_0\rangle$7

For electric-dipole response, the polarizability tensor is obtained as

$|\Phi_0\rangle$8

An important implementation choice is that the present symmetric LR formulation requires only perturbed $|\Phi_0\rangle$9-amplitudes. The alternative non-symmetric formulation involving perturbed $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$0-amplitudes $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$1 is explicitly noted as not yet implemented. This keeps the current response machinery aligned with the existing CCSD, $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$2-CCSD, and tensor-intermediate infrastructure (Caricato et al., 2 Aug 2025).

3. Periodic boundary conditions and k-space treatment

The periodic implementation targets 1D crystals: infinite chains with translation vector $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$3 along one spatial direction. The underlying mean-field data are generated with Gaussian-type orbitals and Bloch orbitals under PBCs. GAUSSIAN provides the crystal-orbital coefficients $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$4, AO-basis Fock and overlap matrices, AO-basis two-electron repulsion integrals, AO dipole integrals, and the real-space and reciprocal-space parameters $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$5 and $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$6. CCResPy reconstructs the k-point mesh internally from the reciprocal-space sampling parameter and transforms all AO quantities to the crystal-orbital/k-point basis (Caricato et al., 2 Aug 2025).

Momentum conservation is enforced through the usual k-selection rules,

$$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$7

with integers $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$8 and $$|\Psi_{\text{CC}}\rangle = e^{\hat T}|\Phi_0\rangle, \qquad \hat T = \hat T_1 + \hat T_2,$$9. In consequence, only $t_I^A$0 independent blocks exist for one-electron quantities and $t_I^A$1 for two-electron quantities. However, CCResPy deliberately stores the tensors with all k-indices present, as if the dimensions were $t_I^A$2 and $t_I^A$3, while setting disallowed elements to zero. This is less efficient in storage and floating-point work, but it permits uniform tensor shapes and a single code path for molecules and periodic systems. The reported implementation treats this as a central design decision: efficiency is sacrificed in favor of correctness, code transparency, and extensibility.

The thermodynamic limit is approached by increasing both the number of replica cells used in real-space integral generation and the number of k-points used to sample the 1D Brillouin zone. No special extrapolation module is built into CCResPy. Instead, the reported study performs direct k-mesh convergence analyses, comparing periodic CCSD and LR-CCSD results with large finite-cluster calculations. A major numerical conclusion is that correlation energies converge much faster with k-point sampling than response properties. In the tested systems, very coarse meshes already yield CCSD correlation energies close to cluster references, whereas polarizabilities, especially along the periodic direction, require substantially denser k-point sampling for semi-quantitative agreement (Caricato et al., 2 Aug 2025).

4. Electric dipole operator and the length-gauge ambiguity

A central conceptual issue in the periodic implementation concerns the electric dipole operator in the length gauge. For finite molecules, the electric dipole operator is a straightforward one-body operator. For periodic solids, the position operator is ill-defined under PBC, and the modern theory of polarization replaces a naive position expectation value by a Berry-phase formulation for static polarization. In the periodic Gaussian framework used here, length-gauge dipole integrals in the Bloch-orbital basis involve k-derivatives of Bloch orbitals,

$t_I^A$4

with $t_I^A$5 analogous to a coupled-perturbed orbital-response matrix (Caricato et al., 2 Aug 2025).

The reported difficulty is that $t_I^A$6 lacks certain diagonal contributions, described as “missing integers,” associated with the arbitrary k-dependent phases of Bloch orbitals. In Hartree-Fock and DFT linear-response polarizabilities in the length gauge, these omitted pieces do not affect the result because only off-diagonal dipole blocks enter. In CCSD linear response, however, the full dipole matrix appears in the working expressions for the property coupling vector and the final response function. Consequently, the ambiguity propagates directly into the CCSD polarizability calculation.

The observed numerical effect depends on the system. For the H$t_I^A$7 chain, the impact is reported to be very small compared with ordinary k-point convergence errors. For LiH chains, deviations of several atomic units in the longitudinal component $t_I^A$8 remain even at large k-meshes and are attributed primarily to the missing-integer issue. The source work emphasizes that this is not a failure of CCSD itself but a problem in the operator definition supplied to the response equations. If a fully phase-consistent $t_I^A$9 were available, the CC equations implemented in CCResPy would not need to change (Caricato et al., 2 Aug 2025).

The velocity gauge is identified as an alternative route, since it expresses the dipole perturbation in terms of the momentum operator and avoids the $t_{IJ}^{AB}$0 construction. At finite basis and finite truncation level, however, length-gauge and velocity-gauge formulations are not numerically identical. A plausible implication is that gauge comparison will be necessary not only as a diagnostic of the operator issue but also as a measure of basis-set and truncation error in future periodic CC response work.

5. Software architecture and numerical strategy

CCResPy is written in Python, with tensor algebra handled through NumPy and einsum contractions. BLAS and LAPACK are used indirectly through NumPy, and shared-memory parallelism is inherited from the linked linear-algebra backend. The implementation is intentionally equation-like: tensor storage and index order are chosen so that the code mirrors the mathematical expressions as closely as possible. Quantities such as $t_{IJ}^{AB}$1, $t_{IJ}^{AB}$2, and $t_{IJ}^{AB}$3 are stored in array layouts that preserve the correspondence between code and working equations (Caricato et al., 2 Aug 2025).

Periodic quantities are generally complex, so the same infrastructure supports real molecular tensors and complex periodic tensors. Reference data are imported rather than generated internally. The workflow uses GAUSSIAN PBC Hartree-Fock calculations together with gauopen to extract binary data into text or NumPy arrays. The imported information includes crystal-orbital coefficients, AO-basis two-electron integrals, AO Fock and overlap matrices, AO dipole integrals, and the PBC parameters required to reconstruct the reciprocal-space description.

Memory management is a prominent part of the implementation. All arrays are kept in memory except the largest $t_{IJ}^{AB}$4 arrays, such as virtual-space two-electron tensors and related $t_{IJ}^{AB}$5-type intermediates. These are written to binary files, loaded only when needed for a contraction, and then released. For particularly large contractions, CCResPy checks available memory and either performs the full contraction or falls back to sliced contractions over an outer index, reducing peak memory at the cost of additional Python-level loops. The code uses psutil, tracemalloc, and, on Linux, resource to query and enforce a global memory cap.

All amplitude equations—ground-state $t_{IJ}^{AB}$6, left-hand $t_{IJ}^{AB}$7, and perturbed $t_{IJ}^{AB}$8—are solved iteratively with Jacobi-style fixed-point updates. Convergence thresholds are reported as $t_{IJ}^{AB}$9 Ha for the energy, together with RMS and MAX thresholds on amplitude updates of $P \equiv p k_p$0 and $P \equiv p k_p$1, respectively. A pseudo-energy analogous to the CCSD energy expression is used to monitor convergence of $P \equiv p k_p$2 and $P \equiv p k_p$3. Pulay DIIS acceleration is employed, with the error vector defined as the difference between successive iterates. In the periodic case, the DIIS overlap matrix is complex Hermitian rather than real symmetric, and old amplitude sets are stored to disk and reloaded when needed. This combination of transparent tensor algebra, explicit iterative solvers, and memory-aware disk staging is characteristic of CCResPy’s role as a research prototype (Caricato et al., 2 Aug 2025).

6. Benchmark systems and numerical behavior

The reported numerical tests use CCResPy on three model classes: an H$P \equiv p k_p$4 chain, a LiH chain with fixed geometry, and an optimized LiH chain. In all cases, the comparisons involve isolated molecules, large finite clusters, and periodic calculations with varying reciprocal-space sampling, allowing separate examination of intermolecular effects, finite-size effects, and k-mesh convergence (Caricato et al., 2 Aug 2025).

For the H$P \equiv p k_p$5 chain, the geometry is $P \equiv p k_p$6 Å with translation vector $P \equiv p k_p$7 Å and a 3-21G basis, giving 2 occupied and 6 virtual orbitals per cell. Only $P \equiv p k_p$8, the polarizability along the chain direction, is nonzero. Moving from the isolated H$P \equiv p k_p$9 molecule to the infinite chain increases the Hartree-Fock energy per cell by about $$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$0 mHa, while correlation partially compensates this increase by roughly half. The CCSD and MP2 correlation energies from PBC converge rapidly with k-mesh; for sufficiently dense meshes the differences from cluster references fall to the $$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$1–$$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$2 Ha range. The polarizability behaves differently: the chain $$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$3 is about $$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$4–$$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$5 a.u. larger than the molecular value, but coarse PBC sampling can qualitatively misrepresent this effect. In particular, PBC(5) underestimates $$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$6 by about $$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$7 a.u., whereas PBC(21) and larger reduce the discrepancy to about $$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$8–$$E_\text{corr} = \frac{1}{N_k}\, t_I^{A} f_{IA} + \frac{1}{4 N_k^3}\, \tau_{IJ}^{AB} \langle IJ || AB \rangle,$$9 a.u.

For the LiH chain with bond length $N_k$0 Å and $N_k$1 Å in STO-3G, all electrons are correlated, yielding 4 occupied and 8 virtual orbitals per cell. Two independent polarizability components occur: $N_k$2 along the chain and $N_k$3 perpendicular to it. The cluster calculations show an inversion relative to the isolated molecule: in the molecule, $N_k$4, whereas in the chain, $N_k$5. Thus the longitudinal polarizability is strongly enhanced by chain formation, while the transverse polarizability decreases modestly. CCResPy reproduces the perpendicular components well even on relatively small PBC meshes, with remaining differences likely due more to finite-cluster artifacts than to PBC deficiencies. The longitudinal component is more problematic: even at PBC(17)–PBC(21), periodic values differ from cluster references by about $N_k$6–$N_k$7 a.u., approximately $N_k$8 of the full molecule-to-chain increase. The source work attributes this persistent discrepancy mainly to the missing-integer problem in the length-gauge dipole operator rather than to inadequacy of the CCSD response formalism itself.

The optimized LiH chain provides a stricter test because the repeat distance shortens to $N_k$9 Å and the Li–H bond length changes to $f_{PQ}$0 Å, increasing intercell interaction. The Hartree-Fock energy per cell is lowered by about $f_{PQ}$1 mHa relative to the previous LiH-chain geometry, while the isolated molecule changes by only about $f_{PQ}$2 mHa. Correlation energies change little, which is reported as evidence that, in this minimal basis, correlation is mostly local. The reference $f_{PQ}$3 becomes $f_{PQ}$4–$f_{PQ}$5 a.u. smaller than in the previous geometry, whereas molecular values remain nearly unchanged. Once again, the perpendicular polarizabilities are reproduced well at sufficiently large k-meshes, but the longitudinal component differs from cluster results by about $f_{PQ}$6–$f_{PQ}$7 a.u. even at PBC(25), reinforcing the interpretation that the dominant limitation is the length-gauge dipole operator rather than geometry-dependent electronic structure.

7. Significance, limitations, and extension paths

The principal significance of CCResPy in its current documented form is that it realizes, within one explicit implementation, periodic CCSD ground state, $f_{PQ}$8-CCSD, perturbed CCSD amplitudes at complex frequencies, and symmetric CC linear response for 1D periodic systems. The source work identifies this as the first implementation of frequency-dependent electric dipole-electric dipole polarizability for 1D periodic systems at the CCSD level with PBCs. In methodological terms, it extends coupled-cluster response theory—long established for molecular optical rotation, Raman optical activity, and higher-order response—into a domain where DFT and TDDFT have largely dominated practical computations (Caricato et al., 2 Aug 2025).

The implementation also makes clear what CCResPy is not yet. It does not exploit restricted closed-shell spin symmetry, reduced k-space tensor representations, local correlation, or low-rank factorization. All tensors are dense and global. Only 1D periodicity is implemented. The code is therefore best understood as a research prototype and development platform rather than a production workhorse for very large systems or very dense k-meshes. The major scientific limitation highlighted by the benchmarks is the electric-dipole operator in the periodic length gauge, especially for longitudinal response in strongly polarizable chains.

The future directions stated for the code are correspondingly method-oriented. They include more efficient periodic implementations with explicit momentum-conserving k-summations, persistent disk-based strategies for large intermediates, a velocity-gauge formulation of polarizability to avoid the missing-integer issue, higher-order response functions and additional operators such as magnetic dipole and electric quadrupole, and extension to 2D and 3D periodicity. Because the present code is organized around collective indices and tensor intermediates rather than separate molecular and periodic engines, these extensions are presented as additions to the intermediate machinery rather than a wholesale redesign of the architecture. This suggests that CCResPy’s lasting importance may lie not only in the specific 1D polarizability implementation already reported, but also in its role as an explicit, extensible reference implementation for periodic coupled-cluster response theory.