一種基于分塊迭代法的高精度原子電荷計(jì)算方法

李焱，湯天寶，金樂(lè)，李冬婷，李國(guó)輝*，許佩軍*

1.中國(guó)科學(xué)院大連化學(xué)物理研究所，分子反應(yīng)動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023

2.遼寧師范大學(xué)，物理與電子技術(shù)學(xué)院，遼寧 大連 116023

一種基于分塊迭代法的高精度原子電荷計(jì)算方法

李焱1，湯天寶2，金樂(lè)2，李冬婷2，李國(guó)輝1*，許佩軍2*

1.中國(guó)科學(xué)院大連化學(xué)物理研究所，分子反應(yīng)動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023

2.遼寧師范大學(xué)，物理與電子技術(shù)學(xué)院，遼寧 大連 116023

原子電荷是研究生物分子的重要物理化學(xué)參數(shù)。本文在系統(tǒng)分塊法的基礎(chǔ)上，通過(guò)迭代的量化計(jì)算，得到多分子體系在特定環(huán)境下的高精度原子電荷，并提出一套完整的多分子體系原子電荷計(jì)算方法，Partition Iterative Quantum(PIQ)。本文采用若干種主要的生物體系作為測(cè)試體系，驗(yàn)證了所提方法的可行性。并以完全的量子化學(xué)計(jì)算得到的原子電荷為標(biāo)準(zhǔn)，比較了不同初始電荷設(shè)置方式對(duì)計(jì)算精度和效率的影響。說(shuō)明了所提方法是可以被應(yīng)用于不同尺度的分子體系的電荷計(jì)算中。

原子電荷；量化計(jì)算；分子力場(chǎng)；分塊法

引言

相對(duì)于傳統(tǒng)的實(shí)驗(yàn)科學(xué)，利用計(jì)算機(jī)研究生物分子結(jié)構(gòu)和功能的計(jì)算化學(xué)可以很大程度上提高研究生物分子的效率和節(jié)省資源。隨著計(jì)算機(jī)硬件的發(fā)展和相關(guān)算法的成熟，利用計(jì)算化學(xué)技術(shù)研究生物大分子的結(jié)構(gòu)和功能變得越來(lái)越不可或缺。特別是在探索生物體系的反應(yīng)機(jī)理，輔助藥物設(shè)計(jì)等方面，計(jì)算化學(xué)已經(jīng)成為常規(guī)手段之一[1]。

基于分子力場(chǎng)(Force Field)[2]的分子動(dòng)力學(xué)模擬是計(jì)算化學(xué)的重要組成部分。通過(guò)迭代求解分子體系內(nèi)的原子位置隨時(shí)間變化的牛頓動(dòng)力學(xué)方程，它可以直接從原子層面上揭示生物大分子結(jié)構(gòu)變化或功能行使的機(jī)制[3]。所謂分子力場(chǎng)，是指用一系列經(jīng)驗(yàn)性的函數(shù)來(lái)描述在體系內(nèi)原子之間的相互作用。一般來(lái)說(shuō)，分子力場(chǎng)包含若干鍵連相互作用項(xiàng)(如鍵伸縮項(xiàng)、鍵角運(yùn)動(dòng)項(xiàng)和平面角運(yùn)動(dòng)項(xiàng)等等)和非鍵相互作用項(xiàng)。其中非鍵相互作用由靜電相互作用和范德華相互作用組成。為了提高計(jì)算效率，人們一般使用位于原子中心的固定點(diǎn)電荷(以下將簡(jiǎn)稱為原子電荷)之間的庫(kù)倫相互作用來(lái)代表靜電項(xiàng)。由于生物大分子的運(yùn)動(dòng)一般是由非鍵相互作用中靜電項(xiàng)來(lái)主導(dǎo)，所以原子電荷的精度直接影響了分子動(dòng)力學(xué)模擬的結(jié)果。因此，得到生物大分子體系在不同的環(huán)境下的精確原子電荷具有非常重要的理論研究意義和實(shí)際應(yīng)用價(jià)值。迄今為止，人們針對(duì)不同的體系或不同的計(jì)算目標(biāo)發(fā)展了很多套分子力場(chǎng)(比如 AMBER[4]，CHARMM[5-8]，OPLSA 等力場(chǎng)[9-11])。一般情況下，在生物大分子的模擬過(guò)程中，原子電荷是依賴于分子力場(chǎng)實(shí)現(xiàn)的。同一原子類型在這些不同力場(chǎng)下對(duì)應(yīng)的原子電荷并不相同，但大都通過(guò)擬合靜電勢(shì)獲得，可以有效的代表特定體系在特定條件下的平均帶點(diǎn)性質(zhì)。然而，在實(shí)際應(yīng)用中，隨著體系構(gòu)象和環(huán)境的變化，屬于原子的有效電荷會(huì)因?yàn)檎T導(dǎo)和極化效應(yīng)而發(fā)生變化。當(dāng)誘導(dǎo)和極化效應(yīng)顯著時(shí)，固定點(diǎn)電荷模型就不再準(zhǔn)確了。為了解決這一問(wèn)題，諸多學(xué)者提出了可極化的靜電模型。Ponder 和 Ren 發(fā)展了 AMOEBA 力場(chǎng),這一力場(chǎng)可以通過(guò)電多極距展開(kāi)和誘導(dǎo)偶極來(lái)描述模擬過(guò)程中的誘導(dǎo)極化效應(yīng)，使得對(duì)原子電荷的描述提升到了一個(gè)更高的精度[12,13]。楊和宮等人提出的ABEEM 相關(guān)方法，通過(guò)原子鍵電負(fù)性均衡理論建立方程組，求解浮動(dòng)的電荷，也得到了很好的效果[14-15]。Benoit 等人發(fā)展了可極化的 Drude 模型[16-17]。另外，還有許多學(xué)者在這方面也作為深入的研究并得到了令人興奮的結(jié)果[18]。

