1.5 提高效率与准确性

发表于 2023-05-09 更新于 2024-09-30 分类于 DFT学习

截断能

截断能设置在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}

#!/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}

#!/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}

#!/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}

#!/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}

#!/bin/sh
root_path='pwd'

#在每个job文件夹里进行静态计算。
for i in {001..1316} #按顺序从job-001到116做静态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

```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}

# -*- 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
# ASE读POSCAR
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] #振动部分 7222到7248行
#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}

#!/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] #振动部分 7222到7248行

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：要替换成的新值。