Question on |g|(qx,qy) distribution

[andreyl]

Dear all,

I am sorry in advance for the long post, but I didn't find "hide" function on the forum interface.

I faced a problem trying to obtain and analyze the q- distribution of electron phonon coupling elements |g| for given k-point and given phonon mode for monolayer graphene.

The motivation of obtaining such plots is: a) verify, that e-p matrix elements reasonably represent the electron-phonon coupling effects of the system b) to obtain the
deformation potential which can be found utilizing linear dependency of electron-phonon coupling from q (for low q's):
M=D*q

(M, can be obtained from |g|)

A number of references exist, where this task was completed for graphene, for example:
https://arxiv.org/pdf/1705.01816.pdf
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.81.121412
In both cases authors successfully obtained both the reasonable distribution of |g|(qx,qy), and deformation potential.
I've been trying to reproduce these results, so far with no success.

With print_gkk functional introduced in the last git version of EPW it seemed like a straightforward task. The averaging over E(k), E(k+q) and Omega(q) are preformed as they are provided in the implementation.
However the distributions I obtain do not even reproduce the symmetry of the system correctly. Moreover, the data seems
noisy and contain strange stripes of q-points with the same value of |g|:


q-resolved electron-phonon matrix elements |g|(qx,qy) [eV], g=K, conduction band. Top row: ZA, TA modes, bottom left LA mode, green crosses mark the high symmetry points on Brillouin zone edge.

At the same time phonon frequencies and K+q energies from the very same run look reasonable, thus I assume, that the interpolation for this quantities was preformed correctly and I don't violate the order of q's and band's when
I parse the run results:


Phonon frequencies for acoustic modes.


E(k+q) [eV].

The electron and phonon spectra in high symmetry directions are available here (this post is already to long to give them here) and look quite OK in my opinion:
https://paper.dropbox.com/doc/Graphene-elph-9QVfbGljBpzkTLb8i7E1R#:uid=365416936615388882424307&h2=Same-from-initial-QE-pw-and-ph
The link also contains all the input files I used.

My questions is:
What can lead to such behavior? I don't believe the parameters of my calculations are unreasonable, at least they are consistent with those in archived paper (authors also used EPW), maybe there is some principle error or some kind of clue that can point on, for example, bad wannierization or inadequate choice of parameters in initial calculations?

The input file for the final calculation is:

Code: Select all

&inputepw
prefix = 'c'
amass(1) = 12.0107,
outdir = '/dev/shm/andrey/epw0/'
iverbosity = 0
system_2d = .true.
eps_acustic = 1.0d0
!
elph = .true.
ep_coupling = .true.
kmaps = .true.
! Coarse blochl grid M_e-ph
epbwrite = .false.
epbread = .false.
! Coarse wannier grid M_e-ph and dynmat
epwwrite = .false.
epwread = .true.

wannierize = .false.
nbndsub = 5
nbndskip = 0

dis_win_min =   -25.000
dis_win_max =   15.000

dis_froz_min =  -25.000
dis_froz_max =    -1.0


num_iter = 3000
iprint = 2

proj(1) = 'f = 0.00000, 0.00000, 0.50000:sp3'
proj(2) = 'f = 0.33333, 0.66667, 0.50000:pz'
wdata(4) = 'Begin Kpoint_Path'
wdata(5) = 'G 0.0000  0.0000 0.0  M 0.5000  0.0000 0.0'
wdata(6) = 'M 0.5000  0.0000 0.0  K 0.3333  0.3333 0.0'
wdata(7) = 'K 0.3333  0.3333 0.0  G 0.0000  0.0000 0.0'
wdata(8)=  'End Kpoint_Path'
wdata(9) = 'bands_plot = .true.'
wdata(10) = 'dis_num_iter = 4000'
wdata(11) = 'guiding_centres = true'
wdata(12) = 'num_bands = 16'

!elecselfen = .true.
!phonselfen = .true.
!a2f = .true.
parallel_k = .true.
parallel_q = .false.
!
fsthick = 3.00000 ! eV
eptemp = 300 !
ngaussw = 0
degaussq = 0.03 ! meV
degaussw = 0.01 ! eV
efermi_read = .true.
fermi_energy = -3.3535 ! just coarse grid value

! elph and grid flags
dvscf_dir = './save'
filukk = './c_ukk'
! More grid control
prtgkk = .true.
band_plot = .true.
filkf = grid_3721.dat
filqf = grid_3721.dat
nk1 = 16
nk2 = 16
nk3 = 1
!
nq1 = 16
nq2 = 16
nq3 = 1
/
30 cartesian
0.000000000 0.000000000 0.000000000 0.0333
0.000000000 0.072168784 0.000000000 0.0333
0.000000000 0.144337567 0.000000000 0.0333
0.000000000 0.216506351 0.000000000 0.0333
0.000000000 0.288675135 0.000000000 0.0333
0.000000000 0.360843918 0.000000000 0.0333
0.000000000 0.433012702 0.000000000 0.0333
0.000000000 0.505181486 0.000000000 0.0333
0.000000000 -0.577350269 0.000000000 0.0333
0.062500000 0.108253175 0.000000000 0.0333
0.062500000 0.180421959 0.000000000 0.0333
0.062500000 0.252590743 0.000000000 0.0333
0.062500000 0.324759526 0.000000000 0.0333
0.062500000 0.396928310 0.000000000 0.0333
0.062500000 0.469097094 0.000000000 0.0333
0.062500000 0.541265877 0.000000000 0.0333
0.125000000 0.216506351 0.000000000 0.0333
0.125000000 0.288675135 0.000000000 0.0333
0.125000000 0.360843918 0.000000000 0.0333
0.125000000 0.433012702 0.000000000 0.0333
0.125000000 0.505181486 0.000000000 0.0333
0.125000000 0.577350269 0.000000000 0.0333
0.187500000 0.324759526 0.000000000 0.0333
0.187500000 0.396928310 0.000000000 0.0333
0.187500000 0.469097094 0.000000000 0.0333
0.187500000 0.541265877 0.000000000 0.0333
0.250000000 0.433012702 0.000000000 0.0333
0.250000000 0.505181486 0.000000000 0.0333
0.250000000 0.577350269 0.000000000 0.0333
0.312500000 0.541265877 0.000000000 0.0333

[carla.verdi]

Hi,

Don't worry for the length, it's good that you put so many details
After having a good look we couldn't spot obvious mistakes, so I would suggest a few checks:
- Can you plot the decays of the Hamiltonian, dynamical matrix, and matrix elements, to check that the decays of the matrix elements are similar to the ones of H and D, i.e. the matrix elements are properly interpolated?
- It seems that the coarse k and q grids used vary a bit in the literature; have you tried using slightly denser ones?
- [just to double check] Are you positive that you extracted from the output precisely the matrix elements corresponding to k=K as the initial state?

Best,
Carla

[andreyl]

Dear Carla

thank you very much for the answer.

- Can you plot the decays of the Hamiltonian, Dynamical matrix, and matrix elements, to check that the decays of the matrix elements are similar to the ones of H and D, i.e. the matrix elements are properly interpolated?

Hamiltonian and dyn matrix decay (H(Re) and Dyn(Rp) respectively)


E-p matrix elements decay (first and second row decay.epmat_wanep:
max_{m,n,nu} |g(m,n,nu;R_e,R_p)| (Re) / (Rp)
respectively, third row g from decay.epwane):


I would say that Hamiltonian and Dyn.mat decay rates are quite comparable with those in EPW CPC paper (http://dx.doi.org/10.1016/j.cpc.2010.08.027)
I am not quite sure how should I interpret what is going on for epmat decays, do they look reasonable?

- It seems that the coarse k and q grids used vary a bit in the literature; have you tried using slightly denser ones?

comparable.

Indeed they do, that is what I am testing now with 36x36x1 k grid for a start.
But, as I mentioned, https://arxiv.org/pdf/1705.01816.pdf used even smaller
coarse grids than mine (10x10x1 k and 5x5x1 q), still their g(qx,qy) looks reasonable.

- [just to double check] Are you positive that you extracted from the output precisely the matrix elements corresponding to k=K as the initial state?

That is how I grep the data from epw output file:

Code: Select all

grep -A 28 'ik =    2481 ' epw_elph.out

And what I obtain starts with

Code: Select all

ik =    2481 coord.:    0.3333333   0.3333333   0.0000000

From that point I parse the e-p matrix elements I need (5, 5, n).
If ik coord. is in crystal (relative/direct) coordinates I think it is correct.
I used the same special points convention for band structure as well.

Best wishes,
Andrey

[sponce]

Hello Andrey,

The reason we asked about your extraction script is because you seem to have vertical black strips (i.e. 0 values).
Such discontinuities is a bit strange. It might be worth taking a look at your script to be sure that you are not missing any values.

Also why do you have this white missing information on the top left edge of your picture?

The decays seems okish to me.

Lets hope a denser grid improves your results or that you made a silly mistake in one of your scripts (happens to everyone).

Best,
Samuel

[andreyl]

Dear Samuel,
thank you for the answer!

Here is the script:

Code: Select all

# band, band, mode indexes
n1=5
n2=5
n3=1

# The order of the loops is:
# {q loop {k loop}}
# ik 2481 yeilds g=K
# the grep call outputs the table of  all |g| for all i,j,nu
# second grep grabs the line, contaning the desired indexes
# epw_elph is an output file
grep -A 28 'ik =    2481 ' epw_elph.out | grep "$n1        $n2        $n3" > $n1$n2$n3\.dat
# thus $n1$n2$n3\.dat for ik = 2481 contains crystal coordinates of q, i,j,nu indexes
# E(k), E(k+q) and omega(nu,q)

# just to be sure that everything is collected
l_nu=`wc -l $n1$n2$n3\.dat | awk '{ print $1 }'`
echo 'points in g array per (band,band,mode):' $l_nu 

# merge the meaningfull data: 
# cart dat contains the grid_3721 points
# transformed to 2pi/a cartesian coordnates
# the same way it is done in QE
paste cart.dat $n1$n2$n3\.dat | awk '{ print $1, $2, $3, $8, $9, $10 }' > d$n1$n2$n3\.dat

(highlighted version is on dropbox paper: https://paper.dropbox.com/doc/Graphene-elph-9QVfbGljBpzkTLb8i7E1R)

And here are the results for coarse k=36x36x1.
It looks just a little bit better, but these strange lines-like structures still persist.

Such discontinuities is a bit strange. It might be worth taking a look at your script to be sure that you are not missing any values.

Yeah, they are the most weird parts of this distributions, not only g=0, but there are some other lines with constant values. I absolutely can't explain this.

Also why do you have this white missing information on the top left edge of your picture?

That is because of conversion from crystal (direct) to 2pi/a Cartesian coordinates. I generate meshes for calculation in direct coordinates but plot in 2pi/a.
I believe I use the same procedure QE uses for this purpose. |g| distributions in direct coordinates look as awful and E and Omega distributions look as good as in Cartesian (same picture as above for one of the modes (TA), q is in direct coordinates):

[sponce]

Hello,

It does look a bit better (but still the stripes).

I have never looked to much at 2D materials yet but I know that some cutoff might be required. This is not implemented at the moment.

If its not too expansive, I would suggest you to try those coarse grid to see if it improves:

Code: Select all

nk1 = 16
nk2 = 16
nk3 = 2
!
nq1 = 16
nq2 = 16
nq3 = 2

If I recall correctly, somebody in our group did that and it helped him for its 2D materials.
In principle you should not need this for a 2D materials but in practice it might required since there is no 2D truncation in the code.

Hope this works,
Best,
Samuel

[andreyl]

Dear Samuel

If I recall correctly, somebody in our group did that and it helped him for its 2D materials.
In principle you should not need this for a 2D materials but in practice it might required since there is no 2D truncation in the code.

This is relatively easy to do, I will indeed try it. Can you please specify, how does it work exactly (i.e. what exactly is to be truncated)?

Also one more question, I found this page
http://epw.org.uk/Main/Troubleshooting

One of the question-answer interested me:

My decay.H file shows no spatial decay

Your Wannier functions are probably very extended. When running the nscf calculation to feed wavefunctions to EPW, did you either use the same number of pools as processors, OR turn the wf_collect flag on? If you have a good match between Wannier-interpolated band structure and DFT band structure but still see no decay in decay.H, please contact us.

Although decays in my calculations seem to be alright, there are wf_collect and pools/processors count mentioned.

On nscf step I have wf_collect=.true. in input file and I run QE in the following manner:

Code: Select all

mpirun.sh -np 16 pw.x  -npool 16 < c.scf.in > c.scf.out
mpirun.sh -np 16 pw.x -npool 16 < c.nscf.in > c.nscf.out

Is it as intended, or it should be any other way?

Best wishes,
Andrey

[sponce]

Hello,

In the past you had to do wf_collect=.false. for the nscf step.

However I have implemented the capability of wf_collect=.true. so it should be fine.
You could try wf_collect=.false. in nscf to test. If it leads to slightly different results, please do contact me as it should be exactly the same results !

For the truncation, the dvscf potential might need to be truncated. See for example https://journals.aps.org/prb/abstract/1 ... .94.085415

Eventually I would like to implement that but it is not yet done.

Another idea from Carla: you could try the silicon case (should be very cheap). Try to do the same plot as for your graphene. This should work and not give stripes.

Best,
Samuel

[andreyl]

Thank you for the link.

Another idea from Carla: you could try the silicon case (should be very cheap). Try to do the same plot as for your graphene. This should work and not give stripes.

Do you mean the bulk silicon or silicene?

Best,
Andrey

[sponce]

Bulk silicon or anything simple and 3D.