空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

ABACUS计算模拟实例 | IX. 表面能的计算

ABACUS

DFT

ABACUSDFT

FermiNoodles

weiqingzhou@whu.edu.cn

发布于 2024-05-21

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

推荐机型 :c32_m64_cpu

数据集

ABACUS-cases(v14)

ABACUS计算模拟实例 | IX. 表面能的计算

代码

文本

©️ Copyright 2023 @ Authors
作者：Jinghang Wang (AISI电子结构团队实习生) 📨
日期：2024-5-25
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的开始连接按钮，选择 ABACUS:3.6.3-user-guide 镜像及 c32_m64_cpu 节点配置，并在连接后点击右上方的切换内核，选择Bash内核稍等片刻即可运行。

代码

文本

🎯 本教程旨在展示使用ABACUS完成表面能的计算，主要介绍表面能计算的流程以及对应的ABACUS输入输出文件。

代码

文本

本节教程目录如下：
1. 表面能的定义
2. 表面能的计算
3.伍尔夫结构定理

代码

文本

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

代码

文本

[1]

cp -r /bohr/ABACUS-cases-nu6t/v14/ABACUS-cases/9 .

代码

文本

[2]

cd 9 && ls

Pt_ONCV_PBE-1.0.upf  Pt_gga_8au_100Ry_4s2p2d1f.orb  relax  scf

代码

文本

relax：包含对称和非对称表面的结构优化计算，以获取优化后的表面体系总能。
scf：包含非对称表面的自洽计算，以获取优化前非对称表面的总能。

代码

文本

1. 表面能的定义

代码

文本

图1 FCC晶体中的晶面示意图（来源: Different planes in FCC ）

代码

文本

表面能是一种度量表面形成时分子间化学键断裂强度的参数。在固体理论中，表面原子处于比物质内部原子能量更高的状态。根据能量最低原理，原子会自发地倾向于聚集在物质内部而非表面。表面能的另一种等价定义是材料表面相对于材料内部所具有的额外能量。通过研究表面能，我们可以预测纳米晶体的形态和可能暴露的表面，以进一步预测相关的晶体性质,探究晶体生长机理，并指导新型晶体材料的设计。
此外，表面能的确定在催化领域具有重要意义，可以用于辅助确认催化剂中的反应活性面，从而加深对催化机理的理解。总之，表面能在纳米材料、晶体生长和催化研究中具有重要作用。

代码

文本

图2 表面的对称模型与非对称模型示意图

(a) 对称表面模型；(b) 非对称表面模型

代码

文本

表面模型分为两种——对称模型和非对称模型，如图1所示。其中对称模型包含更多的原子和更大的真空空间，在表面能计算中需要消耗更多的计算资源。非对称模型包含的原子数更少，但具有较低的对称性，在一些情况下会产生偶极矩，导致收敛困难。在ABACUS中可以通过设置dip_cor_flag和efield_flag参数来实现偶极矫正。

对称模型表面能计算公式： $γ_{s} = \frac{E _{s l ab} - n E _{b u l k}}{2 A} (1)$ 式中： $E_{s l ab}$ 为表面模型的能量； $E_{b u l k}$ 为单个原子的能量； $n$ 为表面模型中的原子数目； $A$ 为表面的面积。
非对称模型表面能计算公式： $γ_{a s} = γ_{c l e a v a g e} + γ_{re l a x a t i o n} = \frac{E _{s l ab}^{u n re l a x} - n E _{b u l k}}{2 A} + \frac{E _{s l ab}^{re l a x} - E _{s l ab}^{u n re l a x}}{A} (2)$ 式中: $E_{s l ab}^{u n re l a x}$ 为未弛豫表面模型的能量; $E_{s l ab}^{re l a x}$ 为结构优化后体系的能量。不同于对称模型,非对称模型由于上、下两端原子层不对称(见图2(b)),对称模型中的表面能计算公式不再适用,因为不能直接简单地将能量除以2。

代码

