I’m unclear on how the CEPA(0) implementation considers singles. Is it analogous to Szabo’s MBPT equations:
-B^\dagger D B = E_{corr}
B = \braket{\Psi_0|H|\Psi_{ij}^{ab}}
D = \braket{\Psi_{kl}^{cd}|H|\Psi_{ij}^{ab}}
but with D including single/single and single/double terms and with B including single terms? The Wennmohs & Neese paper mentions an ambiguity in the treatment of singles, but I believe that ambiguity should go away if all the shifts are 0 anyway, right?
The implementation of CEPA(0) with singles in the FNOCC module should be equivalent to linearized CCSD, which I think is what you’re describing with zero shift for the singles. A long time ago I verified that implementation against the CEPA code in Molpro for CEPA-0,1, and 3, so I believe the singles are handled in a consistent way with that code.
Does that mean that LCCD should be equivalent to using the CEPA_NO_SINGLES flag? (I ask because I’m getting different numbers for these, and I’m not sure if it’s an input issue on my end, a misunderstanding of the theory, or a bug with psi4.) Computed energies:
LCCD:
-460.2612756198385
LCCSD:
-460.26183389044934
CEPA(0)D:
-460.26183389045053
CEPA(0)SD:
-460.26183389045025
Code input:
import numpy as np
molecule = {
0 1
H 0 0 0
Cl 0 0 1
}
psi4.core.be_quiet()
set basis cc-pvdz
set reference rhf
From psi4/psi4/driver/procrouting/proc.py it looks like CEPA_NO_SINGLES is set to false for the method ‘cepa(0)’. I am not sure why - perhaps for consistency with the DCFT and DFOCC modules, which I think also can run cepa0. Maybe @loriab knows why the procedure is defined that way and if it is possible to override that procedure from within the input file?
Ok, I think all is working as expected. Consistent with the table, cepa(0) includes singles and is equivalent to lccsd. lccd is equivalent to cepa(0) without singles. By declaring method ocepa(0), you’re entering into a contract that overrides keywords, so no-singles doesn’t work (@hrgrimsl, you were setting as a kwarg, not a keyword anyways).
So the below is what I get, and I think it’s sane. Yes, use lccd.
molecule {
Ne
}
e = energy('lccsd/cc-pvdz')
print(e)
e = energy('lccd/cc-pvdz')
print(e)
e = energy('cepa(0)/cc-pvdz')
print(e)
This is fine, and I understand what is going on now ,but I had originally intended for CEPA_NO_SINGLES to control whether or not singles were used in any CEPA-n generated through FNOCC. it seems weird that that option will work for CEPA-1/3 but not for CEPA-0.