# Visualize electron density using pyscf

I found the following example of code which uses density functional theory to compute the electron density $$\rho$$:

    #!/usr/bin/env python
#
# Author: Qiming Sun <[email protected]>
#

import numpy
from pyscf import lib
from pyscf.dft import numint, gen_grid

"""
Gaussian cube file format
"""

def density(mol, outfile, dm, nx=80, ny=80, nz=80):
coord = mol.atom_coords()
box = numpy.max(coord, axis=0) - numpy.min(coord, axis=0) + 4
boxorig = numpy.min(coord, axis=0) - 2
xs = numpy.arange(nx) * (box / nx)
ys = numpy.arange(ny) * (box / ny)
zs = numpy.arange(nz) * (box / nz)
coords = lib.cartesian_prod([xs, ys, zs])
coords = numpy.asarray(coords, order="C") - (-boxorig)

nao = mol.nao_nr()
ngrids = nx * ny * nz
blksize = min(200, ngrids)
rho = numpy.empty(ngrids)
for ip0, ip1 in gen_grid.prange(0, ngrids, blksize):
ao = numint.eval_ao(mol, coords[ip0:ip1])
rho[ip0:ip1] = numint.eval_rho(mol, ao, dm)
rho = rho.reshape(nx, ny, nz)

with open(outfile, "w") as f:
f.write("Density in real space\n")
f.write("Comment line\n")
f.write("]" % mol.natm)
f.write(" .8f .8f .8f\n" % tuple(boxorig.tolist()))
f.write("] .8f .8f .8f\n" % (nx, xs, 0, 0))
f.write("] .8f .8f .8f\n" % (ny, 0, ys, 0))
f.write("] .8f .8f .8f\n" % (nz, 0, 0, zs))
for ia in range(mol.natm):
chg = mol.atom_charge(ia)
f.write("%5d %f" % (chg, chg))
f.write(" .8f .8f .8f\n" % tuple(coord[ia]))
fmt = " .8e" * nz + "\n"
for ix in range(nx):
for iy in range(ny):
f.write(fmt % tuple(rho[ix, iy].tolist()))

if __name__ == "__main__":
from pyscf import gto, scf
from pyscf.tools import cubegen

mol = gto.M(atom="H 0 0 0; H 0 0 1")
mf = scf.RHF(mol)
mf.scf()
cubegen.density(mol, "h2.cube", mf.make_rdm1())


I wonder however about the nature of this output rho. So if I was to visualize this electron density how would I plot it? Is anyone familiar with this package? It seems to be an array with three columns $$\rho$$, each of length 80. But what does one point then correspond to?

• Welcome to the site! Just to clarify, is your question about how to the plot the density (ie how to do this in some programming language) or how to interpret the Cube file format?
– Tyberius
Jul 4, 2021 at 17:42
• @Tyberius yes, I wonder if there is a nice way to visualize this density in python. Jul 4, 2021 at 17:44
• Jul 6, 2021 at 15:15
• Multiwfn will visualize electron density as well as perform a bunch of other calculations Jul 11, 2021 at 23:11

This is just the Gaussian cube format. It's essentially a voxel dump of the wave function evaluated on a grid. You can find some documentation at http://paulbourke.net/dataformats/cube/ and https://h5cube-spec.readthedocs.io/en/latest/cubeformat.html .

Most electronic structure programs are able to generate Gaussian cube files. Several molecular viewers can also do so. For instance, I think Avogadro, Jmol and IQmol all support visualizing from Gaussian cube format.

Cube files can be parsed into a 3D scalar grid (i.e. a numpy array) directly. These can then by visualized by applying the marching cubes approach (scikit-image has an implementation) to generate a mesh, and this mesh can be shown in interactive 3d using meshplot.

For example, to parse cube files, in the past I've used:

def parse_cube(filename):
#from: https://github.com/psi4/psi4numpy/blob/6ed03e715689ec82bf96fbb23c1855fbe7835b90/Tutorials/14_Visualization/vizualize.ipynb
""" Parses a cube file, returning a dict of the information contained.
The cubefile itself is stored in a numpy array. """
with open(filename) as fp:
results = {}

# skip over the title

natoms = int(origin)
results['minx'] = minx = float(origin)
results['miny'] = miny = float(origin)
results['minz'] = minz = float(origin)

numx = int(infox)
incx = float(infox)
results['incx'] = incx
results['numx'] = numx
results['maxx'] = minx + incx * numx

numy = int(infoy)
incy = float(infoy)
results['incy'] = incy
results['numy'] = numy
results['maxy'] = miny + incy * numy

numz = int(infoz)
incz = float(infoz)
results['incz'] = incz
results['numz'] = numz
results['maxz'] = minz + incz * numz

atnums = []
coords = []
for atom in range(natoms):
atnums.append(int(coordinfo))
coords.append(list(map(float, coordinfo[2:])))
results['atom_numbers'] = np.array(atnums)
results['atom_coords'] = np.array(coords)

data = np.array([ float(entry) for line in fp for entry in line.split() ])
if len(data) != numx*numy*numz:
raise Exception("Amount of parsed data is inconsistent with header in Cube file!")
results['data'] = data.reshape((numx,numy,numz))

return results


Now load up a cube file:

cube = parse_cube('ESP.cube')
print(cube['data'].min(), cube['data'].max())


I print the min/max because I want to know some range within which to set the cutoff for marching cubes. Here the cutoff is 0.07, and the spacing (i.e. gridsize) is 0.3 angstroms. Note that the spacing is hardcoded into the cube file - one would normally set it when calculating the ESP.

from skimage import measure
vert, faces, norm, values= measure.marching_cubes(cube['data'],
0.07,
spacing=(0.3,0.3,0.3))


finally, draw this mesh with meshplot:

import meshplot as mp
vert = vert
mp.plot(vert, faces)


And the final result looks a bit like this: scikitimage: https://scikit-image.org/