文本

2. 表面能的计算

代码

文本

[3]

cd relax && ls

INPUT       INPUT_s   KPT_slab    STRU_as        STRU_bulk_relax  time.json
INPUT_as    KPT       OUT.ABACUS  STRU_as_relax  STRU_s
INPUT_bulk  KPT_bulk  STRU        STRU_bulk      STRU_s_relax

代码

文本

[4]

cat INPUT_as

INPUT_PARAMETERS
pseudo_dir              ../
orbital_dir             ../

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.05  # ev
# stress_thr            5

basis_type              lcao  # lcao or pw

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

cal_force          1
cal_stress         1
out_stru           1  # print STRU in OUT

#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

代码

文本

[5]

cat INPUT_s

INPUT_PARAMETERS
pseudo_dir              ../
orbital_dir             ../

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.05  # ev
# stress_thr            5

basis_type              lcao  # lcao or pw

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

cal_force          1
cal_stress         1
out_stru           1  # print STRU in OUT

#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

代码

文本

[6]

cat INPUT_bulk

INPUT_PARAMETERS
pseudo_dir              ../
orbital_dir             ../

calculation             relax # scf relax cell-relax md
ecutwfc                 80
scf_thr                 1e-7
scf_nmax                300
relax_nmax              200
relax_method            cg  # cg, bfgs, cg_bfgs, sd, "fire"
force_thr_ev            0.05  # ev
# stress_thr            5

basis_type              lcao  # lcao or pw

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

cal_force          1
cal_stress         1
out_stru           1  # print STRU in OUT

#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

代码

文本

关键参数：

relax_method：bfgs算法更适合表面体系计算；对于bulk，我们选择默认的cg算法。
efield_flag和dip_cor_flag：是否打开偶极矫正。注意，efield_flag需要设为1才能开启偶极矫正，否则会报错 Input warning : dipole correction is not active if efield_flag=false !。然后设置dip_cor_flag=1，同时电场强度设为0（默认值），就是正常的偶极修正了。

代码

文本

[7]

cat KPT_slab

K_POINTS
0
Gamma
5 5 1 0 0 0

代码

文本

[8]

cat KPT_bulk

K_POINTS
0
Gamma
5 5 5 0 0 0

代码

文本

[9]

cat STRU_as

ATOMIC_SPECIES
Pt  195.084       Pt_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
Pt_gga_8au_100Ry_4s2p2d1f.orb

LATTICE_CONSTANT
1.889726

LATTICE_VECTORS
2.7746000000      0.0000000000      0.0000000000      
0.0000000000      2.7746000000      0.0000000000      
0.0000000000      0.0000000000      19.847800000      

ATOMIC_POSITIONS
Direct

Pt
0.0000000000
5
0.0000000000 0.0000000000 0.0000000000 0 0 0 mag 0.0 
0.0000000000 0.0000000000 0.1977000000 0 0 0 mag 0.0 
0.0000000000 0.0000000000 0.3954000000 1 1 1 mag 0.0 
0.5000000000 0.5000000000 0.0988500000 0 0 0 mag 0.0 
0.5000000000 0.5000000000 0.2965500000 1 1 1 mag 0.0

代码

文本

[10]

cat STRU_s

ATOMIC_SPECIES
Pt  195.084       Pt_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
Pt_gga_8au_100Ry_4s2p2d1f.orb

LATTICE_CONSTANT
1.889726

LATTICE_VECTORS
2.7746000000      0.0000000000      0.0000000000      
0.0000000000      2.7746000000      0.0000000000      
0.0000000000      0.0000000000      35.695600000      

ATOMIC_POSITIONS
Direct

