博文
再谈phonopy计算声子谱虚频的处理方法
||
关注：
1） 虚频的物理意义？波矢量？ 特征值？ 动力学矩阵的含义？
2） 什么样的虚频可以调节？ 振动模式的可视化
3） 什么样的不可调节，必须更新结构预测
4） 方法有三个：
a. modulation 调制：得到一系列结构，计算单点能获得能量最低结构
【尚未尝试】
b. 分子动力学处理：扩胞后，分子动力学跑以期得到能量更低结构，找出对称性，优化，扩胞后再计算声子谱
分子动力学跑的过程中可以设计一个脚本把每一步的POSCAR保留下来，看能量在哪一步最低，找到能量最低的构型，拿出来优化即可得到一个较为稳定的构型。
用到的命令：grep 'energy without ' OUTCAR
运行分子动力学处理虚频问题，在某一步得到了最低能量值，请问如何编写一个脚本把该能量值对于的POSCAR提取出来？
XDATCAR vasprun.xml 都有某一步的位型。
可以提取出来。
感觉，前面的晶格矢量都是一样的，改变的只是原子位置，
Direct configuration= 5
.....
Direct configuration= 2000
vaspun 中，我也没注意过。
我一般在XDATCAR里找。
small 0:36:16
嗯，XDATCAR中较直观，谢谢；祝美国工作愉快
c. 类似于软膜的处理：
找出最大虚频如K(1/3,1/3,0)，扩胞成(3,3,1)，即将K点移动至Gamma点；
计算gamma点声子谱(IBRION=6)，找到最大虚频点，将x/y/z+dx/dy/dz，得到笛卡尔坐标形式的POSCAR；
优化后，再计算gamma点声子谱，重复上述步骤，反复进行x/y/z+dx/dy/dz，直至gamma点虚频可忽略；
从而找到能量最低的且动力学稳定的结构，导入castep找出对称性，优化后，扩胞后再计算声子谱 ；如确认消除虚频，则在不同压力点优化，得到焓值，与之前结构的焓值比较，看是否也热力学更稳定
System = ScH3-P63mmc-0G
PREC = Accurate
ENCUT = 400.0
LREAL = .False
ISMEAR = 1
EDIFF = 1E-7
#PSTRESS = 2500
IBRION = 6
LWAVE = .FALSE
LCHARG = .FALSE
ADDGRID = .TRUE
5） 计算声子谱的方法有三种
a. 超胞法
b. DFPT方法
c. 只计算gamma点声子频率，IBRION=6
分子动力学方式处理虚频脚本摘录1：
Thanks to YKT:
SYSTEM = bcc
EDIFF = 2e-5
NSW = 2000
IBRION = 0
POTIM = 1.5
TEBEG =400
ISMEAR = 0
SIGMA=0.05
LCHARG=F
LWAVE=F
NELMIN=6
KBLOCK=1
NBLOCK=5
ISIF=2
SMASS=2.0
modulation 调制处理虚频脚本摘录1：
Thanks to YKT:
for a in 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8
do
cat > INPHON_$a << EOF
ATOM_NAME =N H
DIM= 4 4 4
MODULATION=1 2 1, 0.00 0.50 0.00 1 $a 【 1 $a ； 1的含义仍不清楚】
EOF
phonopy INPHON_$a --tolerance=0.001
cp MPOSCAR-001 MPOSCAR-001-$a
done
脚本摘录2：
Thanks must be give to ZCY:
for a in 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6
do
cat > INPHON_$a << EOF
ATOM_NAME =B Mo
DIM= 2 2 2
MODULATION= 0. 0. 0. 1 1 1,1 $a
【gama点0,0,0；不用扩胞所以是1 1 1； 沿1个方向，按a所代表的数值移动；若两个方向需移动，则写成1 $a， 2 $a】
EOF
phonopy INPHON_$a 【这个命令是什么意思？】
cp MPOSCAR-001 MPOSCAR-001-$a 【得到一系列POSCAR....】
cp MPOSCAR-002 MPOSCAR-002-$a
cp MPOSCAR MPOSCAR-$a
done
下面的脚本表示对上面所得到的MPOSCAR-002-$a做单点计算得到能量最低的结构
vi job
for pressure in 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6
do
cp MPOSCAR-001-$pressure POSCAR
mpirun -np 8 ~/bin/vasp </dev/null > log
cp OUTCAR OUTCAR_$pressure
done
单点计算INCAR设置如下：
cat INCAR
SYSTEM = Various- local optimisation
PREC = Accurate
ENCUT = 500.0
EDIFF = 1e-5
#IBRION = 1
#ISIF =3
#NSW = 100
ISMEAR = 1
SIGMA = 0.2
POTIM = 0.100
#No writing charge density and wavefunction
LCHARG = FALSE
LWAVE = FALSE
#Target Pressure
PSTRESS = 0
#Finer optimization
EDIFFG = 1e-5
#Parallelisation
NPAR = 1
EDIFFG = -0.001
网络摘录：
http://emuch.net/html/201407/7626449.html
1. 在G点存在虚频，这是最好解决的，直接做结构优化消除
2. 在非G点的高对称性点存在虚频，需要做相应的单胞增倍。
如果不是高对称性位置存在虚频，这个处理其实类似于上面2中的处理，比如说(1/3,0,0)这种，变成3*1*1单胞还是可以解决的。【0都弄成1？】
谢谢你的耐心解答，这个就是我计算出来的声子谱。它并不是在某个点处产生虚频，而是一部分的区域都有虚频，虚频最大值处在0 0 0.28处。并不是整数倍，这样说来 是无法解决的了？
_C84V_8G(WBY(JA}4AW4~AL.jpg
...
你可以看看声学简正模对应的原子移动是怎么样的，然后在结构中人为引入这种畸变吧。因为计算软件在结构优化过程中一般是固定初始结构的空间群，所以往往无法产生自发结构对称性的破坏。
摘自phonopy手册
Create modulated structure¶
The MODULATION tag is used to create a crystal structure with displacements along normal modes at q-point in the specified supercell dimension.
Atomic displacement of the j-th atom is created from the real part of the eigenvectors【特征向量】 with amplitudes【振幅】 and phase factors as
![.frac{A} { .sqrt{m_j} } .operatorname{Re} .left[ .exp(i.phi)
.mathbf{e}_j .exp( .mathbf{q} .cdot .mathbf{r}_j ) .right],](http://phonopy.sourceforge.net/_images/math/d8ac653863a6fa27354337d883c5612d14a88542.png)
where 
is the amplitude, 【从哪里可获取这些数值？】
is the phase,
and 
is the mass of the j-th atom,
is the q-point specified,
is the position of the j-th atom and in the l-th unit cell,
and 
is the j-th part of eigenvector.
Convention of eigenvector or dynamical matrix employed in phonopy is shown in Dynamical matrix.
If several modes are specified as shown in the example above, they are overlapped on the structure.
The output filenames are MPOSCAR.... Each modulated structure of a normal mode is written in MPOSCAR-<number> where the numbers correspond to the order of specified sets of modulations.
MPOSCAR is the structure where all the modulations are summed.
MPOSCAR-orig is the structure without containing modulation, but the dimension is the one that is specified.
Some information is written into modulation.yaml.
（1）The first three values correspond to the supercell dimension. 【3 3 1】
（2） The following values are used to describe how the atoms are modulated.
Multiple sets of modulations can be specified by separating by comma ,.
In each set, the first three values give a Q-point in the reduced coordinates in reciprocal space
【倒空间可约坐标？】.
Then the next three values【1/3 1/3 0 1 2】 are the band index from the bottom with ascending 【上升的】order, amplitude, and phase factor 【在此处phase factor没给出】in degrees.
The phase factor is optional. If it is not specified, 0 is used.
Before multiplying user specified phase factor, the phase of the modulation vector is adjusted as the largest absolute value, 
, of element of 3N dimensional modulation vector to be real.
The complex modulation vector is shown in modulation.yaml.
Characters of irreducible representations (IRs) of phonon modes are shown. For this calculation, a primitive cell has to be used. If the input unit cell is a non-primitive cell, it has to be transformed to a primitive cell using PRIMITIVE_AXIS tag.
The first three values gives a q-point in reduced coordinates to be calculated. The degenerated modes are searched only by the closeness of frequencies. The frequency difference to be tolerated is specified by the fourth value in the frequency unit that the user specified.
Only the databases of IRs for a few point group types at the 
point are implemented. If the database is available, the symbols of the IRs and the rotation operations are shown.
Irreducible representations are shown along with character table.
Show irreps of little co-group (point-group of wavevector) instead of little group.
Dynamical matrix¶
In phonopy, a phase convention 【约定，惯例】of dynamical matrix is used as follows:
![D_{.alpha.beta}(jj',.mathbf{q}) = .frac{1}{.sqrt{m_j m_{j'}}}
.sum_{l'}
.Phi_{.alpha.beta}(j0, j'l')
.exp(i.mathbf{q}.cdot[.mathbf{r}(j'l')-.mathbf{r}(j0)]),](http://phonopy.sourceforge.net/_images/math/50e116d4df99a736363b96f4fb08c46b0f2b4172.png)
where 
is the atomic mass and is the wave vector【波矢量？！】. An equation of motion is written as
![.sum_{j'.beta} D_{.alpha.beta}(jj',.mathbf{q}) e_.beta(j', .mathbf{q}.nu) =
m_j [ .omega(.mathbf{q}.nu) ]^2 e_.alpha(j, .mathbf{q}.nu).](http://phonopy.sourceforge.net/_images/math/183f5bf949c284230b2917e7e04d1457c0069145.png)
where the eigenvector of the band index 
at is obtained by the diagonalization of 
:
![.sum_{j .alpha j' .beta}e_.alpha(j',.mathbf{q}.nu)^* D_{.alpha.beta}(jj',.mathbf{q})
e_.beta(j',.mathbf{q}.nu') = [.omega(.mathbf{q}.nu)]^2 .delta_{.nu.nu'}.](http://phonopy.sourceforge.net/_images/math/cb895e7e39b1a3b539163831c7edbf2912d3e549.png)
The atomic displacements 
are given as
![u_.alpha(jl,t) = .left(.frac{.hbar}{2Nm_j}.right)^{.frac{1}{2}}
.sum_{.mathbf{q},.nu}.left[.omega(.mathbf{q}.nu).right]^{-.frac{1}{2}}
.left[.hat{a}(.mathbf{q}.nu).exp(-i.omega(.mathbf{q}.nu)t)+
.hat{a}^.dagger(.mathbf{-q}.nu).exp({i.omega(.mathbf{q}.nu)}t).right]
.exp({i.mathbf{q}.cdot.mathbf{r}(jl)})
e_.alpha(j,.mathbf{q}.nu),](http://phonopy.sourceforge.net/_images/math/f2ff979b7528023acaf17b604b73e84224954e85.png)
where 
and 
are the creation and annihilation operators of phonon, 
is the reduced Planck constant, and 
is the time.
类似于软膜处理方法备忘-OUTCAR中部分内容摘录与解析
Dimension of arrays:
k-points NKPTS = 6 k-points in BZ NKDIM = 6 number of bands NBANDS= 248
number of dos NEDOS = 301 number of ions NIONS = 72
non local maximal LDIM = 8 non local SUM 2l+1 LMDIM = 32
total plane-waves NPLWV = 307200
max r-space proj IRMAX = 1 max aug-charges IRDMAX= 288919
dimension x,y,z NGX = 80 NGY = 80 NGZ = 48
dimension x,y,z NGXF= 160 NGYF= 160 NGZF= 96
support grid NGXF= 320 NGYF= 320 NGZF= 192
ions per type = 54 18
NGX,Y,Z is equivalent to a cutoff of 11.60, 11.60, 11.63 a.u.
NGXF,Y,Z is equivalent to a cutoff of 23.20, 23.20, 23.26 a.u.
I would recommend the setting:
dimension x,y,z NGX = 75 NGY = 75 NGZ = 45 【采用推荐值怎么样？ 是不是不设置该参数，则程序自动采用推荐值？】
***************
energy-cutoff : 400.00
volume of cell : 781.45
direct lattice vectors reciprocal lattice vectors
11.467499972 0.000000000 0.000000000 0.087202965 0.050346655 0.000000000
-5.733749986 9.931146293 0.000000000 0.000000000 0.100693311 0.000000000
0.000000000 0.000000000 6.861700058 0.000000000 0.000000000 0.145736478
length of vectors
11.467499972 11.467499971 6.861700058 0.100693311 0.100693311 0.145736478
k-points in units of 2pi/SCALE and weight: K-Points
0.00000000 0.00000000 0.00000000 0.037
0.02906766 0.01678222 0.00000000 0.222
0.02906766 0.05034666 0.00000000 0.074
0.00000000 0.00000000 0.04857883 0.074
0.02906766 0.01678222 0.04857883 0.444
0.02906766 0.05034666 0.04857883 0.148
k-points in reciprocal lattice and weights: K-Points
0.00000000 0.00000000 0.00000000 0.037
0.33333333 -0.00000000 0.00000000 0.222
0.33333333 0.33333333 0.00000000 0.074
0.00000000 0.00000000 0.33333333 0.074
0.33333333 -0.00000000 0.33333333 0.444
0.33333333 0.33333333 0.33333333 0.148
position of ions in fractional coordinates (direct lattice)
0.00000000 0.00000000 0.75000000
0.33333333 0.00000000 0.75000000
0.66666667 0.00000000 0.75000000
......
0.22222221 0.77777777 0.25000000
0.55555554 0.77777777 0.25000000
0.88888888 0.77777777 0.25000000
position of ions in cartesian coordinates (Angst):
0.00000000 0.00000000 5.14627504
3.82249999 0.00000000 5.14627504
7.64499998 0.00000000 5.14627504
........
-1.91125011 7.72422483 1.71542501
1.91124988 7.72422483 1.71542501
5.73374987 7.72422483 1.71542501
*********
215 f/i= 7.469225 THz 46.930526 2PiTHz 249.146526 cm-1 30.890243 meV
X Y Z 【对应笛卡尔坐标】 dx dy dz
0.000000 0.000000 5.146275 0.000000 -0.000000 -0.184253
3.822500 0.000000 5.146275 -0.000000 -0.000000 0.297550
7.645000 0.000000 5.146275 0.000000 -0.000000 -0.113297
..........
-1.911250 7.724225 1.715425 -0.000000 -0.000000 0.000000
1.911250 7.724225 1.715425 -0.000000 0.000000 0.000000
5.733750 7.724225 1.715425 -0.000000 0.000000 -0.000000
216 f/i= 7.469225 THz 46.930526 2PiTHz 249.146526 cm-1 30.890243 meV
X Y Z dx dy dz
0.000000 0.000000 5.146275 -0.000000 0.000000 -0.237202
3.822500 0.000000 5.146275 0.000000 0.000000 -0.040966
7.645000 0.000000 5.146275 0.000000 0.000000 0.278169
..........
-1.911250 7.724225 1.715425 0.000000 0.000000 -0.000000
1.911250 7.724225 1.715425 0.000000 -0.000000 0.000000
5.733750 7.724225 1.715425 -0.000000 -0.000000 0.000000
ELASTIC MODULI CONTR FROM IONIC RELAXATION (kBar)【与声子谱虚频有何对应关系？】
Direction XX YY ZZ XY YZ ZX
--------------------------------------------------------------------------------
XX -161.3730 152.5077 17.6648 -0.0000 0.0000 -0.0000
YY 152.5077 -161.3730 17.6648 -0.0000 -0.0000 0.0000
ZZ 17.6648 17.6648 -70.3967 0.0000 0.0000 -0.0000
XY -0.0000 -0.0000 0.0000 -156.9404 0.0000 0.0000
YZ 0.0000 -0.0000 0.0000 0.0000 -5.0821 -0.0000
ZX -0.0000 0.0000 -0.0000 0.0000 -0.0000 -5.0822
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------------------------------
LATTYP: Found a hexagonal cell.
ALAT = 11.4674999715
C/A-ratio = 0.5983605908
Lattice vectors:
A1 = ( 11.4674999715, 0.0000000000, 0.0000000000)
A2 = ( -5.7337499856, 9.9311462931, 0.0000000000)
A3 = ( 0.0000000000, 0.0000000000, 6.8617000580)
Analysis of symmetry for initial positions (statically):
=====================================================================
Subroutine PRICEL returns following result:
LATTYP: Found a hexagonal cell.
ALAT = 3.8224999905
C/A-ratio = 1.7950817724
Lattice vectors:
A1 = ( 3.8224999905, 0.0000000000, 0.0000000000)
A2 = ( -1.9112499952, 3.3103820977, 0.0000000000)
A3 = ( 0.0000000000, 0.0000000000, 6.8617000580)
9 primitive cells build up your supercell.
Routine SETGRP: Setting up the symmetry group for a
hexagonal supercell.
Subroutine GETGRP returns: Found 24 space group operations
(whereof 12 operations were pure point group operations)
out of a pool of 24 trial point group operations.
The static configuration has the point symmetry D_3d.
The point group associated with its full space group is D_6h.
Analysis of symmetry for dynamics (positions and initial velocities):
=====================================================================
Subroutine PRICEL returns following result:
LATTYP: Found a hexagonal cell.
ALAT = 3.8224999905
C/A-ratio = 1.7950817724
Lattice vectors: 【对应扩胞前的POSCAR】
A1 = ( 3.8224999905, 0.0000000000, 0.0000000000)
A2 = ( -1.9112499952, 3.3103820977, 0.0000000000)
A3 = ( 0.0000000000, 0.0000000000, 6.8617000580)
9 primitive cells build up your supercell.
Routine SETGRP: Setting up the symmetry group for a
hexagonal supercell.
Subroutine GETGRP returns: Found 24 space group operations
(whereof 12 operations were pure point group operations)
out of a pool of 24 trial point group operations.
The dynamic configuration has the point symmetry D_3d.
The point group associated with its full space group is D_6h.
Analysis of magnetic symmetry:
=====================================================================
Subroutine MAGSYM returns: Found 24 space group operations
(whereof 12 operations were pure point group operations)
out of a pool of 24 trial space group operations
(whereof 12 operations were pure point group operations)
and found also 9 'primitive' translations
The magnetic configuration has the point symmetry D_3d.
The point group associated with its full space group is D_6h.
初始扩胞前 POSCAR
1.00000000000000
3.8224999904999999 0.0000000000000000 0.0000000000000000
-1.9112499951999999 3.3103820976999998 0.0000000000000000
0.0000000000000000 0.0000000000000000 6.8617000580000003
6 2
Direct
0.0000000000000000 0.0000000000000000 0.7500000000000000
0.0000000000000000 0.0000000000000000 0.2500000000000000
0.3333333429999996 0.6666666870000029 0.0946099980000028
0.6666666269999979 0.3333333129999971 0.9053900239999990
0.6666666269999979 0.3333333129999971 0.5946099760000010
0.3333333429999996 0.6666666870000029 0.4053899939999965
0.3333333429999996 0.6666666870000029 0.7500000000000000
0.6666666269999979 0.3333333129999971 0.2500000000000000
https://blog.sciencenet.cn/blog-567091-835071.html
上一篇：结构优化与插值：SED命令使用等
下一篇：Convex hull 图再谈
