背景
在表面科学中,吸附能 $\Delta E_{\mathrm{ads}}$ 定义为吸附前后体系总能量的变化:
$$ \Delta E_{\mathrm{ads}} = E_{\mathrm{total}}^{\mathrm{ads}} - E_{\mathrm{total}}^{\mathrm{clean}} - E_{\mathrm{adatom}} $$
通过计算 adatom 在表面上方不同 $z$ 坐标处的能量,可以得到一条势能曲线,从而判断最有利于吸附的位置和强度。本脚本将多个元素、多个位点的数据统一绘制在同一张图中,便于横向比较。
代码结构概览
脚本 plot_energy_adatom_delta_E_ads_grid.py 整体分为以下几个部分:
- 配置与常量定义 —— 设定数据路径、图形参数、颜色映射等
- 数据读取与平滑 —— 从 CSV 文件读取原始数据并用 PCHIP 插值平滑
- 单面板绘制 —— 绘制一个元素的主图和放大图
- 网格布局组装 —— 将所有面板组合成 $3 \times 2$ 网格并添加图例
核心实现
数据读取
数据存储为 CSV 格式,每行包含 adatom 沿 $z$ 轴的坐标和对应的 $\Delta E_{\mathrm{ads}}$ 值。脚本通过 csv.DictReader 逐行解析:
def read_curve(path: Path) -> tuple[list[float], list[float]]:
with path.open(newline="") as file:
rows = list(csv.DictReader(file))
points = sorted((float(row[X_COL]), float(row[Y_COL])) for row in rows)
return [x for x, _ in points], [y for _, y in points]
数据按 $z$ 坐标升序排列,保证曲线从左到右连续。
PCHIP 平滑插值
原始 DFT 计算数据点数量有限,直接绘制会出现折线。脚本使用 SciPy 的 PchipInterpolator(保形分段三次 Hermite 插值)生成 $500$ 个插值点的平滑曲线:
from scipy.interpolate import PchipInterpolator
def smooth_curve(x_values, y_values):
unique_x, unique_indices = np.unique(x_array, return_index=True)
unique_y = y_array[unique_indices]
spline = PchipInterpolator(unique_x, unique_y)
x_smooth = np.linspace(unique_x.min(), unique_x.max(), SPLINE_POINTS)
y_smooth = spline(x_smooth)
return x_smooth, y_smooth
PCHIP 插值优于普通三次样条之处在于它不会在数据点之间产生过冲(overshoot),保持原始数据的单调性,适合物理量曲线的可视化。
双数据集叠加
脚本同时绘制两套数据:
| 数据集 | 目录 | 线型 | 标签后缀 |
|---|---|---|---|
| TA017(纯净表面) | dat/ta017_energy | 实线 (solid) | 无 |
| TA018(含空位表面) | dat/ta018_energy_vac | 虚线 (dashed) | (vac) |
通过循环遍历 DATASETS 元组,同一吸附位点的两条曲线使用相同颜色、不同线型叠加,直观对比空位对吸附能的影响:
DATASETS = (
(TA017_DIR, "solid", "{site}", "energy_{element}_{site}.csv"),
(TA018_DIR, "dashed", "{site}_vac", "energy_{element}_{site}_vac.csv"),
)网格布局与子图分割
整体图形使用 fig.add_gridspec(nrows=3, ncols=2) 创建 $3 \times 2$ 外层网格,每个外层单元格内部再用 subgridspec 将水平空间按 width_ratios=(1.0, ZOOM_WIDTH_RATIO) 分割为主图(左侧)和放大图(右侧):
inner_grid = outer_grid[row, col].subgridspec(
nrows=1,
ncols=2,
width_ratios=(1.0, ZOOM_WIDTH_RATIO),
wspace=ZOOM_PANEL_WSPACE,
)
main_ax = fig.add_subplot(inner_grid[0, 0])
zoom_ax = fig.add_subplot(inner_grid[0, 1])
放大图的比例由 ZOOM_X_WIDTH / MAIN_X_WIDTH 自动计算,确保放大区域的宽度与主图中矩形选框的宽度对齐。
局部放大区域
每个元素独立配置了 ZOOM_LIMITS,指定放大区域的 $x$ 和 $y$ 范围,聚焦于吸附能最低的过渡态附近区域。主图中用灰色虚线矩形框标出放大区域:
ZOOM_LIMITS = {
"N": {"x": (14.55, 15.75), "y": (-8.40, -2.40)},
"H": {"x": (14.60, 15.80), "y": (-3.20, -0.20)},
# ...
}刻度与标注
- $x$ 轴:使用
FixedLocator按整数间隔放置刻度 - $y$ 轴:使用
MaxNLocator(nbins=4)自动选择 4 个合适刻度,并用FormatStrFormatter("%.1f")保留一位小数 - 表面位置:以垂直虚线标出
SURFACE_Z = 14.059608(Å),标示固体表面的 $z$ 坐标 - 元素标签:每个主图的右上角标注元素符号
图例
图例位于图形底部中央,$8$ 列排列,每个吸附位点包含实线(无空位)和虚线(含空位)两条图例句柄:
fig.legend(
handles=build_legend_handles(),
loc="lower center",
bbox_to_anchor=(0.5, 0.03),
ncol=8,
# ...
)Matplotlib 全局配置
脚本通过 plt.rcParams.update() 统一设置全局样式:
plt.rcParams.update({
"font.family": "serif",
"mathtext.fontset": "cm",
"axes.unicode_minus": False,
"xtick.direction": "in",
"ytick.direction": "in",
"xtick.top": True,
"ytick.right": True,
"savefig.dpi": 400,
})
mathtext.fontset: "cm":使用 Computer Modern 字体渲染数学符号,与 LaTeX 风格一致axes.unicode_minus: False:避免负号显示为连字符- 四边刻度线全部向内,增强专业感
- 输出分辨率 $400$ DPI,保证打印质量
自定义坐标轴范围
脚本提供了 AXIS_LIMITS 字典,允许按元素手动设定 $x$ 和 $y$ 轴范围,设为 None 则回退到 Matplotlib 自动选择:
AXIS_LIMITS = {
"C": {"x": (13.8, 18.0), "y": (-9.0, 3.0)},
"H": {"x": (13.8, 18.0), "y": (-9.0, 3.0)},
# ...
}
这一设计使得六张子图的坐标轴范围保持一致,便于横向对比不同元素的吸附行为。
输出
最终图片保存为 viz/delta_E_ads_eV/energy_all.png,尺寸 $14.4 \times 10.2$ 英寸,$400$ DPI,适合直接用于学术论文插图。
小结
本脚本展示了使用 Matplotlib 制作复杂多面板科学图表的标准实践:
- subgridspec 实现面板内嵌子图
- PCHIP 插值 平滑离散数据
- 统一坐标范围 保证子图可比较性
- 双线型叠加 对比不同表面条件
- 手动刻度和标注 精确控制图表外观
该框架可灵活适配其他元素、吸附位点或表面类型的 DFT 计算数据。
References
- Virtanen, Pauli, et al. "SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python." Nature Methods, vol. 17, 2020, pp. 261–272. link
- Hunter, John D. "Matplotlib: A 2D Graphics Environment." Computing in Science & Engineering, vol. 9, no. 3, 2007, pp. 90–95. link
- Harris, Charles R., et al. "Array Programming with NumPy." Nature, vol. 585, 2020, pp. 357–362. link