One- and two-body integrals from Hartree-Fock

[mtfishman]

Hi @gustavojra,

Thanks for the nice package!

Is there a simple way to obtain the coefficients of the second quantized Hamiltonian (the one- and two-body integrals) from a Hartree-Fock solution?

Thanks,
Matt

[gustavojra]

Hi @mtfishman,

The computation returns the orbital coefficients C. There should be a simple formula relating these coefficients to the elements of this matrix. I just don't remember off the top of my head. My guess is that you are simply looking for the density matrix:

Dμν = Cμi Ciν

If your second quantized Hamiltonian looks like

H = hμν + gμνρσ

then contracting D in there should return the HF energy

Ehf = Dμν*hμν + Dμν*Dρσ*gμνρσ

Just watch out for changes due to spin integration.

If that's so, all you need is to slice the occupied orbitals (Co) and multiply by the transpose of itself (Co')

wfn = @energy rhf
C = wfn.orbitals.C
ndocc = wfn.molecule.Nα
Co = C[:, 1:ndocc]
D = Co * Co'

[mtfishman]

Thanks for the quick response. Apologies if I'm missing something obvious (I'm not so familiar with quantum chemistry notation and conventions), but can I get the second quantized Hamiltonian itself (hμν and gμνρσ) from Fermi.jl (either the original one, or the one in the Hartree-Fock basis)? Just to be clear about the notation, I assume you mean the second quantized Hamiltonian is:

H = ∑_μν hμν a†_μ a_ν + ∑_μνρσ gμνρσ a†_μ a†_ν a_ρ a_σ

Also, what would the convention for the orbital that get summed over: both the spin up and spin down orbitals, or just the spin up orbitals, with the assumption that the spin down orbitals are the same, etc.?

[gustavojra]

1. First about the Hamiltonian - and quantum chemistry lingo

Yes, that is the second quantized Hamiltonian (with the creation and annihilation operators). But to be honest, this formalist is not normally used for Hartree-Fock. That is because the goal of the Hartree Fock method is to obtain a set of independent particle model (IPM) functions (i.e. orbitals) on which second quantized models can be built upon taking the Hartree-Fock determinant as the Fermi vacuum.

Orbitals are linear combinations of the atomic functions in the basis set χ = {μ, ν, ρ, σ, ...}:

|i> =  ∑ cμi χμ

C is the matrix I mentioned above. It can be obtained from Fermi.

The second quantized Hamiltonian you wrote is a little inconvenient because the AO basis χ = {μ, ν, ρ, σ, ...} do not form orthonormal states. If we work with the orbitals from the Hartree-Fock instead, this becomes

H = ∑_pq hpq a†_p a_q + ∑_pqrs gpqrs a†_p a†_q a_r a_s

which is the second quantized Hamiltonian used in correlated methods such as Coupled Cluster.

There is no a† or a objects in Fermi. While these objects are used in deriving working equations they obvious must vanish at the end. Fermi has chewed up equations already, no derivation or fancy Wick's stuff happens here. However, it does return wave function parameters that you can use to build your Hamiltonian.

• hμν and gμνρσ are simply AO integrals, no HF is needed to get it.

hpq and gpqrs depend on RHF orbitals coefficients, they are the AO integrals transformed to MO.

Normally, hpq will be the Fock matrix (effective one-electron operator). In canonical HF, this matrix is diagonal and its values are taken as orbital energies. From the wfn.orbitals object in Fermi, you can get these numbers from the eps field.

wfn = @energy rhf
wfn.orbitals.eps |> Diagonal
7×7 Diagonal{Float64, FermiMDArray{Float64, 1}}:
 -20.2459    ⋅         ⋅          ⋅          ⋅         ⋅         ⋅ 
    ⋅      -1.25223    ⋅          ⋅          ⋅         ⋅         ⋅ 
    ⋅        ⋅       -0.600785    ⋅          ⋅         ⋅         ⋅ 
    ⋅        ⋅         ⋅        -0.449034    ⋅         ⋅         ⋅ 
    ⋅        ⋅         ⋅          ⋅        -0.389376   ⋅         ⋅ 
    ⋅        ⋅         ⋅          ⋅          ⋅        0.573384   ⋅ 
    ⋅        ⋅         ⋅          ⋅          ⋅         ⋅        0.703995

