计算化学公社

标题: 我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment [打印本页]

作者
Author: sobereva 时间: 2021-3-2 02:41
标题: 我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment

我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment

My comment on a cyclo[18]carbon research article in the journal J. Mol. Model. that contains a lot of errors

文/Sobereva@北京科音 2021-Mar-1

1 前言

18碳环（cyclo[18]carbon）具有十分特殊的几何和电子结构，是近期最热门的分子。目前已经有好几十篇关于18碳环的理论研究工作，从不同角度对这个体系及其类似物的特点和性质进行分析探讨，笔者也已经发表了很多相关研究工作，文章汇总以及我写过的大量相关博文和评述见此页面：http://sobereva.com/carbon_ring.html，其中也有所有文章的下载链接。

已发表的18碳环研究文章中有很多优秀的工作，但也有水准不敢恭维的。笔者之前读到一篇名为Theoretical investigation on bond and spectrum of cyclo[18] carbon (C18) with sp-hybridized的理论计算研究18碳环的文章，见J. Mol. Model., 26, 111 (2020) (https://doi.org/10.1007/s00894-020-4344-5)。虽然J. Mol. Model.的影响因子不高，但真是难以想象一个专业的计算化学期刊上居然能出现从头到尾净是严重错误的文章，而且其中很多还是显著的低级错误，文中还有很多明显误导性结论。我非常纳闷审稿人都干嘛去了，就算俩审稿人都是外行或者都懒得认真看，起码编辑基于专业量化知识也不应该让这样的文章发出来才对。我深感不斧正一下的话恐怕会令一些读者对18碳环的特征产生错误的认识，于是连同合作者对此文写了一篇comment一一指出文中的各种问题并予以纠正，近期已发表于J. Mol. Model., 27, 42 (2021)，此文可以在https://rdcu.be/cdOpQ免费阅读。这篇comment文章里的内容不仅对于研究18碳环及类似体系的研究者有益，因为文中提及了此类体系的很多特征、提及了研究时需要注意的问题，研究其它体系的初学者们也值得看一下，可以少犯类似的低级错误。下面我就对我这篇comment文章的内容进行简要介绍，并给出一些额外的相关评论，更多细节请自行阅读comment文章。

2 计算级别的选取问题

此文作者居然用HF/6-31G**计算18碳环（注：对18碳环6-31G**和6-31G*是等同的）！哪怕有一点最基本的量子化学常识，也绝对不会拿这个级别做实际问题的计算。一般情况下用什么级别做计算在《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）和《谈谈量子化学中基组的选择》（http://sobereva.com/336）我都说过了，选择用HF真是匪夷所思，作者居然也不给出理由。需要注意的是，由于18碳环这个体系的几何结构和电子结构都非常特殊，所以切不可凭一般经验和感觉选择计算级别，比如常用的B3LYP算18碳环的结构就是完全错的，会把18碳环的长-短键交替结构算成是所有键长相同的结构，这在《谈谈量子化学研究中什么时候用B3LYP泛函优化几何结构是适当的》（http://sobereva.com/557）的第5节我特意说过。计算18碳环用什么级别合适，笔者在Carbon, 165, 468 (2020)一文中已经专门做了讨论，证明流行的ωB97XD是个很好的选择。在这篇comment文章中笔者又做了更全面系统的展开讨论，不同级别计算的结果汇总在了表1中，请大家查看。

作者用似乎不清楚RHF和UHF有什么区别，大抵也不了解对称破缺单重态的计算，在文章里没有明确提及HF用的是什么形式。笔者重复作者的数据后发现作者在文章中用的HF实际上是RHF。RHF下优化出来的18碳环结构是长-短键交替的，乍看起来与实验相符，但这真的就说明HF能凑合用么？答案是否定的！在comment文章中，我指出作者得到的RHF波函数其实是不稳定的，这在Gaussian里做一个stable任务就会发现。若改成UHF试图计算对称破缺单重态，会得到能量比RHF明显更低的解，这种计算详见《谈谈片段组合波函数与自旋极化单重态》（http://sobereva.com/82）中的介绍。在UHF下优化完的18碳环的无虚频结构是所有键长都相同的，这和实验事实明显不符，并且此时是自旋密度正负交替出现的多自由基状态的电子结构，见comment文中图1。所以作者用HF算出来看似还合理的18碳环结构实际上是典型的good result due to wrong reason！这也提醒初学者们注意，不要以为算出来的结果看上去能和实验基本对上就觉得计算合理、万事大吉了，一定要搞清楚原理和逻辑，以错误的方式算出了看似合理的结果并没有任何意义，碰上内行审稿人一定会被指出问题；而即便文章侥幸发出去了，也很可能之后被严谨、专业的研究者打脸。这里也特意说一下，用某个理论方法算某个体系的基态，不仅需要保证没有虚频，还应同时保证波函数是稳定的，这样的计算在原理上才是合理的，才能把数据拿出来说事，对于电子结构复杂的团簇类体系尤其要注意这点，比如在《揭示各种新奇的碳环体系的振动特征》（http://sobereva.com/578）介绍的笔者使用ωB97XD/def2-TZVP研究各种尺寸碳环的时候都同时检验了有无虚频和波函数稳定性。

对一般体系，6-31G**基组下优化几何结构可以得到定性合理的结果，在《浅谈为什么优化和振动分析不需要用大基组》（http://sobereva.com/387）我也说过几何优化用不着太大基组。但对于诸如18碳环这样的特殊的体系，就不能盲目地对用6-31G**做优化太有信心了。更何况18碳环也不大，文章作者冒着风险用保底的6-31G**是明显不妥的。如comment文章表1所示，6-31G**根本不足以合理描述18碳环的结构，18碳环的正确极小点本应该是D9h对称性的，ωB97XD在6-311G*、def-TZVP、def2-TZVP基组下都能得到这样的结构，但结合6-31G**时对称性会下降到C9h，这是因为所有C-C-C键角不再都是160度，而成了153.6和166.4度交替的结构。comment这篇文章的数据体现出，要想合理描述18碳环，最便宜的级别是ωB97XD/6-311G*。

作者不仅优化结构用的是严重不当的RHF/6-31G**，后面做各种其它计算和分析也都是基于这个级别的波函数，明摆着结果完全不可信。

3 成键分析的问题

作者对18碳环成键方式描述为“单-三键交替”，这是很多研究18碳环的人对这个体系成键方式的明显错误的认识。在18碳环于凝聚相下实验发现后一个月左右，我就写了一篇博文《谈谈18碳环的几何结构和电子结构》（http://sobereva.com/515），初步讨论了其实际的成键方式，后来在Carbon, 165, 468 (2020)一文又从各个角度对成键做了全面透彻的分析。在文章中我指出18碳环的成键方式应当被描述为“长-短键交替”或者“强-弱键交替”。必须正确认识到18碳环的pi电子离域特征，说成单-三键交替等于完全忽视了此体系极为显著的全局pi共轭了。

作者对18碳环做了NBO分析，声称NBO分析支持其所谓的单-三键交替的说法，这实际上是对NBO分析严重的错误理解！作者没说他们是怎么通过NBO得到单-三键交替的结论的，我猜大概率是看NBO轨道搜索所产生的BD型NBO轨道的数目。实际上这种方式判断成键特征极具误导性，因为会这完全忽视体系的pi共轭特征。举个最简单的例子，对苯分子做NBO分析，NBO会告诉你其中三个C-C键各有一个BD型轨道（sigma特征的），三个C-C键各有两个BD型轨道（一个sigma特征一个pi特征），即把苯分子的电子结构描述成了单一Lewis式所示的单-双键交替的情况。如果你就这么不假思索地轻易信了，以为苯分子中真是单双键交替的，那就完全被误导了。对于18碳环这种有高度pi共轭的体系，若你看一下pi特征的BD型轨道的占据数，就会发现它离整数2.0有显著偏离，这体现出当前的BD型轨道对应的Lewis式对此体系的电子结构的描述明显不理想，都算不上定性正确，所以根本不能随便拿它说事。

4 芳香性分析的问题

作者说NMR results show that C18 molecule does not have aromaticity in NICS = 0 position。这句话含糊不明，但大意是说18碳环的NICS(0)计算体现18碳环没有芳香性，这结论实在很不可思议。笔者在Carbon, 165, 468 (2020)中专门对18碳环做过NICS计算，还按照《通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性》（http://sobereva.com/216）做了更直观、对芳香性展现更全面的ICSS分析，还用《使用AICD 2.0绘制磁感应电流图》（http://sobereva.com/294）和《考察分子磁感生电流的程序GIMIC 2.0的使用》（http://sobereva.com/491）介绍的方法考察了18碳环的磁感生环电流，还用《使用Multiwfn计算AV1245指数研究大环的芳香性》（http://sobereva.com/519）介绍的方法通过AV1245指数做了考察，这些结论都充分指出18碳环不仅有芳香性，而且芳香性巨强（有双芳香性），其它研究者关于18碳环芳香性的研究也指出了这一点。为了搞清楚作者怎么在RHF/6-31G*下算出来NICS的结论是此体系无芳香性的，笔者在这个级别特意重复了他们的数据，发现NICS(0)的结果是-3.3 ppm，这不可忽略的负值已经体现出有芳香性了。由于RHF/6-31G*级别太垃圾，所以这个数据根本没什么实际意义，若改用合理得多的ωB97XD/def2-TZVP来算对18碳环更有意义的NICS(0)ZZ，结果将为-25.7 ppm，比苯分子的值-16.1 ppm还要负，这充分展现了18碳环有极强的芳香性。如果不熟悉NICS的话，参看《衡量芳香性的方法以及在Multiwfn中的计算》（http://sobereva.com/176）中的相应部分。

5 电子离域性的分析

作者说HOMO and LUMO also illustrate that there is no delocalization in C18 molecule，即作者通过HOMO和LUMO图形认为18碳环没有电子离域，这是何等**的结论！18碳环没有电子离域的话哪来的那么强的芳香性？哪来的那么显著的感生环电流？根据《在Multiwfn中单独考察pi电子结构特征》（http://sobereva.com/432）中介绍的方法，以及《谈谈18碳环的几何结构和电子结构》（http://sobereva.com/515）所示的具体绘制步骤，笔者在Carbon, 165, 468 (2020)给出了18碳环的LOL-pi图像，图中明确体现出18碳环有非常显著的18中心全局电子离域，而且是平面内和平面外分别有两套高度离域的pi电子。作者光根据HOMO和LUMO图形来讨论电子离域是严重错误的，真不知道作者怎么会认为靠这俩轨道就能讨论电子离域问题。比如LOL-pi对pi电子离域特征的描述，那是所有占据的pi轨道共同所决定的，既不是光由HOMO决定的，而且还和非占据的LUMO轨道毫无关系。顺带一提，作者在给出HOMO和LUMO轨道图形的时候，都只给出了一条轨道，然而18碳环由于其高对称性，前线轨道是简并的，这在笔者研究18碳环的电子光谱和非线性光学特征的Carbon, 165, 461 (2020)一文中专门说了，因此要说HOMO或LUMO的话也应该把简并的轨道一同拿出来说事。在此也提醒初学者，在分析有简并性的轨道的时候，简并的轨道必须一起说，否则会有极大的误导性。有趣的是，如作者原文里图1所示，18碳环的HOMO和LUMO明明是离域在整体的，就算按照作者所认为的HOMO和LUMO能判断电子离域的（错误的）逻辑，那也应该得到18碳环存在离域的pi电子才对，真不知道作者是怎么想的。

6 ECD（电子圆二色谱）的计算

作者对C6, C12, C16, C18, C20这些不同尺寸的碳环，以及18碳环的等电子体B9N9在RHF/6-31G**优化的结构下用TDHF做了ECD计算，并发现这些体系都有ECD信号。作者这里居然用TDHF，这和使用HF一样完全不可思议，而且也没给出任何理由。自从有了TDDFT，精度和CIS一样烂的TDHF就彻底被埋汰了，如今普遍用什么这里说了：《乱谈激发态的计算方法》（http://sobereva.com/265）。我在Carbon, 165, 461 (2020)里都是用ωB97XD做TDDFT研究电子光谱问题。

作者算出来的这些体系的ECD信号完全是虚假的。我不是说作者刻意造假，而是作者对有对称性体系的量子化学计算的基本常识实在太匮乏。由于这些体系极小点的高对称性，它们根本不可能有ECD信号，笔者在作者所用的计算级别下重复了他们的计算也完全证实了这一点。作者是怎么算出来ECD信号的？99.9%的可能性是作者在建模的时候就没有考虑体系的对称性，在GaussView里徒手画一个结构，也不知道做对称化就直接做几何优化了，几何优化完了也不去检验对称性问题。由于Gaussian的几何优化收敛的精度是有限的，作者这样优化完了的结构必定会轻微偏离实际的点群。我估计他们优化后用于ECD计算的结构肯定不是纯平面的，出现（虚假的）ECD信号也就不意外了。在这里提醒初学者，碰到这种原子团簇体系，一开始的结构如果是徒手画的没对称性的无所谓，如果优化后看到最终结构基本上有对称性，那么应当在GaussView里做一下对称化得到有对称性的结构，然后再次优化，这样才可能得到严格处于实际对称性的结构。确保结构处于正确的点群对称性对于原子团簇类体系极为重要。对称性的合理考虑和利用也是笔者每次在讲授北京科音的初级量子化学培训班（http://www.keinsci.com/workshop/KEQC_content.html）和中级量子化学培训班（http://www.keinsci.com/workshop/KBQC_content.html）的课程中都特意强调的。

7 导电性分析

作者在文中说the bandgap of C18 molecule is 9.70 eV，这暴露了作者对gap概念欠缺基本的正确认识，一个孤立体系哪来的band gap？好好看过此文就不会犯这种低级错误了：《正确地认识分子的能隙(gap)、HOMO和LUMO》（http://sobereva.com/543）。我重复了作者的计算，发现他们所谓的bandgap其实是指的HOMO-LUMO gap。作者拿HF方法算的gap说事完全没有意义，因为会严重高估。笔者在comment的文章里给了个具体例子，聚乙炔的实验带隙是1.5 eV，而笔者用HF/6-311G*算出来的带隙却高达5.63 eV，显然用HF的gap讨论没意义，只会产生误导。

8 分子的稳定性分析

作者居然说the binding energy of C18 molecule with HF/6-31G(d,p) is − 680.998 a.u. and − 685.23 a.u. in DFT/B3LYP/Lanl2dz. Obviously, the C18 molecule has good stability in DFT/B3LYP/Lanl2dz method, which is consecutive double bonds。我当时看到这里真是感到无语，很难想象作者搞量化计算到底搞了多久，怎么能写出这样惊悚的文字，而且难道审稿人和编辑都瞎了么？就这么短的一段话，就暴露出很多严重问题：
(1)拿两个不同的计算级别对比算出来的绝对能量根本毫无实际意义。日常用的计算级别算的电子能量都有极大的系统性误差，只有相同计算级别算出来的电子能量在求相对能量时才能把系统性误差绝大部分抵消掉从而得到有实际意义的结果。
(2)明明作者算出来的是电子能量，居然说成是结合能，明显不知道结合能在量化领域是指什么。建议作者看看此文《谈谈该从Gaussian输出文件中的什么地方读电子能量》（http://sobereva.com/488）。
(3)不知道怎么突然就蹦出来一个lanl2dz。lanl2DZ赝势基组对碳根本就没定义，在Gaussian里写lanl2DZ会给碳用完全过时的D95V基组。作者应该好好看看此文：《谈谈赝势基组的选用》（http://sobereva.com/373）。
(4)给出的能量保留的小数位数居然还不相同。而且以a.u.（此处等于Hartree）为单位时，一般都应当保留五位小数。保留两位小数的话可能使得数据记录误差达到13 kJ/mol。
(5)稳定性是由体系自身特征决定的，作者居然根据电子能量大小说18碳环在B3LYP/Lanl2dz下比HF/6-31G**下更稳定，这是何等清奇的逻辑？按照变分原理，基组越大能量就越低，这样的话用越大的基组就能算出来分子稳定性越高了？
(6)稳定性本来就是一个没有明确定义的词，热稳定性、光稳定性、化学稳定性等都是需要通过不同方式来体现的，而作者也没说这里的稳定性具体是指什么。

作者在所谓的B3LYP/Lanl2DZ级别下算出来体系中所有键的键长相等，就认为所有键是双键了（which is consecutive double bonds），这种说法明显不当，作者明显没认识到18碳环具有环平面内和环平面外两套离域pi电子的电子结构特征。作者后面还说了一句：maybe C18 molecule with consecutive double bonds exists and can be obtained one day, which has good stability than that of the C18 molecule with alternating single and triple bonds。这明显毫无道理，就因为B3LYP/lanl2DZ算出来的能量比HF/6-31G*更低，作者居然就认为B3LYP/lanl2DZ给出的结构就是更稳定的，并预测有朝一日能获得。倘若所有C-C键键长均等的结构是真实存在的，2019年Science那篇在NaCl表面STM观测18碳环结构的文章中肯定就直接观测到了。而且笔者在Carbon, 165, 468 (2020)中用可靠的CCSD/def-TZVP计算证实了长-短键交替的结构就是18碳环的极小点的稳定结构，在笔者Chem. Asian J., 16, 56 (2021)中的从头算动力学模拟过程中也进一步体现了这点，体系始终都是处于长短键交替的。而且在其他一些人的文章中还体现出在合理的计算级别下，所有C-C键长度相同的D18h结构实际上是18碳的长-短键反转的异构化过程的过渡态结构，明显不是一个极小点。

9 关于cyclo[24]carbon的结构

作者说cyclo[24]carbon，即环状C24的结构在他们用的HF/6-31G**级别下没有得到，居然什么原因都没说。这明显说不通，于是笔者在RHF和UHF级别下都做了优化，最终都得到了无虚频的结构，在comment文章里给出了，还同时给出了更可靠的在ωB97XD/def2-TZVP级别下优化的结构。

10 总结

不夸张地说，J. Mol. Model., 26, 111 (2020)这篇文章真是从头到尾几乎全是问题，错误的地方比正确的地方还多，反映出作者在量子化学研究方面存在大量知识漏洞。这也告诫初学者，做理论计算一定不能在没有基本常识的情况下瞎算和瞎讨论，否则就算真是特别走运文章发出去了，之后还是很容易被比较严谨认真的研究者写文章批评，在面子上不好看。这篇文章里暴露出的很多问题、漏洞也是许多缺乏经验和理论知识的初学者容易犯的问题、容易忽视的地方，我在本文都一一指出了，希望能令初学者们注意。

笔者的这篇comment文章不仅指出原文本身存在的大量问题，也同时对18碳环的很多特点、计算时应当注意的问题、计算级别的选取做了讨论，欢迎研究18碳环及其类似物的人阅读和引用。

作者
Author: wangxubo 时间: 2021-3-2 05:42
18碳环以及类似的碳环体系，很早就有一些人用很高级的理论方法算过了，现在还拿hf算算居然能发出来，真是恰烂文章恰疯了

作者
Author: wangzhe 时间: 2021-3-2 08:01
好奇这种文章会被撤稿吗

作者
Author: Accelerator 时间: 2021-3-2 10:35

作者
Author: sobereva 时间: 2021-3-2 14:24

wangzhe 发表于 2021-3-2 08:01
好奇这种文章会被撤稿吗

撤稿倒不至于，但我觉得会成为作者的黑历史

作者
Author: 喵星大佬 时间: 2021-3-2 14:38
简直是太惊悚了

作者
Author: niobium 时间: 2021-3-2 15:40
盲猜一作和通讯是夫妻俩。一堆挂名。
好好地材料不去做，要弄篇这种垃圾。计算的门槛还真低呢！

作者
Author: yygong 时间: 2021-3-2 16:42

niobium 发表于 2021-3-2 15:40
盲猜一作和通讯是夫妻俩。一堆挂名。
好好地材料不去做，要弄篇这种垃圾。计算的门槛还真低呢！

瞎算是低；掌握原理，科学解决问题也是挺厉害的，

作者
Author: zjxitcc 时间: 2021-3-2 16:51

niobium 发表于 2021-3-2 15:40
盲猜一作和通讯是夫妻俩。一堆挂名。
好好地材料不去做，要弄篇这种垃圾。计算的门槛还真低呢！

要不要加个多少积分才可以浏览的权限，感觉有点危

作者
Author: niobium 时间: 2021-3-2 17:05

zjxitcc 发表于 2021-3-2 16:51
要不要加个多少积分才可以浏览的权限，感觉有点危

哈哈，没必要

作者
Author: Freeman 时间: 2021-3-2 21:49
“所有C-C键长度相同的D9d结构实际上是18碳的长-短键反转的异构化过程的过渡态结构”
这里是不是应该改为D18h？D9d连平面都不是了。

作者
Author: sobereva 时间: 2021-3-3 01:40

Freeman 发表于 2021-3-2 21:49
“所有C-C键长度相同的D9d结构实际上是18碳的长-短键反转的异构化过程的过渡态结构”
这里是不是应该改为D ...

对。之前gview把结构直接判断成D9d就写上去了，疏忽了
文字已改

作者
Author: Daniel_Arndt 时间: 2021-3-3 08:17
社长关于RHF与UHF那部分的讨论让我想起来Gustavo E. Scuseria的一个讲座，里面提到了C60用RHF优化出的几何构型符合Ih点群，但波函数有UHF不稳定性。如果用UHF去优化C60，得到的几何构型就不符合Ih点群，与实验上C60的NMR只有一个信号相矛盾。他还提到了UHF优化C60后，得到的波函数是GHF不稳定的。用GHF去优化C60，可以得到符合Ih点群的几何构型。他那天的讲座是用spin-projected GHF去研究C60的。

也不知道有没有人用GHF去优化过18碳环的几何构型。（本人不熟悉相对论效应，请勿询问本人

）

作者
Author: 函数与激情 时间: 2021-3-3 19:11
有点好奇，comment需要找审稿人吗

作者
Author: sobereva 时间: 2021-3-4 00:41

函数与激情发表于 2021-3-3 19:11
有点好奇，comment需要找审稿人吗

需要。是两个人，一个给原作者看，另一个是单独找的审稿人评判comment文章内容

作者
Author: sobereva 时间: 2021-3-4 01:00

Daniel_Arndt 发表于 2021-3-3 08:17
社长关于RHF与UHF那部分的讨论让我想起来Gustavo E. Scuseria的一个讲座，里面提到了C60用RHF优化出的几何 ...

在我看过的几十篇文章中目前还没有过。就算GHF能用，但毕竟没考虑动态相关，还是不如直接用合适的泛函做DFT计算。
验证了一下，确实RHF对C60有波函数不稳定性。

作者
Author: shalene 时间: 2021-3-5 17:14
没底线。猜是学生自己做的工作。

作者
Author: niobium 时间: 2021-3-5 17:23

shalene 发表于 2021-3-5 17:14
没底线。猜是学生自己做的工作。

参看我7楼的回复

作者
Author: mfdsrax2 时间: 2022-5-8 15:46
我有一个问题，不管用什么理论方法，都能算出来一个结果，哪个可靠哪个不可靠是怎么判断的，和实验数据比对吗？

作者
Author: sobereva 时间: 2022-5-9 21:38

mfdsrax2 发表于 2022-5-8 15:46
我有一个问题，不管用什么理论方法，都能算出来一个结果，哪个可靠哪个不可靠是怎么判断的，和实验数据比对 ...

此文对你理解这个问题有益
谈谈量子化学研究中什么样的benchmark才有意义
http://sobereva.com/554（http://bbs.keinsci.com/thread-17826-1-1.html）

作者
Author: ionexchangeC 时间: 2023-5-26 01:33
今天发现HF方法对大共轭体系可能是普遍存在文中的问题。我用HF优化18轮烯，也出现了RHF波函数不稳定的情况，换UHF之后会优化成完全的平面结构。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)