用 Wolfram 语言绘制电子轨道

化学研究中可能经常需要绘制电子轨道，来描述原子或分子中电子的波函数。通常，它们是通过电子结构软件（如高斯程序 Gaussian）以多维数据集文件（Cube 文件）形式输出的。这些文件包含三维网格中给定轨道的体数据。

实现多维数据集文件可视化的应用程序有很多（如 VMD 或 GaussView），但我在这里想利用 Mathematica 的功能来轻松地合并图形, 以及使用它的过程自动化能力来高效地创建动画中的帧。

首先，我们需要一个从多维数据集文件中提取数据的函数。在这个过程中, 我们将创建一个 XYZ 文件的文本，这个格式也是由高斯开发的。函数 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}]

用 Wolfram 语言绘制电子轨道相关推荐

Wolfram 语言之父 Stephen Wolfram ：编程的未来
以后说到编程,我们想到的不单单是程序员了.未来的生活与计算机紧密相连,编程的未来与我们息息相关.用计算机完成我们所有感兴趣的事可能成为现实. 作者 |Nick Heath 译者 |弯月,责编 | ma ...
工程系的学生为什么要学 Wolfram 语言而不是 Matlab
Matlab 类似于 Fortran 和 C 代码,估计一般的小朋友是没有兴趣学的.Mathematica 是函数式编程,当然也支持过程式编程,还有对象编程.基于规则等的编程.更绝的是支持自然输入,你 ...
R语言绘制不一样的条形图
绘制条形图的方法有很多,这里介绍如果用R语言绘制一个不一样的条形图准备数据,这里为了方便,我们使用已经存在于gcookbook包中的一个数据集首先需要调用该包,如果该包不存在,可以使用下面的方式安 ...
R语言绘制生存曲线图
R语言绘制生存曲线图 KMunicate是支持按照Morris等人的KMunicate研究推荐的方式生成Kaplan-Meier图. 1958年,Edward L. Kaplan 和Paul Meie ...
R语言绘制Bump Chart
R语言绘制Bump Chart的小示例 # install.packages("tidyverse") # install.packages("ggbump") ...
R语言绘制环形树状图
R语言绘制环形树状图 1.主要用到dendextend和circlize包绘图: library(dendextend) library(circlize)# 距离矩阵 d <- dist(US ...
R语言绘制二维密度图
R语言绘制二维密度图二维密度图显示了两个数值变量之间的关系,一个在x轴上表示,另一个在Y轴上表示,与散点图类似,然后计算二维空间中特定区域内的观测数,并用颜色梯度表示.二维密度图有几种类型,以下主要 ...
R语言绘制带聚类树的堆叠柱形图
R语言绘制带聚类树的堆叠柱形图聚类树与柱形图结合,即可反映样本或分组间的相似性,又能展示样本内的元素组成信息. 例如下图是一个在扩增子测序微生物群落分析中常见的统计图类型,在测序公司给的报告中通常都 ...
R语言绘制线图（line）实战
R语言绘制线图(line)实战目录 R语言绘制线图(line)实战 #仿真数据 #基础线图
R语言绘制空白图实战
R语言绘制空白图实战目录 R语言绘制空白图实战 #绘制空白图1 #绘制空白图2 #绘制空白图3

用 Wolfram 语言绘制电子轨道

用 Wolfram 语言绘制电子轨道相关推荐

最新文章

热门文章