用第一性原理方法获取周期体系中原子的部分电荷

李亚娜 吕 洋 周立川 陈 理 李慎敏

(大连大学,辽宁省生物有机化学重点实验室,辽宁大连 116622)

用第一性原理方法获取周期体系中原子的部分电荷

李亚娜 吕 洋 周立川 陈 理 李慎敏*

(大连大学,辽宁省生物有机化学重点实验室,辽宁大连 116622)

利用量子化学软件包Crystal计算了立方周期性边界条件下液态水体系的静电势(ESP)和静电场(EF).在此基础上,提出了一种由第一性原理方法获取周期体系中原子的部分电荷的快捷方法.该方法把由周期性边界条件引入的平均静电势φmean作为一个拟合参数,通过对第一性原理静电势与Ewald加和法静电势的最小二乘法拟合而实现.值得说明的是,比较静电势与静电场拟合方法,前者的相对拟合误差仅为2%-3%,比后者小一个数量级.考察了四种电荷限制条件下,静电势、静电场拟合的水分子原子部分电荷及偶极矩的分布情况.

原子的部分电荷; 最小二乘拟合; 静电势; 静电场;Ewald加和法

对于化学工作者来说,原子的部分电荷是一个具有重要意义的不可或缺的概念.许多用于描述体系性质的物理、化学概念或物理量,如分子的极性(偶极矩)、反应性(活性)、静电相互作用等都与体系中电子的分布有关,而描述电子分布的简便、有效、经典的方法是原子的部分电荷(以下简称原子电荷).此外,原子电荷也是计算化学方法,如分子模拟、蒙特卡洛模拟技术等所依赖的描述体系电学性质的重要力场参数.然而,原子电荷却又是一个没有明确的量子力学定义,不能被实验直接测量的人为的物理量.这里的“人为”体现在借助量子力学或实验方法获得的原子电荷往往具有一定的随意性和不确定性[1].

随着计算机技术的快速发展,量子力学第一性原理方法已成为获得原子电荷的主要手段.目前常见的方法有:基于分子轨道(Mulliken电荷)[2]、电子密度(Bader电荷)[3]的空间划分方法;基于原子周围静电势(ESP)[4-8]、静电场(EF)[9]的数值拟合方法等. Mulliken电荷的缺点在于其对计算方法、基组的依赖性较大;Bader电荷的基组依赖性虽小,但它只给出Bader体积内的净电荷,一般不用作模拟力场的电荷参数;静电势电荷因其可以较好地描述分子间的相互作用,特别是近年来发展的参数限制的静电势拟合(RESP)电荷的使用[8,10],使得基于第一性原理静电势的拟合电荷已成为目前许多力场的电荷参数,如著名的AMBER力场[11].然而,静电势拟合方法不适合于凝聚相周期体系中原子电荷的获取,原因在于由周期性边界条件所引入的体系的平均静电势是未知的.基于此,Spackman等[9]提出了静电场拟合电荷的方法.由于静电场是静电势梯度的负值,因而该方法避开了平均静电势的求算.不过,静电场的第一性原理计算以及作为矢量的静电场的拟合,其计算工作量要远大于静电势拟合.

对于离子或强极性凝聚相周期体系,体系中原子的电荷明显依赖于分子的几何构象及周围环境的变化,换句话说,体系的运动过程中存在着显著的电荷极化和电荷转移现象,而这是传统的固定电荷力场方法无法正确描述的.完全的第一性原理方法处理工作量巨大,目前还难以实现.对于无限稀溶质-溶剂体系、生物质-水溶液体系等,变通的方法是通过把系统划分为量子-经典两个作用区域,对于所关心的溶质、生物大分子活性区域等进行量子化学处理，而对其他区域进行经典力学处理，即所谓的量子-经典混合方法[12-13].对于不易划分区域的研究体系,如水等纯溶剂体,类似的应用还鲜见报道.

