|
谈谈重复不出来计算化学文献里的数据的可能原因 On the possible reasons why data in computational chemistry literature cannot be reproduced
文/Sobereva@北京科音 2023-Aug-1
0 前言
经常有做计算化学的人在网上问为什么他没法重复出来文献里的计算数据。导致重复不出来的原因实在太多了，然而初学者提问此类问题时几乎总是含糊其辞，提供的有效信息甚少，导致别人根本无从回复，或者需要别人罗列一大堆可能导致差异的原因。我感到每次回复此类问题太费事，于是在这里专门写个文章说说这事。遇到这种问题的读者仔细看完此文后，应该能料想出可能是什么原因导致自己和文献的结果存在差异。就算还搞不明白，至少也应该能明白在网上问别人的时候应当提供什么信息，知道怎么把导致差异的可能原因的范围尽可能缩小，免得别人打很多字罗列诸多当前不需要考虑的可能性，或者别人嫌你提供的有效信息完全不够而根本不想回复。
这里首先要强调一点的是，如果你的目的就是为了精确重复文献数据，那你就要尽最大可能使用和文献完全相同的方式去做计算，而不要管人家的做法合不合理，哪怕你明知道人家的计算方式不当。下文说的内容都是为了让你能尽可能准确重复出文献里的数据，而不是讲怎么合理地计算。如果你的实际目的是想获得相同体系相同研究问题的合理的结果，那你只要确保你自己的计算是合理的就行了，没必要非得把别人的数据重复出来，说不定人家的计算本来就有问题。如果最终验证发现别人的计算有硬伤并因此导致结论错误，或者凭基本知识就能断定别人的计算有问题，那么可以发一篇comment文章指正以免以讹传讹，比如《我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment》（http://sobereva.com/584）里介绍的文章。
1 结构不一样
自己计算用的结构和文献不一致，是导致数据无法重现的最常见原因之一。如果你用的结构就是文献正文或补充材料里给的坐标，那就可以直接排除这个因素了。注意有些文献的SI里键长、笛卡尔坐标单位是Bohr而不是埃，千万别弄错了。
如果你计算的基本结构特征都和文献里不符，也即算的都不是同一个东西，能重现出文献里的数据就怪了。很多体系有不同的形态，比如酚酞在酸性、中性、碱性环境下主要的存在形态明显不一样，氨基酸在真空和水环境下的质子化态明显不一样，姜黄素在水和有机溶剂中一种是烯醇式一种是酮式结构，等等。你必须确保你算的和文献里算的确实是同一种状态的结构，如果文献里给了截图，一定要仔细对照。
做计算要有基本的化学常识，并且要认真仔细，免得搭出来没意义的结构去做计算，自然和文献里对不上。比如文献里给出Lewis结构式时通常是隐去碳上的氢的，如果计算的时候真不画氢，或者少画了个别的氢，自然计算毫无意义。再比如，文献里算一个羧基配位的过渡金属配合物，如果你在计算时羧基上还留着氢，明显也是错的。再比如，文献里画了一个带一个小配体的过渡金属配合物的Lewis式，若你算的是水溶液中的状态并忽略了配位水分子，也自然是错的。再比如，文献里给出一个联苯或者碗烯的二维结构式，看上去是在一个平面上，若你去计算其基态时真的摆成平面来优化，得到的也是大概率是错误的结构（有破坏平面的虚频的结构）。还要注意一些结构的构建不要想当然，比如SiO2表面在水中时表面的O主要是以羟基及其去质子化形态存在的，比例和pH有关，如果都当成是乌秃秃的氧就错了。结构构建的关键性信息在文献里一般都会有具体说明，但也可能一些文献中在描述上有疏漏。
有些体系在计算的时候结构构建方式不唯一，比如《要善用簇模型，不要盲目用ONIOM算蛋白质-小分子相互作用问题》（http://sobereva.com/597）、《18碳环（cyclo[18]carbon）与石墨烯的相互作用：基于簇模型的研究一例》（http://sobereva.com/615）和《使用量子化学程序基于簇模型计算金属表面吸附问题》（http://sobereva.com/540）里提到的簇模型，怎么构建簇明显会影响结果，特别是当簇用得较小的时候。若要重复文献里的数据，就要用尽可能相同的模型，哪怕其模型很昂贵。
很多分子有不止一种构象（势能面极小点），柔性大分子甚至有一万种以上的构象。几何结构优化一般会优化到与初始结构最接近的势能面极小点。而不同构象下计算的属性可能有不小差异，因此你必须确保你用的构象和文献里一致，也即结构优化使用的初猜结构要尽可能像文献里给出的。如果文献没明确给出三维结构信息，按理说文献应当用的是能量最低构象（确切来说是自由能最低构象）做的计算，但也有可能作者忽略了构象搜索并得到了不合理的结果。如果你对构象搜索一点概念都没有，建议看下文了解些常识
《使用molclus程序做团簇构型搜索和分子构象搜索》（http://bbs.keinsci.com/thread-577-1-1.html）
《gentor：扫描方式做分子构象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html）
《将Confab或Frog2与Molclus联用对有机体系做构象搜索》（http://bbs.keinsci.com/thread-20063-1-1.html）
《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）
也不是所有情况都适合用能量最低构象来计算，比如某化学反应可能是从某个能量相对较高的构象发生的，可以对反应的IRC的反应物一端做几何优化得到之。
团簇体系的构型搜索和柔性分子的构象搜索同等重要。哪怕水二聚体这么简单的两个小分子构成的团簇都都有不同的极小点构型，并且对应的结合能明显不同。因此计算这种体系也必须确保和文献用的构型相一致。如果文中没明确交代，一般默认其用的是能量最低构型结构做的计算，但也可能作者忽略了这个问题，导致误用了能量不是最低的构型并因此得到了错误的结果。如果你对构型搜索没有概念，建议看此文：《genmer：生成团簇初始构型结合molclus做团簇结构搜索的超便捷工具》（http://bbs.keinsci.com/thread-2369-1-1.html）。
有的时候两种极小点构型或构象可能结构相差很小。这种情况需要把可视化程序里的视角和文献里的图片摆得尽可能一致再去观察，免得算的其实不是同一个结构。
上面说的是计算极小点结构的情况，自不必说，重复文献里过渡态的计算结果时也必须严格确保和文献里用的结构一致。
对于研究固体表面的情况，要注意文献里的结构是怎么来的。表面是从体相直接切出来的未经优化的（unrelaxed），还是经过优化的（relaxed），还是使用的考虑表面重构的（reconstructed），这需要区分清楚。文献中计算用的是具体哪个晶面也必须要搞清楚。固体表面吸附问题也有类似构象/构型搜索的势能面存在多种极小点问题，一个小分子可能在表面上不同位点以不同方式吸附，实际算的结构必须和文献相一致。文献里算的吸附是物理吸附还是化学吸附也得搞清楚。对于吸附过程有势垒的情况，初猜结构中你把被吸附物和表面摆得相距很近，一般优化完了的是形成新键的化学吸附，如果摆得较远，一般得到的是弱相互作用维持的物理吸附，二者的结构、相互作用能相差甚大。研究表面吸附还有覆盖率的问题，要注意文献中对所用结构的相关描述，研究低覆盖率和饱和覆盖率时用的模型明显不同。
用slab模型研究固体表面还要注意厚度问题。slab太薄的话，得到的功函数、吸附能等性质都不准确。对于二维层状材料的研究，还要注意层数对计算的性质的影响。这些方面都应当和文献严格一致。对表面类体系优化时还经常牵扯到对一部分原子冻结使之固定为体相的结构，之后振动分析也相应地牵扯冻结，冻结方式必须和文献一致。
做某些体系的分子动力学模拟也同样要注意初始结构的差异，否则现象可能和文献不符。比如有文献研究表明水中的蛋白质由于疏水效应会自发进入口径足够大的碳纳米管，如果你一开始把蛋白质放得离碳纳米管老远，由于质量颇大的蛋白质本身整体运动就较慢，再加上碳纳米管的随意旋转（没冻结的话），可能文献里说的那种蛋白质自发进入管内的现象跑挺长时间也出现不了。再比如想模拟结冰过程，哪怕一开始把温度设得显著低于冰点，但没有在水中加入局部的冰的结构的话，跑很长时间也不会看到冰的结构的出现。至于在模拟的时间尺度内一定会出现的现象，初始结构就不重要了。比如水和乙醇的互溶，无论初始结构里把二者混合均匀，还是分成两相，最后肯定都是互溶的。
做周期性体系的模拟还要注意盒子尺寸与文献要一致。用第一性原理计算固体时，假设k点数是固定的，盒子尺寸的不同会影响精度。做动力学模拟时，盒子尺寸还间接影响原子运动行为以及模拟得到的属性。
用Quantum ESPRESSO等平面波程序时，以及CP2K这种将平面波当辅助基函数的程序时，如果算的是某些维度为非周期性体系，还需要注意真空层的厚度，要尽可能和文献一致。因为真空层可能对结果有不可忽视的影响，比如影响非周期性方向上与相邻镜像的作用是否可以忽略。而且比如CP2K里用MT的Psolver时还有非周期性方向上盒子尺寸大于相应方向有电子密度明显分布的两倍以上的要求，显然真空区不能小了。
假设文献用的是实验测定的坐标直接进行性质的计算，还要注意实验来源必须一致。比如文献或晶体数据库里对同一种晶体结构可能提供了不同温度下测定的坐标，会有一定程度的不同（体现热膨胀）。显然压力也直接影响测定的结构。哪怕是相同的测试环境，不同来源的晶体结构也会有差异，特别是氢的位置的确定，无疑这对氢键等问题的计算结果会有直接影响。做量子化学/第一性原理研究的人一般都会至少对氢的坐标重新进行优化来refine之，需要注意文献里是怎么考虑的。蛋白质的高解析度和低解析度测出来的晶体结构里残基侧链的位置也可能相差甚大，甚至明显影响到做酶催化、分子对接之类问题的结果。同样的属性，以不同类型实验测定的结果注定也会有差异，比如气相电子衍射和单晶衍射测的一个是气相结构一个是晶体结构，对于柔性分子差异可能很大，这点参考《实验测定分子结构的方法以及将实验结构用于量子化学计算需要注意的问题》（http://sobereva.com/569）。
上面列举了诸多结构层面影响计算结果的可能性，当然不可能列举得很完备，实际情况多种多样。总之希望读者能意识到严格保证和文献里用的结构的一致对于重现数据有多么重要。
2 计算设置不一样
为了精确重复文献数据，要尽可能仔细阅读文章及其补充材料，了解清楚文献计算细节，使得自己的计算设置能最大程度与文献相一致。本节里谈到的很多可能导致差异的因素在《谈谈不同量子化学程序计算结果的差异问题》（http://sobereva.com/573）里都有详细说明，务必注意阅读。
量子化学和第一性原理计算用的理论方法，以及理论方法中涉及的设定和参数（如范围分离泛函的w参数、色散校正的经验参数、CASSCF和多参考方法的活性空间的选择、TDDFT计算是否用了TDA近似、DFT+U参数、CASPT2计算用的IPEA shift参数、DLPNO方法用的阈值、GW@DFT和DFT-SAPT具体基于什么泛函、sTDA计算用的经验参数，等等）必须和文献精确一致。有的理论方法在同一个程序里甚至有两种定义，比如B3LYP在ORCA里写成B3LYP还是B3LYP/G是不同的，要搞清楚文献实际用的哪个（通常是前者）。还要注意理论方法名和关键词的差异，比如里这里提及的：《Gaussian16里用PBE0关键词计算的实际上是PBE0-DH双杂化泛函》（http://bbs.keinsci.com/thread-13660-1-1.html）。有的文献的写法还可能存在含糊性，比如可能文献就说对某个开壳层体系用了二阶微扰理论，但实际上二阶微扰理论用于开壳层的实现形式有一堆，比如UMP2、PUMP2、ROMP2、ZAPT2、ROMP、OPT2等等，如果文中给了引用的文献，要根据文献判断具体是哪种。
对计算结果可能有明显影响的各种条件必须一致，比如：
• 溶剂效应的考虑。如果用的是隐式溶剂模型，那么用的具体是什么方法？描述的是哪种溶剂？用CP2K的SCCS模型时α、β、γ参数是怎么设的？描述混合溶剂的时是怎么确定溶剂参数的？如果用的是显式溶剂模型，那么溶剂分子排布是怎么考虑的？有没有同时用隐式溶剂模型？
• 计算热力学校正量的模型。比如《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）里介绍的Shermo程序默认用的以及ORCA用的是Grimme的QRRHO模型，而Gaussian用的是较原始的不适合含低频体系的RRHO模型
• 计算热力学量时是否考虑了平动和转动，显然算周期性体系时不应当考虑（在Shermo程序中这由settings.ini里的imode决定）
• 计算转动对热力学量贡献时对转动对称数的考虑，看Shermo手册附录部分了解相关常识
• 计算热力学校正量、振动分析用的原子质量
• k点数目和分布
• 平面波相关计算的截断能
• CP2K中，做杂化泛函时用的库仑截断半径、DFT+U计算用的布居计算方法以及考虑的角动量、Poisson solver的选择
• 对相对论效应的考虑
• 对外场/外势的考虑
• DFT积分格点精度。参看《密度泛函计算中的格点积分方法》（http://sobereva.com/69）
• 是否用了冻核、怎么冻的
• 双电子积分的积分屏蔽设置
• 是否用了RI、DLPNO等加速计算方式
• 是否考虑了偶极校正
• QM/MM计算划分的位置、MM与QM区之间的相互作用的考虑、link原子的设置
等等等等
基组必须和文献严格一致，并注意上述http://sobereva.com/573博文里说的许多细节问题。注意一个基组/赝势名可能有不止一个指代，比如Stuttgart赝势，有不同方式考虑相对论效应的版本、不同价电子数和氧化态的版本，Gaussian内置的和Stuttgart官网上的赝势和赝势基组定义往往又有所不同。再比如，GAMESS-US用户写CCD和CCT关键词分别用的是cc-pV(D+d)Z和cc-pV(T+d)Z，而不是更常见的cc-pVDZ和cc-pVTZ。要弄清楚作者实际用的什么，仔细看文中的描述，并注意作者引的基组原文。
不光一般意义上的基组（叫主基组或者叫轨道基组）要和文献计算时的一致，辅助基组也应当一致。辅助基组种类很多，比如RIJ用的/J辅助基组，RIJK用的/JK辅助基组，电子相关计算用的/C辅助基组，CP2K里加快杂化泛函计算用的ADMM辅助基组，CP2K里做LRIGPW计算用的辅助基组，等等。同一个主基组能够或者适合搭配的某种目的的辅助基组可能不止一种，而且有的时候有含糊性和任意性，比如ORCA里用ma-def2-TZVP主基组时，可能有人用autoaux关键词自动产生辅助基组，可能有的人用标准的def2/J辅助基组。
对于基于经典力场的模拟，力场参数显然必须和文献严格一致。如果参数是从文献里拷来然后用在自己计算中的，注意单位转换，以及1-4相互作用规则和LJ参数的混合规则等细节必须和文章相同。文献给出的LJ参数有的是R0有的是σ需要区分。原子电荷也必须和文献保持一致，如果原子电荷是作者自己算的，要搞清楚计算原子电荷的方法，并且如果是基于波函数算的，那么波函数是在什么几何结构、溶剂环境、计算级别下得到的要弄清楚。这点很重要，比如对于RESP电荷，气相下计算的和水环境下计算的，Hartree-Fock和B3LYP计算的，相差大了去了，参看这些讨论：《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》（http://sobereva.com/441）、《RESP2原子电荷的思想以及在Multiwfn中的计算》（http://sobereva.com/531）。
做动力学模拟的情况，必须确保模拟相关设置和文献一致。比如步长、步数、坐标/速度/受力/能量保存频率、参考温度和压力、控温和控压方法、热浴和压浴的时间常数、可压缩系数、各向同性还是各向异性控压、初速度怎么产生的、范德华和静电作用的计算方式（如简单截断、twin-range、Ewald、PME、SPME等，以及其中涉及的截断半径，是否用了switching function等处理）、水模型是刚性还是柔性的、是否用了键/键角约束、用没用能量/压力色散校正，等等。
如果你用的计算程序和文献里不同，或者和文献里的相同但版本不同，那么由于默认设置不同，以及算法实现的不同，可能对结果带来很大差异，一定要仔细阅读《谈谈不同量子化学程序计算结果的差异问题》（http://sobereva.com/573）。
在量子化学波函数分析领域也有很多细节设置会影响分析结果。这里以非常流行的波函数分析程序Multiwfn（http://sobereva.com/multiwfn）为例：
• 按照《使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析》（http://sobereva.com/179）、《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》（http://sobereva.com/159）、《使用Multiwfn对静电势和范德华势做拓扑分析精确得到极小点位置和数值》（http://sobereva.com/645）等文章介绍的方法进行分析时，格点间距数值越小精度越高，但也越耗时
• 按照《使用Multiwfn绘制态密度(DOS)图考察电子结构》（http://sobereva.com/482）绘制PDOS和OPDOS图时，在Multiwfn里选择的计算轨道成份的方法会直接影响结果
• 按照《使用Multiwfn模拟扫描隧道显微镜(STM)图像》（http://sobereva.com/549）和《使用Multiwfn结合CP2K的波函数模拟周期性体系的隧道扫描显微镜（STM）图像》（http://sobereva.com/671）绘制STM图时，偏压的数值和常高模式扫描的平面的选择都会明显影响结果
• 计算《衡量芳香性的方法以及在Multiwfn中的计算》（http://sobereva.com/176）和《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》（http://sobereva.com/138）里介绍的多中心键级时，基于原本的基函数和基于自然原子轨道（NAO）计算会得到不同的结果
• 计算Hirshfeld原子电荷时需要用到原子自由状态下的密度。用自行提供的原子波函数文件来计算此密度，还是直接用Multiwfn内置的此密度，会有所不同。
本节提到的这些计算细节是理应充分体现在文章的computational details里或者补充材料里的。但由于一些文献的作者为了省篇幅，或者用的很多设置和程序默认的一致，或者由于作者计算水平比较外行、意识不到很多细节对计算非常重要，就在文中没特意说。如果有些要点不知道作者具体怎么考虑的，自己把各种可能性都试试就知道了，文章没特意说的设置就姑且用程序默认的。
3 计算的电子结构状态不一样
净电荷和自旋多重度的设置必须和文献里的相同，否则相当于算的是其它电子态，自然算出来的各种性质都不同。比如计算基态的氧气发生的反应、计算富勒烯的原子化能，如果氧气和碳原子误当成了单重态来算，显然结果是错的。此外，稍微了解点两态反应性的概念就会知道，不少体系在不同自旋态下反应机理都明显不一样，如果自旋多重度或净电荷设错了，自然反应物/产物/中间体/过渡态等无法得到和文献里一致的结果，而且相差悬殊。
如果文献没注明、从文献中看不出相关信息，那么净电荷一般为0。如果被研究的体系的自旋多重度不那么显而易见，或者研究中涉及到了不止一个自旋态，文中一般都会说明自旋多重度。如果确实没说，那就得查阅文献或者数据库，或者自己判断、测试了。
很容易被忽略的是开壳层单重态或者叫自旋极化单重态。用KS-DFT等方法来计算的话应当做对称破缺计算，参看《谈谈片段组合波函数与自旋极化单重态》（http://sobereva.com/82）。比如整体为单重态的反铁磁性耦合体系、双自由基体系、共价键被拉长的均裂状态就是最典型的例子。还有一些特殊的体系也是如此，比如18碳环的D18h点群的过渡态结构、乙烯二面角扭转成90度的结构、六并苯。重复文献里的这些体系的计算时别误当成了闭壳层单重态来计算（文中如果也是做了对称破缺计算，一般都会明确说。如果没说，没准文献也误当成了闭壳层计算）。
对于第一性原理计算磁性材料时必须谨慎，很多情况下不是光设个整体的自旋多重度就完事的，而必须恰当设置初始磁矩才能收敛到真正要考察的磁性状态。为了重复文献（或材料数据库）中的数据，若它们给了收敛的波函数对应的各个原子的自旋磁矩/自旋布居信息的话，最好根据其来设置初猜波函数里的原子自旋磁矩，以尽最大可能收敛到和文献相同的波函数状态。
要注意波函数稳定性问题。即便净电荷、自旋多重度的设置是正确的，也不代表SCF收敛后得到的波函数就是实际想研究的态、和文献里算的状态相一致。如果文献里研究的是基态，你算的结果和文献对不上，而且体系又不是简单有机体系而是电子结构特殊的体系，那么应当做一下波函数稳定性检测，比如Gaussian里可以用stable关键词，发现不稳定的话可以用stable=opt尝试搞出稳定波函数。不稳定的波函数可能意味着这是某个激发态，还可能意味着此波函数是没有任何物理意义的虚假的状态。
4 计算的东西或方式不一样
重复不出文献的数据时还要想清楚到底你算的东西、计算的方式是否和文献一致、有可比性。下面列举一些典型的需要注意的情况。
比如通过能量差E(AB)-E(A)-E(B)考察俩分子A和B间相互作用强度，复合物结构肯定是要优化的，而单分子是直接从复合物优化后的结构里直接取出来，还是也做优化，这在不同文中的做法就不一样了。单分子不优化情况下得到的能量差一般称为（复合物结构下的）相互作用能，单分子也优化情况下得到的能量差一般称为结合能，二者之间相差各个单体的变形能之和。对于诸如《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（http://sobereva.com/513）里面列举的那些刚性小分子，相互作用能和结合能差异甚微，但如果考察的是两个容易变形的分子间的作用强度，比如两个氨基酸分子的结合，那么相互作用能和结合能就必须明确区分了。J. Chem. Phys., 158, 244106 (2023)还专门讨论了什么情况变形能会较大。对于键能，也是类似地有不同计算方式。所以研究此类问题的时候一定要搞清楚文献具体是怎么算的。
再比如做SAPT能量分解计算，会得到很多耦合项，比如交换与色散作用的耦合项，这是归到色散项里，还是归到交换互斥项里，不同文章的做法可能不一样，需要注意文章怎么说的。再比如计算反应的自由能问题，是计算标准反应自由能，还是计算特定浓度下的反应自由能，这是不一样，相差浓度校正项，要搞清楚文献实际算的是什么。再比如电子亲和能、电离能、激发能，都有垂直和绝热之分，文献里计算哪种都有，也必须确保和文献算的东西对应。
文献里有时候对热力学量的标注比较含糊，要分清楚文献算的是什么温度下什么量。比如E，文献里可能指电子能量，也可能指内能，需要结合语境判断。
模拟ECD谱需要转子强度，有长度形式和速度形式两种，诸如Gaussian都会输出，在Multiwfn里按照《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》（http://sobereva.com/224）里说的绘制ECD时可以选采用哪种，一般推荐用速度形式。如果你用的和文献不一致，会导致峰型存在一定差异。
基组的CBS外推的方法不唯一，有不同公式，有的情况外推方式里也涉及到参数，参看《谈谈能量的基组外推》（http://sobereva.com/172）。不光是能量外推，还有比如GW计算的准粒子能量随基组的外推、几何结构随基组的外推等。一些文献对高级别理论方法结合大基组还用一些近似方式估计，比如一种常见的是CCSD(T)/大基组 = CCSD(T)/小基组 + (MP2/大基组-MP2/小基组)，这在《各种后HF方法精度简单横测》（http://sobereva.com/378）里专门说了。涉及到这些高精度数据估算时，需注意保证和文献的做法一致。
还有一些量本来就没有唯一定义，去重复文献的数据更是得弄清楚文献里用什么定义、什么方法算的。比如自由体积，我在《使用Multiwfn计算晶体结构中自由区域的体积、图形化展现自由区域》（http://sobereva.com/617）里介绍了一种合理的计算方式，但显然这不是唯一方式。再比如分子半径，我在《谈谈分子半径的计算和分子形状的描述》（http://sobereva.com/190）就已经列举了多种计算方式。还有诸如原子电荷，计算方法更是多了去了，如ADCH、RESP、Hirshfeld、NPA等，简要介绍见《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）。
为了正确理解文献具体算的是什么，必须把背景知识搞清楚，典型的就是gap这个概念，有HOMO-LUMO gap、band gap、optical gap、fundamental gap等，参见《正确地认识分子的能隙(gap)、HOMO和LUMO》（http://sobereva.com/543）。如果没搞清楚概念和定义方式，可能重复不出文献的数据，或者实际上文献的计算方式或用词不当，你却意识不到这一点。
某些问题的计算中可能涉及到一些实验值，要搞清楚作者用的什么值。比如用量子化学程序计算含碳物质的生成焓，通过构造热力学循环可知其中需要利用实验的（或第一性原理计算的）石墨的升华焓，要搞清楚这个量是从哪来的，具体是多少。再比如， 量子化学计算某分子的标准氧化还原势，其中利用到标准氢电极的电势，其数值有不同来源，从4.05到4.44 V都有报道，Phys. Chem. Chem. Phys., 16, 15068 (2014)建议对理论计算的情况用4.28 V。
还有一些研究涉及到一些额外的参数需要注意，典型的就是频率校正因子，见《谈谈谐振频率校正因子》（http://sobereva.com/221）。是否使用频率校正因子，以及对某一特定级别使用什么目的的、哪个来源的校正因子，都会影响结果。
周期性体系的计算中，轨道能量的绝对数值是没意义的，不同程序用的能级零点也不一样，需要关心的只有相对数值，比如带隙。所以不要问为什么DOS图的形状和文献一致但横坐标和文献不一致。这也正是为什么文献在绘制DOS图时往往一般把价带顶或者费米能级位置设为0点，因为这样才便于公平对比。
还有些量在文献里用的习俗和你用的可能不同。比如算结合能问题，有少数文章把结合过程的能量降低量作为正值记录。再比如NBO给出的E2超共轭作用能，程序输出时把正值作为超共轭作用造成的体系能量降低，但有的文献报道E2值的时候取的是负值。
如以上列举的情况可见，算某个东西的时候必须仔细阅读文献搞清楚作者是怎么算的、用的什么定义、算的到底是什么。
5 数据的可重现性本身就差
有些数据很容易重现，比如在RHF/6-31G**级别下优化水分子的基态结构，不管用什么主流量化程序、怎么优化，得到的都是相同结果，初始结构很随意。有些数据本身就是难以精确重现，典型的就是凝聚相体系的分子动力学。由于分子动力学过程有混沌效应，而且实际计算中往往不可避免地引入随机因素，会导致文献中的动力学模拟的定量结果几乎无法严格精确重现。关于这点我专门写过文章，参看《数值误差对计算化学结果重现性的影响》（http://sobereva.com/88）。更具体来说也看重现文献里的什么动力学模拟数据。比如小分子的液态性质相对来说是容易重复的。比如重复密度值，文献里是801.92 g/L，你算出来801.88 g/L，这就可以认为是很好重复出来了。如果你算出来是810.21 g/L，而且模拟流程和设置是合理的，那就是有明显不可忽视的差异了，已经超过了随机性的影响了，需要根据前面所说的来考虑导致与文献存在差异的可能原因。也有的量的重现相对比较困难，比如磷脂膜里插入了一个分子，去计算它的侧向扩散系数，这就属于相对比较难算准的，因为本身分子在磷脂膜这种拥挤的环境里扩散就困难，而且体系里就这么一个分子，没法通过同类分子的平均来减小统计误差。就算你自己跑的时间非常长以保证统计精度绝对足够高，但文献里给出的数据未必统计精度就够高。像这种问题，对文献里的数据的重现精度就别要求太高了。
还有一种情况是动力学模拟的现象在有必然性的同时本身也有一定随机性，跑出什么现象有运气成分。比如距离蛋白质的口袋一定距离处放一个小分子配体，在水环境下做动力学模拟，考察2 ns模拟时间内配体能否和蛋白质结合，初速度随机生成。像这种模拟的结果就很有随机性，如果跑10次模拟，可能有的模拟刚过几百ps小分子就进入口袋了，也可能有的模拟开始后小分子恰好往蛋白质相反方向游离，跑到2 ns还没进入口袋。对于这种问题的研究必须做大量的模拟来得到统计结果，较严谨的文献普遍都会这么考察这种现象明显有随机性的问题，你若要重复文献也必须以相同的方式去模拟。如果模拟的样本数不够多的话，也不要指望结果特别相符。比如前述问题，可能文献里跑10次发现有5次进入口袋，你模拟10次发现有3次进入口袋，这未必是模拟设置和文献有区别，而纯粹是随机性所致，要以正确的态度来看待这种情况。
顺带一提，量子化学、第一性原理做静态的计算的可重复性比起分子动力学模拟要好得多得多，但也不排除极个别情况前者可能受到随机性的明显影响。比如对势能面很复杂的大分子进行优化，特别是QM/MM计算时可能牵扯到几千原子的情况的优化，有可能由于微小的随机性因素（如使用了并行计算所带来的）导致两次优化分别收敛到了势能面上可能相距挺远的两个极小点，且有肉眼可察觉的差异。
6 对相符程度要求过高
有人其实可能已经合理重复出了文献中的数据，但由于他对重复的精度要求过高，导致误以为没重复出来。
头脑里要有收敛限的概念。比如文献里优化的结构二面角是123.64度，你优化出来的是123.69，这就可以认为是重现出来了，因为这点差异通常都已经小于计算程序的几何优化的默认收敛限了，差异基本来自初始结构的任意性。再比如，文献里DFT计算出的反应能是-8.12 kJ/mol，你算出来的是-8.18 kJ/mol，也完全可以接受了。如果本来你用的计算程序和文献里都不同，在尽可能令设置和文献相同的情况下，你算出来比如是-8.32 kJ/mol也可以算重复出了文献，各种鸡毛蒜皮不值得一提的因素足矣带来这种程度的差异。但如果你算出来是比如-9.8 kJ/mol，那就明显不是不值得一提的差异了。再比如TDDFT计算电子激发能、DFT计算HOMO能级，和文献里相差0.02 eV，这都可以算一致，但如果相差到0.1 eV，这就不算一致了，毕竟TDDFT算激发能的精度都在0.3 eV左右（泛函选择得当的前提下），当差异达到接近方法本身的误差的数量级时明显就算显著差异了。
重复文献数据的时候用的收敛限应当足够严，起码令收敛限显著小于你期望的对文献数据的重现精度。虽然文献作者用的收敛限未必够严，因此光是你把收敛限设严未必有意义，但至少先做自己这里能做的事。
7 自己的计算存在硬伤
重复不出文献数据有可能是自己计算存在硬伤导致的。初学者缺乏常识、经验和敏锐性，计算很容易存在硬伤，例如：
• 关键词没写对或设置不当，导致数据没意义，比如：Gaussian里想要用PBE0泛函但写的关键词是PBE0（此时做的是PBE0-DH泛函计算）、用赝势基组时只定义了基组部分而没定义赝势、自定义基组时漏定义了某些原子或元素、Gaussian里用IOp自定义泛函时用opt freq关键词（IOp没法自动传递到freq这一步导致振动分析不是在与优化相同级别势能面的极小点计算的）、Gaussian里用后HF方法算偶极矩时没写density结果得到的是HF级别的、计算UV-Vis光谱时算的态数太少导致光谱不够完整、CP2K计算时用杂化泛函并且基于密度矩阵做积分屏蔽时没读取纯泛函波函数当初猜（《解决CP2K中SCF不收敛的方法》http://sobereva.com/665文中提到了）、CP2K算磁性体系时没写UKS结果当成了闭壳层算、Gaussian里按照特定方向加电场计算时不知道要写恰当用nosymm导致电场加的方向和期望的不符（nosymm用处见《谈谈Gaussian中的对称性与nosymm关键词的使用》http://sobereva.com/297）、用sobtop（http://sobereva.com/soft/Sobtop）直接产生拓扑文件后没有给里面添加原子电荷、用sobtop产生拓扑文件时误把可旋转二面角以谐振势对待、用Multiwfn做空穴-电子分析时没按照《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（http://sobereva.com/434）说的在Gaussian输入文件里写IOp(9/40=4)，等等
• 缺乏基本知识瞎算，比如：计算小晶胞体系不知道考虑k点、算极化率时不知道要带弥散函数、B3LYP算色散主导的弱相互作用却不加色散校正、CP2K里用MOLOPT基组计算需要CUTOFF很高的Na的时候用的CUTOFF不够、界面体系用了各向同性控压而不知道应该用半各向同性控压、算气-液界面体系在垂直于界面方向开了控压结果导致真空区消失、用Multiwfn基于CP2K产生的molden文件分析前没在里面添加盒子信息（参见《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》http://sobereva.com/651）、做Mulliken布居分析用的基组带了不该带的弥散函数、算溶剂环境下溶质自由能的时候忘了考虑标准态浓度转换问题（参看《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》http://sobereva.com/327）、试图用后HF方法计算分子轨道（却浑然不知程序虽然给出了结果但实际上是HF部分产生的）、用背景电荷表现环境分子的势场用的是对静电势重现性极烂的QEq电荷（相关讨论看《基于背景电荷计算分子在晶体环境中的吸收光谱http://sobereva.com/579），等等
• 计算流程/方式不对，比如：算高精度能量和性质之前结构没经过几何优化、直接把IRC两端的结构当做反应物和产物的准确结构用于计算能量和属性数据、振动分析和IRC计算之前结构没有在严格相同级别下优化过、对动力学轨迹计算某个分子的RMSD衡量内部结构变化的时候不知道先要消除整体平动和转动、算电子发射光谱之前没优化激发态、模拟拉曼光谱时没按照《使用Multiwfn计算特定方向的UV-Vis吸收光谱》（http://sobereva.com/648）里说的将拉曼活性转化成拉曼强度
• 数据读取错误。比如Gaussian里做了双杂化泛函计算，结果读成了其中SCF Done后面的中间量（怎么正确读看《谈谈该从Gaussian输出文件中的什么地方读电子能量》http://sobereva.com/488）。再比如做几何优化，Gaussian会对初始和最终结构都做布居分析，本来想读最终结构下的，结果误读成第一次输出的。再比如做开壳层体系的NBO分析时，NBO会对总密度矩阵、alpha密度矩阵和beta密度矩阵的分析结果依次输出，如果读错地方就糟糕了
• 对计算数据质量缺乏把控。如果数据质量太糟糕，自然可能无法重现文献数据。各个部分的计算都要保证数据质量。比如算稳定状态必须确保波函数稳定、确保无虚频。虽然做波函数稳定性测试和振动分析会增加耗时，但必要时必须做，否则可能得到错误的结果而无法重复文献。不要轻易使用比如Gaussian里opt=loose这种关键词放宽收敛限，否则很可能收敛限太松导致某些情况对文献数据重复性太差，尤其是此时做振动分析得到的频率的准确度很可能很糟糕。用取消SCF收敛情况检查的IOp(5/13=1)这种危险关键词的时候必须谨慎检查最终SCF收敛情况，千万别最后收敛精度很烂却直接使用其数据。再比如对于通过分子动力学计算物质平衡状态属性时，必须保证采样充分使统计误差足够小。再比如考察某种磷脂膜体系的特征时，需要模拟时间足够长以使得体系充分达到平衡，本身这类体系达到平衡所需的时间就偏长，动辄需要100 ns以上
• 没恰当考虑对称性问题。对于有对称性的体系，在Gaussian等支持对称性的量子化学程序里，要尽可能利用对称性，这样计算又好又快，内行的研究者的文章里对于能利用对称性的情况几乎都会利用。不利用对称性甚至可能还会得到和文献不符的结论，比如文献报道某个体系是平面的，但你一开始摆的结构是弯曲的，由于优化的收敛不可能无穷精确，最后你得到的结构可能轻微偏离平面，与文献结论不一致。如果体系本身没有那么高对称性，你一开始当成高对称性结构来计算，最后可能错误地得到了高对称性结构，也和文献不符（这种情况只要做一下振动分析看有没有虚频便知）。
• 其它低级错误：单位转换因子没用对或者用的不准确（比如百毒搜出来的大量转换因子都是错的）、该做构型/构象搜索的情况没做或者做得不当、该考虑构型/构象权重平均的情况没考虑、优化表面吸附问题时把被吸附分子放得太远导致还没优化到吸附状态就满足了收敛判据、处理数据用的Excel表格或者脚本存在问题、科学计数法记录的数据没读对（如漏掉了指数部分），等等
8 文献自身有问题
千万不要只怀疑自己而不怀疑文献，尽信文献不如无文献，很多文献里（包括IF很高的期刊的文章里）的数据是有错误的，甚至是非常明显的错误，上一节提到的各种犯错的例子也都出现在很多文献里。比如《我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment》（http://sobereva.com/584）里我就指出一篇文章里从头到尾的一大堆错误，诸如作者计算C18的时候大概率由于用的是非平面的初始结构，导致优化出的C18不是精确平面的，然后他们发现C18的极小点结构竟然有ECD信号，明显这是虚假的结果。
还有些文献里对数据的解释、标注是错的，是他们对数据理解有问题。比如直接把HOMO、LUMO能量的负值说成是氧化、还原势是很多涉及电化学的文章里常见的事，显然你如果以正确的方式来算，也即计算溶液下的电极反应的自由能变，结果肯定和文献对不上。再比如有不少文献用的费米能级是错的，对有gap的体系直接就把价带顶当做费米能级，这明显不符合正确的确定有限温度下费米能级的方式（正确定义参看wiki相应条目和Multiwfn手册3.300.9节）。
还有些文献里的计算方式从描述上就知道是错的。比如有的文献写G=Eele+ZPE-TS，有本科化学知识的人都知道这明摆着不对，H(T)和U(0)之间的差异上哪去了？就算是作为近似不考虑这部分、假设这部分在求差的时候能极大抵消，至少也应该把等号写成约等号。这种明显不合理的数据就别去重复了，非要重复的话那就暂时性降智故意用错误的公式。
很多文献作者也对数据很不负责任。比如之前我在IF很高的文献里看到补充材料中的Gaussian输入文件里居然全带着IOp(5/13=1)，这文献的数据的可靠性没法令人不放心（虽然作者可能检查了最终收敛情况，但谁也不敢说作者确实这么做了）。现在有很多无良代算机构，没一丁点搞理论计算的人的基本操守，为了能给出令做实验的人满意的与实验相符的计算数据来交差，极有可能在算不出期望的数据的情况下捏造假数据，这样的数据更别指望能被别人重复出来了。甚至有的代算方连计算都不计算，直接凭空胡编乱造数据，建议读者看看《谈谈我对计算化学代算的看法》（http://sobereva.com/505）里提及的情况。
有很多文献对计算的关键要点、直接影响数据可重复性的因素交代得非常粗糙、简略、语焉不详，甚至根本不提，这无疑给读者带来了很多误解、给别人重复计算造成了障碍。重复这些文章的数据只能是考虑各种可能情况去试了。另外，还有可能文章里用的参数、设置本来就不慎写错了，作者写的程序本身存在bug等等。
对于已经恰当考虑了前面说的所有因素，竭尽全力去重复文献，但还是重复不出来（且确实是有不可忽视的差异而不是在复现精度上过度较真），如果没有绝对必要去重复文献，那就别去重复了，极有可能文献的数据就是有问题的，或者没交代关键性计算细节。这种情况只要你能保证自己的计算是合理的就够了。但如果有特殊原因就是非要重复文献不可，显然该做的不是在思想家公社QQ群或者计算化学公社论坛（http://bbs.keinsci.com）等公开场所提问，因为找谁都不如直接找作者，世界上还有谁比作者更清楚数据是怎么来的的人么？世界上还有谁比作者更有回复你的义务？发邮件给通讯作者，向作者索要相应数据的输入输出文件，或者向作者提供你复现文献计算的输入输出文件问作者为什么没重复出来即可。当然很有可能一直没收到作者的回复，可能性有这些：
(1)作者看漏了你的邮件，或者邮件被误归入垃圾邮件，或者正忙着什么事没来及回复你结果后来又忘了。可以过一个礼拜（且可以换个邮箱）再发一次试试。如果通讯作者有多个，也可以给另外的人发
(2)措辞不当，不够礼貌和诚恳，而作者脾气又大或者又比较忙，又觉得回复你也没什么好处，就不回了。可以过一阵重新措辞再发一次。为了让作者有动力花时间给你回复，最好表示你会在未来的文章中会引用他的文献。
(3)文献中的数据是作者找别人代算的，代算方又是非常不负责任的，相关文件也不提供、计算细节也不交代清楚，甚至代算之后就再也联系不上了，因而作者也无能为力。
(4)做计算的是学生，学生毕业后找不到人了，导师又不清楚计算细节、没原始文件。
(5)文献中的计算本身就有严重错误或是假数据，可能收到你的邮件后作者才意识到，又或者可能之前就已经明知道用的是假数据，总之都不希望你知道他的数据存在问题，免得面子上挂不住，甚至怕文章最后被撤稿。
| 
发表于 Post on 2023-8-1 01:13:33
收藏 Add to favorites
发表于 Post on 2023-8-1 10:23:09
楼主 Author