Pt
0.0000000000
9
0.0000000000 0.0000000000 0.0000000000 1 1 1 mag 0.0 
0.0000000000 0.0000000000 0.1099300000 0 0 0 mag 0.0 
0.0000000000 0.0000000000 0.2198500000 0 0 0 mag 0.0 
0.0000000000 0.0000000000 0.3297800000 0 0 0 mag 0.0 
0.0000000000 0.0000000000 0.4397100000 1 1 1 mag 0.0 
0.5000000000 0.5000000000 0.0549600000 1 1 1 mag 0.0 
0.5000000000 0.5000000000 0.1648900000 0 0 0 mag 0.0 
0.5000000000 0.5000000000 0.2748200000 0 0 0 mag 0.0 
0.5000000000 0.5000000000 0.3847400000 1 1 1 mag 0.0

代码

文本

对称模型包含九层原子，我们设置中间五层不参与结构优化，以模拟体相原子；非对称模型包含5层原子，我们设置底部三层不参与结构优化。

代码

文本

[11]

cat STRU_bulk

ATOMIC_SPECIES
Pt  195.084       Pt_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
Pt_gga_8au_100Ry_4s2p2d1f.orb

LATTICE_CONSTANT
1.889726

LATTICE_VECTORS
3.9239000000      0.0000000000      0.0000000000      
0.0000000000      3.9239000000      0.0000000000      
0.0000000000      0.0000000000      3.9239000000      

ATOMIC_POSITIONS
Direct

Pt
0.0000000000
4
0.0000000000 0.0000000000 0.0000000000 1 1 1 mag 0.0 
0.0000000000 0.5000000000 0.5000000000 1 1 1 mag 0.0 
0.5000000000 0.0000000000 0.5000000000 1 1 1 mag 0.0 
0.5000000000 0.5000000000 0.0000000000 1 1 1 mag 0.0

代码

文本

[12]

#bulk

cp INPUT_bulk INPUT && cp KPT_bulk KPT

cp STRU_bulk STRU && mpirun -n 32 abacus | tee run_relax_bulk.log # bulk relaxation

cp OUT.ABACUS/running_relax.log OUT.ABACUS/running_relax_bulk.log # save running log

cp OUT.ABACUS/STRU_ION_D ./STRU_bulk # save STRU

# symmetry slab

cp KPT_slab KPT

cp INPUT_s INPUT

cp STRU_s STRU && mpirun -n 32 abacus | tee run_relax_s.log # symmetry slab ralaxation

cp OUT.ABACUS/running_relax.log OUT.ABACUS/running_relax_s.log # save running log

cp OUT.ABACUS/STRU_ION_D ./STRU_s # save STRU

# asymmetry slab

cp INPUT_as INPUT

cp STRU_as STRU && mpirun -n 32 abacus | tee run_relax_as.log # asymmetry slab relaxation

cp OUT.ABACUS/running_relax.log OUT.ABACUS/running_relax_as.log # save running log

cp OUT.ABACUS/STRU_ION_D ./STRU_as # save STRU

已隐藏输出

代码

文本

双击即可修改

代码

文本

由于这里都进行的是结构优化计算，并且输入文件基本相同，因此在同一个目录中进行计算。每一次计算完成后将OUT.ABACUS中的running_relax.log重命名保存，方便查看能量。

代码

文本

[12]

grep -i -n "FINAL_ETOT" OUT.ABACUS/running_relax_s.log

grep -i -n "FINAL_ETOT" OUT.ABACUS/running_relax_as.log

grep -i -n "FINAL_ETOT" OUT.ABACUS/running_relax_bulk.log

19368: !FINAL_ETOT_IS -29702.5959090686010313 eV
4847: !FINAL_ETOT_IS -16500.7473224414170545 eV
1065: !FINAL_ETOT_IS -13202.2086024597592768 eV

代码

文本

[13]

cd ../scf && ls

INPUT  KPT  OUT.ABACUS  STRU  time.json

代码

文本

[16]

cat INPUT

INPUT_PARAMETERS
pseudo_dir              ../
orbital_dir             ../

calculation             scf # 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.05  # ev
# stress_thr            5

basis_type              lcao  # lcao or pw

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

cal_force          1
cal_stress         1
out_stru           0  # print STRU in OUT

