phonon bandstructure from EPW and matdyn.x don't match

[sridhu88]

Hello,
I'm trying to run EPW to calculate e-ph coupling in perovskite materials. I find that the Wannier-interpolated phonon bandstructure does not match the phonon bandstructure obtained directly from q2r.x and matdyn.x. Here's a plot of the phonon bandstructure from EPW (red) compared to the phonon bandstructure obtained from matdyn.x: https://www.dropbox.com/s/bjg9gd6fdcl26 ... w.png?dl=0

I'm also attaching a plot of the data in 'decay.dynmat': https://www.dropbox.com/s/u09xrkhhem6jw ... t.png?dl=0

As can be seen, the decay in the phonon dynamical matrix in Wannier representation isn't exactly great but I'm just using a 2 x 2 x 2 q-grid for phonon calculations because of the computational expense of these calculations. Also the unit cell size is a cube of dimension ~ 6 Angstrom, so I believe 2 x 2 x 2 is a reasonable start.

My understanding is that phonon dynamical matrix in Wannier representation is the same as inter-atomic force constants (which is what is obtained from q2r.x and matdyn.x). So I don't understand why the phonon dispersion from matdyn.x and EPW are different? I'm especially intrigued why some of the optical modes are flat in matdyn.x but show oscillations in EPW.

I'd appreciate it if someone could explain the possible reason(s) behind this. I have tried both EPW 4.0 and 4.1 with no difference.

Many thanks,
Sridhar

[carla.verdi]

Dear Sridhar,

Such a wrong interpolated phonon dispersion looks very weird indeed. Before thinking about possible problems in the epw interpolation, I would suggest to check if the problem comes from the start - i.e. if the dynamical matrices are being read correctly, and the phonon eigenvalues on the coarse grid are the same as with ph.x. I believe you can do this check without the need to go into the code - if you set iverbosity to 1, the frequencies on the coarse grid should be printed in the output.

Best
Carla

[sridhu88]

Dear Carla,
Thanks for your reply. I set iverbosity=1 and ran the calculation. Indeed, there is a discrepancy in the phonon frequencies from ph.x and the frequencies printed by EPW (see an example below for the 2nd irreducible q-pt). There appears to be a large difference for the low frequency modes while the high frequency modes are fairly close between ph.x and EPW. Also, I should note that since my calculation involves spin-orbit coupling, the dynamical matrix files from ph.x are in XML format.

Does this point to some error while reading the dynamical matrix?
Thanks again.
Sridhar

From ph.x:

Diagonalizing the dynamical matrix

q = ( 0.004197533 0.000109710 -0.502992536 )

**************************************************************************
freq ( 1) = -0.758668 [THz] = -25.306451 [cm-1]
freq ( 2) = 0.162267 [THz] = 5.412657 [cm-1]
freq ( 3) = 0.645132 [THz] = 21.519280 [cm-1]
freq ( 4) = 0.767770 [THz] = 25.610061 [cm-1]
freq ( 5) = 0.874158 [THz] = 29.158784 [cm-1]
freq ( 6) = 0.926805 [THz] = 30.914878 [cm-1]
freq ( 7) = 1.042757 [THz] = 34.782634 [cm-1]
freq ( 8) = 1.228815 [THz] = 40.988849 [cm-1]
freq ( 9) = 1.637458 [THz] = 54.619720 [cm-1]
freq ( 10) = 1.888619 [THz] = 62.997559 [cm-1]
freq ( 11) = 1.931020 [THz] = 64.411882 [cm-1]
freq ( 12) = 2.160524 [THz] = 72.067330 [cm-1]
freq ( 13) = 2.424439 [THz] = 80.870578 [cm-1]
freq ( 14) = 2.718140 [THz] = 90.667386 [cm-1]
freq ( 15) = 2.759219 [THz] = 92.037637 [cm-1]
freq ( 16) = 2.878967 [THz] = 96.031993 [cm-1]
freq ( 17) = 3.980414 [THz] = 132.772304 [cm-1]
freq ( 18) = 4.710575 [THz] = 157.127875 [cm-1]
freq ( 19) = 9.595943 [THz] = 320.086220 [cm-1]
freq ( 20) = 26.883049 [THz] = 896.721992 [cm-1]
freq ( 21) = 27.659692 [THz] = 922.628026 [cm-1]
freq ( 22) = 29.195669 [THz] = 973.862691 [cm-1]
freq ( 23) = 36.980033 [THz] = 1233.521115 [cm-1]
freq ( 24) = 37.537109 [THz] = 1252.103195 [cm-1]
freq ( 25) = 42.223193 [THz] = 1408.414132 [cm-1]
freq ( 26) = 43.354956 [THz] = 1446.165654 [cm-1]
freq ( 27) = 43.548445 [THz] = 1452.619779 [cm-1]
freq ( 28) = 44.842611 [THz] = 1495.788496 [cm-1]
freq ( 29) = 47.207603 [THz] = 1574.676141 [cm-1]
freq ( 30) = 47.886633 [THz] = 1597.326153 [cm-1]
freq ( 31) = 89.921333 [THz] = 2999.452819 [cm-1]
freq ( 32) = 92.855290 [THz] = 3097.319085 [cm-1]
freq ( 33) = 92.970008 [THz] = 3101.145655 [cm-1]
freq ( 34) = 94.174276 [THz] = 3141.315731 [cm-1]
freq ( 35) = 95.104234 [THz] = 3172.335779 [cm-1]
freq ( 36) = 95.148375 [THz] = 3173.808153 [cm-1]
**************************************************************************

