干我们计算材料这行的,估计都逃不过这么一遭:用DFT把某个半导体的带隙算出来,高高兴兴拿去组会汇报,结果导师一句"这个带隙和实验差了0.8 eV,光学吸收谱的第一峰位置也不对",瞬间把人打回原形。问题出在哪儿?DFT的Kohn-Sham本征值本质上是辅助量,不是准粒子激发能,算带隙天然偏低;要算光吸收谱,还得考虑电子-空穴的库仑相互作用,也就是激子效应。这时候,VASP里的GW+BSE组合就成了绕不开的解决方案。
这套方法算下来,带隙精度能拉到0.1 eV量级,吸收谱的激子峰位置也能和实验对上。GW负责把DFT的单粒子能级修正成准粒子能级,BSE负责在准粒子能级基础上解电子-空穴束缚态方程,两者配合,才算把一个材料的光电性质从"定性对"推进到"定量对"。
这篇文章就围绕"用VASP算GW+BSE"这条线,从方法动机、参数设置、实操流程到踩坑记录,把我攒下的经验一次性倒出来。适合刚接触激发态计算、准备给体系上GW+BSE的研究生,也适合被DFT结果搞到怀疑人生、想换方法的老手参考。
1. 为什么DFT算不准带隙之后,我会转向GW+BSE
先说清楚一个容易混淆的点:DFT不是一个"坏"方法,只是它的Kohn-Sham轨道能级在物理上并没有被严格定义成电子激发能。严格来说,Kohn-Sham本征值只是拉格朗日乘子,除了最高占据态在渐近极限下有些意义,其他能级和真实的光电子谱并没有严格对应关系。所以,算出来的带隙偏小不是"误差",而是方法本身的局限。
那GW为什么能修?因为GW通过自能算符Σ对单粒子格林函数做微扰修正,把"电子感受到的其他电子的交换关联效应"重新整理了一遍。自能里面有G(单粒子格林函数)和W(屏蔽库仑相互作用),G负责描述电子-空穴对的传播,W负责描述电子在介质中感受到的动态屏蔽相互作用。把这两个东西卷积在一起,得到的就是准粒子色散关系,准粒子带隙自然比KS带隙更接近实验。
再来说BSE。GW只解决单粒子激发的问题,也就是一个电子被激发到导带之后的能级位置。但光吸收是一个双粒子过程:一个电子从价带跃迁到导带,原来的价带位置留下一个空穴,电子和空穴之间因为有库仑吸引,可能束缚在一起,这个束缚态就是激子。激子能级比GW算出的带边低一个结合能,在吸收谱里表现为带边以下或带边附近的尖锐峰。这些信息单粒子GW给不了,必须解BSE方程。
所以一句话总结我的经验:只关心带隙,GW够用;关心吸收谱形状、激子峰位置、光学响应,必须上BSE。
顺带说一句,什么时候可以不上GW+BSE?如果你只是做结构优化、算形成能、看电荷转移这类基态性质,DFT和DFT+U完全够用。材料带隙不算太小、激子结合能又很弱(比如很多三维无机半导体),用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是常规操作,但绝缘体和半导体用ISMEAR=0(即高斯展宽)没问题。如果是金属,后面算GW会遇到很多麻烦,因为费米面处的占据数在GW里怎么处理都不太舒服,这也是为什么GW+BSE主要应用场景是半导体和绝缘体。
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输出里会有每个能级的自能修正值,你可以对比NBANDS=200和NBANDS=300的结果,差在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点网格约化,NKRED=2表示自能部分只用一半K点数。它会显著降低计算量,但我建议只在测试阶段用,最终精度计算还是设成1。
- ISYM = 0。GW计算建议关闭对称性,因为自能算符和介电函数在对称操作下的变换不如DFT那么直观,关掉省心。代价是计算量上去了,但对体系小的还好。
这一步跑完,VASP会在目录里生成WAVEDER文件。这个文件极其关键——它是BSE计算必需的前置产物,里面存的是占据态与未占据态之间的动量矩阵元。如果这一步用的是旧版本VASP,或者你中途改了NBANDS、K点,WAVEDER就得重新生成,否则后面BSE全都白搭。
3. BSE计算:把激子效应从谱线里"逼"出来
GW算完,文件夹里应该已经有WAVEDER、WAVECAR、CHGCAR这些文件了。BSE的计算就是在保留这些文件的基础上,再一次替换INCAR,VASP会读取已有的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对大多数半导体光吸收谱来说精度足够,而且计算量减半。如果你要算特别精细的线形或者涉及自旋翻转的情况,再考虑ANTIRES=1。
- OMEGAMAX和NOMEGA。和GW里的含义类似,控制频率网格范围和密度。BSE里这个网格决定最终输出的介电函数虚部的能量分辨率。NOMEGA=100配合CSHIFT=0.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. 实操中必然会踩的坑:我的排查记录
说实话,GW+BSE刚上手的时候,我几乎每个体系都要被参数坑一轮。下面这几个问题基本属于"人人都会遇到"级别的,我把排查过程写出来,你遇到的时候可以直接对照。
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收敛性测试。测试方法:分别用NBANDS=200、300、400跑G0W0,对比感兴趣的几根能级的修正值,变化小于0.05 eV才算过。
4.2 BSE报错"WAVEDER not found"或者"Problem in WAVEDER"
这个错我愿称之为BSE新手的第一道门槛。BSE计算前必须有一个完整的WAVEDER文件,这个文件不是DFT默认产生的,必须由GW计算或者LOPTICS=.TRUE.计算生成。如果你从DFT那一步直接切BSE,VASP大概率直接罢工。
我第一次遇到时以为是文件缺失,重新看目录,WAVEDER就在那里。但VASP还是报"Problem in WAVEDER"。后来查手册才明白,WAVEDER和WAVECAR的"配套关系"很严格——WAVEDER里记录的波函数信息必须和当前INCAR里的NBANDS、KPOINTS设置一致。我当时的操作是:GW阶段用NBANDS=300生成WAVEDER,到了BSE那步,为了省内存,我把INCAR里的NBANDS偷偷改成了200,结果WAVEDER和当前波函数不匹配,直接被拒。
解决方式:要么BSE阶段NBANDS保持和GW阶段一样,要么重新生成WAVEDER。千万别自作聪明改NBANDS。后续类似的事情多了,我也学乖了——凡是涉及WAVEDER的操作,所有关键参数一律保持和生成它的那一步完全一致。
4.3 CSHIFT太大把激子峰"抹平"了
CSHIFT作为洛伦兹展宽,设太大确实会改变谱形。但更隐蔽的问题在于:如果我设了CSHIFT=0.3,NOMEGA又不够大,谱线在激子峰附近就特别拉胯,峰高被削得只剩一半,峰位看上去也偏了将近0.1 eV。
为什么?因为BSE输出的谱是离散频率点上的介电函数值,CSHIFT只是给了每个点一个展宽。如果频率点太稀(NOMEGA太小),或者展宽太大,两个邻近点的贡献会互相叠加,把尖锐的激子峰平均成一个矮胖子。
建议:先不管展宽,用CSHIFT=0.05配NOMEGA=200跑一版,把真实的峰位峰形看清楚;等确定物理内容没问题了,再根据你要和实验对比还是画图,调整CSHIFT。我最终常用的组合是CSHIFT=0.1 + NOMEGA=100,既不抹平物理峰,也不会让谱线抖得太厉害。
5. 资源估算与计算成本优化经验
GW+BSE是出了名的"贵",但也贵得有规律。搞清楚计算量的来源,才知道从哪省钱。
5.1 不同体系的资源消耗参考
我在不同规模体系上攒过一组粗略的参考数据,仅供参考,具体数值跟机器配置、编译版本、并行方式都有关系:
| 体系 | 原子数 | K点网格 | G0W0耗时(256核) | BSE耗时(256核) |
|---|---|---|---|---|
| Si 原胞 | 2 | 6×6×6 | 20 min | 15 min |
| MoS₂ 单层 | 3 | 12×12×1 | 2 h | 1 h |
| MAPbI₃ 原胞 | 12 | 4×4×4 | 5 h | 2 h |
| 异质结界面模型 | 40 | 3×3×1 | 24 h | 8 h |
BSE的耗时大头在构造电子-空穴相互作用矩阵元,也就是对每个K点组合计算库仑势和屏蔽交换项的积分。这部分复杂度近似正比于(NBANDSO×NBANDSV×N_k)²,K点越多、带数越多,涨得越快。所以BSE部分的优化策略从来不是靠堆核心,而是靠减少K点或带数。
5.2 几个实用的降成本技巧
第一,ENCUTGW要舍得用。这是VASP官方推荐的做法,省掉的计算量非常可观。我试过ENCUTGW取ENCUT的一半甚至0.4倍,带隙和吸收谱的变化基本可以忽略。
第二,NOMEGA在测试阶段降低。NOMEGA影响的是介电函数频率网格的密度,它对GW自能计算的收敛有一定影响,但在测试阶段你不需要那么高分辨的谱线。先用NOMEGA=20把流程跑通、参数试对,最后正式计算再调回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到GW+BSE,这套方法用到现在,我最深的体会是:GW+BSE不是一个"跑完就出结果"的黑盒子,而是一套需要你不断做收敛性测试的精细工具。参数之间的耦合关系很微妙——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刚结束的那个干净状态。这套方法虽然不能让你少写代码,但确实能省下几个月的机时。
说白了,GW+BSE就是计算材料里那个"门槛高、但过了就通透"的坎。参数调通、流程跑顺之后,你会发现激子物理、准粒子修正这些概念不再只是教科书上的公式,而是你每天都会在OUTCAR和vasprun.xml里看到的具体数字。这篇文章里写的每一条,都是我在反复试错里磨出来的,希望能让后来的人少走几段弯路。