如何构造有限温的力常数¶

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

1. Dynaphopy¶

1.1 Install Dynaphopy¶

:blush: Dynaphopy 可以简单地利用 conda 安装:

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

1.2 How to use Dynaphopy¶

我们需要准备输入文件:

trajectory.lammpstrj \& FORCE_CONSTANTS \& input_file

1.2.1 traj from Lammps¶

首先, 使用原胞 POSCAR, 扩胞为 SPOSCAR, 然后转为 Lammps 的 conf.lmp 文件.

其次, 使用 Lammps 得到 Dynaphopy 所需的 trajectory.lammpstrj. 只需要在 Lammps 的 in.lmp 中写入这样的指令即可.

fix             nvt all nvt temp ${TEMP} ${TEMP} ${TAU_T}
dump            dynaphopy all custom 1 trajectory.lammpstrj x y z
dump_modify     dynaphopy sort id
dump_modify     dynaphopy format line "%16.10f %16.10f %16.10f"   
run             ${1000_ps}

Note: 不应该使用 npt 系综, Dynaphopy 目前 (2025-8-28) 仅能被用于 nvt 系综.

1.2.2 FORCE_CONSTANTS from Phonopy/Phonolammps¶

Dynaphopy 实际上是去拟合 0K 的力常数, 所以需要一个原胞 POSCAR 的 0K 力常数. 这个一般由 Phonopy/Phonolammps 获得.

Note: 可以参考 3. Phonopy 章节来获得这个文件.

1.2.3 input_file¶

Dynaphopy 需要一个 input_file 来指定一些信息.

STRUCTURE FILE POSCAR
./POSCAR # 原胞文件

FORCE CONSTANTS
./FORCE_CONSTANTS # Phonopy 得到的力常数文件

PRIMITIVE MATRIX
1 0 0 
0 1 0
0 0 1 # 和 Phonopy 的 --pa 同理

SUPERCELL MATRIX
10 0 0
0 5 0
0 0 1 # Phonopy 的力常数对应的扩胞大小

BANDS
0.0, 0.0, 0.0         0.5, 0.0, 0.0
0.5, 0.0, 0.0         0.5, 0.5, 0.0
0.5, 0.5, 0.0         0.0, 0.5, 0.0
0.0, 0.5, 0.0         0.0, 0.0, 0.0 # 自己任意指定的声子谱路径

Note: 轨迹文件使用的 SPOSCAR 的扩胞大小, 与这里的二阶力常数的扩胞大小可以不相等. Note: 建议轨迹文件的扩胞大小是二阶力常数扩胞大小的整数倍 (教程上的建议). 我自己轨迹文件的扩胞大小是 6*3*1, 二阶力常数是 10*5*1, 效果也还行.

1.2.4 running command¶

我认为常用的指令如下:

dynaphopy input_file trajectory.lammpstrj -n 1000000 -ts 0.001 -psm 3 -sfc renormalized_force_constants -sdata
# input_file: 刚刚写好的 input_file 文件
# trajectory.lammpstrj: 轨迹文件
# -n 1000000 -ts 0.001: 选择最后 1000000 步的轨迹数据, 时间步长是 0.001 ps
# -psm 3: 选择快速傅立叶变换算法, 建议用这个算法, 又快又稳定
# -sfc renormalized_force_constants: 将拟合好的含温力常数保存为 renormalized_force_constants
# -sdata: 将含温的 linewidth 保存为 quasiparticle_data.yaml

Note: 理论上, 有了含温声子谱, 就有了含温下的声子群速度. 再加上含温的声子线宽. 可以自己写程序直接计算热导率. Note: 拟合的时间长度, 我个人感觉, 1ns 是够用的.

你可以选择使用 shell 脚本, 使用超算来做 Dynaphopy 的计算:

#!/bin/bash -l

conda activate phono3py

dynaphopy input_file trajectory.lammpstrj -n 1000000 -ts 0.001 -psm 3 -sfc renormalized_force_constants -sdata

