表面能的计算

关于表面能

根据VASP官方手册Hands on Session III，表面能的定义如下，大师兄的文章也提到了这个定义( Nanoscale, 2017, 9, 13089–13094)

σ=σunrelax+Erelax=12(Esurf−Natoms⋅Ebulk)+Erelax\sigma=\sigma^{unrelax}+E^{relax}=\dfrac{1}{2}(E_{surf}-N_{atoms}\cdot E_{bulk})+E^{relax}

γs=12A(Esurfunrelax−N⋅Ebulk)+1A(Esurfrelax−Esurfunrelax)\gamma_s=\dfrac{1}{2A}(E_{surf}^{unrelax}-N\cdot E_{bulk})+\dfrac{1}{A}(E_{surf}^{relax}-E_{surf}^{unrelax})

其中，弛豫能量被定义为：

Erelax=Esurfrelax−EsurfunrelaxE^{relax}=E_{surf}^{relax}-E_{surf}^{unrelax}

  • 手册的公式和大师兄文章中公式的区别在于表面积A，表面能应当除以表面积，所以后者更合理一些
  • EbulkE_{bulk}为体相中单个原子的能量，所以N⋅EbulkN\cdot E_{bulk}为这个体相的能量，即为OUTCAR直接读出的能量

关于表面积A的计算

这部分有些难理解，涉及到一点高等数学和线性代数🤕

表面积可通过晶格基矢量的叉积确定，在POSCAR或CONTCAR文件中，晶格基矢量通常以三行的形式定义

对于面积的计算，可设一组不共线的基矢量为a1=(a11,a12,0),a2=(a21,a22,0)\mathbf{a_1}=(a_{11},a_{12},0),\mathbf{a_2}=(a_{21},a_{22},0)

下面这张图应该挺形象

https://zhuanlan.zhihu.com/p/359975221

a⃗×b⃗=(0,0,a11×a22−a12×a21)\vec{a}\times \vec{b}=(0,0,a_{11}\times a_{22}-a_{12}\times a_{21})

面积为外积模乘以缩放银子

实际面积=s2⋅∣a11a22−a12a21∣实际面积=s^2\cdot |a_{11}a_{22}-a_{12}a_{21}|

PS：对于三维的体积

体积=s3⋅∣a1⋅(a2×a3)∣体积=s^3\cdot |\mathbf{a_1}\cdot (\mathbf{a_2\times a_3})|

p(1x1)Cu(111)的表面能

前面已经算得差不多了

slab模型的能量信息

