Visualize electron density using pyscf

Question

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)
    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?

Answer 1

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.

Answer 2

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
        fp.readline()
        fp.readline()

        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()
            atnums.append(int(coordinfo[0]))
            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:

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

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