COHP与PDOS合并图绘制：从数据到化学键的深度分析

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

如果你在计算化学或者材料模拟领域摸爬滚打过一阵子，肯定对COHP（Crystal Orbital Hamilton Population，晶体轨道哈密顿布居）和PDOS（Projected Density of States，投影态密度）这两个分析工具不陌生。前者能定量告诉你原子间成键的强弱，后者则能揭示特定原子轨道对体系电子态的贡献。但很多时候，单独看其中一张图总觉得差点意思：PDOS图能看出轨道有贡献，但说不清这个贡献对成键是“建设性”的还是“破坏性”的；COHP图能告诉你成键强弱，但又看不清具体是哪些轨道在“干活”。于是，“把COHP和PDOS画到同一张图上”就成了一个非常自然且强烈的需求。这不仅仅是简单的图形叠加，而是一种将成键信息（COHP）与轨道成分信息（PDOS）进行空间对齐和关联解读的深度分析方法。

这张合并图的价值在于，它能让你一目了然地看到：在某个特定的能量区间（比如费米能级附近），哪些原子轨道有态密度分布（PDOS），同时这些轨道的相互作用对化学键的形成是起稳定作用（成键，负值COHP）还是 destabilizing 作用（反键，正值COHP）。这对于理解复杂体系中的电荷转移、轨道杂化、以及预测材料的稳定性、磁性、超导等性质至关重要。无论是研究二维材料的层间耦合，还是分析催化反应中活性位点的电子结构，亦或是探索拓扑材料中的奇异能带，这张合并图都是你从“看到现象”走向“理解本质”的关键一步。

接下来，我将以一个典型的案例——比如分析石墨烯中某个C-C键——为线索，拆解从数据准备、工具选型、到绘图美化、再到专业解读的全流程。你会发现，虽然核心步骤可能只需要几行代码，但背后的参数调整、细节处理和解读逻辑，才是真正体现功力的地方。

2. 核心思路与工具链选型：为什么是它们？

在动手画图之前，我们需要先搭建一条高效、可靠的数据处理与可视化流水线。这条流水线的起点是你的第一性原理计算软件（如VASP、Quantum ESPRESSO），终点是一张信息丰富、出版质量的合并图。中间环节的工具选型，直接决定了过程的顺畅度和结果的可信度。

2.1 数据处理核心：LOBSTER 程序

对于COHP分析， LOBSTER 程序是目前事实上的行业标准。它不是一个独立的计算软件，而是一个后处理工具。你需要先通过VASP等软件完成一个精确的电子结构自洽计算，并获得 WAVECAR 和 LOCPOT 等文件。LOBSTER 的作用，就是利用这些文件，并结合你提供的结构信息，通过构造局域原子轨道基组，将平面波基组下的波函数进行投影，从而计算出COHP、COOP（晶体轨道重叠布居）和PDOS等化学键分析数据。

选择LOBSTER的理由很充分：

1. 精度与可靠性 ：它经过了大量文献的验证，其结果与化学直觉和实验观测（如键长、键能）符合得很好。
2. 功能全面 ：除了COHP，还能计算COOP、原子和轨道的布居数（Mulliken / Löwdin），以及 spilling parameter（用于评估投影质量）。
3. 社区与文档 ：拥有相对活跃的用户社区和逐步完善的文档，遇到问题有据可查。

注意 ：运行LOBSTER前，务必确保初始的VASP计算是高质量的。特别是 INCAR 中需要设置 LORBIT = 11 来生成 PROCAR 文件（虽然LOBSTER主要用 WAVECAR ，但某些设置与此相关），并且 KPOINTS 网格要足够密集，以保证能量本征值和波函数的准确性。一个粗糙的 KPOINTS 网格会导致COHP/PDOS图像出现不真实的锯齿和毛刺。

2.2 可视化方案：Python + Matplotlib/Seaborn

得到LOBSTER输出的文本数据（通常是 COHPCAR.lobster 和 DOSCAR.lobster 或 projectedData 目录下的文件）后，我们需要一个灵活强大的工具来绘图。这里首推 Python 生态，特别是 Matplotlib 库，并可以搭配 Seaborn 美化样式。

为什么不直接用Origin或Grace等图形软件？原因在于 批量处理 和 定制化 。一个研究项目往往涉及多个键位、不同掺杂浓度、不同应变条件的对比。用Python写一个脚本，可以轻松实现循环读取、自动对齐费米能级、统一绘图风格、并批量生成数十张高质量的图片，效率是手动操作无法比拟的。此外，对于合并图中复杂的图例、双Y轴、填充颜色、线型样式等细节，Matplotlib提供了像素级的控制能力。

