# Normalization of density in Cubeprop

I was wondering if anybody could tell me how exactly is the electron density in the cube file generated by `psi4.core.CubeProperties` normalized? My tests (see the notebook below) seem to indicate that (i) the density does not integrate out to the number of electrons, and (ii) if I extend the cubic grid (so that the old grid is a subset of the new one), the values at the common grid point change. Am I missing something?

``````%matplotlib notebook
import numpy as np
np.set_printoptions(precision=5, linewidth=200, threshold=2000, suppress=True)
import psi4
import matplotlib.pyplot as plt

# Set Psi4 & NumPy Memory Options
psi4.set_memory('2 GB')
psi4.core.set_output_file('output.dat', False)

geometry = psi4.geometry("""Ne 0.0 0.0 0.0
symmetry c1
""")
psi4.set_options({'basis': 'aug-cc-pvtz',
'df_basis_scf': 'aug-cc-pvtz-ri',
'df_basis_mp2': 'aug-cc-pvtz-ri'})
(e,wfn) = psi4.energy('HF',return_wfn=True)

#generate density cube number 1
psi4.set_options({'cubic_grid_overage': [5.0,5.0,5.0],
'cubic_grid_spacing': [0.2,0.2,0.2]})
cubegen = psi4.core.CubeProperties(wfn)
Dtot = wfn.Da()
cubegen.compute_density(Dtot, "Dtot")

#generate density cube number 2, with increased resolution
psi4.set_options({'cubic_grid_overage': [5.0,5.0,5.0],
'cubic_grid_spacing': [0.1,0.1,0.1]})
cubegen = psi4.core.CubeProperties(wfn)
Dtot = wfn.Da()
cubegen.compute_density(Dtot, "Dtot2")

def integrate_density(cubename):
#reads the complete cube file and performs a simple numerical integration, assuming only 1 atom
l = [line.rstrip() for line in open(cubename,"r")]
(xstart,ystart,zstart) = tuple([float(x) for x in l[2].split()[1:]])
xlen = int(l[3].split()[0])
ylen = int(l[4].split()[0])
zlen = int(l[5].split()[0])
xinc = float(l[3].split()[1])
yinc = float(l[4].split()[2])
zinc = float(l[5].split()[3])
lvalues = []
for s in l[7:]:
lvalues.extend([float(x) for x in s.split()])
integral = np.sum(lvalues) * xinc * yinc * zinc
return integral

Nel = integrate_density("Dtot.cube")
print("Approximate number of electrons 1:",Nel)
Nel = integrate_density("Dtot2.cube")
print("Approximate number of electrons 2:",Nel)

#reads the complete cube file and extracts only the values along the z axis, assuming only 1 atom
l = [line.rstrip() for line in open(cubename,"r")]
(xstart,ystart,zstart) = tuple([float(x) for x in l[2].split()[1:]])
xlen = int(l[3].split()[0])
ylen = int(l[4].split()[0])
zlen = int(l[5].split()[0])
xinc = float(l[3].split()[1])
yinc = float(l[4].split()[2])
zinc = float(l[5].split()[3])
lvalues = []
for s in l[7:]:
lvalues.extend([float(x) for x in s.split()])
lxyz = np.array(lvalues).reshape((xlen,ylen,zlen))
xmid = int((xlen - 1) / 2)
ymid = int((xlen - 1) / 2)
zend = zstart + zinc * (zlen -1)
zs = np.linspace(zstart,zend,zlen)
values = lxyz[xmid,ymid,:]
return (zs,values)

plt.close()
plt.plot(zs,values,'g',linestyle='-',label='Ne Density')
plt.plot(zs2,values2,'b',linestyle='-',label='Ne Density 2')
plt.legend(loc='upper right')
plt.xlim((-5.0,5.0))
plt.ylim((0.0,1.0))
plt.show()
``````

The cube file will only capture a fraction of the density, determined by the keyword `CUBEPROP_ISOCONTOUR_THRESHOLD`, which defaults to 85%.

This function is designed mainly for visualization, where 85% is perfectly reasonable. For numerically accurate integrations, this is not what you want. Depending on your use case, I may be able to recommend another way to get what you want.

Thank you @jmisiewicz for those hints. Yes, please do recommend a more accurate way. My goal is to generate HF and DFT densities on an arbitrary cross section or subset of the 3D space.

I saw this Psi4NumPy tutorial that looks very helpful:
https://github.com/psi4/psi4numpy/blob/master/Tutorials/04_Density_Functional_Theory/4a_Grids.ipynb
I don’t know if it can handle HF as well. `Cubeprop` was my first attempt because it seemed to me that it computes density the same way but already takes care of the 3D grid generation. Looks like I was wrong though…

I’ll be on vacation for a week and won’t be able to look at this until I get back.

The first thing I’d try is setting `CUBEPROP_ISOCONTOUR_THRESHOLD` and seeing what happens.