13
$\begingroup$

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:
    lines = f.read().splitlines()
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.

Improve this question
$\endgroup$
11
  • 2
    $\begingroup$ 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. $\endgroup$
    – Hans Wurst
    Oct 11, 2021 at 9:57
  • $\begingroup$ 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? $\endgroup$
    – johnymm
    Oct 11, 2021 at 16:44
  • 2
    $\begingroup$ 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 ? $\endgroup$
    – Hans Wurst
    Oct 11, 2021 at 19:36
  • 1
    $\begingroup$ @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. $\endgroup$
    – Tyberius ♦
    Oct 26, 2021 at 14:01
  • 3
    $\begingroup$ @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. $\endgroup$
    – johnymm
    Oct 26, 2021 at 17:54

1 Answer 1

Reset to default
5
$\begingroup$

The Python package mayavi uses VTK for 3D visualization and offers thus several poptions for plotting volumetric data. ase has an interface to mayavi and facilitates IO of electronic structure grid formats such as cube files.

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Executing this errors for me with an error at ase/calculators/calculator.py:191: AttributeError: 'NoneType' object has no attribute 'title'. (Running on Ubuntu 20.04) $\endgroup$
    – Raven
    Oct 12, 2022 at 6:58
  • $\begingroup$ Actually, nvm. My file was named file.cub instead of file.cube and thus a different path within the program was taken. Renaming the file to have the proper file extension fixes above error. $\endgroup$
    – Raven
    Oct 12, 2022 at 7:02

You must log in to answer this question.

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