CCSD with Arbitrary MO Coefficients

[hrgrimsl]

Hi all,

I’m trying to use custom MO coefficients within a CCSD calculation, but am obtaining results which are both wrong and inconsistent between different QC modules for a trivial example. This may be related to Post-SCF with custom matrices, but I don’t need it to work for open-shell references. Here is a minimal working example, using pseudocanonicalized MP2-natural orbitals. (I.e. the Fock matrix is only diagonal within the occupied-occupied and virtual-virtual subspaces.)

min_example.py (2.7 KB)

It’s a 2-electron problem, so CCSD should be exact, regardless of orbital choice, but that’s not what I’m seeing:

Exact energy: -1.1400734808025696
Old CCSD energy: -1.1400734808003057
FNOCC CCSD energy with new orbitals: -1.1407206561955716
CCENERGY CCSD energy with new orbitals: -1.1400706133026943

Do I need to be doing something other than wfn.Ca().copy(…) to make CCSD care about the off-diagonal part of F? Is one of the qc_modules preferred for this?

[deprince]

I can’t speak for what CCENERGY does.

With conventional ERIs (without DF/CD), I don’t think FNOCC will respect the fact that you have F_ij != eps_i delta_ij. If you use the DF/CD form of FNOCC, then I think these terms should be included.

If you use CD with a super tight CD threshold (1e-12?), then you should recover the conventional ERI result. Can you try that and let me know the result?

[crawdad]

The ccenergy code does not assume canonical, Hartree-Fock orbitals, even in the spin-restricted variant.

[hrgrimsl]

@deprince CD causes a segmentation fault as soon as DF-CCSD starts:

min_example.py (2.9 KB)

I tried playing with CACHELEVEL, the memory allowance, etc.- am I missing something?

[deprince]

I can’t reproduce the segfault with min_example.py, and I’m getting a different wrong energy than you. I think I need to clone a fresh psi4 and then I’ll revisit this.

[hrgrimsl]

I was using the last 1.3.2 release via a conda install- updating to the dev head removed the seg fault, but I also just got a different wrong energy with a tight cholesky_tolerance:

min_example.py (2.9 KB)

Exact energy: -1.1400734808025708
Old CCSD energy: -1.1400734808029145
FNOCC CCSD energy with new orbitals: -1.1404058717422065

(The energies obtained with PK scf/conventional cc didn’t change.)

[jmisiewicz]

The fnocc and cc problems are probably quite different, so let’s go one at a time.

My hunch for the problem in the cc module is that your Fock matrix is incorrect. The cc module obtains the Fock matrix when executing the relevant function here. Whether RHF or UHF, the Fock matrix (AO basis) is read from the reference wavefunction and transformed into the MO basis.

There’s no mechanism in your script to change the Fock matrix of the wavefunction you send in compared to the Fock matrix of your original Roothaan-Hartree-Fock solution, so the energy should only be correct if the AO basis Fock matrix is invariant to your change of orbitals. It isn’t.

The expression for the Fock matrix for not-necessarily-canonical orbitals is f^p_q = h^p_q + gbar^pi_qi, where the double i denotes summation over occupied orbitals. Different choices of occupied orbitals mean different two-electron integral terms, whatever your basis.

So, is your choice of occupied orbitals different from the SCF choice? Certainly. If your new orbitals preserved the occupied space, the Fock matrix should be block-diagonal in the occupied and virtual spaces, but I gather that it isn’t, if this is “pseudocanonicalized”. Therefore, your Fock matrix is incorrect in the AO basis.

As a demonstration, try the following

tricky_wfn = psi4.energy("hf/cc-pvdz", return_wfn=True)[1]
scf_wfn = psi4.energy('b3lyp/cc-pvdz', return_wfn=True)[1]

scf_wfn.Ca().copy(tricky_wfn.Ca())
scf_wfn.Cb().copy(tricky_wfn.Cb())
trick_ccsd_E = psi4.energy("ccsd/cc-pvdz", ref_wfn = scf_wfn)

Even though the orbitals you send in are perfect, the Fock matrix coming out of DFT isn’t the Fock matrix as defined by CC, so you don’t even get zero contribution to MP2 singles.

I leave it to @crawdad to suggest a fix. I know that the cc mega-module has some orbital rotation machinery, so it may be possible to reconstruct the Fock matrix rather than trusting the wavefunction to be consistent.

[crawdad]

It may be a question of exactly when in the workflow you’re inserting the orbitals. The integral transformation code called within cctransort should be building the necessary Fock operators using the orbitals provided in the ref, should it not? Tagging @andysim who may know more, since he wrote the integral transformation code.

[hrgrimsl]

Manually correcting the AO basis Fock matrix to account for the new orbitals per Jonathon’s comments makes the CCENERGY module give the correct energy. Thank you all for the help.

min_example.py (3.5 KB)

[jmisiewicz]

@crawdad I haven’t seen any part of cc or cctransort that builds the Fock matrix apart from the orbital rotation code. As far as I can tell, they are taken from the wavefunction and assumed to be valid.

@hrgrimsl Are there any other questions, at this point? I think that fnocc assumes canonical orbitals, and there’s nothing you can do about that, short of changing the code.

[hrgrimsl]

No, I think that gives me the tools to do what I’m trying to do- thanks again.

[jmisiewicz]

Marking this as solved, then.