粗?；到y(tǒng)的GROMACS模擬

馮 劍

1 前言

計算機模擬是繼實驗和理論研究方法后的又一種科學研究手段受到極大的關(guān)注。分子動力學模擬作為重要的工具之一也被廣大的化學化工科技工作者所采用?？萍脊ぷ髡叱Ｗ孕芯幹栖浖鉀Q自己的研究問題[1]。隨著科學研究深入和計算技術(shù)的發(fā)展，個人編制完整的計算程序用于各類問題求解變得越來越困難。GROMACS[2]、NAMD、LAMMPS等模擬軟件已被廣泛使用?？萍脊ぷ髡哂袝r常關(guān)注研究問題主要影響因素而拋開原子具體細節(jié)，理論工作者尤為如此，因而也非常關(guān)注粗粒化的模擬工作[3, 4]。

GROMACS特別適用于生物大分子體系的分子動力學模擬，它提供了許多可選的力場，以及很多計算技術(shù)，容易擴展到粗?；M領(lǐng)域。知名的MARTINI粗?；鯷5]就采用了對實際的生物分子中多個原子使用GROMACS的虛擬格點技術(shù)進行粗?；?。使用粗?；Ｐ偷腉ROMACS系統(tǒng)性工作很少報道，本文較完整給出GROMACS粗?；M的一般方法，軟件版本是GROMACS v2018.2。

2 力場擴展和體系構(gòu)建

2.1 擴展OPLS-AA力場

分子動力學模擬必須要選擇具體的力場，參數(shù)化研究體系。利用GROMACS進行粗?；幚硇枰獢U展力場的處理能力。為此，需要熟悉相關(guān)力場以及GROMACS中力場的組織。本文以OPLS-AA力場為例，對此進行粗?；幚?，其它力場方向相同。

OPLS-AA力場包含多項非鍵和成鍵作用，這里僅就其中的非鍵范德華作用進行說明，更詳細的內(nèi)容可參考相關(guān)文獻[6]和GROMACS手冊。GROMACS處理的非鍵范德華作用有Lennard-Jones(LJ)和Buckingham兩種勢能形式，OPLS-AA中使用的是LJ勢能函數(shù)，其原始形式為

在模擬軟件中考慮計算效率，常采用如下變形式

OPLS-AA力場采用后者。不同原子間的相互作用參數(shù)根據(jù)純組分值按以下組合規(guī)則獲得

該組合規(guī)則在GROMACS提供的三種規(guī)則之一，編號為3。在GROMACS提供的OPLS-AA力場中LJ輸入值依然是純組分的ii和ii。OPLS-AA力場的非鍵作用規(guī)則在力場包含文件forcefield.itp中給與表述，該文件有如下內(nèi)容：

[ defaults ]

1 3 yes 0.5 0.5

第一列1表示非鍵范德華作用使用LJ勢能，第二列3表示使用編號3的組合規(guī)則，yes表示為處于化學鍵上1-4位原子產(chǎn)生非鍵作用，后面兩列分別為LJ和靜電1-4作用的乘積因子。

為了避免直接修改GROMACS庫中的原力場造成力場污染，最好的方法是將該力場復制到工作目錄。另外，為了與現(xiàn)有的粗?；鲞M行比較，本文使用MARTINI力場的相關(guān)參數(shù)。修改力場主要有以下幾步。

(1)擴展原子類型。打開atomtypes.atp添加幾個粗?；Ｗ樱Q任意但不要與已有的名稱沖突。作為示范以下僅添加2個，摩爾質(zhì)量為72(四個水分子的粗粒化粒子)

CGA 72.00

CGB 72.00

(2)為新的粒子類型添加相互作用。對于非鍵作用，打開ffnonbonded.itp文件，在最后追加以下幾行數(shù)據(jù)

CGA CGA 72.00 0.0000 A 0 0

CGB CGB 72.00 0.0000 A 0 0

[nonbond_params ]

CGA CGA 1 0.47 5.0

CGB CGB 1 0.47 4.5

CGA CGB 1 0.47 4.0

前面二行分別給定新增兩種粒子類型的具體參數(shù)，第一列是粒子名稱，最后兩列是LJ函數(shù)的兩個參數(shù)，這里設置為0。這是因為它們間的所有作用不是由系統(tǒng)按照組合規(guī)則生成，而是在[ nonbond_params ]語句塊中指定。可以給定這兩個參數(shù)具體的值。如果希望按照組合規(guī)則計算交叉項，則不需要[ nonbond_params ]段。

