Off-diagonal component calculation for optical conductivity

[pereng7180]

Dear all,

I am trying to calculate the optical conductivity matrix through Quantum Espresso with Yambo.
In order to compare the result, calculation by Abinit with Yambo is also implemented according to tutorial in Kerr effect section.
Here is system and computational details.


[System Details]
BCC Fe

[Computational Details]
Abinit 7.4.3 and Yambo 3.4.2
Quantum Espresso(QE) 5.1.2 and Yambo 3.4.2

[Pseudopotential]
Abinit: HGH pseudopotential with semi-core (26fe.16.hgh)
QE: HGH pseudopotential with semi-core (Fe.pbe-sp-hgh.UPF)


Here is the result of two programs.
https://drive.google.com/open?id=0Bz9z8 ... m1vaUZJRkk
Diagonal components show similar trend.
But, off-diagonal component of QE+Yambo is very small and shape is different comparing to Abinit+Yambo result.
What makes the difference between two programs?


Here is the input of Abinit and QE calculation.
[Abinit]

Code: Select all

#Fe in FCC structure

 ndtset  2


# Common variables
####################
# system
 acell   3*5.42
 natom   1
 fband1  2
 ntypat  1
 typat   1
 znucl   26.0
 rprim   -0.5  0.5  0.5
          0.5 -0.5  0.5
          0.5  0.5 -0.5
 xred    0.0  0.0  0.0
 occopt  7
 tsmear  0.01

# convergence parameters
 ecut       50
 toldfe1  1.0d-8
 nstep   100

# kpts
 kptopt     4
 ngkpt      8 8 8
 nshiftk    1
 shiftk     0.5  0.5  0.5
#nshiftk    2
#shiftk     0.25  0.25  0.25
#          -0.25 -0.25 -0.25

# Magnetism and so
 nspinor    2
 nsppol     1
 nspden     4
 spinat     0.0  0.0  5.0 
 so_psp     1


# other
 ixc     1    # xc approximation (should be consistent with psp)
 enunit  1    # in which format eigenvalues are printed 0 (Ha), 1 (eV), 2 (both)

 nsym2       0
 symmorphi2  0

#DATASET 2: WFK GENERATION
 iscf2        -2               # Non self-consistent calculation
 getden2       1               # Read previous density file
 tolwfr2       1.0d-10         # Still get it converged
 nstep2      300
 nband2      34
 nbandkss2   30
 kssform2     3
 nbdbuf2     10

[QE]
1st calculation

Code: Select all

 &control
    prefix='Fe_SCF'
    outdir = './temp/'
    pseudo_dir = './'
    etot_conv_thr = 1.0D-5
 /
 &system    
    ibrav=  3, celldm(1) =5.42, nat=  1, ntyp= 1,
    ecutwfc =300.0, 
    occupations = 'smearing',
    smearing = 'gaussian',
    degauss = 0.0001
    nbnd = 36
    nosym = .true.
    noncolin = .false.
    nspin = 2
    starting_magnetization(1) = 1.0
 /
 &electrons
    electron_maxstep = 100
 /
ATOMIC_SPECIES
 Fe 55.845 Fe.pbe-sp-hgh.UPF

ATOMIC_POSITIONS (alat)
 Fe 0.0 0.0 0.0

K_POINTS (automatic)
  8 8 8 0.5 0.5 0.5

2nd calculation

Code: Select all

 &control
    prefix='Fe_SCF'
    outdir = './temp/'
    pseudo_dir = './'
    wf_collect=.true.
 /
 &system    
    ibrav=  3, celldm(1) =5.42, nat=  1, ntyp= 1,
    ecutwfc =300.0, 
    occupations = 'smearing',
    smearing = 'gaussian',
    degauss = 0.0001
    nbnd = 36
    nosym = .true.
    noncolin = .true.
    starting_spin_angle = .true.
    lspinorb = .true.
    starting_magnetization(1) = 1.0
    constrained_magnetization = 'atomic'
    angle1(1) = 0.0
    angle2(1) = 90.0
    lambda = 5.0
 /
 &electrons
    mixing_beta = 0.5
 /
