|
gromacs教程-7-原子电荷什么是原子电荷？ 原子电荷,即位于原子中心的点电荷,是对化学体系中电荷分布最简单、最直观的描述方式之一.它有很多重要意义,比如帮助化学工作者研究原子在各种化学环境中的状态、考察分子性质、预测反应位点.另外原子电荷模型在计算化学中也有很多实用价值,如作为分子描述符用于药物虚拟筛选、在分子对接和分子动力学/蒙特卡罗模拟中描述静电作用、在量子力学(QM)与分子力学(MM)结合计算中表现QM区域原子对MM区域原子的静电作用势。在其它一些理论方法中也需要借助原子电荷,如Pipek-Mezey轨道定域化方法、电荷自洽的紧束缚密度泛函方法(SCC-DFTB)、Truhlar课题组开发的一系列溶剂模型(SM)等. 原子电荷并不是可观测量,也没有客观、唯一的定义,因此获得原子电荷的方法多种多样.通过一些实验数据可以间接地考察原子电荷,如通过分子的多极矩、红外光谱强度和频率、配体场分裂能、核磁共振(NMR)位移等,但是数据的获取相对不便,与原子电荷的对应关系存在较大经验性,也不能分析不稳定的体系和状态.计算化学的发展使原子电荷方便、快速、可靠的获得成为了可能.自从1955年Mulliken电荷被提出以来,迄今已有不下60种原子电荷计算方法被相继提出,近年仍有许多研究者在改进计算方法. 虽然在许多计算方法的原文和部分文献中都对不同类型原子电荷进行了比较和讨论,但是涉及的方法并不全面,测试分子和测试内容缺乏统一性.因此,很有必要对较为重要的、文献中常涉及的原子电荷计算方法进行综合的比较,使它们的特点、优缺点能够清晰地展现,这将有助于在实际问题研究中选择适合的计算方法.我们选取resp/Gasteiger和电荷均衡(QEq)共3种方法.在文中将首先介绍这些方法的基本原理和特点,之后通过实际体系测试和比较它们的各方面性能,最后讨论计算量的大小,并给出方法选择上的建议. 各种原子电荷计算方法的原理和特点？ 卢天, 陈飞武. 原子电荷计算方法的对比[J]. 物理化学学报, 2012, 28(1):1-18. doi: 10.3866/PKU.WHXB2012281 RESP，全称Restrained ElectroStatic Potential，是由Kollman等人在1993年提出的计算原子电荷的方法。在Amber分子动力学模拟中，RESP是有机小分子计算电荷的标配方法。它需要配合使用量化计算软件（如Gaussian、ORCA）获得静电势，然后进行限制性拟合计算得到。 传统拟合静电势方法存在一些弊端,例如:(1)电荷对构象依赖性较大;(2)单一构象拟合的原子电荷不能体现原子的等价性,例如一氯丙烷中甲基的三个氢是化学等价的,但任何构象下拟合出来的这三个氢的电荷都不是全同的;(3)内部被包埋的原子电荷拟合不准确,这是由于这些原子距离拟合点过远所导致的.这三个问题使得拟合静电势电荷用于分子动力学模拟等问题中可能导致错误的结果.限制性静电势拟合(RESP)方法可以很大程度解决这些弊端,它对内部原子在拟合时施加限制势以降低数值不稳定性,对等价原子在拟合时强制相等,并且将拟合分为两步或者多步以提高拟合质量.然而计算RESP电荷过程相对复杂,而且需要较为准确的量子化学方法计算静电势,用于大批量或者较大分子体系较为耗时. 与之比较接近的电荷是2000-2002年Bayly等人提出的AM1-BCC电荷，在Amber中自带计算程序。它首先在AM1水平下优化结构和计算静电势拟合得到电荷，然后进行键电荷校正（BCC）。AM1-BCC的计算开销比RESP要小，准确性不比RESP更好，但适合于处理大批量的分子，或者没有Gaussian等商业软件版权的情况下使用。 AM1-BCC34是一种简单的近似计算RESP电荷的方法,其电荷是在AM1电荷(即AM1半经验方法获得的波函数产生的Mulliken电荷)基础上进行键电荷校正(BCC)得到的.校正过程只依赖于原子类型和原子间连接关系.校正参数拟合自2755个种类多样的有机分子的HF/6-31G*下的RESP电荷.例如甲醇中的碳原子形成三个“Csp3-single-H”键和一个“Csp3-single-Osp3”键,根据查寻事先拟合的BCC参数就可得知校正电荷应为3×0.0274+1×0.0835= 0.1657(本文电荷单位均为原子单位),加到碳的AM1电荷上就是AM1-BCC电荷. Gasteiger电荷也被称为轨道电负性部分均衡(PEOE)电荷,这种方法利用了电负性均衡概念,即电负性不同的原子成键时,电负性较小的原子附近电子密度会流向电负性较大的原子.在这个过程中原先电负性小的原子电负性会增大.当所有原子间电负性相等时,电子密度的分布就是平衡状态分布. PEOE将原子电负性与原子电荷关系表达为
其中a、b和c为预设参数,对于相同元素但价层轨道处于不同杂化状态的情况参数不同.参数通过实验测定的相应元素的基态和电离态在相应杂化状态下的电离势与电子亲合势来计算.PEOE的迭代过程如同电荷不断在原子间转移的过程,当迭代收敛时电荷分布也就平衡了.在每一步中按照下式计算每一对相连原子间电荷转移量
电荷均衡方法(QEq)得到的是在电负性完全均衡下的原子电荷.电负性表达式为电负性均衡条件要求χ1=χ2=…=χN,它提供了N-1个条件,再将所有原子电荷总和等于分子净电荷作为剩余条件,通过解线性方程即可到得全部原子电荷.式(19)中是广义化的Mulliken电负性为自库仑积分,是原子电离势与电子亲和势的差值.JAB描述了原子间静电作用,若两原子距离较远,则JAB= 1/rAB;若原子间距离较近,电子云静电屏蔽作用不能忽视,此时JAB被定义为两原子间库仑积分.每个原子用单个Slater函数描述电子密度分布,其中指数通过拟合碱金属卤化物实验偶极矩数据确定,对所有原子都相同.普适型力场UFF66使用的正是QEq电荷模型. 值得一提的是QEq方法可以视为另一种使用广泛的方法电负性均衡方法(EEM)的变体.
在分子力学和分子动力学模拟中,计算分子中各原子上的电荷分布是一个十分重要的步骤,特别是在研究极性分子时,电荷分布的准确程度对于模拟结果有很大的影响.计算电荷分布有很多种方法,如Gasteiger方法,Qeq方法,量子力学半经验算法和从头计算法等,也有部分力场自带了电荷分布信息,如OPLS力场.这些方法的计算速度和所得结果的准确程度有较大的差别,Gasteiger方法的速度最快,但准确程度最低,从头计算法最为精确,但耗时也最长.Qeq方法和半经验算法的准确程度介于前两者之间,但在计算速度上,Qeq方法快于半经验算法.从综合的角度来看,在对拥有数百个至数千个原子的体系进行分子力学或分子动力学模拟时,Qeq方法是一种计算耗时和准确性较为平衡的算法,它的计算时间复杂度为T(n)=O(n~2).采用Qeq方法对含有450个原子的聚乙烯链进行了电荷分布的计算,将得到的结果与用Gaussian03(DFT,6-31G)计算得到的结果进行了对比.通过对比可以发现,两者计算得到的电荷分布整体趋势相同,并且Qeq方法得到的单个原子电荷的绝对值约为从头计算法所得结果的4倍,准确程度不及从头计算法,但比Gasteiger法精确得多.在计算时间方面,Qeq方法的计算速度约为从头计算法的30至50倍左右.因此,对于需要实时计算电荷分布的分子动力学模拟,Qeq方法是更为有效的方法（上面就MS会这么弄），建议可以考虑使用RESP和AM1-BCC方法。 |
