ABACUS计算模拟实例 | XII. Pt(111)表面羟基解离的过渡态搜索
©️ Copyright 2023 @ Authors
作者：Jinghang Wang (AISI电子结构团队实习生) 📨
日期：2024-6-21
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：由于本案例需要比较长的计算时间（使用32核需要10h以上），因此本案例不会要求完成具体计算，只涉及预处理和结果分析，故使用c2_m4_cpu即可。如果你想进行具体计算，点击上方的 开始连接 按钮，选择 ABACUS:3.6.3-user-guide 镜像及 c64_m128_cpu 节点配置，并在连接后点击右上方的切换内核，选择Bash内核稍等片刻即可运行。
🎯 本教程旨在帮助初学者学习使用ABACUS及软件包ase-abacus、ATST-Tools完成过渡态初末构型建模，使用AutoNEB方法完成过渡态搜索及后处理，以及对过渡态进行振动分析
- ATST-Tools的使用方法见教程： ATST-Tools （使用ABACUS做过渡态相关计算一定要看！）。
- ATST-Tools官网见：https://github.com/QuantumMisaka/ATST-Tools
本节教程目录如下：
1. 过渡态方法简介
2. Pt(111)表面OH基团构型建模
3.使用ATST-Tools进行AutoNEB计算
4.使用ATST-Tools进行振动计算和自由能较正
在使用ABACUS进行计算时，所有操作都需要在指定的文件夹目录下进行。我们已经准备好了数据集ABACUS-cases，其中文件夹“12”对应本节教程。
autoneb data
可以看到有一共有两个文件夹，data对应Pt(111)表面OH基团构型建模，autoneb对应过渡态及振动频率计算。
1. 过渡态方法及ATST-Tools简介
Arrhenius方程指出，基元反应速率与温度和反应活化能相关。高效催化剂正是通过降低反应活化能来实现净反应速率的极大提升的。在体系由一个稳定态移动至另一个稳定态的动力学变化过程中，存在一条最小能量路径（MEP）。在MEP中，过渡态处的能量最高，反应的活化能等于此能量最高点（过渡态的能量）与反应物能量之差。该能量值对应体系的一个一阶鞍点，沿反应路径方向是极大值，而言其他任何方向均为极小值。该一阶鞍点是沿MEP方向的一个过渡态（transtition state，TS），其振动模式中只有一个频率为虚频。
对于材料老化、演变和化学反应等微观过程，过渡态及其对应的反应势垒的描述非常重要。寻找过渡态的一种常用方法是NEB方法。该方法的主要思路如图1所示，图中的圆圈代表不同的状态坐标（也称为图像）。反应的初始构型initial和最终构型final位于势能面的极小值点上，由一条灰色的线连接。NEB方法通过在其之间插入一系列图像（initial path，初始猜测的反应路径），通过整体优化调整这组图像，使它们沿着MEP方向移动，以在势能面上寻找连接两个相邻极小值的MEP。
不难理解，NEB方法需要插入足够多的中间图像来描述MEP，因此计算成本非常高。此外，NEB方法容易低估势垒。因此，在实际计算中更倾向于使用基于NEB方法发展的攀爬图像微动弹性带（climbing image nudged elastic band,CINEB）方法。CINEB方法考虑了攀爬机制对过渡态定位的影响，通过调整能量最高点的受力，使能量最高点不会因相邻点的弹簧力影响而远离过渡态，同时通过反转该点平行与路径的势能力分量方向，促使能量最高点攀爬至过渡态并且能量升高。
Hammer课题组提出的AutoNEB则结合了IT-NEB（Improved Tangent）和CI-NEB，通过这样的动态NEB过程，新加入的映像将对过渡态和反应路径有更好的解析度，并且显著缓解NEB的计算资源不平衡问题，以及NEB的计算量。本案例使用的就是该方法。
VASP自带NEB方法，使用CINEB方法可以借助Henkelman课题组开发的VTST程序，在这里就不详细赘述。
ABACUS本身并不支持NEB计算，可以通过ase-abacus接口使用ase的过渡态搜索方法。ase-abacus的使用方法见：ase-abacus 。在熟练了ASE和ASE-ABACUS接口的代码逻辑之后，直接写一个python脚本，调用ABACUS等计算组件，基于ASE进行过渡态搜索计算，并不是一件非常困难的事情。但使用ATST-Tools，可以将与过渡态搜索相关的前处理，具体计算与后处理工作被合理模块化与耦合起来，我们只需要进行简单的编辑设置，就可以在服务器等命令行环境下快速开展过渡态搜索工作。ATST-Tools的使用方法见教程： ATST-Tools 。
下面以Pt(111)表面OH基团的解离为例，介绍使用AutoNEB方法搜索过渡态的流程。
2. Pt(111)表面OH基团构型建模
在进行NEB计算前，需要确定两个邻近的能量极小值所对应的构型。通常，这些构型可以通过结构优化获得。对于本案例中的OH解离过程，其初态是OH基团在Pt(111)表面上的一个稳定吸附构型，而末态是H原子和O原子在Pt(111)表面上的一个共吸附构型。对OH基团在Pt(111)表面上的四种吸附位点进行测试，发现它更倾向于结合在两个Pt原子之间的Bri位点（见图2(a)）,因此我们将Bri位点作为OH基团解离的初态构型。Pt(111)位点的不同对H原子吸附能的影响相对较小。因此，在构建本例末态时，可以假定O原子始终占据FCC位点，H原子占据近邻的顶位（见图2(b)），这样的构型可以作为OH基团解离的末态。
本案例建立Pt(111)的slab模型并保留五层原子，在z方向增加12 Å的真空层。
(a) 初态（对应解离前） (b)末态（对应解离后）
FS O_gga_7au_100Ry_2s2p1d.orb H_ONCV_PBE-1.0.upf Pt_ONCV_PBE-1.0.upf H_gga_8au_100Ry_2s1p.orb Pt_gga_8au_100Ry_4s2p2d1f.orb IS cif2STRU.py O_ONCV_PBE-1.0.upf
本案例使用Materials Studio进行建模，将cif文件上传（见IS中的IS.cif以及FS中的FS.cif），并使用cif2STRU.py脚本转为STRU格式文件，最后手动在STRU文件中设置最下面两层原子固定。初态和末态构型的弛豫分别在IS和FS文件夹中进行，这里已经做完了弛豫计算，stdout被重定向到result中用于直接观看。
INPUT_PARAMETERS RUNNING ABACUS-DFT #Parameters (1.General) suffix ABACUS # suffix of OUTPUT DIR nspin 1 # 1/2/4 4 for SOC symmetry 0 # 0/1 1 for open, default # symmetry_autoclose 1 # if symmetry error: set symmetry to 0 # symmetry_prec 1e-5 # default esolver_type ksdft # ksdft, ofdft, sdft, tddft, lj, dp # dft_functional pbe # same as upf file, can be lda/pbe/scan/hf/pbe0/hse ks_solver genelpa #genelpa is default for ksdft-lcao, cg is default for ksdft-pw pseudo_dir ../ orbital_dir ../ #Parameters (2.Iteration) calculation relax # scf relax cell-relax md ecutwfc 80 scf_thr 1e-7 scf_nmax 300 relax_nmax 200 relax_method bfgs # cg, bfgs, cg_bfgs, sd, "fire" force_thr_ev 0.5 # ev # stress_thr 5 #Parameters (3.Basis) basis_type lcao # lcao or pw # kspacing 0.25 # replace KPT # gamma_only 1 # 0/1, replace KPT #Parameters (4.Smearing) smearing_method mp # mp/gaussian/fd/fixed, mp for metal gau for semicon smearing_sigma 0.008 # Rydberg, 0.008 for mp 0.001 for gau #Parameters (5.Mixing) mixing_type broyden # pulay/broyden #mixing_beta 0.2 # for metal: 0.05-0.4, 0 - 0.1 for difficult #mixing_gg0 1.5 # only for transition metal #Parameters (6.Calculation) cal_force 1 cal_stress 1 out_stru 1 # print STRU in OUT out_chg 0 # print CHG or not out_bandgap 0 out_mul 0 # print Mulliken charge and mag of atom in mulliken.txt # restart_save auto # false, auto, other # restart_load false #Parameters (7. Dipole Correction) #efield_flag 1 # open added potential, if 0, all below useless #dip_cor_flag 1 # open dipole correction #efield_dir 2 # direction of dipole correction, 0,1,2 for x,y,z #efield_pos_max 0.0 # max frac-pos of correction , default 0.5. should be in vaccum #efield_pos_dec 0.1 # where the saw-like potential decreases, default 0.1 #efield_amp 0.0 # Amplitude of outer electric field , 0 for only dipole-corr
K_POINTS 0 Gamma 5 5 1 0 0 0
ATOMIC_SPECIES Pt 195.084 Pt_ONCV_PBE-1.0.upf O 15.999 O_ONCV_PBE-1.0.upf H 1.008 H_ONCV_PBE-1.0.upf NUMERICAL_ORBITAL Pt_gga_8au_100Ry_4s2p2d1f.orb O_gga_7au_100Ry_2s2p1d.orb H_gga_8au_100Ry_2s1p.orb LATTICE_CONSTANT 1.889726 LATTICE_VECTORS 5.6158000000 0.0000000000 0.0000000000 -2.807900000 4.8634240000 0.0000000000 0.0000000000 0.0000000000 21.170500000 ATOMIC_POSITIONS Direct Pt 0.0000000000 20 0.0000000000 0.0000000000 0.2027650000 0 0 0 mag 0.0 0.0000000000 0.0000000000 0.5264390000 1 1 1 mag 0.0 0.3333330000 0.1666660000 0.0944710000 0 0 0 mag 0.0 0.3333330000 0.1666660000 0.4171380000 1 1 1 mag 0.0 0.1666660000 0.3333330000 0.3099620000 1 1 1 mag 0.0 0.5000000000 0.0000000000 0.2027650000 0 0 0 mag 0.0 0.5000000000 0.0000000000 0.5264390000 1 1 1 mag 0.0 0.8333330000 0.1666660000 0.0944710000 0 0 0 mag 0.0 0.8333330000 0.1666660000 0.4171380000 1 1 1 mag 0.0 0.6666670000 0.3333330000 0.3099620000 1 1 1 mag 0.0 0.0000000000 0.5000000000 0.2027650000 0 0 0 mag 0.0 0.0000000000 0.5000000000 0.5264390000 1 1 1 mag 0.0 0.3333330000 0.6666670000 0.0944710000 0 0 0 mag 0.0 0.3333330000 0.6666670000 0.4171380000 1 1 1 mag 0.0 0.1666660000 0.8333330000 0.3099620000 1 1 1 mag 0.0 0.5000000000 0.5000000000 0.2027650000 0 0 0 mag 0.0 0.5000000000 0.5000000000 0.5264390000 1 1 1 mag 0.0 0.8333330000 0.6666670000 0.0944710000 0 0 0 mag 0.0 0.8333330000 0.6666670000 0.4171380000 1 1 1 mag 0.0 0.6666670000 0.8333330000 0.3099620000 1 1 1 mag 0.0 O 0.0000000000 1 0.7445320000 0.4927860000 0.6076120000 1 1 1 mag 0.0 H 0.0000000000 1 0.8354300000 0.6741040000 0.6281760000 1 1 1 mag 0.0
Note： ase目前还不支持晶格矢量的初末构型进行neb计算，因此需要保证STRU中LATTICE_VECTORS在初末构型中完全相同，并且INPUT中的calculation参数只能选relax而不能选cell-relax。
3. 使用ATST-Tools进行AutoNEB计算
AutoNEB_iter neb_latest.extxyz run_autoneb000.traj AutoNEBrun-rank0 neb_latest.traj run_autoneb001.traj STRU-final neb_make.py run_autoneb002.traj STRU-init neb_make_pymatgen.py run_autoneb003.traj ase_sort.dat neb_post.py run_autoneb004.traj autoneb_run.py nebplots_all.pdf run_autoneb_init.extxyz init_neb_chain.traj result_nebpost vib_analysis_TS neb_dist.py result_nebrun
ATST-Tools官网 提供了一系列辅助脚本用于NEB计算的预处理，实际计算，后处理等。本镜像已经克隆了其github上的内容到/opt/ATST-Tools下，因此也可以直接在本地获取脚本。
3.1. neb_dist.py确定需要的图像个数
首先，我们可以使用neb_dist.py确定大致需要的图像个数。生成初猜的中间图像是NEB计算的第一步。使用neb_dist.py能够计算初、末态的构型距离，VTST推荐采用“距离值/0.8”取整作为插入图像的数目，这里也可以参考。
==> Distance Between Two Image is 1.9438 Angstrom <==
Note：如果直接使用弛豫完的构型运行上述命令计算距离，很有可能得到过大的距离，这与ase软件本身对周期性的处理有关。例如：对于H原子，(0,0,z_H)/(0,1,z_H)/(1,0,z_H)/(1,1,z_H)这四个相对坐标（弛豫后都可能出现），原本由于周期性边界条件应该是完全等价的，但是ase在计算距离时会将其转换为直角坐标，导致这四种结果不再等价。处在边界上的原子普遍都存在这样的问题。这一问题我们后面还会碰到。
此时一种解决办法是手动微调H原子的位置：在这里我们知道Bri位点的OH基团朝向的Pt原子的坐标是(1,1,z_Pt)，因此可以将H原子的坐标调整为(0.999999,0.999999,z_H)，从而可以得到正确的距离。
3.2.neb_make_pymatgen.py生成初猜图像
初猜图像的生成就是在初末状态之间进行插值，这里我们在初末状态之间插入3个图像。下面的neb_make.py和neb_make_pymatgen.py均默认使用非线性插值：IDPP(image dependent pair potential，相关图像对势)方法。
Reading files: ../data/IS/OUT.ABACUS/running_relax.log and ../data/FS/OUT.ABACUS/running_relax.log Generating path, number of images: 3, sort_tol: 1.0 Optimizing path using IDPP method Fix Atoms below 0.25 in direction z Fix Atoms below 0.25 in direction z Fix Atoms below 0.25 in direction z Writing path: init_neb_chain.traj,Number of images: 5
运行完该命令后会生成init_neb_chain.traj文件，即包含初末状态的初始链。traj这一二进制轨迹格式是ASE官方指定的轨迹存储格式，并能直接被ATST-Tools利用。
Note:
1.还记得我们上面说到的周期性的问题吗？没错，在这里它又出现了，因为ase在这里生成的traj文件时又会将相对坐标转为直角坐标。这会有什么后果呢？想象一下你的原子本来在原点附近，弛豫后从盒子的前面“穿”到了盒子后面（例如相对坐标在弛豫前是(0.000001,0,0），弛豫后是(-0.000001,0,0)），相对来讲几乎没有动，但是转为直角坐标后就成了(0.999999*a,0,0)，相当于移动了整个x方向基矢长度的距离。用这样的初末态计算NEB将会得到十分离谱的结果。所以我们这里给出了两个命令，前者在这里会给出错误的初猜图像，而后者使用了pymatgen模块避免了这个问题。
2.在这里是否可以使用输出文件的STRU_NOW.cif文件进行插值？答案是：目前还不可以。因为调用ase做NEB计算要求IS和FS必须包含能量信息，这就要求必须使用包含能量的running_relax.log文件而不能只使用结构文件。
3.设置约束。在优化初末态构型时，我们固定了底部两层原子不移动。对于中间图像，我们也应该设置相同的约束，否则初末态构型和中间构型将不在同一个势能面上。实现相同约束的方法是在生成初始图像的命令中使用fix参数，这里--fix 0.25:2的含义是：z轴方向相对坐标小于0.25的原子不移动。
3.3. autoneb_run.py进行AutoNEB计算
有了初始链之后，简单修改一下autoneb_run.py中的核数、赝势轨道文件、abacus输入文件相关参数以及与并行NEB相关参数等即可运行。
from ase.optimize import FIRE, BFGS
from ase.io import read, write
from ase.parallel import world, parprint, paropen
#from pathlib import Path
from abacus_autoneb import AbacusAutoNEB
# setting
mpi = 32
omp = 1
neb_optimizer = FIRE # suited for CI-NEB
neb_directory = "AutoNEBrun"
# algorism = 'eb' # default for AutoNEB
algorism = "improvedtangent" # IT-NEB
init_chain = "init_neb_chain.traj"
climb = True
fmax = [0.20, 0.05] # eV / Ang, 2 fmax for others and last CI-NEB in all-images
n_simul = 1 # only for autoneb, number of simultaneous calculation
n_images = 5 # only for autoneb, max number of all image, which should be reached
smooth_curve = False # True to do more neb step to smooth the curve.
k = 0.10 # eV/Ang^2, force constant of spring, 0.05 is from VTST-Tools
# setting for calculator
abacus = "abacus"
lib_dir = "../data"
pseudo_dir = f"{lib_dir}/"
basis_dir = f"{lib_dir}/"
pp = {
'O':'O_ONCV_PBE-1.0.upf',
'H':'H_ONCV_PBE-1.0.upf',
'Pt':'Pt_ONCV_PBE-1.0.upf',
}
basis = {
'O': 'O_gga_7au_100Ry_2s2p1d.orb',
'H': 'H_gga_8au_100Ry_2s1p.orb',
'Pt': 'Pt_gga_8au_100Ry_4s2p2d1f.orb'
,}
kpts = [5, 5, 1]
parameters = {
'calculation': 'scf',
'nspin': 1,
'xc': 'pbe',
'ecutwfc': 80,
'dft_functional': 'pbe',
'ks_solver': 'genelpa',
'symmetry': 0,
'smearing_method': 'gaussian',
'smearing_sigma': 0.0037,
'basis_type': 'lcao',
'mixing_type': 'broyden',
'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': 0,
'out_chg': 0,
}
if __name__ == "__main__":
# running process
# read initial guessed neb chain
init_chain = read(init_chain, index=':')
neb = AbacusAutoNEB(init_chain, parameters, algorism=algorism, directory=neb_directory, k=k, n_simul=n_simul, n_max=n_images, abacus=abacus, mpi=mpi, omp=omp)
neb.run(optimizer=neb_optimizer, climb=climb,
fmax=fmax, smooth_curve=smooth_curve)
⚠️这一步使用32核预计耗时10h，因此在这里不做运行，直接给出输出。
这里注意n_image是图像总数,需要大于n_simul+2。
如果这里提示缺少依赖库，可以使用下列命令下载：
上述计算完成后将会生成AutoNEBrun-rank0：NEB计算中自洽计算的临时文件，每做一步自洽都会刷新，因此最终能看到的就是最后一步的自洽计算（如果你使用了mpi运行gpaw python autoneb_run.py，那么这里会有多个该文件夹，每个文件夹内都是NEB计算中自洽计算的临时文件）；AutoNEB_iter：包含NEB计算的每一步迭代时的力、能量等信息；run_autoneb001.traj至run_autoneb004.traj：NEB计算的每一步迭代对应的轨迹文件。命令行输出被重定向到result_nebrun。
Notice: AutoNEB method is set You manually set n_simul = 1, n_max = 5 ----- Running AutoNEB ----- ----- improvedtangent method is being used ----- The NEB initially has 5 images (including the end-points) Start of evaluation of the initial images Now starting iteration 1 on [0, 3, 4] Now starting iteration 2 on [0, 2, 3] => Notice: Remove existing AutoNEBrun-rank0 for next NEB calculation <= Now starting iteration 3 on [0, 1, 2] => Notice: Remove existing AutoNEBrun-rank0 for next NEB calculation <= Finished initialisation phase. n_max images has been reached ****Now doing the CI-NEB calculation**** Now starting iteration 4 on [1, 2, 3] => Notice: Remove existing AutoNEBrun-rank0 for next NEB calculation <= ----- AutoNEB calculation finished -----
Step Time Energy fmax FIRE: 0 23:07:15 -66455.729609 2.533728 FIRE: 1 23:16:45 -66455.853359 2.267671 FIRE: 2 23:27:40 -66456.046975 1.695272 FIRE: 3 23:34:45 -66456.201581 0.786062 FIRE: 4 23:42:50 -66456.180728 1.596629 FIRE: 5 23:51:45 -66456.193086 1.439411 FIRE: 6 23:59:40 -66456.213447 1.153377 FIRE: 7 00:07:00 -66456.235135 0.784760 FIRE: 8 00:15:05 -66456.252045 0.452501 FIRE: 9 00:23:10 -66456.260827 0.381228 FIRE: 10 00:31:24 -66456.262138 0.540499 FIRE: 11 00:40:09 -66456.262424 0.536357 FIRE: 12 00:48:17 -66456.262991 0.528799 FIRE: 13 00:56:59 -66456.263809 0.517518 FIRE: 14 01:06:24 -66456.264858 0.502790 FIRE: 15 01:15:34 -66456.266091 0.484508 FIRE: 16 01:24:55 -66456.267489 0.464306 FIRE: 17 01:35:05 -66456.269004 0.440861 FIRE: 18 01:43:34 -66456.270726 0.413499 FIRE: 19 01:52:58 -66456.272629 0.381645 FIRE: 20 02:02:09 -66456.274696 0.347420 FIRE: 21 02:11:09 -66456.276876 0.311161 FIRE: 22 02:19:45 -66456.279152 0.276430 FIRE: 23 02:27:53 -66456.281503 0.246155 FIRE: 24 02:37:10 -66456.283957 0.230264 FIRE: 25 02:45:24 -66456.286519 0.202627 FIRE: 26 02:55:05 -66456.288909 0.176002
Step Time Energy fmax FIRE: 0 08:50:38 -66455.524118 0.704231 FIRE: 1 09:00:24 -66455.514961 0.712220 FIRE: 2 09:10:09 -66455.497607 0.709480 FIRE: 3 09:20:05 -66455.475275 0.667487 FIRE: 4 09:28:44 -66455.455390 0.558110 FIRE: 5 09:36:40 -66455.446975 0.389782 FIRE: 6 09:44:25 -66455.455565 0.320521 FIRE: 7 09:50:40 -66455.478108 0.532222 FIRE: 8 09:57:10 -66455.478040 0.527056 FIRE: 9 10:03:45 -66455.477888 0.516362 FIRE: 10 10:10:20 -66455.477621 0.500914 FIRE: 11 10:17:55 -66455.477212 0.479598 FIRE: 12 10:25:07 -66455.476660 0.453041 FIRE: 13 10:32:45 -66455.476008 0.420409 FIRE: 14 10:40:35 -66455.475367 0.381852 FIRE: 15 10:48:40 -66455.474836 0.333177 FIRE: 16 10:56:35 -66455.474795 0.272704 FIRE: 17 11:04:50 -66455.475833 0.201243 FIRE: 18 11:14:15 -66455.478437 0.123002 FIRE: 19 11:22:15 -66455.482812 0.137385 FIRE: 20 11:31:40 -66455.488589 0.155020 FIRE: 21 11:39:40 -66455.494688 0.185833 FIRE: 22 11:47:19 -66455.499336 0.231329 FIRE: 23 11:55:27 -66455.500966 0.273072 FIRE: 24 12:02:05 -66455.499777 0.284775 FIRE: 25 12:09:35 -66455.498683 0.246005 FIRE: 26 12:16:55 -66455.500704 0.172494 FIRE: 27 12:24:50 -66455.504309 0.205741 FIRE: 28 12:33:30 -66455.504410 0.189265 FIRE: 29 12:42:40 -66455.504538 0.158799 FIRE: 30 12:50:20 -66455.504584 0.118618 FIRE: 31 12:59:40 -66455.504539 0.080585 FIRE: 32 13:07:50 -66455.504408 0.064948 FIRE: 33 13:15:40 -66455.504332 0.046122
3.4.neb_post.py后处理
最后，ATST-Tools还提供了反应路径可视化脚本neb_post.py。这里无需修改脚本内容，直接运行。
该命令运行完成后将会生成nebplots_all.pdf：反应路径可视化的pdf文件；neb_latest.traj：NEB计算中最后一次迭代结束后的轨迹文件，即MEP；neb_latest.extxyz：neb_latest.traj对应的.extxyz格式的文件，用于可视化。
同样，命令行输出被重定向到result-nebpost：
=== n_max set to 0, automatically detect the images of chain by NEBTools === Appears to be only one band in the images. SinglePointCalculator(energy=-66456.44354942004, forces=...) num: 0; Energy: -66456.44354942004 (eV) SinglePointCalculator(energy=-66456.381104, forces=...) num: 1; Energy: -66456.381104 (eV) SinglePointCalculator(energy=-66455.504332, forces=..., stress=...) num: 2; Energy: -66455.504332 (eV) SinglePointCalculator(energy=-66456.288909, forces=...) num: 3; Energy: -66456.288909 (eV) SinglePointCalculator(energy=-66456.419796, forces=...) num: 4; Energy: -66456.419796 (eV) Reaction Barrier and Energy Difference: (0.9392174200474979, 0.023753420042339712) (eV) Appears to be only one band in the images. Processing band 0 / 1
所生成的nebplots_all.pdf很好地给出了当前的NEB反应路径对应的势能-反应坐标曲线。
4.使用ATST-Tools进行振动计算和自由能较正
为了验证我们得到的过渡态是否是正确的过渡态结构（即该过渡态是否具有沿反应路径的单一虚频），以及考虑除零温电子能量之外的其他统计热力学信息，比如零点能和适当的振动自由能校正，我们还需要进行振动计算。振动计算在vib_analysis_TS目录下进行。
SCF-rank0 result-vib vib.1.traj vib.4.traj vib.7.traj ase_sort.dat vib vib.2.traj vib.5.traj vib.8.traj neb_latest.traj vib.0.traj vib.3.traj vib.6.traj vib_analysis.py
ATST-Tools中集成了使用ASE-ABACUS基于有限差分进行振动计算的方法。该工作流实际上只使用了一个脚本vib_analysis.py，我们可以直接查看：
# JamesMisaka in 2023-11-30
# Vibrational analysis from finite displacement by using abacus
# part of ATST-Tools scripts
import os
import numpy as np
from ase.vibrations import Vibrations
from ase.thermochemistry import HarmonicThermo
from ase.io import read, write
from ase.calculators.abacus import Abacus, AbacusProfile
from ase.parallel import world, parprint
from ase.mep.neb import NEBTools
from neb2vib import neb2vib
# indices setting for which atoms to be displaced
vib_indices = [1,11,6,16,20,21]
# usage by neb
neb_traj = read('neb_latest.traj', index='2')
atoms, vib_indices = neb2vib(neb_traj)
# traditional
# stru = "STRU"
# atoms = read(stru)
# vib_indices = [0, 1, 37]
T = 523.15 # K
vib_name = 'vib'
delta = 0.01
nfree = 2
abacus = "abacus"
mpi = 32
omp = 1
lib_dir = "../data"
pseudo_dir = f"{lib_dir}/"
basis_dir = f"{lib_dir}/"
pp = {
'O':'O_ONCV_PBE-1.0.upf',
'H':'H_ONCV_PBE-1.0.upf',
'Pt':'Pt_ONCV_PBE-1.0.upf',
}
basis = {
'O': 'O_gga_7au_100Ry_2s2p1d.orb',
'H': 'H_gga_8au_100Ry_2s1p.orb',
'Pt': 'Pt_gga_8au_100Ry_4s2p2d1f.orb'
,}
kpts = [5, 5, 1]
parameters = {
'calculation': 'scf',
'nspin': 1,
'xc': 'pbe',
'ecutwfc': 80,
'dft_functional': 'pbe',
'ks_solver': 'genelpa',
'symmetry': 0,
'smearing_method': 'gaussian',
'smearing_sigma': 0.0037,
'basis_type': 'lcao',
'mixing_type': 'broyden',
'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': 'atomic',
'out_stru': 0,
'out_chg': 0,
}
def set_calculator(abacus, parameters, mpi=1, omp=1) -> Abacus:
"""Set Abacus calculators"""
os.environ['OMP_NUM_THREADS'] = f'{omp}'
profile = AbacusProfile(
argv=['mpirun', '-np', f'{mpi}', abacus])
out_directory = f"SCF-rank{world.rank}"
calc = Abacus(profile=profile, directory=out_directory,
**parameters)
return calc
if __name__ == "__main__":
print("==> Starting Vibrational Analysis <==")
atoms.calc = set_calculator(abacus, parameters, mpi=mpi, omp=omp)
vib = Vibrations(atoms, indices=vib_indices,
name=vib_name, delta=delta, nfree=nfree)
print("==> Running Vibrational Analysis <==")
vib.run()
# post-processing
print("==> Done !!! <==")
print(f"==> All force cache will be in {vib_name} directory <==")
print("==> Vibrational Analysis Summary <==")
vib.summary()
print("==> Writing All Mode Trajectory <==")
vib.write_mode()
# thermochemistry
print("==> Doing Harmonmic Thermodynamic Analysis <==")
vib_energies = vib.get_energies()
#real_vib_energies = np.array([energy for energy in vib.get_energies() if energy.imag == 0 and energy.real > 0], dtype=float)
thermo = HarmonicThermo(vib_energies, ignore_imag_modes=True,)
entropy = thermo.get_entropy(T)
free_energy = thermo.get_helmholtz_energy(T)
print(f"==> Entropy: {entropy:.6e} eV/K <==")
print(f"==> Free Energy: {free_energy:.6f} eV <==")
这里我们需要做的修改是核数、选择需要进行有限差分和计算振动模的图像以及其中的原子（ABACUS输入参数部分可以直接复制之前autoneb_run.py里的）。实际上，严格的振动模计算肯定是需要对晶胞内所有原子进行有限差分的，但由于针对表面催化体系，我们实际上只关心那些直接参与反应的部分原子的振动模式，从而我们只需要选择一部分原子即可。这里我们选择第一层Pt原子以及O、H原子进行振动频率计算。
⚠️这一步使用32核预计耗时2h，因此在这里不做运行，直接给出输出。
该命令完成后将会生成SCF-rank0：振动计算中自洽计算的临时文件，每做一步自洽都会刷新，因此最终能看到的就是最后一步的自洽计算；vib.i.traj：每一个振动模式对应的轨迹文件；vib：以JSON格式存储了计算所得的各个不同结构的能量与原子受力，用于程序复用，在下一次输出振动分析结果，或更换温度等条件进行自由能校正计算时不需要重新调用ABACUS进行计算。
同样，命令行输出被重定向到result-vib：
Appears to be only one band in the images.
=== Normalized displacement vector generated by 1 and 3 images of NEB chain ===
[0.0072698 0.14983481 0.0079008 0.04113835 0.02028012 0.00640004
0.09792615 0.00778042 0.04154566 0.01322169 0.00501568 0.03457091
0.0040057 0.07784615 0.01867339 0.00510257 0.03434381 0.01070796
0.01572058 0.0187964 0.66596701 0.71460546]
=== Locate TS in 2 of 0-4 images ===
=== TS Barrier: 0.9355 (eV) ===
=== TS main moving atoms: [1, 20, 21] ===
==> Starting Vibrational Analysis <==
==> Running Vibrational Analysis <==
==> Done !!! <==
==> All force cache will be in vib directory <==
==> Vibrational Analysis Summary <==
---------------------
# meV cm^-1
---------------------
0 95.0i 766.2i
1 11.3 91.1
2 13.5 108.6
3 13.8 111.5
4 38.3 309.3
5 40.5 326.7
6 60.3 486.7
7 82.1 661.8
8 195.4 1576.0
---------------------
Zero-point energy: 0.228 eV
==> Writing All Mode Trajectory <==
==> Doing Harmonmic Thermodynamic Analysis <==
Entropy components at T = 523.15 K:
=================================================
S T*S
S_harm 0.0009045 eV/K 0.473 eV
-------------------------------------------------
S 0.0009045 eV/K 0.473 eV
=================================================
Internal energy components at T = 523.15 K:
===============================
E_pot 0.000 eV
E_ZPE 0.228 eV
Cv_harm (0->T) 0.213 eV
-------------------------------
U 0.441 eV
===============================
Entropy components at T = 523.15 K:
=================================================
S T*S
S_harm 0.0009045 eV/K 0.473 eV
-------------------------------------------------
S 0.0009045 eV/K 0.473 eV
=================================================
Free energy components at T = 523.15 K:
=======================
U 0.441 eV
-T*S -0.473 eV
-----------------------
F -0.032 eV
=======================
==> Entropy: 9.045226e-04 eV/K <==
==> Free Energy: -0.032390 eV <==
振动计算结果表明，只有一个766.2 cm-1的虚频。
通过查看其振动模式，发现振动矢量的方向正是由初态指向末态，表明成功搜索到了OH基团在Pt(111)表面解离的过渡态。
恭喜你完成了第十二个ABACUS计算模拟实例🎉
如果你还想尝试更多的计算实例，可以点击下方了解ABACUS计算模拟合集链接： ABACUS计算模拟实例 | 概述
kirk0830