Vaspkit的使用

in 从零开始的科研之旅, VASP

安装、配置好vaspkit之后输入vaspkit可启动软件包

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
(base) storm@DESKTOP-HE4FQ8Q:~/my-learn/trash-vaspkit$ vaspkit
            \\\///
           / _  _ \         Hey, you must know what you are doing.
         (| (o)(o) |)       Otherwise you might get wrong results.
 o-----.OOOo--()--oOOO.------------------------------------------o
 |         VASPKIT Standard Edition 1.5.1 (27 Jan. 2024)         |
 |         Lead Developer: Vei WANG (wangvei@icloud.com)         |
 |      Main Contributors: Gang TANG, Nan XU & Jin-Cheng LIU     |
 |  Online Tutorials Available on Website: https://vaspkit.com   |
 o-----.oooO-----------------------------------------------------o
        (   )   Oooo.                          VASPKIT Made Simple
         \ (    (   )
          \_)    ) /
                (_/
 ===================== Structural Utilities ======================
 01) VASP Input-Files Generator    02) Mechanical Properties
 03) K-Path for Band-Structure     04) Structure Editor
 05) Catalysis-ElectroChem Kit     06) Symmetry Analysis
 07) Materials Databases           08) Advanced Structure Models
 ===================== Electronic Utilities ======================
 11) Density-of-States             21) Band-Structure
 23) 3D Band-Structure             25) Hybrid-DFT Band-Structure
 26) Fermi-Surface                 28) Band-Structure Unfolding
 31) Charge-Density Analysis       42) Potential Analysis
 44) Piezoelectric Properties      51) Wave-Function Analysis
 62) Magnetic Analysis             65) Spin-Texture
 68) Transport Properties
 ======================== Misc Utilities =========================
 71) Optical Properties            72) Molecular-Dynamics Kit
 74) User Interface                78) VASP2other Interface
 84) ABACUS Interface              91) Semiconductor Kit
 92) 2D-Material Kit               95) Phonon Analysis
 0)  Quit
 ------------>>

以下是中文翻译对照：

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
===================== 结构工具 ======================
01) VASP 输入文件生成器
02) 力学性质
03) 能带结构的 K 路径
04) 结构编辑器
05) 催化 - 电化学套件
06) 对称性分析
07) 材料数据库
08) 高级结构模型
==================== 电子工具 ======================
11) 密度态
21) 能带结构
23) 三维能带结构
25) 混合 DFT 能带结构
26) 费米面
28) 能带结构展开
31) 电荷密度分析
42) 势分析
44) 压电性质
51) 波函数分析
62) 磁性分析
65) 自旋纹理
68) 输运性质
====================== 杂项工具 ========================
71) 光学性质
72) 分子动力学套件
74) 用户界面
78) VASP 转其他格式接口
84) ABACUS 接口
91) 半导体套件
92) 二维材料套件
95) 声子分析
0) 退出

我目前还是初学者小白，使用得最多的功能还是生成输入文件功能尤其是生成INCAR文件，

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
 ------------>>
01
 ==================== VASP Input Files Options ===================
 101) Customize INCAR File
 102) Generate KPOINTS File for SCF Calculation
 103) Generate POTCAR File with Default Setting
 104) Generate POTCAR File with User Specified Potential
 105) Generate POSCAR File from cif (no fractional occupations)
 106) Generate POSCAR File from Material Studio xsd (retain fixes)
 107) Reformat POSCAR File in Specified Order of Elements
 108) Successive Procedure to Generate VASP Files and Check
 109) Submit Job Queue

 0)   Quit
 9)   Back
 ------------>>
  ------------>>
101
 +---------------------------- Tip ------------------------------+
 |          WARNNING: You MUST know what wou are doing!          |
 |Some Parameters in INCAR file need to be set/adjusted manually.|
 +---------------------------------------------------------------+
 ======================== INCAR Options ==========================
 ST) Static-Calculation            SR) Standard Relaxation
 MG) Magnetic Properties           SO) Spin-Orbit Coupling
 D3) DFT-D3 no-damping Correction  H6) HSE06 Calculation
 PU) DFT+U Calculation             MD) Molecular Dynamics
 GW) GW0 Calculation               BS) BSE Calculation
 DC) Elastic Constant              EL) ELF Calculation
 BD) Bader Charge Analysis         OP) Optical Properties
 EC) Static Dielectric Constant    PC) Decomposed Charge Density
 PH) Phonon-Calculation            PY) Phonon with Phononpy
 NE) Nudged Elastic Band (NEB)     DM) The Dimer Method
 FQ) Frequence Calculation         LR) Lattice Relaxation
 MT) Meta-GGA Calculation          PZ) Piezoelectric Calculation

 0)   Quit
 9)   Back
 ------------>>
 Input Key-Parameters (STH6D3 means HSE06-D3 Static-Calcualtion)

以下是中文对照：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
+---------------------------- 提示 ------------------------------+
|          警告：你必须清楚自己在做什么！          |
|INCAR 文件中的一些参数需要手动设置/调整。|
+---------------------------------------------------------------+
======================= INCAR 选项 ==========================
ST) 静态计算            SR) 标准松弛
MG) 磁性性质            SO) 自旋轨道耦合
D3) DFT-D3 无阻尼校正  H6) HSE06 计算
PU) DFT+U 计算          MD) 分子动力学
GW) GW0 计算            BS) BSE 计算
DC) 弹性常数            EL) ELF 计算
BD) Bader 电荷分析      OP) 光学性质
EC) 静态介电常数        PC) 分解电荷密度
PH) 声子计算            PY) 使用 Phonopy 的声子计算
NE) 强制弹性带（NEB）   DM) 二聚体方法
FQ) 频率计算            LR) 晶格松弛
MT) Meta-GGA 计算       PZ) 压电计算

