Saving CCSD density for read-in

#1

I want to save the CCSD density so that I can read it into a DFT/HF calculation later. It seems (unsurprisingly) that with the input I have below I’m saving the SCF density (at least when I read it in in a new calculation):

molecule {
0 1
H  1.0327262753 0.6722433198 0.0000000000
H  1.7327262753 0.6722433198 0.0000000000
symmetry c1
}

set {
basis                   def2-qzvppd
e_convergence           1e-12
freeze_core true
cc_type df
}

en, wfn = energy('ccsd', return_wfn=True)

set {
basis                   def2-qzvppd
e_convergence           1e-12
guess read
maxiter 0
fail_on_maxiter false
}

energy('scf')

My question is: is there a way to write the CCSD density so that I can read it in an SCF guess as done above?

#2

Do you want to calculate the HF/DFT energy using a CCSD density?
Restarting a HF calculation is technically only possible using orbital coefficients, not a density.

The CCSD density is not available for energy calculations as its not needed and expensive to calculate.
Here you would typically request a gradient computation.

#3

Yes I want to calculate the HF/DFT energy using a CCSD density. Is there a way to do this? Changing “energy” to “gradient” doesn’t seem to have done anything different.

#4

This is not so easily done because our JK object cannot take a density as input only orbitals.
So then you need to compute the integrals and contract them with the density yourself.

This is not overly difficult and there are plenty examples on GitHub - psi4/psi4numpy: Combining Psi4 and Numpy for education and development.
Depending on the system size you could use DF integrals: psi4numpy/density-fitting.ipynb at master · psi4/psi4numpy · GitHub

e, wfn = psi4.gradient('ccsd',return_wfn=True)
D = wfn.Da() # CCSD density

The so obtained density could be plugged into the HF equation.

For DFT the additional complication is to build the potential. Easiest is to do a DFT calculation first and afterwards exchange the density and recalculate the potential and energy. This is described here: Obtain the energy from DFT at arbitrary density matrix - #3 by dgasmith

There might be more elegant solutions that I haven’t though of. Possibly one can reuse the integrals written to disk during the CCSD calculation but for a prototype code the psi4numpy approach will get you far already

#5

This is not the same thing as

and I need to be absolutely clear about which you’re asking about, because they’re both unusual things to ask.

You could be asking for computing the energy delivered by a functional for the CCSD density or using the CCSD density as an initial guess for the iterative solution of the KS-DFT equations. If you want to use CCSD as an initial guess, I have to ask if this is absolutely essential, or you just want a robust “trick” to converge KS-DFT for difficult systems.

There are two annoying technicalities to worry about: first, the “CCSD density” is ambiguous. Some people define the CCSD density excluding the so-called orbital response terms, and some people define the CCSD density including the so-called orbital response terms. If you just want an initial guess, you should be fine without orbital response terms. Otherwise, it is essential that you know the difference and specify which you want.

Second, some functionals may be ill-behaved for non-idempotent densities or even undefined without orbitals to compute the non-density quantities used by modern DFT (example: kinetic energy density).