计算化学公社

标题: Growing String Method (GSM) - pyGSM 安裝、簡介與範例 [打印本页]

作者
Author: ms860309 时间: 2021-5-6 21:25
标题: Growing String Method (GSM) - pyGSM 安裝、簡介與範例
本帖最后由 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
此時資料夾中應該有，並且對應內容為下圖
(, 下载次数 Times of downloads: 71) (, 下载次数 Times of downloads: 75) (, 下载次数 Times of downloads: 74)
四、執行
五、結果
計算後資料夾中的圖片檔為Energy profile，TSNode.xyz為TS initial guess後續再以其他軟體優化TS
至於優化過程可以參考我附檔的scratch/opt_converged_000_000.xyz再丟到VMD看即可(記得去掉開頭的molden format)，這邊我就不做成GIF節省論壇圖片資源
(, 下载次数 Times of downloads: 55)

PART II: Single-Ended Growing String Method (SE-GSM)
與上個例子相同，差別在於不需要給product結構而是給driving coordinate，可以給ADD BREAK bond編號、angle、torsion等資訊，詳情請自行研究
(, 下载次数 Times of downloads: 67) (, 下载次数 Times of downloads: 70) (, 下载次数 Times of downloads: 74) (, 下载次数 Times of downloads: 89)

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

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)反應大概是長這個樣子的：
(, 下载次数 Times of downloads: 61)
那麼這時候你應該將你的結構擺的像是這個樣子：
(, 下载次数 Times of downloads: 67)
而不是這個樣子：
(, 下载次数 Times of downloads: 46)
雖然說下圖pyGSM一樣會試著幫你把分子移動，以致於發生反應，但通常移動過程會跨過一個barrier，如果達到truncate標準，pyGSM會誤認為中間存在一個TS
導致他從這個TS結構開始Climb image
比較好的作法(以DA反應為例)可以參考以下說明：
1.使用avogadro等軟體輸入產物結構
(, 下载次数 Times of downloads: 57)
2.將預計生成的鍵移除
(, 下载次数 Times of downloads: 47)
3.使用force filed優化結構就會自動拉開
(, 下载次数 Times of downloads: 51)
這時候你就可以得到擺放位置不差的初始結構與產物結構

二、
請確保你所給定的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
(, 下载次数 Times of downloads: 51)

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

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

(, 下载次数 Times of downloads: 16)

我的表達並不是非常好，而且對此瞭解也沒有到非常深，如果有任何問題歡迎提出，我很樂意討論！
另外排版可能也非常有問題，說明也可能不到位，有想到什麼後續會再補充！

作者
Author: chenxin199261 时间: 2021-6-20 17:47
请问，Growing String Method 是否可以解决过渡态附近beads数目过少、过渡态很尖锐的问题呢？
我用ASE里面的NEB方法，找到的过渡态非常尖锐。在找Hydrogen atom transfer minimal energy path 时候，这种问题经常出现。
(, 下载次数 Times of downloads: 52)

作者
Author: ms860309 时间: 2021-6-20 18:08

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是比較不適合的

作者
Author: ms860309 时间: 2021-6-20 18:10
除此之外，很sharp會有什麼問題嗎，對應到hessian應該是一個很負的值
應該是不影響TS!?
又或者你是擔心太sharp對於neb的每個node之間是equaly space會有影響？

作者
Author: chenxin199261 时间: 2021-6-21 08:42

ms860309 发表于 2021-6-20 18:10
除此之外，很sharp會有什麼問題嗎，對應到hessian應該是一個很負的值
應該是不影響TS!?
又或者你是擔心太 ...

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

作者
Author: ms860309 时间: 2021-6-21 18:52

chenxin199261 发表于 2021-6-21 08:42
sharp主要是因为，仅有3个beads（我总共用了70个beads）去描述过渡态附近的情况，其余的beads都掉到低能 ...

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

作者
Author: liaojianhan 时间: 2024-4-21 00:50
请问一下各位大佬，我用pygsm的双端gsm方法生成过渡态初猜的结构，但是我不太确定怎么分析那个计算结果，求助（官方那个说明文档感觉啥也没说）。
我就是在想，他每个node的能量算出来都是负的，那我怎么拿能量的最高点作为初猜呢，我不知道我是不是算得有问题

作者
Author: liaojianhan 时间: 2024-4-23 01:12
gsm算法的图大概长这样，想知道这样是不是不太合理，以及我又应该如何去调整呢？

欢迎光临计算化学公社 (http://bbs.keinsci.com/)