关于表面能

根据VASP官方手册Hands on Session III,表面能的定义如下,大师兄的文章也提到了这个定义( Nanoscale, 2017, 9, 13089–13094) $$ \sigma=\sigma^{unrelax}+E^{relax}=\dfrac{1}{2}(E_{surf}-N_{atoms}\cdot E_{bulk})+E^{relax} $$

$$ \gamma_s=\dfrac{1}{2A}(E_{surf}^{unrelax}-N\cdot E_{bulk})+\dfrac{1}{A}(E_{surf}^{relax}-E_{surf}^{unrelax}) $$

其中,弛豫能量被定义为: Erelax = Esurfrelax − Esurfunrelax

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

关于表面积A的计算

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

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

对于面积的计算,可设一组不共线的基矢量为a1 = (a11, a12, 0), a2 = (a21, a22, 0)

下面这张图应该挺形象

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

a⃗ × b⃗ = (0, 0, a11 × a22 − a12 × a21)

面积为外积模乘以缩放银子  = s2 ⋅ |a11a22 − a12a21| PS:对于三维的体积  = s3 ⋅ |a1 ⋅ (a2×a3)|

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.720456642530275 Å2
  • Esurfunrelax = −13.9767 eV
  • Ebulk ⋅ N = −14.9112 eV
  • Esurfacerelax = −13.9773 eV

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

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.8818 Å2
  • Esurfunrelax = -55.9443 eV
  • N ⋅ Ebulk = -59.6452 eV
    • N指的是slab模型中的原子个数而非bulk中的原子个数
  • Esurfrelax = -55.9535 eV

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

这与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
  • Esurfunrelax = -25.97291371 eV
  • Esurfrelax = -25.97620156 eV
  • N ⋅ Ebulk = -21.88952707 eV

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

小结

  • 本节主要是联系表面能的计算
  • Ni(100)表面能的计算选用的赝势为Ni_sv赝势

Comment and share

本文为一个查阅性的记录文章,也可以说是搬运,内容取自侯柱锋老师的使用VASP的个人经验手册

元素 普通赝势 (eV) 高硬赝势 (eV) 软赝势 (eV) d 电子作为半芯态来处理赝势 (eV)
B 318 700 (B_h) 250 (B_s) -
C 400 700 (C_h) 273 (C_s) -
N 400 700 (N_h) 250 (N_s) -
O 400 700 (O_h) 250 (O_s) -
F 400 700 (F_h) 250 (F_s) -
Al 240 295 (Al_h) - -
Si 245 380 (Si_h) - -
P 270 390 (P_h) - -
S 280 402 (S_h) - -
Cl 280 409 (Cl_h) - -
Ga 134 404 (Ga_h) - 282 (Ga_d)
Ge 173 410 (Ge_h) - 287 (Ge_d)
As 208 - - -
Se 211 - - -
Br 216 - - -
In 95 - - 239 (In_d)
Sn 103 - - 241 (Sn_d)
Sb 172 - - -
Te 174 - - -
I 175 - - -
Tl 95 - - 239 (Tl_d)
Pb 98 - - 237 (Pb_d)
Bi 105 - - 242 (Bi_d)
元素 X——普通赝势 (eV) X_h——硬赝势 (eV) 其他赝势 (eV)
H 250 700 (H_h) -
Li 140 - 271 (Li_sv)
Be 300 - 308 (Be_sv)
Na 81 - 300 (Na_pv), 700 (Na_sv)
Mg 210 - 265 (Mg_pv)
K - - 150 (K_pv), 259 (K_sv)
Ca - - 150 (Ca_pv), 290 (Ca_sv)
Rb - - 121 (Rb_pv), 220 (Rb_sv)
Sr - - 226 (Sr_sv)
Cs - - 220 (Cs_sv)
Ba - - 187 (Ba_sv)
元素 普通赝势 (eV) X_pv赝势 (eV) X_sv赝势 (eV)
Sc - - 222 (Sc_sv)
Ti 178 222 (Ti_pv) -
V 192 263 (V_pv) -
Cr 227 265 (Cr_pv) -
Mn 269 269 (Mn_pv) -
Fe 267 293 (Fe_pv) -
Co 267 - -
Ni 269 367 (Ni_pv) -
Cu 273 368 (Cu_pv) -
Zn 276 - -
Y - - 211 (Y_sv)
Zr - - 229 (Zr_sv)
Nb - 207 (Nb_pv) -
Mo 224 224 (Mo_pv) -
Tc 228 228 (Tc_pv) -
Ru 213 230 (Ru_pv) -
Rh 228 271 (Rh_pv) -
Pd 250 350 (Pd_pv) -
Ag 249 - -
Cd 274 - -
Hf 220 220 (Hf_pv) -
Ta 223 223 (Ta_pv) -
W 223 223 (W_pv) -
Re 226 226 (Re_pv) -
Os 228 228 (Os_pv) -
Ir 210 - -
Pt 230 - -
Au 229 - -
Hg 233 - -

Comment and share

本文其实没啥新的内容,只是对上两节的一个总结分析

在表面弛豫前(slab模型优化计算前)

1
2
3
4
5
6
7
8
9
10
11
12
13
SYSTEM = p(1x1)Cu(111)
1.0
2.5701000690 0.0000000000 0.0000000000
-1.2850500345 2.2257719501 0.0000000000
0.0000000000 0.0000000000 21.2954006195
Cu
4
Selective Dynamics
Direct
0.000000000 0.000000000 0.000000000 F F F
0.666670000 0.333330000 0.098540000 F F F
0.333330000 0.666670000 0.197080000 T T T
0.000000000 0.000000000 0.295620000 T T T

