Results of the Computational Chemistry Commune 2024 poll on the popularity of computational chemistry programs and DFT functionals

文/Sobereva@北京科音  2024-May-5

0 前言

2024年4月4号，在北京科音建立的人气最高、专业性最强的综合性计算化学论坛“计算化学公社”（http://bbs.keinsci.com）上开展了为期一个月的五项投票：
你最常用的做 计算的DFT泛函投票（http://bbs.keinsci.com/thread-44387-1-1.html）
你最常用的做第一性原理计算的DFT泛函投票（http://bbs.keinsci.com/thread-44391-1-1.html）
你最常用的 程序投票（http://bbs.keinsci.com/thread-44388-1-1.html）
你最常用的分子动力学程序投票（http://bbs.keinsci.com/thread-44389-1-1.html）
你最常用的第一性原理程序投票（http://bbs.keinsci.com/thread-44390-1-1.html）

现对投票结果进行总结和简单评论。未来预计每三年重新开展一次投票。要强调的是，这个投票只是体现流行程度，和方法/程序的好坏并没直接关系。本投票结果主要反映中国的计算化学领域的专业、内行群体的情况，不反映业余/外行群体的情况。本次投票的结果也有助于业余/外行研究者正确认清什么才是主流，减少被他人忽悠、听信歪曲说辞误入歧途的概率。

上一次的投票于2021年举行，当时的结果和很多相关评论见下文：
2021年计算化学公社论坛“你最常用的计算化学程序和DFT泛函”投票结果统计
//www.umsyar.com/599（http://bbs.keinsci.com/thread-23482-1-1.html）

1 你最常用的做 计算的DFT泛函投票

本次可投的泛函有43种，带不带色散校正算同一种泛函。在2021年的投票条目基础上有所增减，特别是增加了双杂化泛函，明显几乎不会有人用的泛函没纳入可投范围。投票范畴仅限做 计算的情况，不包含第一性原理计算的情况。关系特别近的，比如M05-2X和M06-2X、wB97X和wB97XD和wB97X-D3、SCAN和r2SCAN、revDSD-PBEP86-D3(BJ)和DSD-PBEP86-D3(BJ)等等当做同一个泛函来计。此次投票者共713人。本投票每个人最多选6项，且所投的泛函必须占平时全部研究工作的5%以上。按照得票率（票数除以总投票人数）绘制的图如下。为了避免条目过多，只把得票前30名的列出。此图中诸如某泛函对应50%就代表有50%的人平时较多使用此泛函，后文的统计图同理。

总的来说本年度的投票结果和2021年时没太大变化，前六名的顺次没有改变，还是依次为B3LYP、(M05/M06)-2X、PBE0、wB97X-(/D/D3)、PBE、CAM-B3LYP老几样，各自有各自的用处（参看我对2021年投票结果的评论//www.umsyar.com/599，这里不再赘述）。它们的得票率的相对比例也基本没变，也就是量化领域里用处相对有限的PBE的比例有一定降低，以后肯定还得跌。2021年时候的第7名M06虽然得票率如今还是5%左右，但排名已下滑到第10，被wB97M-V、SCAN/r2-SCAN、PWPB95-D3(BJ)所超过。M06在我来看用处着实不大，虽然计算过渡金属配合物体系有一定用处，但PBE0-D3(BJ)/D4多数情况是更好的选择，而且M06还有后继者MN15可用。wB97M-V的得票率从2018年的3.1%提升到了如今的10%，已经算是增幅很快了，再过3年统计时肯定还会增高。此泛函在国内量化研究者中一定程度的流行，很大程度在于在计算化学公社论坛和思想家公社QQ群的讨论中时常被提及、在《简谈 计算中DFT泛函的选择》（//www.umsyar.com/272）博文中和我在北京科音基础（中级） 培训班（http://www.keinsci.com/workshop/KBQC_content.html）中的推荐、被免费的ORCA程序支持。提出时间不算很长的纯泛函SCAN及其改进版r2SCAN现在得票率能到6%着实不容易，2021年时得票率还不到1%，这主要在于有越来越多的程序已经支持此泛函，而且综合素质整体强于更早的经典泛函PBE，因而在纯泛函范畴内能有重要的位置。

2021年投票的时候没纳入双杂化泛函，这次得票率超过1%的双杂化泛函的排名顺序是PWPB95-D3(BJ)(5.9%) > (rev)DSD-PBEP86-D3(BJ)(3.1%) > B2PLYP-D3(BJ) (2.7%) > ωB97X-2-D3(BJ) (2.0%)。PWPB95-D3(BJ)和(rev)DSD-PBEP86-D3(BJ)能在国内用户中变得流行和上述wB97M-V的情况很类似。本身这俩泛函各有长处，有流行开来的必然性。PWPB95-D3(BJ)比较robust，算过渡金属配合物能量问题较好，能在ORCA里用；而revDSD-PBEP86-D3(BJ)算主族体系反应能、势垒以及弱相互作用能都是双杂化里顶尖的，而且在Gaussian里通过《Gaussian中非内置的理论方法和泛函的用法》（//www.umsyar.com/344）我介绍的方法能直接用。此外，ORCA中DSD-PBEP86适合算TDDFT和NMR目的也是其加分项。这俩泛函现在流行度能超过经典且最早提出的双杂化泛函B2PLYP是其应得的。

BLYP这次的排名降幅很大，从第10名已跌到第22名，本身这个泛函如今鲜有用武之地。BLYP一般也就算算主族体系，在Gaussian里用这个明显不如用B3LYP，耗时也持平，而以前在ORCA里用这个搭配def2-SVP结合RIJ加速做有机体系几何优化速度效率高是个用处，以前我在《大体系弱相互作用计算的解决之道》（//www.umsyar.com/214）里也提过，但如今也不如改用B97-3c来跑。

2 你最常用的做第一性原理计算的DFT泛函投票

可投的泛函有26种，带不带色散校正算同一种泛函。此投票在2021年没有，是本次新加的。此次投票者共442人。本投票每个人最多选6项，且所投的泛函必须占平时全部研究工作的5%以上。结果如下，零票的未显示

96年提出的PBE至今仍稳居第1的位置，如同B3LYP在 领域的地位，而且和第二名相差更悬殊。PBE能有这样的地位是必然的，PBE提出年代早、被程序支持得极为广泛，计算便宜，有丰富的专门为其搞的赝势，几何优化和分子动力学目的大多数时候够用（加色散校正后又拓宽了其普适性），而且在基态能量相关问题方面依然有使用价值而且没被已流行的其它纯泛函甩开特别多（这和B3LYP在 领域的情况截然不同，B3LYP算能量早过时了，很难再发得出去文章，见http://bbs.keinsci.com/thread-12773-1-1.html）。B3LYP在这次投票里得了第2，令我有点意外，大概率是很多人不好好看投票帖子的说明，误把 研究用的泛函也在此进行投票了。PBE0能排第3不意外，需要一个HF成分适当的杂化泛函做TDDFT/TDDFPT算激发态、算强相关体系的问题时经常用得着。HSE06流行度排得上第4主要来自于其计算带隙、能带方面公认很好，以及晶胞参数优化方面表现不错。PBEsol是优化晶体结构、晶胞参数的好把式，而且还是便宜的纯泛函，能排到第5很正常。M06-2X能排第6令我有点意外，一方面必定是有人误当成 计算的情况来投，另一方面是计算主族晶体/液体相关的化学反应、吸附的相关能量问题必定有一些人在用。SCAN/r2SCAN已经有一定流行度，由于在文献中出现频率越来越高，在未来的流行度必定也会逐渐提升，但流行度超越PBE在可预见的未来还不太可能，毕竟对于大量PBE就已经表现得够用的范畴，作为更贵但没带来显著优势的meta-GGA的SCAN/r2SCAN不可能显著侵占这方面的市场。第一性原理领域里用BLYP的人我不很理解是什么心态。revPBE和与之相似的RPBE能有一定流行度在于算化学吸附方面不错，被不少人用。第一性原理方面的泛函投票中CAM-B3LYP显得远不如 领域里来得流行，估计用这个的大部分都是CP2K用户用来算激发态和UV-Vis谱方面，只占投票的少量群体。算化学吸附很好的BEEF-vdW和算物理吸附很好的optB88-vdW能有一定票数不算意外。纯泛函中矮子里拔将军算带隙往往可以接受的HLE17在本次投票中获得了一点流行度，略意外的是算带隙整体更好些的纯泛函mBJ反倒在这次投票中显得无人问津，可能是前者在CP2K里能直接用而后者不能所致。作为PBE后继提出来的知名的TPSS流行度那么低有点令我意外，倒也确实实际用处不太大，现在又有了理论上以及实际整体表现得更好的r2SCAN。PW91虽然在文献里出现得很多，但这次得票率相当低，毕竟实际中有PBE就基本没有更老的PW91能派上用场的时候。B97M-rV能有少量票数，主要在于算热力学量方面在纯泛函里是表现得较突出的。

3 你最常用的 程序投票

可投程序有29种，投票者共679人。本投票每个人最多选三项，且所投的程序必须占平时全部研究工作的10%以上。按照得票率绘制的图如下，只显示了得票前20名的

前三位和2021年投票的结果一样，还是Gaussian > ORCA > xtb，Gaussian依然是约90%的 研究者日常必用的程序，稳稳占据主导位置。而ORCA和xtb的得票率比2021年时都有了约10%增长，这是意料之中的。实际上这三个程序也是我自己最常用的：xtb拿来快速预优化和结合molclus（http://www.keinsci.com/research/molclus.html）做构象搜索的初筛，Gaussian做基于DFT的opt、freq、扫描、IRC等涉及几何变化的任务以及算一些属性（NMR、超极化率等），ORCA算高精度单点。这三个程序的流行度远远甩开了其它程序，一方面是它们比较容易安装和使用，一方面各有各的独特优势，有被大量使用的刚性原因。而且它们把大部分 计算的应用领域都给覆盖了，对于日常应用性研究来说其它程序难以有和它们竞争的显著空间。Dmol3、ADF、Q-Chem、Turbomole这四个商业程序日子愈发不好过。在量化计算方面没有长处还巨贵的Dmol3的用户越来越少，从2021年的6.2%已经进一步萎缩到了4.3%，以后肯定还会明显进一步萎缩。ADF从2021年时候的仅仅2.2%进一步萎缩到了1.5%。Q-Chem从2021年的3.0%萎缩到了1.0%。Turbomole从2021年的1.6%萎缩到了1.0%。以GPU加速为卖点的TeraChem更不幸，2021年时候还有1人投票，今年变成了0人。

4 你最常用的第一性原理程序投票

可投程序有25种，投票者共565人。本投票每个人最多选三项，且所投的程序必须占平时全部研究工作的10%以上。按照得票率绘制的图如下（0票的没显示）

根据这次投票的结果可见，至少在计算化学公社论坛里，CP2K的流行程度已经赶超了VASP。这令我90%程度感到意外，但也有10%程度感觉是在情理之中。由于Multiwfn在2021年开始提供了极其易用的创建CP2K输入文件的功能（//www.umsyar.com/587），我后来又对此功能反复打磨并在Multiwfn中提供了对CP2K绘制DOS、能带、STM、电子激发、成键分析等许多功能，再加上2023年3月、2023年10月和2024年3月开办了三期北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）非常全面系统讲解了CP2K的使用，无疑令中国的CP2K用户猛增。但即便我已预料CP2K的得票率肯定会远高于2021年时候的33.3%，但也没预料到这次居然能达到60.9%，直接翻了将近一倍，甚至把一直霸占流行度榜首的VASP给挤下去了。由于免费且十分高效的CP2K的用户还在不断激增，而且CP2K更新很快，不断完善和发展新功能，Multiwfn在未来还会给其提供更多的相关分析处理功能，CP2K的位置在以后无疑还会更加牢固。相比之下，VASP的流行度从2021年投票时候的65.8%降到了现在的58.1%。由于现在有非常强大的竞争者CP2K，而且CP2K不具备的一些功能还有免费的Quantum ESPRESSO能用来平替VASP，售价较昂贵且算大体系速度远不如CP2K的VASP在未来的前景不乐观。以前很多人一提到第一性原理计算就言必称VASP，以后就不再是如此了。除了CP2K的流行度猛增外，其它程序的流行度都普遍出现了下降，如CASTEP和Dmol3分别从2021年的19.0%和9.3%下降到了13.8%和6.4%。Wien2k今年更是连一票都没有，而2021年时还有3票。那些程序流行度百分比减少，一方面是它们的票数大多数确实有一定减少，用户发现有更理想或免费的程序可用，另一方面原因应当是有很多通过CP2K程序开始入手第一性原理计算的人加入了投票，他们只对CP2K的得票率有贡献而间接地拉低了其它程序的得票率。

5 你最常用的分子动力学程序投票

可投程序有18种，投票者共551人。本投票每个人最多选三项，且所投的程序必须占平时全部研究工作的10%以上。按照得票率绘制的图如下

GROMACS依旧是用户数的龙头老大，而且流行度从2021年投票时的69.3%还进一步略微提升到了71.3%，得票数大约等于其它所有程序用户数之和，和曾经的情况一样。第2、3位依然分别是Lammps和AMBER。Lammps和OpenMM得票率略涨，而AMBER、Forcite和NAMD的流行度都有较多降幅，GULP、DL_POLY和CHARMM更是快跌没了。

Several tools for constructing three-dimensional structures based on nucleic acid sequences

文/Sobereva@北京科音   2023-Dec-20

之前我在《几种基于氨基酸序列构建很简单蛋白质三维结构的工具》（//www.umsyar.com/687）中介绍过一些基于氨基酸序列构建简单蛋白质三维结构的工具，本文将介绍几种基于核酸序列构建DNA/RNA三维结构的工具，可以用于做分子动力学模拟、分子对接等目的。虽然还有很多其它程序也可以构建，如HyperChem等，但本文提供的这些就已经足够用了，且都是免费的。这些工具在产生核酸结构时只需要指定一条链的序列，从5'端到3'端，对于产生双链结构的情况，另一条链的序列总是自动按照规范DNA中标准碱基配对方式自动确定的。这些程序都可以保存成常用的pdb文件格式，并且原子名是规范的。

1 在线工具DNA Sequence to Structure

地址：http://www.scfbio-iitd.res.in/software/drugdesign/bdna.jsp

输入DNA序列以及DNA结构类型，即可立刻返回产生的pdb结构，例如：

返回的结构用VMD查看：

2 在线工具web.x3dna.org

地址：http://web.x3dna.org

进入后，选Rebuilding - combination of A-, B-, or C-form DNA models。之后可以输入DNA序列由几段构成，比如设了3，点next，若三段内容分别按下面这样设，那么DNA序列就是AAACCCCGGG，且其中AAA部分是A-DNA形式、CCCC部分是B-DNA形式、GGG部分是C-DNA形式。

提交之后，过一会儿（有可能时间挺长），看到下图，可以点击链接下载pdb文件

3 AmberTools的NAB

