最近，張英峰博士釋出了國産量子化學程式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,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可以将整個軌迹或者目前結構導出為其他格式。
圖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所示。
圖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的強大功能，以及與其他程式的聯合使用，開展更深入和複雜的探索與研究。