VASP 晶格常数优化
状态方程(EOS)
Birch-Murnaghan1
a relationship between the volume of a body and the pressure to which it is subjected.
Birch proposed this equation based on the work of Francis Dominic Murnaghan of Johns Hopkins University published in 1944, so that the equation is named in honor of both scientists.
$$E(V)=E_{0}+\frac{9 V_{0} B_{0}}{16}\left{\left[\left(\frac{V_{0}}{V}\right)^{\frac{2}{3}}-1\right]^{3} B_{0}^{\prime}+\left[\left(\frac{V_{0}}{V}\right)^{\frac{2}{3}}-1\right]^{2}\left[6-4\left(\frac{V_{0}}{V}\right)^{\frac{2}{3}}\right]\right}$$
vaspkit EOS
[1]: Murnaghan EOS [2]: Birch-Murnaghan 3rd-order EOS [3]: Birch-Murnaghan 4th-order EOS [4]: Mie–Gruneisen [5]: Pack–Evans–James [6]: Poirier-Tarantola EOS [7]: Tait EOS [8]: Vinet EOS [9]: Natural strain 3rd-order EOS [10]: Natural strain 4th-order EOS [11]: Universal EOS [12]: Cubic polynomial in (V-V0)
重复金刚石的例子
实验: A=3.57Å (300K)(3.56683 \AA) 拟合区间影响很大，拟合方程选择不重要，VASP-D3,D2,PBE有不同 同时拟合也没有那么贴切
下图也是把所有的拟合方法都尝试了一下子
当然原图也可以重复出来
迭代法
我们计算资源充足，干脆不拟合了，换迭代法 按照VASP官方指导 我们设置 ISIF=4 The method we have used quite often in the past is to relax the structure (cell shape and internal parameters) for a set of fixed volumes .
标准化POSCAR, 设POSCAR第二行的Scale系数为1
代码
1% matlab bash for run vasp
2warning off;
3tic;
4vasprun = '!mpirun -np 16 vasp_std >log 2>err';
5format long;
6[Rm,sites,Atom_name,Atom_num,~,factor]=POSCAR_readin('POSCAR');
7natom = length(sites);
8V = abs(dot(Rm(:,3),cross(Rm(:,1),Rm(:,2))))*factor;
9fid = fopen('matlabrun.log','w');
10fprintf(fid,'Volume is :%6.5f \\AA ^3 \n',V);
11copyfile('POSCAR','POSCAR.bk0');
12copyfile('POSCAR','POSCAR.bk');
13copyfile('POSCAR','POSCAR.bk_converge');
14% INCAR.RELAX INCAR.static
15
16% mkdir
17mkdir('RELAX');
18mkdir('STATIC');
19WORKDIR = pwd;
20
21DecimalDepth = 3;
22Vnum = 11;
23BeginingCrange = [-0.1,0.1];
24Width = (BeginingCrange(2)-BeginingCrange(1))/2;
25DepthBase = (Vnum-1)/2;
26lastVolume = 1;
27lastVrange = BeginingCrange+lastVolume;
28TOTVlist = zeros(DecimalDepth*Vnum,1);
29TOTElist = zeros(DecimalDepth*Vnum,1);
30
31fid2 =fopen('EvsV.dat','w');
32
33for j = 1:DecimalDepth
34 Vrangej = lastVrange;
35 TOTElistj = zeros(Vnum,1);
36 TOTclistj = zeros(Vnum,1);
37 count = 0;
38 Vrange = linspace(Vrangej(1),Vrangej(2),Vnum);
39 copyfile('POSCAR.bk_converge','POSCAR.bk');
40 Emin = [];
41 for v_i = Vrange
42 count = count + 1;
43 copyfile('POSCAR.bk','POSCAR');
44 evalstring = ['!sed -i "2c ',num2str(v_i),'" POSCAR'];
45 [Rm,~,~,~,~,~] = POSCAR_readin('POSCAR','digits',10);
46 Rstruct = park.Rm2abc(Rm*v_i);
47 eval(evalstring);
48 % cd RELAX
49 cd('RELAX');
50 !cp ../POSCAR ../POTCAR ../KPOINTS .
51 !cp ../INCAR.RELAX ./INCAR
52 eval(vasprun);
53 !grep 'free energy' OUTCAR |tail -n 1|awk '{print $5}' >tmpE.dat
54 E_RELAX = importdata('tmpE.dat')/natom;
55 % cd back
56 cd(WORKDIR);
57 % cd STATIC
58 cd('STATIC');
59 !cp ../RELAX/CONTCAR POSCAR
60 !cp ../POSCAR ../POTCAR ../KPOINTS .
61 !cp ../INCAR.STATIC ./INCAR
62 eval(vasprun);
63 !grep 'free energy' OUTCAR |tail -n 1|awk '{print $5}' >tmpE.dat
64 E_STATIC = importdata('tmpE.dat')/natom;
65 if isempty(Emin)
66 Emin = E_STATIC;
67 !cp POSCAR ../POSCAR.bk_converge
68 else
69 if E_STATIC < Emin
70 Emin = E_STATIC;
71 !cp POSCAR ../POSCAR.bk_converge
72 end
73 end
74 TOTElistj(count) = E_STATIC;
75 [Rm_optimize,~,~,~,~,vfactor]= POSCAR_readin('CONTCAR','digits',10);
76 Rstruct_optimize = park.Rm2abc(Rm_optimize*vfactor);
77 % cd back
78 cd(WORKDIR);
79 fprintf(fid,'--- %s ---%2d-%2d [V:%7.5f E: %13.9f]\n',datetime,j,count,v_i^3,E_STATIC);
80 fprintf(fid,'LatticeC: a b c alpha beta gamma \n');
81 fprintf(fid,'%s %7.5f %7.5f %7.5f %7.5f %7.5f %7.5f\n',' init',Rstruct.a,Rstruct.b,Rstruct.c,Rstruct.alpha,Rstruct.alpha,Rstruct.gamma);
82 fprintf(fid,'%s %7.5f %7.5f %7.5f %7.5f %7.5f %7.5f\n','converge',Rstruct_optimize.a,Rstruct_optimize.b,Rstruct_optimize.c,Rstruct_optimize.alpha,Rstruct_optimize.alpha,Rstruct_optimize.gamma);
83 % write data
84 fprintf(fid2,'%7.5f %11.7f %13.9f ',v_i,v_i^3*V,E_STATIC);
85 fprintf(fid2,'%7.5f %7.5f %7.5f %7.5f %7.5f %7.5f\n',Rstruct_optimize.a,Rstruct_optimize.b,Rstruct_optimize.c,Rstruct_optimize.alpha,Rstruct_optimize.alpha,Rstruct_optimize.gamma);
86 end
87 [Eminj,Eminseq] = min(TOTElistj);
88 Vmin = Vrange(Eminseq);
89 lastVolume = Vmin;
90 VWidth = Width/DepthBase^j;
91 lastVrange = [Vmin-VWidth,Vmin+VWidth];
92 TOTVlist(((j-1)*Vnum+1):j*Vnum) = Vrange;
93 TOTElist(((j-1)*Vnum+1):j*Vnum) = TOTElistj;
94end
95fclose(fid);
96fclose(fid2);
97 TOTVlist_r = (TOTVlist).^3*V/natom;
98save('EvsC.mat');
99toc;
Diamond
迭代法与上树拟合差不太多
Graphite
贝叶斯优化
VASP参数参考
VASP中与晶格常数优化相关的主要参数有 ISIF 以及 IBRION 重要的一点在于 ENCUT的设置
Before you perform relaxations in which the volume or the cell shape is allowed to change you must read and understand the section on energy vs. volume, volume relaxations, and Pulay stress . In general volume changes should be done only with an increased energy cutoff, i.e.,
ENCUT=1.3×max(ENMAX)orPREC=High.
Pulay stress
ISIF
| ISIF | calculate | degrees-of-freedom | |||
|---|---|---|---|---|---|
| forces | Stress tensor | positions | cell shape | cell volume | |
| 0 | yes | no | yes | no | no |
| 1 | yes | trace only | yes | no | no |
| 2 | yes | yes | yes | no | no |
| 3 | yes | yes | yes | yes | yes |
| 4 | yes | yes | yes | yes | no |
| 5 | yes | yes | no | yes | no |
| 6 | yes | yes | no | yes | yes |
| 7 | yes | yes | no | no | yes |
F. Birch, Finite Elastic Strain of Cubic Crystals, Phys. Rev. 71, 809 (1947). ↩︎