AmberTools程序包可以在http://ambermd.org下载，NAB是AmberTools中的组件，AmberTools装好后NAB就可以直接用了。最简单的运行方式为nab test.nab -o test.out，这里test.nab是NAB程序的输入文件（后缀必须是nab）。NAB就像编译器一样会编译出名为test.out的可执行程序，然后运行./test.out即可使里面的指令生效。

NAB可以用于创建DNA和RNA序列。例如创建一个序列为gcgttaacgc的B-DNA结构，就创建一个文本文件比如叫genDNA.nab，里面写以下内容

molecule m;
m = fd_helix("abdna","gcgttaacgc","dna" );
putpdb( "sobDNA.pdb", m );

之后运行nab genDNA.nab -o genDNA，当前目录下就出现了名为genDNA的可执行文件。再输入./genDNA运行之，当前目录下就出现了sobDNA.pdb，是我们要的DNA的结构，DNA的骨架顺着Z轴。

从上面例子可见fd_helix函数里面跟了三个参数，第一个参数控制产生的核酸类型，第二个参数是序列，第三个参数写dna就是生成DNA、写rna就是生成RNA。

4 Gabedit

Gabedit是一个免费的可视化程序，可以在http://gabedit.sourceforge.net下载。启动后点击菜单栏Geometry - Draw，然后点右键选Build - polyNucleic Acid，之后一边点击碱基名字的按钮，三维结构一边不断产生，如下图所示。可见核酸类型和结构形式都可以自己定义。如果选上Add Counter Ion，产生的核酸结构的磷酸基旁边还会自动加上Na+作为抗衡离子。构建好后，在图形窗口上点右键选Save as，就可以选择保存成pdb格式。

5 Avogadro

Avogadro可视化程序可以在http://avogadro.cc免费下载。启动Avogadro后，点击菜单栏的Build - Insert - DNA/RNA，就蹦出了如下窗口。然后一边点击按钮输入核酸序列，一边图形窗口里就可以看到生成的核酸结构。DNA和RNA，单链和双链，结构形式都可以自由选择。

Several tools for constructing three-dimensional structures of very simple proteins based on amino acid sequences

文/Sobereva@北京科音   2023-Oct-15

经常有人问已知很简单蛋白质（小肽）的氨基酸序列，怎么得到其三维立体结构，这里就介绍几种工具。本文并不是对相关程序做完整的收集，介绍的这些已经足够满足各种实际情况的需要了。注意本文的范畴和蛋白质结构预测完全不同，那类程序是用于构造有三级结构的较复杂蛋白质用的，本文涉及的工具顶多只牵扯到二级结构。

需要注意的是，只有其中部分程序可以产生原子名符合IUPAC规范的pdb结构文件。仅对于这种情况，pdb文件才可以载入VMD后正常显示二级结构（至少得骨架重原子名满足IUPAC规范），以及能够直接被比如GROMACS的pdb2gmx所处理来产生拓扑文件。

笔者另外写了篇文章介绍通过输入序列产生DNA或RNA三维结构的工具，见《几种基于核酸序列构建三维结构的工具》（//www.umsyar.com/692）。

1 ProBuilder

这是个在线程序，网址是https://www.ddl.unimi.it/vegaol/probuilder.htm。用户只需要输入氨基酸序列，比如RLCVVPIRLPPLRERT，并且指定要产生的二级结构的phi、psi、omega角度，就可以获得结构文件。能产生的格式很多，包括符合IUPAC命名规范的pdb文件。此程序产生的结构可以要求带和不带侧链，带侧链的情况由于不会自动做结构优化，因此侧链往往有严重重叠，后面介绍的很多程序也都是如此。

2 GaussView

Gaussian的官方界面程序GaussView提供了生物分子片段库，如下所示。可以选择位于链中间、N端、C端的各种氨基酸残基片段，在建模窗口点击一次出现第一个残基后，再一次一次点击肽键末端，就可以逐渐得到满足要求的肽链。

建好后可以保存成结构文件，包括pdb格式。用此程序的优点是可以自行精细修改三维结构，缺点是没法把一长串序列一下子转化为三维结构，而且这么保存的pdb文件里的原子名只是元素名而不是满足IUPAC规范的名字，里面连残基名都没有，而且也没法像ProBuilder那样直接设置对应特定二级结构的phi和psi角度（得手动一个一个调）。

3 Avogadro

免费的可视化程序Avogadro（http://avogadro.cc）的Build - Insert - Peptide界面可以构建有特定二级结构、特定序列的小肽的三维结构，界面如下所示。可以保存出许多格式，保存成pdb的话原子名和残基名都是符合IUPAC规范的。此程序还自己支持GAFF、MMFF94等力场做几何优化，可以用于把建模后侧链的不合理接触消除掉。

4 gabedit

启动免费的可视化程序gabedit（http://gabedit.sourceforge.net）后，点菜单栏的Geometry - Draw，在新窗口里点右键，选Build - PolyPeptide，就可以看到构建多肽的界面，如下所示，也可以类似于Avogadro去产生符合特定二级结构的指定序列的小肽。每点击一次氨基酸的按钮在图形窗口里就会长出来相应的残基。之后点右键选save as，就可以保存成各种格式，包括满足IUPAC原子名和残基名规范的pdb格式。

5 AmberTools中的leap

知名的AmberTools（http://ambermd.org）程序里的leap可以用于根据氨基酸序列产生线型小肽的三维结构并保存成符合IUPAC命名规范的pdb文件。

参考手册或《Amber14安装方法》（//www.umsyar.com/263）安装AmberTools后，在命令行窗口输入诸如tleap -f leaprc.protein.ff14SB即可启动纯文本界面的leap（tleap）且同时载入AMBER14SB力场相关的氨基酸的库文件。在里面可以输入比如如下命令定义一个序列，其中残基的N结尾和C开头分别代表N端和C端残基
TC5b = sequence { NASN LEU TYR ILE CGLN }
之后再运行以下命令就可以把序列以三维结构保存成/sob目录下的TC5b_linear.pdb了
savepdb TC5b /sob/TC5b_linear.pdb

6 PeptideBuilder

PeptideBuilder是个基于Python的轻量级构建小肽的库，网址是https://github.com/clauswilke/PeptideBuilder，具体介绍见原文PeerJ, 1, e80 (2013)。

在Linux下联网状态运行pip install PeptideBuilder即可安装。之后可以基于PeptideBuilder写Python脚本来构造小肽结构，可以指定psi和phi角度、有哪些残基等。此程序只产生重原子坐标，没有氢。产生的pdb文件符合IUPAC命名规范。

以下是个例子，比如保存为createQ6.py，之后运行python createQ6.py，就可以在当前目录下得到含有6个谷氨酰胺的yohane.pdb文件。

from PeptideBuilder import Geometry
import PeptideBuilder

geo = Geometry.geometry("Q")
geo.phi = -60
geo.psi_im1 = -40
structure = PeptideBuilder.initialize_res(geo)
for i in range(5):
    PeptideBuilder.add_residue(structure, geo)
PeptideBuilder.add_terminal_OXT(structure)

import Bio.PDB

out = Bio.PDB.PDBIO()
out.set_structure(structure)
out.save("yohane.pdb")

有人基于PeptideBuilder做了扩展写了额外Python程序，能生成具有指定序列在N端加上Fmoc或ACE（乙酰）封闭、在C端加上NME（甲胺基）封闭的pdb文件。见《生成任意序列的封端短肽pdb的脚本CappedPeptideBuilder.py》（http://bbs.keinsci.com/thread-30278-1-1.html）。一个简单例子：先把PeptideBuilder装上，然后把CappedPeptideBuilder.py放在当前目录下，再运行python CappedPeptideBuilder.py -s KEYIN -p test，此时当前目录下就出现了test目录，其中ACE_KEYIN.pdb就是N端加了乙酰、C端加了甲氨基的序列为KEYIN的小肽结构了（没有氢）。

7 Tinker

动力学模拟程序Tinker有根据序列建立蛋白质三维结构的命令protein，参见计算化学公社论坛的帖子：http://bbs.keinsci.com/forum.php?mod=redirect&goto=findpost&ptid=11478&pid=79011&fromuid=1

Theoretical design of a novel dual-motor supramolecular system based on cyclo[18]carbon

文/Sobereva@北京科音   2023-Aug-14

1 前言

北京科音自然科学研究中心（www.keinsci.com）的卢天和江苏科技大学的刘泽玉等人近期设计了由18碳环（cyclo[18]carbon, C18）和具有两个大环的分子（OPP）构成的非常新颖的双马达超分子体系，并做了详细的分子动力学特征的研究，相关工作已发表在知名的Chemical Communications通讯刊物上，文章信息如下，欢迎阅读和引用：
Zeyu Liu,* Xia Wang, Tian Lu,* et al., Theoretical design of a dual-motor nanorotator composed of all-carboatomic cyclo[18]carbon and a figure-of-eight carbon hoop, Chem. Commun., 59, 9770 (2023) DOI: 10.1039/D3CC02262E
https://pubs.rsc.org/en/content/articlelanding/2023/CC/D3CC02262E

链接：https://pan.baidu.com/s/1H7AcSnDqbjBOE2ajFk-LXA?pwd=tnzv

下面将对此文的主要研究内容和思想进行深入浅出的介绍，以有助于读者更好地理解文章内容。下图是被研究的双马达超分子体系的结构图和运动方式示意

值得一提的是，同作者之前还深入细致地考察了上文研究的OPP对18碳环的吸附和释放的动力学过程及相互作用本质，已发表于Phys. Chem. Chem. Phys., 25, 16707 (2023) DOI: 10.1039/d3cp01896b，并在《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）中做了详细介绍，推荐阅读。同作者之前还对18碳环及衍生物的各方面性质还做了非常广泛的研究，相关工作汇总见//www.umsyar.com/carbon_ring.html。

2 设计思想

18碳环这个奇妙的分子虽然在半个多世纪以前就被理论预测，但直到2019年才首次在凝聚相实验观测到，并迅速引发了巨大关注。它的直径和C60富勒烯很相似，对比见下图

具有较大双环结构的OPP分子在Angew. Chem., Int. Ed., 61, e202113334 (2022)中被首次合成出来，实验上得到了它带有C60和C70富勒烯的晶体，C60@OPP的晶体结构如下所示，可见C60被嵌入在了大环里面

由以上结构特征可以想到两个C18也可以分别内嵌到OPP的两个大环中，成为2C18@OPP。倘若每个C18在OPP的大环里能够由于热的驱动而较容易地发生转动，自然就成了双马达超分子体系。这样的体系在之前的文章中是没有报道的，此体系或类似物在未来有望成为构造复杂分子机器的关键组成部分。本文介绍的Chem. Commun.这篇文章的目的就在于通过 和分子动力学模拟来证实OPP与C18复合作为双马达体系的可能性并考察其特点。

3 基于 的能量角度的研究

文中首先使用 方法从能量角度对C18@OPP和2C18@OPP进行了研究。首先使用Gaussian程序通过ωB97XD泛函优化了复合物结构并做了振动分析。2C18@OPP达到了260个原子，为了节约时间，对占体系大部分原子的OPP部分用了6-31G*基组，而关键的C18部分用了更好的6-311G*基组（在//www.umsyar.com/584中也指出6-31G*无法合理描述C18）。对相应单体也在对应的级别做了优化和振动分析。之后通过ORCA程序使用ωB97X-V/def2-TZVP级别结合counterpoise校正计算了它们的结合能（做法可参考《在ORCA中做counterpoise校正并计算分子间结合能的例子》//www.umsyar.com/542），发现C18 + OPP以及C18 + C18@OPP的结合能分别为-18.4和-18.5 kcal/mol，一方面说明C18与OPP有较强的内在结合能力，另一方面说明C18@OPP的已进入的一个C18几乎不影响OPP对第二个C18的结合能力。之后，文中通过《使用Shermo结合 程序方便地计算分子的各种热力学数据》（//www.umsyar.com/552）介绍的Shermo程序计算了C18、OPP、C18@OPP和2C18@OPP常温下自由能热校正量，并进而得到了常温下的C18 + OPP和C18 + C18@OPP的结合自由能，分别为-6.0和-4.0 kcal/mol，明显为负值说明C18的进入是热力学上可自发的，其大小明显小于结合能是来自于分子间复合带来的熵罚效应。

审稿人问及C18是否有可能与OPP存在其它结合方式。为了说明这一点，此文的补充材料里给出了优化出的其它两个C18与OPP复合物的极小点结构，如下所示，相对于C18平行地嵌入在OPP大环内的结构的能量差也在下图给出了。可见其它两种结构都是能量更高、明显更不稳定的。在《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）介绍的PCCP文章中实际上也已经通过分子动力学模拟证明了C18没有其它的能够与OPP稳定结合的方式。

文中进一步研究了C18的旋转势垒。对C18@OPP和2C18@OPP优化的C18旋转的过渡态结构如下(a)所示，可见虚频数值非常小，暗示旋转势垒肯定很低。图中的箭头体现了虚频方向，可见明显对应的是C18的整体旋转。此图是按照《在VMD中绘制Gaussian计算的分子振动矢量的方法》（//www.umsyar.com/567）介绍的方式绘制的。

进一步，文中对C18在OPP大环中的旋转进行了扫描（计算级别同几何优化），如上图(b)所示，其中横坐标是旋转角度，0处为过渡态结构。可见旋转势垒非常低，只有0.13 kcal/mol，这体现C18在OPP大环中的旋转必然很容易靠热运动来驱动。上图横坐标范围对应一个旋转周期，是40度，这也正对应于18碳环的极小点结构是D9h的特征。图中看上去每20度有一个极大点，这是因为18碳环是18个原子构成的环状体系，如果忽视其长-短键交替特征的话，就相当于20度一个旋转周期。

4 分子动力学模拟C18在OPP中的运动行为

分子动力学模拟研究是本文最关键的部分，模拟出的动力学行为是C18与OPP复合物能否作为分子马达的决定性证据。此文使用GROMACS程序通过经典力场进行动力学模拟，使用灵活的sobtop程序（//www.umsyar.com/soft/Sobtop）构建拓扑文件，主要基于GAFF力场，部分参数利用OPP的Hessian矩阵由sobtop产生，相关细节在《8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究》（//www.umsyar.com/674）里有描述，这里就不多提了。

2C18@OPP处于300 K的平衡状态时的200 ps动画如下所示（是在文章的补充材料里提供的）。为了能让C18的旋转特征看得清楚，其中一个原子以红色来标记。由动画可见，的确C18在OPP中能够非常顺利、自如地转动起来，充分证实了此文提出的双马达超分子体系的设想！

//www.umsyar.com/attach/684/2C18_OPP_rotator.mp4 （点击查看动画）

