天天看點

PyAmesp——Amesp與ASE的聯姻

最近,張英峰博士釋出了國産量子化學程式Amesp。為了進一步降低程式使用的門檻并彌補一些功能上的短闆,我們嘗試為Amesp增加一個接口程式——PyAmesp,通過ASE (Atomic Simulation Environment)調用Amesp進行理論計算,實作Amesp與ASE的內建與“聯姻”。本文将給出PyAmesp的安裝過程,并對ASE做簡要介紹,随後展示利用ASE調用Amesp進行結構的優化與過渡态的計算。(注:本文适合具備一定Python和ASE基礎的讀者,如果您對ASE及其使用方法不熟悉,可以登入ASE官網https://wiki.fysik.dtu.dk/ase/擷取更多資訊。)

一. ASE介紹

1. ASE簡介

ASE是一個旨在原子和電子尺度上建立、控制、可視化和分析計算結果的Python子產品。ASE提供了大量的基礎類(Class),例如“Atoms”存儲了與原子相關的屬性與位置的資訊,可以通過ASE的io子產品實作不同結構檔案之間的格式轉換。同時ASE也提供了大量電子結構計算程式代碼的接口,可将能量、能級、力、應力等諸多資訊導入ASE,進而實作以ASE作為優化器進行結構優化、分子動力學、過渡态搜尋,以及預覽分子/晶體結構、繪制聲子譜和能帶等結果的分析工具。在Python語言強大的文法和生态加持下,通過ASE可以輕松定義和控制分子或晶體結構和計算參數,實作複雜的計算化學工作流,大大提升研究者的工作效率,并在一定程度上避免重複造輪子。

2. 為什麼讓ASE與Amesp聯姻?

将ASE與Amesp進行聯姻,可以利用ASE轉換Amesp的輸入輸出檔案,并調用Amesp使用更多的優化算法進行結構優化以及CI-NEB或Dimer實作過渡态搜尋,即Amesp變為一個Calculator,由此達到從模組化到格式轉換,再到計算、分析與可視化計算結果,實作Amesp計算的“一條龍”服務。

二. PyAmesp安裝

我們假設已經安裝好了Python和ASE(Python版本≥3.6,NumPy版本≥1.11.3,ASE版本≥3.20.0。Python安裝完成後,可通過pip install ase自動安裝ASE)。PyAmesp可在github進行下載下傳:https://github.com/DYang90/PyAmesp再按照如下方式進行安裝(注:“$”表示Linux中的指令提示符,安裝時不用輸入):

$ tar -zxvf PyAmesp.tar.gz              $ cd PyAmesp              $ pip install -e           

注意:使用者應按照Amesp手冊正确設定相關環境變量。

三. PyAmesp包的使用方法

我們以經典的有機化學Diels–Alder反應為例,通過PyAmesp實作ASE調用Amesp完成結構優化和搜尋過渡态的過程。在本例中我們的計算參數是:12核心處理器且每核心記憶體最大占用為1 GB,計算的方法基組為M06-2X/6-31G**,積分格點使用grid5級别,其餘選項采用Amesp中的預設設定。

  1. 1. 反應物及産物結構優化

我們先通過分子結構的模組化程式建構反應物乙烯和1,3-丁二烯以及産物環己烯,将這些結構儲存為xyz或者Gaussian的gjf等一系列格式。然後我們編寫Python程式,通過ASE的io子產品将結構按照Atoms類格式讀入:(注:>>>為Python互動式解釋器的提示符,不用輸入)

>>> from ase import io              >>> label = 'product'              >>> atoms = io.read('%s.xyz' % label)           

按照前文提到的計算參數設定Amesp并将其指定為atoms的計算引擎:

>>> from PyAmesp import Amesp              >>> amesp = Amesp(atoms=atoms,              label=label,              maxcore=1024, npara=12,              charge=0, mult=1,              keywords=['m06-2x', '6-31g**','grid5', 'force']              )              >>> atoms.calc = amesp           

注意:進行結構優化以及過渡态計算需要計算原子受力,在keywords必須手動指定關鍵字force。接下來使用optimize子產品的BFGSLineSearch優化器對結構進行優化直至力收斂到0.02 eV/Å,并儲存優化軌迹:

>>> from ase.optimize import BFGSLineSearch              >>> dyn = BFGSLineSearch(atoms)              >>> traj = io.Trajectory('%s.traj' % label, 'w', atoms)              >>> dyn.attach(traj)              >>> dyn.run(fmax=0.02)           

上述完整程式可以參考PyAmesp/Examples/Ex1/RunAmesp.py。最終程式經過8個BFGS步驟達到收斂,并将優化過程中結構及相應的能量與受力等資訊儲存到product.traj中。通過以下指令:

$ ase gui product.traj           

利用ASE自帶的預覽器檢視結果,程式将會預設顯示優化後的分子結構以及能量相對變化,如圖1所示。在菜單欄中的【視圖】-【簡要資訊】,或是用快捷鍵Ctrl+I來檢視目前結構的能量和受力。在菜單欄中的【檔案】-【儲存】或者通過快捷鍵Ctrl+S可以将整個軌迹或者目前結構導出為其他格式。

PyAmesp——Amesp與ASE的聯姻

圖1. 利用ASE預覽器檢視優化的幾何結構和相對能量

關于優化器的選擇,除了BFGSLineSearch,ASE還提供了如FIRE、GPMin等優化器,可參考https://wiki.fysik.dtu.dk/ase/ase/optimize.html,此處不再贅述。除了ASE所提供的優化器,PyAmesp也支援使用Gaussian作為優化器通過External關鍵字對Amesp進行調用,例如對反應物結構reactant.xyz進行優化,與讀取Amesp參數類似,但從指定計算引擎和設定優化器開始,将會與上面BFGSLineSearch優化産物的例子有些許差别,需要從PyAmesp子產品導入GaussianExternal類來根據Amesp設定的參數自動産生調用腳本gaussian_amesp.py:

>>> from PyAmesp import GaussianExternal              >>> gext = GaussianExternal(label=label, parameters=amesp.parameters)              >>> gext.write_script()           

再通過ASE對Gaussian的接口指定Gaussian使用External關鍵字調用gaussian_amesp.py腳本進行優化:

>>> from ase.calculators.gaussian import Gaussian, GaussianOptimizer              >>> gau = Gaussian(label=label, external='\'python3 gaussian_amesp.py\'')              >>> opt = GaussianOptimizer(atoms, gau)              >>> opt.run(fmax='tight', steps=100, opt='nomicro,gdiis')           

檢視結果時,與使用ASE的optimize子產品不同,優化過程的相關資訊将會記錄在Gaussian的輸出檔案reactant.log中,也可以通過ASE的預覽器檢視log檔案,本例經過30步優化達到收斂。上述案例的完整程式可以參考PyAmesp/Examples/Ex2/RunAmesp.py,但需要注意的是這種使用方法還處于測試階段,暫不确定這種調用形式是否為最終設計。

2. 過渡态搜尋

在獲得反應物與産物的結構之後,可以通過ASE調用Amesp進行CI-NEB或者Dimer計算,進而得到Diels–Alder的反應過渡态。下面的案例将使用CI-NEB搜尋該反應的過渡态(該案例參考了ASE官網NEB的例子https://wiki.fysik.dtu.dk/ase/ase/neb.html)。先将前面優化得到的traj進行讀取,然後建立一個6幀的軌迹,首尾分别為反應物和産物:

>>> from ase import io              >>> ini = io.read('reactant.traj')              >>> fin = io.read('product.traj')              >>> images = [ini.copy() for i in range(6)]              >>> from PyAmesp import Amesp              >>> label = 'neb'              >>> for i,img in enumerate(images):              amesp = Amesp(…) #這裡參數與上面相同,作為示例隻是省略了沒寫              img.calc = amesp              >>> images[-1] = fin           

然後建立NEB類并進行idpp插值擷取初始的chain:

>>> from ase.neb import NEB              >>> neb = NEB(images, climb=True)              >>> neb.interpolate(method='idpp')           

最終通過LBFGS優化直至受力收斂到0.02 eV/Å并儲存最終的chain軌迹:

>>> dyn = LBFGS(neb, trajectory='%s.traj' % label)              >>> dyn.run(fmax=0.02)              >>> io.write('%s.json' % label, images)           

最終經過58步LBFGS優化找到過渡态,其結構如圖2所示。

PyAmesp——Amesp與ASE的聯姻

圖2. ASE預覽器檢視過渡态結構以及軌迹的相對能量

經過頻率計算驗證,得到一個非低頻的虛頻,完整案例可以參考PyAmesp/Examples/Ex3/RunAmesp.py。為了更高效獲得這個過渡态的結構,可以先以受力收斂0.5 eV/Å為标準使用CI-NEB優化,然後用Dimer方法做更細緻的優化。Dimer計算的具體過程這裡不再詳述,可以參考PyAmesp/Examples/Ex4/RunAmesp.py,頻率分析相關的功能目前我們還并沒有內建在PyAmesp當中。

四. 總結

本文主要介紹了PyAmesp的功能、安裝以及使用方法。通過使用Python編寫PyAmesp,使得ASE能夠調用Amesp完成優化以及過渡态等計算任務。通過ASE與Amesp的聯姻可以為量子化學計算提供更多的選擇性與更高的靈活性,進而更加高效地進行理論計算,探索物質的特性和反應機制。

随着Amesp的發展逐漸成熟,其作為Calculator的功能将會進一步擴充。同時Amesp還可以與其他程式進行聯合,實作更複雜的計算任務。相信随着時間的推移,我們将能夠利用Amesp的強大功能,以及與其他程式的聯合使用,開展更深入和複雜的探索與研究。