本專題將介紹一種量子化學(xué)與分子力學(xué)結(jié)合的方法（QM/MM方法），該方法既包括量子化學(xué)的精確性，又利用分子力學(xué)的高效性，其基本思想是用量子力學(xué)處理感興趣的區(qū)域，其余部分用經(jīng)典分子力學(xué)來處理。本章節(jié)以一個典型的沒食子酸分子（Gallic Acid，GA）為例，介紹輸入文件準(zhǔn)備，QM/MM單點計算，QM/MM結(jié)構(gòu)優(yōu)化，QM/MM激發(fā)態(tài)計算。其中，BDF程序主要完成量子化學(xué)計算部分，其余部分由BDF開發(fā)成員修改的pDynamo2.0程序包完成。同時介紹如何讀取數(shù)據(jù)用于結(jié)果分析，幫助用戶深入了解BDF軟件的使用。

輸入文件準(zhǔn)備

一般來說，QM/MM計算之前，需要對目標(biāo)體系進(jìn)行分子動力學(xué)模擬，得到適合的初始構(gòu)象。當(dāng)采用PDB、MOL2或xyz文件作為輸入時，pDynamo2.0程序包僅支持OPLS力場，對于小分子和非標(biāo)準(zhǔn)氨基酸力場參數(shù)不全，不推薦使用。建議優(yōu)先采用Amber程序，通過拓?fù)湮募斎肓鰠?shù)。以Amber為例，從動力學(xué)模擬軌跡提取感興趣的結(jié)構(gòu)存儲于.crd文件中，與對應(yīng)的參數(shù)/拓?fù)湮募?prmtop一起可以作為QM/MM計算的起始點。Python腳本如下：

其中，需要提前安裝好AmberTools，python2.0版本，并正確設(shè)置好AMBERHOME和PDYNAMO環(huán)境變量，關(guān)于如何將GallicAcid.pdb初始結(jié)構(gòu)文件(圖1，晶胞為2*1*1)生成使用AmberTools21程序相對應(yīng)的坐標(biāo)文件GallicAcid.crd和參數(shù)/拓?fù)湮募礼allicAcid.prmop的方法如下：

運行antechamber程序?qū)db文件轉(zhuǎn)化為mol2文件：antechamber -i GallicAcid.pdb -fi pdb -o GallicAcid.mol2 -fo mol2 -j 5 -at amber -dr no

-i 指定輸入文件

-fi 指定輸入文件類型

-o 指定輸出文件

-fo 指定輸出文件類型

-j匹配原子類型和鍵類型

-at定義原子類型

運行parmchk2程序生成對應(yīng)體系的力場參數(shù)文件：parmchk2 -i GallicAcid.mol2 -f mol2 -o GallicAcid.frcmod

運行tleap程序構(gòu)建系統(tǒng)拓?fù)洳榉肿佣x力場參數(shù)步驟如下：

1.使用tleap命令啟動tleap程序

2. 確定并加載體系力場：source leaprc.gaff(此為GAFF力場)

3. 調(diào)入配體mol2文件：GA = loadmol2 GallicAcid.mol2

4. 檢查導(dǎo)入的結(jié)構(gòu)是否準(zhǔn)確或缺失參數(shù)：check GA

5. 調(diào)入體系分子的模板，并補(bǔ)全庫文件中缺失的參數(shù)：loadamberparams GallicAcid.frcmod

6. 準(zhǔn)備生成的Sustiva庫文件：saveoff GA GallicAcid.lib

7. 修改生成的Sustiva庫文件并調(diào)入該文件：loadoff GallicAcid.lib

8. 保存.crd和.prmop文件：saveamberparm GA GallicAcid.prmtop GallicAcid.crd

9. 退出tleap程序：quit

分子動力學(xué)模擬

1.此處采用amber軟件進(jìn)行分子動力學(xué)模擬，首先對體系進(jìn)行能量最小化模擬，輸入文件min.in如下：

Initial minimisation of GallicAcid complex

&cntrl

imin=1, maxcyc=200, ncyc=50,

cut=16, ntb=0, igb=1,

&end

imin=1：運行能量最小化

maxcyc=200：能量最小化的最大循環(huán)數(shù)

ncyc=50：最初的0到ncyc循環(huán)使用最速下降算法, 此后的ncyc到maxcyc循環(huán)切換到共軛梯度算法

cut=16：以埃為單位的非鍵截斷距離

ntb=0：關(guān)閉周期性邊界條件

igb=1：Born模型

