QM/MM結構優化

QM/MM幾何構型優化計算的python腳本如下：

import glob, math, os.path

from pBabel import AmberCrdFile_ToCoordinates3,

AmberTopologyFile_ToSystem ,

SystemGeometryTrajectory ,

AmberCrdFile_FromSystem ,

PDBFile_FromSystem ,

XYZFile_FromSystem

from pCore import Clone, logFile, Selection

from pMolecule import NBModelORCA, QCModelBDF, System

from pMoleculeScripts import ConjugateGradientMinimize_SystemGeometry,

FIREMinimize_SystemGeometry ,

LBFGSMinimize_SystemGeometry ,

SteepestDescentMinimize_SystemGeometry

# 定義結構優化接口

def opt_ConjugateGradientMinimize(molecule, selection):

molecule.DefineFixedAtoms(selection) #固定原子

#定義優化方法

ConjugateGradientMinimize_SystemGeometry(

molecule,

maximumIterations = 40, # 最大優化步數

rmsGradientTolerance = 0.1, #優化收斂控制

trajectories = [(trajectory, 1)]

) # 定義軌跡保存頻率

# 定義能量計算模式

nbModel = NBModelORCA()

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

# 讀取體系坐標和拓撲信息

molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")

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

# 關閉體系對稱性

molecule.DefineSymmetry(crystalClass = None) #QM/MM方法不支持使用周期性邊界條件

#. Define Atoms List

natoms = len(molecule.atoms) # 系統中總原子數

qm_list = range(72, 90) # QM 區原子

activate_list = range(126, 144) + range (144, 162) # MM區活性原子（優化中可以移動）

#定義MM區原子

mm_list = range (natoms)

for i in qm_list:

mm_list.remove(i) # MM 刪除QM原子

mm_inactivate_list = mm_list[:]

for i in activate_list :

mm_inactivate_list.remove(i)

# 輸入QM原子

qmmmtest_qc = Selection.FromIterable(qm_list)

# 定義各選擇區

selection_qm_mm_inactivate = Selection.FromIterable(qm_list + mm_inactivate_list)

selection_mm = Selection.FromIterable(mm_list)

selection_mm_inactivate = Selection.FromIterable(mm_inactivate_list)

# . Define the energy model.

molecule.DefineQCModel(qcModel, qcSelection = qmmmtest_qc)

molecule.DefineNBModel(nbModel)

molecule.Summary()

#計算優化開始時總能量

eStart = molecule.Energy()

#定義輸出文件目錄名

outlabel = 'opt_watbox_bdf'

if os.path.exists(outlabel):

pass

else:

os.mkdir (outlabel)

outlabel = outlabel + '/' + outlabel

# 定義輸出軌跡

trajectory = SystemGeometryTrajectory (outlabel + ".trj" , molecule, mode = "w")

# 開始第一階段優化

# 定義優化兩步

iterations = 2

# 順次固定QM區和MM區進行優化

for i in range(iterations):

opt_ConjugateGradientMinimize(molecule, selection_qm_mm_inactivate) #固定QM區優化

opt_ConjugateGradientMinimize(molecule, selection_mm) #固定MM區優化

# 開始第二階段優化

# QM區和MM區同時優化

opt_ConjugateGradientMinimize(molecule, selection_mm_inactivate)

#輸出優化后總能量

eStop = molecule.Energy()

#保存優化坐標，可以為xyz/crd/pdb等。

XYZFile_FromSystem(outlabel + ".xyz", molecule)

AmberCrdFile_FromSystem(outlabel + ".crd" , molecule)

PDBFile_FromSystem(outlabel + ".pdb" , molecule)

輸出體系收斂信息如下（此處僅展示前20步優化收斂結果）：

----------------------------------------------------------------------------------------------------------------

Iteration Function RMS Gradient Max. |Grad.| RMS Disp. Max. |Disp.|

----------------------------------------------------------------------------------------------------------------

0 I -1696839.69778731 2.46510318 9.94250232 0.00785674 0.03168860

2 L1s -1696839.82030342 1.38615730 5.83254788 0.00043873 0.00126431