需要注意的是，由于原子间频繁碰撞导致的动能交换，单靠热运动驱动的C18在OPP中的旋转并非始终是单向的。300 K下模拟的50000帧（0.2 ps保存一帧，故相当于10 ns）中的2C18@OPP中两个C18的旋转角速度如下所示（数据已通过Savitzky-Golay算法利用相邻1000个数据点做了平滑化处理以消除噪音）。由图可见C18的含符号角速度不断发生变化，时正时负，时大时小，因此C18的旋转不能视为是连贯、无摩擦的。从图上两条曲线也可以看到2C18@OPP的两个C18彼此间的运动并没有什么相关性。由于它们离得远，自然也不会有什么明显的直接的耦合或者藉由OPP的结构产生的间接耦合。

5 分子动力学模拟C18在OPP中转动的统计行为

文中在50、100、200、300、400 K的热浴下都做了100 ns的分子动力学模拟，并通过自写的轨迹分析程序考察了C18在OPP中转动的统计行为，结果如下表所示。在温度不是很低情况的模拟过程中，偶尔C18会倾斜地倚靠在OPP的内侧，此时的体系明显不算是转子。这种情况的帧数占总模拟帧数的百分比是下表的p_out，这些帧不纳入C18转动行为的统计。可以看到从200 K开始，随着温度上升，处于这种亚稳状态的比例随之上升。

上表中ω是平均无符号转动角速度，可见随着温度越高、体系中原子的平均动能越大，C18平均旋转速度也随之上升。但是300 K变化到400 K过程中平均旋转速度增大甚微，这主要在于400 K的时候C18在OPP中的状态已经很不稳定了，过强的热运动很大程度影响了C18转动的有序性。上表中f=ω/(2π)是换算出来的转动频率，2C18@OPP属于亚THz范围的超快转子。

上表中αavg、αSD分别是C18环平面与OPP大环平面之间在模拟过程中的夹角的平均值和标准偏差，d_avg和d_SD分别是C18环中心与OPP环中心之间在模拟过程中的距离的平均值和标准偏差。从这些指标来看，随着模拟温度提升得越高，2C18@OPP偏离其极小点结构（α和d都近乎为0）的程度越高，越缺乏理想转子特征。

下图展现了VMD绘制的2C18@OPP在100 ns动力学过程中的轨迹叠加图，每1 ns绘制一次，根据帧号从前到后按照红-白-蓝着色。可见在100 K的时候由于热运动弱，体系结构基本上就是在极小点结构附近很小范围波动。到了更高的200 K后OPP的运动幅度以及C18的运动范围都显著增加了。在300 K的时候结构波动更大，还明显看到C18都偶尔侧贴在OPP大环内侧了。到了400 K的时候C18在OPP大环里的运动其实已经很大程度算是无序了，甚至仔细看轨迹都会发现C18在大环内出现了整体翻转。可以预期如果温度明显更高的话，C18甚至就能从OPP的大环中跑出去了，比如临时性跑到OPP外侧乃至飞走。

上图的(b)和(c)对不同温度下模拟的C18在OPP中的结构状态做了统计分析，包括C18与OPP大环间的夹角以及中心距离，可见这两个衡量C18在OPP中结构状态的参数的分布都是随着温度增加而显著变宽，也即在模拟过程中C18的朝向和位置的波动程度愈发增大，这和轨迹叠加图展现的信息是一致的。

6 总结

本文介绍了近期发表于Chem. Commun., 59, 9770 (2023)的理论设计的双转子超分子复合物结构，此工作从能量和动力学角度对体系进行了全面的研究，充分证明了此体系作为超快纳米转子的能力，并对其转子性能以及温度依赖性做了系统的考察。由于OPP和C18都已经被合成出来，因此此体系或者其变体有望成为有实际价值的分子机器的组成部分，本文的研究思想和发现对于纳米转子类体系在未来的探索也有显著的启示作用。

Distinguishing tasks and theoretical methods in computational chemistry

文/Sobereva@北京科音

First release: 2023-Aug-5    Last update: 2024-Sep-18

0 前言

在网上答疑时，偶尔看到有初学者搞不清楚任务类型和理论方法，比如今天有人在思想家公社3号群问“md，aimd和qmmm的区别是啥啊？什么需求下会用到这三种呢？”，这三个词明显都不是一个层面的东西。此文就把计算化学中的任务类型和理论方法的关系明确一下，简明扼要地介绍一下相关基本概念，做一个科普，以令计算化学零基础者一次性搞清楚它们的关系。更多的计算化学的总览性的知识在北京科音初级 培训班（http://www.keinsci.com/workshop/KEQC_content.html）里有介绍，此培训也是从零开始入门 研究的极好方式。

1 任务类型

有N个原子的非线型的体系有3N-6个内坐标描述原子间相对位置关系，在Born-Oppenheimer近似下体系的能量是依赖于内坐标的，也即能量是个3N-6维的函数，这被称为势能面（potential energy surface, PES）。计算化学领域有很多常见的任务（task）都是基于势能面做的，即有了势能面信息（能够求得势能面上任意位置的能量及其导数）就能做这些任务。下面罗列一下常见的这种类型的任务：

• 优化极小点：平时说的几何优化（geometry optimization）一般也是指这个。此任务从一个给定的初猜结构开始，根据特定算法去寻找与之最近的势能面上极小点的精确位置。在分子动力学程序如GROMACS里这种任务也被叫做能量极小化（energy minimization, em），只不过实际目的不一样，em的目的主要是在动力学模拟之前释放体系中可能存在的较大斥力（自行构建的初始模型里往往有原子间距离太近、斥力太大），免得动力学模拟一开始由于过大的原子加速度造成过大的速度而导致动力学模拟异常或崩溃。

• 优化过渡态：从一个给定的初猜结构开始，根据特定算法去寻找与之最近的势能面上的一阶鞍点的精确位置。过渡态可以视为反应过程中最有代表性的一个结构，可以由此判断或验证反应机理，利用它和相邻极小点的能量求差可以得到反应势垒，对于讨论反应难易有关键意义，见《谈谈如何通过势垒判断反应是否容易发生》（//www.umsyar.com/506）。优化过渡态的算法介绍见《过渡态、反应路径的计算方法及相关问题》（//www.umsyar.com/44）。

• 产生反应路径：用于把基元反应对应的势能面上两个相邻极小点之间最容易相互转变的路径产生，也相当于得到了一个轨迹并可以观看对应的结构变化过程，过渡态是此路径上能量最高点。在质量权重坐标下产生的反应路径称为IRC（intrinsic reaction coordinate），不考虑质量权重时产生的一般称为MEP（minimum energy path）。具体算法很多，见《过渡态、反应路径的计算方法及相关问题》（//www.umsyar.com/44）。反应路径可认为是实际化学反应过程最有代表性的路径，自然对于理解反应机理有重要意义，还可以用Multiwfn对反应路径上各个点做波函数分析考察反应过程中电子结构的变化以获得更深入的信息，例如《通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程》（//www.umsyar.com/200）。

• 振动分析：通常是在势能面上的驻点（所有原子受力都为0的点）做的，用于得到相应结构下的振动频率，可以用来计算热力学量，公式看《使用Shermo结合 程序方便地计算分子的各种热力学数据》（//www.umsyar.com/552）介绍的Shermo程序手册的附录部分，还可以用来得到振动能级、振动波函数、检验几何优化的准确性（根据虚频判断）、绘制振动光谱（参看《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》//www.umsyar.com/224）。

• 分子动力学（molecular dynamics, MD）：上面的任务都是“静态”的任务，即不含时间这个维度。而分子动力学则引入了时间这个维度，可以模拟体系状态随时间的变化，能研究的问题与那些静态的任务有极大的互补性，在北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/workshop/KGMX_content.html）里面有全面的介绍。分子动力学模拟过程相当于体系在势能面上不断运动，也相当于是对势能面的一种采样方式。虽然其名字叫分子动力学，但研究对象并不限于分子，比如无机固体、金属团簇的动力学模拟也可以叫分子动力学。最常见的分子动力学的实现形式是BOMD（结合 方法时另有CPMD、ADMP等），也相当于根据原子的位置、速度和受力按照经典牛顿运动方程演化原子的坐标和速度。还有其它一些特殊的动力学形式，如朗之万动力学、耗散粒子动力学等，如今用得很有限。

• 蒙特卡罗：和分子动力学并列的另一种常见的对势能面采样的方法，适用场合相对更少，主要是在多孔材料吸附小分子、计算气液共存曲线等问题上用得较多。不需要像MD那样计算受力和速度信息，也没有明确的时间的概念、无法像MD那样严格地考察含时演化。

• 构型/构象搜索：势能面上往往有很多极小点，能量最低的那个是全局最小点，可以视为体系最稳定的构型或构象对应的结构，其它的极小点对应亚稳的构型/构象。前面说的几何优化只能收敛到与初猜结构最近的极小点，显然这未必是最小点。如果想确保得到最小点（全局优化，global optimization），或者能量最低的一批构型/构象，就需要做构型/构象搜索，方法很多，比如免费的molclus（http://www.keinsci.com/research/molclus.html）就是很常用的实现工具。

2 理论方法

上一节介绍的那些任务在进行过程中，至少需要计算体系在不同结构下的能量。有的还需要计算受力（势能面对核坐标一阶导数矢量，或者说梯度，的负矢量），如几何优化、分子动力学。有的还需要计算能量对核坐标的更高阶导数，比如振动分析至少需要算二阶导数（Hessian矩阵由之构成）。怎么获得各种任务中涉及的能量及其对核坐标的导数，就是理论方法层面的问题了。下面简要列举一些常见的理论方法：

• 分子力场（molecular forcefield）：也称为经典力场（classical forcefield）或分子力学（molecular mechanics, MM），或者经验势函数等等。虽然通常带着“分子”俩字，但实际中不限于用于分子体系。这类方法一般使用简单的模型描述体系（通常一个原子作为一个粒子，电子与原子核不分开考虑），并用形式简单的数学函数（势函数）描述原子间相互作用。例如计算原子间静电相互作用时，分子力场普遍是给每个原子核位置分配一个点电荷（数值等于原子电荷），然后基于库仑公式计算。由于分子力场的形式简单，因此计算耗时极低。分子力场的一个关键的局限性是绝大多数分子力场都无法描述化学反应，因为成键关系在计算从始至终是固定不变的，是一开始就定义好的。描述化学反应问题主流的是后面提及的 类型的方法，或者反应力场（远比普通分子力场要贵、支持的程序要少）。分子力场另一个关键局限性是依赖于大量参数，一方面影响模拟精度，一方面决定能描述的体系，寻找和构建合适的参数往往不是易事。

• 机器学习势（machine learning potential）：基本思想是通过机器学习的思想构造依赖于原子坐标的分子描述符（如距离矩阵、能量矩阵等）与能量之间的抽象关系，这种关系比上述传统的分子力场用的势函数形式明显更为复杂。

以下提及方法都是 （quantum chemistry, QC）范畴的，也可以叫基于量子力学（quantum mechanics, QM）的方法。

• 密度泛函理论（DFT）：是目前在 、第一性原理计算领域用得最普遍的理论方法，因为其性价比非常高，交换-相关泛函（影响DFT计算精度的决定性因素）选择得当的话可以得到很不错的结果，足够满足绝大多数研究需要，而且精度显著强于分子力场（除非力场是对某类体系专门训练得特别精妙的力场），普适性强得多得多，但耗时高于分子力场N个数量级，体系越大N越大。B3LYP、ωB97XD、PBE0、TPSS等等都是常见的交换-相关泛函，各种交换-相关泛函的流行程度见《2024年计算化学公社举办的计算化学程序和DFT泛函的流行程度投票结果》（//www.umsyar.com/706）中的得票情况，泛函的选用建议见《简谈 计算中DFT泛函的选择》（//www.umsyar.com/272）。很多泛函如B3LYP描述色散作用糟糕或完全失败，导致这些泛函不适合研究弱相互作用，但可以通过加上色散校正很好地解决，参考《谈谈“计算时是否需要加DFT-D3色散校正？” 》（//www.umsyar.com/413）和《DFT-D色散校正的使用》（//www.umsyar.com/210）。值得一提的是，DFT应用到实际问题中实际上分为Kohn-Sham DFT（KS-DFT）和orbital free DFT两种形式，前者是绝对的主流，计算过程牵扯到轨道波函数，本文以及大家平时说的DFT计算都是用的KS-DFT形式；而后者非主流，不牵扯到轨道，应用的局限性非常大，只有很少数程序支持。

• Hartree-Fock（HF）：在DFT兴起之前用得很普遍，如今已经完全过时。耗时不比DFT显著更低（DFT计算时利用RI技术时反倒还可以明显更快），而精度差得多。

• 后HF：是一大类方法的统称，如CCSD(T)、MP2、CISD等。HF方法精度烂在于同时没考虑到动态相关和静态相关。后HF在HF基础上额外把动态相关在一定程度上考虑了进来，但耗时也高了很多。动态相关、静态相关的相关概念很有学问，在北京科音中级 培训班（http://www.keinsci.com/workshop/KBQC_content.html）里有详细讲授，这里就不多提了。

• MCSCF：主要弥补HF缺乏对静态相关的考虑。但由于几乎没有或很少考虑动态相关，所以结果也说不上好。MCSCF最常见的实现是CASSCF。

• 多参考方法：在MCSCF基础上进一步把动态相关考虑进来，精度整体很好，普适性很强，但使用相对复杂且昂贵。典型的实现如CASPT2、NEVPT2、MRCI。

• 半经验：一大类方法的统称，都是对HF的简化以巨幅降低耗时，耗时只是HF的微小零头，整体精度也有不小降低，而且引入了依赖于元素的参数而只能用于有限的元素。典型的如AM1、PM3、PM6。一些较复杂的机器学习势的耗时已经追平了半经验。

• GFN-xTB：相当于纳入了一定DFT思想的半经验级别的方法，有不同版本，其中较常用的GFN1-xTB和GFN2-xTB耗时和半经验相仿佛，而整体精度和可靠性>=主流的半经验方法，但跟DFT比还有很大差距。GFN-xTB从2017年开始发展，和历史更早的DFTB有极大的共性。

前面介绍的HF、后HF、MCSCF、多参考方法里面都完全没有引入经验参数（极个别的除外，诸如CASPT2可以涉及shift参数），计算中的参数只依赖于物理常数，故它们被称为“从头算（Ab initio）”类型的方法。分子力场、半经验、GFN-xTB由于引入了大量经验参数，因此明显不属于从头算。DFT算不算从头算模棱两可。如今尚未找到的理论上严格的交换-相关泛函是不含任何参数的，但如今只有近似的交换-相关泛函被提出，其中虽包含经验参数，但整体相对较少，因此争论DFT是否属于从头算没意思，只是个说法问题。有的泛函的经验性较小，有的经验性则较高。如今也有的文章将HF、后HF、MCSCF、多参考方法等物理思想完全基于波函数的方法统称为wavefunction theory (WFT)，在讨论时使用WFT一词时DFT就不会和那些方法混淆了。