使用如下命令運行能量最小化:

sander -O -i min.in -o GallicAcid_min.out -p GallicAcid.prmtop -c GallicAcid.crd -r GallicAcid_min.rst &

其中GallicAcid_min.rst為輸出包含坐標(biāo)和速度的重啟文件

2.接下來利用最小化模擬得到的重啟文件升溫系統(tǒng)，從而完成分子動力學(xué)模擬，輸入文件md.in如下：

Initial MD equilibration

&cntrl

imin=0, irest=0,

nstlim=1000,dt=0.001, ntc=1,

ntpr=20, ntwx=20,

cut=16, ntb=0, igb=1,

ntt=3, gamma_ln=1.0,

tempi=0.0, temp0=300.0,

&end

imin=0：進(jìn)行分子動力學(xué)(MD)

irest=0：讀取先前保存的重新啟動文件讀取坐標(biāo)和速度

nstlim=1000：運行的MD步數(shù)

dt=0.001：時間步長（單位：ps）

ntc=1：不啟用SHAKE約束

ntpr=20：每ntpr步輸出能量信息mdout一次

ntwx=20：每ntwx步輸出Amber軌跡文件mdcrd一次

ntt=3：Langevin恒溫器控制溫度

gamma_ln=1.0：Langevin恒溫器的碰撞頻率

tempi=0.0：模擬的初始溫度

temp0=300.0：模擬的最終溫度

使用如下命令運行分子動力學(xué)模擬:

sander -O -i md.in -o md.out -p GallicAcid.prmtop -c GallicAcid_min.rst -r GallicAcid_md.rst -x GallicAcid_md.mdcrd -inf GallicAcid_md.mdinfo

其中GallicAcid_md.mdcrd文件即為MD模擬的軌跡文件，可借助VMD軟件進(jìn)行可視化顯示分子結(jié)構(gòu)，并從動力學(xué)模擬軌跡提取感興趣的結(jié)構(gòu)存儲于.crd文件中。

QM/MM總能量計算

分子動力學(xué)模擬后提取文件為GallicAcid.prmtop， GallicAcid.crd，可對體系進(jìn)行全量子化學(xué)總能量計算，python代碼如下：

import glob, math, os

from pBabel import AmberCrdFile_ToCoordinates3, AmberTopologyFile_ToSystem

from pCore import logFile

from pMolecule import QCModelBDF, System

# 讀取水盒子坐標(biāo)和拓?fù)湫畔?/p>

molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")

molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")

# 定義能量計算模式，此處為全體系密度泛函計算，可以定義方法和基組，分別為GB3LYP和6-31g，

model = QCModelBDF("GB3LYP:6-31g")

molecule.DefineQCModel(model)

molecule.Summary() #輸出體系計算設(shè)置信息

# 計算總能量

energy = molecule.Energy()

除了可以用全量子化學(xué)QM計算體系總能量，也可對感興趣的分子進(jìn)行QM/MM計算（本例為指定第五個分子用QM方法計算），QM/MM組合能量計算python腳本如下：

import glob, math, os

from pBabel import AmberCrdFile_ToCoordinates3, AmberTopologyFile_ToSystem

from pCore import logFile, Selection

from pMolecule import NBModelORCA, QCModelBDF, System

# 定義能量計算模式

nbModel = NBModelORCA() #處理QM和MM區(qū)相互作用

qcModel = QCModelBDF("GB3LYP:6-31g")

# 讀取體系坐標(biāo)和拓?fù)湫畔?/p>

molecule = AmberTopologyFile_ToSystem("GallicAcid.prmtop")

molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")

# 關(guān)閉體系對稱性

molecule.DefineSymmetry(crystalClass = None) #QM/MM方法不支持使用周期性邊界條件，故關(guān)閉周期性邊界條件

# 指定QM區(qū)

qm_area = Selection.FromIterable(range (72, 90)) # 指定第五個分子用QM方法計算，其中(72, 90)指明原子列表索引值為72，73，74…..87,88,89，該值=原子序數(shù)-1

# 定義能量計算模式

molecule.DefineQCModel (qcModel, qcSelection = qm_area)

molecule.DefineNBModel (nbModel)

molecule.Summary()

# 計算總能量

energy = molecule.Energy()

QM/MM模擬的輸出總結(jié)了MM部分，QM部分，QM區(qū)和MM區(qū)相互作用部分的計算細(xì)節(jié)如下：

輸出體系總能量信息以及各部分能量貢獻(xiàn)如下：