为了考察上述研究体系的极化作用,近年来,一种发展较快的浮动电荷力场,即第三代力场受到了越来越多人的关注[14-15].该方法利用Sanderson的电负性均衡原理(EEM)[16-17]快速计算体系中原子、分子的瞬态电荷,因而可以考察体系中分子间极化作用及电荷转移作用.该方法计算电荷的准确性依赖于所研究体系中原子的价态电负性及价态硬度参数,这些参数通常来自于对第一性原理获取的点电荷的拟合,与拟合训练集组成、大小以及第一性原理所获取电荷的准确性密切相关.

水溶剂作为常见的极性溶剂一直备受人们的关注.目前用于分子模拟的水分子模型多为固定电荷模型,如SPC水,TIPnP(n=3-5)系列水等[18-20].最近Yang等利用原子-键电负性均衡原理(ABEEM)方法构建了七点的浮动电荷水模型[21],并将其应用于水团簇及生物大分子水溶液的模拟中[22-24].Yang的七点水模型在处理水的簇合物中得到了与第一性原理相媲美的计算结果[21,25].特别地,当将键电荷回归到原子位点时,原子电荷与量子化学计算的Mulliken电荷具有相当好的线性关系.

本文中,我们将周期性边界条件下的液态水体系作为研究对象,利用第一性原理方法计算体系的网格静电势和静电场;然后,利用点电荷的直接截断加和法(Rcut)、Ewald加和法静电势[26]以及静电场模型，通过电荷限制条件下的最小二乘拟合法拟合第一性原理静电势、静电场的计算数据,进而获取周期水体系的原子部分电荷.我们的目的有二:一是通过第一性原理方法获取的静电势电荷与静电场电荷的比较,寻找适合于凝聚相体系原子电荷获取的快速方法;二是在分子动力学模拟环境(周期性边界条件)下,考察水体系中电荷极化和电荷转移作用,为下一步在EEM理论框架内进行价态电负性及价态硬度的参数化提供可靠的训练集数据.

1 理论方法与计算细节

1.1 第一性原理计算静电势与静电场

目前适用于周期体系第一性原理计算的商用软件包并不多,而直接包含静电场计算的则更少, Dovesi等[27]发展的Crystal软件包是其中的一个.本文中,使用Crystal 06软件包在HF/6-21G基组下,计算了立方体周期性边界条件下液态水体系的网格静电势和静电场.具体步骤如下:首先,利用常规的MD软件包GROMACS[28],选取固定电荷的SPC水模型,在NVT系综下(T=298 K)分别对两个尺寸的液态水盒子进行了500 ps的动力学模拟.模拟中,水的密度 ρ=1.0 g·cm-3,水盒子边长 a为 1.0,1.1 nm,与其对应的水分子数分别为34,45.然后,针对模拟平衡后的两个体系的坐标构象进行数据采集,通过随机取样方法分别选取10个水盒子构象为样本数据.最后,利用Crystal程序计算每个样本构象的单点能,水分子的Mulliken电荷,以及每个原子周围给定网格点的静电势和静电场.其中,第一性原理静电势和静电场的采集坐标(网格点)来自于水盒子(元胞)内每个原子范德华半径以外的均匀分布三维网格点,网格点间距为0.05 nm.对于边长为1.0,1.1 nm的两个体系,元胞内网格点数 Np分别是3388和3483.图1给出了a=1.0 nm体系的计算网格点示意图.

1.2 周期体系的点电荷静电势和静电场

对于孤立分子或分子簇,无穷远处的静电势为零.然而,对于无限重复的周期体系(如凝聚相体系的常用处理方法),由于电子密度分布的周期性,无穷远处的静电势不为零.不过,当元胞内体系中总的净电荷为零时,元胞内静电势φ(r)通常采用如下边界条件[9],

式(1)也是Crystal程序中元胞边界条件的选择方式[29],这与Ewald加和法处理长程静电相互作用的边界条件是一致的[30].为满足式(1)的边界条件,周期体系的静电势可以写成如下两部分:

