Skip to content

如何使用 FourPhonon 计算 晶格热导率

0. POSCAR-D

FourPhonon 最好使用 Direct 坐标的 POSCAR. 如果你手头的 POSCAR 是 Cartesian 坐标, 那么可以用以下方法转换:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
from ase.io import read, write

import sys
from ase.io import read, write

# 直接检查参数数量
if len(sys.argv) != 3:
    print("用法: python Ase.py 输入文件 输出文件")
    sys.exit(1)

atom = read(sys.argv[1], format='vasp')
atom.wrap() # 确保原子位置在 [0,1) 范围内

# 读取和写入文件（无函数封装）
write(sys.argv[2], atom, format='vasp', direct=True, vasp5=True)
print(f"已转换 {sys.argv[1]} -> {sys.argv[2]}")

使用方法是:

1
python Ase.py POSCAR POSCAR-D

1. Install spglib

1
conda create -n fourphonon python=2.7.18

thirdorder.py 和 fourthorder.py 只能在 python2 环境中使用.

1
2
3
4
5
conda activate fourphonon
conda install spglib=1.14.1 numpy=1.16.5 -c conda-forge
pip install scipy==1.2.3
cd ***/envs/fourphonon/include
mkdir spglib; cp ./spglib.h ./spglib/spglib.h

Notes: 在 conda 环境的目录中, 创建一个 spglib 文件夹. 将 spglib.h 拷贝到 spglib 中. Notes: 若 spglib 软件正确安装, 则 include 文件夹中一定有 spglib.h 这个文件.

:relaxed: 现在请记住两个目录地址.

  • 刚刚构建了 spglib 文件夹的 include 目录地址: /Users/mudian/miniconda3/envs/fourphonon/include
  • FourPhonon 环境的 lib 目录地址: /Users/mudian/miniconda3/envs/fourphonon/lib

2. Install thirdorder

2.1 安装 thirdorder

1
2
mkdir ${YourPath}/ThirdOrder
cd ThirdOrder

创建并进入 ThirdOrder, 安装 thirdorder.py.

1
git clone <https://bitbucket.org/sousaw/thirdorder.git>

Note: 建议安装1.1.2版本的, 1.1.3好像安装有问题.

:triumph: 修改 setup.py 这个文件.

1
2
- INCLUDE_DIRS = []
+ INCLUDE_DIRS = ["/Users/mudian/miniconda3/envs/fourphonon/include"]

Note: 将空列表更新为指定conda环境中的 include 文件夹路径.

1
2
- LIBRARY_DIRS = []
+ LIBRARY_DIRS = ["/Users/mudian/miniconda3/envs/fourphonon/lib"]

Note: 将空列表更新为指定conda环境中的 lib 文件夹路径.

1
2
- USE_CYTHON = True
+ USE_CYTHON = False

Note: 不使用cython编译.

:yum: 确保自己在 fourphonon 这个环境下:

1
sh compile.sh

或

1
python setup.py build --build-lib=. --build-platlib=.

Note: 此处是 python2, 而非 python3.

或

1
python setup.py install

2.2 测试 thirdorder 是否安装成功

1
mkdir test; cd test

在 test 目录下新建一个 POSCAR 文件.

1
2
3
python ../thirdorder_vasp.py sow|reap na nb nc cutoff[nm/-integer] # 这里的 python 也还是 python2
python ../thirdorder_vasp.py sow 2 2 2 -3 #扩胞 2\*2\*2, 第三近邻截断
python ../thirdorder_vasp.py sow 2 2 2 0.8 #扩胞 2\*2\*2, 8 埃的截断半径

Note: 产生超胞的指令, 任选其一即可. 可以是指定第三近邻, 也可以是指定截断半径.

成功运行后, 出现:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
.d88888b   .88888.  dP   dP   dP
88.    "' d8'   `8b 88   88   88
`Y88888b. 88     88 88  .8P  .8P
      `8b 88     88 88  d8'  d8'
d8'   .8P Y8.   .8P 88.d8P8.d8P
 Y88888P   `8888P'  8888' Y88'
ooooooooooooooooooooooooooooooooo

Writing undisplaced coordinates to 3RD.SPOSCAR
Writing displaced coordinates to 3RD.POSCAR.*

888888ba   .88888.  888888ba   88888888b
88    `8b d8'   `8b 88    `8b  88
88     88 88     88 88     88 a88aaaa
88     88 88     88 88     88  88
88    .8P Y8.   .8P 88     88  88
8888888P   `8888P'  dP     dP  88888888P
ooooooooooooooooooooooooooooooooooooooooo

Note: 我测试时发现, POSCAR使用 Direct 或 Cartesian 会有超胞数量差异. 建议仅使用 Direct.

2.3 thirdorder 参考链接

  1. VASP环境安装-thirdorder安装 <知乎>
  2. thirdorder安装教程 <个人博客>

3. Install Fourthorder

3.1 安装 Fourthorder

Fourthorder.py 和 thirdorder.py 的安装相似. 为了教程的完整性, 这里也详细说明.

1
2
mkdir FourthOrder; cd FourthOrder
git clone <https://github.com/FourPhonon/Fourthorder.git>

:wink: 同样修改 setup.py 这个文件.

1
2
- INCLUDE_DIRS=["/opt/software/spglib.1.6.4/include"]
+ INCLUDE_DIRS = ["/Users/mudian/miniconda3/envs/fourphonon/include"]

Note: 将空列表更新为指定conda环境中的 include 文件夹路径.

1
2
- LIBRARY_DIRS=["/opt/software/spglib.1.6.4/lib"]
+ LIBRARY_DIRS = ["/Users/mudian/miniconda3/envs/fourphonon/lib"]

Note: 将空列表更新为指定conda环境中的 lib 文件夹路径.

我下载的 compile.sh 中, 执行语句是:

1
~/anaconda/bin/python2.7 setup.py build --build-lib=. --build-platlib=.

如果要用 sh compile.sh 进行安装, 就需要先做修改:

1
2
- ~/anaconda/bin/python2.7 setup.py build --build-lib=. --build-platlib=.
+ python setup.py build --build-lib=. --build-platlib=.

如 thirdorder.py 所示, 随便选择以下一种指令安装即可:

1
2
3
sh compile.sh
python setup.py build --build-lib=. --build-platlib=.
python setup.py install

Note: 我在 MacOS 中, 使用 第3个指令.

3.2 测试 Fourthorder 是否安装成功

1
mkdir test; cd test

在 test 目录下新建一个 POSCAR 文件.

1
2
3
python ../Fourthorder_vasp.py sow|reap na nb nc cutoff[nm/-integer]
python ../Fourthorder_vasp.py sow 2 2 2 -1 #扩胞 2\*2\*2, 第一近邻截断
python ../Fourthorder_vasp.py sow 2 2 2 0.2 #扩胞 2\*2\*2, 2 埃的截断半径

:rainbow: 三阶力常数, 一般 5~10\(\rm \AA\) 的截断. 四阶力常数, 一般 2~4\(\rm {\AA}\) 的截断.

成功运行后, 出现:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
.d88888b   .88888.  dP   dP   dP
88.    "' d8'   `8b 88   88   88
`Y88888b. 88     88 88  .8P  .8P
      `8b 88     88 88  d8'  d8'
