On the new features and changes of ORCA 5.0

文/Sobereva@北京科音

First release: 2021-Jul-2  Last update: 2021-Jul-10

1 前言

2021年7月1日， 程序ORCA 5.0发布了，造成了不小影响。这个程序现在流行度已经相当高，从《2021年计算化学公社论坛“你最常用的计算化学程序和DFT泛函”投票结果统计》（//www.umsyar.com/599）的投票结果上可见一斑。这里就简单说下这个更新对于大多数用户值得一提的新功能和改变，以及我的一些相关看法，便于读者了解新版本都带来了什么。更全面的更新说明见5.0版手册1 ORCA 5.0 Foreword和2 ORCA 5 Changes部分。

5.0版手册比4.2.1扩增了近300页，达到了1300多页的规模。我所知的手册最厚的量化程序是Q-Chem（5.3.2版是1391页）。按照ORCA的势头，据说将在2022年发布的ORCA 6的手册厚度预计将超过Q-Chem。

2 比较重要的改进

5.0的一个很大进步是加入了全新开发的积分代码SHARK。号称对于低角动量和高角动量基函数比以前版本用的Libint电子积分库更快。ORCA自动决定用SHARK还是Libint来达到最好的效率。据说此改变令解析Hessian、TDDFT、磁相关属性速度提升较大。SHARK既适合片段收缩基组，对广义收缩基组也专门做了优化，因此能算得动ANO那种完全广义收缩的基组了。笔者实际测试了一下，ANO这种完全广义收缩的基组确实在ORCA 5.0里能算得动了，而在4.2.1里哪怕用ANO-pVDZ算个很小的分子都极其吃力。我用带有高角动量函数的基组cc-pVQZ也简单测试了一下算一个小体系HF的耗时，5.0比4.2.1耗时低了百分之十几。

DFT积分格点和COSX积分格点利用机器学习技术进行全新的设计。目前版本定义了三种格点的关键词，defgrid1、defgrid2和defgrid3，控制所有涉及到积分格点步骤用的格点（COSX、CPSCF、TDDFT等）。默认是defgrid2，据称没有特殊情况这个档次的格点质量就足够了，而对积分格点要求特高的情况（近乎Gaussian里int=superfine的情况），比如VPT2做非谐振计算的时候，才需要defgrid3。另外，使用对积分格点要求较高的明尼苏达系列泛函的时候程序会提示建议用defgrid3。defgrid1的大小近似相当于之前版本默认的挺糙的格点。由于默认的defgrid2比之前的积分格点明显更好，因此DFT计算的数值精度比以往版本有明显提升，而且这还令DFT的SCF收敛、几何优化收敛等也往往更为容易。据说即便把之前版本的格点数设得和当前版本的差不多，新的积分格点的精度也更好，这是因为格点分布调教得更科学了，期间还利用了GMTKN55数据集进行了训练和测试。defgrid系列格点还有一个特点就是有自适应能力，对于带弥散的基组、带紧s基函数的基组、带核极化函数的基组，新版本在对格点剪裁的时候有相应的特殊考虑来保证数值精度。特别要注意的是，之前版本分别控制DFT和COSX格点的grid和gridx关键词不仅不能用了，而且连出现都不行，这导致原先大量输入文件都不能用了，这点很变态（哪怕自动无视也好）。

新版本里COSX计算解析积分的代码效率得到了优化。结合新的COSX格点，令数值稳健性也显著提升、SCF和几何优化达到收敛的步数有可能更少。由于RIJCOSX已经非常成熟、理想了，因此从5.0开始RIJCOSX成为了杂化和双杂化泛函计算时默认开启的加速方法（以前版本只是纯泛函默认用RIJ）。如果要关闭RIJCOSX的话，写noCOSX。由于RIJCOSX已变得更完美，RIJK可以算是没有任何用武之地了，而且用RIJK时功能限制很大。

提升了CPSCF的效率，因此做freq等需要求解CPSCF方程的任务耗时有所降低。

隐式溶剂模型得到了提升。一个很大的改变是5.0用的CPCM的默认设置是范德华型的溶质孔洞结合Gaussian分布的表面显著电荷，而老版本默认用的设置是gepol_ses型溶质孔洞结合离散的表面点电荷，新版本默认的溶剂模型的数值稳定性更好，解决了原来版本溶剂模型下可能出现的势能面不连续问题，这也使得从头算动力学过程结合溶剂模型的时候更令人放心了。另外，新版本的溶剂模型能用的场合更多了，新版本可以在CPCM溶剂模型下算解析Hessian是个巨大的进步，不过尚不支持SMD下的解析Hessian。并且新版本的CPCM可以结合(DLPNO-)CCSD(T)使用。

由于上述数值方面的各种改进，估计计算Hessian容易出虚频的毛病会得到明显改善。

CPCM（线性响应的形式）正式支持了TDDFT。我发现之前版本的ORCA在溶剂模型下得到的TDDFT激发能不合理，和Gaussian的线性响应溶剂模型相差很大，貌似是因为之前版本没有考虑隐式溶剂对电子激发的响应。我找了个D-pi-A体系做了个简单测试（氨基-苯-硝基），用PBE0/def2-SV(P)结合SMD描述的水环境做TDDFT计算，考察最低激发能和振子强度。G16算出来是3.77 eV，f=0.44；ORCA 5.0算出来是3.76 eV，f=0.45；ORCA 4.2.1算出来是3.90 eV，f=0.36。明显ORCA 5.0比4.2.1合理多了，和G16基本一致了。注：ORCA里的SMD的极性部分用的是CPCM，和Gaussian的SMD的极性部分用的IEFPCM有异，结果注定不可能非常相近，但差异不会太显著。

ORCA作为MPI并行程序，内存消耗量问题一直为人所诟病。5.0在SCF和CPSCF过程中部分程度使用了共享内存方式节约内存消耗量。希望以后能想办法把DLPNO的巨大内存花费降下来。

支持了VV10色散校正的解析梯度，wB97M-V、wB97X-V之类的带VV10的泛函终于都可以用于优化了，以前这是Q-Chem的专利（阿Q一直都活在Gaussian的阴影下，现在免费的ORCA又来势汹汹，把阿Q的卖点一个一个拔掉，在二者的夹缝中日子真不好过）。不过wB97M-V虽然支持了解析梯度优化，但我感觉也没什么实际意义，毕竟还不支持解析Hessian。而优化完了不做个振动分析检验虚频的话很难发文章。而且wB97M-V虽然能量计算精度很好，但并不代表几何优化方面会有什么突出的。实际上B3LYP-D3(BJ)优化就已经挺理想了，这点我在《谈谈 研究中什么时候用B3LYP泛函优化几何结构是适当的》（//www.umsyar.com/557）中已经说过；而且就算遇到B3LYP-D3(BJ)不适合的时候，还有wB97X-D3或wB97M-D4可用，就算真的优化结果比wB97M-V差一些，也就是毫厘之间。

终于支持了meta-GGA、meta-杂化泛函的解析Hessian了，可以用诸如M06-2X在优化完之后做振动分析了。

加入了一些新泛函：B97M-D4、wB97X-D4、wB97M-D4。wB97M-D4的精度略逊于wB97M-V一丁点，但好处是有解析Hessian。不过我实测了一下，发现这泛函结合小基组下对苯分子这么简单的体系做振动分析时CPSCF都特别难收敛，所以在opt freq方面我不觉得这泛函有什么实用性。双杂化加入了wB88PP86、wPBEPP86、RSX-QIDH、RSX-0DH、PBE-QIDH、PBE0-DH，以及它们的为激发态计算优化参数的SCS、SOS形式的变体。不过这些双杂化泛函在我来看并没什么用处，算基态不灵，算激发态的话，对于普通价层激发肯定不如DSD-PBEP86，而算显著的CT态、里德堡态比起之前就有的wB2GP-PLYP也不会有什么优势。

还有一个新加入的双杂化泛函是2009年提出的wB97X-2。这个泛函带DFT-D3(BJ)的版本wB97X-2-D3(BJ)在《简谈 计算中DFT泛函的选择》（//www.umsyar.com/272）里介绍过，对于除了分子间弱相互作用以外，这个泛函通常是最佳选择之一（虽然逊于后继者wB97M(2)，但目前wB97M(2)仅有Q-Chem一个程序支持）。不过我发现ORCA 5.0没有自带这个泛函的DFT-D3(BJ)参数，所以没法直接用wB97X-2-D3(BJ)，而必须自己定义参数才行。为了大家用起来方便，我在《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490）介绍的Multiwfn的相应功能中已经直接加入了产生wB97X-2-D3(BJ)输入文件的选项。

支持了刚提出不久的r2SCAN-3c。其构建思想类似于已经流行开来的B97-3c，而其中的纯泛函部分改用了2020年底才提出的SCAN泛函的第二个修改版，即r2SCAN。虽然r2SCAN-3c耗时比B97-3c高百分之几十，但精度有了全面提升，因此多付出的代价是划得来的。建议r2SCAN-3c能算得动的情况一律用它代替B97-3c。由于其重要性，我在《详谈Multiwfn产生ORCA 程序的输入文件的功能》介绍的功能中也已经加入了产生其输入文件的选项。

当默认的SCF solver收敛失败或者表现不佳时，会默认切换为新加入的augmented Hessian trust radius SCF converger (TRAH)。TRAH比默认的SCF solver要慢，但据说收敛能力很强，号称除非结构不合理、基组有问题、有线性依赖这些严重硬伤外几乎总能收敛，并声称对付收敛普遍偏难的含过渡金属的体系很好。至于实际效果如何笔者没试过，但至少是又多了一个可能解决SCF不收敛的选择。想关闭的话可以用noTrah。

NEB功能得到了改进，并且支持结合TDDFT对激发态势能面找过渡态和反应路径。

支持了共线形式的spin-flip TDDFT以及其解析梯度。这个挺有意义，一些圆锥交叉点的搜索都可以用这个做了。以前spin-flip TDDFT用得比较多的是GAMESS-US和Q-Chem，而前者用着麻烦，后者又收费。

改良了NEVPT2、CASPT2的计算形式，号称对大活性空间情况提升了计算速度。

动力学模块支持了一维和二维的metadynamics获得自由能面。集合变量可以用距离、角度、二面角、配位数，还可以加谐振或高斯限制势。支持了Nose-Hoover chain压浴，支持了velocity-rescale热浴（这个很重要，之前只能用Berendsen热浴，很容易被瞎较真的审稿人吐槽）。新支持了消除质心运动的设置以避免模拟过程中漂移，手册里搜CenterCOM（在之前版本中为了消除漂移，我都是靠质心约束来实现的）。

支持了两层和三层ONIOM，并支持电子嵌入，可以结合到优化、扫描、NEB、IRC、振动分析等各种任务。并且改进了原先的QM/MM。不过我感觉，做ONIOM没有个gview这样的专门开发的界面还是不够方便直观。

DLPNO-MP2和DLPNO-双杂化可以算NMR和极化率了。这对于需要算较大体系高精度NMR者是个福音，有测试表明DLPNO-DSD-PBEP86算NMR很好，不仅胜于普通泛函，比MP2也更好。相比之下，Gaussian的MP2或双杂化算稍大一点体系的NMR都非常吃力。另外值得一提的是，当体系不大时，这种DLPNO形式和RI形式计算耗时相仿佛，而当体系很大时，DLPNO形式耗时显著低于RI形式。

支持了Ionic-Crystal-QMMM和Mol-Crystal-QMMM以簇模型分别计算离子晶体和分子晶体中的局部区域，背景电荷可以自洽地自动确定（对于离子晶体的情况，中心部分和背景电荷之间还有一层capped ECP）。利用这种模型可以算晶体中的电子光谱、NMR，乃至算离子晶体的带隙，还可以算表面吸附问题（这比起第一性原理程序里能用的理论方法丰富得多，还有准确得多得多的方法可用）。对Ionic-Crystal-QMMM方式的计算，ORCA中新加入了专门的orca_crystalprep工具用来构造模型。

在workflow控制方面做了加强，可以在输入文件里控制多个任务的执行、数据提取（计算过程中涉及的大量信息都可以通过变量直接访问）、输出、数学运算，还可以用循环、条件判断，相当于在输入文件里可以直接用脚本语言进行非常灵活的计算流程的控制，也因此可以方便地实现热力学组合方法，而不必依赖于外部的Shell、Python等脚本。甚至还可以靠脚本实现自动循环各个振动频率，发现有虚频则自动按虚频模式调节结构，之后重新做优化，反复如此直到没有虚频为止。

支持了通过简单的关键词实现很多组合方法。比如写compound[W2-2]就可以做W2.2高精度热力学组合方法的计算。支持的组合方法关键词看ORCA 5.0手册9.47.3 List of known Simple input commands部分。

3 其它一些值得一提的改进或变化

新版本默认对双杂化启用RI，但做纯粹的MP2则不默认用。可以noRI关闭RI。

DLPNO支持了multi-level，对不同片段内部和片段间的电子相关用不同的PNO阈值考虑，从而在不明显损失精度的情况下节约耗时。

SARC ZORA/DKH基组支持了Rb-Xe，这是J. Comput. Chem., 41, 1842 (2020)中提出的。

已经可以直接使用月份基组了，现在ORCA的基组库已经做得挺完善了。

加入了J. Chem. Phys., 152, 214110 (2020)中提出的CASPT2-K方法，入侵态问题比CASPT2更小，和IPEA  shift一样对CASPT2结果大有改进，而不需要IPEA shift那种经验性的参数，因此从本质上更好。此方法算激发能不如CASPT2-IPEA和NEVPT2，不过算Cr2的解离曲线则比CASPT2加不加IPEA都好。此方法只能算是对CASPT2自身的一种改良，相对于NEVPT2并没什么实际优势。

Iterative configuration expansion (ICE)这个功能之前就有，属于selected CI的一种，利用微扰方法筛选出对Full CI贡献较大的组态用于计算，从而在精度牺牲很小的前提下大幅减少传统做Full CI要考虑的组态数，因此比Full CI能算更大一些的体系。5.0对这个功能进行了强化（近期Neese在JCC和JCTC上发了相关文章）。做超高精度计算的人可以尝试。

算激发态支持了SCS-CIS(D)和SOS-CIS(D)，不过适用场合很有限。当各种泛函下TDDFT结果都不理想，而DLPNO-STEOM-CCSD又算不动的话可以试一把。

支持了DLPNO-NEVPT2-F12。支持了完全内收缩多参考耦合簇（FIC-MRCC）。

DLPNO系列方法支持了做PNO外推。DLPNO加速方法受TCutPNO参数影响，此外推是在一个较大（较便宜）和一个较小（较昂贵）的TCutPNO下分别计算能量，再根据简单的公式经验外推到TCutPNO=0的极限得到更精确的能量。写compound[EXTRAPOLATE-PNO]关键词的话会在DLPNO-CCSD(T1)/cc-pVTZ下做自动这种外推。

支持了VPT2非谐振模型得到振动平均的NMR化学位移和EPR超精细耦合常数，由此还可以体现不同同位素的差异。

加入了ANISO程序用来计算自旋哈密顿参数和属性。

ab initio ligand field analysis (AILFT)做了改进。

支持了化学位移的定域轨道的分解分析，并给NBO的NCS（自然化学屏蔽）分析留了接口。

4 一些吐槽和局限性

windows版只有静态库版本，解压完了竟然占24.8GB，而4.2.1版则只有9.2GB而已。

PWPB95、PW6B95、DSD-PBEB95没解析梯度，这个之前版本的局限性到现在还没解决。

可惜没有直接支持比DSD-PBEP86-D3(BJ)精度更好的后继者revDSD-PBEP86-D3(BJ)。它在现有双杂化泛函中精度是名列前茅的。

双杂化泛函依然不支持解析Hessian。

meta-泛函做TDDFT时没有解析梯度。

算拉曼活性依然还得用数值Hessian，是个硬伤，算中等体系都超级耗时，和Gaussian完全没法比。

TDDFT的解析Hessian什么时候ORCA才能支持是个未知数，不过我相信ORCA早晚会支持的。这对于一个普适型、大众型，而且又对激发态计算特别重视的程序来说太重要了。

新版本振动分析的输出格式发生了变化，迫使我今日更新了《OfakeG：使GaussView能够可视化ORCA输出文件的工具》（//www.umsyar.com/498）里介绍的OfakeG程序。除了这个外，新版本还不允许出现grid和gridx关键词，也迫使我今日更新了Multiwfn的相关功能。据说2022年可能要出ORCA 6，但愿出下个大版本的时候别再折腾输出格式，并且考虑对以往输入文件的兼容性（之前ORCA 4出来的时候相对于3.x在基组写法方面就变化甚巨）。老折腾的话对于写辅助程序、写相关教程、维护input library的人是很大负担。

Several methods to straighten molecular structures using Gaussian or ORCA program

文/Sobereva@北京科音

 First release: 2021-May-4   Last update: 2021-Jun-15

1 前言

有人问怎么把一个分子从已有的较为卷曲的构象转化成比较平直的构象。这有一些实际意义，比如通过genmixmem程序构造磷脂双层膜、用packmol程序构造表面活性剂的囊泡，都需要用户提供分子比较平直的构象的结构文件。一个非常简单、直观、效果又好的做法是使用SAMSON可视化程序里的twist工具，笔者在《生成混合组分的磷脂双层膜结构文件的工具genmixmem》（//www.umsyar.com/245）里还提供了我录的一段操作视频。但可惜如今SAMSON程序免费版当中已经没有twist工具了，需要买收费版才行（虽然收费版也可以申请免费试用）。还有一个办法是在GaussView等程序里一个一个改二面角，但无疑这非常费时费力，特别是分子可旋转的键很多的时候。

在此文，笔者介绍如何利用Gaussian和ORCA 程序以不同方式实现把分子从卷曲的结构拉直，里面用到的做法对于读者研究其它一些问题可能也会有启发。本文使用的例子是《将Confab或Frog2与Molclus联用对有机体系做构象搜索》（http://bbs.keinsci.com/thread-20063-1-1.html）中的例子，即Actos分子。通过Molclus程序（http://www.keinsci.com/research/molclus.html）搜索出的此体系的能量最低构象如下，我们后文将要把此结构拉直：

本文涉及的所有文件可以在//www.umsyar.com/attach/594/file.rar下载，上面的结构是文件包里的Actos.xyz。Gaussian使用的是G16 B.01版，ORCA是4.2.1版。

2 通过ORCA拉直分子：给两个原子间加外力

笔者写过不少与ORCA的安装和使用有关的文章，参看//www.umsyar.com/category/ORCA/，相关常识性问题这里就不再说了。在ORCA中优化时可以给两个原子在连线方向上加外力驱使它们远离，很适合用于拉直分子。对于上面图中展示的Actos分子，显然我们应当把最末端的两个重原子，即8和25号之间拉远。为此，写一个ORCA输入文件：

! AM1 opt nopop
%geom
POTENTIALS
{ C 7 24 3.0 }
end
end
* xyz   0   1
[Actos.xyz里的坐标]
 *

这里POTENTIALS下面花括号里的C代表constant force。注意原子序号是从0开始算的。3.0是把俩原子间拉远的力的大小，单位是nN。数值大小不是很重要，一般个位数大小就可以。数值太小的话拉不开，而太大的话会最终把一些化学键过度拉长、键角被过度拉大（不过键伸缩、键角弯曲相对于二面角扭转来说是很刚的，只要外力不是特别大这就不是明显问题）。这里用的AM1半经验方法是ORCA里能直接用的几乎最便宜的方法（虽然ORCA也支持力场计算，但还得搞参数，这里不考虑），注意此方法不支持并行计算。如果嫌AM1太low，也可以用稍微好一点的HF-3c（Hartree-Fock结合额外校正项），对不太大的体系也非常便宜。

用ORCA运行上面的输入文件，三分钟就算完了，最终得到的结构如下（即file文件包里ORCA目录下的Actos.xyz），可见完美地实现了我们的目的，分子被拉得很直

file文件包里的Actos_trj.xyz是优化轨迹，感兴趣者可以将之拖到VMD程序（http://www.ks.uiuc.edu/Research/vmd/）里看看优化过程。顺带一提，考察加外力对化学体系的影响专门有人做过大量研究，参见Chem. Rev., 112, 5412 (2012)、Angew. Chem. Int. Ed., 58, 5232 (2019)。

3 通过Gaussian拉直分子

在Gaussian程序中不支持直接给原子加外力，但有其它的方法也可以实现把分子拉直，将在下面几节介绍。

3.1 方法1：柔性扫描

