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 <osirpt.sun@gmail.com>

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.