用VASP优化TiO2_anatase晶格常数（a, c）

Posted on | In

对一系列固定体积的结构进行弛豫，然后根据birch-murnaghan拟合E-V曲线，得到基态的平衡体积$V_0$。得到平衡体积之和再优化形状，得到a，c的值。

  • INCAR.relax

    1
    2
    3
    4
    5
    6
    ENCUT = 520 eV 
    PREC = High
    LREAL = .FALSE
    NSW = 999
    IBRION = 2 
    ISIF = 4

  • INCAR.static

    1
    2
    3
    4
    5
    6
    7
    PREC = High
    ISTART = 0 
    ICHARG = 2
    ISMEAR = -5
    ENCUT = 520 eV

  • KPOINTS

1
2
3
4
5
Automatic mesh
  0
Gamma
  5  5  5
  0  0  0
  • vasp.pbs
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
#!/bin/bash
#$ -S /bin/bash
#$ -cwd
#$ -j y
#$ -N test_opt
#$ -pe make 12
source /share/apps/intel/Compiler/11.1/073/bin/iccvars.sh intel64
source /share/apps/intel/Compiler/11.1/073/bin/ifortvars.sh intel64
source /share/apps/intel/impi/3.2.0.011/bin64/mpivars.sh
for a in 3.74 3.75 3.76 3.78 3.79 3.80 3.81 3.82 3.83 3.84
do 
cat > POSCAR <<!
Anatase
1.0
        $a         0.0000000000         0.0000000000
        0.0000000000         $a         0.0000000000
        0.0000000000         0.0000000000         9.737
   Ti    O
    4    8
Direct
     0.000000000         0.000000000         0.000000000
     0.500000000         0.500000000         0.500000000
     0.000000000         0.500000000         0.250000000
     0.500000000         0.000000000         0.750000000
     0.000000000         0.000000000         0.208059996
     0.500000000         0.500000000         0.708060026
     0.000000000         0.500000000         0.458059996
     0.500000000         0.000000000         0.958060026
     0.500000000         0.000000000         0.541939974
     0.000000000         0.500000000         0.041940004
     0.500000000         0.500000000         0.291940004
     0.000000000         0.000000000         0.791939974
!
cp INCAR.relax INCAR
echo "a = $a"; mpirun -r ssh -np 12 ~/bin/vasp5.3.3
cp CONTCAR POSCAR
cp INCAR.static INCAR
echo "a = $a"; mpirun -r ssh -np 12 ~/bin/vasp5.3.3
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f \n", $5 }'`
V=`grep "volume" OUTCAR | tail -1 | awk '{printf "%12.4f \n" , $5}'`
echo $a $V $E >> SUMMARY.dat
done

结果

  • SUMMARY.dat
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    3.74 136.2000 -105.833996
    3.75 136.9300 -105.865334
    3.76 137.6600 -105.892136
    3.78 139.1300 -105.928546
    3.79 139.8600 -105.944722
    3.80 140.6000 -105.955402
    3.81 141.3400 -105.960574
    3.82 142.0900 -105.960563
    3.83 142.8300 -105.954437
    3.84 143.5800 -105.949877

拟合E-V曲线

  • fit_os.py
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
#!/usr/bin/python
# coding=utf-8
import sys
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def read_data(filename):
    """ 读取能量和体积，文件的形式如下，第一列为晶格常数，第二列为体积，第三列为能量
        3.74 136.2000 -105.833996
        3.75 136.9300 -105.865334
        3.76 137.6600 -105.892136
        3.78 139.1300 -105.928546
        3.79 139.8600 -105.944722
        3.80 140.6000 -105.955402
        3.81 141.3400 -105.960574
        3.82 142.0900 -105.960563
        3.83 142.8300 -105.954437
        3.84 143.5800 -105.949877
    """
    data = np.loadtxt(filename)     
    return data[:,1], data[:,2] 
def eos_murnaghan(vol, E0, B0, BP, V0):
    """ 输入体积和方程参数，返回能量
    Birch-Murnaghan方程：Phys. Rev. B 1983, 28 (10), 5480–5486.
    """
    return E0 + B0*vol/BP*(((V0/vol)**BP)/(BP-1)+1) - V0*B0/(BP-1)
def fit_murnaghan(volume, energy):
    """ 拟合Murnaghan状态方程，返回最优参数 
    """
    # 用一元二次方程拟合，得到初步猜测参数
    p_coefs = np.polyfit(volume, energy, 2)
    # 抛物线的最低点 dE/dV = 0 ( p_coefs = [c,b,a] ) V(min) = -b/2a
    p_min = - p_coefs[1]/(2.*p_coefs[0])
    # warn if min volume not in result range 
    if (p_min < volume.min() or p_min > volume.max()):
        print "Warning: minimum volume not in range of results"
    # 从抛物线最低点估计基态能量    
    E0 = np.polyval(p_coefs, p_min)
    # 体积模量估计
    B0 = 2.*p_coefs[2]*p_min
    # 初步猜测参数 (BP通常很小，取4)
    init_par = [E0, B0, 4, p_min]
    print "guess parameters:"
    print " V0     =  {:1.4f} A^3 ".format(init_par[3])
    print " E0     =  {:1.4f} eV  ".format(init_par[0])
    print " B(V0)  =  {:1.4f} eV/A^3".format(init_par[1])
    print " B'(VO) =  {:1.4f} ".format(init_par[2])
    best_par, cov_matrix = curve_fit(eos_murnaghan, volume, energy, p0 = init_par)
    return best_par