From EPW:
===================================================================
irreducible q point # 2
===================================================================

Symmetries of small group of q: 1
in addition sym. q -> -q+G:

Number of q in the star = 1
List of q in the star:
1 0.004197533 0.000109710 -0.502992536
Frequencies of the matrix for the current q in the star
8.59405 13.82365 24.21344 25.82667 27.87061 29.07170
35.29590 41.91616 54.74756 65.38392 73.16179 74.78191
77.90095 88.84540 92.22573 107.29754 135.97740 149.31190
324.69575 896.63909 923.22700 973.89063 1233.09737 1252.78050
1408.93323 1447.46753 1453.30021 1495.74709 1574.70795 1598.35319
2999.24158 3097.26702 3100.79366 3140.56491 3171.27867 3173.49009

[carla.verdi]

Dear Sridhar,

As a matter of fact I think the differences you are seeing on that q point come from the acoustic sum rule. In fact, in EPW the asr is imposed immediately on the coarse grid, whereas in the phonon code it is imposed in real space (matdyn does it before interpolating back the phonon dispersion).
A couple of extra questions regarding your interpolated phonon dispersion with EPW:
- EPW gives you the phonon dispersion in meV, did you convert this in cm-1? I am sure you did but just checking.
- Did you put 'lpolar=.true.' in the epw.in? Otherwise the LO-TO splitting is not taken into account.

Best
Carla

[sridhu88]

Dear Carla,

But the frequencies I sent in my previous post were for a non-Gamma point. So I'm a bit confused on how ASR would play a role for that q-pt?

About your questions, yes, I did convert from meV to cm^-1 when plotting phband.freq. Also, I set lpolar = .true. in my calculation.

Once again, many thanks for your help.
Sridhar

[sridhu88]

I'm just sending the frequencies from ph.x and EPW for a couple of additional q-pts (see below). I also compared phonon frequencies between EPW and ph.x for a much simpler system (cubic BN) where the interpolation works fine. I saw that although the frequencies are almost the same, they are not exactly the same -- for example the frequency from ph.x is 1281.821 cm-1 while that from epw.out is 1281.7179 cm-1. I believe that's numerical noise due to differences in diagonalization routines and also probably because EPW loses some precision in the dynamical matrix when reading from the file. Is it possible that the small errors in the dynamical matrix show exaggerated errors in the soft low-frequency modes of the perovskite (MAPbI3) crystal I'm studying here?

Just a thought and I may be completely wrong.
Best,
Sridhar

From ph.x

Diagonalizing the dynamical matrix

q = ( 0.000168668 -0.508381476 0.000098126 )

