Electron-phonon matrix elements with spin orbit coupling (SO

Post here questions linked with issue while running the EPW code

Moderator: stiwari

Post Reply
wwwennie

Electron-phonon matrix elements with spin orbit coupling (SO

Post by wwwennie »

Dear EPW developers,

I would like to ask a question concerning the calculation of electron-phonon matrix elements with spin-orbit coupling (SOC).

I have calculated with and without SOC (denoted as SOC and NSOC) the electron-phonon matrix elements interpolated along the high-symmetry k-point path for the cubic symmetry.
The matrix elements for NSOC diverge as expected for q -> 0 for the long range component.
However, this was not the case for matrix elements with SOC. Instead, for q -> 0, the matrix elements are suppressed to zero, as if the long range component was not included.

Please find a plot with the magnitude of the electron-phonon matrix elements at the k = Gamma point for averaged over transitions involving the 3-fold degenerate bands here:
Image

I have separately checked the norm-conserving pseudopotentials used against other pseudopotentials and other codes, and found them to be consistent with each other in terms of effective band masses and the SO slitting of the bands.
I have also checked the wannierization for NSOC and SOC and find the band structure to agree well with the original QE calculations.

Please find below the input files for running scf, nscf, and epw.
Could you advise on where a possible error might be (e.g., input) or if this is a known issue?

Thank you for your time.
Sincerely,
Wennie
------------------------------------------------------------------------------------------------------
epw.in
------------------------------------------------------------------------------------------------------

Code: Select all

--
&inputepw
  prefix      = 'WO3'
  amass(1)    = 183.84
  amass(2)    = 15.999
  outdir      = './'

  elph        = .true.
  ep_coupling = .true.
  kmaps       = .true.
  epbwrite    = .false.
  epbread     = .false.

  etf_mem     = .false.
  lpolar      = .true.
 
  !! ww add for printing out matrix elements
  lgelph      = .true.
  lfmt        = .true.

  epwwrite    = .false.
  epwread     = .true.

  !ephwrite    = .true.  ! write el-ph elements of fine mesh to file, unformatted

  nbndsub     =  42
  nbndskip    =  0

  wannierize  = .false.
  num_iter    = 500
  dis_win_min =-75.000
  dis_win_max = 20.000
  dis_froz_min= -5.000
  dis_froz_max = 10.000
  !dis_win_max = 11
  !dis_froz_max= 11
  proj(1)     = 'W:s,p,d'
  proj(2)     = 'O:s,p'
  wdata(1)    = 'bands_plot=.true.'
  wdata(2)    = 'spinors=.true.'
  wdata(3)    = 'begin kpoint_path' 
  wdata(4)    = 'G 0.0000  0.0000  0.0000     X 0.000  0.5000  0.0000'   
  wdata(5)    = 'X 0.0000  0.5000  0.0000     M 0.500  0.5000  0.0000'
  wdata(6)    = 'M 0.5000  0.5000  0.0000     G 0.000  0.0000  0.0000'
  wdata(7)    = 'G 0.0000  0.0000  0.0000     R 0.500  0.5000  0.5000'
  wdata(8)    = 'R 0.5000  0.5000  0.5000     X 0.000  0.5000  0.0000'
  wdata(9)    = 'end kpoint_path'

  elecselfen  = .false.  ! scattering rate -> ~ tau
  phonselfen  = .false.
  a2f         = .false.

  vme         = .false. ! .true. does not work

  parallel_k  = .true.
  parallel_q  = .false.

  fsthick     = 5.0
  eptemp      = 300.0
  degaussw    = 0.1
  efermi_read = .true.
  fermi_energy= 6.7 !approx mid-gap;  VBM = 6.1500, CBM =  6.8526

  dvscf_dir   = './save'
  filqf       = 'testqpt'
  filkf       = 'testkpt'

  nk1         = 6
  nk2         = 6
  nk3         = 6
  nq1         = 6
  nq2         = 6
  nq3         = 6
 /
20 cartesian
   0.000000000   0.000000000   0.000000000   1
 ...



------------------------------------------------------------------------------------------------------
scf.in
------------------------------------------------------------------------------------------------------

Code: Select all

&CONTROL
    calculation     = 'scf'
    restart_mode    = 'from_scratch',
    prefix          = 'WO3',
    tstress         = .true.
    tprnfor         = .true.
    pseudo_dir      = './',
    outdir          = './'
    wf_collect      = .true.
    verbosity       = 'high'   
/
&SYSTEM
    ibrav           =  1,
    celldm(1)       =  7.1749995265,
    nat             =  4,   
    ntyp            =  2,     
    ecutwfc         =  125.00,
    nbnd            =  48,
    lspinorb        = .true.,
    noncolin        = .true.
    la2f            = .true.
/
&ELECTRONS
    diagonalization = 'cg'
    mixing_mode     = 'plain'
    mixing_beta     = 0.3
    conv_thr        = 1.0d-10
    diago_full_acc  = .true.
/
&IONS
    ion_dynamics = 'damp'
/
&CELL
    cell_dynamics = 'bfgs'
/                           
ATOMIC_SPECIES
 W   183.84  W_frl.upf
 O   15.999  O_frl.upf
ATOMIC_POSITIONS (crystal)
 W    0.000000000         0.000000000         0.000000000
 O    0.500000000         0.000000000         0.000000000
 O    0.000000000         0.500000000         0.000000000
 O    0.000000000         0.000000000         0.500000000
K_POINTS automatic
6 6 6 0 0 0

sponce
Site Admin
Posts: 616
Joined: Wed Jan 13, 2016 7:25 pm
Affiliation: EPFL

Re: Electron-phonon matrix elements with spin orbit coupling

Post by sponce »

Dear Wennie,

Thank you for your interest in the EPW code.

Could you clarify exactly what you did to get the figure ?

The black dots are QE/PH calculations (so adding write statements in the phonon code) and the blue line is the EPW interpolation. Is this correct?

This is the way it was done in the EPW paper for GaN.

Providing that this is correct, then are you saying that there is an issue in the Phonon code with SOC ? Is this correct?

Edit: I have implemented the possiblity to print the vertex |g| inside EPW (see viewtopic.php?f=6&t=353&p=1256#p1256).
I'm guessing you hack the code yourself to do this and it looks OK but you might be interested by double checking.

Best,
Samuel
Prof. Samuel Poncé
Chercheur qualifié F.R.S.-FNRS / Professeur UCLouvain
Institute of Condensed Matter and Nanosciences
UCLouvain, Belgium
Web: https://www.samuelponce.com
wwwennie

Re: Electron-phonon matrix elements with spin orbit coupling

Post by wwwennie »

Dear Samuel,

Thank you for your prompt reply! Please find below answers to your questions.

To produce the attached figure:
- I hacked the ephwann_shuffle.f90 routine with a modified version of what was described in a previous post (see http://epwforum.uk/viewtopic.php?f=3&t=34&p=147&hilit=print+electron+phonon#p147)
- The hacked version prints the real and imaginary parts of g_elph
- The figure consists of averaging over degenerate k-states; the phonon mode is singly degenerate.

Apologies for the lack of clarity:
All black dots are taken from "interpolating" along the high-symmetry path for a cubic structure using EPW from a calculation using a 6x6x6 scf and 6x6x6 phonon grid.
The solid line for NSOC corresponds to the same thing and is used as a guide to the eye to help distinguish between the qpts chosen for NSOC and SOC cases.

Rephrasing the question:
The g_elph matrix elements for NSOC diverge as q-> 0 with EPW. I would have expected to see the same for SOC since the material is polar, but do not.
When selecting for only the long-range or short-range part, I observe that the long-range part is essentially zero.
In looking through the EPW code, it seems the spinors are taken into account at the beginning in memory allocation, which would funnel into the calculation of g_elph for both the long- and short-range contributions.
However, my understanding of how SOC is handled in the Wannier representation is limited, so I am wondering if there is something I am missing, either conceptually or practically in running the code.

Thank you for the heads up on the added prtgkk input flag. I will look into this and compare.

Sincerely,
Wennie
sponce
Site Admin
Posts: 616
Joined: Wed Jan 13, 2016 7:25 pm
Affiliation: EPFL

Re: Electron-phonon matrix elements with spin orbit coupling

Post by sponce »

Dear Wennie,

I can offer the following suggestions:

1) Do the same figures for all modes. You calculation with SOC might change the phonon numbering order and you might be printing a non-polar mode.
2) Make sure to check your Born effective charges with SOC and that they are correct
3) Print the full phonon dispersion with and without SOC and make sure you have the correct LO-TO splitting in both case
4) SOC doubles the number of electronic bands. Make sure you are using the one you expect during the calculations