Gaussian中做柔性扫描在《详谈使用Gaussian做势能面扫描》（//www.umsyar.com/474）里笔者详细介绍过。显然，做两端的两个原子距离逐渐拉长的柔性扫描最终就可以得到被拉得平直的结构。写一个这种任务的Gaussian输入文件，如下所示

%nprocs=1
#T UFF opt(loose,nomicro,modredundant)

Generated by Multiwfn

  0  1
[Actos.xyz里的坐标]
 
8 25 S 50 0.5

这里有几个需要注意的地方：
(1)为了耗时尽可能低，这里使用了Gaussian支持的最便宜的方法，即分子力场。这里用的是较粗糙但支持元素最广、用着最省事（不需要指定原子类型）的UFF力场。由于本例只需要得到一个较为平直的构象，对质量什么没要求，所以UFF就够了。
(2)柔性扫描每一步相当于一次限制性优化，为了节约时间，用了loose收敛限降低每一步扫描的耗时。
(3)opt里必须写nomicro。因为分子力场的优化在Gaussian里有专用的代码，不受modredundant设置的控制。而加了nomicro后，就会用通用的几何优化代码，虽然远没有分子力场专用的优化代码快，但此时可以设置柔性扫描。
(4)分子力场计算时，至少对于这个体系，并行计算会令耗时不降反升，所以这里靠%nprocs=1刻意要求不用并行计算。
(5)扫描的步长太小的话把分子拉直需要的步数太多，显然不行，而步长太大则容易导致中途报错（此时还没充分拉直），根据我的经验步长设0.5埃左右比较合适。由于事先往往不好估计拉直时两个原子相距多远，所以可以把扫描步数设得比较大（比如本例设50步），当分子已经被拉得很直后，继续扫下一步时程序通常会报错，此时用GaussView载入输出文件查看扫描轨迹，若最后一帧结构就是想要的平直的结构的话就取最后一帧结构即可。如果步数一开始设小了，导致最后没有完全拉直也没关系，拿最后一步的结构继续做拉长的扫描即可。
(6)当前体系是之前我用DFT优化过的结构，已经很合理了，所以没有明确靠geom=connectivity指定原子间连接关系，而直接让Gaussian根据当前结构猜连接关系并用于力场计算。如果初始结构没有优化过，最好让GaussView保存出带有恰当连接关系的输入文件以免自动判断的连接关系有误。

使用Gaussian对上面的输入文件进行计算，算了1分钟左右扫描到第31步后报错，报错前最后一步结构如下所示，可见是我们想要的平直的结构。

3.2 方法2：利用外电场加外力

Gaussian里可以加外电场，外电场会给带电粒子施加额外的力，在《一篇文章深入揭示外电场对18碳环的超强调控作用》（//www.umsyar.com/570）里介绍的笔者的文章中有很多介绍和讨论。因此，可以在UFF力场计算时，把25号原子坐标冻结住，而让8号原子带+1.0原子电荷，然后按照下图所示向X轴正方向加足够强度的电场，这样8号原子就会受电场力往X正方向运动，从而实现拉直分子的目的。

此例输入文件如下

%nprocs=1
# UFF opt(loose,nomicro,modredundant) field=x-300 nosymm

Generated by Multiwfn

  0  1
N      2.96734200   -1.21688000    2.24577900
C      2.32051000   -1.62618600    3.34057500
C      1.63462400   -2.84168400    3.45382100
C      1.64984500   -3.67567100    2.33272400
C      2.32098400   -3.26708600    1.18384000
C      2.96534100   -2.02795600    1.17397000
C      0.96378000   -3.24289700    4.74461900
C--1.0      1.95041300   -3.88295500    5.73536500
C      3.62426800   -1.48867000   -0.07067000
[略]
 
25 F

此例C--1.0代表对这个8号原子不手动指定原子类型（两个横杠之间没写内容），原子电荷为+1.0。没设原子电荷的原子的原子电荷默认为0。Gaussian里电场方向和习俗相反，因此field=x后面是负号。电场强度用300（0.03 a.u.）一般就可以，如果发现最终结果怪异，如有些键被拉得太长，结构都扭曲了，应当适当减小电场以减小外力；如果分子拉得还不够直，可以稍微增大电场来加大外力。nosymm必须写，要不然结构被自动弄到标准朝向下后外电场相对于分子的方向就和期望不符了，关于nosymm更多信息见《谈谈Gaussian中的对称性与nosymm关键词的使用》（//www.umsyar.com/297）。nomicro还是要写，要不然冻结设定不生效，其电场设置也不会对优化起作用。

用Gaussian计算上面的输入文件，68步后收敛，耗时才十几秒。最终结构如下，可见是我们想要的

3.3 方法3：令两端的原子彼此间受到静电斥力

这个做法超级简单，而且超级快速（直接用分子力场的专用优化代码），比前述两种方法通常更值得推荐。原理很简单，即给两端的原子都设数值较大的相同符号的原子电荷，二者通过强烈的静电斥力自然就会彼此远离，从而把分子拉直。输入文件如下，只有8号和25号原子设了原子电荷，都设为了+10.0。设多大合适应当看实际情况，比如如果分子特别长，那么原子电荷应设大一些，要不然静电斥力不够强。此例输入文件如下

%nprocs=1
# UFF opt

Generated by Multiwfn

  0  1
[略]
C      0.96378000   -3.24289700    4.74461900
C--10.0      1.95041300   -3.88295500    5.73536500
C      3.62426800   -1.48867000   -0.07067000
[略]
S     -2.51218200   -3.05838400    3.26014600
O--10.0     -4.37962400   -3.91727700   -0.05061200
O     -1.13156500   -5.32519500    2.86775900
[略]

用Gaussian计算，仅花了一秒钟就算完了！结果如下所示，非常平直，特别理想！

网友还给我发过一个含有200多个原子的弯曲的很长的分子，如下所示

用这一节的做法也很顺利地优化成了平直的构象，如下所示。计算仅仅花了5秒钟。

4 总结

本文介绍了笔者想到的四种把分子从弯曲结构拉成平直的方法，思路各有不同。实际当中会遇上各种特殊情况，比如可能需要同时拉直两条链（如磷脂的两条尾巴），应当在理解这些方法思想的基础上根据实际情况恰当选择、灵活运用。

本文的做法适合一般中、小分子体系。如果是大分子体系，比如蛋白质，比较适合GROMACS等专门的基于分子力场计算的程序。在GROMACS里面可以用pull相关的选项做拉伸动力学（见比如《在Gromacs中模拟金纳米线拉伸过程（含视频）》//www.umsyar.com/153），也可以将体系一端的原子定义为冻结组（freezegrps）而另一端定义为受常外力的组（acc-grps，结合accelerate设置）从而拉开。

值得一提的是，本文介绍的虽然是把分子拉直的方法，但实际上也可以将这些方法应用到研究环状、笼状、簇状等化学体系对外力造成的拉伸和压缩问题的研究上。

补充：使用Avogadro程序拉直分子的做法

通过Avogadro可视化程序也可以拉直分子，此程序可以在http://avogadro.cc免费下载。在此程序中拉直分子比较直观，直接按住鼠标拖拽即可，这和SAMSON可视化程序里的twist工具非常类似。启动Avogadro后，先载入pdb之类的结构文件，点击下图上方箭头所指的按钮，然后再点击下图左侧箭头所指的Start按钮，之后按住鼠标左键拖动特定原子即可，按照下图的绿色箭头分别把两端的原子往两头拉就能最终拉直。

Avogadro的这个功能实际上是持续通过指定的力场优化体系几何结构，因此当你拖动某个原子时，其它原子的位置会自发地弛豫来降低能量，因此键长、键角、二面角始终能保持比较合理的状态。但我发现Avogadro的这个功能拉直小分子不错，而拉直第3节最后那个特别长的聚合物则很难实现，而且特别卡顿。

A simple example of performing ab initio dynamics (AIMD) using ORCA

Sobereva@北京科音

First release: 2020-Dec-13  Last update: 2021-Jul-7

1 前言

基于 方法的动力学一般称为从头算动力学（Ab initio molecular dynamics, AIMD），相比于基于一般的经典力场的动力学，其关键优势在于精度高、普适性强、能够描述化学反应，代价是耗时相差N个数量级。ORCA 程序有不错的做BOMD形式的从头算动力学的功能，使用很方便，而且本身ORCA做DFT的效率又高，是做孤立体系AIMD的首选程序之一。虽然有些特性不支持，比如没法像Gaussian的BOMD那样直接做准经典动力学，不能根据原子距离等标准判断什么时候自动结束任务等等，但都不是大问题。对于跑跑普通的AIMD来说，笔者感觉ORCA明显比Gaussian的ADMP或BOMD更好用（Gaussian的AIMD输入文件较为抽象，手册相应部分写得很烂，而且连个像样的热浴都没有），而且速度也明显更快。

笔者发表的18碳环（cyclo[18]carbon）的研究论文中，其单体和二聚体做的分子动力学就是用ORCA跑的，分别见《揭示各种新奇的碳环体系的振动特征》（//www.umsyar.com/578）对中Chem. Asian J., 16, 56 (2021)和《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）中对Carbon, 171, 514 (2021)一文的介绍。在《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）中介绍的笔者的Inorg. Chem., 62, 19986 (2023)文章中还用ORCA跑了B6C6N6的高温动力学以证明其稳定性。在《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（//www.umsyar.com/696）中介绍的笔者的ChemPhysChem, 25, e202400377 (2024)文章中还用ORCA跑了18氮环的动力学以考察其热分解行为。这些文章可以作为ORCA跑AIMD的典型范例，很推荐阅读和引用。

笔者在北京科音高级 培训班（http://www.keinsci.com/workshop/KAQC_content.html）中会用多达两百多页的ppt专门深入详细讲AIMD的模拟，其中也包括ORCA做AIMD的各种细节、大量技巧以及诸多实例，并且此培训中还会系统讲授ORCA程序的使用，因此是使用ORCA专门做AIMD的读者一定不可错过的培训。而本文只是举一个简单的例子，帮助读者快速了解ORCA如何做最简单的AIMD。注意ORCA只适合跑孤立体系的从头算动力学，如果是做周期性体系的从头算动力学，CP2K是最佳的选择，在笔者讲授的北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）中有极为全面、系统、详细的讲解并给了大量例子。

本文内容适用于Multiwfn最新版本、VMD 1.9.3、ORCA 4.2.x及以后版本的情况。

2021-Jul-7 针对ORCA 5.0的补充说明：对于2021-Jul-7及以后更新的Multiwfn，进入本文所述的Multiwfn产生输入文件的界面后，可以通过选项-11选择适合的ORCA版本，默认为ORCA 5.0及以后版本，而下文内容对应的是4.2.x版的情况。对于>=5.0版，Multiwfn自动设的热浴是CSVR，比之前版本唯一支持的Berendsen热浴更好，而且同样普适。而且在run后面多加了CenterCOM，这是从5.0版本开始支持的消除整体质心运动的选项，而下文里提到的constraint add center这一行就没有了。

2 实例：[Al(H2O)6]3+与NH3之间的质子转移

ORCA的MD模块的开发者网站上有个真空中[Al(H2O)6]3+与NH3之间的质子转移的动画，如下所示

可见NH3向带正电的[Al(H2O)6]3+逐渐靠近，水上的一个质子转移到了氨气分子上，然后由于静电互斥，NH4+就逐渐远离[Al(H2O)5(OH)]2+了。这里我们试图重现这个过程。下面提到的文件都可以在//www.umsyar.com/attach/576/file.zip中获得。

我们首先获得[Al(H2O)6]3+的基本合理的结构。当然用ORCA优化也可以，这里笔者习惯性地用Gaussian来优化。在GaussView里搭建Al(H2O)6，保存为Al(H2O)6_optfreq.gjf，将关键词改为# B3LYP/6-31G* opt freq，将电荷改为3，然后用Gaussian运行之，就得到了优化后的[Al(H2O)6]3+结构。再用GaussView打开输出文件Al(H2O)6_optfreq.out，把一个氨气分子画在一个水的旁边，如下所示，然后保存为complex.gjf。

Multiwfn程序（//www.umsyar.com/multiwfn）有很便利的生成ORCA常见任务的输入文件的功能，见《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490），这里我们用Multiwfn生成ORCA的AIMD任务的标准输入文件。

启动Multiwfn，载入complex.gjf，然后输入
oi  // 生成ORCA输入文件
0  // 选择任务类型
6  // 分子动力学
1  // 计算级别用B97-3c
在当前目录下就得到了ORCA的AIMD任务的标准输入文件complex.inp。

在complex.inp里面将相应几行改成下面这样：
%maxcore  10000
%pal nprocs   10 end
timestep 1.0_fs
initvel 298.15_K
* xyz   3   1

现在的complex.inp的完整内容如下。以#为开头的行代表后面的设置被注释掉了，不会生效，想启用可以去掉#。此文件里也有大量Multiwfn自动添加的注释，只要一看注释马上就明白相应的行是什么含义，巨贴心，都省得查手册了。

! B97-3c noautostart miniprint nopop
%maxcore  10000
%pal nprocs   10 end
%md
#restart ifexists  # Continue MD by reading [basename].mdrestart if it exists. In this case "initvel" should be commented
#minimize  # Do minimization prior to MD simulation
 timestep 1.0_fs  # This stepsize is safe at several hundreds of Kelvin
 initvel 298.15_K no_overwrite # Assign velocity according to temperature for atoms whose velocities are not available
 thermostat berendsen 298.15_K timecon 30.0_fs  # Target temperature and coupling time constant
 dump position stride 1 format xyz filename "pos.xyz"  # Dump position every "stride" steps
#dump force stride 1 format xyz filename "force.xyz"  # Dump force every "stride" steps
#dump velocity stride 1 format xyz filename "vel.xyz"  # Dump velocity every "stride" steps
#dump gbw stride 20 filename "wfn"  # Dump wavefunction to "wfn[step].gbw" files every "stride" steps
 constraint add center 0..22  # Fix center of mass at initial position
 run 2000  # Number of MD steps
end
* xyz   3   1
[坐标部分]
*

下面简单说一下complex.inp里这些设置的含义、为什么这么设。
• B97-3c是一个又便宜又不错的计算级别，在ORCA里还自动会启用RIJ加速，速度很快，因此很适合跑AIMD，描述当前体系没有问题。B97-3c的介绍见《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490）的相应部分。但B97-3c也不是什么时候都能用，比如18碳环用纯泛函描述皆失败，见笔者在Carbon, 165, 468 (2020) 里的讨论，显然就不能用这个了，笔者跑涉及18碳环的AIMD的时候都用的是ωB97X-D3。
• %pal nprocs   10 end代表用10核。笔者当前任务实际上是在双路E5-2696v3共36个物理核心的机子上跑的，但却故意用了10核。这是因为根据笔者以前的测试，发现ORCA做DFT的AIMD的并行效率不理想，尤其是对于小体系，用核数太多甚至反倒速度更慢（大家可以对实际情况短暂跑比如5步实测一下设多少核速度最快）。一般来说就设10核就行了，当机子里有明显更多核的时候，可以跑多个AIMD任务来充分利用计算资源，但应当对CPU内核进行绑定，否则AIMD计算速度可能显著降低，见《通过设置CPU内核绑定降低ORCA同时做多任务的耗时》（//www.umsyar.com/553）。
• %maxcore  10000代表每个ORCA进程最多用10000MB。其实完全没必要这么大，普通泛函下的AIMD不怎么耗内存，给1000都绝对够了。
• restart ifexists这句被注释掉了。如果你的AIMD想续跑，且当前目录下有之前跑出来的与当前任务同名的后缀为.mdrestart的文件，可以取消注释，任务就会延续之前AIMD最后的状态续跑，新轨迹会在原有轨迹文件后面续写。
• minimize这句被注释掉了。如果取消注释的话，动力学开始前会自动在当前级别下做几何优化。
• timestep 1.0_fs代表动力学步长用1 fs。对于此例常温下的模拟，1 fs步长是合适的，改大的话可能造成动力学不稳定，而改小的话跑同样的时间长度需要更多步数导致需要更多耗时。如果追求绝对稳妥，或者是在明显更高温度下模拟，可以用0.5 fs。
• initvel 298.15_K代表根据298.15 K温度通过Maxwell分布随机初始化原子速度。no_overwrite代表如果之前已经有初速度信息了（比如可以是续跑时从之前的mdrestart文件里继承来的），就不再产生新的初速度。
• thermostat berendsen 298.15_K timecon 30.0_fs代表用Berendsen热浴控温在298.15 K（此例笔者用的ORCA 4.2.1只能用这个热浴，据悉从ORCA 5.0开始可以用更好的velocity rescale热浴），时间常数为30 fs，通常这个时间常数是合适的。
• dump position stride 1 format xyz filename "pos.xyz"代表每一步都往当前目录下的pos.xyz文件里写入当前的坐标，因此pos.xyz是多帧xyz格式的轨迹文件。此格式的介绍见《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477）。
• dump force和dump velocity开头的行都被注释掉了，这两行用于把模拟过程中的原子受力和原子速度分别保存到相应xyz文件里。dump gbw开头的行也被注释掉了，它可以在MD过程中每隔指定步数保存一次gbw文件，之后可以通过写脚本调用Multiwfn对它们进行批量分析，从而考察动力学过程中电子结构变化，得到丰富的有化学意义的信息（如动力学中的电荷转移情况、成键变化情况、电子定域性变化等等，在北京科音高级 培训班里笔者会给出很多这种例子）
• constraint add center 0..22代表将当前整个体系（0号到22号原子）的质心约束在初始位置。之所以这么做，是因为尽管ORCA在产生初速度时已经去掉了整体平移的速度分量，但实际模拟过程中由于数值问题，仍然往往会看到有明显的整体漂移的现象，因而有碍观察（在VMD里观看轨迹时还得做一下align才能消掉）。因此直接把质心位置约束住就没有这个问题了。
• run 2000代表总共跑2000步，当前步长是1 fs，即最多跑2 ps。当前这个模拟的目的是观察到质子转移，跑多长时间合适并不好预估，所以可以一次先跑2000步看看，若不够到时候再续跑。值得一提的是，至少对于笔者现在用的ORCA 4.2.1，我发现如果一次跑的步数很多，达到2000步左右之后，之后每一步的耗时会有逐渐的上升趋势，原因不明。因此我建议每次跑最多不宜超过3000步，如果需要跑更长的话，最好分多段来跑。

现在用ORCA运行complex.inp，模拟过程中可以看到每一步的实时情况：

Step是当前的步数，Time是当前的时间，t_SCF和t_Grad分别是算这一步的SCF过程和计算受力的耗时，二者加和就是这一步的总耗时。可见每一步耗时约6~7秒，乘以要跑的步数就可以估计跑完整个任务的耗时。后面还可以看到每一步的温度、动能、势能等信息。由于当前体系原子数很少，温度相对于热浴的参考温度波动大是很正常的事。而且由于用了热浴，所以总能量E_tot也有明显波动。

VMD是观看动力学轨迹最好的程序，可以在http://www.ks.uiuc.edu/Research/vmd/免费下载。建议在模拟过程中，隔一阵子就用VMD把pos.xyz载入进去，看看当前的动态行为如何、跑成什么结构了。跑到289步的时候，笔者看了一眼pos.xyz，发现质子不仅已经完全转移，而且NH4+都已经跑走一定距离了，所以就没必要继续跑了，故把ORCA直接杀掉了。

模拟过程中ORCA还产生其它一些文件。complex.md.log是ORCA的MD模块输出的信息，相当于整个输出文件中的子集。complex.mdrestart是用来续跑的文件，每一步都会往里面写入当前步的时间、坐标、速度、受力等信息。complex-md-ener.csv是把每一步的时间、耗时、温度、能量等信息以csv格式保存的文件。还有其它一些零碎的文件，通常不是一般用户关心的，这里就不说了，都可以放心删掉，留着也没用。

现在我们来看模拟轨迹。建议大家根据《VMD初始化文件(vmd.rc)我的推荐设置》（//www.umsyar.com/545）里说的修改VMD的初始化文件，从而添加自定义命令bt，这样在播放轨迹的时候对每一帧会自动重新判断成键关系。

启动VMD，载入pos.xyz（在本文文件包里已提供，共290帧）。然后在文本窗口输入bt，按回车，使得每帧都更新连接关系。然后再输入bw，按回车使得背景变为白色。选Graphics - Representation，将Drawing Method改为CPK。然后点击VMD界面右下角的三角播放动画，会看到如下结果。如果想在VMD图形窗口中显示出帧号或时间，看《使VMD播放轨迹时同步显示帧号》（//www.umsyar.com/13）。

