Saving CCSD density for read-in

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?

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.

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.

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

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

Sorry for the very-long wait and subsequent revive of this topic…

I was interested in this because of a recent Burke paper “Density-Corrected DFT Explained: Questions and Answers” where they show that using the Hartree-Fock density instead of the self-consistent KS-DFT density provides significantly better results in many cases. I wanted to play around with using running DFT on highly accurate electron densities, which I assume CCSD could provide.

I’ll look into the method proposed by Hokru from dgasmith.

Because I don’t know your background, I do need to warn you about something you may already know:

DFT relies on a density, which is a function of coordinates corresponding to a single particle. To enforce that you only optimize to densities that can be produced by some physical wavefunction, Kohn-Sham DFT represents a density by a set of orbitals.

Now, methods like SCF and CCSD provide a one-electron reduced density matrix (1RDM), which we also call a density, but this is not the same object as in the DFT case. You can freely interchange between an SCF 1RDM and a space of occupied orbitals, and you can convert from either of them to a real-space electron density, but you cannot in general go backwards. The technical speak for a 1RDM convertible to a set of orbitals is idempotent; I’m taking Hermiticity for granted.

Although it’s theoretically possible to apply a DFT functional to the electron density delivered by a 1RDM that isn’t idempotent, I have no idea how approximate DFT functionals will react (you don’t automatically have orbitals) and no idea if that would break Psi’s DFT infrastructure in important ways.