4 L1s -1696839.90971371 1.41241184 5.29242524 0.00067556 0.00172485

6 L0s -1696840.01109863 1.41344485 4.70119338 0.00090773 0.00265969

8 L1s -1696840.09635696 1.44964059 5.72496661 0.00108731 0.00328490

10 L1s -1696840.17289698 1.28607709 4.73666387 0.00108469 0.00354577

12 L1s -1696840.23841524 1.03217304 3.00441004 0.00081945 0.00267931

14 L1s -1696840.30741088 1.40349698 5.22220965 0.00162080 0.00519590

16 L1s -1696840.43546466 1.32604042 4.51175225 0.00158796 0.00455431

18 L0s -1696840.52547251 1.27123125 4.20616166 0.00158796 0.00428040

20 L0s -1696840.60265453 1.08553355 3.12355616 0.00158796 0.00470223

----------------------------------------------------------------------------------------------------------------

輸出體系總能量信息如下：

注：QM/MM幾何構型優化一般不容易收斂，在實際操作中需要的技巧較多。常見的有，固定MM區，優化QM區；然后固定QM區優化MM區。如此往復循環幾次后，再同時優化QM區和MM區。優化是否收斂，和QM區的選擇及QM/MM邊界是否有帶電較多的原子等關系很大。為了加速優化，可以在計算時固定MM區，僅選擇離QM區較近的合適區域，作為活性區域，在優化中坐標可以變化。

QM/MM激發態計算

基于上一步的QM/MM幾何構型優化，繼而即可將MM區活性原子添加到QM區進行QM/MM-TDDFT計算，完整的代碼如下:

import glob, math, os.path

from pBabel import AmberCrdFile_ToCoordinates3,

AmberTopologyFile_ToSystem ,

SystemGeometryTrajectory ,

AmberCrdFile_FromSystem ,

PDBFile_FromSystem ,

XYZFile_FromSystem

from pCore import Clone, logFile, Selection

from pMolecule import NBModelORCA, QCModelBDF, System

from pMoleculeScripts import ConjugateGradientMinimize_SystemGeometry,

FIREMinimize_SystemGeometry ,

LBFGSMinimize_SystemGeometry ,

SteepestDescentMinimize_SystemGeometry

# 定義結構優化接口

def opt_ConjugateGradientMinimize(molecule, selection):

molecule.DefineFixedAtoms(selection) #固定原子

#定義優化方法

ConjugateGradientMinimize_SystemGeometry(

molecule,

maximumIterations = 40, # 最大優化步數

rmsGradientTolerance = 0.1, #優化收斂控制

trajectories = [(trajectory, 1)]

) # 定義軌跡保存頻率

# 定義能量計算模式

nbModel = NBModelORCA()

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

# 讀取體系坐標和拓撲信息

molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")

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

# 關閉體系對稱性

molecule.DefineSymmetry(crystalClass = None) #QM/MM方法不支持使用周期性邊界條件

#. Define Atoms List

natoms = len(molecule.atoms) # 系統中總原子數

qm_list = range(72, 90) # QM 區原子

activate_list = range(126, 144) + range (144, 162) # MM區活性原子（優化中可以移動）

#定義MM區原子

mm_list = range (natoms)

for i in qm_list:

mm_list.remove(i) # MM 刪除QM原子

mm_inactivate_list = mm_list[:]

for i in activate_list :

mm_inactivate_list.remove(i)

# 輸入QM原子

qmmmtest_qc = Selection.FromIterable(qm_list)

# 定義各選擇區

selection_qm_mm_inactivate = Selection.FromIterable(qm_list + mm_inactivate_list)

selection_mm = Selection.FromIterable(mm_list)

selection_mm_inactivate = Selection.FromIterable(mm_inactivate_list)

# . Define the energy model.

molecule.DefineQCModel(qcModel, qcSelection = qmmmtest_qc)

molecule.DefineNBModel(nbModel)

molecule.Summary()

#計算優化開始時總能量

eStart = molecule.Energy()

