计算化学公社

 找回密码 Forget password
 注册 Register
Views: 15711|回复 Reply: 5

[辅助/分析程序] Growing String Method (GSM) - pyGSM 安裝、簡介與範例

[复制链接 Copy URL]

31

帖子

1

威望

366

eV
积分
417

Level 3 能力者

发表于 Post on 2021-5-6 21:25:14 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 ms860309 于 2021-5-6 21:33 编辑

#####
#簡介#
#####
傳統上對於瞭解一個化學反應,非常重要的一點是尋找到Transition State(TS),而TS通常難以收斂又或者需要多加嘗試不同的排列組合才能得到一個TS。
一般來講尋找TS最重要的是給一個足夠好的initial guess,常見的方法有對2D PES做掃描(如Angew.Chem.Int.Ed.2021,60,4266–4274)又或者透過NEB(J. Chem. Phys. 113, 9901 (2000))
詳情可以參考:过渡态、反应路径的计算方法及相关问题 - 思想家公社的门口:量子化学·分子模拟·二次元 (sobereva.com)
各些方法其實各有優缺點,例如說NEB彈簧間彼此是等距的,因此常常會在高能量區域sample太多的點,通常高能量的結構難以描述,導致容易發生難以收斂或者高計算資源的狀況。
以下我會介紹String method的安裝方法,以及給一個範例,至於原理請各位參考Paul Zimmerman的論文:
J. Chem. Phys. 120, 7877 (2004)
J. Chem. Theory Comput. 2013, 9, 7, 3043–3050
J. Chem. Phys. 138, 184102 (2013)
Phys. Chem. Chem. Phys., 2018, 20, 27394


我在這裡給一個較為形象化的比喻,你假設目前身處在一個低谷(local minimum)中,我告訴你一個方向(driving coordinate),你朝著那個方向看過去發現了一個山峰(TS guess)。
這時候你就開始往那個方向爬,每隔一段距離(step size)重新抬頭找山峰,當你發現你已經身處最高點時(即環顧四周都是下坡),那麼你就已經在一個TS附近。(你還是需要OptTS)
目前GSM還會結合Climb image來優化TS initial 以致於得到一個更好的initial guess

#####
#流程#
#####

1. 往driving coordinate方向插值走一步
2. 計算gradient並且優化結構,直到達到收斂條件
3. 長node
4. 重複2-3直到string中包含ts
5. optimize string
6. climb image
7. 當TS node 達到收斂條件時結束

driving coordinate對於de-gsm而言是reactant與product的插值,但是對於SE-GSM則是isomer來作為driving coordinate
如果只有1-4則稱為Freezing string method(FSM),通常GSM會得到較優的reaction path(其實沒用因為後續還是要跑TS+IRC)與TS,耗時較長
而FSM則是少去優化整條string的過程,但是其實也能給到不差的TS guess
至於哪一種比較好,其實並沒有一定,沒有任何一種string method能找到所有反應 (J Am Chem Soc. 2018 Jan 24;140(3):1035-1048.)

#####
#安裝#
#####

PART I:安裝pyGSM
1. 利用conda創建一個python3的環境,例如gsm_env
2. conda activate gsm_env
3. git clone https://github.com/ZimmermanGroup/pyGSM.git
4. cd pyGSM
5. python setup install
PART II:安裝計算軟體
目前pyGSM支援大部分主流的量子力學計算軟體列表如下
如果你需要Gaussian可以使用C++版本的molecularGSM
["QChem","Orca","Molpro","PyTC","TeraChemCloud","OpenMM","DFTB","TeraChem","BAGEL","xTB_lot"]

#####
#範例#
#####

PART I: Double-Ended Growing String Method (DE-GSM)
一、選定你的Reactant以及對應的Product
二、優化Reactant與Product結構,通常會與GSM的level of theory相同,以下用xTB舉例
三、準備好你的job file
此時資料夾中應該有,並且對應內容為下圖
202105062047563717..png 202105062048378566..png 202105062048553789..png
四、執行
五、結果
計算後資料夾中的圖片檔為Energy profile,TSNode.xyz為TS initial guess後續再以其他軟體優化TS
至於優化過程可以參考我附檔的scratch/opt_converged_000_000.xyz再丟到VMD看即可(記得去掉開頭的molden format),這邊我就不做成GIF節省論壇圖片資源
0000_string.png

PART II: Single-Ended Growing String Method (SE-GSM)
與上個例子相同,差別在於不需要給product結構而是給driving coordinate,可以給ADD BREAK bond編號、angle、torsion等資訊,詳情請自行研究
202105062101002113..png 202105062101464418..png 202105062102019559..png 202105062102251395..png

#####
#補充#
#####

1.有關於coordinate_type的選擇預設為TRIC(J. Chem, Phys. 144, 214108),你也可以選擇DLC,通常情況下差異不大,但是如果你的系統是金屬cluster又或者你的體系常常有鍵角在180度附近變動時
請改用TRIC

2.-num_nodes 通常使用預設即可,要記得DEGSM node數量選擇單數9 or 11

3.-linesearch 可以試著打開backtrack

4.-ADD_NODE_TOL  你的node優化到多好才添加下一個node,可以使用0.001來使node優化的更好