1
2
3
4
5
6
[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' OUTCAR
  energy  without entropy=      -13.96941183  energy(sigma->0) =      -13.97673168
  energy  without entropy=      -13.96945528  energy(sigma->0) =      -13.97676958
  energy  without entropy=      -13.97006413  energy(sigma->0) =      -13.97726241
  energy  without entropy=      -13.97020643  energy(sigma->0) =      -13.97735554
  energy  without entropy=      -13.97017969  energy(sigma->0) =      -13.97735013

slab模型优化后的结构信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
CONTCAR\(1\1\1)
   1.00000000000000
     2.5701000690000000    0.0000000000000000    0.0000000000000000
    -1.2850500345000000    2.2257719501000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
               4
Selective dynamics
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000   F   F   F
  0.6666700000000034  0.3333299999999966  0.0985399999999998   F   F   F
  0.3333311322281505  0.6666688677718494  0.1972481653092927   T   T   T
 -0.0000004145140272  0.0000004145140272  0.2951613337135311   T   T   T

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

bulk模型的能量信息

1
2
3
4
[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' ../bulk-calc/OUTCAR
  energy  without entropy=      -14.90258332  energy(sigma->0) =      -14.90537746
  energy  without entropy=      -14.90568039  energy(sigma->0) =      -14.90844732
  energy  without entropy=      -14.90859467  energy(sigma->0) =      -14.91129762
  • A=2.570100069∗2.2257719501−0.0000000000000000∗(−1.2850500345)=5.720456642530275A=2.570100069*2.2257719501-0.0000000000000000*(-1.2850500345)=5.720456642530275 Å2^2
  • Esurfunrelax=−13.9767E_{surf}^{unrelax}=-13.9767 eV
  • Ebulk⋅N=−14.9112E_{bulk}\cdot N=-14.9112 eV
  • Esurfacerelax=−13.9773E_{surface}^{relax}=-13.9773 eV

所以表面能γs\gamma_s = 0.08158 eV/Å2^2 ,又1 eV/Å2^2=16.0218 J/m2^2，所以表面能γs=\gamma_s= 1.3070 J/m2^2

c(1x1)Cu(111)的表面能

之前大师兄提到过，对于面心立方最密堆积的形式，slab模型的建立可以使用原胞也可以使用单胞，那么接下来看看用单胞计算表面能吧

c(1x1)Cu(111)优化计算的能量信息：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
[ctan@baifq-hpc141 conventional-slab-opt]$ grep '  without' OUTCAR
  energy  without entropy=      -55.94972046  energy(sigma->0) =      -55.94434499
  energy  without entropy=      -55.94997517  energy(sigma->0) =      -55.94458201
  energy  without entropy=      -55.95396000  energy(sigma->0) =      -55.94820485
  energy  without entropy=      -55.95519341  energy(sigma->0) =      -55.94928561
  energy  without entropy=      -55.95554732  energy(sigma->0) =      -55.94959634
  energy  without entropy=      -55.95713778  energy(sigma->0) =      -55.95095507
  energy  without entropy=      -55.95788514  energy(sigma->0) =      -55.95158140
  energy  without entropy=      -55.95893404  energy(sigma->0) =      -55.95242747
  energy  without entropy=      -55.95963672  energy(sigma->0) =      -55.95296199
  energy  without entropy=      -55.96031694  energy(sigma->0) =      -55.95340564
  energy  without entropy=      -55.96050719  energy(sigma->0) =      -55.95349958
  energy  without entropy=      -55.96053048  energy(sigma->0) =      -55.95316488
  energy  without entropy=      -55.96067995  energy(sigma->0) =      -55.95346716
  energy  without entropy=      -55.96069855  energy(sigma->0) =      -55.95356576
  energy  without entropy=      -55.96051601  energy(sigma->0) =      -55.95360396
  energy  without entropy=      -55.96050802  energy(sigma->0) =      -55.95361124
  energy  without entropy=      -55.96041586  energy(sigma->0) =      -55.95360050
  energy  without entropy=      -55.96018279  energy(sigma->0) =      -55.95347398

结构信息：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
[ctan@baifq-hpc141 conventional-slab-opt]$ cat CONTCAR
CONTCAR\(1\1\1)
   1.00000000000000
     5.1402001381000000    0.0000000000000000    0.0000000000000000
    -2.5701000690000000    4.4515439000999999    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
              16
Direct
  0.0000006270807907 -0.0000006270807907  0.0011656342599333
  0.6666672048402060  0.3333327951597896  0.0988195745854847
  0.3333327951597896  0.6666672048402060  0.1968004254145135
 -0.0000006270807907  0.0000006270807907  0.2944543657400676
  0.0000006270807907  0.4999993729192098  0.0011656342599333
  0.6666672048402060  0.8333327951597940  0.0988195745854847
  0.3333327951597896  0.1666672048402087  0.1968004254145135
 -0.0000006270807907  0.5000006270807924  0.2944543657400676
  0.5000006270807924  0.4999993729192098  0.0011656342599333
  0.1666672048402087  0.8333327951597940  0.0988195745854847
  0.8333327951597940  0.1666672048402087  0.1968004254145135
  0.4999993729192098  0.5000006270807924  0.2944543657400676
  0.5000006270807924 -0.0000006270807907  0.0011656342599333
  0.1666672048402087  0.3333327951597896  0.0988195745854847
  0.8333327951597940  0.6666672048402060  0.1968004254145135
  0.4999993729192098  0.0000006270807907  0.2944543657400676

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

Cu(111)的体相能量：

1
2
3
4
[ctan@baifq-hpc141 conventional-slab-opt]$ grep '  without' ../bulk-calc/OUTCAR
  energy  without entropy=      -14.90258332  energy(sigma->0) =      -14.90537746
  energy  without entropy=      -14.90568039  energy(sigma->0) =      -14.90844732
  energy  without entropy=      -14.90859467  energy(sigma->0) =      -14.91129762
  • 表面积A=5.1402001381000000∗4.4515439000999999=22.88185.1402001381000000*4.4515439000999999=22.8818 Å2^2
  • EsurfunrelaxE_{surf}^{unrelax} = -55.9443 eV
  • N⋅EbulkN \cdot E_{bulk} = -59.6452 eV
    • N指的是slab模型中的原子个数而非bulk中的原子个数
  • EsurfrelaxE_{surf}^{relax} = -55.9535 eV

计算得表面能γs\gamma_s=0.08047 eV/Å2^2 = 1.2892 J/m2^2

这与p(1x1)Cu(111)计算所得的表面能有一定的差距，目前现阶段我还无法解决这个问题🤔😵‍💫

p(1x1)Ni(100)的表面能

这个例子来自于vasp官方的手册Hands on Session III

1
2
3
4
5
6
7
8
9
10
11
12
general:
SYSTEM = clean Ni(100) surface
ISTART = 0 ; ICHARG=2
ENCUT = 270
ISMEAR = 2 ; SIGMA = 0.2
spin:
ISPIN=2
MAGMOM = 5*1
dynamic:
IBRION = 1
NSW = 100
POTIM = 0.2

有个疑惑的点是手册中的INCAR没有加上偶极矫正，我在自己算的时候加上了偶极矫正

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
System = Nickel-slab-opt
ISMEAR = 1
ISTART = 0
ICHARG = 2
SIGMA = 0.05
ENCUT = 360
ALGO = Normal
ISPIN = 2
MAGMOM = 5*1
NSW = 200
NELM = 200
EDIFF = 1E-6
EDIFFG = -0.01
IDIPOL = 3
LDIPOL = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
IBRION = 1
POTIM = 0.05

体相能量信息：

1
2
3
4
[ctan@baifq-hpc141 primitive-opt]$ grep '  without' ../bulk/OUTCAR
  energy  without entropy=      -21.88865316  energy(sigma->0) =      -21.88759144
  energy  without entropy=      -21.88930785  energy(sigma->0) =      -21.88825056
  energy  without entropy=      -21.89056178  energy(sigma->0) =      -21.88952707

弛豫与未弛豫的slab模型的能量信息：

1
2
3
4
5
6
7
8
[ctan@baifq-hpc141 primitive-opt]$ grep '  without' OUTCAR
  energy  without entropy=      -25.97155386  energy(sigma->0) =      -25.97291371
  energy  without entropy=      -25.97233802  energy(sigma->0) =      -25.97366730
  energy  without entropy=      -25.97412653  energy(sigma->0) =      -25.97534340
  energy  without entropy=      -25.97511519  energy(sigma->0) =      -25.97602399
  energy  without entropy=      -25.97516432  energy(sigma->0) =      -25.97619350
  energy  without entropy=      -25.97516747  energy(sigma->0) =      -25.97619463
  energy  without entropy=      -25.97517140  energy(sigma->0) =      -25.97620156

优化后的结构：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
CONTCAR\(1\0\0)
   1.00000000000000
     2.4846000671000001    0.0000000000000000    0.0000000000000000
    -1.2423000336000001    2.1517267763999999    0.0000000000000000
     0.0000000000000000    0.0000000000000000   18.7047996521000002
  Ni_pv/368cd815
               5
Direct
  0.6666685456044492  0.3333314543955508  0.0005454007373249
  0.0000003831810720 -0.0000003831810720  0.1078816295897802
  0.3333305745814666  0.6666694254185335  0.2169200598414882
  0.6666681017676376  0.3333318982323624  0.3259932620996842
  0.0000023948653782 -0.0000023948653782  0.4332496477317282

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  • 表面积A=2.4846000671000001*2.1517267763999999=5.346180493024307 eV
  • EsurfunrelaxE_{surf}^{unrelax} = -25.97291371 eV
  • EsurfrelaxE_{surf}^{relax} = -25.97620156 eV
  • N⋅EbulkN \cdot E_{bulk} = -21.88952707 eV

表面能 γs\gamma_s = 0.1293 eV/Å2^2 = 2.0715 J/m2^2

小结

  • 本节主要是联系表面能的计算
  • Ni(100)表面能的计算选用的赝势为Ni_sv赝势
从零开始的科研之旅 > VASP
#化学 #科研
表面能的计算
https://hydrogen1222.xyz/2025/04/29/科研/VASP/Starry Sky/表面能的计算/
作者
Storm Talia
发布于
2025年4月29日
许可协议