# What's the best code for getting 1-particle and 2-particle density matrices with CCSD(T) in Python?

I have been using PySCF to calculate the 1-particle and 2-particle density matrix from ccsd(T) wavefunction using these modules in-built in PySCF make_rdm1() and make_rdm2(). However, the code takes ~2TB of memory for CCSD(T) calculation of only ~175 basis functions. So, I have been thinking to switch to other commercial/free codes like Gaussian, Psi4. Turbomole or orca for CCSD(T) calculations and then convert the output to obtain the density matrices (with python3). Should I have to start with MO coefficients? How can I do it? Also some codes do not provide hessian for CCSD(T).

Edit 1: Added PySCF Code. Edit 2: Modified PySCF code to generate 2PDM at CCSD(T) (no frozen core).

import numpy
import scipy.linalg
from pyscf import gto, scf
from pyscf import lib, tools
from pyscf import cc, ao2mo, grad
from gto import mole
from pyscf.cc import ccsd
from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda
from pyscf.cc import ccsd_t_slow as ccsd_t
from pyscf.cc import ccsd_t_rdm_slow as ccsd_t_rdm
from pyscf.cc.ccsd_t_rdm_slow import make_rdm1,make_rdm2

mol = gto.M(
verbose = 5,
atom = [
['O' ,  0.074642 , -0.000000 , -0.026389],
['H' , -0.103145 ,  0.000000 ,  0.930145],
['H' ,  0.991917 ,  0.000000 , -0.350692],
['H' , -0.664854 , -0.000000 , -0.658618]],
basis = '3-21g',charge=1,
cart=True
)
mf = scf.RHF(mol)
mf.conv_tol = 1e-13
mf.scf()
orbshape = mf.mo_coeff.shape

HFdm1 = mf.make_rdm1()
E1e = numpy.einsum('pq,qp', mf.get_hcore(), HFdm1)

with open('moldenmo', 'w') as f1:
for imo in range(orbshape[1]):
for j in  range(orbshape[1]):
f1.write('\n{}'.format(mf.mo_coeff[j,imo]))
f1.close()
mcc = ccsd.CCSD(mf)
mcc.conv_tol = 1e-12
ecc, t1, t2 = mcc.kernel()
eris = mcc.ao2mo()
# Next we calculate CCSD(T) which takes huge memory, fails most of the time.

e3ref = ccsd_t.kernel(mcc, eris, t1, t2)
l1, l2 = ccsd_t_lambda.kernel(mcc, eris, t1, t2)[1:]
print(ecc, e3ref)

eri_mo = ao2mo.kernel(mf._eri, mf.mo_coeff, compact=False)
nmo = mf.mo_coeff.shape[1]
eri_mo = eri_mo.reshape(nmo,nmo,nmo,nmo)
dm1 = make_rdm1(mcc, t1, t2, l1, l2, eris=eris)

with open('onepdm', 'w') as f1:
for i in range(orbshape[1]):
for j in range(orbshape[1]):
f1.write('\n{}'.format(dm1[i,j]))
f1.close
# Find the nautral orbitals (MO basis) and eigenvalues (i.e. occupation numbers).
e, c = scipy.linalg.eigh(dm1)
with open('occupnumb', 'w') as f1:
for i in range(orbshape[1]):
f1.write('\n{}'.format(e[i]))
f1.close
#
with open('natorbmo', 'w') as f1:
for i in range(orbshape[1]):
for j in range(orbshape[1]):
f1.write('{}\n'.format(c[j,i]))
f1.close

# Determine the CCSD 1electron energy terms based on density matrix
h1 = numpy.einsum('pi,pq,qj->ij', mf.mo_coeff.conj(), mf.get_hcore(), mf.mo_coeff)
E = numpy.einsum('pq,qp', h1, dm1)
print("E-CCSD(T) 1E part only", E)

# Get the nuclear repulsion energy and the number of electrons
enuc = gto.energy_nuc(mol)
print(" Nuclear Repulsion Energy is ", enuc)
ne = mole.tot_electrons(mol)
print(" Total number of electrons is ", ne)

dm2 = make_rdm2(mcc, t1, t2, l1, l2, eris=eris)

# format of the 2PDM

with open('twopdm', 'w') as f1:
for i in range(orbshape[1]):
for j in range(orbshape[1]):
for k in range(orbshape[1]):
for l in range(orbshape[1]):
f1.write('\n{}'.format(dm2[i,j,k,l]))
f1.close()
h1 = reduce(numpy.dot, (mf.mo_coeff.T, mf.get_hcore(), mf.mo_coeff))
e3 =(numpy.einsum('ij,ji->', h1, dm1)
+ numpy.einsum('ijkl,ijkl->', eri_mo, dm2)*.5 + mf.mol.energy_nuc())
print(e3ref, e3-(mf.e_tot+ecc))

• Please provide the PySCF input file and I'll try it in other software. Commented Aug 21, 2023 at 16:41
• Sure. Adding the details.
– Pro
Commented Aug 22, 2023 at 9:31
• Without "ccsd_opt_geom.xyz" we can't do anything. Also, is there a name for that basis set? I'm going to try to run this in MRCC. Commented Aug 22, 2023 at 14:12
• I am sorry, I have added the xyz-coordinate. Basis was a custom one, so I have replaced that with standard 3-21G basis set. Since this is supposed to be a generic problem, change of basis set should not be a problem. Please have a look at the modified input.
– Pro
Commented Aug 22, 2023 at 14:30
• Thanks for adding the geometry data and the basis set name, but now I see that you have an even number of electrons in the geometry followed by a charge of 1, so an odd number of electrons in total. I don't see the spin being specified anywhere though. Which spin multiplicity are you hoping to target? Also, the basis set does matter because your problem was about not having enough RAM to do the calculation. Commented Aug 22, 2023 at 14:31

## 1 Answer

With Psi4, I have been partially successful in obtaining the 1-PDM at the CCSD(t) level. Sharing the code,

import numpy as np
import psi4
import scipy

mol = psi4.geometry("""
0 1
O  0.025976  -0.000000  -0.067996
H  -0.028998  -0.000000  0.892837
H  0.966586  0.000000  -0.271648

symmetry c1
no_reorient
""")
psi4.set_options({'basis': '3-21G', 'freeze_core': 'false','print_mos': 'false','print':2, 'scf_type': 'Direct'})
_,cc_wfn = psi4.energy('ccsd(t)', return_wfn=True, molecule=mol) # CCSD(t) run
CCdm1=np.asarray(cc_wfn.Da())*2 # Obtain 1-particle density matrix as numpy array



For 2-PDM, I have obtained the T1 and T2 arrays but need to have l1 and l2 arrays from CCSD(T) lambda module. Which I can not do after a lot of research.

amps = cc_wfn.get_amplitudes()
TIjAb = amps['tIjAb'].to_array() # This looks like t2 of PySCF
TIA = amps['tIA'].to_array()     # This looks like t1 of PySCF
tau_IjAb = TIjAb + np.einsum("ia,jb -> ijab", TIA, TIA)