可见，我们跑出来的动力学现象和本文一开始的那个官方动画几乎完全一致。都是两个反应物先接近，然后形成复合状态时质子在二者之间震荡，最后NH4+跑掉。

如果想把质子转移情况通过距离随时间的变化曲线方式呈献给读者，可以在点击VMD的三维图形窗口后，按键盘上的2键（之后按r可以恢复为默认的旋转模式），然后点击两个原子正中心，二者之间就会增加Bond label（默认是以白色字显示距离，在黑背景下才看得清楚），这里笔者把N和转移过去的质子之间增加了Bond label。然后进入Graphics - Labels，切换到Bonds，选Graph，点击Show preview复选框，然后点击1:H 1:N这项，就会看到距离变化已经显示在预览窗口了，如下所示，其中红色和蓝色竖条标记的分别是最小距离和最大距离位置和数值。

如果点击Graph按钮，就会把曲线显示在大窗口中，如下所示。可见，在大约60帧，也即60 fs左右，N-H键就算是基本形成了，之后N-H键不断振动。

以类似方式，我们可以标记Al与N的距离，随时间变化如下所示。可见Al与N先接近，质子转移完毕后，二者就逐渐远离了。

在Labels窗口里还可以点击save，把距离变化数据保存到文本文件里，之后可以导入Origin等程序里绘制曲线。

用VMD还可以测量角度、二面角的变化，分别是在图形窗口里按3然后点三个原子、按4然后点四个原子进行标记，之后在Graphics - Labels里观看。

3 总结&其它

通过上面的例子，可以看到ORCA做AIMD是相当容易的，只要把Multiwfn支持的含有结构信息的文件（如pdb/gjf/xyz/mol/mol2/fch等等，见Multiwfn手册2.5节）载入Multiwfn，敲几下键盘产生标准AIMD任务的输入文件，然后根据实际情况稍微改几个设定就可以跑了。

以几十核的一般双路服务器的运算能力，ORCA里用B97-3c跑几十原子有机体系的几十ps的动力学不是特别困难的事。不过，能跑的时间尺度仍远远比不上xtb跑半经验层面DFT的GFN-xTB方法的动力学，xtb跑动力学的粗浅介绍和例子看此文的相应部分：《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）。因此，拿ORCA跑DFT的动力学之前，先拿xtb初步跑跑，找找感觉，大体摸索出自己期望的现象能出现的条件（如温度、初始结构、反应物相对位置和碰撞方式等），然后再用DFT跑通常是比较好的做法，免得做昂贵的DFT的MD试来试去把时间都耽误了。

本文只涉及了VMD一丁点皮毛，VMD对于做动力学的人是必须玩得非常溜的。笔者在北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/workshop/KGMX_content.html）里对VMD有非常深入全面的介绍，包括tcl脚本的编写。利用VMD的tcl脚本可以对轨迹做千变万化的分析，有些分析诸如质心距离变化、平面间夹角变化、某几何变量分布统计、不同结构出现比率等，是必须靠写脚本才能实现的。

光是分析分析动力学过程的能量、结构变化是很肤浅的。利用Multiwfn，可以对ORCA跑的动力学的全过程的电子结构做深入透彻的分析，从而考察化学键、弱相互作用、电荷分布等等在动力学过程中的变化，由此能够从提供深入的视角，使研究文章信息更丰富、明显更上档次。非常建议详细阅读《详谈Multiwfn的命令行方式运行和批量运行的方法》（//www.umsyar.com/612），里面第4节专门讲了怎么做这样的分析，你会发现特别容易也特别有价值。

如果Multiwfn创建ORCA做动力学输入文件的功能对你的实际研究产生了帮助，希望在写文章的时候提及诸如The input file of ab-initio molecular dynamics was prepared with the help of Multiwfn program并引用Multiwfn原文，这是对Multiwfn此功能开发最好的支持。

Reducing time-consuming of multitasking of ORCA by setting CPU core binding

文/Sobereva@北京科音   2020-Jun-1

1 前言

由于任何程序的并行效率都是随着核数增加而降低，当机子核数比较多的时候，比如有好几十核，而且又有许多任务要跑的时候，比起一个一个调用所有核心来跑，同时跑两个或者多个任务但是每个都用较少的核心数（总和不超过物理核心数），总耗时通常会更低。对于某些程序跑某些任务，甚至并行核数较少的时候反倒比核数较多的时候速度还更快。因此在核数较多的机子上，同时跑多个任务是很常见的事情。

然而，同时跑多个任务涉及到资源争抢问题，如果争抢得比较厉害，跑多个任务的效率会大打折扣，甚至可能还不如一个一个用所有核来跑。

如今很多程序在并行计算的时候都支持线程或进程与CPU内核的绑定，从而避免操作系统自动调度导致任务在不同CPU核心上频繁切换执行引起运行效率的下降，恰当绑定也可以有效降低CPU资源争抢问题。ORCA自身没有直接提供内核绑定的设置，但可以通过自行指定OpenMPI的运行选项来实现，本文就说说具体做法和实际效果。

2 CPU内核绑定的实现

2.1 在ORCA运行时使用内核绑定

ORCA在Linux下并行一般结合OpenMPI使用，本文也只说OpenMPI的情况。OpenMPI基于hwloc库来获取硬件信息，而hwloc库又要利用libnuma。只有在编译OpenMPI时提供了libnuma静态库，之后ORCA运行时才能借助OpenMPI的设置实现内核绑定。多数情况下系统里是没有装libnuma静态库的，对于CentOS 7.x，可以通过yum install numactl-devel命令来装上，之后再按照《 程序ORCA的安装方法》（http://sobereva.com/451）里说的照常编译OpenMPI就可以了。

ORCA运行命令中可以加入传递给OpenMPI的mpirun的选项，比如：
orca test.inp "--bind-to core"
此时ORCA利用mpirun来运行其支持并行的模块的时候就会把--bind-to core选项传递给mpirun。OpenMPI 3.1.x版的mpirun支持的选项在这里都能看到：https://www.open-mpi.org/doc/v3.1/man1/mpirun.1.php，其中就包括了内核绑定的设置的说明。

2.2 内核绑定规则的设置

mpirun的选项中--bind-to [值]指定的是绑定的对象什么，而--map-by [值]是设定按照什么顺序循环被绑定物。例如当前机子是双路的，每个CPU有六个核，看以下两种情况：

--bind-to core --map-by core：按照循环各个核的顺序绑定，图例如下。B是代表运行当前进程的CPU核心
进程0：[B . . . . .][. . . . . .]
进程1：[. B . . . .][. . . . . .]
进程2：[. . B . . .][. . . . . .]
进程3：[. . . B . .][. . . . . .]
...略
这表明进程0绑定了第一个CPU的第一个核心，只有这个核心被用来跑这个进程。

--bind-to socket --map-by socket：按照循环各个CPU的顺序绑定（socket此处是指CPU插槽），图例如下
进程0：[B B B B B B][. . . . . .]
进程1：[. . . . . .][B B B B B B]
进程2：[B B B B B B][. . . . . .]
进程3：[. . . . . .][B B B B B B]
...略
可见，此时进程0绑定了第一个CPU的所有核心，这些核心一起执行这个进程，而第二个CPU核心不会跑这个进程。

mpirun加上--report-bindings选项可以在开始并行执行时显示当前是怎么绑定的，用以确认绑定方式合乎自己的要求。

如果想精细控制绑定规则，需要创建rankfile文件。比如rankfile文件是/sob/miku.txt，可以这么运行./violet.x程序：mpiexec -np 3 -rf /sob/miku.txt --report-bindings ./violet.x，此时一共3个进程会按照rankfile里指定的方式来绑定。

比如rankfile文件内容如下
rank 0=Saber110 slot=1:0,1
rank 1=Saber109 slot=0:*
rank 2=Saber108 slot=1:1-3
rank 3=Saber17 slot=0:1,1:0-2
rank 4=Saber109 slot=0:*,1:*
rank 5=Saber109 slot=0-2
就意味着（注意此处进程、CPU编号和核心编号都是从0开始的）
第0个进程与Saber110节点的第1号CPU的0、1两个核心绑定
第1个进程与Saber109节点的第0号CPU的所有核心绑定
第2个进程与Saber108节点的第1号CPU的1、2、3核心绑定
第3个进程与Saber17节点的第0号CPU的1号核心，以及第1号CPU的0~2号核心绑定
第4个进程与Saber109节点的第0、1号CPU的所有核心绑定
第5个进程与Saber109节点的第0、1、2号逻辑核心绑定
如果比如当前是-np 3跑3个进程，则只有前3个规则生效，而rank 3、rank 4、rank 5的设置不会被利用。

2.3 让每个ORCA任务独占不同的CPU核心

如果只跑一个任务的话，设内核绑定没什么用处。而跑多个任务时，给不同任务指定不同的CPU核心可用范围，则有可能降低CPU资源争抢来降低耗时。为此，需要设置不同rankfile。比如双路共36核72线程机子，其逻辑核心的序号是下面的顺序：
第一个CPU的18个物理核心：依次为0~17
第二个CPU的18个物理核心：依次为18~35
第一个CPU的18个虚拟核心：依次为36~53
第二个CPU的18个虚拟核心：依次为54~71
更确切来说，比如0和36号逻辑核心的任务其实都是第一个CPU的实际的0号核心做的，如果就是一个MPI进程要跑，在0还是36号上执行速度都一样，但如果有两个MPI进程，若在0和36号上面同时跑，速度就只有之前的约一半了。详见《正确认识超线程(HT)技术对计算化学运算的影响》（http://sobereva.com/392）。下文就忽略掉由于超线程技术而多出来的36~71号核心了。

如果有两个ORCA任务要跑，每个都设nprocs 18，那么默认的情况下，每个任务会同时在两个CPU上跑，有可能存在资源争抢。如果希望这两个任务分别在第一和第二个CPU上跑，可以设置这以下两个rankfile文件。
CPU1.txt，内容为
    rank 0=2696v3 slot=0:0
    rank 1=2696v3 slot=0:1
    rank 2=2696v3 slot=0:2
    rank 3=2696v3 slot=0:3
    rank 4=2696v3 slot=0:4
    rank 5=2696v3 slot=0:5
    rank 6=2696v3 slot=0:6
    rank 7=2696v3 slot=0:7
    rank 8=2696v3 slot=0:8
    rank 9=2696v3 slot=0:9
    rank 10=2696v3 slot=0:10
    rank 11=2696v3 slot=0:11
    rank 12=2696v3 slot=0:12
    rank 13=2696v3 slot=0:13
    rank 14=2696v3 slot=0:14
    rank 15=2696v3 slot=0:15
    rank 16=2696v3 slot=0:16
    rank 17=2696v3 slot=0:17
CPU2.txt，内容为
    rank 0=2696v3 slot=1:0
    rank 1=2696v3 slot=1:1
    rank 2=2696v3 slot=1:2
    rank 3=2696v3 slot=1:3
    rank 4=2696v3 slot=1:4
    rank 5=2696v3 slot=1:5
    rank 6=2696v3 slot=1:6
    rank 7=2696v3 slot=1:7
    rank 8=2696v3 slot=1:8
    rank 9=2696v3 slot=1:9
    rank 10=2696v3 slot=1:10
    rank 11=2696v3 slot=1:11
    rank 12=2696v3 slot=1:12
    rank 13=2696v3 slot=1:13
    rank 14=2696v3 slot=1:14
    rank 15=2696v3 slot=1:15
    rank 16=2696v3 slot=1:16
    rank 17=2696v3 slot=1:17
这里2696v3是当前的主机名，运行hostname命令时显示什么就填什么（如果显示的是IP号，这里就填IP号）。

假设这两个rankfile文件都在/sob目录下，然后执行以下两条命令分别跑ORCA任务
orca MD1.inp "-rf /sob/CPU1.txt" > MD1.out
orca MD2.inp "-rf /sob/CPU2.txt" > MD2.out
则MD1和MD2任务就分别只在第一个和第二个CPU上跑了

如果你有四个任务要跑，每个任务给9个核，那么可以创建四个rankfile：
CPU1a.txt的内容
    rank 0=2696v3 slot=0:0
    rank 1=2696v3 slot=0:1
    rank 2=2696v3 slot=0:2
    rank 3=2696v3 slot=0:3
    rank 4=2696v3 slot=0:4
    rank 5=2696v3 slot=0:5
    rank 6=2696v3 slot=0:6
    rank 7=2696v3 slot=0:7
    rank 8=2696v3 slot=0:8
CPU1b.txt的内容
    rank 0=2696v3 slot=0:9
    rank 1=2696v3 slot=0:10
    rank 2=2696v3 slot=0:11
    rank 3=2696v3 slot=0:12
    rank 4=2696v3 slot=0:13
    rank 5=2696v3 slot=0:14
    rank 6=2696v3 slot=0:15
    rank 7=2696v3 slot=0:16
    rank 8=2696v3 slot=0:17
CPU2a.txt的内容
    rank 0=2696v3 slot=1:0
    rank 1=2696v3 slot=1:1
    rank 2=2696v3 slot=1:2
    rank 3=2696v3 slot=1:3
    rank 4=2696v3 slot=1:4
    rank 5=2696v3 slot=1:5
    rank 6=2696v3 slot=1:6
    rank 7=2696v3 slot=1:7
    rank 8=2696v3 slot=1:8
CPU2b.txt的内容
    rank 0=2696v3 slot=1:9
    rank 1=2696v3 slot=1:10
    rank 2=2696v3 slot=1:11
    rank 3=2696v3 slot=1:12
    rank 4=2696v3 slot=1:13
    rank 5=2696v3 slot=1:14
    rank 6=2696v3 slot=1:15
    rank 7=2696v3 slot=1:16
    rank 8=2696v3 slot=1:17
之后跑四个ORCA任务的时候分别指定这四个rankfile，则第一、二个任务就分别在第一个CPU的前9个、后9个核心上跑，第三、四个任务就分别在第二个CPU的前9个、后9个核心上跑，这尽可能减小了四个任务彼此的干扰（但干扰是不可能完全杜绝的）。

3 绑定不同CPU核心跑多任务对耗时的影响

为了测试绑定不同CPU核心对于跑多任务的耗时影响，我选了一个碱基对当测试体系，一共29个原子，在B3LYP-D3(BJ)/def2-QZVP结合RIJCOSX级别下，用笔者的XEON 2696v3双路机子（36核72线程）进行测试。结果如下

其中诸如nprocs=18代表不对内核进行绑定设置的情况。可见跑两个18核任务的时候，是否绑定对耗时影响可以忽略。跑两个9核的任务时，分别在两个CPU上跑的耗时明显比让这俩任务分别占同一个CPU的前9个和后9个核心的耗时要低，显然后者的情况会对第一个CPU有一定程度的资源争抢。内核绑定的效果对于同时跑四个任务很明显，直接跑四个9核的任务的耗时比让这四个任务绑定不同核心（分别绑定0-8、9-17、18-26、27-35号核心）的耗时高了1/3有余！

这种绑定的做法对于不同类型任务产生的效果不一样。笔者之前拿ORCA跑了一个环状体系的AIMD模拟，哪怕只是同时跑两条对应于不同温度下的轨迹，每个都给9核，分别绑在两个CPU上跑的速度都显著快于不设绑定直接提交两个任务的情况。还可以观察到一个现象，如果以普通方式运行，当第二个AIMD任务交上去之后，就会发现ORCA显示的第一个AIMD任务的每步的耗时增加了很多；而如果两个任务在运行时分别绑在不同CPU的核心上，则提交第二个任务后，之前正在跑的第一个任务每步的耗时几乎没有变化。

总之，对于ORCA用户，如果经常同时跑多任务的话，不妨留意一下内核绑定问题。对于跑AIMD，我个人总是将多个任务绑定不同内核来跑，哪怕只是同时跑两个任务。

Gaussian从16版开始支持%cpu设置，可直接控制任务在哪些CPU核心上运行，这在《正确认识超线程(HT)技术对计算化学运算的影响》（http://sobereva.com/392）中我专门提过。笔者也测试了一下在Gaussian中把同时提交的任务指定在不同CPU核心上对速度的影响。还是用碱基对体系，在M06-2X/def2-TZVP下做单点计算，耗时如下

可见是否将不同任务绑在不同CPU核心上运行对耗时基本没影响，所以Gaussian用户不需要特意考虑内核绑定问题。如果两个9线程的任务分别绑定在同一个CPU的前9个和后9个核来跑，速度则会明显不如照常提交这两个任务，这和测试ORCA时看到的情况是一致的。

Example of performing counterpoise correction and calculation of intermolecular binding energy in ORCA

文/Sobereva@北京科音

First release: 2020-Mar-13  Last update: 2022-Jan-9

在《谈谈BSSE校正与Gaussian对它的处理》（//www.umsyar.com/46）一文中介绍过BSSE问题与counterpoise校正。Counterpoise校正在高精度弱相互作用计算时通常需要考虑。Gaussian有现成的counterpoise关键词，因此做counterpoise校正很简单。但是ORCA程序（至少对于撰文时最新的4.2.1版来说）没有提供现成的counterpoise关键词，只能通过定义鬼原子的方式手工实现，因此实现起来稍微麻烦。不熟悉ORCA者看//www.umsyar.com右侧ORCA分类里我写过的大量相关文章。

之前在《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490）中笔者专门介绍了Multiwfn产生ORCA 程序的输入文件的功能。为了使得用ORCA做counterpoise任务更方便，笔者2020-Mar-13给Multiwfn的这个功能增加了相应选项，在这里通过一个实例简单演示一下。最新版Multiwfn可以在//www.umsyar.com/multiwfn免费下载，相关知识见《Multiwfn FAQ》（//www.umsyar.com/452）。

注意按照下文例子算出来的实际上是复合物结构中单体间相互作用能。而一般意义的结合能在计算时是要考虑单体形成复合物过程中形变造成的能量改变（变形能）。由于下面的例子里单体是刚性的，变形能可以忽略不计，因此下文仍称作结合能。

本文用的示例体系是下图这个体系，我们要在DLPNO-CCSD(T)/aug-cc-pVTZ下计算其结合能，并考虑counterpoise校正。

此体系的xyz格式的结构文件，后文中涉及的输入输出文件以及要用的shell脚本都可以在这里下载：//www.umsyar.com/attach/542/file.rar。并不是必须得用xyz格式，只要是Multiwfn支持的含有结构信息的格式（如pdb、mol2...）都可以作为输入文件，见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。

启动Multiwfn，输入dimer.xyz的路径载入之。如果你还不知道两个单体的原子序号范围，可以进入Multiwfn的主功能0来查询。进入主功能0后，选择Tools - Select fragment，输入其中一个片段中的任意一个原子的序号比如16，整个片段就被高亮，而且序号也显示出来了，如下所示。

如此这般，把两个片段序号范围都拷出来备用，然后点Return退出主功能0，之后依次输入以下内容：
oi  //产生ORCA的输入文件
[按回车] //输出文件与输入文件同名但用inp作为后缀
0  //选择任务类型
7  //考虑counterpoise校正的结合能
1-15  //片段1的序号
16-32  //片段2的序号
-2  //加弥散函数
8  //用DLPNO-CCSD(T) with NormalPNO

在当前目录下马上就有了dimer.inp，内容如下

%pal nprocs   4 end

! DLPNO-CCSD(T) normalPNO RIJK aug-cc-pVTZ aug-cc-pVTZ/JK aug-cc-pVTZ/C tightSCF noautostart miniprint nopop
%maxcore 1000
* xyz   0   1
C      0.79991408   -1.02205164    0.68773696
H      0.85355588   -1.12205101   -0.39801435
H      1.49140210   -1.74416936    1.11972040
C      1.11688700    0.42495279    1.09966205
H      1.83814230    0.89014504    0.43045256
H      1.55556959    0.43982464    2.09708356
C     -0.24455916    1.16568959    1.10297714
H     -0.25807760    2.00086313    0.40532333
H     -0.44880450    1.57699582    2.09098447
C     -1.29871418    0.10381191    0.73930899
H     -1.47356078    0.10524338   -0.33800545
H     -2.25673428    0.27804118    1.22715843
C     -0.64687993   -1.22006836    1.13630660
H     -1.12443918   -2.08762702    0.68299327
H     -0.68601864   -1.34528332    2.22022006
C      0.04984615    0.09420760    5.61627735
C     -0.04649805   -0.05787837    7.13191782
H      0.94604832   -0.07334458    7.58427505
H     -0.60542282    0.77000613    7.57035274
H     -0.55366275   -0.98654445    7.39726741
C      0.76389939    1.40111272    5.28065247
H      0.84541894    1.53461185    4.20097059
H      0.22042700    2.25580115    5.68615385
H      1.77150393    1.41176313    5.69888547
C     -1.35516567    0.11403225    5.01895782
H     -1.31823408    0.23122219    3.93510886
H     -1.93746520    0.94145581    5.42730374
H     -1.88506873   -0.81375459    5.24028712
C      0.83774596   -1.07927730    5.03893917
H      0.34252564   -2.02626804    5.25918232
H      0.93258913   -0.99209454    3.95580439
H      1.84246405   -1.11668194    5.46268763
 *

