COHP与PDOS合并图绘制指南：从数据对齐到出版级可视化

1. 项目概述：从“两张图”到“一张图”的化学键分析进阶

在计算化学和材料科学的日常研究中，我们经常遇到一个场景：你已经分别计算并绘制了晶体轨道哈密顿布居（COHP）图和投影态密度（PDOS）图。COHP图清晰地告诉你，在某个能量区间内，特定原子对之间的相互作用是成键（负值）还是反键（正值），量化了化学键的强度。PDOS图则展示了特定原子轨道（比如某个原子的d轨道）对体系总态密度的贡献，让你知道在哪些能量处这些轨道有电子占据。这两张图单独看都很有价值，但总感觉少了点什么——你无法直观地将“轨道成分”与“成键/反键特性”在同一个能量坐标下对应起来。这就是“COHP和PDOS合并图”要解决的痛点：它将COHP的成键信息与PDOS的轨道成分信息叠加在同一张图中，让你一眼就能看出，在某个成键或反键区域，主要是由哪些原子的哪些轨道在“贡献”。这个操作对于深入理解复杂体系中的化学键本质，比如催化活性中心的d带中心与吸附质成键的关系、电池材料中过渡金属的氧化还原机制等，至关重要。

画这样一张合并图，听起来只是简单的图形叠加，但实际操作中会遇到不少细节问题：能量零点如何对齐？坐标轴刻度、颜色、图例怎么设置才既美观又专业？不同软件输出的数据格式如何统一处理？本文将基于我处理大量VASP、Lobster等软件输出数据的经验，手把手带你完成从数据提取、处理到最终用Origin或Python（Matplotlib）绘制出高质量出版级合并图的全过程。无论你是刚接触能带分析的研究生，还是希望优化科研绘图流程的资深研究者，这篇指南都能提供直接的“抄作业”方案和避坑心得。

2. 核心思路与数据准备：对齐是成功的基石

在动笔画图之前，我们必须理清核心逻辑：COHP和PDOS数据之所以能合并，是因为它们共享同一个能量轴。这个能量通常是以费米能级（E_F）为参考零点。因此，合并图的第一步，也是最重要的一步，是确保两份数据的能量标度完全一致。

2.1 数据来源与格式解析

通常，COHP数据来自专门的分析软件，如 Lobster （与VASP联用）或 DPT-COHP （与Quantum ESPRESSO等联用）。PDOS数据则可以直接从主流量子化学软件（VASP、Quantum ESPRESSO、CASTEP）的输出文件中获取。

1. COHP数据文件 ：以Lobster输出为例，你会得到类似 COHPCAR.lobster 或 COOPCAR.lobster 的文件。用文本编辑器打开，其结构通常如下：

   # 能量(eV)     COHP值
   -10.0000     0.005432
   -9.9900      0.005567
   ...          ...
   10.0000      -0.002341
   

   第一列是相对于费米能级的能量，第二列是对应的COHP值（正值反键，负值成键）。有时文件会有多列，代表不同原子对或不同轨道的COHP，需要根据你的分析目标选择正确的列。

2. PDOS数据文件 ：以VASP的 DOSCAR 文件为例，在提取了总态密度后，你需要通过 vasppy 脚本或手动处理来获取特定原子的投影态密度。最终你需要的文件格式可能类似于：

   # 能量(eV)    s轨道    p轨道    d轨道    ... (或 px, py, pz等)
   -10.0000     0.000    0.001    0.050
   -9.9900      0.000    0.001    0.051
   ...          ...      ...      ...
   10.0000      0.002    0.005    0.020
   

   同样，第一列是能量，其他列是各轨道的态密度值。

关键检查点 ：务必确认两个文件中的能量范围（最小值、最大值）和能量步长是否一致。如果不一致，直接绘图会导致错位。通常，在计算时设置相同的能量窗口和采样点数可以避免此问题。

2.2 能量零点对齐与数据裁剪

即使能量范围一致，我们仍需确认能量零点是否都已校正到费米能级。在VASP中， DOSCAR 文件中的能量默认就是以费米能为零点的。而Lobster输出的COHP数据，通常也自动完成了这个校正。但为了保险起见，一个简单的验证方法是：观察PDOS在费米能级（0 eV）附近的分布是否符合金属或半导体的特征。

有时我们只关心费米能级附近某个特定能量窗口（例如 -5 eV 到 5 eV）的电子结构。这时就需要对原始数据进行裁剪。用Python的Pandas或简单的Shell命令（如 awk ）可以轻松完成：

# 示例：用awk提取能量在[-5, 5]区间内的COHP数据
awk '$1>=-5 && $1<=5 {print $1, $2}' COHPCAR.lobster > COHP_cropped.dat

对于PDOS数据也进行同样的裁剪操作，确保两份数据在合并前具有完全相同的能量点序列。这是后续所有操作的基础。

3. 绘图工具选型与基础图形绘制

