Normalization of density in Cubeprop

1

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()
2

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.

3

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…

4

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.