8 辅助工具使用教程

备注

cli 是 dspawpy （一个基于 Python ≥ 3.6 的库）的命令行交互程序，提供了十几个选项，用于快速完成常见任务。对于更复杂的任务，可参考本章第二节开始的脚本示例。

安装 dspawpy后，命令行输入 cli 即可使用：

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

正在联网检查dspawpy版本... 使用 -q 参数启动此程序可跳过检查

        成功导入 x.x.x（最新版）dspawpy

********这是dspawpy命令行交互小工具，预祝您使用愉快********

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

8.1 安装与更新

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

source /data/hzwtech/profile/dspawpy.env

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

使用 conda

推荐使用 conda，它在管理 python 包时会全局考虑依赖关系，出现底层依赖库错误的概率较低

conda install dspawpy -c conda-forge

或使用 pip

pip install dspawpy

更新 dspawpy

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

conda update dspawpy

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

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

备注

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

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

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

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

如果当前用户目录下存在 .local/lib 文件夹，Python 将优先从中寻找并导入包，这可能导致包的版本错误（即使已进入conda虚拟环境）！请尝试在运行时增加 -s 选项避免导入 .local/lib 中的包： python -s your-script.py

8.2 structure结构转化

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

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

# -*- coding: utf-8 -*-
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 支持读写的结构文件列表**
infile/outfile（输入/输出文件名模糊匹配）	infmt（输入文件格式）	outfmt（输出文件格式）	兼容性	说明
*.h5	h5	h5	读+	DS-PAW计算完成后保存的hdf5类型文件
*.json	json	json	读+ 写	DS-PAW计算完成后保存的json类型文件
*.pdb	pdb	pdb	读+ 写+	Protein Data Bank
*.as	as	as	读写	DS-PAW记载原子坐标等信息的结构文件
*.hzw	hzw	hzw	读写	DeviceStudio默认的结构文件
*.xyz	无	xyz	写+	extended-xyz类型轨迹文件
*.dump	无	dump	写+	lammps-dump类型轨迹文件
.cif/.mcif	无	cif/mcif	读写	Crystallographic Information File
.vasp/POSCAR/CONTCAR*	无	poscar	读写	VASP结构文件
.cssr	无	cssr	读写	Crystal Structure Standard Representation
.xsf	无	xsf	读写	eXtended Structural Format
rndstr.in/lat.in/bestsqs	无	mcsqs	读写	Monte Carlo Special Quasirandom Structure
.yaml/.yml	无	yaml	读写	YAML Ain’t Markup Language
inp.xml/.in/inp_	无	fleur-inpgen	读写	FLEUR 结构文件，须安装 pymatgen-io-fleur
prismatic	无	prismatic	写	用于STEM模拟的一种文件格式
*.res	无	res	写	ShelX res 结构文件

备注

读取时，infmt参数优先级高于文件名模糊匹配规则；例如，指定infmt=’h5’后，文件名可以是任意值，甚至可以是 a.json
写入时，outfmt参数优先级高于文件名模糊匹配规则；例如，指定outfmt=’h5’后，文件名可以是任意值，甚至可以是 a.json
“infile/outfile”列 * 号表示任意字符，“兼容性”列 + 号表示支持处理多个结构，否则只能处理单个结构：

警告