上式中φmean是一个与位置无关的常数项,φ′(r)是由点电荷模型定义的周期体系中元胞内位置r处的静电势:

其中,ri是元胞内原子i的位置矢量,qi是原子i的部分电荷,l是水盒子的尺寸,N是元胞中原子的个数,k表示重复单元的序号.当元胞内净电荷为零时,式(3)是收敛的,我们称其为静电势的截断加和法,但其收敛速度较慢.本文中,以元胞为中心,沿坐标轴三个方向分别取k为-10到10的21×21×21个重复单元.静电势的Ewald加和法可以很好地克服式(3)收敛慢的缺点,近年来已被广泛地应用于长程静电相互作用能的计算.根据静电势的Ewald加和法,元胞内静电势可改写为:

上式中,V是元胞体积,α是高斯函数的宽度系数,用以描述原子位点电荷密度的分布,k是r在倒易空间中的对应量,erfc是误差函数的补函数.值得注意的是,与静电势的Ewald加和法不同,这里的φ′(r)ES只包含Fourier变换项和短程校正项两项,不包含静电能中的自能项部分.

需要说明的是,由于式(2)中φmean的存在,给利用静电势法精确地拟合体系的原子电荷带来了一定的困难.最近,Spackman等[9]给出了静电场拟合电荷的方案,其中,E(r)是静电势φ(r)梯度的负值,如(5)式:

也就是说,E(r)与φmean的选取无关,从而避开了φmean的求算,减少了拟合的不确定性.不过,由于E(r)是矢量,相应的第一性原理计算与拟合工作量较静电势的计算与拟合要大许多.

1.3 电荷限制的最小二乘法拟合

通过Lapack软件包中DGGLSE程序,利用周期体系点电荷静电势、静电场模型,针对Crystal 06计算给出的水盒子体系10个构象的网格静电势、静电场数据,分别进行了四个电荷限制条件下的最小二乘法拟合.拟合评价判据由下式定义的相对拟合误差的百分比给出,

上式中,A0(ri)是第一性原理计算的网格点ri处的静电势或静电场,Np是元胞内网格点数.χ2由式(7)给出:

其中,A(ri)是式(2)-(5)定义的ri处静电势或静电场. Q是元胞内中的净电荷(本文中Q=0),qi是原子i的部分电荷,λ是拉格朗日乘因子.需要说明的是,对于周期体系,由于φmean是未知的,拟合中为获得静电势拟合电荷,我们设其为一新的拟合参数qN+1,并构造新的拟合判据:

实际计算中我们还通过引入多个拉格朗日乘因子的方式尝试了其它电荷限制条件下的拟合情况.例如,若只允许分子电荷的极化而不允许分子间电荷转移,即限制每个水分子的净电荷为零.此外,我们还仿照Kollman引入限制函数χ200rstr,考察了Mulliken电荷限制下的电荷拟合(RESP)[10],这时,

式(9)中χ20就是式(8)中的χ2;式(10)中,qMui00为Crystal 06计算的原子i的Mulliken电荷.

2 结果与讨论

2.1 不同方法拟合误差的比较

利用Crystal量子化学软件包,首先计算了立方体周期性边界条件下液态水体系的第一性原理网格静电势和静电场.然后,利用式(2)-(5),通过电荷限制的最小二乘拟合获得了元胞中每个水分子的原子部分电荷.依据拟合数据来源划分,本文的工作可归纳为静电势拟合和静电场拟合两种;而依据对周期体系点电荷静电势计算模型的划分,静电势拟合又可归纳为简单的截断加和法(Rcut)和Ewald加和法两种.此外,为了合理地描述体系中电荷极化与转移作用,针对以上三种拟合方法,分别考察了四种电荷限制条件下的拟合情况:(1)元胞中总的净电荷为零;(2)满足条件(1),同时加入Mulliken电荷限制; (3)满足条件(1),同时加入分子的净电荷为零限制; (4)条件(1)、(2)和(3)同时满足的情况.不同方法的拟合判据由10个水盒子构象的平均拟合误差及均方差给出(见表1).较大的体系,即45个水分子元胞中氧、氢原子电荷的拟合结果见表2,拟合电荷的分布情况分别见图2和图3.

