博文

电声耦合常数计算之：手册学习

已有 9259 次阅读 2013-11-9 04:35 |个人分类:声子谱计算|系统分类:科研笔记

关注：有时感觉发文章才是硬道理哈

http://www.quantum-espresso.org/wp-content/uploads/Doc/ph_user_guide/node1.html

The phonon code ph.x calculates normal modes at a given q-vector, starting from data files produced by pw.x with a simple SCF calculation.

NOTE: the alternative procedure in which a band-structure calculation with calculation='phonon'was performed as an intermediate step is no longer implemented since version 4.1. It is also no longer needed to specify lnscf=.true. for 0 .

The output data files appear in the directory specified by the variable outdir, with names specified by the variable prefix.

After the output file(s) has been produced (do not remove any of the files, unless you know which are used and which are not), you can run ph.x.

cat > nbh2.elph.in << EOF

Electron-phonon coefficients for NbH4

&inputph

tr2_ph=1.0d-12,

prefix='nbh2',

amass(1)=92.90638,

amass(2)=1.008,

outdir='/scratch/gg375/espresso/',

fildyn='nbh2.dyn',

fildvscf='nbh2.dv'

elph=.true.,

trans=.true.,

ldisp = .true.,

nq1=8,nq2=8,nq3=8

EOF

mpirun -np $NPROCS ph.x < nbh2.elph.in > nbh2.elph.out

The first input line of ph.x is a job identifier.

At the second line the namelist &INPUTPH starts.

The meaning of the variables in the namelist (most of them having a default value) is described in file Doc/INPUT_PH.*.

Variables outdirand prefix must be the same as in the input data of pw.x.

Presently you can specify amass(i) (a real variable) the atomic mass of atomic type i or you can use the default one deduced from the periodic table on the basis of the element name. If amass(i) is not given as input of ph.x, the one given as input in pw.x is used. When this is 0 the default one is used.

After the namelist you must specify the q-vector of the phonon mode, in Cartesian coordinates and in units of 2/a .

Notice that the dynamical matrix calculated by ph.x at = 0 does not contain the non-analytic term occurring in polar materials, i.e. there is no LO-TO splitting in insulators. 【不理解】

Moreover no Acoustic Sum Rule (ASR) is applied. In order to have the complete dynamical matrix at = 0 including the non-analytic terms, you need to calculate effective charges by specifying option epsil=.true. to ph.x. This is however not possible (because not physical!) for metals (i.e. any system subject to a broadening).

At = 0 , use program dynmat.x to calculate the correct LO-TO splitting, IR cross sections, and to impose various forms of ASR. If ph.x was instructed to calculate Raman coefficients, dynmat.x will also calculate Raman cross sections for a typical experimental setup. Input documentation in the header of PHonon/PH/dynmat.f90.

See Example 01 for a simple phonon calculations in Si, Example 06 for fully-relativistic calculations (LDA) on Pt, Example 07 for fully-relativistic GGA calculations.

4.2 Calculation of interatomic force constants in real space

First, dynamical matrices are calculated and saved for a suitable uniform grid of q-vectors (only those in the Irreducible Brillouin Zone of the crystal are needed). Although this can be done one q-vector at the time, a simpler procedure is to specify variable ldisp=.true. and to set variables nq1, nq2, nq3 to some suitable Monkhorst-Pack grid, that will be automatically generated, centered at = 0 .

Second, code q2r.x reads the dynamical matrices produced in the preceding step and Fourier-transform them, writing a file of Interatomic Force Constants in real space, up to a distance that depends on the size of the grid of q-vectors. Input documentation in the header of PHonon/PH/q2r.f90.

Program matdyn.x may be used to produce phonon modes and frequencies at any q using the Interatomic Force Constants file as input. Input documentation in the header of PHonon/PH/matdyn.f90.

See Example 02 for a complete calculation of phonon dispersions in AlAs.

4.3 Calculation of electron-phonon interaction coefficients

Since v.5.0, there are two ways of calculating electron-phonon coefficients, distinguished according to the value of variable electron_phonon.

The following holds for the case electron_phonon= 'interpolated' (see also Example 03).

The calculation of electron-phonon coefficients in metals is made difficult by the slow convergence of the sum at the Fermi energy.