d8'   .8P Y8.   .8P 88.d8P8.d8P
 Y88888P   `8888P'  8888' Y88'
ooooooooooooooooooooooooooooooooo

Writing undisplaced coordinates to 4TH.SPOSCAR
Writing displaced coordinates to 4TH.POSCAR.*

888888ba   .88888.  888888ba   88888888b
88    `8b d8'   `8b 88    `8b  88
88     88 88     88 88     88 a88aaaa
88     88 88     88 88     88  88
88    .8P Y8.   .8P 88     88  88
8888888P   `8888P'  dP     dP  88888888P
ooooooooooooooooooooooooooooooooooooooooo

Note: thirdorder 或 Fourthorder 在使用后, 会产生很多超胞文件. 如果是 DFT 方法, 就是把超胞拿去算单点能, 再构建高阶力常数. Note: 后面将介绍 5.3.2 赝 VASP 法, 利用 机器学习势场 构建高阶力常数.

4. Install FourPhonon

4.1 安装 FourPhonon

1
2
3
4
FourPhonon 在 github 上有三个版本
main: 含有格式转换工具, 教程案例. 利用 "迭代法" 计算四声子的弛豫时间.
Wigner_Park: 引入了 "WTE" 计算框架, 可以计算相干热导率 (coherence/wave-like term).
sampling_mathod: 使用 "采样法" 计算四声子的弛豫时间, 节省计算时间. 但是目前 (2025-8-21) 不能计算相干热导率.

Note: 湘潭大学的欧阳滔老师 (QQ号: 375635564), 整合了一个 Wigner_Park + sampling_method 的版本. Note: OK, 今天 (2025-9-30) 拿到了他们开源的版本. 点击前往.

这里以安装 FourPhonon (main) 为例.

1
2
3
mkdir FourPhononMain; cd FourPhononMain
git clone <https://github.com/FourPhonon/FourPhonon.git>
cd ./Src; cp ../arch.make ./arch.make #将 arch.make 文件拷贝到 Src 文件夹中

接下来我们需要修改 Src 文件夹中的 arch.make 文件:

1
2
- export FFLAGS=-qopenmp -traceback -debug -O2 -static_intel
+ export FFLAGS= -debug -O2 # 使用动态链接

1
2
- export LDFLAGS=-L/home/user/REPOSITORY/spglib/lib -lsymspg
+ export LDFLAGS=-L/Users/mudian/miniconda3/envs/fourphonon/lib -lsymspg # spglib的lib库

Note: 没错, 这里也用 fourphonon 环境中的 spglib 的 lib 库.

1
2
3
4
arch.make 中有这样一句:
- export MPIFC=mpiifort # 见于 FourPhonon v1.2 (推荐)
或
- export MPIFC=mpif90 # 见于 FourPhonon v1.1

Note: ==西湖超算==上建议使用 mpiifort, 可以用 module 得到 mpiifort 的依赖:

1
2
module load intel/intel2019
module load intelmpi/2019

在组里的 V100 集群上:

1
2
3
module load intel/mkl/2024.2
module load intel/ifort/2024.2.1
module load intel/mpi/2021.13

这里存放一个 arch.make 文件, 我用的是 V100 上的 FourPhononMain.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
#export FFLAGS=-qopenmp -traceback -debug -O2 -static_intel 
export FFLAGS= -debug -O2 # 使用动态链接
#export LDFLAGS=-L/home/user/REPOSITORY/spglib/lib -lsymspg
export LDFLAGS=-L/home/liushiLab/zhouwenjie/.conda/envs/FourPhonon/lib -lsymspg # spglib的lib库
export MPIFC=mpiifort
MKL=$(MKLROOT)/lib/intel64/libmkl_lapack95_lp64.a -Wl,--start-group \
$(MKLROOT)/lib/intel64/libmkl_intel_lp64.a \
 $(MKLROOT)/lib/intel64/libmkl_sequential.a \
 $(MKLROOT)/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm
export LAPACK=$(MKL)
export LIBS=$(LAPACK)

在 Src 文件夹中, 键入:

1
make # 开始编译 arch.make

Note: 如果安装成功, 那么上一级目录中会出现 可执行文件: ShengBTE.

4.2 测试 FourPhonon 是否安装成功

:yum: FourPhononMain 文件夹中有一个 Example-Si, 我们直接进入 Example-Si/input:

1
2
ls # 查看 input 文件夹中有哪些文件
CONTROL FORCE_CONSTANTS_3RD FORCE_CONSTANTS_4TH espresso.ifc2

准备一个 runscript 文件:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
#!/bin/bash

module load intel/intel2019
module load intelmpi/2019

# spglib 的 lib 路径不在环境路径中, 需要指定
FilePath1="/Users/mudian/miniconda3/envs/fourphonon/lib"
export LD_LIBRARY_PATH=${FilePath}:$LD_LIBRARY_PATH

# 执行计算
FilePath2="/Users/mudian/Softwares/FourPhonon/FourPhononMain" # ShengBTE 可执行文件的位置
mpirun -np 20 ${FilePath2}/ShengBTE 2>BTE.err >BTE.out

计算完成后, 与 Example-Si/output 对比一下输出结果即可.

4.3 FourPhonon 参考链接

  1. ShengBTE简单介绍、安装与使用 <知乎>
  2. Linux安装ShengBTE(以本人超算为例）<知乎>
  3. ShengBTE简单介绍、安装与使用 <CSDN>
  4. ShengBTE接口 <PWMAT>

Note: 如果想在本地 (Linux / MacOS) 调用 intel oneAPI, 需要安装 Intel Parallel Studio XE 2019. :sob: Note: 我现在 (2025-8-21) 还不会本地 (Linux / MacOS) 安装并行的 intel 驱动. Note: 适用于 MacOS 的 Parallel Studio XE 2020 安装指南. 太难了, 放弃 (2025-8-23).

5. FourPhonon 计算流程

:kissing_heart: 本章主要介绍怎么准备 输入文件, 至于 输出文件 以及 后处理, 麻烦自行参考 FourPhonon's Manual.

5.1 CONTROL 文件

这里用 Wigner_Park + Sampling 的一个 In2Se3 的 CONTROL 为例. Notes: 不同 FourPhonon 版本的 CONTROL 文件有微小的差异. 建议先看看 example 文件夹中怎么写, 自己跟着模仿. Notes: 我这里只是为了讲解更多的 tag, 才选择这个版本.

5.1.1 allocations

1
2
3
4
5
&allocations
nelements=2 # 两种元素
natoms=10 # 原胞中的原子个数
ngrid(:)=40 20 1 # FourPhonon 计算使用的网格大小. 网格越大, 计算越慢.
&end

Notes: 如果是二维材料, ngrid 选择 40 20 1 来加快计算. z 轴可以只给一个网格点.

5.1.2 crystral

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
&crystal
lfactor=0.1 # FourPhonon 中使用单位为 nm, POSCAR 则是埃, 所以转换系数就是 0.1
lattvec(:,1)=4.09    0.00    0.00
lattvec(:,2)=-0.00    7.09    0.00
lattvec(:,3)=0.00    0.00   30.00
elements="In" "Se"
types=1 1 1 1 2 2 2 2 2 2
positions(:,1)=0.0008857475281388  0.3341202875578863  0.5297273046666667
...
positions(:,10)=0.5008857481540573  0.5007869546068784  0.5737025553333334
scell(:)=10 5 1 #二阶力常数的扩胞大小
&end

Notes: 这里只能写 Direct 坐标, 不能写 Cartesian 坐标. Notes: 所以一开始就推荐, FourPhonon 全程使用 Direct 坐标的 POSCAR.

5.1.3 parameters

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
&parameters
T_min=100
T_max=800
T_step=50 # 如果只想计算300K的结果, 仅保留 T=300 即可.
scalebroad=0.5 #推荐: 0.5~0.8. 默认值且最大值是 1. 数值越大, 计算越慢, 结果越准.
num_sample_process_3ph_phase_space = -1
num_sample_process_3ph = -1 # -1 表示三声子计算不使用 sampling.
num_sample_process_4ph_phase_space = 200000
num_sample_process_4ph = 200000 # 200000 表示 sampling 的采样点.
&end

Notes: 一般遇到三声子计算不收敛, CONTROL 文件可以做的修改是:

  1. 加密 ngrid, 20*20*5 其实就足够大了
  2. 增大 scalebroad, 测试 0.5 \(\rightarrow\) 0.6 \(\rightarrow\) 0.7 \(\rightarrow\) 0.8

5.1.4 flags

1
2
3
4
5
6
&flags
nonanalytic=.TRUE.
nanowires=.FALSE.
convergence=.TRUE. # 使用迭代法
four_phonon=.TRUE. # 计算四声子
&end

:persevere: Notes: 四声子计算非常非常耗时. 所以我目前 (2025-8-22) 计算四声子, 常拜托湘潭大学的汤准韵同学 (QQ号: 1207375975), 用他们编译的 Wigner_Park + Sampling 版本. 如果没有计算相干热导 (coherence term) 的需求, 也可以自行编译 Sampling 版本, 计算包含四声子的粒子热导 (population term).

5.2 FORCE_CONSTANTS_2ND 文件

FourPhonon 使用的二阶力常数, 来自 Phonopy 输出的 FOCE_CONSTANTS. 关于怎么使用 Phonopy 或 PhonoLammps, 这里不再赘述. 将 FORCE_CONSTANTS 文件拷贝到 FourPhonon 计算文件夹, 并将文件名改为 FORCE_CONSTANTS_2ND 即可.

Notes: 使用 Phonopy 时, band.conf 文件中, 写入 FORCE_CONSTANTS = WRITE, 即可得到目标文件.

1
2
3
4
phonopy --dim d1 d2 d3 -c POSCAR -p band.conf -s --full-fc

# band.conf 中: FORCE_CONSTANTS = WRITE
# --full-fc: 产生完整的力常数, 否则 phonopy 产生具有对称性的简约了的力常数

Notes: PhonoLammps 计算完成, 会自动产生一个 FORCE_CONSTANTS. :sleeping: Notes: FORCE_CONSTANTS 不要与 FORCE_SETS 搞混了. 两者并不等同.

5.3 FORCE_CONSTANTS_3RD 与 FORCE_CONSTANTS_4TH 文件

5.3.1 使用 VASP 计算

这里我直接照搬 FourPhonon 的官方 shell 脚本:

1
2
3
4
5
6
7
8
9
for i in 3RD.POSCAR.*;do
   s=$(echo $i|cut -d"." -f3) &&
   d=job-$s &&
   mkdir $d &&
   cp $i $d/POSCAR &&
   cp ~/vaspinputs/INCAR ~/vaspinputs/POTCAR ~/vaspinputs/KPOINTS $d &&
   cp ~/vaspinputs/runvasp.sh $d &&
   (cd $d && qsub runvasp.sh)
done

计算完成后, 你会有很多个 job-* 文件夹. 每个文件夹都是对应编号的 POSCAR 文件的 SCF 计算结果. 然后键入:

1
2
3
4
5
FilePath="/Users/mudian/Softwares/ThirdOrder"
第三近邻:
find job-\* -name vasprun.xml|sort -n|python \${FilePath}/thirdorder_vasp.py reap 2 2 2 -3
截断 8 埃:
find job-\* -name vasprun.xml|sort -n|python \${FilePath}/thirdorder_vasp.py reap 2 2 2 0.8

使用近邻/截断, 需要与构建 3RD.POSCAR.* 时候使用的指令保持一致. Note: 这个指令会读取当前文件夹下的原胞 POSCAR, 记得拷贝一个到当前文件夹下.

:mask: 这里用三阶举例, 四阶同理. 不过是执行文件换一下, 不再赘述.

Notes: 顺利的话, 你就得到了 FORCE_CONSTANTS_3RD 和 FORCE_CONSTANTS_4TH. Notes: 准备好 CONTROL, FORCE_CONSTANTS_2ND, FORCE_CONSTANTS_3RD 和 FORCE_CONSTANTS_4TH. 之后按照 4.2 测试 FourPhonon 是否安装成功 所述, 再准备一个 runscript, 提交计算即可. Notes: 如果不想计算四声子, 那就不用准备 FORCE_CONSTANTS_4TH. 提供正确的 CONTROL 即可继续计算.

5.3.2 使用 DP or NEP 计算 (赝 VASP 法)

5.3.2.1 创建 phono3py 环境

先构建一个用于调用 机器学习势场 (Machine Learning Potential, MLP) 的 python 环境:

1
2
3
4
conda create -n phono3py deepmd-kit=2.2.11 lammps -c conda-forge
conda activate phono3py
conda install phono3py dynaphopy phonolammps -c conda-forge
pip install calorine

Note: 紧急提示, deepmd-kit 尽量安装 2.2.7 以后的版本. Note: calorine 是用来调用 NEP 势场的. dynaphopy 可以计算含温声子谱. :sweat_smile: Note: 装conda环境, 别用 VSCode 打开的 终端. 记得用 MacOS 自带的 终端.

:rainbow: 可以用以下方法测试, phono3py 环境中的 deepmd-kit 与 lammps 是否对接成功:

1
2
3
4
5
6
7
8
9
(phono3py) mudian@mudiandeiMac 热导 % python
Python 3.11.13 | packaged by conda-forge | (main, Jun  4 2025, 14:48:01) [Clang 18.1.8 ] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from lammps import lammps
>>> lmp=lammps()
LAMMPS (2 Aug 2023)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
>>> # 出现这样的文本, 就算是对接成功了. 我这里用 MacOS 做演示. 超算上一样操作即可.

:kissing_closed_eyes: 接下来, 我们将用 python 脚本构建高阶力常数. 为了操作的一致性, 这里将从构建超胞开始讲起.

1
2
3
mkdir FakeVasp; cd FakeVasp
cp ~/input/POSCAR ./. # 确保这个是充分优化后的结构, 声子谱不可以有虚频
cp ~/input/graph.pb ./. # 这个是 DP 的势场
5.3.2.2 构建三阶超胞

在已经存在 POSCAR 的 FakeVasp 文件夹中, 构建一个 1.GenSupercell.sh.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#!/bin/bash -l

conda activate FourPhonon

FilePath="/Users/mudian/Softwares/ThirdOrder"

mkdir Supercell;
cd Supercell;
cp ../POSCAR ./POSCAR:

# 如果是0.8, 就是8A. 如果是-3, 就是第三近邻.
python ${FilePath}/thirdorder_vasp.py sow 4 4 4 0.8
rm ./POSCAR
cd ..
tar -czvf Supercell.tar Supercell # 压缩包, 占用空间小
conda deactivate
5.3.2.3 FakeVasp.py 构建 vasprun.xml

运行这个脚本后, FakeVasp 下就有一个 Supercell 文件夹. 再在刚刚安装好的 phono3py 环境中运行 FakeVasp-BX-3RD.py:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
import os
import numpy as np
from ase.io import read
from deepmd.calculator import DP
from multiprocessing import Pool, cpu_count
import functools
from tqdm import tqdm
import sys
import tarfile

# 全局变量，用于存储计算器
_calculator = None

def init_calculator():
    """初始化计算器，只会在每个子进程开始时调用一次"""
    global _calculator
    _calculator = DP('graph.pb')

def create_vasprun_xml(forces, energy, output_dir="job-0001"):
    """
    生成 VASP 格式的 vasprun.xml 文件
    Args:
        forces (np.ndarray): 原子力数组 (N, 3)
        energy (float): 体系总能量
        output_dir (str): 输出目录（如 "job-0001"）
    """
    os.makedirs(output_dir, exist_ok=True)
    xml_path = os.path.join(output_dir, "vasprun.xml")

    if os.path.isdir(xml_path):
        raise ValueError(f"Path {xml_path} is a directory, cannot write file")

    with open(xml_path, "w") as f:
        f.write('<?xml version="1.0" encoding="ISO-8859-1"?>\n')
        f.write('<modeling>\n')
        f.write(' <calculation>\n')
        f.write('  <energy>\n')
        f.write(f'   <e_fr_energy>{energy:.10f}</e_fr_energy>\n')
        f.write('  </energy>\n')
        f.write('  <varray name="forces">\n')
        for force in forces:
            f.write(f'   <v> {force[0]:12.6f} {force[1]:12.6f} {force[2]:12.6f} </v>\n')
        f.write('  </varray>\n')
        f.write(' </calculation>\n')
        f.write('</modeling>\n')

def process_single_file(args):
    """
    处理单个POSCAR文件的函数，用于并行计算
    Args:
        args: 包含 (poscar_file, supercell_dir) 的元组
    Returns:
        (poscar_file, output_dir, success, error_msg)
    """
    poscar_file, supercell_dir = args
    try:
        file_num = poscar_file.split(".")[-1]
        job_id = file_num.zfill(4)
        output_dir = os.path.join("JOBS", f"job-{job_id}")

        atoms = read(os.path.join(supercell_dir, poscar_file), format="vasp")
        atoms.calc = _calculator  # 使用全局计算器
        energy = atoms.get_potential_energy()
        forces = atoms.get_forces()

        create_vasprun_xml(forces, energy, output_dir)
        return (poscar_file, output_dir, True, None)
    except Exception as e:
        return (poscar_file, None, False, str(e))

def main():
    # 0. 创建JOBS主目录
    os.makedirs("JOBS", exist_ok=True)

    # 1. 读取所有POSCAR文件
    supercell_dir = "Supercell"
    poscar_files = sorted(
        [f for f in os.listdir(supercell_dir) 
         if f.startswith("3RD.POSCAR.") and f.split(".")[-1].isdigit()],
        key=lambda x: int(x.split(".")[-1])
    )

    if not poscar_files:
        raise FileNotFoundError(f"No 3RD.POSCAR.* files found in {supercell_dir}")

    # 2. 准备并行任务参数
    num_cores = min(4, cpu_count())  # 限制最大4个核心
    tasks = [(f, supercell_dir) for f in poscar_files]

    print(f"Starting parallel processing with {min(4,num_cores)} cores...")
    print(f"Total files to process: {len(poscar_files)}")

    # 3. 并行处理带进度条
    # 使用initializer和initargs来初始化每个子进程的计算器
    with Pool(processes=min(4,num_cores), initializer=init_calculator) as pool:
        results = list(tqdm(
            pool.imap_unordered(process_single_file, tasks),
            total=len(poscar_files),
            file=sys.stdout,
            desc="Processing",
            ncols=80
        ))

    # 4. 打印处理结果
    success_count = sum(1 for r in results if r[2])
    print(f"\nProcessing summary:")
    print(f"Successfully processed: {success_count}/{len(poscar_files)}")

    # 打印失败的文件（如果有）
    failed_files = [r[0] for r in results if not r[2]]
    if failed_files:
        print("\nFailed files:")
        for file, _, _, error in [r for r in results if not r[2]]:
            print(f"{file}: {error}")

if __name__ == "__main__":
    main()
    # 处理完成后压缩文件夹
    with tarfile.open("JOBS.tar", "w") as tar:
        tar.add("JOBS", arcname=os.path.basename("JOBS"))
    print("JOBS folder has been compressed to JOBS.tar")

:rainbow: Note: 这个代码的核心部分, 就是 create_vasprun_xml 函数. 是为每一个需要被计算的超胞结构, 构建一个 job 文件夹, 然后写入一个符合 vasprun.xml 文件, 其中包含这个 POSCAR 的能量和力的信息. 之后我们将会利用这些文件, 构建力常数.

5.3.2.4 构建三阶力常数

成功运行刚刚的 python 代码后, FakeVasp 文件夹内会出现一个 JOBS 文件夹和 JOBS.tar. 再在 FakeVasp 文件夹下, 构建一个 2.GenThird.sh.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#!/bin/bash -l

conda activate FourPhonon

FilePath="/Users/mudian/Softwares/ThirdOrder"

cd JOBS
cp ../POSCAR ./POSCAR

# 这句指令, 要和刚刚构建超胞时候的指令, 保持一致性
find job-* -name vasprun.xml|sort -n|python ${FilePath}/thirdorder_vasp.py reap 4 4 4 0.8

cd ..
cp JOBS/FORCE_CONSTANTS_3RD ./FORCE_CONSTANTS_3RD
rm JOBS.tar # 把 FakeVasp-3RD.py 产生的 JOBS.tar 删除
tar -czvf JOBS.tar JOBS # 打包 JOBS 文件夹为 JOBS.tar, 这样新的 JOBS.tar 内就含有力常数等文件, 相当于更新了 JOBS 文件夹

现在我们得到了三阶力常数. 那么对于四阶力常数的处理, 也是一样的.

Note: 如果你使用的不是 DP, 而是 NEP. 那么, python 程序中需要做的改变是:

1
2
3
4
5
- from deepmd.calculator import DP
+ from calorine.calculators import CPUNEP # 引入 NEP 的计算器

- _calculator = DP('graph.pb')
+ _calculator = CPUNEP('nep.txt') # NEP 计算器调用 nep.txt 势函数

5.3.3 赝 VASP 法参考链接

  1. 秦光照老师组的github <Fourthorder-mtp>

6. 文件转换

:sunglasses: FourPhonon, Phono3py, TDEP 之间, 存在输入/输出文件互相转换的可能. 主要是我比较喜欢 phono3py 的 hdf5 文件作为后处理的输入文件. 所以这里也写一点关于格式转换的脚本.

6.1 Phono3py 的 fc2.hdf5 和 fc3.hdf5 转为 FORCE_CONSTANTS_*

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
### 将phono3py的文件，转换为ShengBTE格式的力常数 ###

from hiphive import ForceConstants
from ase.io import read

prim = read('POSCAR', format='vasp')
supercell_2nd = read('SPOSCAR_2ND', format='vasp')
supercell_3rd = read('SPOSCAR_3RD', format='vasp')

### 将二阶力常数转换为FORCECONSTANTS_2nd的格式 ###
fcs = ForceConstants.read_phonopy(supercell=supercell, fname='fc2.hdf5', format='hdf5')
fcs.write_to_phonopy(fname="FORCECONSTANTS_2ND",format="text")

### 将三阶力常数转换为FORCECONSTANTS_3rd的格式 ###
fcs = ForceConstants.read_phono3py(supercell=supercell, fname='fc3.hdf5')
fcs.write_to_shengBTE(fname="FORCECONSTANTS_3RD",prim=prim)

Note: hiphive可以转换的hdf5文件，不能通过以下指令产生:

1
phono3py --cfc --hdf5-compression gzip

Note: 而应该选择:

1
phono3py --fc-symmetry

6.2 ShengBTE 的 FORCE_CONSTANTS_* 转为 fc*.hdf5

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
### 将phono3py的文件，转换为ShengBTE格式的力常数 ###

from hiphive import ForceConstants
from ase.io import read

prim = read('POSCAR', format='vasp')
supercell_2nd = read('SPOSCAR_2ND', format='vasp')
supercell_3rd = read('SPOSCAR_3RD', format='vasp')

### 将二阶力常数转换为fc2.hdf5的格式 ###
fcs = ForceConstants.read_phonopy(supercell=supercell, fname='FORCE_CONSTANTS_2ND', format='text')
fcs.write_to_phonopy(fname="fc2.hdf5",format="hdf5")

### 将三阶力常数转换为fc3.hdf5的格式 ###
fcs = ForceConstants.read_shengBTE(supercell=supercell, fname='FORCE_CONSTANTS_3RD', prim=prim, symprec=1e-05) # symprec 是默认值, 不用改
fcs.write_to_phono3py(fname='fc3.hdf5')

6.3 TDEP 转为 ShengBTE

这里使用的是, FourPhonon 自带的格式转换工具. 在 FourPhononMain/TDEP-convert-scripts 文件夹下先编译执行文件:

1
2
ifort -o tdep2ShengBTE tdep2ShengBTE.f90
ifort -o tdep2phonopy tdep2phonopy.f90

当使用 TDEP 软件的计算指令是:

1
2
3
4
5
6
7
8
mpirun -np 4 extract_forceconstants --temperature 100 -rc2 12 -rc3 8 -rc4 4 -s 2 > extract_fcs.log

# -np 4: 使用 4 核心并行加速计算
# --temperature 100: 拟合 100K 下的力常数
# -rc2 12: 二阶力常数的截断半径为 12A
# -rc3 8: 三阶力常数的截断半径为 8A
# -rc4 4: 四阶力常数的截断半径为 4A
# -s 2: 轨迹中, 每间隔 2 step 的结构, 才用来做拟合. 假如一共有 1000 step 的轨迹结构, 仅有 1000/2=500 个用于拟合

:rainbow: Note: TDEP 一般仅需要 300 个左右的结构, 即可拟合力常数. 结构数大了算得慢.

TDEP 计算完成后, 就能见到如下的文件:

1
2
3
outfile.forceconstant # 二阶力常数输入文件
outfile.forceconstant_thirdorder # 三阶力常数输入文件
outfile.forceconstant_fourthorder # 四阶力常数输入文件

再执行下面这个 shell 脚本:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
#!/bin/bash

module load intel/intel2019
module load intelmpi/2019

# spglib 的 lib 路径不在环境路径中, 需要指定
FilePath1="/Users/mudian/miniconda3/envs/fourphonon/lib"
export LD_LIBRARY_PATH=${FilePath}:$LD_LIBRARY_PATH

# 执行计算
FilePath2="/Users/mudian/Softwares/FourPhonon/FourPhononMain" # ShengBTE 可执行文件的位置

# 6 3 1: supercell of 2nd FC; 10: 10 atoms in unit cell
${FilePath}/TDEP-convert-scripts/tdep2phonopy outfile.forceconstant 6 3 1 10

# 3: 3rd; 10: 10 atoms in unit cell
${FilePath}/TDEP-convert-scripts/tdep2ShengBTE 3 10 outfile.forceconstant_thirdorder

# 4: 4th; 10: 10 atoms in unit cell
#${FilePath}/TDEP-convert-scripts/tdep2ShengBTE 4 10 outfile.forceconstant_fourthorder

最后得到 FORCE_CONSTANTS_*_tdep, 这个可以直接作为 FourPhonon 的输入文件.

6.4 文件转换的参考链接

  1. Force constants IO <hiphive>
  2. Force constant models <hiphive>
  3. FourPhonon 的 TDEP 转换文件 <FourPhonon>
  4. mlfcs: 非谐力常数计算工具套件 <mlfcs>

Note: 这里的 mlfcs 是我在 2026-1-21 从 KALDO 微信群聊了解到的。好像可以产生三阶和四阶的 0K-力常数, 以及有限温度下的二阶力常数, 且软件安装纯 python, 不依赖 spglib 这个库. 作者好像暂时移除了 0K-二阶力常数 的计算功能. 对于非谐效应显著的体系, 建议使用 mlfcs 的 SSCHA 方法得到 0K-二阶力常数. SSCHA 其实是一种得到==有限温度力常数==的方法. Note: 说到 SSCHA, 这个方法整合在 ALAMODE 软件中. 这个软件也支持计算热导率. 但是我没有深入学习过. 读者可以自行学习 ALAMODE 软件.

7. Phono3py 计算力常数

7.0 引言

写这部分的原因是, 我在尝试做力常数转换的时候发现, ShengBTE 构建的力常数经过 hiphive 转换, 再交给 Phono3py 计算. 其结果不合理, 虽然 ShengBTE 自己计算的结果也不合理. 我不知道为什么 (后来发现是原子结构问题), 于是我认识到使用 Phono3py 构建 独立计算流程 的重要性. 此外, 我在 GPUMD 的官方教程里发现了 曾泽柱老师 的 Phono3py+NEP 计算力常数的方法. 当然其实这部分在 Phono3py 官网上的 Phono3py API 部分 亦有记载. Note: 更意外的事情是, 之前我苦于无法在 Phono3py+NEP 中加入 --cutoff-pair 这个 tag, 可我近期 (2026-1-10) 在官网发现开发者的一句话:

1
2
3
For testing, thermal conductivities with respect to --cutoff-pair values are calculated as follows.
Note that if FORCES_FC3 for full fc3 elements exists,
the same FORCES_FC3 file can be used for generating contracted fc3 for each special phono3py_disp.yaml.

意味着, 我们可以先构建完备的 FORCES_FC3, 再结合 phono3py_disp.yaml 生成 fc3.hdf5. 同理也可以产生 fc2.hdf5.

值得注意的是, 一旦 fc2.hdf5 和 fc3.hdf5 产生, 就不再需要 phono3py_disp.yaml 文件. 直接使用:

1
phono3py -c POSCAR --dim 4 4 4 --dim-fc2 6 6 6 --fc2 --fc3 --lbte --wigner --ts 300

这句计算热导率的指令里, 不再需要 --cutoff-pair 这个 tag. 官网上是这样描述的:

1
2
Once fc3.hdf5 and fc2.hdf5 are created,
--cutoff-pair option and the special phono3py_disp.yaml are not needed anymore.

7.0 参考链接

  1. Phono3py cutoff-pair <Phono3py cutoff-pair>
  2. Phono3py API <Phono3py API>
  3. 曾泽柱老师的教程 <Phono3py+NEP>

7.1 机器学习势场, 得到 FORCES_FC2 和 FORCES_FC3

假设当前目录中, 已经有 POSCAR 和 nep.txt 这两个文件. 并且已经在包含 phono3py, calorine 这两个软件的 python 环境中.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
import os
import shutil
from glob import glob

import numpy as np
import h5py

from ase import Atoms
from ase.io import read, write
from calorine.calculators import CPUNEP
from calorine.tools import relax_structure
from matplotlib import pyplot as plt

from phono3py import Phono3py
from phono3py.interface.phono3py_yaml import Phono3pyYaml
from phonopy.interface.calculator import read_crystal_structure
from phono3py.file_IO import write_fc2_to_hdf5, write_fc3_to_hdf5

# step1 设置工作目录, 读取 POSCAR 和 nep 文件
root_dir = os.getcwd()  # here we set the root directory
work_dir = os.path.join(root_dir, 'phono3py_workdir')
shutil.rmtree(work_dir, ignore_errors=True)  # Deletes current working dir
os.mkdir(work_dir)
os.chdir(work_dir)
print(f'Root directory: {root_dir}')
print(f'Working directory: {work_dir}')
potential_filename = os.path.join(root_dir, 'nep.txt')
prim = read(os.path.join(root_dir, 'POSCAR-D'), format='vasp')
calc = CPUNEP(potential_filename)

# step2 使用 phono3py 生成超胞和位移
dim2 = (6, 6, 1); dim3 = (3, 3, 1)
write('POSCAR', prim, format='vasp')
ph3yml = Phono3pyYaml() # use phono3py api to generate the displacement dataset
unitcell, _ = read_crystal_structure("POSCAR", interface_mode='vasp')
ph3 = Phono3py(unitcell, supercell_matrix=dim3, phonon_supercell_matrix=dim2, primitive_matrix='auto')
ph3.generate_displacements(); ph3.generate_fc2_displacements() # 构建二阶和三阶的位移数据集
ph3.save("phono3py_disp.yaml")

# step3 计算力常数
supercells_2nd = ph3.phonon_supercells_with_displacements; print(f'Number of 2nd order supercells: {len(supercells_2nd)}')
supercells_3rd = ph3.supercells_with_displacements; print(f'Number of 3rd order supercells: {len(supercells_3rd)}')

def read_PhonopyAtoms_to_ase_Atoms(fname: str):
    """Convert from PhonopyAtoms to Atoms"""
    structure = Atoms(symbols=fname.symbols,
                  positions=fname.positions,
                  cell=fname.cell,
                  pbc=True)
    return structure
def GenFC(supercells, calc):
    forces_data = []
    print('Calculating forces for supercells...')
    print('-------------------------------------')
    print(f'Number of supercells: {len(supercells)}')
    for it, fname in enumerate(supercells):
        structure = read_PhonopyAtoms_to_ase_Atoms(fname)
        structure.calc = calc
        forces = structure.get_forces()
        forces_data.append(forces)
        if it % 100 == 0:
            print(f'Calculating supercell {it:5d} / {len(supercells)}, f_max= {np.max(np.abs(forces)):8.5f}')
    print('-------------------------------------')
    print('Finished calculating forces for supercells.')
    return np.array(forces_data).reshape(-1, 3)

forces_2nd = GenFC(supercells_2nd, calc); np.savetxt('FORCES_FC2', forces_2nd)
forces_3rd = GenFC(supercells_3rd, calc); np.savetxt('FORCES_FC3', forces_3rd)

当成功构建出 FORCES_FC2 和 FORCES_FC3 文件后, 就可以开始构建 fc2.hdf5 和 fc3.hdf5 文件了.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
'''
下面开始构建 fc2.hdf5 和 fc3.hdf5 文件
'''

# 先进入一个临时目录，避免文件冲突
tmp_dir = os.path.join(work_dir, 'phono3py_calc')
os.mkdir(tmp_dir); os.chdir(tmp_dir)
# 将上一个目录的 FORCES_FC2 和 FORCES_FC3 复制过来
shutil.copy(os.path.join(work_dir, 'FORCES_FC2'), 'FORCES_FC2')
shutil.copy(os.path.join(work_dir, 'FORCES_FC3'), 'FORCES_FC3')
# 重新产生 POSCAR 文件
write('POSCAR_phono3py', prim, format='vasp')
cutoff_distance = 5.0 # 块体的 cutoff 距离, 设置 5A 比较合适

# step4 使用 cmd 的方法构建 phono3py_disp.yaml 文件
cmd = f'phono3py -c POSCAR_phono3py -d --dim {dim3[0]} {dim3[1]} {dim3[2]} --dim-fc2 {dim2[0]} {dim2[1]} {dim2[2]} --cutoff-pair {cutoff_distance} --pa auto'
print(cmd); os.system(cmd)
# 删除 POSCAR-*
fnames = glob('POSCAR-*')
for fname in fnames:
    os.remove(fname)
# 删除 POSCAR_FC2-*
fnames = glob('POSCAR_FC2-*')
for fname in fnames:
    os.remove(fname)

cmd = f'phono3py --fc-symmetry'
print(cmd); os.system(cmd)

# step4 利用 phono3py 的 api 构建 phono3py_disp.yaml 文件
unitcell, _ = read_crystal_structure("POSCAR_phono3py", interface_mode='vasp')
ph3 = Phono3py(unitcell, supercell_matrix=dim3, phonon_supercell_matrix=dim2, primitive_matrix='auto')
ph3.generate_displacements(cutoff_pair_distance=cutoff_distance); ph3.generate_fc2_displacements() # 构建二阶和三阶的位移数据集
ph3.save("phono3py_disp.yaml")

forces_2nd = np.loadtxt('FORCES_FC2').reshape(-1, len(ph3.phonon_supercell), 3)
forces_3rd = np.loadtxt('FORCES_FC3').reshape(-1, len(ph3.supercell), 3)
ph3.phonon_forces = forces_2nd; ph3.forces = forces_3rd

ph3.produce_fc2() #generate fc2.hdf5
write_fc2_to_hdf5(ph3.fc2) #write fc2.hdf5
ph3.produce_fc3() #generate fc3.hdf5
write_fc3_to_hdf5(ph3.fc3) #write fc3.hdf5

有了 fc2.hdf5 和 fc3.hdf5 文件后, 就可以开始计算晶格热导了. 计算的脚本如下:

1
2
3
4
#!/bin/bash

conda activate phono3py
phono3py -c POSCAR --fc2 --fc3 --lbte --mesh 15 15 5 --ts 300

Note: 显然可以注意到, 计算力常数部分, 可以利用 并行算法 进行优化. 这个就交给后来人自己去实现了. Note: BTW, 其实并行优化部分, 参考我 5.3.2 赝 VASP 法 中的代码, 交给 deepseek 分析一下也能得到吧. Note: 计算晶格热导率, 除了 phono3py 软件, 还可以利用 ALAMODE 软件进行计算. 但是 ALAMODE 也无法支持输入 Dynaphopy 的声子线宽文件, 求解 WTE 方程. 所以, 这个教程有一件没有做完的事情, 就是结合 Dynaphopy 的声子线宽文件, 利用 phono3py 软件计算晶格热导率 (2026-1-22).

8. MacOS 导出 Markdown 为 PDF

1
2
3
4
5
6
7
8
9
# MacOS 13.7.7 安装 homebrew
/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"
# 输入密码后, 按回车键
echo 'export PATH="/usr/local/bin:$PATH"' >> ~/.zshrc # intel 芯片, 将 brew 添加到路径
source ~/.zshrc # 激活路径
brew --version # 查看是否安装成功
brew install --cask prince # 安装 prince
prince --version # 验证 prince 是否安装成功
# prince 的 bin 目录和 brew 是同一个. 我就不再添加路径了