空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

ASE-ABACUS | 第三章：单端过渡态搜索

ASE

ABACUS

ATST-Tools

催化计算

Sella

Dimer

ASEABACUSATST-Tools催化计算TSSellaDimer

量子御坂

更新于 2024-10-25

推荐镜像 :atst-tools:1.5.0

推荐机型 :c2_m4_cpu

ASE-ABACUS | 第三章：单端过渡态搜索

1. 过渡态与过渡态搜索

1.1 简介

1.2 单端过渡态搜索基本思路

2. 两类单端过渡态搜索的基本原理

2.1 Dimer法

2.1.1 旋转Dimer

2.1.2 平移Dimer

2.2 P-RFO法

2.2.1 RFO方法

2.2.2 P-RFO方法

2.3 小结

3. Dimer方法使用

3.1 ATST-Tools中的Dimer

3.2 案例：H2在Au111表面的解离

4. 基于P-RFO的Sella方法使用

4.1 Sella简介与使用

4.2 案例一：H2在Au111表面的解离及其IRC

4.3 案例二：环己烷在Pt1@graphene上的解离

4.4 案例三：Au-Pd掺杂ZSM5分子筛内的甲烷活化

5 总结

ASE-ABACUS | 第三章：单端过渡态搜索

代码

文本

【ATST-Tools】，全名【Advanced ASE Transition State Tools for ABACUS and Deep Potential】, 是基于ASE调用ABACUS和DP势函数进行过渡态搜索计算的工作流脚本集，涵盖ASE中包含的NEB，AutoNEB和Dimer方法，已接入高效单端鞍点优化器Sella，并预期未来基于ASE以及ASE-ABACUS接口进一步接入或开发其他的过渡态搜索方法.

在完成这个教程之后，你将能够：

了解单端过渡态搜索的基础知识，以及它与NEB等过渡态搜索方法的异同
掌握在ATST-Tools工作流组织下，使用ASE调用ABACUS进行Dimer, Sella等单端过渡态搜索方法的操作

欢迎大家提出改进意见！

代码

文本

©️ Copyright 2024 @ Authors
作者： 量子御坂（Quantum Misaka） *** 📨
日期：2024-10-24
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的 开始连接按钮，选择 atst-tools:1.5.0 镜像 和任意配置机型即可开始。注：由于具体计算的案例都需要比较长的计算时间，本案例不会要求完成具体计算，只涉及预处理和结果分析，故使用c2_m4_cpu即可。如果你想进行具体计算，可使用c16_m32_cpu, c32_m64_cpu等具有较高核数的机型

代码

文本

代码

文本

1. 过渡态与过渡态搜索

让我们从复习一下过渡态和过渡态搜索的基本概念开始，这一部分与ASE-ABACUS应用第二章的内容重复，已经看过上一个Notebook的同学可以快速跳过。

1.1 简介

化学是一门研究物质组成，结构及其变化的学科，其主要的研究对象都是物质之间的相互反应，这些反应过程往往由一系列基元反应组成，而这一系列基元反应的发生都离不开反应的过渡态。

过渡态由过渡态理论引出，它认为：在A -> B的反应过程中，反应物A会经历一个能量较高的临界状态，并经此亚稳定态转化为产物B，这种临界状态即为过渡态。它是反应能垒的来源，并与反应速率直接相关。

化学反应图示。黑线为能量曲线，「反应物, 生成物」对应极点，「过渡态」对应鞍点

在（基于BO近似的）势能面上，这种过渡态位置表现为马鞍点，其能量对坐标的一阶梯度为0，且二阶梯度仅在反应坐标方向为负，其他位置为正，并且也仅在反应坐标方向上存在一个虚频。反应经由这样一个**最小能量路径（Minimum Energy Path， MEP）**进行，从反应物经过渡态生成产物。

势能面上的极小点（稳定结构），过渡态（马鞍点结构）与最小能量路径

基于统计热力学的方法可以推导得到Eyring过渡态理论，该理论指出，基元反应的活化自由能，即过渡态自由能与反应物自由能之差，可以直接与反应速率常数联系起来，即： $k = κ \frac{k _{B} T}{h} e^{\frac{- Δ G ^{‡}}{RT}}$ 从而基元反应的反应速率可以通过反应物与过渡态的信息计算得到。

过渡态是一种典型的稀有事件，无论是扩散过程、解离过程、还是其他各种均相/多相的化学反应过程，其微观机制均涉及对过渡态的分析，从而过渡态的理论研究与实验检测都是化学研究的重要组成部分。同时，过渡态搜索则是通过理论模拟方法寻找过渡态的关键。过渡态搜索的方法有很多，其中最常用的是NEB和Dimer方法。当然，新的过渡态搜索方法也在不断涌现，尤其是近年来我们很惊喜地看到基于机器学习的生成扩散模型在过渡态搜索与反应生成方面起到了非常惊艳的效果（戳这里）。

在上一章中，我们在回顾过渡态理论与计算概念的同时，介绍了NEB方法的使用，同时引入了使用ASE-ABACUS进行过渡态计算的工作流脚本集ATST-Tools。可戳：ASE-ABACUS应用第二章。而本章的重点将会放在单端过渡态搜索上。

代码

文本

1.2 单端过渡态搜索基本思路

在多相催化模拟中，NEB方法和Dimer方法是最为常用的两个。这二者正好代表着搜索过渡态的两类不同方法，其一是基于已知反应物和产物的信息去搜索过渡态，即双端（Double-End）过渡态方法，其二是已知过渡态的一个初始猜测去搜索过渡态，即单端（Single-End）过渡态方法。

然而，能用于搜索过渡态的方法其实非常之多，北京科音的卢天老师曾经总结过截至19年以来较为有价值的过渡态方法，可参见博文：过渡态、反应路径的计算方法及相关问题。同时，卢天老师的另一篇博文：简谈Gaussian里找过渡态的关键词opt=TS和QST2、3中，分析了量子化学软件Gaussian内的过渡态搜素关键词以及其算法原理，并指出Gaussian中opt=TS关键词是最为建议适用的过渡态方法，这一方法只需要提供过渡态初猜结构，即可进行以过渡态所在一阶鞍点位置为目标的优化过程。

代码

文本

这类单端过渡态方法做法各不相同，但它们总的原理均是类似的，都是在简谐近似下，利用了体系的势能面二阶导，即Hessian矩阵（或者近似Hessian矩阵）的信息，来进行过渡态搜索。这背后蕴藏的原理其实非常直观：过渡态其实就是势能面的一阶鞍点，也即仅在一个方向上二阶导小于零的鞍点（用更数学的语言就是：二阶导矩阵仅有一个负本征值）。而这种鞍点搜索的方法都可以用类似于牛顿法的求解方式完成。 $x_{n + 1} = x_{n} - H (x_{n})^{- 1} F (x_{n})$ 其中， $F (x_{n})$ 为势能面一阶导（即力向量）， $H (x_{n})$ 即为势能面的二阶导矩阵（Hessian矩阵）。当然了，牛顿法虽然是普适的鞍点搜索方法，却更适用于局部极小点（而非一阶鞍点）的搜索，尤其是在Hessian矩阵具有较好正定性的时候。当过渡态初猜位置就在势能面鞍点的二次区域内时，采用牛顿法也能优化到过渡态，但是这一情形是相对罕见的。

从化学的角度来思考Hessian矩阵的话，会发现Hessian矩阵的leftmost eigenvector对应的就是本征值最小的（虚）振动模，也即振子强度最小的振动模，它对应最容易发生化学反应的方向。

注：leftmost eigenmode中文翻译就是最小本征值对应的向量的自由度，这个leftmost来源于最优化学科在对角化Hessian矩阵之后习惯将本征值从小到大排列

这意味着，如果我们要进行高效的过渡态搜索（而非直接应用牛顿法，导致体系朝着局部极小点（而非过渡态）进行搜索的话，需要在进行优化算法的时候，让体系在leftmost eigenmode的方向，朝着梯度（而非负梯度）的方向进行优化，从而保证在这一反应坐标方向，体系是朝着能量上升的方向优化的。

这也就引出了单端过渡态搜索方法的基本流程

求解体系在势能面上的梯度 $g$ 和Hessian矩阵 $H$ 的信息，进而得到体系的最小本征值及其本征矢 (即leftmost eigenvector) 。
利用这些信息进行优化，使leftmost eigenmode方向朝着能量升高的方向演化，其他eigenmode朝着能量降低的方向演化，最终找到过渡态。

代码

文本

在此基础上，本notebook将简单介绍两类最为常见的单端过渡态搜索方法，也是目前ATST-Tools里面接入的两类方法。

Dimer法
P-RFO法

更多的过渡态搜索和反应路径计算可以参考卢天老师的博文：

http://sobereva.com/44

以及对应参考文献

http://bbs.keinsci.com/thread-105-1-1.html

代码

文本

2. 两类单端过渡态搜索的基本原理

代码

文本

2.1 Dimer法

Dimer方法巧妙地绕开了Hessian矩阵的求解，通过定义一个Dimer直接求解leftmost eigenmode，从而实现了高效的过渡态定位。

具体来说，Dimer法在势能面上定义了由两个点 $R_{1}$ 和 $R_{2}$ 组成的一个Dimer，如下图所示：

Dimer定义及各物理量图示

在这个Dimer中，两个点的能量和受力分别定为 $E_{1}$ 、 $F_{1}$ 、 $E_{2}$ 、 $F_{2}$ ，Dimer方向单位向量取为 $\hat{N}$ ，Dimer在势能面上的距离取为 $Δ R = R_{1} = R_{2}$ 。这一值通常很小，小于0.1 Angstrom。Dimer的中间点为 $R$ ，受力 $F (R) = (F_{1} + F_{2}) /2$ 。Dimer的总能量为 $E = (E_{1} + E_{2}) /2$

代码

文本

Dimer的结构演化步分为两个步骤：

旋转Dimer
平移Dimer

其中，旋转Dimer的目的是寻找势能面当前位置处的leftmost eigenmode，而平移Dimer则是基于该结果进行结构演化。

具体来说：

2.1.1 旋转Dimer

该步保持 $R_{1}, R_{2}$ 中点位置 $R$ 不变，以此为轴在势能面上旋转Dimer，直至总能量 $E$ 最小。

因为Dimer的长度 $Δ R$ 很小，应用有限差分方法可得Dimer方向的势能面曲率为：

$C = \hat{N}^{T} H \hat{N} \approx \frac{( F _{2} - F _{1} ) \cdot N ^}{2Δ R} = \frac{E - 2 E _{0}}{( Δ R ) ^{2}}$

代码

文本

也即在旋转过程中，总能量 $E$ 与曲率 $C$ 满足线性关系。最小化 $E$ 的过程就是最小化 $C$ 的过程，所以优化结束后，Dimer就会处于最小曲率方向。当Dimer优化到过渡态时，Dimer方向即为虚频方向。

这一优化过程一般通过基于有限差分的方法完成

将Dimer旋转一个很小的 $d θ$ 角度，通过有限差分的方法求解Dimer的法向受力（旋转力） $F^{'}$ $F^{'} = \frac{d F}{d θ} \approx \frac{F ^{*} \cdot Θ ^ ^{*} - F \cdot Θ ^}{d θ}_{θ = d θ /2}$
基于该旋转力进行Dimer的旋转优化，这一步和普通结构优化一致，最初采用牛顿法，当前的Dimer算法里大多采用CG（共轭梯度）或L-BFGS法，并且在优化的过程中结合线性插值的方法减少势能面计算。