还有一个词叫“第一性原理（first-principles）”，其字义和从头算一样，但“第一性原理计算”如今的主流指代是“ 方法研究周期性体系”，不仅包括严格的从头算方法，也包括DFT，因为这类计算最常用的就是DFT。DFTB虽然依赖很多参数，但DFTB计算周期性体系也往往被不计较地算作第一性原理计算。

由于有了“第一性原理计算”这个专门描述把 方法用于周期性体系的计算的词，因此平时说的“ 计算”、“从头算方法计算”一般是指计算分子、团簇这样孤立（非周期性）体系的情况。那些通过 方法主要研究孤立体系的程序被称为 程序，如Gaussian、ORCA、GAMESS-US等，而通过 方法主要研究周期性体系的程序被称为第一性原理程序，如CP2K、Quantum ESPRESSO等。注意“量子计算”和 计算又不是一码事，前者强调的是利用在量子计算机上实现的量子算法，可以但绝不限于做 问题的研究。

接下来再说很知名的QM/MM。这个方法是指使用 方法以量子力学（QM）的思想描述体系的一部分，用分子力场以分子力学（MM）的思想描述体系的另一部分，并且恰当考虑两部分之间的耦合。对于很大的体系，这显然比起全都用 方法描述要便宜得多得多（MM部分计算耗时相当于QM部分的九猪一毛），但代价是MM区域描述的精度比QM区要差，而且绝大多数力场不能描述这部分发生的成键/断键过程。因此，必须恰当划分QM和MM区，让最需要精确描述的部分或者涉及化学反应的原子纳入到QM区中，而起到较次要作用的“环境”原子纳入到MM区中。

还有一个词是ONIOM，也是划分不同区域，不同方法着重描述不同部分，它可以把任意两种理论方法相组合。如果把 方法和分子力场相组合，特称为ONIOM(QM:MM)，其用处和QM/MM类似但实现不同，描述时绝对不能搞混。一些相关讨论见《要善用簇模型，不要盲目用ONIOM算蛋白质-小分子相互作用问题》（//www.umsyar.com/597）。

再说一下对电子激发态的描述。前面介绍的半经验、HF、后HF、DFT等方法通过delta-SCF方式可以算激发态，但不够普适。 领域有很多专门的方法可以计算激发态，包括计算激发态的能量及其导数，如最常用的TDDFT、过时的CIS、较高精度的EOM-CCSD等，见《乱谈激发态的计算方法》（//www.umsyar.com/265）。几乎所有分子力场都是描述基态电子态的，但如果专门对激发态势能面拟合力场参数，也可以用力场描述激发态，但已有的这样的力场都是针对极其具体的体系的激发态拟合的，不具备普适性。机器学习势同理。QM/MM也可以描述激发态，但仅限于电子激发完全发生在QM区的情况。

再简单说一下和实际计算关系密切的基函数（basis function）的概念。本节介绍的理论方法，凡是基于量子理论思想提出的（也就是前面说的那些里面除了分子力场和机器学习势以外的方法），在实际数值求解的过程中一般都要涉及分子轨道，关于分子轨道的概念详见《正确地认识分子的能隙(gap)、HOMO和LUMO》（//www.umsyar.com/543）。绝大多数计算程序中都是把分子轨道展开成基函数的线性组合来描述的。最常见的基函数一类是原子中心基函数（如Gaussian等大多数 程序以及CP2K等部分第一性原理程序用的高斯型基函数），其中心一般位于原子核，还有一类常见的基函数是平面波基函数，是大多数计算周期性体系为主的第一性原理程序用的，它的分布覆盖整个被计算的晶胞。基组（basis set）是对于原子中心基函数而言的，例如6-31G*、def2-TZVP、cc-pVDZ等都是很常见的基组，它定义了实际计算时对各种元素原子具体用多少、什么参数的基函数。做HF、DFT、后HF、MCSCF、多参考等方法计算时都需要告诉计算程序用什么基组，基组质量越好，也就是越接近于所谓的完备基组极限，这些方法自身的精度发挥得就越充分，但代价就是耗时越高。一个好方法配一个烂基组，以及一个烂方法配一个好基组，结果都不理想，必须好方法配好基组才能得到较准确结果。半经验、GFN-xTB方法计算时也利用基函数，但是用什么基函数是这类方法自身定义的，就不需要而且也不能再指定基组了。北京科音初级 培训班（http://www.keinsci.com/workshop/KEQC_content.html）和北京科音中级 培训班（http://www.keinsci.com/workshop/KBQC_content.html）都专门有一节对基组做详细介绍，后者讲得明显更深、更广。

3 任务类型与理论方法的结合

第2节介绍的各种理论方法原则上都可以用于第1节介绍的各种任务。例如，分子力场、DFT、QM/MM等理论方法全都可以用于几何优化、振动分析、产生反应路径等任务。这些组合在程序实现上一般没有什么需要特殊考虑的，一个程序支持什么理论方法，就能基于它们做什么依赖于势能面的任务。

搞明白了前面介绍的名词，从头算动力学（ab initio molecular dynamics, AIMD）和第一性原理动力学（first-principle molecular dynamics, FPMD）顾名思义自然就知道是描述什么计算了。这两个词其实没有严格的区分，实际中经常混淆使用，但从严格的习俗来说，AIMD和FPMD分别最适用于描述使用 方法对孤立体系和周期性体系做分子动力学模拟。所以北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里面讲CP2K做基于GFN-xTB和DFT的分子动力学的地方标题用的是第一性原理动力学，但若说成AIMD也完全可以，不会造成什么误会。

由于基于分子力场做分子动力学的文章数目和流行程度远高于非常昂贵的AIMD/FPMD，因此平时说的分子动力学模拟默认是指的基于经典力场的分子动力学，这也往往被叫做经典分子动力学。但要注意“经典”这个词在不同语境下有不同含义，有的场合下是特指把原子核当做经典粒子来考虑，此时基于BOMD形式做基于 方法的分子动力学也是经典分子动力学。而与之相对的是“量子动力学”，强调把原子核以量子力学方式描述，原子核有对应的波函数。显然，此处“量子动力学”和“基于 做分子动力学”又不是同一个概念。

QM/MM MD指什么也不言而喻，显然是用的QM/MM理论方法做分子动力学任务，这结合了QM/MM对体系描述的特点和分子动力学研究问题的能力。显然这比所有原子都用QM描述（相当于AIMD/FPMD）要便宜得多得多，当然代价就是MM区域在动力学过程中描述得相对较糙，而且预期不会发生化学反应，同时还得想办法去找力场参数。因此能否用QM/MM MD代替AIMD/FPMD当然得根据被研究的体系和问题判断，盲目瞎用纯粹是自讨苦吃（吐槽一下，很多初学者没有基本理论常识，一看文献里有QM/MM就觉得好像高大上似的，就整天想着QM/MM而不知道其局限性，太naive了）。

文献里的许多用词并不标准，必须结合语境和常识理解说的是什么。比如可能有的文章没明确用QM/MM MD这个词而只说了QM/MM，若实际上在此基础上跑了动力学，当然就是QM/MM MD。再比如有的文章表面上说是做的AIMD，但从计算细节描述来看作者把一部分当做MM来描述以节约时间，显然实际上也是QM/MM MD。理解名词必须有基本的应变能力。

如果用的理论方法描述的是激发态，那么前述各种任务就是对激发态做的。比如用TDDFT做激发态分子动力学，那就是动力学的每一步用的原子受力都是用TDDFT算的激发态的受力。顺带一提，有些文献里声称做的是激发态动力学研究，但其实内容里根本没做分子动力学，只是对激发态计算了势垒，文中这么说只是由于势垒和研究反应速率的反应动力学领域里的反应速率常数密切相关。另外值得注意的是研究激发态动力学往往需要考虑内转换、系间窜越等势能面切换的过程，考虑这些需要做非绝热动力学，需要额外的算法，这方面的坑很深。所以并不是一个程序支持了计算激发态的理论方法和分子动力学模拟后就能完美、严格地做激发态动力学研究相关问题，情况远没那么简单。

On the possible reasons why data in computational chemistry literature cannot be reproduced

文/Sobereva@北京科音   2023-Aug-1

0 前言

经常有做计算化学的人在网上问为什么他没法重复出来文献里的计算数据。导致重复不出来的原因实在太多了，然而初学者提问此类问题时几乎总是含糊其辞，提供的有效信息甚少，导致别人根本无从回复，或者需要别人罗列一大堆可能导致差异的原因。我感到每次回复此类问题太费事，于是在这里专门写个文章说说这事。遇到这种问题的读者仔细看完此文后，应该能料想出可能是什么原因导致自己和文献的结果存在差异。就算还搞不明白，至少也应该能明白在网上问别人的时候应当提供什么信息，知道怎么把导致差异的可能原因的范围尽可能缩小，免得别人打很多字罗列诸多当前不需要考虑的可能性，或者别人嫌你提供的有效信息完全不够而根本不想回复。

这里首先要强调一点的是，如果你的目的就是为了精确重复文献数据，那你就要尽最大可能使用和文献完全相同的方式去做计算，而不要管人家的做法合不合理，哪怕你明知道人家的计算方式不当。下文说的内容都是为了让你能尽可能准确重复出文献里的数据，而不是讲怎么合理地计算。如果你的实际目的是想获得相同体系相同研究问题的合理的结果，那你只要确保你自己的计算是合理的就行了，没必要非得把别人的数据重复出来，说不定人家的计算本来就有问题。如果最终验证发现别人的计算有硬伤并因此导致结论错误，或者凭基本知识就能断定别人的计算有问题，那么可以发一篇comment文章指正以免以讹传讹，比如《我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment》（//www.umsyar.com/584）里介绍的文章。

1 结构不一样

自己计算用的结构和文献不一致，是导致数据无法重现的最常见原因之一。如果你用的结构就是文献正文或补充材料里给的坐标，那就可以直接排除这个因素了。注意有些文献的SI里键长、笛卡尔坐标单位是Bohr而不是埃，千万别弄错了。

如果你计算的基本结构特征都和文献里不符，也即算的都不是同一个东西，能重现出文献里的数据就怪了。很多体系有不同的形态，比如酚酞在酸性、中性、碱性环境下主要的存在形态明显不一样，氨基酸在真空和水环境下的质子化态明显不一样，姜黄素在水和有机溶剂中一种是烯醇式一种是酮式结构，等等。你必须确保你算的和文献里算的确实是同一种状态的结构，如果文献里给了截图，一定要仔细对照。

做计算要有基本的化学常识，并且要认真仔细，免得搭出来没意义的结构去做计算，自然和文献里对不上。比如文献里给出Lewis结构式时通常是隐去碳上的氢的，如果计算的时候真不画氢，或者少画了个别的氢，自然计算毫无意义。再比如，文献里算一个羧基配位的过渡金属配合物，如果你在计算时羧基上还留着氢，明显也是错的。再比如，文献里画了一个带一个小配体的过渡金属配合物的Lewis式，若你算的是水溶液中的状态并忽略了配位水分子，也自然是错的。再比如，文献里给出一个联苯或者碗烯的二维结构式，看上去是在一个平面上，若你去计算其基态时真的摆成平面来优化，得到的也是大概率是错误的结构（有破坏平面的虚频的结构）。还要注意一些结构的构建不要想当然，比如SiO2表面在水中时表面的O主要是以羟基及其去质子化形态存在的，比例和pH有关，如果都当成是乌秃秃的氧就错了。结构构建的关键性信息在文献里一般都会有具体说明，但也可能一些文献中在描述上有疏漏。

有些体系在计算的时候结构构建方式不唯一，比如《要善用簇模型，不要盲目用ONIOM算蛋白质-小分子相互作用问题》（//www.umsyar.com/597）、《18碳环（cyclo[18]carbon）与石墨烯的相互作用：基于簇模型的研究一例》（//www.umsyar.com/615）和《使用 程序基于簇模型计算金属表面吸附问题》（//www.umsyar.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之，需要注意文献里是怎么考虑的。蛋白质的高解析度和低解析度测出来的晶体结构里残基侧链的位置也可能相差甚大，甚至明显影响到做酶催化、分子对接之类问题的结果。同样的属性，以不同类型实验测定的结果注定也会有差异，比如气相电子衍射和单晶衍射测的一个是气相结构一个是晶体结构，对于柔性分子差异可能很大，这点参考《实验测定分子结构的方法以及将实验结构用于 计算需要注意的问题》（//www.umsyar.com/569）。

上面列举了诸多结构层面影响计算结果的可能性，当然不可能列举得很完备，实际情况多种多样。总之希望读者能意识到严格保证和文献里用的结构的一致对于重现数据有多么重要。

2 计算设置不一样

为了精确重复文献数据，要尽可能仔细阅读文章及其补充材料，了解清楚文献计算细节，使得自己的计算设置能最大程度与文献相一致。本节里谈到的很多可能导致差异的因素在《谈谈不同 程序计算结果的差异问题》（//www.umsyar.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结合 程序方便地计算分子的各种热力学数据》（//www.umsyar.com/552）里介绍的Shermo程序默认用的以及ORCA用的是Grimme的QRRHO模型，而Gaussian用的是较原始的不适合含低频体系的RRHO模型
• 计算热力学量时是否考虑了平动和转动，显然算周期性体系时不应当考虑（在Shermo程序中这由settings.ini里的imode决定）
• 计算转动对热力学量贡献时对转动对称数的考虑，看Shermo手册附录部分了解相关常识
• 计算热力学校正量、振动分析用的原子质量
• k点数目和分布
• 平面波相关计算的截断能
• CP2K中，做杂化泛函时用的库仑截断半径、DFT+U计算用的布居计算方法以及考虑的角动量、Poisson solver的选择
• 对相对论效应的考虑
• 对外场/外势的考虑
• DFT积分格点精度。参看《密度泛函计算中的格点积分方法》（//www.umsyar.com/69）
• 是否用了冻核、怎么冻的
• 双电子积分的积分屏蔽设置
• 是否用了RI、DLPNO等加速计算方式
• 是否考虑了偶极校正
• QM/MM计算划分的位置、MM与QM区之间的相互作用的考虑、link原子的设置
等等等等

基组必须和文献严格一致，并注意上述//www.umsyar.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中的计算》（//www.umsyar.com/441）、《RESP2原子电荷的思想以及在Multiwfn中的计算》（//www.umsyar.com/531）。

做动力学模拟的情况，必须确保模拟相关设置和文献一致。比如步长、步数、坐标/速度/受力/能量保存频率、参考温度和压力、控温和控压方法、热浴和压浴的时间常数、可压缩系数、各向同性还是各向异性控压、初速度怎么产生的、范德华和静电作用的计算方式（如简单截断、twin-range、Ewald、PME、SPME等，以及其中涉及的截断半径，是否用了switching function等处理）、水模型是刚性还是柔性的、是否用了键/键角约束、用没用能量/压力色散校正，等等。

如果你用的计算程序和文献里不同，或者和文献里的相同但版本不同，那么由于默认设置不同，以及算法实现的不同，可能对结果带来很大差异，一定要仔细阅读《谈谈不同 程序计算结果的差异问题》（//www.umsyar.com/573）。

