博文

导师的信：How to follow the imaginary phonon vectors?

已有 11331 次阅读 2014-5-17 17:49 |系统分类:科研笔记

关注：

1) How to follow the imaginary phonon vectors

题记：

Let us discuss the GdH₃ structure for ScH₃. The symmetric P6₃/mmcstructure, shown in Fig. 8(a), is not dynamically stable from 1 atm to 25 GPa (see Fig.S13(a)). We have not yet been able to follow the imaginary phonon to a more stable minimum[rh1] . Recall that the experimentally observed lower symmetry structure isf or a nonstoichiometric ScH_2.9phase. The fractional occupation is consistent with motion away from a moresymmetrical structure, but the details remain to be worked out.

[rh1]Have you ever tried to follow the imaginary phonon vectors? If not,write to Wojciech Grochala, a former group member, and ask him how one does that. He’s friendly…

问答：

Dear Prof. Wojciech Grochala,
First of all, I hope everything goes very well with you.
As a research assistant under the direction of Prof. Roald Hoffmann, I am impressed by your outstanding work about hydrides especially after the introduction by Roald. Recently, I met a problem about scandium hydrides when I did phonon calculations:
"I find the GdH₃ structure for ScH₃ is not dynamically stable from 1 atm to 25 GPa, as shown in the attached file"
Could you kindly tell me how to follow the corresponding eigenvector ( or the imaginary phonon vectors) to a more stable minimum? I also tried to get some hints from the following link, but I failed:
http://phonopy.sourceforge.net/setting-tags.html#create-modulated-structure
You can find the 'cif' file and the result of phonon calculations of ScH₃ in the attachment.
Thank you very much for your time and looking forward to your kindly suggestions.
Best wishes,
Sincerely yours,

答：

Dear Xiaoqiu,
yes, I enjoy H-rich systems a lot. we've recently managed to generate a hydride hosting the strongest AFM intteractions to date, but more work is needed to fully understand this system.

I have the following comments to your problem:
1) it is difficult to calculate any phonon spectrum properly, one must rigorously optimize the structure for energy and than for forces, and only then run the phonon calc. Otherwise artifacts such as pseudo-img modes, may be seen
2) the img parts away from gamma (as you've computed) are not so often but they indeed happen. Formally they imply necessity of using a very large supercell to distort the structure, but sometimes, fortunately, this distorted structure may be symmetrized to a much smaller unit cell of a lower symmetry.
3) recipe for following the img phonon mode sound simple.

for each atom one must obtain a distorted structure according to:

x' = x + 0.02*x(img)

where x is one of 3 coordinates (x,y,z), x' is the new sought coordinate, 0.02 is arbitrary weight constant (realistically could be anything from 0.01 up to 0.05, or 0.1 when the distortion is particularly soft) while x(img) is the amplitude of the normal mode in question.

I never do it manually, phonon @ medea does it for me. I know it is possible to do it while having only phonon without medea but I do not know how.
I am not sure if other programs do it.
sorry that I can't help more

Xiaoqiu

2014-5-22

问：

Dear Prof.Wojciech Grohala,
Thank you very much for your reminder.
Yes,the pseudo-ing phonons appeared near K (-0.333 0.667 0.000) in my case.
Could you kindly give me some hints:
1) How large supercells I should use (along K)?
2) What does "the amplitude of the normal mode in question" mean in the equation of x' = x + 0.02*x(img)? Where I can find this equation?
Best wishes,
Xiaoqiu

答：From Prof.Wojciech Grochala

hi
This means that the smallest supercell to use is 331 (since 3, 1.5, 1 is a physical), but it may reduce to some smaller representation
I do not know, regretfully, where your program (Phonopy) holds information about atomic displacements associated with each normal mode at each k point...
If you have ever used Gaussian, you should know what I am talking about
Gausiian gives atomic (say Cartesian) coordinations for each atom:

B 0.00 0.125 0.234

and the same for B's motion within a given normal mode, say:

B 0.12 0.14 -0.12

the equation I have sent you means simply to get new position of B as:

B (0.00+0.02*0.12) (0.125+0.02*0.14) (0.234-0.02*0.12)

I really do not know Phonopy and I cannot help. you must read the manual.
best,
w

