Large calculation ao_eri seg fault

Hello,

I am attempting to do a calculation on a 573 atom water droplet in the cc-pvdz basis. I successfully converged an SCF calculation, but now I am trying to compute overlap integrals. I get a seg fault, however, when I try to compute them with ao_eri.

mints = psi4.core.MintsHelper(orbital_basis)
result = mints.ao_eri()

This gives me the following seg fault:
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
libgdma.so 00007F9F825540DC for__signal_handl Unknown Unknown
libpthread-2.27.s 00007F9F96A268A0 Unknown Unknown Unknown
libchemps2.so.3 00007F9F82E5B2CA __intel_avx_rep_m Unknown Unknown
core.cpython-38-x 00007F9F858FCFED _ZN3psi6linalg6de Unknown Unknown
core.cpython-38-x 00007F9F858B4C1F _ZN3psi6Matrix5al Unknown Unknown
core.cpython-38-x 00007F9F858B58E6 _ZN3psi6MatrixC1E Unknown Unknown
core.cpython-38-x 00007F9F85804B7C _ZN3psi11MintsHel Unknown Unknown
core.cpython-38-x 00007F9F8584B31B _ZN3psi11MintsHel Unknown Unknown
core.cpython-38-x 00007F9F833FE13F Unknown Unknown Unknown
core.cpython-38-x 00007F9F8332A0C2 Unknown Unknown Unknown
core.cpython-38-x 00007F9F83329F9B Unknown Unknown Unknown
core.cpython-38-x 00007F9F8319234E Unknown Unknown Unknown
python 0000556240E9AF76 PyCFunction_Call Unknown Unknown
python 0000556240E5885F _PyObject_MakeTpC Unknown Unknown
python3.8 0000556240EA6FA1 Unknown Unknown Unknown
python3.8 0000556240E1B77F Unknown Unknown Unknown
python 0000556240EA5A92 _PyEval_EvalCodeW Unknown Unknown
python 0000556240EA6754 PyEval_EvalCodeEx Unknown Unknown
python 0000556240F34EDC PyEval_EvalCode Unknown Unknown
python3.8 0000556240F34F84 Unknown Unknown Unknown
python3.8 0000556240F671F4 Unknown Unknown Unknown
python 0000556240E2F6E1 PyRun_FileExFlags Unknown Unknown
python 0000556240E2FAC6 PyRun_SimpleFileE Unknown Unknown
python3.8 0000556240E3098B Unknown Unknown Unknown
python 0000556240F69D19 Py_BytesMain Unknown Unknown
libc-2.27.so 00007F9F96644B97 __libc_start_main Unknown Unknown
python3.8 0000556240EF9E93 Unknown Unknown Unknown

I am using the version 1.4a2.dev1013+fb7f385 from conda.

I assume the crash is related to the very large size, given that the same calculation works fine for small and medium sized systems. Any ideas for a fix or workaround?

Thanks!
Josh

Should add that the seg fault occurs on the mints.ao_eri line. MintsHelper is initialized just fine.

Is this part of a larger script? Are there other calculations done and a clean() is required?
A minimal breaking example would be helpful. Can be without the coordinates.

I also don’t understand why you’re even using ao_eri if you want overlap integrals.

Am I missing something?

Hi @hokru and @jmisiewicz,

Thanks for the quick responses. This is part of a larger script that does density fitting (credit @andysim). Here is a minimal example that reproduces the error.

import numpy as np
import psi4
from scipy.stats import special_ortho_group

# read in args
xyzfile = sys.argv[1]
basisname = sys.argv[2]
auxbasis = sys.argv[3]

psi4.core.set_output_file('output.dat', False)

with open(xyzfile) as f:
    temp = f.readlines()[2:]

molstr = ' '.join(temp) 
molstr = molstr + "\n symmetry c1 \n no_reorient \n"
mol = psi4.geometry(molstr)

psi4.core.set_global_option('df_basis_scf', auxbasis)

orbital_basis = psi4.core.BasisSet.build(mol, "ORBITAL", basisname)
aux_basis = psi4.core.BasisSet.build(mol, "DF_BASIS_SCF", "", "JFIT", auxbasis)

zero_basis = psi4.core.BasisSet.zero_ao_basis_set()
mints = psi4.core.MintsHelper(orbital_basis)

#
# Form 3 center integrals (P|mn)
#
J_Pmn = np.squeeze(mints.ao_eri(
    aux_basis, zero_basis, orbital_basis, orbital_basis))

print("formed 3 center integrals!")
print("DONE")

If I run with this command,
python minimal.py 8mer.txt cc-pvdz cc-pvdz-jkfit
it runs to completion.

If, however, I run with this command,
python minimal.py droplet_191.txt cc-pvdz cc-pvdz-jkfit
it seg faults.

XYZ files are attached.droplet_191.txt (15.8 KB) 8mer.txt (1.3 KB)

A memory issue perhaps?
For your droplet I calculate the following demands (nbf x nbf x naux):

(4584*4584*22156*8)/1024**3 = 3468.7 GiB

Useful to print the basis set info

orbital_basis.print_out()
aux_basis.print_out()

It looks like you’re probably right. Here is the basis set info:

Basis Set: CC-PVDZ
    Blend: CC-PVDZ
    Number of shells: 2292
    Number of basis function: 4584
    Number of Cartesian functions: 4775
    Spherical Harmonics?: true
    Max angular momentum: 2

  Basis Set: CC-PVDZ-JKFIT
    Blend: CC-PVDZ-JKFIT
    Number of shells: 8022
    Number of basis function: 22156
    Number of Cartesian functions: 25021
    Spherical Harmonics?: true
    Max angular momentum: 3