def fit_and_plot(filename):
    """ 读取文件数据，拟合Murnaghan状态方程，返回最优参数和E-V图形
    """
    # 从文件读取数据
    volume, energy = read_data(filename)   
    # 用Murnaghan状态方程拟合数据
    best_par = fit_murnaghan(volume, energy)
    # 输出最优参数   
    print "Fit parameters:"
    print " V0     =  {:1.4f} A^3 ".format(best_par[3])
    print " E0     =  {:1.4f} eV  ".format(best_par[0])
    print " B(V0)  =  {:1.4f} eV/A^3".format(best_par[1])
    print " B'(VO) =  {:1.4f} ".format(best_par[2])
    # 用拟合的参数生成Murnaghan模型
    m_volume = np.linspace(volume.min(), volume.max(), 1000) 
    m_energy = eos_murnaghan(m_volume, *best_par) 
    # 画E-V图
    lines = plt.plot(volume, energy, 'ok', m_volume, m_energy, '--r' )
    plt.xlabel(r"Volume [$\rm{A}^3$]")
    plt.ylabel(r"Energy [$\rm{eV}$]")
    return  best_par, lines
fit_and_plot("SUMMARY.dat")
plt.draw()
plt.show()

$V_0=141.863A^3$，得到平衡体积之和再固定体积优化形状即可得到晶格常数（a, c）。

  • INCAR

    1
    2
    3
    4
    5
    6
    ENCUT = 520 eV 
    PREC = High
    LREAL = .FALSE 
    NSW = 999
    IBRION = 2 
    ISIF = 4

  • POSCAR

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    Anatase
    -141.863
            3.7844998837         0.0000000000         0.0000000000
            0.0000000000         3.7844998837         0.0000000000
            0.0000000000         0.0000000000         9.5143003464
       Ti    O
        4    8
    Direct
         0.000000000         0.000000000         0.000000000
         0.500000000         0.500000000         0.500000000
         0.000000000         0.500000000         0.250000000
         0.500000000         0.000000000         0.750000000
         0.000000000         0.000000000         0.208059996
         0.500000000         0.500000000         0.708060026
         0.000000000         0.500000000         0.458059996
         0.500000000         0.000000000         0.958060026
         0.500000000         0.000000000         0.541939974
         0.000000000         0.500000000         0.041940004
         0.500000000         0.500000000         0.291940004
         0.000000000         0.000000000         0.791939974

  • KPOINTS

    1
    2
    3
    4
    5
    Automatic mesh
      0
    Gamma
      5  5  5
      0  0  0

  • vasp.pbs

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    #!/bin/bash
    #$ -S /bin/bash
    #$ -cwd
    #$ -j y
    #$ -N isif4
    #$ -pe make 12
    source /share/apps/intel/Compiler/11.1/073/bin/iccvars.sh intel64
    source /share/apps/intel/Compiler/11.1/073/bin/ifortvars.sh intel64
    source /share/apps/intel/impi/3.2.0.011/bin64/mpivars.sh
    mpirun -r ssh -np 12 ~/bin/vasp5.3.3

最后得到CONTCAR 

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
Anatase                                 
   1.01350313575973     
     3.7709606715685475    0.0000000000000000    0.0000000000000000
     0.0000000000000000    3.7709606715685475    0.0000000000000000
     0.0000000000000000    0.0000000000000000    9.5827430546962322
   Ti   O 
     4     8
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.5000000000000000  0.5000000000000000  0.5000000000000000
  0.0000000000000000  0.5000000000000000  0.2500000000000000
  0.5000000000000000  0.0000000000000000  0.7500000000000000
  0.0000000000000000  0.0000000000000000  0.2066888233929389
  0.5000000000000000  0.5000000000000000  0.7066888533929343
  0.0000000000000000  0.5000000000000000  0.4566888233929389
  0.5000000000000000  0.0000000000000000  0.9566888533929343
  0.5000000000000000  0.0000000000000000  0.5433111466070657
  0.0000000000000000  0.5000000000000000  0.0433111766070611
  0.5000000000000000  0.5000000000000000  0.2933111766070611
  0.0000000000000000  0.0000000000000000  0.7933111466070657
 
  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

体积是能对上的： $V_0 = 141.863 = 1.01350313575973^3 \times 3.7709606715685475 \times 3.7709606715685475 \times 9.5827430546962322$

所以最后$a=1.013 \times 3.771=3.820, c=1.013 \times 9.583=9.707$

文献计算值a=3.786, c=9.737（Phys. Rev. B 2001, 63 (15), 155409.）

实验值a=3.782, c=9.502（J. Am. Chem. Soc. 1987, 109 (12), 3639–3646.）