ATOMIC_SPECIES
 Fe 55.845 Fe.pbe-sp-hgh.UPF

ATOMIC_POSITIONS (alat)
 Fe 0.0 0.0 0.0

K_POINTS (automatic)
  8 8 8 0.5 0.5 0.5

There was no error during generation of yambo databases.

I used same Yambo input for input databases from both programs.
[Yambo Input]

Code: Select all

#  __  __   ________   ___ __ __    _______   ______        
# /_/\/_/\ /_______/\ /__//_//_/\ /_______/\ /_____/\       
# \ \ \ \ \\::: _  \ \\::\| \| \ \\::: _  \ \\:::_ \ \      
#  \:\_\ \ \\::(_)  \ \\:.      \ \\::(_)  \/_\:\ \ \ \     
#   \::::_\/ \:: __  \ \\:.\-/\  \ \\::  _  \ \\:\ \ \ \    
#     \::\ \  \:.\ \  \ \\. \  \  \ \\::(_)  \ \\:\_\ \ \   
#      \__\/   \__\/\__\/ \__\/ \__\/ \_______\/ \_____\/   
#                                                           
#               Version 3.3.1 Revision 2085                 
#                http://www.yambo-code.org                  
#
optics                       # [R OPT] Optics
bse                          # [R BSE] Bethe Salpeter Equation.
BSEmod= "causal"             # [BSE] resonant/causal/coupling
BSKmod= "IP"                 # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
Gauge= "velocity"            # [BSE] Gauge (length|velocity). In metals length misses An-Hall
EvalKerr                     # [BSE] Compute the Kerr effect
DrudeWBS= ( 4.800000 , 5.650000 )  eV  # [BSE] Drude plasmon
% BEnRange
 0.000000 | 8.000000 | eV    # [BSS] Energy range
%
% BDmRange
 0.300000 | 0.500000 | eV    # [BSS] Damping range
%
BEnSteps=  400               # [BSS] Energy steps
% BLongDir
 1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field
%
% BSEBands
  1 | 30 |                   # [BSK] Bands range
%

Thank you!
Best regards,
KISUNG KANG

[Davide Sangalli]

Dear Kisng Kang,
can you also upload the report files created by yambo?
Both with abinit and with pwscf.

My best guess is that the SOC is not properly described in this pwscf calculation.
You should not specify nspin=2 and noncolin=.false. in the scf step.
Also nspin=2 should be not there in the nscf step.
Just use noncolin=.true. and lspinor=.true. in both cases

Also pay attention that the 0.5 0.5 0.5 shift of QE is not the 0.5 0.5 0.5 shift of abinit.
I think it is more like
nshiftk 2
shiftk 0.25 0.25 0.25
-0.25 -0.25 -0.25
in abinit. You can check the total number of kpts.

Finally you do not need nosym=.true. , just force_symmorphyic=.true. in pwscf

Best,
D.

[pereng7180]

Thank you for response.
I couldn't upload the report files because of large number of characters.
So, I use dropbox link to share them.
Here are the report files created by yambo.
[Abinit]
https://www.dropbox.com/s/df7o2bcudzcec ... t.txt?dl=0
[QE]
https://www.dropbox.com/s/pde6vkxzaggt379/QE.txt?dl=0


Based on your advice, I change the input of QE.
[1st calculation]

Code: Select all

 &control
    prefix='Fe_SCF'
    outdir = './temp/'
    pseudo_dir = './'
    etot_conv_thr = 1.0D-5
 /
 &system    
    ibrav=  3, celldm(1) =5.42, nat=  1, ntyp= 1,
    ecutwfc =300.0, 
    occupations = 'smearing',
    smearing = 'gaussian',
    degauss = 0.0001
    nbnd = 36
    force_symmorphic=.true.
    starting_magnetization(1) = 1.0
 /
 &electrons
    electron_maxstep = 100
 /
ATOMIC_SPECIES
 Fe 55.845 Fe.pbe-sp-hgh.UPF

