QM/MM幾何構(gòu)型優(yōu)化計(jì)算的python腳本
QM/MM結(jié)構(gòu)優(yōu)化
QM/MM幾何構(gòu)型優(yōu)化計(jì)算的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
# 定義結(jié)構(gòu)優(yōu)化接口
def opt_ConjugateGradientMinimize(molecule, selection):
molecule.DefineFixedAtoms(selection) #固定原子
#定義優(yōu)化方法
ConjugateGradientMinimize_SystemGeometry(
molecule,
maximumIterations = 40, # 最大優(yōu)化步數(shù)
rmsGradientTolerance = 0.1, #優(yōu)化收斂控制
trajectories = [(trajectory, 1)]
) # 定義軌跡保存頻率
# 定義能量計(jì)算模式
nbModel = NBModelORCA()
qcModel = QCModelBDF("GB3LYP:6-31g")
# 讀取體系坐標(biāo)和拓?fù)湫畔?/p>
molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")
molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")
# 關(guān)閉體系對(duì)稱性
molecule.DefineSymmetry(crystalClass = None) #QM/MM方法不支持使用周期性邊界條件
#. Define Atoms List
natoms = len(molecule.atoms) # 系統(tǒng)中總原子數(shù)
qm_list = range(72, 90) # QM 區(qū)原子
activate_list = range(126, 144) + range (144, 162) # MM區(qū)活性原子（優(yōu)化中可以移動(dòng)）
#定義MM區(qū)原子
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)
# 定義各選擇區(qū)
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()
#計(jì)算優(yōu)化開始時(shí)總能量
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")
# 開始第一階段優(yōu)化
# 定義優(yōu)化兩步
iterations = 2
# 順次固定QM區(qū)和MM區(qū)進(jìn)行優(yōu)化
for i in range(iterations):
opt_ConjugateGradientMinimize(molecule, selection_qm_mm_inactivate) #固定QM區(qū)優(yōu)化
opt_ConjugateGradientMinimize(molecule, selection_mm) #固定MM區(qū)優(yōu)化
# 開始第二階段優(yōu)化
# QM區(qū)和MM區(qū)同時(shí)優(yōu)化
opt_ConjugateGradientMinimize(molecule, selection_mm_inactivate)
#輸出優(yōu)化后總能量
eStop = molecule.Energy()
#保存優(yōu)化坐標(biāo)，可以為xyz/crd/pdb等。
XYZFile_FromSystem(outlabel + ".xyz", molecule)
AmberCrdFile_FromSystem(outlabel + ".crd" , molecule)
PDBFile_FromSystem(outlabel + ".pdb" , molecule)
輸出體系收斂信息如下（此處僅展示前20步優(yōu)化收斂結(jié)果）：
----------------------------------------------------------------------------------------------------------------
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幾何構(gòu)型優(yōu)化一般不容易收斂，在實(shí)際操作中需要的技巧較多。常見的有，固定MM區(qū)，優(yōu)化QM區(qū)；然后固定QM區(qū)優(yōu)化MM區(qū)。如此往復(fù)循環(huán)幾次后，再同時(shí)優(yōu)化QM區(qū)和MM區(qū)。優(yōu)化是否收斂，和QM區(qū)的選擇及QM/MM邊界是否有帶電較多的原子等關(guān)系很大。為了加速優(yōu)化，可以在計(jì)算時(shí)固定MM區(qū)，僅選擇離QM區(qū)較近的合適區(qū)域，作為活性區(qū)域，在優(yōu)化中坐標(biāo)可以變化。
QM/MM激發(fā)態(tài)計(jì)算
基于上一步的QM/MM幾何構(gòu)型優(yōu)化，繼而即可將MM區(qū)活性原子添加到QM區(qū)進(jìn)行QM/MM-TDDFT計(jì)算，完整的代碼如下: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
# 定義結(jié)構(gòu)優(yōu)化接口
def opt_ConjugateGradientMinimize(molecule, selection):
molecule.DefineFixedAtoms(selection) #固定原子
#定義優(yōu)化方法
ConjugateGradientMinimize_SystemGeometry(
molecule,
maximumIterations = 40, # 最大優(yōu)化步數(shù)
rmsGradientTolerance = 0.1, #優(yōu)化收斂控制
trajectories = [(trajectory, 1)]
) # 定義軌跡保存頻率
# 定義能量計(jì)算模式
nbModel = NBModelORCA()
qcModel = QCModelBDF("GB3LYP:6-31g")
# 讀取體系坐標(biāo)和拓?fù)湫畔?/p>
molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")
molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")
# 關(guān)閉體系對(duì)稱性
molecule.DefineSymmetry(crystalClass = None) #QM/MM方法不支持使用周期性邊界條件
#. Define Atoms List
natoms = len(molecule.atoms) # 系統(tǒng)中總原子數(shù)
qm_list = range(72, 90) # QM 區(qū)原子
activate_list = range(126, 144) + range (144, 162) # MM區(qū)活性原子（優(yōu)化中可以移動(dòng)）
#定義MM區(qū)原子
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)
# 定義各選擇區(qū)
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()
#計(jì)算優(yōu)化開始時(shí)總能量
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")
# 開始第一階段優(yōu)化
# 定義優(yōu)化兩步
iterations = 2
# 順次固定QM區(qū)和MM區(qū)進(jìn)行優(yōu)化
for i in range(iterations):
opt_ConjugateGradientMinimize(molecule, selection_qm_mm_inactivate) #固定QM區(qū)優(yōu)化
opt_ConjugateGradientMinimize(molecule, selection_mm) #固定MM區(qū)優(yōu)化
# 開始第二階段優(yōu)化
# QM區(qū)和MM區(qū)同時(shí)優(yōu)化
opt_ConjugateGradientMinimize(molecule, selection_mm_inactivate)
#輸出優(yōu)化后總能量
eStop = molecule.Energy()
#保存優(yōu)化坐標(biāo)，可以為xyz/crd/pdb等。
XYZFile_FromSystem(outlabel + ".xyz", molecule)
AmberCrdFile_FromSystem(outlabel + ".crd" , molecule)
PDBFile_FromSystem(outlabel + ".pdb" , molecule)
# TDDFT計(jì)算
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 ( )
輸出體系總能量信息如下：
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
吸收光譜分析
對(duì)于激發(fā)態(tài)我們往往需要得到理論預(yù)測(cè)的吸收譜峰形，也就是將每個(gè)激發(fā)態(tài)的吸收按一定的半峰寬進(jìn)行高斯展寬。在TDDFT計(jì)算正常結(jié)束后，我們需要進(jìn)入終端用命令調(diào)用BDF安裝路徑下的plotspec.py腳本執(zhí)行計(jì)算。若用戶使用鴻之微云算力資源，進(jìn)入命令端方式請(qǐng)查閱鴻之微云指南，此文不做贅述。
進(jìn)入終端后，在目錄下運(yùn)行$BDFHOME/sbin/plotspec.py bdf.out，會(huì)產(chǎn)生兩個(gè)文件，分別為bdf.stick.csv和bdf.spec.csv，前者包含所有激發(fā)態(tài)的吸收波長(zhǎng)和摩爾消光系數(shù)，可以用來(lái)作棒狀圖，后者包含高斯展寬后的吸收譜（默認(rèn)的展寬FWHM為0.5 eV），分別對(duì)比真空環(huán)境以及溶劑化效應(yīng)下高斯展寬后的吸收譜情況，并用Excel、Origin等作圖軟件作圖如下：
-
光譜分析
+關(guān)注
關(guān)注
0文章
41瀏覽量
10380 -
python
+關(guān)注
關(guān)注
57文章
4876瀏覽量
90022
原文標(biāo)題：鴻之微BDF軟件計(jì)算賞析｜采用BDF的QM/MM多尺度計(jì)算方法研究晶體的光物理性質(zhì)（二）
文章出處：【微信號(hào)：hzwtech，微信公眾號(hào)：鴻之微】歡迎添加關(guān)注！文章轉(zhuǎn)載請(qǐng)注明出處。
發(fā)布評(píng)論請(qǐng)先 登錄
Simlab+CAD+Hyperstudy聯(lián)合優(yōu)化幾何參數(shù)
官方新品 | 虹科PCAN-Explorer 7發(fā)布：帶來(lái)Python腳本與靈活授權(quán)新體驗(yàn)
光子精密QM系列閃測(cè)儀優(yōu)化VR眼鏡關(guān)鍵制程質(zhì)量控制
用 Python 給 Amazon 做“全身 CT”——可量產(chǎn)、可擴(kuò)展的商品詳情爬蟲實(shí)戰(zhàn)
Pico Technology發(fā)布Python軟件包pyPicoSDK
termux調(diào)試python猜數(shù)字游戲
termux如何搭建python游戲
Python腳本實(shí)現(xiàn)運(yùn)維工作自動(dòng)化案例
怎么導(dǎo)出python邊緣計(jì)算中的APP，想進(jìn)行修改又找不到源碼？
FRED應(yīng)用：MTF的計(jì)算
請(qǐng)問如何在imx8qm上將PCIe和SMMU用于顯卡？
零基礎(chǔ)入門：如何在樹莓派上編寫和運(yùn)行Python程序？
【MWC 2025】羅德與施瓦茨推出CMX500 AI腳本助手，AI驅(qū)動(dòng)自動(dòng)化助力移動(dòng)設(shè)備測(cè)試
工商網(wǎng)監(jiān)
評(píng)論