(3)添加鍵長作用。打開ffbonded.itp文件，在最后追加以下幾行數(shù)據(jù)，其中第三列1表示諧振勢，后面兩個參數(shù)分別是平衡鍵長和力常數(shù)。

[bondtypes ]

CGA CGA 1 0.47 1250

CGA CGB 1 0.47 1250

CGB CGB 1 0.47 1250

(4)擴展殘基數(shù)據(jù)庫。GROMACS的一個強大的功能是利用pdb2gmx模塊為生物分子產(chǎn)生拓撲文件。要想使用pdb2gmx為自定義的模型分子產(chǎn)生拓撲就需要擴展殘基數(shù)據(jù)庫。新建rtp文件或打開aminoacids.rtp，在最后添加

[ CGA ]

[ atoms ]

CGA CGA 0.000 0

[ CGB ]

[ atoms ]

CGB CGB 0.000 0

如果為多個粒子構(gòu)成的鏈狀分子生成拓撲，則[atoms]還需要列出分子上所有粒子，同時給出[ bonds ]、[ angles ]和[ dihedrals ]等成鍵作用項。由于不同的分子均要如此處理，最直接的辦法是為每種分子構(gòu)建一個拓撲包含文件，使用時將該文件包含在系統(tǒng)拓撲文件中。這里給出處理由這兩者粒子組成鏈長為2的分子(殘基)，成鍵作用僅有一對鍵長作用。

[ CGM ]

[ atoms ]

CGA CGA 0.000 0

CGB CGB 0.000 0

[ bonds ]

CGA CGB

2.2 生成系統(tǒng)拓撲和構(gòu)建分子拓撲

為具體分子構(gòu)建拓撲首先需要有該分子的結(jié)構(gòu)文件。在GROMACS中常用的有兩種格式，一種是pdb格式文件，即蛋白質(zhì)數(shù)據(jù)庫文件[7]，生物大分子的結(jié)構(gòu)一般由該格式提供。還有一種是GROMACS默認的gro文件。對于粗?；肿?，可以使用手工或使用Chem3D等3D分子繪圖軟件制作。對于單原子分子，使用作圖軟件插入一個單原子，保存為pdb文件(命名為atom.pdb)并進行適當修改；或使用一個文本編輯器，直接輸入相關(guān)內(nèi)容制作atom.pdb文件。該文件僅需2行，如：

ATOM 1 CGA CGA 1 0.000 0.000 0.000

END

注意pdb文件有固定格式，其中7-11列是原子編號，13-16列是原子名稱，18-20列是殘基名稱(不是必須的)。本文還處理雙原子分子(命名為mole.pdb)，只需要添加一行并做一些改變，內(nèi)容如：

ATOM 1 CGA CGM 1 0.000 0.000 0.000

ATOM 2 CGB CGM 1 0.000 0.000 0.470

END

這里做了幾個改變，首先對分子上的粒子從1進行編號。將分子(殘基)命名為CGM，與殘基數(shù)據(jù)庫中的[ CGM ]塊對應。對第二個粒子z坐標進行修改，距第一個粒子0.47nm。

力場修改完成后就可以為這些分子生成拓撲，執(zhí)行如下命令即可

gmx pdb2gmx -f atom.pdb -o atom.gro -p topol.top -ff oplsaa -water none

上面命令生成的系統(tǒng)拓撲文件topol.top內(nèi)容非常簡單，其核心內(nèi)容是

#include "./oplsaa.ff/forcefield.itp"

[moleculetype ]

Other 3

[ atoms ]

1 CGA 1 CGA CGA 1 0 72

[ system ]

Protein

[ molecules ]

Other 1

該系統(tǒng)拓撲文件不僅包含了模擬系統(tǒng)的分子拓撲(前面[ moleculetype ]、[ atoms ]等部分)也包含給系統(tǒng)命名和組成系統(tǒng)的各類分子及其數(shù)目等信息(后面[ system]和[molecules]部分)。在GROMACS模擬中，最好為自行制作的分子建立單獨的拓撲文件，在使用時添加#include語句包含在系統(tǒng)中。即將最后的[ system]和[molecules]段刪除，并刪除第一行力場包含文件，將這些內(nèi)容另存為mol_a.itp，它就是A分子的分子拓撲文件。