$new_job
! DLPNO-CCSD(T) normalPNO RIJK aug-cc-pVTZ aug-cc-pVTZ/JK aug-cc-pVTZ/C tightSCF noautostart miniprint nopop Pmodel
%maxcore 1000
* xyz   0   1
C      0.79991408   -1.02205164    0.68773696
H      0.85355588   -1.12205101   -0.39801435
H      1.49140210   -1.74416936    1.11972040
C      1.11688700    0.42495279    1.09966205
H      1.83814230    0.89014504    0.43045256
H      1.55556959    0.43982464    2.09708356
C     -0.24455916    1.16568959    1.10297714
H     -0.25807760    2.00086313    0.40532333
H     -0.44880450    1.57699582    2.09098447
C     -1.29871418    0.10381191    0.73930899
H     -1.47356078    0.10524338   -0.33800545
H     -2.25673428    0.27804118    1.22715843
C     -0.64687993   -1.22006836    1.13630660
H     -1.12443918   -2.08762702    0.68299327
H     -0.68601864   -1.34528332    2.22022006
C:     0.04984615    0.09420760    5.61627735
C:    -0.04649805   -0.05787837    7.13191782
H:     0.94604832   -0.07334458    7.58427505
H:    -0.60542282    0.77000613    7.57035274
H:    -0.55366275   -0.98654445    7.39726741
C:     0.76389939    1.40111272    5.28065247
H:     0.84541894    1.53461185    4.20097059
H:     0.22042700    2.25580115    5.68615385
H:     1.77150393    1.41176313    5.69888547
C:    -1.35516567    0.11403225    5.01895782
H:    -1.31823408    0.23122219    3.93510886
H:    -1.93746520    0.94145581    5.42730374
H:    -1.88506873   -0.81375459    5.24028712
C:     0.83774596   -1.07927730    5.03893917
H:     0.34252564   -2.02626804    5.25918232
H:     0.93258913   -0.99209454    3.95580439
H:     1.84246405   -1.11668194    5.46268763
 *

$new_job
! DLPNO-CCSD(T) normalPNO RIJK aug-cc-pVTZ aug-cc-pVTZ/JK aug-cc-pVTZ/C tightSCF noautostart miniprint nopop Pmodel
%maxcore 1000
* xyz   0   1
C:     0.79991408   -1.02205164    0.68773696
H:     0.85355588   -1.12205101   -0.39801435
H:     1.49140210   -1.74416936    1.11972040
C:     1.11688700    0.42495279    1.09966205
H:     1.83814230    0.89014504    0.43045256
H:     1.55556959    0.43982464    2.09708356
C:    -0.24455916    1.16568959    1.10297714
H:    -0.25807760    2.00086313    0.40532333
H:    -0.44880450    1.57699582    2.09098447
C:    -1.29871418    0.10381191    0.73930899
H:    -1.47356078    0.10524338   -0.33800545
H:    -2.25673428    0.27804118    1.22715843
C:    -0.64687993   -1.22006836    1.13630660
H:    -1.12443918   -2.08762702    0.68299327
H:    -0.68601864   -1.34528332    2.22022006
C      0.04984615    0.09420760    5.61627735
C     -0.04649805   -0.05787837    7.13191782
H      0.94604832   -0.07334458    7.58427505
H     -0.60542282    0.77000613    7.57035274
H     -0.55366275   -0.98654445    7.39726741
C      0.76389939    1.40111272    5.28065247
H      0.84541894    1.53461185    4.20097059
H      0.22042700    2.25580115    5.68615385
H      1.77150393    1.41176313    5.69888547
C     -1.35516567    0.11403225    5.01895782
H     -1.31823408    0.23122219    3.93510886
H     -1.93746520    0.94145581    5.42730374
H     -1.88506873   -0.81375459    5.24028712
C      0.83774596   -1.07927730    5.03893917
H      0.34252564   -2.02626804    5.25918232
H      0.93258913   -0.99209454    3.95580439
H      1.84246405   -1.11668194    5.46268763
 *

这个输入文件包含了DLPNO-CCSD(T)/aug-cc-pVTZ with normalPNO级别下的三个单点任务。第一个子任务是对二聚体做计算，输出的能量我们叫E_AB；第二个子任务是对第1个片段在二聚体基组下做计算，输出的能量我们叫E_A(AB)；第三个子任务是对第2个片段在二聚体基组下做计算，输出的能量我们叫E_B(AB)。原子名后面带冒号的说明这个原子是鬼原子。计算前别忘了把maxcore和nproc都设为恰当情况。其中nproc只要在开头定义一次就够了。

在笔者的2*Intel E5-2696v3（36核）机子上，通过36核并行，每个核给6000MB内存，花了41分钟跑完。在输出文件中搜索FINAL SINGLE POINT ENERGY，总共会搜索到三个：
FINAL SINGLE POINT ENERGY      -393.613585294163
FINAL SINGLE POINT ENERGY      -196.197449643018
FINAL SINGLE POINT ENERGY      -197.412738820804
依次为E_AB、E_A(AB)、E_B(AB)。然后按照计算E_AB - E_A(AB) - E_B(AB)就得到考虑counterpoise校正的结合能了。如果不理解这点，看我的培训班（http://www.keinsci.com/workshop/KBQC_content.html）里的一页ppt

对于当前例子，counterpoise校正的结合能因此为627.51*(-393.613585294163+196.197449643018+197.412738820804)= -2.13 kcal/mol。

本文的这个体系其实是著名的S66弱相互作用测试集里的，用的结构也是其补充材料里的结构。在S66文中通过CCSD(T)/CBS计算出来的结果是-2.40 kcal/mol，可见本文的结果和金标准符合得非常好。

为了让结合能计算更省事，我还专门写了个脚本，是本文文件包里的getEbind.sh。此脚本在运行后会处理当前目录下的所有out文件，自动提取整体和两个片段的能量，并求差给出不同单位的结合能。例如当前目录下有三个.out文件，都是上文的方法用Multiwfn产生的任务的输出文件，getEbind.sh输出的信息如下，可见提取结合能非常方便。

H2O_ext.out
Etot = -761.79953946  E1 = -685.35948656  E2 = -76.43844639 Hartree
-0.0016065 Hartree
    -1.008 kcal/mol
    -4.218 kJ/mol

H2O_int.out
Etot = -761.80467188  E1 = -685.35943188  E2 = -76.43928426 Hartree
-0.0059557 Hartree
    -3.737 kcal/mol
   -15.637 kJ/mol

H2.out
Etot = -686.52275741  E1 = -685.35832899  E2 = -1.16228297 Hartree
-0.0021455 Hartree
    -1.346 kcal/mol
    -5.633 kJ/mol

2022-Jan-9补充

从2022-Jan-9更新的Multiwfn开始，在产生ORCA输入文件的界面里选择任务类型时可以选择-7，用法同前。此时产生的ORCA输入文件会计算5个能量，依次为：整体的能量、整体基函数下片段1的能量、整体基函数下片段2的能量、片段1的能量、片段2的能量。执行本文文件包里的ORCA_CP.sh就可以对当前目录下所有这种任务的输出文件进行处理，此脚本对每个文件会输出这些信息：(1)5个能量的具体数值 (2)原始的相互作用能 (3)BSSE校正后的相互作用能 (4)BSSE校正能 (5)BSSE校正后的整体的能量。这种任务所做的计算、脚本给出的信息，和Gaussian程序的counterpoise关键词完全相同。

本文文件包里的waterdimer.out是一个例子文件，ORCA_CP.sh对它进行处理的输出如下：

File: waterdimer.out
 E(AB) = -152.054002577 Hartree
 E(A)_AB = -76.022655799    E(A) = -76.022397044 Hartree
 E(B)_AB = -76.024059371    E(B) = -76.022540157 Hartree
 Raw interaction energy:
-0.0090654 Hartree
    -5.689 kcal/mol
   -23.801 kJ/mol
 BSSE corrected interaction energy:
-0.0072874 Hartree
    -4.573 kcal/mol
   -19.133 kJ/mol
 BSSE correction energy:
 0.0017780 Hartree
     1.116 kcal/mol
     4.668 kJ/mol
 BSSE corrected complex energy:
-152.0522246 Hartree

Calculation of metal surface adsorption problems based on cluster models using quantum chemistry programs

文/Sobereva@北京科音

 First release: 2020-Mar-8   Last update: 2020-Mar-10

简介：本文首先对 程序基于簇模型计算表面反应或吸附问题做简单介绍，然后基于ORCA 程序，简单演示如何计算苯分子在Ag(111)表面上的吸附能。本文目的是令读者认识到 程序研究表面问题的可行性，并了解实际计算中需要值得注意的地方。

笔者还另有一篇博文章《18碳环（cyclo[18]carbon）与石墨烯的相互作用：基于簇模型的研究一例》（//www.umsyar.com/615），也和本文主题密切相关，建议读者一看。

1 关于簇模型

一说到基于量子力学方法的固体表面反应/吸附的计算，很多人就首先想到用第一性原理程序诸如Quantum Espresso、CP2K去做，甚至以为以计算孤立体系为主的Gaussian、ORCA等 程序完全不适合，实际上这是不对的。很多固体表面问题的计算都可以通过构建簇模型（cluster model）来使用 程序计算。所谓簇模型就是把反应/吸附位点附近区域的原子挖出来当做孤立体系计算。

用 程序通过簇模型来算和使用第一性程序通过周期性方式来算这类问题各有利弊。用簇模型的优点在于：
• 整体来说， 程序支持的理论方法远比第一性原理程序丰富得多得多，选择余地明显更大，还可以做到更高精度（如双杂化、DLPNO-CCSD(T)等）。用杂化泛函相对于纯泛函的耗时增加程度远低于基于平面波的第一性原理程序。
• 在计算设置上更方便，也不需要考虑k点、真空层、分子与相邻镜像作用、偶极校正之类乱七八糟的事。
• 支持的任务类型更丰富，关键词写着省事，可视化程序也多。例如Gaussian vs. VASP
• 结合Multiwfn可以做非常丰富的波函数分析，使文章的分析讨论部分明显更充实、更上档次。参考比如《Multiwfn支持的分析化学键的方法一览》（//www.umsyar.com/471）、《Multiwfn支持的弱相互作用的分析方法概览》（//www.umsyar.com/252）、《Multiwfn支持的电子激发分析方法一览》（//www.umsyar.com/437）等大量文章，在《Multiwfn入门tips》（//www.umsyar.com/167）里有所有Multiwfn做分析的博文汇总。而使用第一性原理程序的话，能做的分析少得多，不过《使用Multiwfn结合CP2K通过NCI和IGM方法图形化考察固体和表面的弱相互作用》（//www.umsyar.com/588）

用簇模型的弊端在于：
• 需要合理考虑边界效应。边界效应此处指的是使用有限的簇模型作为周期性体系的近似导致的对被研究问题的计算精度的影响。边界效应取决于两点：(1)簇的尺寸 (2)对边界的处理。簇的尺寸越大、让边界原子所处的环境越像体相，则边界效应越小。
簇的大小选取有很大任意性。固体部分截的原子太少的话边界效应太强，结果不准；而截得大的话，则计算耗时太高。除非团簇截得极大，否则都要考虑怎么恰当处理边界。比如有机体系（如石墨烯），边界往往用氢饱和，见比如《18碳环（cyclo[18]carbon）与石墨烯的相互作用：基于簇模型的研究一例》（//www.umsyar.com/615）；如果是无机晶体，有用赝氢、capped ECP、用大量背景电荷表现更长程区域原子的静电势等做法；如果是纯金属的话倒是不需要做特别的考虑。
• 对于原子致密排布的块状固体材料（不是那种稀疏的、有孔洞的诸如MOF等，或者单层材料），用簇模型的话耗时比起用第一性原理程序当周期性算往往更高。但这点具体看用什么程序、什么基组、数值方面的设定等等，不能一概而论。
• 用簇模型算含有过渡金属的无机晶体有时会在SCF收敛上遇到麻烦，有时需要大量折腾。但也不是绝对解决不了，可结合《解决SCF不收敛问题的方法》（//www.umsyar.com/61）里的说明恰当尝试。

总的来说，用第一性原理程序做周期性计算研究表面问题是比较主流的做法，文章非常多，但簇模型计算这类问题也绝非很稀罕，相关研究文章也不少。还有的研究将两种方式相结合，取二者长处，比如ACS Catal., 8, 3825 (2018)研究石墨烯上催化二氧化硫与氧气反应形成硫酸盐，就先用了VASP做了周期性计算优化了极小点和过渡态，再用Gausisan基于VASP优化的结构用簇模型做了单点计算，将得到的波函数用Multiwfn分析了催化反应机理以及石墨烯表面与被催化物质的弱相互作用。

2 簇模型计算表面问题的实际例子：苯分子在Ag(111)表面物理吸附的结合能的计算

下面就通过一个简单例子，说明如何基于簇模型计算物理吸附的结合能。

2.1 相关说明

本例涉及的各种文件都可以在这里下载：//www.umsyar.com/attach/540/file.rar。

这个体系在DFT-D3原文J. Chem. Phys., 132, 154104 (2010)中被研究过，苯与银表面间距的刚性扫描图如下所示

文中也说了，实验测定的苯吸附到Ag(111)上的焓变是-13 kcal/mol。

下面我们来算苯分子在Ag(111)表面物理吸附对应的结合能。这里说的结合能是习俗上定义的，即只考虑电子能量的变化，而且片段能量用的是复合物中的结构算的。算这个问题分为以下几步
(1)构建金属簇
(2)把苯分子放进去得到复合物初猜结构
(3)优化复合物
(4)计算各个片段的单点能
(5)求差得到结合能
之后还可以再用Multiwfn做一些波函数分析。

严格来说，就算哪怕不计算热力学量，优化完复合物也应当做一下振动分析来看看虚频情况，如果苯分子有沿表面运动的虚频模式，应当恰当消掉，基本策略参考《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》（//www.umsyar.com/278）。但由于当前体系不小，做振动分析耗时很高，本文就忽略这点了。

本例使用ORCA程序，一方面是因为免费，另一方面是支持RIJCOSX，可以令这种大团簇体系在单点计算和优化时耗时降低极多，明显快过Gaussian。不过ORCA做振动分析慢，也是本文没做振动分析的原因。ORCA的基本常识看《 程序ORCA的安装方法》（//www.umsyar.com/451）和《基于ORCA 程序对分子做优化、振动分析、观看红外光谱、观看轨道的简单演示》（https://www.bilibili.com/video/av59599938）。本文用的是ORCA 4.2.1。

2.2 建模

本文用的对应Ag(111)表面的团簇结构就按照上一节文献里的图来构建就行了，用什么程序随意。比如可以去ICSD找到Ag的晶体结构，用VESTA、GaussView、Avogadro之类复制延展成足够大的复晶胞，然后一点点删原子就可以得到。你也可以用某贵得离谱的有切表面功能的可视化程序先载入自带的库里的Ag原胞，用切表面的选项切出(111)表面，然后再复制延展成面积足够大、至少有三层原子厚度的表面，之后再一点点删多余的原子。为了不让边界效应太显著，对于簇模型研究表面问题，一般建议至少三层原子，而且固体团簇要比被吸附的分子明显大一圈。像上面图片里构造的Ag58团簇对于研究苯的吸附就是很得当的，不大不小。再大一圈的话耗时就忒高了，再小一圈的话边界效应就较为显著了。而且为了节约原子来降低耗时，可见上图中的团簇是上面宽下面窄，比起上下一样宽时原子数明显更少。由于苯分子竖着看起来是近乎圆形的，因此团簇的表面最好也接近圆形。显然，分子如果是长条形的，簇的表面也应当是长条形。

笔者建立好的簇如下所示，为了看得清楚每层用了不同颜色。在gview里可能默认没有把Ag-Ag之间判断为成键，导致构建簇模型的时候不好判断原子位置，可以根据《谈谈原子间是否成键的判断问题》（//www.umsyar.com/414）文中所述自定义成键判断阈值，笔者设的是小于3埃就判断成成键。

然后用把苯分子放在金属簇的上方正中央位置。距离不要太远，要不然优化到平衡位置需要多花很多步；也不要太近，本来过渡金属团簇就是SCF难收敛的典型，如果被吸附分子放得太近导致电子结构变得更复杂的话会雪上加霜。本文把苯分子放在团簇表面上方3埃的位置。从DFT-D3原文里那张势能面扫描图来看，TPSS-D3的极小点位置也差不多就是3埃。

当前的结构不能直接就拿去优化，因为簇模型里边界的原子没有感受到真实体相环境的束缚，如果不把边界的原子进行固定就优化的话，可能导致边界的原子位置大幅改变，和真实情况明显不符。因此我们把最靠边的Ag都冻结，而中间的，即与苯分子密切接触的Ag原子不设冻结，这是为了让它们的位置自发地调整来响应相互作用造成的实际的结构变化。

下图就是笔者在GaussView里最终构建的复合物初始结构，为了看得清楚给了两个视角，黄色是被冻结的原子。在GaussView的Tools - Atom Selection界面里可以直接查询到这些被选成黄色的原子的序号，即1-11,14,17,19-23,26,29,32,34-38,41,44,46-50,52,54-58。

这个初始结构用GaussView保存出来的gjf文件是本文文件包里的Ag58_benzene_freeze.gjf。

值得一提的是，但凡若有可能，都最好让簇模型处于闭壳层状态，这样算起来比开壳层的情况快得多。对于当前的体系，苯分子是闭壳层的，Ag是奇数个电子而且靠s电子成键（成键情况比较简单），故只要Ag原子是偶数个，就可以让整个体系是闭壳层。当前Ag是58个，所以整个体系是闭壳层。如果你用的银团簇是比如Ag59或者Ag57，那么整体就是开壳层了，算起来就费劲多了。所以在设计金属团簇的时候应当注意考虑这点。

2.3 几何优化

下面用Multiwfn产生ORCA程序对上面的体系做限制性优化的输入文件。相关信息看《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490）。本文用的是2020-Mar-5更新的Multiwfn，可在官网//www.umsyar.com/multiwfn下载。启动Multiwfn，然后输入
Ag58_benzene_freeze.gjf  //gjf文件的实际路径，Multiwfn会从中读取结构
oi  //产生ORCA输入文件
complex_opt.inp  //产生的文件的名字
0  //设置任务类型
2  //优化
-3  //设置原子冻结
1-11,14,17,19-23,26,29,32,34-38,41,44,46-50,52,54-58  //被冻结的原子序号
3  //用RI-B3LYP-D3(BJ)/def2-TZVP(-f)

然后我们就在当前目录下得到了complex_opt.inp文件。把其中的B3LYP改成PBE0，把def2-TZVP(-f)改为def2-SV(P)，把CPU核数和内存分配量设为实际情况，其它内容不变。此时这个文件对应于PBE0-D3(BJ)/def2-SV(P)下进行限制性优化，是我们实际要用的。

这里解释下为什么用这个级别做优化。前面我的博文里提到过，在ORCA里开RI时，纯泛函的耗时远低于杂化泛函，而且纯泛函PBE对纯金属的描述也不错，而之所以这里非要用杂化泛函PBE0，是因为用纯泛函时SCF收敛难度比杂化泛函通常明显更大，尤其是当前这种过渡金属团簇类型的体系本来就是SCF难收敛的体系，尝试过PBE就会发现收敛情况就像噩梦，所以用了杂化泛函。之所以把原本输入文件里的B3LYP关键词给改成了PBE0，这是因为J. Chem. Phys., 127, 024103 (2007)中专门说了B3LYP算金属很烂，而 领域里另一个很常用的杂化泛函PBE0则在这方面可靠得多，也有不少文章用PBE0算过Ag表面，如J. Phys. Chem. C, 117, 5075 (2013)。对于PBE0这种描述色散作用垃圾的泛函，在计算物理吸附时DFT-D校正显然是需要的，见《谈谈“计算时是否需要加DFT-D3色散校正？”》（//www.umsyar.com/413）、《乱谈DFT-D》（//www.umsyar.com/83），要不然结果就搞笑了。虽然目前ORCA也支持DFT-D4了，见《DFT-D4色散校正的简介与使用》（//www.umsyar.com/464），但这里还是用D3，一方面是D3目前已经被极为广泛地使用，比D4更放心一些，另一方面是D4比D3强的地方关键在于能考虑原子带电状态对色散校正的影响，而当前体系里金属基本不带电荷，因此用D4也不会带来什么明显好处（改进较大的是离子、过渡金属配合物类型的体系）。

