Basis functions evaluation on a user-defined grid

I am trying to carry out the following task using Psi4 Python API: import a predefined grid (points and weights) and the electronic density mapped onto the grid. Basically, I need to evaluate the basis functions on the grid points. I have found this link and the 4b_LDA_kernel.ipynb tutorial very useful, and my code produces an (almost) reasonable PHI matrix. Eventually, I would like to “convert” the numerical density into its matrix representation with respect to the orbital basis. But the density matrix I obtain features a wrong number of electrons. Probably I have made some mistake when I evaluate the basis functions at grid points.

import numpy as np
import psi4
basis_set = 'cc-pvdz'
mol_desc = ("""
  O     0.00000000    -0.00001441    -0.34824012
  H    -0.00000000     0.76001092    -0.93285191
  H     0.00000000    -0.75999650    -0.93290797
 symmetry c1
p4mol = psi4.geometry(mol_desc)

#read the density from file

list = []
with open("density.txt","r") as f:
 for line in f:
  raw = line.split()
  for i in raw:
rho = np.array(list)

#load the grid and weights
pts = np.loadtxt('pts.txt')

xs = psi4.core.Vector.from_array(pts[0:-6, 0])
ys = psi4.core.Vector.from_array(pts[0:-6, 1])
zs = psi4.core.Vector.from_array(pts[0:-6, 2])
ws = psi4.core.Vector.from_array(pts[0:-6, 3])

#check if dimensions match
if (pts[0:-6,0].shape[0] == rho.shape[0]):
   print("Dim match")

#delta param
delta = 0.00001

basis_obj =, 'ORBITAL', basis_set)
basis_extents = psi4.core.BasisExtents(basis_obj,delta)

blockopoints = psi4.core.BlockOPoints(xs, ys, zs, ws,basis_extents)
npoints = blockopoints.npoints()
print("n points: %i" % npoints)
lpos = np.array(blockopoints.functions_local_to_global())

print("Local basis function mapping")

#print some info

funcs = psi4.core.BasisFunctions(basis_obj,npoints,basis_obj.nbf())
phi = np.array(funcs.basis_values()["PHI"])[:npoints, :lpos.shape[0]]

print("phi matrix")
np.set_printoptions(formatter={'float': '{: 0.12f}'.format})
Da = np.zeros((basis_obj.nbf(),basis_obj.nbf()),dtype=np.float_)
#compute Da mat
for i in range(basis_obj.nbf()):
  for j in range(basis_obj.nbf()):
     for p in range(npoints):

In attachment you can find the points and densitydensity.txt (502.8 KB)
pts.txt (1.6 MB). Thank you for help and support

1 Like