Overwriting Fock Matrix In Wavefunction (FCI)

Hi all,

Is there a way to do an FCI calculation on a wavefunction with a customized Fock matrix? I’m running several FCI calculations using a low-level reference wavefunction (calculated using HF), and I need to apply a potential to specific orbitals within the active space. Overwriting the Fock matrix directly would be ideal. I’ve tried setting Fa and Fb using

low_level_wfn.Fa().copy(new_fock_matrix)
low_level_wfn.Fb().copy(new_fock_matrix)

and then passing the updated wavefunction to FCI with

energy(‘FCI’, ref_wfn=low_level_wfn, molecule=x, return_wfn=True)

However, this doesn’t seem to be working; regardless of the values in new_fock_matrix, the calculation always gives the same result. Is there some other way of overwriting it, either Python-side or C+±side?

Thank you!

Just to be sure, have you tried the clone() function. Something like:

low_level_wfn.Fa().copy(new_fock_matrix.clone())
low_level_wfn.Fb().copy(new_fock_matrix.clone())

I’ve just tried using clone(), and unfortunately, there is still no change in the result. (Thank you for the suggestion!)

Ok. When I use this simple example, I see that I could overwrite the wfn’s fock matrix. Here is my input.

molecule h2o{
O
H 1 R
H 1 R 2 A

R = 0.957
A = 104.5
}

set basis sto-3g
Escf, scf_wfn = energy(‘SCF’, molecule=h2o, return_wfn=True)

print(scf_wfn.Fa().get(0,0,0))
a = scf_wfn.Fa()
a.set(0,0,0, 500)
print(a.get(0,0,0))
print(scf_wfn.Fa().get(0,0,0))

Efci, fci_wfn = energy(‘FCI’, molecule=h2o, ref_wfn=scf_wfn, return_wfn=True)
print(fci_wfn.Fa().get(0,0,0))

Yes, this seems to return the correct values when looking at the Fock matrix… However, changing an element of the Fock matrix ought to affect the calculated energy and orbital occupations, right? If I change the a.set(0,0,0, 500) to a.set(0,0,0, 0), Efci doesn’t change from the original case at all, so it seems like it’s still using the same Fock matrix within the FCI energy calculation. The output files for the two cases are completely identical besides the lines for input file, process ID, and timestamps.

Also, I noticed this line in the output file:

Transforming the one-electron integrals and constructing Fock matrices

Could it be reconstructing the Fock matrices, and if so, is there a way to prevent this?

Yes, you are totally right about this. You need to make changes in the C coefficients as the fock matrix is reconstructed from the SCF’s C coefficients in the detci module, responsible for doing FCI calculations. So,
In the above example, just replace Fa() with Ca() and you could see it does overwrite the Fock matrix.

Perfect! I changed Ca(), and the energy and Fock matrix changed accordingly. Thank you!

Just to chime in here real quick, post-SCF methods do not use the Fock matrix as far as I know. The only two things that they pull from is the one-electron hamiltonian and the orbital coefficient matrix. There are some special cases, but these are the only two to depend without intimate knowledge of the code.

Thank you! Is there a way to directly overwrite the one-particle Hamiltonian or apply an arbitrary potential to it, then? I tried overwriting the Hamiltonian instead of the Fock matrix based on Ashutosh’s input above,

molecule h2o{
O
H 1 R
H 1 R 2 A

R = 0.957
A = 104.5
}

set basis sto-3g
Escf, scf_wfn = energy(‘SCF’, molecule=h2o, return_wfn=True)

print(scf_wfn.H().get(0,0,0))
a = scf_wfn.H()
a.set(0,0,0, 5000)
print(a.get(0,0,0))
print(scf_wfn.H().get(0,0,0))

Efci, fci_wfn = energy(‘FCI’, molecule=h2o, ref_wfn=scf_wfn, return_wfn=True)
print(fci_wfn.H().get(0,0,0))
print(Efci)

Printing the Hamiltonian shows that index (0,0,0) has been changed to the appropriate value, but the FCI calculation always ends up yielding the same energy, regardless of the value I use in the set() function.

I believe and @dgasmith can correct me if I am wrong here that the one particle hamiltonian in post-scf methods are read from binary PSI4 files instead of the reference wavefunction. I don’t think there is a python interface available for writing to the PSIF_OEI file where the one particle hamiltonian matrix is stored. For now, the best option might be to create a psi4 plugin and do these things on the C++ side.

It’s more inconsistent than I thought. We should probably modify libtrans to use the incoming wavefunctions Hamiltonian. This should simplify quite a few operations.

Strangely enough the above will work for CASCF, but not FCI. A quirk of some underlying coding forms.

Should I go ahead and make a small pull request to fix this inconsistency?

@ashutosh Go for it, it may be slightly involved as you probably need at add the one-electron hamiltonian to the libtrans signature and then change the libtrans one-electron transformation to use this passed in H. This should require touching every project that uses libtrans.

Thanks to both of you for the help! I’ll keep an eye out for the pull request, if you end up making one, but in the meantime, I’ll go ahead and make a C++ plugin for my implementation.

Just wanted to let you know that now you can pass your one electron hamiltonian to the wavefunction successfully from the python interface. The pull request just got merged.

1 Like

Excellent. Thank you so much!