优化计算后

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
SYSTEM = p(1x1)Cu(111)
1.0
2.5701000690 0.0000000000 0.0000000000
-1.2850500345 2.2257719501 0.0000000000
0.0000000000 0.0000000000 21.2954006195
Cu
4
Selective Dynamics
Direct
0.000000000 0.000000000 0.000000000 F F F
0.666670000 0.333330000 0.098540000 F F F
0.333330000 0.666670000 0.197080000 T T T
0.000000000 0.000000000 0.295620000 T T T
[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat CONTCAR
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

工作站离线状态下numpy on python2地安装

这两个文件中的坐标都是分数坐标,不太利于分析位置变化,需要将分数坐标转化为笛卡尔坐标,使用大师兄的脚本来进行转换,

用的是python2

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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#Convert direc coordiation to cartesian Writen By Qiang and Xu Nan
"""
The right input way to run this script is :

python dire2cart.py CONTCAR 3

CONTCAR will be converted from Direct to Cartesian and 3 layers will be fixed.

If you do not want to fix layers and keep the same as before, run command:

python dire2cat.py CONTCAR

"""

import sys
import os
import numpy as np

print 'Please read the head part of this script and get more information!'
print """
###################################
# #
#for VASP 5.2 or higher versions #
# #
###################################
"""

if len(sys.argv) <= 1:
print('\n' + ' Warning! ' * 3 + '\n')
print('You did not select the inputfile to be converted. \n By defalut, we are going to convert your CONTCAR.\n')

if not os.path.isfile("POSCAR") and not os.path.isfile("CONTCAR"):
print("Error:" * 3 + "\n Can not find neither POSCAR nor CONTCAR!!\n")
exit()
else:
if os.path.isfile("CONTCAR"):
script = sys.argv[0]
file_to_be_converted = "CONTCAR"
print "\n Convertion starts........\n"

else:
script = sys.argv[0]
file_to_be_converted = "POSCAR"
print "\n There is no CONTCAR in your directory. \n \n POSCAR Convertion starts........\n"

if len(sys.argv) == 2:
print("\n%s Convertion starts......" %sys.argv[1])
script, file_to_be_converted = sys.argv[:2]
else:
print("\n%s Convertion starts......\n" %sys.argv[1])
script, file_to_be_converted, fixedlayer = sys.argv[:3]
fixedlayer = int(fixedlayer)

def get_infor():
f = open(file_to_be_converted, 'r')
lines = f.readlines()
f.close()
num_atoms = sum([int(x) for x in lines[6].split()])
if lines[7][0] == 'S' or lines[7][0] == 's': # # With Selected T T T, coordination starts from line 9
start_num = 9
if lines[8][0] == 'D' or lines[8][0] == 'd':
is_direct = True
else:
is_direct = False
else:
start_num = 8
print '----------------------------------------------------'
print 'Pay Attetion! There is no TTT in %s ' %(file_to_be_converted)
print '----------------------------------------------------'
if lines[7][0] == 'D' or lines[7][0] == 'd':
is_direct = True
else:
is_direct = False
a = []
b = []
c = []
if is_direct:
for i in np.arange(2,5):
line = [float(i) for i in lines[i].strip().split()]
a.append(line[0])
b.append(line[1])
c.append(line[2])
vector = np.array([a,b,c])
else:
vector = np.array([[1, 0 , 0], [0, 1, 0], [0, 0, 1]])
return vector, lines, start_num, num_atoms, is_direct

def determinelayers(z_cartesian):
threshold = 0.5
seq = sorted(z_cartesian)
min = seq[0]
layerscount = {}
sets = [min]
for j in range(len(seq)):
if abs(seq[j]-min) >= threshold:
min = seq[j]
sets.append(min)

for i in range(1,len(sets)+1):
layerscount[i] = []
for k in range(len(z_cartesian)):
if abs(z_cartesian[k]-sets[i-1]) <= threshold:
layerscount[i].append(k)

return layerscount

def convert():
x_cartesian = []
y_cartesian = []
z_cartesian = []
tf = []
for i in range(start_num, num_atoms + start_num):
line_data = [float(ele) for ele in lines[i].split()[0:3]]
line_data = np.array([line_data])
x, y, z = [sum(k) for k in line_data * vector ]
x_cartesian.append(x)
y_cartesian.append(y)
z_cartesian.append(z)

if start_num == 9 : # if T T T exist, the start_num will be 9
tf.append((lines[i].split()[3:]))
else:
tf.append(' ') # if there is no T T T, use space instead.

layerscount =determinelayers(z_cartesian)

file_out = open(file_to_be_converted+'_C', 'w')
for i in range(0,7):
file_out.write(lines[i].rstrip() + '\n') # first 7 lines are kept the same

print '\n Find %s layers!' * 3 %(len(layerscount), len(layerscount), len(layerscount))

if len(sys.argv) >= 3: # This means that the nuber for fixing layers is given by the user.
file_out.write('Selective\n')
file_out.write('Cartesian' + '\n')
for i in range(1,len(layerscount)+1):
if i <= fixedlayer:
for j in layerscount[i]:
tf[j] = ['F','F','F']
else:
for k in layerscount[i]:
tf[k] = ['T','T','T']
else:
if start_num == 9: # 9 means there are T T T or F F F in the file
file_out.write('Selective\n')
file_out.write('Cartesian' + '\n')
else:
file_out.write('Cartesian' + '\n')

for i in range(0,len(x_cartesian)):
file_out.write("\t%+-3.10f %+-3.10f %+-3.10f %s\n" %(x_cartesian[i], y_cartesian[i], z_cartesian[i], ' '.join(tf[i])))

file_out.close()

vector, lines, start_num, num_atoms, is_direct = get_infor()

if is_direct :
print "\n%s has Direct Coordinates, Contersion starts.... " %(file_to_be_converted)
convert()
else:
print "\n%s has Cartesian Coordinates Already! We are going to fix layers only." %(file_to_be_converted)
convert()

print '-----------------------------------------------------\n'
print '\n %s with Cartesian Coordiates is named as %s_C\n' %(file_to_be_converted, file_to_be_converted)
print '-----------------------------------------------------\n'

脚本会调用到numpy库,然而集群上并没有连接互联网,所以没有办法通过pip2 install numpy来安装,只能手动上传.whl文件或者源码编译等。

进入python2环境,输入import pip;print(pip.pep425tags.get_supported())来查看支持的.whl文件版本

1
2
3
4
5
6
[ctan@baifq-hpc141 primitive-slab-opt-1]$ python2
Python 2.7.18 (default, Oct 23 2023, 20:01:03)
[GCC 8.5.0 20210514 (Red Hat 8.5.0-18)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import pip;print(pip.pep425tags.get_supported())
[('cp27', 'cp27mu', 'manylinux1_x86_64'), ('cp27', 'cp27mu', 'linux_x86_64'), ('cp27', 'none', 'manylinux1_x86_64'), ('cp27', 'none', 'linux_x86_64'), ('py2', 'none', 'manylinux1_x86_64'), ('py2', 'none', 'linux_x86_64'), ('cp27', 'none', 'any'), ('cp2', 'none', 'any'), ('py27', 'none', 'any'), ('py2', 'none', 'any'), ('py26', 'none', 'any'), ('py25', 'none', 'any'), ('py24', 'none', 'any'), ('py23', 'none', 'any'), ('py22', 'none', 'any'), ('py21', 'none', 'any'), ('py20', 'none', 'any')]

这里下载的是numpy-1.16.5-cp27-cp27mu-manylinux1_x86_64.whl(不可下载过高版本的numpy,它支持python2的最后一个版本似乎是1.16.6

然后安装,由于我不在sudoer中,没有root权限,所以需要手动指定安装目录

1
pip2 install --prefix=~/python2/numpy_installation numpy-1.16.5-cp27-cp27mu-manylinux1_x86_64.whl

然后添加环境变量,注意numpy库在/lib64/python2.7/site-packages目录中的

1
export PYTHONPATH=/home/ctan/python2/numpy_installation/lib64/python2.7/site-packages:$PYTHONPATH

ok啦,以上是题外话,安装完numpy之后就可以运行脚本了,将POSCARCONTCAR的坐标从分数坐标转换为笛卡尔坐标

POSCAR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat POSCAR_C
CONTCAR\(1\1\1)
1.0
2.5701000690 0.0000000000 0.0000000000
-1.2850500345 2.2257719501 0.0000000000
0.0000000000 0.0000000000 21.2954006195
Cu
4
Selective
Cartesian
+0.0000000000 +0.0000000000 +0.0000000000 F F F
+1.2850628850 +0.7419165641 +2.0984487770 F F F
-0.0000128505 +1.4838553860 +4.1968975541 T T T
+0.0000000000 +0.0000000000 +6.2953463311 T T T

CONTCAR_C

1
2
3
4
5
6
7
8
9
10
11
12
13
14
[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat CONTCAR_C
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
Cartesian
+0.0000000000 +0.0000000000 +0.0000000000 F F F
+1.2850628850 +0.7419165641 +2.0984487770 F F F
-0.0000084856 +1.4838528659 +4.2004787017 T T T
-0.0000015980 +0.0000009226 +6.2855788488 T T T

Alright, I’m back. 突然想起昨天球磨的样品今天结束了,连忙赶回实验室取出来装进样品管里等待封管。

By the way,我真是和鱼一样,记忆只有七秒、、、😅。师兄教的球磨罐和球珠的洗涤方法忘记下来了,今天不出意外地给忘了😂。汲取教训,每次出实验室后都更加努力地总结和复盘,仔细回想每一个实验细节

优化前后对比

Alright,扯远了,分析一下顶层和亚顶层优化前后地结构变化吧:

  • POSCAR: (+6.2953463311) - (+4.1968975541) = 2.0984487769999998
  • CONTCAR: (+6.2855788488) - (+4.2004787017) = 2.0851001471000004

前后变化了:(2.0851001471000004-2.0984487769999998)/2.0984487769999998 = -0.006361189296735172

可以看到,表层原子向体相收缩

能量变化

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

能量从-13.97673168变化为-13.97735013,降低了0.000618449999999271 eV,即弛豫过程是放热,刚切开的表面不稳定,表层原子向体相收缩后,体系能量降低,变得更稳定

小Tips:可以使用python环境便捷地进行基本运算

python环境下进行基本运算

Comment and share

POSCAR的修改

与单点计算不同,优化计算可能涉及到离子步迭代、晶格常数的改变。但Cu(111)晶面的优化计算并==不只是==在INCAR里添加参数IBRION,优化计算确实允许原子弛豫,但在slab模型中,只允许表面的两层原子允许弛豫,而下面的那些原子则直接固定住

这样做的物理意义在于:在真实的环境下,催化剂体相看做是不变的,只有表面的原子参与催化反应

前面在乙醇分子的振动频率分析里面已经学过了,添加以下选择性动力学标记Slective Dynamics即可

PS:看分数坐标即可知道哪些是表层哪些是体相

1
2
3
4
5
6
7
8
9
10
11
12
13
CONTCAR\(1\1\1)
1.0
2.5701000690 0.0000000000 0.0000000000
-1.2850500345 2.2257719501 0.0000000000
0.0000000000 0.0000000000 21.2954006195
Cu
4
Selective Dynamics
Direct
0.000000000 0.000000000 0.000000000 F F F
0.666670000 0.333330000 0.098540000 F F F
0.333330000 0.666670000 0.197080000 T T T
0.000000000 0.000000000 0.295620000 T T T

INCAR的修改

现在是优化计算,就得加上算法了,大师兄用的IBRION = 2,那我就对比一下,用IBRION = 1试一下吧

IBRION = 1

官方wiki说准牛顿迭代法对POTIM比较敏感,那我就调小一点吧

1
2
3
4
5
6
7
8
9
10
11
12
13
SYSTEM = Cu-concentional-slab-statistic
ISMEAR = 1
SIGMA = 0.5
ALGO = Normal
ENCUT = 700
EDIFF = 1E-6
LDIPOL = .TRUE.
IDIPOL = 3
NCORE = 20
NSW = 500
NELM = 200
IBRION = 1
POTIM = 0.05

成功收敛,迭代情况

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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
[ctan@baifq-hpc141 primitive-slab-opt]$ cat OSZICAR
N E dE d eps ncg rms rms(c)
DAV: 1 0.380094340720E+03 0.38009E+03 -0.23566E+04 2548 0.269E+03
DAV: 2 0.163549324393E+02 -0.36374E+03 -0.35588E+03 2548 0.438E+02
DAV: 3 -0.154087935417E+02 -0.31764E+02 -0.30239E+02 2816 0.186E+02
DAV: 4 -0.171245615883E+02 -0.17158E+01 -0.16227E+01 2856 0.463E+01
DAV: 5 -0.171794690835E+02 -0.54907E-01 -0.54864E-01 2954 0.690E+00 0.146E+01
DAV: 6 -0.156158003013E+02 0.15637E+01 -0.42034E+01 2712 0.133E+02 0.829E+00
DAV: 7 -0.140078760986E+02 0.16079E+01 -0.79564E+00 2992 0.596E+01 0.790E-01
DAV: 8 -0.140167702165E+02 -0.88941E-02 -0.53048E-01 3968 0.226E+01 0.122E+00
DAV: 9 -0.140083290175E+02 0.84412E-02 -0.45257E-01 3848 0.291E+01 0.104E+00
DAV: 10 -0.139869448329E+02 0.21384E-01 -0.20408E-01 4800 0.187E+01 0.615E-01
DAV: 11 -0.139810722283E+02 0.58726E-02 -0.60166E-02 4016 0.780E+00 0.113E-01
DAV: 12 -0.139806486192E+02 0.42361E-03 -0.10641E-02 3786 0.272E+00 0.980E-02
DAV: 13 -0.139810968626E+02 -0.44824E-03 -0.87951E-03 4868 0.367E+00 0.844E-02
DAV: 14 -0.139803135502E+02 0.78331E-03 -0.42901E-03 4920 0.238E+00 0.308E-02
DAV: 15 -0.139802543403E+02 0.59210E-04 -0.32281E-04 3540 0.324E-01 0.182E-02
DAV: 16 -0.139803209356E+02 -0.66595E-04 -0.83429E-04 4558 0.937E-01 0.301E-02
DAV: 17 -0.139803228923E+02 -0.19567E-05 -0.15361E-04 3074 0.264E-01 0.168E-02
DAV: 18 -0.139803417113E+02 -0.18819E-04 -0.36368E-04 4750 0.744E-01 0.954E-03
DAV: 19 -0.139803412500E+02 0.46126E-06 -0.10659E-04 2728 0.404E-01 0.517E-03
DAV: 20 -0.139803676358E+02 -0.26386E-04 -0.80023E-05 3638 0.357E-01 0.818E-03
DAV: 21 -0.139803787022E+02 -0.11066E-04 -0.13797E-05 3336 0.920E-02 0.101E-02
DAV: 22 -0.139803765945E+02 0.21077E-05 -0.39937E-05 2614 0.265E-01 0.222E-03
DAV: 23 -0.139803888096E+02 -0.12215E-04 -0.60877E-05 2528 0.297E-01 0.725E-03
DAV: 24 -0.139803907613E+02 -0.19517E-05 -0.71066E-06 3156 0.553E-02 0.611E-03
DAV: 25 -0.139803916036E+02 -0.84239E-06 -0.13171E-06 1634 0.198E-02
1 F= -.13980392E+02 E0= -.13976732E+02 d E =-.139804E+02
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139807173282E+02 -0.32657E-03 0.14112E+02 2602 0.612E+00 0.467E-03
DAV: 2 -0.139804390721E+02 0.27826E-03 -0.16468E-04 4156 0.597E-01 0.165E-02
DAV: 3 -0.139804261532E+02 0.12919E-04 -0.10604E-04 3756 0.481E-01 0.167E-03
DAV: 4 -0.139804267349E+02 -0.58173E-06 -0.46237E-06 2500 0.225E-02
2 F= -.13980427E+02 E0= -.13976770E+02 d E =-.351312E-04
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139804330907E+02 -0.69375E-05 0.16389E+02 2632 0.107E+01 0.452E-02
DAV: 2 -0.139809280760E+02 -0.49499E-03 -0.57833E-03 3892 0.164E+00 0.586E-02
DAV: 3 -0.139808976764E+02 0.30400E-04 -0.10914E-03 4308 0.154E+00 0.348E-02
DAV: 4 -0.139808745784E+02 0.23098E-04 -0.21714E-04 4572 0.607E-01 0.278E-02
DAV: 5 -0.139808671461E+02 0.74323E-05 -0.56092E-05 4164 0.315E-01 0.187E-02
DAV: 6 -0.139808630630E+02 0.40832E-05 -0.44161E-05 4470 0.289E-01 0.698E-03
DAV: 7 -0.139808615417E+02 0.15213E-05 -0.76178E-06 1946 0.121E-01 0.290E-03
DAV: 8 -0.139808615453E+02 -0.35725E-08 -0.25204E-06 1800 0.234E-02
3 F= -.13980862E+02 E0= -.13977262E+02 d E =-.434810E-03
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139808810490E+02 -0.19507E-04 0.54240E+01 2632 0.376E+00 0.136E-02
DAV: 2 -0.139809351461E+02 -0.54097E-04 -0.54142E-04 3890 0.509E-01 0.173E-02
DAV: 3 -0.139809327068E+02 0.24393E-05 -0.10930E-04 4292 0.484E-01 0.104E-02
DAV: 4 -0.139809303964E+02 0.23104E-05 -0.18479E-05 2078 0.186E-01 0.828E-03
DAV: 5 -0.139809300877E+02 0.30864E-06 -0.87511E-06 2780 0.122E-01
4 F= -.13980930E+02 E0= -.13977356E+02 d E =-.685425E-04
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139808203045E+02 0.11009E-03 0.13819E+02 3808 0.576E+00 0.330E-03
DAV: 2 -0.139808743367E+02 -0.54032E-04 0.77025E+00 2180 0.218E+00 0.818E-03
DAV: 3 -0.139809353212E+02 -0.60985E-04 -0.25088E-05 2014 0.239E-01 0.103E-03
DAV: 4 -0.139809353447E+02 -0.23522E-07 -0.62565E-07 1516 0.240E-02
5 F= -.13980935E+02 E0= -.13977350E+02 d E =-.525697E-05

能量:

1
2
3
4
5
6
[ctan@baifq-hpc141 primitive-slab-opt]$ 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

计算时间:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
General timing and accounting informations for this job:
========================================================

Total CPU time used (sec): 53.453
User time (sec): 49.556
System time (sec): 3.896
Elapsed time (sec): 55.466

Maximum memory used (kb): 134456.
Average memory used (kb): N/A

Minor page faults: 12550
Major page faults: 686
Voluntary context switches: 1587

IBRION =2

1
2
3
4
5
6
7
8
9
10
11
12
13
SYSTEM = Cu-concentional-slab-statistic
ISMEAR = 1
SIGMA = 0.5
ALGO = Normal
ENCUT = 700
EDIFF = 1E-6
LDIPOL = .TRUE.
IDIPOL = 3
NCORE = 20
NSW = 500
NELM = 200
IBRION = 2
POTIM = 0.1

成功收敛,迭代情况:

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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV: 1 0.380094340720E+03 0.38009E+03 -0.23566E+04 2548 0.269E+03
DAV: 2 0.163549324393E+02 -0.36374E+03 -0.35588E+03 2548 0.438E+02
DAV: 3 -0.154087935417E+02 -0.31764E+02 -0.30239E+02 2816 0.186E+02
DAV: 4 -0.171245615883E+02 -0.17158E+01 -0.16227E+01 2856 0.463E+01
DAV: 5 -0.171794690835E+02 -0.54907E-01 -0.54864E-01 2954 0.690E+00 0.146E+01
DAV: 6 -0.156158003013E+02 0.15637E+01 -0.42034E+01 2712 0.133E+02 0.829E+00
DAV: 7 -0.140078760986E+02 0.16079E+01 -0.79564E+00 2992 0.596E+01 0.790E-01
DAV: 8 -0.140167702165E+02 -0.88941E-02 -0.53048E-01 3968 0.226E+01 0.122E+00
DAV: 9 -0.140083290175E+02 0.84412E-02 -0.45257E-01 3848 0.291E+01 0.104E+00
DAV: 10 -0.139869448329E+02 0.21384E-01 -0.20408E-01 4800 0.187E+01 0.615E-01
DAV: 11 -0.139810722283E+02 0.58726E-02 -0.60166E-02 4016 0.780E+00 0.113E-01
DAV: 12 -0.139806486192E+02 0.42361E-03 -0.10641E-02 3786 0.272E+00 0.980E-02
DAV: 13 -0.139810968626E+02 -0.44824E-03 -0.87951E-03 4868 0.367E+00 0.844E-02
DAV: 14 -0.139803135502E+02 0.78331E-03 -0.42901E-03 4920 0.238E+00 0.308E-02
DAV: 15 -0.139802543403E+02 0.59210E-04 -0.32281E-04 3540 0.324E-01 0.182E-02
DAV: 16 -0.139803209356E+02 -0.66595E-04 -0.83429E-04 4558 0.937E-01 0.301E-02
DAV: 17 -0.139803228923E+02 -0.19567E-05 -0.15361E-04 3074 0.264E-01 0.168E-02
DAV: 18 -0.139803417113E+02 -0.18819E-04 -0.36368E-04 4750 0.744E-01 0.954E-03
DAV: 19 -0.139803412500E+02 0.46126E-06 -0.10659E-04 2728 0.404E-01 0.517E-03
DAV: 20 -0.139803676358E+02 -0.26386E-04 -0.80023E-05 3638 0.357E-01 0.818E-03
DAV: 21 -0.139803787022E+02 -0.11066E-04 -0.13797E-05 3336 0.920E-02 0.101E-02
DAV: 22 -0.139803765945E+02 0.21077E-05 -0.39937E-05 2614 0.265E-01 0.222E-03
DAV: 23 -0.139803888096E+02 -0.12215E-04 -0.60877E-05 2528 0.297E-01 0.725E-03
DAV: 24 -0.139803907613E+02 -0.19517E-05 -0.71066E-06 3156 0.553E-02 0.611E-03
DAV: 25 -0.139803916036E+02 -0.84239E-06 -0.13171E-06 1634 0.198E-02
1 F= -.13980392E+02 E0= -.13976732E+02 d E =-.139804E+02
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139807447433E+02 -0.35398E-03 0.14111E+02 2590 0.632E+00 0.528E-03
DAV: 2 -0.139804678790E+02 0.27686E-03 -0.17889E-04 4964 0.563E-01 0.155E-02
DAV: 3 -0.139804572701E+02 0.10609E-04 -0.99595E-05 4260 0.455E-01 0.347E-03
DAV: 4 -0.139804576704E+02 -0.40025E-06 -0.38794E-06 2444 0.252E-02
2 F= -.13980458E+02 E0= -.13976804E+02 d E =-.660667E-04
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139806303531E+02 -0.17308E-03 0.10590E+01 2646 0.165E+00 0.137E-02
DAV: 2 -0.139806282211E+02 0.21320E-05 -0.51586E-04 3892 0.441E-01 0.161E-02
DAV: 3 -0.139806267938E+02 0.14273E-05 -0.86921E-05 4092 0.433E-01 0.112E-02
DAV: 4 -0.139806245851E+02 0.22087E-05 -0.18757E-05 2156 0.193E-01 0.871E-03
DAV: 5 -0.139806244737E+02 0.11138E-06 -0.11052E-05 3344 0.135E-01 0.639E-03
DAV: 6 -0.139806241801E+02 0.29361E-06 -0.44497E-06 1850 0.106E-01
3 F= -.13980624E+02 E0= -.13976988E+02 d E =-.232576E-03
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139806399034E+02 -0.15430E-04 0.10673E+02 2632 0.655E+00 0.273E-02
DAV: 2 -0.139808756284E+02 -0.23573E-03 -0.22093E-03 3886 0.102E+00 0.353E-02
DAV: 3 -0.139808651090E+02 0.10519E-04 -0.43844E-04 4364 0.978E-01 0.208E-02
DAV: 4 -0.139808558492E+02 0.92598E-05 -0.82890E-05 3866 0.380E-01 0.163E-02
DAV: 5 -0.139808537696E+02 0.20796E-05 -0.22746E-05 3796 0.195E-01 0.113E-02
DAV: 6 -0.139808525623E+02 0.12073E-05 -0.15366E-05 3678 0.171E-01 0.403E-03
DAV: 7 -0.139808519218E+02 0.64052E-06 -0.31561E-06 1528 0.772E-02
4 F= -.13980852E+02 E0= -.13977253E+02 d E =-.460318E-03
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139800417341E+02 0.81083E-03 0.23560E+01 2630 0.588E+00 0.573E-02
DAV: 2 -0.139809885034E+02 -0.94677E-03 -0.91190E-03 3890 0.208E+00 0.745E-02
DAV: 3 -0.139809389586E+02 0.49545E-04 -0.17575E-03 4372 0.196E+00 0.428E-02
DAV: 4 -0.139809007116E+02 0.38247E-04 -0.35205E-04 4684 0.779E-01 0.339E-02
DAV: 5 -0.139808904102E+02 0.10301E-04 -0.81050E-05 4124 0.377E-01 0.239E-02
DAV: 6 -0.139808824251E+02 0.79850E-05 -0.67930E-05 4464 0.342E-01 0.861E-03
DAV: 7 -0.139808804209E+02 0.20042E-05 -0.13063E-05 2516 0.156E-01 0.364E-03
DAV: 8 -0.139808803022E+02 0.11871E-06 -0.29028E-06 1874 0.341E-02
5 F= -.13980880E+02 E0= -.13977356E+02 d E =-.488699E-03
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139807672910E+02 0.11313E-03 0.15468E+02 2630 0.858E+00 0.205E-02
DAV: 2 -0.139809519135E+02 -0.18462E-03 -0.16290E-03 3860 0.606E-01 0.224E-02
DAV: 3 -0.139809495045E+02 0.24089E-05 -0.16101E-04 3760 0.572E-01 0.195E-02
DAV: 4 -0.139809446950E+02 0.48095E-05 -0.34715E-05 3338 0.266E-01 0.149E-02
DAV: 5 -0.139809410051E+02 0.36899E-05 -0.17275E-05 3888 0.130E-01 0.678E-03
DAV: 6 -0.139809401472E+02 0.85792E-06 -0.69054E-06 2716 0.164E-01
6 F= -.13980940E+02 E0= -.13977383E+02 d E =-.548544E-03
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139808959432E+02 0.45062E-04 0.81868E+01 2640 0.491E+00 0.921E-03
DAV: 2 -0.139810234583E+02 -0.12752E-03 -0.40347E-04 3868 0.335E-01 0.124E-02
DAV: 3 -0.139810222522E+02 0.12061E-05 -0.37434E-05 3890 0.257E-01 0.934E-03
DAV: 4 -0.139810216716E+02 0.58055E-06 -0.56497E-06 1910 0.105E-01
7 F= -.13981022E+02 E0= -.13977454E+02 d E =-.815244E-04
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139809168510E+02 0.10540E-03 0.16544E+02 2640 0.100E+01 0.295E-02
DAV: 2 -0.139812498696E+02 -0.33302E-03 -0.35344E-03 3874 0.942E-01 0.353E-02
DAV: 3 -0.139812400259E+02 0.98437E-05 -0.31428E-04 3728 0.768E-01 0.285E-02
DAV: 4 -0.139812339348E+02 0.60911E-05 -0.66313E-05 3954 0.364E-01 0.240E-02
DAV: 5 -0.139812229797E+02 0.10955E-04 -0.53085E-05 4306 0.249E-01 0.134E-02
DAV: 6 -0.139812204144E+02 0.25653E-05 -0.20199E-05 3228 0.260E-01 0.363E-03
DAV: 7 -0.139812202373E+02 0.17704E-06 -0.19687E-06 1482 0.408E-02
8 F= -.13981220E+02 E0= -.13977623E+02 d E =-.280090E-03
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139799374972E+02 0.12829E-02 0.13429E+02 2638 0.119E+01 0.609E-02
DAV: 2 -0.139815536956E+02 -0.16162E-02 -0.15321E-02 3874 0.237E+00 0.834E-02
DAV: 3 -0.139814840776E+02 0.69618E-04 -0.19373E-03 4340 0.194E+00 0.553E-02
DAV: 4 -0.139814504163E+02 0.33661E-04 -0.33985E-04 4524 0.775E-01 0.477E-02
DAV: 5 -0.139814168287E+02 0.33588E-04 -0.15051E-04 4446 0.480E-01 0.325E-02
DAV: 6 -0.139813975448E+02 0.19284E-04 -0.12765E-04 4566 0.480E-01 0.853E-03
DAV: 7 -0.139813962132E+02 0.13316E-05 -0.15668E-05 3464 0.152E-01 0.507E-03
DAV: 8 -0.139813960169E+02 0.19635E-06 -0.22934E-06 1728 0.311E-02
9 F= -.13981396E+02 E0= -.13977735E+02 d E =-.455870E-03
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139815084243E+02 -0.11221E-03 0.83969E+01 2670 0.421E+00 0.434E-03
DAV: 2 -0.139813984003E+02 0.11002E-03 -0.73291E-05 3840 0.104E-01 0.437E-03
DAV: 3 -0.139813981891E+02 0.21124E-06 -0.11427E-06 1510 0.490E-02
10 F= -.13981398E+02 E0= -.13977729E+02 d E =-.217219E-05
N E dE d eps ncg rms rms(c)
DAV: 1 -0.139815087603E+02 -0.11036E-03 0.15888E+02 2670 0.681E+00 0.323E-03
DAV: 2 -0.139813984696E+02 0.11029E-03 -0.26471E-05 2428 0.210E-01 0.588E-03
DAV: 3 -0.139813976417E+02 0.82788E-06 -0.27049E-05 3780 0.204E-01 0.141E-03
DAV: 4 -0.139813978006E+02 -0.15885E-06 -0.17146E-06 1590 0.262E-02
11 F= -.13981398E+02 E0= -.13977732E+02 d E =-.178370E-05

计算时间:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
General timing and accounting informations for this job:
========================================================

Total CPU time used (sec): 93.113
User time (sec): 87.076
System time (sec): 6.036
Elapsed time (sec): 94.073

Maximum memory used (kb): 126136.
Average memory used (kb): N/A

Minor page faults: 11497
Major page faults: 687
Voluntary context switches: 1456

能量信息:

1
2
3
4
5
6
7
8
9
10
11
12
[ctan@baifq-hpc141 primitive-slab-opt-2]$ grep '  without' OUTCAR
energy without entropy= -13.96941183 energy(sigma->0) = -13.97673168
energy without entropy= -13.96949582 energy(sigma->0) = -13.97680372
energy without entropy= -13.96971572 energy(sigma->0) = -13.97698803
energy without entropy= -13.97005607 energy(sigma->0) = -13.97725330
energy without entropy= -13.97030617 energy(sigma->0) = -13.97735559
energy without entropy= -13.97026927 energy(sigma->0) = -13.97738319
energy without entropy= -13.97031853 energy(sigma->0) = -13.97745396
energy without entropy= -13.97042910 energy(sigma->0) = -13.97762319
energy without entropy= -13.97041285 energy(sigma->0) = -13.97773496
energy without entropy= -13.97039032 energy(sigma->0) = -13.97772890
energy without entropy= -13.97040171 energy(sigma->0) = -13.97773244

对比一下发现,正如wiki所说,IBRION = 1适合初始结构较好的情况,IBRION - 2的适用性更广,更可靠,尽管IBRION = 2POTIM更大,但所用时间更长。二者最终的能量相差不大

Comment and share

首先对Cu单胞进行结构优化(bulk计算),然后使用优化的结构进行切割,再进行计算(由于晶胞已经优化过,所以本次计算其实是单点计算)。大师兄提到:

Cu(111) slab 模型(POSCAR)的制作流程(复习一下)。注意:前面我们用的是Conventional cell,下面用的是Primitive Cell。 一般来说用Conventional cell, FCC的金属可以用Primitive Cell 但对于其他体系,通过Primitive Cell切出来表面模型有问题

我准备两种都验证一下。

使用Conventional Cell建立slab模型

切割结构

将bulk计算完毕的CONTCAR文件用Vesta打开,再转换为.cif文件使用Materials Studio打开

  • Step1:Build-Surface-Cleave Surface

Angstrom是单位,需要切割的晶面是(111),相应的平面方程为x + y + z − 1 = 0,使用点到直线的距离公式即可求出晶面间距离,,即厚度thickness

PS:之前111晶面是啥都给整忘了,算个点到直线的距离也能算错,真该抽自己两巴掌回去恶补高等数学😭

以下是来自维基百科的关于米勒指数以及111晶面的介绍,很形象了

米勒指数以及111晶面
GeoGebra计算器

平面ABC为:x + y + z − 1 = 0,点D为(1,1,0)

算一下厚度: $$ d=\dfrac{|Ax_0+By_0+Cz_0|}{\sqrt{A^2+B^2+C^2}}=\dfrac{1}{\sqrt{3}}a\approx \dfrac{1}{\sqrt{3}}\times 3.63468 \approx 2.0984 $$ Cleave Surface in Materials Studio

在切割晶面时,thickness设置为4.0 (4.0 = 1.0 × 4),即四倍的晶胞参数、四层原子,按一下TAB,软件软件会自动计算四层原子的真是厚度(单位:Angstrom埃),大约是8.394 ≈ 4 × 2.0984,这与我们手动计算的晶面间距离2.0984的四倍相符

  • Step2:添加真空层,Build-Crystal-Build vacuum slab-Vacuum orientation:C
添加真空层
  • Step3:右键-lattice parameters
设置Z轴方向为真空层

导出为.cif,vesta打开转换为POSCAR,

INCAR如下

1
2
3
4
5
6
7
8
9
10
11
12
System = Cu-conventional-slab
ISMEAR = 1
SIGMA = 0.1
ALGO = Fast
ENCUT = 450 # the system is too huge to adopt huge ENCUT
EDIFF = 1E-5
LDIPOL = .TRUE.
IDIPOL = 3
NWRITE = 0
LWAVE = .FALSE.
LCHARG = .FALSE.
NCORE = 20

Slab模型虽然是代表表面,但是实际上在z方向是固体-真空-固体-真空-…的交替。如果我们建立的Slab模型在z方向是非对称的,模型就会产生一个沿z方向的偶极。偶极会产生静电势,静电势接着会影响模型的镜像(周期性边界条件)。最后算出来的模型的总能量和力与真实情况是不相符的。因此我们需要方法去矫正这种虚假的静电影响。

一种常用的方法就是偶极矫正,在真空部分加入一个超窄的但是方向相反的偶极。这样一来,固体模型产生的偶极和真空中的偶极就会相互抵消。模型和其镜像之间的静电势影响就会抵消

微调KPOINTS

1
2
3
4
5
K-POINTS
0
Gamma
13 13 1
0 0 0

计算结果

1
2
[ctan@baifq-hpc141 Cu-conventional-slab]$ grep '  without' OUTCAR
energy without entropy= -55.91176854 energy(sigma->0) = -55.91397643

PS:最近学迷糊了,关于收敛脑海中有一个极其不正且的理解

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
[ctan@baifq-hpc141 bulk-calc]$ grep accuracy OUTCAR
reached required accuracy - stopping structural energy minimisation
[ctan@baifq-hpc141 bulk-calc]$ cat OSZICAR
N E dE d eps ncg rms rms(c)
DAV: 1 0.103476848532E+03 0.10348E+03 -0.27847E+04 6720 0.246E+03
DAV: 2 -0.180387697818E+02 -0.12152E+03 -0.12091E+03 7280 0.155E+02
DAV: 3 -0.185227503579E+02 -0.48398E+00 -0.48396E+00 10160 0.132E+01
DAV: 4 -0.185240376489E+02 -0.12873E-02 -0.12873E-02 8200 0.464E-01
DAV: 5 -0.185240414620E+02 -0.38131E-05 -0.38134E-05 10480 0.143E-02 0.154E+01
DAV: 6 -0.164394189736E+02 0.20846E+01 -0.72731E+01 6800 0.881E+01 0.783E+00
DAV: 7 -0.149577606472E+02 0.14817E+01 -0.11165E+01 6720 0.344E+01 0.111E+00
DAV: 8 -0.149148685722E+02 0.42892E-01 -0.18431E-02 7280 0.139E+00 0.620E-01
DAV: 9 -0.149081819010E+02 0.66867E-02 -0.51836E-02 8440 0.220E+00 0.636E-02
DAV: 10 -0.149081753161E+02 0.65849E-05 -0.22202E-04 9000 0.149E-01 0.115E-02
DAV: 11 -0.149081715196E+02 0.37964E-05 -0.16862E-06 8520 0.968E-03 0.844E-04
DAV: 12 -0.149081716052E+02 -0.85577E-07 -0.88627E-08 4640 0.248E-03
1 F= -.14908172E+02 E0= -.14905377E+02 d E =-.149082E+02
N E dE d eps ncg rms rms(c)
DAV: 1 -0.149363274123E+02 -0.28156E-01 -0.51946E-01 6720 0.830E+00 0.131E+00
DAV: 2 -0.149199449776E+02 0.16382E-01 -0.41291E-01 6720 0.646E+00 0.538E-01
DAV: 3 -0.149119014518E+02 0.80435E-02 -0.68998E-02 6720 0.271E+00 0.192E-01
DAV: 4 -0.149112107699E+02 0.69068E-03 -0.11896E-03 9800 0.216E-01 0.117E-02
DAV: 5 -0.149112123403E+02 -0.15704E-05 -0.80430E-06 7320 0.267E-02 0.594E-03
DAV: 6 -0.149112142891E+02 -0.19488E-05 -0.13229E-06 4760 0.997E-03 0.128E-04
DAV: 7 -0.149112142406E+02 0.48522E-07 -0.22610E-08 3560 0.101E-03
2 F= -.14911214E+02 E0= -.14908447E+02 d E =-.304264E-02
N E dE d eps ncg rms rms(c)
DAV: 1 -0.150389301673E+02 -0.12772E+00 -0.26022E+00 6720 0.187E+01 0.290E+00
DAV: 2 -0.149566578541E+02 0.82272E-01 -0.20265E+00 6720 0.143E+01 0.119E+00
DAV: 3 -0.149175257623E+02 0.39132E-01 -0.34317E-01 6720 0.604E+00 0.437E-01
DAV: 4 -0.149140056079E+02 0.35202E-02 -0.72028E-03 10080 0.565E-01 0.378E-02
DAV: 5 -0.149139977907E+02 0.78172E-05 -0.11528E-04 7840 0.108E-01 0.134E-02
DAV: 6 -0.149140005388E+02 -0.27480E-05 -0.66656E-06 7680 0.220E-02 0.371E-04
DAV: 7 -0.149140005655E+02 -0.26792E-07 -0.64306E-08 3960 0.196E-03
3 F= -.14914001E+02 E0= -.14911298E+02 d E =-.582896E-02

这是上一节里的关于bulk计算(结构优化)的部分内容,可以看到,reached required accuracy - stopping structural energy minimisation——离子步已经收敛,在之后的OSZICAR里,电子步迭代最多7步,并未超过设置的最大电子步数目。因此对于bulk计算,==电子步和离子步都收敛==

而对于slab模型的单点计算,我也傻傻的这么做😢

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
[ctan@baifq-hpc141 conventional-slab-static]$ grep accuracy OUTCAR
[ctan@baifq-hpc141 conventional-slab-static]$ cat OSZICAR
N E dE d eps ncg rms rms(c)
DAV: 1 0.178654188479E+04 0.17865E+04 -0.88940E+04 10388 0.264E+03
DAV: 2 0.200572973408E+03 -0.15860E+04 -0.15355E+04 10388 0.519E+02
DAV: 3 -0.535557174826E+02 -0.25413E+03 -0.23428E+03 10438 0.249E+02
DAV: 4 -0.683234131101E+02 -0.14768E+02 -0.14274E+02 11572 0.815E+01
DAV: 5 -0.687193157556E+02 -0.39590E+00 -0.39405E+00 12400 0.117E+01 0.292E+01
RMM: 6 -0.611494529970E+02 0.75699E+01 -0.15131E+02 11054 0.104E+02 0.161E+01
RMM: 7 -0.558685188476E+02 0.52809E+01 -0.32335E+01 10772 0.530E+01 0.226E+00
RMM: 8 -0.559296519539E+02 -0.61133E-01 -0.22789E+00 11552 0.614E+00 0.111E+00
RMM: 9 -0.559488617332E+02 -0.19210E-01 -0.37622E-01 11539 0.361E+00 0.509E-01
RMM: 10 -0.559454147061E+02 0.34470E-02 -0.13294E-01 12844 0.193E+00 0.214E-01
RMM: 11 -0.559437301097E+02 0.16846E-02 -0.22655E-02 12875 0.691E-01 0.208E-01
RMM: 12 -0.559413428432E+02 0.23873E-02 -0.55145E-03 11960 0.354E-01 0.586E-02
RMM: 13 -0.559414638149E+02 -0.12097E-03 -0.17465E-03 13023 0.187E-01 0.323E-02
RMM: 14 -0.559411451635E+02 0.31865E-03 -0.83089E-04 11756 0.134E-01 0.373E-02
RMM: 15 -0.559412390394E+02 -0.93876E-04 -0.47990E-04 11372 0.168E-01 0.284E-02
RMM: 16 -0.559414153722E+02 -0.17633E-03 -0.21175E-04 11475 0.107E-01 0.129E-02
RMM: 17 -0.559414836424E+02 -0.68270E-04 -0.10022E-04 11751 0.526E-02 0.889E-03
RMM: 18 -0.559415856898E+02 -0.10205E-03 -0.51724E-05 11898 0.366E-02 0.682E-03
RMM: 19 -0.559416259149E+02 -0.40225E-04 -0.26297E-05 11504 0.350E-02 0.861E-03
RMM: 20 -0.559416406670E+02 -0.14752E-04 -0.10748E-05 11263 0.208E-02 0.367E-03
RMM: 21 -0.559416498741E+02 -0.92071E-05 -0.47420E-06 8911 0.102E-02 0.119E-03
RMM: 22 -0.559416543655E+02 -0.44914E-05 -0.22030E-06 7155 0.977E-03 0.193E-03
RMM: 23 -0.559416573370E+02 -0.29715E-05 -0.14726E-06 6702 0.894E-03 0.135E-03
RMM: 24 -0.559416575122E+02 -0.17515E-06 -0.44209E-07 6318 0.507E-03
1 F= -.55941658E+02 E0= -.55944345E+02 d E =0.806296E-02

不管怎么修改ALGOENCUT等参数,grep accuracy OUTCAR始终不返回reached required accuracy - stopping structural energy minimisation,后来我才意识到,这是==单点计算==,NSW = 0,并不会进行离子步迭代呀,所以我们只需关注电子步收敛即可,很明显电子步收敛,OUTCAR里的信息也可说明

1
2
[ctan@baifq-hpc141 conventional-slab-static]$ grep abort OUTCAR
------------------------ aborting loop because EDIFF is reached ----------------------------------------

This is definitely a stupid, stupid and stupid mistake😭

就写到这,吃饭去了,待会还要回实验室

使用Primitive Cell

使用bulk计算且收敛的CONTCAR进行slab模型的建立,操作还是类似的

值得注意的是,需要先==find symmetry==,然后才能转换为primitive cell

降低对称性之前
find symmetry

然后impose symmetry,接着再转换为primitive cell,然后建立slab模型

POSCAR会有些不同

1
2
3
4
5
6
7
8
9
10
11
12
CONTCAR\(1\1\1)
1.0
2.5701000690 0.0000000000 0.0000000000
-1.2850500345 2.2257719501 0.0000000000
0.0000000000 0.0000000000 21.2954006195
Cu
4
Direct
0.000000000 0.000000000 0.000000000
0.666670000 0.333330000 0.098540000
0.333330000 0.666670000 0.197080000
0.000000000 0.000000000 0.295620000

计算后的能量也有所不同:

1
2
(base) storm@DESKTOP-HE4FQ8Q:~/my-learn/ex42/Cu-primitive-slab$ grep '  without' OUTCAR
energy without entropy= -13.97087892 energy(sigma->0) = -13.97120814

值得一提的是,在白老师的计算集群上运行时报错了(vasp.6.4.2)

如果截断能ENCUT过小可能报错Error EDDDAV: Call to ZHEGV failed. Returncode = 7 1 8

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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
[ctan@baifq-hpc141 Cu-primmitive-slab]$ ls
INCAR KPOINTS POSCAR POTCAR
[ctan@baifq-hpc141 Cu-primmitive-slab]$ clear
[ctan@baifq-hpc141 Cu-primmitive-slab]$ sbatch vasp
Submitted batch job 4321
[ctan@baifq-hpc141 Cu-primmitive-slab]$ tail -f OUTCAR
--------------------------------------- Ionic step 1 -------------------------------------------




--------------------------------------- Iteration 1( 1) ---------------------------------------


POTLOK: cpu time 0.0215: real time 0.0258
SETDIJ: cpu time 0.0109: real time 0.0118
EDDAV: cpu time 0.7451: real time 0.7613
DOS: cpu time 0.0052: real time 0.0060
--------------------------------------------
LOOP: cpu time 0.7826: real time 0.8050

eigenvalue-minimisations : 2548
total energy-change (2. order) : 0.3445715E+03 (-0.1837296E+04)
number of electron 44.0000000 magnetization
augmentation part 44.0000000 magnetization

DIPCOR: dipole corrections for dipol
direction 3 min pos 155,
dipolmoment 0.000000 0.000000 0.000026 electrons x Angstroem
Tr[quadrupol] -31.061780

energy correction for charged system 0.000000 eV
dipol+quadrupol energy correction -0.000000 eV
added-field ion interaction 0.001451 eV (added to PSCEN)


Free energy of the ion-electron system (eV)
---------------------------------------------------
alpha Z PSCENC = 96.14159556
Ewald energy TEWEN = 16384.30326894
-Hartree energ DENC = -21746.47963692
-exchange EXHF = 0.00000000
-V(xc)+E(xc) XCENC = 182.01436818
PAW double counting = 4912.82663951 -5336.19386850
entropy T*S EENTRO = 0.00034874
eigenvalues EBANDS = 288.07434016
atomic energy EATOM = 5563.88447798
Solvation Ediel_sol = 0.00000000
---------------------------------------------------
free energy TOTEN = 344.57153365 eV

energy without entropy = 344.57118491 energy(sigma->0) = 344.57141740


--------------------------------------------------------------------------------------------------------




--------------------------------------- Iteration 1( 2) ---------------------------------------


-----------------------------------------------------------------------------
| |
| EEEEEEE RRRRRR RRRRRR OOOOOOO RRRRRR ### ### ### |
| E R R R R O O R R ### ### ### |
| E R R R R O O R R ### ### ### |
| EEEEE RRRRRR RRRRRR O O RRRRRR # # # |
| E R R R R O O R R |
| E R R R R O O R R ### ### ### |
| EEEEEEE R R R R OOOOOOO R R ### ### ### |
| |
| Error EDDDAV: Call to ZHEGV failed. Returncode = 7 1 8 |
| |
| ----> I REFUSE TO CONTINUE WITH THIS SICK JOB ... BYE!!! <---- |
| |
-----------------------------------------------------------------------------

小结

  • 此计算是在bulk计算的基础上(最优结构)进行的单点计算,计算时不能进行离子弛豫,不能进行离子步迭代,默认NSW=0
  • 不进行离子弛豫,那么自然不能设置离子步收敛条件,不能设置EDIFFG
  • 如果不报错运行但结果不收敛可试着增大电子步数NSW以及截断能ENCUT

Comment and share

bulk计算,直译为“体相计算”,指的是模拟理想的三维无限大晶体(无表面、界面或缺陷,不模拟则体系太大算不动),利用 周期性边界条件 描述材料的体相性质,主要用于研究材料的本征物理化学性质(如晶格结构、电子能带、热力学稳定性等)

准备POSCAR

打开materials studio-File-Import-Structures-Metals-pure metals-Cu.xsd,然后export为cif文件

用vesta打开Cu.cif文件,再export data,导出保存为POSCAR文件

准备KPOINTS

1
2
3
4
5
KPOINTS for Cu-concentional-cell-bulk-calc
0
Gamma
13 13 13
0 0 0

准备INCAR

1
2
3
4
5
6
7
8
9
10
11
12
SYSTEM = Cu-conventional-cell-bulk-calculation
ISMEAR = 0
SIGMA = 0.1
ENCUT = 700
ISIF = 3
NSW = 200
NELM = 200
IBRION = 2
POTIM = 0.1
EDIFFG = -0.01
ALGO = Normal
EDIFF = 1E-6

提交任务

第一次使用超算平台,还算顺利,稍微折腾了一下终于能成功提交任务了

1
2
3
[ctan@baifq-hpc141 Cu-statistic]$ sbatch vasp  #vasp is the script
[ctan@baifq-hpc141 Cu-statistic]$ grep accuracy OUTCAR
reached required accuracy - stopping structural energy minimisation #calculation finished normally

优化后的结构CONTCAR如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Cu
1.00000000000000
3.6346781666386048 0.0000000000000000 0.0000000000000000
0.0000000000000000 3.6346781666386048 0.0000000000000000
-0.0000000000000000 0.0000000000000000 3.6346781666386048
Cu/0e71e558f37
4
Direct
0.0000000000000000 0.0000000000000000 -0.0000000000000000
-0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000

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

其他信息:

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
[ctan@baifq-hpc141 bulk-calc]$ cat OSZICAR
N E dE d eps ncg rms rms(c)
DAV: 1 0.103476848532E+03 0.10348E+03 -0.27847E+04 6720 0.246E+03
DAV: 2 -0.180387697818E+02 -0.12152E+03 -0.12091E+03 7280 0.155E+02
DAV: 3 -0.185227503579E+02 -0.48398E+00 -0.48396E+00 10160 0.132E+01
DAV: 4 -0.185240376489E+02 -0.12873E-02 -0.12873E-02 8200 0.464E-01
DAV: 5 -0.185240414620E+02 -0.38131E-05 -0.38134E-05 10480 0.143E-02 0.154E+01
DAV: 6 -0.164394189736E+02 0.20846E+01 -0.72731E+01 6800 0.881E+01 0.783E+00
DAV: 7 -0.149577606472E+02 0.14817E+01 -0.11165E+01 6720 0.344E+01 0.111E+00
DAV: 8 -0.149148685722E+02 0.42892E-01 -0.18431E-02 7280 0.139E+00 0.620E-01
DAV: 9 -0.149081819010E+02 0.66867E-02 -0.51836E-02 8440 0.220E+00 0.636E-02
DAV: 10 -0.149081753161E+02 0.65849E-05 -0.22202E-04 9000 0.149E-01 0.115E-02
DAV: 11 -0.149081715196E+02 0.37964E-05 -0.16862E-06 8520 0.968E-03 0.844E-04
DAV: 12 -0.149081716052E+02 -0.85577E-07 -0.88627E-08 4640 0.248E-03
1 F= -.14908172E+02 E0= -.14905377E+02 d E =-.149082E+02
N E dE d eps ncg rms rms(c)
DAV: 1 -0.149363274123E+02 -0.28156E-01 -0.51946E-01 6720 0.830E+00 0.131E+00
DAV: 2 -0.149199449776E+02 0.16382E-01 -0.41291E-01 6720 0.646E+00 0.538E-01
DAV: 3 -0.149119014518E+02 0.80435E-02 -0.68998E-02 6720 0.271E+00 0.192E-01
DAV: 4 -0.149112107699E+02 0.69068E-03 -0.11896E-03 9800 0.216E-01 0.117E-02
DAV: 5 -0.149112123403E+02 -0.15704E-05 -0.80430E-06 7320 0.267E-02 0.594E-03
DAV: 6 -0.149112142891E+02 -0.19488E-05 -0.13229E-06 4760 0.997E-03 0.128E-04
DAV: 7 -0.149112142406E+02 0.48522E-07 -0.22610E-08 3560 0.101E-03
2 F= -.14911214E+02 E0= -.14908447E+02 d E =-.304264E-02
N E dE d eps ncg rms rms(c)
DAV: 1 -0.150389301673E+02 -0.12772E+00 -0.26022E+00 6720 0.187E+01 0.290E+00
DAV: 2 -0.149566578541E+02 0.82272E-01 -0.20265E+00 6720 0.143E+01 0.119E+00
DAV: 3 -0.149175257623E+02 0.39132E-01 -0.34317E-01 6720 0.604E+00 0.437E-01
DAV: 4 -0.149140056079E+02 0.35202E-02 -0.72028E-03 10080 0.565E-01 0.378E-02
DAV: 5 -0.149139977907E+02 0.78172E-05 -0.11528E-04 7840 0.108E-01 0.134E-02
DAV: 6 -0.149140005388E+02 -0.27480E-05 -0.66656E-06 7680 0.220E-02 0.371E-04
DAV: 7 -0.149140005655E+02 -0.26792E-07 -0.64306E-08 3960 0.196E-03
3 F= -.14914001E+02 E0= -.14911298E+02 d E =-.582896E-02
[ctan@baifq-hpc141 bulk-calc]$ grep accuracy OUTCAR
reached required accuracy - stopping structural energy minimisation
[ctan@baifq-hpc141 bulk-calc]$ grep ' without' 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

Comment and share

乙醇分子的结构优化

INCAR

1
2
3
4
5
6
7
SYSTEM = ethanol
ISMEAR = 0 #Gaussian ismear, for insulator
SIGMA = 0.01
NSW = 200

IBRION = 2 #Optimizarion
POTIM = 0.01 #a little conservative step

KPOINTS

1
2
3
4
5
K-POINTS
0
Gamma
1 1 1
0 0 0

POSCAR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Title
1.0
4.8082499504 0.0000000000 0.0000000000
0.0000000000 4.8082499504 0.0000000000
0.0000000000 0.0000000000 4.8082499504
O C H
1 2 6
Direct
0.747990456 0.438959082 0.500000000
0.508714210 0.577741385 0.500000000
0.269437963 0.438959082 0.500000000
0.900000025 0.568403268 0.500000000
0.508714210 0.689008479 0.307289552
0.508714210 0.689008479 0.692710448
0.100000016 0.583231948 0.500000000
0.257895302 0.310991525 0.318291476
0.257895302 0.310991525 0.681708524

分子结构可以使用ChemDraw绘制,也可以去RSC的网站上下载https://www.chemspider.com/、

但是不管是下载的mol结构还是ChemDraw导出的mol结构,导入进Vesta里面发现都少了氢原子,笨比博主(😭)至今没找到解决方案,只能笨拙地把mol结构导入到GaussView里面,GaussView会自动补加上氢原子,然后再导入到Vesta中,再导出为POSCAR文件,当然理论上可以手动改POSCAR,不一定要使用可视化软件

POTCAR

使用基础的不带后缀的PBE赝势

1
cat ..potpaw_PBE.64/O/POTCAR ..potpaw_PBE.64/C/POTCAR ..potpaw_PBE.64/H/POTCAR >> POTCAR

对比优化前后

简单对比一下键长吧

  • 优化前
    • C-C bond: 1.33001(0) Å
    • O-H bond: 0.96000(0) Å
O-H bond
C-C bond
  • 优化后
    • C-C bond: 1.53893(0) Å
    • O-H bond: 0.96454(0) Å
O-H bond
C-C bond

频率计算

频率计算的作用:

  • 确定结构是否稳定
  • 看振动方式和大小,用来和实验对比
  • 反应热,反应能垒,吸附能等的零点能矫正
  • 确认过渡态(有一个振动的虚频)
  • 热力学中计算entropy,用于计算化学势,微观动力学中的指前因子和反应能垒。

POSCAR

使用先前优化的CONTCAR作为POSCAR

INCAR

修改一下IBRION,计算频率(IBRION = 5)而非结构优化(IBRION = 2)

1
2
3
4
5
6
7
8
SYSTEM = ethanol
ISMEAR = 0 #Gaussian ismear, for insulator
SIGMA = 0.01
NSW = 1
IBRION = 5 #freq calculation
POTIM = 0.01 #a little conservative step
EDIFF = 1E-8
NFREE = 2

Jmol分析频率

1
jmol OUTCAR

Tools-AtomSetChooser-Frequencies

Jmol分析频率

Comment and share

单点能计算又叫做静态计算,计算时不会优化结构

INCAR

1
2
3
4
5
6
7
8
SYSTEM = O
ISMEAR = 0 #绝缘体系ISMEAR = 0
SIGMA = 0.01
ISPIN = 2 #氧气是顺磁性分子,ISPIN = 2开启自旋
MAGMON = 2*2 #设置磁矩

#IBRION = 2 确定计算过程中晶体结构如何变化的参数,不设置使默认IBRION = -1,IBRION =2为共轭梯度算法
#NSW = 10 离子步10步

POSCAR

1
2
3
4
5
6
7
8
9
10
O
1.0
7.5 0.0 0.0
0.0 8.0 0.0
0.0 0.0 8.9 #三斜的格子打破对称性,symmetry broken solution!!!
O
2
Cartesian
0.0 0.0 0.0
0.0 0.0 1.2074

KPOINTS

1
2
3
4
5
K-POINTS
0
Gamma
1 1 1
0 0 0

运行mpirun -n 8 vasp,正常结束,OUTCAR的部分信息如下,显然,单点计算的结果并不可靠

α电子部分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
spin component 1

k-point 1 : 0.0000 0.0000 0.0000
band No. band energies occupation
1 -32.9044 1.00000 #1σ轨道
2 -20.2474 1.00000 #2σ轨道
3 -13.3174 1.00000 #3σ轨道,能量应比1π和2π轨道更低,不合理
4 -13.3171 1.00000 #1π轨道
5 -13.2940 1.00000 #2π轨道,能量应与1π轨道相同,不合理
6 -6.5560 1.00000 #3π轨道
7 -6.5524 1.00000 #4π轨道
8 -0.3909 0.00000
9 1.1220 0.00000
10 1.5493 0.00000
11 1.7604 0.00000
12 2.0536 0.00000
13 2.0592 0.00000
14 2.4672 0.00000
15 3.1984 0.00000
16 3.6244 0.00000
Fermi energy: -6.4848556277

β电子,正常

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
spin component 2

k-point 1 : 0.0000 0.0000 0.0000
band No. band energies occupation
1 -31.6977 1.00000
2 -18.4541 1.00000
3 -12.3485 1.00000
4 -11.4763 1.00000
5 -11.4760 1.00000
6 -4.2772 0.00000
7 -4.2742 0.00000
8 -0.2440 0.00000
9 1.4049 0.00000
10 1.6215 0.00000
11 1.9473 0.00000
12 2.1710 0.00000
13 2.2499 0.00000
14 2.5860 0.00000
15 3.5816 0.00000
16 3.8831 0.00000
Molecule orbital of oxygen

对氧气进行结构优化,更改INCAR

1
2
3
4
5
6
7
8
SYSTEM = O
ISMEAR = 0 #绝缘体系ISMEAR = 0
SIGMA = 0.01
ISPIN = 2 #氧气是顺磁性分子,ISPIN = 2开启自旋
MAGMON = 2*2 #设置磁矩

IBRION = 2 #确定计算过程中晶体结构如何变化的参数,不设置使默认IBRION = -1,IBRION =2为共轭梯度算法
NSW = 10 #离子步10步

计算后发现α电子部分已经修正

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
spin component 1

k-point 1 : 0.0000 0.0000 0.0000
band No. band energies occupation
1 -32.2878 1.00000 #1σ轨道
2 -20.4859 1.00000 #2σ轨道
3 -13.2238 1.00000 #3σ轨道,能量合理
4 -13.0591 1.00000
5 -13.0588 1.00000
6 -6.7814 1.00000
7 -6.7778 1.00000
8 -0.3857 0.00000
9 1.0213 0.00000
10 1.6177 0.00000
11 1.7818 0.00000
12 2.0634 0.00000
13 2.0746 0.00000
14 2.4643 0.00000
15 2.8035 0.00000
16 3.6409 0.00000
Fermi energy: -6.7122607632

用VESTA打开优化好的结构(CONTCAR)可以查看键长

优化完的氧气分子的键长

Comment and share

2025年4月AMD平台下VASP的编译安装

WSL2安装

本次安装选择的Ubuntu作为发行版,使用人数较多,解决方案较多

Ubuntu 24.04.1 on AMD 7840HS

更新源

1
sudo apt update

安装编译器

AOCC

下载aocc-compiler-5.0.0.tar,解压

1
2
3
4
tar -xvf aocc-compiler-5.0.0.tar #解压
cd aocc-compiler-3.2.0 #进入目录
install.sh #如果没有可执行权限请赋予可执行权限, chmod +x install.sh
source .setenv_AOCC.sh #这一行可以写入~/.bashrc,然后source ~/.bashrc,刷新一下环境变量

AOCL

下载AOCL 5.0 binary packages compiled with AOCC 5.0

1
2
3
tar -zxvf XXX.tar.gz #解压tar.gz包,如果报错请查询Linux下的解压命令,up记不太清了、、、  
cd XXX
install.sh -t /home/XXX #安装,可以指定安装目录,up安装在/home目录下的某个文件夹里的,这样不需要管理员权限、

安装OpenMPI

在安装OpenMPI之前请确保AOCC、AOCL以及必备的依赖已安装完全

1
2
3
4
5
6
7
8
which clang #检查clang
which clang++ #检查clang++
which flang #检查flang
以上三个来自于AOCC和AOCL
which c #检查有没有C
which c++ #检查有没有C++
sudo apt install g++
sudo apt install gcc

下载OpenMPI的稳定版本,然后解压,进入目录

1
2
3
configure CC=clang CXX=clang++ FC=flang --prefix=/xxxx #可以手动指定/xxx目录
make -j4 #以4核心编译,大小可调,比如16核处理器就改成make -j16
make install #编译安装

OpenMPI安装时报错很可能是由于依赖问题造成的,建议安装之前sudo apt install g++ cmake gcc

~/.bashrc中添加环境变量

1
2
3
4
MPI_HOME=/xxx/; make #这一步容易出错,记得添加对目录
export PATH=${MPI_HOME}/bin:$PATH
export LD_LIBRARY_PATH=${MPI_HOME}/lib:$LD_LIBRARY_PATH
export MANPATH=${MPI_HOME}/share/man:$MANPATH

安装完后检查一下which mpirun

编译VASP

不提供源码包,可以尝试去各种论坛,微信公众号搜索,仅供学习使用

本次使用的源码包是vasp.6.4.2.tgz

编译之前首先要选择模板,也就是makefile.include,将/vasp.x.x.x/arch文件夹下的makefile.include.aocc_ompi_aocl复制一份到/vasp.x.x.x根目录,因为本次平台是AMD ZEN4,安装了AOCC,AOCL和OpenMPI,所以选择此模板,接下来修改模板

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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
# Default precompiler options
CPP_OPTIONS = -DHOST=\"LinuxGNU\" \
-DMPI -DMPI_BLOCK=8000 -Duse_collective \
-DscaLAPACK \
-DCACHE_SIZE=4000 \
-Davoidalloc \
-Dvasp6 \
-Duse_bse_te \
-Dtbdyn \
-Dfock_dblbuf

CPP = flang -E -ffree-form -C -w $*$(FUFFIX) >$*$(SUFFIX) $(CPP_OPTIONS)

FC = mpif90
FCL = mpif90

FREE = -ffree-form -ffree-line-length-none

FFLAGS = -w -fno-fortran-main -Mbackslash

OFLAG = -O2
OFLAG_IN = $(OFLAG)
DEBUG = -O0

OBJECTS = fftmpiw.o fftmpi_map.o fftw3d.o fft3dlib.o
OBJECTS_O1 += fftw3d.o fftmpi.o fftmpiw.o
OBJECTS_O2 += fft3dlib.o

# For what used to be vasp.5.lib
CPP_LIB = $(CPP)
FC_LIB = $(FC)
CC_LIB = clang
CFLAGS_LIB = -O
FFLAGS_LIB = -O1
FREE_LIB = $(FREE)

OBJECTS_LIB = linpack_double.o

# For the parser library
CXX_PARS = clang++
LLIBS = -lstdc++

##
## Customize as of this point! Of course you may change the preceding
## part of this file as well if you like, but it should rarely be
## necessary ...
##

# When compiling on the target machine itself, change this to the
# relevant target when cross-compiling for another architecture
VASP_TARGET_CPU ?= -march=znver4 #7840HS belongs to ZEN4 platform
FFLAGS += $(VASP_TARGET_CPU)

# BLAS (mandatory)
AMDBLIS_ROOT ?= /home/storm/AOCL_INSTALL/5.0.0/aocc #Installation category of AOCL
BLAS = -L${AMDBLIS_ROOT}/lib -lblis

# LAPACK (mandatory)
AMDLIBFLAME_ROOT ?= /home/storm/AOCL_INSTALL/5.0.0/aocc #Installation category of AOCL
LAPACK = -L${AMDLIBFLAME_ROOT}/lib -lflame

# scaLAPACK (mandatory)
AMDSCALAPACK_ROOT ?= /home/storm/AOCL_INSTALL/5.0.0/aocc #Installation category of AOCL
SCALAPACK = -L${AMDSCALAPACK_ROOT}/lib -lscalapack

LLIBS += $(SCALAPACK) $(LAPACK) $(BLAS)

# FFTW (mandatory)
AMDFFTW_ROOT ?= /home/storm/AOCL_INSTALL/5.0.0/aocc #Installation category of AOCL
LLIBS += -L$(AMDFFTW_ROOT)/lib -lfftw3
INCS += -I$(AMDFFTW_ROOT)/include

# HDF5-support (optional but strongly recommended)
#CPP_OPTIONS+= -DVASP_HDF5
#HDF5_ROOT ?= /path/to/your/hdf5/installation
#LLIBS += -L$(HDF5_ROOT)/lib -lhdf5_fortran
#INCS += -I$(HDF5_ROOT)/include

# For the VASP-2-Wannier90 interface (optional)
#CPP_OPTIONS += -DVASP2WANNIER90
#WANNIER90_ROOT ?= /path/to/your/wannier90/installation
#LLIBS += -L$(WANNIER90_ROOT)/lib -lwannier

然后运行make all即可,官网给的命令是并行安装make DEPS=1 -Jn <target>,我自己尝试并行安装后测试make test会报错,所以make all,这样去测试全部通过不报错

测试安装

为了测试是否正确安装,可以进行测试

1
make test

添加到环境变量

官方提供的是源码包,通过编译,我们得到了CPU可执行的二进制包,再二进制包的目录下(bin/目录下),我可以同在终端中直接输入二进制包的名字来直接执行(VASP编译会产生三个二进制包:vasp_std, vasp_gam, vasp_ncl),比如在bin/目录下直接输入vasp_std就可以看到VASP的标志,但是为了随时随地在任意目录下也能执行计算,因此需要将bin/目录添加到环境变量里(前提是make test不报错)

1
2
export PATH=$PATH:/path/to/vasp.x.x.x/bin
source ~/.bashrc

Comment and share

使用for循环批量创建文件夹

1
for i in {2..9}; do cp 0.01 0.0$i -r ; done

使用sed命令不打开文件而进行替换

1
sed '3s/0.01/0.02/g' ICNAR

对第三行(3)出现的所有0.01(g,全局替换)进行替换(s表示替换操作),替换为0.02,只输出替换后的结果

sed ‘3/0.01/0.02/g’ INCAR > INCAR最后什么也没有

加上-i参数可以直接进行编辑

结合for循环和sed批量命名

1
2
3
4
5
6
7
8
9
10
11
[storm@cachyos-x8664 ex03]$ for i in *; do sed -i "3s/0.05/$i/g" $i/INCAR ; done
[storm@cachyos-x8664 ex03]$ grep SIGMA */INCAR
0.01/INCAR:SIGMA = 0.01
0.02/INCAR:SIGMA = 0.02
0.03/INCAR:SIGMA = 0.03
0.04/INCAR:SIGMA = 0.04
0.05/INCAR:SIGMA = 0.05
0.06/INCAR:SIGMA = 0.06
0.07/INCAR:SIGMA = 0.07
0.08/INCAR:SIGMA = 0.08
0.09/INCAR:SIGMA = 0.09

使用双引号以读取变量的值

注意这里使用的是英文括号而不是花括号

1
for i in $(seq 8 2 16); do cp 888/POSCAR ${i}${i}${i}/POSCAR; done

“提交任务的命令”

1
yhbatch -p gsc -N 1 -J test job_sub 

找能量的命令

1
2
3
grep  without OUTCAR | tail -n 1
grep ' without' OUTCAR | tail -n 1 # 本人常用的是这个
grep sigma OUTCAR | tail -n 1

提交任务,多个文件夹

1
for i in *; do cd $i ; vasp1; cd $OLDPWD; done

alias vasp1=‘mpirun -n 8 vasp’

输出时间信息

1
for i in *0; do echo -e  $i "\t" $(grep User $i/OUTCAR | awk '{print $4}'); done

绘图脚本

1
2
3
4
5
6
7
8
9
import matplotlib.pyplot as plt
import numpy as np

x,y = np.loadtxt('data.dat', delimiter = ',',
usecols=(0,1), unpack=True)
plt.xlabel('ENCUT / eV')
plt.ylabel('Ttme / S')
plt.plot(x,y, 'rs-', linewidth=2.0)
plt.show()

在vim中进行替换

1
: 10,30s/$/T T T/g

$表示每一行的末尾

提取能量(用制表符TAB进行分隔)

1
2
3
4
5
for i in [0-9]*/; do
dir=${i%/}
energy=$(grep ' without' "$dir/OUTCAR" | tail -n 1 | awk '{print $7}')
printf "%s\t%s\n" "$dir" "$energy"
done > data

Comment and share

John Doe

author.bio


author.job


Changchun, China