It is convenient to use a coarse k-point grid to calculate phonons on a suitable wavevector grid;

a dense k-point grid to calculate the sum at the Fermi energy. The calculation proceeds in this way:

a scf calculation for the dense -point grid (or a scf calculation followed by a non-scf one on the dense -point grid); specify option la2f=.true. to pw.x in order to save a file with the eigenvalues on the dense k-point grid. The latter MUST contain all and + grid points used in the subsequent electron-phonon calculation. All grids MUST be unshifted, i.e. include = 0 .
a normal scf + phonon dispersion calculation on the coarse k-point grid, specifying option electron_phonon='interpolated', and the file name where the self-consistent first-order variation of the potential is to be stored: variablefildvscf). The electron-phonon coefficients are calculated using several values of Gaussian broadening (see PHonon/PH/elphon.f90) because this quickly shows whether results are converged or not with respect to the k-point grid and Gaussian broadening.
Finally, you can use matdyn.x and lambda.x (input documentation in the header of PHonon/PH/lambda.f90) to get the $.alpha^{2}_{}$ F() function, the electron-phonon coefficient , and an estimate of the critical temperature T_c .

See the appendix for the relevant formulae. Important notice: the q 0 limit of the contribution to the electron-phonon coefficient diverges for optical modes! please be very careful, consult the relevant literature.

6 Troubleshooting

6.0.0.1 ph.x stops with error reading fileThe data file produced by pw.x is bad or incomplete or produced by an incompatible version of the code. In parallel execution: if you did not set wf_collect=.true., the number of processors and pools for the phonon run should be the same as for the self-consistent run; all files must be visible to all processors.

6.0.0.2 ph.x mumbles something like cannot recover or error reading recover fileYou have a bad restart file from a preceding failed execution. Remove all files recover* in outdir.