合并图的绘制，主流有两种选择：专业的科学绘图软件（如Origin, SigmaPlot）和编程语言库（如Python的Matplotlib）。我将分别介绍，你可以根据熟悉程度选择。

3.1 方案一：使用OriginPro绘制（推荐新手/快速出图）

Origin以其强大的图形控制和友好的GUI著称，非常适合不常编程的研究者。

1. 导入数据 ：将处理好的 COHP_cropped.dat 和 PDOS_cropped.dat 分别导入Origin的两个工作表（Workbook）中。假设COHP数据有两列（能量，COHP值），PDOS数据有四列（能量，s, p, d）。

2. 绘制基础COHP图 ：

   • 选中COHP数据的工作表，点击菜单 Plot -> Basic -> Line 绘制COHP曲线。你会得到一张以能量为X轴，COHP值为Y轴的折线图。
   • 关键美化 ：双击坐标轴，将X轴标题设为“Energy (eV)”，Y轴标题设为“COHP (states/eV)”。由于COHP负值代表成键，为更直观，通常将成键区域填充颜色。右键点击曲线，选择“Plot Details”。在“Line”选项卡，将连接线设为较粗的实线。在“Fill Area Under Curve”选项卡，选择“Fill to bottom of the graph”，并为填充区域选择一个颜色（如浅蓝色代表成键）。

3. 叠加PDOS图 ：

   • 在当前的COHP图窗口，点击菜单 Graph -> Add Plot to Layer -> Line 。在弹出的对话框中选择PDOS数据工作表，并选择能量列作为X，d轨道列（或其他你关心的轨道）作为Y。点击“OK”，这样就在同一图层添加了一条PDOS曲线。
   • 新添加的曲线会与COHP共享X轴，但Y轴在左侧。我们需要为PDOS创建一个右侧Y轴。双击新添加的PDOS曲线，在“Plot Details”的“Line”选项卡设置其样式（如红色实线）。然后进入“Layer”管理，为这个图层添加一个右侧Y轴（ Add Right Y Axis ）。将这条PDOS曲线关联到右侧Y轴，并设置其标题为“PDOS (states/eV)”。

4. 合并与调整 ：

   • 现在你有了左侧Y轴的COHP（带填充）和右侧Y轴的PDOS曲线。调整两条曲线的颜色、线型，使其对比鲜明。
   • 核心技巧：填充方向 ：COHP的填充代表成键，我们希望它向下填充（负值区域）。而有时为了美观，我们希望成键的填充色在X轴下方，反键区域（如果需要显示）用另一种颜色填充在X轴上方。这需要在“Plot Details”中仔细设置填充的基线（Baseline），通常基线设为0。
   • 添加图例，清楚标明哪条线是COHP，哪条是PDOS（以及是什么轨道）。

3.2 方案二：使用Python Matplotlib绘制（推荐批量处理/高度定制）

对于需要批量生成图片或进行复杂自定义的情况，Python脚本是更强大的工具。

import numpy as np
import matplotlib.pyplot as plt

# 1. 加载数据
cohp_data = np.loadtxt('COHP_cropped.dat')
pdos_data = np.loadtxt('PDOS_cropped.dat')

energy_cohp = cohp_data[:, 0]
cohp = cohp_data[:, 1] # 假设第二列是总COHP或特定键的COHP

energy_pdos = pdos_data[:, 0]
pdos_d = pdos_data[:, 3] # 假设第四列是d轨道PDOS

# 2. 创建图形和左侧Y轴
fig, ax1 = plt.subplots(figsize=(8, 6))

# 绘制COHP曲线和填充
# 先画一条透明的线作为填充的边界
ax1.plot(energy_cohp, cohp, color='blue', linewidth=1.5, label='COHP')
# 填充成键区域（cohp<0的部分）
ax1.fill_between(energy_cohp, cohp, where=(cohp < 0), interpolate=True, color='skyblue', alpha=0.6, label='Bonding')
# 如果需要，也可以填充反键区域
# ax1.fill_between(energy_cohp, cohp, where=(cohp > 0), interpolate=True, color='lightcoral', alpha=0.6, label='Anti-bonding')

ax1.set_xlabel('Energy (eV)', fontsize=14)
ax1.set_ylabel('COHP (states/eV)', color='blue', fontsize=14)
ax1.tick_params(axis='y', labelcolor='blue')
ax1.axhline(y=0, color='black', linestyle='-', linewidth=0.8) # 画y=0的参考线
ax1.axvline(x=0, color='grey', linestyle='--', linewidth=1.0) # 画费米能级参考线

# 3. 创建右侧Y轴并绘制PDOS
ax2 = ax1.twinx()
ax2.plot(energy_pdos, pdos_d, color='red', linewidth=1.5, linestyle='-', label='d-PDOS')
ax2.set_ylabel('PDOS (states/eV)', color='red', fontsize=14)
ax2.tick_params(axis='y', labelcolor='red')

# 4. 合并图例（一个小技巧）
lines_1, labels_1 = ax1.get_legend_handles_labels()
lines_2, labels_2 = ax2.get_legend_handles_labels()
ax1.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper left', frameon=False)

