半经验与近似DFT.docx

组里面有个老兄专搞SCC-DFTB。有一次我用SCC-DFTB算了个东西，作报告时称SCC-DFTB为半经验DFT理论，他不同意，说应该是近似DFT理 论， approximated DFT。到底“半经验”和“近似”有什么区别，今天就来梳理一下量子化学里都有哪些简化的招数，特别是SCC-DFTB是怎么回事。纯属个人理解，如有错误，请予指正。1. 什么是从头算？从头算包含了哪些近似？量子化学包括哪些简化招数呢？还得先从HF（ab initio从头算理论）说起。首先，量子理论的本质是态叠加理论，或不确定性。即只能知道它有些确定的能级状态（可观测到的态，又叫本征态），但不知道具体在哪个状态，只 知道它是这些状态的某种叠加（混合态）。量子化学计算的是电子在原子核外的排布状态。电子是量子，因此首先把电子所有的确定性状态，即轨道，都找出来，然后根据叠加态原理线性组合成原子轨 道，原子轨道再叠加成分子轨道。用来表示轨道的数学函数，称为基函数，有两个要求，正交性和完备性。正交性是指状态之间无关（确保是本征态），完备性是说要尽 可能找到所有状态，越完备其叠加态描述越准确。（QM计算中基函数用Slater函数，每个slater基函数再由多个gaussian基函数表示，后者又称基 组。基组越大越精确就是指完备性）。然而，实际计算中，不可能包括所有的轨道，一般只考虑占据轨道；对于每个占据轨道，也只使用有限个基函数来表示。这就是状态完备性近似。除了无法使用完全基组造成的“状态完备性不足”这个近似，还有一些其它近似。第一个近似，是假定原子核固定，原子核的量子效应不考虑。这就是Born-Oppenheimer近似。这个容易理解，因为原子核比电子重的多啊，运动也比电子慢的多，所以假定原子核每动一下，电子都可以迅速找到能量最低的排布方式，紧贴原子核的基态势能面运动，而不会被抛到高能轨道态。所以BO近似又叫绝热近似。但是，在实际运动中，电子并不总是能跟上原子核的运动，有个滞后或能级波动，所以这是一个近似。第二个近似，是从头算中的HF中心场近似。这是个从头算理论误差的关键。本来，为了表示电子-电子之间的相互作用对轨道的影响，就需要依次计算某个电子与其它电子之间的两两作用（量子交换作用和库仑排斥作用），但为了节约时间，实际直接计算某个电子与所有其他电子的总作用，相当于该电子在其它电子形成的平均场下的排斥作用。这样一来，我们就能把电子的单个哈密顿一一分开 来，得到每个电子的波函数，并乘起来得到总波函数。然后自洽求解薛定谔方程，最后得到对应能量最低的系数矩阵，即态叠加系数。这个近似把复杂的两两作用转 化成单体和环境的作用，减少了计算量。然而，中心场近似是平均效应，忽略了“在电子运动的任一瞬时，不允许两个电子同时出现在同一位置”这一禁忌，称为电子相关效应，所以纯HF计算得到的比实际能 量要高(每电子20kcal/mol)。电子相关能占到总能量up to 5%左右，因此需要用各种方法修正。怎么修正呢？要么，你去算电子两两作用，这种计算对于小的体系也很慢的；要么，用变分法或微扰法。比如，因为HF忽略了 瞬时排斥，电子必然排的密了，那我构造一些比较松弛的态，把这些松弛的态的可能贡献包括进来（松弛态中排斥能较低嘛）。怎么构造松弛态呢？就是某个电子从当前占据态跳到外层某个未占据的激发态的状态。这就是变分法的思想。包括进来的松弛态越多越精确。另外一个电子电子相关方法是微扰，通过能量对电子状态变化的泰勒展开进行修正。我们知道HF状态并非能量最低点，因此，通过对当前势能对电子状态变化的多级展开，可以从数学上近似确定一个附近的能量最小点（电子状态还是通过构造那些松弛态得到）。二级微扰展开是MP2,三级是MP3.比较而言，变分法是一种物理方法，能保证最终得到的能量不低于真实态能量，而微扰法为数学方法，结果依赖于你当前HF状态的准 确性，即如果当前HF状态描述的能量-状态的导函数性质不正确，可能得到糟糕的结果。变分和微扰都是对近似计算的补救措施，目的是为了获得电子相关。这些计算极为缓慢和昂贵。2. 什么是DFT？DFT包含哪些近似？我们知道DFT计算的跟HF差不多快，却能大致包括电子相关，为什么会有这样的好事儿呢？首先说，HF与DFT本质一样是量子 理论，因此都需要构建本征态，并进行叠加计算。HF是所有电子的尽可能完备的各轨道状态的波函数，DFT呢，叫Kohn-Sham轨道。无论何种轨道，都遵循越 完备越精确的道理；不同的是，它们的哈密顿（动能+势能）求法却不一样： HF的势能是从电子-平均场作用一一计算而来，可后续处理相关；而DFT的电子间势能作用却是直接定义为电子密度的某个关系，这样在自洽中，直接由密度一步跳到包括 了排斥能-交换相关能的结果。那么它的准确性，就主要维系于这个“密度-能量”的关系了。这个是个什么样的能量-密度关系呢？ 除了动能，和中心场近似下的轨道间作用势能，DFT还包括一个交换相关泛函XC，它是计算相关效应的关键。但是这个“能量-交换相关作用”的关系没人能说的清，因此，一般把它处理成各点电子密度的多级展开，即把能量展开为电子密度的函数（以0密度为参考态），由电子密度得到动能，排斥能，交换相关能。交换相关能对密度的一级展开，是LDA，二级展开是GGA, 三级展开是meta-GGA， 展开级数越高，越能描述能量对电子密度变化的复杂关系。其实这个“密度-能量”关系有点像个黑箱。这个黑箱的构造可以由高等级的量化计算去经验的寻找。不幸的是，这种关系往往因为体系的不同而变化。比如固体物理学家爱用的PBE算固体行，算气态分子就不行，生成焓要差上几十个千卡每摩尔，所以搞化学的要用BLYP或XLYP(Xin Xu, 2004)，这是因为因为不同体系中存在着复杂的效应嘛。此外，不同泛函构造思路以及参数化都影响它的质量。3. 什么是半经验？半经验有哪些简化的招数和提高准确度的招数？在介绍DFT的简化之前，我们先谈谈从头算的简化招数，即半经验方法。这种方法在生物分子模拟领域日渐广泛应用，我们旁边一个课题组就是用半经验方法做蛋白模拟，以研究离子通道中存在的可能的化学反应(QM/MM)。流行的半经验主要有AM1, PM3和PM6，特别是PM6, 准确性有较大提高。这些半经验的算法，简化途径基本有三条：第一条，是只让原子价层轨道参与分子轨道自洽计算，内层轨道直接用一个势函数来表示；第二条，原子价层电子往往用最小基组表示。本来在正常 的HF计算中，一个s价层电子轨道，可以选择2个slater函数（double zeta），或者3个slater函数描述(triple zeta)，每个salter再有多达16个gaussian基函数拟合。但在PM6中，每个价层电子只有一个slater基函数描述，每个slater 基函数只有13个gaussian基函数拟合，这叫STO-3G最小基组近似。第三条近似途径，是把用STO-3G表示的外层原子轨道冻结，即不计算轨道间的作用，没有作用则原子轨道保持不变，相当于冻结了。这些轨道间的作用一般是通过重叠积分计算的。同一原子上的轨道间作用用单原子重叠积分，不同原子的轨道间作用叫双原子重叠积分。忽略掉双原子积分(Neglect of Differential Diatomic Overlap), 叫NDDO；忽略掉所有的重叠积分，叫CNDO(C - Complete Neglect DO), 部分忽略叫INDO(I - intermediate)， 通过拟合分子性质提高准确性叫MNDO, (modified NDO); 通过拟合光谱等实验数据改进氢键等分子间作用性质叫AM1(=MNDO加1到7个Gaussian functions, 改进吸引和排斥)；如果是MNDO只加2个gaussian functions，叫PM3（但PM3同时还改进参数化过程)； PM3参数进一步对生物和过渡金属体系优化叫PM6。上述这些近似都是某种程度地减少变量，节约计算时间。如果需要提高准确性，则研究缺少 了哪些效应，然后增加一些函数（变量）来拟合完整的势能面。一减一加，倒腾。半经验算法背离了“状态越完备，结果越精确”的量子叠加计算原则，当然要付出代价，所以准确性较低，甚至对某些体系比如金属完全无效。不过，因为半经验算法中有可调参数，可以拟合实验数据或从头算数据提高准确性，比如限制个键长，修正个键能啥的。 对于验证过的体系，大多能给出比较理想的结果。4. 什么是SCC-DFTB? SCC-DFTB做了哪些简化？SCC-DFTB全称self-consistent charge density functional tight-binding method，电荷自洽的密度泛函-紧束缚方法。这个方法近几年在德国Marcus Elstner课题组的领导下发展的很快，根本上是源于生物和材料领域对于大规模分子动力学研究反应，固体电子结构及其对外界扰动的响应，输运性质的渴求。SCC-DFTB计算速度很快，我的经验是比DFT快一个数量级以上，用SCC-DFTB轻松跑 300多个原子的MD，几个小时跑上几个ps完全无压力（4核openMP并行）。但同样体系用平面波DFT,32核并行一天只能跑1个ps，而且还用了 赝势。那么，SCC-DFTB到底是怎样一种模型呢？那首先要问，到底什么是紧束缚tight-binding？这是个从固体物理中的一种模型。与固体物理中流行的与原子无关的平面波模型不同，紧束缚借鉴原子轨道组成分子轨道的思想，假定固体中的电子是束缚在原子轨道上的，这样可以用原子轨道的线性组合表示固体中的离域的电子轨道。由于组合后的分子轨道叠加了不同原子的轨道，可以认为该轨道上的电子，可以在不同原子间跳来跳去，所以又称hoping model。在紧束缚的思想中，原子轨道是一组近似轨道，而且不会随轨道间作用而变化，极其类似于半经验方法中的轨道冻结。DFTB对某些CHNO体系的结果还是非常好的，因为这些体系原子的轨道在成键作用中的确变化不大。然而，跟从头算的半经验类似，冻结轨道会对某些强相关 体系(轨道作用较强)造成较大误差，也需要校正。不过，这里并非像AM1/PM3那样添加gaussian function来校正，而是类似于DFT，把能量对电荷密度变化进行泰勒级数展开，这样，如果能得到准确的一阶导，二阶导，三阶导的话，就能准确描述能量随电荷密度变化的情况。这里的电荷密度只考虑原子的patial charge，而不是像正常DFT那样对空间各处电荷密度进行计算。即把能量在中性原子附近展开为原子局部电荷的密度变化q的函数（以中心原子为参考态）：E(q) = E(0) + E(0)*q + 1/2*E(0)*(q)2 + 1/6*E(0)*(q)3 + .展开到二级称为SCC-DFTB2, 展开到三级称为SCC-DFTB3. 在SCC-DFTB的实现中，电荷波动展开包括同一原子的多级展开和不同原子间的非对角项展开，比如二级展开：ESCC-DFTB2 = 轨道能 + 排斥能 + 1/2*sum(q*q*) 其中第二项即为双原子(和)电荷波动展开到二级，求和遍及所有原子和原子对。 电荷变化q怎么计算呢？对价层电子波函数按TB理论做SCF自洽计算，收敛后做Mulliken分析得到原子的partial charge电荷，即为q。这就是SCC的含义，self-consistent charge, 即猜测冻结的原子轨道组成分子轨道的系数，得到Mulliken电荷，然后得到能量，如不收敛，则重新给定系数，计算能量，。直到收敛。其中能量E对电荷的多阶导E(0), E(0), E(0)必须预先确定。怎么确定呢？通过拟合DFT计算的势能面获得。这个多阶导反应了原子电荷在原子间作用下波动的倾向。一阶导是电负性，二阶导称化学硬度。如果拟合并不能完全重现DFT的势能面，怎么办呢？把拟合剩余的能量归到一个排斥势函数中，作为参数放到参数 文件中。除了冻结原子轨道，为了加快计算速度，SCC-DFTB也采用半经验中那些简化手段：第一，价层电子的DFT-SCF计算采用最小基组近似。对于分子轨道，它还通过一个二次函数(r/r0)2到能量项中对核外电子云的分布做一个适当的收缩，以取得更好的结果。这里包括一个可调参数r0，决定收缩的范围（经验参数）。此外，自从SCC-DFTB2以来，还存在一个gamma函数，用来系统 修正原子间排斥过小的问题，它的作用是抬升整个势能面；氢键作用也有一个专门的修正，是对gamma函数做进一步变形和参数化得到的。在文献中看到 SCC-DFTB3-h, h就表示氢键修正。如果有g(或)就表示包括gamma函数修正。这里是增加变量以便能通过拟合提高准确性。第二，分子轨道自洽要计算不同原子对间的原子轨道积分，但是因为原子轨道是冻结的，它们的积分可以预先计算然后保存到参数文件中。如果你想知道你的体系能 不能用SCC-DFTB，那你就需要先看看参数文件，有没有你的体系的所有元素的两两参数。在标准的DFTB+网站，可供下载的参数文件十分有限，目前主 要包括一些CHNO等常用元素，目前是制约SCC-DFTB应用的主要瓶颈。5. SCC-DFTB准确性怎么样？首先，作为一种含展开系数的模型，其准确性必然依赖于拟合的过程，比如，拟合训练基的选取，DFT势能面计算的精度，等等。我使