由表1我们看到,在静电势拟合法中Ewald加和法的相对误差及标准偏差明显低于截断加和法.在较小的34个水分子体系中,Ewald加和法最大误差出现在限制条件4的情况,只有3.75%,比截断加和法最小误差5.05%(限制条件1)还小.另一方面,静电场拟合法虽然不涉及平均静电势φmean,但由静电场拟合电荷给出的网格静电场与第一性原理方法计算得到的网格数值相比,相对误差达到了30%-40%,远大于静电势拟合相对误差.我们认为在相同的网格点处,静电势与静电场数值值域分布的差异是导致上述情况的主要原因.在当前的水体系中,由于φmean的存在,静电势值域在-0.25--0.05 a.u.,其峰值在φmean附近;而作为矢量的静电场，其值域落在以原点对称的-0.10-0.10 a.u.区间内,其峰值在0.0 a.u.处.在较多的网格点上静电场数值趋近于零的事实使得式(6)中分母项数值较小,表现为静电势拟合出现较大的相对误差.事实上,Spackman等[9]在利用静电场法拟合尿素晶体时也发现该方法存在较大的相对拟合误差.值得一提的是,通过截断加和法静电势拟合的电荷,再利用式(5)而得到的网格静电场数值与第一性原理计算的静电场比较,相对误差仅比静电场拟合法略大.这暗示着即使是对静电场计算有较高要求的模拟体系,考虑到静电势拟合计算效率较静电场拟合高很多,通过静电势拟合获取电荷的方法也是值得考虑的.这一点也可以从静电场拟合电荷通过式(3)所计算的静电势与第一性原理比较具有较大的偏差而进一步被确认.

表1 两个水体系四种电荷限制条件下拟合评价Table 1 Fitting evaluation of the four charge restrained fits on the two water systems

此外,无论是静电势拟合还是静电场拟合,限制条件的加入会使拟合误差相应地增加.综合考虑各种拟合情况容易发现,加入限制条件3,即不考虑分子间电荷转移的情况,拟合误差变化较小,而加入Mulliken电荷限制(RESP限制)后,相对的拟合误差变化较大.

表1还给出了由静电势拟合得到的45个水分子的平均静电势.可以看出,不同限制条件下, Ewald加和法得到的样本均方根偏差均很小,说明平均静电势与体系的构象相关性不大,这与我们预期的结果相一致.

需要说明的是,利用Ewald加和法拟合电荷,不仅静电势的相对误差和标准偏差小、计算效率高、收敛快;而且,比较两个水分子体系的计算结果可知,拟合误差与体系尺寸的依赖性几乎可以忽略不计.可以说,对于周期性体系,利用Ewald加和法通过拟合第一性原理静电势获取力场中原子部分电荷的方法是首选方法.

最近,Woo等[31]通过引入改进的误差函数方法,提出了一种利用静电势计算周期体系原子电荷的方法,并成功应用于金属有机物多孔材料体系中.

2.2 不同方法电荷分布的比较

