Coupled-cluster theory


The MP2 and coupled-cluster functionalities of PySCF are similar. See also Second-order Møller–Plesset perturbation theory (MP2).


Coupled-cluster (CC) theory is a post-Hartree-Fock method capable of describing electron correlation in the ground state. It is size extensive but not variational. PySCF has extensive support for CC calculations with single and double excitations (CCSD). It can also include a perturbative treatment of triple excitations (CCSD(T)), which is a very accurate method for single-reference quantum chemistry. CC calculations can be performed with or without density fitting, depending on the initial SCF calculation. Correlated excited states are accessible through the equation-of-motion (EOM) CCSD framework, described below.

A minimal example of a CCSD and CCSD(T) calculation is as follows:

from pyscf import gto, scf, cc
mol = gto.M(
    atom = 'H 0 0 0; F 0 0 1.1',  # in Angstrom
    basis = 'ccpvdz',
    symmetry = True,
mf = scf.HF(mol).run()
# Note that the line following these comments could be replaced by
# mycc = cc.CCSD(mf)
# mycc.kernel()
mycc = cc.CCSD(mf).run()
print('CCSD total energy', mycc.e_tot)
et = mycc.ccsd_t()
print('CCSD(T) total energy', mycc.e_tot + et)

Spin symmetry

The CC module in PySCF supports a number of reference wave functions with broken spin symmetry. In particular, CC can be performed with a spin-restricted, spin-unrestricted, and general (spin-mixed) Hartree-Fock solution, leading to the RCCSD, UCCSD, and GCCSD methods.

The module-level cc.CCSD(mf) constructor can infer the correct method based on the level of symmetry-breaking in the mean-field argument. For more explicit control or inspection, the respective classes and functions can be found in (restricted with real orbitals), (restricted with potentially complex orbitals), (unrestricted), and (general).

For example, a spin-unrestricted calculation on triplet oxygen can be performed as follows:

from pyscf import gto, scf, cc
mol = gto.M(
    atom = 'O 0 0 0; O 0 0 1.2',  # in Angstrom
    basis = 'ccpvdz',
    spin = 2
mf = scf.HF(mol).run() # this is UHF
mycc = cc.CCSD(mf).run() # this is UCCSD
print('UCCSD total energy = ', mycc.e_tot)


A number of properties are available at the CCSD level.

Unrelaxed 1- and 2-electron reduced density matrices can be calculated. They are returned in the MO basis:

dm1 = mycc.make_rdm1()
dm2 = mycc.make_rdm2()

Analytical nuclear gradients can be calculated:

mygrad = mycc.nuc_grad_method().run()

The CCSD Lambda equations can be solved:

l1, l2 = mycc.solve_lambda()

Frozen orbitals

By default, CCSD calculations in PySCF correlate all electrons in all available orbitals. To freeze the lowest-energy core orbitals, use the frozen keyword argument:

mycc = cc.CCSD(mf, frozen=2).run()

To freeze occupied and/or unoccupied orbitals with finer control, a list of 0-based orbital indices can be provided as the frozen keyword argument:

# freeze 2 core orbitals
mycc = cc.CCSD(mf, frozen=[0,1]).run()
# freeze 2 core orbitals and 3 unoccupied orbitals
mycc = cc.CCSD(mf, frozen=[0,1,16,17,18]).run()

Equation-of-motion coupled-cluster theory

EOM-CCSD can be used to calculate neutral excitation energies (EE-EOM-CCSD), spin-flip excitations (SF-EOM-CCSD), or charged excitations, i.e. ionization potentials (IP-EOM-CCSD) or electron affinities (EA-EOM-CCSD). The EOM functions return the requested number of eigenvalues and right-hand eigenvectors. For example:

e_ip, c_ip = mycc.ipccsd(nroots=1)
e_ea, c_ea = mycc.eaccsd(nroots=1)
e_ee, c_ee = mycc.eeccsd(nroots=1)
e_sf, c_sf = mycc.eomsf_ccsd(nroots=1)

The eecsd() function returns neutral excitations with all possible spin multiplicities. For closed-shell calculations (RHF and RCCSD), singlet and triplet excitations can be requested explicitly:

e_s, c_s = mycc.eomee_ccsd_singlet(nroots=1)
e_t, c_t = mycc.eomee_ccsd_triplet(nroots=1)

By default, PySCF calculates the nroots eigenvalues with the lowest energy, which may include states with dominant double-excitation character. To only calculate states with dominant single-excitation character, use the koopmans keyword argument:

e, c = mycc.eeccsd(nroots=3, koopmans=True)

An initial guess wavefunction may be provided, in which case PySCF will try to find the most similar EOM solution vector:

from import amplitudes_to_vector_ee
r1 = np.zeros((nocc,nvir))
r2 = np.zeros((nocc,nocc,nvir,nvir))
r1[occ_index,vir_index] = 1.0
myguess = amplitudes_to_vector_ee(r1,r2)
e_s, c_s = mycc.eomee_ccsd_singlet(nroots=1, guess=myguess)

Job control

Saving and restarting

To allow for future restarts, the SCF information and the CCSD DIIS information must be saved:

mf = scf.HF(mol)
mf.chkfile = 'hf.chk'

mycc = cc.CCSD(mf)
mycc.diis_file = 'ccdiis.h5'

To restart a CCSD calculation, first the molecule and SCF information must be restored:

mol = lib.chkfile.load_mol('hf.chk')
mf = scf.HF(mol)
mf.__dict__.update(lib.chkfile.load('hf.chk', 'scf'))

Next, the CCSD calculation can be restarted by using the previous CCSD amplitudes as the initial guess:

mycc = cc.CCSD(mf)
mycc.kernel(mycc.t1, mycc.t2)

Modifying DIIS

The parameters of the DIIS algorithm can be tuned in cases where convergence is difficult. To increase the size of the DIIS space:

mycc = cc.CCSD(mf)
mycc.diis_space = 10

By default, DIIS is activated on the first CCSD iteration. Sometimes it can be helpful to postpone the use of DIIS:

mycc = cc.CCSD(mf)
mycc.diis_start_cycle = 4

Integral-direct CCSD

In order to avoid large memory requirements, the default behavior in CCSD calculations is to store most two-electron integral tensors on disk. This leads to a potential I/O bottleneck. For medium-sized molecules, an integral-direct AO-driven implementation can be more efficient. The user must manually request an integral-direct CCSD calculation:

mycc = cc.CCSD(mf) = True
e_corr, t1, t2 = mycc.kernel()