# Plotting electronic orbitals with Wolfram Language

GROUPS:
 Jason Biggs 18 Votes 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,     headerTxt,bohr2angstrom},   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]] &;   Do[    {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];   Close[moltxt];   moltxt = OpenRead[cubeFileName];   headerTxt = Read[moltxt, Table[String, {2 + 4 + nAtoms}]];   Close[moltxt];   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_] :=   Module[{stream},    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}];  Close[stream];]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"];     Show[xyzplot,       ListContourPlot3D[Transpose[cub, {3, 2, 1}],        Evaluate[FilterRules[{plotopts}, Options[ListContourPlot3D]]],        Contours -> {-.02, .02}, ContourStyle -> {Blue, Red},        DataRange -> datarange3D, MeshStyle -> Gray,        Lighting -> {{"Ambient", White}}],        Evaluate[        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 viaCubePlot[{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 functionvp = {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 ListContourPlot3DCubePlot[{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.
Answer
3 years ago
11 Replies
 Robert Collyer 2 Votes Have you considered turning ReadCube and WriteCube into an Import/Export filter?
Answer
3 years ago
 Jason Biggs 2 Votes 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?
Answer
3 years ago
 Robert Collyer 1 Vote Ah, I was using the wrong terminology, Import/Export converter. That should get you started.
Answer
3 years ago
 Jason, thank you for posting this, - this is some nice work. I have a question. When I execute this lines:ExtractArchive["cys-MO35cube.gz"];{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?
Answer
3 years ago
 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
Answer
3 years ago
 Jason Biggs 1 Vote 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"];         ImportExportRegisterImport["TestFormat", TestFormatTestFormatImport];         TestFormatTestFormatImport[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]= \$FailedBut 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"];         ImportExportRegisterImport["TestFormat", TestFormatTestFormatImport];         TestFormatTestFormatImport[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]= 1Out[5]= 6Out[6]= 9Any help would be most appreciated.
Answer
3 years ago
 Robert Collyer 1 Vote 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 viaRebuildPacletData[]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.
Answer
3 years ago
 Tal Einav 3 Votes Super cool! I really like this plotting functionality. I have a couple of questions: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?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"}]]}}]
Answer
3 years ago
 Vitaliy Kaurov 3 Votes 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.
Answer
3 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.
Answer
8 months ago
 Jason Biggs 3 Votes 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  COORDINATES OF ALL ATOMS ARE (ANGS) 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, ImportString[ ExportString[ ImportString[" COORDINATES OF ALL ATOMS ARE (ANGS) 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,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).
Answer
8 months 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 (http://paulbourke.net/dataformats/cube/). 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
Answer
8 months ago
 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.
Answer
8 months ago
 I found something that might be equivalent to a cube file export for GAMESS, https://www.youtube.com/watch?v=7n45abDIQ3Q. Thank you again. Wishing you all joy and happiness.
Answer
8 months ago
 Vitaliy Kaurov 2 Votes - 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.
Answer
8 months ago
 Jason Biggs 3 Votes 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 )
Answer
8 months ago
 Hi,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: http://www.chemcraftprog.com/index.html 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.Juan
Answer
8 months ago
 Hi, 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
Answer
1 month ago