干我们计算材料这行的估计都逃不过这么一遭用DFT把某个半导体的带隙算出来高高兴兴拿去组会汇报结果导师一句这个带隙和实验差了0.8 eV光学吸收谱的第一峰位置也不对瞬间把人打回原形。问题出在哪儿DFT的Kohn-Sham本征值本质上是辅助量不是准粒子激发能算带隙天然偏低要算光吸收谱还得考虑电子-空穴的库仑相互作用也就是激子效应。这时候VASP里的GWBSE组合就成了绕不开的解决方案。这套方法算下来带隙精度能拉到0.1 eV量级吸收谱的激子峰位置也能和实验对上。GW负责把DFT的单粒子能级修正成准粒子能级BSE负责在准粒子能级基础上解电子-空穴束缚态方程两者配合才算把一个材料的光电性质从定性对推进到定量对。这篇文章就围绕用VASP算GWBSE这条线从方法动机、参数设置、实操流程到踩坑记录把我攒下的经验一次性倒出来。适合刚接触激发态计算、准备给体系上GWBSE的研究生也适合被DFT结果搞到怀疑人生、想换方法的老手参考。1. 为什么DFT算不准带隙之后我会转向GWBSE先说清楚一个容易混淆的点DFT不是一个坏方法只是它的Kohn-Sham轨道能级在物理上并没有被严格定义成电子激发能。严格来说Kohn-Sham本征值只是拉格朗日乘子除了最高占据态在渐近极限下有些意义其他能级和真实的光电子谱并没有严格对应关系。所以算出来的带隙偏小不是误差而是方法本身的局限。那GW为什么能修因为GW通过自能算符Σ对单粒子格林函数做微扰修正把电子感受到的其他电子的交换关联效应重新整理了一遍。自能里面有G单粒子格林函数和W屏蔽库仑相互作用G负责描述电子-空穴对的传播W负责描述电子在介质中感受到的动态屏蔽相互作用。把这两个东西卷积在一起得到的就是准粒子色散关系准粒子带隙自然比KS带隙更接近实验。再来说BSE。GW只解决单粒子激发的问题也就是一个电子被激发到导带之后的能级位置。但光吸收是一个双粒子过程一个电子从价带跃迁到导带原来的价带位置留下一个空穴电子和空穴之间因为有库仑吸引可能束缚在一起这个束缚态就是激子。激子能级比GW算出的带边低一个结合能在吸收谱里表现为带边以下或带边附近的尖锐峰。这些信息单粒子GW给不了必须解BSE方程。所以一句话总结我的经验只关心带隙GW够用关心吸收谱形状、激子峰位置、光学响应必须上BSE。顺带说一句什么时候可以不上GWBSE如果你只是做结构优化、算形成能、看电荷转移这类基态性质DFT和DFTU完全够用。材料带隙不算太小、激子结合能又很弱比如很多三维无机半导体用G0W0算个带隙也够了。真正非BSE不可的是二维材料、有机半导体、钙钛矿这类激子效应显著、光学谱里有明显激子峰的体系。2. G0W0计算全流程从DFT基态到准粒子修正VASP里跑GW我个人的习惯是分两步走。第一步先把DFT基态做扎实第二步在同一个计算目录下换INCAR直接续算GW。也有人用一条INCAR让VASP自己从DFT无缝切到GW但那样一旦DFT部分没收敛后面全白算排查起来很痛苦。分步走的好处是每一段都能单独验证。2.1 先把DFT基态做扎实这一步没什么花活但有几个细节会影响后续GW质量结构先优化到位。GW非常贵你不可能用GW做结构优化。所以POSCAR里的原子位置必须是DFT优化后的结果力收敛标准我一般开到EDIFFG-0.01甚至-0.005比默认值严一个量级。KPOINTS用Gamma-centered的网格。这一点后面BSE还会强调。GW和BSE都要求K点网格包含Gamma点因为光学跃迁的联合态密度在Gamma点附近贡献最大。千万别图省事用Monkhorst-Pack生成偏移网格。ENCUT要测收敛。对于GW来说ENCUT不仅影响平面波基组还影响介电函数和自能的收敛。我做常见半导体比如Si、GaAs、MAPbI3这类ENCUT取到平面波基组默认截断能的1.3倍左右才开始稳定。INCAR里可以顺手加上PREC Accurate。DFT这一步的典型INCAR长这样SYSTEM dft ground state PREC Accurate ENCUT 400 ISMEAR 0 SIGMA 0.05 EDIFF 1E-6 EDIFFG -0.01 IBRION 2 ISIF 3 NSW 100 LORBIT 11 LCHARG .TRUE. LWAVE .TRUE.注意SIGMA用0.05是常规操作但绝缘体和半导体用ISMEAR0即高斯展宽没问题。如果是金属后面算GW会遇到很多麻烦因为费米面处的占据数在GW里怎么处理都不太舒服这也是为什么GWBSE主要应用场景是半导体和绝缘体。2.2 换INCAR进入GW计算DFT跑完CHGCAR和WAVECAR都在了这时候把INCAR替换成GW版本保持POSCAR、KPOINTS、POTCAR不变直接继续跑。我的G0W0 INCAR模板SYSTEM G0W0 PREC Accurate ENCUT 400 ISMEAR 0 SIGMA 0.05 EDIFF 1E-6 NELM 200 ALGO GW0 NBANDS 300 ENCUTGW 200 OMEGAMAX 50 NOMEGA 50 LSPECTRAL .TRUE. LMAXFOCKAE .TRUE. NKRED 1 ISYM 0 LOPTICS .TRUE.逐个说说关键参数ALGO GW0。VASP一共提供好几种GW级别最便宜的是单次G0W0往上还有EVGW0本征值自洽、scGW完全自洽。BSE的输入是需要一套准粒子能级和波函数G0W0的精度通常已经够用。我自己算过几个体系G0W0和EVGW0的带隙差一般在0.05 eV以内但EVGW0的计算量可能是G0W0的两到三倍性价比不高。NELM 200。GW循环的迭代上限其实G0W0的NELM主要影响自洽迭代过程中的微扰修正次数。放心绝大多数情况跑不满200步。NBANDS 300。这是决定GW精度的核心参数。GW里自能对未占据态的收敛非常慢需要把你体系里所有关心的导带底附近能级对应的未占据态都包进去。经验法则是NBANDS至少是价带数你关心的导带数×3但更直接的方式是看自能实部在感兴趣的能量窗口里有没有收敛。VASP输出里会有每个能级的自能修正值你可以对比NBANDS200和NBANDS300的结果差在0.05 eV以内就说明基本收敛。ENCUTGW 200。这个是把GW中关联能部分的平面波截断单独降下来。GW对基组的收敛比DFT快一般取ENCUT的一半就够能省一半以上计算量。这也是VASP官方推荐的省资源手段。OMEGAMAX 50。频率网格的最大值单位是eV。这个值至少要比你关心的能量范围大一些。算带隙的话50 eV绰绰有余要算到深能级或者X射线范围的谱就需要调大。LSPECTRAL .TRUE.。使用谱方法实现GW中的频率卷积比直接频率积分快得多。几乎所有现代GW计算都开这个。LMAXFOCKAE .TRUE.。在GW中确保HF交换部分使用精确的AE基组处理对d/f电子体系尤其重要。过渡金属、稀土化合物建议必须开。NKRED 1。这个参数是把GW自能计算时的K点网格约化NKRED2表示自能部分只用一半K点数。它会显著降低计算量但我建议只在测试阶段用最终精度计算还是设成1。ISYM 0。GW计算建议关闭对称性因为自能算符和介电函数在对称操作下的变换不如DFT那么直观关掉省心。代价是计算量上去了但对体系小的还好。这一步跑完VASP会在目录里生成WAVEDER文件。这个文件极其关键——它是BSE计算必需的前置产物里面存的是占据态与未占据态之间的动量矩阵元。如果这一步用的是旧版本VASP或者你中途改了NBANDS、K点WAVEDER就得重新生成否则后面BSE全都白搭。3. BSE计算把激子效应从谱线里逼出来GW算完文件夹里应该已经有WAVEDER、WAVECAR、CHGCAR这些文件了。BSE的计算就是在保留这些文件的基础上再一次替换INCARVASP会读取已有的WAVEDER构造电子-空穴相互作用矩阵元。3.1 BSE阶段的INCAR设置我的BSE模板长这样SYSTEM BSE PREC Accurate ENCUT 400 ISMEAR 0 SIGMA 0.05 EDIFF 1E-6 ALGO BSE ANTIRES 0 NBANDS 300 NBANDSO 8 NBANDSV 32 OMEGAMAX 20 NOMEGA 100 CSHIFT 0.1 LOPTICS .TRUE. ISYM 0各个参数的逻辑ALGO BSE。直接在已有波函数的基础上解BSE方程。最好确保WAVEDER是从GW那一步生成的而不是从DFT的LOPTICS那一步生成的。两者在物理上都能用但GW之后WAVEDER里的单粒子能级是准粒子能级激子峰位置会更准。NBANDSO和NBANDSV。这两个决定激子哈密顿量的尺寸。NBANDSO是占据带数NBANDSV是未占据带数。别忘了BSE里的价带空穴态和导带电子态都要在这个空间里展开。选多少看你的能量窗口——激子波函数主要由带边附近的能带贡献。对于大多数体系NBANDSO取3~8条占据带、NBANDSV取10~30条空带就够用了。注意这里说的N是每条K点上的带数不是总带数。粗算时可以小一点精算时一定要做收敛性测试连续增大NBANDSV直到吸收谱第一峰位置变化小于0.05 eV。ANTIRES 0。这个参数控制是否包含反共振项。0表示使用Tamm-Dancoff近似只考虑电子从价带到导带的激发正向跃迁不考虑反过程。TDA对大多数半导体光吸收谱来说精度足够而且计算量减半。如果你要算特别精细的线形或者涉及自旋翻转的情况再考虑ANTIRES1。OMEGAMAX和NOMEGA。和GW里的含义类似控制频率网格范围和密度。BSE里这个网格决定最终输出的介电函数虚部的能量分辨率。NOMEGA100配合CSHIFT0.1一般能画出挺光滑的谱线了。我这里OMEGAMAX只到20 eV因为光吸收谱最关心的带边和激子峰都在这个范围以内。CSHIFT 0.1。严格说是洛伦兹展宽的半宽单位eV。它的作用是给谱线一个自然的展宽模拟有限寿命效应和温度效应。展宽太小谱线会剧烈振荡太大又会把激子峰抹平。我建议先取0.1试算看谱形再调。3.2 BSE结果的提取BSE跑完后最重要的结果在vasprun.xml和OUTCAR里。VASP会输出宏观介电函数ε(ω)的实部和虚部频率网格就是你设定的NOMEGA那两百个点为什么是两倍因为参数里隐含了虚频和实频两条网格后面你会看到输出的点多于NOMEGA这是正常的。手写脚本从vasprun.xml里提取介电函数当然可以但我一般直接上VASPKIT。它能直接读取OUTCAR或者vasprun.xml输出光学吸收相关的数据包括实部和虚部介电函数、吸收系数、反射率、折射率这些省去自己写awk的麻烦。具体到功能编号VASPKIT把光学性质相关的模块分成DFT、GW、BSE三个入口选BSE就行注意别和DFT那套弄混。如果你自己处理数据核心逻辑是介电函数虚部ε2(ω)直接对应吸收谱激子峰在ε2(ω)里表现为一个尖锐的强峰能量位置在准粒子带隙以下。我拿二维MoS₂算过一次BSE的ε2第一峰比GW带隙低了大概0.3 eV这个差值就是激子结合能。看到这样的结果基本可以确定BSE这套流程跑通了。3.3 一个容易忽略的细节K点网格与BSE的兼容性BSE计算里有一个隐形的杀手——K点设置和对称性。VASP文档里明确说了BSE必须在以Gamma为中心的K点网格上算且不支持某些降低对称性的操作。我之前用一套在DFT阶段Monkhorst-Pack偏移网格跑得好好的KPOINTS直接接到BSE结果VASP直接报错退出。排查了半天问题就是K点网格不兼容。所以前后期的KPOINTS必须保持一致都用Gamma-centered网格。好在GW阶段我就统一用Gamma网格了所以到了BSE这一步反而没出这个问题。从DFT就开始用Gamma网格可以避免一个隐藏的坑。4. 实操中必然会踩的坑我的排查记录说实话GWBSE刚上手的时候我几乎每个体系都要被参数坑一轮。下面这几个问题基本属于人人都会遇到级别的我把排查过程写出来你遇到的时候可以直接对照。4.1 NBANDS不足导致的准粒子能量发散第一次跑G0W0那会儿我NBANDS开得特小觉得反正我只关心价带顶和导带底要那么多空带干嘛。结果GW跑完一检查OUTCAR里的准粒子修正值发现导带底的修正值有1.5 eV——那根本不是正常值正常半导体G0W0修正一般不会超过1 eV。更离谱的是某个能量更低的未占据态修正值随迭代步数一直在涨完全没有收敛趋势。排查链路是这样的先怀疑是ALGO设置问题把GW0换成严格的G0W0跑了一遍结果一样然后怀疑是不是LSPECTRAL没有开检查INCAR确认开着最后跟学长讨论对方一句话点醒你是不是空带没给够GW自能里那个库仑散射矩阵对高能态特别敏感空带数不够高能态都没有库仑散射压根没法收敛。验证方法也简单把NBANDS从150加到300其他不动重跑一遍。准粒子修正值直接降到0.8 eV稳了。从此我学乖了每次换新体系必做NBANDS收敛性测试。测试方法分别用NBANDS200、300、400跑G0W0对比感兴趣的几根能级的修正值变化小于0.05 eV才算过。4.2 BSE报错WAVEDER not found或者Problem in WAVEDER这个错我愿称之为BSE新手的第一道门槛。BSE计算前必须有一个完整的WAVEDER文件这个文件不是DFT默认产生的必须由GW计算或者LOPTICS.TRUE.计算生成。如果你从DFT那一步直接切BSEVASP大概率直接罢工。我第一次遇到时以为是文件缺失重新看目录WAVEDER就在那里。但VASP还是报Problem in WAVEDER。后来查手册才明白WAVEDER和WAVECAR的配套关系很严格——WAVEDER里记录的波函数信息必须和当前INCAR里的NBANDS、KPOINTS设置一致。我当时的操作是GW阶段用NBANDS300生成WAVEDER到了BSE那步为了省内存我把INCAR里的NBANDS偷偷改成了200结果WAVEDER和当前波函数不匹配直接被拒。解决方式要么BSE阶段NBANDS保持和GW阶段一样要么重新生成WAVEDER。千万别自作聪明改NBANDS。后续类似的事情多了我也学乖了——凡是涉及WAVEDER的操作所有关键参数一律保持和生成它的那一步完全一致。4.3 CSHIFT太大把激子峰抹平了CSHIFT作为洛伦兹展宽设太大确实会改变谱形。但更隐蔽的问题在于如果我设了CSHIFT0.3NOMEGA又不够大谱线在激子峰附近就特别拉胯峰高被削得只剩一半峰位看上去也偏了将近0.1 eV。为什么因为BSE输出的谱是离散频率点上的介电函数值CSHIFT只是给了每个点一个展宽。如果频率点太稀NOMEGA太小或者展宽太大两个邻近点的贡献会互相叠加把尖锐的激子峰平均成一个矮胖子。建议先不管展宽用CSHIFT0.05配NOMEGA200跑一版把真实的峰位峰形看清楚等确定物理内容没问题了再根据你要和实验对比还是画图调整CSHIFT。我最终常用的组合是CSHIFT0.1 NOMEGA100既不抹平物理峰也不会让谱线抖得太厉害。5. 资源估算与计算成本优化经验GWBSE是出了名的贵但也贵得有规律。搞清楚计算量的来源才知道从哪省钱。5.1 不同体系的资源消耗参考我在不同规模体系上攒过一组粗略的参考数据仅供参考具体数值跟机器配置、编译版本、并行方式都有关系体系原子数K点网格G0W0耗时(256核)BSE耗时(256核)Si 原胞26×6×620 min15 minMoS₂ 单层312×12×12 h1 hMAPbI₃ 原胞124×4×45 h2 h异质结界面模型403×3×124 h8 hBSE的耗时大头在构造电子-空穴相互作用矩阵元也就是对每个K点组合计算库仑势和屏蔽交换项的积分。这部分复杂度近似正比于(NBANDSO×NBANDSV×N_k)²K点越多、带数越多涨得越快。所以BSE部分的优化策略从来不是靠堆核心而是靠减少K点或带数。5.2 几个实用的降成本技巧第一ENCUTGW要舍得用。这是VASP官方推荐的做法省掉的计算量非常可观。我试过ENCUTGW取ENCUT的一半甚至0.4倍带隙和吸收谱的变化基本可以忽略。第二NOMEGA在测试阶段降低。NOMEGA影响的是介电函数频率网格的密度它对GW自能计算的收敛有一定影响但在测试阶段你不需要那么高分辨的谱线。先用NOMEGA20把流程跑通、参数试对最后正式计算再调回50甚至80能省不少时间。第三并行设置不要无脑堆核心。VASP在GW部分的并行效率不是线性的核心数超过一定规模后会因为通信开销而性能下降。我用的原则是NCORE设成每个计算节点物理核数的一半到全部然后KPAR根据K点数调整。小K点网格上开太多KPAR反而会拖慢。第四也是我最想强调的一点先在小体系、粗K点网格上跑通全流程再上正式规模。我每次接新体系都先用一个很小的模型把INCAR调通确认能出合理的带隙和激子峰位置再铺开大规模计算。这套流程跑下来能避免很多算了三周才发现INCAR里某个参数写错的悲剧。如果要自己装VASP建议编译时就用Intel编译器加MKL同时把FFTW和MPI的路径配置好。VASP 6.x对GW的实现比5.4.4完善不少尤其在高并行下稳定性更好新算例直接用6.x别在旧版本上花时间。6. 最后掏心窝子的几条经验从DFT到GWBSE这套方法用到现在我最深的体会是GWBSE不是一个跑完就出结果的黑盒子而是一套需要你不断做收敛性测试的精细工具。参数之间的耦合关系很微妙——NBANDS影响GW准粒子能级准粒子能级又直接影响BSE激子峰的绝对位置K点密度既影响GW自能也影响BSE激子波函数在倒空间里的描述。想一次跑通就拿到准确的激子结合能基本不现实。每换一个新体系老老实实把NBANDS、NBANDSO/NBANDSV、K点密度、CSHIFT这几组参数都测试一遍得到的数值才有底气写进论文。另外一个小技巧算BSE之前先看一眼GW那一步产出的准粒子带隙心里有个预期——激子峰一定落在这个带隙以下如果BSE峰跑到带隙以上去了要么是NBANDSO选少了要么是WAVEDER没配对好。这个预期校验能帮你在一分钟内发现大多数流程错误。最后再分享一个连我自己都踩过好几次的坑中间文件的备份。WAVEDER、WAVECAR、CHGCAR这几个文件动辄几百MB到几个GB计算中途硬盘满了或者误删了某个文件从头再算一遍是真的崩溃。我现在都是算完GW立刻把WAVEDER单独备份一份后面BSE就算把INCAR调了几百遍也随时能回到GW刚结束的那个干净状态。这套方法虽然不能让你少写代码但确实能省下几个月的机时。说白了GWBSE就是计算材料里那个门槛高、但过了就通透的坎。参数调通、流程跑顺之后你会发现激子物理、准粒子修正这些概念不再只是教科书上的公式而是你每天都会在OUTCAR和vasprun.xml里看到的具体数字。这篇文章里写的每一条都是我在反复试错里磨出来的希望能让后来的人少走几段弯路。