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.

5

It took me a while because the Psi4NumPy tutorial mentioned above did the trick and produced the correct density. Therefore, I didn’t have an urgent need to figure out the Cubeprop issue anymore :slight_smile:

So here’s what I found, at last:

  • the density does integrate out to the correct number of electrons in the limit of zero grid spacing. If I don’t get the correct electron number, I haven’t reached this limit yet :wink:
  • the different density values at common grid points were a result of a bug in my code from the original post. Specifically, the computation of density alters wfn.Da() so my computed density is correct the first time and twice as large the second time. The problem can obviously be solved by computing the density once and reusing it, but if you really want to calculate the density twice from the same wavefunction, one way not to mess things up is by playing with the psi4.core.Matrix <-> numpy.array conversions, as in the updated code below.
%matplotlib widget
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': [4.0,4.0,4.0],
                  'cubic_grid_spacing': [0.1,0.1,0.1]})
cubegen = psi4.core.CubeProperties(wfn)
Dtot = wfn.Da().np + wfn.Db().np
cubegen.compute_density(psi4.core.Matrix.from_array(Dtot), "Dtot")

#generate density cube number 2, with increased resolution
psi4.set_options({'cubic_grid_overage': [4.0,4.0,4.0],
                  'cubic_grid_spacing': [0.025,0.025,0.025]})
cubegen = psi4.core.CubeProperties(wfn)
Dtot = wfn.Da().np + wfn.Db().np
cubegen.compute_density(psi4.core.Matrix.from_array(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()

By the way, the value of CUBEPROP_ISOCONTOUR_THRESHOLD does not change the resulting density cube at all. The only change is in the comment line of the cube file, which reports the isosurface value required for the desired percentage of total density.

Hope this helps!