答：from Prasad

Dear Xiaoqiu,

Nice to hear from you. In general when we find an imaginary frequency we would generate a supercell along the imaginary frequency mode direction. Please visualize the imaginary mode and create a supercell (typically we can double the cell along that direction).

This supercell should switch off the imaginary mode, some times the geometry optimzation of the supercell may generate a new structure with a different space group, but may be close to the initial geometry.

I will go through your attachments this week-end, will write you more. Sure, we will see each other, soon.

Best regards,

Prasad

网络摘录：

Re: [Phonopy-users] Imaginary eigenvector

From: Atsushi Togo <atz.togo@gm...> - 2012-09-07 14:06:53

Hi,

You can use the MODULATION tag.
http://phonopy.sourceforge.net/setting-tags.html#create-modulated-structure

Togo

Good afternoon
>
> Having a phonon mode with imaginary frequency, I want to follow the
> corresponding eigenvector to lower the total energy of the structure. The
> eigenvector is given with both real and imaginary parts like the following
>
> - # atom 9
> - [ -0.15557262927497, -0.20071249623661 ]
> - [ -0.00000000000000, 0.00000000000000 ]
> - [ 0.00000000000000, 0.00000000000001 ]
>
> The left column gives the real and the right column give the imaginary parts
> of three vector elements of the eigenvector. I do not know how to "follow"
> the imaginary part of the vector. Can anyone give me a hint?
>
> Thank you
> Ryan
>
> ======================
> With my best regards
> Ryan Nguyen
> Faculty of Physics
> Hanoi University of Science
> 334 Nguyen Trai, Thanh Xuan

Create modulated structure¶

MODULATION¶

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],$

where $A$ is the amplitude, $.phi$ is the phase, and $m_j$ is the mass of the j-th atom, $.mathbf{q}$ is the q-point specified, $.mathbf{r}_{jl}$ is the position of the j-th atom and in the l-th unit cell, and $.mathbf{e}_j$ 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.

Usage¶

The first three values correspond to the supercell dimension. 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 are the band index from the bottom with ascending order, amplitude, and 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, $.left|.mathbf{e}_j.right|/.sqrt{m_j}$ , of element of 3N dimensional modulation vector to be real. The complex modulation vector is shown in modulation.yaml.

MODULATION = 3 3 1, 1/3 1/3 0 1 2, 1/3 1/3 2 3.5

MODULATION = 3 3 1, 1/3 1/3 0 1 2, 1/3 0 0 2 2

MODULATION = 3 3 1, 1/3 1/3 0 1 1 0, 1/3 1/3 0 1 1 90

网络摘录：

http://emuch.net/html/201109/3611259.html

大家好，用siesta算声子振动谱，出现虚频，如附件图所示，怎么消除虚频呢，我看到很多人讨论过，说提高收敛精度，还有就是用虚频对应的简正坐标加到体系的坐标上去，重新优化结构，再计算声子。我精度设到10-8，已经很高，但结果还是不好，所以我想试试用虚频对应简正坐标加到体系坐标上去重新优化试试，但是这个具体操作过程我不是很了解，一个虚频的振动模式不是体系所有的原子的振动结果么？虚频的简正坐标是根据对应gamma点的虚频的频率值，解附件中的常微分方程得到的么，得到的简正坐标又加到体系的哪个原子上呢，是x,y,z的哪一个分量上呢，很多东西不是很明白，还望大家不腻赐教，谢谢!

QQQ:有虚频说明你的结构动力学不稳定，原则上按照虚频的模式的原子振动方式移动原子（冷冻声子方法）可以找到动力学稳定的结构，但是你这个结构虚的太离谱，基本类似能带了，无法用冷冻声子计算

AAA:

谢谢您的回复，DFT优化是在初始模型的基础上小范围移动原子找到局域内稳定结构，可能并不是体系的最稳定结构，所以我就想按照虚频的原子振动方式移动原子，改变初始构型，优化找到动力学最稳定的结构，所以用这种方法消除虚频应该还是合理的吧。但是也就是怎样去移动原子这点我比较模糊，虚频是个别原子造成的么？siesta虚频的振动方向可以可视化么，希望您能够讲述移动原子的细节。非常感谢您，我算的是纳米管的周期性结构。

