1.5 提高效率与准确性

截断能

截断能设置在INCAR文件的ENCUT参数。截断能(energy cutoff)指定了展开波函数用到的平面波基组(planewave basis set)所对应的截断能量。
参考: https://zhuanlan.zhihu.com/p/348826693
截断能越大,代表用来描述波函数的平面波基组越多,计算精度越高,同时计算耗时也越长。
如果没有指定ENCUT,程序将默认ENCUT=ENMAX(注意:当PREC = High时, ENCUT = 1.3*ENMAX)。扩胞并不会影响ENCUT的设置。

  • 测试标准
    ENMAX在保证定性正确的基础上大幅度降低计算耗时,如果想要定量计算某些性质就需要增大截断能了。所以ENCUT测试通常从ENMAX开始测试,一般截止到 1.5*ENMAX 即可(注意:比ENMAX小的测试可能会得到错误的能量)。当体系中每个原子的能量差收敛至0.001 eV/atom时,此ENCUT即为最终取值。

简单测试:
for i in {1..8}; do cp 400 $((400+$i*50)); sed -i "s/400/$((400+$i*50))/g" $((400+$i*50))/INCAR ; done

encut.sh脚本: 参考 https://zhuanlan.zhihu.com/p/115235135

{.line-numbers}
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
#!/bin/bash

#生成KPOINTS文件,KPOINTS不能与INCAR中的KSPACING参数共用,二选一
cat > KPOINTS <<!
A
0
M
9 9 9
0 0 0
!

#产生计算所需POSCAR
cat > POSCAR<<!
Si8
1.0000000000
5.4687280000 0.0000000000 0.0000000000
0.0000000000 5.4687280000 0.0000000000
0.0000000000 0.0000000000 5.4687280000
Si
8
Direct
0.0000000000 0.0000000000 0.0000000000
0.2500000000 0.7500000000 0.7500000000
0.5000000000 0.0000000000 0.5000000000
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.5000000000 0.0000000000
0.7500000000 0.2500000000 0.7500000000
0.7500000000 0.7500000000 0.2500000000
0.2500000000 0.2500000000 0.2500000000
!

for i in $(seq 400 50 700) #截断能从450-800,步数为50
do

#生成vasp静态计算的INCAR
cat > INCAR_static <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = .FALSE. (Projection operators: automatic)
ENCUT = $i (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)

Electronic Relaxation
ISMEAR = -5 (Gaussian smearing; metals:1)
#SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 60 (Max electronic SCF steps)
NELMIN = 4 (Min electronic SCF steps)
EDIFF = 1E-06 (SCF energy convergence; in eV)
GGA = PE (PBEsol exchange-correlation)

Ionic Relaxation
ISIF = 2 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V,4-Shape/Ions)
EDIFFG = -0.001 (Ionic convergence; eV/AA)
#KSPACING = 0.10
!

cp INCAR_static INCAR
echo "ENCUT = $i eV" ; time mpirun -np 16 vasp_std #vasp并行运行命令,根据系统自行修改

rm INCAR_static
rm INCAR
#提取计算得到的能量
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $E >>encut_energy.out
done

K点测试

是倒空间(动量空间)的基本构成点,只取在一个倒空间晶格向量的范围内来描述 k 。总能量计算是对布里渊区内的波函数进行积分,而积分是通过对部分特殊K点的求和完成的。因此,计算一个体系时,需要确定所用KPOINTS的大小。
k.sh如下:

{.line-numbers}
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
#!/bin/bash
#生成vasp静态计算的INCAR
cat > INCAR_static <<!
Global Parameters
ISTART = 0 (Read existing wavefunction; if there)
ICHARG = 2 (Non-self-consistent: GGA/LDA band structures)
LREAL = .FALSE. (Projection operators: automatic)
ENCUT = 500 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level)
LWAVE = .FALSE. (Write WAVECAR or not)
LCHARG = .FALSE. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid; helps GGA convergence)

Electronic Relaxation
ISMEAR = -5 (Gaussian smearing; metals:1)
#SIGMA = 0.05 (Smearing value in eV; metals:0.2)
NELM = 60 (Max electronic SCF steps)
NELMIN = 4 (Min electronic SCF steps)
EDIFF = 1E-06 (SCF energy convergence; in eV)
GGA = PE (PBEsol exchange-correlation)

Ionic Relaxation
ISIF = 2 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V,4-Shape/Ions)
EDIFFG = -0.001 (Ionic convergence; eV/AA)
#KSPACING = 0.10
!

cp INCAR_static INCAR

#产生计算所需POSCAR
cat > POSCAR<<!
Si8
1.0000000000
5.4687280000 0.0000000000 0.0000000000
0.0000000000 5.4687280000 0.0000000000
0.0000000000 0.0000000000 5.4687280000
Si
8
Direct
0.0000000000 0.0000000000 0.0000000000
0.2500000000 0.7500000000 0.7500000000
0.5000000000 0.0000000000 0.5000000000
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.5000000000 0.0000000000
0.7500000000 0.2500000000 0.7500000000
0.7500000000 0.7500000000 0.2500000000
0.2500000000 0.2500000000 0.2500000000
!


for i in $(seq 6 3 15) #K点从6-24,步数为3
do

#生成KPOINTS文件,K点类型:Monkhorst-Pack (M) 或者 Gamma (G)
cat > KPOINTS <<!
Automatic mesh
0
Monkhorst-Pack
$i $i $i
0 0 0
!

echo "KPOINTS = $i" ; time mpirun -n 16 vasp_std #vasp并行运行命令,根据系统自行修改

rm KPOINTS
#提取计算得到的能量
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $E >>kpoints_energy.out
done

显示收敛情况的脚本

参考:https://blog.shishiruqi.com/2019/08/19/paili-energyforce/

  • 简单版:

grep F= OSZICAR |awk '{print $1,$5}'
grep F= OSZICAR |awk '{print $1,$5}' > conv.dat

  • 小脚本:
    {.line-numbers}
    1
    2
    3
    4
    5
    6
    #!/usr/bin/env bash
    for i in *; do
    if [ -e $i/OUTCAR ]; then
    echo -e $i "\t" $(grep ' without' $i/OUTCAR |tail -n 1 | awk '{print $7}')
    fi
    done
  • 复杂:
    ediff.sh如下:(安装gnuplot)
    {.line-numbers}
    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
    #!/bin/sh
    #lipai@mail.ustc.edu.cn

    begin=$1
    if [ $begin =='' ]; then begin=0; fi
    awk -v begin=$begin '/E0/{if ( i<begin ) i++;else print $0 }' OSZICAR >temp.e

    awk '/POSITION/,/drift/{
    if(NF==6) print $4,$5,$6;
    else if($1=="total") print $1 }' OUTCAR >temp.f

    awk '{if($4=="F"||$4=="T") print $4,$5,$6}' CONTCAR >temp.fix
    flag=`wc temp.fix|awk '{print $1}'`
    steps=`grep E0 OSZICAR |tail -1 |awk '{print $1}'`
    if [ flag != '0' ] ; then
    if [ -f temp.fixx ] ; then rm temp.fixx ; fi
    for i in `seq $steps`;do
    cat temp.fix >>temp.fixx
    echo >>temp.fixx
    done
    paste temp.f temp.fixx >temp.ff
    fi

    awk '{ if($1=="total") {print ++i,a;a=0}
    else {
    if($4=="F") x=0; else x=$1;
    if($5=="F") y=0; else y=$2;
    if($6=="F") z=0; else z=$3;
    force=sqrt(x^2+y^2+z^2);
    if(a<force) a=force} }' temp.ff >force.conv

    gnuplot <<EOF
    set term dumb
    set title 'Energy of each ionic step'
    set xlabel 'Ionic steps'
    set ylabel 'Energy(eV)'
    plot 'temp.e' u 1:5 w l t "Energy in eV"
    set title 'Max Force of each ionic step'
    set xlabel 'Ionic steps'
    set ylabel 'Force (eV/Angstrom)'
    plot 'force.conv' w l t "Force in eV/Angstrom"
    EOF

    if [ $steps -gt 8 ] ; then
    tail -5 force.conv >temp.fff
    tail -5 temp.e >temp.ee
    gnuplot <<EOF
    set term dumb
    set title 'Energy of each ionic step for the last few steps'
    set xlabel 'Ionic steps'
    set ylabel 'Energy(eV)'
    plot 'temp.ee' u 1:5 w l t "Energy in eV"
    set title 'Max Force of each ionic step for the last few steps'
    set xlabel 'Ionic steps'
    set ylabel 'Force (eV/Angstrom)'
    plot 'temp.fff' w l t "Force in eV/Angstrom"
    EOF
    rm temp.fff temp.ee
    fi

    rm temp.e temp.f temp.ff temp.fix temp.fixx