接下来,我们讨论不同方法下,拟合得到的水分子中氧原子、氢原子的部分电荷及其分布情况.对于我们考察的两个水盒子体系,由于体系尺寸对拟合的结果影响不大,为节省篇幅,这里我们只给出较大的45个水分子体系的计算结果(见表2及图2和图3).容易发现,无论哪种拟合方法,在引入Mulliken电荷限制后(限制条件(2)与(4)),氧、氢拟合电荷的平均值,相对于限制条件(1)与(3),其绝对值均呈减小的趋势,同时,电荷分布范围也明显变窄.而通过限制条件(3)拟合的氧、氢电荷,不考虑电荷转移,具有最大的涨落,以截断法静电势拟合为例,氧的最大原子电荷为0.327e,而最小值为-1.423e.由图2和图3可见,在限制条件(1)与(3)下,静电势截断法拟合的电荷呈现最大的分布范围,而Ewald加和法与静电场拟合的电荷呈现相似强度的分布;与限制条件(2)与(4)比较,无论是静电势拟合还是静电场拟合,氧原子的电荷具有较大的负电荷,而氢原子具有较大的正电荷.值得强调的是引入 Mulliken限制的Ewald加和法拟合的水分子电荷与Crystal程序计算的Mulliken电荷分布具有良好的一致性.

表2 45个水分子体系中不同拟合方法获得的原子电荷和分子电荷涨落Table 2 Atomic charges by different fitting methods and its charge fluctuation for each molecule of the 45 water systems

下面,我们讨论一下分子的电荷转移问题.在限制条件(1)、(2)下,利用对Crystal程序计算的静电势或静电场拟合所得到的水分子原子电荷,计算了体系中水分子的平均净电荷涨落,ΔQ(见表2)用以衡量水分子中电荷转移能力.我们发现,静电势截断加和法拟合的电荷具有最大的电荷转移,单分子中最大的电荷涨落为0.410e;Ewald加和法拟合的电荷转移最小,单分子中最大的电荷涨落仅为0.084e.

2.3 不同方法水分子偶极矩及其分布的比较

由于分子偶极矩是一个实验可观测量,因此,计算分子偶极矩及其分布已成为评价凝聚相力场好坏的重要一环.本文中,由于所选体系的Crystal几何全优化工作量巨大,难以实施,但考虑到我们所研究的水体系中原子的部分电荷可能是决定偶极矩大小关键因素,因此,对于不同拟合方法得到的原子电荷,分别计算了45个水分子体系的平均偶极矩并给出了动力学模拟中10帧水构象450个水分子的偶极矩分布,计算结果分别见表3和图4.

表3列出了元胞中45个水分子10个构象的平均偶极矩、均方根偏差及样本中单个水分子偶极矩的最大、最小值.可以看到,在Mulliken电荷限制条件(2)和(4)下,静电势和静电场拟合电荷均给出较小的分子偶极矩.偶极矩减小的这种变化趋势,使得其更趋向于固定电荷的SPC刚性水的偶极矩值,7.84× 10-30C·m[31],以及第一性原理Mulliken电荷的偶极矩值,7.71×10-30C·m,而与体相水的实验观测值, 8.34×10-30C·m[32],有一定的偏差.我们认为导致计算与实验观测差异的主要原因是本文中水的几何结构为SPC刚性水构型,同时在第一性原理计算时只进行了单点能计算,没有进行几何结构的优化.

需要说明的是由10帧构象450个水的偶极矩的标准偏差来看,与上节讨论的氧、氢原子的电荷分布一致,不考虑电荷转移体系的(限制条件(4))偶极矩,其涨落最小,这一点可以从图3偶极矩分布图中清晰看到.有趣的是,在Mulliken电荷限制条件下,尽管静电场拟合电荷与Ewald加和法静电势拟合电荷存在一定的差异,但偶极矩的分布,包括峰宽和峰强,却呈现出很好的一致性.由于偶极矩作用决定了许多体系的静电性质,这表明两种方法拟合的电荷所产生的静电性质的相似性.

表3 不同拟合方法获得电荷计算的45个水分子体系中水分子的平均偶极矩Table 3 Average dipole moment(μ)of the 45 water system calculated by the derived charges from different fitting methods

3 结 论