过程大致如下图所示：

旋转Dimer图示

代码

文本

不难发现，在做旋转Dimer之前定义Dimer的时候，需要定义一个初猜的Dimer方向，这一方向通常由两种方式

用户自定义（一般用于已有振动分析或NEB轨迹等信息时）
采用高斯分布随机生成

代码

文本

2.1.2 平移Dimer

使用旋转Dimer收敛之后的Dimer方向 $\hat{N}$ 和Dimer受力 $F (R)$ 进行平移。平移操作也是普通结构优化任务，在做平移时，Dimer的有效受力定义如下：

$F_{e ff} = {- F^{∥} = - (F (R)^{†} \hat{N}) \hat{N} F (R) - 2 F^{∥} if C > 0 if C < 0$ 其中 $F^{∥}$ 指平行于 $\hat{N}$ 方向的 $F$ 。

这一有效受力的目的如下：

当曲率 $C > 0$ 时，认为Dimer位于鞍点附近区域以外，不存在负值最小本征值的本征模，此时使总的有效力为平行于Dimer方向的力的负值。通过这一有效受力促使Dimer离开该区域。
当曲率 $C < 0$ 时，认为Dimer位于鞍点区内，将平行于Dimer方向的力反号。使得在演化步进时，平行于Dimer方向的演化朝能量升高方向，垂直于Dimer方向的演化朝能量降低方向，使Dimer中心逐步靠近鞍点。

平移Dimer图示

在具体执行步进操作的时候，可以采用CG或Quick-Min等优化方法步进方式。

代码

文本

Dimer法实际上也源于牛顿法，但它通过（基于一阶导数有限差分的）旋转Dimer方法代替了计算Hessian矩阵最小本征模的过程，节省了计算成本。Dimer对初始位置要求较为宽松，并不要求在过渡态二次区域内。诸如Dimer这样的方法可以用于NEB计算后的过渡态精优化，也可以做鞍点搜索。

代码

文本

2.2 P-RFO法

Partitional Rational Function Optimization (分部有理函数优化)方法是基于牛顿法原理改进而来的另一种常用的一阶（以及更高阶的）鞍点优化方法。

2.2.1 RFO方法

有理函数优化(Rational Functional Optimization)方法从形式上来说就是对牛顿法的公式做了个简单改动，将牛顿法的基础——目标函数（即势能面）的二阶Taylor展开： $ε = E (x) - E_{0} = g^{+} x + \frac{1}{2} x^{+} Hx$ 变换为了如下形式： $ε = E (x) - E_{0} = \frac{g ^{+} x + \frac{1}{2} x ^{+} Hx}{1 + x ^{+} Sx}$ 其中 $S$ 为一实对称矩阵（和 $H$ 一样）。定义 $λ = 2 ε$ 可以得到一个(n+1)维的本征方程： $(H g^{+} g 0) (1 x) = λ (s 0 01) (1 x)$ 基于这一本征方程，可以求解出(n+1)个 $λ$ 本征值。

代码

文本

同时，上述方程可以变换为两个方程的方程组： $(H + λ S) x + g = 0$

$g^{+} x = λ$

基于这一方程组，若Hessian矩阵 $H$ 已对角化，其各个本征向量方向的本征值为 $h_{i}$ ，可以得到RFO方法优化时的步进公式： $x_{i} = \frac{g _{i}}{λ - h _{i}}$ 对比相同情况下牛顿法的步进公式： $x_{i} = - \frac{g _{i}}{h _{i}}$ 可以发现RFO方法就优化算法而言起到的作用是在牛顿法的分母处加了一个较正因子，这也是牛顿法步进公式的一个主要改进，可以防止Hessian矩阵本征值接近于0的情况。

代码

文本

求解上述本征方程得到的(n+1)个 $λ$ 本征值，将它们从小到大排列，它与Hessian矩阵本征值 $h_{i}$ 可以被证明满足如下关系： $λ_{i} \leq h_{i} \leq λ_{i + 1}$ 于是可以通过求解得到的 $λ_{i}$ 值，可以将Hessian矩阵的各个本征值分开。

此时，如果选择求解得到的 $λ$ 值中最小的 $λ$ ，即可确保RFO步进公式中各个本征模方向 $h_{i} - λ$ 的值均为正，保证每个本征模上的势能面演化都朝着能量下降方向。

而针对鞍点搜索任务，可以选择各 $λ$ 值中，大于第 $m$ 个Hessian本征值的 $λ$ ，此时势能面演化将在本征值最低的 $m$ 个本征模方向上沿能量升高方向演化，在其余 $N - m$ 个方向上朝着能量降低方向演化，从而最终优化到 $m$ 阶的鞍点。 $m = 1$ 时即为一阶鞍点，即选取倒数第二小的 $λ$ 即可开展过渡态搜索

代码

文本

2.2.2 P-RFO方法

代码

文本

上述RFO方法可用于极小化和鞍点优化，但一般还是用于极小化任务居多（据说Gaussian软件做结构优化厉害的地方之一就在于他们的RFO算法写的很好）

以下讨论的P-RFO方法则专用于过渡态优化。传统RFO对所有简正模方向的优化步进均使用同一个 $λ$ ，而P-RFO则显式将过渡态方向和其他方向分开，针对过渡态方向通过求解如下本征方程的一维简化形式： $λ_{TS} = \frac{g _{TS}^{2}}{h _{TS} - λ _{TS}}$ 取该一元二次方程最大的解，就能保证在这个方向上进行步进时朝着能量升高的方向。其余方向的 $λ$ 确定和RFO方法一致。

代码

文本

多数时候，这里所谓的“指向过渡态的方向”即leftmost eigenmode的方向，对应最低的本征值，在RFO方法和Dimer方法中均是如此假设。但有时也会是其它的非最低的模式。P-RFO也可以选择这样的模式作为指向过渡态的模式，相比RFO更加灵活，且在实际应用中也更加高效。

代码

文本

2.3 小结

可以看出Dimer法和P-RFO方法本质都源于牛顿法，且都聚焦于将对应过渡态的简正模与其它的简正模区分开来。不同的是Dimer采用旋转Dimer的方法回避了Hessian矩阵的直接求解，使得Dimer方法一度是相对更高效的过渡态搜索方法，被更广泛地应用在过渡态搜索，尤其是表面反应过渡态的搜索任务中。

但这并不代表P-RFO方法就无人问津，在Hessian矩阵易于求解的场合，P-RFO方法一般都比Dimer方法更高效。在Gaussian软件最常用的Berny过渡态搜索方法（即opt=TS）中，就采用了基于内坐标的P-RFO方法。并且，由于Dimer的旋转Dimer和平移Dimer步都依赖于有限差分，能选用的步进优化方法并不多（一般被限制在BFGS和CG），并且存在依赖初始Dimer方向，迭代步数过长等问题。

近年来，随着Hessian矩阵求解算法和所需算力的迭代，P-RFO方法在周期性体系的应用成为可能。下述的Sella方法是其中的一个代表。

代码

文本

3. Dimer方法使用

代码

文本

3.1 ATST-Tools中的Dimer

ATST-Tools内集成了对ASE内置Dimer方法的调用，它可以通过如下链接访问： https://github.com/QuantumMisaka/ATST-Tools 。其README中有详细的环境配置和使用方法介绍。

在本notebook的推荐镜像中已经装有ATST-Tools及其依赖环境。我们可以进入相关目录查看

代码

文本

[1]

cd /opt/ATST-Tools

/opt/ATST-Tools

代码

文本

[2]

README.md  dimer/     img/  relax/  source/
ase-dp/    examples/  neb/  sella/  vibration/

代码

文本

其中Dimer和Sella文件夹内的执行脚本，以及example文件夹里面的相关案例，为本期Notebook的介绍重点。

ATST-Tools的介绍，配置以及NEB计算可见上期Notebook: ASE-ABACUS应用第二章。在使用ATST-Tools之前，强烈建议认真阅读README相关内容

我们可以看看Dimer的运行脚本集

代码

文本

[1]

cd /opt/ATST-Tools/dimer

/opt/ATST-Tools/dimer

代码

文本

[11]

dimer_run.py  dimer_submit.sh  neb2dimer.py  neb2dimer_abacus.py

代码

文本

可以看到其中有四个脚本，这四个脚本的功能分别如下：

dimer_run.py: 调用ABACUS开展Dimer计算的主要执行脚本
dimer_submit.sh：将Dimer计算主脚本提交到Slurm等服务器队列的队列脚本
neb2dimer.py：通过解析NEB轨迹生成Dimer初猜的脚本
neb2dimer_abacus.py：将NEB和Dimer联动，调用ABACUS用于过渡态计算的脚本

本次教程主要解析如何开展Dimer计算。可以先来看看dimer_run.py脚本的具体内容

代码

文本

[2]

cat dimer_run.py

# JamesMisaka in 2024-09-16
# Run Dimer calculation based on initial guess by ASE-ABACUS
# part of ATST-Tools scripts

from ase.io import Trajectory, read, write
import numpy as np
from abacus_dimer import AbacusDimer

fmax = 0.05
dimer_input_file = 'dimer_init.traj' # "STRU"
init_eigenmode_method = 'displacement'
displacement_input = 'displacement_vector.npy'

# setting for calculator
# abacus calculator setting
abacus = "abacus"
mpi=16
omp=4
moving_atoms_ind = None
lib_dir = "/lustre/home/2201110432/example/abacus"
#lib_dir = "/data/home/liuzq/example"
pseudo_dir = f"{lib_dir}/PP"
basis_dir = f"{lib_dir}/ORB"
pp = {
        "H":"H_ONCV_PBE-1.0.upf",
        "C":"C_ONCV_PBE-1.0.upf",
        "O":"O_ONCV_PBE-1.0.upf",
        "Fe":"Fe_ONCV_PBE-1.0.upf",
      }
basis = {
        "H":"H_gga_6au_100Ry_2s1p.orb",
        "C":"C_gga_7au_100Ry_2s2p1d.orb",
        "O":"O_gga_7au_100Ry_2s2p1d.orb",
        "Fe":"Fe_gga_8au_100Ry_4s2p2d1f.orb",
         }
kpts = [3, 1, 2]
parameters = {
    'calculation': 'scf',
    'nspin': 2,
    'xc': 'pbe',
    'ecutwfc': 100,
    'ks_solver': 'genelpa',
    'symmetry': 0,
    'vdw_method': 'none',
    'smearing_method': 'mp',
    'smearing_sigma': 0.002,
    'basis_type': 'lcao',
    'mixing_type': 'broyden',
    'mixing_ndim': 20,
    'scf_thr': 1e-7,
    'scf_nmax': 100,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'cal_force': 1,
    'cal_stress': 1,
    'init_wfc': 'atomic',
    'init_chg': 'auto',
    'out_stru': 1,
    'out_chg': 0,
    'out_mul': 1,
    'out_wfc_lcao': 0,
    'out_bandgap': 0,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
}

        
if __name__ == "__main__":
# running process
# read initial guessed neb chain
    dimer_init = read(dimer_input_file)
    displacement_vector = np.load(displacement_input)
    dimer = AbacusDimer(dimer_init, parameters, abacus=abacus,
                        mpi=mpi, omp=omp, 
                        init_eigenmode_method=init_eigenmode_method,
                        displacement_vector=displacement_vector)
    dimer.run(fmax=fmax, moving_atoms_ind=moving_atoms_ind)

    # output TS stru file
    write("TS_dimer.stru", dimer_init, format="abacus")
    write("TS_dimer.cif", dimer_init, format="cif")

代码

文本

以if __name__ =="__main__"部分作为分界线，脚本前半部分主要是ABACUS计算相关设置和Dimer计算相关设置，其中：

dimer_input_file 为Dimer计算初始化时的初始结构路径，可以为任何ASE能解析的文件格式
init_eignemode_method 为Dimer计算初始化时的Dimer方向选取方式，displacement模式下，Dimer方向将由外置npy文件提供，即下一个displacement_input参数对应的路径。另一个可选模式是gauss，此时ASE将以高斯分布的方式随机生成一个Dimer方向。
fmax 为Dimer收敛时体系应有的最大受力
moving_atom_ind 可以在Dimer计算初始化

脚本后半部分为具体计算工作流，先读入初始化信息，再调用source/abacus_dimer.py里面的AbacusDimer库进行Dimer计算，最后输出计算得到的过渡态结构。能量信息会在Dimer运行中的标准输出内给出。

代码

文本

我们可以来看看案例，以下是一个简单且被多次复现过的案例：

3.2 案例：H2在Au111表面的解离

我们可以看看examples文件夹里面的案例：

代码

文本

[13]

cd /opt/ATST-Tools/examples

/opt/ATST-Tools/examples

代码

文本

[14]

 CH4-HOAuHOPd-ZSM5/   CO-Pt111/  'Cy-Pt@graphene'/   H2-Au111/   Li-diffu-Si/

代码

文本

打开H2-Au111文件夹，查看里面包含的计算方法案例

代码

文本

[20]

cd /opt/ATST-Tools/examples/H2-Au111

/opt/ATST-Tools/examples/H2-Au111

代码

文本

[22]

autoneb/  dimer/  neb/        neb2sella/  sella_IRC/
data/     dyneb/  neb2dimer/  sella/

代码

文本

可以看到这个案例文件夹内包含几乎所有ATST-Tools支持的计算方法案例。我们先聚焦于Dimer

代码

文本

[10]

cd /opt/ATST-Tools/examples/H2-Au111/dimer

/opt/ATST-Tools/examples/H2-Au111/dimer

代码

文本

[24]

STRU             displacement_vector.npy  running_dimer.out
dimer_init.traj  run_dimer.traj

代码

文本

这个例子里已经提供了dimer计算所需的各个input,即dimer_init.traj和displacement_vector.py，它们来自于NEB计算，可以通过python neb2dimer.py neb_latest.traj通过NEB结果得到。

复制dimer_run.py脚本到这个目录下，稍加编辑之后即可运行这一案例。

代码

文本

[26]

%%writefile dimer_run.py

from ase.io import Trajectory, read, write

import numpy as np

from abacus_dimer import AbacusDimer

fmax = 0.05

dimer_input_file = 'dimer_init.traj'

init_eigenmode_method = 'displacement'

displacement_input = 'displacement_vector.npy'

# setting for calculator

# abacus calculator setting

abacus = "abacus"

mpi=16

omp=1

moving_atoms_ind = None

lib_dir = "../data

pseudo_dir = f"{lib_dir}/"

basis_dir = f"{lib_dir}/"

pp = {

"H":"H_ONCV_PBE-1.0.upf",

"Au":"Au_ONCV_PBE-1.0.upf",

}

basis = {

"H":"H_gga_6au_100Ry_2s1p.orb",

"Au":"Au_gga_7au_100Ry_4s2p2d1f.orb",

}

kpts = [3, 1, 3]

parameters = {

'calculation': 'scf',

'nspin': 2,

'xc': 'pbe',

'ecutwfc': 100,

'ks_solver': 'genelpa',

'symmetry': 0,

'vdw_method': 'd3_bj',

'smearing_method': 'gaussian',

'smearing_sigma': 0.002,

'basis_type': 'lcao',

'mixing_type': 'broyden',

'mixing_ndim': 20,

'scf_thr': 1e-6,

'scf_nmax': 100,

'kpts': kpts,

'pp': pp,

'basis': basis,

'pseudo_dir': pseudo_dir,

'basis_dir': basis_dir,

'cal_force': 1,

'cal_stress': 1,

'init_wfc': 'atomic',

'init_chg': 'atomic',

'out_stru': 1,

'out_chg': -1,

'out_mul': 1,

'out_wfc_lcao': 0,

'out_bandgap': 0,

'efield_flag': 1,

'dip_cor_flag': 1,

'efield_dir': 1,

}

if __name__ == "__main__":

# running process

# read initial guessed neb chain

dimer_init = read(dimer_init)

displacement_vector = np.load(displacement_input)

dimer = AbacusDimer(dimer_init, parameters, abacus=abacus,

mpi=mpi, omp=omp,

init_eigenmode_method=init_eigenmode_method,

displacement_vector=displacement_vector)

dimer.run(fmax=fmax, moving_atoms_ind=moving_atoms_ind)

# output TS stru file

write("TS_dimer.stru", dimer_init, format="abacus")

write("TS_dimer.cif", dimer_init, format="cif")

Writing dimer_run.py

代码

文本

[ ]

# 去掉注释可以运行上述脚本，或者把上面的writefile去掉直接运行

# ! python dimer_run.py

代码

文本

这一Dimer的实际运行需要较多的核数和较长的时间，我们直接来看计算结果。以下计算采用ABACUS 3.7.5版本，在64核的Intel 8358 CPU上完成。在Dimer达到收敛时，有两个需要被满足的阈值

fmax小于阈值
Dimer对应的曲率为负

代码

文本

[27]

cat running_dimer.out

MINMODE:METHOD: Dimer
DIM:CONTROL: Search Parameters:
DIM:CONTROL: ------------------
DIM:CONTROL: eigenmode_method = dimer
DIM:CONTROL: f_rot_min = 0.1
DIM:CONTROL: f_rot_max = 1.0
DIM:CONTROL: max_num_rot = 1
DIM:CONTROL: trial_angle = 0.7853981633974483
DIM:CONTROL: trial_trans_step = 0.001
DIM:CONTROL: maximum_translation = 0.1
DIM:CONTROL: cg_translation = True
DIM:CONTROL: use_central_forces = True
DIM:CONTROL: dimer_separation = 0.0001
DIM:CONTROL: initial_eigenmode_method = displacement
DIM:CONTROL: extrapolate_forces = False
DIM:CONTROL: displacement_method = vector
DIM:CONTROL: gauss_std = 0.1
DIM:CONTROL: order = 1
DIM:CONTROL: mask = [True, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]
DIM:CONTROL: displacement_center = None
DIM:CONTROL: displacement_radius = None
DIM:CONTROL: number_of_displacement_atoms = None
DIM:CONTROL: ------------------
DIM:ROT: OPT-STEP ROT-STEP CURVATURE ROT-ANGLE ROT-FORCE
MINMODE:DISP:       0     -0.00413662      0.00119250     -0.00553278
MINMODE:DISP:       1      0.00280987      0.00103319      0.00633761
MINMODE:DISP:       2      0.00000000      0.00000000      0.00000000
MINMODE:DISP:       3      0.00000000      0.00000000      0.00000000
MINMODE:DISP:       4      0.00000000      0.00000000      0.00000000
MINMODE:DISP:       5      0.00000000      0.00000000      0.00000000
MINMODE:DISP:       6      0.00000000      0.00000000      0.00000000
MINMODE:DISP:       7      0.00000000      0.00000000      0.00000000
MINMODE:DISP:       8      0.00000000      0.00000000      0.00000000
MINMODE:DISP:       9      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      10      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      11      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      12      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      13      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      14      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      15      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      16      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      17      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      18      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      19      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      20      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      21      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      22      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      23      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      24      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      25      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      26      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      27      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      28      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      29      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      30      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      31      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      32      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      33      0.00000000      0.00000000      0.00000000
MINMODE:DISP:      34     -0.00003862     -0.00003197     -0.00010226
MINMODE:DISP:      35      0.00005203      0.00007139      0.00004057
MINMODE:DISP:      36      0.00000466      0.00002195     -0.00002134
MINMODE:DISP:      37     -0.00003895      0.00003978      0.00002570
MINMODE:DISP:      38      0.00002282      0.00009588      0.00005252
MINMODE:DISP:      39     -0.00002487      0.00027001      0.00007477
MINMODE:DISP:      40     -0.00009639      0.00000744      0.00005492
MINMODE:DISP:      41      0.00001900      0.00001689      0.00000790
MINMODE:DISP:      42      0.00009716      0.00009982      0.00004110
MINMODE:DISP:      43      0.00001521      0.00016641      0.00000025
MINMODE:DISP:      44     -0.00008011      0.00026788     -0.00001045
MINMODE:DISP:      45     -0.00000746      0.00006844     -0.00006528
MINMODE:DISP:      46     -0.00000367      0.00008507      0.00000481
MINMODE:DISP:      47      0.00003131      0.00012350     -0.00010242
MINMODE:DISP:      48     -0.00004767      0.00010900     -0.00004354
MINMODE:DISP:      49      0.00006021     -0.00003258      0.00007823
MINMODE:DISP:      50     -0.00001437      0.00011392     -0.00000110
MINMODE:DISP:      51     -0.00003346      0.00009181      0.00005044
MINMODE:DISP:      52      0.00007854      0.00006127     -0.00015137
MINMODE:DISP:      53      0.00000766      0.00016675     -0.00010759
MINMODE:DISP:      54      0.00009892      0.00016595      0.00003957
MINMODE:DISP:      55      0.00002488     -0.00021964      0.00004482
MINMODE:DISP:      56     -0.00018386      0.00018577      0.00021560
MINMODE:DISP:      57     -0.00000604     -0.00004903      0.00000596
MINMODE:DISP:      58      0.00014823      0.00006514      0.00001189
MINMODE:DISP:      59     -0.00012589      0.00027576     -0.00018437
MINMODE:DISP:      60     -0.00027350      0.00061817      0.00013959
MINMODE:DISP:      61     -0.00028025      0.00018537      0.00004952
MINMODE:DISP:      62     -0.00007085      0.00009471     -0.00000201
MINMODE:DISP:      63      0.00023090      0.00006368     -0.00012927
MINMODE:DISP:      64      0.00013885      0.00041092      0.00014208
MINMODE:DISP:      65     -0.00002702     -0.00021591     -0.00005241
MinModeTranslate: STEP      TIME          ENERGY    MAX-FORCE     STEPSIZE    CURVATURE  ROT-STEPS
DIM:ROT:       0         1   -6.1854   24.9082    5.6953
MinModeTranslate:    0  22:02:37  -239254.830167       0.6718     --------    -6.185415          1
MinModeTranslate:    0  22:06:32  -239254.830167       0.6718     0.100000    -6.185415          1
DIM:ROT:       1         1   -7.8175   13.1388    3.2966
MinModeTranslate:    1  22:12:26  -239254.917191       0.4940     --------    -7.817477          1
MinModeTranslate:    1  22:16:22  -239254.917191       0.4940     0.100000    -7.817477          1
DIM:ROT:       2         1   -7.0323   12.2143    2.7513
MinModeTranslate:    2  22:22:15  -239254.976216       0.2836     --------    -7.032276          1
MinModeTranslate:    2  22:26:12  -239254.976216       0.2836     0.100000    -7.032276          1
DIM:ROT:       3         1   -7.1448   14.5687    4.4094
MinModeTranslate:    3  22:32:05  -239255.018172       0.2862     --------    -7.144755          1
MinModeTranslate:    3  22:36:01  -239255.018172       0.2862     0.100000    -7.144755          1
DIM:ROT:       4         1   -5.8353   10.9834    2.9562
MinModeTranslate:    4  22:41:55  -239255.045788       0.3741     --------    -5.835290          1
MinModeTranslate:    4  22:45:51  -239255.045788       0.3741     0.100000    -5.835290          1
DIM:ROT:       5         1   -5.8347   11.8986    3.7521
MinModeTranslate:    5  22:51:45  -239255.074815       0.1507     --------    -5.834693          1
MinModeTranslate:    5  22:55:40  -239255.074815       0.1507     0.079530    -5.834693          1
DIM:ROT:       6         1   -5.9865   12.7154    3.5773
MinModeTranslate:    6  23:01:35  -239255.087236       0.1817     --------    -5.986482          1
MinModeTranslate:    6  23:05:30  -239255.087236       0.1817     0.100000    -5.986482          1
DIM:ROT:       7         1   -6.2198   13.6302    3.8854
MinModeTranslate:    7  23:11:12  -239255.108901       0.1238     --------    -6.219802          1
MinModeTranslate:    7  23:14:59  -239255.108901       0.1238     0.051491    -6.219802          1
DIM:ROT:       8         1   -5.5465   18.3001    5.5101
MinModeTranslate:    8  23:20:40  -239255.107851       0.1596     --------    -5.546527          1
MinModeTranslate:    8  23:24:27  -239255.107851       0.1596     0.053081    -5.546527          1
DIM:ROT:       9         1   -5.5420   15.9728    3.9643
MinModeTranslate:    9  23:30:07  -239255.113742       0.1357     --------    -5.541979          1
MinModeTranslate:    9  23:33:55  -239255.113742       0.1357     0.036917    -5.541979          1
DIM:ROT:      10         1   -6.2439   22.1303    5.1771
MinModeTranslate:   10  23:39:35  -239255.117778       0.0795     --------    -6.243911          1
MinModeTranslate:   10  23:43:23  -239255.117778       0.0795     0.044007    -6.243911          1
DIM:ROT:      11         1   -6.5080   22.0600    5.2399
MinModeTranslate:   11  23:49:04  -239255.120616       0.0671     --------    -6.508041          1
MinModeTranslate:   11  23:52:51  -239255.120616       0.0671     0.043712    -6.508041          1
DIM:ROT:      12         1   -5.0832   16.3283    3.7101
MinModeTranslate:   12  23:58:32  -239255.123520       0.0720     --------    -5.083194          1
MinModeTranslate:   12  00:02:19  -239255.123520       0.0720     0.047643    -5.083194          1
DIM:ROT:      13         1   -5.0205   16.1546    4.2298
MinModeTranslate:   13  00:08:00  -239255.125296       0.0515     --------    -5.020493          1
MinModeTranslate:   13  00:11:47  -239255.125296       0.0515     0.037861    -5.020493          1
DIM:ROT:      14         1   -5.3158   17.5391    3.9117
MinModeTranslate:   14  00:17:28  -239255.126049       0.0316     --------    -5.315827          1

代码

文本

作为对比，我们可以看看NEB的计算结果

代码

文本

[12]

cat ../neb/running_neb.out

Loading mkl version 2023.0.0
Loading tbb version 2021.8.0
Loading compiler-rt version 2023.0.0
Loading mpi version 2021.8.0
Loading compiler version 2023.0.0
Loading oclfpga version 2023.0.0
  Load "debugger" to debug DPC++ applications with the gdb-oneapi debugger.
  Load "dpl" for additional DPC++ APIs: https://github.com/oneapi-src/oneDPL

Loading abacus/3.4.1-dev-icx
  Loading requirement: tbb/latest compiler-rt/latest mkl/latest mpi/latest
    oclfpga/latest compiler/latest
NMAX is 8
 ===== Make Initial NEB Guess =====
--- Successfully make guessed image chain by idpp method ! ---
===== Running NEB =====
----- Running ASE-NEB -----
----- improvedtangent method is being used -----
----- Parallel calculation is being used -----
      Step     Time          Energy          fmax
FIRE:    0 22:22:56  -239253.377660         4.103469
FIRE:    1 22:25:08  -239253.437260         4.096040
FIRE:    2 22:27:25  -239253.561360         3.714696
FIRE:    3 22:29:34  -239253.708860         3.112537
FIRE:    4 22:31:50  -239253.887400         4.614592
FIRE:    5 22:34:04  -239254.094930         6.213053
FIRE:    6 22:36:18  -239254.325040         8.318300
FIRE:    7 22:38:35  -239254.534350         2.047596
FIRE:    8 22:40:46  -239254.713670         5.559694
FIRE:    9 22:42:59  -239254.710950         4.785889
FIRE:   10 22:45:15  -239254.706880         3.275704
FIRE:   11 22:47:27  -239254.704080         1.582266
FIRE:   12 22:49:37  -239254.706430         1.412709
FIRE:   13 22:51:47  -239254.718600         1.379108
FIRE:   14 22:54:04  -239254.719640         1.291100
FIRE:   15 22:56:14  -239254.721750         1.207613
FIRE:   16 22:58:27  -239254.724890         1.193499
FIRE:   17 23:00:39  -239254.729130         1.174229
FIRE:   18 23:02:52  -239254.734450         1.150686
FIRE:   19 23:05:05  -239254.740880         1.122713
FIRE:   20 23:07:16  -239254.748460         1.090173
FIRE:   21 23:09:29  -239254.758150         1.051320
FIRE:   22 23:11:43  -239254.770370         1.003003
FIRE:   23 23:13:54  -239254.785520         0.945088
FIRE:   24 23:16:07  -239254.804180         0.881546
FIRE:   25 23:18:20  -239254.826770         0.827074
FIRE:   26 23:20:34  -239254.853820         0.894126
FIRE:   27 23:22:45  -239254.885600         0.685146
FIRE:   28 23:24:56  -239254.921440         0.610309
FIRE:   29 23:27:10  -239254.959180         0.515415
FIRE:   30 23:29:21  -239254.994930         0.626494
FIRE:   31 23:31:33  -239255.024290         0.632831
FIRE:   32 23:33:43  -239255.044510         0.494768
FIRE:   33 23:35:53  -239255.056240         0.483825
FIRE:   34 23:38:04  -239255.063560         0.673756
FIRE:   35 23:40:17  -239255.070100         0.364123
FIRE:   36 23:42:27  -239255.078380         0.592516
FIRE:   37 23:44:38  -239255.086920         0.365760
FIRE:   38 23:46:49  -239255.093550         0.656566
FIRE:   39 23:49:01  -239255.094940         0.402994
FIRE:   40 23:51:10  -239255.097590         0.306738
FIRE:   41 23:53:21  -239255.101240         0.311911
FIRE:   42 23:55:31  -239255.105460         0.455725
FIRE:   43 23:57:43  -239255.109870         0.416953
FIRE:   44 23:59:53  -239255.114190         0.213299
FIRE:   45 00:02:03  -239255.118500         0.136535
FIRE:   46 00:04:15  -239255.123190         0.359869
FIRE:   47 00:06:28  -239255.128030         0.296227
FIRE:   48 00:08:38  -239255.132870         0.124782
FIRE:   49 00:10:48  -239255.137380         0.328013
FIRE:   50 00:12:58  -239255.140920         0.198620
FIRE:   51 00:15:06  -239255.142930         0.269909
FIRE:   52 00:17:18  -239255.143150         0.175201
FIRE:   53 00:19:28  -239255.141680         0.338259
FIRE:   54 00:21:38  -239255.139480         0.131160
FIRE:   55 00:23:51  -239255.138720         0.406918
FIRE:   56 00:25:59  -239255.139240         0.169032
FIRE:   57 00:28:09  -239255.140240         0.147014
FIRE:   58 00:30:19  -239255.141600         0.313577
FIRE:   59 00:32:31  -239255.143040         0.284662
FIRE:   60 00:34:44  -239255.144490         0.086161
FIRE:   61 00:36:57  -239255.145820         0.189046
FIRE:   62 00:39:06  -239255.146930         0.273623
FIRE:   63 00:41:16  -239255.147780         0.094373
FIRE:   64 00:43:25  -239255.148290         0.205654
FIRE:   65 00:45:36  -239255.148380         0.228877
FIRE:   66 00:47:45  -239255.147940         0.116809
FIRE:   67 00:50:02  -239255.147960         0.078823
FIRE:   68 00:52:13  -239255.148000         0.045560
===== NEB Process Done ! =====
===== Running Post-Processing =====
=== n_max set to 0, automatically detect the images of chain by NEBTools ===
Number of images per band guessed to be 10.
num: 0; Energy: -239256.26878 (eV)
num: 1; Energy: -239256.24673 (eV)
num: 2; Energy: -239256.15888 (eV)
num: 3; Energy: -239255.88496 (eV)
num: 4; Energy: -239255.148 (eV)
num: 5; Energy: -239255.2698 (eV)
num: 6; Energy: -239255.50369 (eV)
num: 7; Energy: -239255.67735 (eV)
num: 8; Energy: -239255.83488 (eV)
num: 9; Energy: -239255.88577 (eV)
Reaction Barrier and Energy Difference: (1.1207800000265422, 0.3830100000195671) (eV)
Number of images per band guessed to be 10.
Processing band         68 /         69
Appears to be only one band in the images.
Processing band          0 /          1
===== Done ! =====

代码

文本

可以发现NEB得到的过渡态绝对能量为 -239255.148 eV，和Dimer结果 -239255.126049 eV是相匹配的。且Dimer计算所得能变 0.38 eV和能垒 1.10 eV （以NEB初末态为参考）也与H2在Au111表面解离过程的已有研究相匹配：参考：Hammer, Norskov. Nature, 1995

文献中，H2在不同金属的111表面解离的最小反应路径计算结果，此处采用PW91泛函

代码

文本

Dimer过程中的优化轨迹run_dimer.traj可以通过ASE可视化，所得过渡态结果的可视化同理。在Bohrium Notebook上进行可视化之前，需要使用ngl可视化器，并pip安装这一可视化器

代码

文本

[3]

%pip install nglview # magic pip is recommended for NB

已隐藏输出

代码

文本

[5]

# in dimer directory

from ase.io import read

from ase.visualize import view

atoms = read("/opt/ATST-Tools/examples/H2-Au111/dimer/run_dimer.traj", ":")

view(atoms, viewer='ngl')

代码

文本

4. 基于P-RFO的Sella方法使用

代码

文本

4.1 Sella简介与使用

ASE提供了一个通用的原子模拟编程平台，除了ASE内置的NEB和Dimer等过渡态方法以外，近年来也有许多新创的过渡态方法基于ASE实现，可以基于ASE的代码逻辑快速调用。

Sella是由E.Hermes, J.Zador等人开发的，基于P-RFO的过渡态搜索方法。该方法的主要思路如下：

通过迭代对角化（iteratively diagonalization)方法，以及TS-BFGS的更新方式，进行Hessian矩阵逆矩阵 $B$ 的求解。这一部分包括Hessian矩阵初始化的过程，它可以通过一些有限差分过程，或者是读入其他软件计算的Hessian矩阵完成
求得逆矩阵 $B$ 之后，即可对它进行本征方程求解，并应用P-RFO方法进行结构演化，在leftmost eigenvector方向朝梯度方向演化，其余方向朝负梯度方向演化。实际上Sella中应用的是restricted step P-RFO(RS-PRFO)方法，该方法对位移矢量的模进行了扫描，使结构演化更为高效。
针对新结构，以上一步结构对应的 $B$ 为预处理器，进行迭代对角化求解该结构的 $B$ ，重复上述过程直至得到目标过渡态。

