Cannot access MO ERI of heavy atoms in large basis set

1

Hi all,

I’m trying to access MO ERI for the list of atoms [Be, Ne, Ca, Zn, Xe] at HF/UGBS.
I use the following script

memory 2000 mb
import numpy as np

molecule at {
  Ca
}

set basis ugbs
energy('scf')

scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)

#Matrix of MO coefficients and AO ERI tensor

Matrix_C = scf_wfn.Ca_subset("AO", "ALL")
mints = psi4.core.MintsHelper(scf_wfn.basisset())
Matrix_I = mints.ao_eri()

#Transforming AO ERI to MO ERI manually

mo_I = np.asarray(mints.mo_transform(Matrix_I, Matrix_C, Matrix_C, Matrix_C, Matrix_C))

#Calculating the total electron-electron repulsion manually

ndocc = scf_wfn.nalpha()
Vee = 0
for i in range(ndocc):
	for j in range(ndocc):
		Vee += 2*mo_I[i][i][j][j] - mo_I[i][j][j][i]

print(Vee)
print()

#Total electron-electron repulsion calculated by the Psi4 code

print(list(scf_wfn.variables().values())[-1])

It works for Be, Ne, and Ca, but crashes without referring to any particular problem for Zn and Xe.
I suspect that the AO ERI tensor for Zn and heavier atoms is too large to be stored in RAM (my machine is freezing and the program stops at the point of creating Matrix_I).

Does anyone know how to solve this problem?
Thanks in advance

P.S. The electron-electron repulsion is calculated by the Psi4 code without any issues (the last string of the code), but I cannot access it manually. Isn’t the similar procedure of transferring AO ERI to MO ERI utilized by the Psi4 code?

That’s the basis set I used
ugbs.txt (26.8 KB)

2

No, this procedure is not used by Psi4. In fact, the “ao_eri” function is specifically labelled as “not recommended for large systems.”

Since I don’t know your background, I’m going to explain this at a level suitable for a first-year graduate student.

We’re discussing the Hartree-Fock method. Hartree-Fock requires building the Fock matrix and computing the current energy. Computing the Fock matrix is a harder problem, so let’s start there. According to the Psi4Numpy tutorial, the two-electron part of the Fock matrix looks like the tensor contraction

    J = np.einsum('pqrs,rs->pq', I, D, optimize=True)
    K = np.einsum('prqs,rs->pq', I, D, optimize=True)

…and all of Hartree-Fock’s two-electron integral needs can be expressed in terms of these two intermediates.

That’s a way to compute the J and K, but it’s not the most efficient way. There’s been a great deal of quantum chemistry research devoted to finding better ways to build the J and K intermediates, so we can do Hartree-Fock better and with less memory. What you’re proposing is the conceptually simplest way, but also the most memory intensive.

For example, you can use the density fitting approximation, where the two-electron integrals pqrs are written as pqQ, rsQ, where Q is the index of an “auxiliary orbital.” In the limit that the auxiliary orbitals span the whole space, this is identical to the normal integrals. In more realistic cases, the energy isn’t exactly the same, but it’s close enough for 99% of non-benchmark purposes, and it greatly reduces the memory cost, because you’re down from norb^4 to norb^2*naux, and naux < norb^2. Psi uses this approximation by default.

For a very slow but memory efficient approach, you can just re-compute the two-electron integrals as needed. In principle, you only need one integral at a time in this approach, but in practice, we compute a small cluster of integrals at once. This is called the integral direct algorithm.

There are other very sophisticated ways to compute these quantities, some more approximate, some less approximate. This is one of the areas of the code that I don’t know very well, so I’d need to consult another developer for details on how J and K are constructed.

What exactly are you trying to accomplish? There’s likely a way to do it using the much more efficient machinery that Psi actually uses to compute J and K. If your needs go beyond Hartree-Fock, then all bets are off.

1 Like
3

Thank you very much for your very detailed response!

What I actually need is to get values of self-repulsion Coulomb integrals [i][i][i][i] in MO.
So I guess the integral direct algorithm is the one needed?

4

This kind of two-electron integral manipulation is outside my area of expertise. I’ve passed this onto other developers, who may know this part of the code better than I do.

1 Like