[moleculetype ]語句塊給定了自行制作分子的名稱，該分子系統(tǒng)不能識別，因而給了“Other”名稱，它可以修改為任何有意義的名稱。第二列3表示不計算連續(xù)三個鍵上的粒子之間的非鍵作用，粗粒化模擬中該值一般設為1。單粒子分子沒有化學鍵，該項不起作用。修改后為

[moleculetype ]

Mol_A 3

后面的[ molecules ]也要做相應的修改，即

[ molecules ]

Mol_A 1

該行的第二列數(shù)據(jù)是系統(tǒng)中Mol_A分子的個數(shù)，后期需要根據(jù)系統(tǒng)的實際情況進行修改。

前面產(chǎn)生的atom.gro結(jié)構(gòu)文件也是具有固定格式的。該文件包含多行內(nèi)容，其中第一行是任意的說明字符串，第二行是總原子數(shù)，最后一行是模擬盒子邊長，單位是nm，這些是自由格式。其它各行是原子信息，是固定格式的。這里只有一行，內(nèi)容如下

1CGA CGA 1 0.000 0.000 0.000

第一列是殘基的編號和名稱(二項)，第二列是原子名稱，第三列是原子編號，每項均為5個字符寬度。后面三列是原子坐標，8個字符寬度其中小數(shù)位寬度為3。

再次為雙粒子分子生成拓撲文件topol2.top：

gmx pdb2gmx -f mole.pdb -o mole.gro -p topol2.top -ff oplsaa -water none

其中該分子拓撲內(nèi)容如下

[moleculetype ]

Other 3

[ atoms ]

1 CGA 1 CGM CGA 1 0 72

2 CGB 1 CGM CGB 1 0 72

[ bonds ]

1 2 1

同樣可以將這部分內(nèi)容保存為單獨的拓撲文件，如moleb.itp，在使用時包含它。在理論研究中的粗?；Ｐ土黧w一般都不復雜，可直接根據(jù)分子結(jié)構(gòu)依次給出分子中所有粒子([ atoms ])，所有鍵長作用([ bonds ])和鍵角作用[ angles ]等。

2.3 構(gòu)建模擬體系

構(gòu)建包括不同類型分子的體系是比較復雜的工作，可以借助工具如packmol或自行編程解決。GROMACS也提供了一些模塊構(gòu)建簡單系統(tǒng)，如insert-molecules和genconf等。

使用insert-molecules在邊長為4.2nm的立方型盒子中插入1000個A粒子，粒子半徑為0.47nm。

gmx insert-molecules -ci atom.pdb -o atom_box.gro -box 4.2 4.2 4.2 -nmol 1000 -radius 0.47

再將topol.top中[ molecules ]語句塊Mol_A的分子數(shù)從1修改為1000。這樣就獲得模擬體系的結(jié)構(gòu)文件和拓撲文件。

對于雙原子分子，構(gòu)建500個分子的體系可以類似處理，這里給出genconf方法。首先將分子放置一個長方體盒子中并位于盒子中心。

gmx editconf -f mole.pdb -o mole_c.gro -box 0.42 0.42 0.84 -c

將這個盒子在xy維堆積10倍，z維堆積5倍，產(chǎn)生500個分子系統(tǒng)，即

gmx genconf -f mole_c.gro -o mole_box.gro -nbox 10 10 5 -rot

將topol2.top中[ molecules ]語句塊Mol_B(原Other修改為Mol_B)的分子數(shù)從1修改為500即可。注意對于較長分子鏈長使用insert-molecules未必能產(chǎn)生足夠高的濃度，而genconf則會產(chǎn)生更加有序的系統(tǒng)。

有了系統(tǒng)拓撲和結(jié)構(gòu)文件，就可以按照常規(guī)的GROMACS模擬步驟執(zhí)行模擬，一般有能量最小化、NVT(約束)平衡、NPT(約束)平衡和MD模擬等步驟。本文研究分子不帶電，因而在上面各步提供的模擬參數(shù)文件(.mdp)中不需要設置靜電作用項，rcoulomb和rvdw設置為1.2，tc-grps=System、tcoupl=V-rescale以及pcoup=Parrinello-Rahman等。

3 結(jié)果與討論