Best,
Samuel
Prof. Samuel Poncé
Chercheur qualifié F.R.S.-FNRS / Professeur UCLouvain
Institute of Condensed Matter and Nanosciences
UCLouvain, Belgium
Web: https://www.samuelponce.com
wwwennie

Re: Electron-phonon matrix elements with spin orbit coupling

Post by wwwennie »

Dear Samuel,

Thank you for the suggestions and patience.

1) I have checked the other modes. The mode plotted is closest to corresponding to the polar mode and consistently corresponds to the polar mode.
The mode of interest is well separated from any other modes by at least 250 cm^-1.

2) I have checked the Born effective charges with SOC and find they are comparable to NSOC and other
non-collinear calculations with other pseudopotentials. For reference, I've included them below.

3) In checking the phonon dispersion, I find the NSOC and SOC dispersion to have the appropriate LO-TO
splitting and frequencies differing by at most a few cm^-1. This is in QE, i.e., with matdyn.

I have also checked if phonon modes are properly being read in, using the iverbosity = 1 tag in the input file,
and find the frequencies to be comparable.
However, I do find that while the output file says it imposes the ASR, the frequencies outputted are unchanged, and so the expected LO-TO splitting is not observed, example output below.
Previously discussed http://epwforum.uk/viewtopic.php?f=3&t=137