LR晶格弛豫就对应着结构优化，当然生成的INCAR文件可能仍需要进行一些修改，比如ENCUT可能需要调整为赝势中ENMAX的1.3倍

生成KPOINTS文件

对于非能带的计算，选择102) Generate KPOINTS File for SCF Calculation使用程序自动撒点即可，但是需要用户选择撒点方式和撒点密度，选择Gamma Scheme，使用0.04精度

关于生成 KPOINTS 文件

  • 在 VASP 计算中，KPOINTS 文件用于指定 k 点网格的设置。对于非能带计算，比如进行能量计算、优化结构等，通常不需要像计算能带结构那样沿特定的高对称路径采样 k 点，而是可以在整个布里渊区均匀地分布 k 点，这就是所说的 “自动撒点”。

撒点方式

  • Monkhorst-Pack 方式 ：这是最常用的自动撒点方式。它通过指定一组三个整数（nx, ny, nz），在三个倒易晶格方向上均匀地分布 k 点，形成的网格可以系统地逼近布里渊区的积分。例如，对于简单晶格，若指定 k 点网格为 3×3×3，程序就会按照 Monkhorst-Pack 方式在这三个方向各取 3 个点，生成一个均匀分布的 k 点网格。
  • Gamma 方式 ：这种方式会在 k 点网格中包含 Γ 点（布里渊区的中心点）。对于金属等具有高对称性的材料，在 Γ 点附近采样很重要，这时候 Gamma 方式撒点更合适。比如指定 k 点网格为 4×4×4 并采用 Gamma 方式撒点，生成的网格会以 Γ 点为中心，向周围均匀分布其他 k 点。

K 点密度

  • K 点密度指的是在布里渊区单位体积内分布的 k 点数量，它决定了 k 点网格的疏密程度。较高的 K 点密度意味着更精细地采样布里渊区，能够更精确地计算系统的物理性质，但也会增加计算量和资源消耗。
  • 例如，对于一个计算资源充足的大型超算任务，为了获得高精度结果，可以选择较大的 K 点密度，如 8×8×8 甚至更高；而对于初步探索或计算资源有限的情况，可以适当降低 K 点密度，如 4×4×4。

生成KPOINTS的同时，POTCAR也会被自动生成，前提是在vaspkit的配置文件中正确设置了赝势的路径和赝势稳健的种类（默认type为PBE）

Comment and share

截断能ENCUT的选取

in 从零开始的科研之旅, VASP

ENCUT是用来指定平面波基组能量截止值的参数

以上是wiki对ENCUT的描述，POTCAR文件中已经包含了元素的ENMAX和ENMIN，这也是互联网上对ENCUT的选择方案对通常是1.3倍的ENMAX

侯柱峰老师提供了另外一种方法并指出通过该方法选择的ENCUT通常能满足1.3倍ENMAX，供参考使用

截断能指定了用于波函数展开的平面波基组的截断能量，此能量越大则用来描述波函数的平面波基组越多，精度越高，但计算也越耗时。

可以通过计算测试来选择合适的ENCUT。

以下是用于测试的脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#!/bin/sh
rm WAVECAR
for i in 150 200 250 300 350 400
do
cat > INCAR <<!
SYSTEM = Si-Diamond
ENCUT = $i
ISTART = 0 ; ICHARG = 2
ISMEAR = -5
PREC = Accurate
!
echo "ENCUT = $i eV" ; time vasp
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f \n", $5 }'`
echo $i $E >>comment
done

计算后会得到comment文件，它列出了在不同ENCUT下计算得到的总能量

1
2
3
4
5
6
150 -11.900655
200 -11.938864
250 -11.944599
300 -11.945248
350 -11.945503
400 -11.945622

==总能量变化在0.001 eV即可==，因此在这个例子中ENCUT可以选择250

Comment and share

表面能的计算

in 从零开始的科研之旅, VASP

关于表面能

根据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

POTCAR赝势的选取

in 从零开始的科研之旅, VASP

本文为一个查阅性的记录文章，也可以说是搬运，内容取自侯柱锋老师的使用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

以p(1x1)Cu(111)为例分析表面弛豫的结果

in 从零开始的科研之旅, VASP

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

在表面弛豫前（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之后就可以运行脚本了，将POSCAR和CONTCAR的坐标从分数坐标转换为笛卡尔坐标

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

p(1x1)Cu(111)的优化计算

in 从零开始的科研之旅, VASP

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 = 2 的POTIM更大，但所用时间更长。二者最终的能量相差不大

Comment and share

Cu(111) slab模型的单点计算

in 从零开始的科研之旅, VASP

首先对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

不管怎么修改ALGO、ENCUT等参数，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

Cu单胞的bulk计算

in 从零开始的科研之旅, VASP

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

乙醇分子的气相优化和振动频率计算

in 从零开始的科研之旅, VASP

乙醇分子的结构优化

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

氧气分子单点能的计算及结构优化

in 从零开始的科研之旅, VASP

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

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

John Doe

author.bio


author.job


Changchun, China