Fa_subset("MO") Incorrect

[JonathonMisiewicz]

Run:

molecule h2o {
0 1
o
h 1 0.958
h 1 0.958 2 104.4776 
}

wfn = energy('scf/cc-pvdz', return_wfn=True)[1]
print(wfn.Fa_subset("MO").nph)

The result MO basis Fock matrix is not diagonal, as expected.

The correct way to transform the SO basis Fock matrix into the MO basis matrix is contraction against the C coefficient matrices. This is not what the Fa_subset function does. It instead goes through a series of matrix mulitplications involving the overlap matrix. This is the correct transformation rule for the density matrix. Not-so-coincidentally, Da_subset goes through exactly the same function, where this behavior is correct. If you try to use the same logic for both helpers, you're going to get one of them wrong, guaranteed.

While the mechanism of the bug is clear to me, I'm not sure what the best fix should be. Perhaps we should intercept "MO" as a special case in Fa_subset so it never sees the matrix_subset_helper? Under the hood, the code would return Fa_->transform(Ca_).

[susilehtola]

He in STO-3G-decon is maybe an even simpler example.

[CDSherrill]

Interesting. I have examined these functions in the past, but it's been a while. Can you remind me what is the original basis, AO or SO?

It's not super obvious to me why behavior for D would be correct while behavior for F would be wrong, if they both start in the same original basis.

[JonathonMisiewicz]

The original basis is SO.

The behavior for D and F are different because "AO basis quantities" mean different things for those two cases. For the F matrix, "AO basis" means that you have a formula in terms of one-electron functions, and you just throw those in. To change basis, you just linearity. F_pq = C_mu,p C_nu,q F_mu,nu

For D, "AO basis" means "if you contract this against an AO basis integral, you'll get the same result as if you had contracted an MO basis quantity against MO basis integral." In that case, D_mu,nu F_mu,nu = D_p,q F_p,q = D_p,q C_mu,p C_nu,q F_mu,nu, so we end with D_mu,nu = C_mu,p C_nu,q D_p,q.

Note that it differs whether you need to contract your C matrices against the AO or the MO basis quantity.

[gpwood]

Hi All,

I see that this issue is still open. I was wondering if there is a current work around? I tried some simple experiments to see if it was still producing the incorrect output:
> psi4 --version
1.7
input file contents:

mol = psi4.geometry("""
H   0.000   0.000   0.000
H   0.000   0.000   0.740
Symmetry c1
""")
psi4.set_options({'basis': 'sto-3g'})
scf_e, wfn = psi4.energy('SCF', return_wfn=True)
print("MO coefficients")
print(wfn.Ca().to_array())
print("Fock in the AO basis")
print(wfn.Fa_subset("AO").to_array())
print("Fock in the MO basis")
print(wfn.Fa_subset("MO").to_array())

this gives the following output:

MO coefficients 
[[ 0.54884228  1.21245192]
 [ 0.54884228 -1.21245192]]
Fock in the AO basis
[[-0.36607883 -0.59428702]
 [-0.59428702 -0.36607883]]
Fock in the MO basis
[[-1.59408547e+00  8.04278118e-16]
 [ 7.39285035e-16  7.76197397e-02]]

and in the output file

    Orbital Energies [Eh]
    ---------------------
    Doubly Occupied:                                                      
       1A     -0.578578  
    Virtual:                                                              
       2A      0.670950 
@DF-RHF Final Energy:    -1.11678331788308

the orbital energies in the output file are correct, which should be the diagonal elements of the Fock matrix in the MO basis, but these are different.

I check the AO and MO coefficients produced by two other QM programs (pyscf is below) to see where the inconsistency might be, they both produced the "correct" Fock matrix (i.e. the diagnoal elements in the MO basis are the orbital energies) with the following MO coefficients:
pyscf input

import numpy
from pyscf import gto, scf, lo

geometry = '''
    H   0.000   0.000   0.000
    H   0.000   0.000   0.740
 ''' 

mol = gto.M(atom=geometry,
           basis='STO-3G') 
mf = scf.RHF(mol)
mf.kernel()
print(mf.kernel())
Fao = mf.get_fock()
print("MO coefficients")
print(mf.mo_coeff)
print("Fock in AO basis")
print(Fao)
Fmo = mf.mo_coeff.T @ Fao @ mf.mo_coeff
print("Fock in MO basis")
print(Fmo)

converged SCF energy = -1.11675930739643
-1.1167593073964255
MO coefficients
[[ 0.54884228 -1.21245192]
 [ 0.54884228  1.21245192]]
Fock in AO basis
[[-0.36602603 -0.59429997]
 [-0.59429997 -0.36602603]]
Fock in MO basis
[[-5.78553860e-01 -2.43968808e-16]
 [-2.12179326e-16  6.71143492e-01]]

[JonathonMisiewicz]

How soon do you need this? This bug is an easy fix, and I can include it in the Psi4 1.8 release coming out in a couple weeks. If you need this sooner, the "workaround" is simple.

[gpwood]

@JonathonMisiewicz that's great news! I'm in no rush, so I can wait. The purpose of my calculations is to look at Fock, Coulomb and Exchange matrices in different orbital representations (i.e. from some localisation procedure). For your tests have you looked at the J and K matrices as well as orbital rotations? Thanks, Geoff.

[JonathonMisiewicz]

Your questions may make sense in the context of your research, but they don't make sense in the context of the actual Psi4 code.

• As far as I'm aware, Psi has no built-in function to return J or K in an MO basis, only in the SO basis.
• Psi stores the Fock matrix in the SO basis and transforms it to an MO basis. The problem right now is that it's using the wrong transformation rule. If we can get right in the canonical basis, then it'll be right in any orthogonalized basis.

[gpwood]

Completely understand. Just wanted to make sure that if you make a call for K, J, H(core) and density matrices it wouldn't have the same transformation issues as for the Fock. Obviously as long as you can get the rotation matrix/MO coefficients I can make the transformations on J and K myself. Thanks.

[JonathonMisiewicz]

Oh, no, those would be just fine. The issue is that Psi isn't using the formula mf.mo_coeff.T @ Fao @ mf.mo_coeff to transform Fock, but a more complicated formula which describes how the density matrix transforms. If you can get the SO Fock matrix and do the transformation yourself, it'd work.