代码

文本

在使用Sella库之前，你需要先安装Sella, 你的镜像应该帮你装好了。如果没有，可执行如下命令：

代码

文本

[ ]

%pip install sella

代码

文本

Sella的使用相对简单，在ATST-Tools内也能以独立脚本形式实现：

代码

文本

[6]

cd /opt/ATST-Tools/sella

/opt/ATST-Tools/sella

代码

文本

[7]

neb2sella_abacus.py  sella_IRC.py  sella_submit.sh
neb2sella_submit.sh  sella_run.py

代码

文本

可以看到其中有五个脚本，这四个脚本的功能分别如下：

sella_run.py: 调用ABACUS开展Sella计算的主要执行脚本
sella_submit.sh：将Sella计算主脚本提交到Slurm等服务器队列的队列脚本
sella_IRC.py：针对所得过渡态调用Sella包和ABACUS进行IRC计算的脚本
neb2sella_abacus.py：将NEB和Sella联动，调用ABACUS用于过渡态计算的执行脚本
neb2sella_submit.py：将NEB和Sella联动，调用ABACUS用于过渡态计算的提交脚本

本次教程主要解析如何开展Sella计算及其IRC方法使用。可以先来看看sella_run.py脚本的具体内容

代码

文本

[8]

cat sella_run.py