至于优化此体系用的基组的选择，考虑到当前体系不仅原子数多，而且绝大多数都是重原子，本来耗时就肯定比较高，再加上几何优化对基组不敏感（见《浅谈为什么优化和振动分析不需要用大基组》//www.umsyar.com/387），所以基组用个便宜的2-zeta档次的就够了。lanl2DZ虽然便宜且也能基本满足几何优化目的，但没有标配的辅助基组，如果用autoaux关键词自动产生辅助基组的话不划算。def2-SVP是ORCA用户非常常用的中等偏小的基组，有标配的辅助基组，且对Ag是赝势基组，耗时不会很高，用于优化比较合适。由于此基组对Ag有f极化函数，它会造成耗时增加不少，而对优化来说影响不大，因此我们实际用的是def2-SV(P)，它把f极化函数给砍掉了，同时也砍掉了对氢的p极化，这对优化影响也甚微。

复合物优化任务的输入文件是本文文件包里的complex_opt.inp，在双路XEON E5-2696v3（共36核）服务器上花了10小时算完。产生的最终结构文件是本文文件包里的complex_opt.xyz。为了观看优化过程的收敛趋势，笔者用《OfakeG：使GaussView能够可视化ORCA输出文件的工具》（//www.umsyar.com/498）里的OfakeG工具将输出文件转成了赝Gaussian输出文件，即本文文件包里的complex_opt_fake.out。用GaussView打开，收敛曲线如下所示，可见收敛还蛮顺利的

大家通过gview看结构变化也可以发现，由于复合物初始结构本来摆得就比较理想，所以优化过程中只是苯分子位置稍微调节了一下而已，整体没有明显变化，和苯接触的银原子也没怎么动。

2.4 单点计算

现在产生单点任务的输入文件。启动Multiwfn，载入complex_opt.xyz，然后输入
oi  //产生ORCA输入文件
complex_SP.inp  //产生的文件的名字
4  //选择计算级别

如果你对结合能计算精度要求不是特别高，或者计算条件不是很好，将产生的complex_SP.inp里的B3LYP改成PBE0就可以开始算了，此时对应于PBE0-D3(BJ)/def2-TZVP级别。

不过，有常识的人都知道，精确计算弱相互作用能需要考虑BSSE问题，见《谈谈BSSE校正与Gaussian对它的处理》（//www.umsyar.com/46）。在def2-TZVP下虽然BSSE问题已经不太大了，但对于弱相互作用能计算的影响还是不可忽视的，有四个解决办法解决BSSE问题：
(1)做counterpoise校正。但这需要多花两倍左右时间。若想要用这种做法，参看《在ORCA中做counterpoise校正并计算分子间结合能的例子》（//www.umsyar.com/542）。
(2)用gCP方法以免费方式近似解决BSSE问题。gCP校正的简介见《大体系弱相互作用计算的解决之道》（//www.umsyar.com/214）。gCP对def2-TZVP有现成的参数，在ORCA里要用这种校正就写上gCP(DFT/TZ)即可。而且gCP校正有解析梯度，完全可以用于前面的几何优化过程中，在其输入文件里写上gCP(DFT/SV(P))即可。然而，为什么我们对此例不用gCP？这是因为gCP原本只对前四周期元素拟合了参数，对于之后的元素，会用同族的第四周期元素来凑合，比如Ag的gCP参数会用Cu的参数来凑合。这种凑合的做法，对于第四周期之后元素的原子数较少（比如过渡金属配合物）的时候是可以接受的，但当前体系有一大堆Ag，其参数是影响结合能计算的重点，显然这么凑合是不靠谱的近似。
(3)给基组加弥散函数，即改成ma-def2-TZVP，这样可有效改进弱相互作用计算精度并减小BSSE问题。但这种做法我不推荐用于当前计算，因为此时没有标配的辅助基组，用autoaux自动产生的不是特别划算。而且本来当前体系SCF就比较难收敛，而众所周知加弥散函数会导致SCF收敛难度暴增，对于当前这种原子堆得紧密的情况更是会导致线性依赖问题极度显著，所以加弥散是馊主意。
(4)用更大的def2-QZVP基组。这个基组对于普通泛函的DFT计算而言已经基本达到了完备基组极限了。此时BSSE问题已可以忽略，而且比def2-TZVP时对电子结构描述又有少量改进，而且DFT-D3参数原本又是在def2-QZVP下搞的，因此在这个层面上显得很完美。def2-QZVP的确相当贵，但在ORCA中利用RIJCOSX技术，用这么大的基组对当前体系进行单点计算，耗时依然是完全可接受的。因此，本例使用PBE0-D3(BJ)/def2-QZVP来做单点计算。

还有需要注意的是，通常基组越大SCF越不容易收敛。笔者尝试过，在优化后的复合物结构下用def2-TZVP可以直接收敛（虽然SCF花了多达46圈），而用def2-QZVP则迭代几十次之后依然没有收敛的迹象，即便读def2-TZVP收敛的波函数当初猜也不行。经尝试，我发现把incremental Fock加速技术关闭后可以直接解决，即在输入文件里加上%scf DirectResetFreq=1 end。

（补充：如果你是用ORCA 5.0及以后版本，请忽略下面grid4 gridx4部分的说明，输入文件里也不要加上grid4 gridx4，否则会报错）
另外，建议在关键词里加上grid4 gridx4关键词，使用比默认更好的交换相关泛函积分格点和COSX格点。之所以这样，是因为当前我们希望算的结果比较准，故提升积分格点精度有好处。而且，当前体系大量原子致密堆积，默认积分格点下的积分精度差点意思。比如在def2-TZVP下用默认格点，在SCF收敛后自动做Final grid步骤（提升积分格点精度算最终能量）的时候提示如下
Change in XC energy                          ...    -0.003667837
Integrated number of electrons               ...  1144.003054396
Previous integrated no of electrons          ...  1143.933959260
可见迭代过程用的积分格点对电子数的积分值1143.9339和整数相差不小。即便在final grid时，积分值1144.0030偏离整数仍不可忽略。反之，如果用grid4 gridx4，则输出变为
Change in XC energy                          ...     0.001404285
Integrated number of electrons               ...  1144.000217645
Previous integrated no of electrons          ...  1144.003114364
可见积分精度改进不小，明显更接近于整数。提升积分格点精度还有一个好处是对于某些体系可能令SCF收敛更容易。当前这种特征的体系结合大基组本来SCF就很难收敛，显然多花点耗时来提升SCF收敛的可能性是划算的。总之，对本文这种类型的体系在算单点的时候，都建议用grid4 gridx4，一举多得，而对耗时的增加完全可以接受，而且还可能因为达到收敛限所需步数变少而反倒令耗时更低（在def2-TZVP下确实如此，默认格点下用了46轮收敛，grid4 gridx4下用了37轮，导致后者反倒耗时比前者还低了几分钟）。

考虑到上面提到的问题，最终我们对复合物用的单点任务的输入文件的关键词部分如下
! PBE0 D3 grid4 gridx4 def2-QZVP def2/J RIJCOSX noautostart miniprint nopop
%scf DirectResetFreq=1 end
计算用了40轮收敛，在2*E5-2696 v3 36核机子上花了5个小时。注意这类体系本身就是比较难收敛的，所以切勿看到跑了比如30轮左右还没收敛，而且收敛情况还有些波动就把任务给断了，一定要沉得住气。

复合物算完后，我们把complex_SP.inp复制为benene_SP.inp并把其中银原子都去掉，就成了苯分子的单点任务文件，计算之。再把complex_SP.inp复制为Ag58_SP.inp并把其中苯分子部分去掉，就成了银团簇的单点任务文件，计算之。算完的同名out文件在本文文件包里都提供了。

2.5 计算结合能

从单点任务输出文件中提取单点能，按照E_bind = E(Ag58_Benzene) - E(benzene) - E(Ag58)即可计算出结合能，即627.51*(-8760.730415916696+232.068063781533+8528.63
2371853557)=-18.8 kcal/mol。

前面提到，苯吸附过程的实验焓变是-13 kcal/mol，我们算的和实验定性吻合（比DFT-D3原文里势能面扫描图体现的结果还更好）。由于这个结合能已经算是很大了（苯二聚体结合能才-2.8 kcal/mol），而且这又不是容易算的体系，能算到这种精度就已经不容易了。我们的计算结果和实验值的差异来自于几个方面：
(1)边界效应仍不可忽略，用更大的簇可能能改进结果
(2)PBE0-D3(BJ)理论方法的误差。可尝试其它泛函或结合其它的色散校正方式
(3)实验数据对应焓变，本身和我们算的结合能就缺乏可比性。尤其是，形成复合物后，被吸附物被基底束缚，从而会诞生新的振动模式，这会对复合物的焓产生明显贡献（使之变得更正），忽略了这点会导致算出来的结合能偏大（即overbinding、结合能过负）
(4)忽略了从无限远离到被吸附状态过程中的分子和金属表面结构变化对能量的影响。但由于当前的体系比较刚，这基本可以忽略
(5)实验数据本身就有误差

有一点非常值得一提，也是绝大多数人都忽略的，也就是DFT-D3的三体校正项，这在《DFT-D色散校正的使用》（//www.umsyar.com/210）中专门说过。平时我们用Gaussian、ORCA做计算的时候加D3校正都没考虑三体校正项，因为这项对结合能的影响通常很小。然而，对于原子数很多，尤其是当前体系这样原子还是紧密堆积的情况，三体校正项对结合能的影响可能是不容忽视的。

如上面的博文所述，ORCA里写上ABC关键词就可以在计算D3校正时也考虑三体校正项，但我们前面计算的时候没写这个。为了快速得到考虑三体校正项之后的结果，我们用Grimme独立的dftd3程序来算。把分子坐标弄成xyz格式，运行dftd3 mol.xyz -func pbe0 -bj -abc，输出的就是mol.xyz里的体系在PBE0泛函下的带三体校正项的DFT-D3(BJ)色散修正能了。下面的表格里把ORCA输出文件里的PBE0能量、dftd3给出的带和不带三体校正项(ABC)的DFT-D3(BJ)色散校正能都汇总到了一起

由上表可见，如果不考虑色散校正，结合能是正值，这是因为PBE0几乎完全不能描述色散作用（这是苯与Ag表面吸引作用的最主要来源），考虑DFT-D3(BJ)校正后才定性正确。如果再把三体校正项考虑进去，算出来的结合能-16.9 kcal/mol比不考虑三体校正项时明显更接近实验的吸附过程的焓变，改进了约2 kcal/mol。如果再考虑到把焓的热校正考虑进去还能进一步改进结果，且实验数据本来就有误差，可以认为当前计算模型的精度已经相当不错了。

我建议大家用ORCA基于簇模型计算表面问题时都带着ABC关键词。几何优化、振动分析带不带无所谓（PS：三体校正项没有解析Hessian），关键是单点任务时应当带上，反正对耗时没有增加，还有可能对结果有不可忽视的改进。

2.6 弱相互作用分析

顺带一提，用簇模型算完之后可以结合Multiwfn做各种各样的波函数分析。比如《使用Multiwfn图形化研究弱相互作用》（//www.umsyar.com/68）和《使用Multiwfn做NCI分析展现分子内和分子间弱相互作用》（https://www.bilibili.com/video/av71561024）中介绍的NCI分析是非常流行的图形化考察弱相互作用的方法，这里我们也将之用于苯与银表面的物理吸附问题上。

是使用单点计算产生的波函数还是几何优化产生的波函数来做NCI分析？我建议用后者，因为后者是在def2-SV(P)下产生的，前者是在def2-QZVP下产生的，用后者耗时显然低得多得多，而且NCI分析对于波函数质量并不敏感，基组基本像样即可。故基于优化任务产生的gbw文件，我们就可以做分析了。如《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）所述，要么用orca_2mkl将gbw转成molden输入文件再载入Multiwfn，要么通过设置Multiwfn的settings.ini里的orca_2mklpath使得multiwfn能直接载入gbw文件了。对当前体系做NCI分析值得注意的一点是，为了节约时间，在设置格点时，应当让空间范围只涵盖苯与金属作用的区域，而不要把整个体系都纳入。具体来说，在NCI分析功能设置格点的步骤中应选择选项10 Set box of grid data visually using a GUI window，然后把盒子范围调整成下面这样

虽然当前体系很大，但由于盒子小，在基本能保证图像精细的0.25 Bohr格点间距下，哪怕用4核机子也能很快就算出来结果。如果想要让图像更光滑，可在上图的界面的右下角把格点间距再改小点，比如0.15 Bohr。下图是0.19 Bohr时用VMD绘制出来的结果，色彩刻度用的-0.04~0.03 a.u.

可见苯和银表面之间有大面积的绿色区域，直观展现出了苯的pi电子与银表面电子海洋的大范围色散相互作用，这有点类似于pi-pi堆积的景象（当然，银没有pi电子）。另外，在银与银之间也有等值面出现，我们不用管它，这体现的是Ag-Ag之间的作用，如果想屏蔽掉也可以，做法参看《用Multiwfn+VMD做RDG分析时的一些要点和常见问题》（//www.umsyar.com/291）。

3 总结&其它

本文通过苯在Ag(111)表面吸附的一个简单例子，示例了怎么用簇模型研究表面问题，并且针对相关知识、涉及到的一些细节和要点做了讨论，希望读者举一反三。本文的情况其实还算比较简单的，不牵扯对边界进行特殊处理的问题。

像当前这种不易SCF收敛的体系，可能光是在优化阶段一开始就已经出现SCF不收敛了（这和复合物初始结构也有关。比如把苯放到距离银表面2.7埃的距离，SCF就比3.0埃的时候难收敛，因为其更偏离平衡距离）。这种情况可以适度降低SCF收敛限。ORCA对优化默认用tightSCF收敛限比较严，如果发现SCF差一点就能收敛，但死活就是达不到收敛限，不妨用scfconv7关键词，即把SCF收敛阈值设成单点任务默认阈值与tightSCF之间，这样优化出来的结果一般也足够精确。如果scfconv7都收敛不了，但是用scfconv6能收敛，在优化完了之后建议再用默认的tightSCF进一步优化确保结构准确。如果是初始结构太离谱导致SCF难收敛，自己又不知道大概怎么摆合适，用xtb程序做个预优化也可以，相关信息见《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（//www.umsyar.com/421）。

对于单元素团簇类型体系，SCF收敛难度通常是：稀土>d族>ds族>主族。本文的Ag团簇的情况虽然已经不是很容易SCF收敛了，但还不是真正难搞的。诸如Fe团簇等，收敛难度远高于Ag团簇，而且还有较高可能性出现波函数不稳定等情况，如果是初学者没经验的话不建议轻易去用簇模型做这些金属的表面问题的计算。

值得一提的是，当前体系是C3v对称性的。如果你用Gaussian来算的话，且在计算之前用gview对称化为C3v，计算会快非常多，收敛也会比较容易。

笔者讲授的北京科音中级 培训班（http://www.keinsci.com/workshop/KBQC_content.html）里的过渡态搜索部分有个实例也是基于簇模型的，如下，大家作为练习可以自行做一下。除了要冻结边缘原子外，搜索过渡态和走IRC的过程和普通分子体系无异。

簇模型不仅限于研究表面问题。比如计算酶催化反应，整个蛋白酶动辄几千原子，直接用 算算不动，此时可以把反应位点区域挖出来，即保留底物以及周围一圈跟它发生反应或者有显著弱相互作用的残基的侧链和alpha碳，并且把alpha碳冻结并加氢饱和。这种研究方法可以参看这篇综述：J. Am. Chem. Soc., 139, 6780 (2017)。笔者在北京科音高级 培训班里会讲这种计算的具体细节。虽然也有很多人拿ONIOM(QM:MM)算这种问题，但补力场参数相当麻烦，好多任务和分析都做不了、容易出现乱七八糟问题，非常折腾，从结果上看一般也不比用簇模型有多大好处，所以没有特殊情况的话我很不建议用ONIOM(QM:MM)。笔者专门有一篇文章说这事：《要善用簇模型，不要盲目用ONIOM算蛋白质-小分子相互作用问题》（//www.umsyar.com/597）。

将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法

The method of using the converged wavefunction of Gaussian or other program as the initial guess wavefunction of ORCA

文/Sobereva@北京科音

First release: 2019-Oct-11  Last update: 2021-Oct-1

1 说明

ORCA是非常强大的 程序，笔者之前也写过不少相关文章，开发了不少辅助工具，见www.umsyar.com右侧的ORCA分类。ORCA相对于最常用的 程序Gaussian来说，有一个缺点是SCF收敛做得不够好，很多Gaussian能收敛的情况ORCA难以收敛，而且Gaussian的SCF不收敛的解决方案比较成熟，见《解决SCF不收敛问题的方法》（//www.umsyar.com/61）。另外，ORCA在产生初猜波函数方面也没Gaussian灵活，比如Gaussian用户可以用GaussView构建片段组合波函数作为初猜，而且利用stable=opt关键词还可以自动优化出稳定的波函数，这对于双自由基、反铁磁性耦合体系的研究十分重要，参看《谈谈片段组合波函数与自旋极化单重态》（//www.umsyar.com/82）。虽说ORCA也不是没法算这类自旋极化单重态体系，利用%SCF FlipSpin也可以实现，但往往明显不如Gaussian方便。

显然，如果能让其它 程序，特别是波函数初猜以及SCF迭代方面做得很好的Gaussian产生的波函数作为ORCA的初猜波函数，就能令ORCA取长补短，从很多恼人的问题中解脱。

Multiwfn（主页&下载地址：//www.umsyar.com/multiwfn）程序支持载入.fch、.molden、GAMESS-US/Firefly输出文件这些含有基组定义以及轨道信息的格式（在Multiwfn里称为“基函数信息”），并且可以导出各种含有波函数信息的格式，支持的文件格式详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。Multiwfn支持绝大多数主流量化程序产生的波函数，并可以作为波函数文件格式的转换器来用。从2019-Oct-10更新的Multiwfn开始，Multiwfn的导出文件功能更是新增了产生.mkl文件的功能。.mkl文件是老版本Molekel的输入文件，可以被ORCA自带的orca_2mkl工具转换成.gbw文件，ORCA计算时可以从.gbw文件中读取轨道作为初猜波函数。显然，利用Multiwfn可以直接实现将其它量化程序产生的波函数作为ORCA的初猜波函数的目的。甚至可以说，只要存在一种量化程序有办法得到当前体系当前级别下收敛的波函数，用ORCA计算也一定能收敛。

有一点要提及的是ORCA是基于球谐型高斯函数做的计算，因此用Multiwfn实现这个目的，其它程序计算时也应当用球谐型高斯函数。不了解这方面的话看《谈谈5d、6d型d壳层基函数与它们在Gaussian中的标识》（//www.umsyar.com/51）。

2 例子：丁烷双自由基

这里举个例子，计算丁烷双自由基C4H8，用以说明如何将Gaussian收敛的波函数作为ORCA的初猜波函数，请大家举一反三。本例涉及的文件都在//www.umsyar.com/attach/517/file.rar里。

本例使用2020-Jul-1更新的Multiwfn，Gaussian使用G16 A.03版，ORCA使用4.2版。

首先用Gaussian运行C4H8.gjf，内容如下。任务是对C4H8在B3LYP/def2-SVP级别下产生最稳定波函数，对此体系也是对称破缺波函数。使用def2-SVP的时候Gaussian默认用的是球谐型基函数。

%chk=C:\C4H8.chk
# B3LYP/def2SVP stable=opt
[空行]
ub3lyp/6-31g(d) opted
[空行]
0 1
 C                 -0.74400100    1.78566400    0.00000000
 H                 -0.60282700    2.33865300    0.92499500
 H                 -0.60282700    2.33865300   -0.92499500
 C                 -0.74400100    0.30988100    0.00000000
 H                 -1.25452600   -0.08746700    0.88463900
 H                 -1.25452600   -0.08746700   -0.88463900
 C                  0.74400100   -0.30988100    0.00000000
 H                  1.25452600    0.08746700   -0.88463900
 H                  1.25452600    0.08746700    0.88463900
 C                  0.74400100   -1.78566400    0.00000000
 H                  0.60282700   -2.33865300   -0.92499500
 H                  0.60282700   -2.33865300    0.92499500