针对凝聚相体系,提出了一种利用第一性原理静电势获取原子部分电荷的快捷方法.该方法把周期性边界条件引入的平均静电势φmean作为一个拟合参数,通过对第一性原理静电势与截断加和法或Ewald加和法点电荷静电势进行数值拟合,进而获得周期体系力场中原子的部分电荷.该方法相对拟合误差与相应的静电场拟合电荷所给出的相对误差要小得多.特别地,Ewald加和法的静电势拟合相对误差仅为2%-3%,比静电场相对误差小一个数量级.此外,与我们期望的一致,Ewald加和法拟合的φmean随体系几何构象、尺寸的选取变化不大.考虑到分子间的静电相互作用主要与静电势有关,并且静电势拟合的计算效率要远高于静电场拟合,我们认为Ewald加和法是一种通过第一性原理计算获取周期性边界条件下原子电荷的有效、快捷方法.

在不明显增加拟合相对误差的前提下,我们讨论了四种电荷限制条件下的静电势和静电场电荷拟合情况.结果表明,无论是静电势还是静电场拟合,引入Mulliken电荷限制的拟合电荷,其平均值较小,分布较窄.另一方面,在只限制元胞净电荷为零时(限制条件(1)),分子间不仅存在电荷极化还存在电荷转移.其中,电荷转移的大小与拟合的方法存在较大的依赖性.静电势截断加和法拟合的电荷,分子平均的电荷转移最大,为0.410e,而Ewald加和法电荷转移最小,仅为0.084e.为了避免电荷转移被人为地夸大,我们还考察了加入分子电中性的两个电荷限制条件的拟合(限制条件(3)和(4)).需要指出的是,在分子电中性限制下,分子中原子部分电荷的分布范围明显增大,显然,这种拟合电荷呈现了较强的电荷极化作用.

利用获得的水体系中原子的部分电荷,我们计算了水分子的平均偶极矩,与第一性原理计算的偶极矩比较,引入Mulliken电荷限制的体系呈现了更好的趋势.此外,利用第一性原理快速获得体系的原子部分电荷,也将为EEM方法中对于原子价态电负性、价态硬度参数的拟合提供可靠的训练集.

1 Verstraelen,T.;Speybroeck,V.V.;Waroquier,M.J.Chem.Phys., 2009,131:044127

2 Mulliken,R.S.J.Chem.Phys.,1955,23:1833

3 Bader,R.F.W.;Matta,C.F.J.Phys.Chem.A,2004,108:8385

4 Breneman,C.M.;Wiberg,K.B.J.Comput.Chem.,1990,11:361

5 Arroyo,S.T.;Martin,J.A.S.;Carcia,A.H.Chem.Phys.Lett., 2002,357:279

6 Besler,B.H.;Merz,K.M.;Kollman Jr.,P.A.J.Comput.Chem., 1990,11:431

7 Singh,U.C.;Kollman,P.A.J.Comput.Chem.,1984,5:129

8 Wang,J.;Wolf,R.M.;Caklwell,J.W.;Kollman,P.A.;Case,D. A.J.Comput.Chem.,2004,25:1157

9 Whitten,A.E.;Mckinnon,J.J.;Spackman,M.A.J.Comput. Chem.,2006,27:1063

10 Wang,J.;Cieplak,P.;Kollman,P.A.J.Comput.Chem.,2000,21: 1049

11 Case,D.A.;Pearlman,D.A.;Caldwell,J.W.;Cheatham III,T.E.; Wang,J.;Ross,W.S.;Simmerling,C.;Darden,T.;Merz,K.M.; Stanton,R.V.;Cheng,A.;Vincent,J.J.;Crowley,M.;Tsui,V.; Gohlke,H.;Radmer,R.;Duan,Y.;Pitera,J.;Massova,I.;Seibel,G. L.;Singh,U.C.;Weiner,P.;Kollman,P.A.AMBER 7 user′s Manual.California:University of California,2002

12 Gao,J.;Xia,X.Science,1992,258:631

13 Field,M.J.;Bash,P.A.;Karplus,M.J.Comput.Chem.,1990,11: 700

14 Patel,S.;Brooks,C.L.J.Comput.Chem.,2004,25:1

