Unexpected phonon frequency shift at Γ when including Quadrupole correction in EPW

[Hust_wty]

I am currently performing electron–phonon interaction calculations for wurtzite GaN using EPW, and I have encountered an issue related to the Quadrupole correction in the phonon self-energy calculation.

As shown in Figure 1

1.png (95.33 KiB) Viewed 2779 times

， when the Quadrupole term is not included, the phonon self-energy calculation gives a reasonable result:
at the Γ (0 0 0) point, the phonon frequency ω = 0 meV, and both lambda and gamma___ are 0, which is physically correct.

However, as shown in Figure 2,

2.png (98.03 KiB) Viewed 2779 times

when the Quadrupole correction is included, the phonon frequency at Γ starts from 18.3356 meV, which is clearly unphysical.
The Quadrupole tensor itself seems correct — it matches the literature and previously published results (see attached tensor tables).
Could this frequency shift be related to an incorrect phonon renormalization or normalization after including the Quadrupole correction?
the follow is epw input
--
&inputepw
prefix = 'gan'
amass(1) = 69.723
amass(2) = 14.0067
outdir = './output'

iverbosity = 3

elph = .true.
!epbwrite = .true.
!epbread = .false.
epwwrite = .fasle.
epwread = .true.
etf_mem = 1
vme = 'wannier'
!selecqread = .true.


lpolar = .true.
lifc = .false.
asr_typ = 'simple'
use_ws = .true.
lphase = .true.

nbndsub = 14
bands_skipped = 'exclude_bands = 1-12'

wannierize = .false.
num_iter = 1500
iprint = 3
dis_win_max = 24.0660
dis_win_min = 1
dis_froz_max= 11.03265
proj(1) = 'Ga:sp3'
proj(2) = 'N:p'
wdata(1) = 'bands_plot = .true.'
wdata(2) = 'begin kpoint_path'
wdata(3) = 'G 0.0000 0.0000 0.0000 M 0.5000 0.0000 0.0000'
wdata(4) = 'M 0.5000 0.0000 0.0000 K 0.3333 0.3333 0.0000'
wdata(5) = 'K 0.3333 0.3333 0.0000 G 0.0000 0.0000 0.0000'
wdata(6) = 'G 0.0000 0.0000 0.0000 A 0.0000 0.0000 0.5000'
wdata(7) = 'A 0.0000 0.0000 0.5000 L 0.5000 0.0000 0.5000'
wdata(8) = 'L 0.5000 0.0000 0.5000 H 0.3333 0.3333 0.5000'
wdata(9) = 'end kpoint_path'
wdata(10) = 'bands_plot_format = gnuplot'
wdata(11) = 'guiding_centres = .true.'
wdata(12) = 'dis_num_iter = 1000'
wdata(13) = 'num_print_cycles = 10'
wdata(14) = 'dis_mix_ratio = 1.0'
wdata(15) = 'conv_tol = 1E-8'
wdata(16) = 'conv_window = 4'

delta_approx= .false.
elecselfen = .false.
phonselfen = .true.
a2f = .false.

fsthick = 0.2 ! eV
temps = 300 ! K
degaussw = 0.02 ! eV

scissor = 1.5579

dvscf_dir = './save/'

efermi_read = .true
fermi_energy= 10.0328
!band_plot = .true.

!filkf = './eband.check'
filqf = './BTE.qpoints_36'

nq1 = 6
nq2 = 6
nq3 = 6
!nqf1 = 12
!nqf2 = 12
!nqf3 = 12

nk1 = 12
nk2 = 12
nk3 = 12
nkf1 = 60
nkf2 = 60
nkf3 = 60
!mp_mesh_k = .true.
/

28 cartesian
0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00
0.000000000000000E+00 0.000000000000000E+00 0.102310586375119E+00
0.000000000000000E+00 0.000000000000000E+00 0.204621172750238E+00
0.000000000000000E+00 0.000000000000000E+00 -0.306931759125357E+00
0.000000000000000E+00 0.192450089729862E+00 0.000000000000000E+00
0.000000000000000E+00 0.192450089729862E+00 0.102310586375119E+00
0.000000000000000E+00 0.192450089729862E+00 0.204621172750238E+00
0.000000000000000E+00 0.192450089729862E+00 -0.306931759125357E+00
0.000000000000000E+00 0.384900179459723E+00 0.000000000000000E+00
0.000000000000000E+00 0.384900179459723E+00 0.102310586375119E+00
0.000000000000000E+00 0.384900179459723E+00 0.204621172750238E+00
0.000000000000000E+00 0.384900179459723E+00 -0.306931759125357E+00
0.000000000000000E+00 -0.577350269189585E+00 0.000000000000000E+00
0.000000000000000E+00 -0.577350269189585E+00 0.102310586375119E+00
0.000000000000000E+00 -0.577350269189585E+00 0.204621172750238E+00
0.000000000000000E+00 -0.577350269189585E+00 -0.306931759125357E+00
0.166666666666667E+00 0.288675134594792E+00 0.000000000000000E+00
0.166666666666667E+00 0.288675134594792E+00 0.102310586375119E+00
0.166666666666667E+00 0.288675134594792E+00 0.204621172750238E+00
0.166666666666667E+00 0.288675134594792E+00 -0.306931759125357E+00
0.166666666666667E+00 0.481125224324654E+00 0.000000000000000E+00
0.166666666666667E+00 0.481125224324654E+00 0.102310586375119E+00
0.166666666666667E+00 0.481125224324654E+00 0.204621172750238E+00
0.166666666666667E+00 0.481125224324654E+00 -0.306931759125357E+00
0.333333333333333E+00 0.577350269189585E+00 0.000000000000000E+00
0.333333333333333E+00 0.577350269189585E+00 0.102310586375119E+00
0.333333333333333E+00 0.577350269189585E+00 0.204621172750238E+00
0.333333333333333E+00 0.577350269189585E+00 -0.306931759125357E+00

