空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

ABACUS计算模拟实例 | XIII. Pt表面的ORR催化路径

ABACUS

计算材料学

ABACUS计算材料学

FermiNoodles

weiqingzhou@whu.edu.cn

更新于 2024-07-05

推荐镜像 :abacus:3.6.3-user-guide

推荐机型 :c2_m4_cpu

数据集

ABACUS-cases(v12)

ABACUS计算模拟实例 | XIII. Pt表面的ORR催化路径

代码

文本

©️ Copyright 2023 @ Authors
作者：Jinghang Wang (AISI电子结构团队实习生) 📨
日期：2024-5-25
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：由于本案例需要比较长的计算时间（使用64核需要10h以上），因此本案例不会要求完成具体计算，只涉及预处理、部分运行命令展示以及结果分析，故使用c2_m4_cpu即可。如果你想进行具体计算，点击上方的开始连接按钮，选择 ABACUS:3.6.3-userguide 镜像及 c64_m128_cpu 节点配置，并在连接后点击右上方的切换内核，选择Bash内核稍等片刻即可运行。

代码

文本

🎯 本教程旨在帮助初学者学习使用ABACUS使用完成完整的Pt表面的ORR催化路径的计算

代码

文本

本节教程目录如下：
1. 计算氢电极(CHE)模型
2. 基于 CHE 模型的氧还原反应
3. ORR反应中的自由能台阶图计算

代码

文本

在使用ABACUS进行计算时，所有操作都需要在指定的文件夹目录下进行。我们已经准备好了数据集ABACUS-cases，，其中文件夹“13”对应本节教程。

代码

文本

首先我们将数据集中的文件复制到根目录下（否则可能碰到文件没有执行权限的情况）：

代码

文本

[1]

cp -r /bohr/ABACUS-cases-nu6t/v12/ABACUS-cases/13 .

代码

文本

[2]

cd 13 && ls

ZPE  relax-OH  relax-OOH  relax-molecule

代码

文本

relax-OH：包含OH在Pt(111)表面吸附的四种位点的吸附能计算。
relax-OOH：包含OH在Pt(111)表面吸附的四种位点的吸附能计算。
（O在Pt(111)表面吸附的四种位点的吸附能计算、零点能矫正以及热力学矫正我们在上上节教程中已经完成）
relax-molecule：包含H₂和H₂O分子的弛豫计算，从而得到每个分子的总能。
（O₂分子的总能在上上节教程中已经完成）
ZPE：包含H₂, O₂, H₂O分子以及O, OH, OOH吸附质的零点能矫正以及热力学矫正。

代码

文本

电解水制氢和氢氧燃料电池等氢能利用技术已成为解决未来能源问题的重要途径，相应地，氢析出反应 (HER)、氢氧化反应 (HOR)、氧析出反应 (OER) 和氧还原反应 (ORR) 等电化学反应也成为研究的热点。而质子交换膜燃料电池 (PEMFC) 由于零污染、高效能和易轻量化等独特优势，被视为一种具有极大潜力的动力来源。发展低贵金属含量、高电催化活性和稳定性的 ORR 催化剂则是 PEMEC 技术成为一门规模化实用技术的关键所在。迄今为止，Pt及其合金仍被认为是最具实用价值的 ORR 催化剂。本节就以 Pt 电极为例，介绍如何通过 DFT 计算研究 Pt表面的ORR 催化路径。

代码

文本

1. 计算氢电极(CHE)模型

代码

文本

电催化反应是一种非常复杂的系统过程，包含众多的影响因素，如溶剂化效应、电极表面结构和固液界面的电荷分布等。这些因素的相互作用，使得完全考虑所有影响因素是十分困难的。Norskov 等人于 2004 年提出了一种基于吸附能的密度泛函理论模型，用于模拟电催化反应过程，并于2010年将其命名为计算氢电极 (computational hydrogen electrode，CHE) 模型。该模型基于密度泛函理论，结合吸附能的概念，能够模拟电催化反应的基本过程，如吸附、解离和反应等。CHE 模型不仅可以真实地反映电极表面催化活性，还可以为电极材料的设计提供指导。因为它提供了一种定量和预测性的方法，可以帮助我们深入了解原子层面上催化剂的行为。它的应用已经带来了在能量储存、燃料电池和可再生能源等各种领域中开发出新型和改进型催化剂的许多突破。以下对其进行简要介绍。

