VASP DOS计算中eg与t2g轨道贡献的提取与分析求助

1. 解析DOSCAR和PROCAR文件的工具与操作示例 #

常用工具 #

• pymatgen（Python库）
  适合灵活的编程处理，安装命令：pip install pymatgen
  示例代码：

  from pymatgen.io.vasp.outputs import Doscar, Procar
  
  # 解析DOSCAR
  dos = Doscar("DOSCAR")
  energy_points = dos.energies  # 所有能量点
  total_dos = dos.total_dos.densities  # 总DOS密度（自旋向上/向下）
  element_dos = dos.get_element_dos()  # 按元素拆分的分波DOS
  
  # 解析PROCAR
  procar = Procar("PROCAR")
  # 提取第2个原子（索引从0开始）的所有轨道DOS
  atom_orbital_dos = procar.get_atom_dos(atom_index=1)
  

  操作说明：运行脚本后，可直接调用变量获取对应数据，后续可结合numpy、matplotlib做进一步处理。

• VASP自带脚本

  • dosplot.pl：VASP工具目录下的Perl脚本，直接在终端运行：perl dosplot.pl，会自动读取DOSCAR生成gnuplot绘图脚本，执行gnuplot dosplot.gp即可得到DOS图。
  • procar2dos.pl：用于从PROCAR提取分波DOS，命令：perl procar2dos.pl -f PROCAR，会生成每个原子、每个轨道的独立DOS数据文件，方便后续手动合并。

• 自定义Python脚本
  适合针对性提取数据，示例读取DOSCAR核心内容：

  import numpy as np
  
  # 跳过前6行系统信息，读取能量与DOS数据
  data = np.loadtxt("DOSCAR", skiprows=6)
  energy = data[:, 0]
  total_dos_up = data[:, 1]  # 自旋向上总DOS
  # 非自旋极化计算无下行，自旋极化时total_dos_down = data[:,2]
  

2. 分离eg与t2g轨道贡献 #

首先确认PROCAR中d轨道的顺序：通常d轨道按dxy、dyz、dxz、dx²−y²、dz²排列（可查看PROCAR表头的轨道标识确认），其中t2g包含前3个，eg包含后2个。

方法一：用pymatgen快速拆分 #

# 承接之前的procar对象，提取目标原子的d轨道DOS
atom_dos = procar.get_atom_dos(atom_index=1)
d_orbs = atom_dos["d"]  # 5个d轨道的DOS数组

# 合并t2g和eg轨道的DOS
t2g_dos = d_orbs[0] + d_orbs[1] + d_orbs[2]
eg_dos = d_orbs[3] + d_orbs[4]

方法二：基于procar2dos.pl的手动合并 #

1. 先用procar2dos.pl生成单个轨道的DOS文件，如atom_1_dxy.dat、atom_1_dyz.dat等；
2. 用awk命令合并t2g/eg轨道：

   # 生成t2g总DOS文件
   awk 'NR==FNR{a[$1]=$2;next}{a[$1]+=$2}END{for(k in a) print k, a[k]}' atom_1_dxy.dat atom_1_dyz.dat atom_1_dxz.dat > t2g_total.dat
   # 生成eg总DOS文件
   awk 'NR==FNR{a[$1]=$2;next}{a[$1]+=$2}END{for(k in a) print k, a[k]}' atom_1_dx2y2.dat atom_1_dz2.dat > eg_total.dat
   

绘图示例（matplotlib） #

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 6))
plt.plot(energy, total_dos_up, label="Total DOS", color="black", linewidth=1)
plt.plot(energy, t2g_dos, label="t2g DOS", color="blue", linestyle="--")
plt.plot(energy, eg_dos, label="eg DOS", color="red", linestyle=":")
plt.xlabel("Energy (eV)", fontsize=12)
plt.ylabel("DOS (states/eV)", fontsize=12)
plt.legend(fontsize=10)
plt.grid(alpha=0.3)
plt.show()

内容的提问来源于stack exchange，提问作者Farah Shehzadi