Skip to content

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

作者: 周文杰 (:postbox: 3170102144@zju.edu.cn) 单位: 西湖大学-刘仕课题组 时间: 2025-8

1. Install spglib

conda create -n fourphonon python=2.7.18

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

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

mkdir ${YourPath}/ThirdOrder cd ThirdOrder

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

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 这个环境下:

sh compile.sh

或

python setup.py build --build-lib=. --build-platlib=. # 此处是 python2, 而非 python3.

或

python setup.py install

2.2 测试 thirdorder 是否安装成功

mkdir test; cd test

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

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 \(\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 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 的安装相似. 为了教程的完整性, 这里也详细说明.

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=.

最后, 随便选择以下一种指令安装即可:

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

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

3.2 测试 Fourthorder 是否安装成功

mkdir test; cd test

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

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 \(\rm \AA\) 的截断半径

:rainbow: 三阶力常数, 一般 8~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

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) 为例.

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 的依赖:

module load intel/intel2019 module load intelmpi/2019

在组里的 V100 集群上:

module load intel/mkl/2024.2 module load intel/ifort/2024.2.1 module load intel/mpi/2021.13

在 Src 文件夹中, 键入:

make # 开始编译 arch.make

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

4.2 测试 FourPhonon 是否安装成功

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

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 则是 $\rm \AA$
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, 30*30*30 其实就非常非常大了
  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.* 时候使用的指令保持一致, 这个指令会读取当前文件夹下的原胞 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.6 lammps -c conda-forge
conda activate phono3py
conda install phono3py dynaphopy phonolammps -c conda-forge
pip install calorine

Note: 我的 DP 势场是 deepmd-kit=2.2.7 训练的, 所以可以用 deepmd-kit=2.2.6 来调用. 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 脚本构建高阶力常数. 为了操作的一致性, 这里将从构建超胞开始讲起.

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
### 将phono3py的文件，转换为ShengBTE格式的力常数 ###

from hiphive import ForceConstants
from ase.io import read

prim = read('POSCAR', format='vasp')
supercell = read('SPOSCAR', 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文件，不能通过以下指令产生:

phono3py --cfc --hdf5-compression gzip

Note: 而应该选择:

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
### 将phono3py的文件，转换为ShengBTE格式的力常数 ###

from hiphive import ForceConstants
from ase.io import read

prim = read('POSCAR', format='vasp')
supercell = read('SPOSCAR', 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 计算完成后, 就能见到如下的文件:

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>

7. 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 是同一个. 我就不再添加路径了