量子化學(xué)是量子力學(xué)在生物化學(xué)領(lǐng)域的應(yīng)用。它是由Heisenberg 和 Schr?dinger等人在20世紀(jì)20年代建立的知識(shí)體系。幾十年來(lái)，量子化學(xué)已經(jīng)得到了長(zhǎng)足的發(fā)展[19]。量子化學(xué)基于薛定諤方程，在非相對(duì)論近似，玻恩–奧本海默近似和單電子近似的條件下，利用變分原理求解由基函數(shù)線性組合而成的電子波函數(shù)，或者基于 Hohenberg-Kohn定理，通過(guò)密度泛函理論求解電子密度。雖然使用量子化學(xué)手段可以計(jì)算特定構(gòu)象的電子密度，進(jìn)而通過(guò)布居分析可以得到比前述分子力場(chǎng)固定點(diǎn)電荷更為可靠的原子電荷，但是巨大的計(jì)算量嚴(yán)重阻礙了其在生物大分子領(lǐng)域的應(yīng)用。

分塊法通過(guò)將生物大分子劃分成相對(duì)較小的片段，分別由量化方法計(jì)算其高精度的物理化學(xué)性質(zhì)，并通過(guò)特定的組合方式得到整體的物理化學(xué)性質(zhì)。由于計(jì)算小片段比計(jì)算整個(gè)體系要更易實(shí)現(xiàn)并且方便通過(guò)并行計(jì)算提高計(jì)算效率，分塊法已經(jīng)發(fā)展出了眾多的實(shí)現(xiàn)方式，并應(yīng)用于研究生物大分子的動(dòng)態(tài)演化過(guò)程中[20]。其中主要的方法有 Divide and Conquer approach[21,22]，Molecular Tailoring approach[23]，和Molecular Fractionatiion with Conjugate Caps(MFCC)[24-26]以及在MFCC的基礎(chǔ)上應(yīng)用于蛋白質(zhì)和核酸能量計(jì)算的改進(jìn)方法[27-30]，另外，還有高加力的 X-POL方法[31-33]等。

本文基于系統(tǒng)分塊法的思想，提出一種通過(guò)迭代量化計(jì)算得到多分子體系中原子電荷的方法。以完全的量化所得到的電荷為標(biāo)準(zhǔn)，在相對(duì)較少的時(shí)間內(nèi)，不需要消耗更多的計(jì)算資源，得到較高精度的原子電荷。為以量子化學(xué)計(jì)算精度實(shí)現(xiàn)分子動(dòng)力學(xué)模擬提供了一種可能的實(shí)現(xiàn)思路。

在第一節(jié)中，詳細(xì)闡述了上述方法的計(jì)算思路與計(jì)算步驟，并分別以雙分子和多分子體系加以說(shuō)明。在第二節(jié)中，本文采用生物大分子的主要組成部分和存在環(huán)境的單元所組成若干體系作為測(cè)試體系，驗(yàn)證所提出方法的可行性。在第三節(jié)中，以完全的量子化學(xué)計(jì)算為標(biāo)準(zhǔn)，討論了不同的原子電荷初始方式對(duì)計(jì)算效率和計(jì)算精度的影響。在之后的結(jié)語(yǔ)部分，概括了本方法的優(yōu)勢(shì)與不足以及可能改進(jìn)的方向。

1 分塊迭代量化計(jì)算方法(PIQ)

1.1 PIQ計(jì)算方法概述

通過(guò)第一性原理來(lái)研究分子(或原子)的物理化學(xué)性質(zhì)，是生物化學(xué)領(lǐng)域的主要研究方法之一，由于受到計(jì)算時(shí)間和計(jì)算資源的約束，量子力學(xué)方法不容易實(shí)現(xiàn)計(jì)算大分子體系中原子的精確物理化學(xué)性質(zhì)。為了能夠高效的得到可以響應(yīng)環(huán)境變化的原子電荷，本文發(fā)展了一種基于體系分塊和對(duì)子體系進(jìn)行量子化學(xué)計(jì)算的迭代計(jì)算方法，用于計(jì)算多分子體系中不同構(gòu)象的原子電荷，稱之為分塊迭代量化計(jì)算方法Partition Iterative Quantum(PIQ)。在由多個(gè)分子組成的生物體系中，按分子個(gè)體將多分子體系劃分為不同的分子片段。隨機(jī)選擇其中一個(gè)分子片段為目標(biāo)分子，以其他分子片段的原子電荷作為背景電荷，通過(guò)高精度的量子化學(xué)計(jì)算得到目標(biāo)分子的原子電荷，完成一次分塊的迭代計(jì)算。重新選擇體系中的一個(gè)分子片段，以之前得到的原子電荷和非目標(biāo)分子片段的原子電荷作為背景電荷，計(jì)算新的目標(biāo)分子的原子電荷。如此往復(fù)，進(jìn)行迭代計(jì)算，直到多分子體系中計(jì)算得到的所有分子片段中的原子電荷不再改變，迭代計(jì)算完成。通過(guò)將所有分子片段的原子電荷統(tǒng)計(jì)，得到整個(gè)多分子體系的原子電荷，下面，以二分子體系和三分子體系進(jìn)行PIQ方法的說(shuō)明。

