//www.umsyar.com/category/第一性原理 //www.umsyar.com/719 2024-07-12T14:52:00+08:00 使用Multiwfn对周期性体系做键级分析和NAdO分析考察成键特征 Using Multiwfn to perform bond order analysis and NAdO analysis to study bonding characters for periodic systems 文/Sobereva@北京科音  2024-Jul-12 0 前言 键级是考察化学键特征的一类重要方法，非常强大的波函数分析程序Multiwfn可以方便地计算很多类型的键级，在《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）的键级部分有具体的说明，Multiwfn手册4.9节有一些实际计算例子，在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里做了特别全面详细的讲解和演示。另外，Multiwfn还支持NAdO方法，可以将原子间的模糊键级以轨道形式图形化展现，对于了解成键本质极其有用，介绍和实例看《使用键级密度(BOD)和自然适应性轨道(NAdO)图形化研究化学键》（//www.umsyar.com/535）。 以上博文都是举的孤立体系的例子，而如今Multiwfn已经将上述方法扩展到了对周期性体系的分析上。将Multiwfn与CP2K做周期性DFT计算产生的波函数相结合，可以很容易地用上述方法考察周期性体系的成键情况，在本文将通过一个层状的共价有机框架化合物（COF）体系，以及一个Pd表面吸附苯的体系，对具体做法进行演示，希望读者举一反三将本文介绍的方法运用到实际研究中。使用本文的方法计算的结果若用于发表，记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。 读者务必使用2024-Jun-24及以后更新的Multiwfn版本，否则情况与本文所述不符。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载。不了解Multiwfn者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。本文要用到非常流行、高效且免费的第一性原理程序CP2K，笔者假定读者已具备了相关常识，推荐不熟悉者通过“北京科音CP2K第一性原理计算培训班”（http://www.keinsci.com/workshop/KFP_content.html）系统性学习。本文要使用Multiwfn创建CP2K输入文件，这在《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里有简要介绍。本文用的CP2K是2024.1版。 下面的例子涉及到的所有文件都可以在//www.umsyar.com/attach/719/file.rar里得到。 1 单层COF体系 这部分示例的COF体系的cif文件是本文文件包里的COF_16371N2.cif，用GaussView显示的结构如下，可见晶胞里有两层。 1.1 准备输入文件 此体系两层之间的相互作用主要是pi-pi堆积，对电子结构影响微乎其微，因此只计算一层就够了，可以节约时间，尤其是能显著节约NAdO分析过程中计算原子重叠矩阵（AOM）的时间。于是在GaussView里删除一层，然后保存为COF.gjf。由于当前体系不是导体且晶胞边长约15埃已经不小了，所以可以不用扩胞就直接做gamma点计算。注意CP2K导出molden文件只支持gamma点的情况，这在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里已明确说了。 启动Multiwfn，然后输入 COF.gjf cp2k  //进入产生CP2K输入文件的界面 [回车]  //输入文件名用默认的COF.inp -2  //要求产生molden文件 -7  //设置周期性 XY  //仅在平行于COF的方向考虑周期性，而垂直于COF的方向用非周期性 -9  //其它设置 16  //要求将Kohn-Sham矩阵导出到.csr文件里。此文件在Multiwfn计算NAdO轨道能量时会用到 0  //返回 0  //产生输入文件 现在当前目录下就有了COF.inp，对应PBE/DZVP-MOLOPT-SR-GTH级别的单点计算。用CP2K运行之，当前目录下会产生记录波函数信息的COF-MOS-1_0.molden和记录KS矩阵的COF-KS_SPIN_1-1_0.csr。按照《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》所述，在molden文件开头插入以下内容定义晶胞信息和当前赝势基组描述的各元素的价电子数  [Cell]  14.94670000     0.00000000     0.00000000  -7.47335000    12.94422190     0.00000000   0.00000000     0.00000000     8.00000000  [Nval]  C 4  N 5  H 1 之后就可以基于这个molden格式的波函数文件做各种周期性体系的分析了。 1.2 计算Mayer键级 首先计算非常常用的Mayer键级，它衡量原子间等效的共享电子的对数。启动Multiwfn，载入COF-MOS-1_0.molden，然后输入 9  //键级分析 1  //Mayer键级 Multiwfn首先计算重叠矩阵，然后立刻输出了Mayer键级，只有大于阈值的（由settings.ini里的bndordthres参数控制）才直接显示了出来。要想看所有原子间的Mayer键级可以再选择y导出键级矩阵然后查看相应的矩阵元。为了便于对照，下图显示了当前体系各个原子的序号，本文主要关注其中绿色虚线圈住的那些原子。 Multiwfn输出的一些有代表性的键的Mayer键级值如下 11(C )   19(C )    1.12454138 11(C )   22(N )    1.48268711 19(C )   48(C )    1.27625346 20(C )   50(C )    1.39369508 18(C )   48(C )    1.24722702 可见这些键的Mayer键级都不同程度地明显大于1.0，因此可以推测它们都不仅仅形成了sigma键，还有一定pi成份。在后文会计算pi键级和使用NAdO分析更进一步确认这一点。上面列出的C-C键键级可以体现出体系中哪种C-C键的强度更强。可见连接萘单元和C3N3单元的C11-C19的键是相对最弱的。 1.3 计算模糊键级（fuzzy bond order） Mayer键级怕弥散函数，而模糊键级则没有这个问题，用于任意基组都可以。Mayer键级用于MOLOPT系列基组都是没问题的，不过这里作为例子，也计算一下模糊键级。由于计算模糊键级需要先构造AOM，对较大的体系很耗时，所以对于Mayer键级适用的情况建议优先用Mayer键级。注意对孤立体系Multiwfn计算模糊键级默认是基于Becke划分的，而对于周期性体系计算模糊键级默认是基于Hirshfeld划分的，因此默认设置下周期性体系的计算结果和孤立体系的没有严格的可比性。 在键级分析主功能的菜单里选择7，然后Multiwfn就会开始计算AOM，之后立刻给出模糊键级的计算结果，前述的那些键的模糊键级如下 11(C )   19(C )    1.05202861 11(C )   22(N )    1.50866214 19(C )   48(C )    1.22313834 20(C )   50(C )    1.38494183 18(C )   48(C )    1.17984661 虽然模糊键级与Mayer键级在定量上有一些差别，但不同的键的键级的大小顺序是完全一致的，没有结论上的差异。 1.4 计算pi电子贡献的Mayer键级 我在《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）中专门介绍了pi键级的计算方法，也即pi占据轨道对Mayer键级的贡献，没读过的话务必先阅读。这一节对单层COF这个周期性体系也做这个计算，看看各个键的pi作用程度的差异。由于当前体系是精确平行于XY平面的，所以让Multiwfn自动指认pi轨道前不需要先做轨道定域化。 启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入 100  //其它功能（Part 1） 22  //检测pi轨道 0  //当前的轨道都是离域形式（如正则分子轨道） 2  //设置其它轨道占据数为0 0  //返回 9  //键级分析 1  //Mayer键级 前面提到的那些键的键级计算结果如下，这对应的就是pi键级。可见C3N3六元环中的C-N键的pi共享电子作用很强，比这里列出的C-C键的还要更显著。 11(C )   19(C )    0.14583409 11(C )   22(N )    0.41775145 19(C )   48(C )    0.28202057 20(C )   50(C )    0.38270736 18(C )   48(C )    0.27414145 使用此做法发表文章时，除了引用Multiwfn原文外还建议同时引用介绍Multiwfn做pi电子结构分析的文章Theor. Chem. Acc., 139, 25 (2020)。 1.5 NAdO分析 利用前述的《使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键》中介绍的NAdO分析，可以将sigma和pi作用分别以轨道形式呈现出来，非常有价值。例如在我的研究18氮环的ChemPhysChem (2024) https://doi.org/10.1002/cphc.202400377文中就用了NAdO分析非常清楚、严格地展示了较短的N-N键的sigma+pi作用特征。本节我们用NAdO分析展示一下C11-C19的sigma和pi轨道相互作用情况。 启动Multiwfn并载入COF-MOS-1_0.molden后，依次输入 15  //模糊空间分析 3  //计算AOM。对周期性体系默认是基于Hirshfeld原子空间计算的 [回车]  //用默认的0.2 Bohr格点间距。此设置对当前体系耗时不高，但对于明显更大体系为了节约耗时可以适当用更大的格点间距如0.35 Bohr，但会多多少少牺牲精度 计算完成后当前目录下出现了记录AOM的AOM.txt。在Multiwfn里接着输入 0  //返回 200  //其它功能（Part 2） 20  //BOD和NAdO分析 -1  //要求计算NAdO的能量 2  //从外部文件读取Fock矩阵 COF-KS_SPIN_1-1_0.csr  //输入此文件的实际路径 1  //基于AOM计算原子间的相互作用 [回车]  //读取当前目录下的AOM.txt 11,19  //分析11和19号原子间的相互作用 马上NAdO分析就做完了，屏幕上显示以下内容。这说明所有NAdO的本征值加和为1.052，和上一节得到的模糊键级一致。并且当前只有前两个NAdO的本征值明显大于0，它们是最值得考察的，其它的基本可以忽略  Eigenvalues of NAdOs: (sum=   1.05204 )    0.78955   0.19252   0.06784   0.02133   0.00807   0.00651   0.00127    0.00096   0.00031   0.00011   0.00010   0.00005   0.00004   0.00000 ...略 所有的NAdO轨道现在都已经导出到了当前目录下的NAdOs.mwfn中。现在输入y让Multiwfn载入此文件，之后退回到Multiwfn主菜单，进入主功能0查看NAdO轨道。不熟悉Multiwfn看轨道功能的话参考《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）。在图形界面里分别选择1号和2号轨道，可分别看到下面的图像，同时在文本窗口显示的NAdO本征值和轨道能量也标上了。这里等值面数值用的是0.08。 由上可见，C11-C19同时具有sigma和pi作用特征，前者远比后者显著得多。而且由于C19处在pi共轭区域，所以NAdO方法产生的主要描述C11-C19的pi作用的轨道并不完全定域在C11和C19上，也同时一定程度离域到了周围原子上。通过NAdO轨道能量还可以看到，C11和C19之间共享的sigma电子的能量明显比pi电子要低，这十分符合一般化学常识。 2 Pd(100)表面吸附苯 本节通过Pd(100)晶面吸附苯的体系，演示一下苯与Pd(100)基底这两个片段间Mayer键级的计算，以及对于像这样很大的体系怎么尽可能节约整个NAdO分析的时间。此体系的CP2K做结构优化的输入和输出文件，以及优化任务跑完后产生的restart文件都在本文文件包里提供了。优化完的结构如下，可见吸附作用很强，以至于Pd-C的作用都令苯产生弯曲了。 由于此例子的molden文件、FOM.txt和NAdOs.mwfn文件太大，本文的文件包里就没提供了。 2.1 准备输入文件 用Multiwfn载入本文文件包里的Pd+ben_opt-1.restart读取优化完的坐标和晶胞信息，然后输入以下命令创建一个做单点且产生molden文件的任务 cp2k  //产生CP2K输入文件 Pd+ben_SP.inp  //产生的输入文件名 -7  //设置周期性 XY -2  //要求产生molden文件 6  //由于Pd基底是导体，开启smearing 0  //产生输入文件 为了帮助SCF收敛，把Pd+ben_SP.inp里的ALPHA设为0.15，NBROYDEN设为12。用CP2K计算Pd+ben_SP.inp，得到Pd+ben_SP-MOS-1_0.molden。然后在里面开头部分加入以下内容。  [Cell]  23.33880000     0.00000000     0.00000000  0.00000000    23.33880000     0.00000000  0.00000000     0.00000000    24.79750184  [Nval]  Pd 18  C 4  H 1 2.2 片段间Mayer键级 这一节将苯和Pd(100)基底分别定义为一个片段计算它们之间的Mayer键级。启动Multiwfn，载入Pd+ben_SP-MOS-1_0.molden，然后输入 9  //键级分析 -1  //定义片段 289-300  //苯的部分 1-288  //Pd(100)基底部分 1  //计算键级 Multiwfn开始产生重叠矩阵，之后给出了原子间键级，以下是其中一部分，可见苯中的C和有的Pd之间的Mayer键级较大，即共价作用很显著。 ...略  # 1412:       246(Pd)  270(Pd)    0.21887832  # 1413:       246(Pd)  271(Pd)    0.21771061  # 1414:       246(Pd)  289(C )    0.08497491  # 1415:       246(Pd)  293(C )    0.09391700  # 1416:       246(Pd)  294(C )    0.55236636  # 1417:       247(Pd)  271(Pd)    0.23611928  # 1418:       247(Pd)  272(Pd)    0.23841126 ...略 老有人问怎么判断固体表面对分子是物理吸附还是化学吸附，判断方法很多，如结合能的数量级、原子间距离和原子半径的关系、电子密度差、IGMH（//www.umsyar.com/621）等等，而用Mayer或模糊键级是最能说明问题的。以上Mayer键级体现的C和Pd之间的明显的共价作用意味着二者之间形成了明显的化学键作用，显然Pd(100)对苯是化学吸附。如果键级只有比如零点零几的程度，那一般可以说是物理吸附，但也不排除形成典型离子键的可能，这通过原子电荷可以判断，计算方式可参考《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）。 最后Multiwfn还给出了片段间键级： The bond order between fragment 1 and 2:    4.309784 即苯分子和Pd(100)之间等效共享电子对数多达4.3，显然太共价了，如果这都不叫化学吸附... 2.3 NAdO分析 这一节我们将Pd(100)与苯之间显著的共价作用通过NAdO方法以轨道形式直观展现。本节的做法和1.5节示例的常规的NAdO分析流程有两个明显的不同： (1)对当前这样原子数又多、轨道数又多的大体系计算所有原子的AOM是相当耗时间的事情，而且导出、载入AOM.txt的耗时巨高而且巨占硬盘（正比于占据轨道数的平方乘以原子数）。为了极大程度节约时间，可以使用Multiwfn提供的一个特殊的策略，也就是在主功能15里用选项33计算感兴趣的两个片段间的片段重叠矩阵（FOM）并导出为FOM.txt，然后对两个片段之间做NAdO分析时直接从中读取要用的FOM，这样不仅省得导出和载入巨大的AOM.txt，而且产生FOM的过程中不需要计算片段里没有涉及的原子的AOM，当片段里的原子数占整个体系的原子数比例较少的时候这可以大幅节约时间。对于当前的例子，苯可以作为片段1，与苯相互作用最密切的那部分Pd原子适合作为片段2，这比起把所有Pd都定义为片段2能节约巨量耗时。但用户很难判断哪些Pd原子与苯作用密切，此时可以利用Multiwfn特意提供的一种定义片段2的方法，即片段1以外的原子中，与片段1任何一个原子之间Mayer键级大于指定阈值的原子都被定义为片段2。利用这个做法可以使得与苯的轨道相互作用最显著的Pd原子被定义为片段2，其数目远少于总的Pd数目。 (2)在计算FOM和做NAdO分析之前，需要先用主功能6的子功能38令电子根据轨道能量由低到高重新排列，使得轨道占据数都为整数。这是因为当前CP2K计算用了smearing，导致费米能级附近的轨道处于分数占据状态，这不是单行列式波函数应具有的特征，而Multiwfn的NAdO分析目前只支持单行列式波函数。重新排列后，当前的波函数就成了标准的单行列式波函数形式，从而可以做NAdO分析。另一方面，这么重排占据数后在计算AOM和FOM时可以只对占据轨道之间计算，也只有这部分是NAdO分析所实际需要的，只计算占据轨道之间的AOM/FOM的耗时、对内存的占用远远低于考虑所有轨道的情况， 启动Multiwfn，载入Pd+ben_SP-MOS-1_0.molden，然后输入 6  //修改波函数 38  //按照轨道能量由低到高重排占据数 -1  //返回 15  //模糊空间分析 33  //计算FOM 3  //定义两个片段并计算FOM，不属于片段1的原子若与片段1的任意原子的Mayer键级大于指定值则被定义为片段2 289-300  //苯的部分 0.001  //Mayer键级的阈值 之后看到以下信息，告诉了你哪些原子被作为了片段2，并且可见片段1和2之间的Mayer键级为4.114，和上一节看到的苯与整个Pd(100)表面的片段间Mayer键级4.310相差很小，这说明当前自动定义的片段2是合理的。如果你把片段2的原子在Multiwfn主功能0里用菜单栏的Other settings - Set atom highlighting功能高亮显示，会看到它们都是离苯不远的原子。  Atoms in fragment 2:  76,77,80-83,86,87,97,101-103,106-108,112,195-197,200-202,205-207,220,221,224-22  6,230,231,241,245-247,250-252,256,265,266,269-272,274-277,280,281  Mayer bond order between fragments 1 and 2:   4.114 之后输入0.35，Multiwfn就开始用0.35 Bohr间距的立方格点计算AOM，然后构造出FOM并导出到当前目录下的FOM.txt。在双路7R32机子上96核并行计算，整个过程花了26分钟（耗时与格点间距的三次方呈反比，耗时太高的话可以用比如0.5 Bohr格点间距，误差也还可以接受，此时只需要不到10分钟）。然后输入 0  //返回 200  //其它功能（Part 2） 20  //BOD/NAdO分析 4  //直接从当前目录下的FOM.txt中读取NAdO分析要用的片段1和2的FOM 之后开始了NAdO的计算，结果如下  Eigenvalues of NAdOs: (sum=   6.95512 )    0.79476   0.77393   0.71745   0.63686   0.62698   0.35965   0.22325    0.16243   0.14986   0.14655   0.14545   0.13818   0.12618   0.12058    0.11388   0.11078   0.10617   0.09335   0.09007   0.08979   0.08423    0.08224   0.07942   0.07434   0.07391   0.07171   0.06475   0.05751    0.05507   0.05478   0.05246   0.04964   0.04669   0.04030   0.03771 ...略 当前显示的所有NAdO本征的加和6.955对应于苯和Pd(100)之间在Hirshfeld原子空间划分下算的模糊键级，和前面看到的片段间Mayer键级4.31有一定差异，这很正常，毕竟对原子空间的定义截然不同。由以上数据可见有不少NAdO轨道对模糊键级的贡献都显著，尤其是前六个。现在输入y载入新产生的NAdOs.mwfn，在主功能0里观看NAdO轨道，其中本征值最大的6个NAdO轨道如下所示，等值面数值用的0.03 从上图中可以看到这些轨道在苯与Pd(100)之间都是相位相同方式叠加的，必然都是起到明显成键作用的，所以本征值都不小。从轨道图形上可以看出这些轨道来自于苯上C原子的垂直于苯环的p轨道与Pd的原子轨道混合，而且有的图里直接就能清楚看出Pd用的是d原子轨道。例如NAdO(2)中和苯接触的Pd明显用的是dz2轨道，从Pd区域的轨道等值面形状就能看出这一点。感兴趣的话可以进一步按照《谈谈轨道成份的计算方法》（//www.umsyar.com/131）介绍的用SCPA方法做一下轨道成分分析。在Multiwfn主菜单里输入 8  //轨道成分分析 3  //SCPA方法 2  //2号NAdO轨道 Multiwfn马上就输出了轨道成份，下面列出的是Multiwfn返回的各个原子上各个角动量基函数产生的贡献。明显可以看到Pd主要都是用D基函数，而C用的都是P基函数。因此对NAdO(2)的轨道成份的分析体现了Pd的d原子轨道与C的p原子轨道的显著混合对Pd吸附苯有关键性贡献。可见以这种方式讨论成键特征能分析得巨清楚透彻！ Composition of each shell Shell  1370 Type: D    in atom  196(Pd) :     0.50345 % Shell  1405 Type: D    in atom  201(Pd) :     0.88530 % Shell  1440 Type: D    in atom  206(Pd) :     0.50251 % Shell  1538 Type: D    in atom  220(Pd) :     0.80192 % Shell  1545 Type: D    in atom  221(Pd) :     0.76004 % Shell  1570 Type: S    in atom  225(Pd) :     0.53572 % Shell  1573 Type: D    in atom  225(Pd) :    10.40186 % Shell  1577 Type: S    in atom  226(Pd) :     0.56745 % Shell  1580 Type: D    in atom  226(Pd) :    10.21554 % Shell  1608 Type: D    in atom  230(Pd) :     0.80255 % Shell  1615 Type: D    in atom  231(Pd) :     0.76058 % Shell  1713 Type: D    in atom  245(Pd) :     0.53415 % Shell  1717 Type: S    in atom  246(Pd) :     0.96339 % Shell  1720 Type: D    in atom  246(Pd) :    21.79157 % Shell  1748 Type: D    in atom  250(Pd) :     0.53493 % Shell  1752 Type: S    in atom  251(Pd) :     0.96056 % Shell  1755 Type: D    in atom  251(Pd) :    21.80729 % Shell  2019 Type: P    in atom  289(C ) :     1.48815 % Shell  2024 Type: P    in atom  290(C ) :     1.47584 % Shell  2029 Type: P    in atom  291(C ) :     6.18539 % Shell  2034 Type: P    in atom  292(C ) :     1.44475 % Shell  2039 Type: P    in atom  293(C ) :     1.41651 % Shell  2044 Type: P    in atom  294(C ) :     6.20361 % 感兴趣的读者还可以用CDA方法从苯的分子轨道与Pd表面的晶体轨道的混合角度来考察Pd对苯吸附造成的电子转移和轨道相互作用的细节。详见《使用Multiwfn结合CP2K对周期性体系做电荷分解分析（CDA）》（//www.umsyar.com/716）。 3 总结 本文通过COF二维层状体系和Pd(100)表面吸附苯作为例子，演示了如何用Multiwfn程序计算原子间和片段间的Mayer键级和模糊键级，还演示了如何计算pi键级单独考察pi电子的贡献，以及讲解了如何做NAdO分析以轨道形式清楚直观地考察原子间或片段间的共价作用的内在特征。可见使用Multiwfn可以把成键特征情况展现得超级清楚透彻。本文的做法对其它情况的周期性体系，诸如三维原子晶体、过渡态结构等等，也都是完全适用的。更多细节和相关知识请参看本文提到的相关博文。 Multiwfn对周期性体系能算的键级不止本文示例的这些。对周期性体系Multiwfn还能算Mulliken键级和Wiberg键级，还能把Mulliken和Mayer键级分解成不同轨道的贡献，后者称为轨道占据数扰动的Mayer键级。具体介绍见《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471），做法见Multiwfn手册4.8节的相关例子，本文就不特意演示了。除了本文介绍的这些外Multiwfn还有很多其它考察成键的方法对周期性体系都可以用，比如AIM拓扑分析，见《使用Multiwfn结合CP2K做周期性体系的atom-in-molecules (AIM)拓扑分析》（//www.umsyar.com/717）。 使用Multiwfn计算的结果若用于发表，记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。给别人代算的话也必须明确告知对方这一点。 //www.umsyar.com/717 2024-07-03T11:37:00+08:00 使用Multiwfn做周期性体系的atom-in-molecules (AIM)拓扑分析 Using Multiwfn to perform atom-in-molecules (AIM) topology analysis for periodic systems 文/Sobereva@北京科音   2024-Jul-3 0 前言 Bader提出的非常知名的atom-in-molecules (AIM)理论在分析孤立体系和周期性体系的电子结构方面有非常广泛的应用，相关资料汇总见《AIM学习资料和重要文献合集》（http://bbs.keinsci.com/thread-362-1-1.html）。AIM分析框架中最常见的分析手段之一是电子密度的拓扑分析，也常被称为AIM拓扑分析，例如根据这种拓扑分析得到的键临界点（BCP）的位置可以用来讨论成键和弱相互作用特征，参考《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）和《Multiwfn支持的弱相互作用的分析方法概览》（//www.umsyar.com/252）里的相关介绍。非常强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）很早以前就已经支持了结合Gaussian、ORCA等 程序产生的波函数文件对孤立体系做AIM拓扑分析，见《使用Multiwfn做拓扑分析以及计算孤对电子角度》（//www.umsyar.com/108）和Multiwfn手册4.2节的大量例子。如今Multiwfn也已支持了对周期性体系做AIM拓扑分析，本文的目的是对此进行演示，以使得读者可以轻松地举一反三充分将这种手段应用于实际研究中。使用Multiwfn做AIM拓扑分析在发表文章时记得需要按照Multiwfn启动时的提示恰当引用Multiwfn的原文。 读者务必使用2024-Jun-16及以后更新的Multiwfn版本，否则情况与本文所述不符。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载。不了解Multiwfn者建议阅读《Multiwfn入门tips》（//www.umsyar.com/167）和《Multiwfn FAQ》（//www.umsyar.com/452）。本文要用到非常流行、高效且免费的第一性原理程序CP2K，不熟悉者推荐通过“北京科音CP2K第一性原理计算培训班”（http://www.keinsci.com/workshop/KFP_content.html）系统性学习。本文要使用Multiwfn创建CP2K输入文件，这在《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）里有简要介绍。本文用的CP2K是2024.1版。 下面的例子涉及到的所有文件都可以在//www.umsyar.com/attach/717/file.rar里得到。 1 实例1：NaCl晶体 这一节对NaCl做AIM拓扑分析，要使用CP2K产生的molden文件作为Multiwfn的输入文件。molden文件的产生方法在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里有详细说明。由于考虑k点时CP2K没法产生molden文件，因此必须使用足够大的超胞模型并只考虑gamma点。本文文件包里NaCl.cif是NaCl单胞的晶体结构文件，基于它我们创建CP2K算单点任务的超胞的输入文件，同时得到molden文件。 用Multiwfn载入NaCl.cif，然后输入 cp2k  //创建CP2K输入文件的界面 [回车]  //产生的输入文件名用默认的NaCl.inp -2  //要求产生molden文件 -11  //几何操作界面 19  //构造超胞 2  //第1个方向平移复制为原先的2倍 2  //第2个方向平移复制为原先的2倍 2  //第3个方向平移复制为原先的2倍 -10  //返回 0  //产生输入文件 现在当前目录下就有了NaCl.inp，对应的是PBE/DZVP-MOLOPT-SR-GTH计算级别，对此体系是合适的。当前晶胞边长为11.28埃，对于绝缘体来说只考虑gamma点基本够用。用CP2K运行NaCl.inp，得到NaCl-MOS-1_0.molden。用文本编辑器打开NaCl.inp，将以下内容插入文件开头来定义晶胞信息从而使得Multiwfn在分析时能考虑周期性，并且定义Na和Cl的价电子数分别为9和7（不了解为什么这么设置的人看前述的《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》）。  [Cell] 11.28112000     0.00000000     0.00000000  0.00000000    11.28112000     0.00000000  0.00000000     0.00000000    11.28112000  [Nval]  Na 9  Cl 7 下面开始做AIM拓扑分析。启动Multiwfn，载入NaCl-MOS-1_0.molden，然后输入 2  //拓扑分析（默认是对电子密度做拓扑分析） 2  //以每个原子核位置为初猜搜索临界点 当前从屏幕上会看到找到了64个(3,-3)临界点，分别对应64个原子核位置的电子密度极大点。之所以虽然当前用的是赝势基组，在原子核位置还能找到电子密度极大点，是因为Multiwfn自动对这些原子指认了合适的EDF信息来描述内核电子，详见《在赝势下做波函数分析的一些说明》（//www.umsyar.com/156）。 然后再输入 3  //以每两个原子间连线中点作为搜索临界点的初猜位置 4  //以每三个原子的中心位置作为搜索临界点的初猜位置 0  //观看结果 此时文本窗口显示了各种临界点的数目，如下所示，并且提示Poincare-Hopf关系已经满足了，因此大概率所有临界点都已经找全了。  The number of critical points of each type:  (3,-3):    64,   (3,-1):   384,   (3,+1):   384,   (3,+3):    64  Poincare-Hopf verification:    64  -   384  +   384  -    64  =     0  Fine, Poincare-Hopf relationship is satisfied, all CPs may have been found 在Multiwfn的图形窗口可以看到下图，棕红色大球是Na，绿色大球是Cl，橙色、黄色、绿色小球分别是电子密度的(3,-1)、(3,+1)和(3,+3)临界点，而(3,-3)临界点都被原子球挡住了。在Multiwfn手册3.14节对Multiwfn的拓扑分析功能有完整、非常详细的说明，建议阅读。 下面再把AIM拓扑分析里常涉及的键径产生出来。在图形界面右上角点Return关闭之，然后选8产生(3,-3)和(3,-1)之间的拓扑路径，也即键径。在笔者的i9-13980HX CPU上24核并行很快就产生完了。之后再选0进入图形界面观看，看到下图，可见所有键径（棕色细线）都产生了出来，既连接相邻的Na和Cl，也连接每一对临近的两个Cl。 默认情况下Multiwfn把处于晶胞边缘的所有镜像原子和临界点也显示了出来，如果不想显示它们的话（也即只显示晶胞中唯一的那些），在图形界面顶端依次选择Other settings - Toggle showing all boundary atoms和Toggle showing all boundary CPs and paths。在文本窗口也会看到当前设置的状态。 下面来考察一下Na和Cl之间的BCP的属性。用选项0的图形界面右边的选项把其它临界点关闭显示而只显示(3,-1)临界点，选择CP labels复选框要求显示临界点序号，关闭键径显示并要求显示原子，然后就能看到下图，这是当前图像的一个边缘区域（点击save picture按钮可在当前目录下得到高分辨率图像文件，下图是剪切出的一部分） 可见189、512、238等序号的临界点都对应Na-Cl的BCP，考察哪个都可以，都是等价的。点击Return按钮关闭图形窗口，然后选择7进入临界点属性计算界面，然后再输入临界点序号189，此时就看到了这个临界点处各种函数的值，如下所示。这些函数在Multiwfn手册2.6节都有简要介绍，且在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有十分全面深入详细的讲解。 CP type: (3,-1)  Density of all electrons:  0.1144132659E-01  Density of Alpha electrons:  0.5720663296E-02  Density of Beta electrons:  0.5720663296E-02  Spin density of electrons:  0.0000000000E+00  Lagrangian kinetic energy G(r):  0.1168592222E-01  G(r) in X,Y,Z:  0.7296255010E-03  0.1022667133E-01  0.7296253875E-03  Hamiltonian kinetic energy K(r): -0.3025433987E-02  Potential energy density V(r): -0.8660488230E-02  Energy density E(r) or H(r):  0.3025433987E-02  Laplacian of electron density:  0.5884598270E-01  Electron localization function (ELF):  0.1993303578E-01  Localized orbital locator (LOL):  0.1249064363E+00 ...略 如我在全面提到的《Multiwfn支持的分析化学键的方法一览》文中所述，BCP位置上电子密度拉普拉斯函数和能量密度都为正是离子键的典型特征。当前BCP这两个值分别为0.0588 a.u.和0.00302 a.u.，都为正值，完全体现出Na-Cl是离子键的实事。此外，当前BCP的电子密度仅为0.01144 a.u.、ELF仅为0.0199，BCP处很小的电子密度和ELF值也是典型的离子键的普遍共性。 我在《使用Multiwfn+VMD快速地绘制高质量AIM拓扑分析图（含视频演示）》（//www.umsyar.com/445）介绍了一种极好的将Multiwfn和VMD联用快速、理想地绘制AIM拓扑分析图的做法，比在Multiwfn里直接显示的效果更好，而且视角可以自由旋转，没看过的话务必先看一下。对于当前体系，也可以像此文一样用Multiwfn文件包里的examples\scripts\AIM.bat和AIM.txt来作图，得到的图像如下。 上图效果明显不好，显得残缺不全，这是因为此时显示的原子、临界点、键径都只有晶胞里唯一的那些，而处于棱和面上的边界原子、临界点、键径的周期镜像都没有显示出来。为了避免这个问题，对于当前体系有原子、临界点、键径处于棱和面上的情况，应当改用examples\scripts\目录下的AIM+boundary.bat和AIM+boundary.txt来绘制，其用法与AIM.bat和AIM.txt的组合完全一致。只不过前者传递给Multiwfn的指令要求Multiwfn导出mol.pdb、CPs.pdb和paths.pdb文件时把边界的镜像也都纳入进去。利用AIM+boundary.bat和AIM+boundary.txt来绘制的结果如下，可见此时的效果就很好了，和前面在Multiwfn里看到的完全相同。注意此时这三个pdb文件里只有序号靠前的那些才是晶胞中唯一的那些。 2 层状COF 下面再举个例子，用Multiwfn对一个层状的共价有机框架化合物（COF）体系做AIM拓扑分析。此体系的cif文件是本文文件包里的COF-12000N2.cif，其结构图如下所示。 由于此体系的晶胞尺寸不小，计算时可以直接对原胞做只考虑gamma点的计算来得到molden文件。如《实验测定分子结构的方法以及将实验结构用于 计算需要注意的问题》（//www.umsyar.com/569）所说，由于X光衍射测的晶体结构中氢的原子位置通常不准确，最好做分析前先冻结重原子而优化一下氢原子的坐标，但此例为了省事就跳过这一步了。如果你要用AIM拓扑分析考察氢键，由于氢的位置是重中之重，那显然是绝对不能略过这一步的。 当前体系用PBE结合CP2K的GTH赝势基组计算完全可以，而此例作为演示，这回使用PBE结合pob-TZVP-rev2全电子基组计算。启动Multiwfn，载入COF-12000N2.cif，然后输入 cp2k  //创建CP2K输入文件的界面 [回车]  //产生的输入文件名用默认的COF-12000N2.inp -2  //要求产生molden文件 2  //选择基组 15  //pob-TZVP-rev2 4  //开OT，对当前用的基组通常比对角化明显更容易收敛 0  //产生输入文件 用CP2K计算COF-12000N2.inp，然后将以下晶胞设置信息插入到计算得到的COF-12000N2-MOS-1_0.molden文件的开头  [Cell]  22.55600000     0.00000000     0.00000000 -11.27800000    19.53406901     0.00000000   0.00000000     0.00000000     6.80000000 启动Multiwfn，载入COF-12000N2-MOS-1_0.molden，然后输入 2  //拓扑分析（默认是对电子密度做拓扑分析） 2  //以每个原子核位置为初猜位置搜索临界点 3  //以每两个原子间连线中点作为搜索临界点的初猜位置 4  //以每三个原子的中心位置作为搜索临界点的初猜位置 8  //产生(3,-3)到(3,-1)的拓扑路径 0  //观看结果 现在看到下图。可见临界点和键径都成功地得到了，不仅有化学键对应的BCP和键径，COF内氢键N-H...O作用以及层之间相互作用相对应的BCP和键径都有。而且不仅有晶胞内两层COF之间相互作用的BCP和键径，当前晶胞里的COF和相邻镜像晶胞里的COF相互作用对应的BCP和键径也都有，显然周期性完全充分考虑了。 虽然Multiwfn的文本窗口提示当前并未符合Poincare-Hopf关系，但由于感兴趣的临界点和键径都有了，所以就没必要进一步搜索了。 再利用前述名为AIM+boundary.bat的Windows批处理脚本结合AIM+boundary.txt里记录的指令，用Multiwfn+VMD绘制体系结构+临界点+键径图。然后再在VMD的文本窗口里输入pbc box命令显示盒子，VMD main界面里选择Display - Orthographic用正交视角。此时VMD显示的图像如下，可见非常理想！ 顺带一提，对于当前体系要想完整展现层间相互作用的话，最理想的做法是用我提出的IGMH，见《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（//www.umsyar.com/621）和《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（//www.umsyar.com/667）。 3 对Ag晶胞基于电子密度格点数据做AIM拓扑分析 Multiwfn不仅可以基于CP2K产生的波函数文件做AIM拓扑分析，还可以基于从cube等文件里载入的电子密度的格点数据做拓扑分析，此时各个位置的电子密度、梯度和Hessian矩阵基于格点数据以B-spline方式插值得到。由于这个特性，Multiwfn做周期性体系的AIM拓扑分析不限于CP2K用户，诸如从VASP产生的CHGCAR文件里读取电子密度格点数据做AIM拓扑分析也可以。Multiwfn支持哪些格点数据文件格式见手册2.5节的说明。而且CP2K用户还可以在考虑k点的情况下对原胞计算得到电子密度cube文件结合Multiwfn做AIM拓扑分析，而不必在晶胞较小时非得弄成超胞模型再算。但基于电子密度格点数据做AIM拓扑分析的关键不足是对临界点只能得到电子密度及其导数信息以及一些与之直接相关的量，而依赖于波函数才能计算的如动能密度、ELF等函数就没法计算了。另外，如果电子密度格点数据是在赝势下计算产生的，基于它插值出的电子密度函数中，原子核处将没有电子密度极大点，也因此无法产生键径。 此例演示用Multiwfn对Ag的单胞基于电子密度cube文件做AIM拓扑分析的过程。使用的级别是PBE/TZVPP-MOLOPT-PBE-GTH-q19，这个基组是CP2K 2024.1的data目录下的BASIS_MOLOPT_UZH基组文件里的。之所以不用常用的BASIS_MOLOPT里的Ag的基组，在于那里面Ag的基组都是-q11的，即只描述11个价电子。使用赝势的情况下做AIM拓扑分析时，被描述的价电子越多，即使用越小核赝势，结果越准确、越可能接近全电子计算的结果。 启动Multiwfn，载入本文文件包里的Ag单胞的晶体结构文件Ag.cif，然后输入 cp2k  //创建CP2K输入文件的界面 [回车]  //产生的输入文件名用默认的Ag.inp -3  //要求产生cube文件 1  //电子密度 6  //开smearing 8  //设置k点 10,10,10 0  //产生输入文件 把Ag.inp里BASIS_SET_FILE_NAME设为BASIS_MOLOPT_UZH，POTENTIAL_FILE_NAME设为POTENTIAL_UZH。Ag的BASIS_SET设为TZVPP-MOLOPT-PBE-GTH-q19，POTENTIAL设为GTH-PBE-q19。然后用CP2K计算，之后会得到记录晶胞内价电子密度的Ag-ELECTRON_DENSITY-1_0.cube。如果不了解cube格式的话，看《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）。用文本编辑器打开这个cube文件，在第一行里写上box2cell。这样Multiwfn载入这个cube文件时，一旦发现第一行的内容包含box2cell，就会把格点数据的盒子信息转化为晶胞信息，从而在计算的时候能够考虑周期性（也可以不修改cube文件，而是载入cube文件后用隐藏的主功能1000里的子功能18来实现盒子信息向晶胞信息的转化）。 把Multiwfn的settings.ini里的iuserfunc设为-3，使得user-defined function对应基于载入的格点数据通过B-spline插值的函数。然后启动Multiwfn，载入Ag-ELECTRON_DENSITY-1_0.cube，之后输入 2  //拓扑分析 -11  //修改被分析的函数 100  //User-defined function，当前对应于插值方式产生的价层电子密度 3  //将每两个原子之间的中点作为搜索临界点的初猜 然后进入选项0，并把原子标签显示出来、把键设细，就会看到下图。可见临近原子之间的BCP都已经找到了 下面再把更多的临界点也找出来。当前情况不适合选择选项4、5分别用3、4个原子中心作为初猜点来搜索临界点。选项3用每两个原子之间中点当做初猜点时是考虑循环相邻晶胞中的镜像原子的，而选项4、5则只会循环当前晶胞里那些唯一的原子，因此此例用选项4、5搜索后还是会遗漏许多重要的临界点。当前最适合的方式是选择选项6里的子选项-1，此时会在每个原子附近特定半径内撒一批初猜点来搜索临界点。这么选过之后，会发现找到了一大堆新的临界点，进入选项0后看到的图如下所示。下图左侧把所有临界点都显示了出来，可以看到离原子核较近的区域临界点有不少，而且各种类型的都有。这是因为当前用了赝势，由于电子密度不是从原子核开始向四周单调下降的，因此会在价层区域出现各种各样的临界点，这些不是我们感兴趣的。把原子球设大掩盖这些没什么意义的临界点，并且旋转体系使得视角与晶轴平行后，就看到了下面的右图。可见在原子间相互作用区域找到的临界点的分布整齐、对称，因此可以认为有意义的临界点都找到了。 有一个概念叫平面性指数，在“ 波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有一页幻灯片做了介绍，如下所示。我们这里来对Ag晶体计算一下这个平面性指数，看看是否和文章里说的情况吻合。注意这里说的平面性指数和《使用Multiwfn定量化和图形化考察分子的平面性（planarity）》（//www.umsyar.com/618）里介绍的我提出的衡量体系几何结构平面性的参数是两码事。 首先我们得确定要考察的笼临界点也即(3,+3)的临界点编号。下图左边只显示了笼临界点，可见当前体系有两类笼临界点，一种是诸如晶胞中心的那个（15号），它处于八面体中心。另一种是处于四面体中心的，比如63号临界点。下图右边只显示了所有BCP，此体系中所有的BCP都是等价的，我们考虑其中任意一个，比如1号。利用当前拓扑分析界面里的选项7依次对这三个临界点的属性进行考察，可以得到的这三个位置的基于格点数据插值产生的电子密度，也即User-defined real space function后面的值，结果如下（注意不要读Multiwfn直接输出的Density of electrons后面的值，那个是用孤立状态原子的电子密度叠加得到的准分子密度） CP(15): 0.1572533472E-01 a.u. CP(63): 0.2212960297E-01 a.u. CP(1): 0.2706526462E-01 a.u. 可见两类笼临界点里电子密度最小值是0.01572 a.u.，BCP处的电子密度为0.02706 a.u.，因此平面性指数为0.01572/0.02706=0.58，在0.5左右，符合幻灯片里说的“其它金属及合金”的情况。 4 总结&其它 本文通过三个有代表性的例子详细介绍了在Multiwfn中对周期性体系做电子密度的拓扑分析，也即一般说的AIM拓扑分析的完整过程，所有要注意的细节都做了解释。通过例子可见，Multiwfn做周期性体系的拓扑分析的过程很简单，并不比分析孤立体系更复杂。读者如果实际操作一遍，还会感受到Multiwfn做这些分析的速度相当快。AIM拓扑分析对于考察周期性体系的电子结构有重要意义，推荐读者在实际研究中运用，使文章增光添彩。 很值得一提的是，Multiwfn的拓扑分析功能极其普适，绝不仅限于电子密度的拓扑分析。对于Multiwfn自身可以直接计算的函数，以及载入的格点数据里记录的函数，Multiwfn都可以做拓扑分析，对分子体系的例子参看比如《使用Multiwfn对静电势和范德华势做拓扑分析精确得到极小点位置和数值》（//www.umsyar.com/645）对静电势、范德华势的拓扑分析，以及《在Multiwfn中单独考察pi电子结构特征》（//www.umsyar.com/432）里面对ELF-pi函数的拓扑分析。对周期性体系，Multiwfn也完全可以以相同的方式对任意函数做拓扑分析。 使用本文的做法做AIM拓扑分析在发表文章时必须按照Multiwfn启动时的提示恰当引用Multiwfn，如果用于给别人代算也要明确告知对方这一点。 //www.umsyar.com/716 2024-07-01T22:17:00+08:00 使用Multiwfn结合CP2K对周期性体系做电荷分解分析（CDA） Using Multiwfn in combination with CP2K to perform charge decomposition analysis (CDA) for periodic systems 文/Sobereva@北京科音  2024-Jul-1 1 前言 笔者之前写过《使用Multiwfn做电荷分解分析(CDA)、绘制轨道相互作用图》（//www.umsyar.com/166）介绍了怎么用Multiwfn程序基于Gaussian等 程序产生的波函数文件对分子体系做电荷分解分析（CDA）。这种方法可以从片段轨道间相互作用的角度深入了解片段间是如何转移电子的。由于CDA的很高价值，Multiwfn做CDA分析已经得到非常广泛的使用。如果读者没看过此//www.umsyar.com/166的话一定要先仔细看一下，否则无法理解下文的内容。 从2024-Jun-5更新的Multiwfn开始，CDA模块已经支持了基于CP2K产生的molden文件做周期性体系的CDA分析，下面将通过一个简单的例子进行演示。Multiwfn可以在主页//www.umsyar.com/multiwfn免费下载，不了解此程序者参看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。我假定读者已经具备了CP2K的基本使用知识，不了解者推荐通过北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）完整系统地学习。 本文的例子是计算NaCl板吸附CO分子的体系，用CDA方法考察CO与NaCl板之间电子转移情况。这个体系优化后的结构如下，C和O分别是1和2号原子，3-110号原子是NaCl板。对这个体系，之前我在《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（//www.umsyar.com/638）中通过密度差图的方式定性考察了电子转移，在《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）里通过片段电荷的方式定量考察了电子转移。而下面用CDA分析，可以以轨道角度在明显更深层次理解电子转移细节。 本文例子涉及的所有文件都在//www.umsyar.com/attach/716/file.7z里。CP2K用的版本是2024.1，Multiwfn是2024-Jun-5更新的版本，读者绝对不要用Multiwfn的更老版本。 2 产生输入文件 首先我们要用CP2K计算NaCl+CO、NaCl、CO各自的molden格式的波函数文件。如果你不知道怎么用CP2K产生molden文件的话，先仔细阅读《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651），我假定读者已经看过本文了。这里有几个要点： (1)整体和片段的molden文件里的原子坐标必须对应 (2)整体和片段的计算必须使用严格相同的级别和设置，用的盒子也应当一致 (3)必须要求CP2K计算所有空轨道 (4)整体的molden文件里需要加入[Cell]字段定义晶胞信息，而片段的molden文件里加不加入不影响CDA结果。[Nval]字段可加入可不加入，也不影响CDA结果 本文文件包里的NaCl_CO.cif是以前我在PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH下优化的NaCl板吸附CO的结构文件。用Multiwfn载入它，然后输入 cp2k  //创建CP2K输入文件 NaCl_CO.inp  //产生的输入文件名 -2  //要求产生molden文件 -9  //其它设置 12  //设置计算的空轨道数 -1  //计算所有空轨道 0  //返回 0  //产生输入文件 用CP2K运行NaCl_CO.inp，得到NaCl_CO-MOS-1_0.molden。然后将以下字段插入到molden文件的开头 [Cell] 16.92168000     0.00000000     0.00000000  0.00000000    16.92168000     0.00000000  0.00000000     0.00000000    25.00000000 将NaCl_CO.inp复制为NaCl.inp，再将NaCl.inp中&COORD里的CO部分删掉、PROJECT NaCl_CO改名为PROJECT NaCl，然后用CP2K运行之，得到NaCl-MOS-1_0.molden。将上面的[Cell]加入其中（这和CDA分析无关，加入这个的目的是为了之后在Multiwfn中能正确观看NaCl板的周期性的晶体轨道）。 将NaCl_CO.inp复制为CO.inp，再将CO.inp中&COORD里的NaCl部分删掉、PROJECT NaCl_CO改名为PROJECT CO，然后用CP2K运行之，得到CO-MOS-1_0.molden。这个可以不用刻意加上[Cell]，因为当前体系中CO离盒子边界很远，显示轨道时是否考虑周期性对看到的图像没影响。 接下来就可以进行CDA分析了。 3 进行CDA分析 启动Multiwfn，载入NaCl_CO-MOS-1_0.molden，然后依次输入 16  //CDA分析 2  //两个片段 CO-MOS-1_0.molden  //CO的波函数文件。注意要按整个体系里原子序号顺序载入片段 NaCl-MOS-1_0.molden  //NaCl板的波函数文件 现在屏幕上输出了一大堆复合物轨道的CDA结果，同时还输出了ECDA结果： The net electrons obtained by frag. 2 = CT( 1-> 2) - CT( 2-> 1) =    0.1580 这个ECDA结果告诉你CO（片段1）往NaCl板（片段2）转移了0.158个电子，这和《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）里通过Mulliken方法计算的CO的片段电荷严格相同。 当前CDA分析输出的复合物轨道太多，一共437个（对应整个体系的占据轨道数），而绝大多数对d、b、r项的贡献微乎其微或完全为0，没有考察的意义。为了简化输出，在当前界面里输入 -3  //设置CDA输出时d、b、r项的阈值，它们之一的绝对值大于这个阈值的复合物轨道才会输出 0.005 0  //显示CDA和ECDA结果 当前结果如下    Orb.      Occ.          d           b        d - b          r     272    2.000000    0.006220   -0.000012    0.006232    0.000488     275    2.000000    0.125862   -0.001874    0.127736    0.051155     297    2.000000    0.001098    0.000305    0.000792   -0.005373     318    2.000000    0.001053    0.000326    0.000727   -0.005720     341    2.000000    0.000781    0.000200    0.000581   -0.006708     352    2.000000    0.002291    0.000731    0.001560   -0.015705     360    2.000000    0.001101    0.000298    0.000803   -0.009297     414    2.000000    0.000635    0.000214    0.000421   -0.005229     427    2.000000    0.000739    0.000219    0.000520   -0.005731 ------------------------------------------------------------------- Sum:     874.000000    0.153805    0.017732    0.136074   -0.037831 Note: The "Sum" includes all terms including those not printed above CDA的总结果是CO向NaCl转移了d=0.154个电子，NaCl向CO反馈了b=0.018个电子，NaCl净获得d-b=0.136个电子。b项基本可以忽略不计，而d的主要贡献是由275号复合物轨道造成的。这个复合物轨道的贡献显然值得进一步探究，看看是由CO和NaCl板的哪些片段轨道混合所造成的。为此，输入 5  //将特定复合物轨道对CDA项的贡献进行分解 275  //要考察的复合物轨道序号 0.005  //输出阈值 现在看到如下结果 Occupation number of orbital   275 of the complex:  2.00000000 FragA Orb(Occ.)  FragB Orb(Occ.)      d           b        d - b          r    5( 2.0000)     271( 2.0000)    0.000000    0.000000    0.000000    0.005183    5( 2.0000)     289( 2.0000)    0.000000    0.000000    0.000000    0.005642    5( 2.0000)     336( 2.0000)    0.000000    0.000000    0.000000    0.005127    5( 2.0000)     341( 2.0000)    0.000000    0.000000    0.000000    0.008960    5( 2.0000)     435( 0.0000)    0.009358    0.000000    0.009358    0.000000    5( 2.0000)     439( 0.0000)    0.011934    0.000000    0.011934    0.000000    5( 2.0000)     443( 0.0000)    0.010929    0.000000    0.010929    0.000000    5( 2.0000)     447( 0.0000)    0.008032    0.000000    0.008032    0.000000    5( 2.0000)     453( 0.0000)    0.009177    0.000000    0.009177    0.000000    5( 2.0000)     465( 0.0000)    0.010755    0.000000    0.010755    0.000000    5( 2.0000)     469( 0.0000)    0.013735    0.000000    0.013735    0.000000    5( 2.0000)     473( 0.0000)    0.011282    0.000000    0.011282    0.000000 Sum of above terms:               0.085201    0.000000    0.085201    0.024911 可以看到CO向NaCl转移电子靠的都是CO的5号占据轨道和NaCl空轨道相互作用。NaCl板接收CO转移来的电子分散在其大量空轨道上，比如NaCl板的443号空轨道接受了0.011个电子。在当前的输出阈值下，上面输出的这些项的所有的d项加和为0.085，明显和275号复合物轨道的总d值0.153有很大差距，说明CO还向上面没输出的很多其它的NaCl板的空轨道转移了电子。 下图把275号复合物轨道、CO的5号占据轨道，以及NaCl接收CO贡献来的电子最多的两个轨道（439、469）展示在了一起，便于大家直观了解情况，等值面数值用的是0.02。复合物和片段的轨道可以用Multiwfn载入相应的波函数文件后用主功能0按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）说的观看。从此图可以清晰直观地认识到CO的孤对电子往吸附CO的那个Na位点的空轨道上转移了电子。还有大量也具有这种特征的NaCl板的空轨道也接收了CO的孤对电子，这里没全都画出来。 下面再看一下上图中275号复合物轨道是怎么由片段轨道构成的。按0返回到CDA菜单，然后进入选项2，输入275，看到以下信息  Note: Only the fragment orbitals with contribution >  1.0 % will be shown below, the threshold can be changed by "compthresCDA" in settings.ini  Occupation number of orbital   275 of the complex:  2.00000000  Orbital     5 of fragment  1, Occ: 2.00000    Contribution:   87.42 %  Sum of values shown above:     87.42 % 这告诉你275号复合物轨道由CO的5号占据轨道贡献了87.4%。默认情况下只有贡献大于1%的片段轨道才会输出，虽然NaCl板的一大堆非占据轨道都通过与CO的5号轨道混合而参与了这个复合物轨道的构成，但贡献都小于输出阈值，所以上面没看到。从当前信息可以认识到275号复合物轨道主要体现CO的5号轨道特征，这是为什么上面图中275号复合物轨道和CO的5号轨道看起来十分相似。 4 总结 此文通过一个简单的表面吸附的例子演示了怎么用Multiwfn的CDA功能分析周期性体系的电荷转移本质。讨论电子转移常见的套路就是绘制电子密度差、算算片段电荷，而从本文的例子看到Multiwfn可以从片段轨道相互作用对电子转移的本质进行深入的考察，掌握了本文介绍的方法明显可以给分析增光添彩。 使用Multiwfn做本文介绍的分析请记得在发表文章时引用Multiwfn程序启动时提示的Multiwfn原文，以及引用介绍CDA/GCDA方法的文章（Multiwfn用的实质上是我对CDA广义化后的GCDA形式），这在进入Multiwfn的CDA功能时屏幕上明确提示了。 //www.umsyar.com/712 2024-06-08T15:43:00+08:00 使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷 Using Multiwfn to calculate Hirshfeld(-I), CM5 and MBIS atomic charges for periodic systems 文/Sobereva@北京科音   2024-Jun-8 1 前言 原子电荷（atomic charge）即原子所带的净电荷，对于讨论化学体系中原子的带电状态有重要意义，作为原子的一种描述符它也和原子的很多属性有密切联系。强烈建议阅读《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（//www.umsyar.com/714）里提到的笔者写的原子电荷的综述以全面了解原子电荷。 强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）支持非常丰富的原子电荷计算方法，其中多数已经支持周期性体系。Multiwfn支持的方法中有一类是基于模糊式原子空间划分，也就是每个原子对应一个平滑的权重函数，权重函数与体系的电子密度的乘积的全空间积分就是这个原子带的电子数（原子的布居数），原子的核电荷数减去它就是原子电荷。Multiwfn支持的这类方法包括Hirshfeld、Hirshfeld-I、MBIS、Becke，以及对Hirshfeld电荷做后校正得到的ADCH和CM5电荷。从2024-May-25更新的Multiwfn版本开始，Hirshfeld、CM5、Hirshfeld-I、MBIS都已经支持了周期性体系，本文将具体介绍怎么用Multiwfn与非常流行、高效且免费的第一性原理程序CP2K相结合非常方便快速地计算这些原子电荷。Multiwfn计算这些电荷的功能是普适的，并不限于结合CP2K，也可以基于其它第一性原理程序产生的体系的电子密度的cube文件来算，还支持VASP产生的记录电子密度的CHGCAR文件。 注：Multiwfn更老的一些版本也支持基于周期性波函数计算Hirshfeld和CM5电荷，但那时候用的算法对于周期性体系效率极低，本文介绍的做法与之有天壤之别。 下文第2节简要介绍一下Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷的基本特点，第3节介绍用Multiwfn+CP2K计算它们的具体方法，第4节给出具体例子，同时还会介绍怎么计算片段电荷。使用Multiwfn计算原子电荷在发表文章时记得需要按Multiwfn启动时的提示恰当引用Multiwfn。 Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载，如果对Multiwfn不了解的话建议看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。如果你用CP2K创建给Multiwfn算原子电荷用的输入文件的话，笔者假定你对CP2K的使用已经有基本且正确的了解，不具备这些知识的话非常建议通过北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）完整、系统地学习一遍。本文的例子利用到Multiwfn创建CP2K的输入文件，相关介绍见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）。 2 Hirshfeld、Hirshfeld-I、CM5和MBIS等原子电荷简介 Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷在Multiwfn手册3.9节的相应小节里都有很详细的介绍，因此这里不介绍细节，只是简单说一下它们的特点。 Theor. Chim. Acta (Berl.), 44, 129 (1977)中提出的Hirshfeld电荷是被广为使用的原子电荷，计算方式简单，物理意义清楚，在 研究孤立体系和第一性原理研究周期性体系的领域中都用得很多。Hirshfeld电荷为人诟病的地方是原子电荷数值整体严重偏小，对偶极矩、静电势重现性很差，这在我写的《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）里的对比测试中有充分体现。 JCP, 126, 144111 (2007)提出的Hirshfeld-I是在Hirshfeld基础上的改进，通过迭代方式不断更新原子权重函数直到收敛，这使得原子权重函数所描述的原子空间可以响应实际化学环境。其原子电荷数值显著大于Hirshfeld，对静电势的重现性也有明显改进。代价是Hirshfeld-I需要做迭代，往往几十轮，计算耗时明显高于Hirshfeld，而且需要事先提供体系中各种元素的各个氧化态的原子径向密度信息，实现起来较麻烦（为了实现Hirshfeld-I，我计算了360多个原子的径向密度文件并在Multiwfn中自带）。 一定要注意，虽然CP2K程序自己也支持Hirshfeld-I电荷的计算，但我测试发现其结果完全不靠谱，我检查了其源代码也确认了其实现方式明显不符合Hirshfeld-I的标准定义，因此切勿用CP2K自身的Hirshfeld-I电荷计算功能。CP2K算的Hirshfeld电荷在结合SHAPE_FUNCTION DENSITY选项时的结果不那么离谱，但结果和Multiwfn算的也存在不可忽视的差异，应当以Multiwfn的结果为准。 JCTC, 8, 527 (2012)中提出的CM5电荷和我在J. Theor. Comput. Chem., 11, 163 (2012)中提出的ADCH电荷一样都是对Hirshfeld电荷的后校正，对大部分情况都能使得它的原子电荷变得更大、更符合化学直觉，对静电势的重现性也明显变得更好。ADCH没有CM5那么强的经验性、不牵扯一堆经验参数，因而原理更为理想、更推荐使用。CM5最初是用于分子体系的，如今在应用于晶体方面也已经有了一定探索，例如JCTC, 16, 5884 (2020)将CM5与其它一些原子电荷对比考察了它们对一些固体的电荷分布的描述。CM5还被用于一些计算化学领域的方法，比如《比SMD算溶解自由能更好的溶剂模型uESE的使用》（//www.umsyar.com/593）里介绍的uESE溶剂模型依赖于CM5电荷，CM5乘以1.2后适合结合OPLS-AA力场做动力学模拟，见《计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本》（//www.umsyar.com/585）。 ADCH电荷在分子体系的研究方面已被广为使用，之前有人问我是否会把ADCH电荷扩展到周期性体系。我目前没这个打算，因为ADCH电荷的思想对周期性体系往往不适用。ADCH方法将原子偶极矩展开为周围原子的校正电荷，但对于诸如NaCl这样的体系，每个原子所处的环境是中心对称的，因此原子偶极矩为0、校正电荷都为0，故ADCH电荷会与Hirshfeld电荷完全相同。 MBIS (Minimal Basis Iterative Stockholder)电荷于JCTC, 12, 3894 (2016)中提出。类似于Hirshfeld-I，此方法也是通过迭代方式（因此耗时较高）优化原子空间直到收敛。MBIS比Hirshfeld-I的一大好处是不需要利用各个元素各个氧化态的原子径向密度信息，因此实现起来简单得多，也更优雅。原文里的测试体现出MBIS比Hirshfeld-I整体更具有优势。但MBIS经常收敛较慢，对于某些体系甚至需要200多轮才能充分收敛。目前Multiwfn里MBIS最高支持到Rn元素。 根据我的测试，对于大多数固体体系，原子电荷大小的关系是MBIS和Hirshfeld-I最大、谁大不一定，CM5大小介中，Hirshfeld电荷最小。如果你对数值的绝对大小不那么在乎，只是想体现相对大小关系的话，最“朴素”且很便宜的Hirshfeld也可以用。但如果对定量数值关心的话，建议用其它的。CM5、Hirshfeld-I和MBIS用于晶体哪个最理想没有定论，都可以算一算看看，如果倾向于得到较大数值的结果可以优先考虑Hirshfeld-I和MBIS。如果其中有的原子电荷明显不合理那就弃掉换别的方法，比如O的电荷如果算出来-2.34那就不应该用，因为显然O最多只能带-2电荷。 特别要强调的是，绝对不要把原子电荷与氧化态搞混，也切勿用原子电荷判断氧化态！氧化态怎么正确地判断看《使用Multiwfn结合CP2K计算晶体中原子的氧化态》（//www.umsyar.com/711）和《使用Multiwfn通过LOBA方法计算氧化态》（//www.umsyar.com/362）。 还有很多其它原子电荷计算方法，诸如Mulliken电荷、Lowdin电荷、SCPA电荷、AIM电荷（被个别文章称为Bader电荷）等，Multiwfn也都支持将它们用于周期性体系，这不属于本文的范畴。Mulliken电荷结合CP2K里常用的MOLOPT系列和pob系列基组用于周期性体系也说得过去。虽然它是基于希尔伯特空间划分的原子电荷中最原始、思想最朴素的一种，但实际结果一般还算定性正确，而且计算耗时极低。被用得很多的AIM电荷实际上很糟糕，往往出现违背化学常识的结果，这在前述的《原子电荷计算方法的对比》里已经充分体现了，不建议用。 3 Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷在Multiwfn中的计算方法 3.1 CP2K用户的情况 首先说对于CP2K用户怎么计算。Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷都是完全基于电子密度计算的，不直接牵扯到波函数，在Multiwfn中计算它们有两种可以用的输入文件： (1)用CP2K产生的体系的电子密度cube文件作为输入文件。这种情况在CP2K计算时可以考虑k点，是我最推荐的，耗时非常低。如果不了解cube文件的格式的话，参看《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）。注意CP2K有bug（起码直到我撰文时最新的CP2K 2024.1为止），产生的cube文件里的原子有效电荷都为0，这会导致Multiwfn无法正确计算原子电荷。这里说的原子的有效核电荷是指当前基组描述的原子的电子数，对于赝势基组它对应价电子数。为此，要么用文本编辑器手动改cube文件中每个原子的有效核电荷，要么把cube文件第一行写成各个元素的有效核电荷，例如N 5 B 3 Ti 12，这代表当前体系中的N、B、Ti的有效核电荷数分别为5、3、12 (2)用CP2K产生的记录了体系波函数的molden文件作为输入文件。详见《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）。这个方法的缺点在于没法考虑k点，对于小晶胞体系必须先扩胞到足够大，从而在计算时只考虑gamma点，这显然使得在CP2K中和Multiwfn中的计算耗时都会比用(1)的方式高很多 用Multiwfn载入上述输入文件其中一种后，进入主功能7，然后选择相应选项即可计算Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷。对于molden文件当输入文件时，Multiwfn还会让你输入计算用的格点间距，格点间距越小耗时越高而计算精度越高，通常用0.2 Bohr就可以，是精度和耗时的较好权衡。而使用cube文件做为输入文件时，计算原子电荷用的格点间距与此文件的格点间距直接对应，CP2K的CUTOFF设得越大则产生的cube文件里的格点间距会越小。 由于Multiwfn计算上述电荷用的是均匀分布的格点，不可能准确积分较重原子的内核区域电子密度，因此上述输入文件应当是CP2K做GPW计算得到的，也需要用赝势基组，此时cube文件记录的电子密度或者molden文件里记录的波函数只对应于价电子。 顺带一提，如果你要计算Mulliken电荷，用Multiwfn载入上述的molden文件，然后进主功能7，依次选择5、1即可得迅速到结果。Multiwfn中计算AIM电荷的方法笔者以后会另文说明。 3.2 其它程序用户的情况 如果你是VASP用户，可以用记录晶胞中价电子密度的CHGCAR文件给Multiwfn用于计算前述原子电荷，等同于3.1节说的第2类输入文件。只要文件名里包含CHGCAR字样就会被Multiwfn当做是CHGCAR文件来载入。注意为了让Multiwfn得知各类原子的有效核电荷数，CHGCAR文件的第一行需要以Nval为开头，后面依次写上各类元素的名字和有效核电荷，并以空格分隔。比如写为Nval Na 1 Cl 7，这代表Na和Cl的有效核电荷数分别为1和7。 如果你是其它第一性原理程序的用户，若程序可以产生记录晶胞中价电子密度的cube文件（即等同于3.1节说的第2类输入文件），也都可以用于给Multiwfn计算前述的原子电荷。 之后在Multiwfn中的计算过程和CP2K用户没任何区别。 4 实例 下面给出Multiwfn结合CP2K计算Hirshfeld、Hirshfeld-I、CM5和MBIS原子电荷的具体例子。相关的文件都可以在//www.umsyar.com/attach/712/file.zip中得到（4.2、4.3和4.5节的cube/molden文件除外，因为文件过大）。 4.1 TiO2晶体 此例对TiO2的金红石晶型计算前述的原子电荷。它的实验晶体结构的cif文件是本文文件包里的TiO2-Rutile.cif，结构图如下所示 先用CP2K计算出这个晶胞的电子密度的cube文件，用于之后Multiwfn做原子电荷计算。首先用Multiwfn创建CP2K输入文件，启动Multiwfn并载入TiO2-Rutile.cif，然后输入 cp2k  //创建CP2K的输入文件 [回车]  //产生的文件名用默认的TiO2-Rutile.inp -3  //产生cube文件 1  //电子密度 8  //设置k点 4,4,6  //与各个晶胞边长相乘皆为18埃左右，够用了 0  //产生输入文件，计算级别为默认的PBE/DZVP-MOLOPT-SR-GTH 用CP2K运行刚得到的TiO2-Rutile.inp，马上就算完了，得到了TiO2-Rutile-ELECTRON_DENSITY-1_0.cube文件，在本文的文件包里已经提供了。用文本编辑器打开此文件，把第一行改为Ti 12 O 6并保存。这代表Ti和O的有效核电荷分别为12和6，直接对应于输入文件里它们分别用的DZVP-MOLOPT-SR-GTH-q12和DZVP-MOLOPT-SR-GTH-q6基组名里面的q后面的值。 启动Multiwfn，输入 TiO2-Rutile-ELECTRON_DENSITY-1_0.cube  //写此文件的实际路径 y  //从第一行载入各个元素有效核电荷 7  //布居分析与原子电荷计算 1  //Hirshfeld电荷 立刻看到如下Hirshfeld电荷计算结果。由于格点间距不是无穷小，因此等价的原子算出来的原子电荷也会有轻微的数值差异，可以自己手动取平均。如果用更大的CUTOFF，则cube文件的格点间距会更小，不等价性会更低。 Final atomic charges: Atom    1(Ti):     0.59349459 Atom    2(O ):    -0.29461903 Atom    3(O ):    -0.29472601 Atom    4(Ti):     0.58363165 Atom    5(O ):    -0.29385236 Atom    6(O ):    -0.29385236 现在Multiwfn还问你是否把原子电荷导出为chg文件，这里选n不导出。 接下来计算CM5电荷。在Multiwfn主功能7里选择16，马上就得到了结果                       ---------- CM5 charges ---------- Atom:    1Ti  CM5 charge:    1.348393  Hirshfeld charge:    0.593495 Atom:    2O   CM5 charge:   -0.672093  Hirshfeld charge:   -0.294619 Atom:    3O   CM5 charge:   -0.672191  Hirshfeld charge:   -0.294726 Atom:    4Ti  CM5 charge:    1.338580  Hirshfeld charge:    0.583632 Atom:    5O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852 Atom:    6O   CM5 charge:   -0.671306  Hirshfeld charge:   -0.293852 Summing up all CM5 charges:     0.00007649 可见CM5电荷的大小明显大于Hirshfeld电荷，看起来更合理，而Hirshfeld电荷值明显偏小。 下面再演示计算Hirshfeld-I电荷。把Multiwfn自带的examples目录下的atmrad目录挪到当前目录下，这样Multiwfn计算Hirshfeld-I电荷时就会自动用此目录下的各种元素各个价态的径向电子密度。之后在主功能7里选15，再选1用默认的设置进行Hirshfeld-I电荷计算，经过16轮迭代得到了结果： Final atomic charges: Atom    1(Ti):     3.01789464 Atom    2(O ):    -1.50888824 Atom    3(O ):    -1.50886131 Atom    4(Ti):     3.01787620 Atom    5(O ):    -1.50897240 Atom    6(O ):    -1.50897240 可见Hirshfeld-I电荷的大小不仅显著大于Hirshfeld电荷，比CM5电荷还要大很多。虽然数值很大，但完全在合理范围内，毕竟都没有大过原子的氧化态。 接下来再计算MBIS电荷。在主功能7里选择20，再选1用默认的设置进行MBIS电荷计算，经过19轮迭代得到了结果： Final atomic charges: Atom    1(Ti):     1.73473636 Atom    2(O ):    -0.86733169 Atom    3(O ):    -0.86731253 Atom    4(Ti):     1.73467041 Atom    5(O ):    -0.86734303 Atom    6(O ):    -0.86734303 对于当前例子，MBIS的数值介于CM5和Hirshfeld-I之间。但也有很多反例，MBIS也经常比Hirshfeld-I电荷还大不少，后文有例子。 以上原子电荷的计算顺序是随意的，不要误以为必须按上述顺序计算，你想算哪个就算哪个即可。 4.2 AlN晶体 这一节计算AlN晶体的原子电荷。完全可以按照和上一节相同的做法基于电子密度的cube文件来计算，但本节演示一下基于超胞的molden文件的计算过程。AlN原胞的cif文件是本文文件包里的AlN.cif，将之载入Multiwfn，然后输入以下命令创建CP2K输入文件 cp2k  //创建CP2K输入文件 AlN_443.inp  //产生的输入文件名 -11  //进入几何操作界面 19  //构造超胞 4  //a方向扩胞成原本的4倍 4  //b方向扩胞成原本的4倍 3  //c方向扩胞成原本的3倍 -10  //返回 -2  //要求产生molden文件 0  //产生输入文件 用CP2K运行刚产生的AlN_443.inp，得到AlN_443-MOS-1_0.molden。在此文件开头插入以下内容定义晶胞和各个元素的价电子数 [Cell] 12.44400000     0.00000000     0.00000000 -6.22200000    10.77682012     0.00000000  0.00000000     0.00000000    14.93400000 [Nval] Al 3 N 5 启动Multiwfn，载入AlN_443-MOS-1_0.molden，然后输入 7  //布居分析与原子电荷计算 16  //CM5电荷 [回车]  //使用均匀分布的格点计算 [回车]  //使用默认的0.2 Bohr格点间距。格点间距越小结果越精确但计算越耗时 在i9-13980HX CPU上24核并行，50秒就算完了，结果为                       ---------- CM5 charges ---------- Atom:    1Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302813 Atom:    2N   CM5 charge:   -0.467945  Hirshfeld charge:   -0.302812 Atom:    3Al  CM5 charge:    0.467945  Hirshfeld charge:    0.302812 Atom:    4N   CM5 charge:   -0.467946  Hirshfeld charge:   -0.302814 ...略 Al的电负性很小，N的电负性很大，因此二者的原子电荷理应相差悬殊。当前的Hirshfeld电荷明显偏小。CM5电荷的数值虽然更大，但还是显得有些偏低。 下面再计算Hirshfeld-I电荷。在主功能7里输入 15  //Hirshfeld-I 1  //开始计算 [回车]  //使用默认的0.2 Bohr格点间距 经过37轮收敛，结果如下。可见Hirshfeld-I电荷的数量级比较合理，体现出AlN中成键的极强离子性  Atom    1(Al):     1.80317384  Atom    2(N ):    -1.80317397  Atom    3(Al):     1.80317257  Atom    4(N ):    -1.80317539 ...略 最后再计算下MBIS电荷。在主功能7里输入 20  //MBIS 1  //开始计算 [回车]  //使用默认的0.2 Bohr格点间距 对此体系MBIS收敛很慢，经过190轮才达到收敛，在i9-13980HX上24核并行算了5分钟，在双路7R32 96核并行下算了不到一分钟。结果如下。可见对此体系，MBIS电荷的数值很大，和氧化态几乎正好是一样的。所以，如前所述，若你就是想算出来比较大的原子电荷，可以先考虑MBIS，但也有一些例外，后面会提到。 Atom    1(Al):     2.99599978 Atom    2(N ):    -2.99601488 Atom    3(Al):     2.99600129 Atom    4(N ):    -2.99599935 4.3 MOF-5晶体 MOF-5晶体含有Zn、O、C、H。其cif文件是本文文件包里的MOF-5.cif。由于此体系的晶胞边长约26埃，因此计算电子密度格点数据时不用考虑k点。按照3.1节的做法，基于PBE/DZVP-MOLOPT-SR-GTH级别的电子密度cube文件（第一行改为Zn 12 O 6 C 4 H 1），用Multiwfn计算各种原子电荷，结果如下。此体系中O有两类，分别给出，C则处于较多不同的化学环境，其电荷就不给出了。 Hirshfeld： Zn=0.372 O=-0.332/-0.200 H=0.056 CM5： Zn=0.644 O=-0.584/-0.305 H=0.116 Hirshfeld-I（36轮收敛，2*7R32机子上96核并行4分多钟算完） Zn=1.568 O=-1.917/-0.746 H=0.139 MBIS（39轮收敛，2*7R32机子上96核并行3分多钟算完） Zn=1.546 O=-1.662/-0.782 H=0.190 对这个体系，MBIS和Hirshfeld-I相符得很好，而且Zn的电荷大小很合理。Hirshfeld电荷再次数值最小、CM5比之大一些。 4.4 其它：MoS2、BaTiO3、MgF2、NaCl 注意MBIS电荷并非总是较大。例如《使用Multiwfn结合CP2K计算晶体中原子的氧化态》（//www.umsyar.com/711）里的单层MoS2的例子，在PBE/DZVP-MOLOPT-SR-GTH级别的波函数下，Mo的电荷的计算结果是Hirshfeld=0.330，CM5=0.827，Hirshfeld-I=1.987，MBIS=-0.059，Mulliken=-0.771。对这个体系来说CM5和Hirshfeld-I显得靠谱，MBIS电荷接近0有点说不通，Mulliken电荷明显有误导性。 BaTiO3在PBE/DZVP-MOLOPT-SR-GTH级别下的原子电荷计算结果如下。Hirshfeld还是严重偏小。Hirshfeld-I大得离谱，Ba的原子电荷都超过氧化态+2了。相较之下CM5和MBIS的结果比较靠谱，正好二者算的O的原子电荷差不多，差异主要在Ba和Ti的电荷的相对大小上。 Atom:    1Ba  CM5 charge:    1.103253  Hirshfeld charge:    0.404357 Atom:    2Ti  CM5 charge:    1.200467  Hirshfeld charge:    0.537212 Atom:    3O   CM5 charge:   -0.767963  Hirshfeld charge:   -0.310211 Atom:    4O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693 Atom:    5O   CM5 charge:   -0.767892  Hirshfeld charge:   -0.315693 Hirshfeld-I Atom    1(Ba):     2.37169239 Atom    2(Ti):     3.02878350 Atom    3(O ):    -1.78799416 Atom    4(O ):    -1.80625420 Atom    5(O ):    -1.80625434 MBIS： Atom    1(Ba):     0.57391985 Atom    2(Ti):     1.69416825 Atom    3(O ):    -0.74678425 Atom    4(O ):    -0.76066530 Atom    5(O ):    -0.76066537 再看两个典型离子化合物的情况，都在PBE/DZVP-MOLOPT-SR-GTH级别下计算。对NaCl计算的Na的原子电荷为：Hirshfeld=0.217，CM5=0.443，Hirshfeld-I：1.053，MBIS：1.043，Mulliken：0.596 对MgF2计算的Mg的原子电荷为： Hirshfeld=0.414，CM5=0.828，Hirshfeld-I：2.048，MBIS：1.888，Mulliken：1.439 可见对于强离子性化合物，Hirshfeld-I和MBIS电荷都在氧化态附近，并且由于方法定义的缺陷，往往电荷的大小还轻微大于氧化态的大小（但如果用大核赝势，对Na和Mg只描述最外层的3s价电子，则不会有这个问题）。其它原子电荷的大小关系为Mulliken>CM5>Hirshfeld。 以MgF2为例，这里顺带一提CP2K算的Hirshfeld-I电荷明显不对，千万不要用。在CP2K输入文件的&DFT/&PRINT里加上以下内容要求输出Hirshfeld-I电荷后，CP2K直接给出的Mg的Hirshfeld-I电荷为0.522，和Multiwfn算的相距甚远。 &HIRSHFELD   SHAPE_FUNCTION DENSITY   SELF_CONSISTENT T &END HIRSHFELD 4.5 NaCl板吸附CO的片段电荷 NaCl板吸附CO是《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（//www.umsyar.com/638）里的一个例子体系，这里基于这个模型计算一下原子电荷和CO片段的电荷。CP2K优化这个体系后产生的NaCl_CO-1.restart文件在本文的文件包里提供了，将之载入Multiwfn，然后输入 cp2k  //创建CP2K输入文件 NaCl_CO.inp  //产生的输入文件名 -7  //设置周期性 XY -3  //要求产生cube文件 1  //电子密度 -2  //产生molden文件（用于之后算Mulliken电荷的目的） 0  //产生输入文件 计算完毕后，把NaCl_CO-ELECTRON_DENSITY-1_0.cube文件第一行改为Na 9 Cl 7 C 4 O 6。先用Hirshfeld方法各个原子的电荷和CO的片段电荷。启动Multiwfn，载入此cube文件，然后输入 7  //布居分析与原子电荷计算 -1  //定义片段 109,110  //CO的原子序号 1  //Hirshfeld电荷 结果为 ...略  Atom  109(C ):     0.13622804  Atom  110(O ):    -0.03514734 Fragment charge:    0.10108070 Fragment population:    9.89891930 即CO的Hirshfeld片段电荷为0.101。 之后再选择CM5方法，结果为 ...略  Atom:  109C   CM5 charge:    0.147306  Hirshfeld charge:    0.136228  Atom:  110O   CM5 charge:   -0.081972  Hirshfeld charge:   -0.035147  Summing up all CM5 charges:     0.00510030  Fragment charge:    0.06533361 之后再选择Hirshfeld-I方法，结果为  Atom  109(C ):     0.09790054  Atom  110(O ):    -0.08744136  Fragment charge:    0.01045918 之后再选择MBIS方法，结果为 ...略  Atom  109(C ):     0.19642171  Atom  110(O ):    -0.15937966  Fragment charge:    0.03704205 再来算Mulliken电荷。用Multiwfn载入加入了恰当的[Cell]和[Nval]信息的CP2K计算后产生的NaCl_CO-MOS-1_0.molden文件，输入 7  //布居分析与原子电荷计算 -1  //定义片段 109,110  //CO的原子序号 5  //Mulliken电荷 1  //在屏幕上输出结果 结果为 ...略 Atom   109(C )    Population:  3.60671601    Net charge:  0.39328399 Atom   110(O )    Population:  6.23484007    Net charge: -0.23484007 Total net charge:   -0.00000086 Fragment charge:    0.15844392 可见不同方法计算的结果在定量上有一定差异，但都有共性，即CO整体带少量正电，即电子从CO往NaCl板上转移了一些。并且C和O分别带正电和负电，且原子电荷的绝对值都是C大于O。 5 总结 本文简要介绍了四种知名的基于模糊式原子空间划分计算原子电荷的方法的特点，包括Hirshfeld、Hirshfeld-I、CM5和MBIS，并介绍了它们在Multiwfn中的计算方法，并且给了很多具体例子，以使得读者能了解如何使用Multiwfn结合CP2K或其它第一性原理程序很容易地计算它们。同时还说明了怎么直接计算片段电荷，这是定量讨论体系中两部分间电子转移的最严格的方式。推荐将这些原子电荷应用于实际问题的研究当中以考察固体与表面体系的电荷分布情况。记得届时应恰当引用Multiwfn启动时提示的Multiwfn原文以及相应的原子电荷的原文。 通过本文的一些体系的对比可以明显看到，Hirshfeld电荷原理上最简单，但明显缺点是整体严重偏小，CM5通常比它大不少，而且电荷的稳定性不错。虽然结合化学直觉来看，CM5电荷的绝对大小仍然往往能偏小，但至少也不会明显离谱，总比Hirshfeld电荷更建议使用。Hirshfeld-I和MBIS电荷普遍比较大，谁更大不一定，至少几乎总比CM5还大不少，从绝对大小来看往往比CM5更合理，但稳健性比CM5弱一些，对极个别体系可能会表现得不符合常识。对于离子性特征很强的化合物，Hirshfeld-I和MBIS电荷往往接近氧化态，由于方法并不完美，有时候原子电荷甚至比氧化态还要大一点点。当然，没有任何原子电荷计算方法是绝对完美的，本文介绍的这些原子电荷在实际当中都可以用，可以根据它们的特点和实际结果决定选用，在横向对比时方法必须统一。 顺带一提，对于计算固体表面或者多孔物质的原子电荷以用于基于力场的动力学模拟的目的，我目前最推荐CP2K算REPEAT电荷，北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里讲“电子结构的分析”的部分做了专门的介绍。但REPEAT电荷需要在真空区定义拟合点，因而无法用于诸如本文所示的AlN、TiO2等致密的固体，而且对于远离真空区的原子也没法拟合出质量靠谱的原子电荷。所以REPEAT电荷有特定用处，但不是像本文介绍的那些原子电荷一样属于普适性的原子电荷计算方法。 //www.umsyar.com/711 2024-05-23T16:05:00+08:00 使用Multiwfn结合CP2K计算晶体中原子的氧化态 Using Multiwfn in combination with CP2K to calculate oxidation state for atoms in crystals 文/Sobereva@北京科音  2024-May-23 1 前言 波函数分析程序Multiwfn可以基于电子波函数计算原子和片段的氧化态，这在笔者之前写的《使用Multiwfn通过LOBA方法计算氧化态》（//www.umsyar.com/362）一文中已做了详细介绍，没读过者一定要先阅读此文，否则无法理解下文的内容。此文介绍的做法已经被广泛用来准确地考察分子中的氧化态。如今Multiwfn已经支持将这种方法用于CP2K计算的周期性体系的波函数，从而考察晶体中原子的氧化态，这非常有实际价值，具体做法将在本文介绍。注意氧化态和原子电荷完全是两码事，周期性体系的原子电荷的计算看《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（//www.umsyar.com/712）。 Multiwfn可以在官网//www.umsyar.com/multiwfn免费下载，本文的读者必须使用2024年5月23日及以后更新的版本，否则没有本文介绍的特性。对Multiwfn不了解者看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。 2 用法 在Multiwfn中对周期性体系做LOBA或mLOBA分析氧化态的具体操作流程和对孤立体系做别无二致，都是先用主功能19产生定域化轨道（介绍见《Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比》//www.umsyar.com/380），然后进入主功能8（轨道成分分析），选择LOBA/mLOBA分析，之后输入一个LOBA方法判断定域化轨道电子归属的阈值，或者输入m让Multiwfn按照mLOBA的方式判断归属即可。定域化轨道的电子归属是什么意思，在前述的《使用Multiwfn通过LOBA方法计算氧化态》里已经专门说过了，这里不再累述。之后各个原子的氧化态就会立马输出出来。用户也可以自定义片段来得到片段的氧化态。 对周期性体系做LOBA/mLOBA分析的输入文件必须是CP2K产生的molden文件，并且需要自己在里面加入盒子信息。具体方法在《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）里有详细介绍。由于考虑k点时不能产生molden文件，因此原本晶胞较小的话需要扩胞成足够大的超胞，从而对其计算时只需要考虑gamma点。Multiwfn不支持Quantum ESPRESSO、VASP等其它第一性原理程序产生的波函数，也没有任何支持的必要，因为对于绝大多数情况，用Multiwfn创建CP2K单点输入文件并用CP2K计算是非常简单且快速的事，可参考《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）。 对于当前目的，CP2K计算时用赝势基组也行，用全电子基组也行，没有什么特别的讲究。CP2K计算的时候不要开smearing，否则轨道占据数不再精确为整数、不再是标准的单行列式波函数，此时也就没法做轨道定域化。如果有特殊原因非开smearing不可，可以在Multiwfn载入molden文件后进入主功能6，选择选项38，此时Multiwfn会按照轨道能量由低到高重排电子确定占据数，轨道占据数就都为整数了，之后退回到主菜单，再照常按前述步骤做LOBA/mLOBA分析即可。 3 例子 下面给出一系列通过Multiwfn使用LOBA/mLOBA计算固体的氧化态的例子。CP2K的输入输出文件在//www.umsyar.com/attach/711/file.zip里都给了。只有3.1节例子对应的molden文件我在里面直接给了，如果所有例子的molden文件都给的话压缩包就太大了。CP2K输入文件里的各种设置我这里就不细说了，在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里有非常全面系统的讲解。 3.1 单层氮化硼 氮化硼（BN）是知名的二维材料，这里完整介绍单层BN的氧化态计算的流程。本文文件包里的BN_bulk.cif是氮化硼晶体的结构文件，用GaussView打开它，可看到晶胞里有两层BN。删除其中的一层，然后保存为BN_ml.gjf。之后启动Multiwfn，载入BN_ml.gjf，输入 cp2k   //创建CP2K输入文件 [回车]  //产生的输入文件用默认的名字BN_ml.inp -2  //要求导出molden文件 -11  //进入几何操作界面 19  //构造超胞 6  //第1个平行于层的方向扩胞成原先的6倍 6  //第2个平行于层的方向扩胞成原先的6倍 1   //垂直于层的方向不扩胞 -10  //返回 -7  //设置周期性 XY 0  //产生输入文件（计算级别为默认的PBE/DZVP-MOLOPT-SR-GTH） 用CP2K运行Multiwfn产生的BN_ml.inp，我这里用双路7R32 96核机子4秒钟就算完了。现在当前目录下有了BN_ml-MOS-1_0.molden，按照《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）所述，把以下信息插入进去，以定义晶胞信息和原子价电子数信息。改好的文件已经在本文的文件包里提供了，是BN_ml-MOS-1_0.molden。  [Cell] 14.98944000     0.00000000     0.00000000 -7.49472000    12.98123580     0.00000000  0.00000000     0.00000000     8.00000000  [Nval] B 3 N 5 启动Multiwfn，载入BN_ml-MOS-1_0.molden，然后输入 19  //轨道定域化 1  //只定域化占据轨道，用默认的Pipek-Mezey方法 8  //轨道成分分析 100  //LOBA/mLOBA分析 输入一个阈值，比如50，马上看到各个原子的以LOBA方法计算的氧化态，如下所示。可见B的氧化态为3、N为-3，整个体系所有原子氧化态加和为0，非常符合常识。这里用的阈值50%对于成键方式比较典型的情况是适合的，它的意思是如果有某原子对某个定域化轨道贡献超过50%（Multiwfn以Hirshfeld方法计算贡献值），则这个轨道的电子就完全归属于这个原子。  Oxidation state of atom   1(B ) :  3  Oxidation state of atom   2(N ) : -3  Oxidation state of atom   3(B ) :  3  Oxidation state of atom   4(N ) : -3  Oxidation state of atom   5(B ) :  3  Oxidation state of atom   6(N ) : -3 ...略  The sum of oxidation states:     0 我提出的modified LOBA (mLOBA)方法是对LOBA的重要改进，它将每个定域化轨道的电子完全归属到对它贡献最大的那个原子上，可以确保所有原子的氧化态加和总为0，而且也没有阈值选择的任意性（有的时候LOBA的结果对阈值敏感），比LOBA在原理上明显更好。在当前界面里输入m即得到mLOBA方法算的氧化态的结果，也是B和N的氧化态分别为3和-3。 3.2 BaTiO3 BaTiO3晶体的cif文件是本文文件包里的BaTiO3.cif。用Multiwfn创建其4*4*4超胞在PBE/DZVP-MOLOPT-SR-GTH下做单点的CP2K计算的输入文件。当前的超胞在每个方向有16埃，足够大了。CP2K计算后得到molden文件，在里面加入恰当的[Cell]和[Nval]字段后，启动Multiwfn并将之载入，然后输入 19  //轨道定域化 -9  //设置产生定域化轨道后自动计算轨道成份的方法 0  //不输出轨道成份。这些信息不是当前感兴趣的，当前体系又大，不输出可以节约不少时间 1  //只定域化占据轨道，用默认的Pipek-Mezey方法 由于当前体系大、轨道数很多，因此轨道定域化过程耗时较高。在EPYC 7R32机子上花了约8分钟。之后输入 8  //轨道成分分析 100  //LOBA/mLOBA分析 m  //显示mLOBA方法计算的原子氧化态 马上看到以下结果。可见Ba、Ti、O的氧化态都和常识完全一致  Oxidation state of atom   1(Ba) :  2  Oxidation state of atom   2(Ti) :  4  Oxidation state of atom   3(O ) : -2  Oxidation state of atom   4(O ) : -2  Oxidation state of atom   5(O ) : -2 ...略 3.3 MOF-5 MOF-5是知名的金属有机框架化合物，介绍见https://en.wikipedia.org/wiki/MOF-5。其中含有Zn，这里算算Zn的氧化态是多少。MOF-5晶体的cif文件在本文的文件包里提供了，是立方晶胞，边长25.86埃。由于边长已经大到可以只考虑gamma点做计算，因此CP2K计算时不需要扩胞，直接对原胞算单点得到molden文件即可。相应的CP2K输入文件已提供在了本文的文件包里。氧化态计算过程同前，7R32机子上几分钟可以算完整个过程。氧化态结果为Zn=2、O=-2、H=1、C=2/0/2。可见Zn、O、H的氧化态都完全符合预期。当前体系中C的氧化态有三种，和所处的化学环境差异有关。 3.4 单层MoS2 MoS2是知名的二维材料体系，它的体相结构文件是本文文件包里的MoS2_bulk.cif。完全效仿3.1节氮化硼的例子，取单层并基于5*5*1的超胞模型做它的氧化态的计算。如果以LOBA方法计算氧化态并输入50作为阈值，结果为  Oxidation state of atom   1(Mo) :  6  Oxidation state of atom   2(S ) : -2  Oxidation state of atom   3(Mo) :  6  Oxidation state of atom   4(S ) : -2  Oxidation state of atom   5(S ) : -2  Oxidation state of atom   6(S ) : -2  Oxidation state of atom   7(Mo) :  6 ...略  The sum of oxidation states:   50 虽然Mo的价电子有6个，看似Mo的氧化态为6也不是不可能，但实际上明显是错的，因为当前所有原子的氧化态加和为100。考虑到S的氧化态不可能比-2更负，明显当前的结果高估了Mo的氧化态。即便尝试各种阈值，也始终得不到能令LOBA方法算的氧化态加和为0的情况。因此LOBA对此体系完全失败。 输入m使用mLOBA计算氧化态，结果为  Oxidation state of atom   1(Mo) :  4  Oxidation state of atom   2(S ) : -2  Oxidation state of atom   3(S ) : -2  Oxidation state of atom   4(Mo) :  4  Oxidation state of atom   5(S ) : -2  Oxidation state of atom   6(S ) : -2  Oxidation state of atom   7(Mo) :  2  Oxidation state of atom   8(S ) : -2  Oxidation state of atom   9(S ) : -2 ...略  The sum of oxidation states:     0 可见mLOBA方法算的所有原子氧化态加和为期望的0，但是Mo有的氧化态是2有的是4，和当前体系中所有Mo等价的事实不符，原因后面会说。解决方法是输入-1定义片段，然后输入所有Mo对应的序号。序号的快速查询方法是在Multiwfn主功能0的菜单里选择Tools - Get atom indices of a given element，然后输入Mo，点OK，就会返回下面的内容，将之粘贴到Multiwfn定义片段的窗口里即可。 1,4,7,10,13,16,19,22,25,28,31,34,37,40,43,46,49,52,55,58,61,64,67,70,73 片段定义好后，再次输入m，可以看到Oxidation state of the fragment: 100，即这个片段总的氧化态为100。由于当前体系里有25个完全等价的Mo，因此Mo的氧化态为100/25=4，是完全正确的值。 为什么mLOBA算出来的明明等价的Mo会有不同的氧化态？这可以通过考察定域化轨道的特征结合mLOBA的原理予以了解。默认情况下，Multiwfn的主功能19做完轨道定域化后会自动输出各个定域化轨道的成份，可以看到除了单中心、双中心的定域化轨道外，还有大量三中心的：  More delocalized LMOs: (Three largest contributions are printed)   301:   64(Mo) 24.8%   46(Mo) 24.8%   61(Mo) 24.8%   302:   46(Mo) 24.8%   49(Mo) 24.8%   31(Mo) 24.8%   303:   16(Mo) 24.8%   34(Mo) 24.8%   31(Mo) 24.8% ...略 从轨道成份上明显可知这些三中心轨道都是等价的。轨道定域化做完后按照《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）说的，在Multiwfn的主功能0里随便看一个三中心轨道，如下所示 可见这个轨道体现的是三个Mo之间的三中心作用。按照mLOBA的思想，这个轨道上的两个电子会完全归属到对它贡献最大的原子上，但由于三个Mo的贡献本质上相同，实际的贡献量由于数值计算误差会有极其轻微的差异，所以这个轨道上的电子归属到哪个Mo是有随机性的，这是为什么Mo的氧化态有的为2有的为4。可见，只要把轨道定域化和mLOBA的思想都理解了，遇到可疑的结果自然就能通过具体分析搞明白是怎么回事、知道怎么恰当解决。 4 总结 本文详细介绍了怎么用Multiwfn基于CP2K对固体做DFT计算得到的波函数通过LOBA和mLOBA方法得到原子的氧化态，给了很多例子。可见本文的方法操作简单，结果非常合理，这对于讨论原子在固体中的状态、对体系进行分类很有价值。作为练习，读者可以对NaCl、MgO、AlN、SiO2晶体也都算算氧化态。值得一提的，如本文的例子所体现的，如有了我后来提出的mLOBA方法，原本的LOBA就没什么使用价值了，我建议总是默认用mLOBA。 使用本文的方法计算在发表结果时务必记得按照Multiwfn程序启动时的提示对Multiwfn进行正确的引用。 //www.umsyar.com/706 2024-05-05T14:31:00+08:00 2024年计算化学公社举办的计算化学程序和DFT泛函的流行程度投票结果 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更是快跌没了。 //www.umsyar.com/690 2023-11-07T17:53:00+08:00 CP2K做杂化泛函计算的关键要点和简单例子 Key points and simple examples of hybrid functional calculations with CP2K 文/Sobereva@北京科音  First release: 2023-Nov-7   Last update: 2024-Mar-9 0 前言 对周期性体系，杂化泛函的计算耗时众所周知远高于纯泛函。著名的第一性原理程序CP2K不光纯泛函的计算超级快，对大体系做杂化泛函计算也非常快（尽管还是远比纯泛函贵），远远快于Quantum ESPRESSO (QE)、VASP等纯粹基于平面波的程序。CP2K结合像样的基组做杂化泛函的单点计算，在一般双路服务器上算到上千原子都没问题。然而，使用CP2K做高效的杂化泛函计算有很多非常关键性的要点必须知道，远远不是像Gaussian程序那样光写个杂化泛函名字就了事的。然而，笔者经常在网上回答CP2K问题时看到很多人完全不具备这些最基本常识，胡用瞎用CP2K做杂化泛函计算，导致耗时巨高、始终出不来结果、任务崩溃、SCF完全不收敛等各种问题。因此笔者觉得很有必要写个小文将CP2K做杂化泛函计算的最关键常识和要点挨个提一下。Multiwfn程序（//www.umsyar.com/multiwfn）具有极其便利的创建CP2K输入文件的功能，在本文的最后会还用Multiwfn产生C60晶体的CP2K杂化泛函计算的输入文件作为实际例子。 本文的内容至少对CP2K 2024.1，以及2023-Nov-7及以后更新的Multiwfn是适用的。 笔者在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）的“能量的计算及相关问题”部分里讲杂化泛函计算远比本文全面、深入得多，并列举诸多经验技巧和例子，还详细介绍HSE06、CAM-B3LYP、wB97X等各种范围分离形式的杂化泛函以及revDSD-PBEP86-D3(BJ)等各种双杂化泛函的用法和要点，使用杂化泛函结合k点计算的要点，欢迎想学习更多知识者关注和参加。 1 CP2K杂化泛函计算的特性 --------- 2024-Feb-18更新：以下之前写的关于k点的说法已经过时。从CP2K 2024.1版开始支持了杂化泛函计算时考虑k点，称为RI-HFXk算法，并且可以无缝结合ADMM降低耗时。此时需要以RI近似计算HF交换部分。不仅能算能量还能做优化（包括变胞）和能带计算。从第3届的“北京科音CP2K第一性原理计算培训班”开始已经加入了RI-HFXk的全面深入的介绍和具体例子，参考http://bbs.keinsci.com/thread-43683-1-1.html。 首先要知道CP2K做杂化泛函计算的一个关键不足是不支持考虑k点，因此如果被计算的是周期性体系，晶胞必须大到只考虑gamma点就足够。如果你要计算的体系的晶胞原本比较小，显然必须先扩胞到足够大才行。如果不知道扩到多大够用，严格来说需要做你感兴趣的属性的计算结果随扩胞倍数的收敛性测试。虽然CP2K能够对含有非常多原子的超胞做得动杂化泛函计算，但耗时会明显高于QE等程序对原胞在考虑k点的情况下做杂化泛函计算。 最能充分体现CP2K杂化泛函计算效率优势的情况是算一些较大分子晶体、MOF、COF等一些大晶胞体系，以及做第一性原理动力学等情况。对小晶胞体系如上所述经过扩胞后也能算，只不过不体现CP2K的长处。 CP2K显然没法用杂化泛函算能带，因为不支持k点，这种目的只能用QE等程序代替。必须用CP2K算的话，应当用CP2K支持的纯泛函里算带隙相对比较好的比如HLE17。 --------- 由于杂化泛函远比普通泛函昂贵，尤其是用于几何优化、NEB、动力学等涉及到结构变化的任务。因此如果你算的问题用纯泛函就足以描述得不错、用杂化泛函不会带来明显好处的话，显然应该用纯泛函。记住，非必要甭考虑杂化泛函，纯泛函的耗时仅有杂化泛函的一个微不足道的零头。对几何优化、振动分析、NEB等高耗时任务用纯泛函，算势垒、反应能的时候用更好的杂化泛函是可以接受且常见的策略。 周期性体系的杂化泛函计算不仅耗时远高于纯泛函计算，对内存要求更是远远高于纯泛函计算。经常用CP2K做杂化泛函任务的服务器的内存容量绝对是多多益善，选择服务器的时候要注意这一点。也别天真地试图拿个人计算机将就着跑CP2K的杂化泛函计算。 CP2K做杂化泛函计算的前提是编译时必须带着libint电子积分库，否则没法算杂化泛函中HF交换项中涉及的双电子积分。CP2K的编译方法见//www.umsyar.com/586。 杂化泛函用的赝势基组和赝势就用常见的给PBE等纯泛函优化的就可以。从CP2K 9.1开始，data目录下的BASIS_MOLOPT_UZH和POTENTIAL_UZH文件里还分别给出了专门给PBE0优化的赝势基组和赝势，如C的DZVP-MOLOPT-PBE0-GTH-q4和GTH-PBE0-q4，原理上来说它们用于杂化泛函计算会更理想一些（但实际未必有什么优势）。 2 库仑截断 杂化泛函计算之所以比纯泛函昂贵得多，在于计算杂化泛函的HF交换项时要计算大量的双电子积分（two-electron integral, ERI），尤其是对于周期性体系来说，不光要计算晶胞内基函数之间的ERI，还要计算中心晶胞与周围镜像晶胞基函数之间的ERI。再加上由于只能算gamma点而晶胞必然较大，原子数、基函数通常非常多，无疑要算的ERI是巨量的。 由于周期边界条件，原理上来说要计算的ERI是无穷的，显然这会导致计算无法进行。因此实际计算时必须使用库仑截断（Coulomb truncation），也就是将ERI中描述电子库仑相互作用的1/r算符截断在某个距离（截断半径），超过截断半径时1/r视为0，相应地ERI就不需要计算了，这个策略使得ERI的计算数目是有限的。如果你计算的是非周期性体系，由于本身ERI数目就是有限的，因此虽然也可以用库仑截断，但这就不是必须的了。 杂化泛函的库仑截断的半径设置显然既影响耗时也影响结果。截断半径越小，要计算的ERI越少，耗时就越低，但结果越偏离泛函原本的精确结果。显然不能为了一味地节约时间而把截断半径设得太小，否则结果会很垃圾。6埃通常是耗时和精度的较好权衡，设得明显更小会令精度有明显损失，设得更大并不会令精度有明显提升却会增加耗时不少。如果对精度要求很高也可以用8埃，设得更大就完全没必要了。 要注意像6、8埃这种程度的截断半径其实还没有大到令绝对能量充分收敛，但能量在相互求差时，由截断半径带来的系统性误差能很大程度抵消。可以认为用不同的截断半径相当于用不同的理论方法，直接对比能量的情况必须保持截断半径的设置相统一！ 值得一提的是，不一定截断半径小结果就一定差，有可能截断半径造成的对ERI忽略的误差和泛函本身的误差相抵消导致恰好对某些体系某些问题的计算比用更大截断半径时结果还更好。 截断半径应当小于晶胞最短尺寸的一半以免导致得到缺乏物理意义的结果。这里所谓的最短尺寸，是(001)晶面间距、(010)晶面间距、(100)晶面间距三者中的那个最小值。也因此为了能够使用6埃截断半径，晶胞最短尺寸应当大于12埃，如果晶胞不够大则应当扩胞（即便抛开库仑截断问题不谈，单从只能计算gamma点这点来说，晶胞也不能过小）。 Multiwfn创建的杂化泛函输入文件中，若周期性设成了NONE，则不使用库仑截断，若是周期性，Multiwfn会自动加上库仑截断的设置。如果晶胞最短尺寸够大，Multiwfn会自动把截断半径设为6埃，如果不够大，则会根据当前晶胞最短尺寸自动设成能设的最大半径。 以下是库仑截断的设置方式，加入到&DFT/&XC里面。POTENTIAL_TYPE TRUNCATED代表用库仑截断的1/r算符，CUTOFF_RADIUS设置截断半径（埃）。T_C_G_DATA t_c_g.dat是默认的可以不写，这指定的是库仑截断用的gamma函数的数据文件t_c_g.dat，在CP2K的data目录下。 &HF   &INTERACTION_POTENTIAL     POTENTIAL_TYPE TRUNCATED     CUTOFF_RADIUS 6.0     T_C_G_DATA t_c_g.dat   &END INTERACTION_POTENTIAL &END HF 3 积分屏蔽 积分屏蔽（integral screening）是指尽可能在不明显牺牲精度的前提下减少要计算的ERI的量，这对于CP2K杂化泛函计算的耗时和内存使用量有至关重要的影响！上一节的库仑截断本身就是一种积分屏蔽策略，CP2K里还有另外两种非常重要的积分屏蔽方法，也是普遍应用在Gaussian等 程序中的： (1)利用Schwarz不等式做积分屏蔽：利用Schwarz不等式关系，可以根据数目较少的二指标ERI估计出数目甚巨的三、四指标ERI中哪些数值肯定小于阈值从而可以忽略掉，这可以大幅减少ERI要算的数目（原本总共要算的ERI中占比最大的就是四个基函数序号不同的那些ERI）。阈值通过&DFT/&XC/&HF/&SCREENING中EPS_SCHWARZ设置，此值越大，被忽略的ERI数目越多、计算精度越差。默认的1E-10虽然足够保证精度，但此时太贵，没必要。通常设为1E-6就可以了，是精度和耗时的较好权衡，想要精度更好点用1E-8。 (2)结合密度矩阵做积分屏蔽：ERI对能量的贡献取决于它与特定密度矩阵元的乘积。因此哪怕ERI不是很小，但与之相乘的密度矩阵元很小，那么这样的ERI也完全不需要计算。CP2K定义了特定方法快速估计一个ERI与相应密度矩阵元乘积值的上限，如果小于上述的EPS_SCHWARZ，则此ERI就会忽略掉。这种依赖于密度矩阵的积分屏蔽默认是关闭的，如果开启则把&DFT/&XC/&HF/&SCREENING里的SCREEN_ON_INITIAL_P设为T。由于这种积分屏蔽能令耗时巨幅降低，因此我强烈建议总是开启，也因此Multiwfn产生的杂化泛函的输入文件里都自动把SCREEN_ON_INITIAL_P设成T。 基于密度矩阵做积分屏蔽有一个极其关键的要点是初猜的密度矩阵质量不能太差！然而很多人对此一无所知！默认情况下初始的密度矩阵是CP2K自动做初猜产生的，由于它和SCF收敛时的密度矩阵相差很大、质量很糙，因此此时积分屏蔽通常做得极烂，不仅可以忽略的ERI可能没忽略，还会把一些重要的ERI给忽略掉，由此可能导致KS矩阵质量很烂、SCF极难收敛。而且忽略哪些ERI是在SCF过程一开始决定好的，即便之后随着SCF迭代密度矩阵质量不断变好，也不会由此使得积分屏蔽逐步变得合理，这导致即便最后SCF收敛了，得到的能量的精度也可能很烂。因此当SCREEN_ON_INITIAL_P设T时，必须读取一个基本合理的波函数当初猜。通常是在杂化泛函计算前先用相同的基组做一个便宜的纯泛函计算，比如PBE0计算前先用PBE做一个单点计算，得到记录了PBE收敛的波函数信息的wfn文件，然后PBE0计算时用SCF_GUESS RESTART，并把WFN_RESTART_FILE_NAME指定成那个wfn文件，以读取其作为初猜。这不仅使得基于密度矩阵的积分屏蔽可以合理进行，还有一个额外的好处是由于初猜波函数质量较好（至少比默认自动产生的好得多），使得杂化泛函计算时达到SCF收敛所需的迭代次数会比较少，这进一步节约了时间（注意即便抛开ERI计算耗时不谈，杂化泛函在构造KS矩阵时花的时间也明显高于纯泛函，因为涉及到密度矩阵元与ERI的大量的乘加操作）。 一般的 程序虽然也基于密度矩阵做积分屏蔽，但不是非得像CP2K这样得读取一个像样的初猜波函数，这是因为它们会根据当前最新的密度矩阵判断这一轮SCF用的ERI哪些是能被忽略的。CP2K的做法和它们不同，主要是因为CP2K倾向于让用户做incore形式的SCF而非一般量化程序一般用的direct形式的SCF，见后文。 当SCREEN_ON_INITIAL_P设T时，对于杂化泛函的几何优化、分子动力学等任务的续算，要注意光有.restart文件还不行，还必须提供上一步的.wfn文件当初猜波函数，否则也是会由于基于自动初猜的密度矩阵做积分屏蔽非常恶心而令任务无法正常跑下去。 如果你嫌麻烦而不想先做纯泛函计算产生收敛的波函数再做杂化泛函计算，那就不要写SCREEN_ON_INITIAL_P T。对于小体系、小基组特别是非周期性的计算，不基于密度矩阵做积分屏蔽时耗时也不高。而对于计算量较大的杂化泛函计算，则强烈建议不要偷这个懒而多花巨量不必要的时间。 4 ADMM CP2K还支持其开发者独创的ADMM (Auxiliary Density Matrix Methods)方法减少要算的ERI数目从而进一步显著加快杂化泛函的计算，其原理见以下北京科音CP2K第一性原理计算培训班的ppt。 ppt里提到的主基组就是指平时实际用的基组，诸如DZVP-MOLOPT-SR-GTH。对于不同的主基组，有一些专门构造的与之匹配的ADMM辅助基组，比如MOLOPT系列基组适合搭配的是CP2K的data目录下的BASIS_ADMM_UZH文件里的那些辅助基组，从小到大依次是admm-dz、admm-dzp、admm-tzp、admm-tz2p。admm-dzp是2-zeta结合极化函数的辅助基组，通常算是保底质量。DZVP-MOLOPT-SR-GTH主基组适合搭配admm-dzp，更好的主基组如TZVP-MOLOPT-GTH结合admm-dzp也可以接受，若结合admm-tzp精度会更好些，这都比起不用ADMM能巨幅节约计算时间和内存。辅助基组的大小与主基组越接近，ADMM给杂化泛函HF交换项带来的误差越小，但需要计算越多ERI因而越耗时。如果辅助基组和主基组设为相同的，那么ADMM不会造成任何误差，相应地也不会节约任何耗时（反倒还会增加一些耗时）。 data目录下的BASIS_ADMM文件中的FIT3、pFIT3等辅助基组都比较过时了，对元素涵盖也远不如BASIS_ADMM_UZH里的完整，因此如今不再建议使用（它们出现在不少官方例子文件和文献中，如今不要效仿）。目前版本的Multiwfn产生利用ADMM的杂化泛函的输入文件时默认设的是admm-dzp。 注意前述这些辅助基组都是对GTH赝势基组用的，千万别做全电子计算时也用它们。要想全电子计算时用ADMM，可以用(aug-)pcseg系列基组作为主基组，在《在线基组和赝势数据库一览》（//www.umsyar.com/309）中介绍的BSE基组数据库中有和它匹配的(aug-)admm辅助基组，诸如pcseg-2适合搭配admm-2。pcseg是一类构造得很好的对DFT计算非常划算的基组，《谈谈 中基组的选择》（//www.umsyar.com/336）里简要提及了，在北京科音中级 培训班（http://www.keinsci.com/workshop/KBQC_content.html）里讲基组的地方做了具体介绍。 用MOLOPT系列基组做周期性体系的杂化泛函计算几乎是一定要结合ADMM的，否则极难算得动！因为MOLOPT是完全广义收缩基组，而CP2K利用的计算源高斯函数（PGTF）之间ERI的libint库又没有专门考虑广义收缩的情况，因此造成等效的PGTF数目甚巨，ERI的数目又与PGTF数目的四次方近似成正比，数目显然多得恐怖。如果你用片段收缩基组，比如def2-SVP、6-31G**、pcseg-1、pob-DZVP-rev2之类去做杂化泛函计算，由于需要考虑的PGTF数目远比档次与之相仿佛的DZVP-MOLOPT-SR-GTH要少，因此ADMM方法不是非用不可，但不用的话耗时会显著高于DZVP-MOLOPT-SR-GTH搭配admm-dzp的情况。 不能有的原子用ADMM加速计算而有的原子不用。因此如果你用混合基组，其中有的原子用的主基组有ADMM辅助基组，有的没有，你又想用ADMM，那么没有ADMM辅助基组的那些原子的辅助基组就只能设定成其主基组。虽然这些原子没法被ADMM技术加速，但至少其它那部分原子能被加速。 前面提到的各种节约耗时的策略是彼此完全兼容的，同时利用可以最大程度节约时间。即周期性体系杂化泛函计算时一般要用库仑截断、利用Schwarz不等式做积分屏蔽、结合密度矩阵做积分屏蔽、开启ADMM。所有这些设定都会对能量产生影响，因此写文章的时候最好注明，横向对比时应统一。 启用ADMM需要在&KIND里加上BASIS_SET AUX_FIT [辅助基组名]，同时增加一行BASIS_SET_FILE_NAME [辅助基组文件名]。并且在&DFT里加入以下内容 &AUXILIARY_DENSITY_MATRIX_METHOD &END AUXILIARY_DENSITY_MATRIX_METHOD 如果用的是对角化而非OT，还需要在以上字段里加入ADMM_PURIFICATION_METHOD NONE要求不做纯化。 5 ERI的储存和内存使用量的控制 SCF有两种常见形式，incore SCF是指SCF计算一开始将所有要用到的ERI全都算出来并储存在内存里，之后SCF每一轮在构造KS矩阵时会直接从内存里读取ERI而不再重新计算。direct SCF是指SCF每一轮都重新计算要用到的ERI，这也叫on-the-fly方式计算ERI。显然，incore比direct总耗时低得多，但代价是机子的必须内存非常大，因为要存的ERI通常非常多。在CP2K的杂化泛函计算中，如果你给的内存足够存得下所有ERI，那么只有SCF第一轮的时候由于要计算所有ERI所以耗时很高，而之后SCF每一轮的耗时就很低了。如果给的内存不够储存所有ERI，那么能在内存里存多少ERI就存多少，存的那部分在之后每一轮SCF中会直接读取，存不下的ERI在每一轮SCF中会重新算，相当于介于incore和direct之间的情况。显然，做大体系杂化泛函计算时应当尽量在大内存机子上做、分配给CP2K尽可能多的内存，以让尽可能多的ERI能存到内存中、尽可能减少每轮SCF过程中要on-the-fly计算的ERI的量来节约耗时。 常有人问为什么他做杂化泛函计算时老是卡在SCF刚开始的位置不动，这是因为SCF第一轮肯定是要计算所有需要考虑的ERI，耗时显然往往很高，对大体系、大基组的情况花一、两个小时都是常事，看上去好像一直停在这里似的。如果你给的内存较少而导致大部分ERI都得on-the-fly方式去算，那么之后SCF每一轮都得花将近这么多时间，整个任务跑完可能得半天甚至一天时间。如果PRINT_LEVEL设的是MEDIUM，在SCF第一轮算完后就会输出Number of sph. ERI's calculated on the fly和Number of sph. ERI's stored in-core，前者是之后每一轮SCF都需要重新算的ERI数，后者是已存到内存中而之后不需要SCF每一轮重新算的ERI数。显然前者相对于后者越小，之后每一轮SCF比第一轮时间少得越多。 CP2K用户通常默认用popt版本，每个MPI进程在计算HF交换项时最大内存使用量由&DFT/&XC/&HF/&MEMORY里的MAX_MEMORY控制，单位为MB，扣除掉储存零碎数据的占用量外，其余都用来储存ERI。MAX_MEMORY与MPI并行进程数的乘积必须小于当前机子空余物理内存量，显然应当尽量给大以尽可能减少on-the-fly方式算的ERI量，但也不要太顶着头分配，因为进入HF交换项计算模块之前CP2K还会占一定量的内存，故要有适当余量。在SCF第一轮不断往内存里写入ERI、内存占用量不断增大的过程中，若空余内存已用完时CP2K还在继续往里存ERI，CP2K就会马上崩溃，甚至还可能导致计算机暂时失去响应。所以MAX_MEMORY应当设得恰到好处。 如果空余物理内存很有限，为了能让MAX_MEMORY能设得比较大，可以减小MPI进程数。为了不因此浪费CPU运算能力，此时可以用psmp版，此时每个MPI进程下属的OpenMP线程都可以共享这个MPI进程储存的ERI。 6 杂化泛函的定义方式 这一节简要说一下杂化泛函怎么正确在输入文件里定义。对于非周期性体系，比如用PBE0，直接在&DFT里写上以下内容即可。此时不用库仑截断而用完整的1/r算符，EPS_SCHWARZ为默认的1E-10、SCREEN_ON_INITIAL_P为默认的F。 &XC   &XC_FUNCTIONAL PBE0   &END XC_FUNCTIONAL &END XC 显然对周期性体系绝对不能简单写成上面那样，至少也得自己去指定为使用库仑截断。结合前面的讲解，周期性体系PBE0计算的泛函定义部分一般写为下面这样。PBE0泛函由75%的PBE交换项+25%的HF交换项+100%的PBE相关项构成，所以需要恰当地定义&PBE里面的参数，并在&HF里用FRACTION指定HF成份。 &XC   &XC_FUNCTIONAL     &PBE       SCALE_X 0.75       SCALE_C 1.0     &END PBE   &END XC_FUNCTIONAL   &HF     FRACTION 0.25     &SCREENING       EPS_SCHWARZ 1E-6       SCREEN_ON_INITIAL_P T     &END SCREENING     &INTERACTION_POTENTIAL       POTENTIAL_TYPE TRUNCATED       CUTOFF_RADIUS 6.0     &END INTERACTION_POTENTIAL     &MEMORY       MAX_MEMORY 3000     &END MEMORY   &END HF &END XC 有的泛函在&XC_FUNCTIONAL里有特定的字段，可以省得自己去定义细节参数，但HF成份依然必须自己定义，否则相当于只做了残缺的纯泛函计算，结果毫无意义。比如HF成分为10%的TPSSh泛函可以写为 &XC   &XC_FUNCTIONAL        &HYB_MGGA_XC_TPSSH     &END HYB_MGGA_XC_TPSSH   &END XC_FUNCTIONAL   &HF     FRACTION 0.1     ...同上   &END HF &END XC 更多的泛函定义方式可以直接看Multiwfn产生的相应杂化泛函的输入文件里的写法，都是严格正确的。 7 一个常见警告：The Kohn Sham matrix is not 100% occupied 做杂化泛函计算，特别是用了ADMM时，很容易在SCF开始后看到The Kohn Sham matrix is not 100% occupied警告。具体原因在https://www.cp2k.org/faq:hfx_eps_warning有介绍，这里就不多说了。对于实际计算来说，看到这个警告时，先把&QS里的EPS_PGF_ORB设为比默认明显更小的1E-12，如果计算正常，即便还有这个警告也可直接无视。EPS_PGF_ORB设得越小，越无法利用重叠矩阵和KS矩阵的稀疏性节约时间，但这对耗时的影响相对于杂化泛函计算的总耗时来说微不足道。目前版本的Multiwfn产生的杂化泛函的输入文件里自动就把EPS_PGF_ORB设成了1E-12。 对于晶胞较小的情况，本身矩阵稀疏性就差，开发者建议上来就在&QS里设MIN_PAIR_LIST_RADIUS -1从而完全不利用矩阵的稀疏性节约时间，此时绝对不会再出现那个警告，但比起用EPS_PGF_ORB 1E-12时耗时会高一些。 8 简单例子：C60分子晶体 最后，这里以C60分子晶体的单点任务作为例子演示一个标准的周期性杂化泛函的计算。如果你还不了解Multiwfn创建CP2K输入文件的功能，建议先阅读《CP2K第一性原理程序在CentOS中的简易安装方法》（//www.umsyar.com/586）。 此体系的cif文件和下面涉及的输入文件都可以在//www.umsyar.com/attach/690/file.zip中找到。此体系结构如下所示，一共240个C原子。晶胞边长是14.07埃，对于分子晶体来说不扩胞的情况下只考虑gamma点可以接受。 启动Multiwfn，输入C60.cif的路径载入之，然后输入 cp2k PBE.inp   //要产生的PBE计算的输入文件名 0   //基于默认设置产生输入文件 现在当前目录下就有了PBE.inp，对应PBE/DZVP-MOLOPT-SR-GTH单点任务。然后接着在Multiwfn里输入 cp2k PBE0.inp   //要产生的PBE0计算的输入文件名 1   //选择理论方法 -6   //PBE0结合ADMM 0   //产生输入文件 由PBE0.inp内容可见，Multiwfn自动指定的计算设置完全合适，对应的是PBE0/DZVP-MOLOPT-SR-GTH结合admm-dzp辅助基组、6埃库仑截断、EPS_SCHWARZ 1E-6、SCREEN_ON_INITIAL_P T、开启OT。 手动把PBE0.inp中的WFN_RESTART_FILE_NAME前头的#去掉并设为PBE-RESTART.wfn，把SCF_GUESS RESTART开头的#去掉，从而使得PBE0计算时读取PBE收敛的波函数当初猜。然后根据你的机子的空余物理内存量和要用的MPI进程数设置合适的MAX_MEMORY。 把PBE.inp和PBE0.inp放到同一个目录下，先运行PBE.inp，算完后当前目录下就出现了PBE-RESTART.wfn，然后再运行PBE0.inp。在《淘宝上购买的双路EPYC 7R32 96核服务器的使用感受和杂谈》（//www.umsyar.com/653）介绍的笔者的2*7R32 512GB双路服务器上用96核并行、MAX_MEMORY设5000，用CP2K 2023.2 popt版，PBE计算花了8秒，PBE0计算花了178秒。PBE0任务的SCF迭代过程如下所示，可见第一轮耗时较高，由于当前给的内存足够储存所有ERI，因而ERI没有on-the-fly计算的，所以之后每一轮SCF的耗时远远低于第一轮。 9 总结 本文把强大的CP2K高效地做杂化泛函计算的各方面要点进行了介绍，以避免初学者对各种关键信息毫不知情而凭感觉胡用瞎用。本文的内容也充分体现出用Multiwfn产生CP2K输入文件非常容易，令研究者从复杂的输入文件编写中充分解放，但这绝不意味着CP2K的使用者就可以对关键性技术细节毫不知情，否则会各种踩坑、白浪费时间、算出无意义的结果。CP2K用户必须掌握的知识、要领远多于大部分其它计算程序。笔者开设的北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）可以使得学员完整系统透彻掌握全部这些知识。 //www.umsyar.com/678 2023-07-28T19:53:00+08:00 谈谈重复不出来计算化学文献里的数据的可能原因 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)文献中的计算本身就有严重错误或是假数据，可能收到你的邮件后作者才意识到，又或者可能之前就已经明知道用的是假数据，总之都不希望你知道他的数据存在问题，免得面子上挂不住，甚至怕文章最后被撤稿。 //www.umsyar.com/671 2023-06-04T16:01:00+08:00 使用Multiwfn结合CP2K的波函数模拟周期性体系的隧道扫描显微镜（STM）图像 Simulating scanning tunneling microscope (STM) images of periodic systems using Multiwfn in combination with wavefunction produced by CP2K 文/Sobereva@北京科音    2023-Jun-4 1 前言 强大的波函数分析程序Multiwfn（//www.umsyar.com/multiwfn）有非常方便好用的模拟隧道扫描显微镜（STM）图像的功能，模拟的原理以及Multiwfn中此功能的用法我在《使用Multiwfn模拟扫描隧道显微镜(STM)图像》（//www.umsyar.com/549）里有专门说明，但那篇博文是以孤立体系为例。Multiwfn支持对CP2K产生的周期性体系波函数进行大量的分析，在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里我做了很多介绍，其中包括模拟周期性体系的STM图。本文的目的是举一个完整的例子，能令甚至之前没用过CP2K的人都能轻易使用Multiwfn模拟周期性体系的STM图。没读过//www.umsyar.com/549的读者必须要先仔细阅读此文。不了解Multiwfn的话建议看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。值得一提的是，CP2K自己也有模拟STM的功能，但远没有Multiwfn来得好用和方便，所以这里就不多说了。 2D Mater., 6, 015005 (2019)文中给出了几层厚度的黑磷薄片表面的0.7 V偏压下的常电流模式的STM图，如下所示，本文将对黑磷计算同样条件下的STM图并与之对照。 本文为了节约时间只用单层黑磷，而且图省事直接从三维晶体中截出来一层结构直接做单点计算得到波函数。严格来说应该用多层模型，并且固定最下层黑磷并令上层的自发弛豫。 本文的例子涉及的文件都可以在//www.umsyar.com/attach/671/file.rar中得到。读者务必使用2023-Jun-4及以后更新的Multiwfn版本，否则情况与本文所述不同。 2 得到单层黑磷的molden文件 如《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（//www.umsyar.com/651）所述，用Multiwfn分析CP2K的波函数必须让CP2K产生记录波函数信息的molden文件。这一节我们就对单层黑磷做一个单点计算来得到此文件。 本文文件包里的Phosphorus-black.cif是黑磷的晶体结构，用GaussView打开，只保留中间的单层部分而删除其它原子，如下所示，然后保存为Phosphorus-black.gjf。 Multiwfn有极其便利的创建CP2K输入文件的功能，见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587），这里我们用Multiwfn创建产生molden文件的单点任务文件。由于CP2K没法产生考虑k点时的molden文件，因此必须构造超胞用gamma点计算。当前的晶胞参数是a=3.31埃、b=4.38埃，在a方向复制为7倍、b方向复制为5倍后两个方向都有>20埃，此时只考虑gamma点是可以接受的。当前模拟STM用的偏压为正，此时能量最低的一个或多个空轨道对STM会产生贡献，因此必须要求CP2K求解出一些空轨道才行，几十个就够了。模拟STM用的偏压越大应当计算的空轨道越多，因为越高的空轨道会被涉及。 启动Multiwfn，载入Phosphorus-black.gjf，然后输入 cp2k   //产生CP2K输入文件 SP.inp  //将要产生的CP2K输入文件名 -7  //设置周期性 XY  //XY二维周期性 -2  //产生molden文件 -11  //进入结构编辑界面 19  //产生超胞 7  //a方向复制的倍数 5  //b方向复制的倍数 1  //c方向（垂直于表面的方向）不变 -10  //返回 -9  //其它设置 12   //计算空轨道 20  //算最低20个空轨道 0  //返回 0  //产生输入文件 此时SP.inp就产生在当前目录下了，在本文的文件包里提供了，计算级别是Multiwfn默认的PBE/DZVP-MOLOPT-SR-GTH。用CP2K运行SP.inp，在我的2*EPYC 7R32双路服务器上不到10秒钟就算完了。之后按照《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）所述，把晶胞信息和价电子信息手动加入.molden文件，即在开头插入以下字段。本文文件包里的SP-MOS-1_0.molden是已经改好的。  [Cell]  23.17000000     0.00000000     0.00000000   0.00000000    21.90000000     0.00000000   0.00000000     0.00000000    12.11600000  [Nval]  P 5 3 用Multiwfn模拟黑磷常距离模式的STM 首先模拟一下常距离模式的STM。启动Multiwfn，载入SP-MOS-1_0.molden，然后输入 300  //其它功能 Part 3 4  //模拟STM 2  //设置偏压 0.7  //0.7V，和文献里的相同 0   //开始计算 默认计算的是Z最大值原子上方0.7埃的XY平面，X和Y方向计算范围和当前晶胞的X和Y方向的跨度一致。瞬间就算完了，在新出现的作图设置菜单里输入0绘图，图像如下所示 可以再做一些调整，只让最上层的P用粗体字标注，下层的用细体字显示以作区分，并且用和文献里相似的配色，并且调整坐标轴让刻度整齐。为此，把图像关闭后，接着输入 4  //修改显示原子标签的距离阈值 0.71 A //在图上显示距离作图平面0.71埃以内的原子的标签 y  //将距离作图平面更远的原子的标签以细体字显示 9  //修改色彩变化方式 13   //黑-橙-黄 7  //修改标签间隔 3,3,0.0005 8  //修改色彩刻度 0,0.002  //下限和上限 1  //保存图像 此时的图像如下所示 此图的效果已经很好了。和实验STM图像展现出的信息近似一致。 4 用Multiwfn模拟黑磷常电流模式的STM 这次模拟常电流模式的STM，这与文献里的STM图对应。启动Multiwfn，载入SP-MOS-1_0.molden，然后输入 300  //其它功能 Part 3 4  //模拟STM 2  //设置偏压 0.7  //0.7V，和文献里的相同 1  //切换模式为常电流模式 0   //开始计算 3   //绘制STM图像 0.0009   //常电流值（反复尝试，取一个图像效果较好、和实验图像较接近的即可） 4  //修改显示原子标签的距离阈值 2.6 A  //在图上显示距离被计算区域顶端2.6埃以内的原子的标签（计算的区域顶端的z坐标是计算开始之前界面里7 Set range in Z direction选项后面显示的第二个值。当前设2.6埃可以只把顶层磷原子纳入进来） y  //将距离更远的原子的标签以细体字显示 9  //修改色彩变化方式 13   //黑-橙-黄 7  //修改标签间隔 3,3,0.1 1  //保存图像 现在看到下图，和常高模式的STM图展现出的信息基本一致 5 其它 本文的例子充分体现出用Multiwfn+CP2K绘制周期性体系表面的STM真是又快又方便，鼓励大家应用在实际当中。请别忘了恰当引用Multiwfn。作为练习，大家可以绘制一下单层MoS2和石墨炔的STM图。 对于无gap的体系，如金属表面（包括吸附小分子后），绘制STM有一点需要注意。这种情况一般都是开smearing的，此时费米能级附近的轨道占据数不为整数，导出的molden文件里的轨道占据数也是如此，此时Multiwfn没法直接绘制STM（Multiwfn会以为是记录自然轨道的多组态波函数）。解决方法是先进入主功能300的子功能9（计算费米能级功能），此功能会将电子按照轨道能量由低到高以整数方式填充来修改轨道占据数，相当于0 K的占据方式。此外，在此功能里输入当前的温度，程序会给你费米能级值，将之记下来。退出此功能后，就可以照常进入STM绘制界面，届时再把费米能级值输入进选项3 Set Fermi level即可。另外值得一提的是，用Multiwfn创建CP2K输入文件时开启smearing的话，产生的输入文件里自动就会有计算空轨道的设置，就不用自己再单独设置了。 //www.umsyar.com/665 2023-05-14T13:25:00+08:00 CP2K中遇到SCF难收敛时的解决方法 Solution when SCF is difficult to converge in CP2K 文/Sobereva@北京科音 First release: 2023-May-15   Last update: 2024-Feb-22 0 前言 经常有人问CP2K中遇到SCF难收敛怎么解决，这是很家常便饭的问题。这里就基于笔者讲授的“北京科音CP2K第一性原理培训班”（http://www.keinsci.com/workshop/KFP_content.html）里面相关部分的内容专门写个博文，供CP2K用户参考。这里假定读者已经了解CP2K计算的基本的背景知识（如k点、CUTOFF、对角化、OT、&BS、smearing等等），如果不了解相关名词、算法的话，十分推荐参加这个培训全面、系统学习一遍。 笔者之前专门写过Gaussian用户遇到SCF不收敛的解决方法的文章《解决SCF不收敛问题的方法》（//www.umsyar.com/61），本文和那篇文章一样，已经把所有可能有用的方法和要注意的问题说得尽可能全面、充分了，读者遇到难收敛/不收敛时需要结合实际情况一个个尝试。切勿问笔者“这文章我看了，但还是解决不了”之类的，我没有任何可额外补充的（至多是给出当前计算的尽可能详细的信息，比如给出输入输出文件，然后我能回答哪些问题需要优先考虑、哪些方法值得优先尝试）。 下面第1节先说说遇到SCF不收敛时一般情况需要考虑的问题和解决策略，然后第2、3节再分别说使用对角化和OT时分别可以进一步考虑的。最后再说一下几何优化、分子动力学时与SCF收敛有关的问题。在此之前先看培训班ppt了解一下CP2K是怎么判断SCF收敛的： 1 一般情况时考虑的 • 检查结构合理性。SCF是否容易收敛和结构的合理性关系极大，化学意义越强的结构通常越容易收敛。因此遇到SCF很难收敛时，需要结合结构化学、固体物理常识检查结构是否靠谱，如模型是否有意义、是否漏了原子、是否某些原子所处的化学环境和实际明显不符（如过渡金属本该与配体成键的地方却光秃秃地悬着）等等。且要结合图像肉眼检查，确保不存在不该出现的结构严重扭曲、盒子边缘原子和周期镜像原子发生过近接触等问题。若有这些问题，即便SCF收敛了，结果显然也毫无意义。 • 确认关键词、设置没明显不合理的。CP2K的关键词很复杂，对使用者的知识水平要求很高，若用了明显不当的关键词和设置，SCF很难收敛或收敛到离谱结果是极其正常的事，然而很多要点是初学者难以注意到的。比如计算自旋极化体系的时候必须得在&DFT里写UKS（或等价的LSD），否则即便自旋多重度设的是>1，程序也会愣当成是闭壳层算。再比如用wavelet的Poisson solver算孤立体系的时候，如果体系有明显电子密度分布的区域跨越了盒子，结果就无意义。 CP2K属于那种是完全手动挡还很个性的车。如果经验不够丰富的话，笔者建议总是用Multiwfn（//www.umsyar.com/multiwfn）产生其输入文件，见《使用Multiwfn非常便利地创建CP2K程序的输入文件》（//www.umsyar.com/587）。Multiwfn的CP2K输入文件创建功能将CP2K尽最大可能包装成黑箱，不仅用起来极其方便省事，对初学者还能大大减小关键词使用不当的可能性。比如将自旋多重度设成>1时自动就会添加UKS、用wavelet做非周期性计算时会自动设置足够大又不太浪费的盒子尺寸并且加入自动让体系在盒子里居中的设置。 • 确保体系的净电荷、自旋多重度的设置合理。这俩如果设置不当的话，相当于计算了无意义或激发态的电子结构，SCF极难收敛也是很可能出现的事。特别是对于磁性体系，自旋多重度的设置对SCF收敛往往影响极大。 • 如果SCF有收敛趋势，可增大SCF迭代次数上限再试，即&SCF里的MAX_SCF，其默认的50对于很多情况明显太小。Multiwfn产生的CP2K输入文件里默认用128（128是为了优雅，2的7次方。和Gaussian程序一致），绝大多数情况都够了，但也可能碰上极个别体系仍需要加大。 笔者超级反感一碰到SCF不收敛，也不观察收敛趋势，上来就增加迭代次数上限的做法，因为这是最典型的菜鸟思维，绝大多数情况毫无用处而白浪费时间。因为如果SCF过程都发生显著震荡了、没收敛的苗头，再怎么迭代下去都是白搭，纯粹是毫无意义地浪费时间和计算资源。 • 如果发现用其它泛函或基组下计算能收敛的话，可以读取其收敛的波函数作为初猜，这至少比自动产生的初猜波函数要好、更容易收敛。也即&SCF里用SCF_GUESS RESTART，且&DFT里WFN_RESTART_FILE_NAME指定成之前任务的.wfn文件（若考虑了k点是.kp文件）。注意两次的基组特征必须定性一致，比如之前计算用的大核赝势，而现在用小核赝势或全电子基组，这时候读之前的波函数就没意义了。 通常较小的基组下计算比大基组下计算更容易收敛，比如直接用TZV2P档次的MOLOPT基组不收敛的话，可以尝试用DZVP档次的算，若收敛了，读它收敛的波函数为初猜。 往往HF成份越大的泛函越容易收敛，很大程度是因为其HOMO-LUMO gap更大、SCF过程中占据和空轨道混合程度相对更缓和所致。因此遇上SCF难收敛时可以尝试用HF成份比实际要用的泛函更大的泛函，若能收敛，读取其收敛的波函数当初猜。HF成份可以参考《不同DFT泛函的HF成份一览》（//www.umsyar.com/282）。 如果原本用的是纯泛函，想尝试杂化泛函但用不动，也可以用DFT+U试试，耗时和纯泛函基本一致，但此时不支持k点，此外怎么获得+U用的Ueff参数也是个问题。如果实际计算不打算+U，也可以尝试3或4 eV的+U看看能否令SCF收敛，如果能，则不+U时读取其收敛的波函数当初猜。 • 做杂化泛函计算时（特别是对于周期性体系）通常会在&DFT / &XC / &HF / &SCREENING里开启SCREEN_ON_INITIAL_P从而基于初始的密度矩阵做双电子积分的屏蔽来巨幅节约耗时，Multiwfn产生的杂化泛函计算的输入文件里默认也开了这个。若此时没有读取之前收敛的波函数当初猜（通常是便宜得多的纯泛函计算的），则CP2K会基于自动初猜的波函数对应的密度矩阵做积分屏蔽，此时结果通常无意义，SCF能收敛的概率也非常低！然而这个杂化泛函计算的关键要点却经常被初学者所忽略！笔者强烈建议在做杂化泛函计算之前仔细把此文认真读一遍：《CP2K做杂化泛函计算的关键要点和简单例子》（//www.umsyar.com/690）。 • 确保恰当考虑了k点、k点数够多。如果晶胞边长不够大，显然必须考虑k点而不能只用默认的Gamma点。而且k点设得还得够多，否则SCF难收敛极为正常。k点如果设置不当，就算SCF能收敛了也没意义。众所周知，算能量问题，原则上k点应当大到足以令能量随k点数增加而基本收敛的程度。 • 尝试增大&DFT / &MGRID里的CUTOFF、REL_CUTOFF。这俩截断能设置控制格点精度，不够大的话不仅结果不准，甚至可能造成SCF难收敛。如果当前明显已经够大了，就没必要再盲目尝试再增加了。注：还有种说法是尝试&DFT / &XC / &XC_GRID里写USE_FINER_GRID T将计算交换-相关泛函部分的格点精度提升为CUTOFF的四倍，但根据我的经验这并没有什么实际用处，而耗时暴增。 • CP2K默认的产生初猜波函数的方法是孤立原子密度的叠加，在整个任务一开始时会对每个KIND对应的原子在孤立状态下做计算，用PRINT_LEVEL MEDIUM时可以看到具体细节。计算孤立原子所用的电子组态可以通过&KIND里的MAGNETIZATION（定义alpha电子数减beta电子数）或者&BS（可以具体控制各个主量子数的各个角量子数上的alpha和beta电子数）来设置。初猜波函数中的孤立原子的组态越接近整个体系收敛的波函数中原子的实际组态则SCF越容易收敛、SCF不收敛的可能性越低。 如果体系有磁性，遇到SCF不收敛时一定别忘了尝试用MAGNETIZATION或&BS尽可能令初猜波函数中原子组态与实际情况接近，这对SCF是否能够顺利收敛有关键性影响。如果当前体系在文献（可用google搜索来找）或高通量数据库（如Material Project）中已经被计算过，可以参考那里面报道的原子自旋磁矩或原子组态来设。有时同一种元素在体系里处于不止一种化学环境，差异很大时（比如明显处于不同价态）可分成多个&KIND分别设置。 即便体系没有磁性，利用&BS使初猜波函数中的原子带电状态使之尽可能符合实际也往往是有益的。比如对于NaCl这样高度离子性体系计算时SCF不收敛的话，可以尝试令初猜中Na和Cl的原子组态分别对应Na+和Cl-的状态（默认情况下分别对应中性的Na和Cl）。曾经我用PBEsol/pob-TZVP-rev2计算10层的NaCl板时通过这种做法成功解决了SCF不收敛问题。 • 尝试将&QS / EPS_DEFAULT减小几个数量级，默认为1E-10。EPS_DEFAULT这个参数控制一堆涉及数值计算精度的参数，如果设得太大，有可能造成SCF难收敛。若这个数值当前已经很小了，比如都1E-12了，尝试继续减小则没必要。 • 改用其它泛函、基组。如果原先的计算级别就是极难收敛，而有其它级别性价比、精度相仿佛，且发现能收敛，那就别死磕当前的了。 • 在GPW与GAPW之间切换。如果之前用的GPW结合GTH赝势基组不收敛，尝试用GPAW结合一个全电子基组计算试试，反之亦然。 • 修改默认产生波函数的方式。参见手册里&SCF / SCF_GUESS设置。但其实这没什么用，除了默认的ATOMIC（即前述的基于孤立原子密度叠加产生的）以及RESTART（读取之前的波函数当初猜），其它的一个能打的都没有。 用ATOMIC产生初猜时，CP2K计算每个KIND对应的孤立原子也是需要SCF迭代的，上限是固定的50轮（没法靠关键词改），在PRINT_LEVEL MEDIUM时其迭代过程会在一开始输出。有时候我发现到了50轮时收敛精度还很差，故有可能在某些情况下由于这个问题导致初猜质量烂，不利于实际体系的SCF收敛。我发现可以通过atom_kind_orbitals.F里atom%optimization%max_iter修改迭代次数上限，比如改成300然后重新编译，一般孤立原子的SCF就都能收敛了。 • 如果加了外电场时SCF很难收敛，可能是外电场过大所致，可尝试更小的数值。如果就是要用较大的外电场，可以试试读取较小电场或不带电场时收敛的波函数当初猜。有的体系在很强外电场下本来就不可能稳定存在，这时候SCF不收敛自然很正常，是物理上决定的，也别盲目去计算。 • 实在没辙时干脆放宽SCF收敛限，即加大&SCF / EPS_SCF，默认是1E-5（这其实已经挺宽松了），这是走投无路时的做法。或者换程序，如免费且和CP2K同样流行的Quantum ESPRESSO，或者尝试CP2K里的纯粹基于平面波的SIRIUS模块。某些体系用CP2K的quickstep模块就是很难收敛，如某些d族纯金属的slab，这和算法、基函数特征有关。 2 用对角化的情况可专门考虑的 • 对无gap（如金属、石墨烯）或gap很小的体系（非常窄带半导体）应当使用smearing。即加入&SCF / &SMEAR / METHOD FERMI_DIRAC，且在&SCF里用ADDED_MOS要求求解一些空轨道从而能被热激发的电子所占据。ADDED_MOS数目够用就行，设太大了会增加耗时而且也没用处。Multiwfn产生CP2K输入文件时如果开了smearing会自动做ADDED_MOS设置，并根据体系原子数估计一个通常够用的值。用较小电子温度（如200K、300K）时若不收敛可尝试再提升温度，即&SMEAR / ELECTRONIC_TEMPERATURE设的。温度设高时ADDED_MOS往往需要相应地增加，否则热激发的电子没法按Fermi-Dirac分布占该占的空轨道（此时计算时会有警告）。 我老在网上看到有人不理解smearing原理是什么就胡用、乱试smearing以试图解决SCF不收敛，所以这里强调一下。一定要注意，对于gap不是很小的体系（比如氧化锂等绝缘体），根本没必要尝试smearing，本来这在原理上就不是smearing能起到帮助SCF收敛的情况，而且本文还有其它一堆方法可以尝试用于帮助SCF收敛，完全轮不到考虑用smearing。 用smearing还会对结果产生一定影响，结果会依赖于电子温度T的设置，用Fermi-Dirac函数做smearing时最终给出的能量相当于是T温度下电子的自由能。如果你就是要得到基态（0K）的电子能量，原则上最好做T→0的外推。如果你不做外推，在能量横向对比时至少也应当用相同的smearing设置。关于smearing的物理意义和实际意义，以及做外推的方法，我在北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里有专门详细的讲解。 • 对gap较小的体系可尝试使用能级移动，即&SCF / LEVEL_SHIFT，设比如0.3、0.4左右，默认单位是Hartree。这使得SCF过程中人为加大空轨道与占据轨道之间的能级差，从而减轻彼此间混合的剧烈程度而有时令SCF更容易收敛。 • 修改&SCF / &MIXING / METHOD设置SCF过程中将旧密度矩阵与新密度矩阵混合得到实际每一轮密度矩阵的方法。默认的DIRECT_P_MIXING往往表现极差。通常用BROYDEN_MIXING，这也是Multiwfn产生的输入文件里默认的。亦可尝试PULAY_MIXING，但一般不比BROYDEN_MIXING好。 • 对于使用BROYDEN_MIXING的情况，尝试修改&MIXING中的ALPHA（新密度矩阵混入旧的的比例，一般用默认的0.4，即混入40%。难收敛时可尝试诸如0.1、0.2、0.3）和NBROYDEN（默认的4明显太小，改大到8通常很有助于收敛，这是Multiwfn产生的输入文件默认的。也可以尝试进一步改大到比如12） • 对于使用PULAY_MIXING的情况，尝试修改&MIXING中的NPULAY （默认的4偏小，改大往往很有助于收敛）。注：NPULAY、NBROYDEN、NBUFFER三个关键词实际上等价。 • 对于使用DIRECT_P_MIXING的情况，收敛到&SCF / EPS_DIIS值后会切换到加速收敛的DIIS算法。但默认的EPS_DIIS为0.1，有点偏小，可以尝试改大以早点启用DIIS（但借助DIIS通常仍不如BROYDEN_MIXING）。还可以把控制最大DIIS矢量的&SCF / MAX_DIIS从默认的4改大到诸如7。 • 在对角化框架下各种方法都不行的话，尝试改用OT。如果OT能收敛的话，而且当前由于特殊目的又必须用对角化，也可以用OT收敛的波函数当对角化计算的初猜。由于OT不支持k点，因此如果之前对角化是结合k点算的较小晶胞情况，用OT时就只能靠超胞模型了。OT也不支持用ADDED_MOS设置在SCF计算过程中求解空轨道，因此必须用smearing的导体体系就别指望OT了。如果是用GFN1-xTB，总应当用OT，这也是Multiwfn默认的。 3 用OT的情况可专门考虑的 • 如果之前没用outer SCF则开启之，这对解决SCF不收敛非常有效，也即加入&SCF / &OUTER_SCF字段。由于其重要性，Multiwfn产生的输入文件里开OT后默认就是启用outer SCF的。此时每次到达&SCF / MAX_SCF所设的inner SCF的步数上限时，就会进入outer SCF过程更新OT的preconditioner，然后进入下一次的inner SCF迭代过程。每次更新preconditioner后经常会看到比之前以更快的速度趋近收敛限。原理上来说，OT每一步都更新preconditioner是最理想、SCF最容易收敛的，但由于这样太昂贵，因此默认情况下（即不用outer SCF时）只会在最开始计算一次preconditioner。 开启outer SCF时inner SCF的迭代次数上限应当设得比平时更小（否则经历很多次SCF迭代才做一次outer SCF，无法有效达到用outer SCF的目的），15至35是比较合适的范围。outer SCF迭代次数上限通过&SCF / &OUTER_SCF / MAX_SCF控制，总SCF迭代次数上限相当于inner SCF次数上限和(outer SCF次数上限+1)的乘积。 &SCF / EPS_SCF设的是inner SCF收敛限，&SCF / &OUTER_SCF / EPS_SCF设的是outer SCF收敛限，后者应当小于或等于前者（一般设为相同）。 • 尝试其它的preconditioner。FULL_ALL通常最好，但最贵。大体系可以尝试明显更便宜的FULL_SINGLE_INVERSE和FULL_KINETIC。个别时候FULL_KINETIC反倒比FULL_ALL更容易收敛。GAPW计算总应当用FULL_ALL。 • 修改OT算法，把&OT里的ALGORITHM从默认的STRICT改为IRAC，SCF经常会收敛得更稳健，非常值得尝试！特别是我发现用pob系列全电子基组的时候，用IRAC往往容易收敛得多，因此Multiwfn产生基于pob系列基组的输入文件时若用户选了OT，则默认用IRAC。而且用约束性DFT（CDFT）时，我发现用IRAC也会令SCF容易收敛得多得多。 • &OT / MINIMIZER控制OT用的极小化方法，若之前用的是较快的DIIS（也是Multiwfn产生的用OT的输入文件里默认的），可改为更稳健但耗时更高的CG（共轭梯度法），也可以尝试BROYDEN。 • MINIMIZER用CG时，尝试用&OT / LINESEARCH选项把线搜索方式设为比默认2PNT更贵但也更稳的3PNT。 4 几何优化和分子动力学过程中的SCF收敛问题 在几何优化过程中，有可能中途偶尔出现几步SCF不收敛。此时程序还会继续算下去。只要最终几何优化收敛时SCF也是收敛的状态就行了。（注：从CP2K 2024.1开始，优化中途SCF不收敛会导致任务自动停止，&SCF里写IGNORE_CONVERGENCE_FAILURE才会在SCF不收敛的情况下也继续优化） 有时候在几何优化初期SCF一直不收敛，不要急于把任务停掉，可以先凑合跑跑看看情况。初始几何结构相对于实际的极小点结构而言肯定是存在多多少少不合理性的，此时SCF也肯定比极小点结构处更难收敛。在SCF不收敛、受力也较糙的情况下凑合跑一定步数后，结构一般会变得比一开始更为合理，此时SCF也往往变得更容易收敛。 在几何优化、分子动力学过程中，程序会自动利用之前一些结构的波函数外推出当前步的初猜波函数，这比起直接产生的初猜波函数质量通常好得多。外推方式可以通过&DFT / &QS / EXTRAPOLATION来具体设置，默认的ASPC通常是最好的选择。由于这一点，当几何优化、分子动力学跑一定步数之后，SCF达到收敛限所需的迭代次数明显减小、出现SCF不收敛的可能性也大为降低。 注意当使用k点时CP2K不支持外推产生初猜波函数，因此会在几何优化、分子动力学等涉及结构变化的任务中每一步重新产生初猜波函数，即EXTRAPOLATION USE_GUESS的情况。我一般建议改为EXTRAPOLATION USE_PREV_P，代表直接使用上一步SCF收敛的波函数当初猜。 如果想减小MD中途出现SCF不收敛的可能，可以尝试减小步长（&MD / TIMESTEP）。如果想减小几何优化中途出现SCF不收敛的可能，可减小优化的步长上限，即置信半径（如通过&MOTION / &GEO_OPT / &BFGS / TRUST_RADIUS）。这样做之后由于相邻两个结构间差异变小了，ASPC或USE_PREV_P给出的初始波函数对于当前结构而言就越理想，自然SCF收敛也会更容易，且有利于让SCF能收敛的状态保持下去。