-
Notifications
You must be signed in to change notification settings - Fork 83
/
current.jl
21 lines (21 loc) · 935 Bytes
/
current.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
"""
Computes the *probability* (not charge) current, ∑ fn Im(ψn* ∇ψn)
"""
function compute_current(basis::PlaneWaveBasis, ψ, occupation)
@assert all(symop -> length(symop) == 1, basis.ksymops) == 1 # TODO lift this
current = [zeros(eltype(basis), basis.fft_size) for α = 1:3]
for (ik, kpt) in enumerate(basis.kpoints)
for (n, ψnk) in enumerate(eachcol(ψ[ik]))
ψnk_real = G_to_r(basis, kpt, ψnk)
for α = 1:3
dαψnk = [im*(G+kpt.coordinate_cart)[α] for G in G_vectors_cart(kpt)] .* ψnk
dαψnk_real = G_to_r(basis, kpt, dαψnk)
current[α] .+= @. basis.kweights[ik] *
occupation[ik][n] *
imag(conj(ψnk_real) * dαψnk_real)
end
end
end
mpi_sum!.(current, Ref(basis.comm_kpts))
[from_real(basis, current[α]) for α = 1:3]
end