4) The previously attached plot does include the change and doubling in electronic band indices as well as lowering of degeneracy from NSOC to SOC, and adjusts the average of g_elph accordingly to take that into account.

I have also tested the code with prntgkk added functionality using the same input files as before.
I was missing averaging over k+q. However, I still run into the same problem as explained before.
That is, the long-range component of g_elph in SO calculations all output to zero and the initial
divergence as q->0 as seen in the NSOC calculations is no longer there.
This seems to be most likely related to the lack of LO-TO splitting mentioned earlier.

My best guess is to look at the read-in routines for the dynamical matrices.
Do you have any further suggestions?

Sincerely,
Wennie
- - - - - - - - - - - - - - - - - - - - - -
excerpt from EPW output

Code: Select all

  ===================================================================
   irreducible q point #    1
   ===================================================================

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

   Number of q in the star =    1
   List of q in the star:
        1   0.000000000   0.000000000   0.000000000
      Read dielectric tensor and effective charges
      Imposing acoustic sum rule on the dynamical matrix
   Frequencies of the matrix for the current q in the star
-295.33519  -295.33519  -295.33519    -0.00001     0.00000     0.00001
 224.69599   224.69599   224.69599   565.68021   565.68021   565.68021


excerpt from matdyn

Code: Select all

     diagonalizing the dynamical matrix ...

 q =       0.0000      0.0000      0.0000
 **************************************************************************
     freq (    1) =      -8.763603 [THz] =    -292.322334 [cm-1]
 (  0.148420  -0.000000     0.000000   0.000000     0.000111  -0.000000   )
 ( -0.643546   0.000000     0.000000  -0.000000    -0.000396  -0.000000   )
 ( -0.530951   0.000000     0.000000   0.000000    -0.000396   0.000000   )
 ( -0.530951   0.000000    -0.000000   0.000000    -0.000480   0.000000   )
     freq (    2) =      -8.763603 [THz] =    -292.322334 [cm-1]
 (  0.000111  -0.000000    -0.000000  -0.000000    -0.148420  -0.000000   )
 ( -0.000480   0.000000     0.000000  -0.000000     0.530951  -0.000000   )
 ( -0.000396   0.000000    -0.000000  -0.000000     0.530951   0.000000   )
 ( -0.000396   0.000000    -0.000000   0.000000     0.643546   0.000000   )
     freq (    3) =      -0.000000 [THz] =      -0.000008 [cm-1]
 (  0.006518  -0.000000     0.000922  -0.000000    -0.499957  -0.000000   )
 (  0.006518  -0.000000     0.000922  -0.000000    -0.499957   0.000000   )
 (  0.006518  -0.000000     0.000922  -0.000000    -0.499957  -0.000000   )
 (  0.006518  -0.000000     0.000922  -0.000000    -0.499957   0.000000   )
     freq (    4) =      -0.000000 [THz] =      -0.000004 [cm-1]
 (  0.469462   0.000000    -0.171964  -0.000000     0.005803   0.000000   )
 (  0.469462   0.000000    -0.171964  -0.000000     0.005803  -0.000000   )
 (  0.469462   0.000000    -0.171964  -0.000000     0.005803   0.000000   )
 (  0.469462   0.000000    -0.171964  -0.000000     0.005803   0.000000   )
     freq (    5) =       0.000000 [THz] =       0.000004 [cm-1]
 (  0.171939  -0.000000     0.469497  -0.000000     0.003107   0.000000   )
 (  0.171939  -0.000000     0.469497  -0.000000     0.003107   0.000000   )
 (  0.171939  -0.000000     0.469497  -0.000000     0.003107  -0.000000   )
 (  0.171939  -0.000000     0.469497  -0.000000     0.003107   0.000000   )
     freq (    6) =       7.060748 [THz] =     235.521214 [cm-1]
 (  0.000000  -0.000000     0.000000   0.000000    -0.000000  -0.000000   )
 ( -0.000000   0.000000    -0.078310   0.000000     0.702411  -0.000000   )
 ( -0.022067  -0.000000     0.000000   0.000000    -0.702411   0.000000   )
 (  0.022067   0.000000     0.078310  -0.000000    -0.000000   0.000000   )
     freq (    7) =       7.060748 [THz] =     235.521214 [cm-1]
 ( -0.000000   0.000000     0.000000   0.000000     0.000000  -0.000000   )
 ( -0.000000  -0.000000    -0.678679   0.000000    -0.069828  -0.000000   )
 (  0.185793  -0.000000    -0.000000   0.000000     0.069828   0.000000   )
 ( -0.185793   0.000000     0.678679  -0.000000     0.000000   0.000000   )
     freq (    8) =       7.060748 [THz] =     235.521214 [cm-1]
 (  0.000000  -0.000000     0.000000  -0.000000    -0.000000   0.000000   )
 (  0.000000  -0.000000    -0.182380   0.000000    -0.041756   0.000000   )
 ( -0.681905   0.000000     0.000000  -0.000000     0.041756  -0.000000   )
 (  0.681905  -0.000000     0.182380  -0.000000     0.000000   0.000000   )
     freq (    9) =      11.609088 [THz] =     387.237495 [cm-1]
 (  0.000000   0.000000     0.116782  -0.000000     0.000000   0.000000   )
 (  0.000000   0.000000    -0.700983   0.000000     0.000000   0.000000   )
 ( -0.000000  -0.000000     0.060059  -0.000000    -0.000000  -0.000000   )
 ( -0.000000  -0.000000    -0.700983   0.000000     0.000000   0.000000   )
     freq (   10) =      17.090891 [THz] =     570.090774 [cm-1]
 ( -0.000002  -0.000000     0.000000  -0.000000     0.011129   0.000000   )
 ( -0.000173  -0.000000    -0.000000   0.000000    -0.449735   0.000000   )
 (  0.000101   0.000000     0.000000   0.000000    -0.449735  -0.000000   )
 (  0.000101   0.000000    -0.000000   0.000000     0.771591   0.000000   )
     freq (   11) =      17.090891 [THz] =     570.090774 [cm-1]
 (  0.011129   0.000000    -0.000000   0.000000     0.000002   0.000000   )
 (  0.771591   0.000000     0.000000  -0.000000    -0.000101   0.000000   )
 ( -0.449735  -0.000000    -0.000000   0.000000    -0.000101   0.000000   )
 ( -0.449735  -0.000000     0.000000  -0.000000     0.000173   0.000000   )
     freq (   12) =      26.747517 [THz] =     892.201145 [cm-1]
 ( -0.000000   0.000000    -0.080630  -0.000000     0.000000   0.000000   )
 ( -0.000000   0.000000    -0.034527  -0.000000     0.000000  -0.000000   )
 ( -0.000000  -0.000000     0.995547   0.000000    -0.000000   0.000000   )
 (  0.000000  -0.000000    -0.034527  -0.000000    -0.000000   0.000000   )
 **************************************************************************