代码

文本

对于一个标准氢电极 (standard hydrogen electrode，SHE)，氧化还原半反应式可写为 $\frac{1}{2} H_{2} (g) ⇌ H^{+} + e^{-} (1)$ 当氢电极处于pH=0、氢气压力为1 bar、温度为 298.15 K 的环境中且外加电压 U=0 V时，氢电极氧化还原半反应的自由能变化为 $Δ G = 0$ eV。因此，质子电子对的自由能可写为 $G (H^{+} + e^{-}) = \frac{1}{2} G (H_{2} (g)) (2)$ 若有外加电压，则式(1)需加上电压修正项-eU,其中e为基本电荷，U 为外加电压。此外，在电催化过程中经常通过调整 pH 值来促进反应的进行，因此也需对 pH 值进行修正。pH的修正项可以表示为 $Δ G (p H) = - k_{B} T l n [H^{+}] = k_{B} T l n (10) \times p H (3)$ 式中：k_B为玻尔兹曼常数。若以eV 为单位代入各项值，则 G(pH)=0.0592 pH。pH值可视为质子活性。综上所，在外加电场和 pH 值均不为0的情况下,质子电子对的自由能可写为 $G (H^{+} + e^{-}) = \frac{1}{2} G (H_{2} (g)) - e U - 0.0592 p H (4)$ 进一步，对于一个还原的吸热反应,有 $A^{*} + (H^{+} + e^{-}) = A H^{*} (5)$ 吸附物种 A(A*)经还原产生对应的氢化物(AH*)，该反应的自由能变化为 $Δ G = G (A H^{*}) - G (A^{*}) - G (H^{+} + e^{-}) = G (A H^{*}) - (G (A^{*}) + \frac{1}{2} G (H_{2} (g)) - e U - 0.0592 p H) (6)$ 当 pH=0 且有外加电压时，为使反应达到平衡状态，所需外加电压为 $U = \frac{G ( A H ^{*} ) - G ( A ^{*} ) - G ( H ^{+} + e ^{-} )}{e} (7)$ 这个电压值即为平衡电位(equilibrium potential)。至于各反应物种的自由能，则可通过DFT计算获得 $G = E_{D FT} + G_{cor} = E_{D FT} + (E_{ZPE} + ν \sum \frac{h ν}{e ^{h ν / k_{B} T} - 1} - TS) (8)$ Norskov通过这样一个相对的单的模型将表面科学程架与复杂的电化学环境成功联系起来，为电催化的反应机理研究和催化剂的理性设计等提供了巧妙、实用的方案。

代码

文本

2. 基于 CHE 模型的氧还原反应

代码

文本

如图 1 所示，氧还原反应是一种复杂的多电子、多步骤反应，涉及多种反应中间体，并伴随着质子迁移过程。氢燃料电池的整体反应式为 $2 H_{2} (g) + O_{2} (g) \to 2 H_{2} O (l) + 反应热 (9)$

图1 PEMFC中的氧还原反应示意图（来源: FuelCellsWorks ）

这个反应式表明，在氢燃料电池中，氢气和氧气在催化剂的作用下将发生氧化还原反应，生成水和热能。反应热则是在 298.15 K、1 atm 的条件下反应释放的能量，等于 4.92 eV。该反应通常在酸性电解液中进行，在 SHE 电极发生的半反应为

2 H_{2} (g) \to 4 H^{+} + 4 e^{-} (10)

而在酸性电解液中，ORR 电极上的氧还原通常存在两种反应路径：
(1) 直接四电子转移路径，即发生四个电子和四个质子的转移：

O_{2} (g) + 4 H^{+} + 4 e^{-} \to 2 H_{2} O (l) (11)

(2) 连续两电子转移路径,即发生电子的连续转移:
$O_{2} (g) + 2 H^{+} + 2 e^{-} \to H_{2} O (l) (12)$ $H_{2} O_{2} + 2 H^{+} + 2 e^{-} \to 2 H_{2} O (l) (13)$

