1. 可视化优化步骤

(* 定义文件路径 *)
filepath = "http://pic.wxyh.top/c3.log";

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] (*显示分子模型和静电势图*)