15 Varekova,R.S.;Koca,J.J.Comput.Chem.,2006,27:396

16 Sanderson,R.T.Chemical bond and bond energies.New York: Academic Press,1976

17 Sanderson,R.T.Polar covalence.New York:Academic Press, 1983

18 Berendsen,H.J.C.;Grigera,J.R.;Straatsma,T.P.J.Phys.Chem., 1987,91:6269

19 Jorgensen,W.L.;Madura,J.D.J.Am.Chem.Soc.,1983,105: 1407

20 Mahoney,M.W.;Jorgensen,W.L.J.Chem.Phys.,2000,112: 8910

21 Yang,Z.Z.;Wu,Y.;Zhao,D.X.J.Chem.Phys.,2004,120:2541

22 Zhang,Q.;Yang,Z.Z.Chem.Phys.Lett.,2005,403:242

23 Yang,Z.Z.;Zhang,Q.J.Comput.Chem.,2006,27:1

24 Yang,Z.Z.;Qian,P.J.Chem.Phys.,2006,125:064311

25 Wu,Y.;Yang,Z.Z.J.Phys.Chem.A,2004,108:7563

26 Ewald,P.P.Ann.Phys.,1921,64:253

27 Dovesi,R.;Saunders,V.R.;Roetti,C.;Orlando,R.;Zicovich-Wilson,C.M.;Pascale,F.;Civalleri,B.;Doll,K.;Harrison,N.M.; Bush,I.J.;D′Arco,P.;Llunell,M.Crystal 06 user′s manual. Torino:University of Torino,2006

28 Spoel,D.v.d.;Lindahl,E.;Hess,B.;Groenhof,G.;Mark,A.E.; Berendsen,H.J.C.J.Comput.Chem.,2005,26:1701

29 Saunders,V.R.;Freyria-Fava,C.;Dovesi,R.;Salasco,L.;Roetti, C.Mol.Phys.,1992,77:629

30 Stewart,R.F.God Jugosl Cent Kristalogr,1982,17:1

31 Campañá,C.;Mussard,B.;Woo,T.K.J.Chem.Theory Comput., 2009,5:2866

32 Shirono,K.;Daiguji,H.Chem.Phys.Lett.,2006,417:251

Atomic Partial Charges for Periodic Systems from First-Principles Calculations

LI Ya-Na LÜ Yang ZHOU Li-Chuan CHEN Li LI Shen-Min*
(Liaoning Key Laboratory of Bio-Organic Chemistry,Dalian University,Dalian 116622,Liaoning Province,P.R.China)

We calculated the electrostatic potential(ESP)and electric field(EF)of periodic liquid water systems using the quantum chemistry software package,Crystal.We propose a method to obtain atomic partial charges rapidly for periodic systems based on first-principles calculations.In this method,the average electrostatic potential φmean, which is introduced to meet the periodic boundary condition,is taken as a parameter during the least squares fitting of the ESP from first-principles calculations and used in the Ewald summation.A comparison of the two methods,i.e., ESP and EF fitting,reveals that the relative root mean-square deviation(RMS)of the former method is only 2%-3%, which is one order of magnitude smaller than that of the latter method.In addition,the distribution of the derived atomic partial charges and dipole moments for the water system are discussed using four charge restrained fits.

Atomic partial charge; Least square fitting; Electrostatic potential; Electrostatic field; Ewald summation

O641

Received:April 24,2010;Revised:June 21,2010;Published on Web:August 27,2010.

*Corresponding author.Email:shenmin@dl.cn;Tel:+86-411-87402384.

The project was supported by the National Natural Science Foundation of China(20973027,20633050)and Program for Liaoning Excellent Talents in University,China(2007R02).

国家自然科学基金(20973027,20633050)和辽宁省高等学校优秀人才支持计划项目(2007R02)资助

ⒸEditorial office of Acta Physico-Chimica Sinica