基础工具链 ：

• Python 3.8+ ： 编程环境。
• NumPy ： 用于高效读取和处理文本数据。
• Matplotlib ： 核心绘图库。我们主要使用其 pyplot 接口。
• Seaborn (可选)： 提供更美观的绘图样式（如 sns.set_style(“whitegrid”) ），让科研图表瞬间提升一个档次。
• Pandas (可选)： 如果数据整理复杂，可以用它来管理DataFrame。

3. 数据准备与核心参数解析

拿到LOBSTER的输出文件后，别急着画图。花几分钟理解文件结构和关键参数，能避免很多低级错误。

3.1 理解LOBSTER输出文件

一次典型的LOBSTER成功运行后，你会得到一系列文件，我们需要重点关注这几个：

1. COHPCAR.lobster ： 这是核心文件。它通常包含多列数据：

   • 第一列：能量（相对于费米能级，单位通常是eV）。
   • 第二列：总COHP（所有选定原子对和轨道的总和）。
   • 第三列及以后：特定原子对、特定轨道组合的COHP。例如，“C1-2s_C2-2pz”表示原子C1的2s轨道与原子C2的2pz轨道的COHP。
   • 文件头部会有注释行，说明每一列对应的具体轨道。 务必仔细阅读头部信息 ，这是你正确选择绘图数据列的唯一依据。

2. DOSCAR.lobster 或 projectedData/DOS 开头的文件： 这些是投影态密度数据。 DOSCAR.lobster 可能包含总DOS和分轨道PDOS。更常见的是，在 projectedData 文件夹下，会有 DOS_{AtomIndex}.{Orbital}.dat 这样的文件，例如 DOS_1.s.dat 代表第1个原子的s轨道PDOS。这种分文件存储的结构更清晰，便于程序化读取。

3. ICOHPLIST.lobster ： 这个文件列出了所有分析的原子对及其积分COHP（ICOHP，即COHP曲线在费米能级以下的面积）。ICOHP的负值越大，表示成键越强。在画图前，可以先查看这个文件，筛选出成键作用最强的几对原子进行重点分析。

3.2 关键参数与预处理

1. 费米能级对齐 ： LOBSTER输出的能量零点默认是计算设置的费米能级。但为了与实验对比（如ARPES）或不同计算之间对比，我们通常需要将费米能级对齐到0 eV处。检查你的 COHPCAR.lobster 和 DOSCAR.lobster 第一列数据的范围，如果费米能级不在0点，需要在读取数据后，将所有能量值减去费米能级 E_fermi 。这个 E_fermi 通常可以在VASP的 OUTCAR 或 vasprun.xml 中找到。

2. 能量范围选取 ： 绘图时不需要展示全能量区间的数据。通常关注费米能级附近[-10, 5] eV的范围就足够了，这个区间包含了价带顶和导带底的关键信息。设置合理的 xlim 可以让图像重点突出。

3. 轨道选择与归类 ： 这是合并图信息组织的关键。你不需要把所有轨道都画上去，那会导致图面杂乱。例如，分析过渡金属氧化物的金属-氧键时，你可能会将过渡金属的d轨道（ dxy, dyz, dz2, dxz, dx2-y2 ）的PDOS求和，画成一条“总d-PDOS”曲线，同时将金属d轨道与氧p轨道的COHP也画出来。这样就能清晰看到d-p杂化对成键的贡献。在LOBSTER中，可能需要手动对多个数据列进行求和来得到归类的PDOS或COHP。

4. 实操绘图：一步步生成合并图

假设我们已经通过LOBSTER分析了一对C原子（索引为1和2）之间的键，并且得到了分轨道的PDOS数据。现在，我们用Python脚本将C1原子的 2pz 轨道PDOS与C1-C2原子间的 pz-pz 轨道COHP画在同一张图上。

4.1 环境准备与数据读取

首先，确保安装了必要的库，并读取数据。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 设置Seaborn样式，让图表更美观
sns.set_style("whitegrid")
sns.set_context("notebook", font_scale=1.2)

# 1. 读取COHP数据 (假设COHPCAR.lobster中，第4列是C1-pz和C2-pz的COHP)
# 跳过文件头部的注释行（具体行数需要根据你的文件查看）
cohp_data = np.loadtxt('COHPCAR.lobster', skiprows=5) # skiprows需要调整
energy_cohp = cohp_data[:, 0] - e_fermi  # 第一列是能量，对齐费米能级
cohp_pz_pz = cohp_data[:, 3]  # 假设第4列是目标COHP数据，单位是eV