6.0.0.3 ph.x says occupation numbers probably wrong and continuesYou have a metallic or spin-polarized system but occupations are not set to `smearing'.

6.0.0.4 ph.x does not yield acoustic modes with zero frequency at = 0This may not be an error: the Acoustic Sum Rule (ASR) is never exactly verified, because the system is never exactly translationally invariant as it should be. The calculated frequency of the acoustic mode is typically less than 10 cm^-1 , but in some cases it may be much higher, up to 100 cm^-1 . The ultimate test is to diagonalize the dynamical matrix with program dynmat.x, imposing the ASR. If you obtain an acoustic mode with a much smaller (let us say< 1cm^-1 ) with all other modes virtually unchanged, you can trust your results.

``The problem is [...] in the fact that the XC energy is computed in real space on a discrete grid and hence the total energy is invariant (...) only for translation in the FFT grid. Increasing the charge density cutoff increases the grid density thus making the integral more exact thus reducing the problem, unfortunately rather slowly...This problem is usually more severe for GGA than with LDA because the GGA functionals have functional forms that vary more strongly with the position; particularly so for isolated molecules or system with significant portions of ``vacuum'' because in the exponential tail of the charge density a) the finite cutoff (hence there is an effect due to cutoff) induces oscillations in rho and b) the reduced gradient is diverging.''(info by Stefano de Gironcoli, June 2008)

6.0.0.5 ph.x yields really lousy phonons, with bad or negative frequencies or wrong symmetries or gross ASR violationsPossible reasons

if this happens only for acoustic modes at = 0 that should have = 0 : Acoustic Sum Rule violation, see the item before this one.
wrong data file read.
wrong atomic masses given in input will yield wrong frequencies (but the content of file fildyn should be valid, since the force constants, not the dynamical matrix, are written to file).
convergence threshold for either SCF (conv_thr) or phonon calculation (tr2_ph) too large: try to reduce them.
maybe your system does have negative or strange phonon frequencies, with the approximations you used. A negative frequency signals a mechanical instability of the chosen structure. Check that the structure is reasonable, and check the following parameters:
- The cutoff for wavefunctions, ecutwfc
- For USPP and PAW: the cutoff for the charge density, ecutrho
- The k-point grid, especially for metallic systems.

Note that ``negative'' frequencies are actually imaginary: the negative sign flags eigenvalues of the dynamical matrix for which $.omega^{2}_{}$ < 0 .

6.0.0.6 Wrong degeneracy error in star_qVerify the q-vector for which you are calculating phonons. In order to check whether a symmetry operation belongs to the small group of , the code compares and the rotated , with an acceptance tolerance of 10^-5 (set in routine PW/eqvect.f90). You may run into trouble if your q-vector differs from a high-symmetry point by an amount in that order of magnitude.

A. Appendix: Electron-phonon coefficients

The electron-phonon coefficients g are defined as

, i, j) = $.displaystyle .left(.vphantom{{.hbar.over 2M.omega_{{.bf q}.nu}}}.right.$ $.displaystyle {.hbar.over 2M.omega_{{.bf q}.nu}}$ $.displaystyle .left..vphantom{{.hbar.over 2M.omega_{{.bf q}.nu}}}.right)^{{1/2}}_{}$

$.displaystyle .psi_{{i,{.bf k}}}^{}$ | $.displaystyle {dV_{SCF}.over d {.hat u}_{{.bf q}.nu}}$ ^. $.displaystyle .hat{.epsilon}_{{{.bf q}.nu}}^{}$ | $.displaystyle .psi_{{j,{.bf k}+{.bf q}}}^{}$

(1)

The phonon linewidth $.gamma_{{{.bf q}.nu}}^{}$ is defined by

$.displaystyle .gamma_{{{.bf q}.nu}}^{}$ = 2

$.displaystyle .omega_{{{.bf q}.nu}}^{}$ $.displaystyle .sum_{{ij}}^{}$

$.displaystyle {d^3k.over .Omega_{BZ}}$ | g(

, i, j)|²

(e_{, i} - e_F)

(e_{+, j} - e_F),

(2)

while the electron-phonon coupling constant $.lambda_{{{.bf q}.nu}}^{}$ for mode at wavevector is defined as

$.displaystyle .lambda_{{{.bf q}.nu}}^{}$ = $.displaystyle {.gamma_{{.bf q}.nu} .over .pi.hbar N(e_F).omega^2_{{.bf q}.nu}}$

(3)

where N(e_F) is the DOS at the Fermi level. The spectral function is defined as

$.displaystyle .alpha^{2}_{}$ F(

) = $.displaystyle {1.over 2.pi N(e_F)}$ $.displaystyle .sum_{{{.bf q}.nu}}^{}$

(

- $.displaystyle .omega_{{{.bf q}.nu}}^{}$ ) $.displaystyle {.gamma_{{.bf q}.nu}.over.hbar.omega_{{.bf q}.nu}}$ .

(4)

The electron-phonon mass enhancement parameter can also be defined as the first reciprocal momentum of the spectral function:

= $.displaystyle .sum_{{{.bf q}.nu}}^{}$ $.displaystyle .lambda_{{{.bf q}.nu}}^{}$ = 2

$.displaystyle {.alpha^2F(.omega) .over .omega}$ d

(5)

Note that a factor M^-1/2 is hidden in the definition of normal modes as used in the code.

McMillan:

T_c = $.displaystyle {.Theta_D .over 1.45}$ exp $.displaystyle .left[.vphantom{ {-1.04(1+.lambda).over .lambda(1-0.62.mu^*)-.mu^*}}.right.$ $.displaystyle {-1.04(1+.lambda).over .lambda(1-0.62.mu^*)-.mu^*}$ $.displaystyle .left..vphantom{ {-1.04(1+.lambda).over .lambda(1-0.62.mu^*)-.mu^*}}.right]$

(6)

or (better?)

T_c = $.displaystyle {.omega_{log}.over 1.2}$ exp $.displaystyle .left[.vphantom{ {-1.04(1+.lambda).over .lambda(1-0.62.mu^*)-.mu^*}}.right.$ $.displaystyle {-1.04(1+.lambda).over .lambda(1-0.62.mu^*)-.mu^*}$ $.displaystyle .left..vphantom{ {-1.04(1+.lambda).over .lambda(1-0.62.mu^*)-.mu^*}}.right]$

(7)

where

$.displaystyle .omega_{{log}}^{}$ = exp $.displaystyle .left[.vphantom{ {2.over.lambda} .int {d.omega.over.omega} .alpha^2F(.omega) .mbox{log}.omega }.right.$