- - - - - - - - - - - - - - - - - - - - - -

srl-ph-cubic
---------------------------------------

Code: Select all

    Program PHONON v.6.0 (svn rev. 13079) starts on 26Jan2017 at 23:36:57
         Dielectric constant in cartesian axis

    (       9.831707195       0.000000000       0.000000000 )
    (       0.000000000       9.831707195      -0.000000000 )
    (       0.000000000       0.000000000       9.831707195 )

          Effective charges (d Force / dE) in cartesian axis

           atom      1   W     
      Ex  (       13.72323        0.00000       -0.00000 )
      Ey  (       -0.00000       13.72323        0.00000 )
      Ez  (       -0.00000        0.00000       13.72323 )
           atom      2   O     
      Ex  (      -10.12893        0.00000        0.00000 )
      Ey  (        0.00000       -1.86620        0.00000 )
      Ez  (        0.00000       -0.00000       -1.86620 )
           atom      3   O     
      Ex  (       -1.86620       -0.00000       -0.00000 )
      Ey  (        0.00000      -10.12893       -0.00000 )
      Ez  (        0.00000       -0.00000       -1.86620 )
           atom      4   O     
      Ex  (       -1.86620        0.00000        0.00000 )
      Ey  (        0.00000       -1.86620        0.00000 )
      Ez  (       -0.00000        0.00000      -10.12893 )