代码

文本

[16]

mpirun -n 32 abacus | tee run_scf_as.log # asymmetry slab scf

已隐藏输出

代码

文本

[17]

grep -i -n "FINAL_ETOT" OUT.ABACUS/running_scf.log

1479: !FINAL_ETOT_IS -16500.7469489816721762 eV

代码

文本

根据计算结果，我们可以得到：
对称模型：

$E_{b u l k} = - 3, 300.5521506149484594 e V$
$A = 7.6984 \overset{˚}{A}^{2}$
$n = 9$
$E_{s l ab} = - 29702.5927456220779277 e V$

非对称模型：

$n = 5$
$E_{s l ab}^{u n re l a x} = - 16500.7469491951305827 e V$
$E_{s l ab}^{re l a x} = - 16500.7472156202093174 e V$

代入公式（1）和（2）中得：

$γ_{s} = 0.134 e V / \overset{˚}{A}^{2} = 2.147 J / m^{2}$
$γ_{a s} = 0.131 e V / \overset{˚}{A}^{2} = 2.099 J / m^{2}$

其中 $1 e V / \overset{˚}{A}^{2} = 16.02 J / m^{2}$ 。可以看到，非对称模型的计算结果与九层原子对称模型的结果相仿。需要注意的是，在计算表面能时，我们需要考虑表面模型中真空层厚度的影响。当真空层厚度较小时，上表面原子和下表面原子之间可能会存在一定的相互作用，此时有必要进行测试，观察表面模型能量随真空厚度的变化，确保真空层的厚度收敛。

代码

文本

3. 伍尔夫结构定理

代码

文本

伍尔夫结构定理（Wulff construction rule）是一项确定纳米晶体平衡形状的理论，而纳米晶体是指在1~100 nm尺寸内具有典型形状的晶体。理想的纳米晶体通常呈现出规则的几何多面体形态，为了使表面能尽可能地小，它们通常倾向于暴露具有较低表面能的晶面，或者调整整体形状至接近球形以降低表面积。这两种作用共同决定了纳米晶体的形状。基于晶体的表面能，我们可以预测纳米晶体的形态和形状，进而探索其生长机制，为新晶体材料的设计提供指导。
纳米颗粒在生长过程中总是倾向于使不同晶面的表面能与晶面到原点的距离的比值趋向于一个常量，即伍尔夫结构定律，用公式可表示为 $\frac{σ _{1}}{d _{1}} = \frac{σ _{2}}{d _{2}} = ... = \frac{σ _{i}}{d _{i}} = C (3)$ 式中： $σ_{i}$ 为 $i$ 面的表面能； $d$ 为晶面到原点的距离； $C$ 为一常量。

代码

文本

使用开源软件VESTA绘制Pt纳米颗粒的伍尔夫结构的的教程可见单斌老师的视频：【计算材料学-从算法原理到代码实现】视频教程 | 4.12_基于表面能的纳米颗粒Wulff构建。我们还计算了 (110) 和 (111) 面的表面能分别为 $γ_{110} = 2.172 J / m^{2}$ 和 $γ_{111} = 1.801 J / m^{2}$ 。代入公式（3），并在VESTA中构建模型，得到结果如图3所示。可以看到，能量最高的(110)面在纳米颗粒中并不会出现。

图3 Pt 纳米颗粒的伍尔夫结构示意图。其中紫色代表(111)面，绿色代表(100)面。

代码

文本

恭喜你完成了第九个ABACUS计算模拟实例🎉

代码

文本

如果你还想尝试更多的计算实例，可以点击下方了解ABACUS计算模拟合集链接： ABACUS计算模拟实例 | 概述

代码

文本

ABACUS

DFT

ABACUSDFT

点个赞吧

本文被以下合集收录

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

weiqingzhou@whu.edu.cn

更新于 2024-08-09

13 篇8 人关注

ABACUS计算模拟实例

FermiNoodles

更新于 2024-07-15

13 篇1 人关注