from ase.calculators.abacus import Abacus, AbacusProfile
from ase.io import read, write
from ase.visualize import view
from sella import Sella, Constraints
import os

ts = read("STRU", format="abacus")
# abacus calculator setting
abacus = "abacus"
mpi=16
omp=4
lib_dir = "/lustre/home/2201110432/example/abacus"
#lib_dir = "/data/home/liuzq/example"
pseudo_dir = f"{lib_dir}/PP"
basis_dir = f"{lib_dir}/ORB"
pp = {
        "H":"H_ONCV_PBE-1.0.upf",
        "C":"C_ONCV_PBE-1.0.upf",
        "O":"O_ONCV_PBE-1.0.upf",
        "Fe":"Fe_ONCV_PBE-1.0.upf",
      }
basis = {
        "H":"H_gga_6au_100Ry_2s1p.orb",
        "C":"C_gga_7au_100Ry_2s2p1d.orb",
        "O":"O_gga_7au_100Ry_2s2p1d.orb",
        "Fe":"Fe_gga_8au_100Ry_4s2p2d1f.orb",
         }
kpts = [3, 1, 2]
parameters = {
    'calculation': 'scf',
    'nspin': 2,
    'xc': 'pbe',
    'ecutwfc': 100,
    'ks_solver': 'genelpa',
    'symmetry': 0,
    'vdw_method': 'none',
    'smearing_method': 'mp',
    'smearing_sigma': 0.002,
    'basis_type': 'lcao',
    'mixing_type': 'broyden',
    'mixing_ndim': 20,
    'scf_thr': 1e-7,
    'scf_nmax': 300,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'cal_force': 1,
    'cal_stress': 1,
    'init_wfc': 'atomic',
    'init_chg': 'atomic',
    'out_stru': 1,
    'out_chg': -1,
    'out_mul': 1,
    'out_wfc_lcao': 0,
    'out_bandgap': 0,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
}

