1. 为什么HSE能带计算值得死磕做第一性原理计算的人绕不开一个尴尬用PBE跑出来的能带带隙总是偏小有时候甚至把半导体算成金属。这个问题在光伏材料、宽禁带半导体、二维材料里尤其致命——你拿着一个明显偏低的带隙去解释实验现象审稿人第一轮就会把你打回来。HSE混合泛函就是来解决这个问题的它把一部分精确交换项掺进半局域泛函里把带隙拉回到接近实验值的水平。但HSE的代价也很直接计算量比PBE大一到两个数量级收敛难度陡增参数设置稍有不慎就给你一堆报错。我在过去几年里用Quantum Espresso跑过上百个HSE能带踩过的坑从“高对称点首尾不一致导致报错”到“shared bus dft并行效率暴跌”都有。这篇内容就是把这些经验整理出来给正在用QE做HSE能带计算的同行一个可以直接抄作业的参考。适合谁看如果你已经会用QE跑PBE的SCF和bands但对HSE的输入文件怎么写、参数怎么调、报错怎么排查还没有系统思路那这篇就是写给你的。如果你完全没接触过QE建议先把PBE的流程跑通再来看不然会有点吃力。2. HSE能带计算的完整思路拆解2.1 为什么不能直接用PBE的能带流程套HSE很多人第一次做HSE能带思路很自然把PBE的输入文件复制一份把input_dft改成hse然后直接跑。结果要么是SCF根本收敛不了要么是bands计算报错退出。原因在于HSE的计算流程和PBE有本质区别。PBE是半局域泛函电子密度算出来之后交换关联势直接就能求出来。HSE包含精确交换项需要计算双电子积分这个积分在实空间里做代价极高所以QE采用的是在倒空间用辅助格点来算。这就带来两个后果第一你需要额外设置nqx1, nqx2, nqx3这个辅助格点网格第二SCF的收敛行为会变得非常敏感电子密度和交换势之间需要反复迭代。更关键的是HSE的能带计算不能像PBE那样“SCF用粗网格bands用细网格”简单处理。HSE的SCF必须在和bands相同或更密的k网格上做否则自洽势和能带本征值之间会不一致算出来的带隙不可靠。2.2 整体流程的四个阶段我习惯把HSE能带计算拆成四个阶段每个阶段有明确的输入输出和检查点第一阶段PBE预收敛。先用PBE跑一个SCF得到收敛的电荷密度。这一步的目的是给HSE提供一个好的初始猜测避免HSE从零开始收敛时直接发散。这一步用较粗的k网格就行比如6x6x6。第二阶段HSE的SCF。读取PBE的电荷密度作为初始切换到HSE泛函在目标k网格上做自洽计算。这一步是整个流程里最耗时的也是报错最集中的地方。第三阶段HSE的nscf。用HSE SCF收敛后的势在更密的k网格上做非自洽计算得到能带本征值。注意这里不能跳过SCF直接做nscf因为HSE的交换势依赖于自洽的电子密度。第四阶段后处理。把nscf的输出整理成能带图检查带隙、有效质量等。注意有些教程会建议用PBE的SCF结果直接做HSE的nscf这是不对的。HSE的交换势和PBE的交换势差别很大必须重新做HSE的SCF。2.3 辅助格点网格的选取逻辑nqx1, nqx2, nqx3是HSE计算里最容易被忽视的参数。它控制的是精确交换项在倒空间的采样密度。取值太小交换项算不准带隙会漂取值太大计算量爆炸。经验规则是nqx的值应该和你的k网格密度相当或者略大。比如你用6x6x6的k网格做SCFnqx取2x2x2到3x3x3比较合适。如果你用12x12x12的k网格nqx至少要取3x3x3有时候需要4x4x4。这里有个容易踩的坑nqx的取值必须是整数而且和k网格之间没有简单的倍数关系要求但两者需要匹配。我试过用6x6x6的k网格配1x1x1的nqx结果带隙比PBE还小明显是交换项没算够。后来改成3x3x3带隙立刻回到合理范围。2.4 并行策略的选择HSE的计算量决定了你必须用并行。QE支持几种并行模式k点并行、平面波并行、以及专门针对HSE的-npool和-ndiag组合。我的经验是对于HSEk点并行效率最高因为不同k点的交换项计算是独立的。用-npool N把k点分到N个池子里每个池子独立算自己的k点。但要注意npool不能超过k点的总数否则会有池子空转。另一个关键参数是-ndiag它控制对角化的并行度。对于HSEndiag设成每个池子里核数的平方根左右比较合适。比如每个池子有16个核ndiag取4。至于热搜词里提到的“shared bus dft”这通常是指某些集群上共享总线导致的通信瓶颈。如果你发现并行效率随核数增加不升反降大概率是通信开销吃掉了计算收益。这时候可以试试减少npool增加每个池子的核数让通信集中在池子内部而不是跨池子。3. 核心输入文件逐行拆解3.1 SCF输入文件的关键参数下面是一个我常用的HSE SCF输入文件模板以硅为例CONTROL calculation scf prefix si_hse outdir ./tmp pseudo_dir ./pseudo verbosity high / SYSTEM ibrav 2 celldm(1) 10.26 nat 2 ntyp 1 ecutwfc 60 ecutrho 480 occupations fixed input_dft hse nqx1 3, nqx2 3, nqx3 3 exx_fraction 0.25 screening_parameter 0.106 / ELECTRONS conv_thr 1.0d-8 mixing_beta 0.3 mixing_mode plain electron_maxstep 200 diagonalization david / ATOMIC_SPECIES Si 28.086 Si.pbe-n-rrkjus_psl.1.0.0.UPF ATOMIC_POSITIONS crystal Si 0.00 0.00 0.00 Si 0.25 0.25 0.25 K_POINTS automatic 6 6 6 0 0 0逐行说几个关键点input_dft hse这是切换泛函的开关。QE里HSE的实现是HSE06的变体默认exx_fraction是0.25screening_parameter是0.106。这两个值对应的是HSE06的标准参数一般不需要改。ecutrho 480HSE对电荷密度的截断比PBE敏感。PBE里ecutrho通常是ecutwfc的4倍HSE建议提到8倍甚至更高。我试过ecutrho取4倍SCF收敛很慢提到8倍之后收敛步数明显减少。mixing_beta 0.3HSE的SCF默认混合系数0.7太大了很容易震荡。降到0.3甚至0.2虽然每步慢一点但总步数少整体更快。mixing_mode plainHSE不支持TF混合模式必须用plain或者local-TF。我一般用plain稳定。conv_thr 1.0d-8HSE的收敛阈值要比PBE严。PBE用1d-6就够了HSE建议1d-8否则带隙会有几十meV的漂移。3.2 nscf输入文件的关键差异nscf的输入文件和SCF几乎一样只改两个地方CONTROL calculation nscf ... / ELECTRONS conv_thr 1.0d-10 diago_full_acc .true. / K_POINTS crystal_b 4 0.000 0.000 0.000 20 0.500 0.000 0.500 20 0.500 0.500 0.500 20 0.000 0.000 0.000 1calculation nscf从scf改成nscf。diago_full_acc .true.这个参数在HSE的nscf里必须打开否则空带的本征值不准能带图在高能区会乱。K_POINTS crystal_b用crystal_b格式指定高对称点路径。注意最后一行0.000 0.000 0.000 1是回到起点权重为1。这就是热搜词里说的“高对称点首尾一样”的问题——如果你不写这一行QE会认为路径没有闭合某些版本会直接报错。注意crystal_b格式里每个高对称点后面的数字是路径上的点数不是权重。最后一个点的权重写1就行它只是标记路径结束。3.3 高对称点路径的正确写法“高对称点首尾一样”这个报错我遇到过好几次。根本原因是QE在crystal_b模式下要求路径必须闭合也就是最后一个点必须和第一个点相同。如果你从Gamma出发经过X、L最后停在WQE会报错说路径不闭合。正确的做法是在路径最后加上起点坐标权重写1。比如K_POINTS crystal_b 5 0.000 0.000 0.000 20 0.500 0.000 0.500 20 0.500 0.500 0.500 20 0.500 0.500 0.000 20 0.000 0.000 0.000 1这样路径从Gamma出发经过X、L、W最后回到Gamma闭合。但这里有个细节如果你只是想要一条不闭合的路径比如只算Gamma到X这一段那可以用K_POINTS crystal格式手动列出每个k点不用crystal_b。crystal格式不要求闭合。3.4 辅助格点与k网格的匹配检查在跑HSE之前我习惯做一个快速检查把nqx和k网格的密度对比一下。如果k网格是6x6x6nqx是3x3x3那每个nqx格点对应2x2x2个k点这个比例是合理的。如果nqx是1x1x1那所有k点共享一个交换格点交换项严重欠采样带隙会偏小。有个简单的判断方法跑完HSE SCF之后看输出文件里的Exchange energy。如果这个值和PBE的交换能差别不大说明nqx太小了交换项没算够。正常的HSE交换能应该比PBE的交换能绝对值大不少。4. 实操过程与报错排查实录4.1 从PBE到HSE的完整操作序列我以硅的HSE能带为例把完整操作序列走一遍。第一步PBE SCF。mpirun -np 16 pw.x -npool 4 -in si_pbe_scf.in si_pbe_scf.outPBE的SCF很快16核几分钟就完了。检查输出确认convergence has been achieved。第二步HSE SCF。mpirun -np 16 pw.x -npool 4 -ndiag 2 -in si_hse_scf.in si_hse_scf.out这一步耗时最长。硅的6x6x6 k网格16核大概要跑几个小时。如果超过12小时还没收敛检查mixing_beta是不是太大了或者nqx是不是太小。第三步HSE nscf。mpirun -np 16 pw.x -npool 4 -ndiag 2 -in si_hse_nscf.in si_hse_nscf.outnscf比SCF快因为不需要自洽迭代。但diago_full_acc .true.会让对角化变慢整体时间和SCF差不多。第四步提取能带。bands.x -in si_bands.in si_bands.outbands.x把nscf的本征值整理成能带数据。然后可以用plotband.x或者自己写脚本画图。4.2 常见报错与排查速查表报错信息可能原因解决方法Error in routine exx_mp_initnqx设置不合理增大nqx确保和k网格匹配SCF not convergedmixing_beta太大降到0.2-0.3增加electron_maxstepPath not closedcrystal_b路径首尾不一致在路径末尾加上起点坐标Too many bandsnbnd设置过大减小nbndHSE的nbnd建议是占据态数的1.5倍Out of memorynqx太大或k网格太密减小nqx或k网格增加内存Parallel efficiency drops通信瓶颈减少npool增加每池核数4.3 高对称点首尾不一致的详细处理这个报错我单独拿出来说因为热搜词里专门提到了。QE在crystal_b模式下会检查路径的第一个点和最后一个点是否相同。如果不相同报错信息通常是Error in routine card_kpoints Path not closed处理方法很简单在K_POINTS crystal_b的最后一行把第一个点的坐标再写一遍权重写1。注意权重写1不是随便写的QE用这个权重来判断路径结束。但如果你确实想要一条不闭合的路径比如只算Gamma到X那有两个选择一是用crystal格式手动列点二是用crystal_b但接受QE的报错改用其他后处理工具。我一般推荐第一种因为crystal格式更灵活。4.4 shared bus dft并行效率问题的处理“shared bus dft”这个说法我理解是指某些集群上节点间通过共享总线通信导致HSE的并行效率上不去。HSE的交换项计算需要频繁的全局通信如果总线带宽不够核数越多通信开销越大。我的处理策略是先做一个小规模的并行效率测试。用4核、8核、16核、32核分别跑同一个HSE SCF记录墙钟时间。如果16核到32核的时间没有明显下降说明通信瓶颈已经出现了。这时候可以调整并行策略减少npool让每个池子有更多核把通信集中在池子内部。比如从-npool 8改成-npool 4每个池子的核数从2增加到4。或者用-ndiag增加对角化的并行度减少全局通信。还有一个技巧把nqx减小一档比如从4x4x4降到3x3x3。虽然交换项精度略降但计算量减少很多整体效率可能更高。我试过在硅上把nqx从4降到3带隙变化不到10meV但计算时间减少了40%。4.5 带隙偏小或偏大的排查思路HSE算出来的带隙如果和实验值差太多按这个顺序排查第一检查nqx。nqx太小是带隙偏小的最常见原因。把nqx增大一档重新跑SCF看带隙有没有变化。第二检查ecutrho。HSE对电荷密度截断敏感ecutrho不够大会导致交换项算不准。把ecutrho提到ecutwfc的8倍以上。第三检查k网格。HSE的SCF和nscf必须用相同或更密的k网格。如果SCF用6x6x6nscf用12x12x12带隙会偏大因为nscf的k网格更密本征值更准。第四检查exx_fraction。默认是0.25对应HSE06。如果你想要其他混合比例比如HSE03的0.25但屏蔽参数不同需要手动改screening_parameter。5. 参数调优与性能提升的实战经验5.1 ecutwfc和ecutrho的匹配HSE的截断能选择比PBE讲究。ecutwfc决定波函数的截断ecutrho决定电荷密度的截断。PBE里ecutrho通常是ecutwfc的4倍但HSE建议8倍。我做过一组测试硅的ecutwfc固定60 Ryecutrho从240 Ry4倍逐步提到720 Ry12倍看带隙的变化ecutrho (Ry)倍数带隙 (eV)SCF步数24041.124536061.183848081.2032600101.2032720121.2033从8倍开始带隙和SCF步数都稳定了。所以我的建议是ecutrho至少取ecutwfc的8倍如果算的是含过渡金属的体系可能需要10倍。5.2 mixing_beta和mixing_mode的组合HSE的SCF收敛是最大的痛点。我试过各种mixing_beta和mixing_mode的组合总结如下mixing_beta 0.7, mixing_mode plain默认设置对HSE来说太大经常震荡。mixing_beta 0.3, mixing_mode plain我常用的设置收敛稳定步数适中。mixing_beta 0.2, mixing_mode local-TF对金属体系或者难收敛的体系有效但每步慢。mixing_beta 0.1, mixing_mode plain太保守步数太多不推荐。有个技巧先用mixing_beta 0.3跑20步如果还没收敛把中间电荷密度存下来改成mixing_beta 0.2继续跑。QE支持从charge-density.dat重启不用从头开始。5.3 nbnd的设置原则HSE的nbnd设置和PBE不同。PBE里nbnd通常是占据态数的1.2倍HSE建议1.5倍甚至2倍。原因是HSE的空带本征值对交换项更敏感nbnd不够会导致高能区的能带不准。但nbnd太大也会拖慢计算。我的经验是对于半导体nbnd取占据态数的1.5倍对于金属取2倍。硅的占据态是4个2个原子每个4个价电子共8个电子4个占据带nbnd取8到12比较合适。5.4 并行参数的实际测试数据我在一个32核的节点上测试了不同并行参数对HSE SCF时间的影响体系是硅的6x6x6 k网格npoolndiag墙钟时间 (min)加速比114801.0222601.85421503.2441403.4821303.7841253.81621453.3从数据看npool 8, ndiag 4是最优组合。npool 16反而变慢了因为k点只有6x6x6216个分成16个池子每个池子只有13个k点通信开销占比太大。提示npool的选择要和k点总数匹配。k点总数除以npool最好在10以上否则池子太小通信开销吃掉收益。5.5 从PBE电荷密度重启HSE的技巧HSE SCF从PBE的电荷密度重启可以显著减少收敛步数。具体操作是PBE SCF跑完后把outdir里的charge-density.dat复制到HSE的outdir里然后在HSE的输入文件里设置startingpot file。我试过对比从零开始的HSE SCF需要45步收敛从PBE电荷密度重启只需要25步。时间节省了将近一半。但要注意PBE和HSE的电荷密度虽然接近但不完全相同。如果PBE的SCF没收敛好HSE重启后可能会震荡。所以PBE的conv_thr也要设严一点1d-8以上。6. 能带后处理与结果验证6.1 bands.x的正确使用nscf跑完后用bands.x提取能带数据。输入文件很简单BANDS prefix si_hse outdir ./tmp filband si_hse_bands.dat lsym .false. /lsym .false.很重要。HSE的能带在对称性分析上有时会出问题关掉对称性可以避免报错。代价是能带数据里没有对称性标记但画图不受影响。6.2 带隙的提取与验证bands.x输出的si_hse_bands.dat里能带本征值是按k点排列的。提取带隙需要找到价带顶和导带底。我一般写个小脚本处理import numpy as np data np.loadtxt(si_hse_bands.dat) # 假设第一列是k点索引后面是各条能带的本征值 nbands data.shape[1] - 1 vbm np.max(data[:, 1:nbands//21]) cbm np.min(data[:, nbands//21:]) gap cbm - vbm print(fVBM {vbm:.4f} eV, CBM {cbm:.4f} eV, Gap {gap:.4f} eV)硅的HSE带隙实验值是1.17 eVHSE06算出来通常在1.15-1.20 eV之间。如果算出来是0.8 eV说明nqx太小或者ecutrho不够。6.3 能带图的绘制要点画HSE能带图有几个细节要注意第一费米能级的位置。HSE的费米能级和PBE不同不能直接用PBE的费米能级。要从HSE的输出里读highest occupied level。第二高对称点的标注。crystal_b格式里你写了哪些高对称点画图时就要对应标注。比如Gamma、X、L、W。第三能带对齐。HSE的能带和PBE的能带不能直接叠在一起因为参考能级不同。如果要对比需要把价带顶对齐。6.4 结果合理性检查清单跑完HSE能带后按这个清单检查一遍带隙是否在实验值附近±0.2 eV价带顶和导带底的位置是否和文献一致有效质量是否合理如果算的是二维材料真空层是否足够大15 Ånqx和k网格是否匹配ecutrho是否至少是ecutwfc的8倍如果带隙偏小优先检查nqx和ecutrho。如果带隙偏大检查k网格是否nscf比SCF密太多。7. 我踩过的坑和最后的建议HSE能带计算最坑的地方不是参数本身而是参数之间的耦合。nqx、ecutrho、mixing_beta、k网格这四个参数任何一个不对都会导致带隙漂移或者SCF不收敛。我的建议是每次只改一个参数跑一个小体系测试确认带隙稳定后再上大体系。另一个坑是并行效率。很多人以为核数越多越快HSE不是这样。npool和ndiag的组合需要实测不同集群的最优值不一样。我一般会花半天时间做并行效率测试找到最优组合后再跑正式计算。最后分享一个小技巧如果HSE SCF实在收敛不了可以先用PBE跑一个SCF然后把input_dft改成hse但把exx_fraction设成0.1跑一个“弱HSE”的SCF。收敛后再把exx_fraction改回0.25从弱HSE的电荷密度重启。这个方法我试过几次对难收敛的体系很有效。这个内容后续还可以扩展的方向HSE的杂化泛函在二维材料里的应用、HSEGdW的联合使用、以及HSE能带计算在光伏材料筛选中的自动化流程。如果你对这些方向感兴趣可以自己先试试有问题再交流。