用 Wolfram 语言绘制电子轨道
发布时间
阅读量：
阅读量
化学研究中常常需要描绘原子或分子内电子的波函数行为。这些描述通常借助于电子结构软件（例如高斯程序 Gaussian），它们会以多维数据集文件（Cube files）的形式输出结果。每个 Cube 文件都包含了特定轨道在三维网格中的数据表示。
多种多维数据集可视化的应用程序已知存在（如VMD或GaussView），但基于Mathematica的强大功能体系，在这里我计划运用其特性来便捷地整合图形数据，并借助其自动化处理功能实现动画帧的高效生成。
为了从多维数据集文件中读取数据并进行分析, 我们将编写一个名为 OutForm 的函数. 在这一过程中, 我们将生成一种XYZ格式的数据描述文档, 这种格式最初由Gaussian公司设计. OutForm 功能类似于其他编程语言中的printf指令.
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}
];如果需要创建多维数据集文件，可以使用以下函数：
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];]接下来，我们需要用该函数来绘制轨道：
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
{cubedata,xg,yg,zg,xyz,header}= ReadCube["cys-MO35.cube"];然后通过下式绘图：
CubePlot[{cubedata, xg, yg, zg, xyz}]在制作一个动画文件的过程中, 自然要求所有图像具备一致的视图参数, 包括统一的视角（ViewAngle）、观察点（ViewPoint）以及观察中心（ViewCenter）。一旦我们将这些参数配置好并传递给CubePlot函数后, Show函数会自动获取这些设置。
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 的任何选项。
CubePlot[{cubedata, xg, yg, zg, xyz}, vp,
ContourStyle -> {Texture[ExampleData[{"ColorTexture", "Vavona"}]],
Texture[ExampleData[{"ColorTexture", "Amboyna"}]]},
Contours -> {-.015, .015}]
全部评论 (0)
还没有任何评论哟~