8 辅助工具使用教程

备注

用DSPAW完成DFT计算后，想快速分析结果、画图，或者完成一些常见的数据处理任务？

dspawpy （ Python ≥ 3.8 ）就是这样一款辅助工具，它可以编程调用（参考下文示例脚本），也提供了一个命令行交互程序

参考下方教程安装后，命令行输入 dspawpy 后回车即可使用交互程序：

 ... loading dspawpy cli ...

 ********这是dspawpy命令行交互小工具，预祝您使用愉快********
    ( )
   _| |  ___  _ _      _ _  _   _   _  _ _    _   _
 /'_` |/',__)( '_`\  /'_` )( ) ( ) ( )( '_`\ ( ) ( )
( (_| |\__, \| (_) )( (_| || \_/ \_/ || (_) )| (_) |
 \__,_)(____/| ,__/'`\__,_) \___x___/ | ,__/  \__, |
             | |                      | |    ( )_| |
             (_)                      (_)     \___/

 [版本号]: [安装路径]

 ======================================
 |  1: update更新
 |  2: structure结构转化
 |  3: volumetricData数据处理
 |  4: band能带数据处理
 |  5: dos态密度数据处理
 |  6: bandDos能带和态密度共同显示
 |  7: optical光学性质数据处理
 |  8: neb过渡态计算数据处理
 |  9: phonon声子计算数据处理
 |  10: aimd分子动力学模拟数据处理
 |  11: Polarization铁电极化数据处理
 |  12: ZPE零点振动能数据处理
 |  13: TS的热校正能
 |
 |  q: 退出
 ======================================
 --> 输入数字后回车选择功能：

亮点：

自动补全：按Tab键即可生效，有助于快速且正确地输入程序所需要的参数
多线程懒加载：在等待用户输入时后台加载模块，大幅缩短等待时间；仅加载必要模块，减少内存占用

注意：

在远程服务器上使用时，由于机器的磁盘读写性能不佳，可能需要较长的启动时间，极端情况下可能长达半分钟（与服务器当时的磁盘读写性能直接相关）。如果无法忍受，请在自己的电脑上安装dspawpy使用
输入 cli 回车后，python会先加载内建模块，完成后显示 … loading dspawpy cli … 提示语句，表明进入第二阶段（开始加载第三方依赖库）
等待第二阶段完成后，会显示欢迎界面，此时dspwapy已完成前期加载，进入第三阶段。后续将动态根据选择的功能模块加载相应的依赖库，从而最大限度减少等待时间

8.1 安装与更新

在HZW机器上，已提前安装 dspawpy，使用以下命令激活虚拟环境即可使用：

source /data/hzwtech/profile/dspawpy.env

在其他机器上，请自行安装 dspawpy（下列两种方式二选一）：

使用 mamba 或 conda ，安装包可前往 https://conda-forge.org/download/ 下载

mamba install dspawpy -c conda-forge
#conda install dspawpy -c conda-forge

或使用 pip3 (部分操作系统中没有pip3可执行程序，此时请尝试pip)

pip3 install dspawpy

关于如何设置pip和conda镜像地址以加速安装过程，可参考 https://mirrors.tuna.tsinghua.edu.cn/help/pypi/ 和 https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/

如果安装依然失败，请尝试上面的 mamba/conda 安装法。

警告

在集群上，由于权限问题，公共路径中的pip很可能不支持全局安装python库，必须在 pip install 后面追加 --user 选项安装到家目录 ~/.local/lib/python3.x/site-packages/ 下，其中3.x表示python解释器版本，x可能是8-12之间的任意整数

python将优先读取家目录下的dspawpy，即使公共环境中的dspawpy的版本比家目录中的dspawpy的版本新! 因此，如果以前用 --user 安装过dspawpy，又忘记手动更新家目录下的老版本dspawpy，即使source了公共环境，也无法调用公共环境中的dspawpy，依旧会使用老版本dspawpy，导致一些BUG。

因此，鉴于HZW集群每周自动更新dspawpy，建议不要在家目录下重复安装，已安装的建议删除 。如果在其他集群上，请确保及时手动更新家目录下的 dspawpy

如果不想删除家目录中的dspawpy，又不想更新它，可以在运行python脚本时，尝试 -s 选项避免导入家目录中的 dspawpy： python -s your-script.py。cli 不支持这个命令行选项

更新 dspawpy

如果 dspawpy 是通过 mamba/conda 安装的，使用以下命令更新：

mamba update dspawpy
#conda update dspawpy

如果 dspawpy 是通过 pip 安装的，那么：

pip install dspawpy -U # -U 表示安装最新版

备注

如果 pip 使用了国内的镜像网站，可能由于镜像网站尚未同步最新版 dspawpy 而导致无法顺利升级，请使用如下命令告诉 pip 从pypi官网下载安装：

pip install dspawpy -i https://pypi.org/simple --user -U # -i 指定下载地址，--user 表示仅为当前用户安装，-U 表示安装最新版

如果运行时出现 dspawpy 相关 错误信息，请先检查是否已正确导入最新版本的 dspawpy 并检查安装路径：

$ python3 # 或 python
>>> import dspawpy
>>> dspawpy.__version__ # 将输出版本号
>>> dspawpy.__file__ # 将输出安装路径

8.2 structure结构转化

读取结构信息使用 read 函数；将结构信息写入文件，使用 write 函数；快速结构转化，使用 convert 函数：

可参考 2conversion.py 脚本完成转化：

from dspawpy.io.structure import convert

convert(
    infile="dspawpy_proj/dspawpy_tests/inputs/2.1/relax.h5",  # 待转化结构，如果在当前路径，可以只写文件名
    si=None,  # 筛选构型编号，如果不指定，默认读取全部
    ele=None,  # 筛选元素符号，默认读取所有元素的原子信息
    ai=None,  # 筛选原子编号，从 1 开始，默认读取所有原子信息
    infmt=None,  # 输入结构文件类型，例如 'h5'，如果为None，将根据文件名规则模糊匹配
    task="relax",  # 任务类型，此参数仅当 infile 为文件夹而非文件名时有效
    outfile="dspawpy_proj/dspawpy_tests/outputs/us/relaxed.xyz",  # 结构文件名
    outfmt=None,  # 输出结构文件类型，例如 'xyz'，如果为None，将根据文件名规则模糊匹配
    coords_are_cartesian=True,  # 默认以笛卡尔坐标写入
)

convert 函数的几个关键参数设置规则见下表：

**dspawpy 支持读写的结构文件列表**
infmt（输入文件格式）	infile（输入文件名模糊匹配）	outfmt（输出文件格式）	outfile（输出文件名模糊匹配）	说明
h5	*.h5	X	X	DS-PAW计算完成后保存的hdf5类型文件
json	*.json	json	*.json	DS-PAW计算完成后保存的json类型文件
pdb	*.pdb	pdb	*.pdb	Protein Data Bank
as	*.as	as	*.as	DS-PAW记载原子坐标等信息的结构文件
hzw	*.hzw	hzw	*.hzw	DeviceStudio默认的结构文件
xyz	*.xyz	xyz	*.xyz	读取时仅支持分子结构的单构型，写入时则是包含晶胞的extended-xyz类型轨迹文件
X	X	dump	*.dump	lammps-dump类型轨迹文件
X	.cif/.mcif	cif/mcif	.cif/.mcif	Crystallographic Information File
X	POSCAR/CONTCAR/.vasp/CHGCAR/LOCPOT/vasprun.xml*	poscar	POSCAR	VASP文件
X	.cssr	cssr	.cssr	Crystal Structure Standard Representation
X	.yaml/.yml	yaml/yml	.yaml/.yml	YAML Ain’t Markup Language
X	.xsf	xsf	.xsf	eXtended Structural Format
X	rndstr.in/lat.in/bestsqs	mcsqs	rndstr.in/lat.in/bestsqs	Monte Carlo Special Quasirandom Structure
X	inp.xml/.in/inp_	fleur-inpgen	.in	FLEUR 结构文件，须额外安装 pymatgen-io-fleur 库
X	*.res	res	*.res	ShelX res 结构文件
X	.config/.pwmat	pwmat	.config/.pwmat	PWmat 文件
X	X	prismatic	prismatic	用于STEM模拟的一种文件格式
X	CTRL*	X	X	Stuttgart LMTO-ASA 文件

备注

上述表格中 * 号表示任意字符，X 表示不支持
h5, json, pdb, xyz, dump, CONTCAR等格式支持多个结构组成的轨迹信息（常见于结构优化或者NEB或者AIMD任务）
in(out)fmt 参数优先级高于文件名模糊匹配规则；例如，指定 in(out)fmt=’h5’ 后，文件名可以是任意值，甚至可以是 a.json
将结构信息写为 json 格式时，仅支持可视化 NEB 链任务，详见观察NEB链节相关说明
DS-PAW 输出的 neb.h5、phonon.h5、phonon.json、neb.json以及wannier.json 暂时无法读取结构信息

8.3 volumetricData数据处理

volumetricData可视化

参考 3vis_vol.py ：

from dspawpy.io.write import write_VESTA

# 读取数据文件（h5或json格式），处理后输出到cube文件
write_VESTA(
    in_filename="dspawpy_proj/dspawpy_tests/inputs/2.2/rho.h5",  # 包含电子体系信息的json或h5文件路径
    data_type="rho",  # 数据类型，支持 "rho", "potential", "elf", "pcharge", "rhoBound"
    out_filename="dspawpy_proj/dspawpy_tests/outputs/us/DS-PAW_rho.cube",  # 输出文件路径
    gridsize=(10, 10, 10),  # 指定插值网格大小
    format="cube",  # 支持cube, vesta和txt格式（xyz格点坐标+数值）
)

将转换后的文件 DS-PAW_rho.cube 拖入 VESTA 软件中显示：

晶体硅单胞的电荷密度分布示意图

差分volumetricData可视化

参考 3dvol.py ：

from dspawpy.io.write import write_delta_rho_vesta

# 读取数据文件（h5或json格式），处理后输出到cube文件，可直接使用vesta打开且体积小
write_delta_rho_vesta(
    total="dspawpy_proj/dspawpy_tests/inputs/supplement/AB.h5",  # 包含所有组元的体系的数据文件
    individuals=[
        "dspawpy_proj/dspawpy_tests/inputs/supplement/A.h5",
        "dspawpy_proj/dspawpy_tests/inputs/supplement/B.h5",
    ],  # 包含各组元的体系的数据文件
    output="dspawpy_proj/dspawpy_tests/outputs/us/3delta_rho.cube",  # 输出文件路径
)

上述脚本支持处理多元体系的电荷密度差，示例以二元体系为例，得到了 AB.h5 与 A.h5 和 B.h5 之间的电荷密度差文件 delta_rho.cube ，该文件可直接使用 VESTA 打开。

volumetricData面平均

参考 3planar_ave.py ：

from dspawpy.plot import average_along_axis

axes = ["2"]  # “0”, "1", "2" 分别对应xyz轴向，选择要沿着哪些轴平均
axes_indices = [int(i) for i in axes]
for ai in axes_indices:
    plt = average_along_axis(
        datafile="dspawpy_proj/dspawpy_tests/inputs/3.3/scf.h5",  # 数据文件路径
        task="potential",  # 任务名称，可以是 'rho', 'potential', 'elf', 'pcharge', 'rhoBound'
        axis=ai,  # 沿着哪个轴向绘制势能曲线
        smooth=False,  # 是否平滑
        smooth_frac=0.8,  # 平滑系数
        subtype=None,  # 用于指定task数据子类，目前仅用于 Potential
        label=f"axis{ai}",  # 图例标签
    )
if len(axes_indices) > 1:
    plt.legend()

plt.xlabel("Grid Index")
plt.ylabel("TotalElectrostaticPotential (eV)")
plt.savefig("dspawpy_proj/dspawpy_tests/outputs/us/3pot_ave.png", dpi=300)  # 图片名称

处理应用案例 3.3 小节所得静电势文件，可得真空方向势函数曲线如下所示：

AuAl势函数示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

8.4 band能带数据处理

知识点：

脚本中调用get_band_data()读取数据，在读取数据时设置efermi=XX可将能量零点修改为指定值；设置zero_to_efermi=True, 可将能量零点修改问所读取的文件中的费米能级处
脚本调用pymatgen的BSPlotter.get_plot()画图，在画图时可设置zero_to_efermi=True，将坐标轴能量零点修改到费米能级处。由于pymatgen在2023.8.17的一处关键更新将绘图函数返回对象从plt改成了axes，导致后续脚本可能出现不兼容。因此用户脚本相关部分已添加一个判断语句加以处理。
能带两步算需从第一步自洽计算获取准确费米能级（从自洽获得的system.json），若获取失败，用户可在调用get_band_data读取数据时，利用参数efermi修改能量零点。例如：band_data=get_band_data(‘band.h5’,efermi=-1.5)
脚本中画图调用pymatgen中的BSPlotter.get_plot, 当判断体系为非金属时，设置zero_to_efermi会认为VBM为费米能级能量而非文件读取时的费米能级。故当体系为非金属时，在数据读取时设置zero_to_efermi=True和在画图时设置zero_to_efermi=True得到的图会有区别

执行本节所列的python脚本，程序会判断是否为金属体系。若为非金属体系，将要求选择是否平移费米能级到能量零点，请按提示操作。

普通能带处理

参考 4bandplot.py ：

import os

import matplotlib.pyplot as plt
from dspawpy.io.read import get_band_data
from pymatgen.electronic_structure.plotter import BSPlotter

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # 指定数据文件路径
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=True,  # 非金属体系应将零点能移动到费米能级
)

bsp = BSPlotter(band_data)
axes_or_plt = bsp.get_plot(
    zero_to_efermi=False,  # 上面读取数据时已经平移，此处应关闭
    ylim=[-10, 10],  # 能带图的纵坐标范围
    smooth=False,  # 是否对能带图进行平滑处理
    vbm_cbm_marker=False,  # 是否在能带图中标记导带和价带
    smooth_tol=0,  # 平滑处理的阈值
    smooth_k=3,  # 平滑处理的阶数
    smooth_np=100,  # 平滑处理的点数
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# 添加能量零点参考线
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/4bandplot.png"  # 输出的能带图文件名
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

知识点：

能带两步算需从自洽计算获取准确的费米能级（从system.json获取），若获取失败，用户可在 get_band_data 函数中指定 efermi 参数修改

执行代码可以得到类似以下能带图：

二硫化钼能带示意图

将能带投影到每一种元素分别作图，数据点大小表示该元素对该轨道的贡献

参考 4bandplot_elt.py ：

import os

import matplotlib.pyplot as plt
import numpy as np
from dspawpy.io.read import get_band_data
from pymatgen.electronic_structure.plotter import BSPlotterProjected

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # 指定数据文件路径
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=True,  # 非金属体系应将零点能移动到费米能级
)

bsp = BSPlotterProjected(bs=band_data)  # 初始化BSPlotterProjected类
axes_or_plt = bsp.get_elt_projected_plots(
    zero_to_efermi=False,  # 上面读取数据时已经平移，此处应关闭
    ylim=[-8, 5],  # 设置能量范围
    vbm_cbm_marker=False,  # 是否标记导带底和价带顶
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# 添加能量零点参考线
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = (
    "dspawpy_proj/dspawpy_tests/outputs/us/4bandplot_elt.png"  # 输出的能带图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

知识点：

用户如果需要绘制能带投影的数据，此时需要使用 BSPlotterProjected模块
使用 BSPlotterProjected模块中 get_elt_projected_plots 函数能够绘制每种元素对轨道贡献的能带图

执行代码可以得到类似以下能带图：

二硫化钼元素投影能带示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

能带投影到不同元素的不同轨道

参考 4bandplot_elt_orbit.py ：

import os

import matplotlib.pyplot as plt
import numpy as np
from dspawpy.io.read import get_band_data
from pymatgen.electronic_structure.plotter import BSPlotterProjected

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # 指定数据文件路径
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=True,  # 非金属体系应将零点能移动到费米能级
)

bsp = BSPlotterProjected(bs=band_data)  # 初始化BSPlotterProjected类
# ! 选定元素与轨道，创建字典
dict_elem_orbit = {"Mo": ["d"], "S": ["s"]}

axes_or_plt = bsp.get_projected_plots_dots(
    dictio=dict_elem_orbit,
    zero_to_efermi=False,  # 上面读取数据时已经平移，此处应关闭
    ylim=[-8, 5],  # 设置能量范围
    vbm_cbm_marker=False,  # 是否标记导带底和价带顶
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# 添加能量零点参考线
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/4bandplot_elt_orbit.png"  # 输出的能带图文件名
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

知识点：

使用 BSPlotterProjected模块中 get_projected_plots_dots可以让用户来自定义需要绘制的某种元素某种轨道(L)的能带图
例如 get_projected_plots_dots ({‘Mo’:[‘d’],’S’:[‘s’]})就是绘制Mo的d轨道和S的s轨道

执行代码可以得到类似以下能带图：

二硫化钼元素轨道投影能带示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

将能带投影到不同原子的不同轨道

参考 4bandplot_patom_porbit.py ：

import os

import matplotlib.pyplot as plt
import numpy as np
from dspawpy.io.read import get_band_data
from pymatgen.electronic_structure.plotter import BSPlotterProjected

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # 指定数据文件路径
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=True,  # 非金属体系应将零点能移动到费米能级
)

bsp = BSPlotterProjected(bs=band_data)
# ! 指定元素和轨道和原子序号
dict_elem_orbit = {"Mo": ["px", "py", "pz"]}
dict_elem_index = {"Mo": [1]}

axes_or_plt = bsp.get_projected_plots_dots_patom_pmorb(
    dictio=dict_elem_orbit,  # 指定元素-轨道字典
    dictpa=dict_elem_index,  # 指定元素-原子序号字典
    sum_atoms=None,  # 是否对原子求和
    sum_morbs=None,  # 是否对轨道求和
    zero_to_efermi=False,  # 上面读取数据时已经平移，此处应关闭
    ylim=None,  # 设置能量范围
    vbm_cbm_marker=False,  # 是否标记导带底和价带顶
    selected_branches=None,  # 指定绘制的能带分支
    w_h_size=(12, 8),  # 设置图片宽高
    num_column=None,  # 设置每行显示的图片数量
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# 添加能量零点参考线
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = (
    "dspawpy_proj/dspawpy_tests/outputs/us/4band_patom_porbit.png"  # 输出的能带图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

知识点：

使用 BSPlotterProjected模块中 get_projected_plots_dots_patom_pmorb 的自由度更高，可以让用户来自定义需要绘制的某些原子某些轨道的能带图
dictpa指定原子，dictio 指定该原子的轨道
如果要将一些原子或一些轨道的投影分量叠加起来，请根据 get_projected_plots_dots_patom_pmorb 函数文档指定 sum_atoms 或 sum_morbs 参数

警告

如果只选择单个轨道，且轨道名称不止一个字母（例如px、dxy、dxz），get_projected_plots_dots_patom_pmorb 函数将报错，详情见此处

执行代码可以得到类似以下能带图：

二硫化钼原子投影能带示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

能带反折叠处理

参考 4bandunfolding.py ：

import os

from dspawpy.plot import plot_bandunfolding

plt = plot_bandunfolding(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.22.1/band.h5",  # 读取数据
    ef=None,  # 费米能级，从文件中读取
    de=0.05,  # 能带宽度，默认0.05
    dele=0.06,  # 能带间隔，默认0.06
)

plt.ylim(-15, 10)
figname = (
    "dspawpy_proj/dspawpy_tests/outputs/us/4bandunfolding.png"  # 输出的能带图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
plt.savefig(figname, dpi=300)
# plt.show()

执行代码可以得到类似以下能带图：

Cu3Au能带反折叠示意图

警告

此功能暂不支持将非金属材料的费米能级设为能量零点（默认价带顶为能量零点）。

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

band-compare能带对比图处理

将普通能带和wannier能带绘制在同一张图上

参考 4bandcompare.py ：

import os

from dspawpy.io.read import get_band_data
from pymatgen.electronic_structure.plotter import BSPlotter

band_data = get_band_data(
    band_dir="dspawpy_proj/dspawpy_tests/inputs/2.30/wannier.h5",  # 瓦尼尔能带文件路径
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=False,  # 是否将零点能移动到费米能级
)
bsp = BSPlotter(bs=band_data)
band_data = get_band_data(
    band_dir="dspawpy_proj/dspawpy_tests/inputs/2.3/band.h5",  # 读取DFT能带
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=False,  # 是否将零点能移动到费米能级
)

bsp2 = BSPlotter(bs=band_data)
bsp.add_bs(bsp2._bs)
axes_or_plt = bsp.get_plot(
    zero_to_efermi=True,  # 将零点能移动到费米能级
    ylim=[-10, 10],  # 能带图的纵坐标范围
    smooth=False,  # 是否对能带图进行平滑处理
    vbm_cbm_marker=False,  # 是否在能带图中标记导带和价带
    smooth_tol=0,  # 平滑处理的阈值
    smooth_k=3,  # 平滑处理的阶数
    smooth_np=100,  # 平滑处理的点数
    bs_labels=["wannier interpolated", "DFT"],  # 能带图中的标签
)

import matplotlib.pyplot as plt  # noqa: E402

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# 添加能量零点参考线
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/4wanierBand.png"  # 输出的能带图文件名
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

执行代码可得到类似如下能带对比曲线：

晶体硅瓦尼尔能带与DFT能带示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

8.5 dos态密度数据处理

总的态密度

参考 5dosplot_total.py ：

import os

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos
from pymatgen.electronic_structure.plotter import DosPlotter

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # 读取投影态密度数据
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # 是否将费米能级作为零点
    stack=False,  # True表示绘制面积图
    sigma=None,  # 高斯展宽，None表示不进行平滑处理
)
dos_plotter.add_dos(label="total dos", dos=dos_data)  # 设置态密度图例  # 传入态密度数据

ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # 设置能量范围
    ylim=[-15, 15],  # 设置态密度范围
)
ax.axhline(0, lw=2, ls="-.", color="gray")

filename = (
    "dspawpy_proj/dspawpy_tests/outputs/us/5dos_total.png"  # 输出的态密度图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(filename, dpi=300)

知识点：

使用 get_dos_data 函数可以将DS-PAW计算得到的 dos.h5 文件转化为 pymatgen 支持的格式
使用 DosPlotter模块获取到DS-PAW计算的 dos.h5 的数据
DosPlotter函数可以传递参数：stack参数表示画态密度是否加阴影， zero_at_efermi 表示是否在态密度图中进行将费米能量置零，这里设置 stack=False , zero_at_efermi=False
使用 DosPlotter 模块中 add_dos 获取态密度的数据
DosPlotter模块中 get_plot函数绘制态密度图

执行代码可以得到类似以下态密度图：

NiO态密度示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

将态密度投影到不同的轨道上

参考 5dosplot_spd.py ：

import os

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos
from pymatgen.electronic_structure.plotter import DosPlotter

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # 读取投影态密度数据
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # 是否将费米能级作为零点
    stack=False,  # True表示绘制面积图
    sigma=None,  # 高斯展宽，None表示不进行平滑处理
)
dos_plotter.add_dos_dict(
    dos_dict=dos_data.get_spd_dos(),
    key_sort_func=None,  # 轨道投影  # 指定排序函数
)
ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # 设置能量范围
    ylim=None,  # 设置态密度范围
)

ax.axhline(0, lw=2, ls="-.", color="gray")

filename = "dspawpy_proj/dspawpy_tests/outputs/us/5dos_spd.png"  # 输出的态密度图文件名
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(filename, dpi=300)

知识点：

使用 DosPlotter模块中 add_dos_dict 函数获取投影态密度的数据，之后使用 get_spd_dos 将投影信息按照 spd 轨道投影输出

执行代码可以得到类似以下态密度图：

NiO轨道投影态密度示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

将态密度投影到不同的元素上

参考 5dosplot_elt.py ：

import os

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos
from pymatgen.electronic_structure.plotter import DosPlotter

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # 读取投影态密度数据
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # 是否将费米能级作为零点
    stack=False,  # True表示绘制面积图
    sigma=None,  # 高斯展宽，None表示不进行平滑处理
)
dos_plotter.add_dos_dict(
    dos_dict=dos_data.get_element_dos(),
    key_sort_func=None,  # 元素投影态密度  # 指定排序函数
)

ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # 设置能量范围
    ylim=None,  # 设置态密度范围
)
ax.axhline(0, lw=2, ls="-.", color="gray")

filename = "dspawpy_proj/dspawpy_tests/outputs/us/5dos_elt.png"  # 输出的态密度图文件名
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(filename, dpi=300)

知识点：

使用 DosPlotter模块中 add_dos_dict 函数获取投影态密度的数据，之后使用 get_element_dos 将投影信息按照不同元素投影输出

执行代码可以得到类似以下态密度图：

NiO元素投影态密度示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

将态密度投影到不同原子的不同轨道上

参考 5dosplot_atom_orbit.py ：

import os

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos
from pymatgen.electronic_structure.core import Orbital
from pymatgen.electronic_structure.plotter import DosPlotter

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # 读取投影态密度数据
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # 是否将费米能级作为零点
    stack=False,  # True表示绘制面积图
    sigma=None,  # 高斯展宽，None表示不进行平滑处理
)