**************************************************************************
freq ( 1) = 0.505328 [THz] = 16.855917 [cm-1]
freq ( 2) = 0.516492 [THz] = 17.228326 [cm-1]
freq ( 3) = 0.572471 [THz] = 19.095591 [cm-1]
freq ( 4) = 0.831062 [THz] = 27.721242 [cm-1]
freq ( 5) = 0.927202 [THz] = 30.928130 [cm-1]
freq ( 6) = 1.038924 [THz] = 34.654783 [cm-1]
freq ( 7) = 1.114696 [THz] = 37.182252 [cm-1]
freq ( 8) = 1.343627 [THz] = 44.818564 [cm-1]
freq ( 9) = 1.628899 [THz] = 54.334208 [cm-1]
freq ( 10) = 1.827189 [THz] = 60.948471 [cm-1]
freq ( 11) = 2.040067 [THz] = 68.049312 [cm-1]
freq ( 12) = 2.253948 [THz] = 75.183617 [cm-1]
freq ( 13) = 2.421498 [THz] = 80.772490 [cm-1]
freq ( 14) = 2.579582 [THz] = 86.045580 [cm-1]
freq ( 15) = 2.856184 [THz] = 95.272047 [cm-1]
freq ( 16) = 3.231171 [THz] = 107.780277 [cm-1]
freq ( 17) = 4.075960 [THz] = 135.959391 [cm-1]
freq ( 18) = 4.520714 [THz] = 150.794772 [cm-1]
freq ( 19) = 9.672628 [THz] = 322.644131 [cm-1]
freq ( 20) = 26.845737 [THz] = 895.477407 [cm-1]
freq ( 21) = 27.753159 [THz] = 925.745740 [cm-1]
freq ( 22) = 29.190972 [THz] = 973.706011 [cm-1]
freq ( 23) = 36.952509 [THz] = 1232.603032 [cm-1]
freq ( 24) = 37.573356 [THz] = 1253.312241 [cm-1]
freq ( 25) = 42.203298 [THz] = 1407.750506 [cm-1]
freq ( 26) = 43.370977 [THz] = 1446.700056 [cm-1]
freq ( 27) = 43.527810 [THz] = 1451.931462 [cm-1]
freq ( 28) = 44.772502 [THz] = 1493.449913 [cm-1]
freq ( 29) = 47.102587 [THz] = 1571.173179 [cm-1]
freq ( 30) = 47.984900 [THz] = 1600.603961 [cm-1]
freq ( 31) = 89.904799 [THz] = 2998.901277 [cm-1]
freq ( 32) = 92.856884 [THz] = 3097.372258 [cm-1]
freq ( 33) = 92.969222 [THz] = 3101.119432 [cm-1]
freq ( 34) = 94.110756 [THz] = 3139.196913 [cm-1]
freq ( 35) = 94.881408 [THz] = 3164.903091 [cm-1]
freq ( 36) = 95.658657 [THz] = 3190.829350 [cm-1]
**************************************************************************

From EPW:
===================================================================
irreducible q point # 3
===================================================================

Symmetries of small group of q: 1
in addition sym. q -> -q+G:

Number of q in the star = 1
List of q in the star:
1 0.000168668 -0.508381476 0.000098126
Frequencies of the matrix for the current q in the star
13.00647 15.86393 18.27632 23.00463 28.85729 35.71825
37.39739 41.79450 51.01952 66.65791 75.35999 75.95179
79.98413 85.35635 97.55271 121.60487 128.31661 153.40686
327.11242 895.38935 926.34116 973.73406 1232.19449 1253.99304
1408.26002 1447.99955 1452.62126 1493.39151 1571.20749 1601.63069
2998.68882 3097.32073 3100.76643 3138.39476 3163.88902 3190.51484


From ph.x:
Diagonalizing the dynamical matrix

q = ( 0.004366201 -0.508271765 -0.502894410 )