1.3 Dynaphopy 参考链接¶

Dynaphopy的功能介绍 <Dynaphopy's Features description>
Dynaphopy和Lammps <Example of Si>

2. TDEP¶

2.1 Install TDEP¶

2.1.1 conda 环境¶

首先构建一个 conda 环境:

conda create -n TDEP python=3.10
conda activate TDEP
conda install gfortran openmpi-mpifort scalapack lapack hdf5 fftw spglib=2.4.0 -c conda-forge
pip install https://github.com/flokno/tools.tdep/archive/refs/tags/v0.0.5.zip # 这个是 TDEP 的一个工具, 和 TDEP 的本体安装无关

Note: tdeptools 太老了, 依赖 spglib=2.4.0. 且对 python 版本也有要求.

which python
> /home/liushiLab/zhouwenjie/.conda/envs/TDEP/bin/python
realpath /home/liushiLab/zhouwenjie/.conda/envs/TDEP/bin
> /home/liushiLab/zhouwenjie/.conda/envs/TDEP

Note: 这里记一下 TDEP 环境的位置: /home/liushiLab/zhouwenjie/.conda/envs/TDEP

我这里选择安装 TDEP 的并行版本, 使用的是 mpifort. 检查一下 gcc 和 mpifort 是否安装成功.

# 查看 gcc
which gcc
> ~/.conda/envs/TDEP/bin/gcc
gcc -v
# 查看 mpifort
which mpifort
> ~/.conda/envs/TDEP/bin/mpifort
mpifort -v

2.1.2 TDEP 安装¶

接下来, 直接 git clone 获得 TDEP 的安装包:

mkdir TDEP; cd TDEP
git clone https://github.com/tdep-developers/tdep.git

# 拷贝编译文件
cp examples/build/important_settings.conda ./important_settings

接下来要对这个 important_settings 做一些修改:

- PREFIX=/path/to/your/conda/environment
+ PREFIX=/home/liushiLab/zhouwenjie/.conda/envs/TDEP # 刚刚记下的 conda 环境路径

- FORTRAN_COMPILER="$PREFIX/bin/gfortran"
+ FORTRAN_COMPILER="$PREFIX/bin/mpifort" # 2025-8-29 的默认是 mpifort, 也就是说自己不用做这个修改

# 你需要检查你的 conda 环境的 lib 库, 如果其中有 mpifh, 就不改. 如果你安装的是 mpich-mpifort, 就得改了
# 下面的mpi_mpifh 和 mpifort 都是 mpifort 的一种
- MPI_LIBS="-lmpi_mpifh -lmpi"
+ MPI_LIBS="-lmpifort -lmpi"

# 以下内容也是需要检查你的 conda 环境中到底有什么. 如果有 scalapack, 就不改. 如果你安装的是 blas, 就得改了
- BLASLAPACK_LIBS="-llapack -lscalapack"
+ BLASLAPACK_LIBS="-llapack -lblas"

Note: TDEP 比较依赖软件, 建议多看几遍这里的软件选择. 不懂的地方可以问问 AI. Note: gfortran 和 blas, 似乎只能支持单线程. 我安装的是 openmpi-mpifort 和 scalapack, 所以能使用并行计算.

修改完成后, 在终端执行:

bash build_things.sh --nthreads_make 4

如果安装成功, 终端会出现:

Printing bashrc_tdep, append these lines to your .bashrc for stuff to work nicely
MANPATH=$MANPATH:/home/liushiLab/zhouwenjie/software/TDEP/man
PATH=$PATH:/home/liushiLab/zhouwenjie/software/TDEP/bin
TDEP_BIN_DIR=/home/liushiLab/zhouwenjie/software/TDEP/bin
export MANPATH
export PATH
alias gnuplot='gnuplot -persist'

Everything should be set up and good to go!

这里备份一下修改好的 important_settings:

#!/bin/bash
# A central place to put all the important paths. You probably have to modify this to make things work.

# prefix of your conda environment folder (in which you have ./bin, ./lib, ./include, ... directories)
#PREFIX=/path/to/your/conda/environment
PREFIX=/home/liushiLab/zhouwenjie/.conda/envs/TDEP

# the fortran compiler
FORTRAN_COMPILER="$PREFIX/bin/mpifort"
# required compiler flags
FCFLAGS="-ffree-line-length-none -std=f2008 -cpp -fallow-argument-mismatch"
# extra flags, for debugging and such
FCFLAGS_EXTRA=""
# FCFLAGS_EXTRA="-fbacktrace -fcheck=all -finit-real=nan -finit-derived -fmax-errors=10 --pedantic --warn-all"
# FCFLAGS_EXTRA="-g -fbacktrace -fcheck=all -finit-real=nan -finit-derived -fmax-errors=10 --pedantic -Wall -Wextra -pedantic -Wcast-align -Wdisabled-optimization -Wmissing-include-dirs -Wshadow -Wunused -fdiagnostics-show-option -fcheck=all -Wstrict-overflow=0 -Wrealloc-lhs"

# optimization stuff. Go all in, sometimes
# OPTIMIZATION_LEVEL="-Ofast"
OPTIMIZATION_LEVEL="-O0"
OPTIMIZATION_SENSITIVE="-O0"

# the flag that sets the default real to a double.
DOUBLE_FLAG=  # "-fdefault-real-8"
# The flag that tells the compiler where to put .o and .mod files.
MODULE_FLAG="-J"

# the header to put in python scripts.
PYTHONHEADER="#!/usr/bin/env python"

# Which gnuplot terminal to use by default.
# Choices: aqua, qt, wxt
GNUPLOTTERMINAL="qt"

# Precompiler flags. Selecting default gnuplot terminal, and make the progressbars work.
PRECOMPILER_FLAGS="-DGP${GNUPLOTTERMINAL} -Dclusterprogressbar -DAGRESSIVE_SANITY"
PRECOMPILER_FLAGS="-DGP${GNUPLOTTERMINAL} -Dclusterprogressbar"

# These are the BLAS/LAPACK libraries. On OSX with gfortran, use the built-in 'framework accelerate'
PATH_TO_BLASLAPACK_LIB="-L/$PREFIX/lib"
PATH_TO_BLASLAPACK_INC="-I/$PREFIX/include"
BLASLAPACK_LIBS="-llapack -lscalapack"

# I use fftw for Fourier transforms.
PATH_TO_FFTW_LIB="-L/$PREFIX/lib"
PATH_TO_FFTW_INC="-I/$PREFIX/include"
FFTW_LIBS="-lfftw3"

# Also need MPI
PATH_TO_MPI_LIB="-L/$PREFIX/lib"
PATH_TO_MPI_INC="-I/$PREFIX/include"
MPI_LIBS="-lmpi_mpifh -lmpi"

# I also use HDF5 every now and then
PATH_TO_HDF5_LIB="-L/$PREFIX/lib"
PATH_TO_HDF5_INC="-I/$PREFIX/include"
HDF5_LIBS="-lhdf5 -lhdf5_fortran"

2.1.3 安装报错整理¶

2.1.3.1 git 包缺失¶

如果出现:

fatal: 不是 git 仓库（或者直至挂载点 /home 的任何父目录停止在文件系统边界（未设置 GIT_DISCOVERY_ACROSS_FILESYSTEM）。

这是因为安装包内没有 git 的版本信息和哈希值. 是缺失了两个文件夹 .git 和 github. 解决办法如下 (任选其一):

从官网下载之前版本的 TDEP 的 .zip 文件, 在里面找到这两个文件夹, 然后上传到安装目录.
找到 build_things.sh 的约第55-56行:

- gitbranch=`git rev-parse --abbrev-ref HEAD`
- gitrevision=`git rev-parse HEAD`
+ gitbranch="unknown"
+ gitrevision="unknown" # 相当于不进行版本号的检测