# developers only
sella_eta = 0.005 # 0.002 ~ 0.01

# set calculator
def set_abacus_calc(abacus, parameters, directory, mpi, omp) -> Abacus:
    """Set Abacus calculators"""
    os.environ['OMP_NUM_THREADS'] = f'{omp}'
    profile = AbacusProfile(f"mpirun -np {mpi} {abacus}")
    out_directory = directory
    calc = Abacus(profile=profile, directory=out_directory,
                **parameters)
    return calc

if __name__ == "__main__":

    ts.calc = set_abacus_calc(abacus, parameters, f"ABACUS", mpi, omp)

    # run from TS stru, but fmax more required
    # set cons is optional
    #cons = Constraints(ts_neb)
    #cons.fix_translation(ts_neb._get_constraints()[0].get_indices())
    dyn = Sella(
        ts,
        trajectory='run_sella.traj',
        eta = sella_eta,
    )
    dyn.run(fmax=0.05)

    # output TS stru file
    write("TS_sella.stru", ts, format="abacus")
    write("TS_sella.cif", ts, format="cif")

代码

文本

同样的，前半段为计算设置，后半段为调用Sella库开展计算的工作流。由于Sella本身的封装已经非常好，可以像调用一个普通结构优化器一样直接使用，故这一脚本并不依赖source文件夹内脚本，是一个可直接使用的独立脚本。

这里有一个值得注意的参数是sella_eta，它决定了Sella在通过有限差分方法生成初始Hessian矩阵时的步长 $η$ （类似于Dimer长度），该值需要是一个小量，但不能过小，因为在有限差分过程中过小的差分步长会放大计算器本身的数值误差，导致有限差分结果出现问题。在使用力场作为势能面计算器时这一问题尚不显著，但使用第一性原理方法作为势能面计算器时这一问题将变得明显。对此感兴趣的朋友可以参考如下链接：

Sella这一默认 $η$ 值为1e-4（即0.0001），这一值是非常小的，适用于传统力场的使用，但并不适用于调用ABACUS等DFT计算器。经本人测试（Notes: https://github.com/QuantumMisaka/ATST-Tools/issues/17）：

在使用ABACUS通过Sella开展表面反应过渡态搜索时，推荐的 $η$ 值为5e-3，
针对特定情况，2e-3也会有更好的表现。
可选择的 $η$ 区间为1e-2 $\to$ 1e-3
有时 $η$ 到1e-4量级也能得到结果，但计算效率明显降低
对于DP势（以DPA2势为例），将 $η$ 适当调大到2e-4能提升Sella计算性能

代码

文本

IRC（Internal Reaction Coordinate，内秉反应坐标）为从过渡态出发，连接反应初末态产物的（质量加权）最小能量路径。IRC的计算一般通过从过渡态位置出发，沿虚频方向进行固定步长扫描与限制性优化，直至体系能量不再下降来完成。

通过计算IRC，我们可以从过渡态出发确定单个基元反应通道的初末态，并确定当前过渡态确实是我们要找的过渡态。我们再来看看在ATST-Tools中如何基于Sella展开IRC计算。

代码

文本

[9]

cat sella_IRC.py

from ase.calculators.abacus import Abacus, AbacusProfile
from ase.io import read, write, Trajectory
from ase.visualize import view
from sella import IRC
import os

ts_opt = read("STRU", format="abacus")
# abacus calculator setting
abacus = "abacus"
mpi=16
omp=4
lib_dir = "/lustre/home/2201110432/example/abacus"
#lib_dir = "/data/home/liuzq/example"
pseudo_dir = f"{lib_dir}/PP"
basis_dir = f"{lib_dir}/ORB"
pp = {
        "H":"H_ONCV_PBE-1.0.upf",
        "C":"C_ONCV_PBE-1.0.upf",
        "O":"O_ONCV_PBE-1.0.upf",
        "Fe":"Fe_ONCV_PBE-1.0.upf",
      }
basis = {
        "H":"H_gga_6au_100Ry_2s1p.orb",
        "C":"C_gga_7au_100Ry_2s2p1d.orb",
        "O":"O_gga_7au_100Ry_2s2p1d.orb",
        "Fe":"Fe_gga_8au_100Ry_4s2p2d1f.orb",
         }
kpts = [3, 1, 2]
irc_log = "irc_log.traj"
dx = 0.1
fmax = 0.05
steps = 1000
parameters = {
    'calculation': 'scf',
    'nspin': 2,
    'xc': 'pbe',
    'ecutwfc': 100,
    'ks_solver': 'genelpa',
    'symmetry': 0,
    'vdw_method': 'none',
    'smearing_method': 'mp',
    'smearing_sigma': 0.002,
    'basis_type': 'lcao',
    'mixing_type': 'broyden',
    'mixing_ndim': 20,
    'scf_thr': 1e-7,
    'scf_nmax': 300,
    'kpts': kpts,
    'pp': pp,
    'basis': basis,
    'pseudo_dir': pseudo_dir,
    'basis_dir': basis_dir,
    'cal_force': 1,
    'cal_stress': 0,
    'init_wfc': 'atomic',
    'init_chg': 'atomic',
    'out_stru': 1,
    'out_chg': -1,
    'out_mul': 1,
    'out_wfc_lcao': 0,
    'out_bandgap': 0,
    'efield_flag': 1,
    'dip_cor_flag': 1,
    'efield_dir': 1,
}

# developers only
sella_eta = 0.002


# set calculator
def set_abacus_calc(abacus, parameters, directory, mpi, omp) -> Abacus:
    """Set Abacus calculators"""
    os.environ['OMP_NUM_THREADS'] = f'{omp}'
    profile = AbacusProfile(f"mpirun -np {mpi} {abacus}")
    out_directory = directory
    calc = Abacus(profile=profile, directory=out_directory,
                **parameters)
    return calc

if __name__ == "__main__":
    ts_opt.calc = set_abacus_calc(abacus, parameters, f"ABACUS", mpi, omp)
    # run from TS stru, but fmax more required
    # set cons is optional
    #cons = Constraints(ts_neb)
    #cons.fix_translation(ts_neb._get_constraints()[0].get_indices())
    irc_traj = Trajectory(irc_log, 'w')
    irc = IRC(ts_opt, trajectory=irc_traj, dx=dx, eta=sella_eta)
    irc.run(fmax, steps=steps, direction='forward')
    irc.run(fmax, steps=steps, direction='reverse')
    # normalize the trajectory
    irc_log_norm = []
    irc_origin = read(irc_log,":")
    ene_last = -99999999
    for stru in irc_origin[::-1]:
        if stru.get_potential_energy() > ene_last:
            irc_log_norm.append(stru)
            ene_last = stru.get_potential_energy()
        else:
            break
    ene_last = 99999999
    for stru in irc_origin:
        if stru.get_potential_energy() < ene_last:
            irc_log_norm.append(stru)
            ene_last = stru.get_potential_energy()
        else:
            break
    write(f"norm_{irc_log}", irc_log_norm, format="traj")

代码

文本

IRC计算有如下值得注意的参数：

dx 为IRC每步扫描的步长，需要在一个合适值，默认为0.
fmax 为IRC沿特定方向扫描收敛的受力阈值
irc_log 为IRC过程的原始轨迹
steps 为最大IRC计算步数

