1、聊点八卦
Dalhousie在加拿大东海岸的“大”城市Halifax，街上冷冷清清，中国餐馆绝不超过三家。那里有加拿大最大的海洋博物馆（Maritime Museum of the Atlantic），里面有大量Titanic的珍贵文物。当然，那里还有加拿大海军屈指可数的盾舰，没准就是刚刚从博物馆里拉出来的。走在海边的步道上，一眼望去就是大西洋。大洋的澎湃绝不是在大连、青岛或者三亚的海边走走所能感受到的。在那附近还有 称加拿大最上镜小渔村的Piggy’s Cove以及最美灯塔。怎么说呢，可能我去的时候又刮风又下雨的确不容易让人体会到美感，或许是海鲜味道还不错就忽略了美景？
Dalhousie的确是个神奇的学校，按说Halifax这种放在中国只能算18线小县城的地方哪里容得下大神，有Becke在就是个奇迹，然而，他们还有Jeff Dahn。做锂电的小朋友一定知道这位大叔，Tesla给大叔投钱弄了个Tesla Lab，最近大叔说Tesla的电池寿命可以达到100万英里了。换句话说，车开废了电池还没事，还可以继续发光发热做Powerwall、Megapack。这里再多说一句加拿大这个神奇的国家，HSE里面的E则是在加拿大Université deMontréal (UdeM)的Matthias Ernzerhof，UdeM在蒙特利尔的位置大概相当于景山在北京的位置。如果再看看深度神经 络的大牛们都在哪的话，大家会发现加拿大虽然是个人口跟北京差不多的国家，但人家的科技成果还是挺耀眼的。
2、Becke-Johnson解决了什么问题？
这段主要聊聊背景，已经知道的朋友可以把这段也跳过去。大家从学DFT开始就总会听到有人在不停的念叨：DFT会低估半导体/绝缘体的能隙、低估能隙、低估能隙。能隙计算是个复杂的问题，这里只能草草说下。其实DFT不应该一个人背上低估带隙的锅。严格的DFT交换关联泛函应该是依赖于电子数变化所导致的阶跃变化的电荷密度的，本质上是不连续的，而我们通常把这玩意做成连续的了，就是所谓的derivative discontinuity问题。当然，这里还有自相互作用误差(self-interaction error, SIE)的锅，简而言之DFT存在过度离域化的占据态而Hartree-Fock (HF) 存在过度离域化空态的问题。这样的结果就是他们俩一个推高了占据态能量，而另一个推高了空态的能量，相当于一个低估、一个高估了能隙。实践中大家会用HSE泛函，引入了25%的HF exchange和DFT混合在一起来“对付”这个问题。想想这俩方法一个高估、一个低估能隙，混在一起就完事啦。这里的25%是从PBE0继承过来的，HSE的无屏蔽极限就是PBE0。PBE0里这个25%则是John Perdew大神（没错就是PBE的P）通过细致的观察和敏锐的洞察力给出来的。他把HF的占比写为1/n的形式，而n可以跟多体微扰修正阶数联系起来。对于大部分闭壳层分子而言，MP4（四阶多体微扰修正）已经足够准确了，就有了25%(0.25)这个系数。当然，实在算不对的时候也是有人调这个比例系数的，例如调成20%（相当于MP5
)。本质上，这种混合是让占据态和空态都有了接近的SIE问题，一起推高了占据态和空态，得到了相对准确的能隙。
这事看起来挺美好的，但现实总是那么残酷——HSE计算真的太费劲了。我的一位博士后做了个12个原子的HSE-SOC计算，程序两天才做了10步迭代，着实令人绝望。这里主要是HF的exact exchange要做遍历所有电子的两体积分，这个积分还是自旋和轨道依赖的。
当然，Slater老爷爷大刀阔斧地把轨道依赖扔掉了，提出了个exchange hole和exchange electron的图像，把轨道依赖的交换势用个局域的有效势代替了。于是，前面看着就吓人的轨道依赖的积分变成了下面这样的：
活生生少了一重积分，这事让人梦见都能笑醒了。慢着，先别急着笑，虽然这个有效势对于球对称的原子轨道是可以严格解的，但是Slater老爷爷并没有说这事情在凝聚态体系里到底怎么弄，这就是所谓的OEP问题(optimized effective potential for exchange)。在过去很长的一段日子里，这事情基本是一地鸡毛，反正解决的不怎么样，大多数的方法还是没绕开两体积分。
然而，奇迹出现了！重点重点，敲黑板敲到碎：2006年Becke大叔和Johnson小姐姐发现，其实真的不需要做两体积分，用一个特别简单的形式就能近似的给出交换泛函的形式，而且看起来还挺普适的。
这里的Tao_sigma是个动能密度，就是所有波函数梯度模的平方求和，Rho_sigma就是电荷密度啦，C_deltaV是个常数，他们写了篇论文扔去了JCP——又见JCP。这么个好东东发表了之后貌似没引起太大的动静，或许是大家一直在蠢蠢欲动，反正看起来故事似乎要结束了。
他具体是这么干，先定义了mBJ的交换势如下：
从形式上看跟原始BJ交换势还是很像的，这里面的BR交换势是Becke-Roussel交换势用来描述交换空穴的库伦势的。当然，这不是重点，重点是他们定义了一个跟体系电荷密度Rho和它的梯度相关的系数c。
这里的c其实是个无量纲数，由一个系数alpha = -0.012和一个系数beta = 1.023 bohr^1/2和一个与Rho相关的积分乘起来决定。这样一来是不是显得更有道理了，是不是普适性更强了，嘿嘿。
你们是不是想问我alpha和beta这两个这么奇怪的数是怎么来的。我想说，你们看了上面这张图，能不能把alpha和beta拟合出来？我想你们一定是可以的，Blaha大叔和Tran小哥也做到了，更重要的是他们把这玩意写到了程序里。虽然这部分程序不难写，但掌握着DFT程序资源的组就那么几家，而且个别组（例如VASP）已经进入了非常美妙的良性循环：全世界的用户帮他们试错、找bug、开发新功能。想想20多年前，VASP就是个在March Meeting会场免费送的小软盘，现在已经成了个庞然大物，还卖的挺贵。从这个角度看，咱们也应该多支持下国产的DFT程序，不然咱们连试错的机会都没有，发展就更难了。
我们再来看看两篇论文的引用情况—-。BJ的论文2006发表出来到2010年每年的引用还是个位数，到了2012年勉强到了20次，一路跌跌撞撞在去年达到了120次。而Blaha大叔的论文呢，2010年从22次起跳，一路飙升到去年的400次。这个事情告诉我们，一个好的方法很重要，但降低使用门槛、提高用户满意度对引用数这个KPI似乎更加重要。
3、VASP实践：这里我要由衷的感谢人民大学刘凯教授、麦吉尔大学孔祥华博士、天津大学胡智鑫教授、教育部的贾茜博士和新加坡国立大学的乔婧思博士。他们在将近十年前详细阅读、讨论了(m)BJ相关的原始文献，测试了大量泛函与mBJ组合的优劣和效果，优化了计算流程，这才有了下面的步骤1-3，以及注意事项1-3，而注意事项4是北京大学蒋鸿教授的贡献。
1) 用GGA/LDA/VDW-DF等做结构优化：这里具体想用什么泛函做结构优化其实不是那么重要，但对一些窄带隙的体系就不一定了。从我们之前在二维材料里的实践看，optB86b-vdW这个泛函得到的结构相关性质都非常好，这在黑磷、TMD、Te等等体系里都反复验证过。可如果用这个泛函得到的结构去做mBJ计算，最终得到的能隙却真是一言难尽，黑磷块体材料说好的0.3 eV的能隙被活生生算没有了(如下图)。测试了一堆泛函之后，孔博士和乔博士发现：如果用optB88-vdW弛豫得到的结构，再用mBJ去算个能隙，这种组合得到的结果(0.31 eV)似乎跟实验符合的挺不错的。这么看来optB88-vdW这个泛函还是很美好的，果然是Becke一家亲，当然，我们要”祈祷”他（们）不出现收敛问题。
2) 做个LDA/GGA总能计算，写下波函数：同样，这里用LDA还是GGA泛函其实并不是很重要，但是我们为了偷懒这里可能会做些难以描述的事情，例如把LDA扔一边，都用PBE。总之，这个总能计算的目的是为了输出波函数，用于做Tao_sigma的动能修正，所以需要把会在能带结构图里出现的每个k点的波函数都算出来。这就导致选取k点的时候有个权重的问题，既要保证MP k 格的优雅度，又要照顾到所有能带结构图里会出现的k点。于是，我们可以这么干(感谢乔博士准备的KPOINTS文件)：
前面是带有正常权重的均匀撒点k-mesh（优雅且符合处女座的审美），后面是所有能带结构图里要出现的且0权重的k点。你说这里是不是有重复的点，的确有的。但我想说的是：程序计算多花的时间大概率比你把重复的点挑出来再在画图的时候把重复点的数据拼回去的时间短。所以，不要纠结这么多了，虽然我们组的处女座程序猿表示这就是两个Python脚本的事
。
3) 读前面的波函数，做个自洽的mBJ计算：这里注意啊，是读波函数，不是电荷密度。K点的设置不要改，直接沿用就可以了。在INCAR里写上下面几个tag就好：
这里面GGA的标签是设置为CA(LDA)的，是否可以设置成其他泛函形式，例如PBE(PE)之类的，就需要大家自己去测试了。
4) 注意事项：
这里的CMBJ就是我们前面说依赖于体系电荷密度的和它的梯度的系数c。在三维体系的计算里，这个由程序自洽得到的c应该是可靠的。但是在低维体系里，程序自洽算出来的c就一言难尽了。我们的建议是，用低维体系的体相结构算出来的c，直接在计算里指定CMBJ = 1.xxx。
如果SCF收敛困难的话（要以平常心面对mBJ计算中的一两百步迭代），可以尝试换下mixing，例如用linear (straight) mixing (IMIX = 1，BMIX = 0.0001) 或者试试用其他算法，例如Damped (ALGO = D)或者CG (ALGO = C or A)。如果折腾一通之后的结果还是令人沮丧的话，请参考“4、未完待续”
。
大家也许都注意到了不论是BJ的原始形式还是mBJ的形式，我们只能看到交换势的形式但看不到交换能的形式。这意味着“(m)BJ的总能到底是个啥”看起来的确是个好问题，作为连带后果，大家就别想着用(m)BJ做弛豫了——当然，用它做弛豫是图什么呢？
同样，(m)BJ并没有很好的考虑derivative discontinuity (DD)的问题，它也会在描述f和d电子的时候遇到各种各样的困难，好在+U本质上考虑了DD的问题，则(m)BJ+U可以帮助(m)BJ更好地描述d和f电子。
4、未完待续：晚上9点半，我把这篇文章的草稿发给了几位朋友和同事，请他们审阅。北京大学蒋鸿教授反应敏捷，秒回我的消息，扔过来一篇JCP——再见神刊JCP。
这篇文章里，他系统比较了几十种材料能隙的自洽mBJ、微扰mBJ、PBE、mBJ+U@PBE+U、GW计算和实验测量值。这里的微扰mBJ是个看起来有点陌生的新名词，一会告诉大家这是什么意思。总之，比较了一通之后，他发现这个微扰mBJ对于小能隙体系的计算值比自洽mBJ更符合实验测量值；而到大能隙的时候，自洽mBJ才有了更好的表现（如下图）。那么大家关系的这个“大能隙”是多大呢，大概是6 eV，我想这大过了绝大多数我们关心的材料能隙。当然，我们这里说的只是微扰mBJ的计算值更接近实验值，并不代表这种方法更严格、更高阶、更优美。不论如何，用它来“凑”数的方法看起来的确无可挑剔。
好了，我们接着说这个微扰mBJ是什么意思。微扰是在谁的基础上微扰呢，没看错，就是PBE。蒋教授组里的孙怀洋博士大概率是被动辄200步的mBJ自洽迭代过程整惨了，这种看不到希望的心理感受的确不太好。于是，他索性不做mBJ自洽了，而是在一个PBE自洽计算基础上跑一步mBJ自洽迭代让程序算个c (CMBJ)，再把这个c做输入，基于PBE轨道做一步mBJ的严格对角化计算就完事了。具体可以参见蒋教授提供的如下PPT截图。
严格对角化要算多久的确不好说，但至少我们确定它只要算一步就完事了。跟动辄跑200步还不知道什么时候能结束的自洽迭代过程相比，孙博士和蒋教授的做法有着巨大的心理优势，想必也一定是有利于博士生们的身心健康的。
这里再做个广告啊，招博士后、博士后、博士后：后疫情时代的博士后需求
声明：本站部分文章及图片源自用户投稿，如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢！