# How to plot a Gaussian cube file with python/matplotlib?

I know there are programs that can automatically visualize .cube files (like VMD), but I'm trying to do it myself in python/matplotlib (or equivalent) to better understand what the numbers mean.

I understand that there are nx*ny*nz grid points in the file (and can be saved in an array fo dimensions [nx,ny,nz]), but I can't quite figure out how to visualize them or what kind of plot they correspond to (with respect to the modules available on matplotlib).

Part of a script that extracts info from cube file :

filename1 = 'density.cube'
with open(filename1, 'r') as f:
no_of_atoms, _, _, _ = lines[2].split()
no_of_atoms = int(no_of_atoms)
xdim, _, _, _ = lines[3].split()
xdim = int(xdim)
ydim, _, _, _ = lines[4].split()
ydim = int(ydim)
zdim, _, _, _ = lines[5].split()
zdim = int(zdim)

elements = [None] * no_of_atoms
atoms_x_coords_in_density = [None] * no_of_atoms
atoms_y_coords_in_density = [None] * no_of_atoms
atoms_z_coords_in_density = [None] * no_of_atoms

for _ in range(no_of_atoms):
(elements[_], __,
atoms_x_coords_in_density[_],
atoms_y_coords_in_density[_],
atoms_z_coords_in_density[_]) = lines[6 + _].split()
elements[_] = int(elements[_])
atoms_x_coords_in_density[_] = float(atoms_x_coords_in_density[_])
atoms_y_coords_in_density[_] = float(atoms_y_coords_in_density[_])
atoms_z_coords_in_density[_] = float(atoms_z_coords_in_density[_])

# Load the data, we need to remove the first 8 lines and the space after
str = ' '.join(file(filename1).readlines()[(6 + no_of_atoms):])
data = np.fromstring(str, sep=' ')
data.shape = (xdim, ydim, zdim)


So, the grid points are in the data array, but I'm not sure how to plot them. I've seen an example using the mayavi package, but that's also a black box solution (plus I couldn't get it installed for some reason -- the conda version seems broken). I'm trying to understand about what's happening under the hood.

Any help is appreciated.

• The data is electron density on a discretized volume. Each data point is the value of the electron density inside a cube of size dxdydz. The position of this cube is encoded in the arrangement of the data points. I.e. you can map the linear index of each data point to a point in 3 dimensional space. Three dimensional densities are often plotted as volume plots or as iso-surface plots. But 3D plots of this kind are not the strong suite of matplotlib and I wouldn't expect particularly good results, especially not straight out-of-the-box. Oct 11, 2021 at 9:57
• I think the plotly package can do isosurface plots. I don't quite understand what you mean by: "The position of this cube is encoded in the arrangement of the data points. I.e. you can map the linear index of each data point to a point in 3 dimensional space." Could you elaborate? Or point me to a suitable reference? Oct 11, 2021 at 16:44
• A good explanation of the cube format is given here, paulbourke.net/dataformats/cube, I hope that link makes it clear what I meant. The Gaussian website also provides a description, gaussian.com/cubegen. I do not clearly understand what your question is. Is all that you need a clear explanation of the cube format so that you can parse it ? Or do you not understand what kind of physical property the cube file represents ? Oct 11, 2021 at 19:36
• @HansWurst I might try to add an answer here when I find the time, but from the comments it looks like you already have most of whats needed for an answer, so I'll give you the first crack at it if you want it. It would also essentially address this prior question, so it seems like an explicit answer could be useful to a decent number of users.
– Tyberius
Oct 26, 2021 at 14:01
• @Tyberius I was able to figure it out. The plotly "volume plot" function essentially does it: plotly.com/python/3d-volume-plots .When I have the time I'll try to write up an answer. Oct 26, 2021 at 17:54

python3 -m ase.visualize.mlab file.cube

will visualize the grid data from file.cube. To plot from within your own python code, use the method ase.visualize.mlab.plot, which takes an ase atoms object, the grid data, and a list of floats to plot iso-surfaces at as parameters.