以广受认可的四电子转移路径为例，ORR 反应过程可进一步根据电子转移拆解为四个放热反应步骤： $* + O_{2} (g) + (H^{+} + e^{-}) \to OO H^{*} OO H^{*} + (H^{+} + e^{-}) \to O^{*} + H_{2} O (l) O^{*} + (H^{+} + e^{-}) \to O H^{*} O H^{*} + (H^{+} + e^{-}) \to * + H_{2} O (l) (14)$

其中 * 表示洁净表面。根据 CHE 模型，pH=0且无外加电压时，每一步的自由能变化为： $Δ G (1) = G (OO H^{*}) - \frac{1}{2} G (H_{2}) - G (*) - G (O_{2}) (15.1)$ $Δ G (2) = G (O^{*}) - \frac{1}{2} G (H_{2}) - G (OO H^{*}) + G (H_{2} O) (15.2)$ $Δ G (3) = G (O H^{*}) - \frac{1}{2} G (H_{2}) - G (O^{*}) (15.3)$ $Δ G (1) = G (*) - \frac{1}{2} G (H_{2}) - G (O H^{*}) + G (H_{2} O) (15.4)$

代码

文本

其中,H₂(g)、H₂O(l)和 O₂(g)分子的自由能可通过ABACUS计算结合ase热力学修正得到。

代码

文本

前面我们已经介绍了使用ATST-Tools脚本完成零点能矫正。ATST-Tools同样提供了同时完成零点能和热力学矫正的脚本。在推荐镜像中ATST-Tools被放在了/opt文件夹下，因此可以直接拷贝脚本到计算文件夹（这里我们直接提供修改好了的脚本）。
热力学矫正使用的是ase.thermochemistry模块，包含了四种情况：理想气体限制，只考虑平动和转动自由度，适用于气体分子；谐振限制，只考虑振动自由度，适用于吸附模型；受限平动/受限转动模型，考虑两个自由度为平动，一个自由度为转动，剩余3N-3个自由度为振动，同样适用于吸附模型；晶体固体模型，在这种模型中，一个含有N个原子的晶格被视为一个包含3N个独立谐振子的系统，适用于固体。
在下面的计算中，对于分子体系，我们都采取理想气体限制；对于吸附质，我们都采取谐振限制。

代码

文本

[3]

cd ZPE && ls # ZPE

H2  H2O  O2  OH  OOH

代码

文本

[4]

cd H2

代码

文本

[5]

cat STRU

ATOMIC_SPECIES
H   1.008         H_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
H_gga_8au_100Ry_2s1p.orb

LATTICE_CONSTANT
1.889726

LATTICE_VECTORS
10.000 0.000 0.000
0.000 10.000 0.000
0.000 0.000 10.000

ATOMIC_POSITIONS
Cartesian                          # Cartesian or Direct coordinate.
H                                  # Element type
0.0                                # magnetism
2                                  # number of atoms
0.000 0.000 0.000 0 0 0            # the position of atoms and other parameter specify by key word
0.000 0.000 0.742 0 0 0

代码

文本

[6]

cat idealgas_analysis.py # script for ZPE and thermodynamic correction

# JamesMisaka in 2024-01-16
# Thermodynamic analysis for molecular (in a box)
# part of ATST-Tools scripts

import os
import numpy as np
from ase.vibrations import Vibrations
from ase.thermochemistry import IdealGasThermo
from ase.io import read, write
from ase.calculators.abacus import Abacus, AbacusProfile
from ase.parallel import world, parprint
from ase import units

# indices setting for which atoms to be displaced
atoms = read('STRU', format='abacus')
vib_indices = None
geometry = 'linear'
symmetrynumber = 1

T = 298.15 # K
P = 1E5 # Pa

vib_name = 'vib'
delta = 0.01
nfree = 2

abacus = "abacus"
mpi = 10
omp = 1
lib_dir = "/abacus-develop/tests/PP_ORB"
pseudo_dir = f"{lib_dir}/"
basis_dir = f"{lib_dir}/"
pp = {
        'H':'H_ONCV_PBE-1.0.upf',
      }
