新版phonopy2.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.xzxz phonopy_params.yaml如果需要测试不同数据计算量对声子谱的影响可在phonopy_params.yaml.xz文件后加数据范围比如重命名为phonopy_params.yaml-120.xz。本文教程使用全部1000个计算数据进行示例所以省去区分。使用已压缩好的数据集进行训练生成phonopy-load phonopy_params.yaml.xz --pypolymlp --mlp-paramsntrain100, ntest20将产生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.pdf0K和300K下声子谱对比不同温度下声子谱对比批处理脚本#!/bin/bash# 文件名: run_sscha_temperature.sh# 功能: 计算300-1000K每50K的温度依赖力常数并导出能带数据# 检查必要文件if [ ! -f polymlp.yaml ]; thenecho 错误未找到 polymlp.yaml请先完成MLP训练exit 1fiif [ ! -f phonopy_params.yaml.xz ]; thenecho 错误未找到 phonopy_params.yaml.xzexit 1fi# 创建结果汇总目录mkdir -p resultsecho 开始计算温度范围300K 到 1000K步长 50Kecho # 循环温度for T in $(seq 300 50 1000); doecho 正在处理: ${T}K ...# 为每个温度创建独立工作目录WORKDIRT${T}Kmkdir -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 21# 检查是否成功生成最终力常数if [ ! -f phonopy_sscha_fc_10.yaml.xz ]; thenecho 警告: ${T}K 计算可能失败检查 ${WORKDIR}/sscha_${T}K.logcd ..continuefi# 2. 为最终力常数计算能带结构-s保存-p绘图for i in {1..10}; doif [ -f phonopy_sscha_fc_${i}.yaml.xz ]; thenphonopy-load phonopy_sscha_fc_${i}.yaml.xz --band auto --band-points 101 -s -p pdf_${T}_${i}K.log 21phonopy-bandplot --gnuplot phono.txtcp band.yaml band_${T}K_iter${i}.yamlcp band.pdf band_${T}K_iter${i}.pdfcp phono.txt band_${T}K_iter${i}.txtfidone# 保存所有结果到汇总目录文件名带温度标记if [ -f phonopy_sscha_fc_10.yaml.xz ]; thencp phonopy_sscha_fc_10.yaml.xz ../results/phonopy_sscha_fc_${T}K.yaml.xz# cp band.yaml ../results/band_${T}K.yamlphonopy-bandplot band-{1..10}.yaml --legend --gnuplot phoon_${T}.txt[ -f phoon_${T}.txt ] cp phoon_${T}.txt ../results/phono_${T}K.txtphonopy-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.pdffi# 返回上级目录cd ..echo ${T}K 完成: results/phonopy_sscha_fc_${T}K.yaml.xzdoneecho echo 全部完成结果保存在 ./results/ 目录:ls -lh results/phonopy-bandplot 0k/band.yaml aresults/band_*.yaml --legend -o ../phono-T.pdf