# 2. 读取PDOS数据 (假设从projectedData/DOS_1.pz.dat读取)
pdos_data = np.loadtxt('projectedData/DOS_1.pz.dat')
energy_pdos = pdos_data[:, 0] - e_fermi  # 同样对齐费米能级
pdos_pz = pdos_data[:, 1]  # 假设第二列是态密度值，单位是states/eV

# 检查能量网格是否一致。通常LOBSTER输出时，COHP和PDOS的能量网格是相同的。
# 如果不一致，可能需要插值。这里假设它们一致。
assert np.allclose(energy_cohp, energy_pdos), “能量网格不一致，需要插值处理！”

4.2 构建双Y轴合并图

这是最关键的一步。我们将使用Matplotlib的 twinx() 方法来创建共享X轴（能量轴）的双Y轴图。

# 创建图形和第一个Y轴（用于PDOS）
fig, ax1 = plt.subplots(figsize=(8, 6))

# 在ax1上绘制PDOS (填充曲线，更直观显示轨道占据)
# 选择一种颜色，比如蓝色，代表轨道贡献
color_pdos = ‘tab:blue’
ax1.fill_between(energy_pdos, 0, pdos_pz, where=(pdos_pz>0),
                  facecolor=color_pdos, alpha=0.3, label=‘C1-2pz PDOS’)
ax1.plot(energy_pdos, pdos_pz, color=color_pdos, linewidth=2)
ax1.set_xlabel(‘Energy (eV)’, fontsize=14)
ax1.set_ylabel(‘Projected DOS (states/eV)’, color=color_pdos, fontsize=14)
ax1.tick_params(axis=‘y’, labelcolor=color_pdos)
ax1.set_xlim([-10, 5])  # 设置合理的能量范围
ax1.axvline(x=0, color=‘k’, linestyle=‘--’, linewidth=1, alpha=0.5) # 画出费米能级线

# 创建第二个Y轴，与ax1共享X轴，用于绘制COHP
ax2 = ax1.twinx()
# 选择另一种颜色，比如红色，代表成键信息
color_cohp = ‘tab:red’
# COHP的绘图惯例：成键区域（负值）通常在X轴上方，反键区域（正值）在下方。
# 我们可以通过取反cohp_pz_pz来实现，这样负的COHP（成键）在图上就是正的填充。
cohp_to_plot = -cohp_pz_pz  # 注意：这里取反了！
ax2.fill_between(energy_cohp, 0, cohp_to_plot, where=(cohp_to_plot>0),
                  facecolor=color_cohp, alpha=0.3, label=‘-COHP(C1-pz, C2-pz)’)
ax2.plot(energy_cohp, cohp_to_plot, color=color_cohp, linewidth=2)
ax2.set_ylabel(‘-COHP (eV)’, color=color_cohp, fontsize=14)
ax2.tick_params(axis=‘y’, labelcolor=color_cohp)
# 设置COHP Y轴的范围，可能需要根据数据调整
# ax2.set_ylim([-5, 10])

# 合并图例是一个小技巧。分别从两个轴上获取句柄和标签
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=True)

# 添加标题
plt.title(‘C-C Bond Analysis: PDOS vs. COHP’, fontsize=16, pad=20)

# 调整布局，防止标签重叠
fig.tight_layout()

# 保存图片，设置高DPI以保证出版质量
plt.savefig(‘C_C_bond_pdos_cohp.png’, dpi=300, bbox_inches=‘tight’)
plt.show()

这段代码生成了核心的合并图。左侧Y轴（蓝色）表示C1原子2pz轨道的态密度，右侧Y轴（红色）表示C1与C2原子pz轨道间相互作用的负COHP。填充区域使正负值更加直观。

4.3 图表美化与出版级调整

上面的代码生成了可用的图，但要达到投稿或报告的水平，还需要精细调整。

1. 线型与颜色 ：

   • 使用颜色盲友好的配色方案（如 tab:blue , tab:orange , tab:green ）。避免同时使用红色和绿色，以免色盲读者无法区分。
   • 对于多条PDOS曲线（如s, p_x, p_y, p_z），可以使用不同线型（实线、虚线、点划线）配合颜色来区分。

2. 坐标轴与刻度 ：

   • 统一刻度线的朝向和宽度： ax1.tick_params(axis=‘both’, which=‘major’, length=6, width=1.5)
   • 设置科学计数法或调整刻度密度，避免数字重叠。