IRC需要朝初态和末态两端计算，其原始轨迹并不是接在一起的，而是两段从TS出发的轨迹。 ATST-Tools的IRC计算工作流会生成一个normalized的轨迹，这一轨迹会将两段IRC轨迹拼接成一个完整的IRC结果。

了解了相关计算脚本，我们就可以看案例了。由于实际案例运行时间较长，这里建议大家在自己的服务器上调整好相关参数之后运行这些案例。

代码

文本

4.2 案例一：H2在Au111表面的解离及其IRC

代码

文本

同样是这个简单的案例，让我们看看Sella的计算结果

代码

文本

[13]

cd /opt/ATST-Tools/examples/H2-Au111/sella

/opt/ATST-Tools/examples/H2-Au111/sella

代码

文本

[14]

STRU  run_sella.traj  running_sella.out

代码

文本

[15]

cat running_sella.out

     Step     Time          Energy         fmax         cmax       rtrust          rho
Sella   0 21:16:11  -239254.830921       0.6455       0.0000       0.1000       1.0000
Sella   1 21:29:55  -239254.965985       0.3227       0.0000       0.1000       1.0397
Sella   2 21:31:55  -239255.045439       0.2091       0.0000       0.1000       1.0399
Sella   3 21:33:55  -239255.089664       0.1260       0.0000       0.1150       1.0093
Sella   4 21:35:50  -239255.112886       0.1094       0.0000       0.1150       0.9399
Sella   5 21:37:44  -239255.116946       0.0788       0.0000       0.1150       0.6098
Sella   6 21:39:39  -239255.120211       0.0561       0.0000       0.1150       1.2665
Sella   7 21:41:34  -239255.122160       0.0483       0.0000       0.1150       0.9641

代码

文本

输出信息中，rtrust代表优化任务中的信赖半径(trust radius)，而rho代表该点的势能面信息以二阶Taylor展开近似描述的可信度，该值接近于1时，rtrust可以逐步增大以加速结构优化，但该值过大，接近于0或为负时，二次展开的近似描述存在问题，rtrust将随之减小。

这里Sella用了和Dimer同等水平的计算资源（Intel 8358 64core）和ABACUS版本（3.7.5），以及完全一样的初始结构。不难发现Sella得到和Dimer一致的结果，并且Sella计算效率（用时：25 min）相比Dimer计算时间（用时：2h 30min）有质的飞跃。

如果你想在Bohrium Notebook里重复该案例，可以调用16核以上机器，运行如下代码：

代码

文本

[46]

cd /opt/ATST-Tools/examples/H2-Au111/

/opt/ATST-Tools/examples/H2-Au111

代码

文本

[31]