Any ideas for a workaround for this?

Not sure… This is quite a large calculation even with psi4 itself. Actually disk/memory demands will be smaller with screening and permutational symmetry.

Maybe using the DFTensor class directly is an option? http://psicode.org/psi4manual/master/psi4api.html#psi4.core.DFTensor
An example can be found the psi4numpy MP2: https://github.com/psi4/psi4numpy/blob/master/Moller-Plesset/DF-MP2.py#L66
You want the Qso tensor, but even then you wont be able to fit it into memory.
Storing the files on disk and reading batches is the typical way. Or a Direct algorithm.

In any case giving the script more memory than the default is a good idea.
psi4.set_memory('X GB')

I have written a plugin borrows the concept from the original code of mints.ao_eri to save these integrals to disk and access them by other scripts, maybe it would help. Also, I’m not familiar with C++, so they probably have a lot of place can be improved. Note, although this can help you avoid the memory problem, your basis set is still too large, and it will cost a lot of time and disk space.

.cc is not allowed to be uploaded here, so I change it to .txt
aointegrals.txt (7.3 KB)

@joshrackers what’s your ultimate goal? Where do you want to consume these integrals? Problems of this size are usually formulated in “direct” mode, where the integrals aren’t ever in memory or disk simultaneously.

@lhes30412, thanks for sharing your code! I will have a look and see if it might be helpful.

Hi @susilehtola,

I want to consume these integrals along with the auxiliary metric matrix to fit the 1 particle density matrix to an auxiliary basis set. The full code is below. There definitely should be some way to formulate this in a “direct” way, but it’s not obvious to me. If you can come up with a good way, that would be great!

(again, code credit: @andysim)


import sys
import tempfile
import numpy as np
import psi4
from scipy.stats import special_ortho_group
import time


# read in args
xyzfile = sys.argv[1]
basisname = sys.argv[2]
auxbasis = sys.argv[3]


xyzprefix = xyzfile.split('.')[0]

psi4.set_memory('30 GB')
psi4.set_num_threads(8)
psi4.core.set_output_file('output.dat', False)

ang2bohr = 1.88973
bohr2ang = 1/ang2bohr

#necessary to skip the first two lines of standard xyz file format
with open(xyzfile) as f:
    temp = f.readlines()[2:]

molstr = ' '.join(temp) 
molstr = molstr + "\n symmetry c1 \n no_reorient \n"
mol = psi4.geometry(molstr)

E, wfn = psi4.energy('PBE0/{}'.format(basisname), return_wfn=True)
print("PBE0 energy: ", E)

psi4.core.set_global_option('df_basis_scf', auxbasis)

orbital_basis = wfn.basisset()
aux_basis = psi4.core.BasisSet.build(mol, "DF_BASIS_SCF", "", "JFIT", auxbasis)

zero_basis = psi4.core.BasisSet.zero_ao_basis_set()
mints = psi4.core.MintsHelper(orbital_basis)

#
# Check normalization of the aux basis
#
Saux = np.array(mints.ao_overlap(aux_basis, aux_basis))

#
# Form 3 center integrals (P|mn)
#
# THIS IS THE CURRENT PROBLEM
###
J_Pmn = np.squeeze(mints.ao_eri(
    aux_basis, zero_basis, orbital_basis, orbital_basis))
###

#
# Form metric (P|Q) and invert, filtering out small eigenvalues for stability
#
J_PQ = np.squeeze(mints.ao_eri(aux_basis, zero_basis, aux_basis, zero_basis))
evals, evecs = np.linalg.eigh(J_PQ)
evals = np.where(evals < 1e-10, 0.0, 1.0/evals)
J_PQinv = np.einsum('ik,k,jk->ij', evecs, evals, evecs)

#
# Finally, compute and print the fit coefficients.  From the density matrix, D, the
# coefficients of the vector of basis aux basis funcions |P) is given by
#
# D_P = Sum_mnQ D_mn (mn|Q) PQinv[P,Q]
#

D = np.array(wfn.Da()) + np.array(wfn.Db())

D_P = np.einsum('mn,Pmn,PQ->Q', D, J_Pmn, J_PQinv, optimize=True)

print("-----------")
print(D_P)

@joshrackers so… you just want to do an RI… And since you’re talking about the 1PDM, I assume you’re interested in the RI-J.

So, brief intro: the Coulomb matrix is J(uv) = (uv|st) P(st), where P is the density matrix. In the RI approximation the ERI (uv|st) is approximated as sum_{AB} (uv|A) (A|B)^{-1} (B|st). If you plug this in, you see that the Coulomb matrix becomes

J(uv) = (uv|st) P(st) = sum_{AB} (uv|A) (A|B)^{-1} (B|st) P(st).

Now, the trick is to notice that the sum over s and t is only in the last two terms, so you can build the vector gamma(B) = (B|st) P(st) by computing the (B|st) integrals one triplet at a time, and incrementing gamma.

Next, you can multiply gamma with the inverse overlap matrix
delta(A) = (A|B)^-1 gamma(B)
and at the last step you build the Coulomb matrix as
J(uv) = (uv|A) delta(A)

You can do a similar approach for exchange, but it’s a bit more tricky since you’re not contracting over all indices straight away. The algorithms for both cases have been presented by e.g. Weigend in https://doi.org/10.1039/B204199P.

Unfortunately, neither algorithm is not implemented in Psi4 at the moment. It shouldn’t be too hard to add them, though…