辨析计算化学中的任务类型和理论方法
辨析计算化学中的任务类型和理论方法
Distinguishing tasks and theoretical methods in computational chemistry
文/Sobereva@北京科音
First release: 2023-Aug-5 Last update: 2024-Feb-23
0 前言
在网上答疑时，偶尔看到有初学者搞不清楚任务类型和理论方法，比如今天有人在思想家公社3号群问“md，aimd和qmmm的区别是啥啊？什么需求下会用到这三种呢？”，这三个词明显都不是一个层面的东西。此文就把计算化学中的任务类型和理论方法的关系明确一下，简明扼要地介绍一下相关基本概念，做一个科普，以令计算化学零基础者一次性搞清楚它们的关系。更多的计算化学的总览性的知识在《北京科音初级量子化学培训班》（http://www.keinsci.com/workshop/KEQC_content.html）里有介绍，此培训也是从零开始入门量子化学研究的最好方式。
1 任务类型
有N个原子的非线型的体系有3N-6个内坐标描述原子间相对位置关系，在Born-Oppenheimer近似下体系的能量是依赖于内坐标的，也即能量是个3N-6维的函数，这被称为势能面（potential energy surface, PES）。计算化学领域有很多常见的任务（task）都是基于势能面做的，即有了势能面信息（能够求得势能面上任意位置的能量及其导数）就能做这些任务。下面罗列一下常见的这种类型的任务：
• 优化极小点：平时说的几何优化（geometry optimization）一般也是指这个。此任务从一个给定的初猜结构开始，根据特定算法去寻找与之最近的势能面上极小点的精确位置。在分子动力学程序如GROMACS里这种任务也被叫做能量极小化（energy minimization, em），只不过实际目的不一样，em的目的主要是在动力学模拟之前释放体系中可能存在的较大斥力（自行构建的初始模型里往往有原子间距离太近、斥力太大），免得动力学模拟一开始由于过大的原子加速度造成过大的速度而导致动力学模拟异常或崩溃。
• 优化过渡态：从一个给定的初猜结构开始，根据特定算法去寻找与之最近的势能面上的一阶鞍点的精确位置。过渡态可以视为反应过程中最有代表性的一个结构，可以由此判断或验证反应机理，利用它和相邻极小点的能量求差可以得到反应势垒，对于讨论反应难易有关键意义，见《谈谈如何通过势垒判断反应是否容易发生》（http://sobereva.com/506）。优化过渡态的算法介绍见《过渡态、反应路径的计算方法及相关问题》（http://sobereva.com/44）。
• 产生反应路径：用于把基元反应对应的势能面上两个相邻极小点之间最容易相互转变的路径产生，也相当于得到了一个轨迹并可以观看对应的结构变化过程，过渡态是此路径上能量最高点。在质量权重坐标下产生的反应路径称为IRC（intrinsic reaction coordinate），不考虑质量权重时产生的一般称为MEP（minimum energy path）。具体算法很多，见《过渡态、反应路径的计算方法及相关问题》（http://sobereva.com/44）。反应路径可认为是实际化学反应过程最有代表性的路径，自然对于理解反应机理有重要意义，还可以用Multiwfn对反应路径上各个点做波函数分析考察反应过程中电子结构的变化以获得更深入的信息，例如《通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程》（http://sobereva.com/200）。
• 振动分析：通常是在势能面上的驻点（所有原子受力都为0的点）做的，用于得到相应结构下的振动频率，可以用来计算热力学量，公式看《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）介绍的Shermo程序手册的附录部分，还可以用来得到振动能级、振动波函数、检验几何优化的准确性（根据虚频判断）、绘制振动光谱（参看《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》http://sobereva.com/224）。
• 分子动力学（molecular dynamics, MD）：上面的任务都是“静态”的任务，即不含时间这个维度。而分子动力学则引入了时间这个维度，可以模拟体系状态随时间的变化，能研究的问题与那些静态的任务有极大的互补性，在北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/workshop/KGMX_content.html）里面有全面的介绍。分子动力学模拟过程相当于体系在势能面上不断运动，也相当于是对势能面的一种采样方式。虽然其名字叫分子动力学，但研究对象并不限于分子，比如无机固体、金属团簇的动力学模拟也可以叫分子动力学。最常见的分子动力学的实现形式是BOMD（结合量子化学方法时另有CPMD、ADMP等），也相当于根据原子的位置、速度和受力按照经典牛顿运动方程演化原子的坐标和速度。还有其它一些特殊的动力学形式，如朗之万动力学、耗散粒子动力学等，如今用得很有限。
• 蒙特卡罗：和分子动力学并列的另一种常见的对势能面采样的方法，适用场合相对更少，主要是在多孔材料吸附小分子、计算气液共存曲线等问题上用得较多。不需要像MD那样计算受力和速度信息，也没有明确的时间的概念、无法像MD那样严格地考察含时演化。
• 构型/构象搜索：势能面上往往有很多极小点，能量最低的那个是全局最小点，可以视为体系最稳定的构型或构象对应的结构，其它的极小点对应亚稳的构型/构象。前面说的几何优化只能收敛到与初猜结构最近的极小点，显然这未必是最小点。如果想确保得到最小点（全局优化，global optimization），或者能量最低的一批构型/构象，就需要做构型/构象搜索，方法很多，比如免费的molclus（http://www.keinsci.com/research/molclus.html）就是很常用的实现工具。
2 理论方法
上一节介绍的那些任务在进行过程中，至少需要计算体系在不同结构下的能量。有的还需要计算受力（势能面对核坐标一阶导数矢量，或者说梯度，的负矢量），如几何优化、分子动力学。有的还需要计算能量对核坐标的更高阶导数，比如振动分析至少需要算二阶导数（Hessian矩阵由之构成）。怎么获得各种任务中涉及的能量及其对核坐标的导数，就是理论方法层面的问题了。下面简要列举一些常见的理论方法：
• 分子力场（molecular forcefield）：也称为经典力场（classical forcefield）或分子力学（molecular mechanics, MM），或者经验势函数等等。虽然通常带着“分子”俩字，但实际中不限于用于分子体系。这类方法一般使用简单的模型描述体系（通常一个原子作为一个粒子，电子与原子核不分开考虑），并用形式简单的数学函数（势函数）描述原子间相互作用。例如计算原子间静电相互作用时，分子力场普遍是给每个原子核位置分配一个点电荷（数值等于原子电荷），然后基于库仑公式计算。由于分子力场的形式简单，因此计算耗时极低。分子力场的一个关键的局限性是绝大多数分子力场都无法描述化学反应，因为成键关系在计算从始至终是固定不变的，是一开始就定义好的。描述化学反应问题主流的是后面提及的量子化学类型的方法，或者反应力场（远比普通分子力场要贵、支持的程序要少）。分子力场另一个关键局限性是依赖于大量参数，一方面影响模拟精度，一方面决定能描述的体系，寻找和构建合适的参数往往不是易事。
• 机器学习势（machine learning potential）：基本思想是通过机器学习的思想构造依赖于原子坐标的分子描述符（如距离矩阵、能量矩阵等）与能量之间的抽象关系，这种关系比上述传统的分子力场用的势函数形式明显更为复杂。
以下提及方法都是量子化学（quantum chemistry, QC）范畴的，也可以叫基于量子力学（quantum mechanism, QM）的方法。
• 密度泛函理论（DFT）：是目前在量子化学、第一性原理计算领域用得最普遍的理论方法，因为其性价比非常高，交换-相关泛函（影响DFT计算精度的决定性因素）选择得当的话可以得到很不错的结果，足够满足绝大多数研究需要，而且精度显著强于分子力场（除非力场是对某类体系专门训练得特别精妙的力场），普适性强得多得多，但耗时高于分子力场N个数量级，体系越大N越大。DFT应用到实际问题中实际上分为Kohn-Sham DFT（KS-DFT）和orbital free DFT两种形式，前者是绝对的主流，计算过程牵扯到轨道波函数，本文以及大家平时说的DFT计算都是用的KS-DFT形式；而后者非主流，不牵扯到轨道，应用的局限性非常大，只有很少数程序支持。
• Hartree-Fock（HF）：在DFT兴起之前用得很普遍，如今已经完全过时。耗时不比DFT显著更低（DFT计算时利用RI技术时反倒还可以明显更快），而精度差得多。
• 后HF：是一大类方法的统称，如CCSD(T)、MP2、CISD等。HF方法精度烂在于同时没考虑到动态相关和静态相关。后HF在HF基础上额外把动态相关在一定程度上考虑了进来，但耗时也高了很多。动态相关、静态相关的相关概念很有学问，在北京科音基础量子化学培训班（http://www.keinsci.com/workshop/KBQC_content.html）里有详细讲授，这里就不多提了。
• MCSCF：主要弥补HF缺乏对静态相关的考虑。但由于几乎没有或很少考虑动态相关，所以结果也说不上好。MCSCF最常见的实现是CASSCF。
• 多参考方法：在MCSCF基础上进一步把动态相关考虑进来，精度整体很好，普适性很强，但使用相对复杂且昂贵。典型的实现如CASPT2、NEVPT2、MRCI。
• 半经验：一大类方法的统称，都是对HF的简化以巨幅降低耗时，耗时只是HF的微小零头，整体精度也有不小降低，而且引入了依赖于元素的参数而只能用于有限的元素。典型的如AM1、PM3、PM6。一些较复杂的机器学习势的耗时已经追平了半经验。
• GFN-xTB：相当于纳入了一定DFT思想的半经验级别的方法，耗时和半经验相仿佛，而整体精度和可靠性>=主流的半经验方法，但跟DFT比还有很大差距。GFN-xTB从2017年开始发展，和历史更早的DFTB有极大的共性。
前面介绍的HF、后HF、MCSCF、多参考方法里面都完全没有引入经验参数（极个别的除外，诸如CASPT2可以涉及shift参数），计算中的参数只依赖于物理常数，故它们被称为“从头算（Ab initio）”类型的方法。分子力场、半经验、GFN-xTB由于引入了大量经验参数，因此明显不属于从头算。DFT算不算从头算模棱两可。如今尚未找到的理论上严格的交换-相关泛函是不含任何参数的，但如今只有近似的交换-相关泛函被提出，其中虽包含经验参数，但整体相对较少，因此争论DFT是否属于从头算没意思，只是个说法问题。有的泛函的经验性较小，有的经验性则较高。如今也有的文章将HF、后HF、MCSCF、多参考方法等物理思想完全基于波函数的方法统称为wavefunction theory (WFT)，在讨论时使用WFT一词时DFT就不会和那些方法混淆了。
还有一个词叫“第一性原理（first-principle）”，其字义和从头算一样，但“第一性原理计算”如今的主流指代是“量子化学方法研究周期性体系”，不仅包括严格的从头算方法，也包括DFT，因为这类计算最常用的就是DFT。DFTB虽然依赖很多参数，但DFTB计算周期性体系也往往被不计较地算作第一性原理计算。
由于有了“第一性原理计算”这个专门描述把量子化学方法用于周期性体系的计算的词，因此平时说的“量子化学计算”、“从头算方法计算”一般是指计算分子、团簇这样孤立（非周期性）体系的情况。那些通过量子化学方法主要研究孤立体系的程序被称为量子化学程序，如Gaussian、ORCA、GAMESS-US等，而通过量子化学方法主要研究周期性体系的程序被称为第一性原理程序，如CP2K、Quantum ESPRESSO等。注意“量子计算”和量子化学计算又不是一码事，前者强调的是利用在量子计算机上实现的量子算法，可以但绝不限于做量子化学问题的研究。
接下来再说很知名的QM/MM。这个方法是指使用量子化学方法以量子力学（QM）的思想描述体系的一部分，用分子力场以分子力学（MM）的思想描述体系的另一部分，并且恰当考虑两部分之间的耦合。对于很大的体系，这显然比起全都用量子化学方法描述要便宜得多得多（MM部分计算耗时相当于QM部分的九猪一毛），但代价是MM区域描述的精度比QM区要差，而且绝大多数力场不能描述这部分发生的成键/断键过程。因此，必须恰当划分QM和MM区，让最需要精确描述的部分或者涉及化学反应的原子纳入到QM区中，而起到较次要作用的“环境”原子纳入到MM区中。
还有一个词是ONIOM，也是划分不同区域，不同方法着重描述不同部分，它可以把任意两种理论方法相组合。如果把量子化学方法和分子力场相组合，特称为ONIOM(QM:MM)，其用处和QM/MM类似但实现不同，描述时绝对不能搞混。一些相关讨论见《要善用簇模型，不要盲目用ONIOM算蛋白质-小分子相互作用问题》（http://sobereva.com/597）。
再说一下对电子激发态的描述。前面介绍的半经验、HF、后HF、DFT等方法通过delta-SCF方式可以算激发态，但不够普适。量子化学领域有很多专门的方法可以计算激发态，包括计算激发态的能量及其导数，如最常用的TDDFT、过时的CIS、较高精度的EOM-CCSD等，见《乱谈激发态的计算方法》（http://sobereva.com/265）。几乎所有分子力场都是描述基态电子态的，但如果专门对激发态势能面拟合力场参数，也可以用力场描述激发态，但已有的这样的力场都是针对极其具体的体系的激发态拟合的，不具备普适性。机器学习势同理。QM/MM也可以描述激发态，但仅限于电子激发完全发生在QM区的情况。
再简单说一下和实际计算关系密切的基函数（basis function）的概念。本节介绍的理论方法，凡是基于量子理论思想提出的（也就是前面说的那些里面除了分子力场和机器学习势以外的方法），在实际数值求解的过程中一般都要涉及分子轨道，关于分子轨道的概念详见《正确地认识分子的能隙(gap)、HOMO和LUMO》（http://sobereva.com/543）。绝大多数计算程序中都是把分子轨道展开成基函数的线性组合来描述的。最常见的基函数一类是原子中心基函数（如Gaussian等大多数量子化学程序以及CP2K等部分第一性原理程序用的高斯型基函数），其中心一般位于原子核，还有一类常见的基函数是平面波基函数，是大多数计算周期性体系为主的第一性原理程序用的，它的分布覆盖整个被计算的晶胞。基组（basis set）是对于原子中心基函数而言的，例如6-31G*、def2-TZVP、cc-pVDZ等都是很常见的基组，它定义了实际计算时对各种元素原子具体用多少、什么参数的基函数。做HF、DFT、后HF、MCSCF、多参考等方法计算时都需要告诉计算程序用什么基组，基组质量越好，也就是越接近于所谓的完备基组极限，这些方法自身的精度发挥得就越充分，但代价就是耗时越高。一个好方法配一个烂基组，以及一个烂方法配一个好基组，结果都不理想，必须好方法配好基组才能得到较准确结果。半经验、GFN-xTB方法计算时也利用基函数，但是用什么基函数是这类方法自身定义的，就不需要而且也不能再指定基组了。北京科音初级量子化学培训班（http://www.keinsci.com/workshop/KEQC_content.html）和北京科音基础（中级）量子化学培训班（http://www.keinsci.com/workshop/KBQC_content.html）都专门有一节对基组做详细介绍，后者讲得明显更深、更广。
3 任务类型与理论方法的结合
第2节介绍的各种理论方法原则上都可以用于第1节介绍的各种任务。例如，分子力场、DFT、QM/MM等理论方法全都可以用于几何优化、振动分析、产生反应路径等任务。这些组合在程序实现上一般没有什么需要特殊考虑的，一个程序支持什么理论方法，就能基于它们做什么依赖于势能面的任务。
搞明白了前面介绍的名词，从头算动力学（ab initio molecular dynamics, AIMD）和第一性原理动力学（first-principle molecular dynamics, FPMD）顾名思义自然就知道是描述什么计算了。这两个词其实没有严格的区分，实际中经常混淆使用，但从严格的习俗来说，AIMD和FPMD分别最适用于描述使用量子化学方法对孤立体系和周期性体系做分子动力学模拟。所以北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）里面讲CP2K做基于GFN-xTB和DFT的分子动力学的地方标题用的是第一性原理动力学，但若说成AIMD也完全可以，不会造成什么误会。
由于基于分子力场做分子动力学的文章数目和流行程度远高于非常昂贵的AIMD/FPMD，因此平时说的分子动力学模拟默认是指的基于经典力场的分子动力学，这也往往被叫做经典分子动力学。但要注意“经典”这个词在不同语境下有不同含义，有的场合下是特指把原子核当做经典粒子来考虑，此时基于BOMD形式做基于量子化学方法的分子动力学也是经典分子动力学。而与之相对的是“量子动力学”，强调把原子核以量子力学方式描述，原子核有对应的波函数。显然，此处“量子动力学”和“基于量子化学做分子动力学”又不是同一个概念。
QM/MM MD指什么也不言而喻，显然是用的QM/MM理论方法做分子动力学任务，这结合了QM/MM对体系描述的特点和分子动力学研究问题的能力。显然这比所有原子都用QM描述（相当于AIMD/FPMD）要便宜得多得多，当然代价就是MM区域在动力学过程中描述得相对较糙，而且预期不会发生化学反应，同时还得想办法去找力场参数。因此能否用QM/MM MD代替AIMD/FPMD当然得根据被研究的体系和问题判断，盲目瞎用纯粹是自讨苦吃（吐槽一下，很多初学者没有基本理论常识，一看文献里有QM/MM就觉得好像高大上似的，就整天想着QM/MM而不知道其局限性，太naive了）。
文献里的许多用词并不标准，必须结合语境和常识理解说的是什么。比如可能有的文章没明确用QM/MM MD这个词而只说了QM/MM，若实际上在此基础上跑了动力学，当然就是QM/MM MD。再比如有的文章表面上说是做的AIMD，但从计算细节描述来看作者把一部分当做MM来描述以节约时间，显然实际上也是QM/MM MD。理解名词必须有基本的应变能力。
如果用的理论方法描述的是激发态，那么前述各种任务就是对激发态做的。比如用TDDFT做激发态分子动力学，那就是动力学的每一步用的原子受力都是用TDDFT算的激发态的受力。顺带一提，有些文献里声称做的是激发态动力学研究，但其实内容里根本没做分子动力学，只是对激发态计算了势垒，文中这么说只是由于势垒和研究反应速率的反应动力学领域里的反应速率常数密切相关。另外值得注意的是研究激发态动力学往往需要考虑内转换、系间窜越等势能面切换的过程，考虑这些需要做非绝热动力学，需要额外的算法，这方面的坑很深。所以并不是一个程序支持了计算激发态的理论方法和分子动力学模拟后就能完美、严格地做激发态动力学研究相关问题，情况远没那么简单。