博文
VASP-DFPT+phonopy 计算声子谱：最好的老师是英文手册
|||
1. VASP-INCAR文件-用于计算声子谱前的优化
System = Li-diamondoid_from_LiNH2_SG19_20GPa
PREC = Accurate
ENCUT = 650.0
LREAL = .False
ISMEAR = 1
SIGMA = 0.05
EDIFF = 1E-7
EDIFFG = -1E-3
NSW = 50
ISIF = 3
IBRION = 2
PSTRESS = 300
2. VASP INCAR文件： DFPT in VASP- used to calculated phonon spectrum
System = exp-LiN3_SG12
PREC = Accurate
ENCUT = 400.0
LREAL = .False
ISMEAR = 1
EDIFF = 1E-7
IBRION = 8
LWAVE = .FALSE
LCHARG = .FALSE
ADDGRID = .TRUE
3. 采用Phonopy扩胞得到超晶胞用于第二步的自洽计算
利用(1)中的INCAR文件精细优化得到结构后，采用Phonopy软件，扩胞，一般应包含接近100个原子，超胞的晶格尺寸应达到10A。具体步骤如下：
(1) copy CONTCAR POSCAR；拷贝第一步优化得到的CONTCAR文件为POSCAR，并修改POSCAR的格式，使第一行符合新版Phonopy软件的要求，即，第一行为体系中所包含元素对应的符号，顺序与POTCAR中一致。
(2) phonopy -d --dim ="4 4 2"；最好记下扩胞的大小4 4 2。得到disp.yamal、各POSCAR-00x及SPOSCAR文件；
(3)cp SPOSCAR POSCAR采用2中的INCAR进行优化。此时，应确保设置的KPOINTS产生的NKPTS个数达到20个左右。【NKPT保持20个左右时，kmesh可取在0.04~0.05之间的某一数值，对计算结果影响不大，但相比0.03取值，可显著提高计算速度！！】
(4) 优化完毕，check OUTCAR文件中的频率，如果频率行包含i，则表明存在虚频；需修改参数继续优化，以检查是否是由于参数设置不当引起的虚频。如：
grep 'cm-1' OUTCAR
(5) 如没有虚频，则可利用Phonopy提供的工具对OUTCAR进行处理，得到声子谱及态密度。
a. 得到态密度所需文件及操作
所需文件：
FORCE_CONSTANTS mesh.conf POSCAR( 扩胞前的unitcell)
其中mesh.conf内容设置如下：
ATOM_NAME =H Sc
DIM = 2 3 2
MP = 8 8 8
FORCE_CONSTANTS = READ
操作：
phonopy -p --factor=521.471 mesh.conf
b. 得到声子谱所需文件及操作
所需文件：band.conf FORCE_CONSTANTS POSCAR(扩胞前的unitcell)
其中 band.conf 文件内容设置如下：
ATOM_NAME = H Sc
DIM = 2 3 2
FORCE_CONSTANTS = READ
BAND= 0.000 0.000 0.000 0.000 0.000 0.500 -0.333 0.667 0.500 -0.333 0.667 0.000 0.000 0.000 0.000 0.000 0.500 0.000 0.000 0.500 0.500 -0.333 0.667 0.500
操作：1) phonopy -p --factor=521.471 band.conf 得到band.yaml
2) bandplot --gnuplot >PhononBAND.dat 将band.yaml转化为可用orgin绘图的数据文件
4. 手册解读：VASP-DFPT & phonopy calculation
VASP can calculate force constants in real space using DFPT. The procedure to calculate phonon properties may be as
follows:
（1）. Prepare unit cell structure named, e.g., POSCAR-unitcell. The following structure is a conventional unit
cell of NaCl.
Na Cl
1.00000000000000
5.6903014761756712 0.0000000000000000 0.0000000000000000
0.0000000000000000 5.6903014761756712 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.6903014761756712
4 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
(2). Prepare a perfect supercell structure from POSCAR-unitcell, e.g.,
% phonopy -d --dim="2 2 2" -c POSCAR-unitcell
(3). Rename SPOSCAR created in (2) to POSCAR (POSCAR-{number} and disp.yaml files will never be used.)
% mv SPOSCAR POSCAR
(4). Calculate force constants of the perfect supercell by running VASP with IBRION = 8 and NSW = 1.
An example of INCAR for insulator may be such like (just an example!):
PREC = Accurate
ENCUT = 500
IBRION = 8
EDIFF = 1.0e-08
IALGO = 38
ISMEAR = 0; SIGMA = 0.1 #对于金属ISMEAR=1,SIGMA=0.2
LREAL = .FALSE.
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
(5). After finishing the VASP calculation, confirm vasprun.xml contains hessian elements, and then create FORCE_CONSTANTS:
% phonopy --fc vasprun.xml
(6). Run phonopy with the original unit cell POSCAR-unitcell and setting tag FORCE_CONSTANTS = READ
or --readfc option, e.g., as found in example/NaCl-VASPdfpt
% phonopy --dim="2 2 2" -c POSCAR-unitcell band.conf
_
_ __ | |__ ___ _ __ ___ _ __ _ _
| ’_ | ’_ / _ | ’_ / _ | ’_ | | | |
| |_) | | | | (_) | | | | (_) || |_) | |_| |
| .__/|_| |_|___/|_| |_|___(_) .__/ __, |
|_| |_| |___/
1.1
Band structure mode
Settings:
Force constants: read
Supercell: [2 2 2]
Primitive axis:
[ 0. 0.5 0.5]
[ 0.5 0. 0.5]
[ 0.5 0.5 0. ]
Spacegroup: Fm-3m (225)
Paths in reciprocal reduced coordinates:
[ 0.00 0.00 0.00] --> [ 0.50 0.00 0.00]
[ 0.50 0.00 0.00] --> [ 0.50 0.50 0.00]
[ 0.50 0.50 0.00] --> [-0.00 -0.00 0.00]
[ 0.00 0.00 0.00] --> [ 0.50 0.50 0.50]
5. INCAR中部分参数注释-最好的老师是英文手册
(1) IBRION:
Setting IBRION=8 or IBRION=7 selects the calculation of the interatomic force constants using density functional perturbation theory.
For IBRION=8, symmetry is taken into account, whereas IBRION=7 neglects symmetry considerations and is thus usually significantly more expensive
(2) LREAL=false 表示投影操作在倒空间进行
(3)ADDGRID = .TRUE. The accuracy of forces can be further improved by specifying ADDGRID = .TRUE.
(4)
NSW = number of ionic steps
ISIF = what to relax
IBRION = which method to use for the relaxation
POTIM = size of trial step for ions
IBRION=-1表示
No update; ions are not moved, but NSW outer loops are performed. In each outer loop the electronic degrees of freedom are re-optimized (for NSW>0 this obviously does not make much sense, except for test purposes). If no ionic update is required use NSW=0 instead.
(5) ISMEAR(对于半导体和绝缘体：ISMEAR= 0，SIGMA = 0.05；ISMEAR= 0 表示采用Gaussian smearing 方法； 对于金属：ISMEAR= 1，SIGMA = 0.2;ISMEAR= 1 表示采用1阶smearing方法
（1）ISMEAR = -5 | -4 | -3 | -2 | 0 | N
SIGMA = width of the smearing in eV
Default
ISMEAR = 1
SIGMA = 0.2
ISMEAR determines how the partial occupancies fnk are set for each wavefunction. For the finite temperature LDA SIGMA determines the width of the smearing in eV.
ISMEAR:
-5 tetrahedron method with Bl¨ochl corrections
For the calculation of the total energy in bulk materials we recommend the tetrahedron method with Bl¨ochl corrections (ISMEAR=-5). This method also gives a good account for the electronic density of states (DOS).
The only drawback is that the methods is not variational with respect to the partial occupancies. Therefore the calculated forces and the stress tensor can be wrong by up to 5 to 10 % for metals.
For the calculation of phonon frequencies based on forces we recommend the method of Methfessel-Paxton (ISMEAR>0). For semiconductors and insulators the forces are correct, because partial occupancies do not vary and are zero or one.
For the calculations of the DOS and very accurate total energy calculations (no relaxation in metals) use the tetrahedron method (ISMEAR=-5).
For semiconductors or insulators use the tetrahedron method (ISMEAR=-5), if the cell is too large (or if you use only a single or two k-points) use ISMEAR=0 in combination with a small SIGMA=0.05.
For relaxations in metals always use ISMEAR=1 or ISMEAR=2 and an appropriated SIGMA value (the entropy term should be less than 1 meV per atom). Mind: Avoid to use ISMEAR>0 for semiconductors and insulators, since it might cause problems. For metals a sensible value is usually SIGMA= 0.2 (which is the default).
（2）RWIGS or LORBIT
If RWIGS or LORBIT (Wigner Seitz radii, see section 6.326.33) is set in the INCAR file, a lm- and site-projected DOS is calculated and also written to the file DOSCAR. One set of data is written for each ion, each set of data holds NDOS lines with the following data
energy s-DOS p-DOS d-DOS
or
energy s-DOS(up) s-DOS(down) p-DOS(up) p-DOS(dwn) d-DOS(up) d-DOS(dwn)
for the non spin-polarized and spin polarized case respectively. As before the written densities are understood as the difference of the integrated DOS between two pins.
(3) DOSCAR file
Mind: For relaxations, the DOSCAR is usually useless. If you want to get an accurate DOS for the final configuration, first copy CONTCARto POSCAR and continue with one static (ISTART=1; NSW=0) calculation.【静态自洽计算？】
（4）LORBIT
The default for LORBIT is .FALSE. (respectively 0).
This flag determines, together with an appropriate RWIGS (see section 6.32), whether the PROCAR or PROOUT files (see section 5.21) arewritten.
The file PROCAR contains the spd- and site projected wave function character of each band.
The wave function character is calculated, either by projecting the wavefunctions onto spherical harmonics that are non zero within spheres of a radius RWIGS around each ion (LORIT=1, 2), or using a quick projection scheme relying that works only for the PAW method (LORBIT=10,11,12, see below). If the LORBIT flag is not equal zero, the site and l-projected density of states is also calculated.
If the projector augmented wave method is used, LORBIT can also be set to 10, 11 or 12. This alternative setting selects a quick method for the determination of the spd- and site projected wave function character and does not require the specification of a Wigner-Seitz radius in the INCAR file (the RWIGS line is neglected in this case). The method works only for PAW POTCAR files and not for ultrasoft or norm conserving pseudopotentials. 即对每个能带的波函数进行s pd和site分解或投影；LORBIT = 10，不设RWIGS，针对PAW势的计算。
LORBIT = 10表示对波函数按原子site进行s 、p、d轨道分解或投影；并输出投影波函数到PROCAR or PROOUT 文件中，详见上面英文注释部分：This alternative setting selects a quick method for the determination of the spd- and site projected wave function character。什么样的quick method?
超胞法处理声子谱的步骤
1） 计算INCAR
#########################
for p in 001 002 003 004
do
cat > INCAR << EOF
PREC = Accurate
IBRION = -1
ENCUT = 500
EDIFF = 1.0e-08
#LSDA-plus-U
ISPIN=2
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = -1 3
LDAUU = 0.00 4.10
LDAUJ = 0.0 0.1
LDAUPRINT = 0
LMAXMIX = 6
#PSTRESS = 2450
#ISMEAR = 1; SIGMA = 0.2
ISMEAR = 0; SIGMA = 0.01
IALGO = 38
LREAL = .FALSE.
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
EOF
cp POSCAR-$p POSCAR
# start calculation
mpirun -np 8 mpivasp > log 2>/dev/null
cp OUTCAR OUTCAR_$p
cp vasprun.xml vasprun.xml-$p
cp STDOUT-$p
done
2） 后处理
In the post-process
0. Then create FORCE_SETS file using VASP interface:
phonopy -f vasprun.xml-00*
1. Force constants are calculated from the sets of forces
2. A part of dynamical matrix is built from the force constants
3. Phonon frequencies and eigenvetors are calculated from the dynamical matrices with the specified q-points.
For mesh sampling calculation, prepare the following setting file named, e
.g., mesh.conf:
ATOM_NAME = Si O
DIM = 2 2 3
MP = 8 8 8
The density of states (DOS) is plotted by:
% phonopy -p mesh.conf
Thermal properties are calculated with the sampling mesh by:
% phonopy -t mesh.conf
You should check the convergece with respect to the mesh numbers.
Thermal properties can be plotted by:
% phonopy -t -p mesh.conf
Projected DOS is calcualted by the following setting file named, e.g., pdos.conf:
ATOM_NAME = Si O
DIM = 2 2 3
MP = 8 8 8
PDOS = 1 2, 3 4 5 6
and plotted by:
% phonopy -p pdos.conf
Band structure is calculated with the following setting file named, e.g., band.conf by:
ATOM_NAME = Si O
DIM = 2 2 3
BAND = 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
The band structure is plotted by:
% phonopy -p band.conf
In either case, by setting the -s option, the plot is going to be saved in the PDF format. If you don’t need to plot DOS,
the (partial) DOS is just calculated using the --dos option.
Details
Following files are required in your working directory.
• POSCAR, and FORCE_SETS or FORCE_CONSTANTS
• disp.yaml is required to create FORCE_SETS.
In the case of finite difference approach, there are three steps.
1. Create supercells and introduce atomic displacements. Each supercell contains one atomic displacement. It is done by using -d option with --dim option that specifies supercell dimension. The files of supercells with
atomic displacements like as POSCAR-001, POSCAR-002, ..., are created in current directory (the file format
and names are different in WIEN2k mode.) by running phonopy. The files disp.yaml and SPOSCAR are also
created. The file SPOSCAR is the perfect supercell that contains no atomic displacement. This file is not usually
used.
2. Calculate forces on atoms of the supercells with atomic displacements. Currently phonopy has VASP and
WIEN2k interfaces to create FORCE_SETS. After obtaining forces on atoms that calculated by some calculator
(it’s out of phonopy), the forces are summarized in FORCE_SETS file following the format.
3. Calculate phonon related properties. See Features.
If you already have force constants, the first and second steps can be omitted. However your force constants have to be
converted to the format that phonopy can read. The VASP interface to convert force constants is prepared in phonopy.
不扩胞利用DFPT计算声子谱文件解读
dim111-from-gaotao]$ vi INCAR
System=XH3
PREC=Accurate
GGA=PE
ISMEAR=0
SIGMA=0.05
ENCUT=550
EDIFF=0.1E-5
EDIFFG=-0.01
LREAL=.FALSE.
IBRION=8
IALGO=38
LORBIT=11
LASPH=.TRUE.
#MAGMOM= 0 0 0 0 0 0 2 2
#KSPACING=.251327412287
#KGAMMA=.TRUE.
参数解读
Please mind, that the VASP.4.5 default is IALGO = 38 (a Davidson block iteration scheme). IALGO = 8 is not supported for copyright reasons in VASP.4.5, but IALGO = 38 is roughly 2 times faster for large systems than IALGO = 8 and at least as stable. You can select the algorithm also by setting ALGO= Normal | Fast | Very_Fast in the INCAR file (see Sec. 6.46).
This algorithm is the default in VASP.4.6 and VASP.5.X. It optimizes a subset of NSIMbands simultaneously (Sec. 6.48). The optimized bands are kept orthogonal to all other bands. If problems are encountered with the algorithm, try to decrease NSIM. Such problems are encountered, if linear dependencies develop in the search space. By reducing NSIM the rank of the search space is decreased.
| 11 | not read | DOSCAR and lm decomposed PROCAR file |
If the projector augmented wave method is used, LORBIT can also be set to 10, 11 or 12. This alternative setting selects a quick method for the determination of the spd- and site projected wave function character and does not require the specification of a Wigner-Seitz radius in the INCAR file (the RWIGS line is neglected in this case). The method works only for PAW POTCAR files and not for ultrasoft or norm conserving pseudopotentials.
| Default | ||
| LASPH | = | .FALSE. |
Usually VASP calculates only the spherical contribution to the gradient corrections inside the PAW spheres (non-sperical contributions for the LDA part of the potential and the Hartree potential are always included).
Using LASPH = .TRUE., VASP also includes non-spherical contributions from the gradient corrections inside the PAW spheres. For VASP.4.6, these contributions are only included in the total energy, after self-consistency has been reached disregarding the aspherical contributions in the gradient corrections.
For VASP.5.X the aspherical contributions are properly accounted for in the Kohn-Sham potential as well. This is essential for accurate total energies and band structure calculations for f-elements (e.g. ceria), all 3d-elements (transition metal oxides), and magnetic atoms in the 2nd row (B-F atom), in particular if LDA+U or hybrid functionals or meta-GGAs are used, since these functionals often result in aspherical charge densities.
https://blog.sciencenet.cn/blog-567091-675178.html
下一篇：【整理自好友lpf文章】用VASP计算能量态密度（DOS）和能带
