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[0] / nx)
    ys = numpy.arange(ny) * (box[1] / ny)
    zs = numpy.arange(nz) * (box[2] / 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[1], 0, 0))
        f.write("] .8f .8f .8f\n" % (ny, 0, ys[1], 0))
        f.write("] .8f .8f .8f\n" % (nz, 0, 0, zs[1]))
        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)
    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?

  • 1
    $\begingroup$ 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? $\endgroup$
    – Tyberius
    Commented Jul 4, 2021 at 17:42
  • 1
    $\begingroup$ @Tyberius yes, I wonder if there is a nice way to visualize this density in python. $\endgroup$ Commented Jul 4, 2021 at 17:44
  • 1
    $\begingroup$ See also: mattermodeling.stackexchange.com/q/6235/30 $\endgroup$ Commented Jul 6, 2021 at 15:15
  • 1
    $\begingroup$ Multiwfn will visualize electron density as well as perform a bunch of other calculations $\endgroup$
    – B. Kelly
    Commented Jul 11, 2021 at 23:11

2 Answers 2


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

        origin = fp.readline().split()
        natoms = int(origin[0])
        results['minx'] = minx = float(origin[1])
        results['miny'] = miny = float(origin[2])
        results['minz'] = minz = float(origin[3])

        infox = fp.readline().split()
        numx = int(infox[0])
        incx = float(infox[1])
        results['incx'] = incx
        results['numx'] = numx
        results['maxx'] = minx + incx * numx

        infoy = fp.readline().split()
        numy = int(infoy[0])
        incy = float(infoy[2])
        results['incy'] = incy
        results['numy'] = numy
        results['maxy'] = miny + incy * numy

        infoz = fp.readline().split()
        numz = int(infoz[0])
        incz = float(infoz[3])
        results['incz'] = incz
        results['numz'] = numz
        results['maxz'] = minz + incz * numz

        atnums = []
        coords = []
        for atom in range(natoms):
            coordinfo = fp.readline().split()
            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'], 

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:

enter image description here

meshplot: https://skoch9.github.io/meshplot/tutorial/

scikitimage: https://scikit-image.org/


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .