diff --git a/note/note.tm b/note/note.tm
index 7bbc3b0..554f3ab 100644
--- a/note/note.tm
+++ b/note/note.tm
@@ -1,4 +1,4 @@
-<TeXmacs|2.1.1>
+<TeXmacs|2.1.2>
 
 <style|<tuple|generic|invisible-multiply|reduced-margins|old-lengths|number-long-article>>
 
@@ -310,6 +310,24 @@
   <\eqnarray*>
     <tformat|<table|<row|<cell|\<Lambda\><rsup|\<alpha\>\<beta\>\<gamma\>>>|<cell|=>|<cell|-<frac|e<rsup|3>|\<hbar\>><big|sum><rsub|n><big|int><rsub|FS<rsub|n>><frac|\<mathd\>\<sigma\>|<around*|(|2\<mathpi\>|)><rsup|3>><around*|(|\<partial\><rsub|\<alpha\>>\<varepsilon\><rsub|n>|)><around*|(|\<partial\><rsub|\<beta\>>\<partial\><rsub|\<gamma\>>\<varepsilon\><rsub|n>|)><frac|1|<around*|\||\<nabla\>\<varepsilon\><rsub|n>|\|>>.<eq-number>>>>>
   </eqnarray*>
+
+  <section|Injection current>
+
+  The injection current conductivity is
+
+  <\eqnarray*>
+    <tformat|<table|<row|<cell|\<sigma\><rsup|\<alpha\>\<beta\>\<gamma\>><around*|(|\<omega\>|)>>|<cell|=>|<cell|<frac|\<mathpi\>e<rsup|3>|\<eta\>><big|sum><rsub|n,m><big|int><frac|\<mathd\><rsup|3>\<b-k\>|<around*|(|2\<mathpi\>|)><rsup|3>><around*|(|v<rsup|\<alpha\>><rsub|n>-v<rsub|m><rsup|\<alpha\>>|)>A<rsup|\<beta\>><rsub|n\<nocomma\>m>A<rsup|\<gamma\>><rsub|m\<nocomma\>n>f<rsub|n\<nocomma\>m>\<delta\><around*|(|\<hbar\>\<omega\>-\<varepsilon\><rsub|m\<nocomma\>n>|)>.>>>>
+  </eqnarray*>
+
+  HopTB.jl calculates
+
+  <\eqnarray*>
+    <tformat|<table|<row|<cell|\<Lambda\><rsup|\<alpha\>\<beta\>\<gamma\>><around*|(|\<omega\>|)>>|<cell|=>|<cell|\<mathpi\>e<rsup|3><big|sum><rsub|n,m><big|int><frac|\<mathd\><rsup|3>\<b-k\>|<around*|(|2\<mathpi\>|)><rsup|3>><around*|(|v<rsup|\<alpha\>><rsub|n>-v<rsub|m><rsup|\<alpha\>>|)>A<rsup|\<beta\>><rsub|n\<nocomma\>m>A<rsup|\<gamma\>><rsub|m\<nocomma\>n>f<rsub|n\<nocomma\>m>\<delta\><around*|(|\<hbar\>\<omega\>-\<varepsilon\><rsub|m\<nocomma\>n>|)>.>>>>
+  </eqnarray*>
+
+  For nonmagnetic materials, <math|\<Lambda\><around*|(|\<omega\>|)>> is
+  purely imaginary. In addition, for nonmagnetic materials,
+  <math|\<Lambda\><rsup|\<alpha\>\<beta\>\<gamma\>><around*|(|\<omega\>|)>=-\<Lambda\><rsup|\<alpha\>\<gamma\>\<beta\>><around*|(|\<omega\>|)>>.
 </body>
 
 <\initial>
diff --git a/src/optics.jl b/src/optics.jl
index 9ff6640..5658fd6 100644
--- a/src/optics.jl
+++ b/src/optics.jl
@@ -286,7 +286,6 @@ The returned matrix berryconnection is
 `berryconnection[n, m, ikpt] = |A[n, m]|^2` at `ikpt`.
 """
 function cltberryconnection(atm::AbstractTBModel, α::Int64, kpts::AbstractMatrix{Float64})
-    nkpts = size(kpts, 2)
     jobs = Vector{Future}()
     berryconnection = zeros((atm.norbits, atm.norbits, 0))
     kptslist = HopTB.Utilities.splitkpts(kpts, nworkers())
@@ -340,7 +339,6 @@ R^{α,β}_{nm} = ∂_αϕ_{mn}^β-A_{mm}^α+A_{nn}^α.
 ```
 """
 function cltshiftvector(atm::AbstractTBModel, α::Int64, β::Int64, kpts::AbstractMatrix{Float64})
