VASP+PHONOPY+pypolymlpj计算不同温度下声子谱，附批处理脚本-开发者社区

新版phonopy（>2.31）集成了与pypolymlp和 symfc结合进行机器学习势-自洽谐波近似（MLP-SSCHA）计算的功能，可通过少量DFT计算进行势函数生成并产生有限温度下的力常数，并以此获得有限温度下的声子谱。

教程参考：

http://phonopy.github.io/phonopy/mlp-sscha.html

软件安装

#建议新建环境conda install pypolymlp -c conda-forge -yconda activate polymlpconda install -c conda-forge symfc phonopy seekpath -y

准备工作

准备原胞结构，并进行扩胞与结构产生，其中 --rd 1000 为将产生的目标超胞数量，这里1000仅为测试，如需准确结果可增加。

phonopy --pa auto --rd 1000 -c POSCAR --dim 222 --amin 0.03 --amax 1.5

对产生的1000个结构进行VASP精度统一的自洽计算，计算完成后，使用phonopy读取计算所产生的vasprun.xml文件

其中，{001..120}为生成势所选择的计算文件，可通过人为手动更改使用数据集的大小，最好依据此来测试不同计算数据量下生成的有限温度力常数的声子谱的收敛性。

phonopy --sp -f vasprun_xmls/vasprun-{001..120}.xml

处理完成后，将产生phonopy_params.yaml文件，将其压缩为phonopy_params.yaml.xz

xz phonopy_params.yaml

如果需要测试不同数据计算量对声子谱的影响，可在phonopy_params.yaml.xz文件后加数据范围比如重命名为phonopy_params.yaml-120.xz。

本文教程使用全部1000个计算数据进行示例，所以省去区分。

使用已压缩好的数据集进行训练生成

phonopy-load phonopy_params.yaml.xz --pypolymlp --mlp-params="ntrain=100, ntest=20"

将产生polymlp.yaml文件

有限温度力常数

对于指定温度如300K，生成该温度下的力常数命令为，其中--sscha 10 对应迭代次数，--rd-temperature 300 对应温度，--rd 1000对应数据量大小。

phonopy-load phonopy_params.yaml.xz --pypolymlp --sscha 10 --rd-temperature 300 --rd 1000

将包含十次迭代和十个力常数的产生

需注意SSCHA free energy是否收敛

对于生成的十次迭代所产生的力常数，可依次计算其声子谱

for i in {1..10}; do phonopy-load phonopy_sscha_fc_$i.yaml.xz --band auto --band-points 101 -s -p; mv band.yaml band-$i.yaml; done

最后将每次迭代的声子谱汇总绘制在一起

phonopy-bandplot band-{1..10}.yaml --legend -o band_cov_300K.pdf

0K和300K下声子谱对比

不同温度下声子谱对比

批处理脚本

#!/bin/bash
# 文件名: run_sscha_temperature.sh
# 功能: 计算300-1000K每50K的温度依赖力常数并导出能带数据
# 检查必要文件
if [ ! -f "polymlp.yaml" ]; then
echo "错误：未找到 polymlp.yaml，请先完成MLP训练！"
exit 1
fi
if [ ! -f "phonopy_params.yaml.xz" ]; then
echo "错误：未找到 phonopy_params.yaml.xz！"
exit 1
fi
# 创建结果汇总目录
mkdir -p results
echo "开始计算温度范围：300K 到 1000K，步长 50K"
echo "============================================"
# 循环温度
for T in $(seq 300 50 1000); do
echo "正在处理: ${T}K ..."

# 为每个温度创建独立工作目录
WORKDIR="T${T}K"
mkdir -p ${WORKDIR}

# 复制必要文件进去
cp phonopy_params.yaml.xz ${WORKDIR}/
cp polymlp.yaml ${WORKDIR}/
[ -f "phonopy_disp.yaml" ] && cp phonopy_disp.yaml ${WORKDIR}/

# 进入工作目录
cd ${WORKDIR}

# 1. 运行 SSCHA 计算（10次迭代，1000个随机位移）
phonopy-load phonopy_params.yaml.xz \
--pypolymlp \
--sscha 10 \
--rd-temperature ${T} \
--rd 1000 > sscha_${T}K.log 2>&1

# 检查是否成功生成最终力常数
if [ ! -f "phonopy_sscha_fc_10.yaml.xz" ]; then
echo " 警告: ${T}K 计算可能失败，检查 ${WORKDIR}/sscha_${T}K.log"
cd ..
continue
fi

# 2. 为最终力常数计算能带结构（-s保存，-p绘图）
for i in {1..10}; do
if [ -f "phonopy_sscha_fc_${i}.yaml.xz" ]; then
phonopy-load phonopy_sscha_fc_${i}.yaml.xz --band auto --band-points 101 -s -p > pdf_${T}_${i}K.log 2>&1
phonopy-bandplot --gnuplot > phono.txt
cp band.yaml band_${T}K_iter${i}.yaml
cp band.pdf band_${T}K_iter${i}.pdf
cp phono.txt band_${T}K_iter${i}.txt
fi
done

# 保存所有结果到汇总目录（文件名带温度标记）
if [ -f "phonopy_sscha_fc_10.yaml.xz" ]; then
cp phonopy_sscha_fc_10.yaml.xz ../results/phonopy_sscha_fc_${T}K.yaml.xz
# cp band.yaml ../results/band_${T}K.yaml
phonopy-bandplot band-{1..10}.yaml --legend --gnuplot >phoon_${T}.txt
[ -f "phoon_${T}.txt" ] && cp phoon_${T}.txt ../results/phono_${T}K.txt
phonopy-bandplot band-{1..10}.yaml --legend -o band_${T}.pdf
[ -f "band_${T}.pdf" ] && cp band_${T}.pdf ../results/band_${T}K_sscha_cov.pdf
[ -f "band.pdf" ] && cp band.pdf ../results/band_${T}K.pdf
fi

# 返回上级目录
cd ..
echo " ${T}K 完成: results/phonopy_sscha_fc_${T}K.yaml.xz"
done
echo "============================================"
echo "全部完成！结果保存在 ./results/ 目录:"
ls -lh results/
phonopy-bandplot 0k/band.yaml aresults/band_*.yaml --legend -o ../phono-T.pdf