2.1.4 TDEP 安装参考链接¶

如何安装TDEP? <如何安装TDEP?>
TDEP 官方安装教程 <官方安装>

2.2 How to use TDEP¶

2.2.1 TDEP 的轨迹采样文件¶

Note: TDEP 对 POSCAR 的对称性要求非常高. Note: 如果你有一个 Direct 坐标的 POSCAR, 你需要先使用 tdeptools 的 ase_geometry_refine 得到一个高对称性的结构, 然后再扩胞做 MD 计算对应温度下的轨迹.

ase_geometry_refine POSCAR -t 1e-10 # 可以增大 -t, 使得软件识别出对称性

这样你就会得到一个 POSCAR.standardized, 这个就是高对称性的结构了. 可以扩胞后基于 GPUMD 计算 MD 轨迹.

ensemble    nvt_bdp 300 300 100
time_step    1
dump_thermo    100
dump_exyz 100 1 1 1 # GPUMD 的 run.in 文件需要这个, 输出: 位置, 速度, 力, 能量. 得到 dump.xyz
run    1000000

Note: TDEP 并不像 Dynaphopy, 需要一个连续的轨迹. TDEP 一般仅需要 300 个左右的轨迹采样结构. Note: TDEP 当然也可以使用 Lammps 产生轨迹采样文件, 但是我尝试学习的时候发现, 比 GPUMD 麻烦很多.

2.2.2 TDEP 计算脚本¶

#!/bin/bash

# 激活 Conda 环境
conda activate TDEP

MANPATH=$MANPATH:/home/liushiLab/zhouwenjie/software/TDEP/man
PATH=$PATH:/home/liushiLab/zhouwenjie/software/TDEP/bin
TDEP_BIN_DIR=/home/liushiLab/zhouwenjie/software/TDEP/bin
export MANPATH
export PATH
alias gnuplot='gnuplot -persist'

# 启用 CUDA 感知 MPI
export OMPI_MCA_opal_cuda_support=true
export UCX_MEMTYPE_CACHE=0

# 处理dum.xyz文件
tdep_parse_output --temperature 100 dump.xyz # 将 dump.xyz 拟合为 100K
# 等待
echo "当前任务 tdep_parse_output 已经计算完成"

# 构建力常数
mpirun -np 4 extract_forceconstants --temperature 100 -rc2 12 -rc3 8 -rc4 4 -s 2 > extract_fcs.log
# extract_forceconstants --temperature 100 -rc2 12 -rc3 8 -rc4 4 -s 2 > extract_fcs.log # 单线程

# 这句指令是我认为最常用的. 你可以不选择拟合四阶力常数. 注释如下:
mpirun -np 4 extract_forceconstants --temperature 100 -rc2 12 -rc3 8 -rc4 4 -s 2 > extract_fcs.log

--temperature 100: 拟合温度为 100K -rc2 12: 拟合二阶力常数, 使用 12\(\rm \AA\) 的截断半径 -rc 8: 拟合三阶力常数, 使用 8\(\rm \AA\) 的截断半径 -rc4 4: 拟合四阶力常数, 使用 4\(\rm \AA\) 的截断半径 -s 2: 对于给定的结构文件, 每2个结构, 使用其中一个. 减少使用到的结构数.

2.2.3 TDEP 使用参考¶

TDEP 的官方说明书 <Extract forceconstants>
TDEP 的官方教程 <tdep-tutorials>
陈泽坤同学的使用教程 <Silicon_project>

3. Phonopy¶

写 Phonopy 和 Phonolammps 仅是出于教程的完整性. 如果有不懂的地方, 请先移步别人的脚本, 这一部分更像是我写给自己的说明书.

3.1 Phonopy 计算声子谱的常规流程¶