在 波函数分析领域也有很多细节设置会影响分析结果。这里以非常流行的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）为例：
• 按照《使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析》（//www.umsyar.com/179）、《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》（//www.umsyar.com/159）、《使用Multiwfn对静电势和范德华势做拓扑分析精确得到极小点位置和数值》（//www.umsyar.com/645）等文章介绍的方法进行分析时，格点间距数值越小精度越高，但也越耗时
• 按照《使用Multiwfn绘制态密度(DOS)图考察电子结构》（//www.umsyar.com/482）绘制PDOS和OPDOS图时，在Multiwfn里选择的计算轨道成份的方法会直接影响结果
• 按照《使用Multiwfn模拟扫描隧道显微镜(STM)图像》（//www.umsyar.com/549）和《使用Multiwfn结合CP2K的波函数模拟周期性体系的隧道扫描显微镜（STM）图像》（//www.umsyar.com/671）绘制STM图时，偏压的数值和常高模式扫描的平面的选择都会明显影响结果
• 计算《衡量芳香性的方法以及在Multiwfn中的计算》（//www.umsyar.com/176）和《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》（//www.umsyar.com/138）里介绍的多中心键级时，基于原本的基函数和基于自然原子轨道（NAO）计算会得到不同的结果
• 计算Hirshfeld原子电荷时需要用到原子自由状态下的密度。用自行提供的原子波函数文件来计算此密度，还是直接用Multiwfn内置的此密度，会有所不同。

本节提到的这些计算细节是理应充分体现在文章的computational details里或者补充材料里的。但由于一些文献的作者为了省篇幅，或者用的很多设置和程序默认的一致，或者由于作者计算水平比较外行、意识不到很多细节对计算非常重要，就在文中没特意说。如果有些要点不知道作者具体怎么考虑的，自己把各种可能性都试试就知道了，文章没特意说的设置就姑且用程序默认的。

3 计算的电子结构状态不一样

净电荷和自旋多重度的设置必须和文献里的相同，否则相当于算的是其它电子态，自然算出来的各种性质都不同。比如计算基态的氧气发生的反应、计算富勒烯的原子化能，如果氧气和碳原子误当成了单重态来算，显然结果是错的。此外，稍微了解点两态反应性的概念就会知道，不少体系在不同自旋态下反应机理都明显不一样，如果自旋多重度或净电荷设错了，自然反应物/产物/中间体/过渡态等无法得到和文献里一致的结果，而且相差悬殊。

如果文献没注明、从文献中看不出相关信息，那么净电荷一般为0。如果被研究的体系的自旋多重度不那么显而易见，或者研究中涉及到了不止一个自旋态，文中一般都会说明自旋多重度。如果确实没说，那就得查阅文献或者数据库，或者自己判断、测试了。

很容易被忽略的是开壳层单重态或者叫自旋极化单重态。用KS-DFT等方法来计算的话应当做对称破缺计算，参看《谈谈片段组合波函数与自旋极化单重态》（//www.umsyar.com/82）。比如整体为单重态的反铁磁性耦合体系、双自由基体系、共价键被拉长的均裂状态就是最典型的例子。还有一些特殊的体系也是如此，比如18碳环的D18h点群的过渡态结构、乙烯二面角扭转成90度的结构、六并苯。重复文献里的这些体系的计算时别误当成了闭壳层单重态来计算（文中如果也是做了对称破缺计算，一般都会明确说。如果没说，没准文献也误当成了闭壳层计算）。

对于第一性原理计算磁性材料时必须谨慎，很多情况下不是光设个整体的自旋多重度就完事的，而必须恰当设置初始磁矩才能收敛到真正要考察的磁性状态。为了重复文献（或材料数据库）中的数据，若它们给了收敛的波函数对应的各个原子的自旋磁矩/自旋布居信息的话，最好根据其来设置初猜波函数里的原子自旋磁矩，以尽最大可能收敛到和文献相同的波函数状态。

要注意波函数稳定性问题。即便净电荷、自旋多重度的设置是正确的，也不代表SCF收敛后得到的波函数就是实际想研究的态、和文献里算的状态相一致。如果文献里研究的是基态，你算的结果和文献对不上，而且体系又不是简单有机体系而是电子结构特殊的体系，那么应当做一下波函数稳定性检测，比如Gaussian里可以用stable关键词，发现不稳定的话可以用stable=opt尝试搞出稳定波函数。不稳定的波函数可能意味着这是某个激发态，还可能意味着此波函数是没有任何物理意义的虚假的状态。

4 计算的东西或方式不一样

重复不出文献的数据时还要想清楚到底你算的东西、计算的方式是否和文献一致、有可比性。下面列举一些典型的需要注意的情况。

比如通过能量差E(AB)-E(A)-E(B)考察俩分子A和B间相互作用强度，复合物结构肯定是要优化的，而单分子是直接从复合物优化后的结构里直接取出来，还是也做优化，这在不同文中的做法就不一样了。单分子不优化情况下得到的能量差一般称为（复合物结构下的）相互作用能，单分子也优化情况下得到的能量差一般称为结合能，二者之间相差各个单体的变形能之和。对于诸如《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（//www.umsyar.com/513）里面列举的那些刚性小分子，相互作用能和结合能差异甚微，但如果考察的是两个容易变形的分子间的作用强度，比如两个氨基酸分子的结合，那么相互作用能和结合能就必须明确区分了。J. Chem. Phys., 158, 244106 (2023)还专门讨论了什么情况变形能会较大。对于键能，也是类似地有不同计算方式。所以研究此类问题的时候一定要搞清楚文献具体是怎么算的。

再比如做SAPT能量分解计算，会得到很多耦合项，比如交换与色散作用的耦合项，这是归到色散项里，还是归到交换互斥项里，不同文章的做法可能不一样，需要注意文章怎么说的。再比如计算反应的自由能问题，是计算标准反应自由能，还是计算特定浓度下的反应自由能，这是不一样，相差浓度校正项，要搞清楚文献实际算的是什么。再比如电子亲和能、电离能、激发能，都有垂直和绝热之分，文献里计算哪种都有，也必须确保和文献算的东西对应。

文献里有时候对热力学量的标注比较含糊，要分清楚文献算的是什么温度下什么量。比如E，文献里可能指电子能量，也可能指内能，需要结合语境判断。

模拟ECD谱需要转子强度，有长度形式和速度形式两种，诸如Gaussian都会输出，在Multiwfn里按照《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》（//www.umsyar.com/224）里说的绘制ECD时可以选采用哪种，一般推荐用速度形式。如果你用的和文献不一致，会导致峰型存在一定差异。

基组的CBS外推的方法不唯一，有不同公式，有的情况外推方式里也涉及到参数，参看《谈谈能量的基组外推》（//www.umsyar.com/172）。不光是能量外推，还有比如GW计算的准粒子能量随基组的外推、几何结构随基组的外推等。一些文献对高级别理论方法结合大基组还用一些近似方式估计，比如一种常见的是CCSD(T)/大基组 = CCSD(T)/小基组 + (MP2/大基组-MP2/小基组)，这在《各种后HF方法精度简单横测》（//www.umsyar.com/378）里专门说了。涉及到这些高精度数据估算时，需注意保证和文献的做法一致。

还有一些量本来就没有唯一定义，去重复文献的数据更是得弄清楚文献里用什么定义、什么方法算的。比如自由体积，我在《使用Multiwfn计算晶体结构中自由区域的体积、图形化展现自由区域》（//www.umsyar.com/617）里介绍了一种合理的计算方式，但显然这不是唯一方式。再比如分子半径，我在《谈谈分子半径的计算和分子形状的描述》（//www.umsyar.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》（//www.umsyar.com/543）。如果没搞清楚概念和定义方式，可能重复不出文献的数据，或者实际上文献的计算方式或用词不当，你却意识不到这一点。

某些问题的计算中可能涉及到一些实验值，要搞清楚作者用的什么值。比如用 程序计算含碳物质的生成焓，通过构造热力学循环可知其中需要利用实验的（或第一性原理计算的）石墨的升华焓，要搞清楚这个量是从哪来的，具体是多少。再比如， 计算某分子的标准氧化还原势，其中利用到标准氢电极的电势，其数值有不同来源，从4.05到4.44 V都有报道，Phys. Chem. Chem. Phys., 16, 15068 (2014)建议对理论计算的情况用4.28 V。

还有一些研究涉及到一些额外的参数需要注意，典型的就是频率校正因子，见《谈谈谐振频率校正因子》（//www.umsyar.com/221）。是否使用频率校正因子，以及对某一特定级别使用什么目的的、哪个来源的校正因子，都会影响结果。

周期性体系的计算中，轨道能量的绝对数值是没意义的，不同程序用的能级零点也不一样，需要关心的只有相对数值，比如带隙。所以不要问为什么DOS图的形状和文献一致但横坐标和文献不一致。这也正是为什么文献在绘制DOS图时往往一般把价带顶或者费米能级位置设为0点，因为这样才便于公平对比。

还有些量在文献里用的习俗和你用的可能不同。比如算结合能问题，有少数文章把结合过程的能量降低量作为正值记录。再比如NBO给出的E2超共轭作用能，程序输出时把正值作为超共轭作用造成的体系能量降低，但有的文献报道E2值的时候取的是负值。

如以上列举的情况可见，算某个东西的时候必须仔细阅读文献搞清楚作者是怎么算的、用的什么定义、算的到底是什么。

5 数据的可重现性本身就差

有些数据很容易重现，比如在RHF/6-31G**级别下优化水分子的基态结构，不管用什么主流量化程序、怎么优化，得到的都是相同结果，初始结构很随意。有些数据本身就是难以精确重现，典型的就是凝聚相体系的分子动力学。由于分子动力学过程有混沌效应，而且实际计算中往往不可避免地引入随机因素，会导致文献中的动力学模拟的定量结果几乎无法严格精确重现。关于这点我专门写过文章，参看《数值误差对计算化学结果重现性的影响》（//www.umsyar.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不收敛的方法》//www.umsyar.com/665文中提到了）、CP2K算磁性体系时没写UKS结果当成了闭壳层算、Gaussian里按照特定方向加电场计算时不知道要写恰当用nosymm导致电场加的方向和期望的不符（nosymm用处见《谈谈Gaussian中的对称性与nosymm关键词的使用》//www.umsyar.com/297）、用sobtop（//www.umsyar.com/soft/Sobtop）直接产生拓扑文件后没有给里面添加原子电荷、用sobtop产生拓扑文件时误把可旋转二面角以谐振势对待、用Multiwfn做空穴-电子分析时没按照《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（//www.umsyar.com/434）说的在Gaussian输入文件里写IOp(9/40=4)，等等
• 缺乏基本知识瞎算，比如：计算小晶胞体系不知道考虑k点、算极化率时不知道要带弥散函数、B3LYP算色散主导的弱相互作用却不加色散校正、CP2K里用MOLOPT基组计算需要CUTOFF很高的Na的时候用的CUTOFF不够、界面体系用了各向同性控压而不知道应该用半各向同性控压、算气-液界面体系在垂直于界面方向开了控压结果导致真空区消失、用Multiwfn基于CP2K产生的molden文件分析前没在里面添加盒子信息（参见《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》//www.umsyar.com/651）、做Mulliken布居分析用的基组带了不该带的弥散函数、算溶剂环境下溶质自由能的时候忘了考虑标准态浓度转换问题（参看《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》//www.umsyar.com/327）、试图用后HF方法计算分子轨道（却浑然不知程序虽然给出了结果但实际上是HF部分产生的）、用背景电荷表现环境分子的势场用的是对静电势重现性极烂的QEq电荷（相关讨论看《基于背景电荷计算分子在晶体环境中的吸收光谱//www.umsyar.com/579），等等
• 计算流程/方式不对，比如：算高精度能量和性质之前结构没经过几何优化、直接把IRC两端的结构当做反应物和产物的准确结构用于计算能量和属性数据、振动分析和IRC计算之前结构没有在严格相同级别下优化过、对动力学轨迹计算某个分子的RMSD衡量内部结构变化的时候不知道先要消除整体平动和转动、算电子发射光谱之前没优化激发态、模拟拉曼光谱时没按照《使用Multiwfn计算特定方向的UV-Vis吸收光谱》（//www.umsyar.com/648）里说的将拉曼活性转化成拉曼强度
• 数据读取错误。比如Gaussian里做了双杂化泛函计算，结果读成了其中SCF Done后面的中间量（怎么正确读看《谈谈该从Gaussian输出文件中的什么地方读电子能量》//www.umsyar.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》（//www.umsyar.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)，这文献的数据的可靠性没法令人不放心（虽然作者可能检查了最终收敛情况，但谁也不敢说作者确实这么做了）。现在有很多无良代算机构，没一丁点搞理论计算的人的基本操守，为了能给出令做实验的人满意的与实验相符的计算数据来交差，极有可能在算不出期望的数据的情况下捏造假数据，这样的数据更别指望能被别人重复出来了。甚至有的代算方连计算都不计算，直接凭空胡编乱造数据，建议读者看看《谈谈我对计算化学代算的看法》（//www.umsyar.com/505）里提及的情况。

有很多文献对计算的关键要点、直接影响数据可重复性的因素交代得非常粗糙、简略、语焉不详，甚至根本不提，这无疑给读者带来了很多误解、给别人重复计算造成了障碍。重复这些文章的数据只能是考虑各种可能情况去试了。另外，还有可能文章里用的参数、设置本来就不慎写错了，作者写的程序本身存在bug等等。

对于已经恰当考虑了前面说的所有因素，竭尽全力去重复文献，但还是重复不出来（且确实是有不可忽视的差异而不是在复现精度上过度较真），如果没有绝对必要去重复文献，那就别去重复了，极有可能文献的数据就是有问题的，或者没交代关键性计算细节。这种情况只要你能保证自己的计算是合理的就够了。但如果有特殊原因就是非要重复文献不可，显然该做的不是在思想家公社QQ群或者计算化学公社论坛（http://bbs.keinsci.com）等公开场所提问，因为找谁都不如直接找作者，世界上还有谁比作者更清楚数据是怎么来的的人么？世界上还有谁比作者更有回复你的义务？发邮件给通讯作者，向作者索要相应数据的输入输出文件，或者向作者提供你复现文献计算的输入输出文件问作者为什么没重复出来即可。当然很有可能一直没收到作者的回复，可能性有这些：
(1)作者看漏了你的邮件，或者邮件被误归入垃圾邮件，或者正忙着什么事没来及回复你结果后来又忘了。可以过一个礼拜（且可以换个邮箱）再发一次试试。如果通讯作者有多个，也可以给另外的人发
(2)措辞不当，不够礼貌和诚恳，而作者脾气又大或者又比较忙，又觉得回复你也没什么好处，就不回了。可以过一阵重新措辞再发一次。为了让作者有动力花时间给你回复，最好表示你会在未来的文章中会引用他的文献。
(3)文献中的数据是作者找别人代算的，代算方又是非常不负责任的，相关文件也不提供、计算细节也不交代清楚，甚至代算之后就再也联系不上了，因而作者也无能为力。
(4)做计算的是学生，学生毕业后找不到人了，导师又不清楚计算细节、没原始文件。
(5)文献中的计算本身就有严重错误或是假数据，可能收到你的邮件后作者才意识到，又或者可能之前就已经明知道用的是假数据，总之都不希望你知道他的数据存在问题，免得面子上挂不住，甚至怕文章最后被撤稿。