将结构信息写入文件时，所有列出的格式(outfmt)都支持，而文件名模糊匹配规则不一定完全适用
*.cif*/*.mcif* 及以下的匹配规则已简化，更多参数细节，请查阅 pymatgen 的 from_file() 和 to() 两个函数的源代码。
将结构信息写为 json 格式时，仅支持可视化 NEB 链任务，详见观察NEB链节相关说明
DS-PAW 输出的 neb.h5、phonon.h5、phonon.json、neb.json以及wannier.json 暂时无法读取结构信息。

8.3 volumetricData数据处理

volumetricData可视化

参考 3vis_vol.py ：

# -*- coding: utf-8 -*-
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 ：

# -*- coding: utf-8 -*-
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 ：

# -*- coding: utf-8 -*-
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势函数示意图

API: write_VESTA(), write_delta_rho_vesta(), average_along_axis()

write_VESTA 函数负责处理volumetricData可视化：

dspawpy.io.write.write_VESTA(in_filename: str, data_type: str, out_filename: str = 'DS-PAW.cube', subtype: Optional[str] = None, format: str = 'cube', compact: bool = False, inorm: bool = False, gridsize: Optional[Tuple[int]] = None)

从包含电子体系信息的json或h5文件中读取数据并写入VESTA格式的文件中

参数

in_filename (str) – 包含电子体系信息的json或h5文件路径
data_type (str) – 数据类型，支持 “rho”, “potential”, “elf”, “pcharge”, “rhoBound”
out_filename (str) – 输出文件路径, 默认 “DS-PAW.cube”
subtype (str) – 用于指定data_type的数据子类型，默认为None，将读取 potential 的 TotalElectrostaticPotential 数据
format (str) – 输出的数据格式，支持 “cube” 和 “vesta” （”vasp”），默认为 “cube”，大小写不敏感
compact (bool) – 每个格点的数据都换行，通过减少空格数量减小文件体积（不会影响VESTA软件的解析），默认为False
inorm (bool) – 是否归一化体积数据，使其总和为1，默认为False
gridsize (tuple) – 重新定义的格点数，格式为 (ngx, ngy, ngz)，默认为None，即使用原始格点数

返回

out_filename – VESTA格式的文件

返回类型

file

示例

>>> from dspawpy.io.write import write_VESTA
>>> write_VESTA("dspawpy_proj/dspawpy_tests/inputs/2.2/rho.json", "rho", out_filename='dspawpy_proj/dspawpy_tests/inputs/doctest_out/rho.cube') 
Reading ...dspawpy_proj/dspawpy_tests/inputs/2.2/rho.json...
==> ...dspawpy_proj/dspawpy_tests/inputs/doctest_out/rho.cube

>>> from dspawpy.io.write import write_VESTA
>>> write_VESTA(
...     in_filename="dspawpy_proj/dspawpy_tests/inputs/2.7/potential.h5",
...     data_type="potential",
...     out_filename="dspawpy_proj/dspawpy_tests/inputs/doctest_out/my_potential.txt",
...     subtype='TotalElectrostaticPotential', # or 'TotalLocalPotential'
...     format='txt',
...     gridsize=(50,50,50), # all integer, can be larger or less than the original gridsize
... ) 
Reading .../potential.h5...
Interpolating volumetric data...
volumetric data interpolated
==> .../my_potential.txt

>>> with open("dspawpy_proj/dspawpy_tests/inputs/doctest_out/my_potential.txt") as t:
...     contents = t.readlines()
...     for line in contents[:10]:
...         print(line.strip())
# 2 atoms
# 50 50 50 grid size
# x y z value
0.000 0.000 0.000      0.3279418
0.055 0.055 0.000     -0.0740864
0.110 0.110 0.000     -0.8811763
0.165 0.165 0.000     -2.1283865
0.220 0.220 0.000     -4.0559145
0.275 0.275 0.000     -6.8291030
0.330 0.330 0.000    -10.1550909

volumetricData 指的是随空间位置变化的物理量，如电荷密度 rho，势能函数 potential，局域电荷密度 elf，部分电荷密度 pcharge，溶剂束缚电荷密度 rhoBound 等。这些数据在 DS-PAW 中以 volumetricData 类型保存。

write_delta_rho_vesta 函数负责处理差分volumetricData可视化：
dspawpy.io.write.write_delta_rho_vesta(total, individuals, output='delta_rho.cube', format='cube', compact=False, inorm=False, gridsize: Optional[tuple] = None)
电荷密度差分可视化

DeviceStudio暂不支持大文件，临时写成可以用VESTA打开的格式
参数
total (str) – 体系总电荷密度文件路径，可以是h5或json格式

individuals (list of str) – 体系各组分电荷密度文件路径，可以是h5或json格式

output (str) – 输出文件路径，默认 “delta_rho.cube”

format (str) – 输出的数据格式，支持 “cube” 和 “vasp”，默认为 “cube”

compact (bool) – 每个格点的数据都换行，通过减少空格数量减小文件体积（不会影响VESTA软件的解析），默认为False

inorm (bool) – 是否归一化体积数据，使其总和为1，默认为False

gridsize (tuple) – 重新定义的格点数，格式为 (ngx, ngy, ngz)，默认为None，即使用原始格点数
返回

output – 电荷差分（total-individual1-individual2-…）后的电荷密度文件，

返回类型

file
示例
>>> from dspawpy.io.write import write_delta_rho_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/inputs/doctest_out/delta_rho.cube') Reading ...dspawpy_proj/dspawpy_tests/inputs/supplement/AB.h5... Reading ...dspawpy_proj/dspawpy_tests/inputs/supplement/A.h5... Reading ...dspawpy_proj/dspawpy_tests/inputs/supplement/B.h5... ==> ...dspawpy_proj/dspawpy_tests/inputs/doctest_out/delta_rho.cube
average_along_axis 函数负责处理volumetricData面平均数据：
dspawpy.plot.average_along_axis(datafile='potential.h5', task: str = 'potential', axis=2, smooth=False, smooth_frac=0.8, raw=False, subtype: Optional[str] = None, **kwargs)
绘制沿着某个轴向的物理量平均值曲线
参数
datafile (str or list or np.ndarray) – h5或json文件路径或包含任意这些文件的文件夹，默认 ‘potential.h5’

task (str) – 任务类型，可以是 ‘rho’, ‘potential’, ‘elf’, ‘pcharge’, ‘rhoBound’

axis (int) – 沿着哪个轴向绘制势能曲线, 默认2

smooth (bool) – 是否平滑, 默认False

smooth_frac (float) – 平滑系数, 默认0.8

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

subtype (str) – 用于指定task数据子类，默认None，代表绘制 Potential/TotalElectrostaticPotential

**kwargs (dict) – 其他参数, 传递给 matplotlib.pyplot.plot
返回

axes – 可传递给其他函数进行进一步处理

返回类型

matplotlib.axes._subplots.AxesSubplot
示例
>>> from dspawpy.plot import average_along_axis
读取 potential.h5 文件中的数据，绘图并保存原始绘图数据到csv文件
>>> average_along_axis(datafile='dspawpy_proj/dspawpy_tests/inputs/2.7/potential.h5', task='potential', axis=2, smooth=True, smooth_frac=0.8, raw=True) Reading .../dspawpy_proj/dspawpy_tests/inputs/2.7/potential.h5... <module 'matplotlib.pyplot' from '.../matplotlib/pyplot.py'>

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 ：

# -*- coding: utf-8 -*-
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 ：

# -*- coding: utf-8 -*-
import os
import numpy as np
import matplotlib.pyplot as plt
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 函数能够绘制每种元素对轨道贡献的能带图

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

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

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

参考 4bandplot_elt_orbit.py ：

# -*- coding: utf-8 -*-
import os
import numpy as np
import matplotlib.pyplot as plt
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轨道

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

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

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

参考 4bandplot_patom_porbit.py ：

# -*- coding: utf-8 -*-
import os
import numpy as np
import matplotlib.pyplot as plt
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 函数将报错，详情见此处

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

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

能带反折叠处理

参考 4bandunfolding.py ：

# -*- coding: utf-8 -*-
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能带反折叠示意图

警告

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

band-compare能带对比图处理

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

参考 4bandcompare.py ：

# -*- coding: utf-8 -*-
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能带示意图

8.5 dos态密度数据处理

总的态密度

参考 5dosplot_total.py ：

# -*- coding: utf-8 -*-
import os
from dspawpy.io.read import get_dos_data
from pymatgen.electronic_structure.plotter import DosPlotter
from dspawpy.plot import plot_dos

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态密度示意图

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

参考 5dosplot_spd.py ：

# -*- coding: utf-8 -*-
import os
from dspawpy.plot import plot_dos
from dspawpy.io.read import get_dos_data
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轨道投影态密度示意图

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

参考 5dosplot_elt.py ：

# -*- coding: utf-8 -*-
import os
from dspawpy.plot import plot_dos
from dspawpy.io.read import get_dos_data
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元素投影态密度示意图

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

参考 5dosplot_atom_orbit.py ：

# -*- coding: utf-8 -*-
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原子轨道态密度示意图

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

参考 5dosplot_t2g_eg.py ：

# -*- coding: utf-8 -*-
import os
from dspawpy.plot import plot_dos
from dspawpy.io.read import get_dos_data
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轨道，会画成空白图片

d-带中心分析

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

参考 5center_dband.py ：

# -*- coding: utf-8 -*-
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 ：

# -*- coding: utf-8 -*-
import os
import numpy as np
from matplotlib.axes import Axes
from dspawpy.io.read import get_band_data, get_dos_data
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="elements",  # 能带的投影方式，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)

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

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

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

参考 6bandPdosplot.py ：

# -*- coding: utf-8 -*-
import os
import numpy as np
from matplotlib.axes import Axes
from dspawpy.io.read import get_band_data, get_dos_data
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.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 文件的原因。

8.7 optical光学性质数据处理

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

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

# -*- coding: utf-8 -*-
import os

import matplotlib.pyplot as plt  # noqa: E402
from dspawpy.plot import plot_optical

# 读取数据，处理反射率数据
plot_optical(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.12/scf.h5",
    keys=["ExtinctionCoefficient", "Reflectance"],
    raw=True,
)
# 绘图细节
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)
# plt.show()

知识点：

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

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

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

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

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

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

参考 8neb_interpolate_structures.py ：

# -*- coding: utf-8 -*-
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为线性插值方法，若使用idpp插值请注释该行放开下一行

绘制能垒图

neb.iniFin = true/false

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

参考 8neb_barrier_CubicSpline.py ：

# -*- coding: utf-8 -*-
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 脚本，对比三种算法插值绘制出的曲线：

# -*- coding: utf-8 -*-
import os

import matplotlib.pyplot as plt  # noqa: E402
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 ：

# -*- coding: utf-8 -*-
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）

过渡态计算数据处理

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

参考 8neb_check_results.py ：

# -*- coding: utf-8 -*-
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构型为例）

观察NEB链

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

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

# -*- coding: utf-8 -*-
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 ：

# -*- coding: utf-8 -*-
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 ：

# -*- coding: utf-8 -*-
from dspawpy.diffusion.nebtools import restart
from shutil import copytree
import os

if not os.path.exists("dspawpy_proj/dspawpy_tests/inputs/supplement/neb4bk"):
    copytree(
        "dspawpy_proj/dspawpy_tests/inputs/2.15",
        "dspawpy_proj/dspawpy_tests/inputs/supplement/neb4bk",
    )
restart(
    directory="dspawpy_proj/dspawpy_tests/inputs/supplement/neb4bk",
    output="dspawpy_proj/dspawpy_tests/outputs/us/8neb_restart",  # NEB任务路径  # 备份到哪
)

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

8.9 phonon声子计算数据处理

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

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

声子能带数据处理

参考 9phonon_bandplot.py ：

# -*- coding: utf-8 -*-
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)

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

声子态密度数据处理

参考 9phonon_dosplot.py ：

# -*- coding: utf-8 -*-
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)

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

单晶硅声子态密度示意图

声子热力学数据处理

可以参考 9phonon_thermal.py ：

# -*- coding: utf-8 -*-
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)
读取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") # 读取声子能带

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

dspawpy.io.read.get_phonon_dos_data(phonon_dos_dir: str)

读取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.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=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/inputs/out/phonon_thermal.png', show=False) Reading .../dspawpy_proj/dspawpy_tests/inputs/2.26/phonon.h5... ==> .../dspawpy_proj/dspawpy_tests/inputs/out/phonon_thermal.png

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 ：

# -*- coding: utf-8 -*-
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 ：

# -*- coding: utf-8 -*-
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. 体积

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

均方位移（MSD）分析

参考 10aimd_msd.py ：

# -*- coding: utf-8 -*-
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）示意图

均方根偏差（RMSD）分析

参考 10aimd_rmsd.py ：

# -*- coding: utf-8 -*-
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）示意图

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

参考 10aimd_rdf.py ：

# -*- coding: utf-8 -*-
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

8.11 Polarization铁电极化数据处理

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

参考 11Ferri.py ：

# -*- coding: utf-8 -*-
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组结构对应极化数值（带首尾构型数值）

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 ：

# -*- coding: utf-8 -*-
from dspawpy.io.utils import getZPE

# 导入频率计算得到的frequency.txt文件
getZPE(
    fretxt="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt",
    outfile="dspawpy_proj/dspawpy_tests/outputs/us/12ZPE.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 ：

# -*- coding: utf-8 -*-
from dspawpy.io.utils import getTSads

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

执行代码结果文件将保存到 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]\]

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

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

参考 13getTSgas.py 脚本处理：

# -*- coding: utf-8 -*-
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,
    outfile="dspawpy_proj/dspawpy_tests/outputs/us/13TSgas.txt",
)
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)

API: getTSads(), getTSgas()

getTSads 函数负责计算吸附质熵变对能量的贡献：
dspawpy.io.utils.getTSads(fretxt: str = 'frequency.txt', T: float = 298.15, outfile: str = 'TSads.dat')
从fretext中读取数据，计算ZPE和TS

将另外保存结果到 TSads.dat 中
参数
fretxt (str) – 记录频率信息的文件所在路径, 默认当前路径下的’frequency.txt’

T (float) – 温度，单位K, 默认298.15

outfile (str) – 保存结果的文件名，默认为’TSads.dat’
返回

TS – 熵校正

返回类型

float
示例
>>> from dspawpy.io.utils import getTSads >>> TS = getTSads(fretxt='dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt', T=298.15, outfile='dspawpy_proj/dspawpy_tests/outputs/doctest_out/TSads.dat') Reading ...dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt ... Frequency (THz) 0 68.873993 ==> .../doctest_out/TSads.dat --> T*S (eV): 4.7566997225177686e-06
getTSgas 函数负责计算理想气体熵变对能量的贡献：
dspawpy.io.utils.getTSgas(fretxt='frequency.txt', datafile='.', potentialenergy=0, elements=None, geometry='linear', positions=None, symmetrynumber=1, spin=1, temperature=298.15, pressure=101325, outfile='TSgas.txt')
理想气体近似下，计算熵的能量贡献
参数
fretxt (str) – 记录频率信息的文件所在路径, 默认当前路径下的’frequency.txt’

datafile (str) – 计算结果json或h5文件或包含它们的文件夹路径, 默认当前路径；如果设置为None，则需要传入elements和positions参数

potentialenergy (float) – 势能，单位eV

elements (list) – 元素列表，如果

geometry (str) – 分子几何，monatomic, linear, nonlinear

positions (list) – 原子坐标，单位Angstrom

symmetrynumber (int) – 对称数

spin (int) – 自旋数

temperature (float) – 温度，单位K

pressure (float) – 压强，单位Pa

outfile (str) – 输出文件路径，默认当前路径下的’TSgas.txt’
返回

TSgas – 理想气体近似下，计算熵的能量贡献，单位eV

返回类型

float
示例
>>> from dspawpy.io.utils import getTSgas >>> 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, outfile='dspawpy_proj/dspawpy_tests/outputs/doctest_out/TSgas.txt') Reading ...dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.h5... ==> .../doctest_out/TSgas.txt >>> print(TSgas) 0.8515317035550232 >>> 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, outfile='dspawpy_proj/dspawpy_tests/outputs/doctest_out/TSgas.txt') Reading ...dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.h5... ==> .../doctest_out/TSgas.txt >>> TSgas 0.8515317035550232

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续算

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）分析