conda activate phon3py # 激活 phonopy 的使用环境
InputPath="/home/liushiLab/zhouwenjie/3.In2Se3/0.InputFiles/2.phonopy"
WorkPath="/home/liushiLab/zhouwenjie/3.In2Se3/2.Phonopy/4-4-1"

cd ${WorkPath}
mkdir -p supercells; cd supercells
cp ../POSCAR ./.
phonopy -c POSCAR -d --dim 4 4 1 --pa auto --tolerence 0.0001 # --tolerence 控制对称性精度

Note: POSCAR 应该先经过充分的 vcrelax, 确保软件能识别到对称性. 对称性越高, 需要计算的超胞越少. 因为对称性可以减少计算量.

之后 supercells 文件夹中会有一堆超胞结构. 将这些结构拿来计算 scf.

InputPath="/home/liushiLab/zhouwenjie/3.In2Se3/0.InputFiles/2.phonopy"
WorkPath="/home/liushiLab/zhouwenjie/3.In2Se3/2.Phonopy/4-4-1"

cd ${WorkPath}
mkdir -p scfs; cd scfs
for poscar_file in ../supercells/POSCAR-*
do
    echo "处理文件: $poscar_file"
    if [[ -f "$poscar_file" ]]; then
        # 提取序号
        index="${poscar_file##*-}"
        # 在scfs目录下创建对应的scf-*目录
        cd ${WorkPath}/scfs
        scf_dir="scfs/scf-${index}"
        mkdir -p "$scf_dir"
        # 复制input文件到目标目录
        cp "$poscar_file" "${scf_dir}/POSCAR"
        cp ${InputPath}/POTCAR "${scf_dir}/POTCAR"
        cp ${InputPath}/vasp.sbatch "${scf_dir}/vasp.sbatch"
        cp ${InputPath}/INCAR "${scf_dir}/INCAR"
        cp ${InputPath}/KPOINTS "${scf_dir}/KPOINTS"
        echo "已创建: ${scf_dir}/POSCAR (来自 ${poscar_file})"
        # 提交任务
        cd "${scf_dir}"
        sbatch -J "scf-${index}" vasp.sbatch
    fi
done

scf 计算完成后, 再执行以下脚本:

InputPath="/home/liushiLab/zhouwenjie/3.In2Se3/0.InputFiles/2.phonopy"
WorkPath="/home/liushiLab/zhouwenjie/3.In2Se3/2.Phonopy/4-4-1"

