可快速计算水团簇三体作用强度的新方法

引用本文

李晓蕾, 王长生. 可快速计算水团簇三体作用强度的新方法[J]. 高等学校化学学报, 2019,40(5): 950-958.
LI Xiaolei, WANG Changsheng. Rapid Calculation of the Three-body Interaction Energies in Water Clusters^† [J]. Chemical Journal of Chinese Universities, 2019,40(5): 950-958.
Doi:10.7503/cjcu20190038 复制到剪切板

Permissions

可快速计算水团簇三体作用强度的新方法

李晓蕾, 王长生

辽宁师范大学化学化工学院, 大连 116029

联系人简介: 王长生, 男, 博士, 教授, 博士生导师, 主要从事理论与计算化学研究. E-mail: chwangcs@lnnu.edu.cn

收稿日期: 2019-01-14

基金资助: 国家自然科学基金(批准号: 21773102)资助.

摘要

将水分子视为由2个O—H键偶极构成, 再将水分子间的三体作用视为长程诱导作用和短程校正之和, 使用Thole模型计算长程诱导作用, 通过同时考虑不同水分子间的置换和同一个水分子中2个键偶极间的置换计算短程校正, 从而提出了一个可快速计算水团簇三体作用强度的新方法. 根据已报道的12347个水三聚体的结构和CCSD(T)三体作用能, 确定了该方法所需参数. 将该方法和所确定的参数应用于67个水团簇体系, 计算这些体系的三体作用能, 并与CCSD(T), MP2, M06-2X方法的计算结果进行比较. 结果表明, 相对于CCSD(T)方法的总三体作用能, 本文方法的均方根偏差(RMSD)仅为3.32 kJ/mol, 平均相对偏差(MRD)仅为2.43%; 对较大水团簇体系, 该方法计算精度稍优于MP2方法, 明显优于M06-2X方法, 并且更快捷高效.

关键词: 水团簇; 三体作用; 键偶极

中图分类号:O641 文献标志码:A

Rapid Calculation of the Three-body Interaction Energies in Water Clusters^†

LI Xiaolei, WANG Changsheng^*

School of Chemistry and Chemical Engineering, Liaoning Normal University, Dalian 116029, China

Fund:† Supported by the National Natural Science Foundation of China(No.21773102).

Abstract

A method was proposed to calculate the three-body interaction energies in water clusters. The O—H bond of a water molecule was regarded as a bond-dipole, and the three-body interaction energy were expressed as the sum of long-range induction and short-range correction. The long-range induction was treated by using Thole model, and the short-range correction was derived by considering both permutations among different water molecules and permutations between the two dipoles within each molecule. The parameters needed in the method were determined by fitting to the CCSD(T)/aug-cc-pVTZ three-body interaction energies of 12347 water trimers taken from the reference. The accuracy of the method is validated through comparing the three-body interaction energies of 67 water clusters ranging from trimer to 30-mer obtained by CCSD(T) method, MP2 method, and M06-2X method. The method can reproduce the CCSD(T)/aug-cc-pVDZ total three-body interaction energies, with the RMSD of 3.32 kJ/mol and the MRD of 2.43%, slightly better than the MP2 method, and outperforms the M06-2X method for large water clusters. The method is helpful in simulating the dynamic properties of water and other biomolecules.

Keyword: Water cluster; Three-body interaction; Bond-dipole

文章图片

蛋白质核酸等生物分子体系中以及较大水团簇体系中存在着大量的水分子, 这些水分子间存在着大量的多体作用. 研究发现, 二体作用对水团簇稳定性的贡献最大, 三体作用对水团簇稳定性的贡献也达20%以上^[1~4], 因此仅仅准确描述分子间的二体作用不足以保证对生物分子体系稳定性的准确预测, 准确计算水分子间存在的诸多三体作用强度十分重要. Paesani等^[5]发现, 对于水团簇体系, 使用CCSD(T)/aug-cc-pVTZ方法结合三体基函数(Trimer basis)计算得到的三体作用能与CCSD(T)/CBS的计算结果相近. Bettens等^[6]发现, 使用MP2/aug-cc-pVDZ方法结合三体基函数可以得到与MP2/CBS方法精度相当的三体作用能. 但是当水团簇较大时, 团簇所含的三体作用数目急剧增加(如1个由40个水分子构成的水团簇含有40× 39× 38/6=9880个三体作用), 使用高精度从头算量子化学方法计算水团簇的三体作用能会变得十分昂贵. 因此建立一种可快速准确预测水团簇三体作用能的方法十分必要.

本文提出了一种可快速准确计算水团簇三体作用强度的新方法. 将水团簇中的每个三体作用能表述为长程诱导作用和短程校正之和, 长程诱导作用由Thole模型计算获得, 短程校正项由关于键偶极和水分子的置换不变量构建. 计算结果表明, 该方法不但具有较高的计算精度, 并且非常高效快捷, 有望在生物分子体系的动力学模拟中发挥作用.

1 模型和计算方法

将水分子中的每个O— H键视作一个键偶极, 并将水团簇中的每个三体作用能V_3B(kJ/mol)视为长程诱导作用 $V_{3 B}^{ind}$ (kJ/mol)和短程校正 $V_{3 B}^{short}$ (kJ/mol)之和, 即

$V_{3 B} = V_{3 B}^{ind} + V_{3 B}^{short}$ (1)

1.1 长程诱导作用$V_{3B}^{ind}$的计算

仿照多体展开理论对三体作用能的定义^{[6, 7]}, 由水团簇中第m, n, l个水分子组成的三聚体的长程诱导作用 $V_{3 B}^{ind}$ 的计算如下:

