CCSD calculation between two orbitals

I am trying to calculate the CCSD energies between two orbitals. Suppose I have a system,

from pyscf import gto, scf, ao2mo, cc
from pyscf.gto import mole

mol = gto.Mole()
mol.atom = '''
O  0.0000000, 0.0000000, 0.0000000
H  0.757606,  0.586006,  0.0000000
H -0.757606,  0.586006,  0.0000000'''
mol.basis = 'sto-3g'
mol.charge = 0
mol.spin = 0
mol.build()

# Perform SCF calculation
mf = scf.RHF(mol)
mf.kernel()

# Get MO coefficient matrix
mo_coeff = mf.mo_coeff

# Obtain the CCSD energy from canonical
mcc = cc.CCSD(mf)
ecc,t1,t2 = mcc.kernel()
print("CCSD correlation energy = ",mcc.e_corr)
print("Total CCSD energy = ",mcc.e_tot)


The correlation energy is -0.0494284459751351 Hartree. Now, this calculation takes all HF orbitals and occupancies as arguments to the CCSD function. If I know correctly, I can calculate the CCSD correlation energies between each pair of MO (of course, one orbital should be doubly occupied, and the other should be vacant) and sum them up to obtain the total correlation energy. To do this, I take a pair of orbitals (one from doubly occupied and another from vacant MO) and run the CCSD calculation for each pair. So, for water molecules with a sto-3g basis, the MO coefficient matrix is 7*7, and the occupancy is a numpy array [2.0,2.0,2.0,2.0,2.0,0.0,0.0]. So I obtained the MO columns using this (just for testing and without any beautiful coding)

import numpy as np
coeff_0 = mf.mo_coeff[:,0]; coeff_1 = mf.mo_coeff[:,1]; coeff_2 = mf.mo_coeff[:,2]; coeff_3 = mf.mo_coeff[:,3]
coeff_4 = mf.mo_coeff[:,4]; coeff_5 = mf.mo_coeff[:,5]; coeff_6 = mf.mo_coeff[:,6]


Then, I merged the pairs for occupied and virtual MOs. In the case of the MO coefficient matrix, columns 0-4 are occupied, and columns 5-6 are virtual. So, for each occupied column, there will be two virtual columns. I merged the columns as,

mo_coeff05 = np.column_stack((coeff_0,coeff_5));mo_coeff06 = np.column_stack((coeff_0,coeff_6))
mo_coeff15 = np.column_stack((coeff_1,coeff_5));mo_coeff16 = np.column_stack((coeff_1,coeff_6))
mo_coeff25 = np.column_stack((coeff_2,coeff_5));mo_coeff26 = np.column_stack((coeff_2,coeff_6))
mo_coeff35 = np.column_stack((coeff_3,coeff_5));mo_coeff36 = np.column_stack((coeff_3,coeff_6))
mo_coeff45 = np.column_stack((coeff_4,coeff_5));mo_coeff46 = np.column_stack((coeff_4,coeff_6))


And at last, I run the CCSD for each merged MO coefficient matrix.

mcc.mo_occ = np.array([2.0,0.0])
mcc.mo_coeff = mo_coeff05
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff06
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff15
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff16
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff25
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff26
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff35
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff36
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff45
mcc.kernel()[0]
mcc.mo_coeff = mo_coeff46
mcc.kernel()[0]


What I expected that if I sum all the correlation energies, I will get the total CCSD correlation energy, but some of the correlation energy are even higher than the total system.

E(CCSD) = -51.47188187443665  E_corr = -0.0003962543997611398
E(CCSD) = -51.47149476545477  E_corr = -9.145417883697846e-06
E(CCSD) = -6.483145580303841  E_corr = -1.059192232873223
E(CCSD) = -5.424905601367101  E_corr = -0.0009522539364820224
E(CCSD) = -2.90788145525101  E_corr = -0.002489550153102168
E(CCSD) = -5.339446749364194  E_corr = -2.434054844266286
E(CCSD) = -4.937559682283515  E_corr = -0.9363818461883874
E(CCSD) = -4.003003233876766  E_corr = -0.001825397781638726
E(CCSD) = -4.845351104114307  E_corr = -0.0003760041376774681
E(CCSD) = -4.84514723864436  E_corr = -0.0001721386677304669


Also, the SCF + CCSD energy defined by this, E(CCSD) = is very different.

Is my approach wrong?

• This question has one significant problem. Instead of using one pair of occupied and virtual orbitals, one must use two pairs of occupied and virtual orbitals. As the author of this question, can I modify the question significantly?
– Pro
Feb 16 at 17:04