致敬经典：一篇引用超过16万次的论文
1 一篇引用超过16万次的论文
科研论文的引用率，某种程度上是一篇论文影响力的一种反映，尤其是当一篇论文能在数十年里始终保持极高的引用。今天老司机为您介绍这样一篇论文，其引用高达16万次，就是1996年发表于Phys. Rev. Lett. 的论文。这是从事密度泛函理论计算的朋友们最常引用的一个工作，如下图：
这篇论文是纯理论推导，在已有的广义梯度近似（GGA）框架下提出了一种新的泛函。不同于当时已经建立起来的PW系列泛函（如PW86， PW91），这种新的泛函是非经验的（non-empirical），这就是后来被广泛使用的PBE泛函，其名称以三位作者名字命名。
PBE泛函一经推出，就迅速获得广泛使用，随后进一步发展为系列泛函，如RPBE，PBEsol, PBE0等，俗称PBE家族。如今使用广泛的杂化泛函HSE系列，也与PBE系列泛函密切相关。在1996年这篇PRL经典论文中，作者展示该泛函可以较为准确地计算出系列分子结构对应的原子化能，如下图所示：
从表中可以看到，PBE计算平均误差与PW91接近，显著优于同时期的LSD和UHF方法。
2 老司机致敬经典，重算H2的原子化能
老司机在2003年初接触密度泛函理论时，一位年长的引路人把这篇文章推荐给我，并让我重复上述表格中的数据，并鼓励我寻找新的方法，看看能不能讲表中数据计算的更为准确。20年后，我以相同的模式来训练计算老司机训练营的小白，也是鼓励他们用最常用的软件，试图复刻这篇经典。以下是老司机本人在工作中用Materials Studio中CASTEP，计算H2的原子化能的结果。
第一步：计算氢原子的能量
对于很多零基础学员，这一步就容易错，三个细节：A 单个H原子需要放入一个周期性盒子里，注意盒子不能太小；B 单个H原子实际上没有周期性，因此不需要设置Gamma点以外的K点；C 单个H是奇数电子，必然存在自旋极化（注意偶数电子的体系也可能存在极化），因此需要把自旋极化从一开始就打开。当然只有一个原子，不需要做结构优化，可以直接计算能量，如下图：
第二步：计算氢分子的结构和能量
要获得原子化能，必须准确优化H2分子并进一步获取能量。严谨起见，我们不默认H2是反铁磁的，而是采用测试的原则，即分别考虑两个H原子具有相同和相反的自旋，进行了优化和计算，如下图对应的反铁磁(+1, -1)和铁磁（+1，+1）两种初始自旋设置。结果是前者能优化出与实验接近的键长，并且能量显著低于铁磁情形，所以铁磁构型不是基态，反铁磁才是；
对于是否测试磁矩，我学习之初是有困惑的，也反问了教我计算的老师。他告诉我：测试是一种习惯，想当然则是计算人员的大敌，不是每一个体系我们都提前知道答案，科学研究毕竟是开放且充满未知的。我谨记在心；
基于测试结果，计算获得了原子化能，数值为：104.46 kca/mol, 与PRL原文数值105 kcal/mol极为接近，略低于实验结果（109kcal/mol）。二十年前，当我得到这个结果时，整整一个晚上兴奋地没睡着，盘算着第二天把余下的其它分子全部重复一遍！
3 致敬经典的最佳方式是超越
当我因为重复出这篇经典论文的H2数据时，有些沾沾自喜，而我的引路人问我：你能算得更准确一些吗？你觉得还有哪些地方可以提高？我瞬间有些羞愧，赶紧重新思考，后来采用量化方法得到了更准确结果。
就在最近，我指导训练营的零基础学员，依然采用相同的软件，如VASP和CASTEP等大众软件，如何进一步提高。一个性价比极高的处理手段就是改用更准确的MetaGGA泛函，如SCAN泛函（对相关泛函和密度泛函理论不熟悉的同学，请移步观看老司机的公开课：https://www.bilibili.com/video/BV11B4y1C7Fk/?vd_source=aca19dff321fb0979ec90bf26485b141）。
SCAN泛函在CASTEP中是经过升级的版本，即RSCAN，计算结果如下：原子化能达到106.24 kcal/mol, 提高了1.78 kal/mol, 更接近实验测量值（109 kcal/mol）。
后记1 ：大约6年前，我偶然得知当年教我理论计算的引路人，已经不再从事计算模拟，而是完全转向实验研究。
后记2：需要上述测试PPT的同学，请加微(jisuanlaosiji)注明“H2原子化能”，索取资料，如果需要老司机进一步解说PPT中的测试参数，请在文后留言。