5.-DQMAG_MAX  插值的ratio,對於小系統0.4 大系統0.8 一般情況 0.6

6. 待補 ...

#####
#備註#
#####

一、
GSM對於initial geometry有著較為嚴格的要求(較為敏感),因此將結構排序好是非常重要的一件事情。
我這邊可以舉個例子,對於一個最簡單的Diels-Alder(DA)反應大概是長這個樣子的:
202105061959048759..png
那麼這時候你應該將你的結構擺的像是這個樣子:
da.png
而不是這個樣子:
da2.png
雖然說下圖pyGSM一樣會試著幫你把分子移動,以致於發生反應,但通常移動過程會跨過一個barrier,如果達到truncate標準,pyGSM會誤認為中間存在一個TS
導致他從這個TS結構開始Climb image
比較好的作法(以DA反應為例)可以參考以下說明:
1.使用avogadro等軟體輸入產物結構
da3.png
2.將預計生成的鍵移除
da4.png
3.使用force filed優化結構就會自動拉開
da5.png
這時候你就可以得到擺放位置不差的初始結構與產物結構

二、
請確保你所給定的xyz,reactant部分與product部分atom mapping是一一對應的

三、
使用Q-Chem、Orca等軟體記得額外準備level of theory file並且加在command line
也就是你要算force的lot,在ORCA稱為!Engrad
例如:
gsm -xyzfile ./diels_alder.xyz -mode DE_GSM -package QChem -lot_inp_file qstart > status.log 2>&1
202105062105137633..png

四、
或許你會覺得奇怪為什麼上述DE-GSM的energy profile長得很奇怪,原因是我是拿pyGSM給的範例沒有優化過結構直接執行。
如果全程使用DFT level計算準確對會好很多!

#####
#附件#
#####

degsm.7z (324.79 KB, 下载次数 Times of downloads: 2)

评分 Rate

参与人数
Participants 6
威望 +1 eV +22 收起 理由
Reason
ggdh + 5 精品内容
ABetaCarw + 5 <font style="vertical-align: inh
zsu007 + 5 赞!
asdf + 2 赞!
sobereva + 1
hebrewsnabla + 5 好物!

查看全部评分 View all ratings

98

帖子

1

威望

540

eV
积分
658

Level 4 (黑子)

发表于 Post on 2021-6-20 17:47:32 | 显示全部楼层 Show all
请问,Growing String Method 是否可以解决过渡态附近beads数目过少、过渡态很尖锐的问题呢?
我用ASE里面的NEB方法,找到的过渡态非常尖锐。在找Hydrogen atom transfer minimal energy path 时候,这种问题经常出现。

energy profile

energy profile





31

帖子

1

威望

366

eV
积分
417

Level 3 能力者

 楼主 Author| 发表于 Post on 2021-6-20 18:08:47 | 显示全部楼层 Show all
chenxin199261 发表于 2021-6-20 17:47
请问,Growing String Method 是否可以解决过渡态附近beads数目过少、过渡态很尖锐的问题呢?
我用ASE里面 ...

我想你說的應該是node過少,在string method通常會用node
一般來說string method其中一個優點是不像NEB在高能量區域sample太多nodes
因此如果照你說的NEB算出來很sharp在string method"可能"會更嚴重,或許可以考慮減少插值ratio
-DQMAG_MAX 0.4
並且每個node優化的好一些再添加下個node
-ADD_NODE_TOL 0.005
另外是hydrogen transfer在string method時如果是用GFN-xTB必須用GFN1-xTB因為GFN2-xTB描述hydrogen transfer是比較不適合的

31

帖子

1

威望

366

eV
积分
417

Level 3 能力者

 楼主 Author| 发表于 Post on 2021-6-20 18:10:39 | 显示全部楼层 Show all
除此之外,很sharp會有什麼問題嗎,對應到hessian應該是一個很負的值
應該是不影響TS!?
又或者你是擔心太sharp對於neb的每個node之間是equaly space會有影響?

98

帖子

1

威望

540

eV
积分
658

Level 4 (黑子)

发表于 Post on 2021-6-21 08:42:01 | 显示全部楼层 Show all
ms860309 发表于 2021-6-20 18:10
除此之外,很sharp會有什麼問題嗎,對應到hessian應該是一個很負的值
應該是不影響TS!?
又或者你是擔心太 ...

sharp主要是因为,仅有3个beads(我总共用了70个beads)去描述过渡态附近的情况,其余的beads都掉到低能量区了。这样的话,最高点并不是过渡态。

31

帖子

1

威望

366

eV
积分
417

Level 3 能力者

 楼主 Author| 发表于 Post on 2021-6-21 18:52:21 | 显示全部楼层 Show all
chenxin199261 发表于 2021-6-21 08:42
sharp主要是因为,仅有3个beads(我总共用了70个beads)去描述过渡态附近的情况,其余的beads都掉到低能 ...

你可以嘗試GSM 但是不保證成功
事實上應該沒有任何一種找TS的initial guess方法可以保證找到TS
而且TS算法本身就是不好收斂的

评分 Rate

参与人数
Participants 1
eV +1 收起 理由
Reason
卡开发发 + 1 我很赞同

查看全部评分 View all ratings

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 03:03 , Processed in 0.214325 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list