To validate my setup, I also tested the official SiC Quadrupole example provided with EPW. the SiC calculation without Quadrupole gives ω = 0 meV at Γ (correct), and the computed scattering rates are reasonable.
When the Quadrupole correction is added, however, ω becomes –10.9810 meV at Γ, again unphysical.

I have checked Gaussian broadening, and all relevant parameters, but the results remain the same.

I have also attached a comparison of my calculated GaN electron–phonon linewidths with and without Quadrupole corrections.

10.png (82.38 KiB) Viewed 2779 times

Could you please advise whether this issue could be due to a normalization problem in the phonon renormalization procedure with Quadrupoles, or if there is another known reason for this behavior?

[Nina Girotto]

Hello,

Correct me if I am wrong, but it seems that the issue is related with the interpolated phonon frequencies, and not with the phonon self-energy calculation.
Did you run the quadrupole EPW calculation as the initial calculation, with flags:

epwwrite = .true.
wannierize = .true.

or did you run it only as a restart calculation?
It is important to include the quadrupole contributions also in the interpolation procedure.

Best,

Nina.

[Hust_wty]

Hello,

Correct me if I am wrong, but it seems that the issue is related with the interpolated phonon frequencies, and not with the phonon self-energy calculation.
Did you run the quadrupole EPW calculation as the initial calculation, with flags:

epwwrite = .true.
wannierize = .true.

or did you run it only as a restart calculation?
It is important to include the quadrupole contributions also in the interpolation procedure.

Best,

Nina.

Hello,

Correct me if I am wrong, but it seems that the issue is related with the interpolated phonon frequencies, and not with the phonon self-energy calculation.
Did you run the quadrupole EPW calculation as the initial calculation, with flags:

epwwrite = .true.
wannierize = .true.

or did you run it only as a restart calculation?
It is important to include the quadrupole contributions also in the interpolation procedure.

Best,

Nina.

Dear Nina,

Thanks for your reply. My run was a restart (I didn’t redo the interpolation with epwwrite=.true. and wannierize=.true.).

I computed the quadrupole tensor for w-GaN (3m) with ABINIT (per the EPW tutorial); the values agree with PRL 125, 136602 (2020) in figure 1.

1.png (104.22 KiB) Viewed 423 times

After placing them in quadrupole.fmt in figure 2,

2.png (145.83 KiB) Viewed 423 times

EPW fails during interpolation if Qxy ≠ 0 (“Unexpected phonon frequency”); setting Qxy = 0 avoids the error.

Two questions:

Basis mapping: The paper/ABINIT use crystal-basis components (Q11, Q22, …) while quadrupole.fmt uses Cartesian (Qxx, Qyy, Qxy, …). Should I rotate via 𝑄cart=𝑅 𝑄crys 𝑅T Qcart=R Qcrys RT
using the QE lattice? Could a wrong mapping be turning a symmetry-allowed crystal component into a spurious Cartesian Qxy, triggering the failure?

Symmetry: For w-GaN (3m), should Cartesian Qxy be symmetry-forbidden (≈0), and does EPW expect the input to already satisfy that?

I’ll rerun from scratch with quadrupole included in the initial interpolation.

[Shashi]

Hi Hust,

When you want to use Quadrupole corrections, you need to keep the .fmt files from scratch. But if you read quadrupole.fmt file only from the restart step, there might be an issue. As in the Bloch-Wannier electron-phonon matrix uses dipole correction, while in the Wannier-Bloch step, you are using quadrupole correction. This might be a cause of failure. Please start from scratch and let us know if the issue still persists.

The Quadrupole file uses Cartesian format. Symmetry is not not issue with the Quadrupole calculation. It works with symmetry.

regards,
Shashi