# ! 指定原子序号和轨道
dict_index_orbit = {0: ["dxy"], 2: ["s"]}

print("正在绘图...")
for index in dict_index_orbit:
    _os = dict_index_orbit[index]
    _e = str(dos_data.structure.sites[index].species)
    for _orb in _os:
        dos_plotter.add_dos(
            f"{_e}(atom-{index}) {_orb}",  # label
            dos_data.get_site_orbital_dos(
                dos_data.structure[index], getattr(Orbital, _orb),
            ),
        )

ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # 设置能量范围
    ylim=None,  # 设置态密度范围
)
ax.axhline(0, lw=2, ls="-.", color="gray")

figname = (
    "dspawpy_proj/dspawpy_tests/outputs/us/5dos_atom_orbit.png"  # 输出的态密度图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(figname, dpi=300)

知识点：

使用 get_site_orbital_dos函数提取dos数据中特定原子特定轨道的贡献，dos_data.structure[0],Orbital(4) 表示获取第1个原子的dxy轨道的态密度，get_site_orbital_dos函数中序号从0开始
运行此脚本，根据提示选择元素和轨道，可以得到相应的态密度图

执行代码可以得到类似以下态密度图：

NiO原子轨道态密度示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

将态密度投影到不同原子的分裂d轨道(t2g, eg)上

参考 5dosplot_t2g_eg.py ：

import os

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos
from pymatgen.electronic_structure.plotter import DosPlotter

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # 读取投影态密度数据
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # 是否将费米能级作为零点
    stack=False,  # True表示绘制面积图
    sigma=None,  # 高斯展宽，None表示不进行平滑处理
)
# print(dos_data.structure)