圖1是使用粗粒化的MARTINI水分子力場參數(shù)進行模擬獲得298.15K、308.15K至338.15K等五個溫度下水的密度，并與美國地質(zhì)勘探局水科學學院(https://water.usgs.gov/edu/density.html)的水密度數(shù)據(jù)進行了比較。模擬結(jié)果普遍比實驗結(jié)果大。該力場的創(chuàng)立者也發(fā)現(xiàn)這個問題，建議在水中增加少量防凍結(jié)粒子，其直徑增大到0.57nm[5]。

圖1 不同溫度下水密度的模擬值與實驗值比較

圖2 不同溫度下水擴散系數(shù)的模擬值與實驗值比較

MARTINI粗?；Ｐ瞳@得的密度與實驗值有10%左右的誤差，而擴散系數(shù)差別更大。圖2是相應溫度下模擬的水擴散系數(shù)與實驗值的比較[8]。模擬值普遍比實驗數(shù)據(jù)小，除了在318.15K值與實驗值數(shù)量級一致外，其它溫度下有1-2個數(shù)量級差別。造成這種現(xiàn)象，一方面是實驗水分子和粗?；乃肿哟笮『唾|(zhì)量不同，需要對模擬值進行標度；另一方面是這個粗?；乃Ｐ蜔o法精確描述擴散系數(shù)之類物理量。為實際體系構(gòu)筑粗粒化模型本身就是一個艱巨的難題，這些粗?；哪Ｐ屯荒芙毯妹枋錾贁?shù)性質(zhì)，一般無法準確描述眾多性質(zhì)。

圖3是兩個粒子構(gòu)成的分子體系各粒子間徑向分布函數(shù)的模擬。在該模擬中所有粒子間的作用勢均設為5.0。相同粒子對CGA-CGA和CGB-CGB間的徑向分布函數(shù)值在全距離范圍內(nèi)幾乎一致。不同粒子對CGA-CGB與它們之間的值大約在小于0.49nm時才表現(xiàn)處明顯的差異。從0.472nm始徑向分布函數(shù)值大于1，這與CGA間CGB存在化學鍵，鍵長為0.47nm是一致的。由該圖可以得出，盡管在構(gòu)建該分子體系是有較大的有序性。但經(jīng)過MD模擬，這種有序被打破，獲得了應有的無序狀態(tài)。

圖3 粒子間的徑向分布函數(shù)

任何MD模擬都有自己的單位集，在GROMACS中長度r單位為1nm，質(zhì)量單位為原子質(zhì)量(C質(zhì)量的1/12)，即1u，時間為ps，能量為kJ/mol，力為kJ/mol/nm等。而在粗粒化模擬中質(zhì)量為m，長度為，時間為，能量為，力為/等。這些量直接設為1，即所謂的約化單位或?qū)Ρ葐挝?。由于導出單位之間彼此關(guān)聯(lián)，一旦選定基礎單位后，其它單位就不能再隨意變動。因此，在使用GROMACS進行粗?；M時還需要考慮約化單位與GROMACS單位間的對應關(guān)系。如將對應于GROMACS中的0.47nm，則由GROMACS得出的長度值只要除以0.47nm，則得到粗?；到y(tǒng)的長度量。換句話說，將粗?；瘑挝环謩e對應于GROMACS中單位的具體值，將GROMACS中獲得的相應物理量的值除以其粗?；瘑挝换卮膶导纯?。如上面318.15K模擬的擴散系數(shù)為1.0310-9m2/s，擴散系數(shù)單位是2/。對應關(guān)系—0.47nm、—1ps，則粗?；M的約化擴散系數(shù)為4.6610-3。該溫度下密度為1051.07kg/m3，粗?；瘑挝籱/3，對應關(guān)系為—0.47nm、m—72u，則約化密度為0.91。

4 結(jié)論

本文給出了通過GROMACS擴展系統(tǒng)力場，強化軟件能力，構(gòu)建粗?；w系的方法，并以現(xiàn)有的MARTINI力場參數(shù)為例，對兩個分子體系進行模擬，獲得了滿意的結(jié)果。同時也介紹了GROMACS體系與用戶選擇的粗粒化系統(tǒng)間的映射關(guān)系。本文工作表明功能強大的GROMACS分子動力學軟件具有較大的靈活性和較強的擴展性，可以用于研究更廣泛的系統(tǒng)和更復雜的問題。