cd ${WorkPath}
# 构建一个 band.conf 文件
mkdir -p vaspkit; cd vaspkit
cp ../supercells/POSCAR ./POSCAR
echo -e "3\n305\n3\n" | vaspkit # 你需要自己设定 vaspkit 的输入参数, 我这里是 2D 结构的声子谱 K-Path
cd ..
# 构建 band 文件夹
mkdir -p band; cd band
cp ../supercells/POSCAR ./POSCAR
cp ../supercells/phonopy_disp.yaml ./.
phonopy -f ../scfs/*/vasprun.xml
cp ../vaspkit/KPATH.phonopy ./band.conf
# 1. 注释掉 DIM 这一行（在行首添加 #）
sed -i 's/^DIM =/#DIM =/' band.conf
# 2. 将 FORCE_CONSTANTS = READ 改为 FORCE_CONSTANTS = WRITE
sed -i 's/FORCE_CONSTANTS = READ/FORCE_CONSTANTS = WRITE/' band.conf

# 启动 phonopy 绘制声子谱图
phonopy --dim 4 4 1 -c POSCAR -p band.conf --pa auto --tolerence 0.001 -s --full-fc;
phonopy-bandplot --gnuplot band.yaml > band.dat; # 将声子谱信息保存为 txt 文件格式
rm mesh.yaml; # 删除掉占用存储空间的 mesh.yaml

Note: phonopy 想得到可用于 ShengBTE 的二阶力常数. 需要使用 --full-fc 这个 tag.

3.2 沿虚频位移法¶

实际上, 对于一个初始的 POSCAR, 出现虚频是正常的. 这里提供一种沿着虚频位移原子的方法. 对于 band.conf 文件:

# 下面三个必须注释掉
#MP = 30 30 30
#TETRAHEDRON = .TRUE.
#PDOS = 1 2, 3 4 5
EIGENVECTOR = .TRUE.
MODULATION = 2 2 1, 0 0 0 1 -2
# 2 2 1: 扩胞 2*2*1
# 0 0 0: 选择声子谱中, k = (0, 0, 0) 的点
# 1: 选择第 1 条声子能带
# -2: 沿着虚频位移的大小的参数, 绝对值越大, 位移量越大

Note: 修改好 band.conf, 就继续使用 phonopy 的绘制声子谱的指令即可. Note: 之后会产生一个 MPOSCAR, 拿这个新的结构去做 vcrelax, 重算一遍声子谱.

4. Phonolammps¶

4.1 Phonolammps 计算声子谱的常规流程¶

4.1.1 Min¶

in_min.lmp 如下:

units           metal
boundary        p p p
atom_style      atomic

neighbor        2 bin
neigh_modify    every 1 delay 0 check yes

box             tilt large
read_data       conf.lmp
change_box      all triclinic

pair_style      nep nep.txt
pair_coeff      * *
mass            1 114.818
mass            2 78.96

thermo 1000
thermo_style custom step pe lx ly lz

# 二维, 不对 z 轴做优化了
#fix 1 all box/relax  x 0.0 y 0.0 xy 0.0 vmax 0.001

# 块体， 直接上 iso 优化
#fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 0.001
fix 1 all box/relax iso 0.0 vmax 0.001

# fix 结束, 保存一下结构
write_data fix_done.data

# cg 算法跑不通的话, 就试一下 hftn 算法
#min_style hftn
#minimize 1.0e-12 1.0e-12 100000 100000
min_style cg
minimize 0 1e-10 1000000 1000000

# minimize 结束, 保存一下结构
write_data min.data

minimize 的 shell 脚本如下:

#!/bin/bash

module load intelmpi/oneapi-2022.2;
FilePath="/home/liushiLab/zhouwenjie/software/lammps/lammps-stable_23Jun2022_update4"
export PYTHONPATH="${FilePath}/python:$PYTHONPATH"
export LD_LIBRARY_PATH="${FilePath}/src${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}"


FilePath="/home/liushiLab/zhouwenjie/software/lammps/lammps-stable_23Jun2022_update4/src"
mpirun -n 1 ${FilePath}/lmp_mpi -i in_min.lmp > log.min

4.1.2 Phonolammps¶

先把刚刚 minimize 好的结构, 转为 CONTCAR, 方便后面绘制声子谱.

lmp2pos.py min.data;
sed  -i  "s/Type_0 Type_1 /In Se /g" CONTCAR;

in_phono.lmp 如下:

units           metal
boundary        p p p
atom_style      atomic

neighbor        2 bin
neigh_modify    every 1 delay 0 check yes

box             tilt large

read_data       min.data

change_box      all triclinic

pair_style      nep nep.txt
pair_coeff      * *
mass            1 114.818
mass            2 78.96

Phonolammps 的 shell 脚本如下:

#!/bin/bash

FilePath="/home/liushiLab/zhouwenjie/software/lammps/lammps-stable_23Jun2022_update4"
export PYTHONPATH="${FilePath}/python:$PYTHONPATH"
export LD_LIBRARY_PATH="${FilePath}/src${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}"

conda activate phono3py;

phonolammps in_phono.lmp --dim 6 6 1 > log.phonolammps
phonopy -c CONTCAR -s -p band.conf --dim 6 6 1
phonopy-bandplot --gnuplot band.yaml > band.dat;

rm mesh.yaml;

Note: 如果使用的是 NEP 势场, 还涉及到将 NEP 安装进 lammps, 以及将 lammps 的 python 和 lib 接入 python 环境的事情.