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,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的强大功能，以及与其他程序的联合使用，开展更深入和复杂的探索与研究。