Mathematica 处理 Gaussian 输出文件
type
status
date
slug
summary
tags
category
icon
password
Property
Apr 25, 2023 02:25 PM
URL
1. 可视化优化步骤
只做了 SCF 过程的能量,其他的懒得搞了,绘图范围也有点小问题,以后在搞吧
(* 定义文件路径 *) filepath = "http://pic.wxyh.top/c3.log";
上面的 filepath 是自己的文件的路径,需要自己自行修改,这里我改成我的 CDN 地址,这样联网就能用了
1.1 将 Gaussian 计算后的 log 文件的坐标提取出来
log2xyz[filepath_] := Module[ {log, indices, coordsList, dotline}, (*导入log文件*) log = StringSplit[Import[filepath, "Text"], "\n"]; (*找到"Standard orientation"所在行*) indices = Flatten[Position[log, s_String /; StringContainsQ[s, "Standard orientation"]]] + 5; dotline = "---------------------------------------------------------------------\ "; (*找到所有坐标*) coordsList = Map[ Function[ ind, Module[ {lines, data, coords}, lines = TakeWhile[log[[ind ;;]], Not@StringContainsQ[#, dotline] &]; data = N@Table[ToExpression /@ StringCases[line, RegularExpression[ "[-+]?[0-9]*\\.?[0-9]+(?:[eE][-+]?[0-9]+)?"]], {line, lines}]; coords = data[[All, {2, 4, 5, 6}]]; coords] ], indices]; coordsList ] (* coordsList 是从 Gaussian log 文件中提取的坐标 *) coordsList = log2xyz@filepath;
1.2 将找到的坐标处理成 mathematica 中的列表
molStringFromCoords[coords_List] := Module[ {symbol, coordStrings}, (*将原子编号转换为原子符号*) symbol[number_] := Entity["Element", number]["Abbreviation"]; (*将坐标转换为字符串*) coordStrings = StringRiffle[{symbol[Round@#1], " ", #2, " ", #3, " ", #4}, " "] & @@@ coords; (*创建分子字符串*) molString = ToString[Length[coords]] <> "\n\n" <> StringRiffle[coordStrings, "\n"]; (*导入字符串为 Molecule 对象*) mol = ImportString[molString, "XYZ"]; (*返回 Molecule 对象*) mol]
1.3 获取所有的能量
log = StringSplit[Import[filepath, "Text"], "\n"]; (*获取所有的能量*) pattern = RegularExpression["SCF Done: E\\([^)]+\\) =\\s+(-?\\d+\\.\\d+)"]; energies = ToExpression@StringCases[log, pattern -> "$1", 1] // Flatten;
1.4 绘图
(* 这是绘制 *) molPlots = Map[MoleculePlot3D[molStringFromCoords@#, ViewPoint -> {2.251449039826572`, 0.29525999152285654`, -2.5036371059859954`}] &, coordsList][[2 ;;]]; molPlotsThumbnails = Rasterize[#, ImageSize -> {90, 90}, Background -> None] & /@ molPlots; (*绘图样式*) plotOptions = { PlotRange -> {{0, 8}, {Min[energies] - 0.007, Max[energies] + 0.005}}, PlotStyle -> {Thickness[0.003], Red}, PlotMarkers -> "OpenMarkers", Frame -> True, FrameTicks -> {{All, None}, {Range[Length@energies], None}}, FrameStyle -> Directive[FontSize -> 14, FontFamily -> "Times", Black], Axes -> False, FrameLabel -> {"优化步骤", "能量(Hatree)"}, BaseStyle -> {FontFamily -> "Times", FontColor -> Black, FontSize -> 12}, GridLines -> {Range[Length@energies], Automatic}, GridLinesStyle -> Directive[Gray, Dashed], ImageSize -> Large}; energyPlot = ListLinePlot[ energies, Epilog -> Table[Inset[molPlotsThumbnails[[i]], {i, energies[[i]]}, {90, 90}, Scaled[0.2]], {i, 1, Length[energies]}], Sequence @@ plotOptions ]

2. 等值面静电势图
(* === 1. 定义文件路径 === *) path1 = "http://pic.wxyh.top/C3H4O2-esp.cube"; (*导入esp.cube文件的路径*) path2 = "http://pic.wxyh.top/C3H4O2-Potential.cube"; (*导入Potential.cube文件的路径*) (* === 2. 导入体积数据 === *) espCube = Import[path1, "VolumetricData"]; (*导入esp.cube文件*) potentialCube = Import[path2, "VolumetricData"]; (*导入Potential.cube文件*) (* === 3. 获取数据信息 === *) espData = espCube["Data"]; (*从esp.cube文件中提取数据*) potentialData = potentialCube["Data"]; (*从Potential.cube文件中提取数据*) dataRange = espCube["DataRange"]; (*获取数据的范围*) (* === 4. 处理数据 === *) espDataNormal = Normal[espData]; (*将espData数据标准化*) potentialDataNormal = Normal[potentialData]; (*将potentialData数据标准化*) (* === 5. 绘图 === *) mole = Import[path1, "Molecule"]; (*导入分子数据*) moleplot = MoleculePlot3D@mole; (*生成分子模型*) surface = ListContourPlot3D[espData[[1]], Contours -> 0.03 {-1, 1}, DataRange -> dataRange] // DiscretizeGraphics; (*生成等值面*) potentialPlot = ListSliceDensityPlot3D[potentialData[[1]], surface, DataRange -> dataRange, ColorFunction -> ({Opacity[0.8], Blend["M10DefaultDensityGradient", #1]} &), PlotLegends -> Automatic]; (*生成静电势图*) Show[moleplot, potentialPlot] (*显示分子模型和静电势图*)
