Group Abstract Group Abstract

Message Boards Message Boards

Plotting electronic orbitals using Mathematica

GROUPS:
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 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
Answer
11 months ago
Have you considered turning ReadCube and WriteCube into an Import/Export filter?
POSTED BY: Robert Collyer
Answer
11 months ago
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
Answer
11 months ago
Ah, I was using the wrong terminology, Import/Export converter. That should get you started.
POSTED BY: Robert Collyer
Answer
11 months 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?
POSTED BY: Vitaliy Kaurov
Answer
11 months 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
POSTED BY: Jason Biggs
Answer
11 months ago
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
Answer
11 months ago
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
RebuildPacletData[]
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
Answer
11 months 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
Answer
11 months ago
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
Answer
11 months ago