-    nkpts = size(kpts, 2)
     jobs = Vector{Future}()
     shiftvector = zeros((atm.norbits, atm.norbits, 0))
     kptslist = HopTB.Utilities.splitkpts(kpts, nworkers())
@@ -598,4 +596,106 @@ function get_Drude_weight(
     return result * bzvol / nks
 end
 
+
+@doc raw"""
+```julia
+get_injection_conductivity_k!(
+    σs::AbstractVector{ComplexF64},
+    tm::AbstractTBModel,
+    α::Int64,
+    β::Int64,
+    γ::Int64,
+    ωs::AbstractVector{Float64},
+    μ::Float64,
+    k::Vector{Float64};
+    ϵ::Float64=0.1
+)
+```
+
+Calculate injection conductivity at `k` point and add the result to `σs`.
+"""
+function get_injection_conductivity_k!(
+    σs::AbstractVector{ComplexF64},
+    tm::AbstractTBModel,
+    α::Int64,
+    β::Int64,
+    γ::Int64,
+    ωs::AbstractVector{Float64},
+    μ::Float64,
+    k::Vector{Float64};
+    ϵ::Float64=0.1
+)
+    nωs = length(ωs)
+    vα = getvelocity(tm, α, k)
+    rβ = getA(tm, β, k)
+    rγ = getA(tm, γ, k)
+    Es = geteig(tm, k).values
+    for n in 1:tm.norbits, m in 1:tm.norbits
+        En = Es[n]
+        Em = Es[m]
+        if abs(En - Em) > only(HopTB.DEGEN_THRESH)
+            fn = (En < μ) ? 1 : 0
+            fm = (Em < μ) ? 1 : 0
+            for iω in 1:nωs
+                ω = ωs[iω]
+                delta = exp(-(Em - En - ω)^2 / ϵ^2) / ϵ / √π
+                σs[iω] += (vα[n, n] - vα[m, m]) * rβ[n, m] * rγ[m, n] * (fn - fm) * delta * 3.082867745843085 # π * e^2 / (ħ * (2π)^3)
+            end
+        end
+    end
+    return nothing
+end
+
+function get_injection_conductivity_k(
+    tm::AbstractTBModel,
+    α::Int64,
+    β::Int64,
+    γ::Int64,
+    ωs::AbstractVector{Float64},
+    μ::Float64,
+    k::Vector{Float64};
+    ϵ::Float64=0.1
+)
+    σs = zeros(ComplexF64, length(ωs))
+    get_injection_conductivity_k!(σs, tm, α, β, γ, ωs, μ, k; ϵ=ϵ)
+    return σs
+end
+
+@doc raw"""
+```julia
+get_injection_conductivity(
+    tm::AbstractTBModel,
+    α::Int64,
+    β::Int64,
+    γ::Int64,
+    ωs::AbstractVector{Float64},
+    μ::Float64,
+    meshsize::Vector{Int64};
+    ϵ::Float64=0.1,
+    batchsize::Int64=1
+)
+```
+
+This function returns injection conductivity in μA * eV / V^2.
+"""
+function get_injection_conductivity(
+    tm::AbstractTBModel,
+    α::Int64,
+    β::Int64,
+    γ::Int64,
+    ωs::AbstractVector{Float64},
+    μ::Float64,
+    meshsize::Vector{Int64};
+    ϵ::Float64=0.1,
+    batchsize::Int64=1
+)
+    nks = prod(meshsize)
+    nωs = length(ωs)
+    σs = [SharedArray{ComplexF64}(nωs) for _ in 1:nprocs()]
+    parallel_do(k -> get_injection_conductivity_k!(σs[myid()], tm, α, β, γ, ωs, μ, k; ϵ=ϵ),
+        UniformMesh(meshsize), batchsize=batchsize)
+    bzvol = abs(det(tm.rlat))
+    return sum(σs) * bzvol / nks
+end
+
 end