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

Question

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

Answer 1

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)