8字形双环分子对18碳环的独特吸附行为的 、波函数分析与分子动力学研究

Quantum chemistry, wavefunction analysis and molecular dynamics study of the unique adsorption behavior of 8-shaped bicyclic molecule to cyclo[18]carbon

文/Sobereva@北京科音   2023-Jun-30

0 前言

18碳环于2019年首次在凝聚相实验中被观测到后，其独特的几何和电子结构引发了学术界的巨大关注，笔者对18碳环及其衍生物做了大量、系统的理论研究工作，论文汇总以及深入浅出的介绍文章见//www.umsyar.com/carbon_ring.html。2022年Angew. Chem., 134, e202113334 (2022)新合成出一种具有两个大环的整体看上去像8字形的分子，每个大环都由一串苯环相连构成，此体系以下简称oligoparaphenylene (OPP)，如下图最下侧的结构所示。此文实验发现OPP可以稳定结合两个C60或C70富勒烯。另外，由下图所示意的，11个苯环构成的环状分子[11]CPP也已知可以吸附C70。

考虑到18碳环的直径和C70富勒烯相近（如上图所标注的），OPP是否也有可能吸附18碳环分子来实现18碳环的富集？此外，18碳环本身不稳定、易反应，若利用OPP吸附来保护住的话，还有可能实现18碳环的稳定化，无疑这很有实际意义。受到这个想法的启发，笔者和江苏科技大学的刘泽玉等人通过 、波函数分析和分子动力学对OPP吸附18碳环的可能性从各个角度做了全面深入的探究，并在近期发表了通讯文章，非常欢迎读者们阅读和引用：

Molecular assembly with a figure-of-eight nanohoop as a strategy for the collection and stabilization of cyclo[18]carbon, Phys. Chem. Chem. Phys., 25, 16707 (2023) https://doi.org/10.1039/d3cp01896b

在下文中，笔者将对这篇研究工作的主要内容和研究思想进行深入浅出的介绍，还将同时还做许多扩展讨论和细节说明，图片来自于上文以及其补充材料（有的图略有差异，版本不同）。希望本研究能对读者研究类似问题有所启发、将本文的研究手段应用于相似问题的研究上。如果大家不了解18碳环的基本特征的话，建议阅读《谈谈18碳环的几何结构和电子结构》（//www.umsyar.com/515）和笔者之前发表的Carbon, 165, 468 (2020) https://doi.org/10.1016/j.carbon.2020.04.099中关于18碳环电子结构的讨论以了解这方面的背景知识。

1 吸附形成的结构

文中首先研究OPP吸附一个和两个18碳环后的结构。几何优化和振动分析使用Gaussian 16在wB97XD泛函下进行，因为在Carbon, 165, 468 (2020)文中证明了此泛函可以很好地描述18碳环的结构，而且此泛函又可以合理地描述弱相互作用，故拿来优化2C18@OPP复合物很适合。2C18@OPP多达260个原子，而优化+振动分析任务又很昂贵，因此这个环节用的基组不能太贵，文中对包含224个原子的OPP部分使用6-31G*，对更重要但占较小部分的C18使用更大的6-311G*。顺带一提，对C18不能图便宜也用6-31G*，因为如《我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment》（//www.umsyar.com/584）里介绍的论文所示，6-31G*都不足以合理描述C18的几何结构。优化后的带一个和两个18碳环的OPP结构如下所示

从上图可见由于尺寸匹配，18碳环能非常理想地嵌入进OPP的大环中，C18环平面和大环的平面平行。此外，C18的吸附基本不改变大环原有的形状。

2 吸附的热力学

为了从热力学角度探究吸附的可行性，文中对OPP吸附一个18碳环的结合能ΔG进行了计算，并且为了更进一步探究焓和熵所产生的贡献，根据ΔG=ΔH-TΔS，文中将焓变ΔH（结合焓）和熵增对结合自由能的贡献（-TΔS）项分别给了出来。同时为了考察温度对吸附的影响，文中利用《使用Shermo结合 程序方便地计算分子的各种热力学数据》（//www.umsyar.com/552）中介绍的非常方便的Shermo程序一次性得到了各个温度下的热力学量，从而绘制了下图。

从上图可见，结合焓受温度影响相对不大，是个明显的负值，这是C18能够与OPP结合的驱动力。而分子间结合会造成熵的显著降低，因此-TΔS总为正，且温度越高这项越正，越不利于结合，是C18与OPP结合的阻碍。这两部分共同决定结合自由能。在404 K以下，结合自由能为负，C18能够与OPP结合，而在此温度以上则热力学上无法结合；即便之前C18已经吸附了，也会发生脱附。由于这一点，可以通过温度来实现OPP对C18的可控的吸附和释放。即低温下吸附，高温下释放。

以上数据的计算有两点值得说明：

(1)文中使用Shermo基于Grimme提出的准RRHO模型计算热力学校正量中的熵，这一点非常重要，Gaussian用的RRHO模型会严重高估诸如C18@OPP这种含有很低频体系的熵。所以哪怕不靠Shermo考察热力学量随温度的影响，也不要直接用Gaussian给出的熵和自由能热校正量而应当用Shermo来计算它们。

(2)为了准确地计算结合自由能、结合焓，以上数据中涉及的电子能量是使用ORCA在wB97M-V/def2-QZVPP下开着RIJCOSX计算的，这个级别算弱相互作用相当准确，而且是算当前这种较大体系最理想的选择。同时由于wB97M-V是长程100% HF成份的长程校正泛函，自相互作用误差较小，因此可以和具有同等特征的wB97XD一样正确描述18碳环。由于def2-QZVPP很昂贵，故直接算C18@OPP这样两百多原子的体系算不动，因此文中使用了半边模型， 如下所示，是在C18@OPP结构基础上把和吸附无关的那一头的大环去掉并用氢封闭的结构。由于另一头的大环距离吸附的C18甚远，因此这种处理不会影响计算结合能的精度。

3 吸附作用的本质

什么样的机制导致了OPP能够吸附C18？为了考察这一点，本文绘制了如今已被广泛使用的笔者所推广的分子表面静电势穿透图，利用Multiwfn结合VMD可以按照《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图（含视频演示）》（//www.umsyar.com/443）中所示的过程非常容易地绘制，结果如下，是在wB97XD/6-311G*级别绘制的，色彩刻度条单位为kcal/mol。由图可见，18碳环与OPP的范德华表面存在一定程度的穿透，这是两个分子结合后在平衡结构下会出现的典型特征。此图也体现出18碳环的形状着实和OPP的大环吻合得很好，匹配得很完美。由图还可以看出，18碳环表面的静电势非常小，这点在《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）介绍的文章中我也专门指出过。OPP表面静电势虽然比C18更大，但也仅在较小的数值区间内。通过表面静电势的特征就可以看出C18与OPP之间的静电作用是甚微的，不可能对结合起主导作用。

下图是按//www.umsyar.com/443博文说的做法用Multiwfn+VMD绘制的OPP分子表面极值点数值（kcal/mol），可以把OPP表面静电势分布情况展现得更清楚。青色小球是表面静电势极小点，橙色是极大点。

然后我们再关注OPP与18碳环之间可能存在的范德华作用。笔者提出的范德华势是研究这种问题绝佳的手段，介绍见《谈谈范德华势以及在Multiwfn中的计算、分析和绘制》（//www.umsyar.com/551）、使用《Multiwfn对静电势和范德华势做拓扑分析精确得到极小点位置和数值》（//www.umsyar.com/645）。以碳原子作为探针原子（对应18碳环上的原子），OPP的-1.2 kcal/mol的范德华势等值面如下图所示。可见OPP产生的范德华势最负的区域，也即它的范德华作用对碳原子吸引最强的部分，正好与C18原子在吸附后出现的区域相重叠。这充分体现出OPP对18碳环必定有显著的范德华吸引作用，这应当是OPP对18碳环吸附的最主要本质。当前体系是范德华势分析的很好的应用实例。

文中进一步使用Multiwfn结合VMD使用笔者提出的IGMH方法展现了OPP与18碳环的相互作用，对2C18@OPP绘制的0.002 a.u.的分子间的IGMH等值面如下图所示。IGMH的介绍见《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）和《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）。IGMH的等值面非常清晰地把OPP与18碳环之间的相互作用展现了出来，而且其颜色完全为绿色，说明作用区域电子密度很低，而且等值面又很扁且广阔，这明显是色散吸引主导的相互作用的特征。再加上18碳环具有in-plane pi电子，在吸附进OPP后正好与OPP大环的苯环的pi电子分布区域对着，因此可以明确指出18碳环与OPP与形成了显著的pi-pi堆积作用。值得强调的是，靠in-plane电子形成pi-pi堆积，这是碳环体系所独有的特征，非常具有个性。

本文还使用笔者提出的基于力场的能量分解方法EDA-FF考察了各个原子对C18与OPP结合产生的贡献。EDA-FF的介绍见《使用Multiwfn做基于分子力场的能量分解分析》（//www.umsyar.com/442）。此文的EDA-FF完全基于GAFF力场实现。由于C18的原子电荷均为0，所以此分析中没有静电作用项而只有范德华作用部分。EDA-FF计算的平衡结构下C18与OPP的相互作用能为-26.0 kcal/mol，这与高精度wB97M-V/def2-QZVPP计算的-20.3 kcal/mol相差并不很大，因此EDA-FF分析结果完全能说明问题。在上图中，各个原子对C18-OPP结合能的贡献通过原子着色予以了展示。由于C18和OPP上的原子对相互作用贡献的数量级差异较明显，因此用了不同的色彩刻度。由上图可清楚看出，C18上每个原子对相互作用能的贡献在-0.7至-1.0 kcal/mol的范围，离OPP中心（pi-linker部分）越远的原子由于与OPP大环上的原子帖得越近，因而贡献也越大。相应地，OPP大环上的原子与C18离得越近的那些原子对相互作用贡献也越大。而在OPP两个大环之间的pi-linker附近的原子则对结合的贡献基本为0。通过此例可见通过EDA-FF能非常清楚地展现出各个原子对结合所起到的作用，极具实用性。

4 吸附的分子动力学

为了准确地从动态角度考察完整的吸附过程，此文使用GROMACS程序基于经典力场做了吸附过程的分子动力学模拟。在300K下C18进入OPP的总共2000 ps的模拟轨迹的完整视频可在此观看，情况一目了然（其中18碳环的一个原子用红色高亮以便于观察18碳环的旋转情况）：//www.umsyar.com/attach/674/C18_enter_OPP_2ns_300K.mp4

下面说明一下研究C18进入过程的具体细节。

由于OPP是结构很特殊的体系，绝对不能直接用诸如GAFF力场去模拟OPP，否则实测结构会严重变形，两个大环区域会几乎成为严格的圆形而不是液滴状。为了能通过力场合理地描述OPP，就必须利用笔者开发的sobtop程序（//www.umsyar.com/soft/Sobtop）在拓扑文件构建方面提供的灵活性，即允许一部分用刚性参数，一部分用柔性参数。本研究将两个大环之间的相对刚性的pi-linker部分用 计算的Hessian矩阵通过mSeminario方法计算出的力常数和平衡参数来描述，而大环部分则用GAFF力场参数描述，这样可以允许其中的苯环部分发生应有的旋转。另外，模拟C18用的C-C键的参数也都是sobtop基于Hessian矩阵严格算出来的。实测发现，使用这种方式构建的拓扑文件通过GROMACS做能量极小化得到的C18@OPP构型与 优化得到的相当吻合，见下图的对比，蓝色是wB97XD优化出来的，棕色是基于分子力场做能量极小化得到的。

基于以上方式得到的参数，在模拟时，笔者将C18@OPP四周扩了一定区域设立了一个大的周期性盒子，将两个18碳环随机插入在其中，然后在300 K下进行了模拟。下图将碳环进入过程的一些关键时间区间通过多帧叠加方式予以了展现，每个时间区间内碳环按照白到黄的方式着色来清楚区分不同时刻所在的位置，图像通过VMD绘制。下面三幅子图只显示了感兴趣的OPP右半边的情况。第一个C18在模拟才到200 ps时就已经进入OPP了，下图就只关注进入得晚一些的第2个C18的情况。由图可见，852-884 ps时间内C18从OPP外部擦着它的外沿顺利滑入大环空腔，而之前就已经进入另一侧大环中的C18还在原处晃悠没受什么影响。在884-1448 ps期间，刚进入的C18自发贴到大环的内壁倾斜地呆了一阵子。在1448-1460 ps期间，原本贴着大环内壁的亚稳的吸附构型发生了改变，C18自发变得平行于大环。在1460ps到模拟结束（2000 ps）期间内，在大环里躺平的C18就一直稳定地呆在那了，也说明这是在300 K下能够稳定存在的吸附结构。

下图还对C18进入过程做了细致的几何参数分析，包括C18与OPP大环的中心之间的距离，即下图(a)，以及C18与OPP大环之间的夹角，即下图(b)。图中对两个C18的情况都做了统计，分别是红线和蓝线。结合这些几何参数的变化和轨迹，文中将OPP吸附C18的过程分为三个阶段：(1)wander，即C18在OPP周围漫游的过程 (2)adsorption，即C18进入OPP并处于亚稳态构型的过程 (3)equilibrium，即C18在OPP大环中处于完全平衡状态。拿第2个C18为例，下图可见当它处于真空区中游荡而没撞上OPP的时间内，它与OPP大环的距离和夹角都在反复大幅波动，在刚吸附到大环里成为亚稳态时，其夹角和中心间的距离和最终稳定状态都存在一定差异。第1个C18进入过程也是分为这么三个阶段，只不过恰好阶段2非常短暂，即这个C18刚吸附进大环后迅速就躺平达到稳定构型了。这些阶段出现的时刻、长短，都有一定随机性，直接受到初始结构里C18所处的位置、初速度的影响。笔者通过大量重复模拟，都验证了以上描述的C18进入OPP的过程是普遍现象。

下图是模拟过程中C18与OPP的相互作用能的变化，可见吸附时相互作用能有大幅增加（变得更负），并且从C18在大环里倾斜的亚稳态构型变成躺平的稳定构型过程中相互作用能又有所增加。明显体现出这些过程都是能量降低所驱动的。

5 吸附受电子激发的影响