ATOMIC_POSITIONS (alat)
 Fe 0.0 0.0 0.0

K_POINTS (automatic)
  8 8 8 0.5 0.5 0.5

[2nd calculation]

Code: Select all

 &control
    prefix='Fe_SCF'
    outdir = './temp/'
    pseudo_dir = './'
    wf_collect=.true.
 /
 &system    
    ibrav=  3, celldm(1) =5.42, nat=  1, ntyp= 1,
    ecutwfc =300.0, 
    occupations = 'smearing',
    smearing = 'gaussian',
    degauss = 0.0001
    nbnd = 36
    force_symmorphic= .true.
    noncolin = .true.
    starting_spin_angle = .true.
    lspinorb = .true.
    starting_magnetization(1) = 1.0
    constrained_magnetization = 'atomic'
    angle1(1) = 0.0
    angle2(1) = 90.0
    lambda = 5.0
 /
 &electrons
    mixing_beta = 0.5
 /
ATOMIC_SPECIES
 Fe 55.845 Fe.pbe-sp-hgh.UPF

ATOMIC_POSITIONS (alat)
 Fe 0.0 0.0 0.0

K_POINTS (automatic)
  8 8 8 0.5 0.5 0.5

After those calculations, I try to convert wavefunction file by using 'p2y'.
But it couldn't and it shows error message as following.

Code: Select all

 ___ __  _____  __ __  _____   _____
|   Y  ||  _  ||  Y  ||  _  \ |  _  |
|   |  ||. |  ||.    ||. |  / |. |  |
 \_  _/ |. _  ||.\_/ ||. _  \ |. |  |
  |: |  |: |  ||: |  ||: |   \|: |  |
  |::|  |:.|:.||:.|:.||::.   /|::.  |
  `--"  `-- --"`-- --"`-----" `-----"

 <---> P(W) 2 Y(ambo) Ver. 5.0
 <---> DBs path set to .
 <---> Index file set to data-file.xml
 <---> Header/K-points/Energies...done
 <---> Cell data...done
 <---> Atomic data...done
 <---> Symmetries...[SI yes]...[I yes]...[-I yes]...[TR yes]
 <---> [MAG-SYMS] T-rev with (A) Magnetic field or (B) SOC and magnetization
 <---> XC functional...Perdew, Burke & Ernzerhof(X)+Perdew, Burke & Ernzerhof(C)
 <---> K-points mesh...done
 <---> RL vectors...

[ERROR] STOP signal received while in :
[ERROR]Error in qexml_read_planewaves IOTK error ierr:  1

How can I deal with this?
Thanks you

Best regards,
KISUNG KANG

[Davide Sangalli]

It maybe a problem with yambo 3.4
Could you try the last patch of yambo 4.1 and let us know if it works?

Meanwhile I'll also try to test your system.

Best,
D.

P.S.: One thing I noticed just now.
Is this input working for you ?

Code: Select all

K_POINTS (automatic)
  8 8 8 0.5 0.5 0.5

I think it should be

Code: Select all

K_POINTS (automatic)
  8 8 8  1 1 1

[Davide Sangalli]

One more observation.
The idea of doing two consecutive calculations with pwscf is that (i) first the scf density is computed and (ii) then the empty states are be computed via an nscf calculation. The two steps should be done with the same "approximations" and in the second step the nscf mode of pwscf can be set.

Code: Select all

    calculation='nscf'

In your previous post I see instead the nscf variable was not set-up in the second step.

Best,
D.

[pereng7180]

Based on your advice, I use Yambo 4.1.2 version.
Before QE+Yambo calculation, I compare the results from different version of Yambo.
https://drive.google.com/open?id=0B3eXI ... 3B3OFEwUEE

With 3.4.2, it shows similar results in tutorial. But, 4.1.2 shows different pattern.
I used same KSS file from abinit and same input for Yambo.
Only difference is version of a2y and yambo_kerr.

Here is input of abinit and yambo.
[Abinit]

Code: Select all

#Fe in FCC structure

 ndtset  2