然后用formchk将chk转换为fch，载入Multiwfn后依次输入
100  //其它功能 Part 1
2   //导出文件
9   //导出mkl文件
C4H8.mkl
y   //表明产生的mkl文件是给ORCA当初猜用，程序会做特殊处理
现在当前目录下就有了C4H8.mkl。

在当前目录下运行orca_2mkl C4H8 -gbw，就把C4H8.mkl转换为了C4H8.gbw。下面，我们将C4H8.gbw里的轨道作为初猜波函数，用ORCA也在B3LYP/def2-SVP下跑一下这个双自由基的单点。输入文件如下，将之命名为C4H8.inp，把它和C4H8.gbw都放在当前目录下，然后进行计算。运行时ORCA会自动从C4H8.gbw中读取轨道作为初猜波函数。由于C4H8.gbw里的波函数是非限制性波函数，因此写了UKS关键词。由于ORCA的B3LYP和Gaussian的定义不同，所以加了/G来和Gaussian一致。此文件里的结构和上面Gaussian输入文件里的结构精确一致。

! UKS B3LYP/G def2-SVP miniprint nopop
* xyz   0   1
 C                 -0.74400100    1.78566400    0.00000000
 H                 -0.60282700    2.33865300    0.92499500
 H                 -0.60282700    2.33865300   -0.92499500
 C                 -0.74400100    0.30988100    0.00000000
 H                 -1.25452600   -0.08746700    0.88463900
 H                 -1.25452600   -0.08746700   -0.88463900
 C                  0.74400100   -0.30988100    0.00000000
 H                  1.25452600    0.08746700   -0.88463900
 H                  1.25452600    0.08746700    0.88463900
 C                  0.74400100   -1.78566400    0.00000000
 H                  0.60282700   -2.33865300   -0.92499500
 H                  0.60282700   -2.33865300    0.92499500
 *

迭代过程信息如下
ITER       Energy         Delta-E        Max-DP      RMS-DP      [F,P]     Damp
               ***  Starting incremental Fock matrix formation  ***
  0   -157.0020361024   0.000000000000 0.00099750  0.00002892  0.0004319 0.7000
  1   -157.0020389471  -0.000002844760 0.00101224  0.00002810  0.0003439 0.7000
                               ***Turning on DIIS***
  2   -157.0020411727  -0.000002225540 0.00299633  0.00007975  0.0002632 0.0000
  3   -157.0020467355  -0.000005562801 0.00029108  0.00000721  0.0000572 0.0000
  4   -157.0020467759  -0.000000040442 0.00006147  0.00000180  0.0000542 0.0000
                 **** Energy Check signals convergence ****

               *****************************************************
               *                     SUCCESS                       *
               *           SCF CONVERGED AFTER   5 CYCLES          *
               *****************************************************

可见由于用了Gaussian收敛的波函数当做初猜，收敛非常快和顺利，第一轮迭代时的能量和密度矩阵变化就已经极小了，之后很快就收敛了。之所以不是一轮就收敛，是因为Gaussian和ORCA用的DFT积分格点有异。如果二者都用的是HF方法的话，ORCA里仅仅一轮就能收敛。

如果gbw和输入文件不同名，为了能从gbw中读取初猜波函数，应当写上moread关键词，然后加上一行诸如%moinp "/sob/Azusa_Nakano.gbw"来指明读取的gbw文件位置。

为了让Gaussian收敛的波函数放到ORCA里能够收敛的概率尽可能高，有以下一些注意事项和建议
(1)确保Gaussian里用的基组和ORCA里精确一致。比如Gaussian里用6-31G系列基组时，默认是用笛卡尔型d基函数，而ORCA总是用球谐型基函数，因此Gaussian计算时要写5d关键词（不过ORCA用户极少会用Pople系列基组，所以这无所谓）
(2)让Gaussian计算时实际坐标与ORCA一致。Gaussian在计算时会自动将输入朝向摆到标准朝向，因此chk文件最后转出来的gbw里的笛卡尔坐标（标准朝向的）可能和ORCA计算时的不一样，因此要么ORCA输入文件里的坐标就用标准朝向的，要么Gaussian计算时候加上nosymm关键词避免摆到标准朝向下，详见《谈谈Gaussian中的对称性与nosymm关键词的使用》（//www.umsyar.com/297）。
(3)Gaussian为了节约电子积分计算的耗时，用Dunning相关一致性基组等（部分）广义收缩的基组在计算的时候会自动把某些指数重复的primitive GTF给去掉。程序还会把收缩系数很小的primitive GTF给去掉。由于ORCA不自动做这种处理，可能导致个别情况下Gaussian里收敛的波函数拿到ORCA里还是没法收敛。为避免这个问题，可以在Gaussian计算时加上int=NoBasisTransform关键词避免去掉任何primitive GTF，即基组原始怎么定义的就怎么用。如果你没有经验的话，我强烈建议总是带上这个关键词！
(4)建议带上IOp(3/32=2)关键词（尤其是带弥散函数时）避免Gaussian自动去除线性依赖的基函数，以确保和ORCA计算使用的基函数的一致性。
(5)对于DFT计算，如果上面问题都已经考虑了，但放到ORCA里还是不能收敛，可以让Gaussian和ORCA在计算时都用更好的DFT积分格点精度，并且在ORCA里不开RI。

顺带一提，利用Multiwfn还能将其它量化程序产生的波函数作为GAMESS-US的初猜波函数，因为如果载入的波函数文件含有基函数信息，Multiwfn产生的GAMESS-US输入文件里可以直接带初猜信息。Multiwfn也可以将其它量化程序产生的波函数当Gaussian的初猜，因为Multiwfn可以导出fch格式，用Gaussian的unfchk转成chk后，就可以用guess=read来从中读取波函数当初猜了。若有不清楚的，参看前述的《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》。

3 chk到gbw的批量转换脚本

为了让大家可以更方便地把chk转换成gbw，笔者写了一个批量自动转换的bash小脚本chk2gbw.sh，可以在这里下载：//www.umsyar.com/attach/517/chk2gbw.sh。运行前需要先把Gaussian、Multiwfn、ORCA都安装好，运行过程中会自动调用formchk、Multiwfn和orca_2mkl工具。

运行此脚本后，脚本会对当前目录下每个chk文件进行转换。比如有个文件叫popipa.chk，运行此脚本后会得到popipa_Gau.gbw，还会得到popipa.inp，这是ORCA输入文件，打开一看就会看到里面已经写了moread关键词和%moinp "popipa_Gau.gbw"，坐标和chk里的直接对应，因此用户现在只需要把定义计算级别的关键词改成实际情况即可开始ORCA计算。可见此脚本方便至极！

运行时输出信息示例：
Converting B2H6.chk to B2H6.gbw ... (1 of 3)
Converting C2H5F.chk to C2H5F.gbw ... (2 of 3)
Converting CCl4.chk to CCl4.gbw ... (3 of 3)

4 gbw到chk的批量转换脚本

从2020年8月21日更新的Multiwfn开始，在其examples\scripts目录下提供了gbw2chk.sh文件，只要运行，就会把当前目录下的所有gbw文件转化为同名的chk文件，以便通过Gaussian使用guess=read关键词读取其中的波函数当初猜。运行前需要先把Gaussian、Multiwfn、ORCA都安装好。注：如果波函数是UHF/UKS计算得到的，必须用2021年10月1日及以后更新的Multiwfn。

OfakeG: A tool that enables GaussView to visualize ORCA output files

文/Sobereva@北京科音

First release: 2019-Jul-17  Last update: 2024-Oct-30

1 前言

程序ORCA用的人越来越多，功能很强大而且免费，用户数在所有量化程序中已经是第二高（虽然跟Gaussian比还遥不可及）。但至少在笔者撰写此文时，相对于用户数占绝对主导地位的Gaussian程序而言，仍有一个不足之处是没有像GaussView那样的很理想的图形界面。虽然也有Avogadro、Chemcraft、Gabedit等程序能支持ORCA，但都没GaussView用着舒服。

在产生输入文件方面，Multiwfn已经提供了产生ORCA常见任务的输入文件的功能，见《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490），用户只需要用GaussView画好结构，保存为gjf/mol/mol2/pdb格式，就可以用Multiwfn很方便地得到ORCA输入文件，所以在建模、产生输入文件方面，对ORCA用户没什么困难的。在观看ORCA产生的轨道、做波函数分析方面，Multiwfn也都提供了极其丰富的功能，相关信息见《使用Multiwfn观看分子轨道》（//www.umsyar.com/269）、《Multiwfn FAQ》（//www.umsyar.com/452）等文章，因此ORCA用户在后处理分析方面也没任何压力。另外，Multiwfn可以基于ORCA的输出文件绘制各类光谱图，所以ORCA用户在光谱研究方面也已经很方便了，相关信息见《Simulating UV-Vis and ECD spectra using ORCA and Multiwfn》（//www.umsyar.com/485）、《Simulating UV-Vis and ECD spectra using ORCA and Multiwfn》（//www.umsyar.com/485）。笔者在北京科音高级 培训班（http://www.keinsci.com/workshop/KAQC_content.html）里还对ORCA的各方面使用做了相当全面系统的讲解，可以很快上手并用得游刃有余。

虽然Multiwfn已经解决了ORCA用户在使用方面的大量障碍，但几何优化轨迹/收敛情况的考察以及振动模式的观看不属于Multiwfn的范畴，而目前却没有理想的解决办法。虽然Chemcraft也能观看ORCA的优化轨迹，但终究没有常用的GaussView用着舒服，而且还是收费程序；Avogadro虽然也观看ORCA的振动分析对应的振动动画，但显示效果不理想，在Windows下容易崩溃、出现异常（至少在笔者的Win7-64bit机子上是如此）。尽管也可以让Gaussian与ORCA挂接，这样可以使ORCA做计算但是输出的是Gaussian格式的信息，从而等效实现让GaussView观看ORCA计算的结果的目的，见《将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析》（//www.umsyar.com/422），但这样做稍显麻烦、在Windows下也没法用。

显然，如果能开发个程序把ORCA的优化、振动分析、优化+振动分析的输出文件“伪造”成Gaussian的，这样就可以令GaussView直接支持读取ORCA的输出文件了，使得ORCA对于常见问题的研究用起来方便得多，也明显便于Gaussian用户同时掌握ORCA程序。笔者开发的OfakeG程序就是实现这个目的，下面介绍一下。如果你想一睹为快这个程序的实际效果，可以看这个视频：《基于ORCA 程序对分子做优化、振动分析、观看红外光谱、观看轨道的简单演示》（https://www.bilibili.com/video/av59599938），其中用到了此程序。

OfakeG的学术合理性声明：本程序的开发灵感来自于Grimme的xtb程序。xtb程序做振动分析的时候会自动输出一个伪造Gaussian的振动分析输出文件，目的是为了让用户看振动模式方便；xtb程序为了兼容GSM也官方支持伪造ORCA的输出文件。大牛Grimme直接用的就是fake这个词。显然令A程序输出B程序的格式在学术界是非常正常的事情。Gaussian的输出格式是公开的而非加密的，GaussView能读入的格式也相当于是公开的，本文的OfakeG程序亦没有对GaussView本身做任何篡改，明显从各个角度上本程序的开发是完全学术正当的。本程序所做的事仅仅是将ORCA的输出信息转化成Gaussian的格式而已，文中所谓“伪造”只不过是常规的文件格式转换而已，和数据层面的“造假”有天壤之别，转换出的文件里也根本没有任何文字体现这文件是靠Gaussian程序算出来的。此程序的开发目的是给广大科研工作者提供个便利，开发/使用此程序不涉及任何侵权和学术不端（除非你利用OfakeG之后，把ORCA算的结果说成是Gaussian算的）。此程序愿意用就用，不爱用的、缺乏对计算化学领域程序状况基本认识的、怀有恶意的、有特殊利益驱动的人以及杠精，不要强词夺理在学术合理性上乱喷此程序。

2 OfakeG程序

2.1 介绍+使用方法

OfakeG程序可以在官方页面下载：//www.umsyar.com/soft/OfakeG
其中带.exe后缀的是Windows版，不带后缀的是Linux版。

此程序目前支持处理ORCA的opt、freq和opt freq任务的输出文件（不支持单点任务文件，因为根本没有任何转换的意义！）。此程序对ORCA 4.x、5.x、6.x版经测试可用，对于其它版本可能兼容也可能不兼容，请读者自行尝试。等ORCA以后出新版本，并且笔者发现和OfakeG不兼容时，预计笔者会更新此程序并更新本文。

OfakeG使用非常简单。启动此程序后，把上述任务的ORCA输出文件路径输入进去（对于Windows也可以直接把文件拖进去，路径会直接显示出来），一按回车，就会在当前目录下产生伪造的Gaussian输出文件。如果输入文件名字是yuri.out，则输出文件将是yuri_fake.out。这个输出文件可以直接载入GaussView，对opt或opt freq任务可以播放优化过程的动画、用results - Optimization观看优化过程的收敛情况，对freq或opt freq任务可以用results - Vibrations观看振动模式。

在Windows下还有更省事的运行方式，即可以直接将ORCA输出文件拖到OfakeG.exe图标上，此时在ORCA输出文件的目录下会出现文件名带_fake的伪造的Gaussian输出文件。

OfakeG也可以通过命令行方式使用，比如在Linux下可以在OfakeG所在目录下运行./OfakeG Aika.out，将在当前目录下得到Aika_fake.out。显然，你也可以自写shell脚本用这个程序大批量转换ORCA输出文件。

OfakeG文件包里的.out文件是一些ORCA的示例输出文件。如果你的输出文件转换不成功，请尝试通过对照这些示例文件搞清楚是怎么回事。目前OfakeG名义上只支持HF/DFT的输出文件，其它理论方法不一定能支持。对于加了乱七八糟复杂关键词的情况，OfakeG也不一定能处理。

如果OfakeG处理你的文件时崩溃，且得到的_fake后缀的文件里只有几行信息，很有可能是因为你的ORCA输出文件的编码是UTF16造成的，OfakeG是处理不了前者的情况的，是什么编码和你用的终端有关系。比如Windows的cmd终端重定向输出的文件是ASCII编码的，而PowerShell是UTF16编码的。对于UTF16编码的输出文件，你可以用比如Ultraedit打开，选另存为，把编码改成Unicode或UTF8，之后再用OFakeG处理。如果你平时习惯用PowerShell且希望重定向出的文件直接就ASCII编码，可以用诸如这样的命令运行test.inp得到test.out：D:\study\orca\orca test.inp | out-file test.out -encoding ascii
如果是Win10，还可以直接指定默认的重定向的编码，详见https://stackoverflow.com/questions/40098771/changing-powershells-default-output-encoding-to-utf-8。

如果你怎么也搞不清楚为什么你的ORCA输出文件无法转化成功，或者可判定OfakeG程序有bug，请在http://bbs.keinsci.com/thread-13952-1-1.html贴子里发回帖，把文件压缩后上传。

2.2 OfakeG的几个细节

以下内容建议留意一下，以更好地理解OfakeG的细节，但初学者不看也可以。

OfakeG给出的是简化到不能再简化的能令GaussView正常读取的伪造的Gaussian输出文件，因此如果你写类似工具把其它程序的输出文件也伪造成类似格式，就也可以令GaussView读取。

GaussView要求输出文件里必须有basis functions、alpha electrons、beta electrons信息，但ORCA输出文件里不直接体现，而且这仨对于观看优化和振动分析没有意义，因此在伪造的输出文件开头有这仨信息，但数值都为0。

OfakeG从ORCA输出文件里读能量的时候读的是FINAL SINGLE POINT ENERGY，即当前计算级别下的最终能量。而产生伪造的Gaussian输出文件时，为了省事和统一，是以SCF Done标签来输出的。

Gaussian做优化任务的时候，对每一步，输出次序是[结构i]-[结构i的能量]-[结构i的受力]-[结构i的收敛情况]，所以结构、能量、收敛信息都是一一对应的。而对于ORCA，输出也是这样的顺序，但最后在第i步时发现已经满足收敛限了，之后还会根据第i步的信息再预测出第i+1步的结构，并且计算这个结构下的能量（也顺带得到波函数），而这个i+1结构就不再计算受力了，也因此对这个结构也不再输出收敛判断信息。所以OfakeG产生的伪造的Gaussian输出文件中，第1步到第i+1步的结构、能量、收敛情况都会给出，但最后一次输出的收敛情况信息里当前值全都被设成0来占位。

优化过程中除了像Gaussian一样用受力/位移的最大/RMS值作为判断标准外，ORCA还用能量变化作为判断标准。为了体现这点，在伪造的Gaussian输出文件中也在收敛判断部分添加了这项，但这项不会被GaussView所读取，大家可以自行考察。

对于振动分析，由于ORCA不会给出约化质量和振动模式的力常数，所以伪造的Gaussian输出文件里也没这项，这不影响一般的分析。由于ORCA不给出振动模式的不可约表示，所以OfakeG把不可约表示都一律输出为A。

OfakeG把ORCA振动分析输出的热力学数据也都转化为了Gaussian的输出形式，对于用惯了Gaussian的人来说读起来方便不少，并且还顺带多显示了一项Electronic energy=，后面是振动分析对应的结构的电子能量。

OfakeG以后版本也有可能支持处理ORCA的IRC任务的输出文件，但目前没有打算支持。因为笔者撰文时最新的ORCA 4.1.2版的IRC功能非常弱、速度慢，甚至就连反应坐标都不给出来，原理上没法转换成Gaussian的格式。另外OfakeG也不会去支持转换ORCA的TDDFT等电子激发任务的输出文件，因为做这个转换没有任何实际意义。Multiwfn直接就能基于ORCA的TDDFT输出信息绘制各种电子光谱和做电子激发分析（后者我都有现成的例子，见//www.umsyar.com/485。Multiwfn绘制光谱的更详细介绍见//www.umsyar.com/224），而且ORCA目前版本给出的是TDDFT组态函数的贡献而不是系数，原理上也不可能转换为Gaussian形式的输出。

OfakeG是100%纯Fortran写的，没有利用任何库和其它任何编程语言（或许有的人能猜到我为什么刻意彰显这点）。

详谈Multiwfn产生ORCA 程序的输入文件的功能

On the function of generating ORCA input file in Multiwfn

文/Sobereva@北京科音

First release: 2019-Jun-23  Last update: 2024-Oct-17

1 简介

ORCA 程序如今的用户越来越多，用户数仅次于Gaussian，不仅免费，而且有很多Gaussian远不能及的优点，支持很多Gaussian不支持的重要方法，并且由于充分利用了RI技术，使得DFT、TDDFT的效率极高，这点在《大体系弱相互作用计算的解决之道》（//www.umsyar.com/214）、《使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱》（//www.umsyar.com/462）我都已经说过。ORCA的安装方法见《 程序ORCA的安装方法》（//www.umsyar.com/451）。

波函数分析程序Multiwfn可以产生ORCA常见任务的输入文件，这个功能实用性极高。产生的关键词对于ORCA最新版都是合适的。随着ORCA的更新、用法的变动，Multiwfn也会相应地更新。最新版Multiwfn可以从主页//www.umsyar.com/multiwfn免费获取。

Multiwfn的这个功能并不是那种类似GaussView产生Gaussian输入文件的界面那样可以充分指定各种各样的选项，我觉得这种创建输入文件的工具对于ORCA用户而言并没什么实际意义。实际上有很大一批人用ORCA只是冲着ORCA的那些长处、特色功能去的，特别常用的任务类型就那么些，他们没精力也不需要把ORCA的所有细节都摸清楚。Multiwfn的这个产生ORCA输入文件的功能更像是把当前载入的体系的坐标套入特定类型任务的关键词模板里，之后用户只需把内存使用量、核数等参数简单地改一下，马上就可以跑，这显著拉低了ORCA的使用门槛，也使得ORCA用起来方便不少。另外，ORCA很多功能的速度之所以远胜于Gaussian关键是靠着RI，然而如何正确地在ORCA里使用RI，给初学者可能三言两语说不明白，而使用Multiwfn产生的ORCA输入文件中与RI相关的设定总是最恰当的，使得初学者也可以轻松享受ORCA带来的巨大优越。