本文研究的OPP的位于中央的pi-linker已知在激发态下具有Baird芳香性，因此电子激发会导致结构发生明显的变形（开口更大、更接近平面），无疑这可能明显影响大环的形状以及对C18的吸附行为。为了探究这一点，本文对OPP的最低单重激发态S1使用TDDFT结合wB97XD泛函做了几何优化，S1结构如下图粉色所示，基态S0结构如下图蓝色所示。可见确实电子激发明显影响了OPP的结构，令pi-linker开口更大了、大环部分变得更圆了。

电子激发计算发现S0-S1跃迁由HOMO到LUMO+8的跃迁所主导，贡献达到86.8%，这两个轨道如下图所示，通过Multiwfn结合VMD按照《用VMD绘制艺术级轨道等值面图的方法（含演示视频）》（//www.umsyar.com/449）介绍的方法绘制。可见都在pi-linker部分，是高度定域化的激发，而大环部分几乎没受影响。这也体现出为什么电子激发能对pi-linker的结构有显著影响。

本文还对2C18@OPP复合物用TDDFT做了几何优化，之后计算了结合能。S0态的结合能为-18.8 kcal/mol，而S1态为-17.1 kcal/mol，说明S1态下OPP对18碳环的吸附变弱了。虽然结合能的改变幅度不算很大，但毕竟吸附的平衡常数受结合自由能影响很大，因此在OPP能够吸附C18的临界温度（前述的404 K）附近，OPP对C18的吸附在一定程度上会受到光控制。

对S1态2C18@OPP复合物，文中也用IGMH方法直观展现了相互作用情况，和前面S0态的图一样也是用的0.002 a.u.绘制。相比之下，可以看出在S1态下等值面往OPP pi-linker部分的延伸程度有所减小，体现出S1态下由于pi-linker开口更大，使得靠近pi-linker的大环原子与C18之间的相互作用明显变弱，这解释了为什么S1态下OPP与C18的结合能的大小变小了。

6 吸附对红外和UV-Vis光谱的影响

文中还考察了OPP吸附18碳环如何影响红外和UV-Vis光谱，如下图所示。可见OPP吸附一个18碳环后会在379.6和2122.0 cm-1处出现新的峰，前者对应碳环的C-C-C键角弯曲振动，后者对应于C-C键伸缩振动。而再吸附第二个C18之后，虽然峰的波数没有进一步变化，但18碳环所引入的峰明显变得更强了。关于碳环的振动特征，强烈建议阅读《揭示各种新奇的碳环体系的振动特征》（//www.umsyar.com/578）里的介绍。

从上图的UV-Vis吸收光谱中可以看到，OPP吸附越多的C18分子，吸收峰越红移，吸收强度越低。

根据以上所示的OPP吸附不同数目的C18所对应的明显不同的红外和UV-Vis光谱，实验化学家可以通过光谱技术检测C18是否吸附到了OPP，以及吸附的量如何。

7 总结

本文介绍了近来笔者在Phys. Chem. Chem. Phys., 25, 16707 (2023)上发表的很有意思的具有8字形双环结构的OPP与新颖独特的18碳环相吸附形成1:2 主-客体复合物的研究。文中利用 计算、波函数分析和分子动力学模拟，从各个角度全面揭示了吸附的本质和特点，还考察了光激发对吸附的影响。本文也是波函数分析与Multiwfn程序研究实际问题的很好范例。相信本文的研究工作的思想和手段能对其他人通过理论化学研究吸附问题产生启发作用。

Introduction to the mol2 format for recording chemical structures

文/Sobereva@北京科音  2022-Dec-31

1 mol2文件简介

mol2文件是计算化学领域非常常用的记录分子结构的格式，被很多程序所支持和利用。例如VMD、GaussView、Chem3D、OpenBabel、AmberTools里的Antechamber等程序都可以导出mol2文件，笔者开发的Multiwfn（//www.umsyar.com/multiwfn）可以基于mol2文件计算EEM原子电荷，笔者开发的Sobtop（//www.umsyar.com/soft/Sobtop）可以基于mol2文件产生GROMACS的拓扑文件，等等。

相对于非常常见的在《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477）中介绍的xyz格式，mol2格式关键优点在于可以记录成键信息，即谁与谁成什么形式的键，这对于判断原子所处的化学环境非常重要，比如Sobtop需要有这样的成键信息才能自动指认GAFF原子类型，Multiwfn计算EEM原子电荷时需要有这样的信息才能判断各个原子要用的EEM参数。另外mol2文件还可以记录其它诸多信息（距离约束、可旋转的键、原子类型和电荷等等），对于分子模拟、QSAR、化学信息学等一些方面有特殊的意义。

mol2格式略微复杂，不同程序产生的mol2文件有所出入，有的程序产生的mol2文件还不规矩，导致经常有人由于用的mol2文件有问题而无法用Sobtop和Multiwfn等程序正常处理，甚至导致程序崩溃。我遂觉得有必要写一篇文章介绍一下mol2格式，便于读者了解怎么读取mol2文件的信息、判断自己手头的mol2文件是否规范，以及拿到不标准的mol2文件时怎么修改。

mol2文件是文本格式，包含大量的字段，每个字段各有各的用处和定义规范。mol2文件的详细说明可以下载此文档查看：//www.umsyar.com/attach/655/Tripos_Mol2_File_Format.pdf。这些字段并不是全都需要出现的，常见的字段只有几个而已，每个字段涉及的信息中往往也只有其中少部分会经常涉及到，将在下文进行介绍，若想更全面详细了解mol2格式可阅读上述文档，共54页。

本文提到的程序的版本是VMD 1.9.3、Multiwfn 3.8(dev, 2022-Dec-18)、Sobtop 1.0(dev3.1)、GaussView 6.0.16、OpenBabel 3.1.1。其它版本可能与本文相同也可能不同。

2 mol2文件的例子和解读

OpenBabel程序产生的mol2格式相对来说是属于比较规矩的，这里结合OpenBabel程序产生的苯甲醛的mol2文件的内容进行讲解。我先用一个可视化程序画了一个苯甲醛的结构，保存为了benzaldehyde.pdb文件，然后用obabel benzaldehyde.pdb -O benzaldehyde.mol2命令用OpenBabel把pdb格式转成了mol2格式的benzaldehyde.mol2文件，其内容如下，GaussView载入此文件后显示的结构图也给出了以便于对照

@<TRIPOS>MOLECULE
benzaldehyde.pdb
 14 14 0 0 0
SMALL
GASTEIGER

@<TRIPOS>ATOM
      1  C          1.9920    0.4700   -0.0000 C.2     0  UNK0        0.1502
      2  C          0.5340    0.2150   -0.0000 C.ar    0  UNK0        0.0142
      3  C         -0.3610    1.2920   -0.0000 C.ar    0  UNK0       -0.0515
      4  C         -1.7360    1.0600    0.0000 C.ar    0  UNK0       -0.0611
      5  C         -2.2160   -0.2520    0.0000 C.ar    0  UNK0       -0.0617
      6  C         -1.3250   -1.3310   -0.0000 C.ar    0  UNK0       -0.0611
      7  C          0.0460   -1.1010   -0.0000 C.ar    0  UNK0       -0.0515
      8  O          2.8450   -0.3960    0.0000 O.2     0  UNK0       -0.2957
      9  H          2.2730    1.5470    0.0000 H       0  UNK0        0.1081
     10  H          0.0250    2.3090   -0.0000 H       0  UNK0        0.0624
     11  H         -2.4300    1.8950    0.0000 H       0  UNK0        0.0618
     12  H         -3.2860   -0.4350    0.0000 H       0  UNK0        0.0618
     13  H         -1.7060   -2.3480   -0.0000 H       0  UNK0        0.0618
     14  H          0.7610   -1.9170   -0.0000 H       0  UNK0        0.0624
@<TRIPOS>BOND
     1     1     2    1
     2     1     8    2
     3     1     9    1
     4     2     3   ar
     5     2     7   ar
     6     3     4   ar
     7     3    10    1
     8     4     5   ar
     9     4    11    1
    10     5     6   ar
    11     5    12    1
    12     6     7   ar
    13     6    13    1
    14     7    14    1

首先要知道mol2文件里以#作为第一列的是注释行，空行也被完全无视。mol2文件是自由格式，因此空格数目完全随意。

第一列以@开头的叫做字段，从上面的benzaldehyde.mol2可见，当前文件有@<TRIPOS>MOLECULE、@<TRIPOS>ATOM、@<TRIPOS>BOND三个字段。一般来说这三个字段都是必须出现的，一起提供了描述一个分子最起码的信息。

2.1 MOLECULE字段

@<TRIPOS>MOLECULE字段记录了体系的基本信息，包括：
第1行：体系的名字。可见OpenBabel把转换出mol2文件的源文件的名字benzaldehyde.pdb当做了当前体系的名字
第2行：五个数字分别是体系中的原子数、化学键数、子结构数、特征数、set数。对于单纯记录体系结构信息，只要提供前两者就够了，后三个可以省略。所谓的子结构是指体系中的一个部分，比如每个分子、每个残基、蛋白质的每条链等等都可以在@<TRIPOS>SUBSTRUCTURE字段里定义为一个子结构。所谓的set是指基于体系中的一些原子/键/子结构根据特定规则和需要定义的集合，可以在@<TRIPOS>SET里具体定义。
第3行：体系的类型。可以为SMALL（小分子）、BIOPOLYMER、PROTEIN、NUCLEIC_ACID、SACCHARIDE
第4行：原子电荷。如果mol2文件没记录原子电荷信息这里就为NO_CHARGES。而在产生当前benzaldehyde.mol2文件的时候OpenBabel自动计算了Gasteiger电荷，因此此处写的是GASTEIGER。还可以为MULLIKEN_CHARGES（Mulliken电荷）、MMFF94_CHARGES（MMFF94力场定义的电荷）等等，不同种类电荷都有固定名字。如果记录的原子电荷是比如Multiwfn算的ADCH、RESP、1.2*CM5等电荷，在mol2格式规范中没有对应的名字，则这里应当写USER_CHARGES。

2.2 ATOM字段

@<TRIPOS>ATOM字段每一行定义一个原子的信息，每一列记录的信息为：
(1)原子序号（整数）
(2)原子名（字符串）
(3)X坐标（埃）
(4)Y坐标（埃）
(5)Z坐标（埃）
(6)原子类型（atom type。字符串）
(7)原子所属的子结构序号（整数），可省略
(8)原子所述的子结构名字（字符串），可省略
(9)原子电荷（浮点数），可省略

原子名部分可以为比如C2、H4等等，完全随意。记录生物分子结构时通常用IUPAC定义的各种残基中的原子名。

原子类型部分可以记录做分子模拟用的力场中此原子实际对应的原子类型。mol2格式自己也有一套原子类型定义，见前述的Tripos_Mol2_File_Format.pdf文档的末尾，比如sp3杂化的碳的原子类型是C.3，C.ar代表芳香碳，Any代表任意，Hal泛指卤素，Cl代表氯，Ca代表钙，H代表氢，H.spc特指SPC水模型的氢，LP代表孤对电子（lone pair），Du代表虚原子（dummy），等等。

一定要特别注意，mol2格式虽然定义了一大堆字段，但（居然）没有一个地方是专门用来记录元素的，这在我来看是mol2格式的严重不足！！！mol2记录的原子名和原子类型信息可以与元素名相同也可以不同，不同程序产生的mol2文件的情况各有不同。例如如此例可见，OpenBabel产生的这个mol2文件里原子名恰等于元素名，原子类型是根据mol2格式自己定义的原子类型指认的。而在GaussView产生的mol2文件中，原子名是给元素名后面加了数字（因此不会有重名的原子），而原子类型恰等于元素名。由于情况混乱，所以一个程序在读取mol2文件的时候并没有严格的办法能准确判断元素，只能靠猜。Multiwfn和Sobtop在读取mol2文件时是根据原子类型的字符串判断元素的：如果字符串中没有.，就直接将之当做元素名来判断元素；如果有.，比如是C.3，就把.前面的内容当做元素名判断元素。因此，读者应该知道Multiwfn和Sobtop没判断对元素时该怎么办了，最简单的做法就是手动修改mol2文件以让@<TRIPOS>ATOM字段每一行的第6列对应元素名。

由于子结构信息在原子电荷前面，因此即便你不想定义原子所属的子结构信息而只想定义原子电荷，也必须随便写上子结构序号和子结构名字来占位，比如此例用0  UNK0来占位。

2.3 BOND字段

@<TRIPOS>BOND字段每一行定义一个键的信息，其每一列记录的信息为：
(1)键的序号（整数）
(2)第1个原子的序号
(3)第2个原子的序号
(4)键的类型

键的类型有以下这些
• 1 = 单键
• 2 = 双键
• 3 = 三键
• am = 酰胺的N-C键（这种键有一定pi共轭作用，这是为什么mol2格式里特意用am来与单键区分）
• ar = 芳香环（aromatic）上的键，以下简称芳香键
• du = 虚键
• un = 未知/无法判断
• nc = 不相连（俩原子不成键就没必要在BOND字段出现，但可以靠nc强调某两个原子间就是没成键）
绝大多数程序产生的mol2文件里没有du、un、nc。

有的程序产生的键的类型名不规矩。比如GaussView对于芳香环上的键（具体来说，是图形窗口里看到一个实线+一个虚线的那种键）用的类型是Ar，但按照mol2规范应当是ar，因此这会导致一些程序无法识别（而Sobtop在读取时已经考虑到了GaussView的这个bug，因此用户不用自己做替换）。GaussView还自行给mol2格式做了扩展，把图形窗口里看到是一个虚键的那种键记录为Wk代表Weak，而这可能导致很多程序无法正常识别和载入。

不同程序对成键的指认也往往有很大不同。比如甲酰胺，Avogadro产生的它的mol2里把N-C键记录为am，严格符合mol2格式的要求。而GaussView则无法将之记录为am，而是可能记录为单键也可能记录为Ar（取决于当前图形窗口里显示的成键方式）。另外，我之前在《谈谈原子间是否成键的判断问题》（//www.umsyar.com/414）中说过GaussView是根据原子半径和原子间距离判断成键形式的，导致很有可能判断出的成键方式很不“经典”，甚至很违背化学常识，比如可能显示一个碳原子连着一个双键和两个单+虚键。而如果把结构保存成pdb然后再用OpenBabel转成mol2格式，则成键方式就很经典了，因为OpenBabel能够自动让成键方式满足经典Lewis式且同时识别芳香区域。

VMD程序里如果载入了xyz、gro、pdb等不含成键关系信息的文件（虽然pdb有CONECT字段，但VMD不利用），保存出的mol2文件将没有BOND字段，明显是不符合规范的。在同时载入拓扑文件以提供拓扑信息后，保存出的mol2文件才是有BOND字段的（与此同时也有了原子类型、原子电荷、原子所属的残基号和残基名信息）。VMD如果载入的是mol2文件，也可以从中获取成键信息，使得保存出的mol2文件里也有BOND部分。但是VMD并不能记录芳香键，而且它自己也没有像OpenBabel那样根据几何结构和元素就能判断出芳香区域的能力，因此即便载入的是本例的benzaldehyde.mol2，保存出的mol2里苯环上的键也都会简单地记录成单键。