# Common variables
####################
# system
 acell   3*5.42
 natom   1
 fband1  2
 ntypat  1
 typat   1
 znucl   26.0
 rprim   -0.5  0.5  0.5
          0.5 -0.5  0.5
          0.5  0.5 -0.5
 xred    0.0  0.0  0.0
 occopt  7
 tsmear  0.01

# convergence parameters
 ecut       50
 toldfe1  1.0d-8
 nstep   100

# kpts
 kptopt     4
 ngkpt      8 8 8
 nshiftk    1
 shiftk     0.5  0.5  0.5
#nshiftk    2
#shiftk     0.25  0.25  0.25
#          -0.25 -0.25 -0.25

# Magnetism and so
 nspinor    2
 nsppol     1
 nspden     4
 spinat     0.0  0.0  5.0 
 so_psp     1


# other
 ixc     1    # xc approximation (should be consistent with psp)
 enunit  1    # in which format eigenvalues are printed 0 (Ha), 1 (eV), 2 (both)

 nsym2       0
 symmorphi2  0

#DATASET 2: WFK GENERATION
 iscf2        -2               # Non self-consistent calculation
 getden2       1               # Read previous density file
 tolwfr2       1.0d-10         # Still get it converged
 nstep2      300
 nband2      34
 nbandkss2   30
 kssform2     3
 nbdbuf2     10

[Yambo]

Code: Select all

#  __  __   ________   ___ __ __    _______   ______        
# /_/\/_/\ /_______/\ /__//_//_/\ /_______/\ /_____/\       
# \ \ \ \ \\::: _  \ \\::\| \| \ \\::: _  \ \\:::_ \ \      
#  \:\_\ \ \\::(_)  \ \\:.      \ \\::(_)  \/_\:\ \ \ \     
#   \::::_\/ \:: __  \ \\:.\-/\  \ \\::  _  \ \\:\ \ \ \    
#     \::\ \  \:.\ \  \ \\. \  \  \ \\::(_)  \ \\:\_\ \ \   
#      \__\/   \__\/\__\/ \__\/ \__\/ \_______\/ \_____\/   
#                                                           
#               Version 3.3.1 Revision 2085                 
#                http://www.yambo-code.org                  
#
optics                       # [R OPT] Optics
bse                          # [R BSE] Bethe Salpeter Equation.
BSEmod= "causal"             # [BSE] resonant/causal/coupling
BSKmod= "IP"                 # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
Gauge= "velocity"            # [BSE] Gauge (length|velocity). In metals length misses An-Hall
EvalKerr                     # [BSE] Compute the Kerr effect
DrudeWBS= ( 4.800000 , 5.650000 )  eV  # [BSE] Drude plasmon
% BEnRange
 0.000000 | 8.000000 | eV    # [BSS] Energy range
%
% BDmRange
 0.300000 | 0.500000 | eV    # [BSS] Damping range
%
BEnSteps=  400               # [BSS] Energy steps
% BLongDir
 1.000000 | 0.000000 | 0.000000 |        # [BSS] [cc] Electric Field
%
% BSEBands
  1 | 30 |                   # [BSK] Bands range
%

For QE calculation, I try to recalculate it with new input.
I understand that I have to do nscf calculation.
But, I feel confused about detail process.
Do you mean that another nscf calculation should be done after two calculation I mentioned before?
Or Should 2nd calculation contain "calculation='nscf'" term?

Thanks a lot.
Best regards,
KISUNG KANG

[Davide Sangalli]

For pwscf.
Yes, the second calculation should be nscf.
Unless you were doing two calculations for some other reason

For the difference you are experiencing in between the two versions of yambo.
Can you attach the two reports ?
You should be able to attach files, instead of putting them in the text. (Maybe you need to rename them.)

Best,
D.

[pereng7180]

I finally find the attachment button. Thanks
Report files from different version of Yambo with Abinit are attached.


For PWscf,
I calculate with new input based on your advice (force_symmorphic, Kpoint etc.)
[1st calculation]

Code: Select all

 &control
    prefix='Fe_SCF'
    outdir = './temp/'
    pseudo_dir = './'
    etot_conv_thr = 1.0D-5
 /
 &system    
    ibrav=  3, celldm(1) =5.42, nat=  1, ntyp= 1,
    ecutwfc =300.0, 
    occupations = 'smearing',
    smearing = 'gaussian',
    degauss = 0.0001
    nbnd = 36
    force_symmorphic=.true.
    starting_magnetization(1) = 1.0
 /
 &electrons
    electron_maxstep = 100
 /
ATOMIC_SPECIES
 Fe 55.845 Fe.pbe-sp-hgh.UPF

ATOMIC_POSITIONS (alat)
 Fe 0.0 0.0 0.0

K_POINTS (automatic)
  8 8 8 1 1 1

After 1st calculation, I change the input for 2nd calculation with "calculation='nscf'" and 1st calculation's wavefunction file.
[2nd calculation with nscf input]

Code: Select all

 &control
    calculation='nscf'
    prefix='Fe_SCF'
    outdir = './temp/'
    pseudo_dir = './'
    wf_collect=.true.
 /
 &system    
    ibrav=  3, celldm(1) =5.42, nat=  1, ntyp= 1,
    ecutwfc =300.0, 
    occupations = 'smearing',
    smearing = 'gaussian',
    degauss = 0.0001
    nbnd = 36
    force_symmorphic= .true.
    noncolin = .true.
    starting_spin_angle = .true.
    lspinorb = .true.
    starting_magnetization(1) = 1.0
    constrained_magnetization = 'atomic'
    angle1(1) = 0.0
    angle2(1) = 90.0
    lambda = 5.0
 /
 &electrons
    mixing_beta = 0.5
 /
ATOMIC_SPECIES
 Fe 55.845 Fe.pbe-sp-hgh.UPF

ATOMIC_POSITIONS (alat)
 Fe 0.0 0.0 0.0

K_POINTS (automatic)
  8 8 8 1 1 1

But, it shows error message and calculation is terminated.

Code: Select all

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     task #         7
     from read_rho_xml : error #        10
     searching for ./temp/Fe_SCF.save/magnetization.x.xml
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Thus, I calculate it again without "calculation='nscf'".
There is no error message but it did scf calculation.
Then, with wavefunction file from this scf calculation, I just add "calculation='nscf'" in input file and do nscf calculation.
From this calculation, there is no error message with nscf calculation.
I attach figure for result and report from different version of Yambo.
The result still shows huge difference between Abinit+Yambo and PWscf+Yambo.

Thank you
Best regards,
KISUNG KANG

[pereng7180]

Here are report files from QWscf+Yambo.

[Davide Sangalli]

Ok. This what I think is happening.

1) For the problem with the nscf calculation in QE the reason is that the scf and the ncf stap are not done at the same level of approximations.
I think in one case (scf) you are not doing the full non collinear calculation with SOC included while in the nscf the parameters are there

2) For the fact that QE+yambo 4.1.2 shows almost zero kerr signal I suspect the reason is that in the pwscf calculation the SOC is not fully accounted for.
This should also explain the difference in the eps_xx

3) eps_xx abinit + yambo 3.4 vs 4.1 are on top. I think the difference you instead experience in the eps_xy is due to a different description of the smearing.
In the first implementation (3.4) we had a non exact description of the smearing parameter in the current based approach.
This is why here http://journals.aps.org/prb/abstract/10 ... .86.125139 we had to onsider the zero smearing limit in figure to prove the equivalence between length and velocity gauge.
In the new implementation we fixed this limit of the velocity gauge. Indeed the solution was to replace the 1/w^2 term (which probably gives the divergence you see at low omega) with 1/(w+i\eta^2), plus other details.
In any case try to plot sigma_xy "propto" w*eps_xy if you want to see better the off diagonal component.
http://www.yambo-code.org/tutorials/KER ... /index.php
Also try to compute the effect in the length gauge adding on top the AnHall effect. Thi is usually more stable than the velocity gauge

4) Also pay attention that you are using two different energy cutoff in the two calculations

Ok, this is more or less all I noticed.
Hope it helps.

Best,
D.