如果你的研究中使用了Multiwfn创建了ORCA输入文件并给你带来了便利，希望在写文章的时候提及诸如The input files of ORCA program were prepared with the help of Multiwfn code并引用Multiwfn程序在刚启动时显示的Multiwfn原文，这是对Multiwfn此功能开发最好的支持。

注意，本文介绍的功能仅仅能让你用ORCA做一些最简单、常见的计算，以及在使用上提供便利。而正经用ORCA做计算的话是绝对要花时间专门具体学习ORCA的用法和相关理论知识的。非常推荐参加北京科音高级 培训班（http://www.keinsci.com/workshop/KAQC_content.html），里面专门有一节用大约300页幻灯片对ORCA的使用进行非常系统的讲解并给出大量使用经验，并且还有一节对于RI/密度拟合和COSX这个关键知识有特别详细的介绍，具备了这些知识才可能把ORCA用得得心应手。

2 产生ORCA输入文件的界面

可以用任意Multiwfn支持的含有结构信息的文件作为Multiwfn的这个功能的输入文件，比如fch、molden、wfn、mol、mol2、pdb、xyz、gjf等。然后进入主功能100的子功能2中的选项12（更快速进入的方法是在主菜单里直接输入oi），输入要产生的ORCA输入文件的文件名，之后就会看到以下界面。注意随着Multiwfn的版本更新，这个界面也可能会有所变化，届时笔者也会相应地更新本文。

-100 Use template input file provided by user to generate new input file
-11 Choose ORCA version compatibility, current: >ORCA 5.0
-10 Set computational resources, core:   8 memory/core:  1000 MB
-2 Toggle adding diffuse functions, current: No
-1 Toggle employing implicit solvation model, current: No
 0 Select task, current: Single point
 1 B97-3c      1b r2SCAN-3c
 2 RI-BLYP-D3(BJ)/def2-TZVP
 3 RI-B3LYP-D3(BJ)/def2-TZVP(-f)     4 RI-B3LYP-D3(BJ)/def2-TZVP
 5 RI-wB97M-V/def2-TZVP
 6 RI-PWPB95-D4/def2-TZVPP       7 RI-PWPB95-D4/def2-QZVPP
 6b RI-wB97X-2-D3(BJ)/def2-TZVPP     7b RI-wB97X-2-D3(BJ)/def2-QZVPP
 6c RI-revDSD-PBEP86-D4/def2-TZVPP   7c RI-revDSD-PBEP86-D4/def2-QZVPP
 8 DLPNO-CCSD(T)/def2-TZVPP with normalPNO and RIJCOSX
 9 DLPNO-CCSD(T)/def2-TZVPP with tightPNO and RIJCOSX
 10 CCSD(T)/cc-pVTZ
 11 CCSD(T)-F12/cc-pVDZ-F12 with RI
 12 Approximated CCSD(T)/CBS with help of MP2 (cc-pVTZ->QZ extrapolation)
 13 DLPNO-CCSD(T)/CBS with RIJCOSX & tightPNO (def2-TZVPP->QZVPP extrapolation)
 14 CCSD(T)/CBS (cc-pVTZ->QZ extrapolation)
 20 sTD-DFT based on RI-wB97X-D3/def2-SV(P) orbitals
 22 TDDFT RI-PBE0/def2-SV(P)
 23 TDDFT RI-RSX-QIDH/def2-TZVP     231 TDDFT RI-DSD-PBEP86/def2-TZVP
 24 EOM-CCSD/cc-pVTZ                 25 STEOM-DLPNO-CCSD/def2-TZVP(-f)

上述选项里，从1开始都是最常用的ORCA计算级别，20之前的都是基态计算任务，从20开始的都是激发态计算任务，都是序号越大精度越高但耗时也越高。选了之后就会生成ORCA输入文件，关键词是相应级别计算时最恰当的。之后用户可以再根据实际需要把输入文件里的基组、核数（%pal nprocs后面的值）、内存使用量（%maxcore后面的值）改成实际情况。

由于ORCA不同大版本的情况不同，某些关键词写法相差甚大，为了兼容性考虑，界面里专门留了个选项-11可以用来切换ORCA版本。默认对应的是最新的ORCA版本。

选择计算级别前，可以先用选项0选择任务类型。默认是单点，也可以选优化极小点、振动分析、优化极小点+振动分析、优化过渡态+振动分析、分子动力学、考虑counterpoise校正算结合能。对于振动分析任务，关键词会自动加上tightSCF来使用更严格的SCF收敛限（而对于优化任务ORCA默认就是用tightSCF）。分子动力学任务默认设置对应的是控温在298.15 K下用0.5 fs步长模拟500步。

默认不使用溶剂，如果想用SMD溶剂模型表现溶剂环境，可以选择选项-1，然后输入溶剂名，之后导出的输入文件里就带着溶剂模型关键词了。ORCA直接支持的SMD溶剂名详见手册（一搜1-HEXANOL就可以跳到对应的页）。若在这个选项里输入c，也可以用CPCM溶剂模型、自定义溶剂。

如《谈谈弥散函数和“月份”基组》（//www.umsyar.com/119）中所述，很多情况用弥散函数是必须的，比如阴离子体系的单点等。如果你想把基组改成带弥散函数的版本，就选择一次选项-2将其状态切换为Yes。此时cc-pVnZ系列基组都会被改为aug-cc-pVnZ，def2系列基组都会被改为ma-def2系列，不熟悉此基组的话见《给ahlrichs的def2系列基组加弥散的方法》（//www.umsyar.com/340）。与此同时，辅助基组也会改成恰当的。对于ma-def2系列，由于没有标配的辅助基组，所以会自动用autoaux关键词来自动构建。仅B97-3c、r2SCAN-3c和CCSD(T)-F12/cc-pVDZ with RI没法加弥散函数。

如果你载入进Multiwfn的文件是比如fch、wfn、molden这样带有波函数信息的，那么原本波函数对应的电荷和自旋多重度是多少，则产生的ORCA输入文件里的值也是多少。如果载入的是诸如xyz、pdb这样只含有结构信息的文件，那么别忘了需要根据实际情况自己改一下电荷和自旋多重度。

注意有些方法和任务类型之间不兼容，比如对于ORCA 5.0版而言，PWPB95-D3(BJ)没有解析梯度而无法直接用于几何优化任务。也有些方法本身可以用于优化，比如TDDFT，但使用了SMD溶剂模型后就不再支持。这种情况请大家自行手动调整输入文件里的关键词，比如需要改用数值梯度、数值Hessian的自己分别写上numgrad、numfreq

输入文件里的%maxcore控制ORCA并行计算时每个进程用的内存量（MB）的上限。因此比如你想8核并行，机子有32GB物理内存，若扣除操作系统、后台任务占的内存，假设剩30GB（约30000MB），那么maxcore不应超过30000/8=3750，但鉴于ORCA计算时进程实际使用的内存量上限往往会超过maxcore，因此为保险起见此时maxcore建议设为3000。计算使用的CPU核数、maxcore设置既可以自己编辑Multiwfn产生的ORCA输入文件，也可以直接在当前界面里用选项-10修改，默认的核数与Multiwfn的settings.ini里的nthreads参数等同。

Multiwfn产生的ORCA输入文件里都带着noautostart miniprint关键词。ORCA会自动试图从当前目录下与当前任务同名的gbw里读取结构和波函数作为初猜，noautostart代表要求不这样干。miniprint代表避免输出一些对普通用户没什么意义的中间信息，且也不自动做布居分析，因为ORCA默认做的布居分析的输出会导致输出文件的信息量巨大，而这些信息没太大意义而且格式还不好读。真想要做布居分析讨论电子结构，随时可以用Multiwfn基于ORCA产生的molden文件来实现，又快又方便信息又容易读，见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）以及Multiwfn手册4.7节的布居分析例子。

此界面里有个选项-100，选了之后用户要提供一个ORCA模板文件，其中包含实际计算要用的所有设置，只不过坐标部分写成[geometry]，Multiwfn会将这部分替换为当前体系的坐标来产生ORCA输入文件。通过这个功能可以完全自定义地创建ORCA输入文件，还可以通过自写简单的脚本实现诸如把一大批xyz文件转化为带有特殊关键词的ORCA输入文件。一个模板文件的例子如下，用来产生DSD-BLYP双杂化泛函级别的自然轨道（算完后产生的以.mp2nat为后缀的文件可自行改成.gbw后缀，然后就可以用orca_2mkl转成.molden输入文件，之后可以载入Multiwfn做波函数分析）。
! DSD-BLYP def2-TZVP tightscf miniprint
%maxcore  6000
%pal nprocs   10 end
%mp2
density relaxed
NatOrbs true
end
* xyz   0   1
[geometry]
*

3 使用Multiwfn产生ORCA输入文件的简单例子

此例我们要用ORCA对甲醇在B97-3c级别下优化，之后用wB97M-V/def2-TZVP算单点，算单点时使用SMD模型表现乙醇环境。

首先用GaussView（或其它建模程序）画个甲醇，保存为gjf或mol或mol2或pdb文件，然后载入启动Multiwfn，载入此文件，依次输入
100  //其它功能，Part 1
2   //导出文件、产生量化程序输入文件
12  //产生ORCA输入文件
opt.inp   //输出的文件名
0  //选择任务类型
2  //优化
1  //B97-3c
现在当前目录下出现了opt.inp，恰当设置里面的maxcore和nprocs，然后用ORCA运行之。

算完后当前目录下出现了opt.xyz，此为优化后的结构文件，将之载入Multiwfn，依次输入
100  //其它功能，Part 1
2   //导出文件、产生量化程序输入文件
12  //产生ORCA输入文件
SP.inp   //文件名
-1  //启用SMD溶剂模型
ethanol  //溶剂名
5  //RI-wB97M-V/def2-TZVP
然后将当前目录下得到的SP.inp再次用ORCA运算即可。

在《Simulating UV-Vis and ECD spectra using ORCA and Multiwfn》（//www.umsyar.com/485）一文中，笔者还演示了使用了Multiwfn创建ORCA输入文件，然后通过ORCA做TDDFT计算，最后用Multiwfn产生电子光谱的全过程。在《基于ORCA 程序对分子做优化、振动分析、观看红外光谱、观看轨道的简单演示》（https://www.bilibili.com/video/av59599938）视频中，笔者演示了利用ORCA结合Multiwfn等程序，对有机小分子做一系列最常见的计算任务的基本操作。

4 计算级别简介和关键词解读

下面解释一下前述计算级别的特点以及解读一下相关的关键词，以便于读者准确地认识这些计算级别的特点、什么时候适合用什么级别、为什么笔者令Multiwfn支持这些级别，以及理解为什么关键词是那样定义的。

注意采用的基组没有一个是Pople系列基组，因为如《谈谈 中基组的选择》（//www.umsyar.com/336）提到的，这套基组对于较好精度的计算没有一个是划算的，而对于便宜的级别，这套基组又没有标配的给RI用的辅助基组，因此Pople系列基组在ORCA中几乎是摆设。DFT相关的计算笔者配的都是def2系列，这非常适合DFT，而对于后HF类型计算用的都是Dunning相关一致性基组，但实际上改用档次相当的def2系列基组也完全可以。

从ORCA 5.0开始，对于杂化泛函默认就会用RIJCOSX，因此下面关键词里面RIJCOSX对于ORCA >=5.0版做杂化泛函计算其实是多余的，但为了与老版本的兼容性考虑还是留着了。

4.1 基态计算

• B97-3c。关键词：B97-3c
这是Grimme提出来的一个又便宜又快的组合式方法，用的是纯泛函，基组是方法直接内定的，还带了DFT-D3、SRB校正项。对于主族和过渡金属体系都适用。在ORCA里的耗时仅略高于RI-BLYP/def2-SVP一点点，但结果肯定整体更好。更详细介绍见《盘点Grimme迄今对理论化学的贡献》（//www.umsyar.com/388）。我建议使用此方法去优化、计算大体系能量，或者粗略计算小体系用于初步筛选的目的，比如笔者在《gentor：扫描方式做分子构象搜索的便捷工具》（http://bbs.keinsci.com/thread-2388-1-1.html）的实例中就是这样做的。根据笔者自己的一些研究，用Intel 36核机子用B97-3c在真空下优化一个160原子的弱相互作用显著的有机体系，用了75步收敛，总共耗时才不到三个小时，比起在Gaussian下用算这种体系常用的B3LYP-D3(BJ)/6-311G*耗时低得多得多。之后在同样条件下用B97-3c对这个体系做振动分析，耗时三个半小时，也完全可以接受。使用B97-3c时不需要写与RI相关的关键词和指定辅助基组，因为默认就会用RIJ和恰当的辅助基组加速计算。

• r2SCAN-3c。关键词：r2SCAN-3c
这是B97-3c的后继者，基于2020年末提出的r2SCAN泛函弄的-3c方法。耗时比之高了百分之几十，但精度全面提升了不少。因此如果不是特别穷的话，建议总是用r2SCAN-3c代替B97-3c。

• RI-BLYP-D3(BJ)/def2-TZVP。关键词：BLYP D3 def2-TZVP def2/J
BLYP泛函平时用得不多，毕竟算有机体系精度远不如杂化泛函，但是BLYP这样的纯泛函在ORCA下利用RIJ方式加速库仑部分的计算后（记为RI-BLYP），对于稍大点的体系，速度能比Gaussian快一个数量级以上（振动分析除外）！。因此以BLYP为典型代表的纯泛函还是很有存在意义的。关键词里D3代表加上DFT-D3(BJ)色散校正，详见《DFT-D色散校正的使用》（//www.umsyar.com/210）。对于算弱相互作用，根据GMTKN55测试集的测试，BLYP加上DFT-D3(BJ)后，在所有的“纯泛函-D3”里是顶尖的。def2-TZVP虽然是一个中等偏上的基组，在Gaussian里算是偏贵的，但是在ORCA里利用RIJ加速后其实挺便宜的。使用RIJ技术加速需要指定辅助基组，关键词中def2/J就代表使用对def2系列基组通用的用于RIJ目的的辅助基组。不需要特意写RIJ关键词，因为对纯泛函默认会使用RIJ。注意BLYP-D3(BJ)/def2-TZVP的耗时比B97-3c稍高，但精度未必比B97-3c好。特别是从ORCA 5.0开始支持了r2SCAN-3c后，RI-BLYP-D3(BJ)/def2-TZVP这个组合其实就没什么任何使用价值了。

• RI-B3LYP-D3(BJ)/def2-TZVP(-f)。关键词：B3LYP D3 def2-TZVP(-f) def2/J RIJCOSX
由于B3LYP便宜、稳健、被支持广泛，虽然算热化学方面的精度已经严重落伍，但即便到现在还是被使用得最多的泛函，结合DFT-D色散校正后又令它的生命周期进一步延长，相关信息看《简谈 计算中DFT泛函的选择》（//www.umsyar.com/272）。对于B3LYP这样的杂化泛函，ORCA里可以用RIJK或RIJCOSX来加速，后者在ORCA里能用的特征更多，而且对于二、三十个原子以上的情况，通常RIJCOSX比RIJK明显更快，RIJCOSX关键词代表启用之。def2-TZVP(-f)是在def2-TZVP基组的基础之上把f极化函数阉割掉的结果，这可以显著节约时间，而对精度损失不算特别严重，此时尺寸与6-311G(2d,p)相仿佛。杂化泛函结合RIJCOSX的耗时显著高于纯泛函结合RIJ，因此RI-B3LYP即便结合def2-TZVP(-f)，耗时也比RI-BLYP/def2-TZVP高得多，且速度远胜于在Gaussian里的计算速度（Gaussian对杂化泛函不支持RI）。对于主族体系各方面计算以及弱相互作用计算而言，B3LYP-D3(BJ)比上述的BLYP-D3(BJ)好是肯定的，因此是否适合使用B3LYP-D3(BJ)应当自行掂量。

• RI-B3LYP-D3(BJ)/def2-TZVP。关键词：B3LYP D3 def2-TZVP def2/J RIJCOSX
同上，只不过基组变成了完整的def2-TZVP

• RI-wB97M-V/def2-TZVP。关键词：wB97M-V def2-TZVP def2/J RIJCOSX strongSCF
根据Phys. Chem. Chem. Phys., 19, 32184 (2017)、Mol. Phys., 115, 2315 (2017)、J. Chem. Theory Comput., 15, 3610 (2019)等测试，ωB97M-V在除了双杂化泛函以外的泛函中，无论算主族还是含过渡金属的体系，无论算反应能、势垒还是弱相互作用，性能都是顶级的，在《简谈 计算中DFT泛函的选择》我也提到了。虽然从ORCA 5.0开始支持了ωB97M-V的解析梯度，但没有解析Hessian，因此不便于做振动分析检验虚频情况，所以我不建议用这个泛函做优化，而且本身这个泛函也比B3LYP等更贵一些。另外，众所周知，几何优化耗时远高于单点计算，而几何优化用的级别可以比单点要低，因为几何优化结果对计算级别敏感度远低于能量计算，这点在《浅谈为什么优化和振动分析不需要用大基组》（//www.umsyar.com/387）专门提了，因此推荐大家用相对便宜的B97-3c或比它更好也更贵的RI-B3LYP-D3(BJ)/def2-TZVP(-f)来优化结构，最后用RI-wB97M-V/def2-TZVP来算单点。对于ωB97M-V这样的非双杂化泛函而言，对于一般问题用def2-TZVP就已经够用了，要求更高的话也可以考虑def2-QZVPP。如果你要算弱相互作用能的话也推荐用def2-QZVPP，若结合counterpoise校正还能更好一点。由于ORCA的默认的SCF收敛限相对于wB97M-V的精度水平来说过于宽松了些，因此同时用了strongSCF来让能量收敛到更高精度。

值得一提的是虽然M06-2X很流行，很适合算主族体系，但在上面只字未提。这是因为M06-2X做几何优化方面比B3LYP-D3(BJ)优势并不明显，而耗时明显更高，对积分格点要求也高，还更难收敛。虽然M06-2X用于算主族体系的能量很不错，但是和ωB97M-V比又明显逊色，在速度上也没优势，故M06-2X在当前的ORCA中没多大用武之地。

• RI-PWPB95-D4/def2-TZVPP。关键词：PWPB95 D4 def2-TZVPP def2/J def2-TZVPP/C RIJCOSX tightSCF
PWPPWPB95-D4在《简谈 计算中DFT泛函的选择》中专门说了，是很稳健且精度很好的双杂化泛函，精度介于ωB97M-V和CCSD(T)档次之间。当前关键词令PWPB95-D4泛函的杂化泛函部分做SCF的过程中通过RIJCOSX加速来降低耗时，而在之后计算类似MP2部分的时候RI来显著降低耗时（只要用了RIJCOSX关键词，默认就会在MP2部分也用RI），这一步使用def2-TZVPP/C辅助基组。为了让SCF部分尽量准确从而得到数值层面较准确的PWPB95-D4的结果，因此用了tightSCF。笔者用Intel 36核服务器做过测试，对于168个原子的有机体系的单点，这样的计算级别耗时不到三个小时，耗硬盘也不多，而如果用Gaussian的话根本没戏，差不多80个原子就封顶了，还特别耗硬盘。对于66个原子的有机体系的单点，ORCA下用这个级别仅需六分钟就能算完。总的来说，RI-双杂化/def2-TZVPP用Intel三、四十核的服务器跑200原子以内的单点毫无压力。另外，如果你希望耗时更低的话，可以手动加上float关键词，代表通过单精度变量而非默认的双精度变量储存中间数据，耗时能节约一小半，而且还能省一倍硬盘。aussian的话根本没戏，差不多80个原子就封顶了，还特别耗硬盘。对于66个原子的有机体系的单点，ORCA下用这个级别仅需六分钟就能算完。总的来说，RI-双杂化/def2-TZVPP用Intel三、四十核的服务器跑200原子以内的单点毫无压力。另外，如果你希望耗时更低的话，可以手动加上float关键词，代表通过单精度变量而非默认的双精度变量储存中间数据，耗时能节约一小半，而且还能省一倍硬盘。

• RI-PWPB95-D4/def2-QZVPP。关键词：PWPB95 D4 def2-QZVPP def2/J def2-QZVPP/C RIJCOSX tightSCF
我之前在《谈谈 中基组的选择》（//www.umsyar.com/336）中说过，后HF、双杂化泛函对基组的要求显著高于普通泛函，因此为了充分发挥PWPB95-D4的潜力，如果你能接受更大计算量的话，推荐结合def2-QZVPP。对于66个原子的有机体系的单点，在Intel 36核服务器下ORCA下用这个级别用了30分钟算完。

