ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

VASP中GW+BSE方法计算材料激发态光谱:从原理到实操

VASP中GW+BSE方法计算材料激发态光谱:从原理到实操 1. 项目概览与适用场景做材料计算这几年我遇到最多的问题不是“怎么跑程序”而是“跑完DFT之后带隙和实验对不上”。尤其碰到半导体和绝缘体体系用PBE算出来的带隙普遍偏低有的甚至直接给成了金属这种情况做能带分析、缺陷态评估、光学性质预测基本等于白算。后来接触到GWBSE这套组合才算是把“基态算得差不多、激发态完全没谱”这个老大难问题真正打通了。VASP是目前第一性原理计算里用得最顺手的平面波软件之一GWBSE是它内置的高精度激发态计算方法。简单说GW负责修正准粒子能级把DFT那套被低估的带隙拉回实验值附近BSE负责在GW基础上处理电子-空穴相互作用也就是激子效应最后给出可以和实验吸收光谱、损失函数直接对照的数据。这两个方法叠加是当前计算材料光谱性质比较主流且结果可靠的一整套方案。这篇文章面向的是已经能跑通VASP常规DFT计算、但对GWBSE流程还不够熟悉的从业者和研究生。我会从原理、参数、实操到报错排查完整走一遍尽量把每个参数为什么这么设、哪个步骤容易翻车、哪个文件必须保留这类细节讲透。跟着这套流程走你至少能少踩一半的坑。2. 核心原理深挖GW和BSE分别解决什么问题2.1 GW近似把DFT的带隙误差拉回来要理解GW先得理解DFT的尴尬。Kohn-Sham DFT的能带本征值严格来说不是准粒子激发能它的带隙误差主要来源于交换关联泛函对自相互作用的错误处理以及缺少频率依赖的屏蔽效应。PBE这类泛函把电子间的排斥作用平均化了导致电子觉得“自己没有被自己排斥”费米能级附近的占据态被抬得不够高、空态被放得过低带隙自然就偏小。GW近似做的事情是把自能算符写成一个动态屏蔽的交换作用形式也就是Green函数G和屏蔽库仑作用W的乘积在里面显式包含频率依赖的屏蔽效应。准粒子方程里面用这个自能算符替换DFT的交换关联势解出来的本征值才是真正的电子加空穴激发能。实际VASP计算里最常见的是G0W0也叫一次性GW也就是用DFT的波函数和本征值固定不动用它们构造G和W只算一遍自能修正。G0W0对大多数半导体和绝缘体已经足够准带隙修正效果非常好而且比自洽GW便宜很多是我个人默认选择。不过G0W0依赖初始的DFT波函数质量如果体系有强关联电子比如含有过渡金属d电子或稀土f电子的材料那么PBE初始波函数本身就不太对一次性GW修正常常不收敛或者偏离实验很大。这种体系一般需要先做DFTU或者改用杂化泛函提供更好的初始波函数实在不行再考虑evGW或QSGW这类自洽方案。VASP里ALGOGW0对应的是部分自洽GW0它对本征值做了自洽但固定了WALGOGW则是W也一起自洽计算量会成倍上涨普通体系性价比不高。2.2 BSE处理电子-空穴相互作用GW算完准粒子能级准了但还缺一个关键物理光吸收过程中电子被激发到导带之后它和价带留下的空穴之间是有库仑吸引的这个被束缚的电子-空穴对就是激子。GW描述的是单粒子激发BSE才是处理双粒子电子空穴激发的方程两者物理图像完全不同。BSE的核心思路是把电子-空穴格林函数的运动方程写成类似两点散射问题的形式这里面包含三部分电子和空穴各自传播的准粒子能级项、电子-空穴交换作用对应亮激子/暗激子的选择规则以及电子-空穴吸引项主要由屏蔽库仑作用W构成。VASP求解BSE的时候会先用GW步骤得到W和准粒子波函数然后把激子哈密顿量显式对角化得到激子本征态和对应的光学跃迁强度最终输出介电函数虚部也就是吸收光谱。为什么不能用DFT的波函数直接做BSE因为DFT带隙本身就低估了激子束缚能在很多有机半导体里能到0.5甚至1 eV以上两个误差叠在一起光谱整体位置会明显偏离实验。先做GW把带隙和波函数修正了BSE里的单粒子项才是在一个可靠平台上这样才能确保算出来的激子束缚能和吸收峰位置接近实验值。我算过几种有机-无机杂化钙钛矿如果用PBE拿到的能级直接算光学性质吸收峰普遍偏低0.5 eV以上GWBSE之后基本能和实验光谱对到误差0.1 eV以内。这套组合值钱就值钱在这里。3. 环境准备与VASP编译常见问题3.1 Ubuntu下VASP编译的前置依赖VASP收费且License管理严格这是绕不开的客观条件我这里默认你已经拿到了合法授权只聊编译和运行环境的坑。我自己的主力编译环境是Ubuntu 22.04 LTS Intel oneAPI这套组合比GCCOpenBLAS更容易发挥VASP的峰值性能尤其GWBSE这种对BLAS库和FFT性能极其敏感的任务快慢能差出30%以上。编译VASP之前有几个前置工具必须备齐Intel oneAPI Base Kit和HPC Kit里面包含ifort/ifx编译器、MKL数学库和MPI库或者是GCC OpenMPI OpenBLAS FFTW的组合性能略低但在纯AMD平台上反而更稳定GNU Make和CMake取决于你用的VASP版本老版本5.4.4一般用make新版本6.3也可以走CMake实际编译时最容易翻车的地方是MKL库接口不匹配。VASP 5.4.4默认用ifort编译但在新版oneAPI里ifort已经被ifx替代很多旧makefile里的MKL链接写法在新版本编译器中会提示找不到mkl_blas95_ilp64这类库。我的解决方法是直接改用6.x系列的CMake构建流程它会自动探测MKL接口比手工改makefile省心得多。3.2 新增配置与硬件估算如果你用的是VASP 6.x编译前需要确认把-Dtbdyn、-Dscalapack这些选项打开GW和BSE本身不需要额外特殊的宏但必须在编译时开启了MPI和SCALAPACK支持否则并行跑大规模GW任务会非常吃力。我自己遇到过一个情况是集群管理员只编译了串行版VASPBSE计算时内存直接爆炸因为BSE激子哈密顿矩阵规模是(空带数 x k点数)的平方级别串行下根本装不下。硬件方面GWBSE和普通DFT对资源的需求差别很大。普通DFT结构优化32核加64GB内存就能跑得像模像样。到了GW步骤需要存储大量空带波函数和频率依赖的屏蔽函数通常建议至少128GB内存起步BSE步骤如果体系稍大比如超胞包含50个原子以上、k点6x6x6左右256GB内存是及格线对存储IO的吞吐要求也非常高。我见过不少人在工作站上跑GWBSE最后不是让CPU等计算而是让计算等IO所以强烈建议用NVMe SSD或者并行文件系统来存放WAVEDER和WAVECAR这类大文件。4. 完整实操流程从结构优化到BSE光谱4.1 第零步结构优化和基态DFT计算GWBSE是一个计算链不是单步任务。我在正式跑之前一定会先把前两步基础工作做扎实否则后面全白费。第一步是结构优化。对于晶体体系用PBE泛函做原子位置和晶胞参数的弛豫基本够用但需要注意加适当的k点密度和截断能。我的经验是体系的ENMAX来自POTCAR文件但结构优化时ENCUT要取ENMAX的1.2到1.3倍这样做是为了避免由于截断能不足导致的基态波函数人为畸变。晶胞优化要用ISIF3原子位置优化用ISIF2在GW之前晶胞参数最好用实验值或者更高精度的泛函校正过因为GW步骤对晶格常数的敏感性远大于普通DFT。第二步是基态静态计算这一步的输出是GW的输入基础。有几个关键点使用与后续GW一致的ENCUT不要中途改动采用较密的k点网格GW和BSE的k点必须比一般的静态计算更密通常建议至少保证单胞4x4x4以上必须打开LWAVE.TRUE.保存WAVECAR文件建议写入WAVEDER文件这个文件包含了占据态和空态之间的动量矩阵元GW步骤计算介电函数时离不开它如果你在旧版本VASP里跑过可能会注意到WAVEDER通常是和BSE强绑定实际上GW也需要它它描述的是能带间光学跃迁矩阵元在倒空间的变化行为。这个文件生成一次要花不少时间一定记得用LOPTICS.TRUE.提前生成并保留好。4.2 第一步GW0计算的核心参数基态算完之后进入GW0步骤。INCAR里最核心的设置是下面这些ALGOGW0表示采用GW0近似即自能计算中更新准粒子本征值但不自洽更新屏蔽作用W比全自洽GW便宜且数值稳定NELM200GW自洽循环的外层电子迭代步数上限GW0一般100到200步足够NBANDSGW对空带数量极其敏感需要从默认值往上加通常至少取“价带数导带数”的2到3倍具体需要在后续收敛性测试中确认NOMEGA50频率网格点数默认值往往不够尤其计算能量损失谱时建议增加到80左右很多人第一次跑GW0会踩同一个坑直接在静态计算的INCAR文件上改ALGOGW0但忘了检查NBANDS。DFT静态计算默认的空带数极少对于半导体大约只有十几个导带GW需要几十甚至上百个空带才能把屏蔽作用W描述到收敛。NBANDS不足时GW自能高频部分严重缺失带隙结果会偏低且不随空带数稳定。判断标准很简单把NBANDS翻倍如果带隙变化超过0.1 eV说明还没收敛继续加直到变化量小于0.05 eV为止。GW0跑完后OUTCAR里会输出每个k点、每个能带的准粒子修正量同时会在vasprun.xml里写入修正后的能带。我一般习惯把GW0和BSE分两步跑而不是用ALGOBSE一步到底这样可以在中间检查GW的准粒子带隙是否合理如果这一步已经严重偏离实验预期再往下跑BSE就没意义了先得回去检查初始波函数。4.3 第二步BSE计算核心参数GW0收敛后进入真正吃硬件的BSE步骤。保持已有的WAVECAR、WAVEDER文件不变INCAR改成BSE模式ALGOBSE让VASP进入BSE求解流程NBANDS这是BSE激发态空带数BSE需要比GW更多的空带因为它要建激子哈密顿量空带太少会把激子束缚在高能态压制吸收峰NBSE指定BSE实际使用的能带数量通常和NBANDS一致但可以从NBANDS里截取LOPTICS.TRUE.开启光学性质计算在BSE里必须配合WAVEDERNEDOS2000介电函数能量网格点数光谱想要平滑就设大点CSHIFT0.1展宽系数单位是eV用于给光谱加洛伦兹展宽取值太小光谱会像梳子一样离散太大则过度平滑模糊了峰位我一般用0.1到0.2BSE计算可调参数里很容易被忽略的是LSPECTRAL.TRUE.。当你设置ALGOBSE时打开LSPECTRAL可以让VASP先进行一次无相互作用的虚频介电函数扫描这样得到的无相互作用光谱可以当作参考线便于分析激子束缚能即相互作用光谱峰值和无相互作用光谱峰值的能量差。这个量在讨论激子效应强弱时非常有用强烈建议打开。4.4 文件保留与中间产物管理整套流程里会产生几个超大的文件WAVECAR动辄几十到数百GBWAVEDER也不小中途如果磁盘紧张很多人会手滑删掉WAVEDER。这个文件是BSE输入的必需品删了就必须重跑GW0重新生成代价极大。我的建议是给每个步骤建立独立目录保留上一层的WAVCAR和WAVEDER作为下一层的输入并做软链接而不是把文件复制来复制去。还有个经验是GW0和BSE都建议在同一个k点网格下运行。很少有人会注意到VASP对BSE计算中有无空带截断与k点对称性的要求BSE对k点的对称性处理会利用晶格点群如果k点网格非均匀或者用了太粗糙的网格BSE激子色散会被严重扭曲。从GW到BSEk点网格必须保持一致这一点比带宽收敛还重要。5. 关键参数表与收敛性测试策略5.1 一套可直接套用的INCAR模板下面给出一套我在实际项目里验证过、可以直接修改使用的INCAR参数模板。这里以单胞半导体体系为例超胞或分子体系需按需调整NBANDS和k点。# 基态静态计算 SYSTEM my_system ENCUT 1.3 * ENMAX ISMEAR 0 SIGMA 0.05 LREAL .FALSE. LWAVE .TRUE. LCHARG .TRUE. LOPTICS .TRUE. NEDOS 2000 # GW0步骤 ALGO GW0 NELM 200 NBANDS 128 NOMEGA 60 LSPECTRAL .FALSE. # BSE步骤 ALGO BSE NBANDS 160 NBSE 120 NEDOS 2000 CSHIFT 0.1 LOPTICS .TRUE. LSPECTRAL .TRUE.NBANDS128和NBSE120这个设定是针对大概20个原子以下的半导体单胞价带数量在10个左右。如果你的体系原子更多、电子更多128这个数远远不够需要按收敛测试结果来调整别直接复制套用。5.2 收敛性测试的正确顺序GWBSE的计算成本很高收敛性测试不可能每个参数都做一遍完整的矩阵扫描得有优先级。我自己的测试顺序是先固定k点和ENCUT扫描NBANDS。这是最重要的一步因为空带数对GW带隙和BSE光谱位置的影响最大NBANDS收敛后再测k点密度是否需要增加。这一步成本也很高但至少不用担心空带不够造成的假象最后再测NOMEGA和CSHIFT这类相对次要的参数实际操作中还有一个经常被忽略的步骤用较便宜的设置比如较少的k点、较大的展宽CSHIFT先跑一遍全流程确认输入文件和流程没毛病再上完整参数全面跑。我第一次跑BSE就吃了这个亏直接上全参数结果静算一步就卡了两天发现INCAR里忘了关掉结构优化相关设置白白浪费计算资源。5.3 不同体系的计算成本分级GWBSE的计算成本与体系尺寸、维度、k点数目关系很大这里给一个粗略的成本分级参考体系类型例子计算资源配置建议GW0耗时参考BSE耗时参考小单胞半导体硅、金刚石、GaAs单胞32核128GB内存数小时到1天数小时到2天中等双胞/三胞半导体2-3倍超胞、表面模型64核256GB内存1-3天2-5天有机分子晶体50原子以上分子晶体128核512GB内存以上3-7天5-15天强关联体系含稀土f电子氧化物建议先用DFTU或杂化泛函评估可能不收敛物理解释需谨慎这个表只是经验值实际耗时受硬件、版本和参数影响很大。但有一点是共通的如果你发现BSE比GW0还便宜很多多半是哪里出了问题——BSE要显式对角化一个巨大的激子哈密顿矩阵计算量只增不减。6. 常见报错与排查技巧实录6.1 报错速查表我在VASP的GWBSE实战中收集了不少报错和问题这里挑出频率最高、最影响后续运行的几个现象可能原因解决方案报错WAVEDER file not foundBSE步骤没找到跃迁矩阵元文件回到GW0步骤或静态计算打开LOPTICS.TRUE.重新生成WAVEDERGW0不收敛准粒子能级震荡初始波函数不适合GWNELM不够或空带不足加NELM到300检查NBANDS是否收敛考虑用HSE或DFTU提供初始波函数BSE内存溢出OOMNBSE设置过大k点过多或内存配置不足减少空带数减少k点网格增加计算节点内存检查SCALAPACK是否启用吸收光谱出现明显的负吸收或谱线不光滑NEDOS太少CSHIFT太小或LSPECTRAL与BSE结果混用NEDOS增到3000以上CSHIFT加到0.2试一下单独输出无相互作用光谱对照GW和BSE算出的带隙差异极大BSE的能带基础没对准可能是GW和BSE之间NBANDS不同GW和BSE使用一致的NBANDS上限和相同的KPOINTS文件输出光谱整体红移峰位低于实验0.3eV以上激子束缚能或带隙误差叠加常见于基态设置不佳先检查GW带隙是否合理再检查是否缺少激子效应导致的红移用LSPECTRAL光谱做参考线判断激子束缚能6.2 排查思路与心得遇到问题时我的第一反应不是改参数而是回到中间文件检视整个计算链。GW步骤看OUTCAR中的准粒子修正量是否合理特别关注VBM和CBM处的自能修正大小。如果VBM修正量异常地大超过2 eV基本可以断定是初始基态波函数问题或者空带不足。BSE步骤出问题则优先检查WAVEDER有没有正确加载、KPOINTS是否和GW一致。还有一个细节VASP不同版本之间GW和BSE的内部实现有差异。比如5.4.4版本里BSE的NBANDS不能超过NBANDSGW6.x版本则引入NBSE来显式控制BSE的能带选取。升级VASP版本后如果沿用旧INCAR轻则警告重则算错结果。我的习惯是每换一个新版本先用一个已知基准体系跑一遍全流程对比输出光谱和已知文献数据确认没问题再上实际项目。6.3 关于数值噪声与物理判据计算光谱很容易陷入“数值上跑出来了物理上不知道对不对”的状态。我自己的判据有两条一是GW带隙与实验值差距是否在0.1-0.2 eV以内二是BSE吸收峰与实验光谱峰位差距是否在0.3 eV以内。如果都满足说明参数收敛、物理合理如果GW已经准了但BSE峰位还是偏那大概率要审视激子束缚能讨论束缚能本质是相互作用谱和无相互作用谱的峰位差用LSPECTRAL.TRUE.输出两条曲线差值就是激子束缚能和实验激子结合能对比一下偏差太大就要警惕是不是屏蔽W算得不准、空带截断不够或者k点太稀。还有一个老生常谈但值得再提的点不要把CSHIFT黑影当成物理峰。展宽过大的光谱会把两个本应分离的激子峰糊成一个峰有时候还会把连续谱边缘抬起来被误认为新峰。对比不同CSHIFT取值下的峰位是否稳定是一个低成本的验证办法。7. 实操中的几个个人经验与扩展建议GWBSE这套流程我在实战里反复调、反复踩坑积累了几个百试不爽的小经验这里一并分享。第一个是关于WAVEDER和WAVECAR的保存策略。磁盘空间吃紧时很多人会在GW0结束后删除WAVEDER因为还没到BSE步骤等BSE步骤报缺文件才追悔莫及。我现在的做法是跑完基态静态计算后立刻把WAVEDER做一次压缩归档放在独立的存储目录然后建立软链接给后续步骤用。这套习惯救了我很多次尤其是中途有人误清临时目录的场景。第二个是关于BSE内存用量的预估。BSE激子哈密顿矩阵的存储规模大约是(NSEED^2)级别NSEED是BSE实际使用的价带空带组合数。粗略估算时可以这样算假如你设置NBSE120价带空带各选了一部分最终的矩阵规模可能在几十万乘几十万存储量轻轻松松超过几十GB。如果内存不足但不想砍k点还可以考虑把部分空带态用能带投影方法降维这个在VASP 6.3之后有更多的控制参数支持。第三个是后续扩展。GWBSE不只是能做吸收光谱还能扩展到圆二色谱、反射率、能量损失谱、介电函数实部虚部、甚至非线性光学。VASP的BSE模块输出的是介电函数张量只要把不同极化方向写清楚就能分析光学各向异性。我最近在算一个低对称性有机晶体的偏振吸收谱BSE出来的三组对角分量差异非常明显这类信息对实验上的偏振光吸收测量有直接参考价值。算完一套GWBSE我的原则是永远保留完整的输入输出档案和版本信息。这套计算非常贵算完不把数据整理好、不把中间过程记录清楚之后写论文或者复现结果时会非常痛苦。用表格记录每步的INCAR、版本、时间和关键结果长期来看是性价比最高的一笔投资。
返回列表