(T) in frozen-core BCCD(T) calculations

I think there is a problem in the (T) component of frozen-core BCCD(T) calculations, probably both RHF and ROHF. Upon comparison to Molpro or CFour, everything agrees as along as one correlates all the electrons (which is the only option for CFour). However, in the frozen core case, the BCCD energies of Molpro and PSI4 agree very well, but the (T) contributions differ. I believe this is due to the semi-canonicalization done for the (T) - probably PSI4 uses all the orbitals instead of just the correlated ones. I don’t think this is the desired effect since then one cannot unambiguously compare CCSD(T) with BCCD(T). Here are results for the CO molecule (cc-pVQZ basis set with r=1.6 Å).

         HF          CCSD        (T)          BCCD        (T)

MOLPRO -112.600216 -113.034548 -0.041877 -113.026871 -0.046423
PSI4 -112.600216 -113.034548 -0.041878 -113.026868 -0.050950

I forgot to note that these energies correspond to the X2C Hamiltonian with cc-pVQZ-DK basis sets. The analogous situation occurs for the non-relativistic case as well however.

Kirk,

Is this happening for both RHF and ROHF? I know you said “probably”, but I’d like to confirm, because our semi-canonicalization procedure should certainly not include frozen orbitals. On the other hand, ROHF orbital definitions vary from code to code, which could yield different frozen orbitals.

-Daniel

I can only test the frozen-core BCCD(T) against Molpro since CFour requires all-electrons to be correlated. Molpro only has closed shell BCCD(T). That’s why the example I posted was just for the closed-shell CO molecule which already exhibits the problem. I’m just assuming the same thing will happen in the ROHF case.

-Kirk

Could you either send me an input file or try setting “RUN_CCTRANSORT = false” (or both)?

-TDC

Simple input below for stretched CO. The corresponding molpro energies are:

RHF: -112.5333186
CCSD: -112.9674459
CCSD(T): -113.0092946
BCCD: -112.9597779
BCCD(T): -113.0061622

#! Sample

memory 450 mb

molecule CO {
0 1
C
O 1 r
r=1.6
}

set {
reference rhf
basis cc-pvqz
maxiter 150
docc [5,0,1,1]
frozen_docc [2,0,0,0]
units angstrom
MOM_START 10
SCF_TYPE PK
}

energy(‘scf’)
energy(‘ccsd(t)’)
energy(‘bccd(t)’)

As you probably know by now, setting RUN_CCTRANSORT = false does not change the situation.
-Kirk

Just wondering if this issue has gotten any escalation?

@Kirk Could you send me the full output files from molpro for the streched CO test case? I am trying to track down the source of differences you are seeing, but since I don’t have access to molpro I am “flying blind”. This is being looked into though and thank you for bringing it to our attention.

Attached are both Molpro and Psi4 outputs for a frozen-core cc-pVDZ calculation on closed-shell CO. Note that the CCSD(T) and BCCD results match very closely. It’s only the (T) contribution to BCCD(T) that differs.
co_mol.txt (23.7 KB)

co.txt (172.6 KB)

Oh, and thanks very much for looking into this!

@Kirk Could I also see the output file for the same test case in molpro with all electrons correlated?

We have done some tinkering in psi4 and found that if we explicitly ignore contributions from the off-diagonal elements of fock matrix to the (T) correction our frozen core BCCD(T) energy matches Molpro dead on. These terms are zero for canonical Hartree Fock orbitals but they should be nonzero for semi-canonical orbitals used in BCCD(T).

If Molpro was ignoring these terms for BCCD(T) that would lead to differences between psi4 and Molpro when all electrons are correlated as well, and not just for frozen core calculations.

Strange, since Molpro explicitly states in the output that it is block diagonalizing the Fock operator before it starts the triples calculation. I just ran the all-electron calculations and the difference between the two codes persists. I could have swore I tested this before but I guess I only looked at psi4 compared to cfour and g09. I just ran cfour for this same job (all-electrons correlated) and it agrees with psi4.

After all of this, I’m afraid perhaps it is Molpro that has a problem and not psi4. My apologies for time wasted! If I find out something more, I’ll let you know.

To clarify AJ’s post a bit, block diagonalization (semicanonicalization) of the Fock operator is not the source of the discrepancy; indeed, it’s necessary for the non-iterative (T) correction to be valid and both codes apparently do this correctly. However, PSI4 also includes contributions to the (T) correction involving the occupied-virtual block of the Fock operator, which are of the same perturbational order as the disconnected singles terms. If we leave these out, we match the Molpro B-CCD(T) energy for the CO test case.

Aha, thanks Daniel. That makes perfect sense. Are these additional contributions well known and published somewhere? I’d like to educate myself a bit more on this.

Yes, these terms are well known. A rather complete discussion can be found in Watts et al., J. Chem. Phys. 98, 8178-8733 (1993). Note, in particular, the E[4]_DT term in Eq. (24).

Thanks, I’ve been tracking them down. Looks like you also discuss this is a bit in your chapter with Fritz in Rev. Comp. Chem.

Molpro does this correctly for the open-shell ROHF case, so perhaps the fix will not be so difficult.