ABACUS计算模拟实例 | XI. Pt表面简单物种的吸附能计算
©️ Copyright 2023 @ Authors
作者：Jinghang Wang (AISI电子结构团队实习生) 📨
日期：2024-6-20
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：由于本案例需要比较长的计算时间，因此本案例不会要求完成具体计算，只涉及预处理和结果分析，故使用c2_m4_cpu及ABACUS:3.6.3-user-guide 镜像，点击上方的 开始连接 按钮并在连接后点击右上方的切换内核，选择Bash内核稍等片刻即可运行。
🎯 本教程旨在帮助初学者学习使用ABACUS及软件包ase-abacus、ATST-Tools完成过渡态初末构型建模，使用AutoNEB方法完成过渡态搜索及后处理，以及对过渡态进行振动分析
本节教程目录如下：
1. Pt(111)表面简单物种的吸附行为
2. Pt(111)表面吸附位点分析与吸附能提取
3. 吸附能的零点能矫正
在使用ABACUS进行计算时，所有操作都需要在指定的文件夹目录下进行。我们已经准备好了数据集ABACUS-cases，其中文件夹“11”对应本节教程。
ZPE relax
relax和ZPE两个文件夹分别对应本节教程中的吸附能和零点能计算。
1.Pt(111)表面简单物种的吸附行为
在催化反应中，有效的吸附是发生高效催化反应的必要条件。当气体与固体表面接触时，表面上气体浓度高于气相本体浓度的现象即为吸附现象。固体表面上气体浓度随时间增加而增大，则所发生的为吸附过程，反之则为脱附过程。吸附气体的固体物质为吸附剂，被吸附的气体为吸附质，吸附质在吸附剂表面吸附后的状态为吸附态。
根据吸附分子在固体表面吸附时的结合力（例如分子间作用力、静电相互作用力）不同，可将吸附分为物理吸附和化学吸附，其吸附能曲线如下图所示。其中：蓝色曲线代表物理吸附。物理吸附依靠较弱的分子间作用力，通常对吸附分子结构的影响不大。气体分子接近固体表面时，因范德瓦尔斯相互作用而产生引力，导致一个能量极小值——吸附热 出现；红色曲线代表化学吸附。化学吸附能够改变吸附分子的键合状态，使吸附中心和吸附质之间发生电子的重新分配，一般包含电子共享或转移，而非简单的微扰或弱极化作用。
吸附是催化反应中最为重要的环节之一。在实际应用中，固体催化剂是最为常见的催化剂之一，其中应用最广泛的是Pt催化剂。采用Pt催化剂时，简单物种在其表面的吸附行为和反应路径对催化效率具有重要影响。
第一性原理计算中，吸附能Ead(X)常被用于描述吸附作用的强弱，对于确定固体表面反应机理至关重要。吸附能通常由形如下面的公式所定义： 式中：EX为吸附质的能量；Eslab为干净吸附剂表面模型（slab）的能量；Eslab+X为吸附体系的总能量。下面以O原子在Pt（111）表面的吸附为例，介绍吸附能的典型计算流程。
2. Pt(111)表面吸附位点分析与吸附能提取
首先我们建模Pt(111)表面，包含4个原子，之后设置厚度Thickness为5（即5层Pt原子）并设置12Å的真空层，并且手动设置底部两层原子在计算时位置固定（直接在STRU文件中的原子坐标后添加“0 0 0”即可，无须像VASP中设置Selective dynamics）。slab优化后的能量为-66006.58 eV.
对于氧原子，如图2所示，在Pt(111)表面有四种可能的吸附位点：顶位 (Top),桥位 (Bri),和两种间隙位 (FCC和HCP位点)，我们需要测试在这些可能的位点上的吸附能以确定哪种吸附位点最稳定。在建立吸附模型时，需要注意吸附质原子或分子与表面原子之间应该保持适当的距离（可以参考相关键长），否则，在结构优化过程中可能会出现吸附质远离表面等未能有效吸附的情况。
Bri HCP O_ONCV_PBE-1.0.upf Pt_ONCV_PBE-1.0.upf Top FCC O2 O_gga_7au_100Ry_2s2p1d.orb Pt_gga_8au_100Ry_4s2p2d1f.orb relax.sh
输入输出文件如下：
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
Note:
1.所有结构优化的INPUT中必须设置相同的ecutwfc值，否则吸附能没有对比意义。
2.固定的那部分原子只是不参与结构优化，而并非不参与自洽计算，他们的能量对于体系总能量仍然有贡献，需要被计算，所以即使固定了部分原子，不论是这里的弛豫计算还是下面的振动频率计算，仍然需要大量的计算时间。
ATOMIC_SPECIES Pt 195.084 Pt_ONCV_PBE-1.0.upf O 15.999 O_ONCV_PBE-1.0.upf NUMERICAL_ORBITAL Pt_gga_8au_100Ry_4s2p2d1f.orb O_gga_7au_100Ry_2s2p1d.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.6656666667 0.3343333333 0.5750000000 1 1 1 mag 0.0
对于初始构型，我们设置氧原子在间隙位时距离表面Pt原子1Å，在顶位时2Å，在桥位时1.5Å.
K_POINTS 0 Gamma 5 5 1 0 0 0
这里我们直接查看优化后的能量：
./OUT.ABACUS/running_relax.log: !FINAL_ETOT_IS -66440.1328772621200187 eV ../Bri/OUT.ABACUS/running_relax.log: !FINAL_ETOT_IS -66440.1232593575550709 eV ../HCP/OUT.ABACUS/running_relax.log: !FINAL_ETOT_IS -66439.7028961999749299 eV ../Top/OUT.ABACUS/running_relax.log: !FINAL_ETOT_IS -66438.8006687207089271 eV
可以看到能量最低的位点为FCC，是最稳定的构型，其次是Bri、HCP，最后是Top，并且Bri和FCC的能量几乎一样。但实际O原子的Bri位点并不是一个稳定的位点，在优化后会落到附近的间隙位中。所以我们只要看一下优化后的结构文件就可以发现，位于Bri位点的氧原子优化后会落到FCC位点（如图3(a) 所示）。所以在这里（包括后续教程中）我们只考虑其余三个位点即可。
我们还计算了气态O2分子的能量，用其一半作为吸附质O原子的能量。
!FINAL_ETOT_IS -864.9422011353332209 eV
将FCC位点的slab+O，slab以及O的能量带入吸附能公式得到O原子在Pt(111)表面的最稳定吸附能为Eads(O) = -1.08 eV.
3. 吸附能的零点能矫正
以经典力学预测，当温度为绝对零度时分子/原子的动能为0，因此不会在平衡位置附近振动。根据量子力学的谐振子模型则可以发现其振动能具有E = (n+1/2)hv的形式，因此最低振动能级也具有非零的能量贡献。将所有简正振动模加和，得到对电子能量的零点能矫正量。 式中：h为普朗克常量；为谐振子的振动频率。该最低能量E和经典力学最小能量E0之间的差值为零点能（zero-point energy, ZPE）EZPE。若原子合集与其对应的正态振型可分别独立的决定每个振型的零点能，则该原子合集的最小能量可表示为式中：E0为DFT计算得到的能量，即EDFT；i为正态振型的频率。
O2 Oads
我们先来计算一下氧气分子的零点能矫正。气体分子的振动频率计算我们在之前的教程中已经介绍过：ABACUS计算模拟实例 | II. C2H5OH的振动模式与频率计算 。这里我们只需要准备优化后的氧气分子的结构文件STRU（即复制优化后的OUT文件夹中的STRU_ION_D文件过来，类似VASP中的CONTCAR）以及振动频率计算文件vib.py作为输入文件即可运行。
OUT STRU ase_sort.dat result vib vib.py
ATOMIC_SPECIES
O 15.999 O_ONCV_PBE-1.0.upf upf201
NUMERICAL_ORBITAL
O_gga_7au_100Ry_2s2p1d.orb
LATTICE_CONSTANT
1.8897259886
LATTICE_VECTORS
10.0000000000 0.0000000000 0.0000000000 #latvec1
0.0000000000 10.0000000000 0.0000000000 #latvec2
0.0000000000 0.0000000000 10.0000000000 #latvec3
ATOMIC_POSITIONS
Direct
O #label
1 #magnetism
2 #number of atoms
0.0000000000 0.0000000000 0.0000000000 m 0 0 0 mag 1.0000000000
0.0000000000 0.0000000000 0.1208000000 m 0 0 0 mag 1.0000000000
import os
import subprocess
import shutil
from ase.optimize import QuasiNewton
from ase.vibrations import Vibrations
from ase.thermochemistry import HarmonicThermo
from ase.io import read, write
from ase.io import Trajectory
from ase.visualize import view
from ase.calculators.abacus import Abacus, AbacusProfile
# env setting
os.environ['ABACUS_PP_PATH'] = '/abacus-develop/tests/PP_ORB'
os.environ['ABACUS_ORBITAL_PATH'] = '/abacus-develop/tests/PP_ORB'
stru = read('STRU', format='abacus')
pseudo_dir = os.environ['ABACUS_PP_PATH']
basis_dir = os.environ['ABACUS_ORBITAL_PATH']
out_dir = "OUT"
properties = ["energy", "forces", "stress"]
pp = {
'O':'O_ONCV_PBE-1.0.upf',
}
basis = {
'O': 'O_gga_7au_100Ry_2s2p1d.orb',
}
kpts = [1, 1, 1]
# INPUT setting
parameters = {
'suffix': 'Oxygen',
'calculation': 'scf',
'basis_type': 'lcao',
'ks_solver': 'genelpa',
'nspin' : 2,
'xc': 'pbe',
'symmetry': '0',
'ecutwfc': 80,
'scf_thr': 1e-7,
'kpts': kpts,
'pp': pp,
'basis': basis,
'pseudo_dir': pseudo_dir,
'basis_dir': basis_dir,
'smearing_method': 'gauss',
'smearing_sigma': 0.0037,
'cal_force': 1,
'cal_stress': 1,
'out_stru': 0,
'out_chg': 0,
}
# attch calculator to the model
profile = AbacusProfile(argv=['abacus'])
stru.calc = Abacus(profile=profile, directory=out_dir, **parameters)
# vibration calculation
vib = Vibrations(stru, nfree=2)
vib.run()
vib.summary()
--------------------- # meV cm^-1 --------------------- 0 38.6i 311.4i 1 38.6i 311.4i 2 0.1i 0.8i 3 0.6 5.1 4 0.6 5.1 5 208.6 1682.1 --------------------- Zero-point energy: 0.105 eV
可以看到输出中一共有六种振动模式，其中频率最高的模式对应两个原子向相反方向的振动，可视为氧气分子的拉伸；另外两个实的振动频率对应的模式为氧气分子在xy和yz平面的旋转；最后三个虚的振动频率对应的模式为氧气分子分别沿x、y、z方向的平移。
Note：
理论上，对于线性分子，振动模式的数量等于3N-5，因此这里的氧气分子应该只有一种振动模式（另外两个实频对应旋转运动，三个虚频对应平移运动）。DFT计算出现的其他小频率振动模式是有限差分方法普遍存在的误差。因此我们这里只取最大的振动频率对应的能量0.1043 eV作为氧气分子的零点能。
接下来我们计算O原子在Pt(111)表面的FCC位点时的零点能矫正。零点能计算中，我们经常忽略对零点能修正贡献不大的原子，比如slab原子，因为吸附前后它们几乎不移动，因此在振动计算时固定这部分原子，只计算吸附质的正态振型。
类似上面，这里我们只需要准备优化后的吸附构型文件STRU以及振动频率计算文件vib_analysis.py作为输入文件即可运行。这里的频率计算文件来自于ATST-Tools（ATST-Tools官网见：https://github.com/QuantumMisaka/ATST-Tools， 教程见： ATST-Tools ），是基于ASE调用ABACUS进行过渡态搜索计算的工作流脚本集，涵盖ASE中包含的NEB，AutoNEB和Dimer方法。在ABACUS:3.6.0-userguide镜像中我们将源文件克隆到了/opt下，因此可以直接复制脚本到工作文件夹，然后简单修改即可使用（这里我们给出已经修改好的脚本）。
INPUT SCF-rank0 vib.0.traj O_ONCV_PBE-1.0.upf STRU vib.1.traj O_gga_7au_100Ry_2s2p1d.orb ase_sort.dat vib.2.traj Pt_ONCV_PBE-1.0.upf result vib_analysis.py Pt_gga_8au_100Ry_4s2p2d1f.orb vib
# 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
atoms = read('STRU', format = 'abacus')
vib_indices = [20]
T = 298.15 # K
vib_name = 'vib'
delta = 0.01
nfree = 2
abacus = "abacus"
mpi = 32
omp = 1
lib_dir = "./"
pseudo_dir = f"{lib_dir}/"
basis_dir = f"{lib_dir}/"
pp = {
'O': 'O_ONCV_PBE-1.0.upf',
'Pt': 'Pt_ONCV_PBE-1.0.upf',
}
basis = {
'O': 'O_gga_7au_100Ry_2s2p1d.orb',
'Pt': 'Pt_gga_8au_100Ry_4s2p2d1f.orb',
}
kpts = [5, 5, 1]
parameters = {
'calculation': 'scf',
'basis_type': 'lcao',
'ks_solver': 'genelpa',
'nspin' : 1,
'xc': 'pbe',
'symmetry': '0',
'ecutwfc': 80,
'scf_thr': 1e-7,
'kpts': kpts,
'pp': pp,
'basis': basis,
'pseudo_dir': pseudo_dir,
'basis_dir': basis_dir,
'smearing_method': 'mp',
'smearing_sigma': 0.008,
'cal_force': 1,
'cal_stress': 1,
'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 <==")
==> Starting Vibrational Analysis <==
==> Running Vibrational Analysis <==
==> Done !!! <==
==> All force cache will be in vib directory <==
==> Vibrational Analysis Summary <==
---------------------
# meV cm^-1
---------------------
0 47.5 382.9
1 47.7 384.9
2 54.0 435.6
---------------------
Zero-point energy: 0.075 eV
==> Writing All Mode Trajectory <==
==> Doing Harmonmic Thermodynamic Analysis <==
Entropy components at T = 298.15 K:
=================================================
S T*S
S_harm 0.0001253 eV/K 0.037 eV
-------------------------------------------------
S 0.0001253 eV/K 0.037 eV
=================================================
Internal energy components at T = 298.15 K:
===============================
E_pot 0.000 eV
E_ZPE 0.075 eV
Cv_harm (0->T) 0.025 eV
-------------------------------
U 0.100 eV
===============================
Entropy components at T = 298.15 K:
=================================================
S T*S
S_harm 0.0001253 eV/K 0.037 eV
-------------------------------------------------
S 0.0001253 eV/K 0.037 eV
=================================================
Free energy components at T = 298.15 K:
=======================
U 0.100 eV
-T*S -0.037 eV
-----------------------
F 0.062 eV
=======================
==> Entropy: 1.252518e-04 eV/K <==
==> Free Energy: 0.062480 eV <==
这里氧原子的三种振动频率对应的振型均为平移，方向分别沿三个基矢的方向，其中振动频率最大的振型为沿z轴的平移。
Note：
使用abacus做振动计算时，在STRU中固定原子是没有用的，因为振动计算使用的是ase的ase.vibrations.Vibrations模块，abacus在这里只是充当计算器。因此需要在振动计算脚本中设置进行振动的原子的index，并在创建该类的对象时作为初始化参数传递。
由此，结合吸附能和零点能的定义，我们可以得到以下经零点能修正的吸附能计算公式
恭喜你完成了第十一个ABACUS计算模拟实例🎉
如果你还想尝试更多的计算实例，可以点击下方了解ABACUS计算模拟合集链接： ABACUS计算模拟实例 | 概述
kirk0830