Message Boards Message Boards

Plotting electronic orbitals with Wolfram Language

Posted 11 years ago
I am reposting this here from the Stackexchange Mathematica blog so that more people might see it.  I'd be very happy to get some feedback on this plotting function.  If anyone can use the function, let me know how it works out for you, and if you'd recommend any changes.  If so, I can edit this post to have to most up-to-date version.

As a chemist it is often useful to plot electronic orbitals. These are used to describe the wave function of electrons in atoms or molecules. Typically, these are output from electronic structure software in the form of a cube file, first developed by Gaussian. These files contain volumetric data for a given orbital on a three-dimensional grid.

There exist many applications to visualize cube files, such as VMD or GaussView, but I wanted to take advantage of Mathematica‘s capability to easily combine graphics, as well as the ability to automate the process in order to efficiently create frames for a movie.

First off, we need a function to extract the data from the cube file. In the process, we will create the text for an XYZ file, a format also developed by Gaussian. The function OutForm is used here to mimic the printf function found in other programming languages.

 OutForm[num_?NumericQ, width_Integer, ndig_Integer,
    OptionsPattern[]] :=
   Module[{mant, exp, val},
    {mant, exp} = MantissaExponent[num];
    mant = ToString[NumberForm[mant, {ndig, ndig}]];
    exp = If[Sign[exp] == -1, "-", "+"] <> IntegerString[exp, 10, 2];
    val = mant <> "E" <> exp;
    StringJoin@PadLeft[Characters[val], width, " "]

ReadCube[cubeFileName_?StringQ] := Module[
   {moltxt, nAtoms, lowerCorner, nx, ny, nz, xstep, ystep, zstep,
    atoms, desc1, desc2, xyzText, cubeDat, xgrid, ygrid, zgrid,
    dummy1, dummy2, atomicNumber, atomx, atomy, atomz, tmpString,
   bohr2angstrom = 0.529177249;
   moltxt = OpenRead[cubeFileName];
   desc1 = Read[moltxt, String];
   desc2 = Read[moltxt, String];
   lowerCorner = {0, 0, 0};
   {nAtoms, lowerCorner[[1]], lowerCorner[[2]], lowerCorner[[3]]} =
    Read[moltxt, String] // ImportString[#, "Table"][[1]] &;
   xyzText = ToString[nAtoms] <> "\n";
   xyzText = xyzText <> desc1 <> desc2 <> "\n";
   {nx, xstep, dummy1, dummy2} =
    Read[moltxt, String] // ImportString[#, "Table"][[1]] &;
   {ny, dummy1, ystep, dummy2} =
    Read[moltxt, String] // ImportString[#, "Table"][[1]] &;
   {nz, dummy1, dummy2, zstep} =
    Read[moltxt, String] // ImportString[#, "Table"][[1]] &;
    {atomicNumber, dummy1, atomx, atomy, atomz} =
     Read[moltxt, String] // ImportString[#, "Table"][[1]] &;
    xyzText = If[Sign[lowerCorner[[1]]] == 1,
      xyzText <> ElementData[atomicNumber, "Abbreviation"] <>
       OutForm[atomx, 17, 7] <> OutForm[atomy, 17, 7] <>
       OutForm[atomz, 17, 7] <> "\n",
      xyzText <> ElementData[atomicNumber, "Abbreviation"] <>
       OutForm[bohr2angstrom atomx, 17, 7] <>
       OutForm[bohr2angstrom atomy, 17, 7] <>
       OutForm[bohr2angstrom atomz, 17, 7] <> "\n"];
    , {nAtoms}];
   cubeDat =
    Partition[Partition[ReadList[moltxt, Number, nx ny nz], nz], ny];
   moltxt = OpenRead[cubeFileName];
   headerTxt = Read[moltxt, Table[String, {2 + 4 + nAtoms}]];
   headerTxt = StringJoin@Riffle[headerTxt, "\n"];
   xgrid =
    Range[lowerCorner[[1]], lowerCorner[[1]] + xstep (nx - 1), xstep];
   ygrid =
    Range[lowerCorner[[2]], lowerCorner[[2]] + ystep (ny - 1), ystep];
   zgrid =
    Range[lowerCorner[[3]], lowerCorner[[3]] + zstep (nz - 1), zstep];
   {cubeDat, xgrid, ygrid, zgrid, xyzText, headerTxt}

If you need to create a cube file, then the following function can be used:
 WriteCube[cubeFileName_?StringQ, headerTxt_?StringQ, cubeData_] :=
   stream = OpenWrite[cubeFileName, FormatType -> FortranForm];
   WriteString[stream, headerTxt, "\n"];
   Map[WriteString[stream, ##, "\n"] & @@
      Riffle[ScientificForm[#, {3, 4},
          NumberFormat -> (Row[{#1, "E", If[#3 == "", "+00", #3],
               "\t"}] &), NumberPadding -> {"", "0"},
          NumberSigns -> {"-", " "}] & /@ #, "\n", {7, -1, 7}] &,
   cubeData, {2}];
Next we need the function to plot the orbital,
 CubePlot[{cub_, xg_, yg_, zg_, xyz_}, plotopts : OptionsPattern[]] :=
    Module[{xyzplot, bohr2picometer, datarange3D, pr},
     bohr2picometer = 52.9177249;
     datarange3D =
       bohr2picometer {{xg[[1]], xg[[-1]]}, {yg[[1]],
          yg[[-1]]}, {zg[[1]], zg[[-1]]}};
     xyzplot = ImportString[xyz, "XYZ"];
      ListContourPlot3D[Transpose[cub, {3, 2, 1}],
       Evaluate[FilterRules[{plotopts}, Options[ListContourPlot3D]]],
       Contours -> {-.02, .02}, ContourStyle -> {Blue, Red},
       DataRange -> datarange3D, MeshStyle -> Gray,
       Lighting -> {{"Ambient", White}}],
        FilterRules[{plotopts}, {ViewPoint, ViewVertical, ImageSize}]]]
Let’s look at an example. First we need to read in a cube file, download this cube file and place it in your base directory:  cys-MO35cube
{cubedata,xg,yg,zg,xyz,header}= ReadCube["cys-MO35.cube"];
Then plot it via
CubePlot[{cubedata, xg, yg, zg, xyz}]

When I want to create a movie file, I want all the images to have exactly the same ViewAngle, ViewPoint, and ViewCenter. When you give these options to CubePlot, it feeds them directly to the Show function
vp = {ViewCenter -> {0.5, 0.5, 0.5},
   ViewPoint -> {1.072, 0.665, -3.13},
   ViewVertical -> {0.443, 0.2477, 1.527}};

CubePlot[{cubedata, xg, yg, zg, xyz}, vp]

Finally, you can also give any options that normally go to ListContourPlot3D
CubePlot[{cubedata, xg, yg, zg, xyz}, vp,
ContourStyle -> {Texture[ExampleData[{"ColorTexture", "Vavona"}]],
   Texture[ExampleData[{"ColorTexture", "Amboyna"}]]},
Contours -> {-.015, .015}]

Many thanks to Daniel Healion for the ReadCube and WriteCube functions.
POSTED BY: Jason Biggs
21 Replies

Hi, Jason. I attended your talk in last WTC2020. It was about Molecule Fingerprint and Visualization. It was a good talk; it inspired me to use Mathematica for my research. Some speakers recommended visiting Wolfram Community site, and I found this interesting post since I also use Gaussian software in my research.

I tried your code, which I believe you wrote with a version before 12.1. I just copied, pasted, and run. Everything was good until I reached the following line

CubePlot[{cubedata, xg, yg, zg, xyz}]    

The line gave me the following message:

Show::gcomb: Could not combine the graphics objects in Show[<<1>>].

I hope this message is something to do with Version 12.1 that I am using. If it is, could you renew your code for Version 12.1? If not, could you give me some advice so that I can use it for my work?

Another thing. In his opening speech, Stephan Wolfram encouraged users to cite Mathematica in their work. I am willing to do so. Therefore, please advise me on how to cite your code when I plot the HOMO-LUMO wave function in my future manuscripts.

Thank you very much.

POSTED BY: Febdian Rusydi

Hello Febdian, I'm glad that you enjoyed the talk.

As of version 12.1 the functionality from CubePlot is available as part of the system via Import. See the documentation for the cube format for examples.

To reproduce the plots above you would use something like

	"", "Graphics3D",
	Contours -> {-0.02, 0.02},
	ContourStyle -> {Blue, Red},
	MeshStyle -> Gray, Mesh -> Automatic

enter image description here

It's great that you are going to use Mathematica for your research. I recommend this site as well as the Mathematica stack exchange as resources for help. Don't hesitate to ask if you can't figure something out - or if you find a bug.

See this support article for information about citing Wolfram products in your research.

POSTED BY: Jason Biggs

That is amazing. Thank you, Jason. It really inspires me.

And, thanks a lot for your offer. Sure, I will ask you in the future if I stumble some challenges that I can't overcome.

POSTED BY: Febdian Rusydi
Posted 8 years ago


I was wondering if you know a method how to manipulate the data in the cube file to obtain lets say the gradient or the laplacian of the density?

Regards, N

Posted 8 years ago


for All out there trying to plot Molecular or atomic orbitals with Wolfram Mathematica. You can use output files from Ab Intio Simulators GAMESS, Gaussian, NWChem, ADF, Molpro, Dalton, Jaguar, Orca, QChem and generate a Cube file using the Chemcraft software: from there, you can use the collection of functions that Jason Biggs and Daniel Healion kindly created and published here. Of course you can also write your own function in Mathematica or adjust the function that Jason and Daniel wrote. The latter is what I did, just before knowing I could convert with Chemcraft, more time investment, but great learning of the "Read" and "ImportString" functions in Wolfram Mathematica.

Thank you All involved, whishing you joy and happiness.


POSTED BY: Juan Tamara

enter image description here - you earned "Featured Contributor" badge, congratulations !

Dear @Jason Biggs, this is a great post and it has been selected for the curated Staff Picks group. Your profile is now distinguished by a "Featured Contributor" badge and displayed on the "Featured Contributor" board.

POSTED BY: Vitaliy Kaurov

Dear Vitaly - I'm glad it is appreciated, I really enjoy using Mathematica, for teaching and scientific purposes. I'd love to spend more time integrating it with various quantum chemistry packages. I've used it in the past to create input files for Gaussian and molpro.

(on an unrelated note, just last week I put in an application with Wolfram, so right now my fingers are crossed )

POSTED BY: Jason Biggs
Posted 8 years ago

Hi guys, I need some help here since I am starting with this new branch of Quantum chemistry. I do not understand the functions you have here, can you add more explanation please, like: - Do I have to have a cube File to use the functions? - How do I create a XYZ from other files, like a GAMESS-US Output?, can I use the function directly on a XYZ file generated with other software?. Thanks in advance for your help and time, wishing you joy and happiness.

POSTED BY: Juan Tamara

Juan, What specifically do you want to do? There are many programs that can create xyz files and convert between different quantum chemistry formats, or you can write a bash (or even Mathematica) script to do it yourself.

For example, if you look at the GAMESS-US output on this page, you find the final optimized geometry of the compound written

   ATOM   CHARGE       X              Y              Z
 SI         14.0   0.0000000000   0.0000000000  -0.3383854758
 H           1.0   0.6541232567  -1.1329747151   0.4374726899
 H           1.0   0.6541232567   1.1329747151   0.4374726899
 H           1.0  -1.3082465135   0.0000000000   0.4374726899

To make an XYZ file that mathematica can recognize, we need to get rid of the first few rows, then get rid of the column with the atomic masses, to make it look like this:

     SI            0.0000000000   0.0000000000  -0.3383854758
     H             0.6541232567  -1.1329747151   0.4374726899
     H             0.6541232567   1.1329747151   0.4374726899
     H            -1.3082465135   0.0000000000   0.4374726899

You could do this using awk pretty easily, and then import the results to Mathematica, or you could use this (inelegant) solution,

           ATOM   CHARGE       X              Y              Z
         SI         14.0   0.0000000000   0.0000000000  -0.3383854758
         H           1.0   0.6541232567  -1.1329747151   0.4374726899
         H           1.0   0.6541232567   1.1329747151   0.4374726899
         H           1.0  -1.3082465135   0.0000000000   0.4374726899","Table"][[4;;,{1,3,4,5}]],"Table"],"XYZ"]

which will give you this output,

enter image description here

a 3D modely you can combine with other Graphics3D objects as you choose.

For example, if you want to combine it with an atomic orbital, you need that in some plottable format (like a cube file, or expressions for the orbitals).

POSTED BY: Jason Biggs
Posted 8 years ago

Dear @Jason Biggs, I am trying to do the graphic of a molecule and Molecular orbitals together, just the way you did in the post, but with other starting files (GAMESS Output). Searching on the internet I just learned, what you said: The cube file contains the information of one atomic/molecular orbital and position of atoms ( The XYZ file contains the positions of the atoms in the molecule.

The point is that I am not using Gaussian and I don´t have it available to extract the information of the atomic/molecular orbitals from my simulations output in GAMESS and pack it in the cube file. I would need some kind of converter (to cube files), or an equivalent file to get the molecular orbital information, do you know of any?. I know some softwares like Molden or Avogadro, and only know a few about their plotting capabilities but I don´t know if I can get a cube file or similar from them. Thank you for your time and courtesy. Wishing you joy and happiness

POSTED BY: Juan Tamara

Dear Juan, Given time, we could probably work out how to plot the orbitals directly from the information in the GAMESS output file, but that could take a while. Basically if you know the basis set and you know the Hartree-Fock basis occupation numbers then you could bootstrap your way there.

But I just googled "gamess cube file" and the first couple of links look promising. Once you get the cube file, if you have trouble using the code here just let me know.

POSTED BY: Jason Biggs
Posted 8 years ago

I found something that might be equivalent to a cube file export for GAMESS, Thank you again. Wishing you all joy and happiness.

POSTED BY: Juan Tamara
What a coincidence, we just have published a related Demonstration:

Some Examples of Molecule Orbitals 

by Guenther Gsaller and Norbert Mueller from Institute of Organic Chemistry & Johannes Kepler University in Austria. Quoting:

Molecule orbitals are formed when atomic orbitals overlap. Mathematically, this is represented by a linear combination of atomic orbitals (as in the LCAO-MO method). Some possible classifications of orbitals are bonding and antibonding and Pi- and Sigma-symmetry.
This Demonstration shows the basic characteristics for a set of six molecules: the label, the description, the number of electrons in the chosen molecular orbital, and a 3D view of the probability density (with boundary surface, phase-coloring included) and also a ball and stick model for each example.

POSTED BY: Vitaliy Kaurov
Posted 11 years ago
Super cool! I really like this plotting functionality. I have a couple of questions:
  1. Can you explain a little more about the format for the cube file (what the columns mean and what the data represents)? Also, what are the specifications for the orbital in "cys-MO35.cube", and what molecule is it?
  2. Can you use your WriteCube function to write a cube file for some molecule? For instance, could you reproduce how you would go about writing "cys-MO35.cube"? If this enables you to visualize the orbitals of any molecule that would be insanely awesome!
Also, one style-oriented suggestion I can offer is that if you make the orbitals slightly transparent, you will always be able to see the underlying molecule in place throughout your movie:
CubePlot[{cubedata, xg, yg, zg, xyz}, vp,
ContourStyle -> {{Opacity[0.5],
    Texture[ExampleData[{"ColorTexture", "Vavona"}]]}, {Opacity[0.5],
    Texture[ExampleData[{"ColorTexture", "Amboyna"}]]}}]
POSTED BY: Tal Einav
Robert Collyer, I'm trying to implement your suggestion but I've run into a snag.  It seems that if I make a call to `ElementData` during the `Import` convertor program, then I cannot use `ElementData` in another context afterwards.   I get an `InputStream` error.  Here I've got it down to the simplest example that reproduces the error.

Here is where it fails:
 In[1]:= Export["testfile.txt", "hello"];
         ImportExport`RegisterImport["TestFormat", TestFormat`TestFormatImport];
         TestFormat`TestFormatImport[filename_String, options___] :=
           Module[{temp},(*Do some stuff,all works just fine*)
            temp = ElementData["C", "AtomicNumber"];
            {"Test" -> temp}];
         Import["testfile.txt", {"TestFormat", "Test"}]
         ElementData["F", "AtomicNumber"]
Out[4]= 6

        During evaluation of In[1]:= General::openx: InputStream[C:\Users\myusername\AppData\Roaming\Mathematica\Paclets\Repository\ElementData-7.0.28\Data\ElementData.wdx,126] is not open. >>

        During evaluation of In[1]:= BinaryRead::openx: InputStream[C:\Users\myusername\AppData\Roaming\Mathematica\Paclets\Repository\ElementData-7.0.28\Data\ElementData.wdx,126] is not open. >>

Out[5]= $Failed

But if, after I've defined the convertor program but before I invoke it, I call `ElementData` then I get no error message:
 In[1]:= Export["testfile.txt", "hello"];
         ImportExport`RegisterImport["TestFormat", TestFormat`TestFormatImport];
         TestFormat`TestFormatImport[filename_String, options___] :=
           Module[{temp},(*Do some stuff,all works just fine*)
            temp = ElementData["C", "AtomicNumber"];
            {"Test" -> temp}];
         ElementData["H", "AtomicNumber"]
         Import["testfile.txt", {"TestFormat", "Test"}]
         ElementData["F", "AtomicNumber"]

Out[4]= 1

Out[5]= 6

Out[6]= 9

Any help would be most appreciated.
POSTED BY: Jason Biggs
I can't seem to reproduce that effect. First, have you tried this on a fresh kernel?  Presuming you have, you could try rebuilding the paclets via
That likely isn't the source of the issue, but it may not hurt to try. Lastly, you could pre-build an index of atomic numbers and use that instead.
POSTED BY: Robert Collyer
Jason, thank you for posting this, - this is some nice work. I have a question. When I execute this lines:
{cubedata,xg,yg,zg,xyz,header}= ReadCube["cys-MO35.cube"];
I get some error messages. I do save linked file in the base directory of course preliminary. Can you confirm that these lines run fine on your system?
POSTED BY: Vitaliy Kaurov
It wasn't working because I don't understand how Mathematica extracts archive files, it seems to change the file name.  Since I can't seem to find a better place to host this 3MB file, I'm just linking to my dropbox.

It should work now
POSTED BY: Jason Biggs
Ah, I was using the wrong terminology, Import/Export converter. That should get you started.
POSTED BY: Robert Collyer
Not only have I not considered it, I didn't realize it was an option.  How does one go about writing custom Import and Export filters?
POSTED BY: Jason Biggs
Have you considered turning ReadCube and WriteCube into an Import/Export filter?
POSTED BY: Robert Collyer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract