EPW to calculate anisotropic superconducting gap of MgB2

[Cai Cheng]

Dear Dr. Samuel,

I want to calculate the anisotropic superconducting gap of MgB2 according to website:http://epw.org.uk/Documentation/MgB2

First, I try to run phonon dispersion relation to obtain the .dyn* and .dvscf*
there are no problem.

Then, I calculate the epw.x

Code: Select all

[color=#BF0080]epw.in[/color]
--
&inputepw
  prefix      = 'MgB2',
  amass(1)    = 24.305,
  amass(2)    = 10.811
  outdir      = './'

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

  epwwrite = .true.
  epwread  = .false.

  etf_mem     = .true.

  nbndsub     =  5
  nbndskip    =  0

  wannierize  = .true.
  num_iter    = 500
  dis_froz_max= 8.8
  proj(1)     = 'B:pz'
  proj(2)     = 'f=0.5,1.0,0.5:s'
  proj(3)     = 'f=0.0,0.5,0.5:s'
  proj(4)     = 'f=0.5,0.5,0.5:s'

  iverbosity  = 2

  elinterp    = .true.
  phinterp    = .true.

  tshuffle2   = .true.
  tphases     = .false.

  parallel_k  = .true.
  parallel_q  = .false.

  eps_acustic = 5.0    ! Lowest boundary for the 
  ephwrite    = .true. ! Writes .ephmat files used when wliasberg = .true.

  fsthick     = 0.4  ! eV
  eptemp      = 300  ! K
  degaussw    = 0.10 ! eV
  nsmear      = 1
  delta_smear = 0.04 ! eV

  degaussq     = 0.5 ! meV
  nqstep       = 500

  eliashberg  = .true.

  laniso = .true.
  limag = .true.
  lpade = .true.

  conv_thr_iaxis = 1.0d-2

  wscut = 0.525   ! eV   Upper limit over frequency integration/summation in the Elisashberg eq

  nstemp   = 1
  tempsmin = 25.00
  tempsmax = 30.00

  nsiter   = 1000

  muc     = 0.16

  dvscf_dir   = '../phonons/save'
  
  nk1         = 6
  nk2         = 6
  nk3         = 6

  nq1         = 6
  nq2         = 6
  nq3         = 6

  mp_mesh_k = .true.
  nkf1 = 20
  nkf2 = 20
  nkf3 = 20

  nqf1 = 20
  nqf2 = 20
  nqf3 = 20
 /
  28 cartesian
  0.000000000   0.000000000   0.000000000  0.004629630
  0.000000000   0.000000000   0.145933920  0.009259259
  0.000000000   0.000000000   0.291867841  0.009259259
  0.000000000   0.000000000  -0.437801761  0.004629630
  0.000000000   0.192450090   0.000000000  0.027777778
  0.000000000   0.192450090   0.145933920  0.055555556
  0.000000000   0.192450090   0.291867841  0.055555556
  0.000000000   0.192450090  -0.437801761  0.027777778
  0.000000000   0.384900179   0.000000000  0.027777778
  0.000000000   0.384900179   0.145933920  0.055555556
  0.000000000   0.384900179   0.291867841  0.055555556
  0.000000000   0.384900179  -0.437801761  0.027777778
  0.000000000  -0.577350269   0.000000000  0.013888889
  0.000000000  -0.577350269   0.145933920  0.027777778
  0.000000000  -0.577350269   0.291867841  0.027777778
  0.000000000  -0.577350269  -0.437801761  0.013888889
  0.166666667   0.288675135   0.000000000  0.027777778
  0.166666667   0.288675135   0.145933920  0.055555556
  0.166666667   0.288675135   0.291867841  0.055555556
  0.166666667   0.288675135  -0.437801761  0.027777778
  0.166666667   0.481125224   0.000000000  0.055555556
  0.166666667   0.481125224   0.145933920  0.111111111
  0.166666667   0.481125224   0.291867841  0.111111111
  0.166666667   0.481125224  -0.437801761  0.055555556
  0.333333333   0.577350269   0.000000000  0.009259259
  0.333333333   0.577350269   0.145933920  0.018518519
  0.333333333   0.577350269   0.291867841  0.018518519
  0.333333333   0.577350269  -0.437801761  0.009259259

some of the output is incorrenct

Code: Select all

   
[color=#FF0000]epw.out[/color]

  Symmetries of small group of q: 12
 
     Number of q in the star =    2
     List of q in the star:
          1   0.333333333   0.577350269  -0.437801761
          2  -0.333333333  -0.577350269  -0.437801761
 
        q(  215 ) = (   0.3333333   0.5773503  -0.4378018 )
        q(  216 ) = (  -0.3333333  -0.5773503  -0.4378018 )

     Writing epmatq on .epb files


              band disentanglement is used:  nbndsub =    5

     Writing Hamiltonian, Dynamical matrix and EP vertex in Wann rep to file


     Reading Hamiltonian, Dynamical matrix and EP vertex in Wann rep from file


     Finished reading Wann rep data from file

     Using uniform q-mesh:   20  20  20
     Size of q point mesh for interpolation:       8000
     Using uniform MP k-mesh:   20  20  20
     Size of k point mesh for interpolation:        968
     Max number of k points per pool:               42

     Fermi energy coarse grid =   7.567624 eV

     Fermi energy is calculated from the fine k-mesh: Ef =   7.900198 eV

     Warning: check if difference with Fermi level fine grid makes sense

     ===================================================================

              ibndmin =     3  ebndmin =     0.551
              ibndmax =     5  ebndmax =     0.610


     Number of ep-matrix elements per pool :         1701 ~=   13.29 Kb (@ 8 bytes/ DP)

     Nr. of irreducible k-points on the uniform grid:       484


     Finished writing .ikmap file 


     Finished mapping k+sign*q onto the fine irreducibe k-mesh

     Nr irreducible k-points within the Fermi shell =       104 out of       484
                  Fermi level (eV) =    0.790019830659852D+01
     DOS(states/spin/eV/Unit Cell) =    0.219570491047075D+00
            Electron smearing (eV) =    0.100000000000000D+00
                 Fermi window (eV) =    0.400000000000000D+00
      
     Finished writing .ephmat files

     ===================================================================
     Solve anisotropic Eliashberg equations
     ===================================================================


     Finish reading .freq file 

                  Fermi level (eV) =     7.9001983068E+00
     DOS(states/spin/eV/Unit Cell) =     2.1957049093E-01
            Electron smearing (eV) =     1.0000000000E-01
                 Fermi window (eV) =     4.0000000000E-01
     Nr irreducible k-points within the Fermi shell =       104 out of       484
           2 bands within the Fermi window


     Finish reading .egnv file 


     Max nr of q-points =      1599


     Finish reading .ikmap files


     Start reading .ephmat files


     Finish reading .ephmat files 

     lambda_max =           221.5756907   lambda_k_max =             1.1007955
 
     Electron-phonon coupling strength =    0.1415293
 
     Estimated Allen-Dynes Tc = *************** K for muc =    0.16000
  
     Estimated BCS superconducting gap = *************** eV
  
     temp(  1) =  25.0000 K 
  
     Solve anisotropic Eliashberg equations on imaginary-axis 
  
     Total number of frequency points nsiw (      1 ) =     39
     Cutoff frequency wscut =     0.5347
  

     Size of allocated memory per pool : ~=    0.0255 Gb
     iter =      1   relerr =    0.557210204D+16   abserr =    0.236761053D+15   Znormi(1) =    0.100392116D+01   Deltai(1) =   -0.392920088D-01
     iter =      2   relerr =    0.111878749D+00   abserr =    0.427544549D-02   Znormi(1) =    0.100390563D+01   Deltai(1) =   -0.347549316D-01
     iter =      3   relerr =    0.338365872D-01   abserr =    0.125074400D-02   Znormi(1) =    0.100389517D+01   Deltai(1) =   -0.334261973D-01
     iter =      4   relerr =    0.309152553D-02   abserr =    0.114630298D-03   Znormi(1) =    0.100389632D+01   Deltai(1) =   -0.335476482D-01
     Convergence was reached in nsiter =      4
  
     iaxis_imag   :      2.43s CPU      2.60s WALL (       1 calls)
  
     Pade approximant of anisotropic Eliashberg equations from imaginary-axis to real-axis
     Cutoff frequency wscut =     0.5250
  
     pade =     36   error =    0.129066073D-01   Re[Znorm(1)] =    0.112393012D+01   Re[Delta(1)] =   -0.335820500D-01
  
     raxis_pade   :      0.85s CPU      2.27s WALL (       1 calls)
  
     itemp =   1   total cpu time :              4.87 secs
 
     Unfolding on the coarse grid
     elphon_wrap  :    682.33s CPU    814.38s WALL (       1 calls)
 
     INITIALIZATION: 
 
     set_drhoc    :      5.82s CPU      5.84s WALL (     217 calls)
     init_vloc    :      0.33s CPU      0.36s WALL (     218 calls)
     init_us_1    :      1.47s CPU      1.55s WALL (     218 calls)
 
 
     Electron-Phonon interpolation
     ephwann      :    228.62s CPU    262.77s WALL (       1 calls)
     ep-interp    :    206.87s CPU    232.76s WALL (    8000 calls)
 
     Ham: step 1  :      0.00s CPU      0.00s WALL (       1 calls)
     Ham: step 2  :      0.02s CPU      0.03s WALL (       1 calls)
     ep: step 1   :      0.01s CPU      0.04s WALL (    1944 calls)
     ep: step 2   :      3.12s CPU      9.14s WALL (    1944 calls)
     DynW2B       :      0.16s CPU      1.16s WALL (    8000 calls)
     HamW2B       :     14.12s CPU     14.39s WALL (  336042 calls)
     ephW2Bp      :    155.31s CPU    158.38s WALL (    8000 calls)
 
     ELIASHBERG   :     72.32s CPU     75.03s WALL (       1 calls)
 
     Total program execution
     EPW          : 16m36.32s CPU    19m33.79s WALL

I find that the lambda_max = 221.5756907 lambda_k_max = 1.1007955

Electron-phonon coupling strength = 0.1415293

Estimated Allen-Dynes Tc = *************** K for muc = 0.16000

Estimated BCS superconducting gap = *************** eV

Is there some unreasonable papameters in the input file (epw.in)?
It's pleasure to get your advice

Best wishes

Cai Cheng

Email: ccheng@iphy.ac.cn
TEl: +86-10-82648045 (office)

Institute of Physics, Chinese Academy of Sciences

[sponce]

Dear Cai,

Your input seems ok. Off course nkf and nqf grids needs to be increase to achieve convergence but you should get acceptable number already.
Maybe try to increase

Code: Select all

wscut = 1.0

Also, when the numbers gets very large, it can create the ***** issue (depending on compiler).

You can make the following modification in eliashberg_setup.f90 line 589:

Code: Select all

WRITE(stdout,'(5x,a,1E15.7,a,f10.5)') 'Estimated Allen-Dynes Tc = ', tc, ' K for muc = ', muc

The same can be done for the BCS gap.

You can then recompile.

Best,

Samuel

[Cai Cheng]

Thank you for you clearly reply

Then, according your advice, I have change some parameters in epw.in

Code: Select all

--
&inputepw
  prefix      = 'MgB2',
  amass(1)    = 24.305,
  amass(2)    = 10.811
  outdir      = './'

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

  epwwrite = .true.
  epwread  = .false.

  etf_mem     = .true.

  nbndsub     =  5
  nbndskip    =  0

  wannierize  = .true.
  num_iter    = 500
  dis_froz_max= 8.8
  proj(1)     = 'B:pz'
  proj(2)     = 'f=0.5,1.0,0.5:s'
  proj(3)     = 'f=0.0,0.5,0.5:s'
  proj(4)     = 'f=0.5,0.5,0.5:s'

  iverbosity  = 2

  elinterp    = .true.
  phinterp    = .true.

  tshuffle2   = .true.
  tphases     = .false.

  parallel_k  = .true.
  parallel_q  = .false.

  eps_acustic = 5.0    ! Lowest boundary for the 
  ephwrite    = .true. ! Writes .ephmat files used when wliasberg = .true.

  fsthick     = 1.0  ! eV
  eptemp      = 300  ! K
  degaussw    = 0.10 ! eV
  nsmear      = 1
  delta_smear = 0.04 ! eV

  degaussq     = 0.5 ! meV
  nqstep       = 500

  eliashberg  = .true.

  laniso = .true.
  limag = .true.
  lpade = .true.

  conv_thr_iaxis = 1.0d-4

  wscut = 0.4   ! eV   Upper limit over frequency integration/summation in the Elisashberg eq

  nstemp   = 1
  tempsmin = 10.00
  tempsmax = 30.00

  nsiter   = 1000

  muc     = 0.16

  dvscf_dir   = '../phonons/save'
  
  nk1         = 6
  nk2         = 6
  nk3         = 6

  nq1         = 6
  nq2         = 6
  nq3         = 6

  mp_mesh_k = .true.
  nkf1 = 60
  nkf2 = 60
  nkf3 = 60

  nqf1 = 30
  nqf2 = 30
  nqf3 = 30
 /
  28 cartesian
  0.000000000   0.000000000   0.000000000  0.004629630
  0.000000000   0.000000000   0.145933920  0.009259259
  .................................................................................
  0.333333333   0.577350269   0.291867841  0.018518519
  0.333333333   0.577350269  -0.437801761  0.009259259

Code: Select all

...................................................................
Parallel version (MPI), running on    24 processors
     K-points division:     npool     =      24
.............................................................
     G cutoff =  206.3461  (  12303 G-vectors)     FFT grid: ( 30, 30, 36)
     number of k points=  216  gaussian broad. (Ry)=  0.0200     ngauss =   1
.......................................................................................................

     lambda_max =          1401.2509107   lambda_k_max =             3.0757448
 
     Electron-phonon coupling strength =    0.1000114
 
     Estimated Allen-Dynes Tc = *************** K for muc =    0.16000
  
     Estimated BCS superconducting gap =  334621.0054370 eV
  
     temp(  1) =  10.0000 K 
  
     Solve anisotropic Eliashberg equations on imaginary-axis 
.......................................................................................
     AKeri is calculated on the fly since its size exceedes max_memlt

     iter =      1   relerr =    0.976601019D+07   abserr =    0.334621040D+06   Znormi(1) =    0.100209383D+01   Deltai(1) =   -0.327169173D-01
     iter =      2   relerr =    0.669194853D-01   abserr =    0.214910194D-02   Znormi(1) =    0.100208606D+01   Deltai(1) =   -0.304518325D-01
     iter =      3   relerr =    0.251548732D-01   abserr =    0.788019704D-03   Znormi(1) =    0.100326317D+01   Deltai(1) =   -0.296066170D-01
     iter =      4   relerr =    0.311340827D+01   abserr =    0.461495654D-01   Znormi(1) =    0.105772670D+01   Deltai(1) =    0.129652805D-01
     iter =      5   relerr =    0.307403365D+01   abserr =    0.219697118D-01   Znormi(1) =    0.107362333D+01   Deltai(1) =   -0.569268894D-02
     iter =      6   relerr =    0.164905987D+01   abserr =    0.181579747D-01   Znormi(1) =    0.105306935D+01   Deltai(1) =    0.964044120D-02
     iter =      7   relerr =    0.314748025D+01   abserr =    0.835621621D-02   Znormi(1) =    0.108483893D+01   Deltai(1) =    0.224175688D-02
     iter =      8   relerr =    0.217644033D+01   abserr =    0.491160595D-02   Znormi(1) =    0.108922821D+01   Deltai(1) =   -0.165581935D-02
     iter =      9   relerr =    0.353303267D+01   abserr =    0.314762956D-02   Znormi(1) =    0.109258294D+01   Deltai(1) =    0.685569360D-03
     iter =     10   relerr =    0.147102307D+01   abserr =    0.530369627D-03   Znormi(1) =    0.109375557D+01   Deltai(1) =    0.262745202D-03
     iter =     11   relerr =    0.208503722D+01   abserr =    0.692832637D-03   Znormi(1) =    0.109
..............................................................................................................................................

MgB2.a2f

Code: Select all

 

...........................................................................................................................................................
Integrated el-ph coupling
  #            0.1000114   0.1001133   0.1002300   0.1003635   0.1005162   0.1006900   0.1008856   0.1011024   0.1013385   0.1015910
 Phonon smearing (meV)
  #            0.5000000   0.5500000   0.6000000   0.6500000   0.7000000   0.7500000   0.8000000   0.8500000   0.9000000   0.9500000
 Electron smearing (eV)   0.1000000
 Fermi window (eV)   1.0000000
 Summed el-ph coupling    0.0995704

To summary,I have several questions,

1. What the difference about the code modified?

Code: Select all

WRITE(stdout,'(5x,a,1E15.7,a,f10.5)') 'Estimated Allen-Dynes Tc = ', tc, ' K for muc = ', muc

2. Why obtain the small lambda about 0.1 (it should be about 0.75)?

Code: Select all

Electron smearing (eV)   0.1000000
Fermi window (eV)   1.0000000
Summed el-ph coupling    0.0995704

Thank you for your reply again.

[sponce]

Dear Cai,

I previously suggested that you increase your value to

Code: Select all

wscut = 1.0

but it seems like you decreased it to

Code: Select all

wscut = 0.4

.

To answer your questions:
1) The difference is just the exponential format (E) instead of the float format (f). If the numbers get very big, the exponential format is more suited.

2) As reported in the technical paper, I got a lambda of 0.75 using a fine k grid of 60x60x60 and a fine q grid of 30x30x30 with an ecut (for scf/nscf run of 60Ry) and a wscut of 1.0 eV and a degaussw of 0.1 eV.

Best,

Samuel

[roxana]

Dear Cai,

I would advise that you try to reproduce the results from the test directory EPW/tests/Refs/t05/paral
The input parameters won't give you converged results, but you can directly check your output files with the references given.

The input files are located in EPW/tests/Inputs/t05

Best,
Roxana