Note that this is water sto-3g. There should be 7 alpha + 7 beta orbitals. The one-electron matrix is always the same for both cases, thus we only really solve for the alpha case. That is, we reduce the 14x14 problem into a 7x7 problem.

gpqrs will be the repulstion integral (or a combination of permutations of itself). It can be obtained with the IntegralHelper

wfn = @energy rhf
moints = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals)
points["ERI"] 
7×7×7×7 FermiMDArray{Float64, 4}:
[:, :, 1, 1] =
...

Summing up
While we don't have annihilation or creation formalism built into Fermi, we provide wave function parameters that allows you to reconstruct the Hamiltonian you want using second quantization (though you may need to work out how it looks like for you, since there are many equivalent ways to do that)

Sorry for the verbose, I am trying to make sure we are on the same page

2. On spin

As a mentioned for the Fock matrix, we normally do not solve a giant problem for all possible spin combinations. We separate the problem into alpha and beta. For closed shell systems, only one of those needs to be solved for. The density matrix that you get from RHF is defined for station orbitals only, that is, spin has been integrated out.

It is possible to reconstruct thing in a spin-orbital format. It is not computationally advantageous though. But let me show you an example

@molecule {
       H  0.0 0.0 0.0
       H  0.76 0.0 0.0
       }
@set {
basis sto-3g
diis false
}
wfn = @energy rhf
Cspatial = wfn.orbitals.C
2×2 FermiMDArray{Float64, 2}:
 -0.551016  -1.18981
 -0.551016   1.18981

This coefficient matrix is only 2x2, if we were working with spin-orbitals it should be 4x4. We can reconstruct the Cspinorbital as

z = zeros(size(Cspatial))
Cspinorbital = [Cspatial z ; z Cspatial]
4×4 Matrix{Float64}:
 -0.551016  -1.18981   0.0             0.0
 -0.551016   1.18981   0.0             0.0
  0.0             0.0          -0.551016  -1.18981
  0.0             0.0          -0.551016   1.18981

To construct the density matrix you sum over occupied orbitals only.

When you run a RHF calculations you get N orbitals back, N is the size of your basis set. Thus, you always get more orbitals than electrons. The extra orbitals are called virtual and do not contribute to the energy, therefore are sliced out of the density matrix

Co = Cspatial[:, 1:1]
D = Co * Co'
2×2 FermiMDArray{Float64, 2}:
 0.303618  0.303618
 0.303618  0.303618

In spin-orbital format

D = 
4×4 Matrix{Float64}:
 0.303619  0.303619  0.0       0.0
 0.303619  0.303619  0.0       0.0
 0.0       0.0       0.303619  0.303619
 0.0       0.0       0.303619  0.303619

Again sorry about the lengthy response.

I think ultimately what you want is the integrals (AO or MO ones) and the coefficient matrix C. However, you need to figure out how to piece those things together. I assure you that those arrays will contain all the information one could need from Hartree-Fock.

[mtfishman]

Thanks a lot for taking the time to give a thorough response, that is very helpful and I think answers all of my questions.

Perhaps I should have given more context from the start, but my purpose is to obtain the coefficients hpq and gpqrs for the full many-body quantum chemistry Hamiltonian:

H = ∑_pq hpq a†_p a_q + ∑_pqrs gpqrs a†_p a†_q a_r a_s

since I want to form that many-body Hamiltonian as an MPO in ITensors.jl, for use with DMRG and other tensor network applications. In practice, I want to form the Hamiltonian in terms of the spin-orbitals (so the sums would be over spin and spatial orbitals), since we want to get dominant eigenvectors of the full many-body quantum chemistry Hamiltonian. The role of using Hartree-Fock is just to have a particular single particle basis which might be useful for DMRG (if I'm using the terminology correctly, in the molecular orbital (MO)/Hartree-Fock basis, the MPS would be a product state in the non-interacting limit, though the Hamiltonian is very nonlocal and away from that limit there are better bases to use, which is an active area of research).

If I'm following your description correctly, to get hpq and gpqrs all I need to do is:

using Fermi
using LinearAlgebra

@molecule {
  O        1.2091536548      1.7664118189     -0.0171613972
  H        2.1984800075      1.7977100627      0.0121161719
  H        0.9197881882      2.4580185570      0.6297938830
}

@set basis sto-3g

wfn = @energy rhf

hα = Diagonal(wfn.orbitals.eps)
gα = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals)["ERI"]

