# Basis functions evaluation on a user-defined grid

Hi,
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 https://github.com/psi4/psi4/issues/1051 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
no_reorient
no_com
""")
p4mol = psi4.geometry(mol_desc)

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

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 == rho.shape):
print("Dim match")

#delta param
delta = 0.00001

basis_obj = psi4.core.BasisSet.build(p4mol, '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(lpos)

#print some info
blockopoints.print_out('b_info.txt')

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

print("phi matrix")
np.set_printoptions(formatter={'float': '{: 0.12f}'.format})
print(phi)
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):
Da[i,j]+=phi[p,i]*rho[p]*phi[p,j]*np.array(ws)[p]
``````

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