截断能 截断能设置在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 mesh0 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 ]; thenecho -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 ] ; thentail -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-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 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 npfrom ase.io import read, writeimport 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} 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 npimport 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:要替换成的新值。