1.2 二分子體系的PIQ計(jì)算

圖1 中的二分子體系是由蛋白質(zhì)中的氨基酸片段酪氨酸(Tyrosine，TYR)和核酸中的脫氧核苷酸 DNA中的腺嘌呤單鏈 DA 片段組成。TYR 是一種不帶電的極化氨基酸片段[34]，在之后的驗(yàn)證體系介紹中，有更詳細(xì)的說(shuō)明。

圖1 由TYR片段和DA片段組成的雙分子體系的PIQ計(jì)算示意圖。Fig.1 The illustration ofPIQmethod to calculate the atom charges in two molecular systems composed of amino acide TYR fragement and DA fragement.

在第一次迭代時(shí)，我們以TYR中的原子電荷為背景電荷(在圖1 中，以 TYR 原子周圍的電勢(shì)表面表示)，TYR 分子片段和 DA 片段均采用球棍的方式顯示，通過(guò)Gaussian View 05軟件產(chǎn)生[35]。通過(guò)量子化學(xué)計(jì)算軟件Gaussian，計(jì)算DA片段中的各原子的電荷，具體的計(jì)算方法與計(jì)算參數(shù)設(shè)置在之后章節(jié)有描述。在第二次迭代時(shí)，以之前計(jì)算得到的DA片段中的各原子電荷作為背景電荷，計(jì)算TYR片段的原子電荷，通過(guò)兩次迭代過(guò)程中，分子體系中所有原子電荷的均方根誤差 Root Mean Square Deviation(RMSD)判斷本次雙分子體系與上一次迭代時(shí)雙分子體系中所有的原子電荷的差異。如果差異小于預(yù)先設(shè)定的閾值，則退出迭代計(jì)算并輸出體系的原子電荷。否則，繼續(xù)迭代，直至差異小于閾值，完成PIQ對(duì)兩分子體系的迭代電荷計(jì)算。

1.3 多分子體系的PIQ計(jì)算

因?yàn)榻^大多數(shù)的生物分子反應(yīng)是在溶液中進(jìn)行的，所以，對(duì)于溶液，特別是水分子的原子電荷的研究，對(duì)準(zhǔn)確的刻畫(huà)和表達(dá)溶劑的作用有非常重要的意義。所以為了說(shuō)明如何通過(guò)PIQ方法計(jì)算多分子體系的原子電荷，本文以相對(duì)容易實(shí)現(xiàn)的三個(gè)水分子為例，依據(jù)圖2所示進(jìn)行說(shuō)明。

圖2 三個(gè)水分子組成的多分子體系的PIQ示意圖。Fig.2 The illustratiion ofPIQmethod to calculate the atom charges in multi-molecular systems composed of three water moleculers.

PIQ 方法具體的計(jì)算過(guò)程與二分子體系基本相同，在每次迭代過(guò)程中，隨機(jī)選擇一個(gè)水分子作為目標(biāo)分子片段，以剩余兩個(gè)水分子的原子電荷作為背景電荷，在新的背景電荷的環(huán)境下，重新選擇目標(biāo)水分子。通過(guò)量子力學(xué)計(jì)算得到目標(biāo)水分子的原子電荷。通過(guò)迭代，不斷完成量子力學(xué)的電荷計(jì)算過(guò)程，直到多個(gè)水分子的電荷不再改變。

1.4 PIQ 計(jì)算流程圖

雙分子體系和三分子體系只是在每次迭代時(shí)，選擇的目標(biāo)片段以及背景電荷有所差異，本質(zhì)上都是通過(guò)劃分系統(tǒng)，以達(dá)到可以通過(guò)相對(duì)少的計(jì)算資源快速得出精度較高的原子電荷的目標(biāo)，將上述過(guò)程歸納，給出PIQ方法的計(jì)算流程圖，具體的，如圖3所示。

2 驗(yàn)證體系與計(jì)算設(shè)置

因?yàn)樘岢鯬IQ方法的初衷是通過(guò)高精度的原子電荷參與到計(jì)算生物大分子的動(dòng)力學(xué)模擬過(guò)程中，所以，我們盡可能的通過(guò)生物大分子的組成單元和其所處環(huán)境所組成的體系來(lái)驗(yàn)證所提出的PIQ方法的適用性。

圖3 PIQ方法的計(jì)算流程圖Fig.3The fl owchart ofPIQmethod

2.1 驗(yàn)證體系組成單元

蛋白質(zhì)是生物分子的必要組成成分，它參與了幾乎細(xì)胞生命活動(dòng)的絕大多數(shù)進(jìn)程。蛋白質(zhì)是由20種主要的氨基酸線性組成[36]，我們選取其中的若干種有代表性的作為測(cè)試體系的組成單元。將所選擇的氨基酸片段按照極性與非極性以及氨基酸是否帶電荷進(jìn)行分類，具體如表1所示。

核酸，作為另外一類組成生物大分子的重要單元，負(fù)責(zé)生物體遺傳信息的攜帶和傳遞，是生物體遺傳信息的載體。金屬離子廣泛存在于生物體中，它們參與生物體的物質(zhì)運(yùn)送，能量轉(zhuǎn)換，信息傳遞等眾多生命過(guò)程，是生命體系維持正?；顒?dòng)的重要因素。水幾乎是生物分子所有生命過(guò)程的參與者，是重要的溶劑和介質(zhì)環(huán)境。上述三種測(cè)試體系的組成單元，如表2所示。