================================================
frl-ph-cubic

Code: Select all

Program PHONON v.6.0 (svn rev. 13079) starts on 30Jan2017 at 20:50: 2

          Dielectric constant in cartesian axis

          (      10.136153744       0.000000000       0.000000000 )
          (       0.000000000      10.136153744       0.000000000 )
          (       0.000000000       0.000000000      10.136153744 )

          Effective charges (d Force / dE) in cartesian axis

           atom      1   W
      Ex  (       13.87320        0.00000       -0.00000 )
      Ey  (       -0.00000       13.87320       -0.00000 )
      Ez  (       -0.00000       -0.00000       13.87320 )
           atom      2   O
      Ex  (      -10.21654        0.00000       -0.00000 )
      Ey  (        0.00000       -1.90022       -0.00000 )
      Ez  (        0.00000       -0.00000       -1.90022 )
           atom      3   O
      Ex  (       -1.90022       -0.00000        0.00000 )
      Ey  (       -0.00000      -10.21654        0.00000 )
      Ez  (       -0.00000       -0.00000       -1.90022 )
           atom      4   O
      Ex  (       -1.90022       -0.00000        0.00000 )
      Ey  (       -0.00000       -1.90022       -0.00000 )
      Ez  (       -0.00000       -0.00000      -10.21654 )

sponce
Site Admin
Posts: 616
Joined: Wed Jan 13, 2016 7:25 pm
Affiliation: EPFL

Re: Electron-phonon matrix elements with spin orbit coupling

Post by sponce »

Hello,

Yes it's definitely the problem. So,

1) EPW relies on DFPT. Perturbation theory works only close to the minimum of the vibration energy landscape. Therefore those negative/soft mode -295.33519 -295.33519 -295.33519 are really problematic. You should fix this. You structure might not be the stable one or you are not converged well enough. It is ok to have small negative fre (like -10 cm^-1) but not so big.

2) However I agree there seems to be an issue with the LO-TO splitting. So natively EPW implement the "simple" ASR.
However you can change this and use the "crystal" sum rule.