3. 图例 ：

   • 将图例放在图外（ bbox_to_anchor ）可以节省图内空间，使画面更整洁。

   ax1.legend(lines_1 + lines_2, labels_1 + labels_2,
              loc=‘upper left’, bbox_to_anchor=(1.05, 1),
              borderaxespad=0., frameon=True)
   

4. 网格线 ：

   • 保留 sns.set_style(“whitegrid”) 提供的浅灰色网格线，有助于读者追踪数据点。

5. 字体 ：

   • 将字体改为常用的衬线字体（如Times New Roman）或无衬线字体（如Arial），以符合出版要求。

   plt.rcParams[‘font.family’] = ‘sans-serif’
   plt.rcParams[‘font.sans-serif’] = [‘Arial’]
   

5. 深度解读：从图中读出化学故事

图画好了，但真正的功夫在于解读。面对合并图，我们应该按什么逻辑去分析？

5.1 分区域解读法

以我们生成的C-C键 pz-pz 相互作用图为例：

1. 深能级区域（-10 eV 到 ~-5 eV） ：

   • 观察 ：PDOS曲线（蓝色）可能有峰值，-COHP曲线（红色）在此区域也有显著的 正 填充（记住我们画的是-COHP）。
   • 解读 ：PDOS峰值表示C原子的2pz轨道在这个能量区间有电子态分布。同时，正的-COHP（即负的COHP）意味着此处的轨道相互作用是 成键性 的。这说明在深能级处，两个C原子的pz轨道形成了强成键态（σ键或强杂化键），这些态被电子完全占据，是化学键的主要贡献者之一。

2. 费米能级附近区域（-5 eV 到 2 eV） ：

   • 观察 ：PDOS在费米能级处可能有一个有限的值（对于金属或半金属），而-COHP曲线穿过零点，从正变负。
   • 解读 ：在费米能级处有PDOS，说明该轨道对费米面的电子态有贡献。COHP从正（成键）变为负（反键），意味着在费米能级附近，两个pz轨道的相互作用从成键态过渡到了反键态。对于石墨烯，这对应着π键的成键态（价带）和反键态（导带），狄拉克点就在它们之间。如果费米能级穿过的是反键态区域，说明这个键可能因为电子占据反键态而被削弱。

3. 高能反键区域（2 eV 以上） ：

   • 观察 ：-COHP显示为 负 填充区域（即原始的COHP为正值）。
   • 解读 ：这是明确的反键态区域。如果这些反键态被电子占据（比如通过掺杂），将会严重削弱原子间的键强，可能导致结构失稳。

5.2 积分COHP（ICOHP）的辅助判断

光看曲线形状还不够定量。这时需要回头查看 ICOHPLIST.lobster 文件。找到对应键的ICOHP值（例如 ICOHP: -5.32 eV ）。这个负值表示从能量最低态积分到费米能级，总的成键相互作用强度为5.32 eV。绝对值越大，键越强。你可以对比不同键的ICOHP，快速判断键的相对强弱。

5.3 多轨道合并分析实战

实际分析中，我们经常需要合并多个轨道。例如，分析过渡金属（TM）与配体（L）的键：

• 步骤一 ：分别计算TM的 d 轨道总PDOS（将 dxy, dyz, dz2, dxz, dx2-y2 五个文件的PDOS相加）和配体L的 p 轨道总PDOS。
• 步骤二 ：计算TM-d与L-p轨道之间的总COHP（LOBSTER可以输出分轨道的COHP，你需要把TM-d和L-p所有可能的轨道组合对的COHP相加，或者直接使用LOBSTER输出的“总COHP”，但需明确它包含了所有轨道）。
• 步骤三 ：将TM-d的PDOS和TM-d与L-p的总COHP画在同一张合并图上。
• 解读 ：你可以清晰地看到，在哪些能量区间，TM-d轨道有电子态（PDOS峰值），并且这些态与配体p轨道的相互作用是强成键（负COHP）还是强反键（正COHP）。这直接揭示了配体场的强弱和共价键的程度。

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

在实际操作中，你一定会遇到各种问题。下面是我踩过的一些坑和解决方案。

6.1 LOBSTER计算相关

问题1：Spilling Parameter 过高（>10%）

• 现象 ：LOBSTER运行日志中提示spilling parameter过大。
• 原因 ：这意味着你选择的局域原子轨道基组无法很好地展开平面波波函数，投影质量差，结果不可信。
• 解决 ：
  1. 检查并尝试使用LOBSTER内置的不同基组（如 pbeVaspFit2015 或 pbeVaspFit2012 ）。
  2. 确保VASP计算时 ENCUT 值足够高（通常比默认值高1.3倍）。
  3. 对于包含重元素的体系，考虑在 lobsterin 文件中使用 basisfunctions 关键词手动调整基组，例如为过渡金属添加更灵活的极化函数。