由于所選的測(cè)試體系單元較多，將所有的組合用于驗(yàn)證過(guò)程并不現(xiàn)實(shí)，因此，我們?cè)谒念悳y(cè)試單元中選取有代表性的單元個(gè)體組成雙分子體系。具體的如圖4中的下三角矩陣圖形所示，相同顏色為組成體系所采用的大類相同。在之后的分析過(guò)程中，我們將DNA和 RNA 作為兩部分，分別進(jìn)行討論。

2.2 分子相對(duì)位置

對(duì)于由蛋白質(zhì)-蛋白質(zhì)組成的雙分子體系。我們采用了不同的空間相對(duì)位置，以檢驗(yàn)不同相對(duì)位置對(duì)PIQ 結(jié)果的影響，兩種不同的相對(duì)位置分別是平行放置和前后放置，如圖5的A所示，在平行放置時(shí)，基本保持氨基酸的側(cè)鏈指向相同方向，前后放置時(shí)，基本保持主鏈在相同方向上。對(duì)于其它的雙分子測(cè)試體系，為了考察分子相對(duì)位置對(duì)計(jì)算結(jié)果的影響，將其中一個(gè)分子片段A固定，將另一個(gè)分子片段B在3-7 埃的最短距離范圍內(nèi)，隨機(jī)擺放，生成5組測(cè)試體系，如圖5的B所示，用于驗(yàn)證。其中，兩個(gè)分子片段之間的距離Dis(fragA，fragB)被定義為兩個(gè)分子片段中，原子之間的距離集合的最小值，如(1)式所示。

表1 組成測(cè)試體系的氨基酸單元Table 1 The amino acid units to compose the benchmark of two molecular system

表2 組成測(cè)試體系的核酸單元、離子和水分子Table 2 The nucleic acid,ion and water composed the benchmark of two and three moleculars system

其中，dis(·)是歐氏距離，atomi是分子片段 fragA中的原子，atomj是分子片段 fragB中的原子。

對(duì)于三分子體系，固定一個(gè)水分子，將另一個(gè)水分子向距離另兩個(gè)水分子更遠(yuǎn)的方向移動(dòng)0.5埃，共計(jì)移動(dòng)10次，在每次移動(dòng)過(guò)程中，將第3個(gè)水分子向距離其它兩個(gè)水分子更遠(yuǎn)的方向移動(dòng)9次。如圖5的C所示，對(duì)于每一組相對(duì)位置的測(cè)試體系，均進(jìn)行PIQ的迭代計(jì)算。

圖4 驗(yàn)證體系的組成以及用于驗(yàn)證的雙分子和三分子體系Fig.4 The componds of benchmark systems and two molecular and three molecular systems which are applied as test cases

圖5A: 蛋白質(zhì)-蛋白質(zhì)體系的不同擺位置。B: 生成雙分子體系的相對(duì)位置。C: 生成多分子體系的相對(duì)位置Fig.5 The position of protein-protein system.The different positions of two molecular system.The different positions of three molecular system

2.3 計(jì)算參數(shù)設(shè)置

本論文中的所有量子力學(xué)計(jì)算均在基于 Linux 的CentOs 5.6系統(tǒng)下通過(guò)量化計(jì)算軟件 Gaussian 09[37]實(shí)現(xiàn)。計(jì)算中，方法選擇 MP2，基組為 6-311++G**。核酸相關(guān)的計(jì)算使用12核進(jìn)行并行計(jì)算，其它的測(cè)試體系均勻采用8核實(shí)現(xiàn)并行計(jì)算，CPU 為Intel Xeon E5-2650。內(nèi)存根據(jù)計(jì)算的測(cè)試體系中原子個(gè)數(shù)的不同在8G-32G之間設(shè)置。電荷從計(jì)算結(jié)果中提取Mulinken 電荷，判斷兩次計(jì)算的差異的閾值，設(shè)置為0.0001。

2.4 原子電荷初始設(shè)置

在PIQ方法的計(jì)算過(guò)程中，對(duì)分子片段的初始原子電荷設(shè)置是非常重要的步驟。我們以三種不同的方式設(shè)置初始背景電荷討論如何設(shè)置初始原子電荷才能使迭代次數(shù)盡可能的少，能盡快達(dá)到需要的精度。

01.Gaussian 初始電荷：對(duì)于作為非目標(biāo)分子的分子片段，通過(guò)量子力學(xué)單獨(dú)計(jì)算其原子電荷，計(jì)算所采用的基組和方法如 2.3 節(jié)所述。

02.CHARMM 初始電荷：通過(guò)已有的力場(chǎng)中的電荷參數(shù)，設(shè)置非目標(biāo)分子片段的原子電荷。因?yàn)镃HARMM力場(chǎng)對(duì)從小分子到溶劑化的生物大分子均有很好的物理化學(xué)性質(zhì)的描述，所以本文采用CHARMM22CMP力場(chǎng)[38]。

03.Random 初始電荷：在滿足一定約束條件的前提下，隨機(jī)設(shè)置非目標(biāo)分子片段的原子電荷，其滿足的約束條件如(2-4)式。

其中，是第個(gè)原子的初始電荷，并且滿足

所有的原子電荷的加和值的符號(hào)與分子片段的總帶電荷數(shù)的符號(hào)一致。