mkdir sella_test && cp data/* sella_test

cp: -r not specified; omitting directory 'data/FINAL'
cp: -r not specified; omitting directory 'data/INIT'

代码

文本

[47]

cd sella_test

/opt/ATST-Tools/examples/H2-Au111/sella_test

代码

文本

[38]

cp ../sella/STRU .

代码

文本

[51]

%%writefile sella_run.py

# modifired sella_run.py for H2-Au111

from ase.calculators.abacus import Abacus, AbacusProfile

from ase.io import read, write

from ase.visualize import view

from sella import Sella, Constraints

import os

ts = read("STRU", format="abacus")

# abacus calculator setting

abacus = "abacus"

mpi = 16

omp = 1

abacus = "abacus"

lib_dir = "."

pseudo_dir = f"{lib_dir}"

basis_dir = f"{lib_dir}"

pp = {

'H':'H_ONCV_PBE-1.0.upf',

'Au':'Au_ONCV_PBE-1.0.upf'

}

basis = {

'H':'H_gga_6au_100Ry_2s1p.orb',

'Au': 'Au_gga_7au_100Ry_4s2p2d1f.orb'

}

kpts = [3, 1, 3]

# INPUT setting

parameters = {

'calculation': 'scf',

'basis_type': 'lcao',

'xc': 'pbe',

'ks_solver': 'genelpa',

'vdw_method': 'd3_bj',

'nspin' : 2,

'ecutwfc': 100,

'pp': pp,

'basis': basis,

'pseudo_dir': pseudo_dir,

'basis_dir': basis_dir,

'smearing_method': 'mp',

'smearing_sigma': 0.002,

'mixing_type': 'broyden',

'mixing_ndim': 20,

'scf_thr': 1e-7,

'scf_nmax': 300,

'kpts': kpts,

'init_wfc': 'atomic',

'init_chg': 'atomic',

'cal_force': 1,

'cal_stress': 0,

'out_stru': 1,

'out_chg': -1,

'out_mul': 1,

'out_bandgap': 1,

'out_wfc_lcao': 0,

'efield_flag': 1,

'dip_cor_flag': 1,

'efield_dir': 1,

}

# developers only

sella_eta = 0.005 # 0.002 ~ 0.01

# set calculator

def set_abacus_calc(abacus, parameters, directory, mpi, omp) -> Abacus:

"""Set Abacus calculators"""

os.environ['OMP_NUM_THREADS'] = f'{omp}'

profile = AbacusProfile(f"mpirun -np {mpi} {abacus}")

out_directory = directory

calc = Abacus(profile=profile, directory=out_directory,

**parameters)

return calc

if __name__ == "__main__":

ts.calc = set_abacus_calc(abacus, parameters, f"ABACUS", mpi, omp)

# run from TS stru, but fmax more required

# set cons is optional

#cons = Constraints(ts_neb)

#cons.fix_translation(ts_neb._get_constraints()[0].get_indices())

dyn = Sella(

ts,

trajectory='run_sella.traj',

eta = sella_eta,

)

dyn.run(fmax=0.05)

# output TS stru file

write("TS_sella.stru", ts, format="abacus")

write("TS_sella.cif", ts, format="cif")

Overwriting sella_run.py

代码

文本

[49]

Au_ONCV_PBE-1.0.upf            H_ONCV_PBE-1.0.upf        STRU
Au_gga_7au_100Ry_4s2p2d1f.orb  H_gga_6au_100Ry_2s1p.orb  sella_run.py

代码

文本

确认计算所需文件之后即可去掉以下注释，运行Sella计算：（或者直接去掉上面的%%writefile标签，运行python程序）。在c32_m64_cpu这一16核机器上相关计算预期需2-3小时完成

代码

文本

[50]

# ! python sella_run.py

代码

文本

Sella的轨迹同样可以可视化，我们不在该案例中重复展示，但IRC轨迹是需要通过该案例展示的

代码

文本

[19]

cd /opt/ATST-Tools/examples/H2-Au111/sella_IRC

/opt/ATST-Tools/examples/H2-Au111/sella_IRC

代码

文本

[17]

STRU          norm_irc_log.traj  running_sella.out  sella_submit.sh
irc_log.traj  running_sella.err  sella_IRC.py

代码

文本

[20]

cat running_sella.out # Sella IRC运行过程

Loading mkl version 2023.0.0
Loading tbb version 2021.8.0
Loading compiler-rt version 2023.0.0
Loading mpi version 2021.8.0
Loading compiler version 2023.0.0
Loading oclfpga version 2023.0.0
  Load "debugger" to debug DPC++ applications with the gdb-oneapi debugger.
  Load "dpl" for additional DPC++ APIs: https://github.com/oneapi-src/oneDPL

Loading abacus/3.7.5-icx
  Loading requirement: tbb/latest compiler-rt/latest mkl/latest mpi/latest
    oclfpga/latest compiler/latest
     Step     Time          Energy          fmax
IRC:    0 11:19:10  -239255.122800        0.048260
IRC:    1 11:25:02  -239255.139780        0.267745
IRC:    2 11:31:00  -239255.184990        0.435348
IRC:    3 11:36:54  -239255.244350        0.483288
IRC:    4 11:40:51  -239255.306200        0.494103
IRC:    5 11:44:50  -239255.365800        0.520505
IRC:    6 11:48:49  -239255.422090        0.531321
IRC:    7 11:54:40  -239255.475620        0.524178
IRC:    8 11:58:34  -239255.526430        0.503247
IRC:    9 12:00:32  -239255.573880        0.469298
IRC:   10 12:04:28  -239255.617880        0.428430
IRC:   11 12:06:26  -239255.659000        0.383776
IRC:   12 12:08:24  -239255.697570        0.372479
IRC:   13 12:12:19  -239255.732720        0.393580
IRC:   14 12:16:14  -239255.762580        0.410280
IRC:   15 12:20:08  -239255.785260        0.420515
IRC:   16 12:26:01  -239255.799380        0.417349
IRC:   17 12:29:56  -239255.806820        0.399448
IRC:   18 12:31:53  -239255.812380        0.376896
IRC:   19 12:33:51  -239255.817500        0.353711
IRC:   20 12:35:48  -239255.822240        0.330837
IRC:   21 12:37:46  -239255.826680        0.308420
IRC:   22 12:39:43  -239255.830820        0.286773
IRC:   23 12:41:41  -239255.834650        0.265522
IRC:   24 12:43:38  -239255.838240        0.245192
IRC:   25 12:45:35  -239255.841580        0.225531
IRC:   26 12:47:33  -239255.844690        0.206864
IRC:   27 12:49:29  -239255.847560        0.188967
IRC:   28 12:51:27  -239255.850230        0.172293
IRC:   29 12:53:23  -239255.852710        0.156395
IRC:   30 12:55:21  -239255.855000        0.142305
IRC:   31 12:57:18  -239255.857130        0.128483
IRC:   32 12:59:15  -239255.859130        0.116161
IRC:   33 13:01:12  -239255.860960        0.104142
IRC:   34 13:03:09  -239255.862670        0.093549
IRC:   35 13:05:06  -239255.864240        0.083213
IRC:   36 13:07:04  -239255.865710        0.074138
IRC:   37 13:09:01  -239255.867040        0.065556
IRC:   38 13:11:01  -239255.868270        0.057697
IRC:   39 13:13:00  -239255.869390        0.050679
IRC:   40 13:14:57  -239255.870430        0.045884
IRC:   41 13:22:42  -239255.140220        0.229534
IRC:   42 13:28:34  -239255.179620        0.350761
IRC:   43 13:34:26  -239255.227100        0.434114
IRC:   44 13:40:32  -239255.278790        0.473188
IRC:   45 13:44:34  -239255.333790        0.495775
IRC:   46 13:48:38  -239255.391270        0.512511
IRC:   47 13:52:42  -239255.450270        0.526034
IRC:   48 13:56:45  -239255.510050        0.537831
IRC:   49 14:00:48  -239255.569870        0.548641
IRC:   50 14:04:50  -239255.629020        0.557362
IRC:   51 14:08:53  -239255.686660        0.561305
IRC:   52 14:10:55  -239255.742000        0.559586
IRC:   53 14:14:58  -239255.794180        0.553824
IRC:   54 14:19:04  -239255.842700        0.547568
IRC:   55 14:21:06  -239255.887380        0.542633
IRC:   56 14:25:08  -239255.928260        0.538278
IRC:   57 14:27:10  -239255.965400        0.532962
IRC:   58 14:29:11  -239255.998670        0.526186
IRC:   59 14:31:13  -239256.028020        0.519117
IRC:   60 14:33:14  -239256.053650        0.512395
IRC:   61 14:37:16  -239256.075940        0.506088
IRC:   62 14:41:20  -239256.095370        0.499515
IRC:   63 14:43:21  -239256.112290        0.491805
IRC:   64 14:47:24  -239256.126900        0.482904
IRC:   65 14:49:30  -239256.139510        0.472734
IRC:   66 14:53:41  -239256.150460        0.461447
IRC:   67 14:55:46  -239256.160090        0.448696
IRC:   68 14:59:56  -239256.168700        0.435023
IRC:   69 15:04:08  -239256.176520        0.420484
IRC:   70 15:08:20  -239256.183760        0.405403
IRC:   71 15:12:31  -239256.190520        0.389726
IRC:   72 15:16:44  -239256.196890        0.373885
IRC:   73 15:20:56  -239256.202960        0.357954
IRC:   74 15:27:16  -239256.208810        0.342676
IRC:   75 15:33:34  -239256.214470        0.327193
IRC:   76 15:42:00  -239256.219910        0.312339
IRC:   77 15:50:23  -239256.225060        0.297482
IRC:   78 15:58:49  -239256.229890        0.282315
IRC:   79 16:05:07  -239256.234360        0.266492
IRC:   80 16:11:26  -239256.238560        0.251117
IRC:   81 16:17:45  -239256.242420        0.235966
IRC:   82 16:21:58  -239256.245970        0.220192
IRC:   83 16:30:23  -239256.249250        0.205184
IRC:   84 16:34:35  -239256.252230        0.189905
IRC:   85 16:38:47  -239256.254960        0.174869
IRC:   86 16:43:00  -239256.257470        0.160229
IRC:   87 16:45:07  -239256.259770        0.145797
IRC:   88 16:47:14  -239256.261900        0.132037
IRC:   89 16:49:20  -239256.263810        0.118672
IRC:   90 16:55:39  -239256.265560        0.105969
IRC:   91 16:57:43  -239256.267120        0.093700
IRC:   92 16:59:49  -239256.268530        0.082637
IRC:   93 17:01:57  -239256.269820        0.073134
IRC:   94 17:04:03  -239256.270990        0.064534
IRC:   95 17:06:07  -239256.272020        0.056720
IRC:   96 17:08:13  -239256.272920        0.049736
IRC:   97 17:10:18  -239256.273670        0.043403
IRC:   98 17:12:25  -239256.274330        0.037683
IRC:   99 17:14:30  -239256.274900        0.032519
IRC:  100 17:16:36  -239256.275390        0.027619
IRC:  101 17:18:43  -239256.275770        0.023177
IRC:  102 17:20:51  -239256.276110        0.018749
IRC:  103 17:22:57  -239256.276370        0.015436
IRC:  104 17:25:03  -239256.276610        0.010885
IRC:  105 17:27:08  -239256.276790        0.007886
IRC:  106 17:29:15  -239256.276940        0.005296
IRC:  107 17:31:21  -239256.277050        0.004214
IRC:  108 17:33:27  -239256.277120        0.003564
IRC:  109 17:35:36  -239256.277150        0.003929
IRC:  110 17:41:58  -239256.277170        0.003843
===== Sella Calculation Done at Sun Sep 22 17:42:03 CST 2024! =====

代码

文本

IRC的轨迹也可以通过ASE可视化出来,这里我们直接看normalization之后的轨迹。

H2在Au111表面解离的IRC曲线，横轴为计算步数，纵轴为能量，顶点代表解离过渡态，左端为解离初态，右端为解离末态

代码

文本

[22]

# 通过ase -T gui norm_irc_log.traj或如下方法可视化，在NB中用ngl可视化器

from ase.io import read

from ase.visualize import view

atoms = read("/opt/ATST-Tools/examples/H2-Au111/sella_IRC/norm_irc_log.traj", ":")

view(atoms, viewer='ngl')

代码

文本

4.3 案例二：环己烷在Pt1@graphene上的解离

代码

文本

这一案例曾在ASE-ABACUS教程第二章的NEB教学中有详细讨论。这里我们仅展示Sella计算结果，并说明Sella优化器可以收敛到正确的过渡态

代码

文本

[27]

cd /opt/ATST-Tools/examples/Cy-Pt@graphene/sella

/opt/ATST-Tools/examples/Cy-Pt@graphene/sella

代码

文本

[28]

STRU  TS_sella.stru  run_sella.traj  running_sella.out

代码

文本

[29]

cat running_sella.out

     Step     Time          Energy         fmax         cmax       rtrust          rho
Sella   0 02:07:14   -11865.114363       0.8749       0.0000       0.1000       1.0000
Sella   1 02:26:15   -11865.241812       0.3090       0.0000       0.1000       0.9207
Sella   2 02:27:36   -11865.282593       0.4343       0.0000       0.1000       0.7929
Sella   3 02:29:00   -11865.325022       0.2571       0.0000       0.1000       0.8962
Sella   4 02:30:23   -11865.350766       0.4631       0.0000       0.1000       0.7874
Sella   5 02:31:15   -11865.344625       0.5155       0.0000       0.0317      -0.6285
Sella   6 02:32:52   -11865.351992       0.4045       0.0000       0.0317       0.8885
Sella   7 02:34:13   -11865.361207       0.2864       0.0000       0.0365       0.9978
Sella   8 02:35:33   -11865.378145       0.3228       0.0000       0.0365       0.7531
Sella   9 02:36:56   -11865.380040       0.2679       0.0000       0.0365       0.3669
Sella  10 02:38:09   -11865.386524       0.2469       0.0000       0.0365       0.4448
Sella  11 02:39:28   -11865.404985       0.2301       0.0000       0.0365       1.2390
Sella  12 02:40:47   -11865.406920       0.4033       0.0000       0.0237       0.1032
Sella  13 02:42:03   -11865.426712       0.1402       0.0000       0.0237       1.0519
Sella  14 02:43:28   -11865.435073       0.2045       0.0000       0.0237       0.7395
Sella  15 02:44:49   -11865.447161       0.1197       0.0000       0.0237       1.1546
Sella  16 02:46:16   -11865.455603       0.2031       0.0000       0.0237       0.7800
Sella  17 02:47:43   -11865.466259       0.1208       0.0000       0.0237       1.0911
Sella  18 02:48:59   -11865.474991       0.1603       0.0000       0.0237       0.8265
Sella  19 02:50:20   -11865.482908       0.1299       0.0000       0.0237       0.9513
Sella  20 02:51:46   -11865.487931       0.1006       0.0000       0.0273       0.9944
Sella  21 02:53:12   -11865.493055       0.1492       0.0000       0.0314       0.9716
Sella  22 02:54:25   -11865.498289       0.1130       0.0000       0.0314       0.8624
Sella  23 02:55:17   -11865.501364       0.1056       0.0000       0.0314       0.8332
Sella  24 02:56:05   -11865.504961       0.0733       0.0000       0.0314       1.3274
Sella  25 02:56:54   -11865.507132       0.2218       0.0000       0.0314       0.3862
Sella  26 02:57:42   -11865.511491       0.0789       0.0000       0.0314       1.1675
Sella  27 02:58:35   -11865.513616       0.1197       0.0000       0.0314       0.7995
Sella  28 02:59:28   -11865.516915       0.0791       0.0000       0.0314       1.2920
Sella  29 03:00:45   -11865.520128       0.1844       0.0000       0.0314       0.6938
Sella  30 03:01:37   -11865.522811       0.0870       0.0000       0.0314       1.0552
Sella  31 03:02:27   -11865.524934       0.0939       0.0000       0.0314       0.8749
Sella  32 03:03:19   -11865.527394       0.0804       0.0000       0.0314       1.2669
Sella  33 03:04:12   -11865.529410       0.0900       0.0000       0.0314       0.6949
Sella  34 03:05:04   -11865.530798       0.0836       0.0000       0.0314       0.8858
Sella  35 03:05:53   -11865.532060       0.0518       0.0000       0.0314       1.1885
Sella  36 03:06:42   -11865.533617       0.0661       0.0000       0.0314       1.0092
Sella  37 03:07:28   -11865.535115       0.0754       0.0000       0.0314       1.0745
Sella  38 03:08:18   -11865.536663       0.0600       0.0000       0.0314       1.0550
Sella  39 03:09:10   -11865.538748       0.0737       0.0000       0.0314       1.1875
Sella  40 03:09:58   -11865.540976       0.0937       0.0000       0.0314       0.8653
Sella  41 03:10:49   -11865.543111       0.0796       0.0000       0.0314       1.1595
Sella  42 03:11:40   -11865.545712       0.0941       0.0000       0.0314       1.1118
Sella  43 03:12:31   -11865.548326       0.1476       0.0000       0.0314       0.9011
Sella  44 03:13:21   -11865.550817       0.0875       0.0000       0.0314       1.0633
Sella  45 03:14:10   -11865.552783       0.0919       0.0000       0.0314       0.8785
Sella  46 03:14:57   -11865.554390       0.0773       0.0000       0.0314       0.9764
Sella  47 03:15:45   -11865.556025       0.0674       0.0000       0.0314       1.1775
Sella  48 03:16:36   -11865.557428       0.0592       0.0000       0.0314       0.9154
Sella  49 03:17:25   -11865.558034       0.0583       0.0000       0.0314       0.6801
Sella  50 03:18:17   -11865.558644       0.0462       0.0000       0.0314       1.0556

代码

文本

所得能量，结构与上期NEB讲解中AutoNEB结果一致。可以可视化观看

代码

文本

[42]

from ase.io import read

from ase.visualize import view

# 可视化过渡态结构

atoms = read("/opt/ATST-Tools/examples/Cy-Pt@graphene/sella/run_sella.traj") # 加一个":"选项即可读入和可视化轨迹

view(atoms, viewer='ngl')

代码

文本

4.4 案例三：Au-Pd掺杂ZSM5分子筛内的甲烷活化

代码

文本

这一案例旨在体现ABACUS针对大体系计算的高效，以及Sella方法对大体系催化问题的过渡态计算的可行性。

以ZSM-5为代表的沸石分子筛体系具有多孔结构，是非常优秀的催化材料，如果在其中加入适量的金属掺杂，其催化活性将进一步提高。

以下计算考察双原子Au-Pd掺杂的ZSM5分子筛在过氧化氢/水存在下中活化甲烷的过程。这一体系具有近三百个原子，其第一性原理计算和过渡态搜索都相对耗时，而Sella与ABACUS的结合显著加速了这一体系的过渡态搜索过程。

代码

文本

[52]

cd /opt/ATST-Tools/examples/CH4-HOAuHOPd-ZSM5

/opt/ATST-Tools/examples/CH4-HOAuHOPd-ZSM5

代码

文本

[53]

TS-CH4-HOAu/  react-info/

代码

文本

在react-info里面可以查看本次计算所得的各吸附态和过渡态结构与能量，它们以ASE等软件广泛支持的extended-xyz(extxyz)格式存储

代码

文本

[55]

cd react-info

/opt/ATST-Tools/examples/CH4-HOAuHOPd-ZSM5/react-info

代码

文本

[56]

CH3-Au-HOPd-H2O-ZSM5.extxyz  CH4-HOAuHOPd-ZSM5.extxyz
CH3-H2O-AuHOPd-ZSM5.extxyz   TS-CH4-HOAuHOPd-ZSM5.extxyz

代码

文本

我们可以可视化看看这些结构长啥样，以及有多大

代码

文本

[61]

from ase.io import read

from ase.visualize import view

atoms = read("CH4-HOAuHOPd-ZSM5.extxyz",)

view(atoms, viewer='ngl')

代码

文本

[63]

from ase.io import read

from ase.visualize import view

atoms = read("TS-CH4-HOAuHOPd-ZSM5.extxyz",)

view(atoms, viewer='ngl')

代码

文本

然后看看计算结果

代码

文本

[58]

! grep energy *.extxyz # check head message from extxyz

CH3-Au-HOPd-H2O-ZSM5.extxyz:Lattice="20.531 0.0 0.0 0.0 20.3034 0.0 0.0 0.0 13.6297" Properties=species:S:1:pos:R:3:momenta:R:3:forces:R:3 energy=-105977.52643714681 free_energy=-105977.52643714681 stress="0.005796891191598145 0.0006643679896967325 -0.0007623203718265117 0.0006643679896967325 0.00690745729498539 -0.0009217585102837123 -0.0007623203718265117 -0.0009217585102837123 0.005693940242833433" magmom=1.0 pbc="T T T"
CH3-H2O-AuHOPd-ZSM5.extxyz:Lattice="20.531 0.0 0.0 0.0 20.3034 0.0 0.0 0.0 13.6297" Properties=species:S:1:pos:R:3:momenta:R:3:forces:R:3 energy=-105975.75705447442 free_energy=-105975.75705447442 stress="0.004959613838661751 0.0005083688559837459 -0.0004968469567996332 0.0005083688559837459 0.0063455661804149585 -0.0006640122737958755 -0.0004968469567996332 -0.0006640122737958755 0.005581196517856466" magmom=1.0 pbc="T T T"
CH4-HOAuHOPd-ZSM5.extxyz:Lattice="20.531 0.0 0.0 0.0 20.3034 0.0 0.0 0.0 13.6297" Properties=species:S:1:pos:R:3:momenta:R:3:forces:R:3 energy=-105975.99869130056 free_energy=-105975.99869130056 stress="0.004790770014773642 0.000322342410627691 -0.0005387575034449036 0.000322342410627691 0.006680892664415042 -0.0005854518379076326 -0.0005387575034449036 -0.0005854518379076326 0.00531209580760848" magmom=1.0 pbc="T T T"
TS-CH4-HOAuHOPd-ZSM5.extxyz:Lattice="20.531028632838 0.0 0.0 0.0 20.30339864799589 0.0 0.0 0.0 13.629729092395806" Properties=species:S:1:pos:R:3:initial_magmoms:R:1:forces:R:3 energy=-105975.5048 stress="0.004901927855855528 0.00047342690003881 -0.00038130596918519214 0.00047342690003881 0.006487933980468179 -0.0005988968718198388 -0.00038130596918519214 -0.0005988968718198388 0.005511617366686268" magmom=1.0 free_energy=-105975.5048 pbc="T T T"

代码

文本

可以发现这四个结构经由ABACUS优化计算所得的绝对能量：

CH3-Au-HOPd-H2O-ZSM5: -105977.5264 eV （CH4解离且CH3吸附在Au上的末态，稳态）
CH3-H2O-AuHOPd-ZSM5: -105975.7571 eV （CH4解离之后CH3自由基以氢键连接到活性位点的末态，亚稳态）
CH4-HOAuHOPd-ZSM5: -105975.9987 eV （初态）
TS-CH4-HOAuHOPd-ZSM5: -105975.5048 eV （过渡态）

可得

初态能垒：0.49 eV
到亚稳态的能变：0.24 eV
到稳态的能变：-1.53 eV

它们与VASP计算所得的结果是可以相匹配的。

代码

文本

相关的过渡态计算过程可以看这里

代码

文本

[64]

cd /opt/ATST-Tools/examples/CH4-HOAuHOPd-ZSM5/TS-CH4-HOAu

/opt/ATST-Tools/examples/CH4-HOAuHOPd-ZSM5/TS-CH4-HOAu

代码

文本

[65]

JobRun.state  TS_sella.stru  run_sella.traj     sella_run.py
STRU          abacus.slurm   running_sella.err  sella_submit.sh
TS_sella.cif  ase_sort.dat   running_sella.out

代码

文本

[67]

cat sella_submit.sh # 查看计算资源

#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=64
#SBATCH --cpus-per-task=1
#SBATCH -J Sella-CH4HOAuPd
#SBATCH -o running_sella.out
#SBATCH -e running_sella.err
#SBATCH -p C064M0256G
#SBATCH --qos=low

# JamesMisaka in 2023-11-14
# workflow of ase-abacus-sella method

# in developer's PKU-WM2 server
source /lustre/home/2201110432/apps/miniconda3/bin/activate
conda activate ase
module load abacus/3.7.5-icx

# if one just done neb calculation and done neb_post.py
# python neb2dimer.py neb_latest.traj # or neb.traj

# Job state 
echo $SLURM_JOB_ID > JobRun.state
echo "Start at $(date)" >> JobRun.state

# just do dimer by run !
# do check your ntasks(mpi) and cpu(omp) setting is right
python sella_run.py

echo "===== Sella Calculation Done at $(date)! ====="

# Job State
echo "End at $(date)" >> JobRun.state

代码

文本

[66]

cat running_sella.out # 查看计算过程

Loading mkl version 2023.0.0
Loading tbb version 2021.8.0
Loading compiler-rt version 2023.0.0
Loading mpi version 2021.8.0
Loading compiler version 2023.0.0
Loading oclfpga version 2023.0.0
  Load "debugger" to debug DPC++ applications with the gdb-oneapi debugger.
  Load "dpl" for additional DPC++ APIs: https://github.com/oneapi-src/oneDPL

Loading abacus/3.7.5-icx
  Loading requirement: tbb/latest compiler-rt/latest mkl/latest mpi/latest
    oclfpga/latest compiler/latest
     Step     Time          Energy         fmax         cmax       rtrust          rho
Sella   0 08:59:29  -105973.795560       0.3720       0.0000       0.1000       1.0000
Sella   1 12:04:26  -105974.915420       0.2707       0.0000       0.1000       0.8563
Sella   2 12:19:47  -105975.009900       0.3921       0.0000       0.1000       0.3381
Sella   3 12:33:49  -105975.066130       0.4098       0.0000       0.1000       0.5912
Sella   4 12:47:30  -105975.054470       0.3545       0.0000       0.0514      -0.1113
Sella   5 13:01:44  -105975.137490       0.4391       0.0000       0.0514       0.8994
Sella   6 13:16:02  -105975.158080       0.3707       0.0000       0.0514       0.3737
Sella   7 13:30:11  -105975.197880       0.1459       0.0000       0.0514       1.2184
Sella   8 13:44:23  -105975.215080       0.5218       0.0000       0.0514       0.4145
Sella   9 13:58:29  -105975.251800       0.1729       0.0000       0.0514       1.2390
Sella  10 14:12:36  -105975.260010       0.4079       0.0000       0.0514       0.2472
Sella  11 14:26:43  -105975.288100       0.1205       0.0000       0.0514       1.1130
Sella  12 14:40:57  -105975.294910       0.1987       0.0000       0.0514       0.3992
Sella  13 14:54:45  -105975.311230       0.0920       0.0000       0.0514       1.2828
Sella  14 15:08:32  -105975.317280       0.2054       0.0000       0.0514       0.3058
Sella  15 15:22:27  -105975.334820       0.1000       0.0000       0.0514       1.1888
Sella  16 15:36:22  -105975.339120       0.1742       0.0000       0.0514       0.3395
Sella  17 15:50:26  -105975.351920       0.0849       0.0000       0.0514       1.3074
Sella  18 16:04:10  -105975.356090       0.1937       0.0000       0.0514       0.3418
Sella  19 16:18:18  -105975.368930       0.0784       0.0000       0.0514       1.2650
Sella  20 16:32:10  -105975.374850       0.1908       0.0000       0.0514       0.6703
Sella  21 16:46:02  -105975.384550       0.1128       0.0000       0.0514       1.1964
Sella  22 16:59:24  -105975.392420       0.1317       0.0000       0.0514       0.7945
Sella  23 17:12:47  -105975.402330       0.1179       0.0000       0.0514       1.2069
Sella  24 17:26:07  -105975.412310       0.1694       0.0000       0.0514       0.6668
Sella  25 17:39:37  -105975.423460       0.1175       0.0000       0.0514       1.1014
Sella  26 17:53:26  -105975.432190       0.1406       0.0000       0.0514       0.7400
Sella  27 18:07:17  -105975.440910       0.1041       0.0000       0.0514       1.0919
Sella  28 18:20:34  -105975.447860       0.1148       0.0000       0.0514       0.7655
Sella  29 18:34:02  -105975.454790       0.1082       0.0000       0.0514       1.0703
Sella  30 18:47:26  -105975.460770       0.1196       0.0000       0.0514       0.9719
Sella  31 19:00:58  -105975.466610       0.1235       0.0000       0.0514       0.8561
Sella  32 19:15:25  -105975.472490       0.0892       0.0000       0.0514       1.2482
Sella  33 19:28:50  -105975.475070       0.2617       0.0000       0.0514       0.2351
Sella  34 19:42:21  -105975.483360       0.0646       0.0000       0.0514       1.0913
Sella  35 19:55:49  -105975.485980       0.0931       0.0000       0.0514       0.6344
Sella  36 20:09:27  -105975.490380       0.0569       0.0000       0.0514       1.3309
Sella  37 20:25:01  -105975.493640       0.1404       0.0000       0.0514       0.3831
Sella  38 20:38:07  -105975.499100       0.0644       0.0000       0.0514       1.0223
Sella  39 20:51:28  -105975.501380       0.0753       0.0000       0.0514       0.6383
Sella  40 21:04:38  -105975.504800       0.0472       0.0000       0.0514       1.3637
===== Sella Calculation Done at Fri Sep 20 21:04:44 CST 2024! =====

代码

文本

这一体系的Sella计算在ABACUS的加持下，仅采用单个node，64core的计算资源，在16小时内计算完毕。

针对这一体系，ABACUS LCAO的单步SCF计算仅需15分钟，如果放在其他PW基组软件上，相关计算耗时将显著增加，如VASP的单步SCF计算耗时将达到4个小时，这再次体现出ABACUS的高效。同时，这一复杂体系的过渡态搜索本身就有一定难度，但在Sella包的加持下，这一任务也变得容易起来。

代码

文本

5 总结

代码

文本

通过本Notebook，你学到了：

过渡态理论与单端过渡态搜索的基础知识，包括Dimer法和P-RFO方法的理论基础
通过ATST-Tools，基于ASE和ASE-ABACUS接口完成Dimer和Sella过渡态计算

不难发现，这类单端方法是可以和NEB这样的双端方法结合起来的，我们可以先跑一个粗略的NEB，再取出其能量最高点作为单端过渡态的初始结构，将过渡态优化到目标阈值。这一做法相比直接采用NEB能显著加速过渡态搜索过程，并使所得的过渡态更加可靠。

当然了，单端过渡态方法还可以用来做势能面上各个位置的过渡态采样，比NEB灵活许多。

如果你对ASE-ABACUS的应用，以及过渡态和过渡态搜索产生了更多兴趣，或者有更多见解，欢迎一起讨论!

代码

文本

ASE

ABACUS

ATST-Tools

催化计算

Sella

Dimer

ASEABACUSATST-Tools催化计算TSSellaDimer

已赞2