# psi4.core.NumIntHelper - how to use?

Hi all,

I would like to evaluate a scalar function on a DFT grid and get a representation in the AO basis.

In the manual, I found psi4.core.NumIntHelper.potential_integral(), which seemed custom-made for this purpose. However, I don’t understand the expected data type:
List[psi4.core.Vector]

Indeed, I am not even 100% sure whether potential_integral() is what I want:

• Evaluate a function on the DFT grid.
• Pass the vector (array) with function values into the grid-integrate-function or grid object, and have Psi4 work out the AO representation.

An example would be very much appreciated.

### For the grid creation, I follow psi4numpy:

basis = psi4.core.BasisSet.build(mol, “ORBITAL”, psi4.core.get_global_option(“BASIS”))
sup = psi4.driver.dft.build_superfunctional(“PBE”, True)
Vpot = psi4.core.VBase.build(basis, sup, “RV”)
Vpot.initialize()
xs, ys, zs, ws = Vpot.get_np_xyzw()

### Evaluate my potential or other scalar function on the grid.

Vgrid = Eval_potential(xs, ys, zs)

### Get a matrix object nAO-by-nAO

V_AO_rep = ?

Any suggestions?

All the best,
Thomas

I don’t understand the question. If you have a scalar function to integrate on a grid, what would an “AO representation” be?

In KS-DFT, you can have an AO representation of the V matrix because it corresponds to the derivative of the exchange-correlation energy with respect to each matrix element of the electron density in the AO basis. Is something like that happening in your case?

Are you trying to integrate f * phi_mu * phi_nu for each mu, nu?

Yes, you got that right, I would like the integrals f * phi_mu * phi_nu for each mu, nu, where phi_mu are AOs.

So a matrix V having dimension nAOs times nAOs.

Let’s say V(r) is a local potential. Then the values of V(r) at the grid point are a “grid representation”, a bit like a DVR, and V is the AO-representation of V(r), right?

BTW, this has nothing to do with DFT—frankly, I am unfamiliar with the KS intricacies—I just have an external potential I cannot integrate analytically, so I would like to use the DFT grid to do the integrals numerically.

Yes, `psi4.core.NumIntHelper.potential_integral()` is the function you want. Its signature is designed for efficiency, not user-friendliness.

In Psi, points on the integration grid are organized into “blocks”. (See Section 2.2 of the Integration Grid tutorial in Psi4NumPy.) Element i in `List[psi4.core.Vector]` in the method signature is a psi4.core.Vector object that contains all function values for integration grid points in block i.

Does that answer the question, or is there still more to do? I only dabble in numerical integration, so I’ll probably need to call in another developer if you’re concerned about efficiency.

Thanks Jonathan.

I shall take a closer look at that tutorial and report back (may take a day or two).

For the time being efficiency is not an issue since the numbers of atoms, basis functions, and accordingly grid points is pretty small. More like a toy model. But is never hurts to plan for the future.

This is pretty much what the SAP code does in Psi4. The matrix elements are the same as for LDA.

Hello Jonathan and Susi, thanks so much for the hints. I think I am getting closer, but I am not quite there.

I apologize for the long message; just want to be very clear this time.

For the sake of the argument, let’s consider a particle in three-dimensional harmonic oscillator. [Yes, that can be done analytically, so it’s good for testing, and the harmonic oscillator will change eventually.]

I choose a GTO basis; Psi4 supplies the S and T matrix, one function call each.

Let’s pretend I cannot do the V matrix analytically, but I would rather like to use a numerical integral over a DFT grid.

Tutorial 4a is close, but I actually want the DFT grid suitable for integral evaluation, not a Cartesian grid suitable for plotting.
Tutorial 4b [LDA] is also close, but I don’t have a density, I just want the integrals
V_nu_mu = Int chi_nu(r) V(r) chi_mu(r) dr where chi_nu are the AOs.
On the grid that becomes a sum over grid points: Sum chi_nu(r_i) V(r_i) chi_mu(r_i) w_i, right?

As Jonathan pointed out, `psi4.core.NumIntHelper.potential_integral(List[psi4.core.Vector])` should do the trick.

I still wonder how and what `List[psi4.core.Vector]` refers to. After looking over the tutorials, I thought this must be the potential evaluated block-by-block:

``````Lst = []
for i in range(Vpot.nblocks()):
block=grid.blocks()[i]
v = np.ones(grid.npoints())
xbl=block.x().to_array()
ybl=block.y().to_array()
zbl=block.z().to_array()
# compute the function here
# for j in grid points do:
#     v[j] = 0.5*(xbl[j]**2 + ybl[j]**2 + zbl[j]**2)
#
# here we leave it as v = 1 to test against the S matrix
#
Lst.append(psi4.core.Vector.from_array(v))
``````

But that’s not it. Can someone please supply an example on how to use `psi4.core.NumIntHelper.potential_integral()`?

``````import psi4

# Create molecule
h2 = psi4.geometry("""
h 0.0 0.0 -1.0
h 0.0 0.0  1.0
symmetry c1
units bohr
nocom
noreorient
""")
# Create basisset
h2_basisset = psi4.core.Wavefunction.build(h2, "sto-6g").basisset()
# Create DFT grid. Despite the name, this is used by DFT but is not DFT-specific.
# First dict is for int options. Second dict is for string options.
# For demonstration purposes, I've given an awful number of radial points.
# Increase the number and watch how the numerical result converges to the analytical.
grid = psi4.core.DFTGrid.build(h2, h2_basisset, {"DFT_RADIAL_POINTS": 6}, {})
# Create numerical integration object.
numint = psi4.core.NumIntHelper(grid)

# Now we create the function values. We're choosing the identity operator, which should match the overlap integrals.
# function_values[i][j] is the value of the potential at point j of block i
function_values = [psi4.core.Vector.from_list( * block.npoints()) for block in grid.blocks()]
print("Numerical Integration Overlap")
print(numint.potential_integral(function_values).np)
print("Analytical Integration Overlap")
print(psi4.core.MintsHelper(h2_basisset).ao_overlap().np)
``````
2 Likes

First of all, again, thank you so much for all the help.

A few concluding remarks:

• Jonathan’s code solves the problem and does so in the most compact fashion.
• My code (see above) came close; it can be continued as follows:
``````numint = psi4.core.NumIntHelper(Vpot.grid())
Vmat=numint.potential_integral(Lst)
``````
• I used variants of Jonathan’s code to integrate a model potential, and I have versions with and without symmetry. This gets a bit too long for a post, but if anyone is interested, I’ll be happy to share.