• RI-wB97X-2-D3(BJ)/def2-TZVPP（RI-wB97X-2-D3(BJ)/def2-QZVPP与之类似）
关键词：wB97X-2 D3 def2-TZVPP def2/J def2-TZVPP/C RIJCOSX tightSCF
%method
D3S6 0.547
D3A1 3.520
D3S8 0.0
D3A2 7.795
end
ωB97X-2-D3(BJ)在PCCP, 20, 23175 (2018)的双杂化泛函横测当中表现出众（不过在算分子间弱相互作用的测试上表现一般）。ωB97X-2从ORCA 5.0开始支持，但支持的只是2009年当时提出的原版，其搭配的DFT-D3(BJ)色散校正参数在PCCP, 20, 23175 (2018)的补充材料里才给出，没有在ORCA里内置。因此Multiwfn产生的输入文件里如上所示自动用%method ... end字段补充了色散校正参数。

• RI-revDSD-PBEP86-D4/def2-TZVPP（RI-revDSD-PBEP86-D4/def2-QZVPP与之类似）
关键词：revDSD-PBEP86-D4/2021 def2-TZVPP def2/J def2-TZVPP/C RIJCOSX tightSCF
revDSD-PBEP86-D4在计算有机反应能以及弱相互作用能方面在现有的双杂化泛函里都是顶级的，建议对这类问题使用。

• DLPNO-CCSD(T)/def2-TZVPP with normalPNO and RIJCOSX。关键词：DLPNO-CCSD(T) normalPNO RIJCOSX def2-TZVPP def2/J def2-TZVPP/C tightSCF
CCSD(T)普遍被认为是计算静态相关不是特别强的体系的金标准。DLPNO-CCSD(T)是ORCA中的黑科技，是一种对CCSD(T)的数值近似，可以把原本最多只能算得动不超过30原子的CCSD(T)的方法扩展到好几十甚至上百原子，精度和耗时通过此方法中的一些阈值参数来控制。ORCA里有LoosePNO、NormalPNO、TightPNO三种标准，越往后越贵，但结果也越接近CCSD(T)。LoosePNO就太烂了，不建议用，而NormalPNO较有实用价值。用NormalPNO的时候精度就已经显著超过revDSD-PBEP86-D4了，相对于CCSD(T)的误差通常在1 kcal/mol以内。DLPNO-CCSD(T)的HF计算部分可以用RIJCOSX方法来加速，这是为什么关键词里写了RIJCOSX，并且指定了def2/J辅助基组。def2-TZVPP/C辅助基组是用于DLPNO-CCSD(T)的电子相关部分计算的。顺带一提，在北京科音高级 培训班（http://www.keinsci.com/workshop/KAQC_content.html）中专门有一节“低标度耦合簇方法”很详细讲授DLPNO相关方法以及LNO-CCSD(T)等其它低标度方法。

• DLPNO-CCSD(T)/def2-TZVPP with tightPNO and RIJCOSX。关键词：DLPNO-CCSD(T) tightPNO RIJCOSX def2-TZVPP def2/J def2-TZVPP/C tightSCF
DLPNO-CCSD(T)结合tightPNO的时候，根据J. Chem. Theory Comput., 11, 4054 (2015)的测试（注意此文里有很多在耗时方面的严重误导性说法），与CCSD(T)的误差在1 kJ/mol的程度，几乎可认为没有差别。tightPNO的耗时比NormalPNO通常高几倍。在笔者的Intel 36核机子上，用当前级别计算66个原子的有机体系耗时约8小时，耗硬盘最多时候为120GB。如果你还想要更好的精度，建议将基组提升至def2-QZVPP，但也会贵非常多。如果不用RIJCOSX，即去掉RIJCOSX def2/J关键词，精度会有很轻微改进，但对较大的体系，SCF部分的耗时会增加许多。

• CCSD(T)/cc-pVTZ。关键词：CCSD(T) cc-pVTZ tightSCF
这就是原版的CCSD(T)/cc-pVTZ，没什么好说的。

注：上述几个计算级别中，如果有cc-pVTZ不支持的元素，把cc-pVTZ替换为质量差不多的def2-TZVPP。如果提示某些元素没有对应的辅助基组，把cc-pVTZ cc-pVTZ/JK cc-pVTZ/C替换为def2-TZVPP def2/JK def2-TZVPP/C。

• CCSD(T)-F12/cc-pVDZ-F12 with RI。关键词：CCSD(T)-F12/RI cc-pVDZ-F12 cc-pVDZ-F12-CABS cc-pVTZ/C tightSCF
F12是显式相关方法，可以结合到CCSD(T)、MP2等方法上。CCSD(T)-F12结合cc-pVDZ-F12基组时，耗时显著低于CCSD(T)/cc-pVQZ，但精度则与之接近。因此如果你需要CCSD(T)/cc-pVQZ档次的数据但算得很吃力的话，当前级别是个很好的选择。计算过程中利用RI可以进一步显著节约时间，所以写了/RI。顺带一提，北京科音高级 培训班（http://www.keinsci.com/workshop/KAQC_content.html）专门有一节“显式相关计算”用30多页幻灯片非常完整深入讲F12这类计算。

• Approximated CCSD(T)/CBS with help of MP2 (cc-pVTZ->QZ extrapolation)。关键词：ExtrapolateEP2(3/4,cc,MP2) tightSCF
我在《谈谈能量的基组外推》（//www.umsyar.com/172）中专门介绍过基组外推的概念，这也叫CBS外推，即假定外推到了完备基组极限(CBS)。CCSD(T)档次下最常用的外推就是基于cc-pVTZ和cc-pVQZ来外推，可以免费地得到更高一个档次基组的结果，即大约cc-pV5Z的结果，这种计算方式常见于各种高精度小体系的研究文章当中。直接利用CCSD(T)结合cc-pVTZ和cc-pVQZ能量进行外推的话，后者耗时非常高，往往很难承受，因此可以借助MP2来进行“近似的外推”，即利用MP2/cc-pVTZ和QZ先外推出MP2/CBS相关能，然后CCSD(T)/CBS的相关能就近似可以估计为CCSD(T)/cc-pVTZ + MP2/CBS - MP2/cc-pVTZ，这就避免了直接做非常昂贵的CCSD(T)/cc-pVQZ计算了。这个近似精度很好，一般也就带来0.1 kcal/mol左右程度的误差，笔者在《各种后HF方法精度简单横测》（//www.umsyar.com/378）中专门做过测试，有兴趣可以看看。ORCA提供了ExtrapolateEP2关键词，可以直接实现上述CCSD(T)结合MP2的CBS外推。

其实如果将这里的CCSD(T)改为tightPNO的DLPNO-CCSD(T)，能以相近的精度算明显更大体系，但可惜这没法在目前的ORCA里通过一套关键词直接实现。

• DLPNO-CCSD(T)/CBS with RIJCOSX & tightPNO (def2-TZVPP->QZVPP extrapolation)。关键词DLPNO-CCSD(T) tightPNO Extrapolate(3/4,def2) RIJCOSX def2/JK tightSCF
• CCSD(T)/CBS (cc-pVTZ->QZ extrapolation)。关键词：CCSD(T) Extrapolate(3/4,cc) tightSCF
ORCA里可以用Extrapolate关键词对任意后HF方法进行外推，以上两个级别就用了这点。一个是直接用CCSD(T)，算十个原子都已经极度吃力；另一个是便宜得多但结果也糙一些的DLPNO-CCSD(T) with tightPNO，同时为了节约SCF时间而用了RIJCOSX，之所以改用了def2，是因为辅助基组层面的原因

4.2 激发态计算

• sTD-DFT based on RI-wB97X-D3/def2-SV(P) orbitals
关键词：wB97X-D3 def2-SV(P) def2/J RIJCOSX
%tddft
Mode sTDDFT
Ethresh 7.0
PThresh 1e-4
PTLimit 30
maxcore 6000
end
对于好几百个原子的很大体系的电子光谱计算，往往要算几百个态才够覆盖感兴趣的波长范围。用wB97X-D3/def2-SV(P)对这样大小的体系算单点能往往算得动，但是用TDDFT算这么多态往往太困难。此时可以用Grimme提出的sTD-DFT方法，此方法基于DFT计算求解的轨道，通过对TDDFT矩阵元进行高度近似，只需要非常少的计算量就可以算出大量激发态（sTDDFT与TDDFT的耗时关系有点像半经验方法与DFT的耗时关系）。但代价是精度有所降低，不过对于很大体系的粗放式研究来说够用了。sTD-DFT功能已被内置于ORCA中，上面的关键词就是先做RI-wB97X-D3/def2-SV(P)单点计算，再基于其轨道做sTD-DFT，可以算出激发能为7 eV以内的所有态（这是一般感兴趣的电子光谱能量范围）。用wB97X-D3泛函是因为根据Phys.Chem.Chem.Phys.,16,14408(2014)的测试，将它结合sTD-DFT来用对于大体系（通常有显著电子共轭）结果较理想。PThresh和PTLimit用于控制纳入考虑的组态函数的范围，不写它们的话会考虑所有组态函数，这对于大体系可能耗时较高，当前的PThresh和PTLimit设置可以避免这个问题，对精度的牺牲可忽略不计。注意sTD-DFT目前只能对单个结构做电子激发计算，而不能用于激发态优化、振动分析等目的。顺带一提，北京科音高级 培训班（http://www.keinsci.com/workshop/KAQC_content.html）专门有一节“使用sTDA方法快速计算大体系电子光谱”用约40页幻灯片非常全面完整详细讲这类计算。

• TDDFT RI-PBE0/def2-SV(P)
关键词：PBE0 def2-SV(P) def2/J RIJCOSX tightSCF
%tddft
nroots 10
TDA false
end
这是用PBE0结合大小与6-31G*差不多的def2-SV(P)做TDDFT计算。PBE0是TDDFT计算激发态常用的泛函，但碰见大共轭体系建议改为CAM-B3LYP等，详见《乱谈激发态的计算方法》（//www.umsyar.com/265）。def2-SV(P)对于计算中、大体系的电子光谱是很适合的，既不贵，精度也基本够用，如果有余力且希望结果更好建议把def2-SV(P)改为def2-TZVP(-f)或def-TZVP（关键词写为TZVP）。关键词里nroots 10代表计算10个激发态，当然需要根据实际情况进行调节，同Gaussian的TDDFT计算中的nstates，详见《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（//www.umsyar.com/314）。用tightSCF是为了降低数值层面的误差。

• TDDFT RI-DSD-PBEP86/def2-TZVP
关键词：DSD-PBEP86 def2-TZVP def2/J def2-TZVP/C RIJCOSX tightSCF
（用TDDFT RI-RSX-QIDH/def2-TZVP则把上面的DSD-PBEP86改为RSX-QIDH）
%tddft
nroots 10
TDA false
end
ORCA是为数不多的支持双杂化泛函下做TDDFT的程序，此级别精度比普通泛函做TDDFT更高，但代价是耗时也明显更高。计算时应当同时指定/C辅助基组。根据JCTC, 17, 4211 (2021)的测试，在ORCA能用的双杂化泛函中DSD-PBEP86算局域激发和分子内电荷转移激发最好，在JCTC, 18, 1646 (2022)的图5和表2中体现RSX-QIDH是算分子间CT激发的首选，所以Multiwfn给了两个选项。用上述关键词用对一个47原子的有机体系做TDDFT计算，笔者在Intel 2*2696v3 36核机子上用不到20分钟算完，可见耗时毫不夸张，算稍大一些的体系也能算得动。

• EOM-CCSD/cc-pVTZ
关键词：EOM-CCSD cc-pVTZ tightSCF
%mdci
nroots 3
end
EOM-CCSD计算激发态的精度算是很不错了，比TD-双杂化精度更高，也更稳健和普适，但耗时颇高，而且随算的态数增加总耗时迅速增加，因此主要适合的是小体系的低阶激发态的较高精度的研究。RI在EOM-CCSD计算时派不上用场，和Gaussian相比也没有显著优势。

• STEOM-DLPNO-CCSD/def2-TZVP(-f)
关键词：STEOM-DLPNO-CCSD RIJCOSX def2-TZVP(-f) def2/J def2-TZVP/C tightSCF
%mdci
nroots 3
end
STEOM-CCSD方法的初衷是降低EOM-CCSD的耗时，但不应将STEOM-CCSD视为是EOM-CCSD的近似，而应当视为两种不同方法。二者精度差不多，而STEOM-CCSD对CT激发的结果往往比EOM-CCSD更好。EOM-CCSD明显算不动的体系STEOM-CCSD照样算不动，结合像样的基组时一般顶多也就用于二十多个原子。DLPNO技术与STEOM-CCSD的结合诞生的STEOM-DLPNO-CCSD则可以用于大得多的体系，在def2-TZVP(-f)这样的基组下甚至对超过60个原子的体系也能算（但注意超级耗硬盘，没有好几个TB的剩余空间就别指望了）。

Written by Tian Lu (sobereva@sina.com), 2019-May-20

 Beijing Kein Research Center for Natural Sciences (http://www.keinsci.com) 

1 Introduction

In this tutorial, I will show how to very easily using ORCA and Multiwfn to simulate UV-Vis and electronic circular dichroism (ECD) spectra for a typical organic system alanine in water environment. The ORCA is a free quantum chemistry program, it can be downloaded via https://orcaforum.kofo.mpg.de/app.php/portal. Although Multiwfn is a code aiming for wavefunction analysis, it also has a powerful module used to plot various kinds of spectra based on output file of Gaussian, ORCA, xtb, sTDA or plain text file. Multiwfn can be freely obtained at //www.umsyar.com/multiwfn. All the ORCA and Multiwfn supports Windows, Linux and MacOS platforms.

In this tutorial, I assume that you are using Windows platform. The version of ORCA is 4.1.1, the version of Multiwfn is 3.7(dev) updated on 2019-Jul-14 (do not use older version). The overall computational cost should be less than 10 minutes for an ordinary Intel quad-core CPU.

All files involved in this tutorial can be downloaded at //www.umsyar.com/attach/485/file.rar.

2 Optimization of ground state

Both UV-Vis and ECD are absorption spectra, commonly it is assume when a molecule undergo electronic excitation, the geometry is at minimum point of potential energy surface of ground state. Therefore, the ground state geometry should be optimized first. Since the conformational space of alanine is not complicated, construction of the alanine structure is relatively arbitrary. It is important to note that under water environment the alanine is a zwitterion. The initial geometry (initgeom.xyz) I built is shown below:

Now we generate input file of optimization task of ORCA. The easiest way in my opinion is using Multiwfn to do this. Double click icon of Multiwfn.exe to boot up Multiwfn, then input below commands. The texts behind // are comment.
initgeom.xyz   //Other formats are also supported, e.g. mol, mol2, pdb, gjf, fch, molden...
100  //Other functions (Part 1)
2   //Exporting files or generating input file of mainstream quantum chemistry codes
12  //Generate ORCA input file
S0_opt.inp  //The path of the input file to be generated
-1   //Enable using SMD solvation model
[Press ENTER button to use default solvent (water)]
0  //Select type of task
2  //Optimization
1  //B97-3c. This is an economical level but able to provide satisfactory geometry

Now S0_opt.inp has been generated in current folder. Open it by text editor, change the value behind "nprocs" to the actual number of physical CPU cores of your machine, also properly set the value behind "maxcore", which controls the maximal amount of memory can be utilized per CPU core (in MB).

Move the S0_opt.inp to an empty folder, open console window of your system and enter this folder, then input command such as D:\study\orca411\orca S0_opt.inp > S0_opt.out to conduct the calculation, where D:\study\orca411\orca is the absolute path of ORCA executable file in my machine.

After calculation, from the output file S0_opt.out you can find the optimization converged after 11 steps without any error. The S0_opt.xyz in current folder is the optimized geometry, while S0_opt.trj contains every structure of optimization trajectory (If you want to visualize the trajectory, you can manually change the suffix to .xyz, and then drag this file into VMD, which is an excellent visualization program and can be freely obtained via http://www.ks.uiuc.edu/Research/vmd/).

Use Multiwfn to load the S0_opt.xyz, and then enter main function 0 to visualize the geometry, it can be seen that the optimized geometry fully meets our expectation, as shown below

3 Calculation of excited states

Next, we perform excited state calculation using the popular TDDFT method. According to many benchmark articles, for TDDFT calculation of such a small molecule, PBE0 functional is very suitable (however, for large conjugated systems such as most organic dyes, in particular the states with strong charge-transfer character, CAM-B3LYP and wB97XD often work much better)

Return to main menu of Multiwfn, then input
100  //Other functions (Part 1)
2   //Exporting files or generating input file of quantum chemistry codes
12  //Generate ORCA input file
TDDFT.inp  //The path of the input file to be generated
-1   //Enable using SMD solvation model
[Press ENTER button to use default solvent (water)]
22  //TDDFT task at PBE0/def2-SV(P) level with RI acceleration technique

The TDDFT.inp has been generated in current folder, do not forget to properly change the "nprocs" and "maxcore" parameters. Note that in this file the "nroots" is set to 10, namely ten lowest excited states will be calculate. For plotting absorption spectra purpose of small systems like alanine, this setting is adequate, however, if your system consists of much more atoms, "nroots" should also be accordingly increased, otherwise the finally simulated spectra cannot fully cover wavelength range of common interest (e.g. >250nm). In addition, since this system is fairly small, in order to reach higher accuracy, we can use basis set better than the def2-SV(P), therefore we replace the "def2-SV(P)" keyword with "def2-TZVP(-f)", which is a good basis set with size similar to 6-311(2d,p).

Run command like this: D:\study\orca411\orca TDDFT.inp > TDDFT.out

4 Simulating UV-Vis spectrum

From the TDDFT output file, you can find excitation energies and oscillator strengths under "ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS". Based on these data, Multiwfn is able to plot theoretical UV-Vis spectrum. Now we do this.

Boot up Multiwfn and input
TDDFT.out
11  //Plot spectrum
3  //UV-Vis
0  //Show the spectrum on the screen

You can use abundant options in the interface to gradually adjust plotting parameters and then replot the graph until you are satisfied. Please check Section 3.13 of Multiwfn manual on the detailed introduction of the spectrum plotting module and related theoretical background. Section 4.11 provided examples of plotting other kind of spectra.

NOTE: If Multiwfn crashes before entering the plotting interface, it is possible that your ORCA output file is in Unicode encoding, while Multiwfn only supports parsing text file in ASCII encoding. Please check #6 in this post on how to solve this problem: //www.umsyar.com/wfnbbs/viewtopic.php?id=213.

5 Simulating ECD spectrum

Multiwfn is also able to simulate ECD spectrum based on excitation energies and rotatory strengths in the ORCA TDDFT output file. Now we do this.

Input below commands in the Multiwfn window:
-3  //Return to main menu
11  //Plot spectrum
4  //ECD
0  //Show the spectrum on the screen

6 Some worth noting points

There are several points I want to mention.

For most neutral system, geometry optimization can be carried out without solvation model to decrease time cost. However, for systems whose local region shows evident ionic character, such as the alanine in water environment, the solvation model should also be applied in optimization stage. For calculation of excited states, solvation model should always be employed to mimic real environment because its impact on electronic excitation is large, in particular for polar solvents.

ECD spectrum is very sensitive to conformation. The alanine under our present study does not have more than one accessible conformation, therefore we can safely use single structure to simulate the spectra. However, if there may be multiple thermally accessible conformations under present condition (this is usually true for flexible molecules containing rotable bonds at room temperature), commonly you should optimize structure and calculate free energy for each one, then use Boltzmann relationship to estimate occurrence percentage, then perform electron excitation calculation for all conformations having probability higher than e.g. 5%, and finally use Multiwfn to plot conformationally averaged ECD spectrum, as illustrated in Section 4.11.4 of Multiwfn manual.

BTW 1: To obtain all possible conformations of a highly flexible molecule, commonly you need to use conformational search software, there are many choice. Among which, one of best choice is the Molclus code developed by me, it is free and very flexible, the official site is http://www.keinsci.com/research/molclus.html (English version of manual will be available later).

BTW 2: To characterize and understand nature of electron excitations, Multiwfn provide numerous functions, see Section 4.A.12 of the manual for an overview.

BTW 3: There is a video tutorial illustrating how to use ORCA quantum chemistry program in combination with Multiwfn, OfakeG and GaussView to realize very common calculation tasks and analyses for a simple organic molecule methanol, taking a look at it is suggested: https://youtu.be/tiTmTbtbtig.