VASP(1)_参数测试_晶格常数测试

62 阅读 0 评论 41 点赞

我是靠谱客的博主潇洒皮带，最近开发中收集的这篇文章主要介绍VASP(1)_参数测试_晶格常数测试，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

晶格常数测试 (Equation of state method)

必要输入文件

run_a0.sh
POTCAR
INCAR 和 KPOINTS可在run_a0.sh中直接设置，也可以单独给出。
EOS.in

三维立方晶格脚本示例：
Si

#!/bin/sh

cat > INCAR.relax <<!
Global Parameters
  ISTART =  0            (Read existing wavefunction; if there)
  ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
  LREAL  = F             (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) 
  # LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
  # LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
  # NELECT =             (No. of electrons: charged cells; be careful)
  # LPLANE = .TRUE.      (Real space distribution; supercells)
  # NPAR   = 4           (Max is no. nodes; don't set for hybrids)
 
Electronic Relaxation
  ISMEAR =  0            (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-08        (SCF energy convergence; in eV) 
  GGA  =  PE             (PBEsol exchange-correlation)
 
Ionic Relaxation
  NELMIN =  6            (Min electronic SCF steps) 
  NSW    =  60           (Max electronic SCF steps)
  IBRION =  2            (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
  ISIF   =  4            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
  EDIFFG = -0.001        (Ionic convergence; eV/AA)
  # ISYM =  2            (Symmetry: 0=none; 2=GGA; 3=hybrids)
!

cat > INCAR.static <<!
Global Parameters
  ISTART =  0            (Read existing wavefunction; if there)
  ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
  LREAL  = F             (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) 
  # LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
  # LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
  # NELECT =             (No. of electrons: charged cells; be careful)
  # LPLANE = .TRUE.      (Real space distribution; supercells)
  # NPAR   = 4           (Max is no. nodes; don't set for hybrids)
  GGA   = PE
 
Static Calculation
  ISMEAR = -5            (tetrahedron method for DOS) 
  #LORBIT =  11          (PAW radii for projected DOS)
  #NEDOS  =  2001        (DOSCAR points) 
  NELM   =  60           (Max electronic SCF steps)
  EDIFF  =  1E-08        (SCF energy convergence; in eV)
  EDIFFG = -0.001
!

cat > KPOINTS <<!
A
0
M
9  9  9
0  0  0
!

echo 'a0' 'volume' 'free_energy(eV)' >ev.out

for i in $(seq 5.00 0.05 5.90)
do

cat > POSCAR <<!
Si8
1.0000000000
    $i                  0.0000000000        0.0000000000
    0.0000000000        $i                  0.0000000000
    0.0000000000        0.0000000000        $i
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
!

#优化计算
cp INCAR.relax INCAR
echo "a=$i"; time mpirun -n 16 vasp_std
cp CONTCAR POSCAR
rm INCAR

#静态计算
cp INCAR.static INCAR
echo "a=$i"; time mpirun -n 16 vasp_std
V=$(grep "volume" OUTCAR | tail -1 | awk '{printf "%12.9f n", $5 }')
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f n", $5 }')
echo $i $V $E >> ev.out
rm INCAR
done

提交脚本运算后得到ev.out。
利用vaspkit软件进行状态方程拟合。

vaspkit需要准EOS.in文件，文件格式及说明如下：

 cname               : name of crystal up to 256 characters
 natoms              : number of atoms in unit cell
 etype               : equation of state type (see below)
 vplt1, vplt2, nvplt : volume interval over which to plot energy, pressure etc.
                       as well as the number of points in the plot
 nevpt               : number of energy-volume points to be inputted
 vpt(i) ept(i)       : energy-volume points (VASP units)

Note that the input units are VASP default untis (i.e., A^3 and eV).

The equations of state currently implemented are:
 1. Universal EOS (Vinet P et al., J. Phys.: Condens. Matter 1, p1941 (1989))
 2. Murnaghan EOS (Murnaghan F D, Am. J. Math. 49, p235 (1937))
 3. Birch-Murnaghan 3rd-order EOS (Birch F, Phys. Rev. 71, p809 (1947))
 4. Birch-Murnaghan 4th-order EOS
 5. Natural strain 3rd-order EOS (Poirier J-P and Tarantola A, Phys. Earth
    Planet Int. 109, p1 (1998))
 6. Natural strain 4th-order EOS
 7. Cubic polynomial in (V-V0)

参考例子

Si
8
3
124.00 206 500
19
125.000000000 -39.404864520
128.790000000 -40.307996560
132.650000000 -41.074573950
136.590000000 -41.714702140
140.610000000 -42.237893500
144.700000000 -42.653098650
148.880000000 -42.968706720
153.130000000 -43.192650580
157.460000000 -43.332358630
161.880000000 -43.394826990
166.380000000 -43.386635660
170.950000000 -43.313964240
175.620000000 -43.182635890
180.360000000 -42.998104870
185.190000000 -42.765493980
190.110000000 -42.489611560
195.110000000 -42.174957610
200.200000000 -41.825750520
205.380000000 -41.445910010

运行 vaspkit -task 205
得到主要的输出文件 PARAM.out 如下：

Si
 
Birch-Murnaghan 3rd-order EOS
Birch F, Phys. Rev. 71, p809 (1947)
 
(Default VASP units: eV, Angstrom etc.) 
 
 V0 (A^3)          =            163.6206779
 E0 (eV)           =           -43.39626828
 B0                =           0.2982734905E-02
 B0'               =            4.260966253
 
 B0 (GPa)          =            87.75507591

通过平衡体积可以得到晶格常数 $a_0=5.4694803016$

如果不想状态方程拟合，在较高精度下直接使用 $I S I F = 3$ 的参数直接弛豫优化晶胞。
INCAR示例如下：

Global Parameters
  ISTART =  0            (Read existing wavefunction; if there)
  ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
  LREAL  = F             (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) 
  # LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
  # LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
  # NELECT =             (No. of electrons: charged cells; be careful)
  # LPLANE = .TRUE.      (Real space distribution; supercells)
  # NPAR   = 4           (Max is no. nodes; don't set for hybrids)
 
Electronic Relaxation
  ISMEAR =  0            (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-08        (SCF energy convergence; in eV) 
  GGA  =  PE             (PBEsol exchange-correlation)
 
Ionic Relaxation
  NELMIN =  6            (Min electronic SCF steps) 
  NSW    =  60           (Max electronic SCF steps)
  IBRION =  2            (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
  ISIF   =  3            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
  EDIFFG = -0.001        (Ionic convergence; eV/AA)
  # ISYM =  2            (Symmetry: 0=none; 2=GGA; 3=hybrids)

注：在新版本的vaspkit中，给出了更多状态方程及其整个新的拟合流程，具体可参考程序实例。

二维材料的六方晶格示例：
Tl2O


#!/bin/sh

cat > INCAR.relax <<!
Global Parameters
  ISTART =  0            (Read existing wavefunction; if there)
  ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
  LREAL  = F             (Projection operators: automatic)
  ENCUT  =  600          (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) 
  # LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
  # LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
  # NELECT =             (No. of electrons: charged cells; be careful)
  # LPLANE = .TRUE.      (Real space distribution; supercells)
  # NPAR   = 4           (Max is no. nodes; don't set for hybrids)

Electronic Relaxation
  ISMEAR =  0            (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-08        (SCF energy convergence; in eV) 
  GGA  =  PE             (PBEsol exchange-correlation)
 
Ionic Relaxation
  NELMIN =  6            (Min electronic SCF steps) 
  NSW    =  100          (Max electronic SCF steps)
  IBRION =  2            (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
  ISIF   =  2            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
  EDIFFG = -0.001         (Ionic convergence; eV/AA)
  # ISYM =  2            (Symmetry: 0=none; 2=GGA; 3=hybrids)
!

cat > INCAR.static <<!
Global Parameters
  ISTART =  0            (Read existing wavefunction; if there)
  ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
  LREAL  = F             (Projection operators: automatic)
  ENCUT  =  600          (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) 
  # LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
  # LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
  # NELECT =             (No. of electrons: charged cells; be careful)
  # LPLANE = .TRUE.      (Real space distribution; supercells)
  # NPAR   = 4           (Max is no. nodes; don't set for hybrids)
  GGA   = PE

Static Calculation
  ISMEAR = -5           (tetrahedron method for DOS) 
  #SIGMA = 0.05
  #LORBIT =  11           (PAW radii for projected DOS)
  #NEDOS  =  2001         (DOSCAR points) 
  NELM   =  60           (Max electronic SCF steps)
  EDIFF  =  1E-08        (SCF energy convergence; in eV)
  EDIFFG = -0.001

!

cat > KPOINTS <<!
A
0
G
12 12 1
0  0  0
!

for i in $(seq 3.10 0.05 4.10)
do

j=`echo "$i * 0.5000000000"|bc`
k=`echo "$i * 0.8660254038"|bc`

cat > POSCAR <<!
Tl2O
1.0
 $i               0.00000000000    0.00000000000
-$j               $k               0.00000000000
 0.00000000000    0.00000000000   25.00000000000
Tl O
2  1
Direct
  0.3333330099999969  0.6666670059999973  0.5606771909028483
  0.6666669459999994  0.3333329709999973  0.4393227710971557
  0.0000000000000000  0.0000000000000000  0.5000000000000000
!

cp INCAR.relax INCAR
echo "a=$j"; time mpirun -np 16 vasp_std
cp CONTCAR POSCAR

rm INCAR
cp INCAR.static INCAR
echo "a=$j"; time mpirun -np 16 vasp_std

V=$(grep "volume" OUTCAR | tail -1 | awk '{printf "%12.9f n", $5 }')
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f n", $5 }')
echo $i $V $E >> ev.dat
rm POSCAR
done