Shift current

Introduction

In this tutorial we will show how to calculate Shift Current in bulk materials.
We suppose you are already familiar with the non-linear response using the Yambo code. If it is not the case please study the previous tutorials:
Linear response using Dynamical Berry Phase and Real time approach to non-linear response (SHG).
This tutorial was created with the help of Yuncheng Mao.

Setup calculations

In this tutorial we will take as example the two dimensional hBN.
DFT wave-functions and inputs can be downloaded here: hBN-2D-RT.tar.gz.

First of all run the setup, then remove symmetries along the y direction, as explained in the tutorial above.

Real-time setup and calculations

In order to generate input file for shift current you do: yambo_nl -u n -V par -F input.in

nloptics                         # [R] Non-linear spectroscopy
NLogCPUs=0                       # [PARALLEL] Live-timing CPU`s (0 for all)
PAR_def_mode= "balanced"         # [PARALLEL] Default distribution mode ("balanced"/"memory"/"workload"/"KQmemory")
NL_CPU= " 8 1"                       # [PARALLEL] CPUs for each role
NL_ROLEs= " w k"                     # [PARALLEL] CPUs roles (w,k)
DIP_CPU= ""                      # [PARALLEL] CPUs for each role
DIP_ROLEs= ""                    # [PARALLEL] CPUs roles (k,c,v)
OSCLL_CPU= ""                    # [PARALLEL] CPUs for each role
OSCLL_ROLEs= ""                  # [PARALLEL] CPUs roles (k,b)
DIP_Threads=0                    # [OPENMP/X] Number of threads for dipoles
NL_Threads=0                     # [OPENMP/NL] Number of threads for nl-optics
OSCLL_Threads=0                  # [OPENMP/X] Number of threads for Oscillators
% NLBands
  17 |  24 |                           # [NL] Bands range
%
NLverbosity= "high"              # [NL] Verbosity level (low | high)
NLtime=-1.000000           fs    # [NL] Simulation Time
NLintegrator= "INVINT"           # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA"             # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX/LSEX/LHF")
NLLrcAlpha= 0.000000             # [NL] Long Range Correction
% NLEnRange
 1.000000 | 8.000000  |         eV    # [NL] Energy range (for loop on frequencies NLEnSteps/=0
%
NLEnSteps=  24                    # [NL] Energy steps for the loop on frequencies
% NLrotaxis
 0.000000 | 0.000000 | 0.000000 |        # [NL] Rotation axis (for the loop on angles NLAngSteps/=0)
%
NLAngSteps=0                     # [NL] Angular steps (if NLAngSteps/=0 field versor will be ignored)
NLDamping= 0.200000        eV    # [NL] Damping (or dephasing)
RADLifeTime=-1.000000      fs    # [RT] Radiative life-time (if negative RADLifeTime=Phase_LifeTime)
EvalCurrent                    # [NL] Evaluate the current
#FrPolPerdic                   # [DIP] Force periodicity of polarization respect to the external field
% Field1_Freq
 0.100000 | 0.100000 |         eV    # [RT Field1] Frequency
%
Field1_NFreqs= 1                 # [RT Field1] Frequency
Field1_Int=  1000.00       kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000     fs    # [RT Field1] Width
Field1_kind= "SOFTSIN"           # [RT Field1] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field1_pol= "linear"             # [RT Field1] Pol(linear|circular)
% Field1_Dir
 0.000000 | 1.000000 | 0.000000 |        # [RT Field1] Versor 
%
Field1_Tstart= 0.010000    fs    # [RT Field1] Initial Time

Notice that in this input we turned one the evaluation of current EvalCurrent and force the parallelization on the frequencies,
NL_CPU= " 8 1, NL_ROLEs= " w k , that is much more efficient than the one on k-points.

Then you can run simulation by doing: yambo_nl -F input.in

Analysis of the results

In order to analyze the result we use Yambopy: YamboPy,
and the following script for the post-processing:

from yambopy import *
from yambopy.nl.harmonic_analysis import Harmonic_Analysis
from yambopy.units import fs2aut
X_order=4
NLDB=YamboNLDB()
Harmonic_Analysis(NLDB,X_order=X_order,prn_Peff=True)

This script will produce different file containing the response function of the current and polarization respect to the total field. From the second order response of the current at zero frequency, file o.YamboPy-Sigma_probe_order_0 we can extract the shift current coefficient along the 'y' direction that correspond to the 5th column of the file.
Hereafter we report the shift current spectrum compared with the results of Ibañez-Azpiroz et al. .

Then we used this script to plot the result:

import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot()
plt.title("$|\sigma^{(2)}(0,\omega_1-\omega_1$)|",fontsize=20)
fontsize2=14
XHI0 = np.genfromtxt('o.YamboPy-Sigma_probe_order_0',comments="#")
plt.tick_params(\
   axis='x',          # changes apply to the x-axis
   which='both',      # both major and minor ticks are affected
   bottom='off',      # ticks along the bottom edge are off
   top='off',         # ticks along the top edge are off
   labelbottom='off')
plt.xlim(1.0,7.5)
plt.xlabel('Energy [eV]',fontsize=20)
plt.tick_params(axis='y', which='major', labelsize=14)
# plt.plot(XHI0[:,0], XHI0[:,2], linewidth=2.5, linestyle="-",label='$\sigma^{(xyy)}$')
plt.plot(XHI0[:,0], XHI0[:,4], linewidth=2.5, linestyle="-",label='Re $\sigma^{(yyy)}$')
plt.axhline(y=0.0, color='black', linestyle='-')
plt.legend(loc="upper right",fontsize=fontsize2)
plt.xticks(fontsize=fontsize2)
plt.savefig('sigma2.png', format='png')
plt.show()