Skip to content

Commit

Permalink
Merge pull request #21 from jarvist/holstein
Browse files Browse the repository at this point in the history 
Holstein code and updated tests
  • Loading branch information
@Neutrino155
Neutrino155 authored Jun 14, 2024
2 parents fc0b345 + 680bc1f commit 637a7b2
Show file tree
Hide file tree
Showing 15 changed files with 317 additions and 482 deletions.
1 change: 0 additions & 1 deletion Project.toml
100644 → 100755
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"

[compat]
JLD = "0.13"
Expand Down
40 changes: 28 additions & 12 deletions src/FrohlichPolaron.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ Outer constructor for the Polaron type. This function evaluates model data for t
julia> polaron(6, 300, 3, 1.0, 3.6, 2.8)
```
"""
function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesses=3.1, w_guesses=2.9, dims=3, kspace=false, verbose=false)
function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesses=false, w_guesses=false, dims=3, kspace=false, cycleguesses=false, verbose=false)

# v_guesses and w_guesses are initial values for v and w (including many v and w parameters).
# These guesses are generally not needed unless instabilities are found in the minimisation and better initial values improve stability.
Expand Down Expand Up @@ -251,13 +251,21 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse
# Extract the ground-state, athermal polaron properties (energy (enthalpy) and variational parameters v and w).
# w is also the frequency of oscillation of the SHM trial system composed of the bare particle and fictitous mass.
# A, B, C are components of the total energy: A is the bare electron energy, B the electron-phonon interaction energy, C is the energy of the harmonic trial system.



if !cycleguesses && v_guesses == false
w_guesses = 2 .+ tanh.((6 .- α) ./ 3)
v_guesses = αeff < 7 ? 3 .+ α ./ 4 : 4 .* α .^2 / 9π .- 3/2 * (2 * log(2) + 0.5772) .- 3/4
end

athermal_energy(v, w) = !kspace ? frohlich_energy(v, w, α, ω; dims = dims[d], mb = mb) : frohlich_energy_k_space(v, w, α, ω; dims = dims[d])

v_gs, w_gs, F_gs, A_gs, B_gs, C_gs = vw_variation(athermal_energy, v_guesses, w_guesses)

# Update the guesses to keep them close-ish to the true solutions during loops over alphas.
v_guesses, w_guesses = v_gs, w_gs
if cycleguesses && (v_guesses == false)
v_guesses, w_guesses = v_gs, w_gs
end

# Store the athermal data.
p["v0"][d, j, :] .= v_gs
Expand Down Expand Up @@ -352,11 +360,19 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse
# w is also the frequency of oscillation of the SHM trial system composed of the bare particle and fictitous mass.
# A, B, C are components of the total energy: A is the bare electron energy, B the electron-phonon interaction energy, C is the energy of the harmonic trial system.

if !cycleguesses && (v_guesses == false)
w_guesses = w_gs
v_guesses = v_gs
end

thermal_energy(v, w) = !kspace ? frohlich_energy(v, w, α, ω, β; dims = dims[d], mb = mb) : frohlich_energy_k_space(v, w, α, ω, β; dims = dims[d])
v, w, F, A, B, C = vw_variation(thermal_energy, v_guesses, w_guesses)

# Update the guesses to keep them close-ish to the true solutions during loops over temperatures.
v_guesses, w_guesses = v, w
if cycleguesses && (v_guesses == false)
println("true")
v_guesses, w_guesses = v, w
end

# Store thermal data.
p["v"][i, d, j, :] .= v
Expand Down Expand Up @@ -637,45 +653,45 @@ end
"""
Single alpha parameter. polaron() expects alpha parameters to be in a Vector.
"""
frohlichpolaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose)
frohlichpolaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, cycleguesses=false, verbose=false) = frohlichpolaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, cycleguesses=cycleguesses, verbose=verbose)

"""
No frequency input.
"""
frohlichpolaron(αrange, Trange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(αrange, Trange, 0; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose)
frohlichpolaron(αrange, Trange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, cycleguesses=false, verbose=false) = frohlichpolaron(αrange, Trange, 0; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose)

"""
No temperature input => 300 K.
"""
frohlichpolaron(αrange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(αrange, 0, 0; ω=ω, ωeff=ωeff, mb=mb,v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose)
frohlichpolaron(αrange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, cycleguesses=false, verbose=false) = frohlichpolaron(αrange, 0, 0; ω=ω, ωeff=ωeff, mb=mb,v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, cycleguesses=cycleguesses, verbose=verbose)

"""
No input => α = 1
"""
frohlichpolaron(; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(1, 0, 0; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose)
frohlichpolaron(; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, cycleguesses=false, verbose=false) = frohlichpolaron(1, 0, 0; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, cycleguesses=cycleguesses, verbose=verbose)

"""
polaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, verbose=false)
Material specific constructors that use material specific parameters to parameterise the polaron.
Material data is inputted through the `Material` type.
Returns all data in either SI units or other common, suitable units otherwise.
"""
function frohlichpolaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false)
function frohlichpolaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, cycleguesses=false, verbose=false)

# Show material data.
if verbose
display(material)
end

# Extract material data from Material type.
phonon_freqs = material.f
phonon_eff_freq = material.feff
phonon_freqs = material.f .* pustrip(1u"THz2π")
phonon_eff_freq = material.feff .* pustrip(1u"THz2π")
mb = material.mb

TΩrange = length(TΩrange) == 1 ? TΩrange .* pustrip(1u"K") : TΩrange .* (pustrip(1u"K"),1)

# Generate polaron data from the arbitrary model constructor.
p = frohlichpolaron(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose)
p = frohlichpolaron(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, cycleguesses=cycleguesses, verbose=verbose)

# Return material-specific, unitful Polaron type.
return p
Expand Down
6 changes: 3 additions & 3 deletions src/FrohlichTheory.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,7 @@ function frohlich_interaction_energy(v, w, α, ωβ...; dims = 3)
integrand(τ) = phonon_propagator(τ, ωβ...) / sqrt(propagator(τ))
upper_limit = length(ωβ) == 1 ? Inf : ωβ[2] / 2
integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit)
return coupling * ball_surface(dims) / (2π)^dims * sqrt(π / 2) * integral
return coupling * ball_surface(dims) / (2π)^dims * sqrt(π / 2) * integral
end

frohlich_interaction_energy(v, w, α::Vector, ω::Vector; dims = 3) = sum(frohlich_interaction_energy(v, w, α[j], ω[j]; dims = dims) for j in eachindex(α))
Expand Down Expand Up @@ -278,7 +278,7 @@ frohlich_structure_factor_k_space(t, v, w, α::Vector, ω::Vector, β; limits =

function frohlich_memory_function(Ω, v, w, α, ωβ...; dims = 3)
structure_factor(t) = frohlich_structure_factor(t, v, w, α, ωβ...; dims = dims)
return polaron_memory_function(Ω, structure_factor)
return polaron_memory_function(Ω, structure_factor, limits = [0, 1e4])
end

frohlich_memory_function(Ω, v, w, α::Vector, ω::Vector; dims = 3) = sum(frohlich_memory_function(Ω, v, w, α[j], ω[j]; dims = dims) for j in eachindex(α))
Expand Down Expand Up @@ -369,7 +369,7 @@ See also [`polaron_mobility`](@ref), [`polaron_complex_conductivity`](@ref)
"""
function inverse_frohlich_mobility(v, w, α, ω, β; dims = 3)
structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω, β; dims = dims)
return abs(imag(polaron_memory_function(structure_factor)))
return abs(imag(polaron_memory_function(structure_factor, limits = [0, 1e4])))
end

"""
Expand Down
56 changes: 30 additions & 26 deletions src/HolsteinPolaron.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,18 +155,22 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=
dprocess += 1
end

# Update the guesses to keep them close-ish to the true solutions during loops over alphas.
if j > 1
v_guesses = reduce_array(p["v0"][d, j-1, :])
w_guesses = reduce_array(p["w0"][d, j-1, :])
end

# Extract the ground-state, athermal polaron properties (energy (enthalpy) and variational parameters v and w).
# w is also the frequency of oscillation of the SHM trial system composed of the bare particle and fictitous mass.
# A, B, C are components of the total energy: A is the bare electron energy, B the electron-phonon interaction energy, C is the energy of the harmonic trial system.
athermal_energy(v, w) = !kspace ? holstein_energy(v, w, α, ω; dims = dims[d]) : holstein_energy_k_space(v, w, α, ω; dims = dims[d])
athermal_energy(v, w) = !kspace ? holstein_energy(v, w, α, J, ω; a = a, dims = dims[d]) : holstein_energy_k_space(v, w, α, J, ω; a = a, dims = dims[d])
v_gs, w_gs, F_gs, A_gs, B_gs, C_gs = vw_variation(athermal_energy, v_guesses, w_guesses)

# Update the guesses to keep them close-ish to the true solutions during loops over alphas.
v_guesses, w_guesses = v_gs, w_gs

# Store the athermal data.
p["v0"][d, j, :] .= v_gs ./ ωeff
p["w0"][d, j, :] .= w_gs ./ ωeff
p["v0"][d, j, :] .= v_gs
p["w0"][d, j, :] .= w_gs
p["F0"][d, j] = F_gs
p["A0"][d, j] = A_gs
p["B0"][d, j] = B_gs
Expand All @@ -177,8 +181,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=
println(io, "\e[K-----------------------------------------------------------------------")
println(io, "\e[K Zero Temperature Information: ")
println(io, "\e[K-----------------------------------------------------------------------")
println(io, "\e[KVariational parameter | v0 = ", v_gs ./ ωeff)
println(io, "\e[KVariational parameter | w0 = ", w_gs ./ ωeff)
println(io, "\e[KVariational parameter | v0 = ", v_gs)
println(io, "\e[KVariational parameter | w0 = ", w_gs)
println(io, "\e[KEnergy | F0 = ", F_gs)
println(io, "\e[KElectron energy | A0 = ", A_gs)
println(io, "\e[KInteraction energy | B0 = ", B_gs)
Expand All @@ -188,7 +192,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=

# Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal
κ_gs = (v_gs .^ 2 .- w_gs .^ 2)
p["κ0"][d, j, :] .= κ_gs ./ ωeff .^2
p["κ0"][d, j, :] .= κ_gs

# Print athermal fictitious spring constant.
if verbose
Expand Down Expand Up @@ -224,7 +228,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=

# Calculate and store polaron radii. Approximates the polaron wavefunction as a Gaussian and relates the size to the standard deviation. Eqn. (2.4) in Schultz1959. Athermal.
R_gs = sqrt.(3 ./ 2 .* M_reduced_gs .* v_gs)
p["R0"][d, j, :] .= R_gs ./ ωeff .^(1/2)
p["R0"][d, j, :] .= R_gs

# Print athermal polaron radius.
if verbose
Expand Down Expand Up @@ -256,24 +260,24 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=
# Calculate thermal polaron properties (energy (Gibbs free energy) and variational parameters v and w).
# w is also the frequency of oscillation of the SHM trial system composed of the bare particle and fictitous mass.
# A, B, C are components of the total energy: A is the bare electron energy, B the electron-phonon interaction energy, C is the energy of the harmonic trial system.
thermal_energy(v, w) = !kspace ? holstein_energy(v, w, α, ω, β; dims = dims[d]) : holstein_energy_k_space(v, w, α, ω, β; dims = dims[d])
thermal_energy(v, w) = !kspace ? holstein_energy(v, w, α, J, ω, β; a = a, dims = dims[d]) : holstein_energy_k_space(v, w, α, J, ω, β; a = a, dims = dims[d])
v, w, F, A, B, C = vw_variation(thermal_energy, v_guesses, w_guesses)

# Update the guesses to keep them close-ish to the true solutions during loops over temperatures.
v_guesses, w_guesses = v, w

# Store thermal data.
p["v"][i, d, j, :] .= v ./ ωeff
p["w"][i, d, j, :] .= w ./ ωeff
p["v"][i, d, j, :] .= v
p["w"][i, d, j, :] .= w
p["F"][i, d, j] = F
p["A"][i, d, j] = A
p["B"][i, d, j] = B
p["C"][i, d, j] = C

# Print thermal data.
if verbose
println(io, "\e[KVariational parameter | v = ", v ./ ωeff)
println(io, "\e[KVariational parameter | w = ", w ./ ωeff)
println(io, "\e[KVariational parameter | v = ", v)
println(io, "\e[KVariational parameter | w = ", w)
println(io, "\e[KFree energy | F = ", F)
println(io, "\e[KElectron energy | A = ", A)
println(io, "\e[KInteraction energy | B = ", B)
Expand All @@ -282,7 +286,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=

# Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal
κ = (v .^ 2 .- w .^ 2)
p["κ"][i, d, j, :] .= κ ./ ωeff .^2
p["κ"][i, d, j, :] .= κ

# Print spring constants.
if verbose
Expand Down Expand Up @@ -318,7 +322,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=

# Calculate and store polaron radii.
R = sqrt.(3 ./ 2 .* v .* M_reduced)
p["R"][i, d, j, :] .= R ./ ωeff .^(1/2)
p["R"][i, d, j, :] .= R

# Print polaron radius.
if verbose
Expand All @@ -332,7 +336,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=
end

# Calculate and store the DC mobiliies.
μ = !kspace ? holstein_mobility(v, w, α, ω, β; dims = dims[d]) : holstein_mobility_k_space(v, w, α, ω, β; dims = dims[d])
μ = !kspace ? holstein_mobility(v, w, α, J, ω, β; a = a, dims = dims[d]) : holstein_mobility_k_space(v, w, α, J, ω, β; a = a, dims = dims[d])
p["μ"][i, d, j] = μ

# Print DC mobilities.
Expand All @@ -344,8 +348,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=
# If zero temperature.
v, w, β = v_gs, w_gs, Inf
p["β"][i, :] .= β
p["v"][i, d, j, :] .= v_gs ./ ωeff
p["w"][i, d, j, :] .= w_gs ./ ωeff
p["v"][i, d, j, :] .= v_gs
p["w"][i, d, j, :] .= w_gs
p["F"][i, d, j] = F_gs
p["A"][i, d, j] = A_gs
p["B"][i, d, j] = B_gs
Expand Down Expand Up @@ -377,7 +381,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=
end

# Calculate and store polaron memory functions (akin to self energy).
χ = !kspace ? holstein_memory_function(Ω, v, w, α, ω, β; dims = dims[d]) : holstein_memory_function_k_space(Ω, v, w, α, ω, β; dims = dims[d])
χ = !kspace ? holstein_memory_function(Ω, v, w, α, J, ω, β; a = a, dims = dims[d]) : holstein_memory_function_k_space(Ω, v, w, α, J, ω, β; a = a, dims = dims[d])
p["χ"][k, i, d, j] = χ

# Print memory function.
Expand Down Expand Up @@ -512,15 +516,15 @@ function holsteinpolaron(material::Material, TΩrange...; v_guesses=3.11, w_gues
end

# Extract material data from Material type.
ω = material.f
ωeff = material.feff
mb = material.mb
J = material.J
ω = material.f .* pustrip(1Unitful.THz2π)
ωeff = material.feff .* pustrip(1Unitful.THz2π)
mb = material.mb .* pustrip(1Unitful.me)
J = material.J .* pustrip(1Unitful.meV)
γ = material.γ
a = material.a
a = material.a .* pustrip(1Unitful.Å)
dims = material.z ./ 2

TΩrange = length(TΩrange) == 1 ? TΩrange .* phustrip(1u"K") : TΩrange .* (phustrip(1u"K"),1)
TΩrange = length(TΩrange) == 1 ? TΩrange .* pustrip(1u"K") : TΩrange .* (pustrip(1u"K"),1)

# Generate Holstein polaron data from the arbitrary model constructor.
p = holsteinpolaron(material.α', TΩrange...; ω=ω, ωeff=ωeff, γ=γ, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose)
Expand Down
Loading

0 comments on commit 637a7b2

Please sign in to comment.