3 數(shù)值實(shí)驗(yàn)結(jié)果

在驗(yàn)證過(guò)程中，用于比較的真值是通過(guò) Gaussian 09計(jì)算整體多分子體系得到的原子電荷?；M和方法如2.3節(jié)所示。對(duì)于相同組成的測(cè)試體系，用于對(duì)比的是三種不同的初始原子電荷設(shè)置方式，不同的空間相對(duì)位置以及初始分子片段的選擇所對(duì)應(yīng)的測(cè)試?yán)印?/p>

對(duì)于每個(gè)測(cè)試體系，以迭代結(jié)束時(shí)的體系的原子電荷與真值進(jìn)行比較分析。

用于比較的 RMSD 的計(jì)算(5)所示。

其中，m 為體系分子片段的個(gè)數(shù)，ni為第i個(gè)分子片段中的原子個(gè)數(shù)。為第 i個(gè)分子片段的第j個(gè)原子以pat的方式進(jìn)行初始化時(shí)，經(jīng)過(guò)迭代計(jì)算得到的原子電荷。上標(biāo) pat 分別為 Gaussian，CHARMM和 Random 三種方式。

3.1 原子電荷的回歸分析

在眾多的計(jì)算測(cè)試體系中，選取和蛋白質(zhì)相關(guān)的四個(gè)測(cè)試體系。蛋白質(zhì)-蛋白質(zhì)(LEU-LYS，圖6A，蛋白質(zhì)-核酸(TYR-DA，圖6B)，蛋白質(zhì)-離子(ARG-Ca，圖6C)，蛋白質(zhì)-水(GLU-H2O，圖6D)進(jìn)行原子電荷計(jì)算結(jié)果的說(shuō)明。在圖6的四個(gè)子圖中，橫軸為完全由量子力學(xué)計(jì)算得到的整體體系中的原子電荷，縱軸為經(jīng)過(guò)迭代達(dá)到收斂精度時(shí)，整個(gè)體系的原子電荷。不同形狀(顏色)的點(diǎn)分別代表三種不同的原子電荷初始設(shè)置方式。在四個(gè)測(cè)試?yán)又校?jīng)過(guò)若干次的迭代計(jì)算，不同的初始方式的PIQ方法都能得到與完全量子力學(xué)計(jì)算相接近的原子電荷數(shù)值。只是各種初始方式所消耗的迭代次數(shù)不同。

在絕大多數(shù)的測(cè)試體系中，體現(xiàn)的結(jié)果與圖6的子圖的結(jié)果相似。不同體系中的PIQ得到的原子電荷與完全量子力學(xué)得到的原子電荷有非常好的線性相關(guān)性，只有極個(gè)別的原子，電荷的差異值相對(duì)較大，這在RMSD的平均統(tǒng)計(jì)中有所體現(xiàn)。

3.2 電荷RMSD統(tǒng)計(jì)分析

由于每類測(cè)試單元的組合較多，通過(guò) RMSD 的均值來(lái)說(shuō)明PIQ方法所能達(dá)到的總體效果。在圖7中，不同的測(cè)試類型通過(guò)下三角矩陣體現(xiàn)，A，B，C中網(wǎng)格的顏色根據(jù)RMSD的均值體現(xiàn)，數(shù)值顯示在對(duì)應(yīng)的網(wǎng)格中，計(jì)算如(6)式。

圖6 蛋白質(zhì)相關(guān)的測(cè)試體系的原子電荷回歸圖。A: LEULYS; B: TYR-DA; C: ARG-Ca; D: GLU-H2OFig.6 The congression of atom charges of protein related test cases,A: LEU-LYS B: TYR-DA C: ARG-Ca D: GLU-H2O

其中，每類測(cè)試單元的測(cè)試個(gè)數(shù)ncase在圖7 的 D 中顯示。

不同的初始原子電荷設(shè)置方式經(jīng)過(guò)迭代所達(dá)到的精度差異不明顯。除個(gè)別原子的電荷外，絕大部分的原子電荷都能達(dá)到完全量子力學(xué)的計(jì)算精度。

蛋白質(zhì)相關(guān)的PIQ計(jì)算中除了蛋白質(zhì)-離子組合誤差較大。主要是 ARG-Mg,ARG-Zn 測(cè)試體系，在去除這幾個(gè)明顯的誤差點(diǎn)的情況下，RMSD 的平均值由0.193 下降至 0.11。

3.3 迭代完成次數(shù)分析

圖7 不同測(cè)試體系的 RMSD 的平均值的分布圖。A:Gaussian 原子電荷。B: CHARMM 原子電荷。C: Random 原子電荷。D: 每類測(cè)試單元的測(cè)試體系的個(gè)數(shù)Fig.7 The mean(RMSD)of different test castes via three initial atom charge patterns.A: Gaussian atom charge.B: CHARMM atom charge.C: Random atom charge.D: The number of test cases