还值得顺带一提的是有个叫mol的格式，介绍见https://en.wikipedia.org/wiki/Chemical_table_file。不要将它和mol2混淆，二者格式截然不同。mol也能像mol2一样记录键的存在性及其类型。GaussView产生的mol文件中会把芳香键用4来记录，这正是mol标准格式中的芳香键的类型序号，而OpenBabel在产生mol文件时则会把芳香环描述成单双键交替的形式。

3 关于记录晶胞信息

mol2文件通常用来记录孤立体系，但实际上此格式也定义了记录晶胞信息的字段@<TRIPOS>CRYSIN，在其下一行写晶胞的a、b、c三个边长（埃）以及alpha、beta、gamma夹角（度），每个值之间以逗号分隔。例如：
@<TRIPOS>CRYSIN
3.785,3.785,9.514,90,90,90
Sobtop和Multiwfn在载入mol2文件时都会试图载入晶胞信息以考虑周期性。GaussView有很好用的对周期性体系建模的功能，但即便在图形窗口里能看到晶胞边框，代表当前是周期性体系，其保存出的mol2文件里也没有@<TRIPOS>CRYSIN字段，因此需要自行补充。在VMD 1.9.3里，即便当前有晶胞信息，保存的mol2文件里也没有以上字段。VMD 1.9.3载入mol2时也不会从以上字段中载入晶胞信息。

A lazy script for ORCA combined with Multiwfn to calculate RESP, RESP2 and 1.2*CM5 atomic charges

文/Sobereva@北京科音

First release: 2022-Mar-15  Last update: 2022-Aug-6

之前笔者在以下文章中提供了三个Linux shell脚本，分别用来自动调用机子里的Gaussian和Multiwfn程序实现一键计算1.2*CM5、RESP和RESP2原子电荷，它们对于做经典力场的分子动力学非常重要。
计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本
//www.umsyar.com/585（http://bbs.keinsci.com/thread-21462-1-1.html）
计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）
//www.umsyar.com/476（http://bbs.keinsci.com/thread-12858-1-1.html）
RESP2原子电荷的思想以及在Multiwfn中的计算
//www.umsyar.com/531（http://bbs.keinsci.com/thread-16190-1-1.html）

如今用免费的ORCA 程序的人也很多。为了便于那些主要做分子动力学模拟、没买Gaussian又不太懂ORCA程序使用的人也能便利地计算上述原子电荷，笔者写了能自动调用ORCA和Multiwfn的计算那些电荷的Linux shell脚本。这些脚本提供在了Multiwfn程序文件包里（可以在Multiwfn主页//www.umsyar.com/multiwfn免费下载），必须是2022年3月8日及以后更新的Multiwfn版本里才有，包括：
examples\RESP\RESP_ORCA.sh：计算RESP电荷的脚本
examples\RESP\RESP2_ORCA.sh：计算RESP2电荷的脚本
examples\scripts\1.2CM5_ORCA.sh：计算1.2*CM5电荷的脚本

这些脚本的用法和前述帖子里介绍的基于Gaussian的脚本精确一致，需要留意的地方也都相同，只不过脚本中调用Gaussian的地方变成了调用ORCA而已，故不再累述用法。如果不会装ORCA的话看《 程序ORCA的安装方法》（//www.umsyar.com/451）。

这些脚本运行之前记得用文本编辑器打开，把ORCA=和orca_2mkl=后面的内容分别改为当前机子里实际的ORCA和orca_2mkl工具的路径。并把nprocs=后面的值改为计算时要调用的CPU核心数。脚本里maxcore=后面的值是ORCA的每个并行进程的内存使用量（MB）上限，与nprocs的乘积必须明显小于空余物理内存量，其默认值一般是合适的。如果空余物理内存不够则需要适度减小maxcore或nprocs，而如果要算几百个原子的大体系则需要适度加大maxcore，否则可能计算崩溃。

RESP_ORCA.sh和RESP2_ORCA.sh里默认用的优化级别和基于Gaussian的脚本（RESP.sh和RESP2.sh）有所不同，这里用的是ORCA才支持的B97-3c，因为这个级别做优化很快，结果准确度也不错。由于这个差异，以及ORCA和Gaussian在溶剂模型的实现上有所差异，所以基于Gaussian和基于ORCA的脚本得到的RESP或RESP2电荷可能有零点零几的差别，这点没必要在意，都是合理的。

关于使用脚本时哪些溶剂可以直接用、溶剂名怎么写，请在ORCA手册里搜“solvents in the SMD library”查看内置的溶剂名列表。注意ORCA里有些溶剂名是带空格的，对这种情况要把溶剂名用双引号扩住，例如./RESP2_ORCA.sh HF.pdb 0 1 "ETHYL ETHANOATE"。

如果你的输入的结构文件里的结构就已经足够好，不想让脚本自动再做优化浪费时间，可以用examples\RESP\目录下的RESP_ORCA_noopt.sh和RESP2_ORCA_noopt.sh分别代替前述的RESP_ORCA.sh和RESP2_ORCA.sh，它们的用法完全一样，只不过带_noopt后缀的不做优化步骤。

使用这些脚本计算原子电荷发表文章的话请在文中恰当引用Multiwfn和ORCA。

Using Sobtop to create GROMACS topology file for ferrocene super conveniently

文/Sobereva@北京科音   2022-Feb-17

0 前言

Sobtop是笔者开发的主要用于产生各种类型体系的GROMACS拓扑文件的程序，可以在主页//www.umsyar.com/soft/Sobtop免费下载。过渡金属配合物体系的拓扑文件在以往是很难产生的一类情况，因为里面牵扯到GAFF、OPLS-AA、GROMOS等主流有机分子力场不支持的元素，里面涉及到的键/键角/二面角参数在那些力场里也没有，所以acpype、Ligpargen、ATB等程序根本没法搞。而这类体系，用Sobtop则可以超级容易地产生拓扑文件。AmberTools里的MCPB.py程序虽然能产生这类体系的拓扑文件，但是过程麻烦至极，在http://bbs.keinsci.com/thread-27796-1-1.html里有人演示了怎么借助MCPB.py产生二茂铁的GROMACS的拓扑文件，估计大多数人看了后会打退堂鼓。Sobtop的详情以后会另文专门系统性地介绍，本文仅简单演示一下怎么用Sobtop快捷地产生二茂铁的拓扑文件，通过对比大家可以充分认识到用Sobtop有多么方便。

本文牵扯到的所有输入输出文件都可以在//www.umsyar.com/attach/635/file.rar里找到。Sobtop用的是2022-Feb-16发布的预览版1.0(dev)。Multiwfn用的是2022-Feb-17更新的版本，Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载。

1 用Gaussian做优化和振动分析

首先用Gaussian对二茂铁做优化和振动分析。要注意，二茂铁的两个茂环间仅在低温晶体中是交错式的，而在真空中极小点结构是重叠式的（参看Materials, 8, 7723 (2015)和https://en.wikipedia.org/wiki/Ferrocene#Determining_the_structure中的相关信息），因此本例用重叠式（D5h点群）的初始结构。记得Gaussian计算前最好用GaussView做对称化，这样优化得又精准速度又快。

关于计算级别，本例用的TPSSh/def2-TZVP级别对于过渡金属配合物体系通常是很理想的选择。关于泛函的选择看《简谈 计算中DFT泛函的选择》（//www.umsyar.com/272）。如果想明显更便宜，Fe可以用SDD，配体可以用6-311G*。

此例的Gaussian输入文件如下

%chk=Ferrocene.chk
#p tpssh/def2tzvp opt freq
[空行]
Generated by Multiwfn
[空行]
0 1
 C                  0.71056100   -0.97800331    1.65426528
 C                 -0.00000000    1.20887857    1.65426528
 C                 -0.71056100   -0.97800331    1.65426528
 H                  2.17848562    0.70783289    1.62195576
 H                 -2.17848562    0.70783289    1.62195576
 Fe                 0.00000000    0.00000000   -0.00000000
 C                  1.14971184    0.37356402    1.65426528
 C                 -1.14971184    0.37356402    1.65426528
 H                  1.34637816   -1.85313055    1.62195576
 H                 -0.00000000    2.29059533    1.62195576
 H                 -1.34637816   -1.85313055    1.62195576
 C                  1.14971184    0.37356402   -1.65426528
 C                 -1.14971184    0.37356402   -1.65426528
 C                  0.71056100   -0.97800331   -1.65426528
 C                  0.00000000    1.20887857   -1.65426528
 C                 -0.71056100   -0.97800331   -1.65426528
 H                  2.17848562    0.70783289   -1.62195576
 H                 -2.17848562    0.70783289   -1.62195576
 H                  1.34637816   -1.85313055   -1.62195576
 H                  0.00000000    2.29059533   -1.62195576
 H                 -1.34637816   -1.85313055   -1.62195576
[空行]
[空行]

算完了之后，用formchk把当前目录下的Ferrocene.chk转换为Ferrocene.fchk。不知道formchk是什么和怎么用的话看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）里相关部分。

还需要得到二茂铁的mol2文件，并且要求其中记录的键连关系和期望的一致，这决定了Sobtop给出的拓扑文件里都有哪些力场项。很多程序都可以实现这个，此例我们用常用的GaussView打开Ferrocene.fchk，如果看到Fe和各个C之间没显示键的话就手动连上（设几重键都无所谓，Sobtop只看有没有键）。然后保存为Ferrocene.mol2。

2 计算RESP电荷

RESP电荷很适合用于动力学模拟的目的，我们要让拓扑文件里的原子电荷对应RESP电荷。计算RESP电荷非常强大、方便、灵活的程序是Multiwfn，见《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》（//www.umsyar.com/441）。要注意对二茂铁体系不宜直接按照常规方式算RESP电荷，因为拟合点的分布不会恰好满足D5h对称性，因此得到的等价原子的原子电荷会有微小的偏差。为了让所有等价的原子的电荷都严格相同，在Multiwfn中计算RESP电荷时应当设置等价性。

启动Multiwfn，然后输入
Ferrocene.fchk
7  //计算原子电荷
18  //计算RESP电荷
5  //等价性约束设置
11  //根据点群自动判断等价性
a  //根据整个体系的结构判断点群
y  //当前屏幕上显示的点群确实是D5h，列出的原子等价性没错，所以这里输入y确认。然后当前目录下有了等价性设置文件eqvcons_PG.txt
q  //返回
1  //从文本文件里读取等价性设置
eqvcons_PG.txt  //刚才产生的文件
2  //开始一步式RESP电荷拟合
[回车]  //程序提示Fe缺半径，按回车表示用推荐的半径
Multiwfn算RESP电荷效率很高，RESP电荷很快就产生了，从屏幕上给出的结果看确实原子电荷合理，而且满足等价性。按y把产生的原子电荷导出为当前目录下的Ferrocene.chg。

3 用Sobtop产生拓扑文件

启动Sobtop，然后依次输入
Ferrocene.mol2
7  //指定原子电荷
10  //从Multiwfn的chg文件中读取原子电荷
Ferrocene.chg
0  //返回
2  //产生GROMACS的.gro文件
[回车]  //此时当前目录下产生了Ferrocene.gro
1  //产生GROMACS的拓扑文件
2  //优先用GAFF原子类型，缺的自动用UFF力场的。原子类型决定了拓扑文件里的LJ参数
2  //所有成键相关参数基于振动分析产生的fchk文件里的Hessian矩阵生成
Ferrocene.fchk  //如果此文件和Ferrocene.mol2在相同目录下，此处直接按回车可以快捷载入
[回车]  //在当前目录下产生Ferrocene.top
[回车]  //在当前目录下产生Ferrocene.itp

现在拓扑文件在当前目录下出现了，Sobtop可以关了。现在可以打开itp文件看看内容，会看到没有任何问题，文件整齐干净，注释也特别清楚。

4 测试拓扑文件合理性

产生拓扑文件之后，我一般建议在真空下对这个体系跑一下动力学，根据模拟过程中结构变化情况大致判断拓扑文件是否合理。做这个模拟要用的Ferrocene.gro、Ferrocene.top、Ferrocene.itp前面我们都已经产生了。做这个模拟用的md.mdp在本文的文件包里也给了，可见是控温在298.15K下模拟100 ps，1 fs一步，每100 fs往xtc里写入一帧，没有用PBC。笔者此例用的是GROMACS 2018.8。这个mdp不兼容GROMACS较新版本，较新版本用户需要在一个足够大的空盒子里在考虑PBC的情况下跑NVT模拟。

把前述的gro、itp、top、mdp都放到当前目录下，运行以下命令
gmx grompp -f md.mdp -c Ferrocene.gro -p Ferrocene.top -o md.tpr
gmx mdrun -v -deffnm md

然后用VMD载入轨迹，用Extensions - Analysis - RMSD Trajectory Tool工具对system选区做align消除平动转动后，每10帧叠加显示一次，得到下图

可见模拟过程中二茂铁很好地维持了刚性，说明拓扑文件合理。

5 总结&其它

本例充分体现了使用Sobtop结合Multiwfn程序产生过渡金属配合物体系拓扑文件的便利，把复杂度降低到了极限。此文只以极为简单的体系做了示例，对更复杂的情况，诸如多核金属配合物、除了刚性部分还带有可旋转键的柔性部分的配合物，Sobtop也都能很容易地产生拓扑文件。凡是体系全为刚性的情况都可以直接效仿本文的做法；对于既有刚性配位区域也有柔性基团的配合物，在Sobtop问你怎么产生成成键相关(bonded)参数时，建议选择相应选项，让bond和angle参数用特定方法基于Hessian产生，而其它的用预置力场参数，这样可以使得GAFF的参数用于描述柔性部分二面角可旋转特征。

实际中二茂铁的茂环的旋转势垒是非常低的，常温下是可以旋转的，但本文得到拓扑文件对应的是纯刚性二茂铁。本文暂不考虑茂环间可相对旋转的问题，这也不是什么重要问题。

本文的做法比起MCPB.py实在方便百倍，而且MCPB.py用的Seminario方法产生力常数也不如本文的方法合理（当前版本Sobtop默认用的是modified Seminario，产生的键角力常数远好于Seminario）。

做优化和振动分析用免费的ORCA亦可，文中用fchk的地方改成用ORCA振动分析产生的.hess文件即可。

对于体系没有对称性的情况，算RESP电荷时不需要本文这样做对称等价设置，直接在Multiwfn的RESP界面里选择1用常规的两步式拟合即可。更多细节看《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》（//www.umsyar.com/441）。

本文的例子直接从chg文件中载入了RESP电荷。也可以在Sobtop里先不设置电荷，等itp产生出来之后再手动把要用的电荷写入到[ atoms ]部分。

使用Sobtop产生拓扑文件发表文章时记得按照Sobtop主页//www.umsyar.com/soft/Sobtop的说明进行引用。引用方式在未来会有所变化，请记得发文章之前看一下主页上最新的引用说明。