**************************************************************************
freq ( 1) = 0.458527 [THz] = 15.294800 [cm-1]
freq ( 2) = 0.669768 [THz] = 22.341049 [cm-1]
freq ( 3) = 0.810661 [THz] = 27.040756 [cm-1]
freq ( 4) = 0.865478 [THz] = 28.869240 [cm-1]
freq ( 5) = 0.930919 [THz] = 31.052132 [cm-1]
freq ( 6) = 1.042329 [THz] = 34.768343 [cm-1]
freq ( 7) = 1.331478 [THz] = 44.413330 [cm-1]
freq ( 8) = 1.488400 [THz] = 49.647672 [cm-1]
freq ( 9) = 1.836231 [THz] = 61.250074 [cm-1]
freq ( 10) = 2.048938 [THz] = 68.345215 [cm-1]
freq ( 11) = 2.218131 [THz] = 73.988872 [cm-1]
freq ( 12) = 2.361016 [THz] = 78.755005 [cm-1]
freq ( 13) = 2.594236 [THz] = 86.534392 [cm-1]
freq ( 14) = 2.761214 [THz] = 92.104197 [cm-1]
freq ( 15) = 2.897797 [THz] = 96.660104 [cm-1]
freq ( 16) = 3.136809 [THz] = 104.632677 [cm-1]
freq ( 17) = 4.068208 [THz] = 135.700827 [cm-1]
freq ( 18) = 4.529847 [THz] = 151.099446 [cm-1]
freq ( 19) = 9.649192 [THz] = 321.862401 [cm-1]
freq ( 20) = 26.808176 [THz] = 894.224506 [cm-1]
freq ( 21) = 27.698874 [THz] = 923.934983 [cm-1]
freq ( 22) = 29.191629 [THz] = 973.727946 [cm-1]
freq ( 23) = 36.872169 [THz] = 1229.923166 [cm-1]
freq ( 24) = 37.578308 [THz] = 1253.477447 [cm-1]
freq ( 25) = 42.199265 [THz] = 1407.615971 [cm-1]
freq ( 26) = 43.430602 [THz] = 1448.688938 [cm-1]
freq ( 27) = 43.520914 [THz] = 1451.701421 [cm-1]
freq ( 28) = 44.796797 [THz] = 1494.260296 [cm-1]
freq ( 29) = 47.188819 [THz] = 1574.049563 [cm-1]
freq ( 30) = 47.981352 [THz] = 1600.485637 [cm-1]
freq ( 31) = 89.906171 [THz] = 2998.947047 [cm-1]
freq ( 32) = 92.861891 [THz] = 3097.539275 [cm-1]
freq ( 33) = 92.965562 [THz] = 3100.997372 [cm-1]
freq ( 34) = 94.160782 [THz] = 3140.865620 [cm-1]
freq ( 35) = 94.871608 [THz] = 3164.576219 [cm-1]
freq ( 36) = 95.540824 [THz] = 3186.898849 [cm-1]
**************************************************************************

From EPW:
===================================================================
irreducible q point # 4
===================================================================

Symmetries of small group of q: 1
in addition sym. q -> -q+G:

Number of q in the star = 1
List of q in the star:
1 0.004366201 -0.508271765 -0.502894410
Frequencies of the matrix for the current q in the star
18.24692 23.53322 26.61737 26.96061 33.68623 33.92830
42.50112 56.81356 64.48117 70.74882 75.17201 80.55632
84.96034 90.15693 101.04761 113.78587 129.34646 152.58117
326.35203 894.12728 924.53491 973.75548 1229.51312 1254.15680
1408.13077 1449.98991 1452.38397 1494.21109 1574.08044 1601.50976
2998.73594 3097.48779 3100.64262 3140.14908 3163.47862 3186.58516

[carla.verdi]

Dear Sridhar,

The ASR changes the dynamical matrix (hence the frequencies) at every q, not just at gamma, either imposing it in reciprocal space or in real space - see for example Eq.(81) in http://journals.aps.org/prb/pdf/10.1103 ... B.55.10355
At this stage unfortunately I am not sure yet about what is happening for your system. Another important check to do would be to use a fine grid that is the same as the coarse grid, so that you can get the phonon frequencies again on the coarse grid but AFTER the interpolation. These should be the same - if they are not it gives additional information as in what can be wrong in this case.
Ps. if you input the coarse q-points in a 'filqf' file, remember the coordinates should be crystal instead of cartesian as EPW normally works with crystal coordinates after the interpolation.

Best
Carla

[nadeemnatt]

Hi Sridhar and Carla

How did you solve this issue? I am having same issue in my calculations of hydrocarbons.

Regards
Nadeem