雖然不同的初始原子設(shè)置方式在迭代后達(dá)到幾乎相同的精度，但盡可能少的消耗計(jì)算資源和計(jì)算時(shí)間，實(shí)現(xiàn)迭代達(dá)到要求精度，正是PIQ提出的初衷。圖8的ABC三個(gè)子圖分別是 Gaussian，CHARMM，Random 的對(duì)于不同測(cè)試體系的迭代次數(shù)的統(tǒng)計(jì)分析柱狀圖。其中，不同迭代次數(shù)ite所占同類測(cè)試體系總迭代次數(shù)的百分比。因?yàn)閳D8 A 是以單獨(dú)的量子力學(xué)為基礎(chǔ)的(沒(méi)有考慮到環(huán)境中的電荷信息)，其本身就比較接近分子片段中的原子電荷的真值，所以相對(duì)的，迭代次數(shù)的分布以ite=2，3為主。只有三分子體系，因?yàn)榉肿悠蝹€(gè)數(shù)多，所以迭代次數(shù)比較高?；?CHARMM 的原子電荷也是基于大量的量化計(jì)算并依賴于實(shí)驗(yàn)擬合最終得到，但不同的分子片段的原子電荷無(wú)法達(dá)到完全量子力學(xué)計(jì)算得到的精度，迭代次數(shù)的分布以 ite=3，4 為主，少數(shù)的兩次迭代便能達(dá)到精度?；?Random 的原子電荷除分子片段總電荷符號(hào)與分子片段相同，其它完全與用于背景的原子電荷沒(méi)有相關(guān)性，在開(kāi)始的2-3次迭代時(shí)，完全是通過(guò)目標(biāo)分子片段的原子類型在非正確的背景電荷的情況下實(shí)現(xiàn)計(jì)算，所以，在第三種情況下，迭代次數(shù)消耗比較多，以ite=3，4為主要部分。

初始原子電荷的設(shè)置方式不會(huì)影響PIQ計(jì)算的精度，但接近背景分子片段的原子電荷能使迭代向收斂較快的方向進(jìn)行。但初始計(jì)算分子片段的原子電荷也要消耗計(jì)算資源和時(shí)間，所以，當(dāng)無(wú)法得到精確的背景分子片段的原子電荷時(shí)，基于分子力場(chǎng)的作為初始背景電荷，可以平衡計(jì)算精度和計(jì)算資源，是PIQ計(jì)算方法的首選方式。

3.4 PIQ 計(jì)算時(shí)間分析

對(duì)體系整體的量化計(jì)算需要消耗大量的計(jì)算時(shí)間，雖然對(duì)分子片段的計(jì)算所需要時(shí)間相對(duì)比較少，但多次的迭代也會(huì)相應(yīng)的增加PIQ總體的計(jì)算時(shí)間。對(duì)于原子個(gè)數(shù)比較少的體系，特別是對(duì)由水分子組成的三分子體系，經(jīng)過(guò)多次迭代才能收斂，PIQ方法與完全量化計(jì)算在計(jì)算時(shí)間上優(yōu)勢(shì)并不明顯，甚至是要更長(zhǎng)一點(diǎn)。在DA片段和LYS片段的一個(gè)測(cè)試?yán)又?，完全的量化需要消耗大約67個(gè)CPU小時(shí)，而三種初始方式的PIQ方法，在以DA片段中的三種原子電荷作為背景電荷的情況下，分別需要消耗大約8小時(shí)(Gaussian)，8小時(shí)(CHARMM)和10小時(shí)(Random)CPU計(jì)算時(shí)間，PIQ方法的計(jì)算效率顯著提升。所以，隨著體系中原子數(shù)目逐漸增加，PIQ方法的優(yōu)勢(shì)越明顯。

圖8 三種不同的原子電荷初始試，迭代次數(shù)的統(tǒng)計(jì)。A: Gaussian 初始電荷。B: CHARMM 初始電荷。C: Random 初始電荷.D: 迭代次數(shù)對(duì)應(yīng)的 colormap.Fig.8 The statics of three different initial pattern of atom charge.A: Gaussian.B: CHARMM.C: Random.D: Iteration colormap.

4 結(jié)語(yǔ)

PIQ 是通過(guò)迭代分塊計(jì)算分子體系原子電荷的計(jì)算方法。PIQ 基本思路就是將不容易或無(wú)法完全由量子化學(xué)計(jì)算的多分子體系分割成若干相對(duì)較小的容易計(jì)算的子體系，采用背景電荷保持各分子片段之間的聯(lián)系，通過(guò)迭代計(jì)算得到接近完全的量子化學(xué)的精度，這使得將量化精度的電荷應(yīng)用于分子動(dòng)力學(xué)模擬成為可能。本文采用的測(cè)試單元，均是原子個(gè)數(shù)比較少的生物分子片段，對(duì)于大分子或是多分子體系，完全可以通過(guò)相同的計(jì)算策略實(shí)現(xiàn)原子電荷的計(jì)算。對(duì)于更多分子的分子體系，如何通過(guò)并行，進(jìn)一步實(shí)現(xiàn)計(jì)算速度和效率的提升，是一步PIQ方法要嘗試解決的主要問(wèn)題。

致謝

本項(xiàng)工作由國(guó)家自然科學(xué)基金提供資助，基金號(hào)為 21573217，91430110，31370714，21625302。

[1] 張明波,李國(guó)輝,計(jì)算化學(xué)在生命科學(xué)中的應(yīng)用.科研信息化技術(shù)與應(yīng)用,2010.1(1): 60-69.

[2] Guvench,O.,A.D.MacKerell,Comparison of Protein Force Fields for Molecular Dynamics Simulations,in Molecular Modeling of Proteins,A.Kukol,Editor.2008,Humana Press: Totowa,NJ.63-88.

[3] Leach,A.,Molecular Modelling: Principles and Applications(2nd Edition).2 ed.2001: Prentice Hall.