批量提交任务

适用于测试、有限位移法求力常数矩阵等情况。
简单: for i in *; do cd $i ; qsub sub4; cd $OLDPWD; done

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/sh
root_path='pwd'

#在每个job文件夹里进行静态计算。
for i in {001..1316} #按顺序从job-001116做静态vasp计算,根据生成的文件夹设定
do
cd job-$i
echo "job-$i", time
sbatch sub.sh
wait
cd ..
done

批量产生文件夹

i in {2..9}; do cp 0.01 0.0$i -r ; done ``` #批量复制0.01文件夹复制成0.02-0.09
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
```sed '3s/0.01/0.02/g' INCAR > INCAR_new```  #将INCAR中的第三行的0.01全部替换成0.02并输出位INCAR_new文件。
```sed -i '3s/0.01/0.02/g' INCAR``` # -i表示直接对源文件进行编辑,也就是说编辑之后源文件被新文件替换掉。
```for i in *; do sed -i "3s/0.05/$i/g" $i/INCAR ; done``` #将文件夹中所有的0.01替换成了文件夹序号。

我们主要对KPOINTS的文件的第四行进行批量操作,将1 1 1改成 2 2 2, 3 3 3 等。复制文件夹:
```for i in {1..6}; do cp ../ex03/0.01/ ${i}${i}${i} -r ; done```
显示行号:
```cat 333/KPOINTS -n ```
修改文件;
```for i in {1..6}; do sed -i '4s/1 1 1/$i $i $i/g' ${i}${i}${i}/KPOINTS ; done```

脚本:
```javascript {.line-numbers}
for i in 3RD.POSCAR.*
do
s=$(echo $i| cut -d "." -f 3) #代表.的3个字符
d=job-$s
mkdir $d
cp $i $d/POSCAR
cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d
pushd $d
popd
done

节约机时

grep User OUTCAR 后得出的结果被空格分成了4部分,时间信息在第4部分里面。User 或者 Elapsed都可以。

  • 一个合理的初始结构,可以避免很多意外的错误以及快速得到正确的结果。
  • 如果你设置的收敛标准比较严格,精度比较高,则每一个离子步需要更多的电子步数,需要的时间也会随着离子步的增加成线性关系增长。
  • ENCUT和KPOINTS文件需要测试,二者越大时间越久。
  • 这里我们就会需要一些其他的参数,NCORE和NPAR。
    NCORE:控制多少个核同时计算;
    NPAR:如何把计算任务分配到计算资源上面计算。
    它们之间的关系是:NCORE= 计算使用的核数 / NPAR。注意:这两个参数只能选取一个来使用:
    NCORE=单节点的核数,单节点的核数/2,单节点的核数/4…….个人使用经验是:NCORE = 单个节点核数 / 2 的时候,运行最省时间,设置也最方便。
    当然NCORE也需要测试。

消虚频

结构优化消虚频

X, Y, Z为体系中原子的坐标,也就是POSCAR中的内容。 dx, dy,dz为虚频对应的原子振动方向上的位移。我们将(dx, dy,dz)这个位移乘以一个介于0-1之间的校正因子,然后跟POSCAR中的坐标加起来即可。以z方向为例,校正因子设为0.1, 微调后的z方向坐标为: 10.709980 + (-0.676634)0.1。
脚本vib_correct.py如下:https://www.bigbrosci.com/2022/04/29/A29/

