从一篇CEJ文章来看报道计算化学结果时的常见问题
随着计算化学的普及，许多不同领域的科学研究工作者开始尝试使用计算化学的工具。很多情况下，简单的计算即可带来许多有意义的启示。这一方面使得更多人从中受益，另一方面，这些看似门槛较低的计算，使得许多并非专门受过计算化学研究训练的人参与到相关内容的研究、书写甚至审稿中，造成文献中对计算化学结果的不规范或错误报道屡见不鲜。本文以一篇CEJ文章为例，探讨一些报道计算化学结果时的常见问题。
本文标题为Polyarylimide and porphyrin based polymer microspheres for zinc ion hybrid capacitors （DOI: 10.1016/j.cej.2020.127038），报道了一类苝酰亚胺与卟啉共价结合的高分子在锌离子电容器中的应用。此类器件的原理与锌离子电池类似，在发挥作用时主要依赖锌离子与高分子结合/解离并转移电子，从而实现充放电。因此，为了探究电极材料在充放电时的分子行为，作者报道了如下计算结果：
作者首先展示了苝酰亚胺与卟啉连接形成的分子片段的分子表面静电势，从而指出分子中存在多处负电中心，容易结合锌离子。接下来计算了包含一个“卟啉”与2个苝酰亚胺片段的分子与1个和2个“锌离子”结合过程的能量变化，发现“锌离子”结合在“卟啉”环中心表现出极大的结合能（-6.49 eV），导致结合过程不可逆，而在酰胺氧原子上的结合则弱得多（-0.18 eV），从而可以实现可逆的结合和解离。此外作者对比了结合“锌离子”前后的高分子片段的前线轨道，发现在结合“锌离子”后HOMO/LUMO能级差降低，从而可推断呈现出更好的导电性。
01结构的不规范和不清晰
或许你已经注意到了，在介绍上述结果时，凡是提到锌离子和卟啉时都打上了引号，这是为什么呢？
仔细查看作者展示出来的结果，就会发现端倪。在Figure 4a中，作者展示的高分子结构，的确是卟啉；而在Figure 4b中，卟啉上的两个氢已经消失了！在这种情况下，作者声称的“卟啉与锌离子的结合能”，如果采用的真的是锌离子的话，可能对应如下2种反应（取代基略去）：
反应式（1）中，采用卟啉双负离子与锌离子结合，而反应式（2）中，则是直接将卟啉的两个氢去掉，此时已经变成了完全不同的物质。众所周知，卟啉是平面的大环芳香性结构，而如果将两个氢直接去掉，则将变成非芳香性甚至反芳香性的物种，两者将呈现截然不同的性质。作者并没有在文中注明是如何得到结合能的，因此我们只有尝试重复文中计算才能知道。
使用Gaussian 16结合PBE-D3BJ/6-31G(d) (for C, H, O, N)/SDD (for Zn)水平进行构型优化和频率计算，结合PBE-D3BJ/6-311+G(d,p) (for C, H, O, N)/def2-TZVP (for Zn)水平进行单点计算，得到结合过程的电子能量变化。虽然文中使用了CASTEP进行计算，但对于分子体系，Gaussian是好得多的选择，我们与作者同样采用PBE泛函，在可靠的基组水平下进行重复，由此得到反应（2）的能量变化分别为-20.3 eV，可见与文中报道的结果差异巨大。反应（1）涉及+2和-2电荷的结合，能量变化将更大，因此显然文中计算的不是其中任何一个反应。
这种情况下，我们不免怀疑虽然作者声称计算的是结构与锌离子的结合能，实际计算采用的却是锌原子。于是可以尝试计算如下反应的能量变化：
在上述提到的可靠计算水平下，得到上述反应的能量变化为-6.40 eV，与文中展示的-6.49 eV非常接近。因此我们可以得出结论，作者实际计算的，并不是认定合成出来、并展示在Figure 4a中的卟啉分子与文中所声称的锌离子的结合能，而是一个非芳香性的双脱氢卟啉分子与锌原子的结合能！
以下展示了在重复计算过程中优化得到的双脱氢卟啉结构，可见其结构跟卟啉差异很大。
那么作者为什么会采取这种做法呢？
我们假定作者在进行这些计算时是深思熟虑的，代入作者的角度，可能可以想到如下的理由，即卟啉与锌离子结合，是一个配体交换的过程，需要有碱协助脱去质子，又或许伴随着电子转移同时生成氢气，又或许伴随着其他反应：
而作者不想涉及这些反应过程，于是将卟啉的两个氢去掉，希望由此来反映“卟啉N原子”跟“锌元素”的结合。然而，这样带来的结果就是，最终计算的反应，与实际发生的反应，涉及的反应物存在很大的差异，由此带来的结果无法反映真实发生的过程。
事实上，如果想要计算上面列出来的这些反应，都是轻而易举的事情，只要根据实际体系中存在的溶剂环境正确处理质子和锌离子的自由能，就可以得到锌离子与卟啉发生配体交换的自由能变；如果认为反应伴随着电子转移，去计算反应发生的电位，也是很容易的事情。
退一步讲，即使作者认为对结构的选取有自己的道理，也应该在文中注明，以便读者能够正确理解文中展示的结果并进行重复。如果在描述计算结果时声称是锌离子与高分子的结合能，读者将不得不通过上述曲折的尝试才能猜出原文究竟计算了什么内容，对于别人是很大的误导。
作为一个与本文不直接相关但确实存在的问题，关于某物质是当做原子还是当做离子计算的问题，在许多表面吸附的计算中都会存在。这是由于主流第一性原理计算程序对于带电表面的处理存在缺陷，会带来不合理的能量，因此很长时间内在计算周期性表面上带电粒子的吸附现象时，都会默认将吸附质当做电中性粒子处理。而本文讨论的分子体系显然不在此列。另外即使是表面吸附，在描述基于中性粒子的计算结果时，将吸附质称作“离子”也不是好的做法，对于没有了解过此类计算的读者将构成误导。
02术语使用的错误
在正文讨论化合物的轨道信息时，作者将HOMO和LUMO的名字写错了。鉴于文中正确的写法和错误的写法同时出现，这应当属于笔误。虽然笔误难以避免，但对于这种极具标志性的术语应当格外注意，特别是书写论文属于正式场合，更不应该出现这种错误。
03缺失和不恰当的引用
对与本工作直接相关的研究的引用，既是对前人的尊重，也是读者理解文中结论的必要信息。在计算方法部分，往往涉及到许多不同的软件、方法，对它们进行恰当的引用是必须的。除了诸如Hartree-Fock这种属于教科书中的基本常识的方法外，均需要引用方法的原文，否则读者没有办法知道这个方法究竟是从哪里来的。
在本文中，作者使用PBE泛函进行计算，但跟在PBE后面引用的33~35分别是An All-electron numerical method for solving the local density functional for polyatomic molecules (J. Phys. Chem., 92 (1990), pp. 508-517), Density Functional Theory for Battery Materials (Energy Environ. Mater., 2 (4) (2019), pp. 264-279), Projector augmented-wave method (Phys. Rev. B, 50 (1994), pp. 17953-17979)，没有一篇是PBE的原文（Generalized gradient approximation made simple, Phys. Rev. Lett., 77 (1996) 3865-68.），属于无关引用。类似地，段落最后放在VMD程序后面的是一篇2017年的研究文章，同样属于无关引用。
更有甚者，虽然使用了CASTEP、Gaussian 09、Multiwfn、VMD等程序，以及M06-2X泛函和6-311+G(d)基组，但这些相关的引用基本缺失了。不仅如此，泛函和基组的名字还写得不对，M06-2X泛函的名字是其作者Truhlar专门规定的，中间的横杠必不可少；Pople系列基组中的G始终是大写，这些都是学术界约定俗成的规范，需要加强注意。
以CASTEP为例，CASTEP程序官网明确要求，凡是使用了CASTEP的工作，必须按照如下方式进行引用：
其他所有的程序，全部在官网或手册醒目位置对于引用有非常具体的要求。在使用这些程序一定要进行引用。有的程序作者比较佛系，不会追究，而如果遇到严肃认真的作者，轻则被拉入黑名单，重则被起诉、撤稿。以Multiwfn为例，在多种场合特别强调了引用的重要性：
以上是这篇CEJ文中暴露出的若干问题。在众多文献中，认真观察会有很多文章出现文章中提及的相关问题，本文希望给广大科研工作者普及更规范及严谨的计算化学常见问题，能在实际应用中更正确、规范地呈现结果。
以下是这篇推文的计算方法部分所涉及的引文列表：
Gaussian 16, Revision C.01, M. J. Frisch et al., Gaussian, Inc., Wallingford CT, 2016. (Gaussian)
Phys. Rev. Lett., 77 (1996) 3865-68. (PBE)
J. Chem. Phys., 2010, 132, 154104. (DFT-D3)
Theor. Chim. Acta 28, 213-222 (1973) (6-31G(d))
J. Chem. Phys. 56, 2257-2261 (1972) (6-31G(d))
Theor. Chem. Acc. 99, 231 (1998) (SDD)
J. Chem. Phys. 72, 650-654 (1980) (6-311+G(d,p))
J. Comput. Chem. 4, 294-301 (1983) (6-311+G(d,p))
Phys. Chem. Chem. Phys. 7, 3297 (2005) (def2-TZVP)
CYLview, 1.0b; Legault, C. Y., Université de Sherbrooke, 2009 (http://www.cylview.org) (CYLView)
