做材料计算的人基本都会走到声子谱这一步。我在VASP声子谱计算这条路上踩过的坑比想象中多得多——虚频满天飞、力不收敛、超胞选小了力常数还没衰减到零、装完VASP一跑就报错。这篇笔记不打算写成教科书式的流程讲解而是把我自己的实操经验按完整链路整理出来从物理图像、方法选型到Ubuntu环境准备、有限位移法完整实操再到参数对照和问题排查希望给准备跑声子谱的朋友一条能直接走通的路。1. 声子谱到底在算什么先把物理图像建起来1.1 从晶格振动说起声子不是玄学是原子的集体舞晶格里的原子并不是钉死在平衡位置上的温度高了会振动零点能也让它不可能完全静止。原子偏离平衡位置时周围原子会把它拉回来这就形成了振动模式。大量原子耦合在一起振动就不再是孤立的单个原子运动而是整块晶体里周期性排列的集体行为。把这种集体振动量子化就得到声子phonon这个概念。声子谱就是以波矢为横轴、振动频率为纵轴画出来的色散关系。它告诉我们晶格里不同波长的振动模式各自拥有多高的频率以及这些模式是往哪个方向传播的。用VASP算声子谱本质上就是通过第一性原理把原子间的受力关系算准然后换算出晶格振动的频率。声子谱里有两个关键区域必须会看。一个是声学支在Gamma点波矢为零附近的行为——对于非极性材料声学支频率在Gamma点应当趋近于零另一个是是否存在虚频负频率。所谓虚频其实就是频率平方为负从数值上表现为画在零轴下方的曲线。虚频的出现通常意味着结构处于鞍点或亚稳态原子有自发向更低能量构型移动的趋势。判断一个结构是否动力学稳定看虚频基本是一眼定生死。1.2 两条技术路线有限位移法和DFPT怎么选VASP算声子谱主流有两条路。第一条是有限位移法也叫冻声子法。做法很直观在优化的超胞里人为地给某个原子一个很小的位移通常0.01埃然后计算所有原子受到的Hellmann-Feynman力。力对位移求偏导就得到力常数矩阵再通过傅里叶变换得到动力学矩阵最终解出声子频率。因为位移是显式的所以物理图像非常清楚实现也简单配合phonopy这类工具非常顺手。第二条是DFPT也就是密度泛函微扰理论在VASP里通过IBRION7或8调用。它不是在实空间里给一个有限位移而是在倒空间自洽地求解电荷密度对原子位移的线性响应。好处是不需要超胞也能算对极性材料还能很方便地引入非解析项修正LO-TO分裂缺点是内存占用量大对参数的敏感性也更高。我自己做体系时选型逻辑大概是这样的对比维度有限位移法DFPT物理直观性高直接反映力-位移关系低黑箱子感强是否需要超胞需要且越大越稳妥不需要大超胞内存开销低高尤其对大体系极性材料处理需要额外算Born有效电荷和介电常数天然支持配合NAC修正上手难度低中等偏上所以如果你是新手上路或者体系是半导体、绝缘体、简单金属我建议先用有限位移法。等你把流程跑通了再根据体系特性考虑要不要切换到DFPT。对于极性体系比如很多钙钛矿氧化物DFPT或有限位移法NAC修正是必须的否则Gamma点光学支会有明显误差。1.3 为什么非要超胞力常数的实空间局域性有限位移法里绕不开一个概念超胞。原胞太小原子之间的距离近周期性镜像之间会互相干扰算出来的力常数包含假信号。超胞的作用就是把原子的周期镜像推远让真实的力常数在实空间内衰减到可以忽略的程度。打个比方你在操场上推一个人他倒地时会牵动身边的人。如果操场太小四周还是这个人的镜像传递的扰动就会叠加出不真实的效果。把操场扩到足够大扰动传到边界时已经衰减干净测到的才是真实的响应。超胞大小怎么选核心判断标准只有一个力常数是否已经在超胞范围内衰减到接近零。一般来说半导体和绝缘体的力常数衰减得较快3×3×3或4×4×4的超胞常常够用金属体系因为费米面的存在力常数有长程振荡往往需要更大的超胞。后面我会专门讲这个怎么测试。2. 环境准备Ubuntu下把VASP跑起来的几个关键点2.1 编译安装VASP时最容易被卡住的三件事热搜词里有“ubuntu装vasp”看来这是很多人的第一道坎。其实VASP本身编译难度不算高但配套环境不弄好后患无穷。第一件事是编译器。VASP对Intel编译器ifort的兼容性最好性能也最好。建议直接装Intel oneAPI套件里面同时包含了ifort、MKL数学库和MPI库。如果你用gfortran编译时要注意某些算法部分对gfortran的支持不如ifort稳定性能也会有差距。装好之后务必检查一下环境变量是否生效source /opt/intel/oneapi/setvars.sh which ifort which mpiifort第二件事是数学库。VASP大量调用BLAS、LAPACK、FFTW直接用MKL最省心。很多编译失败案例都是因为MKL路径没写对或者库的顺序不对。VASP编译时用的是makefile.include文件需要在源码根目录下把这个文件放好。Intel平台可以参考安装包自带的arch/makefile.include.linux_intel模板然后根据实际路径修改。第三件事是MPI。VASP并行计算依赖MPI用Intel MPI和ifort搭配最稳。如果你机器上同时装了OpenMPI和Intel MPI很容易因为mpif90指向不对导致编译错乱。建议在编译前明确mpif90 --version确保这个命令指向你期望的MPI版本。2.2 编译完了怎么验证先算个小体系VASP编译完成后不要急着直接跑声子谱。先用一个简单体系比如金刚石结构Si做一次自洽计算和能带计算验证编译结果的正确性。你不需要自己准备测试文件VASP官方测试集里就有现成的例子。跑通之后再检查OUTCAR里的总能量是否和已知结果吻合比如Si的基态能量在LDA或PBE下都有公认值。这一步的意义是把你后续排查问题的范围缩小。如果连Si都跑不对那问题大概率出在编译或环境上如果Si没问题声子谱算坏了就可以放心去查计算参数和物理设置。另外提醒一点VASP是商业软件如果你用的是自己的license或机器上的现有版本记得确认版本号。不同版本对某些功能比如IBRION7/8的DFPT支持不完全一致VASP 5.4和6.x之间的默认行为也有差异。3. 有限位移法算声子谱从结构优化到后处理的完整实操3.1 第一步高精度结构优化是声子谱的命根子声子谱的前提是结构必须处于受力平衡的状态。如果原子上有残余力位移之后的力响应就会叠加一个“初始漂移”直接导致力常数包含假信号反映到声子谱上就是虚频或频率偏移。所以结构优化的收敛标准必须远比普通电子结构计算严格。我的习惯是分两步走。第一步常规优化INCAR (第一步优化) PREC Accurate EDIFF 1E-8 EDIFFG -0.02 IBRION 1 ISIF 3 ISMEAR 0 SIGMA 0.05 ENCUT 600其中EDIFFG取负值表示力收敛单位是eV/埃-0.02 eV/埃是比较常规的精度。等这一步跑完接着用CONTCAR替换POSCAR把力的收敛标准提到-1E-6或更严INCAR (第二步精细优化) EDIFF 1E-8 EDIFFG -1E-6 IBRION 1 ISIF 3有人会问ISIF3同时优化晶胞参数和原子位置会不会和声子谱要求的“固定晶格”矛盾其实不矛盾。声子谱的位移是在优化后的平衡结构上做微扰晶格常数当然要用平衡值。先把晶格和原子优化到位后续生成位移构型时晶格矢量保持不变只动原子位置。优化结束后检查OUTCAR里的forces输出。如果最大力分量已经低于1E-6 eV/埃量级就可以进入下一步。如果力始终降不下去不要强行继续先排查是不是对称性设置、k点密度或ENCUT的问题。3.2 第二步用phonopy生成位移构型phonopy是目前主流的声子后处理工具支持VASP的接口。安装很简单pip或conda都能装conda install -c conda-forge phonopy这里提一句phonopy的分析基于群论对称性所以同一套位移构型里会考虑空间群的等效原子只会对不等价的原子方向生成位移构型能省下大量计算量。实际操作时先准备一个POSCAR。这个POSCAR最好来自优化后的CONTCAR并且保持晶格矢量的完整信息。接着设置超胞大小phonopy -d --dim3 3 3 --paF POSCAR--dim3 3 3表示在三个晶矢方向各扩3倍也就是27倍超胞。--paF是把超胞转换到原胞基矢可选但建议加上后面画图时路径还原更方便。phonopy运行完会生成一个SPOSCAR超胞POSCAR以及一堆disp-xxx目录xxx是序号。每个disp目录代表一个位移构型里面包含了被移动原子的坐标。这时你把原始POSCAR复制成每个disp目录里的POSCAR然后配置好INCAR、KPOINTS、POTCAR就可以进入下一步了。有一个细节位移大小。phonopy默认位移是0.01埃对大多数体系是合适的。如果之后发现声子谱对位移大小敏感比如换0.005和0.02结果差异大说明数值噪声过大或体系具有强非谐性要按情况调整。3.3 第三步对每个位移构型算受力这一步是最机械的但也是最能拉开质量差距的地方。每个disp目录里跑一次性静态计算INCAR设置如下INCAR (静态力计算) PREC Accurate EDIFF 1E-8 IBRION -1 NSW 0 ISMEAR 0 SIGMA 0.05注意几个关键点IBRION -1或NSW 0代表不进行离子弛豫只算电子自洽和原子受力。EDIFF 1E-8这是能量收敛标准力是通过Hellmann-Feynman定理得到的电子密度收敛得越准力越可靠。别省这一步。金属体系要注意ISMEAR。半导体或绝缘体用ISMEAR 0高斯展宽配合小的SIGMA没问题金属体系建议用ISMEAR 1Methfessel-Paxton或ISMEAR -5tetrahedron with Blöchl修正。其中-5不适合用于力的计算实际上是可以用-5的但对于金属-5在部分情况下会引入不连续稳妥起见还是用MP方法。k点密度上因为用的是超胞倒空间尺寸已经缩小k点数可以不用太大。比如原胞时是8×8×8扩成3×3×3超胞后k点用4×4×4一般就够了。具体还是测试一下收敛。POTCAR按元素顺序拼接好。然后每个disp目录提交VASP计算。为了监控进度我一般写一个循环脚本#!/bin/bash for dir in disp-*; do cd $dir mpirun -np 16 vasp_std run.log 21 cd .. done注意如果你有SLURM或PBS调度系统改成对应的任务提交脚本即可。这一步最花时间。27倍超胞加上足够的k点Si这种小体系几分钟到十几分钟就能跑完一个位移构型大一点的金属或合金一个构型算几个小时也很正常。3.4 第四步收集力常数并还原声子谱所有disp目录跑完后在父目录就是SPOSCAR所在的目录执行phonopy -f disp-*/vasprun.xmlphonopy会读取每个位移构型里的应变-力响应关系结合位移量拟合出力常数矩阵生成FORCE_CONSTANTS文件。这一步如果报错最常见的原因是某个disp目录的vasprun.xml找不到或计算中断检查后重跑对应构型即可。然后生成力常数phonopy --fc这一步会读取FORCE_CONSTANTS并构建动力学矩阵。在画声子谱之前需要先定义高对称点路径。phonopy里通过band.conf文件指定ATOM_NAME Si DIM 3 3 3 BAND G X W K G L U W L K U X BAND_POINTS 200注意BAND路径里的高对称点写法phonopy支持G或Gamma、X、K等常见标记具体按不同布拉伐格子会有差异。用之前可以先跑一下phonopy --symmetry检查原胞的对称性标签。然后画图phonopy -p band.conf它会输出一张声子色散关系图同时把数据存到band.yaml里。如果你还想要声子态密度再准备一个dos.confATOM_NAME Si DIM 3 3 3 MESH 50 50 50 DOS 1.0执行phonopy -p dos.conf即可。态密度配合色散曲线基本就是一篇声子计算的标准输出内容了。4. 关键参数对照直接抄作业的配置方案4.1 INCAR参数速查表下面这张表是我做有限位移法声子谱计算时常用的参数组合按体系类型区分参数半导体/绝缘体金属备注PRECAccurateAccurate高精度起步EDIFF1E-81E-8能量收敛宁严勿松EDIFFG优化时-1E-8-1E-8声子谱前建议再降一档IBRION静态力计算-1-1只算力不弛豫ISIF优化时33同时优化晶格和原子ISMEAR01金属用MP方法更稳SIGMA0.050.1~0.2金属适当加大避免电子步振荡ENCUT1.3×ENMAX1.3×ENMAX从POTCAR读取ENMAX留30%余量LREAL.FALSE..FALSE.声子计算必须关掉LREAL避免受力近似LWAVE.FALSE..FALSE.静态计算可关掉波函数输出省磁盘LCHARG.FALSE..FALSE.不需要自洽电荷密度文件时关闭这里面有两条特别值得强调。第一LREAL必须设为.FALSE.。LREAL是实空间投影近似主要用于大体系的电子结构计算提升速度但按我的记忆实空间近似下的局域投影会对原子受力引入误差声子这类对力极其敏感的物理量绝不能省这个精度。第二SIGMA在金属体系里要仔细测太大或太小都会造成电子步收敛振动间接恶化受力的数值精度。4.2 超胞大小的经验法则与试算策略超胞大小没有绝对标准完全由力常数的衰减范围决定。按我的经验可以按以下步骤来试先做一个中等超胞比如2×2×2或3×3×3跑完声子谱之后观察Gamma点附近的声学支和力常数在实空间中的衰减行为phonopy的--fc配合FORCE_CONSTANTS文件可以查看。如果在超胞边界处力常数仍大于某个量级比如10^-3 eV/埃²量级说明超胞不够大需要扩大到4×4×4甚至更大。体系类型上的经验值大概是硅、砷化镓、氧化物绝缘体3×3×3到4×4×4过渡金属及其合金4×4×4起步有时需要5×5×5层状材料垂直于平面方向至少2倍比如3×3×1或4×4×2金属体系要特别小心因为费米面的影响力常数在实空间中呈现振荡衰减衰减很慢。如果发现用5×5×5和4×4×4的结果差异仍然明显那基本就是金属长程力常数效应此时要么继续加大超胞要么考虑DFPT方案绕开实空间截断问题。5. 常见的坑虚频、力不收敛和phonopy报错5.1 虚频排查别急着改位移大小虚频是声子谱里最常见的“事故现场”。我看到很多人一出现虚频第一反应是去改位移或增大超胞但虚频的根源往往更基础。按优先级排查第一是结构有没有真正优化到受力为零。很多人第一步优化用的EDIFFG-0.02然后直接拿来算声子这是很危险的。特别是那些存在软模低频振动模式的体系残余力虽然只有0.01 eV/埃量级也足以让软模变成虚频。这时候把EDIFFG降到-1E-6甚至-1E-8重新优化虚频可能自己就消失了。第二是超胞够不够大。如果虚频集中在Gamma点附近而力常数在超胞边界还没衰减干净大概率是超胞尺寸的伪周期效应。第三是位移大小是否合理。位移太小数值噪声对力的影响会被放大位移太大高次非谐项污染力常数。如果这两个方向都调整过仍无改善要考虑是否体系本身确实不稳定。有一个技巧用声子态密度配合虚频出现的频率范围判断虚频涉及的原子种类和方向。结合实空间结构往往能判断出是某个原子在某个方向上的力常数偏弱对应着潜在的相变或结构失稳。这不是bug而是物理。5.2 电子步不收敛金属体系的SIGMA和混合参数在静态力计算里最让人抓狂的报错是电子自洽不收敛。现象一般是SCF循环里总能一直在某个值附近振荡或者根本降不下去。金属体系里常见原因是SIGMA设置不当。SIGMA太小布里渊区积分在费米面附近取样不足电子占据函数跳动导致总能振荡SIGMA太大电子展宽过度总能的误差变大。一般金属体系从SIGMA0.1开始试不行就0.2。同时检查ISMEAR是否用了适合金属的MP方法ISMEAR1。如果SIGMA没问题再检查混合参数。VASP里默认的混合方式对大多数体系够用但有些金属或含d/f电子的体系收敛很慢可以在INCAR里加AMIX 0.2 BMIX 0.0001 AMIX_MAG 0.2 BMIX_MAG 0.0001通过减少电荷密度混合比例来稳定SCF过程。也有人用ALGO VeryFast或ALGO Normal切换但我更推荐先检查k点密度和SIGMA动混合参数是最后手段。5.3 phonopy后处理时的几个经典报错phonopy -f disp-*/vasprun.xml这一步最常见的报错是某个目录里找不到vasprun.xml或者vasprun.xml里缺受力信息。前者多半是VASP没有正常结束后者可能是VASP版本和phonopy的兼容问题。比如VASP 6.x默认输出格式基本兼容但老版本VASP 5.x配合新phonopy偶尔会出现解析失败。VASP版本如果太老建议至少升级到5.4.4。还有一类报错和对称性有关。如果初始POSCAR的对称性标注有问题或者用了--pa参数后原胞选错生成的位移构型可能不符合空间群约束导致最后的力常数矩阵维度对不上。遇到这种情况返回去重新检查POSCAR的对称性和phonopy的对称性识别结果。画图时如果BAND路径不对通常是因为选取的高对称点不在该布拉伐格子的倒空间路径上。phonopy官网有不同晶系的高对称点路径图直接对着查。5.4 金属体系特有的长程振荡问题这里单独提醒一句金属体系里力常数的实空间衰减非常慢这会带来一个有点反直觉的后果——你的超胞可能已经大得离谱了结果仍然和更大超胞明显不同。遇到这种情况有两条路。一条是加大超胞硬算但成本急剧上升另一条是改用DFPTIBRION7/8因为DFPT在倒空间处理天然不受实空间截断半径限制。如果体系不大DFPT反而比有限位移法在金属环境下更干净。另一个和金属相关的点是费米面处理对力的影响。金属的布里渊区积分需要展宽但不同的SIGMA会影响受力的数值结果进而影响声子频率。建议至少用两个不同的SIGMA值测试同一套位移构型确认声子谱结果对SIGMA不敏感。5.5 磁场、自旋极化和声子计算的额外注意事项如果你的体系带磁性声子计算要更谨慎。自旋极化的计算里受力和磁矩的耦合会让问题更复杂。第一步优化必须在磁性自由度上充分收敛否则声子虚频可能纯粹来自磁性构型没有平衡。对于反铁磁或非共线磁性体系VASP声子计算的可靠性和成本都成倍上升建议先做无磁或固定磁矩的测试确认非磁性部分的力学性质是合理的再打开自旋极化。声子谱计算的结果强烈依赖你提供结构的稳定性。如果一个结构本身处于“力学不稳定”状态哪怕流程完全正确声子谱也会给出虚频——这未必是你算错了也可能你这个结构本来就该发生结构相变。判断的标准就是虚频出现在哪个高对称点、哪个原子方向对应着什么样的晶格畸变模式这就是后续研究的起点。6. 后处理进阶声子态密度、热力学性质和不稳定模式分析6.1 从声子谱到自由能零点能、熵和热容声子谱的作用不止是判断稳定性。有了完整的声子态密度就可以进一步得到晶格对热力学势的贡献。理想气体谐振子近似下声子频率直接决定了零点能、晶格熵和定容热容随温度的变化。phonopy自带热力学量计算phonopy -t --mesh50 thermo.conf其中thermo.conf需要指定超胞尺寸和温度范围ATOM_NAME Si DIM 3 3 3 MESH 50 50 50 T_MIN 0 T_MAX 1000 T_STEP 10它会输出自由能、熵、热容随温度的变化。虽然这是谐波近似的结果在高温或强非谐体系里误差会变大但对绝大多数晶体材料已经是相变、同素异构体稳定性对比时的有力参考。6.2 虚频的后续分析找到不稳定模式对应的原子运动当声子谱出现虚频后不要急着删参数查流程而是先分析这个虚频模式对应的原子运动方式。phonopy可以输出虚频处本征矢量的实空间表示phonopy -p band.conf --irreps或者直接看band.yaml里对应虚频频率处的特征向量在VESTA里可视化。你会看到原子沿着特定方向左右摆动这往往就是结构相变的软模对应的位移模式正是从母相过渡到子相的关键路径。这种情况下你不仅没有算错反而发现了更有意思的物理。我自己有一次算一个层状氧化物就在Gamma点出现了一个虚频本来想放弃但顺手画了虚频模式位移发现正好是层间剪切模式的A1g振动顺藤摸瓜发现该结构确实有一个更稳定的层错构型。后来认真看文献这个体系的层间滑移相变恰好是当时研究的热点。所以说虚频不一定是坑也可能是金矿。根据我个人经验做声子谱计算心态和手艺同样重要。算出一套没有虚频的声子谱不代表你行了算出一堆虚频还能冷静地把原因抠出来才算是真正入了门。建议新手拿到题目先花半天时间把一个简单体系Si、MgO这类从优化到出图完整跑一遍感受力常数、超胞、位移这些变量对结果的实际影响再上自己研究的体系会省掉很多弯路。这套流程吃透之后往后无论换哪个材料体系你都会知道该在什么地方较真、什么地方可以偷懒。
VASP声子谱计算实战：从物理图像到虚频排查全流程指南
做材料计算的人基本都会走到声子谱这一步。我在VASP声子谱计算这条路上踩过的坑比想象中多得多——虚频满天飞、力不收敛、超胞选小了力常数还没衰减到零、装完VASP一跑就报错。这篇笔记不打算写成教科书式的流程讲解而是把我自己的实操经验按完整链路整理出来从物理图像、方法选型到Ubuntu环境准备、有限位移法完整实操再到参数对照和问题排查希望给准备跑声子谱的朋友一条能直接走通的路。1. 声子谱到底在算什么先把物理图像建起来1.1 从晶格振动说起声子不是玄学是原子的集体舞晶格里的原子并不是钉死在平衡位置上的温度高了会振动零点能也让它不可能完全静止。原子偏离平衡位置时周围原子会把它拉回来这就形成了振动模式。大量原子耦合在一起振动就不再是孤立的单个原子运动而是整块晶体里周期性排列的集体行为。把这种集体振动量子化就得到声子phonon这个概念。声子谱就是以波矢为横轴、振动频率为纵轴画出来的色散关系。它告诉我们晶格里不同波长的振动模式各自拥有多高的频率以及这些模式是往哪个方向传播的。用VASP算声子谱本质上就是通过第一性原理把原子间的受力关系算准然后换算出晶格振动的频率。声子谱里有两个关键区域必须会看。一个是声学支在Gamma点波矢为零附近的行为——对于非极性材料声学支频率在Gamma点应当趋近于零另一个是是否存在虚频负频率。所谓虚频其实就是频率平方为负从数值上表现为画在零轴下方的曲线。虚频的出现通常意味着结构处于鞍点或亚稳态原子有自发向更低能量构型移动的趋势。判断一个结构是否动力学稳定看虚频基本是一眼定生死。1.2 两条技术路线有限位移法和DFPT怎么选VASP算声子谱主流有两条路。第一条是有限位移法也叫冻声子法。做法很直观在优化的超胞里人为地给某个原子一个很小的位移通常0.01埃然后计算所有原子受到的Hellmann-Feynman力。力对位移求偏导就得到力常数矩阵再通过傅里叶变换得到动力学矩阵最终解出声子频率。因为位移是显式的所以物理图像非常清楚实现也简单配合phonopy这类工具非常顺手。第二条是DFPT也就是密度泛函微扰理论在VASP里通过IBRION7或8调用。它不是在实空间里给一个有限位移而是在倒空间自洽地求解电荷密度对原子位移的线性响应。好处是不需要超胞也能算对极性材料还能很方便地引入非解析项修正LO-TO分裂缺点是内存占用量大对参数的敏感性也更高。我自己做体系时选型逻辑大概是这样的对比维度有限位移法DFPT物理直观性高直接反映力-位移关系低黑箱子感强是否需要超胞需要且越大越稳妥不需要大超胞内存开销低高尤其对大体系极性材料处理需要额外算Born有效电荷和介电常数天然支持配合NAC修正上手难度低中等偏上所以如果你是新手上路或者体系是半导体、绝缘体、简单金属我建议先用有限位移法。等你把流程跑通了再根据体系特性考虑要不要切换到DFPT。对于极性体系比如很多钙钛矿氧化物DFPT或有限位移法NAC修正是必须的否则Gamma点光学支会有明显误差。1.3 为什么非要超胞力常数的实空间局域性有限位移法里绕不开一个概念超胞。原胞太小原子之间的距离近周期性镜像之间会互相干扰算出来的力常数包含假信号。超胞的作用就是把原子的周期镜像推远让真实的力常数在实空间内衰减到可以忽略的程度。打个比方你在操场上推一个人他倒地时会牵动身边的人。如果操场太小四周还是这个人的镜像传递的扰动就会叠加出不真实的效果。把操场扩到足够大扰动传到边界时已经衰减干净测到的才是真实的响应。超胞大小怎么选核心判断标准只有一个力常数是否已经在超胞范围内衰减到接近零。一般来说半导体和绝缘体的力常数衰减得较快3×3×3或4×4×4的超胞常常够用金属体系因为费米面的存在力常数有长程振荡往往需要更大的超胞。后面我会专门讲这个怎么测试。2. 环境准备Ubuntu下把VASP跑起来的几个关键点2.1 编译安装VASP时最容易被卡住的三件事热搜词里有“ubuntu装vasp”看来这是很多人的第一道坎。其实VASP本身编译难度不算高但配套环境不弄好后患无穷。第一件事是编译器。VASP对Intel编译器ifort的兼容性最好性能也最好。建议直接装Intel oneAPI套件里面同时包含了ifort、MKL数学库和MPI库。如果你用gfortran编译时要注意某些算法部分对gfortran的支持不如ifort稳定性能也会有差距。装好之后务必检查一下环境变量是否生效source /opt/intel/oneapi/setvars.sh which ifort which mpiifort第二件事是数学库。VASP大量调用BLAS、LAPACK、FFTW直接用MKL最省心。很多编译失败案例都是因为MKL路径没写对或者库的顺序不对。VASP编译时用的是makefile.include文件需要在源码根目录下把这个文件放好。Intel平台可以参考安装包自带的arch/makefile.include.linux_intel模板然后根据实际路径修改。第三件事是MPI。VASP并行计算依赖MPI用Intel MPI和ifort搭配最稳。如果你机器上同时装了OpenMPI和Intel MPI很容易因为mpif90指向不对导致编译错乱。建议在编译前明确mpif90 --version确保这个命令指向你期望的MPI版本。2.2 编译完了怎么验证先算个小体系VASP编译完成后不要急着直接跑声子谱。先用一个简单体系比如金刚石结构Si做一次自洽计算和能带计算验证编译结果的正确性。你不需要自己准备测试文件VASP官方测试集里就有现成的例子。跑通之后再检查OUTCAR里的总能量是否和已知结果吻合比如Si的基态能量在LDA或PBE下都有公认值。这一步的意义是把你后续排查问题的范围缩小。如果连Si都跑不对那问题大概率出在编译或环境上如果Si没问题声子谱算坏了就可以放心去查计算参数和物理设置。另外提醒一点VASP是商业软件如果你用的是自己的license或机器上的现有版本记得确认版本号。不同版本对某些功能比如IBRION7/8的DFPT支持不完全一致VASP 5.4和6.x之间的默认行为也有差异。3. 有限位移法算声子谱从结构优化到后处理的完整实操3.1 第一步高精度结构优化是声子谱的命根子声子谱的前提是结构必须处于受力平衡的状态。如果原子上有残余力位移之后的力响应就会叠加一个“初始漂移”直接导致力常数包含假信号反映到声子谱上就是虚频或频率偏移。所以结构优化的收敛标准必须远比普通电子结构计算严格。我的习惯是分两步走。第一步常规优化INCAR (第一步优化) PREC Accurate EDIFF 1E-8 EDIFFG -0.02 IBRION 1 ISIF 3 ISMEAR 0 SIGMA 0.05 ENCUT 600其中EDIFFG取负值表示力收敛单位是eV/埃-0.02 eV/埃是比较常规的精度。等这一步跑完接着用CONTCAR替换POSCAR把力的收敛标准提到-1E-6或更严INCAR (第二步精细优化) EDIFF 1E-8 EDIFFG -1E-6 IBRION 1 ISIF 3有人会问ISIF3同时优化晶胞参数和原子位置会不会和声子谱要求的“固定晶格”矛盾其实不矛盾。声子谱的位移是在优化后的平衡结构上做微扰晶格常数当然要用平衡值。先把晶格和原子优化到位后续生成位移构型时晶格矢量保持不变只动原子位置。优化结束后检查OUTCAR里的forces输出。如果最大力分量已经低于1E-6 eV/埃量级就可以进入下一步。如果力始终降不下去不要强行继续先排查是不是对称性设置、k点密度或ENCUT的问题。3.2 第二步用phonopy生成位移构型phonopy是目前主流的声子后处理工具支持VASP的接口。安装很简单pip或conda都能装conda install -c conda-forge phonopy这里提一句phonopy的分析基于群论对称性所以同一套位移构型里会考虑空间群的等效原子只会对不等价的原子方向生成位移构型能省下大量计算量。实际操作时先准备一个POSCAR。这个POSCAR最好来自优化后的CONTCAR并且保持晶格矢量的完整信息。接着设置超胞大小phonopy -d --dim3 3 3 --paF POSCAR--dim3 3 3表示在三个晶矢方向各扩3倍也就是27倍超胞。--paF是把超胞转换到原胞基矢可选但建议加上后面画图时路径还原更方便。phonopy运行完会生成一个SPOSCAR超胞POSCAR以及一堆disp-xxx目录xxx是序号。每个disp目录代表一个位移构型里面包含了被移动原子的坐标。这时你把原始POSCAR复制成每个disp目录里的POSCAR然后配置好INCAR、KPOINTS、POTCAR就可以进入下一步了。有一个细节位移大小。phonopy默认位移是0.01埃对大多数体系是合适的。如果之后发现声子谱对位移大小敏感比如换0.005和0.02结果差异大说明数值噪声过大或体系具有强非谐性要按情况调整。3.3 第三步对每个位移构型算受力这一步是最机械的但也是最能拉开质量差距的地方。每个disp目录里跑一次性静态计算INCAR设置如下INCAR (静态力计算) PREC Accurate EDIFF 1E-8 IBRION -1 NSW 0 ISMEAR 0 SIGMA 0.05注意几个关键点IBRION -1或NSW 0代表不进行离子弛豫只算电子自洽和原子受力。EDIFF 1E-8这是能量收敛标准力是通过Hellmann-Feynman定理得到的电子密度收敛得越准力越可靠。别省这一步。金属体系要注意ISMEAR。半导体或绝缘体用ISMEAR 0高斯展宽配合小的SIGMA没问题金属体系建议用ISMEAR 1Methfessel-Paxton或ISMEAR -5tetrahedron with Blöchl修正。其中-5不适合用于力的计算实际上是可以用-5的但对于金属-5在部分情况下会引入不连续稳妥起见还是用MP方法。k点密度上因为用的是超胞倒空间尺寸已经缩小k点数可以不用太大。比如原胞时是8×8×8扩成3×3×3超胞后k点用4×4×4一般就够了。具体还是测试一下收敛。POTCAR按元素顺序拼接好。然后每个disp目录提交VASP计算。为了监控进度我一般写一个循环脚本#!/bin/bash for dir in disp-*; do cd $dir mpirun -np 16 vasp_std run.log 21 cd .. done注意如果你有SLURM或PBS调度系统改成对应的任务提交脚本即可。这一步最花时间。27倍超胞加上足够的k点Si这种小体系几分钟到十几分钟就能跑完一个位移构型大一点的金属或合金一个构型算几个小时也很正常。3.4 第四步收集力常数并还原声子谱所有disp目录跑完后在父目录就是SPOSCAR所在的目录执行phonopy -f disp-*/vasprun.xmlphonopy会读取每个位移构型里的应变-力响应关系结合位移量拟合出力常数矩阵生成FORCE_CONSTANTS文件。这一步如果报错最常见的原因是某个disp目录的vasprun.xml找不到或计算中断检查后重跑对应构型即可。然后生成力常数phonopy --fc这一步会读取FORCE_CONSTANTS并构建动力学矩阵。在画声子谱之前需要先定义高对称点路径。phonopy里通过band.conf文件指定ATOM_NAME Si DIM 3 3 3 BAND G X W K G L U W L K U X BAND_POINTS 200注意BAND路径里的高对称点写法phonopy支持G或Gamma、X、K等常见标记具体按不同布拉伐格子会有差异。用之前可以先跑一下phonopy --symmetry检查原胞的对称性标签。然后画图phonopy -p band.conf它会输出一张声子色散关系图同时把数据存到band.yaml里。如果你还想要声子态密度再准备一个dos.confATOM_NAME Si DIM 3 3 3 MESH 50 50 50 DOS 1.0执行phonopy -p dos.conf即可。态密度配合色散曲线基本就是一篇声子计算的标准输出内容了。4. 关键参数对照直接抄作业的配置方案4.1 INCAR参数速查表下面这张表是我做有限位移法声子谱计算时常用的参数组合按体系类型区分参数半导体/绝缘体金属备注PRECAccurateAccurate高精度起步EDIFF1E-81E-8能量收敛宁严勿松EDIFFG优化时-1E-8-1E-8声子谱前建议再降一档IBRION静态力计算-1-1只算力不弛豫ISIF优化时33同时优化晶格和原子ISMEAR01金属用MP方法更稳SIGMA0.050.1~0.2金属适当加大避免电子步振荡ENCUT1.3×ENMAX1.3×ENMAX从POTCAR读取ENMAX留30%余量LREAL.FALSE..FALSE.声子计算必须关掉LREAL避免受力近似LWAVE.FALSE..FALSE.静态计算可关掉波函数输出省磁盘LCHARG.FALSE..FALSE.不需要自洽电荷密度文件时关闭这里面有两条特别值得强调。第一LREAL必须设为.FALSE.。LREAL是实空间投影近似主要用于大体系的电子结构计算提升速度但按我的记忆实空间近似下的局域投影会对原子受力引入误差声子这类对力极其敏感的物理量绝不能省这个精度。第二SIGMA在金属体系里要仔细测太大或太小都会造成电子步收敛振动间接恶化受力的数值精度。4.2 超胞大小的经验法则与试算策略超胞大小没有绝对标准完全由力常数的衰减范围决定。按我的经验可以按以下步骤来试先做一个中等超胞比如2×2×2或3×3×3跑完声子谱之后观察Gamma点附近的声学支和力常数在实空间中的衰减行为phonopy的--fc配合FORCE_CONSTANTS文件可以查看。如果在超胞边界处力常数仍大于某个量级比如10^-3 eV/埃²量级说明超胞不够大需要扩大到4×4×4甚至更大。体系类型上的经验值大概是硅、砷化镓、氧化物绝缘体3×3×3到4×4×4过渡金属及其合金4×4×4起步有时需要5×5×5层状材料垂直于平面方向至少2倍比如3×3×1或4×4×2金属体系要特别小心因为费米面的影响力常数在实空间中呈现振荡衰减衰减很慢。如果发现用5×5×5和4×4×4的结果差异仍然明显那基本就是金属长程力常数效应此时要么继续加大超胞要么考虑DFPT方案绕开实空间截断问题。5. 常见的坑虚频、力不收敛和phonopy报错5.1 虚频排查别急着改位移大小虚频是声子谱里最常见的“事故现场”。我看到很多人一出现虚频第一反应是去改位移或增大超胞但虚频的根源往往更基础。按优先级排查第一是结构有没有真正优化到受力为零。很多人第一步优化用的EDIFFG-0.02然后直接拿来算声子这是很危险的。特别是那些存在软模低频振动模式的体系残余力虽然只有0.01 eV/埃量级也足以让软模变成虚频。这时候把EDIFFG降到-1E-6甚至-1E-8重新优化虚频可能自己就消失了。第二是超胞够不够大。如果虚频集中在Gamma点附近而力常数在超胞边界还没衰减干净大概率是超胞尺寸的伪周期效应。第三是位移大小是否合理。位移太小数值噪声对力的影响会被放大位移太大高次非谐项污染力常数。如果这两个方向都调整过仍无改善要考虑是否体系本身确实不稳定。有一个技巧用声子态密度配合虚频出现的频率范围判断虚频涉及的原子种类和方向。结合实空间结构往往能判断出是某个原子在某个方向上的力常数偏弱对应着潜在的相变或结构失稳。这不是bug而是物理。5.2 电子步不收敛金属体系的SIGMA和混合参数在静态力计算里最让人抓狂的报错是电子自洽不收敛。现象一般是SCF循环里总能一直在某个值附近振荡或者根本降不下去。金属体系里常见原因是SIGMA设置不当。SIGMA太小布里渊区积分在费米面附近取样不足电子占据函数跳动导致总能振荡SIGMA太大电子展宽过度总能的误差变大。一般金属体系从SIGMA0.1开始试不行就0.2。同时检查ISMEAR是否用了适合金属的MP方法ISMEAR1。如果SIGMA没问题再检查混合参数。VASP里默认的混合方式对大多数体系够用但有些金属或含d/f电子的体系收敛很慢可以在INCAR里加AMIX 0.2 BMIX 0.0001 AMIX_MAG 0.2 BMIX_MAG 0.0001通过减少电荷密度混合比例来稳定SCF过程。也有人用ALGO VeryFast或ALGO Normal切换但我更推荐先检查k点密度和SIGMA动混合参数是最后手段。5.3 phonopy后处理时的几个经典报错phonopy -f disp-*/vasprun.xml这一步最常见的报错是某个目录里找不到vasprun.xml或者vasprun.xml里缺受力信息。前者多半是VASP没有正常结束后者可能是VASP版本和phonopy的兼容问题。比如VASP 6.x默认输出格式基本兼容但老版本VASP 5.x配合新phonopy偶尔会出现解析失败。VASP版本如果太老建议至少升级到5.4.4。还有一类报错和对称性有关。如果初始POSCAR的对称性标注有问题或者用了--pa参数后原胞选错生成的位移构型可能不符合空间群约束导致最后的力常数矩阵维度对不上。遇到这种情况返回去重新检查POSCAR的对称性和phonopy的对称性识别结果。画图时如果BAND路径不对通常是因为选取的高对称点不在该布拉伐格子的倒空间路径上。phonopy官网有不同晶系的高对称点路径图直接对着查。5.4 金属体系特有的长程振荡问题这里单独提醒一句金属体系里力常数的实空间衰减非常慢这会带来一个有点反直觉的后果——你的超胞可能已经大得离谱了结果仍然和更大超胞明显不同。遇到这种情况有两条路。一条是加大超胞硬算但成本急剧上升另一条是改用DFPTIBRION7/8因为DFPT在倒空间处理天然不受实空间截断半径限制。如果体系不大DFPT反而比有限位移法在金属环境下更干净。另一个和金属相关的点是费米面处理对力的影响。金属的布里渊区积分需要展宽但不同的SIGMA会影响受力的数值结果进而影响声子频率。建议至少用两个不同的SIGMA值测试同一套位移构型确认声子谱结果对SIGMA不敏感。5.5 磁场、自旋极化和声子计算的额外注意事项如果你的体系带磁性声子计算要更谨慎。自旋极化的计算里受力和磁矩的耦合会让问题更复杂。第一步优化必须在磁性自由度上充分收敛否则声子虚频可能纯粹来自磁性构型没有平衡。对于反铁磁或非共线磁性体系VASP声子计算的可靠性和成本都成倍上升建议先做无磁或固定磁矩的测试确认非磁性部分的力学性质是合理的再打开自旋极化。声子谱计算的结果强烈依赖你提供结构的稳定性。如果一个结构本身处于“力学不稳定”状态哪怕流程完全正确声子谱也会给出虚频——这未必是你算错了也可能你这个结构本来就该发生结构相变。判断的标准就是虚频出现在哪个高对称点、哪个原子方向对应着什么样的晶格畸变模式这就是后续研究的起点。6. 后处理进阶声子态密度、热力学性质和不稳定模式分析6.1 从声子谱到自由能零点能、熵和热容声子谱的作用不止是判断稳定性。有了完整的声子态密度就可以进一步得到晶格对热力学势的贡献。理想气体谐振子近似下声子频率直接决定了零点能、晶格熵和定容热容随温度的变化。phonopy自带热力学量计算phonopy -t --mesh50 thermo.conf其中thermo.conf需要指定超胞尺寸和温度范围ATOM_NAME Si DIM 3 3 3 MESH 50 50 50 T_MIN 0 T_MAX 1000 T_STEP 10它会输出自由能、熵、热容随温度的变化。虽然这是谐波近似的结果在高温或强非谐体系里误差会变大但对绝大多数晶体材料已经是相变、同素异构体稳定性对比时的有力参考。6.2 虚频的后续分析找到不稳定模式对应的原子运动当声子谱出现虚频后不要急着删参数查流程而是先分析这个虚频模式对应的原子运动方式。phonopy可以输出虚频处本征矢量的实空间表示phonopy -p band.conf --irreps或者直接看band.yaml里对应虚频频率处的特征向量在VESTA里可视化。你会看到原子沿着特定方向左右摆动这往往就是结构相变的软模对应的位移模式正是从母相过渡到子相的关键路径。这种情况下你不仅没有算错反而发现了更有意思的物理。我自己有一次算一个层状氧化物就在Gamma点出现了一个虚频本来想放弃但顺手画了虚频模式位移发现正好是层间剪切模式的A1g振动顺藤摸瓜发现该结构确实有一个更稳定的层错构型。后来认真看文献这个体系的层间滑移相变恰好是当时研究的热点。所以说虚频不一定是坑也可能是金矿。根据我个人经验做声子谱计算心态和手艺同样重要。算出一套没有虚频的声子谱不代表你行了算出一堆虚频还能冷静地把原因抠出来才算是真正入了门。建议新手拿到题目先花半天时间把一个简单体系Si、MgO这类从优化到出图完整跑一遍感受力常数、超胞、位移这些变量对结果的实际影响再上自己研究的体系会省掉很多弯路。这套流程吃透之后往后无论换哪个材料体系你都会知道该在什么地方较真、什么地方可以偷懒。