问题2：COHP曲线数值异常小或形状奇怪

• 现象 ：COHP曲线几乎是一条接近零的直线，或者在不该有信号的能量区间出现剧烈震荡。
• 原因 ：
  • KPOINTS 网格太稀疏。COHP对k点采样非常敏感。
  • 在 INCAR 中未设置 LORBIT=11 或 LWAVE=.TRUE. 。
  • 在 lobsterin 文件中指定的原子对或轨道不存在或索引错误。
• 解决 ：
  1. 首要检查 ：使用更密的k点网格重新进行VASP计算和LOBSTER分析。对于能带结构复杂的体系， KPOINTS 可能需要达到 15x15x15 甚至更高。
  2. 确认VASP计算正确输出了 WAVECAR 和 LOCPOT 文件。
  3. 仔细核对 lobsterin 文件中的 cohpbetween 和 basisfunctions 设置。

6.2 数据处理与绘图相关

问题3：PDOS和COHP的能量范围或点数不一致

• 现象 ：绘图时报错，或曲线错位。
• 解决 ：使用 numpy.interp 函数将其中一组数据插值到另一组数据的能量网格上。

  # 将cohp数据插值到与pdos相同的能量网格上
  energy_new = energy_pdos # 目标网格
  cohp_interp = np.interp(energy_new, energy_cohp, cohp_pz_pz)
  # 然后用 cohp_interp 和 pdos_pz 进行绘图
  

  注意 ：插值会引入微小误差，尤其是原始数据点很稀疏时。最佳实践是确保LOBSTER计算时使用相同的能量网格参数输出COHP和PDOS。

问题4：合并图图例重叠或显示不全

• 现象 ：双Y轴的图例挤在一起，或第二个Y轴的标签被遮挡。
• 解决 ：
  1. 使用 fig.tight_layout() 或 plt.subplots_adjust() 调整整体布局，给图例留出空间。
  2. 将图例放置在图外，如 ax1.legend(..., bbox_to_anchor=(1.05, 1), loc=‘upper left’) 。
  3. 调整图形尺寸 figsize=(10, 6) 以获得更宽的画布。

问题5：如何绘制“分轨道”贡献的COHP？

• 需求 ：不仅想看总COHP，还想知道是哪些具体的轨道对（如 dxy-px ）对成键贡献最大。
• 方法 ：LOBSTER的 COHPCAR.lobster 文件通常包含了分轨道的数据列。你需要根据文件头部说明，找到对应列。在绘图时，可以将这些分轨道的COHP分别用细线、半透明的方式画在图上，再用一条粗线表示总COHP。

  # 假设第4，5，6列是三个分轨道COHP
  cohp_orb1 = cohp_data[:, 3]
  cohp_orb2 = cohp_data[:, 4]
  cohp_orb3 = cohp_data[:, 5]
  cohp_total = cohp_data[:, 1] # 第二列通常是总COHP
  
  ax2.plot(energy, -cohp_orb1, color=‘red’, alpha=0.4, linewidth=0.8, label=‘orb-pair1’)
  ax2.plot(energy, -cohp_orb2, color=‘green’, alpha=0.4, linewidth=0.8, label=‘orb-pair2’)
  ax2.plot(energy, -cohp_total, color=‘black’, linewidth=2.5, label=‘Total -COHP’)
  

  这样，你就能在总趋势下，清晰地辨别出各个轨道相互作用的细节。

画出一张漂亮的COHP和PDOS合并图只是第一步，更重要的是养成一套系统的分析习惯。我的体会是，永远不要孤立地看待一张图。把它和能带结构、差分电荷密度、声子谱等其它分析手段的结果相互印证。例如，合并图上显示费米能级处有一个反键态峰被部分占据，那么你应该能在能带结构上找到对应的穿过费米能的平带，在差分电荷密度图上看到相关原子间电子密度的减少。这种多角度证据链的闭合，才是计算模拟工作说服力的真正来源。最后一个小技巧，建立一个自己的绘图脚本模板，将数据读取、费米能级对齐、颜色方案、字体设置、保存格式等通用部分封装好。下次分析新体系时，你只需要修改文件路径和轨道索引，就能快速生成风格统一、可直接用于组会或论文的图片，这将极大提升你的研究效率。