[4] Duan,Y.,C.Wu,S.Chowdhury,M.C.Lee,G.Xiong,W.Zhang,R.Yang,P.Cieplak,R.Luo,T.Lee,J.Caldwell,J.Wang,and P.Kollman,A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations.Journal of Computational Chemistry,2003.24(16): 1999-2012.

[5] Jo,S.,X.Cheng,J.Lee,S.Kim,S.-J.Park,D.S.Patel,A.H.Beaven,K.I.Lee,H.Rui,S.Park,H.S.Lee,B.Roux,A.D.MacKerell,J.B.Klauda,Y.Qi,and W.Im,CHARMM-GUI 10 years for biomolecular modeling and simulation.Journal of Computational Chemistry,2017.38(15): 1114-1124.

[6] Zhu,X.,P.E.M.Lopes,A.D.MacKerell,Recent Developments and Applications of the CHARMM force fields.Wiley interdisciplinary reviews.Computational molecular science,2012.2(1): 167-185.

[7] Reiher,W.E.,Theoretical studies of hydrogen bonding.1985.

[8] Nilsson,L.,M.Karplus,Empirical energy functions for energy minimization and dynamics of nucleic acids.Journal of Computational Chemistry,1986.7(5): 591-616.

[9] Jorgensen,W.L.,D.S.M.And,J.Tiradorives,Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids.Journal of the American Chemical Society,2017.118(45): 11225-11236.

[10] Weiner,S.J.,P.A.Kollman,D.T.Nguyen,D.A.Case,An all atom force fi eld for simulations of proteins and nucleic acids.Journal of Computational Chemistry,1986.7(2):230-252.

[11] And,G.A.K.,R.A.Friesner,T.R.And,W.L.Jorgensen,Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides?.Journal of Physical Chemistry B,2001.105(28): 6474-6487.

[12] Shi,Y.,Z.Xia,J.Zhang,R.Best,C.Wu,J.W.Ponder,and P.Ren,Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins.Journal of Chemical Theory and Computation,2013.9(9): 4046-4063.

[13] Ren,P.,C.Wu,J.W.Ponder,Polarizable Atomic Multipole-Based Molecular Mechanics for Organic Molecules.Journal of Chemical Theory and Computation,2011.7(10): 3143-3161.

[14] 宮利東,ABEEM浮動(dòng)電荷分子力場(chǎng)在含離子體系中的發(fā)展和應(yīng)用.中國(guó)科學(xué):化學(xué),2013(05): 519-535.

[15] Yang,Z.-Z.,X.-T.Lin,D.-X.Zhao,Fast Calculation of Molecular Total Energy with ABEEMσπ/MM Method------ For Some Series of Organic Molecules and Peptides.Chemical Physics,2016.

[16] Jana,M.,A.D.MacKerell,CHARMM Drude Polarizable Force Field for Aldopentofuranoses and Methylaldopentofuranosides.The Journal of Physical Chemistry B,2015.119(25): 7846-7859.

[17] Lamoureux,G.,A.D.M.Jr.,B.t.Roux,A simple polarizable model of water based on classical Drude oscillators.The Journal of Chemical Physics,2003.119(10): 5185-5197.

[18] Wang,J.,P.Cieplak,P.A.Kollman,How well does a restrained electrostatic potential(RESP)model perform in calculating conformational energies of organic and biological molecules? Journal of Computational Chemistry,2015.21(12): 1049-1074.

[19] Griffiths,D.J.,Introduction to Quantum Mechanics(2nd Edition).2nd ed.2004: Pearson Prentice Hall.

[20] Gordon,M.S.,D.G.Fedorov,S.R.Pruitt,L.V.Slipchenko,Fragmentation Methods: A Route to Accurate Calculations on Large Systems.Chemical Reviews,2012.112(1): 632-672.

[21] Elliott,P.,K.Burke,M.H.Cohen,A.Wasserman,Partition density-functional theory.Physical Review A,2010.82(2):024501.

[22] Song,G.-L.,Z.H.Li,Z.-P.Liu,X.-M.Cao,W.Wang,K.-N.Fan,Y.Xie,and H.F.Schaefer,Local Hybrid Divideand-Conquer Method for the Computation of Medium and Large Molecules.Journal of Chemical Theory and Computation,2008.4(12): 2049-2056.

[23] Babu,K.,S.R.Gadre,Ab initio quality one-electron properties of large molecules: Development and testing of molecular tailoring approach.Journal of Computational Chemistry,2003.24(4): 484-495.

[24] 梅曄,何曉,季長(zhǎng)鴿,張大為,and 張?jiān)鲚x,大分子體系的量子化學(xué)分塊方法(英文).化學(xué)進(jìn)展,2012(06): 1058-1064.

[25] Zhang,D.W.,Y.Xiang,J.Z.H.Zhang,New Advance in Computational Chemistry:? Full Quantum Mechanical ab Initio Computation of Streptavidin?Biotin Interaction Energy.The Journal of Physical Chemistry B,2003.107(44): 12039-12041.

[26] Zhang,D.W.,J.Z.H.Zhang,Molecular fractionation with conjugate caps for full quantum mechanical calculation of protein–molecule interaction energy.The Journal of Chemical Physics,2003.119(7): 3599-3605.

[27] Wang,X.,J.Liu,J.Z.H.Zhang,X.He,Electrostatically Embedded Generalized Molecular Fractionation with Conjugate Caps Method for Full Quantum Mechanical Calculation of Protein Energy.The Journal of Physical Chemistry A,2013.117(32): 7149–7161.