You can try "asr_typ = 'crystal' " see http://epw.org.uk/Documentation/Inputs#asr_typ
To use this, you need to put to true the flag "lifc" in EPW.
You also need to generate the real-space IFC using q2r.x
The resulting file must be named "ifc.q2r"

Note however, that this real-space ASR does NOT work for SOC at the moment.

Best,
Samuel
Prof. Samuel Poncé
Chercheur qualifié F.R.S.-FNRS / Professeur UCLouvain
Institute of Condensed Matter and Nanosciences
UCLouvain, Belgium
Web: https://www.samuelponce.com
wwwennie

Re: Electron-phonon matrix elements with spin orbit coupling

Post by wwwennie »

Dear Samuel,

Thank you for the helpful information.

1) Yes, we are aware that the structure we are using is not stable for the DFPT calculations being done.
However, could this be a cause for the ASR routines to not work?

2) Thank you for the suggestion.
I have tried other ways to impose ASR.
Interestingly, with iverbosity = 1 for the NSOC calculations, the phonon frequencies being printed out are
similar to those in SOC, i.e., it does not appear that ASR was imposed successfully or that the LO-TO splitting
was properly captured.
I attempted this with both the 'simple' and 'crystal' schemes of ASR.
I noticed that the actual code for these schemes are more or less directly ported from the q2r.f90 and q2trans.f90
routines, so it's curious that the printed phonon frequencies when running EPW
do not match with what QE outputs with q2r.x and matdyn.x.

When you say it does not work for SOC, do you mean that the code is unable to properly handle the *.xml
format, or something else?

Sincerely,
Wennie
sponce
Site Admin
Posts: 616
Joined: Wed Jan 13, 2016 7:25 pm
Affiliation: EPFL

Re: Electron-phonon matrix elements with spin orbit coupling

Post by sponce »

Hello,

Indeed the code is directly taken from q2r routines with some adaptation for EPW.

I was wrong, it does work with SOC (I forgot that it was working). You just need to call the file with .xml.

I updated the info accordingly in the documentation: http://epw.org.uk/Documentation/Inputs#lifc

PS: I would need to check but the iverbosity > 0 sometimes gives debugging information and might not always be what you expect.
band_plot will give you the used el and or ph frequencies.

Best,
Samuel
Prof. Samuel Poncé
Chercheur qualifié F.R.S.-FNRS / Professeur UCLouvain
Institute of Condensed Matter and Nanosciences
UCLouvain, Belgium
Web: https://www.samuelponce.com
wwwennie

Re: Electron-phonon matrix elements with spin orbit coupling

Post by wwwennie »

Dear Samuel,

I have tested different ways to impose the ASR for both the NSOC and SOC cases.

While the NSOC was not sensitive to the choice of ASR algorithm, the SOC calculations were.
In fact, the 'simple' scheme was not sufficiently able to capture the correct LO-TO splitting, whereas the 'crystal' scheme was more suitable.
This was in fact responsible for the strange lack of divergence in the electron-phonon matrix elements as q -> 0.

Understanding why this is case still eludes me.
My cursory guess is that the space spanned by the spinors is not amenable to a simple subtractions scheme and then imposition of index symmetries of the IFCs.
Any additional insights would be appreciated.
Any new findings with printing for iverbosity would be appreciated as well.

Thank you again for your help.
Sincerely,
Wennie
sponce
Site Admin
Posts: 616
Joined: Wed Jan 13, 2016 7:25 pm
Affiliation: EPFL

Re: Electron-phonon matrix elements with spin orbit coupling

Post by sponce »

Hello Wennie,

The iverbosity > 0 is mainly used for debugging or in some special cases to have some additional output files. If those are reported in the doc, then they should be good.

However printing inside some routine might not be what you expect. What I mean by this is that it might be printing stuff on each core (so before an mpi gather or so).

For the reason why the crystal sum rule is better than the simple and why the simple fails in some case, you can read the very good PhD thesis of Nicolas Mounet.

Best,
Samuel
Prof. Samuel Poncé
Chercheur qualifié F.R.S.-FNRS / Professeur UCLouvain
Institute of Condensed Matter and Nanosciences
UCLouvain, Belgium
Web: https://www.samuelponce.com
Post Reply