#定義輸出文件目錄名

outlabel = 'opt_watbox_bdf'

if os.path.exists(outlabel):

pass

else:

os.mkdir (outlabel)

outlabel = outlabel + '/' + outlabel

# 定義輸出軌跡

trajectory = SystemGeometryTrajectory (outlabel + ".trj" , molecule, mode = "w")

# 開始第一階段優化

# 定義優化兩步

iterations = 2

# 順次固定QM區和MM區進行優化

for i in range(iterations):

opt_ConjugateGradientMinimize(molecule, selection_qm_mm_inactivate) #固定QM區優化

opt_ConjugateGradientMinimize(molecule, selection_mm) #固定MM區優化

# 開始第二階段優化

# QM區和MM區同時優化

opt_ConjugateGradientMinimize(molecule, selection_mm_inactivate)

#輸出優化后總能量

eStop = molecule.Energy()

#保存優化坐標，可以為xyz/crd/pdb等。

XYZFile_FromSystem(outlabel + ".xyz", molecule)

AmberCrdFile_FromSystem(outlabel + ".crd" , molecule)

PDBFile_FromSystem(outlabel + ".pdb" , molecule)

# TDDFT計算

qcModel = QCModelBDF_template ( )

qcModel.UseTemplate (template = 'head_bdf_nosymm.inp' )

tdtest = Selection.FromIterable ( qm_list + activate_list )

# . Define the energy model.

molecule.DefineQCModel ( qcModel, qcSelection = tdtest )

molecule.DefineNBModel ( nbModel )

molecule.Summary ( )

# . Calculate

energy = molecule.Energy ( )

輸出體系總能量信息如下：

同時生成.log結果文件，和普通的激發態計算一樣，可以看到振子強度，激發能，激發態的總能量等信息

No. 1 w= 4.7116 eV -1937.8276358207 a.u. f= 0.0217 D= 0.0000 Ova= 0.6704

CV(0): A( 129 )-> A( 135 ) c_i: 0.7254 Per: 52.6% IPA: 7.721 eV Oai: 0.6606

CV(0): A( 129 )-> A( 138 ) c_i: 0.2292 Per: 5.3% IPA: 9.104 eV Oai: 0.8139

CV(0): A( 132 )-> A( 135 ) c_i: 0.4722 Per: 22.3% IPA: 7.562 eV Oai: 0.6924

CV(0): A( 132 )-> A( 138 ) c_i: -0.4062 Per: 16.5% IPA: 8.946 eV Oai: 0.6542

隨后還打印了躍遷偶極矩

*** Ground to excited state Transition electric dipole moments (Au) ***

State X Y Z Osc.

1 0.0959 0.1531 0.3937 0.0217 0.0217

2 0.0632 -0.1286 0.3984 0.0207 0.0207

3 -0.0797 -0.2409 0.4272 0.0287 0.0287

4 0.0384 -0.0172 -0.0189 0.0003 0.0003

5 1.1981 0.8618 -0.1305 0.2751 0.2751

吸收光譜分析

對于激發態我們往往需要得到理論預測的吸收譜峰形，也就是將每個激發態的吸收按一定的半峰寬進行高斯展寬。在TDDFT計算正常結束后，我們需要進入終端用命令調用BDF安裝路徑下的plotspec.py腳本執行計算。若用戶使用鴻之微云算力資源，進入命令端方式請查閱鴻之微云指南，此文不做贅述。

進入終端后，在目錄下運行$BDFHOME/sbin/plotspec.py bdf.out，會產生兩個文件，分別為bdf.stick.csv和bdf.spec.csv，前者包含所有激發態的吸收波長和摩爾消光系數，可以用來作棒狀圖，后者包含高斯展寬后的吸收譜（默認的展寬FWHM為0.5 eV），分別對比真空環境以及溶劑化效應下高斯展寬后的吸收譜情況，并用Excel、Origin等作圖軟件作圖如下：

上圖也可說明由于QM/MM計算考慮溶劑化效應，存在周圍其他分子的相互作用，從而使得吸收光譜發生紅移現象。