# 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 <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 / 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 at 17:42
• @Tyberius yes, I wonder if there is a nice way to visualize this density in python. Jul 4 at 17:44
• Jul 6 at 15:15
• Multiwfn will visualize electron density as well as perform a bunch of other calculations Jul 11 at 23:11