basis = {
        'H': 'H_gga_8au_100Ry_2s1p.orb',
        }
kpts = [1, 1, 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 Ideal-Gas Thermodynamic Analysis <==")
    # initias = (atoms.get_moments_of_inertia() * units._amu /
    #                     (10.0**10)**2)
    # print(f"==> Initial Moments of Inertia for {geometry} {atoms.symbols} molecular: {initias} kg*m^2 <==")
    gasthermo = IdealGasThermo(vib.get_energies(),
                           geometry=geometry,
                           atoms=atoms,
                           ignore_imag_modes=True,
                           symmetrynumber=symmetrynumber,
                           spin=0,)
    entropy = gasthermo.get_entropy(T,P)
    free_energy = gasthermo.get_gibbs_energy(T,P)
    print(f"==> Entropy: {entropy:.6e} eV/K <==")
    print(f"==> Gibbs Free Energy: {free_energy:.6f} eV <==")

代码

文本

脚本中的关键参数：
1.geometry：monatomic,linear or nonlinear
2.symmetrynumber: symmetry number of the molecule. See, for example, Table 10.1 and Appendix B of C. Cramer “Essentials of Computational Chemistry”, 2nd Ed.
3.spin: the total electronic spin. (0 for molecules in which all electrons are paired, 0.5 for a free radical with a single unpaired electron, 1.0 for a triplet with two unpaired electrons, such as O₂.)

代码

文本

输出被保存在result文件中。这里我们没有提供势能信息给ase，因此输出的吉布斯自由能即零点能和热力学矫正值之和。

代码

文本

[7]

cat result

==> Starting Vibrational Analysis <==
==> Running Vibrational Analysis <==
==> Done !!! <==
==> All force cache will be in vib directory <==
==> Vibrational Analysis Summary <==
---------------------
  #    meV     cm^-1
---------------------
  0   70.8i    570.9i
  1   70.8i    570.9i
  2    0.0i      0.0i
  3    0.0       0.0
  4    0.0       0.0
  5  557.5    4496.4
---------------------
Zero-point energy: 0.279 eV
==> Writing All Mode Trajectory <==
==> Doing Ideal-Gas Thermodynamic Analysis <==
Entropy components at T = 298.15 K and P = 100000.0 Pa:
=================================================
                           S               T*S
S_trans (1 bar)    0.0012188 eV/K        0.363 eV
S_rot              0.0001919 eV/K        0.057 eV
S_elec             0.0000000 eV/K        0.000 eV
S_vib              0.0000000 eV/K        0.000 eV
S (1 bar -> P)    -0.0000000 eV/K       -0.000 eV
-------------------------------------------------
S                  0.0014107 eV/K        0.421 eV
=================================================
Enthalpy components at T = 298.15 K:
===============================
E_pot                  0.000 eV
E_ZPE                  0.279 eV
Cv_trans (0->T)        0.039 eV
Cv_rot (0->T)          0.026 eV
Cv_vib (0->T)          0.000 eV
(C_v -> C_p)           0.026 eV
-------------------------------
H                      0.369 eV
===============================

Entropy components at T = 298.15 K and P = 100000.0 Pa:
=================================================
                           S               T*S
S_trans (1 bar)    0.0012188 eV/K        0.363 eV
S_rot              0.0001919 eV/K        0.057 eV
S_elec             0.0000000 eV/K        0.000 eV
S_vib              0.0000000 eV/K        0.000 eV
S (1 bar -> P)    -0.0000000 eV/K       -0.000 eV
-------------------------------------------------
S                  0.0014107 eV/K        0.421 eV
=================================================

Free energy components at T = 298.15 K and P = 100000.0 Pa:
=======================
    H          0.369 eV
 -T*S         -0.421 eV
-----------------------
    G         -0.052 eV
=======================
==> Entropy: 1.410740e-03 eV/K <==
==> Gibbs Free Energy: -0.051945 eV <==

代码

文本

表1列出了标准氢电极条件下H₂，O₂ 和 H₂O 分子的自由能的计算结果。

代码

文本

[8]

%%python

import pandas as pd

from IPython.display import HTML

# data dictionary

data = {

'反应物': ['H<sub>2</sub>(g)', 'H<sub>2</sub>O(l)', 'O<sub>2</sub>(g)'],

'压强/atm': ['1.000', '0.035', '1.000'],

'温度/K': ['298.15']*3,

'E<sub>DFT</sub>/eV': ['-31.689', '-466.511', '—'],

'G<sub>cor</sub>/eV':['-0.052', '0.007', '—'],

'G/eV':['-31.741', '-466.504', '-864.606'],

}

# create DataFrame

df = pd.DataFrame(data)

# CSS properties

style = """

<style>

table, th, td {

border: 1px solid black;

border-collapse: collapse;

font-family: Courier, monospace; /* 使用等宽字体 */

}

th, td {

padding: 10px;

text-align: center; /* 文本居中 */

}

th {

background-color: #f2f2f2;

}

table {

width: 100%; /* 表格宽度为100%，这将导致列等间隔 */

margin-left: auto;

margin-right: auto;

}

</style>

"""

table_title = "<h2>表1 标准氢电极条件下H<sub>2</sub>，O<sub>2</sub> 和 H<sub>2</sub>O 分子的自由能</h2>"

# use HTML

html = HTML(style + table_title + df.to_html(classes='table table-striped', escape=False, index=False))

# show

html

表1 标准氢电极条件下H₂，O₂ 和 H₂O 分子的自由能

反应物	压强/atm	温度/K	E_DFT/eV	G_cor/eV	G/eV
H₂(g)	1.000	298.15	-31.689	-0.052	-31.741
H₂O(l)	0.035	298.15	-466.511	0.007	-466.504
O₂(g)	1.000	298.15	—	—	-864.606

代码

文本

注： $G (O_{2}) = 2 G (H_{2} O) - 2 G (H_{2}) + 4.920$ 。

代码

文本

进一步，针对不同的 ORR 催化剂，计算催化剂纯净表面 (* ) 以及分别吸附三种中间物种 (OH、O和 OOH^*) 时的自由能，代入式（15）即可得到每一步的自由能变化值，最后通过式（16）计算起始电位 $η_{o n se t}$ $η_{o n se t} = - min (Δ (1), Δ (2), Δ (3), Δ (4)) (16)$

代码

文本

如图 2 所示，理想情况下，将总反应的自由能变化平均分配到四个电子反应步骤中，即每一步放热 1.23 eV。若平衡电位与起始电位均为 1.23 V,则过电位为 0.00 V。

代码

文本

代码

文本

下面以 Pt(111)体系为例,介绍其 ORR 四电子转移路径的自由能台阶图绘制方法。

代码

文本

3. ORR反应中的自由能台阶图计算

代码

文本

基于前面的分析，只需通过结构优化计算Pt(111)的洁净slab模型及其分别吸附OH、O和OOH物种的体系总能，并通过计算吸附物种的频率做进一步的热力学校正，将所得结果代入式（15），得到每一步的自由能。
不同物种在 Pt(111)表面不同位点的吸附强度不同，因此，有必要测试三个中间物种在所有可能位点的吸附能，从而判断出各自最稳定的吸附构型用于自由能的计算。（其中OH在Bri位点和O的四种位点的总能我们分别在上一节和上上节已经计算过）。对于本节中的弛豫计算，输入文件除了STRU外完全一致。

代码

文本

[8]

cd ../../relax-OH/ && ls

FCC                       O_ONCV_PBE-1.0.upf             Top
HCP                       O_gga_7au_100Ry_2s2p1d.orb     relax.sh
H_ONCV_PBE-1.0.upf        Pt_ONCV_PBE-1.0.upf
H_gga_8au_100Ry_2s1p.orb  Pt_gga_8au_100Ry_4s2p2d1f.orb

代码

文本

[9]

cat relax.sh # script for doing all relaxation calculations

## relaxation of 3 adsorption sites
cd FCC
echo -e "FCC relax begin!\n"
mpirun -n 32 abacus > result
echo -e "FCC relax complete!\n"
cd ../HCP
echo -e "HCP relax begin!\n"
mpirun -n 32 abacus > result
echo -e "HCP relax complete!\n"
cd ../Top
echo -e "Top relax begin!\n"
mpirun -n 32 abacus > result
echo -e "Top relax complete!\n"

代码

文本

[12]

#. relax.sh #如果你想尝试运行，请切换c32_m64_cpu节点，去掉井号即可运行

代码

文本

[10]

grep -i -n "FINAL_ETOT" FCC/OUT.ABACUS/running_relax.log

grep -i -n "FINAL_ETOT" Top/OUT.ABACUS/running_relax.log

grep -i -n "FINAL_ETOT" HCP/OUT.ABACUS/running_relax.log

34761: !FINAL_ETOT_IS -66456.1120421677187551 eV
33610: !FINAL_ETOT_IS -66456.4965231645182939 eV
34420: !FINAL_ETOT_IS -66455.8749915616062935 eV

代码

文本

[11]

cd FCC && cat INPUT # input files

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

代码

文本

[12]

cat KPT

K_POINTS
0
Gamma
5 5 1 0 0 0

代码

文本

[13]

cat STRU

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.6656666667 0.3343333333 0.5750000000 1 1 1 mag 0.0 

H
0.0000000000
1
0.6656666667 0.3343333333 0.6250000000 1 1 1 mag 0.0

代码

文本

[14]

cd ../../relax-OOH && ls

Bri  H_ONCV_PBE-1.0.upf        O_gga_7au_100Ry_2s2p1d.orb     Top
FCC  H_gga_8au_100Ry_2s1p.orb  Pt_ONCV_PBE-1.0.upf            relax.sh
HCP  O_ONCV_PBE-1.0.upf        Pt_gga_8au_100Ry_4s2p2d1f.orb

代码

文本

[18]

#. relax.sh

代码

文本

[15]

grep -i -n "FINAL_ETOT" FCC/OUT.ABACUS/running_relax.log

grep -i -n "FINAL_ETOT" Top/OUT.ABACUS/running_relax.log

grep -i -n "FINAL_ETOT" HCP/OUT.ABACUS/running_relax.log

grep -i -n "FINAL_ETOT" Bri/OUT.ABACUS/running_relax.log

19427: !FINAL_ETOT_IS -66887.9345830670354189 eV
19968: !FINAL_ETOT_IS -66888.3929906853591092 eV
7638: !FINAL_ETOT_IS -66887.7688679306447739 eV
26516: !FINAL_ETOT_IS -66888.0983206890086876 eV

代码

文本

slab的总能为：-66006.58 eV。我们将物种X的吸附能的计算算结果总结于表2。不难看出,O倾向于吸附在 Pt(111)的 FCC 位点,而 OH和OOH吸附在 Top 位点时更稳定。

代码

文本

[20]

%%python

import pandas as pd

from IPython.display import HTML

# data dictionary

data = {

'Top': ['0.25', '-1.60', '-1.03'],

'Bri': ['—', '-1.52', '-0.73'],

'HCP': ['-0.65', '-0.98', '-0.40'],

'FCC': ['-1.08', '-1.22', '-0.57']

}

# createDataFrame

df = pd.DataFrame(data, index=['O*', 'OH*', 'OOH*'])

# CSS properties

style = """

<style>

table, th, td {

border: 1px solid black;

border-collapse: collapse;

font-family: Courier, monospace; # 使用等宽字体

}

th, td {

padding: 10px;

text-align: right; # 文本右对齐可以更好地对齐数字

width: 60px; # 指定每列的固定宽度

}

th {

background-color: #f2f2f2;

}

table {

width: auto; # 表格宽度自适应内容

margin-left: auto;

margin-right: auto;

}

</style>

"""

table_title = "<h2>表2 O*、OH*和OOH*在Pt(111)表面上不同位点的吸附能</h2>"

# use HTML

html = HTML(style + table_title + df.to_html(classes='table table-striped', index=True))

# show

html

表2 O、OH和OOH*在PT(111)表面上不同位点的吸附能

	Top	Bri	HCP	FCC
O*	0.25	—	-0.65	-1.08
OH*	-1.60	-1.52	-0.98	-1.22
OOH*	-1.03	-0.73	-0.40	-0.57

代码

文本

进而针对每种最稳定的吸附构型,计算相应中间物种的零点能和热力学矫正，得到相应的自由能。对于吸附模型的热力学修正，我们还是使用ATST-Tools中的脚本vib_analysis.py来完成。同样，输出中最后一行的Free Energy即零点能和热力学矫正值之和。

代码

文本

[16]

cd ../ZPE

代码

文本

[17]

cd OH && ls

SCF-rank0  ase_sort.dat  vib         vib.1.traj  vib.3.traj  vib.5.traj
STRU       result        vib.0.traj  vib.2.traj  vib.4.traj  vib_analysis.py

代码

文本

[18]

cat vib_analysis.py # script for ZPE and thermodynamic correction

# 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,21]

T = 298.15 # K

vib_name = 'vib'
delta = 0.01
nfree = 2

abacus = "abacus"
mpi = 32
omp = 1
lib_dir = "../../relax-OH"
pseudo_dir = f"{lib_dir}/"
basis_dir = f"{lib_dir}/"
pp = {
      'H': 'H_ONCV_PBE-1.0.upf',
      'O': 'O_ONCV_PBE-1.0.upf',
      'Pt': 'Pt_ONCV_PBE-1.0.upf',
      }
basis = {
         'H': 'H_gga_8au_100Ry_2s1p.orb',
         '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 <==")

代码

文本

[24]

#python3 vib_analysis.py > result

代码

文本

[19]

cat result

==> Starting Vibrational Analysis <==
==> Running Vibrational Analysis <==
==> Done !!! <==
==> All force cache will be in vib directory <==
==> Vibrational Analysis Summary <==
---------------------
  #    meV     cm^-1
---------------------
  0    5.5      44.1
  1   15.7     126.9
  2   18.5     149.5
  3   64.9     523.8
  4  115.2     929.5
  5  451.8    3644.3
---------------------
Zero-point energy: 0.336 eV
==> Writing All Mode Trajectory <==
==> Doing Harmonmic Thermodynamic Analysis <==
Entropy components at T = 298.15 K:
=================================================
                           S               T*S
S_harm             0.0004971 eV/K        0.148 eV
-------------------------------------------------
S                  0.0004971 eV/K        0.148 eV
=================================================
Internal energy components at T = 298.15 K:
===============================
E_pot                  0.000 eV
E_ZPE                  0.336 eV
Cv_harm (0->T)         0.066 eV
-------------------------------
U                      0.402 eV
===============================

Entropy components at T = 298.15 K:
=================================================
                           S               T*S
S_harm             0.0004971 eV/K        0.148 eV
-------------------------------------------------
S                  0.0004971 eV/K        0.148 eV
=================================================

Free energy components at T = 298.15 K:
=======================
    U          0.402 eV
 -T*S         -0.148 eV
-----------------------
    F          0.254 eV
=======================
==> Entropy: 4.971353e-04 eV/K <==
==> Free Energy: 0.253815 eV <==

代码

文本

整理后的结果如表 3 所示：

代码

文本

[26]

%%python

import pandas as pd

from IPython.display import HTML

# data dictionary

data = {

'物种': ['slab', 'slab+OOH*', 'slab+O*', 'slab+OH*', 'slab'],

'E<sub>DFT</sub>': ['-66006.58', '-66888.393', '-66440.133','-66456.497', '-66006.58'],

'G<sub>cor</sub>':['0', '0.353', '0.062', '0.254', '0'],

'分子与电子质子对':['O<sub>2</sub>+4(H<sup>+</sup>+e<sup>-</sup>)', '3(H<sup>+</sup>+e<sup>-</sup>)', 'H<sub>2</sub>O+2(H<sup>+</sup>+e<sup>-</sup>)', 'H<sub>2</sub>O+(H<sup>+</sup>+e<sup>-</sup>)', '2H<sub>2</sub>O'],

'G<sub>tot</sub>':['-66,934.668', '-66,935.652', '-66938.316', '-66938.618', '-66939.588'],

'G<sub>step</sub><sup>U=0.00eV</sup>':['4.92', '3.94', '1.27', '0.97', '0.00'],

'G<sub>step</sub><sup>U=1.23eV</sup>':['-31.741', '-466.504', '-2.1585', '-0.2595', '0.00']

}

# create DataFrame

df = pd.DataFrame(data)

# CSS properties

style = """

<style>

table, th, td {

border: 1px solid black;

border-collapse: collapse;

font-family: Courier, monospace; /* 使用等宽字体 */

}

th, td {

padding: 10px;

text-align: center; /* 文本居中 */

}

th {

background-color: #f2f2f2;

}

table {

width: 80%;

margin-left: auto;

margin-right: auto;

}

</style>

"""

table_title = "<h2>表3 标准氢电极条件下Pt(111)表面上ORR中每一步的自由能"

# use HTML

html = HTML(style + table_title + df.to_html(classes='table table-striped', escape=False, index=False))

# show

html

表3 标准氢电极条件下Pt(111)表面上ORR中每一步的自由能

物种 E_DFT G_cor 分子与电子质子对 G_tot G_step^U=0.00eV G_step^U=1.23eV

slab -66006.58 0 O₂+4(H⁺+e^-) -66,934.668 4.92 -31.741

slab+OOH* -66888.393 0.353 3(H⁺+e^-) -66,935.652 3.94 -466.504

slab+O* -66440.133 0.062 H₂O+2(H⁺+e^-) -66938.316 1.27 -2.1585

slab+OH* -66456.497 0.254 H₂O+(H⁺+e^-) -66938.618 0.97 -0.2595

slab -66006.58 0 2H₂O -66939.588 0.00 0.00

物种	E_DFT	G_cor	分子与电子质子对	G_tot	G_step^U=0.00eV	G_step^U=1.23eV
slab	-66006.58	0	O₂+4(H⁺+e^-)	-66,934.668	4.92	-31.741
slab+OOH*	-66888.393	0.353	3(H⁺+e^-)	-66,935.652	3.94	-466.504
slab+O*	-66440.133	0.062	H₂O+2(H⁺+e^-)	-66938.316	1.27	-2.1585
slab+OH*	-66456.497	0.254	H₂O+(H⁺+e^-)	-66938.618	0.97	-0.2595
slab	-66006.58	0	2H₂O	-66939.588	0.00	0.00

代码

文本

最后，结合式（15）和表 1 计算出每步的自由能变化，并绘制出在外加电压为 0.00 V和1.23 V时 Pt(111)表面的 ORR 自由能台阶图，如图 3 所示。可以看到，Pt(111) 表面对 OOH* 物种吸附作用较弱，对 O*、OH* 物种吸附作用较强，从而导致每一步的自由能变化分配不均匀。其中，向 OH* 物种转变的这一步在 U=0.00 V 时放热最少，在U=1.23 V 时自由能台阶上升幅度最大，这是最难的一步，也是整个 ORR 过程的决速步骤。由此可计算 $η_{o n se t} = Δ G (3) = - (0.97 - 1.27) = 0.3 V$ 。

代码

文本

代码

文本

恭喜你完成了最后一个ABACUS计算模拟实例🎉

代码

文本

ABACUS

计算材料学

ABACUS计算材料学

点个赞吧

本文被以下合集收录

ABACUS@计算材料学 | 计算模拟实例

weiqingzhou@whu.edu.cn

更新于 2024-08-09

13 篇8 人关注

ABACUS计算模拟实例

FermiNoodles

更新于 2024-07-15

13 篇1 人关注

ABACUS计算模拟实例 | XIII. Pt表面的ORR催化路径

1. 计算氢电极(CHE)模型

2. 基于 CHE 模型的氧还原反应

表1 标准氢电极条件下H2，O2 和 H2O 分子的自由能

3. ORR反应中的自由能台阶图计算

表2 O*、OH*和OOH*在PT(111)表面上不同位点的吸附能

表1 标准氢电极条件下H₂，O₂ 和 H₂O 分子的自由能

表2 O、OH和OOH*在PT(111)表面上不同位点的吸附能