# 5. 设置X轴范围，确保两者一致
ax1.set_xlim([-5, 5])
ax2.set_xlim([-5, 5])

# 6. 调整图形布局并保存
plt.title('COHP and d-PDOS for Metal-O Bond', fontsize=16)
fig.tight_layout()
plt.savefig('COHP_PDOS_merged.png', dpi=300, bbox_inches='tight')
plt.show()

这段代码提供了高度定制化的基础。你可以轻松修改颜色、线宽、填充透明度，添加多条PDOS曲线（如s, p, d分别用不同颜色和线型），或者将多个原子对的COHP绘制在一起进行比较。

4. 高级技巧与实战问题排查

掌握了基础绘制后，要产出真正能用于高水平论文的图表，还需要一些进阶技巧和问题处理经验。

4.1 多轨道PDOS与多键COHP的合并

现实分析中，我们常常需要同时展示多个轨道的PDOS，或者比较多个化学键的COHP。

• 多轨道PDOS ：在Python中，只需为 ax2 添加多条 plot 语句即可，并为每条线设置不同的颜色和标签。注意调整图例，避免过于拥挤。可以使用 linestyle （‘-’， ‘--’， ‘:’）来进一步区分。
• 多键COHP ：如果你想比较A-B键和A-C键的COHP，可以将它们绘制在同一个左侧Y轴上，用不同颜色和填充区分。此时，填充区域可能会重叠，需要通过设置透明度（ alpha ）来确保重叠部分可见。一个更清晰的方案是使用子图（ subplots ），将不同键的COHP-PDOS合并图并排展示。

4.2 能量对齐与插值处理

这是最容易出问题的地方。有时COHP和PDOS数据的能量网格点不完全相同（虽然范围一致，但点数或精确值有微小差异）。直接绘图会导致曲线不对齐或锯齿。

• 解决方案：插值 。使用 numpy.interp 或 scipy.interpolate 将其中一组数据插值到另一组数据的能量网格上。

  from scipy import interpolate
  # 假设我们希望将PDOS插值到COHP的能量网格上
  f_interp = interpolate.interp1d(energy_pdos, pdos_d, kind='linear', bounds_error=False, fill_value=0)
  pdos_d_interpolated = f_interp(energy_cohp)
  # 然后用 energy_cohp 和 pdos_d_interpolated 绘图
  

  bounds_error=False 和 fill_value=0 可以处理插值范围略微超出原数据范围的情况。

4.3 图形美化与出版标准

• 字体与尺寸 ：确保所有轴标签、刻度标签、图例的字体大小足够大（通常不小于10pt），以便在论文中缩小后仍清晰可读。推荐使用无衬线字体如Arial或Helvetica。
• 线宽与颜色 ：印刷时细线可能消失，确保线宽至少1.5pt。颜色选择要兼顾对比度和色盲友好性。避免同时使用红色和绿色作为主要区分色。
• 图例 ：图例应清晰指明每条曲线和每个填充区域的含义。如果图形元素太多，考虑直接在图上用箭头和文字标注关键特征，而不是全部塞进图例。
• 负值填充的视觉优化 ：在COHP图中，成键（负值）填充在X轴下方是惯例。在Matplotlib中， fill_between 的 where 参数和基线设置 y2=0 可以实现：

  ax1.fill_between(energy_cohp, cohp, 0, where=(cohp < 0), interpolate=True, color='skyblue', alpha=0.6)
  

5. 常见问题与排查技巧实录

在实际操作中，你肯定会遇到各种“坑”。下面是我总结的一些典型问题及解决方法。

个人实操心得 ：

• 数据备份 ：在进行任何裁剪、插值操作前，务必保留一份原始数据。所有处理步骤最好通过脚本完成，并记录下参数，确保可重复性。
• 分步验证 ：不要急于画合并图。先分别画出单独的COHP图和PDOS图，检查它们各自是否合理（如PDOS积分是否等于预期电子数，COHP成键反键区域是否符合化学直觉）。确认无误后再进行合并。
• 脚本化是王道 ：即使你更喜欢用Origin点击操作，也建议尝试用Python写一个基本的绘图脚本。当你需要为体系中的10个不同化学键生成10张类似的合并图时，脚本只需循环一下就能搞定，而手动点击会耗费大量时间且容易出错。
• 化学意义优先 ：图形再漂亮，也要服务于分析目的。始终问自己：这张合并图我想说明什么？是证明某个轨道主导了成键？还是比较不同键的强度？围绕核心结论来调整图形的呈现重点，比如突出某个能量区间，或高亮某条曲线。

绘制COHP与PDOS合并图，是将计算数据转化为深刻化学见解的关键一步。它不再是简单的画图任务，而是一个涉及数据校验、工具选择和美学表达的综合过程。希望这份详细的指南能帮你绕过我当年踩过的那些坑，高效地产出清晰、准确、有力的分析图表，让你在论文评审和学术报告中，能用更直观的方式讲述你发现的“电子故事”。