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()
Dtot.add(wfn.Db())
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()
Dtot.add(wfn.Db())
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)
def cube_to_radial(cubename):
#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)
(zs,values) = cube_to_radial("Dtot.cube")
(zs2,values2) = cube_to_radial("Dtot2.cube")
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()
```