█ 本文译自 Wolfram|Alpha 化学组开发人员 Jason.Biggs 在 Wolfram 社区发表的热点文章：Plotting electronic orbitals with Wolfram Language
化学研究中可能经常需要绘制电子轨道，来描述原子或分子中电子的波函数。通常，它们是通过电子结构软件（如高斯程序 Gaussian）以多维数据集文件（Cube 文件）形式输出的。这些文件包含三维网格中给定轨道的体数据。 实现多维数据集文件可视化的应用程序有很多（如 VMD 或 GaussView），但我在这里想利用 Mathematica 的功能来轻松地合并图形, 以及使用它的过程自动化能力来高效地创建动画中的帧。 首先，我们需要一个从多维数据集文件中提取数据的函数。在这个过程中, 我们将创建一个 XYZ 文件的文本，这个格式也是由高斯开发的。函数 OutForm 用于模拟其他编程语言中的 printf 函数。
代码语言：javascript复制 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}
];（滑动屏幕查看全部代码）
如果需要创建多维数据集文件，可以使用以下函数：
代码语言：javascript复制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];]（滑动屏幕查看全部代码）
接下来，我们需要用该函数来绘制轨道：
代码语言：javascript复制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}]]]
]; （滑动屏幕查看全部代码）
让我们来看一个实例。首先，复制并在浏览器中打开此链接： https://dl.dropboxusercontent.com/s/rdsxcnqudn1s76n/cys-MO35.cube ，这是一个多维数据集文件，将该文件保存到你的基目录下：
代码语言：javascript复制{cubedata,xg,yg,zg,xyz,header}= ReadCube["cys-MO35.cube"];然后通过下式绘图：
代码语言：javascript复制CubePlot[{cubedata, xg, yg, zg, xyz}]如果要制作一个动画文件，我们当然希望所有的图像都具有完全相同的视角（ViewAngle）、视点（ViewPoint）和视图中心（ViewCenter）。当将这些选项赋给 CubePlot 时, 它将直接提供给 Show 函数。
代码语言：javascript复制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]
最后，您还可以使用通常用于 ListContourPlot3D 的任何选项。
代码语言：javascript复制CubePlot[{cubedata, xg, yg, zg, xyz}, vp,
ContourStyle -> {Texture[ExampleData[{"ColorTexture", "Vavona"}]],
Texture[ExampleData[{"ColorTexture", "Amboyna"}]]},
Contours -> {-.015, .015}]