$V_{3 B}^{ind} (m, n, l) = E^{ind} (m, n, l) - [E^{ind} (m, n) + E^{ind} (m, l) + E^{ind} (n, l)]$ (2)

式中: $E^{ind} (m, n, l)$ (kJ/mol)和 $E^{ind} (m, n)$ (kJ/mol)分别为三聚体和二聚体的诱导能, 其数值可以通过不同方案计算得到. 本文使用Thole模型计算二聚体和三体的诱导能^{[8, 9, 10, 11, 12]}:

$E^{ind} = - \frac{1}{2} \sum_{i} μ_{i} E_{i}$ (3)

$\begin{array}{l} μ_{i} = α_{i} (\sum_{j} T^{ij} M_{j} + \sum_{j^{'}} T^{i j^{'}} μ_{j^{'}}) \end{array}$ (4)

式中: E_i(N/C)为位点i处的电场; μ _i(C· m)为诱导点偶极; α (C· m²· V^-1)为原子极化率; T(C^-2· J· m^-2)为相互作用矩阵; M(C· m)为固有多极矩矢量.

1.2 短程校正$V_{3B}^{short}$的计算

Bowman等^[13]利用关于原子的和分子的置换不变量构建了多维势函数, 预测甲烷-水-水复合物的三体作用能. 受此启发, 将水分子看成2个O— H键偶极, 构建关于化学键的置换不变量, 进而计算短程校正 $V_{3 B}^{short}$ :

$V_{3 B}^{short} = s \cdot \overset{453}{\sum_{p = 1}} υ_{p} κ_{p}$ (5)

式中: s为开关函数; κ _p为关于键偶极和水分子的置换不变量; υ _p为参数.

1.2.1 置换不变量κ _p的构建将置换不变量κ _p表述为与属于不同分子的2个键偶极间距r有关的e指数函数, 将e指数函数的指数部分包含n个r的置换不变量称为n阶置换不变量. 对一个水的三聚体, 若考虑到六阶置换, 可得2个二阶置换不变量(κ ₁~κ ₂)、 10个三阶置换不变量(κ ₃~κ ₁₂)、 37个四阶置换不变量(κ ₁₃~κ ₄₉)、 106个五阶置换不变量(κ ₅₀~κ ₁₅₅)及298个六阶置换不变量(κ ₁₅₆~κ ₄₅₃), 共计453个置换不变量.

以水团簇中第m, n, l个水分子组成的三聚体为例, 设m₁和m₂为第m个水分子的2个O— H键偶极, n₁和n₂为第n个水分子的2个O— H键偶极, l₁和l₂为第l个水分子的2个O— H键偶极. 仅考虑分子内键偶极间的置换以及3个水分子间的置换, 可得如下所示的一个二阶置换不变量κ ₁:

$κ_{1} = \sum_{\begin{array}{l} a \in m_{1}, m_{2} \\ b \in n_{1}, n_{2} \\ c \in l_{1}, l_{2} \end{array}} [(ε_{ab} + ε_{ac}) e^{- k_{2} (r_{ab} + r_{ac} - r_{2}^{^{0}})} + (ε_{ab} + ε_{bc}) e^{- k_{2} (r_{ab} + r_{bc} - r_{2}^{^{0}})} + (ε_{ac} + ε_{bc}) e^{- k_{2} (r_{ac} + r_{bc} - r_{2}^{^{0}})}]$ (6)

式中: k₂和 $r_{2}^{0}$ 为参数; r_ab为键偶极a, b中心间的距离; ε _ab为键偶极a, b间的固有偶极-偶极作用能. 固有偶极-偶极作用能ε _ab的计算如下^[1421]:

$ε_{ab} = \frac{μ_{0, a} μ_{0, b}}{r_{ab}^{3}} (2 \cos θ_{ab} \cos θ_{ab}^{'} + \sin θ_{ab} \sin θ_{ab}^{'} \cos ζ_{ab})$ (7)

式中: ${μ_{0,}}_{a}$ 和 ${μ_{0,}}_{b}$ 分别为键偶极a和b的固有键偶极矩, 本文中O— H键的固有键偶极矩取实验值 5.04×10^-30C·m^[22]; θ , θ ^'和ζ 为描述2个键偶极相对位置的结构参数.

式(6)中每项均涉及a, b, c 3个键偶极, 它们分属于3个不同的水分子. 若将式(6)中括号内第一项作为水分子置换的基准, 置换水分子后得第二、三项, 在此基础上对分子内键偶极的置换求和, 可得κ ₁. 若定义M^p为水分子的置换算符, 定义D^p为分子内键偶极的置换算符, 则式(6)可表述为

$κ_{1} = D^{p} \{M^{p} [(ε_{m_{1} n_{1}} + ε_{m_{1} l_{1}}) \cdot e^{- k_{2} (r_{m_{1} n_{1}} + r_{m_{1} l_{1}} - r_{2}^{^{0}})}]\}$ (8)

同理, 三阶置换不变量κ ₃、四阶置换不变量κ ₁₃、五阶置换不变量κ ₅₀和六阶置换不变量κ ₁₅₆可分别表述为

${κ_{3} = D}^{p} \{M^{p} [(ε_{m_{1} n_{1}} + ε_{m_{1} l_{1}} + ε_{n_{1} l_{1}}) \cdot e^{- k_{3} (r_{m_{1} n_{1}} + r_{m_{1} l_{1}} + r_{n_{1} l_{1}} - r_{3}^{^{0}})}]\}$ (9)

${κ_{13} = D}^{p} \{M^{p} [(ε_{m_{1} n_{1}} + ε_{m_{1} n_{2}} + ε_{m_{2} n_{2}} + ε_{n_{1} l_{1}}) e^{- k_{4} (r_{m_{1} n_{1}} + r_{m_{1} n_{2}} + r_{m_{2} n_{2}} + r_{n_{1} l_{1}} - r_{4}^{^{0}})}]\}$ (10)

${κ_{50} = D}^{p} \{M^{p} [(ε_{m_{2} n_{2}} + ε_{n_{1} l_{1}} + ε_{n_{1} l_{2}} + ε_{n_{2} l_{1}} + ε_{n_{2} l_{2}}) \cdot e^{- k_{5} (r_{m_{2} n_{2}} + r_{n_{1} l_{1}} + r_{n_{1} l_{2}} + r_{n_{2} l_{1}} + r_{n_{2} l_{2}} - r_{5}^{^{0}})}]\}$ (11)

$κ_{156} = D^{p} \{M^{p} [(ε_{m_{1} l_{2}} + ε_{m_{2} l_{1}} + ε_{n_{1} l_{1}} + ε_{n_{1} l_{2}} + ε_{n_{2} l_{1}} + ε_{n_{2} l_{2}}) e^{- k_{6} (r_{m_{1} l_{2}} + r_{m_{2} l_{1}} + r_{n_{1} l_{1}} + r_{n_{1} l_{2}} + r_{n_{2} l_{1}} + r_{n_{2} l_{2}} - r_{6}^{^{0}})}]\}$ (12)

式中, k_n(n=2~6)和 $r_{n}^{0}$ (n=2~6)为参数.

1.2.2 开关函数s 三体作用强度V_3B随三聚体中任意一个水分子的远离而趋近于0. 在保证计算精度的前提下, 使用合适的开关函数, 可去除与分子间距离较远的体系的计算, 减少计算量. 受Paesani等^[5]和Li等^[23]工作的启发, 将开关函数s设为,

$s = w_{12} \times w_{13} + w_{12} \times w_{23} + w_{13} \times w_{23} (13)$

$\begin{array}{l} w_{ij} = \{\begin{array}{l} 1 x_{ij} < 0 \\ - 6 x_{ij}^{5} + 15 x_{ij}^{4} - 10 x_{ij}^{3} + 1 0 \leq x_{ij} < 1 \\ 0 x_{ij} \geq 1 \end{array} (14) \\ x_{ij} = (R_{ij} - R_{cut}) / Δ R (15) \end{array}$

式中: R_ij为i, j 2个水分子的O原子间距; R_cut和Δ R为参数; x_ij为2个水分子的权函数(Weight function). 开关函数s随三聚体中任意一个水分子的远离而逐渐变小, 当有2个以上R_ij≥ R_cut+Δ R时, s=0, 则短程作用 $V_{3 B}^{short}$ 为0.

1.3 参数的确定

对于任意水三聚体结构, 使用所提出的方法计算三体作用能V_3B, 需已知开关函数s中的参数R_cut和Δ R, 短程校正项中的453个组合系数υ _p(p=1~453), 参数k_n(n=2~6)和 $r_{n}^{0}$ (n=2~6).

选取文献[5]报道的12347个水三聚体结构和CCSD(T)/aug-cc-pVTZ三体作用能作为标准数据, 改变非线性参数R_cut, Δ R, k_n(n=2~6)和 $r_{n}^{0}$ (n=2~6), 解线性函数获得453个线性参数υ _p, 反复迭代, 直至均方根偏差(RMSD)最小, 从而确定本文方法所需参数. 图S1(见本文支持信息)为本文方法和CCSD(T)方法计算的12347个三体作用能的相关性图. 由图S1可知, 对12347个模型分子, 本文方法计算的三体作用能与CCSD(T)方法的结果相符.

2 应用

由N个水分子构成的团簇中共包含 $C_{N}^{3}$ 个三体作用, 团簇的总三体作用能E_3B为团簇中每个三体作用能V_3B之和, 即

$E_{3 B} = \overset{N}{\sum_{m < n < l}} V_{3 B} (m, n, l) (16)$

为了检测所提出方法的准确性, 使用该方法和所确定的参数计算了67个不同水团簇的总三体作用能以及三体作用能随距离和角度变化的能量曲线, 并与使用CCSD(T)/aug-cc-pVDZ方法结合三体基函数的结果进行了对比. 由表S1(见本文支持信息)可见, 使用CCSD(T)/aug-cc-pVDZ方法结合三体基函数、 CCSD(T)/aug-cc-pVTZ方法结合三体基函数计算的总三体作用能差别很小. 为减少计算量, 选择CCSD(T)/aug-cc-pVDZ方法结合三体基函数[简称CCSD(T)方法]的三体作用能作为标准数据. 为了比较, 还给出了使用MP2/aug-cc-pVDZ方法结合三体基函数(简称MP2方法)以及M06-2X/jul-cc-pVTZ方法(jul-cc-pVTZ表示O原子使用aug-cc-pVTZ基组、 H原子使用cc-pVTZ基组^[24])结合三体基函数(简称M06-2X方法)计算得到的总三体作用能. 所有从头算均使用Gaussian 09程序包^[25]完成.

2.1 预测水团簇的总三体作用能

表1表4给出了不同方法计算得到的67个水团簇的总三体作用能. 其中Δ 为各方法计算的总三体作用能与CCSD(T)结果间的绝对误差, δ 为各方法计算的总三体作用能与CCSD(T)结果间的相对误差.

Table 1 Calculated total three-body interaction energies(E_3B) of 21 small water clusters^a

No.	Cluster	Ref.	MP2/aug-cc-pVDZ			M06-2X/jul-cc-pVTZ			This work
No.	Cluster	${E_{3 B}}^{b}$	E_3B	Δ	δ (%)	${E_{3 B}}^{b}$	Δ	δ (%)	E_3B	Δ	δ (%)
1	3w-1a	-5.82^c	-5.77	0.05	0.86	-5.65	0.17	2.92	-5.77	0.05	0.86
2	3w-1b	4.52^c	4.35	-0.17	3.76	5.36	0.84	18.58	4.85	0.33	7.30
3	3w-1c	-10.13^c	-10.29	-0.16	1.58	-11.34	-1.21	11.94	-10.00	0.13	1.28
4	4w-1	-25.82^d	-25.94	-0.12	0.46	-26.53	-0.71	2.75	-25.27	0.55	2.13
5	4w-2	-24.73^d	-24.85	-0.12	0.49	-25.23	-0.50	2.02	-24.60	0.13	0.53
6	4w-3	-15.40^d	-16.11	-0.71	4.61	-19.04	-3.64	23.64	-14.64	0.76	4.94
7	5w-1	-37.95^d	-37.87	0.08	0.21	-36.90	1.05	2.77	-38.41	-0.46	1.21
8	5w-2	-31.46^d	-32.22	-0.76	2.42	-32.38	-0.92	2.92	-30.84	0.62	1.97
9	5w-3	-26.57^d	-27.82	-1.25	4.70	-29.71	-3.14	11.82	-27.07	-0.50	1.88
10	5w-4	-24.85^d	-25.94	-1.09	4.39	-27.95	-3.10	12.47	-24.27	0.58	2.33
11	5w-5	-24.18^d	-25.31	-1.13	4.67	-26.07	-1.89	7.82	-23.30	0.88	3.64
12	5w-6	-26.07^d	-26.53	-0.46	1.76	-26.32	-0.25	0.96	-24.48	1.59	6.10
13	5w-7	-25.94^d	-26.36	-0.42	1.62	-26.48	-0.54	2.08	-25.06	0.88	3.39
14	6w-prism	-36.40^d	-38.07	-1.67	4.59	-37.99	-1.59	4.37	-37.91	-1.51	4.15
15	6w-cage	-37.53^d	-38.99	-1.46	3.89	-40.29	-2.76	7.35	-38.07	-0.54	1.44
16	6w-book1	-43.43^d	-43.76	-0.33	0.76	-42.97	0.46	1.06	-43.22	0.21	0.48
17	6w-book2	-42.30^d	-42.84	-0.54	1.28	-41.21	1.09	2.58	-41.80	0.50	1.18
18	6w-bag	-43.30^d	-44.10	-0.80	1.85	-41.59	1.71	3.95	-41.59	1.71	3.95
19	6w-chair	-49.29^d	-49.08	0.21	0.43	-48.53	0.76	1.54	-51.13	-1.84	3.73
20	6w-boat1	-47.45^d	-47.36	0.09	0.19	-46.28	1.17	2.47	-49.08	-1.63	3.44
21	6w-boat2	-47.45^d	-47.40	0.05	0.11	-46.32	1.13	2.38	-49.41	-1.96	4.13
	RMSD			0.74			1.68			1.02
	MAD			1.67			3.64			1.96
	MRD(%)				2.12			6.11			2.86

a. E_3B, Δ , RMSD, MAD are in kJ/mol; b. calculated with the trimer basis; c. CCSD(T)/CBS three-body interaction energies taken from ref.[26]; d. CCSD(T)/aug-cc-pVTZ three-body interaction energies taken from ref.[27].

Table 1 Calculated total three-body interaction energies(E_3B) of 21 small water clusters^a

2.1.1 较小水团簇表1给出不同方法计算得到的(H₂O)_n(n=3~6)的总三体作用能. 其中, (H₂O)₃的结构和CCSD(T)/CBS三体作用能取自文献[26], (H₂O)_n(n=4~6)的结构和CCSD(T)/aug-cc-pVTZ总三体作用能取自文献[27].

由表1可知, 对这21个较小水团簇, 本文方法计算的总三体作用能与CCSD(T)结果间的最大绝对误差(MAD)在计算六聚体水团簇6w-boat2时产生, 绝对误差为1.96 kJ/mol, 相对误差为4.13%. 对这21个水团簇, 与CCSD(T)方法的结果对比, 本文方法的RMSD为1.02 kJ/mol, 平均相对误差(MRD)为2.86%; MP2方法的RMSD为0.74 kJ/mol, MRD为2.12%; M06-2X方法的RMSD为1.68 kJ/mol, MRD为6.11%; 以上结果表明, 对于这21个较小水团簇, 本文方法与MP2方法计算精度相当, 优于M06-2X方法.

2.1.2 中等水团簇表2给出了不同方法计算的23个中等尺寸水团簇的总三体作用能. 其中, (H₂O)₁₀的结构取自文献[28], 其余水团簇的结构均取自文献[29].

Table 2 Calculated total three-body interaction energies of 23 medium-sized water clusters^a

No.	Cluster	CCSD(T)	MP2/aug-cc-pVDZ			M06-2X/jul-cc-pVTZ			This work
No.	Cluster	${E_{3 B}}^{b}$	${E_{3 B}}^{c}$	Δ	δ (%)	${E_{3 B}}^{c}$	Δ	δ (%)	E_3B	Δ	δ (%)
1	10w-1	-70.58	-73.64	-3.06	4.34	-61.63	8.95	12.68	-70.63	-0.05	0.07
2	10w-2	-71.13	-74.52	-3.39	4.77	-65.81	5.32	7.48	-68.74	2.39	3.36
3	10w-3	-66.94	-70.84	-3.90	5.83	-61.09	5.85	8.74	-66.40	0.54	0.81
4	10w-4	-69.41	-72.97	-3.56	5.13	-62.76	6.65	9.58	-69.79	-0.38	0.55
5	10w-5	-70.84	-73.30	-2.46	3.47	-61.04	9.80	13.83	-75.02	-4.18	5.90
6	10w-6	-68.37	-73.01	-4.64	6.79	-64.43	3.94	5.76	-71.09	-2.72	3.98
No.	Cluster	CCSD(T)	MP2/aug-cc-pVDZ			M06-2X/jul-cc-pVTZ			This work
No.	Cluster	${E_{3 B}}^{b}$	${E_{3 B}}^{c}$	Δ	δ (%)	${E_{3 B}}^{c}$	Δ	δ (%)	E_3B	Δ	δ (%)
7	10w-7	-70.42	-72.59	-2.17	3.08	-61.67	8.75	12.43	-75.27	-4.85	6.89
8	10w-8	-70.08	-72.89	-2.81	4.01	-61.09	8.99	12.83	-73.47	-3.39	4.84
9	10w-9	-69.45	-71.84	-2.39	3.44	-59.20	10.25	14.76	-72.01	-2.56	3.69
10	10w-10	-66.57	-70.17	-3.60	5.41	-61.21	5.36	8.05	-66.44	0.13	0.20
11	11w-1	-79.54	-82.09	-2.55	3.21	-70.00	9.54	11.99	-82.01	-2.47	3.11
12	11w-2	-82.63	-84.89	-2.26	2.74	-73.43	9.20	11.13	-82.30	0.33	0.40
13	11w-3	-78.87	-81.80	-2.93	3.71	-69.54	9.33	11.83	-79.54	-0.67	0.85
14	11w-4	-82.26	-84.22	-1.96	2.38	-69.41	12.85	15.62	-81.92	0.34	0.41
15	11w-5	-81.76	-83.60	-1.84	2.25	-71.84	9.92	12.13	-82.09	-0.33	0.40
16	12w-1	-88.70	-92.47	-3.77	4.25	-77.78	10.92	12.31	-87.91	0.79	0.89
17	13w-1	-99.87	-103.51	-3.64	3.64	-85.77	14.10	14.12	-99.54	0.33	0.33
18	14w-1	-111.55	-115.48	-3.93	3.52	-95.10	16.45	14.75	-108.95	2.60	2.33
19	15w-1	-123.39	-127.90	-4.51	3.66	-101.96	21.43	17.37	-121.92	1.47	1.19
20	17w-1	-142.05	-148.36	-6.31	4.44	-117.32	24.73	17.41	-138.87	3.18	2.24
21	17w-2	-146.15	-150.92	-4.77	3.26	-120.75	25.40	17.38	-142.67	3.48	2.38
22	21w-1	-182.17	-190.54	-8.37	4.59	-141.46	40.71	22.35	-181.71	0.46	0.25
23	21w-2	-184.81	-191.42	-6.61	3.58	-136.86	47.95	25.95	-179.16	5.65	3.06
	RMSD			4.04			17.91			2.50
	MAD			8.37			47.95			5.65
	MRD(%)				3.98			13.50			2.09

a. E_3B, Δ , RMSD, MAD are in kJ/mol; b. calculated by the CCSD(T)/aug-cc-pVDZ method with the trimer basis; c. calculated with the trimer basis.

Table 2 Calculated total three-body interaction energies of 23 medium-sized water clusters^a

由表2可见, 对于这23个水团簇, 本文方法计算的总三体作用能与CCSD(T)结果间的最大绝对误差在计算水二十一聚体第2个构型(21w-2)时产生, 绝对误差为5.65 kJ/mol, 相对误差为3.06%. 由表2还可知, 与CCSD(T)方法的总三体作用能相比, 按3种方法的RMSD, MAD, MRD大小顺序均为本文方法< MP2方法< M06-2X方法. 可见, 本文方法对中等大小水团簇三体作用能的计算精度稍优于MP2方法, 明显优于M06-2X方法.

2.1.3 纳米尺寸水团簇 Truhlar等^[24]将至少有一对核间距大于1.0 nm的水团簇定义为水纳米粒子(Water nanoparticle). 相比小的水团簇, 水纳米粒子的成键性质与大量水体系更为接近^[24]. 为了考察模型对大量水体系的适用性, 计算了23个水纳米粒子的总三体作用能(见表3). (H₂O)₁₆的结构取自文献[30], (H₂O)₂₀和(H₂O)₃₀的结构取自文献[28], 其余水团簇的结构均取自文献[29]. 对23个水纳米粒子, 与CCSD(T)结果相比, 本文方法的最大绝对误差在计算水二十二聚体[(H₂O)₂₂, 22w-1]时产生, 为10.17 kJ/mol, 但相对误差为6.17%.

Table 3 Calculated total three-body interaction energies of 23 water nanoparticles^a

No.	Cluster	CCSD(T)	MP2/aug-cc-pVDZ			M06-2X/jul-cc-pVTZ			This work
No.	Cluster	${E_{3 B}}^{b}$	${E_{3 B}}^{c}$	Δ	δ (%)	${E_{3 B}}^{c}$	Δ	δ (%)	E_3B	Δ	δ (%)
1	16w-1	-135.65	-141.75	-6.10	4.50	-120.16	15.49	11.42	-135.44	0.21	0.15
2	16w-2	-141.00	-147.70	-6.70	4.75	-125.69	15.31	10.86	-139.58	1.42	1.01
3	16w-3	-150.21	-154.60	-4.39	2.92	-129.45	20.76	13.82	-149.91	0.30	0.20
4	16w-4	-150.21	-154.60	-4.39	2.92	-130.46	19.75	13.15	-150.79	-0.58	0.39
5	16w-5	-152.17	-156.82	-4.65	3.06	-132.17	20.00	13.14	-149.83	2.34	1.54
6	17w-3	-132.59	-137.90	-5.31	4.00	-111.67	20.92	15.78	-128.24	4.35	3.28
7	17w-4	-140.88	-147.07	-6.19	4.39	-116.02	24.86	17.65	-139.49	1.39	0.99
8	18w-1	-147.07	-152.80	-5.73	3.90	-123.30	23.77	16.16	-141.21	5.86	3.98
9	19w-1	-164.89	-172.00	-7.11	4.31	-128.57	36.32	22.03	-155.98	8.91	5.40
10	20w-1	-173.05	-180.54	-7.49	4.33	-140.50	32.55	18.81	-179.24	-6.19	3.58
11	20w-2	-170.75	-178.45	-7.70	4.51	-144.01	26.74	15.66	-170.75	0	0
12	20w-3	-166.23	-173.22	-6.99	4.21	-140.25	25.98	15.63	-173.34	-7.11	4.28
13	20w-4	-167.99	-175.98	-7.99	4.76	-133.93	34.06	20.28	-173.76	-5.77	3.43
14	20w-5	-170.87	-178.49	-7.62	4.46	-143.30	27.57	16.14	-175.60	-4.73	2.77
15	20w-6	-166.23	-175.81	-9.58	5.76	-144.10	22.13	13.31	-172.09	-5.86	3.53
16	20w-7	-166.65	-175.39	-8.74	5.24	-136.02	30.63	18.38	-167.86	-1.21	0.73
17	20w-8	-166.86	-174.22	-7.36	4.41	-138.41	28.45	17.05	-171.25	-4.39	2.63
18	20w-9	-168.49	-176.19	-7.70	4.57	-133.97	34.52	20.49	-167.44	1.05	0.62
19	20w-10	-167.82	-175.77	-7.95	4.74	-137.44	30.38	18.10	-172.42	-4.60	2.74
20	22w-1	-164.81	-173.43	-8.62	5.23	-124.73	40.08	24.32	-154.64	10.17	6.17
21	30w-1	-209.87	-221.42	-11.55	5.50	-166.73	43.14	20.56	-217.36	-7.49	3.57
22	30w-2	-196.77	-211.88	-15.11	7.68	-157.19	39.58	20.11	-203.64	-6.87	3.49
23	30w-3	-193.34	-215.60	-22.26	11.51	-164.47	28.87	14.93	-193.30	0.04	0.02
	RMSD			8.98			28.91			4.98
	MAD			22.26			43.14			10.17
	MRD(%)				4.85			16.86			2.37

a. E_3B, Δ , RMSD, MAD are in kJ/mol; b. calculated by the CCSD(T)/aug-cc-pVDZ method with the trimer basis; c. calculated with the trimer basis.

Table 3 Calculated total three-body interaction energies of 23 water nanoparticles^a

与CCSD(T)结果相比, 3种方法的RMSD, MAD, MRD大小顺序均为本文方法< MP2方法< M06-2X方法. 可见, 本文方法对纳米尺寸水团簇三体作用能的计算精度稍优于MP2方法, 明显优于M06-2X方法.

上述67个水团簇中, (H₂O)_n(n=11~15, 17~19, 21~22)的结构取自文献[29], 该文献报道了WHBB方法的总三体作用能. WHBB方法使用H原子和O原子构建二至六阶置换不变量, 共涉及5894项. 本文方法使用O— H键偶极构建二至六阶置换不变量, 仅涉及453项. 因此本文方法大大减少了置换不变量的个数. 表4表明, 与CCSD(T)方法的总三体作用能相比, 本文方法的计算精度明显优于WHBB方法.

Table 4 Calculated total three-body interaction energies of 18 water clusters

No.	Cluster	CCSD(T)	WHBB			This work
No.	Cluster	${E_{3 B}}^{*}$ /(kJ· mol^-1)	E_3B/(kJ· mol^-1)	Δ /(kJ· mol^-1)	δ (%)	E_3B/(kJ· mol^-1)	Δ /(kJ· mol^-1)	δ (%)
1	11w-1	-79.54	-81.17	-1.63	2.05	-82.01	-2.47	3.11
2	11w-2	-82.63	-82.84	-0.21	0.25	-82.30	0.33	0.40
3	11w-3	-78.87	-81.59	-2.72	3.45	-79.54	-0.67	0.85
4	11w-4	-82.26	-83.68	-1.42	1.73	-81.92	0.34	0.41
5	11w-5	-81.76	-83.26	-1.50	1.83	-82.09	-0.33	0.40
6	12w-1	-88.70	-93.30	-4.60	5.19	-87.91	0.79	0.89
7	13w-1	-99.87	-103.76	-3.89	3.90	-99.54	0.33	0.33
8	14w-1	-111.55	-115.06	-3.51	3.15	-108.95	2.60	2.33
9	15w-1	-123.39	-123.01	0.38	0.31	-121.92	1.47	1.19
10	17w-1	-142.05	-148.95	-6.90	4.86	-138.87	3.18	2.24
11	17w-2	-146.15	-144.77	1.38	0.94	-142.67	3.48	2.38
12	17w-3	-132.59	-140.16	-7.57	5.71	-128.24	4.35	3.28
13	17w-4	-140.88	-141.84	-0.96	0.68	-139.49	1.39	0.99
14	18w-1	-147.07	-153.97	-6.90	4.69	-141.21	5.86	3.98
15	19w-1	-164.89	-169.03	-4.14	2.51	-155.98	8.91	5.40
16	21w-1	-182.17	-193.30	-11.13	6.11	-181.71	0.46	0.25
17	21w-2	-184.81	-190.37	-5.56	3.01	-179.16	5.65	3.06
18	22w-1	-164.81	-190.37	-25.56	15.51	-154.64	10.17	6.17
	RMSD/(kJ· mol^-1)			7.61			4.14
	MAD/(kJ· mol^-1)			25.56			10.17
	MRD(%)				3.66			2.09

* Calculated by the CCSD(T)/aug-cc-pVDZ method with the trimer basis.

Table 4 Calculated total three-body interaction energies of 18 water clusters

表5给出了不同方法相对于CCSD(T)总三体作用能的误差统计. 图1为本文方法的总三体作用能与CCSD(T)结果间的相关性图.

Table 5 Deviations of different methods with respect to the CCSD(T) three-body interaction energies of 67 water clusters

	Figure Option View Download New Window
	Fig.1 Correlation plots for the total three-body interaction energies of the 67 water clusters (A) MP2/aug-cc-pVDZ; (B) M06-2X/jul-cc-pVTZ; (C) this work.

表5和图1进一步表明本文方法对水团簇三体作用能的计算精度稍优于MP2方法, 明显优于M06-2X方法.

2.2 三体作用能随水团簇构型的变化

图2为水三聚体三体作用能随2个O原子间的距离 $R_{O \dots O}$ 和O— H…O氢键角度 $θ_{O — H \dots O}$ 变化的能量曲线. 棱柱构型的水六聚体^[27]中包含上下2个三元水环,

	Figure Option View Download New Window
	Fig.2 Three-body interaction energy curves with respect to the intermolecular distances(A— C) and angles(D— F)

保持2个三元水环的构型不变, 改变2个环中心间的距离 $R_{r \dots r}$ , 可获得总三体作用能随两环间距离变化的能量曲线(图3).

	Figure Option View Download New Window
	Fig.3 Total three-body interaction energy curves of the prism hexamer

由图2和图3可见, 本文方法可以较好地描述三体作用强度随结构的变化.

2.3 计算效率

比较了不同方法在配置为16核芯CPU(Intel Xeon E5-2670)、 96 GB内存和4 TB硬盘服务器上的计算效率. 对于由5个水分子构成的水团簇5w-1的总三体作用能计算结果表明, 使用CCSD(T)方法所需CPU时间为16769 s, 使用MP2方法所需CPU时间为1571 s, 使用M06-2X方法所需CPU时间为11133 s, 使用本文方法所需CPU时间仅为3 s. 对于由10个水分子构成的水团簇10w-1的总三体作用能计算结果表明, 使用CCSD(T)方法、 MP2方法、 M06-2X方法及本文方法计算该团簇总三体作用能所需CPU时间依次为200627, 18476, 137629和3 s. 可见, 本文方法的计算效率远高于CCSD(T), MP2和M06-2X方法, 并且体系越大, 本文方法的计算优势越明显.

3 结论

提出了一个可快速准确计算水团簇三体作用能的新方法. 使用化学键偶极构建置换不变量, 与使用原子构建置换不变量的WHBB方法相比, 该方法大大地减少了置换不变量的个数, 同时提高了计算精度. 计算结果表明, 该方法不仅可以准确地预测CCSD(T)方法的三体作用能, 而且能够很好地描述三体作用能随水团簇构型的变化. 对较大水团簇三体作用能的计算结果表明, 该方法的计算精度稍优于MP2方法, 明显优于M06-2X方法, 并且运算速度具有明显优势, 体系越大优势越明显.

支持信息见http://www.cjcu.jlu.edu.cn/CN/10.7503/cjcu20190038.

参考文献

文献列表

[1]	Xantheas S. S. , Chem. Phys. , 2000, 258, 225—231 [本文引用:1]
[2]	Milet A. , Moszynski R. , Wormer P. E. S. , Avoird A. V. D. , J. Phys. Chem. A, 1999, 103(34), 6811—6819 [本文引用:1]
[3]	Li S. S. , Jiang X. N. , Wang C. S. , Chem. J. Chinese Universities, 2014, 35(11), 2403—2409 (李书实, 姜笑楠, 王长生. 高等学校化学学报, 2014, 35(11), 2403—2409) [本文引用:1]
[4]	Yang W. , Li X. L. , Wang C. S. , Acta Phys. Chim. Sin. , 2015, 31(12), 2285—2293 (杨微, 李晓蕾, 王长生. 物理化学学报, 2015, 31(12), 2285—2293) [本文引用:1]
[5]	Babin V. , Medders G. R. , Paesani F. , J. Chem. Theory Comput. , 2014, 10(4), 1599—1607 [本文引用:3]
[6]	Ouyang J. F. , Bettens R. P. A. , J. Chem. Theory Comput. , 2015, 11(11), 5132—5143 [本文引用:2]
[7]	Stone A. J. , The Theory of Intermolecular Forces, Oxford University Press, Oxford, 2013, 3 [本文引用:1]
[8]	Thole B. T. , Chem. Phys. , 1981, 59, 341—350 [本文引用:1]
[9]	Kumar R. , Wang F. F. , Jenness G. R. , Jordan K. D. , J. Chem. Phys. , 2010, 132, 014309 [本文引用:1]
[10]	Burnham C. J. , Li J. , Xantheas S. S. , Leslie M. , J. Chem. Phys. , 1999, 110(9), 4566—4581 [本文引用:1]
[11]	Burnham C. J. , Anick D. J. , Mankoo P. K. , Reiter G. F. , J. Chem. Phys. , 2008, 128(15), 154519 [本文引用:1]
[12]	Ren P. , Ponder J. W. , J. Phys. Chem. B, 2003, 107(24), 5933—5947 [本文引用:1]
[13]	Conte R. , Qu C. , Bowman J. M. , J. Chem. Theory Comput. , 2015, 11(4), 1631—1638 [本文引用:1]
[14]	Sun C. L. , Jiang X. N. , Wang C. S. , J. Comput. Chem. , 2009, 30(15), 2567—2575 [本文引用:1]
[15]	Jiang X. N. , Sun C. L. , Wang C. S. , J. Comput. Chem. , 2010, 31(7), 1410—1420 [本文引用:1]
[16]	Li Y. , Jiang X. N. , Wang C. S. , J. Comput. Chem. , 2011, 32(5), 953—966 [本文引用:1]
[17]	Li Y. , Wang C. S. , J. Comput. Chem. , 2011, 32(13), 2765—2773 [本文引用:1]
[18]	Li S. S. , Huang C. Y. , Hao J. J. , Wang C. S. , J. Comput. Chem. , 2014, 35(6), 415—426 [本文引用:1]
[19]	Hao J. J. , Wang C. S. , RSC Adv. , 2015, 5, 6452—6461 [本文引用:1]
[20]	Gao X. C. , Hao. Q. , Wang C. S. , J. Chem. Theory Comput. , 2017, 13(6), 2730—2741 [本文引用:1]
[21]	Huang C. , Hao Q. , Wang C. , Chem. Res. Chinese Universities, 2017, 33(1), 94—99 [本文引用:1]
[22]	Dean J. A. , Lange's Hand book of Chemistry, McGraw-Hill, New York, 1999 [本文引用:1]
[23]	Chen G. D. , Weng J. , Song G. , Li Z. H. , J. Chem. Theory Comput. , 2017, 13(5), 2010—2020 [本文引用:1]
[24]	Leverentz H. R. , Qi H. W. , Truhlar D. G. , J. Chem. Theory Comput. , 2013, 9(2), 995—1006 [本文引用:3]
[25]	Frisch M. J. , Trucks G. W. , Schlegel H. B. , Scuseria G. E. , Robb M. A. , Cheeseman J. R. , Scalmani G. , Barone V. , Mennucci B. , Petersson G. A. , Nakatsuji H. , Caricato M. , Li X. , Hratchian H. P. , Izmaylov A. F. , Bloino J. , Zheng G. , Sonnenberg J. L. , Hada M. , Ehara M. , Toyota K. , Fukuda R. , Hasegawa J. , Ishida M. , Nakajima T. , Honda Y. , Kitao O. , Nakai H. , Vreven T. , Montgomery J. A. Jr. , Peralta J. E. , Ogliaro F. , Bearpark M. , Heyd J. J. , Brothers E. , Kudin K. N. , Staroverov V. N. , Keith T. , Kobayashi R. , Normand J. , Raghavachari K. , Rendell A. , Burant J. C. , Iyengar S. S. , Tomasi J. , Cossi M. , Rega N. , Millam J. M. , Klene M. , Knox J. E. , Cross J. B. , Bakken V. , Adamo C. , Jaramillo J. , Gomperts R. , Stratmann R. E. , Yazyev O. , Austin A. J. , Cammi R. , Pomelli C. , Ochterski J. W. , Martin R. L. , Morokuma K. , Zakrzewski V. G. , Voth G. A. , Salvador P. , Dannenberg J. J. , Dapprich S. , Daniels A. D. , Farkas O. , Foresman J. B. , Ortiz J. V. , Cioslowski J. , Fox D. J. , Gaussian 09, Revision D. 01, Gaussian Inc. , Wallingford CT, 2013 [本文引用:1]
[26]	R$\check{e}$zá$\check{c}$ J. , Huang Y. , Hobza P. , Beran G. J. O. , J. Chem. Theory Comput. , 2015, 11(7), 3065—3079 [本文引用:1]
[27]	Reddy S. K. , Straight S. C. , Bajaj P. , Pham C. H. , Riera M. , Moberg D. R. , Morales M. A. , Knight C. , Götz A. W. , Paesani F. , J. Chem. Phys. , 2016, 145(19), 194504 [本文引用:2]
[28]	Yuan D. , Shen X. , Li W. , Li S. , Phys. Chem. Chem. Phys. , 2016, 18(24), 16491—16500 [本文引用:2]
[29]	Wang Y. , Huang X. , Shepler B. C. , Braams B. J. , Bowman J. M. , J. Chem. Phys. , 2011, 134, 094509 [本文引用:3]
[30]	Ouyang J. F. , Bettens R. P. A. , J. Chem. Theory Comput. , 2016, 12(12), 5860—5867 [本文引用:1]

2000

0.0

... 研究发现, 二体作用对水团簇稳定性的贡献最大, 三体作用对水团簇稳定性的贡献也达20%以上^[1~4], 因此仅仅准确描述分子间的二体作用不足以保证对生物分子体系稳定性的准确预测, 准确计算水分子间存在的诸多三体作用强度十分重要 ...

1999

0.0

2014

0.0

2015

0.0

2014

0.0

... Paesani等^[5]发现, 对于水团簇体系, 使用CCSD(T)/aug-cc-pVTZ方法结合三体基函数(Trimer basis)计算得到的三体作用能与CCSD(T)/CBS的计算结果相近 ...

... 受Paesani等^[5]和Li等^[23]工作的启发, 将开关函数s设为, ...

... 选取文献[5]报道的12347个水三聚体结构和CCSD(T)/aug-cc-pVTZ三体作用能作为标准数据, 改变非线性参数R_cut, #cod#x00394 ...

2015

0.0

... Bettens等^[6]发现, 使用MP2/aug-cc-pVDZ方法结合三体基函数可以得到与MP2/CBS方法精度相当的三体作用能 ...

... 仿照多体展开理论对三体作用能的定义^[6,7], 由水团簇中第m, n, l个水分子组成的三聚体的长程诱导作用 V3Bind的计算如下: ...

2013

0.0

... 仿照多体展开理论对三体作用能的定义^[6,7], 由水团簇中第m, n, l个水分子组成的三聚体的长程诱导作用 V3Bind的计算如下: ...

1981

0.0