# or interleave, depending on the convention of the
# orbital ordering.
h = cat(hα, hα; dims=(1, 2))
g = cat(gα, gα; dims=(1, 2, 3, 4))

(I'm sure I could have guess that if I knew more about Hartree-Fock and dug into the code a bit more.)

Does that look right to you?

[gustavojra]

There are just a few things I would reinforce

1. Is hpq = fpq? In some Hamiltonian formulations, it is convenient to use the Fock operator as the one-electron part. This matters because the two-electron part will change depending on what you bake into hpq. I think the most general second quantized Hamiltonian looks like this:

If hpq is just the core Hamiltonian, then you get it also from the integral helper

aoints = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals)
hα = aoints["T"] + aoints["V"]

2. The two electron part that you get from Fermi is not antisymmetrized. So you may need to do it yourself if your Hamiltonian is defined in that way. Another thing to notice is that the integrals are in the so called "chemists' notation". If you want to convert it to "physicists' notation" you can simply do

permutedims(g, (1,3,2,4))

Finally, the spin-orbital array for the two-electron operator is not as simple as in the one electron case.

For the integral <pq|rs> (Physicists notation) you have 4 non-zero "blocks"

<αα|αα>
<αβ|αβ>
<βα|βα>
<ββ|ββ>

I think it is easier to construct it piece by piece

gα = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals)["ERI"]
# Convert to physicists notation
permutedims!(gα, (1,3,2,4))
# Create full g
N = size(g,1)
g = zeros(2*N, 2*N, 2*N, 2*N)
α = 1:N
β = (N+1):2*N
g[α, α, α, α] .= gα
g[α, β, α, β] .= gα
g[β, α, β, α] .= gα
g[β, β, β, β] .= gα

Okay, let's say I want the first Hamiltonian from that figure (not the antisymmetrized one) and I want to use physicists' notation. The full code would look like:

using Fermi
using LinearAlgebra
using TensorOperations # Just for a final sanity check

@molecule {
  O        1.2091536548      1.7664118189     -0.0171613972
  H        2.1984800075      1.7977100627      0.0121161719
  H        0.9197881882      2.4580185570      0.6297938830
}

@set basis sto-3g

wfn = @energy rhf
aoints = Fermi.Integrals.IntegralHelper(orbitals=wfn.orbitals)
hα = aoints["T"] + aoints["V"] 
gα = aoints["ERI"]

# Convert g to Physicists notation
permutedims!(gα, (1,3,2,4))

# Construct the full spin-orbital arrays
h = cat(hα, hα; dims=(1, 2))
N = size(gα,1)
g = zeros(2*N, 2*N, 2*N, 2*N)
α = 1:N
β = (N+1):2*N
g[α, α, α, α] .= gα
g[α, β, α, β] .= gα
g[β, α, β, α] .= gα
g[β, β, β, β] .= gα

# Let's make a sanity check by computing HF energy

# Define occupied ranges
αocc = 1:5
βocc = 8:12

# One particle density matrix for HF is just 1.0 for occupied orbitals and zero elsewhere (because we are working in the HF orbitals space)
Γ = [i == j && (i in βocc || i in αocc) ? 1.0 : 0.0 for i = 1:14, j = 1:14]

# The two paricle density can be build from the 1P one
@tensor Γ2[p,q,r,s] := Γ[p,r]*Γ[q,s] - Γ[p,s]*Γ[r,q]

Now contract the densities with the Hamiltonian
@tensor E = h[p,q]*Γ[p,q] + 0.5*g[p,q,r,s]*Γ2[p,q,r,s]

# If you check E + aoints.molecule.Vnuc (nuclear repulsion) you should get the same as wfn.energy -74.96500289468511

[mtfishman]

Wow thanks @gustavojra, that's a great write-up! The story is more involved than I expected, there are all sorts of subtle differences in terminology and conventions between quantum chemists and physicists. I think I will wrap this up in a small package (maybe ITensorFermi.jl) to make the pipeline from a HF solution in Fermi.jl to our second quantized Hamiltonian representation in ITensors.jl easier.

[gustavojra]

Yes, I notice that every time I try to read a DMRG paper haha... Anyways, let me know if you need anything else.