yinwei_li (站内联系TA)

3楼: Originally posted by xhzha at 2011-09-20 09:52:38:
谢谢您的回复，DFT优化是在初始模型的基础上小范围移动原子找到局域内稳定结构，可能并不是体系的最稳定结构，所以我就想按照虚频的原子振动方式移动原子，改变初始构型，优化找到动力学最稳定的结构，所以用这种 ...

通过原子移动来找稳定结构的方法是有一定限制的，如果你的声子谱只是在布里渊区某个特殊点有虚频，找到此点对应的本征矢（本征矢可表征原子如何振动），沿着本征矢逐步移动原子并计算每移动一步体系的能量的变化，你会发现能量的变化会出现一个低谷，那么这个能量低谷的位置就对应新的稳定结构。计算这个稳定结构的声子谱应该是没有虚频的。但是你给出的声子谱在整个布里渊区都有虚频，用这种方法可能就不适应了，因为你没办法同时按照所有布里渊区点的本征矢移动原子。
或许大家有别的方法，这个我就不太清楚了！

Tutorial using VASP as calculator¶

Pre-process¶

The input stureture of POSCAR is supposed to be this.

In the pre-process, supercell structures with (or without) displacements are created from a unit cell fully consiering crystal symmetry.

To obtain supercells ( $2.times 2.times 3$ ) with displacements, run phonopy:

phonopy-d--dim="2 2 3"

You should find the files, SPOSCAR, disp.yaml, and POSCAR-{number} as follows:

% lsdisp.yaml POSCAR POSCAR-001 POSCAR-002 POSCAR-003 SPOSCAR

SPOSCAR is the perfect supercell structure, disp.yaml contains the information on displacements, and POSCAR-{number} are the supercells with atomic displacements. POSCAR-{number} corresponds to the different atomic displacements written in disp.yaml.

Calculation of sets of forces¶

Force constants are calculated using the structure files POSCAR-{number} (from forces on atoms) or using the SPOSCAR file (direct calculation of force constants) by your favorite calculator. See the details.

In the case of VASP, the calculations for the finite displacement method can be proceeded just using the POSCAR-{number} files as POSCAR of VASP calculations. An example of the INCAR is as follows:

   PREC = Accurate IBRION = -1  ENCUT = 500  EDIFF = 1.0e-08 ISMEAR = 0; SIGMA = 0.01  IALGO = 38  LREAL = .FALSE.ADDGRID = .TRUE.  LWAVE = .FALSE. LCHARG = .FALSE.

Be careful not to relax the structures. Then create FORCE_SETS file using VASP interface:

% phonopy -f disp-001/vasprun.xml disp-002/vasprun.xml disp-003/vasprun.xml

% phonopy -f disp-{001..003}/vasprun.xml

If you want to calculate force constants by VASP-DFPT directory, see VASP-DFPT & phonopy calculation.

Post-process¶

In the post-process,

Force constants are calculated from the sets of forces
A part of dynamical matrix is built from the force constants
Phonon frequencies and eigenvectors 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 ODIM = 2 2 3MP = 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 convergence with respect to the mesh numbers. Thermal properties can be plotted by:

% phonopy -t -p mesh.conf

Projected DOS is calculated by the following setting file named, e.g., pdos.conf:

ATOM_NAME = Si ODIM = 2 2 3MP = 8 8 8PDOS = 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 ODIM =  2 2 3BAND = 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.

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.
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.
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.

VASP can calculate force constants in real space using DFPT. The procedure to calculate phonon properties may be as follows:

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 4Direct 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
Prepare a perfect supercell structure from POSCAR-unitcell, e.g.,
% phonopy -d --dim="2 2 2" -c POSCAR-unitcell
Rename SPOSCAR created in (2) to POSCAR (POSCAR-{number} and disp.yaml files will never be used.)
% mv SPOSCAR POSCAR
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 LREAL = .FALSE.ADDGRID = .TRUE. LWAVE = .FALSE. LCHARG = .FALSE.
After finishing the VASP calculation, confirm vasprun.xml contains hessian elements, and then create FORCE_CONSTANTS:
% phonopy --fc vasprun.xml
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.1Band structure modeSettings: 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]