{.line-numbers}
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
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
from ase.io import read, write
import os

#获取虚频开始行
l_position = 0 #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
lines = f_in.readlines()
wave_num = 0.0
for num, line in enumerate(lines):
if 'f/i' in line:
wave_tem = float(line.rstrip().split()[6])
if wave_tem > wave_num: #获取最大的虚频
wave_num = wave_tem
l_position = num+2
# ASEPOSCAR
model = read('POSCAR')
model_positions = model.get_positions()
num_atoms = len(model)
#print(model_positions)

# 获取虚频对应的OUTCAR部分
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 72227248
#print(vib_lines)
vib_dis = []
for line in vib_lines:
#model_positions = [float(i) for i in line.rstrip().split()[:3]]
vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
vib_dis.append(vib_infor)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。

# 微调结构
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。
model.positions = new_positions

###保存结构
write('POSCAR_new', model, vasp5=True) # POSCAR_new是微调后的结构,用于下一步的计算(别忘了把POSCAR_new改成POSCAR)。

过渡态虚频

如果过渡态计算的时候出现好几个虚频,自己直接写个python脚本也能实现。不过使用ASE可以很方便地获取POSCAR中原子数目,保存原子固定的信息。参考如下。

{.line-numbers}
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
import os

#获取虚频开始行
l_position = 0 #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
lines = f_in.readlines()
wave_num = 0.0
for num, line in enumerate(lines):
if 'f/i' in line:
wave_tem = float(line.rstrip().split()[6])
if wave_tem > wave_num: #获取最大的虚频
wave_num = wave_tem
l_position = num+2
# 读POSCAR
with open('POSCAR') as f_pos:
lines_pos = f_pos.readlines()

# 获取虚频对应的OUTCAR部分
num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()])
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 72227248

model_positions = []
vib_dis = []
for line in vib_lines:
position = [float(i) for i in line.rstrip().split()[:3]]
vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
model_positions.append(position)
vib_dis.append(vib_infor)

# 微调结构
model_positions = np.array(model_positions)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。

###保存结构
f_out = open('POSCAR_new','w')
f_out.writelines(lines_pos[:8])
f_out.write('Cartesian\n')
for i in new_positions:
f_out.write(' '.join([str(coord) for coord in i]) + ' F F F\n')
f_out.close()

常用命令行

grep "33 81" -A3 FORCE_CONSTANTS | xargs echo -n| awk '{print $1,$2,$3+$7+$11}' >> trace.dat

计算原子间二阶力常数矩阵的迹

grep 'F=' OSZICAR|awk '{print $1, "\t"$3,"\t"$9}' >AIMD.dat

提取分子动力学温度能量数值

for i in 01 02 03 04 05;do cd $i && cp CONTCAR POSCAR && cd ..;done

CINEB计算重新计算

awk '{print $1,$2,$3+$4+$5,$6+$7+$8+$9+$10,$11}' PDOS > pdos.dat

pdos轨道加和

grep LOOP+ OUTCAR |awk '{sum+=$4} END {print "Average=",sum/NR}'

grep LOOP+ OUTCAR |awk '{sum+=$4} END {print "Sum=",sum}'

grep LOOP+ OUTCAR |wc -l

提取平均离子步时间,总时间,总步数

echo -e "4\n403\n49-96 121-144\n1\n"|vaspkit

vaspkit命令行

for i in *; do mkdir $i && cp INCAR KPOINTS POSCAR POTCAR sub.sh $i && sed -i "3s/0.05/$i/g" $i/INCAR ; done

批量修改INCAR

sed -i '5s/\<0\>/0.5/' filename
请将 filename 替换为您要修改的文件名。
5s/<0>/0.5/:
5:指定要修改的行号(第 5 行)。
s:表示替换操作。
<0>:匹配完整的数字 0,确保只替换单独的 0。
0.5:要替换成的新值。