[28] Liu,J.,T.Zhu,X.Wang,X.He,and J.Z.H.Zhang,Quantum Fragment Basedab InitioMolecular Dynamics for Proteins.Journal of Chemical Theory and Computation,2015.11(12): 5897–5905.

[29] Jin,X.,J.Z.H.Zhang,X.He,Full QM Calculation of RNA Energy Using Electrostatically Embedded Generalized Molecular Fractionation with Conjugate Caps Method.The Journal of Physical Chemistry A,2017.121(12):2503–2514.

[30] Ji,C.,Y.Mei,J.Z.H.Zhang,Developing Polarized Protein-Speci fi c Charges for Protein Dynamics: MD} Free Energy Calculation of pKa Shifts for Asp26/Asp20 in Thioredoxin.Biophysical Journal,2008.95(3): 1080 - 1088.

[31] Wang,Y.J.,C.P.Sosa,A.Cembran,D.G.Truhlar,and J.L.Gao,Multilevel X-Pol: A Fragment-Based Method with Mixed Quantum Mechanical Representations of Different Fragments.Journal of Physical Chemistry B,2012.116(23): 6781-6788.

[32] Gao,J.L.,D.G.Truhlar,Y.J.Wang,M.J.M.Mazack,P.Loffler,M.R.Provorse,and P.Rehak,Explicit Polarization: A Quantum Mechanical Framework for Developing Next Generation Force Fields.Accounts of Chemical Research,2014.47(9): 2837-2845.

[33] Gao,J.,D.G.Truhlar,Y.Wang,M.J.M.Mazack,P.L?ffler,M.R.Provorse,and P.Rehak,Explicit Polarization: A Quantum Mechanical Framework for Developing Next Generation Force Fields.Accounts of Chemical Research,2014.47(9): 2837–2845.

[34] 吳相玉,陳守良,葛明德,陳閱增 普通生物學(xué).2010.

[35] Dennington,R.,T.Keith,J.Millam,GaussView Version 5.

[36] 黃熙泰,于自然,李翠鳳,現(xiàn)代生物化學(xué).2005: 化學(xué)工業(yè)出版社出版.

[37] Frisch,M.J.,G.W.Trucks,H.B.Schlegel,G.E.Scuseria,Gaussian 09.2009: Wallingford,CT.

[38] Mackerell,A.D.,M.Feig,C.L.Brooks,Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations.Journal of Computational Chemistry,2004.25(11): 1400-1415.

2016年10月7日

李 焱：中國(guó)科學(xué)院大連化學(xué)物理研究所,分子反應(yīng)動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，助理研究員，博士，主要研究方向?yàn)槔碚撚?jì)算方法的發(fā)展和算法優(yōu)化。

E-mail: liyan1982@dicp.ac.cn

湯天寶：遼寧師范大學(xué)，物理與電子技術(shù)學(xué)院，在讀碩士研究生，主要研究方向?yàn)槔碚撚?jì)算方法的發(fā)展。

E-mail: 962054297@qq.com

金 樂(lè)：遼寧師范大學(xué)，物理與電子技術(shù)學(xué)院，在讀碩士研究生，主要研究方向?yàn)槔碚撚?jì)算方法的發(fā)展。

E-mail: 870677753@qq.com

李冬婷：遼寧師范大學(xué)，物理與電子技術(shù)學(xué)院，在讀碩士研究生，主要研究方向?yàn)槔碚撚?jì)算方法的發(fā)展。

E-mail: 292850157@qq.com

李國(guó)輝：中國(guó)科學(xué)院大連化學(xué)物理研究所,分子反應(yīng)動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，研究員，博士，主要研究方向?yàn)槔碚撚?jì)算方法的發(fā)展及其在生命科學(xué)、生物醫(yī)藥和生物能源等領(lǐng)域中的應(yīng)用。

E-mail: ghli@dicp.ac.cn

許佩軍：遼寧師范大學(xué)，物理與電子技術(shù)學(xué)院，副教授，碩士，主要研究方向?yàn)槔碚撚?jì)算方法的發(fā)展。

E-mail: xpj631114@sohu.com

A Partition Iterative Quantum Method to Calculate High Precision Atom Charge

Li Yan1, Tang Tianbao2, Jin Le2, Li Dongting2,Li Guohui1＊, Xu Peijun2＊

1.State Key Laboratory of Molecular Reaction Dynamics,Dalian Institute of Chemical Physics,Chinese Academy of Sciences,Dalian,Liaoning 116023,China
2.School of Physics and Electronic Technology,Liaoning Normal University,Dalian,Liaoning 116023,China

The charge is one of the most important physicochemical parameters in studying the biomolecular.Based on the partition,the iterative Quantum Mechanics is introduced to calculate the high precision charge of atoms in multi-molecular system.And a calculating starategy named as Partition Iterative Quantum(PIQ)is proposed.Serval common biomolecular systems are applied as benchmarks to test the vailidation of PIQ.The charges calculated from QM are applied as the real values,and different initial patterns of charges are studied to compare the calculating ef fi ciency and presions.ThePIQis approved to be applied into different scales of biomoleculars charge calculations.

atom charge; quantum mechanics; force fi eld; partition method

10.11871/j.issn.1674-9480.2017.01.011

國(guó)家自然科學(xué)基金項(xiàng)目(21573217，91430110，31370714，21625302)

*通訊作者：李國(guó)輝(ghli@dicp.ac.cn)；許佩軍(xpj631114@sohu.com)