Basis functions evaluation on a user-defined grid

1

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)

#read the density from file

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)

#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 = 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[0]]

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