# ! 指定原子序号，从0开始
ais = [1]

print("正在绘图...")
atom_indices = [int(ai) for ai in ais]
for atom_index in atom_indices:
    dos_plotter.add_dos_dict(
        dos_data.get_site_t2g_eg_resolved_dos(dos_data.structure[atom_index]),
    )

ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # 设置能量范围
    ylim=None,  # 设置态密度范围
)
ax.axhline(0, lw=2, ls="-.", color="gray")

filename = (
    "dspawpy_proj/dspawpy_tests/outputs/us/5dos_t2g_eg.png"  # 输出的态密度图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(filename, dpi=300)

知识点：

使用 get_site_t2g_eg_resolved_dos函数提取dos数据中特定原子的 t2g和 eg轨道的贡献，这是获取第2个原子的t2g和eg轨道的贡献
运行此脚本，根据提示选择原子序号，可以得到相应的态密度图

执行代码可以得到类似以下态密度图：

NiO分裂d轨道原子投影态密度示意图

知识点：

如果元素不含d轨道，会画成空白图片

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

d-带中心分析

以 Pb-slab 体系为例，对 Pt 原子进行 d-band 中心分析：

参考 5center_dband.py ：

from dspawpy.io.read import get_dos_data
from dspawpy.io.utils import d_band

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/supplement/dos.h5",  # 读取投影态密度数据
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
for spin in dos_data.densities:
    print("spin=", spin)
    c = d_band(spin, dos_data)
    print(c)

执行代码可以得到类似以下结果：

spin=1
-1.785319344084034

备注

目前仅支持全部原子平均的 d 轨道中心，不支持元素、原子投影或其他轨道，也不支持选择自旋方向或能量范围。

get_dos_data 函数负责处理态密度数据：

8.6 bandDos能带和态密度共同显示

以应用教程中Si体系为例：

将能带和态密度显示在一张图上

参考 6bandDosplot.py ：

import os

import numpy as np
from dspawpy.io.read import get_band_data, get_dos_data
from matplotlib.axes import Axes
from pymatgen.electronic_structure.plotter import BSDOSPlotter

bandfile = "dspawpy_proj/dspawpy_tests/inputs/2.3/band.h5"  # 普通能带数据
band_data = get_band_data(
    band_dir=bandfile,
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=True,  # 非金属体系应将零点能移动到费米能级
)
dosfile = "dspawpy_proj/dspawpy_tests/inputs/2.5/dos.h5"  # 态密度数据
dos_data = get_dos_data(
    dos_dir=dosfile,
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
bdp = BSDOSPlotter(
    bs_projection=None,  # 能带的投影方式，None表示不投影
    dos_projection=None,  # 态密度的投影方式，None表示不投影
    vb_energy_range=4,  # 价带能量范围
    cb_energy_range=4,  # 导带能量范围
    fixed_cb_energy=False,  # 是否固定导带能量范围
    egrid_interval=1,  # 能量格点间隔
    font="DejaVu Sans",  # 默认Times New Roman，修改成 DejaVu Sans，可避免在linux上因字体未安装而出现的警告
    axis_fontsize=20,  # 坐标轴字体大小
    tick_fontsize=15,  # 刻度字体大小
    legend_fontsize=14,  # 图例字体大小
    bs_legend="best",  # 能带图图例位置
    dos_legend="best",  # 态密度图图例位置
    rgb_legend=True,  # 是否使用彩色图例
    fig_size=(11, 8.5),  # 图片大小
)
axes_or_plt = bdp.get_plot(bs=band_data, dos=dos_data)  # 传入能带数据  # 传入态密度数据

if isinstance(axes_or_plt, Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = (
    "dspawpy_proj/dspawpy_tests/outputs/us/6bandDos.png"  # 输出的能带-态密度图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)

执行代码可以得到类似以下能带态密度图：

单晶硅能带-态密度示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

将能带和投影态密度显示在一张图上

参考 6bandPdosplot.py ：

import os

import numpy as np
from dspawpy.io.read import get_band_data, get_dos_data
from matplotlib.axes import Axes
from pymatgen.electronic_structure.plotter import BSDOSPlotter

bandfile = "dspawpy_proj/dspawpy_tests/inputs/2.4/band.h5"  # 普通能带数据
band_data = get_band_data(
    band_dir=bandfile,
    syst_dir=None,  # system.json文件路径，仅当band_dir为json文件时需要
    efermi=None,  # 用于手动修正费米能级
    zero_to_efermi=True,  # 非金属体系应将零点能移动到费米能级
)
dosfile = "dspawpy_proj/dspawpy_tests/inputs/2.6/dos.h5"  # 投影态密度数据
dos_data = get_dos_data(
    dos_dir=dosfile,
    return_dos=False,  # 如果为False，则统一返回CompleteDos对象（无论计算时是否开了投影）
)
bdp = BSDOSPlotter(
    bs_projection="elements",  # 能带的投影方式，None表示不投影
    dos_projection="elements",  # 投影态密度到元素上
    vb_energy_range=4,  # 价带能量范围
    cb_energy_range=4,  # 导带能量范围
    fixed_cb_energy=False,  # 是否固定导带能量范围
    egrid_interval=1,  # 能量格点间隔
    font="DejaVu Sans",  # 默认Times New Roman，修改成 DejaVu Sans，可避免在linux上因字体未安装而出现的警告
    axis_fontsize=20,  # 坐标轴字体大小
    tick_fontsize=15,  # 刻度字体大小
    legend_fontsize=14,  # 图例字体大小
    bs_legend="best",  # 能带图图例位置
    dos_legend="best",  # 投影态密度图图例位置
    rgb_legend=True,  # 是否使用彩色图例
    fig_size=(11, 8.5),  # 图片大小
)
axes_or_plt = bdp.get_plot(
    bs=band_data, dos=dos_data,
)  # 传入能带数据  # 传入投影态密度数据

if isinstance(axes_or_plt, Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = "dspawpy_proj/dspawpy_tests/outputs/us/6bandPdos.png"  # 输出的能带-投影态密度图文件名
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)

执行代码可以得到类似以下能带态密度图：

单晶硅能带-投影态密度示意图

警告

给定投影能带数据，默认将沿着元素投影；给定普通能带数据（或体系所含元素种类超过4种），则不投影，并输出警告
给定投影态密度数据，默认也将沿着元素投影，可以换成沿着轨道投影，或者不投影；给定普通态密度数据且未关闭态密度投影选项 BSDOSPlotter(dos_projection=None)，pymatgen绘图程序将报错，这就是上面特意准备了一个 6bandDosplot.py 文件的原因。

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

8.7 optical光学性质数据处理

以快速入门Si体系光学性质计算得到的 scf.h5 为例（注：输出文件名与task一致，task = scf，io.optical = true可计算光学性质）：

反射率数据处理，参考 7optical.py ：

from dspawpy.plot import plot_optical

plot_optical(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.12/scf.h5",
    keys=["ExtinctionCoefficient", "Reflectance"],
    axes=["X"],  # ["X", "Y", "Z", "XY", "YZ", "ZX"]
    prefix="",  # 保存到哪个文件夹中，若为空，表示当前文件夹
    save=True,  # 是否工具性质名称保存图片，若为 False，请参考下方脚本自行保存
)
# 以上函数将绘制并分别保存 ExtinctionCoefficient 和 Reflectance 的图像
# 要将多个性质绘制在同一张图中，请取消注释以下代码，并将上方 save 参数设置为 False
# import os
# import matplotlib.pyplot as plt
#
# plt.tick_params(labelsize=16)
# plt.tight_layout()
# filename = (
#     "dspawpy_proj/dspawpy_tests/outputs/us/7optical.png"  # 输出的光学性质图文件名
# )
# os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
# plt.savefig(filename, dpi=300)

知识点：

Reflectance为光学性质中的一种，用户可以根据自己的需求将该关键词修改为“AbsorptionCoefficient”或”ExtinctionCoefficient”或”RefractiveIndex”，分别对应吸收系数、消光系数和折射率

执行代码可以得到类似以下反射率随能量变化的曲线：

单晶硅Reflectance随光子能量变化趋势示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

8.8 neb过渡态计算数据处理

以快速入门 H 在 Pt(100) 表面扩散为例进行介绍：

输入文件之生成中间构型

参考 8neb_interpolate_structures.py ：

from dspawpy.diffusion.neb import NEB, write_neb_structures
from dspawpy.diffusion.nebtools import write_json_chain
from dspawpy.io.structure import read

# 读取初态构型
init_struct = read("dspawpy_proj/dspawpy_tests/inputs/2.15/00/structure00.as")[0]
# 读取末态构型
final_struct = read("dspawpy_proj/dspawpy_tests/inputs/2.15/04/structure04.as")[0]

neb = NEB(
    initial_structure=init_struct,  # 初态构型
    final_structure=final_struct,  # 末态构型
    nimages=8,  # 共8个构型，包括初末态
)
structures = neb.linear_interpolate()  # 线性插值
# structures = neb.idpp_interpolate()   #idpp插值

# 保存 as 结构文件到 dest 路径下
write_neb_structures(
    structures=structures,  # 插入插值后的构型链
    coords_are_cartesian=True,  # 是否以笛卡尔坐标保存
    fmt="as",  # 保存格式，支持'json', 'as', 'hzw', 'pdb', 'xyz', 'dump' 六种
    path="dspawpy_proj/dspawpy_tests/outputs/us/8neb_interpolate_structures",  # 保存路径
    prefix="structure",  # 文件名前缀
)

# 预览初插构型链
write_json_chain(
    preview=True,  # 是否为预览模式
    directory="dspawpy_proj/dspawpy_tests/outputs/us/8neb_interpolate_structures",  # NEB计算目录
    step=-1,  # 默认保存最后一个离子步（最新）的构型链
    dst="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # 保存路径
)
# write_xyz_chain(preview=True, # 是否为预览模式
#                  directory="dspawpy_proj/dspawpy_tests/outputs/us/8neb_interpolate_structures", # NEB计算目录
#                  step=-1, # 默认保存最后一个离子步（最新）的构型链
#                  dst='dspawpy_proj/dspawpy_tests/outputs/us/8neb' # 保存路径
# )

知识点：

用户可以根据需要自行修改插点个数，设置为8得到包含8个结构文件的文件夹，其中中间构型为6个
neb.linear_interpolate为线性插值方法，pbc参数为 True 时将锁定寻找最短扩散路径，默认为 False 以增加用户控制的灵活性，这是因为
举个例子，初态某原子的分数坐标是 0.2，终态变成 0.8。pbc = True 将强制设置扩散路径为 0.2 -> -0.2。pbc = False 则用户可以令程序按照 0.2 -> 0.8 的扩散路径进行插值；如果要用最短路径，手动将 0.8 改为 -0.2 即可，从而确保程序按照使用者的意图完成插值初猜。

绘制能垒图

neb.iniFin = true/false

当设置 neb.iniFin = true/false 时，都可使用读取neb计算的路径进行势垒分析（确保初末态计算文件在neb计算路径下）：

参考 8neb_barrier_CubicSpline.py ：

from dspawpy.diffusion.nebtools import plot_barrier

directory_of_neb_task = (
    "dspawpy_proj/dspawpy_tests/inputs/2.15"  # <-- 请修改成实际NEB路径
)

# 三次样条插值 (CubicSpline) 绘制能垒图
plot_barrier(
    directory=directory_of_neb_task,  # neb任务的路径
    method="CubicSpline",  # 三次样条插值
    figname="dspawpy_proj/dspawpy_tests/outputs/us/8neb_barrier_CubicSpline.png",  # 输出的能垒图文件名
    show=False,  # 是否显示能垒图
)

备注

运行上述脚本后，可得到类似以下三次样条插值的势垒曲线：

对于这个算例，三次样条插值后，曲线会出现不合理的“下坠”区域，这是三次样条插值算法的特性所决定的。

dspawpy 内部调用 scipy 的插值算法，上面这个脚本我们以三次样条插值算法为例，它在 scipy 文档中的定义为：

class scipy.interpolate.CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None)

关键词参数包括 axis, bc_type, extrapolate，具体含义见 scipy.interpolate.CubicSpline 。我们可以往 plot_barrier 函数中指定相应的关键词参数（axis, bc_type, extrapolate），将其传递给 scipy.interpolate.CubicSpline 这个类处理。

下面我们采用 8neb_barrier.py 脚本，对比三种算法插值绘制出的曲线：

import os

import matplotlib.pyplot as plt
from dspawpy.diffusion.nebtools import plot_barrier

# 对比不同插值方法绘制出的能垒曲线的区别，此时show应当设置为False
# 1. interp1d
plot_barrier(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",  # neb计算的路径
    ri=None,  # 初始构型与第二个构型之间的反应坐标，当NEB任务只计算了中间构型时需要
    rf=None,  # 最后一个构型与倒数第二个构型之间的反应坐标，当NEB任务只计算了中间构型时需要
    ei=None,  # 初始构型的能量，当NEB任务只计算了中间构型时需要
    ef=None,  # 最后一个构型的能量，当NEB任务只计算了中间构型时需要
    method="interp1d",  # 插值方法
    figname=None,  # 输出的能垒图文件名
    show=False,  # 是否显示能垒图
    kind="quadratic",  # 插值方法的参数
)
# 2. CubicSpline
plot_barrier(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",
    method="CubicSpline",
    figname=None,
    show=False,
)
# 3. pchip
plot_barrier(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",
    method="pchip",
    figname=None,
    show=False,
)

filename = "dspawpy_proj/dspawpy_tests/outputs/us/8neb_barrier_comparison.png"  # 输出的能垒图文件名
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
plt.savefig(filename, dpi=300)
# plt.show()

知识点：

选择合适的插值算法对于优化最终曲线的呈现效果至关重要
多数情况下，选择 pchip 这种单调三次样条插值算法即可达到较好效果，这也是默认调用的插值算法

neb.iniFin = true

当设置 neb.iniFin = true 时，读取neb计算所得的 neb.h5/neb.json 文件可快速进行势垒分析：

参考 8neb_barrier_CubicSpline.py ：

from dspawpy.diffusion.nebtools import plot_barrier

# 三次样条插值 (CubicSpline) 绘制能垒图
plot_barrier(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.15/neb.h5",  # neb.h5的路径
    method="CubicSpline",  # 三次样条插值
    figname="dspawpy_proj/dspawpy_tests/outputs/us/8neb_barrier_.png",  # 输出的能垒图文件名
    show=False,  # 是否显示能垒图
)

处理得到的势垒图与之前读取路径效果一致。

备注

neb.h5 和 neb.json文件所存能量为TotalEnergy，如需得精确的势垒值，建议采用读取neb计算路径的方法进行处理（取TotalEnergy0）

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

过渡态计算数据处理

NEB计算完成后，一般要画出能垒图观察，并检查各插值构型最终的受力大小，从而确保受力小于指定的阈值。如果结果异常，还需要检查各插值构型在结构优化过程中的受力与能量变化趋势，判断是否真正“收敛”。上述这些操作至少需要三次循环，为简化操作，我们提供了一个一步到位的总结函数 summary ：

参考 8neb_check_results.py ：

from dspawpy.diffusion.nebtools import summary

# 导入neb计算目录，需提供neb计算完成后的完整文件夹
summary(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",
    show_converge=False,  # 是否显示能量和受力收敛过程图
    outdir="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # 能量和受力收敛过程图保存路径
    figname="dspawpy_proj/dspawpy_tests/outputs/us/8neb/neb_barrier_summary.png",  # 能垒图保存路径
)
# 可以进一步设置其他关键字参数，用于绘制能垒图，例如：
# summary(directory = 'dspawpy_proj/dspawpy_tests/inputs/2.15', method='CubicSpline') # 改为三次样条插值CubicSpline

知识点：

此脚本将以表格形式打印各构型的能量和受力、绘制能垒图，并绘制中间构型的能量与受力收敛过程图
若 neb.iniFin = false ，用户必须将自洽计算的结果文件 scf.h5 或 system.json 文件复制到对应的初末态子文件夹中，否则程序无法读取初末态能量和受力信息，将报错退出
默认情况下，能垒图存储在neb计算目录外层，各中间构型的能量与受力收敛图存储在各子文件夹

执行代码可以得到类似如下NEB计算各构型的能量和受力表格：

Image	Force (eV/Å)	Reaction coordinate (Å)	Energy (eV)	Delta energy (eV)
00	0.1803	0.0000	-39637.0984	0.0000
01	0.0263	0.5428	-39637.0186	0.0798
02	0.0248	1.0868	-39636.8801	0.2183
03	0.2344	1.5884	-39636.9984	0.1000
04	0.0141	2.0892	-39637.0900	0.0084

除了可以获得能垒图外，还可以得到各中间构型的能量和受力收敛曲线（以02构型为例）

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

观察NEB链

此处所说NEB链指的是各插值构型（structure00.as, structure01.as, …）之间的几何关系，而不是单个构型在优化过程中的变化。

NEB任务计算量较大，观察NEB链有助于判断NEB计算的收敛速度。另外，通过插值算法生成中间结构后，经常需要预览NEB链。这些需求可以通过 8neb_visualize.py 脚本实现：

from dspawpy.diffusion.nebtools import write_json_chain, write_xyz_chain

# 将NEB计算路径下的构型链转成json格式文件
write_json_chain(
    preview=False,  # NEB计算如果已经完成，则可以不是预览模式
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",  # NEB计算目录
    step=-1,  # 默认保存最后一个离子步（最新）的构型链
    dst="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # 保存路径
    ignorels=False,  # 设为True则忽略 latestStructureXX.as 文件
)

# 将NEB计算路径下的构型链转成xyz格式文件
write_xyz_chain(
    preview=False,  # NEB计算如果已经完成，则可以不是预览模式
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",  # NEB计算目录
    step=-1,  # 默认保存最后一个离子步（最新）的构型链
    dst="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # 保存路径
    ignorels=False,  # 设为True则忽略 latestStructureXX.as 文件
)

知识点：

此脚本生成neb_movie*.json文件后，通过 Device Studio –> Simulator –> DS-PAW –> Analysis Plot 打开 json 文件即可观察
directory指定为NEB计算主路径，需提供neb计算完成后的完整文件夹
该脚本支持处理正在进行中的（即未完成的）neb计算文件，方便用户监测实时轨迹
xyz文件可通过 OVITO 软件打开查看：通过 Device Studio –> Simulator –> OVITO 打开可视化界面，将xyz文件拖入即可
结构信息读取优先级：latestStructureXX.as > h5 > json；设置ignorels为True后，先尝试读取 h5 中的数据，失败则读取 json 中的数据

计算构型间距

参考这个脚本 8calc_dist.py ：

from dspawpy.diffusion.nebtools import get_distance
from dspawpy.io.structure import read

# 请根据实际情况修改 structure01.as 和 structure02.as 结构文件的路径
# 先读取两个构型的分数坐标、元素列表和晶胞信息
s1 = read("dspawpy_proj/dspawpy_tests/inputs/2.15/01/structure01.as")[0]
s2 = read("dspawpy_proj/dspawpy_tests/inputs/2.15/02/structure02.as")[0]
# 计算两个构型的距离，注意这个函数只接受分数坐标
dist = get_distance(
    spo1=s1.frac_coords,
    spo2=s2.frac_coords,
    lat1=s1.lattice.matrix,
    lat2=s2.lattice.matrix,
)
print("两个构型的距离为：", dist, "Angstrom")

neb续算

如果需对neb进行续算，可参考 8neb_restart.py ：

import os
from shutil import copytree, rmtree

from dspawpy.diffusion.nebtools import restart

if os.path.isdir("dspawpy_proj/dspawpy_tests/outputs/us/neb4bk"):
    rmtree("dspawpy_proj/dspawpy_tests/outputs/us/neb4bk")

copytree(
    "dspawpy_proj/dspawpy_tests/inputs/2.15",
    "dspawpy_proj/dspawpy_tests/outputs/us/neb4bk",
)
restart(
    directory="dspawpy_proj/dspawpy_tests/outputs/us/neb4bk",  # NEB任务路径
    output="dspawpy_proj/dspawpy_tests/outputs/us/8neb_restart",  # 备份到哪
)

具体效果见 neb过渡态计算续算说明。

neb计算过程中能量和最大原子受力的变化趋势图

要查看neb计算过程中能量和最大原子受力的变化趋势图，可参考 8neb_energy_force_curves.py ：

from dspawpy.diffusion.nebtools import monitor_force_energy

# 指定NEB计算文件夹路径，运行后会在指定目录生成 Energies.png 和 MaxForce.png 图片
unfinished_neb_folder = "dspawpy_proj/dspawpy_tests/inputs/supplement/neb_unfinished"
monitor_force_energy(
    directory=unfinished_neb_folder,
    outdir="imgs",  # 输出图片路径
)

将生成能量和受力变化趋势图：

8.9 phonon声子计算数据处理

以MgO体系的声子能带态密度计算得到的 phonon.h5 为例：

如果未安装 phonopy，运行下列脚本时会弹出 no module named 'phonopy' 信息，不影响程序正常运行

声子能带数据处理

参考 9phonon_bandplot.py ：

import os

from dspawpy.io.read import get_phonon_band_data
from pymatgen.phonon.plotter import PhononBSPlotter

band_data = get_phonon_band_data(
    "dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.h5",
)  # 读取声子能带
bsp = PhononBSPlotter(band_data)
axes_or_plt = bsp.get_plot(ylim=None, units="thz")  # y轴范围  # 单位
import matplotlib.pyplot as plt  # noqa: E402

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif isinstance(axes_or_plt, tuple):
    fig = axes_or_plt[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = "dspawpy_proj/dspawpy_tests/outputs/us/9phonon_bandplot.png"  # 输出的声子能带图文件名
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)

执行代码可以得到类似以下声子能带曲线：

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

声子态密度数据处理

参考 9phonon_dosplot.py ：

import os

from dspawpy.io.read import get_phonon_dos_data
from pymatgen.phonon.plotter import PhononDosPlotter

dos = get_phonon_dos_data("dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.h5")
dp = PhononDosPlotter(
    stack=False,  # True表示绘制面积图
    sigma=None,  # 高斯模糊参数
)
dp.add_dos(label="Phonon", dos=dos)  # 图例  # 要绘制的声子态密度
axes_or_plt = dp.get_plot(
    xlim=[0, 20],  # x轴范围
    ylim=None,  # y轴范围
    units="thz",  # 单位
)
import matplotlib.pyplot as plt  # noqa: E402

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif isinstance(axes_or_plt, tuple):
    fig = axes_or_plt[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = (
    "dspawpy_proj/dspawpy_tests/outputs/us/9phonon_dosplot.png"  # 输出的能垒图文件名
)
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)

执行代码可以得到类似以下声子态密度曲线：

单晶硅声子态密度示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

声子热力学数据处理

可以参考 9phonon_thermal.py ：

from dspawpy.plot import plot_phonon_thermal

plot_phonon_thermal(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.26/phonon.h5",  # phonon.h5数据文件路径
    figname="dspawpy_proj/dspawpy_tests/outputs/us/9phonon.png",  # 输出的声子热力学图文件名
    show=False,  # 是否显示图片
)

执行代码可以得到类似以下声子热力学曲线：

单晶硅声子热力学性质示意图

API: get_phonon_band_data(), get_phonon_dos_data(), plot_phonon_thermal()

get_phonon_band_data 函数负责读取声子能带：
dspawpy.io.read.get_phonon_band_data(phonon_band_dir: str, verbose: bool = False)
读取h5或json文件中的声子能带数据，构建PhononBandStructureSymmLine对象

参数

phonon_band_dir (str) – 能带文件路径，phonon.h5 / phonon.json 或包含这两个文件的文件夹

返回类型

PhononBandStructureSymmLine

实际案例
>>> from dspawpy.io.read import get_phonon_band_data >>> band_data = get_phonon_band_data("dspawpy_proj/dspawpy_tests/inputs/2.16/phonon.h5") # 读取声子能带 >>> band_data = get_phonon_band_data("dspawpy_proj/dspawpy_tests/inputs/2.16/phonon.json") # 读取声子能带

get_phonon_dos_data 函数负责读取声子态密度：

dspawpy.io.read.get_phonon_dos_data(phonon_dos_dir: str, verbose: bool = False)

读取h5或json文件中的声子态密度数据，构建PhononDos对象

参数: phonon_dos_dir (str) – 声子态密度文件路径，phonon_dos.h5 / phonon_dos.json 或包含这两个文件的文件夹
返回类型: PhononDos

实际案例

>>> from dspawpy.io.read import get_phonon_dos_data
>>> phdos = get_phonon_dos_data(phonon_dos_dir='dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.json')
>>> phdos = get_phonon_dos_data(phonon_dos_dir='dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.h5')
>>> phdos.frequencies
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ,
        1.1,  1.2,  1.3,  1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2. ,  2.1,
        2.2,  2.3,  2.4,  2.5,  2.6,  2.7,  2.8,  2.9,  3. ,  3.1,  3.2,
        3.3,  3.4,  3.5,  3.6,  3.7,  3.8,  3.9,  4. ,  4.1,  4.2,  4.3,
        4.4,  4.5,  4.6,  4.7,  4.8,  4.9,  5. ,  5.1,  5.2,  5.3,  5.4,
        5.5,  5.6,  5.7,  5.8,  5.9,  6. ,  6.1,  6.2,  6.3,  6.4,  6.5,
        6.6,  6.7,  6.8,  6.9,  7. ,  7.1,  7.2,  7.3,  7.4,  7.5,  7.6,
        7.7,  7.8,  7.9,  8. ,  8.1,  8.2,  8.3,  8.4,  8.5,  8.6,  8.7,
        8.8,  8.9,  9. ,  9.1,  9.2,  9.3,  9.4,  9.5,  9.6,  9.7,  9.8,
        9.9, 10. , 10.1, 10.2, 10.3, 10.4, 10.5, 10.6, 10.7, 10.8, 10.9,
       11. , 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7, 11.8, 11.9, 12. ,
       12.1, 12.2, 12.3, 12.4, 12.5, 12.6, 12.7, 12.8, 12.9, 13. , 13.1,
       13.2, 13.3, 13.4, 13.5, 13.6, 13.7, 13.8, 13.9, 14. , 14.1, 14.2,
       14.3, 14.4, 14.5, 14.6, 14.7, 14.8, 14.9, 15. , 15.1, 15.2, 15.3,
       15.4, 15.5, 15.6, 15.7, 15.8, 15.9, 16. , 16.1, 16.2, 16.3, 16.4,
       16.5, 16.6, 16.7, 16.8, 16.9, 17. , 17.1, 17.2, 17.3, 17.4, 17.5,
       17.6, 17.7, 17.8, 17.9, 18. , 18.1, 18.2, 18.3, 18.4, 18.5, 18.6,
       18.7, 18.8, 18.9, 19. , 19.1, 19.2, 19.3, 19.4, 19.5, 19.6, 19.7,
       19.8, 19.9, 20. ])

plot_phonon_thermal 函数负责绘制声子热力学性质图：
dspawpy.plot.plot_phonon_thermal(datafile: str = 'phonon.h5', figname: str = 'phonon.png', show: bool = True, raw: bool = False, verbose: bool = False)
声子热力学计算任务完成后，绘制相关物理量随温度变化曲线

phonon.h5/phonon.json -> phonon.png
参数
datafile (str) – h5或json文件路径或包含任意这些文件的文件夹，默认 ‘phonon.h5’

figname (str) – 保存图片的文件名

show (bool) – 是否弹出交互界面

raw (bool) – 是否保存绘图数据到rawphonon.csv文件
返回

figname – 图片路径，默认 ‘phonon.png’

返回类型

str
实际案例
>>> from dspawpy.plot import plot_phonon_thermal >>> plot_phonon_thermal('dspawpy_proj/dspawpy_tests/inputs/2.26/phonon.h5', figname='dspawpy_proj/dspawpy_tests/outputs/doctest/phonon_thermal_h5.png', show=False) >>> plot_phonon_thermal('dspawpy_proj/dspawpy_tests/inputs/2.26/phonon.json', figname='dspawpy_proj/dspawpy_tests/outputs/doctest/phonon_thermal_json.png', show=False, raw=True)

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

8.10 aimd分子动力学模拟数据处理

以快速入门 \(H_{2}O\) 分子体系分子动力学模拟得到的 aimd.h5 为例：

轨迹文件转换格式为.xyz或.dump

从aimd输出的hdf5文件中读取数据，并生成轨迹文件

生成的.xyz或.dump格式文件，可拖入OVITO观察，通过 Device Studio –> Simulator –> OVITO 打开 OVITO 可视化界面，将xyz文件或dump文件拖入OVITO即可。

参考 10write_aimd_traj.py ：

from dspawpy.io.structure import convert

convert(
    infile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # 待转化结构，如果在当前路径，可以只写文件名
    si=None,  # 筛选构型编号，如果不指定，默认读取全部
    ele=None,  # 筛选元素符号，默认读取所有元素的原子信息
    ai=None,  # 筛选原子编号，从 1 开始，默认读取所有原子信息
    outfile="dspawpy_proj/dspawpy_tests/outputs/us/10aimdTraj.xyz",  # 也可以生成 .dump 文件（精度较低），暂时只支持正交晶胞
)

执行代码将生成.xyz和.dump格式的轨迹文件，该文件可通过OVITO打开。关于结构文件转化的更多细节，请参考 structure结构转化

知识点：

OVITO 与 dspawpy 都不支持将非正交晶胞的体系保存为dump文件

动力学过程中能量、温度等变化曲线

参考 10check_aimd_conv.py ：

from dspawpy.plot import plot_aimd

plot_aimd(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # 数据文件路径
    show=False,  # 是否弹出图像窗口
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10aimd.png",  # 输出的图像文件名
    flags_str="1 2 3 4 5",  # 选择数据类型
)
# 其中flags_str的含义如下
# 1. 动能
# 2. 总能
# 3. 压力
# 4. 温度
# 5. 体积

执行代码将生成如下组合图：

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

均方位移（MSD）分析

参考 10aimd_msd.py ：

from dspawpy.analysis.aimdtools import get_lagtime_msd, plot_msd

# 如果AIMD非一次性完成，可以将多个h5文件路径以列表形式赋值给datafile参数
lagtime, msd = get_lagtime_msd(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # 数据文件路径
    select="all",  # 默认选择所有原子
    msd_type="xyz",  # 默认计算xyz方向的msd
    timestep=None,  # 默认从datafile文件中读取时间步长
)
# 用获取的数据画图并保存
plot_msd(
    lagtime,  # 横坐标
    msd,  # 纵坐标
    xlim=None,  # 设置x轴的显示范围
    ylim=None,  # 设置y轴的显示范围
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10MSD.png",  # 输出的图像文件名
    show=False,  # 是否弹出图像窗口
    #   ax = None, # 可指定子图
)

执行代码将生成类似如下图片：

均方位移（MSD）示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

均方根偏差（RMSD）分析

参考 10aimd_rmsd.py ：

from dspawpy.analysis.aimdtools import get_lagtime_rmsd, plot_rmsd

# 如果AIMD非一次性完成，可以将多个h5文件路径以列表形式赋值给datafile参数
lagtime, rmsd = get_lagtime_rmsd(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",
    timestep=None,  # 数据文件路径  # 默认从datafile文件中读取时间步长
)
plot_rmsd(
    lagtime,  # 横坐标
    rmsd,  # 纵坐标
    xlim=None,  # 设置x轴的显示范围
    ylim=None,  # 设置y轴的显示范围
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10RMSD.png",  # 输出的图像文件名
    show=False,  # 是否弹出图像窗口
    ax=None,  # 可指定子图
)

执行代码将生成类似如下图片：

均方根偏差（RMSD）示意图

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

原子对分布函数或径向分布函数（RDFs）分析

参考 10aimd_rdf.py ：

from dspawpy.analysis.aimdtools import get_rs_rdfs, plot_rdf

# 如果AIMD非一次性完成，可以将多个h5文件路径以列表形式赋值给datafile参数
rs, rdfs = get_rs_rdfs(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # 数据文件路径
    ele1="H",  # 中心元素
    ele2="O",  # 对象元素
    rmin=0.0,  # 最小半径
    rmax=10.0,  # 最大半径
    ngrid=1000,  # 网格数
    sigma=0.1,  # sigma值
)
plot_rdf(
    rs,  # 横坐标
    rdfs,  # 纵坐标
    "H",  # 中心元素
    "O",  # 对象元素
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10RDF.png",  # 图像保存路径
    show=False,  # 是否弹出图像窗口
    ax=None,  # 可指定子图
)

执行代码将生成类似如下图片：

径向分布函数（RDFs）示意图

这部分涉及的统计学计算较复杂，更多细节请参考函数API

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

8.11 Polarization铁电极化数据处理

以快速入门 \(HfO_{2}\) 体系铁电计算得到的系列 scf.h5 为例：

参考 11Ferri.py ：

from dspawpy.plot import plot_polarization_figure

plot_polarization_figure(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.20",  # 铁电极化计算路径
    repetition=2,  # 数据点绘图时重复次数
    figname="dspawpy_proj/dspawpy_tests/outputs/us/11pol.png",  # 输出的极化图文件名
    show=False,  # 是否显示极化图
)  # --> pol.png

执行代码将生成如下组合图：

12组结构对应极化数值

查看首尾构型的铁电极化数值，可以参考如下：

from dspawpy.plot import plot_polarization_figure

plot_polarization_figure(directory='.', annotation=True, annotation_style=1) # 显示首尾构型的铁电极化数值

执行代码将生成如下组合图：

12组结构对应极化数值（带首尾构型数值）

也可以使用第二种批注样式：

from dspawpy.plot import plot_polarization_figure

plot_polarization_figure(directory='.', annotation=True, annotation_style=2) # 显示首尾构型的铁电极化数值

执行代码将生成如下组合图：

12组结构对应极化数值（带首尾构型数值）

警告

如果通过SSH连接到远程服务器执行上述脚本，出现QT相关的报错信息，可能是使用的程序（比如MobaXterm等）和QT库不兼容，要么更换程序（例如VSCode或者系统自带的终端命令行），要么请在python脚本第二行开始添加以下代码：

import matplotlib
matplotlib.use('agg')

8.12 ZPE零点振动能数据处理

以快速入门 \(CO\) 体系频率计算得到的 frequency.txt 文件为例，计算零点振动能，基于以下公式：

\[ZPE=\sum_{i=1}^{3 N} \frac{h \nu_i}{2}\]

其中， \(\nu_i\) 是简正频率， \(h\) 是普朗克常数（ \(6.626\times10^{-34} J\cdot s\)）， \(N\) 是原子数。

参考 12getZPE.py ：

from dspawpy.io.utils import getZPE

# 导入频率计算得到的frequency.txt文件
getZPE(
    fretxt="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt",
)

执行代码结果文件将保存到 ZPE.dat 文件中，文件内容如下：

Data read from D:\quickstart\2.13\frequency.txt
Frequency (meV)
284.840038

--> Zero-point energy,  ZPE (eV): 0.142420019

8.13 TS的热校正能

吸附质的熵变对能量的贡献

计算基于以下公式：

\[\Delta S_{a d s}\left(0 \rightarrow T, P^0\right)=S_{v i b}=\sum_{i=1}^{3 N}\left[\frac{N_{\mathrm{A}} h \nu_i}{T\left(e^{h \nu_i / k_{\mathrm{B}} T}-1\right)}-R \ln \left(1-e^{-h \nu_i / k_{\mathrm{B}} T}\right)\right]\]

其中， \(\Delta S_{a d s}\) 表示吸附质的熵变，根据简谐近似计算。\(S_{v i b}\) 表示振动熵。 \(\nu_i\) 是简正频率，\(N_A\) 是阿伏伽德罗常数（ \(6.022\times10^{23} mol^{-1}\) ）， \(h\) 是普朗克常数（ \(6.626\times10^{-34} J\cdot s\)）， \(k_B\) 是玻尔兹曼常数（ \(1.38\times10^{-23} J\cdot K^{-1}\)）， \(R\) 是理想气体常数（ \(8.314 J\cdot mol^{-1}\cdot K^{-1}\)）， \(T\) 是体系温度，单位 \(K\)。

参考 13getTSads.py ：

from dspawpy.io.utils import getTSads

# 导入频率计算得到的frequency.txt文件，温度可自行修改
getTSads(
    fretxt="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt",
    T=298.15,
)

执行代码结果文件将保存到 TS.dat 文件中，文件内容如下：

Data read from D:\quickstart\2.13\frequency.txt
Frequency (THz)
68.873994

--> Entropy contribution,  T*S (eV): 4.7566990201851275e-06

理想气体的熵变对能量的贡献

计算基于如下公式：

\[\begin{split}\begin{aligned} S(T, P) & =S\left(T, P^{\circ}\right)-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \\ & =S_{\text {trans }}+S_{\text {rot }}+S_{\text {elec }}+S_{\text {vib }}-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \end{aligned}\end{split}\]

其中：

\[S_{\text {trans }}=k_{\mathrm{B}}\left\{\ln \left[\left(\frac{2 \pi M k_{\mathrm{B}} T}{h^2}\right)^{3 / 2} \frac{k_{\mathrm{B}} T}{P^{\circ}}\right]+\frac{5}{2}\right\}\]

\[\begin{split}S_{\mathrm{rot}}= \begin{cases}0 & \text {, monatomic } \\ k_{\mathrm{B}}\left[\ln \left(\frac{8 \pi^2 I k_{\mathrm{B}} T}{\sigma h^2}\right)+1\right] & \text {, linear } \\ k_{\mathrm{B}}\left\{\ln \left[\frac{\sqrt{\pi I_A I_{\mathrm{B}} I_{\mathrm{C}}}}{\sigma}\left(\frac{8 \pi^2 k_B T}{h^2}\right)^{3 / 2}\right]+\frac{3}{2}\right\} & \text {, nonlinear }\end{cases}\end{split}\]

\[S_{\text {elec }}=k_{\mathrm{B}} \ln [2 \times(\text { total spin })+1]\]

\[S_{\text {vib }}=k_{\mathrm{B}} \sum_i^{\text {vib DOF }}\left[\frac{\epsilon_i}{k_{\mathrm{B}} T\left(e^{\epsilon_i / k_{\mathrm{B}} T}-1\right)}-\ln \left(1-e^{-\epsilon_i / k_{\mathrm{B}} T}\right)\right]\]

其中：\(I_A\) 到 \(I_C\) 是非线性分子的三个主惯性矩， \(I\) 是线性分子的简并惯性矩， \(\sigma\) 是分子的对称数。另外，monatomic 表示单原子分子，linear 表示线性分子，nonlinear 表示非线性分子。total spin 是总自旋数。 vib DOF 表示振动自由度。

参考 13getTSgas.py 脚本处理：

from dspawpy.io.utils import getTSgas

# 从计算结果文件json或h5中读取元素和坐标
TSgas = getTSgas(
    fretxt="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt",
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.h5",
    potentialenergy=-0.0,
    geometry="linear",
    symmetrynumber=1,
    spin=1,
    temperature=298.15,
    pressure=101325.0,
)
print("Entropy contribution, T*S (eV)：", TSgas)

# 如果只有频率txt文件，可通过手动指定元素和坐标完成计算
# TSgas = getTSgas(fretxt='dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt', datafile=None, potentialenergy=-0.0, elements=['H', 'H'], geometry='linear', positions=[[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], symmetrynumber=1, spin=1, temperature=298.15, pressure=101325.0)

8.14 附录

快速下载所有脚本，点击 用户脚本.zip
dspawpy更新日志

8 辅助工具使用教程

8.1 安装与更新

更新 dspawpy

8.2 structure结构转化

8.3 volumetricData数据处理

volumetricData可视化

差分volumetricData可视化

volumetricData面平均

8.4 band能带数据处理

普通能带处理

将能带投影到每一种元素分别作图，数据点大小表示该元素对该轨道的贡献

能带投影到不同元素的不同轨道

将能带投影到不同原子的不同轨道

能带反折叠处理

band-compare能带对比图处理

8.5 dos态密度数据处理

总的态密度

将态密度投影到不同的轨道上

将态密度投影到不同的元素上

将态密度投影到不同原子的不同轨道上

将态密度投影到不同原子的分裂d轨道(t2g, eg)上

d-带中心分析

8.6 bandDos能带和态密度共同显示

将能带和态密度显示在一张图上

将能带和投影态密度显示在一张图上

8.7 optical光学性质数据处理

8.8 neb过渡态计算数据处理

输入文件之生成中间构型

绘制能垒图

neb.iniFin = true/false

neb.iniFin = true

过渡态计算数据处理

观察NEB链

计算构型间距

neb续算

neb计算过程中能量和最大原子受力的变化趋势图

8.9 phonon声子计算数据处理

声子能带数据处理

声子态密度数据处理

声子热力学数据处理

8.10 aimd分子动力学模拟数据处理

轨迹文件转换格式为.xyz或.dump

动力学过程中能量、温度等变化曲线

均方位移（MSD） 分析

均方根偏差（RMSD） 分析

原子对分布函数或径向分布函数（RDFs） 分析

8.11 Polarization铁电极化数据处理

8.12 ZPE零点振动能数据处理

8.13 TS的热校正能

吸附质的熵变对能量的贡献

理想气体的熵变对能量的贡献

8.14 附录

均方位移（MSD）分析

均方根偏差（RMSD）分析

原子对分布函数或径向分布函数（RDFs）分析