A function to plot molecular orbitals in Mathematica [closed]

I am making plots of electronic orbitals, and I wrote the following functions to plot them, as well as to read in the common cube file format from the program Gaussian. I know there are already programs to do this, but I find it easier to automate things in Mathematica, which is necessary when I go to make a movie file.

This is the main plotting function

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}]]]
];


These functions are necessary to read from and write to the .cube file format:

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;
];

{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;
lowerCorner = {0, 0, 0}; (* initialize *)
{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 =
Close[moltxt];
Close[moltxt];
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}
];

Module[{stream},
stream = OpenWrite[cubeFileName, FormatType -> FortranForm];
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];]


First read in the data (if anyone can tell me a convenient place on the web to put files that I can link to from here that isn't my own personal dropbox, please let me know and I'll edit this post accordingly),

{cubedata,xg,yg,zg,xyz,header}=


Then plot it via

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


The options for the viewing point or angle are fed 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}]


Any feedback would be most appreciated.

-
Wolfram Community is probably a better environment for these kinds of posts that deviate from specific Q&A form. StackExchange forums are much more rigorously Q&A oriented; flak is more likely, even when it would seem unfair. –  kirma Sep 24 at 21:06
This question appears to be off-topic because it isn't a question –  belisarius Sep 24 at 21:12
@JasonCasamir You can share this here, but the "right" way to do it is to ask a question that goes along with it (Jeopardy style) and then post your solution as an answer to the question. Next time, there's a checkbox when you compose a question that says "Answer my own question". You can click that to do this easily. If you don't do that, you'll have to wait 8 hrs (I think) before posting your own answer since you're a new user. –  rm -rf Sep 24 at 21:20
@belisarius Yes, the blog would be a good place as well. Jason, would you be interested in posting this as a blogpost here? Our last post was almost a year ago, so it can surely use some love... Let me know if you're interested and I'll set you up with an account and stuff. (btw, I remember you from a few other nifty plotting tools you've answered here before... I'm sure some of those would be great on the blog.) –  rm -rf Sep 24 at 21:27
I voted to close for the reason mentioned by Belisarius. I'd like to state for the record, though, that this doesn't reflect my opinion on the quality of your code, which I think is great. –  Sjoerd C. de Vries Sep 24 at 21:58