KE from CCSD density matrix not matching value in output file

Hello all,
I’m running a RHF and CCSD calculation for H2O and trying to reproduce the 1-electron level energies (i.e. everything except the 2-electron repulsion) using the corresponding 1-electron density matrices and the K, V matrices like this.

    grad, ccwfn = psi4.gradient('ccsd', return_wfn=True, ref_wfn=wfn)
    mints = psi4.core.MintsHelper(ccwfn.basisset())
    Da_cc = ccwfn.Da().to_array(copy=False, dense=True)
    K = mints.ao_kinetic().to_array(copy=False, dense=True)
    V = mints.ao_potential().to_array(copy=False, dense=True)
    E_ke_cc = 2 * np.trace(np.matmul(K, Da_cc))
    E_v_cc = 2 * np.trace(np.matmul(V, Da_cc))

The HF energies match the output file numbers. However the KE from the CCSD 1-RDM differs from the one in the output file by 1.69mH.

  1. Why aren’t they matching exactly (till floating point precision)?
  2. Also where can I find the CCSD potential (nuclear-electron) energy in the output file?

Here’s the full script which will print an output like the following.

Any help would be greatly appreciated. Thank you!

Energies from HF 1-RDM
Nuclear: 9.168193296424349
Kinetic: 74.58220348386018
Potential: -196.9063503606068
Coulomb: 47.29435182717615
Exchange: -9.101796533055854
Total : -74.96339828620196

Energies from CCSD 1-RDM
Nuclear: 9.168193296424349
Kinetic: 74.62366469678517
Potential: -196.82537980565303
Coulomb: 47.24918874993827
Exchange: -9.082702767239162
Total : -74.8670358297444

The issue appears to have something to do with whether the OPDM is relaxed or not. I think that the CCSD kinetic energy printed in the output is calculated using the unrelaxed OPDM, but the one you extract in your script is relaxed. I modified your script to ask for the unrelaxed OPDM, and then these values agree. Try

psi4.set_options({'opdm_relax': False})

Awesome, thank you so much!

Can you/someone please help with the second question as well? Where can I find the CCSD potential (nuclear-electron) energy in the output file or as a variable? I tried the following but none except the last one work.

    print(psi4.variable("CC CORRELATION KINETIC ENERGY"))
    print(psi4.variable("CC CORRELATION POTENTIAL ENERGY"))

    print(psi4.variable("CCSD CORRELATION KINETIC ENERGY"))
    print(psi4.variable("CCSD CORRELATION POTENTIAL ENERGY"))

   # Only the following works
    print(psi4.variable("CCSD CORRELATION ENERGY"))

Thanks again!

Whoops, I didn’t notice the second question. From here:

it looks like those variables should be set, but the potential energy will the total potential correlation energy, not just the one-electron part.

If you want the one-electron part, you can modify that function. You can see where the kinetic energy integrals are grabbed and used here:

You would just need to grab the potential energy integrals instead [wfn->mintshelper()->so_potential()] and carry out the same operations.

1 Like