先说一个很多做计算材料的朋友都撞过的南墙你把自己的“本命体系”结构优化好了能带也算了兴冲冲地把带隙和实验值对比——结果GGA/LDA算出来的带隙低了百分之三四十。更别说做吸收光谱时第一吸收峰的位置、强度、线形全对不上。问题出在哪DFT本质上算的是基态性质Kohn-Sham本征值并不严格等于电子激发能而光吸收过程里电子和空穴之间还有很强的关联也就是激子效应。要相对严格地处理这类问题就得请出GW和BSE这两个名字很唬人、实际用起来却很有套路的方法。这篇博文我就用VASP把GWBSE这套组合拳完整拆开先讲清楚它们到底在算什么再说怎么一步步跑通最后把我这几年踩过的参数坑、报错坑、并行坑都列出来。内容偏实操适合已经有DFT计算基础、现在被激发态或光学性质搞到头大的朋友。环境方面如果你正在Ubuntu上装VASP、纠结编译和MPI库配置也有对应的小节可以看。1. GWBSE到底在解决什么问题1.1 DFT的遗憾Kohn-Sham能级不是激发能在DFT框架里Kohn-Sham方程给出的本征值 (\varepsilon_{n\mathbf{k}}) 只是一个数学构造严格来说它并不是电子的激发能或者准粒子能量。LDA/GGA这种局域或半局域泛函更是自带诟病——它会在电子自相互作用上犯错导致价带往上偏、导带往下偏结果就是带隙被系统性低估。举个最常被引用的例子硅的实验带隙约1.17 eVLDA/GGA算出来只有0.6 eV左右误差将近一半。氧化物、氮化物这类体系更夸张强关联的过渡金属氧化物甚至可能被算成金属。普通杂化泛函HSE能改善不少但它毕竟还是泛函修正路线对激发态的描述仍然缺乏严格的物理出发点。如果你只是要一个大致带隙值HSE通常够用。但一旦涉及到“电子被激发之后发生了什么”——比如吸收光谱、激子束缚能、发光过程——那就绕不开多体微扰论了。GW和BSE正是从这一层物理出发的工具。1.2 GW近似给准粒子穿上“自能修正”的外衣GW这个名字乍看像某个软件版本号其实它指的是自能 (\Sigma iGW) 的近似形式。这里 (G) 是单粒子格林函数(W) 是动态屏蔽库仑相互作用。你不需要第一时间吃透整个推导链只需要抓住一个核心思想电子在固体里并不是孤零零地在平均势场里运动它周围的环境会被极化、会被屏蔽电子自己也会在周围产生极化云这套复杂作用统一体现在“自能”里。自能 (\Sigma) 与能量有关、与波函数有关把自能代入准粒子方程求解后得到的能级就是带上多体修正的准粒子能量。G0W0是最常见的一种实现方式——(G) 和 (W) 都是从DFT波函数出发一次性构建的然后算一次自能修正就结束不做自洽。这就像一个用“标准件”估算成本的办法速度尚可大多数体系精度也不错。VASP里ALGOGW0对应的就是G0W0。实际经验是它能把半导体的带隙误差从LDA/GGA的30%~50%压缩到0.1~0.2 eV量级。比如硅的G0W0带隙能到1.2 eV左右已经相当贴近实验值。1.3 BSE把电子和空穴绑在一起看光学响应GW解决的是“准粒子能级”的问题可光吸收过程还有一个关键角色没登场——激子。光激发会同时产生一个导带电子和一个价带空穴它们之间带着库仑吸引可能形成束缚态。这个束缚态的能级落在带隙里面会在吸收谱上产生额外峰甚至改变整个光谱线形。Bethe-Salpeter方程BSE就是描述电子-空穴两粒子关联的方程。它把电子和空穴看成一对“若即若离的舞伴”对角项算的是从准粒子带隙出发的单粒子跃迁非对角项算的是它们之间的直接相互作用和交换作用。把BSE哈密顿量对角化之后你得到的就是包含激子效应的激发态本征值和振子强度再求和就能得到宏观介电函数。形象一点说DFT是一张地形图GW告诉你每座山准确的海拔而BSE告诉你登山绳把两个人绑在一起之后这个“绑定组合”的振动模式是什么。三条信息叠起来你才真正看清体系在光照下的完整行为。2. 用VASP跑GWBSE的整体设计2.1 为什么选VASP做这套计算能实现GW和BSE的代码不少BerkeleyGW、Gaussian、Quantum ESPRESSO加插件、FHI-aims、ABINIT都能做。VASP的优势在于平面波基组下收敛性好赝势库成熟补全GW之后又补了BSE模块一套代码不用来回倒腾格式。更重要的是VASP把“从DFT到GW再到BSE”的文件依赖关系串得很顺。你只需按顺序跑几个单点任务WAVECAR、WAVEDER这些中间文件会自然被后续步骤读取。对于习惯一门代码打天下的人这条路的学习成本最低。VASP另一个实际优点是对并行规模控制比较灵活既能在工作站上小跑也能在集群上扩展。不过要提醒一句GW和BSE的内存消耗和带宽开销超出普通DFT非常多不是拿跑结构优化的节点配置来就能直接撑住的。2.2 我为什么把流程拆成三步走VASP官网手册其实也给出了标准的GWBSE流程但很多新手一上来就卡在“为什么我这一步A参数和那一步不一样”的困惑里。我自己习惯把整个计算拆成三步第一步结构优化和DFT基态自洽。这一步做到结构合理、波函数正确同时把WAVEDER算出来。WAVEDER是后续GW和BSE都要用的跃迁矩阵元文件后面会详细讲。第二步G0W0准粒子修正。读入DFT的WAVECAR和WAVEDER算自能输出修正后的准粒子能级。这一步是整个流程中计算量最大、也是最容易“翻车”的环节。第三步BSE求解。读入GW完成后的波函数和WAVEDER构建并求解BSE哈密顿量输出介电函数和吸收光谱。这一步的时间和内存主要取决于你允许BSE纳入多少条占据带和导带。文件传递关系我后面会集中说但时刻记住这三步的k点网格、NBANDS、ENCUT、POTCAR必须保持一致。不一致轻则文件读不进去重则给你算出一堆看似漂亮但完全错误的结果。2.3 文件传递关系与逻辑顺序VASP里做GWBSE最重要的中间文件就两个WAVECAR和WAVEDER。WAVECAR存的是波函数和能级静态DFT结束之后每一步后续任务都在这个基础上继续。WAVEDER存的是占据态与空带之间的动量矩阵元由LOPTICS.TRUE.触发计算GW和BSE都需要它。实际操作的顺序是这样静态DFT跑完目录里已经有WAVECAR和WAVEDER。把INCAR改成GW参数重新提交VASP会读到WAVECAR基于DFT波函数算GW。算完GW之后这个WAVECAR里已经更新成GW的准粒子能级信息。接着再改INCAR为BSE参数再次提交VASP依旧从WAVECAR出发结合WAVEDER去构建BSE哈密顿量。所以你会发现三个步骤是“接力”而不是“重跑”。这也是为什么我不建议中途清掉WAVECAR或者执着地换个INCAR就从零开始。你只需要严格按流程走VASP会自己完成大部分工作。还有个容易忽略的点WAVEDER文件体积不小它在不同版本VASP中的格式可能不同。如果你换了VASP版本最好把WAVEDER删掉在静态DFT步重新生成否则BSE步骤很可能报版本不兼容的错误。3. 完整实操Si的GWBSE计算3.1 第一步结构优化与DFT基态计算为了让例子简单又典型我用硅的金刚石结构来演示。这是教科书级别的一步面心立方原胞里两个Si原子分数坐标分别放在(0,0,0)和(0.25,0.25,0.25)。先做结构优化。INCAR可以这么写SYSTEM Si bulk optimization ENCUT 400 PREC Accurate EDIFF 1E-6 EDIFFG -0.01 IBRION 2 ISIF 3 NSW 60 ISMEAR 0 SIGMA 0.05 LREAL .FALSE.这里有几个关键选择。ENCUT400 eV对硅的PAW-PBE势来说已经较充裕实际上收敛测试可以从300开始看总能和带隙随截断能的变化。ISMEAR0配合SIGMA0.05适合半导体避免金属展宽带来的虚假占据。LREAL.FALSE.是我特意强调的GW和BSE阶段必须用倒空间投影所以结构优化阶段就养成关掉实空间投影的习惯免得后面忘记。KPOINTS用Gamma-centered的6×6×6网格。硅是间接带隙半导体6×6×6对能带和GW都够起步后面做收敛测试时可以升到8×8×8。优化收敛后把INCAR改成静态DFT因为我们要生成高质量WAVECAR和WAVEDERSYSTEM Si static for GW ENCUT 400 PREC Accurate EDIFF 1E-7 ISMEAR 0 SIGMA 0.05 LREAL .FALSE. LORBIT 11 LOPTICS .TRUE. NELM 200 NBANDS 128注意这里出现了NBANDS128。硅一个原胞两个原子、每个原子4个价电子总共8个电子占据4条带。普通DFT默认NBANDS大约是带数加一小部分空带但GW需要大量空带来计算屏蔽和自能。我直接设成128意思是给体系准备120多条空带这在GW里是比较充足的起点。LOPTICS.TRUE.的目的是生成WAVEDER文件。这里有个常见的坑如果你在静态DFT这一步才想起LOPTICS并且计算的k点很多WAVEDER的生成会额外花不少时间所以有人会在结构优化前就先跑一遍带LOPTICS的静态。不过只要不是频繁改动结构按“优化完再静态”的顺序完全没问题。POSCAR和POTCAR按常规准备。POTCAR用PBE势硅可以用默认的Si势也可以考虑Si_sv版本后者把半芯态放进行价带在GW里往往能改善精度代价是计算量变大。我一般先用默认Si势跑通全流程再回头评估是否需要换更贵的势。3.2 第二步G0W0准粒子修正静态DFT跑完后目录里应该已经有WAVECAR和WAVEDER。现在修改INCARSYSTEM Si G0W0 ENCUT 400 PREC Accurate EDIFF 1E-7 ISMEAR 0 SIGMA 0.05 LREAL .FALSE. ALGO GW0 NBANDS 128 NOMEGA 16 ENCUTGW 200 NELM 200 LOPTICS .TRUE.ALGOGW0是G0W0计算的关键。NBANDS要跟静态DFT保持一致不然后面环节读取WAVECAR会出问题。NOMEGA是频率网格数量VASP默认是16很多体系用16起步问题不大但严谨些要测24、32看准粒子能级是否稳定。ENCUTGW是自能计算时的平面波截断我习惯取ENCUT的一半这里是200 eV然后做一次ENCUTGW300的测试比较结果。NELM设到200是为了避免GW自能迭代过程中达到默认电子步数上限而被截断。这一步跑完你会看到OUTCAR里记录了准粒子能级修正。我的经验是别急着看光谱先盯几件事第一GW总能或者每个带的自能修正是否已经收敛第二带隙是否从PBE的0.6 eV左右提升到1.1~1.2 eV附近。如果你看到的带隙几乎没变化多半是NBANDS不够或者ENCUTGW太低。GW这一步消耗的资源是三个步骤里最大的。对于6×6×6的k点网格128条带硅这个体系不算大普通8核16核的机器也能在合理时间跑完。但如果你换成一个几十个原子的界面模型NBANDS轻松上千内存需求会暴涨后面的并行部分必须提前留心。3.3 第三步BSE求解吸收光谱GW完成后同一目录里已经有带GW信息的WAVECAR和WAVEDER。修改INCAR进入BSESYSTEM Si BSE ENCUT 400 PREC Accurate EDIFF 1E-7 ISMEAR 0 SIGMA 0.05 LREAL .FALSE. ALGO BSE NBANDS 128 NBANDV 4 NBANDC 24 CSHIFT 0.1 LSPECTRAL .TRUE. LOPTICS .TRUE.这里多余的两个参数是NBANDV和NBANDC它们是BSE专属的。NBANDV4表示把所有占据带硅只有4条都纳入电子-空穴对空间NBANDC24表示取24条导带参与BSE计算。导带取数越多激发态描述越全但矩阵规模会急剧扩大因为BSE问题的维度近似正比于NBANDV×NBANDC×k点数。CSHIFT是给介电函数分母加的一个复常数相当于人为展宽。取0.1 eV配合实验光谱常用的洛伦兹展宽很合适。如果你想看清楚精细的激子峰结构可以减到0.05如果体系较大、数值噪声明显可以加到0.2。LSPECTRAL.TRUE.的作用是让VASP输出谱函数方便进一步分析和绘图。BSE计算跑完后宏观介电函数 (\varepsilon_2(\omega)) 会在vasprun.xml里提取方法我会在后处理部分讲。3.4 可直接抄作业的INCAR参数速查表我把三阶段的INCAR和关键参数汇总成表方便对照。任务阶段关键参数说明与理由结构优化IBRION2, ISIF3, EDIFFG-0.01共轭梯度法弛豫优化晶胞与原子静态DFTLOPTICS.TRUE., NBANDS128, EDIFF1E-7生成高质量波函数和WAVEDERG0W0ALGOGW0, NOMEGA16, ENCUTGW200从DFT出发算一次自能修正BSEALGOBSE, NBANDV4, NBANDC24求解激子方程输出介电函数全程ENCUT400, LREAL.FALSE.保证一致性与数值可靠性这张表是“起点”不是“终点”。每个数值背后都有测试空间。比如遇到更宽带隙的氧化物NBANDC可能需要加到30甚至更多遇到窄带隙但仍需激子信息的体系NBANDV要仔细核对包含哪些占据带。宁可先小规模跑通再逐步加参数不要一上来就压上全部资源。4. 常见报错与排查实录4.1 不收敛问题的调整思路先回忆我遇到最多的一类RMM-DIIS不收敛或者GW步骤里出现Sub-Space-Matrix非厄米的提示。这类问题的核心往往不是某一个参数而是“当前任务没给出足够数值空间”。如果发生在DFT静态阶段检查ENCUT是否太低、SIGMA是否和ISMEAR匹配、初始磁矩或电荷密度是否合理。半导体体系用ISMEAR0配SIGMA在0.02~0.1之间一般没问题。实在不收敛可以试试AMIX0.2、BMIX0.0001这种保守混合参数。如果发生在GW阶段第一怀疑对象就是NBANDS。GW对空带数量极其敏感空带不够时自能的高能贡献被截断不仅不收敛算出的带隙也会明显偏低。第二是NOMEGA频率网格太疏会导致积分振荡我遇到过把NOMEGA从16调到24之后整个循环立刻稳定的情况。BSE步骤的收敛问题则更多表现为“光谱随CSHIFT剧烈变化”。这通常是展宽太小叠加数值噪声导致的适当增大CSHIFT就能缓解。还有一个不太起眼的因素k点网格。BSE的激子束缚能对布里渊区采样很敏感6×6×6不够就上8×8×8代价是计算量成倍上涨所以要平衡。4.2 内存与并行配置问题GW和BSE是VASP里最吃内存的计算之一。普通DFT用NCORE8甚至NCORE16都很常见但到了GW阶段同样的设置很容易直接把节点内存榨干。GW计算每个MPI rank需要存自能、格林函数、屏蔽相互作用等大矩阵内存需求大致随NBANDS的增长呈平方甚至更高次方增长。硅这种128条带的小体系还好几百条带以上时动辄几十GB的内存消耗都不稀奇。我一般先估算当前节点总内存再决定进程数。比如64核节点配256GB内存我会减少到8~16个MPI进程配合每进程4~8个OpenMP线程而不是盲目用满64核。另一个常见问题是“Error EDDDAV: Call to ZHEGV failed”。这类对角化报错大概率是内存不足或者进程间分布不均导致。排查的顺序是先降NBANDS试试再看NCORE和NPAR有没有设得过大最后检查是不是多个任务抢同一批节点。并行配置的具体经验是GW阶段别迷信核数BSE阶段也类似。跑大规模BSE前先用小k点和少量导带做一次短测试记录峰值内存用这个数据推算正式任务需要的节点数。实测下来这个“先探后打”的方式比直接提交正式任务然后被OOM杀掉省时省钱得多。4.3 结果异常的排查方向如果你的GW带隙算出来比PBE还低或者BSE光谱首峰位置离谱大概率和参数无关而是某些基础设置被破坏了。第一个要查的是WAVEDER和WAVECAR是否匹配。它们必须来自同一组k点、同一套NBANDS、同一个ENCUT。如果你在中间重新跑过静态DFT但没删掉旧WAVEDERBSE阶段就可能读到过期的跃迁矩阵元产生完全错乱的光谱。处理方式简单粗暴删掉WAVEDER重跑一次带LOPTICS的静态DFT。第二个要查的是POTCAR。换用不同POTCAR比如从Si换成Si_sv后WAVECAR和WAVEDER都必须重新生成因为它们都跟投影无关跟具体的赝势和电子数绑得很紧。偷懒不重算的结果就是带位置漂移、光谱峰错位。第三个要查的是对称性。BSE求解对外界电场、偶极近似下的跃迁选择规则非常敏感如果KPOINTS设置与静态DFT不一致或者用户手工改过WAVEDER生成时的对称性设置报错和异常结果都算轻的。我的习惯是从静态DFT到GW到BSE除了INCAR里的ALGO等必要改动其余涉及k点、NBANDS、POSCAR、POTCAR的字段全部原样保留。4.4 Ubuntu上装VASP的环境坑最后提一句安装和编译。越来越多人在自己的Ubuntu工作站上装VASP官方源码需要合法license这点按下不表只说编译环境的常见坑。VASP编译最常见的问题是MKL和MPI库路径不对。用Intel oneAPI时记得先source setvars.sh然后在makefile.include里正确设置BLAS、LAPACK、FFTW路径。很多报错——比如找不到libmkl_intel_lp64.so、链接阶段报重复符号——都和顺序或版本有关。我建议的顺序是先确认ifort/icx版本再确认mpiifort版本最后确认MKL版本三者尽量来自同一套oneAPI版本否则容易遇到ABI不兼容。编译时用make veryclean清一次旧对象再make能避免不少诡异的“Object file not found”问题。至于VASP版本5.4.4是很多论文的老搭档6.x系列对GW和BSE做了大量优化支持更好。如果刚接触我建议直接用6.x最新版前提是服务器上有匹配的编译器。5. 后处理从输出文件到发表级光谱5.1 从OUTCAR和vasprun.xml读GW结果GW计算完成后准粒子能量的修正信息在OUTCAR里都有记录。搜索quasiparticle或QP-en窄关键字你会看到每个k点、每条带对应的自能修正值以及G0W0修正后的本征能量。很多文章需要报告“G0W0带隙比PBE提高了多少”这一步直接从文本里提取即可。如果你想批量提取或画能带用vaspkit的能带提取功能会更方便。它能直接读取vasprun.xml中的本征值数据再跟GW的准粒子修正合并出带结构图。需要注意GW修正通常是按每个k点每个带单独列出不能简单地整体平移带边。5.2 提取BSE吸收光谱与激子峰BSE的计算结果都汇总在vasprun.xml里。介电函数的实部(\varepsilon_1)和虚部(\varepsilon_2)会作为频率的函数输出。提取时可以直接用Python做一个简单的SAX解析找dielectric function字段也可以看你熟悉的后处理工具是否支持。展宽的处理值得说两句。CSHIFT在计算里起了一个人为洛伦兹展宽的作用你提取出的(\varepsilon_2)本身就是有展宽的曲线直接作图就是一条连续光谱。如果你想得到更接近实验的谱形可以在后处理再做一次卷积。但一定要意识到较大展宽会淹没激子峰的细节特别是束缚能只有几十meV的弱束缚激子小展宽才能揭示真实结构。硅的BSE光谱里最有名的是临界点附近的吸收增强以及第一直接带隙跃迁对应的峰位。如果看到首峰落在带隙以下约0.1~0.2 eV并且强度明显高于普通单粒子跃迁的预测这就是激子效应在起作用。拿这个曲线跟实验椭圆偏振谱对比验证模型是否靠谱。5.3 论文里如何报告这些计算参数审稿人盯GW和BSE盯的往往不是结果本身而是参数是否收敛。我建议在论文方法部分写明这几个数DFT阶段的ENCUT、k点网格、NBANDSGW阶段的NOMEGA、ENCUTGW、自洽方式G0W0还是其他BSE阶段的NBANDV、NBANDC、CSHIFT。更重要的是给出至少一个收敛性测试结果表明这些参数已经可靠。比如固定其他参数把NBANDS从128加到256带隙变化小于0.02 eV把NOMEGA从16加到24光谱基本不变。这些测试放在补充材料里就是审稿人挑不出毛病的底气。如果不做收敛测试直接给一个“好看”的光谱风险在于你自己都不知道这个结果是物理还是数值噪声。GWBSE这种级别的计算重复成本很高出了文章后想补测代价非常大。宁可先花一周做参数扫描也别省这个时间。6. 最后几句经验话我在实际使用中最大的体会是GWBSE这套方法最忌讳“一步到位”的心态。它不像普通DFT那样随便开几个默认参数就能跑出像样的结果NBANDS、NOMEGA、ENCUTGW、k点网格每一步都在互相制约必须耐心做收敛测试。硅这种小体系跑起来轻松可以大胆尝试不同参数组合一旦换到界面、缺陷或大分子体系资源有限的情况下更要先用小模型摸清参数规律再上大任务。还有一个容易忽略的小技巧BSE计算之前花一点时间确认WAVEDER的确是最近生成的。很多人被“WAVEDER not compatible”这种报错折磨一整天最后发现只是忘了在静态DFT阶段打开LOPTICS或者中间换过VASP版本。这类问题排查起来并不难难的是下决心把目录清干净、重跑一次静态计算。老实说重跑一次静态DFT远比在错误结果上反复检查更省时间。另外结果判读时一定把“带隙修正”和“激子效应”分开看。GW管的是单粒子能级的修正BSE管的是电子-空穴关联带来的额外束缚。一个体系的光吸收峰位置是两者共同作用的结果如果你拿着BSE光谱去跟实验峰比对偏差较大先判断是哪个环节出了问题。我以前激进地把CSHIFT调得很小想看清激子峰结果数值噪声被放得很大反而掩盖了真实峰位后来老老实实把展宽调到合理范围结果立刻清晰起来。这套流程后续可以扩展的方向其实很广算缺陷态的吸收谱、异质结的光学响应、有机分子的激发态、溶剂模型下的光谱模拟都可以在同一套VASP框架内做。关键还是一开始把基态、GW、BSE三步之间的文件和参数关系理清楚后面的路就会顺很多。