From ad0dc2b1c0fce82df7512af3890892dac4b13e1d Mon Sep 17 00:00:00 2001 From: Bradley Martin Date: Fri, 8 Dec 2023 21:30:09 +0000 Subject: [PATCH] Many changes. Refactoring package. Reduced absolute functions. Made polaron constructors capable of multiple dimensions and provide the ability to exclude temperature or frequency dependent data generation. --- src/EffectiveMass.jl | 21 - src/FeynmanTheory.jl | 569 ------------------ src/{Polaron.jl => FrohlichPolaron.jl} | 484 +++++++++------- src/FrohlichTheory.jl | 625 ++++++++++++++++++++ src/HolsteinPolaron.jl | 438 ++++++++------ src/HolsteinTheory.jl | 526 +++++------------ src/KSpaceTheory.jl | 771 ------------------------- src/LegacyFunctions.jl | 91 +++ src/Material.jl | 1 + src/MemoryFunction.jl | 39 -- src/PolaronFunctions.jl | 223 +++++++ src/PolaronMobility.jl | 59 +- src/PolaronUnits.jl | 64 +- src/ResponseFunctions.jl | 281 --------- src/TrialPolaron.jl | 310 ++++++++++ 15 files changed, 1970 insertions(+), 2532 deletions(-) delete mode 100644 src/EffectiveMass.jl delete mode 100644 src/FeynmanTheory.jl rename src/{Polaron.jl => FrohlichPolaron.jl} (61%) create mode 100644 src/FrohlichTheory.jl delete mode 100644 src/KSpaceTheory.jl create mode 100644 src/LegacyFunctions.jl delete mode 100644 src/MemoryFunction.jl create mode 100644 src/PolaronFunctions.jl delete mode 100644 src/ResponseFunctions.jl create mode 100644 src/TrialPolaron.jl diff --git a/src/EffectiveMass.jl b/src/EffectiveMass.jl deleted file mode 100644 index bc6af2b..0000000 --- a/src/EffectiveMass.jl +++ /dev/null @@ -1,21 +0,0 @@ -# EffectiveMass.jl - -function polaron_effective_mass(v, w, α, ω, β) - - D(t) = 2 * (v^2 - w^2) * sin(v * t / 2) * sin(v * (t - im * ω * β) / 2) / sinh(v * ω * β / 2) / v^3 - im * w^2 * t * (1 + im * t / β / ω) / v^2 - - # # FHIP1962, page 1009, eqn (36). - S(t) = cos(t - im * β * ω / 2) / sinh(ω * β / 2) / D(t)^(3 / 2) - - # FHIP1962, page 1009, eqn (35a). - integrand(t) = imag(S(t)) * t^2 - - integral = quadgk(t -> integrand(t/(1-t))/(1-t)^2, 0, 1-eps(Float64))[1] - - mass = α * integral / (3 * √π) - - return mass -end - -polaron_effective_mass(v, w, α::Vector, ω::Vector, β) = sum(polaron_effective_mass.(v, w, α, ω, β)) - diff --git a/src/FeynmanTheory.jl b/src/FeynmanTheory.jl deleted file mode 100644 index 295c0ef..0000000 --- a/src/FeynmanTheory.jl +++ /dev/null @@ -1,569 +0,0 @@ -# FeynmanTheory.jl - -""" - frohlichalpha(ε_Inf, ε_S, freq, m_eff) - -Calculates the Frohlich alpha parameter, for a given dielectric constant, frequency (f) of phonon in Hertz, and effective mass (in units of the bare electron mass). - -See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. -""" -function frohlichalpha(ϵ_optic, ϵ_static, freq, m_eff) - ω = freq * 2π * 1e12 # frequency to angular velocity - # Note: we need to add a 4*pi factor to the permitivity of freespace. - # This gives numeric agreement with literature values. This is required as - # the contemporary 1950s and 1960s literature implicitly used atomic units, - # where the electric constant ^-1 has this factor baked in, k_e=1/(4πϵ_0). - α = 1 / 2 / (4 * π * ϵ_0) * # Units: m/F - (1 / ϵ_optic - 1 / ϵ_static) * # Units: none - (eV^2 / (ħ * ω)) * # Units: F - sqrt(2 * m_eff * me * ω / ħ) # Units: 1/m - return α -end - -# Athermal (Feynman 1955) model. -# Set up equations for the polaron free energy, which we will variationally improve upon. - -D_imag(τ, v, w, ω, β) = (v^2 - w^2) / v^3 * (1 - exp(-v * τ)) * (1 - exp(-v * (β * ω - τ))) / (1 - exp(-v * β * ω)) + w^2 / v^2 * τ * (1 - τ / β / ω) + eps(Float64) - -D_imag(τ, v, w) = w^2 * τ / v^2 + (v^2 - w^2) / v^3 * (1 - exp(-v * τ)) + eps(Float64) - -""" - B(v, w, α, ω) - -Integral of Eqn. (31) in Feynman 1955. Part of the overall ground-state energy expression. - -See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. -""" -B(v, w, α, ω; dims = 3) = frohlich_coupling(1, α, ω; dims = dims) * ball_surface(dims) / (2π)^dims * sqrt(π / 2) * quadgk(τ -> phonon_propagator(τ / ω, ω) / sqrt(polaron_propagator(τ, v, w)), 0, Inf)[1] / ω - -B(v, w, α::Vector, ω::Vector; dims = 3) = sum(B.(v, w, α, ω; dims = dims)) - -A(v, w, ω) = -3 * (v - w) / 2 * ω - -A(v, w, ω::Vector) = sum(A.(v, w, ω)) / length(ω) - -C(v, w, ω) = (3 / (4 * v)) * (v^2 - w^2) * ω - -C(v, w, ω::Vector) = sum(C.(v, w, ω)) / length(ω) - -# Hellwarth et al. 1999 PRB - Part IV; T-dep of the Feynman variation parameter - -# In Julia we have 'Multiple dispatch', so let's just construct the free -# energies (temperature-dependent) with the same name as before, but with the thermodynamic beta where required. - -# Define Osaka's free-energies (Hellwarth1999 version) as Julia functions. - -""" - A(v, w, ω, β) - -Hellwarth's A expression from Eqn. (62b) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression. - -See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. -""" -A(v, w, ω, β) = 3 / β * (log(v / w) - 1 / 2 * log(2π * β * ω -) - (v - w) * ω * β / 2 - log(1 - exp(-v * β * ω)) + log(1 - exp(-w * β * ω))) - -A(v, w, ω::Vector, β) = sum(A.(v, w, ω, β)) / length(ω) - -""" - B(v, w, α, ω, β) - -Hellwarth's B expression from Eqn. (62c) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression. - -See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. -""" -B(v, w, α, ω, β; dims=3) = frohlich_coupling(1, α, ω; dims = dims) * ball_surface(dims) / (2π)^dims * sqrt(π / 2) * quadgk(τ -> phonon_propagator(τ / ω, ω, β) / sqrt(polaron_propagator(τ, v, w, β * ω)), 0, β * ω / 2)[1] / ω - -B(v, w, α::Vector, ω::Vector, β) = sum(B.(v, w, α, ω, β)) - -""" - C(v, w, ω, β) - -Hellwarth's C expression from Eqn. (62e) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression. - -See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. -""" -C(v, w, ω, β) = 3 / 4 * (v^2 - w^2) * ω / v * (coth(v * β * ω / 2) - 2 / (v * β * ω)) - -C(v, w, ω::Vector, β) = sum(C.(v, w, ω, β)) / length(ω) - -# Extending the Feynman theory to multiple phonon branches - -# Partial dielectric electron-phonon coupling parameter, decomposed into ionic dielectric contributions from each phonon mode of the material. - -""" - ϵ_ionic_mode(phonon_mode_freq, ir_activity, volume) - -Calculate the ionic contribution to the dielectric function for a given phonon mode. - -# Arguments -- `phonon_mode_freq::Float64`: is the frequency of the mode in THz. -- `ir_activity::Float64`: is the infra-red activity of the mode in e²amu⁻¹. -- `volume::Float64`: is the volume of the unit cell of the material in m³. -""" -function ϵ_ionic_mode(phonon_mode_freq, ir_activity, volume) # single ionic mode - - # Angular phonon frequency for the phonon mode (rad Hz) - ω_j = phonon_mode_freq * 2π * 1e12 - - # Dielectric contribution from a single ionic phonon mode - ϵ_mode = eV^2 * ir_activity / (3 * volume * ω_j^2 * amu) - - # Normalise ionic dielectric contribution with 1 / (4π ϵ_0) (NB: the 4π has been pre-cancelled) - return ϵ_mode / ϵ_0 -end - -""" - ϵ_total(freqs_and_ir_activity, volume) - -Calculate the total ionic contribution to the dielectric function from all phonon modes. - -# Arguments -- `freqs_and_ir_activity::Matrix{Float64}`: is a matrix containeing the phonon mode frequencies (in THz) in the first column and the infra-red activities (in e²amu⁻¹) in the second column. -- `volume::Float64`: is the volume of the unit cell of the material in m^3. -""" -function ϵ_total(freqs_and_ir_activity, volume) # total ionic contribution to dielectric - - # Extract phonon frequencies (THz) - phonon_freqs = freqs_and_ir_activity[:, 1] - - # Extra infra-red activities (e^2 amu^-1) - ir_activity = freqs_and_ir_activity[:, 2] - - # Sum over all ionic contribution from each phonon mode - total_ionic = 0 - - for t in eachindex(phonon_freqs) - total_ionic += ϵ_ionic_mode(phonon_freqs[t], ir_activity[t], volume) - end - - return total_ionic -end - -""" - effective_freqs(freqs_and_ir_activity, num_var_params) - -Generates a matrix of effective phonon modes with frequencies and infra-red activities derived from a larger matrix using the Principal Component Analysis (PCA) method. - -# Arguments -- `freqs_and_ir_activity::Matrix{Float64}`: is a matrix containing the phonon mode frequencies (in THz) in the first column and the infra-red activities (in e²amu⁻¹) in the second column. -- `num_var_params::Integer`: is the number of effective modes required (which needs to be less than the number of modes in `freqs_and_ir_activity``). - -*** POSSIBLY REDUNDANT *** -""" -function effective_freqs(freqs_and_ir_activity, num_var_params) # PCA Algorithm - - # Check that the number of effective modes is less than the number of actual phonon modes. - if num_var_params >= size(freqs_and_ir_activity)[1] - - println("The number of effective phonon modes has to be less than the total number of phonon modes.") - - else - - # Centralise data by subtracting the columnwise mean - standardized_matrix = freqs_and_ir_activity' .- mean(freqs_and_ir_activity', dims=2) - - # Calculate the covariance matrix S' * S. Matrix size is (n - 1) x (n - 1) for number of params (here n = 2) - covariance_matrix = standardized_matrix' * standardized_matrix - - # Extract eigenvectors of the covariance matrix - eigenvectors = eigvecs(covariance_matrix) - - # Project the original data along the covariance matrix eigenvectors and undo the centralisation - reduced_matrix = standardized_matrix[:, 1:num_var_params] * eigenvectors[1:num_var_params, 1:num_var_params] * - eigenvectors[1:num_var_params, 1:num_var_params]' .+ mean(freqs_and_ir_activity', dims=2) - - # Resultant matrix is positive definite and transposed. - return abs.(reduced_matrix') - end -end - -""" - frohlichalpha(ϵ_optic, ϵ_ionic, ϵ_total, phonon_mode_freq, m_eff) - -Calculates the partial dielectric electron-phonon coupling parameter for a given longitudinal optical phonon mode. - -This decomposes the original Frohlich alpha coupling parameter (defined for a single phonon branch) into contributions from multiple phonon -branches. - -# Arguments -- `ϵ_optic::Float64`: is the optical dielectric constant of the material. -- `ϵ_ionic::Float64`: is the ionic dielectric contribution from the phonon mode. -- `ϵ_total::Float64`: is the total ionic dielectric contribution from all phonon modes of the material. -- `phonon_mode_freq::Float64`: is the frequency of the phonon mode (THz). -- `m_eff::Float64` is the band mass of the electron (in units of electron mass m_e). -""" -function frohlichalpha(ϵ_optic, ϵ_ionic, ϵ_total, phonon_mode_freq, m_eff) - - # The Rydberg energy unit - Ry = eV^4 * me / (2 * ħ^2) - - # Angular phonon frequency for the phonon mode (rad Hz). - ω = phonon_mode_freq * 2π * 1e12 - - # The static dielectric constant. Calculated here instead of inputted so that ionic modes are properly normalised. - ϵ_static = ϵ_total + ϵ_optic - - # The contribution to the electron-phonon parameter from the currrent phonon mode. 1 / (4π ϵ_0) is the dielectric normalisation. - α_j = (m_eff * Ry / (ħ * ω))^(1 / 2) * ϵ_ionic / (4π * ϵ_0) / (ϵ_optic * ϵ_static) - - return α_j -end - -# Extending the Feynman theory to multiple variational parameters -# Multiple Parameter Polaron Free Energy -# Calculate the polaron free energy, generalised from Osaka's expression to the case where multiple phonon modes are present in the material. - -""" - κ_i(i, v, w) - -Calculates the spring-constant coupling the electron to the 'ith' fictitious mass that approximates the exact electron-phonon interaction with a harmonic coupling to a massive fictitious particle. - -Required for calculating the polaron free energy. - -Note: Not to be confused with the number of physical phonon branches; many phonon branches could be approximated with one or two etc. fictitious masses for example. The number of fictitious mass does not necessarily need to match the number of phonon branches. - -# Arguments -- `i::Integer`: enumerates the current fictitious mass. -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -""" -function κ_i(i, v::Vector, w::Vector) - κ = v[i]^2 - w[i]^2 - κ *= prod(j != i ? (v[j]^2 - w[i]^2) / (w[j]^2 - w[i]^2) : 1 for j in eachindex(v)) - return κ -end - -""" - h_i(i, v, w) - -Calculates the normal-mode (the eigenmodes) frequency of the coupling between the electron and the `ith' fictitious mass that approximates the exact electron-phonon interaction with a harmonic coupling to a massive fictitious particle. - -Required for calculating the polaron free energy. - -Note: Not to be confused with the number of physical phonon branches; many phonon branches could be approximated with one or two etc. fictitious masses for example. The number of fictitious mass does not necessarily need to match the number of phonon branches. - -# Arguments -- `i::Integer`: enumerates the current fictitious mass. -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -""" -function h_i(i, v::Vector, w::Vector) - h = v[i]^2 - w[i]^2 - h *= prod(j != i ? (w[j]^2 - v[i]^2) / (v[j]^2 - v[i]^2) : 1 for j in eachindex(v)) - return h -end - -""" - C_ij(i, j, v, w) - -Calculates the element to the coupling matrix C_ij (a generalisation of Feynman's `C` coupling variational parameter in Feynman 1955) between the electron and the `ith' and `jth' fictitious masses that approximates the exact electron-phonon interaction with a harmonic coupling to a massive fictitious particle. - -Required for calculating the polaron free energy. - -Note: Not to be confused with the number of physical phonon branches; many phonon branches could be approximated with one or two etc. fictitious masses for example. The number of fictitious mass does not necessarily need to match the number of phonon branches. - -# Arguments -- `i::Integer, j::Integer`: enumerate the current fictitious masses under focus (also the index of the element in the coupling matrix C) -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. - -See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. -""" -function C_ij(i, j, v::Vector, w::Vector) - C = w[i] * κ_i(i, v, w) * h_i(j, v, w) / (4 * (v[j]^2 - w[i]^2)) - return C -end - -""" - D(τ, v, w, β) - -Calculates the recoil function (a generalisation of D(u) in Eqn. (35c) in FHIP 1962) that approximates the exact influence (recoil effects) of the phonon bath on the electron with the influence of the fictitious masses attached by springs to the electron. It appears in the exponent of the intermediate scattering function. - -# Arguments -- `τ::Float64`: is the imaginary time variable. -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. - -See FHIP 1962: https://doi.org/10.1103/PhysRev.127.1004. -""" -function D(τ, v::Vector, w::Vector, ω, β) - return τ * (1 - τ / ω / β) + sum((h_i(i, v, w) / v[i]^2) * ((1 + exp(-v[i] * ω * β) - exp(-v[i] * τ) - exp(v[i] * (τ - ω * β))) / (v[i] * (1 - exp(-v[i] * ω * β))) - τ * (1 - τ / ω / β)) for i in eachindex(v)) -end - -""" - D(τ, v, w) - -Calculates the recoil function at zero-temperature. - -# Arguments -- `τ::Float64`: is the imaginary time variable. -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -""" -function D(τ, v::Vector, w::Vector) - return τ + sum((h_i(i, v, w) / v[i]^2) * ((1 - exp(-v[i] * τ)) / v[i] - τ) for i in eachindex(v)) -end - -""" - B(v, w, α, β) - -Generalisation of the B function from Eqn. (62c) in Hellwarth et al. 1999. This is the expected value of the exact action taken w.r.t trial action, given for the 'jth' phonon mode. - -Required for calculating the polaron free energy. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for the 'jth' phonon mode. -- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. - -See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. -""" -function B(v::Vector, w::Vector, α, ω, β) - B_integrand(τ) = cosh(τ - ω * β / 2) / sqrt(abs(D(τ, v, w, ω, β))) - return α / (√π * sinh(ω * β / 2)) * quadgk(τ -> B_integrand(τ * ω * β / 2), 0.0, 1.0)[1] * ω^2 * β / 2 -end - -B(v::Vector, w::Vector, α::Vector, ω::Vector, β) = sum(B(v, w, α[i], ω[i], β) for i in eachindex(ω)) - -""" - B(v, w, α; rtol = 1e-3) - -Calculates `B(v, w, α, β)` but at zero-temperature, `β = Inf`. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for the 'jth' phonon mode. -""" -function B(v::Vector, w::Vector, α, ω) - B_integrand(τ) = exp(-abs(τ)) / sqrt(abs(D(abs(τ), v, w))) - return α / √π * quadgk(τ -> B_integrand(τ), 0, Inf)[1] * ω -end - -B(v::Vector, w::Vector, α::Vector, ω::Vector) = sum(B(v, w, α[i], ω[i]) for i in eachindex(ω)) - -""" - C(v, w, β) - -Generalisation of the C function from Eqn. (62e) in Hellwarth et al. 1999. This is the expected value of the trial action taken w.r.t trial action. - -Required for calculating the polaron free energy. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. - - -See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. -""" -function C(v::Vector, w::Vector, ω, β) - # Sum over the contributions from each fictitious mass. - s = sum(C_ij(i, j, v, w) / (v[j] * w[i]) * (coth(ω * β * v[j] / 2) - 2 / (ω * β * v[j])) for i in eachindex(v), j in eachindex(w)) - - # Divide by the number of phonon modes to give an average contribution per phonon mode. - return 3 * s * ω -end - -C(v::Vector, w::Vector, ω::Vector, β) = sum(C(v, w, ω[i], β) for i in eachindex(ω)) / length(ω) - -""" - C(v, w) - -Calculates `C(v, w, β)` but at zero-temperature, `β = Inf`. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -""" -function C(v::Vector, w::Vector, ω) - # Sum over the contributions from each fictitious mass. - s = sum(C_ij(i, j, v, w) / (v[j] * w[i]) for i in eachindex(v), j in eachindex(w)) - # Divide by the number of phonon modes to give an average contribution per phonon mode. - return 3 * s * ω -end - -C(v::Vector, w::Vector, ω::Vector) = sum(C(v, w, ω[i]) for i in eachindex(ω)) / length(ω) - -""" - A(v, w, β) - -Generalisation of the A function from Eqn. (62b) in Hellwarth et al. 1999. This is the Helmholtz free energy of the trial model. - -Required for calculating the polaron free energy. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. - -See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. -""" -function A(v::Vector, w::Vector, ω, β) - # Sum over the contributions from each fictitious mass. - s = -log(2π * ω * β) / 2 + sum(v[i] == w[i] ? 0 : - log(v[i]) - log(w[i]) - ω * β / 2 * (v[i] - w[i]) - log(1 - exp(-v[i] * ω * β)) + log(1 - exp(-w[i] * ω * β)) - for i in eachindex(v)) - # Divide by the number of phonon modes to give an average contribution per phonon mode. - 3 / β * s * ω -end - -A(v::Vector, w::Vector, ω::Vector, β) = sum(A(v, w, ω[i], β) for i in eachindex(ω)) / length(ω) - -""" - A(v, w, n) - -Calculates `A(v, w, β)` but at zero-temperature, `β = Inf`. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -""" -function A(v::Vector, w::Vector, ω) - s = sum(v .- w) - return -3 * s / 2 * ω -end - -A(v::Vector, w::Vector, ω::Vector) = sum(A(v, w, ω[i]) for i in eachindex(ω)) / length(ω) - -""" - F(v, w, α, ω, β) - -Calculates the Helmholtz free energy of the polaron for a material with multiple phonon branches. - -This generalises the Osaka 1959 (below Eqn. (22)) and Hellwarth. et al 1999 (Eqn. (62a)) free energy expressions. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for the 'jth' phonon mode. -- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz). -- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. - -See Osaka, Y. (1959): https://doi.org/10.1143/ptp.22.437 and Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. -""" -function frohlich_energy(v, w, α, ω, β; dims = 3) - - # Insurance to avoid breaking the integrals with Infinite beta. - if β == Inf - return frohlich_energy(v, w, α, ω; dims = dims) - end - Ar = A(v, w, ω, β) * dims / 3 - Br = B(v, w, α, ω, β) - Cr = C(v, w, ω, β) * dims / 3 - # Free energy in units of meV - return -(Ar + Br + Cr), Ar, Br, Cr -end - -""" - F(v, w, α, ω) - -Calculates the zero-temperature ground-state energy of the polaron for a material with multiple phonon branches. Similar to `F(v, w, α, ω, β)` but with `β = Inf`. Generalises Eqn. (33) in Feynman 1955. - -# Arguments -- `v::Vector{Float64}`: is a vector of the v variational parameters. -- `w::Vector{Float64}`: is a vector of the w variational parameters. -- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for the 'jth' phonon mode. -- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz). - -See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. -""" -function frohlich_energy(v, w, α, ω; dims = 3) - Ar = A(v, w, ω) * dims / 3 - Br = B(v, w, α, ω; dims = dims) - Cr = C(v, w, ω) * dims / 3 - # Free energy in units of meV - return -(Ar + Br + Cr), Ar, Br, Cr -end - -""" - feynmanvw(v, w, αωβ...; upper_limit = Inf64) - -Minimises the multiple phonon mode free energy function for a set of vₚ and wₚ variational parameters. The variational parameters follow the inequality: v₁ > w₁ > v₂ > w₂ > ... > vₙ > wₙ. Generalises `feynmanvw` to multiple variational parameters. - -# Arguments -- `v::Vector{Float64}`: vector of initial v parameters. -- `w::Vector{Float64}`: vector of initial w parameters. -- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for one or many phonon modes. -- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz) for one or many phonon modes. -- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) for one or many phonon modes. - -See also [`F`](@ref). -""" -function feynmanvw(v::Vector, w::Vector, αωβ...; upper_limit=1e6) - - if length(v) != length(w) - return error("The number of variational parameters v & w must be equal.") - end - - N_params = length(v) - - Δv = v .- w - initial = vcat(Δv .+ repeat([eps(Float64)], N_params), w) - - # Limits of the optimisation. - lower = fill(0.0, 2 * N_params) - upper = fill(upper_limit, 2 * N_params) - - # The multiple phonon mode free energy function to minimise. - f(x) = frohlich_energy([x[2*n-1] for n in 1:N_params] .+ [x[2*n] for n in 1:N_params], [x[2*n] for n in 1:N_params], αωβ...)[1] - - # Use Optim to optimise the free energy function w.r.t the set of v and w parameters. - solution = Optim.optimize( - Optim.OnceDifferentiable(f, initial; autodiff=:forward), - lower, - upper, - initial, - Fminbox(BFGS()) - ) - - # Extract the v and w parameters that minimised the free energy. - var_params = Optim.minimizer(solution) - - # Separate the v and w parameters into one-dimensional arrays (vectors). - Δv = [var_params[2*n-1] for n in 1:N_params] - w = [var_params[2*n] for n in 1:N_params] - E, A, B, C = frohlich_energy(Δv .+ w, w, αωβ...) - - # if Optim.converged(solution) == false - # @warn "Failed to converge. v = $(Δv .+ w), w = $w, E = $E" - # end - - # Return the variational parameters that minimised the free energy. - return Δv .+ w, w, E, A, B, C -end - -function feynmanvw(v::Real, w::Real, αωβ...; upper = [Inf, Inf], lower = [0, 0]) - - Δv = v .- w - initial = [Δv + eps(Float64), w] - - # The multiple phonon mode free energy function to minimise. - f(x) = frohlich_energy(x[1] .+ x[2], x[2], αωβ...)[1] - - # Use Optim to optimise the free energy function w.r.t the set of v and w parameters. - solution = Optim.optimize( - Optim.OnceDifferentiable(f, initial; autodiff=:forward), - lower, - upper, - initial, - Fminbox(BFGS()) - ) - - # Extract the v and w parameters that minimised the free energy. - Δv, w = Optim.minimizer(solution) - E, A, B, C = frohlich_energy(Δv .+ w, w, αωβ...) - - if Optim.converged(solution) == false - @warn "Failed to converge. v = $(Δv .+ w), w = $w, E = $E" - end - - # Return the variational parameters that minimised the free energy. - return Δv .+ w, w, E, A, B, C -end - -feynmanvw(αωβ...) = feynmanvw(3.4, 2.6, αωβ...) - diff --git a/src/Polaron.jl b/src/FrohlichPolaron.jl similarity index 61% rename from src/Polaron.jl rename to src/FrohlichPolaron.jl index 6020576..d76eeee 100644 --- a/src/Polaron.jl +++ b/src/FrohlichPolaron.jl @@ -1,11 +1,12 @@ -# Polaron.jl +# FrohlichPolaron.jl +# A do-it-all function for calculating properties of the Frohlich polaron. Data is stored as a struct. """ Polaron(x...) Type for storing the polaron information, `x...`. """ -mutable struct Polaron +mutable struct FrohlichPolaron α # Fröhlich coupling. αeff # Effective Fröhlich coupling summed for multiple modes. T # Temperature. @@ -54,7 +55,7 @@ mutable struct Polaron z # Complex impedence. σ # Complex conductivity. - function Polaron(x...) + function FrohlichPolaron(x...) new(reduce_array.(x)...) end end @@ -69,7 +70,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 polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.01, w_guesses=2.99, dims=3,verbose=false) +function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.1, w_guesses=2.9, dims=3, kspace=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. @@ -79,6 +80,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses num_T = length(Trange) num_Ω = length(Ωrange) num_ω = length(ω) + num_d = length(dims) ω = reduce_array(ω) @@ -96,45 +98,45 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses "β" => Matrix{Float64}(undef, num_T, num_ω), # Betas. "Ω" => Ωrange, # Photon frequencies. "d" => dims, # Number of spatial dimensions. - "v0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal v params. - "w0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal w params. - "F0" => Vector{Float64}(undef, num_α), # Athermal energies. - "A0" => Vector{Float64}(undef, num_α), # Athermal A parameter. - "B0" => Vector{Float64}(undef, num_α), # Athermal B parameter. - "C0" => Vector{Float64}(undef, num_α), # Athermal C parameter. - "Fs" => Vector{Float64}(undef, num_α), # Small alpha (α→0) approximate energy. - "Fl" => Vector{Float64}(undef, num_α), # Large alpha (α→∞) approximate energy. - "κ0" => Matrix{Float64}(undef, num_α, num_vw), # Fictitious spring constant. - "M0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal fictitious mass. + "v0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal v params. + "w0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal w params. + "F0" => Matrix{Float64}(undef, num_d, num_α), # Athermal energies. + "A0" => Matrix{Float64}(undef, num_d, num_α), # Athermal A parameter. + "B0" => Matrix{Float64}(undef, num_d, num_α), # Athermal B parameter. + "C0" => Matrix{Float64}(undef, num_d, num_α), # Athermal C parameter. + "Fs" => Matrix{Float64}(undef, num_d, num_α), # Small alpha (α→0) approximate energy. + "Fl" => Matrix{Float64}(undef, num_d, num_α), # Large alpha (α→∞) approximate energy. + "κ0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Fictitious spring constant. + "M0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal fictitious mass. "Ms" => Vector{Float64}(undef, num_α), # Small alpha (α→0) approximate fictitious mass. "Ml" => Vector{Float64}(undef, num_α), # Large alpha (α→∞) approximate fictitious mass. - "M0a" => Matrix{Float64}(undef, num_α, num_vw), # Athermal asymptotic approximate fictitious mass (v0/w0). - "M0r" => Matrix{Float64}(undef, num_α, num_vw), # Athermal reduced mass of particle + fictitious particle system. - "R0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal polaron radius (s. d. of a Gaussian wavefunction). + "M0a" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal asymptotic approximate fictitious mass (v0/w0). + "M0r" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal reduced mass of particle + fictitious particle system. + "R0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal polaron radius (s. d. of a Gaussian wavefunction). "Rs" => Vector{Float64}(undef, num_α), # Small alpha (α→0) approximate polaron radius. "Rl" => Vector{Float64}(undef, num_α), # Large alpha (α→∞) approximate polaron radius. "ΩFC" => Vector{Float64}(undef, num_α), # Large alpha (α → ∞)approximate Franck-Condon peak frequency. - "v" => Array{Float64,3}(undef, num_T, num_α, num_vw), # v params. - "w" => Array{Float64,3}(undef, num_T, num_α, num_vw), # w params. - "F" => Matrix{Float64}(undef, num_T, num_α), # Energies. - "A" => Matrix{Float64}(undef, num_T, num_α), # A parameter. - "B" => Matrix{Float64}(undef, num_T, num_α), # B parameter. - "C" => Matrix{Float64}(undef, num_T, num_α), # C parameter. - "κ" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Spring constants. - "M" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Fictitious masses. - "Ma" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Thermal asymptotic approximate fictitious mass (v/w). - "Mr" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Thermal reduced mass of particle + fictitious particle system. - "R" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Polaron radii. - "μ" => Matrix{Float64}(undef, num_T, num_α), # Mobilities. + "v" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # v params. + "w" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # w params. + "F" => Array{Float64,3}(undef, num_T, num_d, num_α), # Energies. + "A" => Array{Float64,3}(undef, num_T, num_d, num_α), # A parameter. + "B" => Array{Float64,3}(undef, num_T, num_d, num_α), # B parameter. + "C" => Array{Float64,3}(undef, num_T, num_d, num_α), # C parameter. + "κ" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Spring constants. + "M" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Fictitious masses. + "Ma" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Thermal asymptotic approximate fictitious mass (v/w). + "Mr" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Thermal reduced mass of particle + fictitious particle system. + "R" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Polaron radii. + "μ" => Array{Float64,3}(undef, num_T, num_d, num_α), # Mobilities. "μFHIP" => Matrix{Float64}(undef, num_T, num_α), # FHIP low-temperature DC mobiltiy approximation. "μD" => Matrix{Float64}(undef, num_T, num_α), # Kadanoff DC mobility, low-temperature approximation using Devreese2016. "μK" => Matrix{Float64}(undef, num_T, num_α), # Kadanoff DC mobility, low-temperature approximation. "τ" => Matrix{Float64}(undef, num_T, num_α), # relaxation time from Kadanoff Boltzmann transport equation. "μH" => Matrix{Float64}(undef, num_T, num_α), # Hellwarth DC mobility. "μH0" => Matrix{Float64}(undef, num_T, num_α), # Hellwarth DC mobility with b=0. - "χ" => Array{ComplexF64,3}(undef, num_Ω, num_T, num_α), # 'Memory function' or self energy. - "z" => Array{ComplexF64,3}(undef, num_Ω, num_T, num_α), # Complex impedences. - "σ" => Array{ComplexF64,3}(undef, num_Ω, num_T, num_α), # Complex conductivities. + "χ" => Array{ComplexF64,4}(undef, num_Ω, num_T, num_d, num_α), # 'Memory function' or self energy. + "z" => Array{ComplexF64,4}(undef, num_Ω, num_T, num_d, num_α), # Complex impedences. + "σ" => Array{ComplexF64,4}(undef, num_Ω, num_T, num_d, num_α), # Complex conductivities. ) # Print the phonon frequencies. @@ -147,6 +149,8 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses for j in axes(αrange, 1) # αrange loop. + dprocess = 1 + # Reduce αrange to either a Vector of values per phonon mode, or a single scalar. α = reduce_array(αrange[j, :]) @@ -169,7 +173,8 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses println(io, "\e[KPhonon frequencies | ωeff = ", ωeff, " | ω = ", join(round.(first(ω, 2), digits=2), ", ")..., " ... ", join(round.(last(ω, 2), digits=2), ", ")...) println(io, "\e[KFröhlich coupling | αeff = ", αeff, " | α = ", join(round.(first(α, 2), digits=3), ", ")..., " ... ", join(round.(last(α, 2), digits=3), ", ")...) end - println(io, "\e[KNumber of dimensions | d = ", dims) + + αprocess += 1 # Increment αrange iteration. end # Small alpha (α → 0) approximate energy. @@ -235,22 +240,30 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses println(io, "\e[KLarge α→∞ FC peak freq. | ΩFC = ", Ω_FC) end + for d in 1:num_d + if verbose + println(io, "\e[KNumber of dimensions", " [", dprocess, " / ", num_d, "]", " | d = ", dims[d]) + dprocess += 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) = frohlich_energy(v, w, α, ω; dims = dims) + + athermal_energy(v, w) = !kspace ? frohlich_energy(v, w, α, ω; dims = dims[d]) : 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 # Store the athermal data. - p["v0"][j, :] .= v_gs - p["w0"][j, :] .= w_gs - p["F0"][j] = F_gs - p["A0"][j] = A_gs - p["B0"][j] = B_gs - p["C0"][j] = C_gs + 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 + p["C0"][d, j] = C_gs # Print athermal variational parameter and energy data. if verbose @@ -264,12 +277,11 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses println(io, "\e[KInteraction energy | B0 = ", B_gs) println(io, "\e[KTrial energy | C0 = ", C_gs) Tprocess = 1 # Counter for Trange. - αprocess += 1 # Increment αrange iteration. end # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal κ_gs = (v_gs .^ 2 .- w_gs .^ 2) - p["κ0"][j, :] .= κ_gs + p["κ0"][d, j, :] .= κ_gs # Print athermal fictitious spring constant. if verbose @@ -278,7 +290,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses # Calculate and store fictitious masses. Athermal. M_gs = κ_gs ./ w_gs .^ 2 - p["M0"][j, :] .= M_gs + p["M0"][d, j, :] .= M_gs # Print athermal fictitious mass. if verbose @@ -287,7 +299,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Athermal M_asymp_gs = v_gs ./ w_gs - p["M0a"][j, :] .= M_asymp_gs + p["M0a"][d, j, :] .= M_asymp_gs # Print athermal asymptotic fictitious mass. if verbose @@ -296,7 +308,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. M_reduced_gs = (v_gs .^2 .- w_gs .^2) / v_gs .^ 2 - p["M0r"][j, :] .= M_reduced_gs + p["M0r"][d, j, :] .= M_reduced_gs # Print athermal reduced mass. if verbose @@ -305,7 +317,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses # 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 .* v_gs ./ (v_gs .^ 2 .- w_gs .^ 2) .^ 2) - p["R0"][j, :] .= R_gs + p["R0"][d, j, :] .= R_gs # Print athermal polaron radius. if verbose @@ -315,152 +327,181 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses for i in eachindex(Trange) # Temperatures loop. T = Trange[i] - # Print temperature. - if verbose - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[K Finite Temperature Information: [$(Tprocess[]) / $(num_T) ($(round(Tprocess[] / (num_T) * 100, digits=1)) %)]") - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[KTemperatures | T = ", T) - end + if !iszero(T) - # Calculate the reduced (unitless) thermodynamic betas for each phonon mode. - β = 1 ./ T .* β0 - p["β"][i, :] .= β + # Print temperature. + if verbose + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[K Finite Temperature Information: [$(Tprocess[]) / $(num_T) ($(round(Tprocess[] / (num_T) * 100, digits=1)) %)]") + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[KTemperatures | T = ", T) + end - # Print thermodynamic betas. - if verbose - println(io, "\e[KReduced thermodynamic | β = ", β) - end + # Calculate the reduced (unitless) thermodynamic betas for each phonon mode. + β = 1 ./ T .* β0 + p["β"][i, :] .= β - # 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) = frohlich_energy(v, w, α, ω, β; dims = dims) - v, w, F, A, B, C = vw_variation(thermal_energy, v_guesses, w_guesses) + # Print thermodynamic betas. + if verbose + println(io, "\e[KReduced thermodynamic | β = ", β) + end - # Update the guesses to keep them close-ish to the true solutions during loops over temperatures. - v_guesses, w_guesses = v, w + # 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. - # Store thermal data. - p["v"][i, j, :] .= v - p["w"][i, j, :] .= w - p["F"][i, j] = F - p["A"][i, j] = A - p["B"][i, j] = B - p["C"][i, j] = C + thermal_energy(v, w) = !kspace ? frohlich_energy(v, w, α, ω, β; dims = dims[d]) : frohlich_energy_k_space(v, w, α, ω, β; dims = dims[d]) + v, w, F, A, B, C = vw_variation(thermal_energy, v_guesses, w_guesses) - # Print thermal data. - if verbose - 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) - println(io, "\e[KTrial energy | C = ", C) - end + # Update the guesses to keep them close-ish to the true solutions during loops over temperatures. + v_guesses, w_guesses = v, w - # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal - κ = (v .^ 2 .- w .^ 2) - p["κ"][i, j, :] .= κ + # Store thermal data. + 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 spring constants. - if verbose - println(io, "\e[KFictitious spring constant | κ = ", κ) - end - - # Calculate and store fictitious masses. Thermal. - M = κ ./ w .^ 2 - p["M"][i, j, :] .= M + # Print thermal data. + if verbose + 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) + println(io, "\e[KTrial energy | C = ", C) + end - # Print masses. - if verbose - println(io, "\e[KFictitious mass | M = ", M) - end + # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal + κ = (v .^ 2 .- w .^ 2) + p["κ"][i, d, j, :] .= κ - # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Thermal - M_asymp = v ./ w - p["Ma"][i, j, :] .= M_asymp + # Print spring constants. + if verbose + println(io, "\e[KFictitious spring constant | κ = ", κ) + end - # Print asymptotic masses. - if verbose - println(io, "\e[KFictitious mass (asymptotic) | Ma = ", M_asymp) - end + # Calculate and store fictitious masses. Thermal. + M = κ ./ w .^ 2 + p["M"][i, d, j, :] .= M - # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. - M_reduced = (v .^2 .- w .^2) / v .^ 2 - p["Mr"][i, j, :] .= M_reduced + # Print masses. + if verbose + println(io, "\e[KFictitious mass | M = ", M) + end - # Print redcued masses. - if verbose - println(io, "\e[KReduced mass | Mr = ", M_reduced) - end + # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Thermal + M_asymp = v ./ w + p["Ma"][i, d, j, :] .= M_asymp - # Calculate and store polaron radii. - R = sqrt.(3 .* v ./ (v .^ 2 .- w .^ 2) .^ 2) - p["R"][i, j, :] .= R + # Print asymptotic masses. + if verbose + println(io, "\e[KFictitious mass (asymptotic) | Ma = ", M_asymp) + end - # Print polaron radius. - if verbose - println(io, "\e[KPolaron radius | R = ", R) - end + # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. + M_reduced = (v .^2 .- w .^2) / v .^ 2 + p["Mr"][i, d, j, :] .= M_reduced + + # Print redcued masses. + if verbose + println(io, "\e[KReduced mass | Mr = ", M_reduced) + end - if verbose - println("\e[K-----------------------------------------------------------------------") - println("\e[K DC Mobility Information: ") - println("\e[K-----------------------------------------------------------------------") - end + # Calculate and store polaron radii. + R = sqrt.(3 .* v ./ (v .^ 2 .- w .^ 2) .^ 2) + p["R"][i, d, j, :] .= R + + # Print polaron radius. + if verbose + println(io, "\e[KPolaron radius | R = ", R) + end - # Calculate and store the DC mobiliies. - μ = frohlich_mobility(v, w, α, ω, β; dims = dims) / mb - p["μ"][i, j] = μ + if verbose + println("\e[K-----------------------------------------------------------------------") + println("\e[K DC Mobility Information: ") + println("\e[K-----------------------------------------------------------------------") + end - # Print DC mobilities. - if verbose - println(io, "\e[KFinite temperature mobility | μ = ", μ) - end + # Calculate and store the DC mobiliies. + μ = !kspace ? frohlich_mobility(v, w, α, ω, β; dims = dims[d]) / mb : frohlich_mobility_k_space(v, w, α, ω, β; dims = dims[d]) / mb + p["μ"][i, d, j] = μ + + # Print DC mobilities. + if verbose + println(io, "\e[KFinite temperature mobility | μ = ", μ) + end - # FHIP low-temperature mobility, final result of Feynman1962. - # Eqn. (1.60) in Devreese2016 page 36; 6th Edition of Frohlich polaron notes (ArXiv). - μ_FHIP = FHIP_mobility_lowT(v, w, α, ω, β) / mb - p["μFHIP"][i, j] = μ_FHIP + # FHIP low-temperature mobility, final result of Feynman1962. + # Eqn. (1.60) in Devreese2016 page 36; 6th Edition of Frohlich polaron notes (ArXiv). + μ_FHIP = FHIP_mobility_lowT(v, w, α, ω, β) / mb + p["μFHIP"][i, j] = μ_FHIP + + # Print low-temperature FHIP mobility. + if verbose + println(io, "\e[KFHIP low-temp. mobility | μFHIP = ", μ_FHIP) + end - # Print low-temperature FHIP mobility. - if verbose - println(io, "\e[KFHIP low-temp. mobility | μFHIP = ", μ_FHIP) - end + # Kadanoff low-temperaure mobility, constructed around Boltzmann equation. + # Adds factor of 3 / (2 * β) c.f. FHIP, correcting phonon emission behaviour. + # Provide also the Kadanoff mobility that is consistent with the + # FHIP, and later statements (Devreese) of the Kadanoff mobility. + # It suggests that Kadanoff used the wrong identy for Nbar in Eqn. (23b) for + # the Γ₀ function, and should have used a version with the -1 to + # account for Bose / phonon statistics! + μ_Kadanoff_Devreese, μ_Kadanoff, rel_time = Kadanoff_mobility_lowT(v, w, α, ω, β) ./ mb + p["μD"][i, j] = μ_Kadanoff_Devreese + p["μK"][i, j] = μ_Kadanoff + p["τ"][i, j] = rel_time + + # Print low-temperature Kadanoff mobility (both original and Devreese-corrected versions). + if verbose + println(io, "\e[KDevreese low-temp. mobility | μD = ", μ_Kadanoff_Devreese) + println(io, "\e[KKadanoff low-temp. mobility | μK = ", μ_Kadanoff) + end - # Kadanoff low-temperaure mobility, constructed around Boltzmann equation. - # Adds factor of 3 / (2 * β) c.f. FHIP, correcting phonon emission behaviour. - # Provide also the Kadanoff mobility that is consistent with the - # FHIP, and later statements (Devreese) of the Kadanoff mobility. - # It suggests that Kadanoff used the wrong identy for Nbar in Eqn. (23b) for - # the Γ₀ function, and should have used a version with the -1 to - # account for Bose / phonon statistics! - μ_Kadanoff_Devreese, μ_Kadanoff, rel_time = Kadanoff_mobility_lowT(v, w, α, ω, β) ./ mb - p["μD"][i, j] = μ_Kadanoff_Devreese - p["μK"][i, j] = μ_Kadanoff - p["τ"][i, j] = rel_time - - # Print low-temperature Kadanoff mobility (both original and Devreese-corrected versions). - if verbose - println(io, "\e[KDevreese low-temp. mobility | μD = ", μ_Kadanoff_Devreese) - println(io, "\e[KKadanoff low-temp. mobility | μK = ", μ_Kadanoff) - end + #C alculates the DC mobility using Hellwarth et al. 1999 Eqn. (2). + # Directly performs contour integration in Feynman1962, for finite temperature DC mobility. + # Eqns. (2) and (1) are going back to the general (pre low-T limit) formulas in Feynman1962. + # To evaluate these, you need to do the explicit contour integration to get the polaron self-energy. + # See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. + μ_Hellwarth, μ_Hellwarth_b0 = Hellwarth_mobility(v, w, α, ω, β) ./ mb + p["μH"][i, j] = μ_Hellwarth + p["μH0"][i, j] = μ_Hellwarth_b0 + + # Print Hellwarth mobility (both with b and b=0) and Kadanoff relaxation time. + if verbose + println(io, "\e[KHellwarth mobility | μH = ", μ_Hellwarth) + println(io, "\e[KHellwarth mobility (b=0) | μH0 = ", μ_Hellwarth_b0) + println(io, "\e[KKadanoff relaxation time | τ = ", rel_time) + end - #C alculates the DC mobility using Hellwarth et al. 1999 Eqn. (2). - # Directly performs contour integration in Feynman1962, for finite temperature DC mobility. - # Eqns. (2) and (1) are going back to the general (pre low-T limit) formulas in Feynman1962. - # To evaluate these, you need to do the explicit contour integration to get the polaron self-energy. - # See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. - μ_Hellwarth, μ_Hellwarth_b0 = Hellwarth_mobility(v, w, α, ω, β) ./ mb - p["μH"][i, j] = μ_Hellwarth - p["μH0"][i, j] = μ_Hellwarth_b0 + else + # If zero temperature. + v, w, β = v_gs, w_gs, Inf + 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 + p["C"][i, d, j] = C_gs + p["κ"][i, d, j, :] .= κ_gs + p["M"][i, d, j, :] .= M_gs + p["μH"][i, j] = Inf + p["μH0"][i, j] = Inf + p["μD"][i, j] = Inf + p["μK"][i, j] = Inf + p["τ"][i, j] = 0 + p["μFHIP"][i, j] = Inf + p["μ"][i, d, j] = Inf + p["R"][i, d, j, :] .= R_gs + p["Mr"][i, d, j, :] .= M_reduced_gs + p["Ma"][i, d, j, :] .= M_asymp_gs + end # End of zero temperature if statement. - # Print Hellwarth mobility (both with b and b=0) and Kadanoff relaxation time. if verbose - println(io, "\e[KHellwarth mobility | μH = ", μ_Hellwarth) - println(io, "\e[KHellwarth mobility (b=0) | μH0 = ", μ_Hellwarth_b0) - println(io, "\e[KKadanoff relaxation time | τ = ", rel_time) Ωprocess = 1 # Counter for Ωrange. Tprocess += 1 # Increment Trange iterator. end @@ -468,50 +509,69 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses for k in eachindex(Ωrange) # E-field frequencies loop. Ω = Ωrange[k] - # Print E-field frequency. - if verbose - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[K Frequency Response Information: [$(Ωprocess[]) / $(num_Ω) ($(round(Ωprocess[] / (num_Ω) * 100, digits=1)) %)]") - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[KElectric field frequency | Ω = ", Ω) - end + if !iszero(Ω) - # Calculate and store polaron memory functions (akin to self energy). - χ = frohlich_memory_function(Ω, v, w, α, ω, β; dims = dims) - p["χ"][k, i, j] = χ + # Print E-field frequency. + if verbose + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[K Frequency Response Information: [$(Ωprocess[]) / $(num_Ω) ($(round(Ωprocess[] / (num_Ω) * 100, digits=1)) %)]") + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[KElectric field frequency | Ω = ", Ω) + end - # Print memory function. - if verbose - println(io, "\e[KMemory function | χ = ", χ) - end + # Calculate and store polaron memory functions (akin to self energy). + χ = !kspace ? frohlich_memory_function(Ω, v, w, α, ω, β; dims = dims[d]) : frohlich_memory_function_k_space(Ω, v, w, α, ω, β; dims = dims[d]) + p["χ"][k, i, d, j] = χ - # Calculate and store polaron complex impedances. + # Print memory function. + if verbose + println(io, "\e[KMemory function | χ = ", χ) + end - z = -(im * Ω + im * χ) .* mb - p["z"][k, i, j] = z + # Calculate and store polaron complex impedances. - # Print complex impedances. - if verbose - println(io, "\e[KComplex impedance | z = ", z) - end + z = -(im * Ω + im * χ) .* mb + p["z"][k, i, d, j] = z + + # Print complex impedances. + if verbose + println(io, "\e[KComplex impedance | z = ", z) + end + + # Calculate and store polaron complex conductivities. + σ = 1 / z + p["σ"][k, i, d, j] = σ + + # Print complex conductivities and show total algorithm progress. + if verbose + println(io, "\e[KComplex conductivity | σ = ", σ) + end + + else - # Calculate and store polaron complex conductivities. - σ = 1 / z - p["σ"][k, i, j] = σ + # If zero frequency. + p["χ"][k, i, d, j] = Inf + 0 * im + p["z"][k, i, d, j] = 0 + im * Inf + p["σ"][k, i, d, j] = 0 + 0 * im + + end # End of zero frequency if statement. - # Print complex conductivities and show total algorithm progress. if verbose - println(io, "\e[KComplex conductivity | σ = ", σ) println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[K[Total Progress: $(process[]) / $(num_α * num_T * num_Ω) ($(round(process[] / (num_α * num_T * num_Ω) * 100, digits=1)) %)]") - print(io, "\e[9F") + println(io, "\e[K[Total Progress: $(process[]) / $(num_α * num_T * num_Ω * num_d) ($(round(process[] / (num_α * num_T * num_Ω * num_d) * 100, digits=1)) %)]") Ωprocess += 1 # Increment Ωrange iterator. process += 1 # Increment total iterator. + print(io, "\e[2F") end + + if verbose && !iszero(Ω) print(io, "\e[7F") end end - if verbose print(io, "\e[26F") end # Move up 26 lines and erase. + if verbose && !iszero(T) print(io, "\e[26F") end # Move up 26 lines and erase. end - if verbose print(io, "\e[27F") end # Move up 26 lines and erase. + if verbose print(io, "\e[15F") end # Move up 15 lines and erase. + end # End dimensions loop. + + if verbose print(io, "\e[12F") end # Move up 12 lines and erase. end if verbose print("\e[?25h") end # Show cursor again. @@ -567,28 +627,28 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses ] # Return Polaron type containing all generated data over any coupling strengths, temperatures and frequencies. - return Polaron(polaron_data...) + return FrohlichPolaron(polaron_data...) end """ Single alpha parameter. polaron() expects alpha parameters to be in a Vector. """ -polaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = polaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +frohlichpolaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No frequency input. """ -polaron(αrange, Trange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = polaron(αrange, Trange, 1; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +frohlichpolaron(αrange, Trange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(αrange, Trange, 0; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No temperature input => 300 K. """ -polaron(αrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = polaron(αrange, 300, 1; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +frohlichpolaron(αrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(αrange, 0, 0; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No input => α = 1 """ -polaron(; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = polaron(1, 300, 1; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +frohlichpolaron(; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(1, 0, 0; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ polaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, verbose=false) @@ -596,7 +656,7 @@ Material specific constructors that use material specific parameters to paramete Material data is inputted through the `Material` type. Returns all data in either SI units or other common, suitable units otherwise. """ -function polaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) +function frohlichpolaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) # Show material data. if verbose @@ -609,14 +669,14 @@ function polaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87 mb = material.mb # Generate polaron data from the arbitrary model constructor. - p = polaron(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, β0=ħ/kB*1e12*2π, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) + p = frohlichpolaron(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, β0=ħ/kB*1e12*2π, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) # Return material-specific, unitful Polaron type. return p end # Broadcast Polaron data. -function Base.show(io::IO, ::MIME"text/plain", x::Polaron) +function Base.show(io::IO, ::MIME"text/plain", x::FrohlichPolaron) io_lim = IOContext(io, :limit => true, :compact => true) @@ -694,13 +754,13 @@ function Base.show(io::IO, ::MIME"text/plain", x::Polaron) end """ - save_polaron(p::Polaron, prefix) + save_polaron(p::FrohlichPolaron, prefix) Saves data from 'polaron' into file "prefix". This is a .jdl file for storing the polaron data whilst preserving types. Allows for saving multidimensional arrays that sometimes arise in the polaron data. Each parameter in the NewPolaron type is saved as a dictionary entry. E.g. NewPolaron.α is saved under JLD.load("prefix.jld")["alpha"]. """ -function save_polaron(polaron::Polaron, prefix) +function save_frohlich_polaron(polaron::FrohlichPolaron, prefix) println("Saving polaron data to $prefix.jld ...") @@ -758,19 +818,17 @@ function save_polaron(polaron::Polaron, prefix) end """ - load_polaron(p::NewPolaron, prefix) + load_polaron(p::FrohlichPolaron, prefix) Loads data from file "polaron_file_path" into a NewPolaron type. """ -function load_polaron(polaron_file_path) +function load_frohlich_polaron(polaron_file_path) println("Loading polaron data from $polaron_file_path ...") data = JLD.load("$polaron_file_path") - println(typeof(data["alpha"])) - - polaron = Polaron( + polaron = FrohlichPolaron( data["alpha"], data["alpha eff"], data["temperature"], diff --git a/src/FrohlichTheory.jl b/src/FrohlichTheory.jl new file mode 100644 index 0000000..e71d94c --- /dev/null +++ b/src/FrohlichTheory.jl @@ -0,0 +1,625 @@ +# FrohlichTheory.jl +# Everything specifically associated with calculating properties of the Frohlich polaron model. + +""" + frohlichalpha(ε_Inf, ε_S, freq, m_eff) + +Calculates the Frohlich alpha parameter, for a given dielectric constant, frequency (f) of phonon in Hertz, and effective mass (in units of the bare electron mass). + +See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. +""" +function frohlichalpha(ϵ_optic, ϵ_static, freq, m_eff) + ω = freq * 2π * 1e12 # frequency to angular velocity + # Note: we need to add a 4*pi factor to the permitivity of freespace. + # This gives numeric agreement with literature values. This is required as + # the contemporary 1950s and 1960s literature implicitly used atomic units, + # where the electric constant ^-1 has this factor baked in, k_e=1/(4πϵ_0). + α = 1 / 2 / (4 * π * ϵ_0) * # Units: m/F + (1 / ϵ_optic - 1 / ϵ_static) * # Units: none + (eV^2 / (ħ * ω)) * # Units: F + sqrt(2 * m_eff * me * ω / ħ) # Units: 1/m + return α +end + +""" + frohlichalpha(ϵ_optic, ϵ_ionic, ϵ_total, phonon_mode_freq, m_eff) + +Calculates the partial dielectric electron-phonon coupling parameter for a given longitudinal optical phonon mode. + +This decomposes the original Frohlich alpha coupling parameter (defined for a single phonon branch) into contributions from multiple phonon +branches. + +# Arguments +- `ϵ_optic::Float64`: is the optical dielectric constant of the material. +- `ϵ_ionic::Float64`: is the ionic dielectric contribution from the phonon mode. +- `ϵ_total::Float64`: is the total ionic dielectric contribution from all phonon modes of the material. +- `phonon_mode_freq::Float64`: is the frequency of the phonon mode (THz). +- `m_eff::Float64` is the band mass of the electron (in units of electron mass m_e). +""" +function frohlichalpha(ϵ_optic, ϵ_ionic, ϵ_total, phonon_mode_freq, m_eff) + + # The Rydberg energy unit + Ry = eV^4 * me / (2 * ħ^2) + + # Angular phonon frequency for the phonon mode (rad Hz). + ω = phonon_mode_freq * 2π * 1e12 + + # The static dielectric constant. Calculated here instead of inputted so that ionic modes are properly normalised. + ϵ_static = ϵ_total + ϵ_optic + + # The contribution to the electron-phonon parameter from the currrent phonon mode. 1 / (4π ϵ_0) is the dielectric normalisation. + α_j = (m_eff * Ry / (ħ * ω))^(1 / 2) * ϵ_ionic / (4π * ϵ_0) / (ϵ_optic * ϵ_static) + + return α_j +end + +# Partial dielectric electron-phonon coupling parameter, decomposed into ionic dielectric contributions from each phonon mode of the material. + +""" + ϵ_ionic_mode(phonon_mode_freq, ir_activity, volume) + +Calculate the ionic contribution to the dielectric function for a given phonon mode. + +# Arguments +- `phonon_mode_freq::Float64`: is the frequency of the mode in THz. +- `ir_activity::Float64`: is the infra-red activity of the mode in e²amu⁻¹. +- `volume::Float64`: is the volume of the unit cell of the material in m³. +""" +function ϵ_ionic_mode(phonon_mode_freq, ir_activity, volume) # single ionic mode + + # Angular phonon frequency for the phonon mode (rad Hz) + ω_j = phonon_mode_freq * 2π * 1e12 + + # Dielectric contribution from a single ionic phonon mode + ϵ_mode = eV^2 * ir_activity / (3 * volume * ω_j^2 * amu) + + # Normalise ionic dielectric contribution with 1 / (4π ϵ_0) (NB: the 4π has been pre-cancelled) + return ϵ_mode / ϵ_0 +end + +""" + ϵ_total(freqs_and_ir_activity, volume) + +Calculate the total ionic contribution to the dielectric function from all phonon modes. + +# Arguments +- `freqs_and_ir_activity::Matrix{Float64}`: is a matrix containeing the phonon mode frequencies (in THz) in the first column and the infra-red activities (in e²amu⁻¹) in the second column. +- `volume::Float64`: is the volume of the unit cell of the material in m^3. +""" +function ϵ_total(freqs_and_ir_activity, volume) # total ionic contribution to dielectric + + # Extract phonon frequencies (THz) + phonon_freqs = freqs_and_ir_activity[:, 1] + + # Extra infra-red activities (e^2 amu^-1) + ir_activity = freqs_and_ir_activity[:, 2] + + # Sum over all ionic contribution from each phonon mode + total_ionic = 0 + + for t in eachindex(phonon_freqs) + total_ionic += ϵ_ionic_mode(phonon_freqs[t], ir_activity[t], volume) + end + + return total_ionic +end + +""" + frohlich_coupling(k, α, ω) + +Calculate the coupling strength for the Frohlich continuum polaron model. + +# Arguments +- `k`: a scalar value representing the k-coordinate in k-space +- `α`: a scalar value representing the coupling constant +- `ω`: a scalar value representing the phonon frequency + +# Returns +The coupling strength for the Frohlich continuum polaron model. + +# Example +```julia +result = frohlich_coupling(2.0, 0.5, 1.0) +println(result) +``` +Expected Output: +`6.0` +""" +function frohlich_coupling(k, α, ω; mb = 1, dims = 3) + r_p = 1 / sqrt(2) + ω^2 * α * r_p * gamma((dims - 1) / 2) * (2√π / k)^(dims - 1) +end + + +""" + frohlich_interaction_energy(v, w, α, ωβ...) + +Integral of Eqn. (31) in Feynman 1955. Part of the overall ground-state energy expression. + +See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. +""" +function frohlich_interaction_energy(v, w, α, ωβ...; dims = 3) + coupling = frohlich_coupling(1, α, ωβ[1]; dims = dims) + integrand(τ) = phonon_propagator(τ, ωβ...) / sqrt(polaron_propagator(τ, v, w, ωβ...)) + upper_limit = length(ωβ) == 1 ? Inf : prod(ωβ) / 2 + integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) + 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, α, ω; dims = dims)) + +frohlich_interaction_energy(v, w, α::Vector, ω::Vector, β; dims = 3) = sum(frohlich_interaction_energy.(v, w, α, ω, β; dims = dims)) + +""" + frohlich_interaction_energy_k_space(v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) + +Calculate the Frohlich polaron interaction energy in k-space at finite temperaure. + +## Arguments +- `τ`: a scalar value representing the imaginary time. +- `v`: a scalar value representing a variational paramater. +- `w`: a scalar value representing a variational paramater. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `a`: an optional scalar value representing the lattice constant. Default value is 1. +- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. + +## Returns +A scalar value representing the Frohlich polaron interaction energy in k-space at finite temperature. +""" +function frohlich_interaction_energy_k_space(v, w, α, ωβ...; limits = [0, Inf], dims = 3) + coupling(k) = frohlich_coupling(k, α, ωβ[1]; dims = dims) + propagator(τ) = polaron_propagator(τ, v, w, ωβ...) + integrand(τ) = phonon_propagator(τ, ωβ...) * spherical_k_integral(coupling, propagator(τ); dims = dims, limits = limits) + upper_limit = length(ωβ) == 1 ? Inf : prod(ωβ) / 2 + integral, _ = quadgk(τ -> frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; limits = limits, dims = dims), 0, upper_limit) + return integral +end + +""" + B(v, w, α, β) + +Generalisation of the B function from Eqn. (62c) in Hellwarth et al. 1999. This is the expected value of the exact action taken w.r.t trial action, given for the 'jth' phonon mode. + +Required for calculating the polaron free energy. + +# Arguments +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for the 'jth' phonon mode. +- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. + +See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. +""" +function B(v::Vector, w::Vector, α, ω, β) + B_integrand(τ) = cosh(τ - ω * β / 2) / sqrt(abs(D(τ, v, w, ω, β))) + return α / (√π * sinh(ω * β / 2)) * quadgk(τ -> B_integrand(τ * ω * β / 2), 0.0, 1.0)[1] * ω^2 * β / 2 +end + +B(v::Vector, w::Vector, α::Vector, ω::Vector, β) = sum(B(v, w, α[i], ω[i], β) for i in eachindex(ω)) + +""" + B(v, w, α; rtol = 1e-3) + +Calculates `B(v, w, α, β)` but at zero-temperature, `β = Inf`. + +# Arguments +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for the 'jth' phonon mode. +""" +function B(v::Vector, w::Vector, α, ω) + B_integrand(τ) = exp(-abs(τ)) / sqrt(abs(D(abs(τ), v, w))) + return α / √π * quadgk(τ -> B_integrand(τ), 0, Inf)[1] * ω +end + +B(v::Vector, w::Vector, α::Vector, ω::Vector) = sum(B(v, w, α[i], ω[i]) for i in eachindex(ω)) + +""" + F(v, w, α, ω, β) + +Calculates the Helmholtz free energy of the polaron for a material with multiple phonon branches. + +This generalises the Osaka 1959 (below Eqn. (22)) and Hellwarth. et al 1999 (Eqn. (62a)) free energy expressions. + +# Arguments +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for the 'jth' phonon mode. +- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz). +- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. + +See Osaka, Y. (1959): https://doi.org/10.1143/ptp.22.437 and Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. +""" +function frohlich_energy(v, w, α, ωβ...; dims = 3) + Ar, Cr = trial_energy(v, w, ωβ...; dims = dims) + Br = frohlich_interaction_energy(v, w, α, ωβ...; dims = dims) + return -(Ar + Br + Cr), Ar, Br, Cr +end + +""" + frohlich_energy_k_space(v, w, α, ωβ...; rₚ = 1, limits = [0, Inf]) + +Calculate the total energy, kinetic energy, and interaction energy of the Frohlich lattice polaron. + +## Arguments +- `v`: a scalar value representing a variational paramater. +- `w`: a scalar value representing a variational paramater. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: (optional) a scalar value representing the inverse temperature. +- `rₚ`: The characteristic polaron radius (default is 1). +- `limits`: The limits of integration for the interaction energy calculation (default is [0, Inf]). + +## Returns +- `total_energy`: The calculated total polaron energy. +- `kinetic_energy`: The calculated polaron kinetic energy. +- `interaction_energy`: The calculated polaron interaction energy. +""" +function frohlich_energy_k_space(v, w, α, ωβ...; limits = [0, Inf], dims = 3) + A, C = electron_energy(v, w, ωβ...) + B = frohlich_interaction_energy_k_space(v, w, α, ωβ...; limits = limits, dims = dims) + return -(A + B + C), A, B, C +end + +""" +---------------------------------------------------------------------- +The Complex Impedence and Conductivity of the Polaron. +---------------------------------------------------------------------- +""" + +# - Polaron optical absorption, complex impedence and complex conductivity +# References: +# [1] R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Mobility of slow electrons in a polar crystal, PhysicalReview127, 1004 (1962) +# [2] Devreese, J. De Sitter, and M. Goovaerts, “Optical absorption of polarons in thefeynman-hellwarth-iddings-platzman approximation,”Phys. Rev. B, vol. 5, pp. 2367–2381, Mar 1972. +# [3] F. Peeters and J. Devreese, “Theory of polaron mobility,” inSolid State Physics,pp. 81–133, Elsevier, 1984. + +function frohlich_structure_factor(t, v, w, α, ωβ...; dims = 3) + coupling = frohlich_coupling(1, α, ωβ[1]; dims = dims) + propagator = polaron_propagator(im * t, v, w, ωβ...) / 2 + integral = ball_surface(dims) / (2π)^dims * √π / 4 / propagator^(3/2) + return 2 / dims * coupling * integral * phonon_propagator(im * t, ωβ...) +end + +""" + frohlich_structure_factor_k_space(t, v, w, α, ωβ...; limits = [0, Inf]) + +Calculate the structure factor in k-space for the Frohlich continuum polaron model at finite temperature. + +## Arguments +- `t`: a scalar value representing the real time. +- `v`: a scalar value representing a variational parameter. +- `w`: a scalar value representing a variational parameter. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. + +## Returns +A scalar value representing the calculated structure factor in k-space for the Frohlich continuum polaron model at finite temperature. +""" +function frohlich_structure_factor_k_space(t, v, w, α, ωβ...; limits = [0, Inf], dims = 3) + coupling(k) = frohlich_coupling(k, α, ωβ[1]; dims = dims) * k^2 + propagator = polaron_propagator(im * t, v, w, ωβ...) + integral = spherical_k_integral(coupling, propagator; limits = limits, dims = dims) + return 2 / dims * phonon_propagator(im * t, ω, β) * integral +end + +function frohlich_memory_function(Ω, v, w, α, ωβ...; dims = 3) + structure_factor(t) = frohlich_structure_factor(t, v, w, α, ωβ...; dims = dims) + return polaron_memory_function(Ω / ωβ[1], structure_factor; limits = [0, 1e5]) +end + +""" + frohlich_memory_function_k_space(Ω, v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) + +Calculate the memory function for the Frohlich model in k-space at finite temperature and frequency. + +## Arguments +- `Ω`: a scalar value representing the frequency at which the memory function is evaluated. +- `v`: a scalar value representing the optimal variational parameter. +- `w`: a scalar value representing the optimal variational parameter. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. + +## Returns +A scalar value representing the memory function of the Frohlich model in k-space at finite temperature evaluated at the frequency `Ω`. +""" +function frohlich_memory_function_k_space(Ω, v, w, α, ωβ...; dims = 3, limits = [0, Inf]) + structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ωβ...; dims = dims, limits = limits) + return polaron_memory_function(Ω, structure_factor) +end + +frohlich_memory_function(Ω, v::Vector, w::Vector, α, ωβ...) = sum(frohlich_memory_function.(Ω, v, w, α, ωβ...)) + +frohlich_memory_function(Ω, v, w, α::Vector, ωβ...) = sum(frohlich_memory_function.(Ω, v, w, α, ωβ...)) + + +""" + frohlich_complex_impedence(Ω, β, α, v, w; rtol = 1e-3, T = nothing, verbose = false) + +Calculate the complex impedence Z(Ω) of the polaron at finite temperatures for a given frequency Ω (Eqn. (41) in FHIP 1962). + +# Arguments +- `Ω::Float64`: is the frequency (2π THz) of applied electric field. +- `β::Float64`: is the reduced thermodynamic betas. +- `α::Float64`: is the Frohlich alpha coupling parameter. +- `v::Float64`: is the 'v' variational parameter. +- `w::Float64`: is the 'w' variational parameter. +- `rtol`: relative tolerance for the accuracy of any quadrature integrations. +- `T`: is a token used by `make_polaron()` to keep track of the temperature for printing during a calculation. Do not alter. +- `verbose`: is used by `make_polaron()` to specify whether or not to print. Ignore. + +See FHIP 1962: https://doi.org/10.1103/PhysRev.127.1004. + +""" +function frohlich_complex_impedence(Ω, v, w, α, ωβ...; dims = 3) + -im * (Ω + frohlich_memory_function(Ω, v, w, α, ωβ...; dims = dims)) +end + +""" + frohlich_complex_conductivity(Ω, β, α, v, w; rtol = 1e-3) + +Calculate the complex conductivity σ(Ω) of the polaron at finite temperatures for a given frequency Ω (equal to 1 / Z(Ω) with Z the complex impedence). + +# Arguments +- `Ω::Float64`: is the frequency (2π THz) of applied electric field. +- `β::Float64`: is the reduced thermodynamic betas. +- `α::Float64`: is the Frohlich alpha coupling parameter. +- `v::Float64`: is the 'v' variational parameter. +- `w::Float64`: is the 'w' variational parameter. +- `rtol`: relative tolerance for the accuracy of any quadrature integrations. + +See also [`polaron_complex_impedence`](@ref) +""" +function frohlich_complex_conductivity(Ω, v, w, α, ωβ...; dims = 3) + return 1 / frohlich_complex_impedence(Ω, v, w, α, ωβ...; dims = dims) +end + +""" + inverse_polaron_mobility(v, w, α, ω, β) + +Calculate the inverse of the dc mobility μ of the polaron at finite temperatues (Eqn. (11.5) in F. Peeters and J. Devreese 1984) for a given frequency Ω. + +# Arguments +- `v::Float64`: is the 'v' variational parameter. +- `w::Float64`: is the 'w' variational parameter. +- `α::Float64`: is the Frohlich alpha coupling parameter. +- `ω::Float64`: is the angular phonon frequency. +- `β::Float64`: is the reduced thermodynamic beta. +- `verbose`: is used by `make_polaron()` to specify whether or not to print. Ignore. + +See F. Peeters and J. Devreese 1984: https://doi.org/10.1016/S0081-1947(08)60312-4. + +See also [`polaron_mobility`](@ref), [`polaron_complex_conductivity`](@ref) +""" +function inverse_frohlich_mobility(v, w, α, ω, β; dims = 3) + if β == Inf + return zero(β) + end + structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω, β; dims = dims) + return abs(imag(polaron_memory_function(structure_factor; limits = [0, 1e5]))) +end + +""" + inverse_polaron_mobility(v, w, α::Vector, ω::Vector, β::Vector) + +inverse of the polaron mobility, but for multiple phonon modes. +""" +inverse_frohlich_mobility(v, w, α::Vector, ω::Vector, β) = sum(inverse_frohlich_mobility.(v, w, α, ω, β)) + +""" + polaron_mobility(v, w, α, ω, β) + +The polaron mobility. + +See also [`inverse_polaron_mobility`](@ref) +""" +frohlich_mobility(v, w, α, ω, β; dims = 3) = 1 ./ inverse_frohlich_mobility(v, w, α, ω, β; dims = dims) + +""" + frohlich_mobility_k_space(v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) + +Calculate the DC mobility in k-space for a Frohlich polaron system at finite temperature. + +## Arguments +- `v`: a scalar value representing the optimal variational parameter. +- `w`: a scalar value representing the optimal variational parameter. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. + + +## Returns +The DC mobility in k-space for the Frohlich polaron system at finite temperature. +""" +function frohlich_mobility_k_space(v, w, α, ω, β; dims = 3, limits = [0, Inf]) + if β == Inf + return Inf + end + structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ω, β; dims = dims, limits = limits) + 1 / imag(polaron_memory_function(structure_factor)) +end + +""" + inverse_FHIP_mobility_lowT(v, w, α, ω, β) + +FHIP low-temperature mobility, final result of Feynman1962. +[1.60] in Devreese2016 page 36; 6th Edition of Frohlich polaron notes (ArXiv). + +See also [`FHIP_mobility_lowT`](@ref) +""" +function inverse_FHIP_mobility_lowT(v, w, α, ω, β) + if β == Inf + return 0 + else + μ = (w / v)^3 * 3 / (4 * ω^2 * α * β) * exp(ω * β) * exp((v^2 - w^2) / (w^2 * v)) + return 1 / μ + end +end + +""" + inverse_FHIP_mobility_lowT(v, w, α::Vector, ω::Vector, β::Vector) + +Inverse FHIP mobility for multiple phonon modes. + +See also [`FHIP_mobility_lowT`](@ref) +""" +inverse_FHIP_mobility_lowT(v, w, α::Vector, ω::Vector, β) = sum(inverse_FHIP_mobility_lowT.(v, w, α, ω, β)) + +""" + FHIP_mobility_lowT(v, w, α, ω, β) + +FHIP mobility in the low-temperature approximation. + +See also [`inverse_FHIP_mobility_lowT`](@ref) +""" +FHIP_mobility_lowT(v, w, α, ω, β) = 1 ./ inverse_FHIP_mobility_lowT(v, w, α, ω, β) + +""" + inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) + +Kadanoff low-temperaure mobility, constructed around Boltzmann equation. +Adds factor of 3 / (2 * β) c.f. FHIP, correcting phonon emission behaviour. + +See also [`Kadanoff_mobility_lowT`](@ref) +""" +function inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) + + if β == Inf + return 0, 0, 0 + else + + # [1.61] in Devreese2016 - Kadanoff's Boltzmann eqn derived mobility. + μ_Devreese2016 = (w / v)^3 / (2 * ω * α) * exp(ω * β) * exp((v^2 - w^2) / (w^2 * v)) + + # From 1963 Kadanoff (particularly on right-hand-side of page 1367), Eqn. (9), + # we define equilibrium number of phonons (just from temperature T and phonon ω): + # N̄ = (exp(β) - 1)^-1. + # But! We find that: + N̄ = exp(-β * ω) + + # is only way to get Kadanoff1963 to be self-consistent with + # FHIP, and later statements (Devreese) of the Kadanoff mobility. + # It suggests that Kadanoff used the wrong identy for Nbar in Eqn. (23b) for + # the Γ₀ function, and should have used a version with the -1 to + # account for Bose / phonon statistics! + + # Between Eqns. (23) and (24) in Kadanoff1963, for small momenta skip intergration + # and sing the fictitious mass M: + M = (v^2 - w^2) / w^2 + + # we get for Γ₀: + Γ₀ = 2 * α * N̄ * (M + 1)^(1 / 2) * exp(-M / v) * ω + + # NB: Kadanoff1963 uses ħ = ω = mb = 1 units. + # Factor of omega to get it as a rate relative to phonon frequency + # Factor of omega*hbar to get it as a rate per energy window. + # Hence, for the mobility from Eqn. (25) in Kadanoff1963 (SI effective mass q / mb): + μ_Kadanoff1963 = 1 / ((M + 1) * Γ₀) + + # # Energy loss is simply energy * rate: + # energy_loss = ω * Γ₀ / 2π + + # # Lifetime is: + # τ = 1 / Γ₀ + + return 1 / μ_Devreese2016, 1 / μ_Kadanoff1963, Γ₀ + end +end + +""" + inverse_Kadanoff_mobility_lowT(v, w, α::Vector, ω::Vector, β::Vector) + +Inverse Kadanoff mobility for multiple phonon modes. + +See also [`Kadanoff_mobility_lowT`](@ref) +""" +inverse_Kadanoff_mobility_lowT(v, w, α::Vector, ω::Vector, β) = map(+, map(x -> 1 ./ x, inverse_Kadanoff_mobility_lowT.(v, w, α, ω, β))...) + +""" + Kadanoff_mobility_lowT(v, w, α, ω, β) + +Kadanoff mobility in the low-temperature approximation. + +See also [`inverse_Kadanoff_mobility_lowT`](@ref) +""" +Kadanoff_mobility_lowT(v, w, α, ω, β) = 1 ./ inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) + +""" + inverse_Hellwarth_mobility(v, w, α, ω, β) + +Calculates the DC mobility using Hellwarth et al. 1999 Eqn. (2). +Directly performs contour integration in Feynman1962, for finite temperature DC mobility. +Eqns. (2) and (1) are going back to the general (pre low-T limit) formulas in Feynman1962. +To evaluate these, you need to do the explicit contour integration to get the polaron self-energy. + +See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. + +See also [`Hellwarth_mobility`](@ref) +""" +function inverse_Hellwarth_mobility(v, w, α, ω, β) + + if β == Inf + return 0, 0 + else + + R = (v^2 - w^2) / (w^2 * v) # inline, page 300 just after Eqn (2) + b = R * ω * β / sinh(ω * β * v / 2) # Feynman1962 version; page 1010, Eqn (47b) + a = sqrt((ω * β / 2)^2 + R * ω * β * coth(ω * β * v / 2)) + k(u, a, b, v) = (u^2 + a^2 - b * cos(v * u) + eps(Float64))^(-3 / 2) * cos(u) # integrand in (2) + K = quadgk(u -> k(u, a, b, v), 0, 1e5)[1] # numerical quadrature integration of (2) + + # Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997 + RHS = α / (3 * sqrt(π)) * (ω * β)^(5 / 2) / sinh(ω * β / 2) * (v^3 / w^3) * K + μ = RHS + + # Hellwarth1999/Biaggio1997, b=0 version... 'Setting b=0 makes less than 0.1% error' + # So let's test this: + K_0 = quadgk(u -> k(u, a, 0, v), 0, 1e5)[1] # Inserted b=0 into k(u, a, b, v). + RHS_0 = α / (3 * sqrt(π)) * (ω * β)^(5 / 2) / sinh(ω * β / 2) * (v^3 / w^3) * K_0 + μ_0 = RHS_0 + + return μ, μ_0 + end +end + +""" + inverse_Hellwarth_mobility(v, w, α::Vector, ω::Vector, β::Vector) + +Inverse Hellwarth mobility for multiple phonon modes. + +See also [`Hellwarth_mobility`](@ref) +""" +inverse_Hellwarth_mobility(v, w, α::Vector, ω::Vector, β) = map(+, map(x -> 1 ./ x, inverse_Hellwarth_mobility.(v, w, α, ω, β))...) + +""" + Hellwarth_mobility(v, w, α, ω, β) + +The Hellwarth polaron mobility. + +See also [`inverse_Hellwarth_mobility`](@ref) +""" +Hellwarth_mobility(v, w, α, ω, β) = 1 ./ inverse_Hellwarth_mobility(v, w, α, ω, β) + +function polaron_effective_mass(v, w, α, ω, β, Ω) + + # FHIP1962, page 1009, eqn (35a). + integrand(t) = (1 - exp(im * Ω * t / ω)) * imag(S(t, v, w, ω, β)) + + integral = real(quadgk(t -> integrand(t), 0, Inf)[1]) + + mass = 1 - 2 * α * ω^2 * integral / (3 * √π * Ω^2) + + return mass +end + +polaron_effective_mass(v, w, α::Vector, ω::Vector, β) = sum(polaron_effective_mass.(v, w, α, ω, β)) \ No newline at end of file diff --git a/src/HolsteinPolaron.jl b/src/HolsteinPolaron.jl index ace0e02..e2ee4a5 100644 --- a/src/HolsteinPolaron.jl +++ b/src/HolsteinPolaron.jl @@ -1,11 +1,12 @@ -# Polaron.jl +# HolsteinPolaron.jl +# A do-it-all function for calculating properties of the Holstein polaron. Data is stored as a struct. """ - Holstein(x...) + HolsteinPolaron(x...) Type for storing the Holstein polaron information, `x...`. """ -mutable struct Holstein +mutable struct HolsteinPolaron α # Holstein coupling. αeff # Effective Holstein coupling summed for multiple modes. T # Temperature. @@ -17,8 +18,9 @@ mutable struct Holstein v0 # Athermal variational parameter v. w0 # Athermal variational parameter w. F0 # Polaron athermal total energy (enthalpy). - K0 # Polaron athermal kinetic energy. - P0 # Polaron athermal potential energy. + A0 # Bare electron enthalpy. + B0 # ⟨S⟩ₜ interaction enthalpy. + C0 # ⟨Sₜ⟩ₜ enthalpy of trial system. κ0 # Fictitious spring constant. M0 # Athermal fictitious mass. M0a # Athermal asymptotic approximate fictitious mass (v0/w0). @@ -27,8 +29,9 @@ mutable struct Holstein v # Thermal variational parameter v. w # Thermal variational parameter w. F # Polaron free energy. - K # Polaron thermal kinetic energy. - P # Polaron thermal potential energy. + A # Bare electron free energy. + B # ⟨S⟩ₜ interaction energy. + C # ⟨Sₜ⟩ₜ free energy of trial system. κ # Fictitious spring constant. M # Thermal fictitious mass. Ma # Thermal asymptotic approximate fictitious mass (v/w). @@ -39,12 +42,12 @@ mutable struct Holstein z # Complex impedence. σ # Complex conductivity. - function Holstein(x...) + function HolsteinPolaron(x...) new(reduce_array.(x)...) end end -function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=4.0, w_guesses=3.9, dims=3, verbose=false) +function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=4.0, w_guesses=3.9, dims=3, kspace=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. @@ -54,6 +57,7 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse num_T = length(Trange) num_Ω = length(Ωrange) num_ω = length(ω) + num_d = length(dims) ω = reduce_array(ω) @@ -71,30 +75,32 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse "β" => Matrix{Float64}(undef, num_T, num_ω), # Betas. "Ω" => Ωrange, # Photon frequencies. "d" => dims, # Number of spatial dimensions. - "v0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal v params. - "w0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal w params. - "F0" => Vector{Float64}(undef, num_α), # Athermal total energies. - "K0" => Vector{Float64}(undef, num_α), # Athermal kinetic energies. - "P0" => Vector{Float64}(undef, num_α), # Athermal potential energies. - "κ0" => Matrix{Float64}(undef, num_α, num_vw), # Fictitious spring constant. - "M0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal fictitious mass. - "M0a" => Matrix{Float64}(undef, num_α, num_vw), # Athermal asymptotic approximate fictitious mass (v0/w0). - "M0r" => Matrix{Float64}(undef, num_α, num_vw), # Athermal reduced mass of particle + fictitious particle system. - "R0" => Matrix{Float64}(undef, num_α, num_vw), # Athermal polaron radius (s. d. of a Gaussian wavefunction). - "v" => Array{Float64,3}(undef, num_T, num_α, num_vw), # v params. - "w" => Array{Float64,3}(undef, num_T, num_α, num_vw), # w params. - "F" => Matrix{Float64}(undef, num_T, num_α), # Total energies. - "K" => Matrix{Float64}(undef, num_T, num_α), # Kinetic energies. - "P" => Matrix{Float64}(undef, num_T, num_α), # Potential energies. - "κ" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Spring constants. - "M" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Fictitious masses. - "Ma" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Thermal asymptotic approximate fictitious mass (v/w). - "Mr" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Thermal reduced mass of particle + fictitious particle system. - "R" => Array{Float64,3}(undef, num_T, num_α, num_vw), # Polaron radii. - "μ" => Matrix{Float64}(undef, num_T, num_α), # Mobilities. - "χ" => Array{ComplexF64,3}(undef, num_Ω, num_T, num_α), # 'Memory function' or self energy. - "z" => Array{ComplexF64,3}(undef, num_Ω, num_T, num_α), # Complex impedences. - "σ" => Array{ComplexF64,3}(undef, num_Ω, num_T, num_α), # Complex conductivities. + "v0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal v params. + "w0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal w params. + "F0" => Matrix{Float64}(undef, num_d, num_α), # Athermal total energies. + "A0" => Matrix{Float64}(undef, num_d, num_α), # Athermal A parameter. + "B0" => Matrix{Float64}(undef, num_d, num_α), # Athermal B parameter. + "C0" => Matrix{Float64}(undef, num_d, num_α), # Athermal C parameter. + "κ0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Fictitious spring constant. + "M0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal fictitious mass. + "M0a" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal asymptotic approximate fictitious mass (v0/w0). + "M0r" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal reduced mass of particle + fictitious particle system. + "R0" => Array{Float64,3}(undef, num_d, num_α, num_vw), # Athermal polaron radius (s. d. of a Gaussian wavefunction). + "v" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # v params. + "w" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # w params. + "F" => Array{Float64,3}(undef, num_T, num_d, num_α), # Total energies. + "A" => Array{Float64,3}(undef, num_T, num_d, num_α), # A parameter. + "B" => Array{Float64,3}(undef, num_T, num_d, num_α), # B parameter. + "C" => Array{Float64,3}(undef, num_T, num_d, num_α), # C parameter. + "κ" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Spring constants. + "M" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Fictitious masses. + "Ma" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Thermal asymptotic approximate fictitious mass (v/w). + "Mr" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Thermal reduced mass of particle + fictitious particle system. + "R" => Array{Float64,4}(undef, num_T, num_d, num_α, num_vw), # Polaron radii. + "μ" => Array{Float64,3}(undef, num_T, num_d, num_α), # Mobilities. + "χ" => Array{ComplexF64,4}(undef, num_Ω, num_T, num_d, num_α), # 'Memory function' or self energy. + "z" => Array{ComplexF64,4}(undef, num_Ω, num_T, num_d, num_α), # Complex impedences. + "σ" => Array{ComplexF64,4}(undef, num_Ω, num_T, num_d, num_α), # Complex conductivities. ) # Print the phonon frequencies. @@ -114,6 +120,8 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse αeff = sum(α) if verbose + dprocess = 1 # Counter for drange iterations. + # Hide cursor whilst printing live data. print(io, "\e[?25l") @@ -121,6 +129,8 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse println(io, "\e[K Polaron Information: [$(αprocess[]) / $(num_α) ($(round(αprocess[] / (num_α) * 100, digits=1)) %)]") println(io, "\e[K-----------------------------------------------------------------------") + αprocess += 1 # Increment αrange iteration. + # Different formatting for single vs multiple frequencies (limits an array to head and tail to limit prints). if num_ω == 1 println(io, "\e[KPhonon frequencies | ωeff = ", ωeff, " | ω = ", ω) @@ -129,24 +139,30 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse println(io, "\e[KPhonon frequencies | ωeff = ", ωeff, " | ω = ", join(round.(first(ω, 2), digits=2), ", ")..., " ... ", join(round.(last(ω, 2), digits=2), ", ")...) println(io, "\e[KHolstein coupling | αeff = ", αeff, " | α = ", join(round.(first(α, 2), digits=3), ", ")..., " ... ", join(round.(last(α, 2), digits=3), ", ")...) end - println(io, "\e[KNumber of dimensions | d = ", dims) end + for d in 1:num_d + if verbose + println(io, "\e[KNumber of dimensions", " [", dprocess, " / ", num_d, "]", " | d = ", dims[d]) + dprocess += 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) = holstein_energy(v, w, α, ω; dims = dims) - v_gs, w_gs, F_gs, K_gs, P_gs = vw_variation(athermal_energy, v_guesses, w_guesses) + athermal_energy(v, w) = !kspace ? holstein_energy(v, w, α, ω; dims = dims[d]) : holstein_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 # Store the athermal data. - p["v0"][j, :] .= v_gs - p["w0"][j, :] .= w_gs - p["F0"][j] = F_gs - p["K0"][j] = K_gs - p["P0"][j] = P_gs + 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 + p["C0"][d, j] = C_gs # Print athermal variational parameter and energy data. if verbose @@ -156,15 +172,15 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse 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[KKinetic energy | K0 = ", K_gs) - println(io, "\e[KPotential energy | P0 = ", P_gs) + println(io, "\e[KElectron energy | A0 = ", A_gs) + println(io, "\e[KInteraction energy | B0 = ", B_gs) + println(io, "\e[KTrial energy | C0 = ", C_gs) Tprocess = 1 # Counter for Trange. - αprocess += 1 # Increment αrange iteration. end # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal κ_gs = (v_gs .^ 2 .- w_gs .^ 2) - p["κ0"][j, :] .= κ_gs + p["κ0"][d, j, :] .= κ_gs # Print athermal fictitious spring constant. if verbose @@ -173,7 +189,7 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse # Calculate and store fictitious masses. Athermal. M_gs = κ_gs ./ w_gs .^ 2 - p["M0"][j, :] .= M_gs + p["M0"][d, j, :] .= M_gs # Print athermal fictitious mass. if verbose @@ -182,7 +198,7 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Athermal M_asymp_gs = v_gs ./ w_gs - p["M0a"][j, :] .= M_asymp_gs + p["M0a"][d, j, :] .= M_asymp_gs # Print athermal asymptotic fictitious mass. if verbose @@ -191,7 +207,7 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. M_reduced_gs = (v_gs .^2 .- w_gs .^2) / v_gs .^ 2 - p["M0r"][j, :] .= M_reduced_gs + p["M0r"][d, j, :] .= M_reduced_gs # Print athermal reduced mass. if verbose @@ -200,7 +216,7 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse # 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 .* v_gs ./ (v_gs .^ 2 .- w_gs .^ 2) .^ 2) - p["R0"][j, :] .= R_gs + p["R0"][d, j, :] .= R_gs # Print athermal polaron radius. if verbose @@ -210,106 +226,131 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse for i in eachindex(Trange) # Temperatures loop. T = Trange[i] - # Print temperature. - if verbose - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[K Finite Temperature Information: [$(Tprocess[]) / $(num_T) ($(round(Tprocess[] / (num_T) * 100, digits=1)) %)]") - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[KTemperatures | T = ", T) - end + if !iszero(T) - # Calculate the reduced (unitless) thermodynamic betas for each phonon mode. - β = 1 ./ T .* β0 - p["β"][i, :] .= β + # Print temperature. + if verbose + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[K Finite Temperature Information: [$(Tprocess[]) / $(num_T) ($(round(Tprocess[] / (num_T) * 100, digits=1)) %)]") + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[KTemperatures | T = ", T) + end - # Print thermodynamic betas. - if verbose - println(io, "\e[KReduced thermodynamic | β = ", β) - end + # Calculate the reduced (unitless) thermodynamic betas for each phonon mode. + β = 1 ./ T .* β0 + p["β"][i, :] .= β - # 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) = holstein_energy(v, w, α, ω, β; dims = dims) - v, w, F, K, P = vw_variation(thermal_energy, v_guesses, w_guesses) + # Print thermodynamic betas. + if verbose + println(io, "\e[KReduced thermodynamic | β = ", β) + end - # Update the guesses to keep them close-ish to the true solutions during loops over temperatures. - v_guesses, w_guesses = v, w + # 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]) + v, w, F, A, B, C = vw_variation(thermal_energy, v_guesses, w_guesses) - # Store thermal data. - p["v"][i, j, :] .= v - p["w"][i, j, :] .= w - p["F"][i, j] = F - p["K"][i, j] = K - p["P"][i, j] = P + # Update the guesses to keep them close-ish to the true solutions during loops over temperatures. + v_guesses, w_guesses = v, w - # Print thermal data. - if verbose - println(io, "\e[KVariational parameter | v = ", v) - println(io, "\e[KVariational parameter | w = ", w) - println(io, "\e[KFree energy | F = ", F) - println(io, "\e[KKinetic energy | K = ", K) - println(io, "\e[KPotential energy | P = ", P) - end + # Store thermal data. + 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) + 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) + println(io, "\e[KTrial energy | C = ", C) + end - # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal - κ = (v .^ 2 .- w .^ 2) - p["κ"][i, j, :] .= κ + # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal + κ = (v .^ 2 .- w .^ 2) + p["κ"][i, d, j, :] .= κ - # Print spring constants. - if verbose - println(io, "\e[KFictitious spring constant | κ = ", κ) - end + # Print spring constants. + if verbose + println(io, "\e[KFictitious spring constant | κ = ", κ) + end - # Calculate and store fictitious masses. Thermal. - M = κ ./ w .^ 2 - p["M"][i, j, :] .= M + # Calculate and store fictitious masses. Thermal. + M = κ ./ w .^ 2 + p["M"][i, d, j, :] .= M - # Print masses. - if verbose - println(io, "\e[KFictitious mass | M = ", M) - end + # Print masses. + if verbose + println(io, "\e[KFictitious mass | M = ", M) + end - # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Thermal - M_asymp = v ./ w - p["Ma"][i, j, :] .= M_asymp + # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Thermal + M_asymp = v ./ w + p["Ma"][i, d, j, :] .= M_asymp - # Print asymptotic masses. - if verbose - println(io, "\e[KFictitious mass (asymptotic) | Ma = ", M_asymp) - end + # Print asymptotic masses. + if verbose + println(io, "\e[KFictitious mass (asymptotic) | Ma = ", M_asymp) + end - # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. - M_reduced = (v .^2 .- w .^2) / v .^ 2 - p["Mr"][i, j, :] .= M_reduced + # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. + M_reduced = (v .^2 .- w .^2) / v .^ 2 + p["Mr"][i, d, j, :] .= M_reduced - # Print redcued masses. - if verbose - println(io, "\e[KReduced mass | Mr = ", M_reduced) - end + # Print redcued masses. + if verbose + println(io, "\e[KReduced mass | Mr = ", M_reduced) + end - # Calculate and store polaron radii. - R = sqrt.(3 .* v ./ (v .^ 2 .- w .^ 2) .^ 2) - p["R"][i, j, :] .= R + # Calculate and store polaron radii. + R = sqrt.(3 .* v ./ (v .^ 2 .- w .^ 2) .^ 2) + p["R"][i, d, j, :] .= R - # Print polaron radius. - if verbose - println(io, "\e[KPolaron radius | R = ", R) - end + # Print polaron radius. + if verbose + println(io, "\e[KPolaron radius | R = ", R) + end - if verbose - println("\e[K-----------------------------------------------------------------------") - println("\e[K DC Mobility Information: ") - println("\e[K-----------------------------------------------------------------------") - end + if verbose + println("\e[K-----------------------------------------------------------------------") + println("\e[K DC Mobility Information: ") + println("\e[K-----------------------------------------------------------------------") + end - # Calculate and store the DC mobiliies. - μ = holstein_mobility(v, w, α, ω, β; dims = dims) / mb - p["μ"][i, j] = μ + # Calculate and store the DC mobiliies. + μ = !kspace ? holstein_mobility(v, w, α, ω, β; dims = dims[d]) / mb : holstein_mobility_k_space(v, w, α, ω, β; dims = dims[d]) / mb + p["μ"][i, d, j] = μ + + # Print DC mobilities. + if verbose + println(io, "\e[KFinite temperature mobility | μ = ", μ) + end + + else + # If zero temperature. + v, w, β = v_gs, w_gs, Inf + p["β"][i, :] .= β + 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 + p["C"][i, d, j] = C_gs + p["κ"][i, d, j, :] .= κ_gs + p["M"][i, d, j, :] .= M_gs + p["μ"][i, d, j] = Inf + p["R"][i, d, j, :] .= R_gs + p["Mr"][i, d, j, :] .= M_reduced_gs + p["Ma"][i, d, j, :] .= M_asymp_gs + end # End of zero temperature if statement. - # Print DC mobilities. if verbose - println(io, "\e[KFinite temperature mobility | μ = ", μ) Ωprocess = 1 # Counter for Ωrange. Tprocess += 1 # Increment Trange iterator. end @@ -317,50 +358,69 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse for k in eachindex(Ωrange) # E-field frequencies loop. Ω = Ωrange[k] - # Print E-field frequency. - if verbose - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[K Frequency Response Information: [$(Ωprocess[]) / $(num_Ω) ($(round(Ωprocess[] / (num_Ω) * 100, digits=1)) %)]") - println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[KElectric field frequency | Ω = ", Ω) - end + if !iszero(Ω) - # Calculate and store polaron memory functions (akin to self energy). - χ = holstein_memory_function(Ω, v, w, α, ω, β; dims = dims) - p["χ"][k, i, j] = χ + # Print E-field frequency. + if verbose + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[K Frequency Response Information: [$(Ωprocess[]) / $(num_Ω) ($(round(Ωprocess[] / (num_Ω) * 100, digits=1)) %)]") + println(io, "\e[K-----------------------------------------------------------------------") + println(io, "\e[KElectric field frequency | Ω = ", Ω) + end - # Print memory function. - if verbose - println(io, "\e[KMemory function | χ = ", χ) - 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]) + p["χ"][k, i, d, j] = χ - # Calculate and store polaron complex impedances. + # Print memory function. + if verbose + println(io, "\e[KMemory function | χ = ", χ) + end - z = -(im * Ω + im * χ) .* mb - p["z"][k, i, j] = z + # Calculate and store polaron complex impedances. - # Print complex impedances. - if verbose - println(io, "\e[KComplex impedance | z = ", z) - end + z = -(im * Ω + im * χ) .* mb + p["z"][k, i, d, j] = z + + # Print complex impedances. + if verbose + println(io, "\e[KComplex impedance | z = ", z) + end + + # Calculate and store polaron complex conductivities. + σ = 1 / z + p["σ"][k, i, d, j] = σ + + # Print complex conductivities and show total algorithm progress. + if verbose + println(io, "\e[KComplex conductivity | σ = ", σ) + end - # Calculate and store polaron complex conductivities. - σ = 1 / z - p["σ"][k, i, j] = σ + else + + # If zero frequency. + p["χ"][k, i, d, j] = Inf + 0 * im + p["z"][k, i, d, j] = 0 + im * Inf + p["σ"][k, i, d, j] = 0 + 0 * im + + end # End of zero frequency if statement. - # Print complex conductivities and show total algorithm progress. if verbose - println(io, "\e[KComplex conductivity | σ = ", σ) println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[K[Total Progress: $(process[]) / $(num_α * num_T * num_Ω) ($(round(process[] / (num_α * num_T * num_Ω) * 100, digits=1)) %)]") - print(io, "\e[9F") + println(io, "\e[K[Total Progress: $(process[]) / $(num_α * num_T * num_Ω * num_d) ($(round(process[] / (num_α * num_T * num_Ω * num_d) * 100, digits=1)) %)]") Ωprocess += 1 # Increment Ωrange iterator. process += 1 # Increment total iterator. + print(io, "\e[2F") end + + if verbose && !iszero(Ω) print(io, "\e[7F") end end - if verbose print(io, "\e[19F")end # Move up 19 lines and erase. + if verbose && !iszero(T) print(io, "\e[20F") end # Move up 20 lines and erase. end - if verbose print(io, "\e[19F") end # Move up 19 lines and erase. + if verbose print(io, "\e[15F") end # Move up 20 lines and erase. + end # End dimensions loop. + + if verbose print(io, "\e[5F") end end if verbose print("\e[?25h") end # Show cursor again. @@ -377,8 +437,9 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse p["v0"], # Athermal v params. p["w0"], # Athermal w params. p["F0"], # Athermal energies. - p["K0"], # Athermal A parameter. - p["P0"], # Athermal B parameter. + p["A0"], # Athermal A parameter. + p["B0"], # Athermal B parameter. + p["C0"], # Athermal C parameter. p["κ0"], # Fictitious spring constant. p["M0"], # Athermal fictitious mass. p["M0a"], # Athermal asymptotic approximate fictitious mass (v0/w0). @@ -387,8 +448,9 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse p["v"], # v params. p["w"], # w params. p["F"], # Energies. - p["K"], # A parameter. - p["P"], # B parameter. + p["A"], # A parameter. + p["B"], # B parameter. + p["C"], # C parameter. p["κ"], # Spring constants. p["M"], # Fictitious masses. p["Ma"], # Thermal asymptotic approximate fictitious mass (v/w). @@ -401,28 +463,28 @@ function holstein(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesse ] # Return Polaron type containing all generated data over any coupling strengths, temperatures and frequencies. - return Holstein(polaron_data...) + return HolsteinPolaron(polaron_data...) end """ Single alpha parameter. holstein() expects alpha parameters to be in a Vector. """ -holstein(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = holstein([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +holsteinpolaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No frequency input => Phonon peak. """ -holstein(αrange, Trange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = holstein(αrange, Trange, ω; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +holsteinpolaron(αrange, Trange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(αrange, Trange, 0; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No temperature input => 300 K. """ -holstein(αrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = holstein(αrange, 300, ω; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +holsteinpolaron(αrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(αrange, 0, 0; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No input => α = 1 """ -holstein(; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) = holstein(1, 300, ω; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) +holsteinpolaron(; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(1, 0, 0; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ holstein(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, verbose=false) @@ -430,7 +492,7 @@ Material specific constructors that use material specific parameters to paramete Material data is inputted through the `Material` type. Returns all data in either SI units or other common, suitable units otherwise. """ -function holstein(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) +function holsteinpolaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) # Show material data. if verbose @@ -443,14 +505,14 @@ function holstein(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.8 mb = material.mb # Generate Holstein polaron data from the arbitrary model constructor. - p = holstein(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, β0=ħ/kB*1e12*2π, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, verbose=verbose) + p = holsteinpolaron(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, β0=ħ/kB*1e12*2π, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) # Return material-specific, unitful Holstein type. return p end # Broadcast Polaron data. -function Base.show(io::IO, ::MIME"text/plain", x::Holstein) +function Base.show(io::IO, ::MIME"text/plain", x::HolsteinPolaron) io_lim = IOContext(io, :limit => true, :compact => true) @@ -469,8 +531,9 @@ function Base.show(io::IO, ::MIME"text/plain", x::Holstein) println(io_lim, "\e[KVariational parameter | v0 = ", x.v0) println(io_lim, "\e[KVariational parameter | w0 = ", x.w0) println(io_lim, "\e[KTotal energy | F0 = ", x.F0) - println(io_lim, "\e[KKinetic energy | K0 = ", x.K0) - println(io_lim, "\e[KPotential energy | P0 = ", x.P0) + println(io_lim, "\e[KElectron energy | A0 = ", x.A0) + println(io_lim, "\e[KInteraction energy | B0 = ", x.B0) + println(io_lim, "\e[KTrial energy | C0 = ", x.C0) println(io_lim, "\e[KFictitious spring constant | κ0 = ", x.κ0) println(io_lim, "\e[KFictitious mass | M0 = ", x.M0) println(io_lim, "\e[KFictitious mass (asymptotic) | M0a = ", x.M0a) @@ -486,8 +549,9 @@ function Base.show(io::IO, ::MIME"text/plain", x::Holstein) println(io_lim, "\e[KVariational parameter | v = ", x.v) println(io_lim, "\e[KVariational parameter | w = ", x.w) println(io_lim, "\e[KFree energy | F = ", x.F) - println(io_lim, "\e[KKinetic energy | K = ", x.K) - println(io_lim, "\e[KPotenital energy | P = ", x.P) + println(io_lim, "\e[KElectron energy | A = ", x.A) + println(io_lim, "\e[KInteraction energy | B = ", x.B) + println(io_lim, "\e[KTrial energy | C = ", x.C) println(io_lim, "\e[KFictitious spring constant | κ = ", x.κ) println(io_lim, "\e[KFictitious mass | M = ", x.M) println(io_lim, "\e[KFictitious mass (asymptotic) | Ma = ", x.Ma) @@ -519,7 +583,7 @@ Saves data from 'polaron' into file "prefix". This is a .jdl file for storing the polaron data whilst preserving types. Allows for saving multidimensional arrays that sometimes arise in the polaron data. Each parameter in the NewPolaron type is saved as a dictionary entry. E.g. NewPolaron.α is saved under JLD.load("prefix.jld")["alpha"]. """ -function save_holstein_polaron(polaron::Holstein, prefix) +function save_holstein_polaron(polaron::HolsteinPolaron, prefix) println("Saving polaron data to $prefix.jld ...") @@ -535,8 +599,9 @@ function save_holstein_polaron(polaron::Holstein, prefix) "v0", polaron.v0, "w0", polaron.w0, "athermal energy", polaron.F0, - "athermal K", polaron.K0, - "athermal P", polaron.P0, + "athermal A", polaron.A0, + "athermal B", polaron.B0, + "athermal C", polaron.C0, "athermal spring", polaron.κ0, "athermal mass", polaron.M0, "athermal asympt mass", polaron.M0a, @@ -545,8 +610,9 @@ function save_holstein_polaron(polaron::Holstein, prefix) "v", polaron.v, "w", polaron.w, "thermal energy", polaron.F, - "thermal K", polaron.K, - "thermal P", polaron.P, + "thermal A", polaron.A, + "thermal B", polaron.B, + "thermal C", polaron.C, "thermal spring", polaron.κ, "thermal mass", polaron.M, "thermal asympt mass", polaron.Ma, @@ -572,7 +638,7 @@ function load_holstein_polaron(polaron_file_path) data = JLD.load("$polaron_file_path") - holstein_polaron = Holstein( + holstein_polaron = HolsteinPolaron( data["alpha"], data["alpha eff"], data["temperature"], @@ -584,8 +650,9 @@ function load_holstein_polaron(polaron_file_path) data["v0"], data["w0"], data["athermal energy"], - data["athermal K"], - data["athermal P"], + data["athermal A"], + data["athermal B"], + data["athermal C"], data["athermal spring"], data["athermal mass"], data["athermal asympt mass"], @@ -594,8 +661,9 @@ function load_holstein_polaron(polaron_file_path) data["v"], data["w"], data["thermal energy"], - data["thermal K"], - data["thermal P"], + data["thermal A"], + data["thermal B"], + data["thermal C"], data["thermal spring"], data["thermal mass"], data["thermal asympt mass"], diff --git a/src/HolsteinTheory.jl b/src/HolsteinTheory.jl index 660e043..33c4960 100644 --- a/src/HolsteinTheory.jl +++ b/src/HolsteinTheory.jl @@ -1,197 +1,34 @@ # HolsteinTheory.jl +# Everything specifically associated with calculating properties of the Holstein polaron model. """ - polaron_propagator(τ, v, w, β) + holstein_coupling(k, α, ω; dims = 1) -Calculate the imaginary time polaron Green's function with temperature dependence. +Calculate the coupling strength for the Holstein lattice polaron model. # Arguments -- `τ`: A scalar representing the imaginary time. -- `v`: A scalar representing a variational parameter. -- `w`: A scalar representing a variational parameter. -- `β`: A scalar representing the inverse temperature. +- `k`: a scalar value representing the k-coordinate in k-space +- `α`: a scalar value representing the coupling constant +- `ω`: a scalar value representing the phonon frequency +- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. # Returns -The value of the polaron propagator. +The coupling strength for the Holstein model. # Example ```julia -τ = 0.5 -v = 0.2 -w = 0.1 -β = 1.0 -result = polaron_propagator(τ, v, w, β) -println(result) -``` -This example calculates the polaron propagator for given values of τ, v, w, and β. The result is then printed. -""" -function polaron_propagator(τ, v, w, β) - (v^2 - w^2) / v^3 * (1 - exp(-v * τ)) * (1 - exp(-v * (β - τ))) / (1 - exp(-v * β)) + w^2 / v^2 * τ * (1 - τ / β) + eps(Float64) -end - -""" - polaron_propagator(τ, v, w) - -Calculate the value of the polaron propagator based on the given inputs. - -# Arguments -- `τ::Number`: A scalar representing the imaginary time. -- `v::Number`: A scalar representing a variational parameter. -- `w::Number`: A scalar representing a variational parameter. - -# Returns -The value of the polaron propagator. - -# Example -```julia -τ = 0.5 -v = 0.2 -w = 0.1 -result = polaron_propagator(τ, v, w) -println(result) -``` -This example calculates the polaron propagator for the given values of τ, v, and w. The result is then printed. -""" -function polaron_propagator(τ, v, w) - w^2 * τ / v^2 + (v^2 - w^2) / v^3 * (1 - exp(-v * τ)) + eps(Float64) -end - -""" - phonon_propagator(τ, ω) - -Calculate the value of the phonon propagator based on the given inputs. - -# Arguments -- `τ`: A scalar representing the imaginary time. -- `ω`: A scalar representing the adiabaticity. - -# Example -```julia -τ = 0.5 -ω = 0.2 -result = phonon_propagator(τ, ω) -println(result) -``` -This example calculates the value of the phonon propagator for the given values of τ and ω. The result is then printed. - -# Returns -The value of the phonon propagator. -""" -function phonon_propagator(τ, ω) - exp(-ω * τ) -end - -""" - phonon_propagator(τ, ω, β) - -Calculate the value of the phonon propagator at a given imaginary time (`τ`), phonon frequency (`ω`), and inverse temperature (`β`). - -# Arguments -- `τ`: A scalar representing the imaginary time. -- `ω`: A scalar representing the phonon frequency. -- `β`: A scalar representing the inverse temperature. - -# Returns -The value of the phonon propagator. - -# Example -```julia -τ = 0.5 -ω = 0.2 -β = 1.0 -result = phonon_propagator(τ, ω, β) +result = holstein_coupling(2.0, 0.5, 1.0; dims = 3) println(result) ``` -This example calculates the value of the phonon propagator for `τ = 0.5`, `ω = 0.2`, and `β = 1.0`. The result is then printed. +Expected Output: +6.0 """ -function phonon_propagator(τ, ω, β) - n = 1 / (exp(β * ω) - 1) - result = n * exp(ω * τ) + (1 + n) * exp(-ω * τ) - if isnan(result) - return phonon_propagator(τ, ω) - else - return result - end -end - -""" - holstein_interaction_energy_integrand(τ, v, w, α, ω, β; dims = 3) - -Calculate the integrand for the Holstein interaction energy. - -## Arguments -- `τ`: A scalar representing the imaginary time. -- `v`: A scalar representing a variational parameter. -- `w`: A scalar representing a variational parameter. -- `α`: A scalar representing the electron-phonon coupling. -- `ω`: A scalar representing the adiabaticity. -- `β`: A scalar representing the inverse temperature. -- `dims`: An optional parameter representing the number of dimensions (default is 3). - -## Returns -The integrand for the Holstein interaction energy. - -## Example -```julia -τ = 0.5 -v = 0.2 -w = 0.1 -α = 0.3 -ω = 0.4 -β = 1.0 -result = holstein_interaction_energy_integrand(τ, v, w, α, ω, β) -println(result) -``` -This example calculates the integrand for the Holstein interaction energy using the given values of `τ`, `v`, `w`, `α`, and `ω`. The result is then printed. -""" -function holstein_interaction_energy_integrand(τ, v, w, α, ω, β; dims = 3) - if β == Inf - return holstein_interaction_energy_integrand(τ, v, w, α, ω; dims = dims) - end - momentum_cutoff = gamma(dims / 2 + 1)^(1/ dims) * (2√π) - coupling = holstein_coupling(1, α, ω; dims = dims) - propagator = polaron_propagator(τ, v, w, β * ω) / 2 - phonon_propagator(τ / ω, ω, β) * coupling * P(dims, momentum_cutoff^2 * propagator) / (4π * propagator)^(dims/2) -end - -""" - holstein_interaction_energy_integrand(τ, v, w, α, ω; dims = 3) - -Calculate the integrand for the Holstein interaction energy. - -## Arguments -- `τ`: A scalar representing the imaginary time. -- `v`: A scalar representing a variational parameter. -- `w`: A scalar representing a variational parameter. -- `α`: A scalar representing the electron-phonon coupling. -- `ω`: A scalar representing the adiabaticity. -- `dims`: An optional parameter representing the number of dimensions (default is 3). - -## Returns -The integrand for the Holstein interaction energy. - -## Example -```julia -τ = 0.5 -v = 0.2 -w = 0.1 -α = 0.3 -ω = 0.4 -result = holstein_interaction_energy_integrand(τ, v, w, α, ω) -println(result) -``` -This example calculates the integrand for the Holstein interaction energy using the given values of `τ`, `v`, `w`, `α`, and `ω`. The result is then printed. -""" - -function holstein_interaction_energy_integrand(τ, v, w, α, ω; dims = 3) - momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) - coupling = holstein_coupling(1, α, ω; dims = dims) / dims - propagator = polaron_propagator(τ, v, w) / 2 - phonon_propagator(τ / ω, ω) * coupling * P(dims, momentum_cutoff^2 * propagator) / (4π * propagator)^(dims/2) +function holstein_coupling(k, α, ω; dims = 3) + 2 * α * ω * dims end """ - holstein_interaction_energy(v, w, α, ω, β; dims = 3) + holstein_interaction_energy(v, w, α, ωβ...; dims = 3) Electron-phonon interaction energy for the Holstein mode at finite temperature. Here the k-space integral is evaluated analytically. @@ -218,42 +55,43 @@ println(result) ``` This example calculates the electron-phonon interaction energy for given values of `v`, `w`, `α`, `ω`, and `β`. The result is then printed. """ -function holstein_interaction_energy(v, w, α, ω, β; dims = 3) - if β == Inf - return holstein_interaction_energy(v, w, α, ω; dims = dims) - end - integral, _ = quadgk(τ -> holstein_interaction_energy_integrand(τ, v, w, α, ω, β; dims = dims), 0, β * ω / 2) - return integral +function holstein_interaction_energy(v, w, α, ωβ...; dims = 3) + momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) + coupling = holstein_coupling(1, α, ωβ[1]; dims = dims) + propagator(τ) = polaron_propagator(τ, v, w, ωβ...) / 2 + integrand(τ) = phonon_propagator(τ, ωβ...) * P(dims, momentum_cutoff^2 * propagator(τ)) / (propagator(τ))^(dims/2) + upper_limit = length(ωβ) == 1 ? Inf : prod(ωβ) / 2 + integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) + return integral * coupling / (4π)^(dims / 2) end """ - holstein_interaction_energy(v, w, α, ω; dims = 3) + holstein_interaction_energy_k_space(v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) -Calculate the interaction energy between electrons and phonons in the Holstein model. +Calculate the Holstein polaron interaction energy in k-space at finite temperaure. -# Arguments -- `v`: A scalar representing a variational parameter. -- `w`: A scalar representing a variational parameter. -- `α`: A scalar representing the electron-phonon coupling. -- `ω`: A scalar representing the adiabaticity. -- `dims`: An optional parameter representing the number of dimensions (default is 3). - -# Returns -- `integral`: The interaction energy between electrons and phonons in the Holstein model. +## Arguments +- `τ`: a scalar value representing the imaginary time. +- `v`: a scalar value representing a variational paramater. +- `w`: a scalar value representing a variational paramater. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `a`: an optional scalar value representing the lattice constant. Default value is 1. +- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. -# Example -```julia -v = 0.2 -w = 0.1 -α = 0.5 -ω = 0.3 -result = holstein_interaction_energy(v, w, α, ω) -println(result) -``` -This example calculates the interaction energy for given values of `v`, `w`, `α`, and `ω`. The result is then printed. +## Returns +A scalar value representing the Holstein polaron interaction energy in k-space at finite temperature. """ -function holstein_interaction_energy(v, w, α, ω; dims = 3) - integral, _ = quadgk(τ -> holstein_interaction_energy_integrand(τ, v, w, α, ω; dims = dims), 0, Inf) +function holstein_interaction_energy_k_space(v, w, α, ωβ...; dims = 3) + momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) + coupling(k) = holstein_coupling(k, α, ωβ[1]; dims = dims) + propagator(τ) = polaron_propagator(τ, v, w, ωβ...) + integrand(τ) = phonon_propagator(τ, ωβ...) * spherical_k_integral(coupling, propagator(τ); dims = dims, limits = [0, momentum_cutoff]) + upper_limit = length(ωβ) == 1 ? Inf : prod(ωβ) / 2 + integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) return integral end @@ -261,137 +99,36 @@ end Total free energy for the Holstein model. Here the k-space integral is evaluated analytically. """ function holstein_energy(v, w, α, ωβ...; dims = 3) - kinetic_energy = electron_energy(v, w, ωβ...; dims = dims) - potential_energy = -holstein_interaction_energy(v, w, α, ωβ...; dims = dims) - return (kinetic_energy + potential_energy), kinetic_energy, potential_energy -end - - -""" - vw_variation(energy, initial_v, initial_w; lower_bounds = [0, 0], upper_bounds = [Inf, Inf]) - -This function optimizes the values of `v` and `w` to minimize the energy function `energy(x[1] + x[2], x[2])[1]`. It uses the `Optim` package to perform the optimization and returns the optimized values of `v` and `w`, as well as the minimized energy, kinetic energy, and potential energy. - -## Arguments -- `energy`: A function that takes two arguments `x` and `y` and returns an array of energy components. -- `initial_v`: The initial value of `v`. -- `initial_w`: The initial value of `w`. -- `lower_bounds`: An optional array representing the lower bounds for `v` and `w` optimization. Default is `[0, 0]`. -- `upper_bounds`: An optional array representing the upper bounds for `v` and `w` optimization. Default is `[Inf, Inf]`. - -## Returns -- `Δv + w`: The optimized value of `v`. -- `w`: The optimized value of `w`. -- `energy_minimized`: The minimized energy. -- `kinetic_energy`: The kinetic energy corresponding to the minimized energy. -- `potential_energy`: The potential energy corresponding to the minimized energy. - -## Example -```julia -energy(x, y) = [x^2 + y^2, x^2, y^2] -initial_v = 0.5 -initial_w = 0.2 -lower_bounds = [0, 0] -upper_bounds = [Inf, Inf] -result = vw_variation(energy, initial_v, initial_w; lower_bounds, upper_bounds) -println(result) -``` -This example demonstrates how to use the `vw_variation` function. It defines an energy function `energy(x, y)` that returns an array of energy components. It then calls `vw_variation` with the initial values of `v` and `w`, as well as lower and upper bounds for the optimization. The function optimizes `v` and `w` to minimize the energy and returns the optimized values, as well as the minimized energy, kinetic energy, and potential energy. The result is then printed. -""" -function vw_variation(energy, initial_v, initial_w; lower = [0, 0], upper = [1e4, 1e4]) - - Δv = initial_v - initial_w # defines a constraint, so that v>w - initial = [Δv + eps(Float64), initial_w] - - # Feynman 1955 athermal action - f(x) = energy(x[1] + x[2], x[2])[1] - - # Use a more efficient optimization algorithm or library to optimise v and w to minimise enthalpy. - solution = Optim.optimize( - Optim.OnceDifferentiable(f, initial; autodiff=:forward), - lower, - upper, - initial, - Fminbox(BFGS()), - # Optim.Options(show_trace = false, g_tol = 1e-12) - ) - - # Get v and w values that minimise the free energy. - Δv, w = Optim.minimizer(solution) # v=Δv+w - energy_minimized = energy(Δv + w, w) - - if !Optim.converged(solution) - @warn "Failed to converge. v = $(Δv .+ w), w = $w, E = $energy_minimized" - end - - # Return variational parameters that minimise the free energy. - return Δv + w, w, energy_minimized... -end - -vw_variation(energy) = vw_variation(energy, 5, 3; lower = [0, 0], upper = [1e4, 1e4]) - -""" - general_memory_function(Ω, structure_factor; limits = [0, Inf]) - -This function calculates the integral of a given structure factor with respect to time using the `quadgk` function from the `QuadGK` package. - -## Arguments -- `Ω`: A scalar representing the frequency. -- `structure_factor`: A function that returns the value of the structure factor at a given time. -- `limits` (optional): A 2-element array specifying the lower and upper limits of integration. Default is `[0, Inf]`. - -## Returns -The integral of the structure factor with respect to time. - -## Example -```julia -# Define a structure factor function -function structure_factor(t) - # Implementation of the structure factor - # ... -end - -# Calculate the memory function for a given frequency and structure factor -Ω = 0.5 -limits = [0, 10] -result = general_memory_function(Ω, structure_factor; limits = limits) -println(result) -``` -This example demonstrates how to use the `general_memory_function` to calculate the memory function for a given frequency `Ω` and structure factor function `structure_factor`. The `limits` argument is optional and specifies the lower and upper limits of integration. The result is then printed. -""" -function general_memory_function(Ω, structure_factor; limits = [0, Inf]) - integral, _ = quadgk(t -> (1 - exp(im * Ω * t)) / Ω * imag(structure_factor(t)), limits[1], limits[2]) - return integral + A, C = trial_energy(v, w, ωβ...; dims = dims) + B = holstein_interaction_energy(v, w, α, ωβ...; dims = dims) + return -(A + B + C) - 2 * dims, A, B, C end """ - general_memory_function(structure_factor; limits = [0, Inf]) + holstein_energy_k_space(v, w, α, ωβ...; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) -Calculate the integral of a given function `structure_factor` using the `quadgk` function in Julia. +Calculate the total energy, kinetic energy, and interaction energy of the Holstein lattice polaron. ## Arguments -- `structure_factor`: A function that takes a single argument `t` and returns a value. -- `limits`: An optional array specifying the lower and upper limits of integration. Default is `[0, Inf]`. +- `v`: a scalar value representing a variational paramater. +- `w`: a scalar value representing a variational paramater. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: (optional) a scalar value representing the inverse temperature. +- `dims`: The number of dimensions of the system (default is 3). +- `rₚ`: The characteristic polaron radius (default is 1). +- `a`: The lattice constant (default is 1). +- `limits`: The limits of integration for the interaction energy calculation (default is [-π, π]). ## Returns -- `integral`: The result of the numerical integration of the function `structure_factor` over the specified limits. - -## Example -```julia -# Define the structure factor function -function structure_factor(t) - return t^2 + 2t + 1 -end - -# Call the general_memory_function with the structure_factor function -result = general_memory_function(structure_factor; limits = [0, 10]) - -println(result) # Output: 383.3333333333333 -``` -""" -function general_memory_function(structure_factor; limits = [0, Inf]) - integral, _ = quadgk(t -> -im * t * imag(structure_factor(t)), limits[1], limits[2]) - return integral +- `total_energy`: The calculated total polaron energy. +- `kinetic_energy`: The calculated polaron kinetic energy. +- `interaction_energy`: The calculated polaron interaction energy. +""" +function holstein_energy_k_space(v, w, α, ωβ...; dims = 3) + A, C = trial_energy(v, w, ωβ...; dims = dims) + B = holstein_interaction_energy_k_space(v, w, α, ωβ...; dims = dims) + return -(A + B + C) - 2 * dims, A, B, C end """ @@ -424,54 +161,42 @@ println(result) ``` This example demonstrates how to use the `holstein_structure_factor` function to calculate the structure factor for a given set of parameters. The function is called with the values of `t`, `v`, `w`, `α`, `ω`, and `β` as arguments, and the result is then printed. """ -function holstein_structure_factor(t, v, w, α, ω, β; dims = 3) - if β == Inf - return holstein_structure_factor(t, v, w, α, ω; dims = dims) - end +function holstein_structure_factor(t, v, w, α, ωβ...; dims = 3) momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) - coupling = holstein_coupling(1, α, ω; dims = dims) - propagator = polaron_propagator(im * t, v, w, β * ω) / 2 + coupling = holstein_coupling(1, α, ωβ[1]; dims = dims) + propagator = polaron_propagator(im * t, v, w, ωβ...) / 2 integral = dims / 2 * ball_surface(dims) / (2π)^3 * P_plus_one(dims, propagator * momentum_cutoff^2) / propagator^(dims/2 + 1) - phonon_propagator(im * t / ω, ω, β) * coupling * integral * 2 / dims / ω + return 2 / dims * phonon_propagator(im * t, ωβ...) * coupling * integral end """ - holstein_structure_factor(t, v, w, α, ω; dims = 3) + holstein_structure_factor_k_space(t, v, w, α, ωβ...; dims = 3) -Calculate the structure factor for the Holstein model. +Calculate the structure factor in k-space for the Holstein lattice polaron model at finite temperature. ## Arguments -- `t`: A scalar representing the time. -- `v`: A scalar representing a parameter. -- `w`: A scalar representing a parameter. -- `α`: A scalar representing a parameter. -- `ω`: A scalar representing a parameter. -- `dims`: An optional integer representing the number of dimensions (default is 3). +- `t`: a scalar value representing the real time. +- `v`: a scalar value representing a variational parameter. +- `w`: a scalar value representing a variational parameter. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `a`: an optional scalar value representing the lattice constant. Default value is 1. +- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. ## Returns -The value of the structure factor. - -## Example -```julia -t = 0.5 -v = 0.2 -w = 0.1 -α = 1.0 -ω = 0.5 -result = holstein_structure_factor(t, v, w, α, ω) -println(result) -``` -This example calculates the structure factor for the given values of `t`, `v`, `w`, `α`, and `ω`. The result is then printed. +A scalar value representing the calculated structure factor in k-space for the Holstein lattice polaron model at finite temperature. """ -function holstein_structure_factor(t, v, w, α, ω; dims = 3) - momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) - coupling = holstein_coupling(1, α, ω; dims = dims) - propagator = polaron_propagator(im * t, v, w) / 2 - integral = dims / 2 * ball_surface(dims) / (2π)^3 * P_plus_one(dims, propagator * momentum_cutoff^2) / propagator^(dims/2 + 1) - phonon_propagator(im * t / ω, ω) * coupling * integral * 2 / ω / dims +function holstein_structure_factor_k_space(t, v, w, α, ωβ...; dims = 3) + momentum_cutoff = gamma(dims / 2 + 1)^(1/ dims) * (2√π) + coupling(k) = holstein_coupling(k, α, ωβ[1]; dims = dims) * k^2 + propagator = polaron_propagator(im * t, v, w, ωβ...) + integral = spherical_k_integral(coupling_one, propagator; dims = dims, limits = [0, momentum_cutoff]) + return 2 / dims * phonon_propagator(im * t, ωβ...) * integral end - """ holstein_memory_function(Ω, v, w, α, ω, β; dims = 3) @@ -491,36 +216,34 @@ result = holstein_memory_function(Ω, v, w, α, ω, β; dims = 3) ``` In this example, the `holstein_memory_function` is called with the input parameters `Ω`, `v`, `w`, `α`, `ω`, and `β`, and the optional parameter `dims` set to 3. The function calculates the memory function using the `general_memory_function` function and returns the result. """ -function holstein_memory_function(Ω, v, w, α, ω, β; dims = 3) - if β == Inf - return holstein_memory_function(Ω, v, w, α, ω; dims = dims) - end - structure_factor(t) = holstein_structure_factor(t, v, w, α, ω, β; dims = dims) - return general_memory_function(Ω / ω, structure_factor) +function holstein_memory_function(Ω, v, w, α, ωβ...; dims = 3) + structure_factor(t) = holstein_structure_factor(t, v, w, α, ωβ...; dims = dims) + return polaron_memory_function(Ω / ωβ[1], structure_factor; limits = [0, 1e5]) end """ - holstein_memory_function(Ω, v, w, α, ω; dims = 3) + holstein_memory_function_k_space(Ω, v, w, α, ωβ...; dims = 3) -Calculate the integral of a structure factor using the `general_memory_function` function. +Calculate the memory function for the Holstein model in k-space at finite temperature and frequency. ## Arguments -- `Ω`: The frequency parameter. -- `v`, `w`, `α`, `ω`: Parameters used in the `holstein_structure_factor` function. -- `dims`: Optional parameter specifying the number of dimensions, defaults to 3. +- `Ω`: a scalar value representing the frequency at which the memory function is evaluated. +- `v`: a scalar value representing the optimal variational parameter. +- `w`: a scalar value representing the optimal variational parameter. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `a`: an optional scalar value representing the lattice constant. Default value is 1. +- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. ## Returns -The integral value of the structure factor. - -## Example -```julia -result = holstein_memory_function(Ω, v, w, α, ω; dims = 3) -``` -In this example, the `holstein_memory_function` is called with the parameters `Ω`, `v`, `w`, `α`, and `ω`. The `dims` parameter is optional and defaults to 3. The function calculates the structure factor using the `holstein_structure_factor` function and then calls the `general_memory_function` to calculate the integral of the structure factor. The result is stored in the `result` variable. +A scalar value representing the memory function of the Holstein model in k-space at finite temperature evaluated at the frequency `Ω`. """ -function holstein_memory_function(Ω, v, w, α, ω; dims = 3) - structure_factor(t) = holstein_structure_factor(t, v, w, α, ω; dims = dims) - return general_memory_function(Ω / ω, structure_factor, limits = [0, 1e6]) +function holstein_memory_function_k_space(Ω, v, w, α, ωβ...; dims = 3) + structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ωβ...; dims = dims) + return general_memory_function(Ω / ωβ[1], structure_factor; limits = [0, 1e5]) end """ @@ -558,7 +281,35 @@ function holstein_mobility(v, w, α, ω, β; dims = 3) return Inf end structure_factor(t) = holstein_structure_factor(t, v, w, α, ω, β; dims = dims) - abs(1 / imag(general_memory_function(structure_factor))) + abs(1 / imag(polaron_memory_function(structure_factor))) +end + +""" + holstein_mobility_k_space(v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) + +Calculate the DC mobility in k-space for a Holstein polaron system at finite temperature. + +## Arguments +- `v`: a scalar value representing the optimal variational parameter. +- `w`: a scalar value representing the optimal variational parameter. +- `α`: a scalar value representing the coupling constant. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. +- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. +- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. +- `a`: an optional scalar value representing the lattice constant. Default value is 1. +- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. + + +## Returns +The DC mobility in k-space for the Holstein polaron system at finite temperature. +""" +function holstein_mobility_k_space(v, w, α, ω, β; dims = 3) + if β == Inf + return Inf + end + structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = dims) + 1 / imag(general_memory_function(structure_factor)) end function holstein_complex_impedence(Ω, v, w, α, ωβ...; dims = 3) @@ -567,4 +318,5 @@ end function holstein_complex_conductivity(Ω, v, w, α, ωβ...; dims = 3) return 1 / holstein_complex_impedence(Ω, v, w, α, ωβ...; dims = dims) -end \ No newline at end of file +end + diff --git a/src/KSpaceTheory.jl b/src/KSpaceTheory.jl deleted file mode 100644 index 8a1fbcc..0000000 --- a/src/KSpaceTheory.jl +++ /dev/null @@ -1,771 +0,0 @@ -# KSpaceTheory.jl - -""" - cartesian_k_integrand(k, coupling, propagator; rₚ = 1) - -Calculate the integrand for a k-space integral in Cartesian coordinates. - -# Arguments -- `k`: a scalar value representing the k-coordinate in k-space -- `coupling`: a function that takes a scalar k value and returns a scalar value representing the coupling strength -- `propagator`: a scalar value representing the propagator -- `rₚ`: an optional scalar value representing the radius - -# Example Usage -```julia -coupling(k) = k^2 # define the coupling function -propagator = 0.5 # define the propagator value -result = cartesian_k_integrand(2.0, coupling, propagator; rₚ = 1) # calculate the integrand for k = 2.0 -println(result) # print the result -``` -Expected Output: -0.5 -""" -function cartesian_k_integrand(k, coupling, propagator) - coupling(k) * exp(-k^2 * propagator / 2) -end - - -""" - spherical_k_integrand(k, coupling, propagator; rₚ = 1) - -Calculate the integrand for a k-space integral in spherical coordinates. - -## Arguments -- `k`: a scalar value representing the k-coordinate in k-space -- `coupling`: a function that takes a scalar k value and returns a scalar value representing the coupling strength -- `propagator`: a scalar value representing the propagator -- `rₚ`: an optional scalar value representing the radius - -## Returns -The calculated integrand as a scalar value. - -## Example -```julia -coupling(k) = k^2 # define the coupling function -propagator = 0.5 # define the propagator value -result = spherical_k_integrand(2.0, coupling, propagator; rₚ = 1) # calculate the integrand for k = 2.0 -println(result) # print the result -``` -Expected Output: -`0.5` -""" -function spherical_k_integrand(k, coupling, propagator; dims = 3) - 2 * π^(dims/2) / gamma(dims/2) * k^(dims-1) * coupling(k) * exp(-k^2 * propagator / 2) -end - - -""" - cartesian_k_integral(coupling, propagator; rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the k-space integral in cartesian coordinates of the integrand `cartesian_k_integrand` over a specified range in k-space. - -# Arguments -- `coupling`: A function that takes a scalar k value and returns a scalar value representing the coupling strength. -- `propagator`: A scalar value representing the propagator. -- `rₚ`: An optional scalar value representing the charactersitc polaron radius. Default value is 1. -- `a`: An optional scalar value representing ther lattice constant. Default value is 1. -- `limits`: An array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. - -# Returns -A scalar value representing the calculated k-space integral over the specified range in cartesian coordinates. - -# Example -```julia -coupling(k) = k^2 # define the coupling function -propagator = 0.5 # define the propagator value -result = cartesian_k_integral(coupling, propagator; rₚ = 1, a = 1, limits = [-π, π]) # calculate the integral -println(result) # print the result -``` -Expected Output: -A scalar value representing the calculated k-space integral over the specified range in cartesian coordinates. -""" -function cartesian_k_integral(coupling, propagator; limits = [-π, π]) - integral, _ = quadgk(k -> cartesian_k_integrand(k, coupling, propagator), limits[1], limits[2]) - integral * a / 2π -end - - -""" - spherical_k_integral(coupling, propagator; rₚ = 1, limits = [0, Inf]) - -Calculate the k-space integral in spherical coordinates of the integrand `spherical_k_integrand` over a specified radius in k-space. - -# Arguments -- `coupling`: A function that takes a scalar k value and returns a scalar value representing the coupling strength. -- `propagator`: A scalar value representing the propagator. -- `rₚ`: An optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: An array of two scalar values representing the lower and upper limits of the integration range in spherical coordinates. Default is [0, Inf] for infinite sphere. - -# Returns -A scalar value representing the calculated integral over the specified range in spherical coordinates. - -# Example -```julia -coupling(k) = k^2 # define the coupling function -propagator = 0.5 # define the propagator value -result = spherical_k_integral(coupling, propagator; rₚ = 1, limits = [0, Inf]) # calculate the integral -println(result) # print the result -``` -Expected Output: -A scalar value representing the calculated k-space integral over the specified range in spherical coordinates. -""" -function spherical_k_integral(coupling, propagator; dims = 3, limits = [0, Inf]) - integral, _ = quadgk(k -> spherical_k_integrand(k, coupling, propagator; dims = dims), limits[1], limits[2]) - integral / (2π)^dims -end - - -""" - holstein_coupling(k, α, ω; dims = 1) - -Calculate the coupling strength for the Holstein lattice polaron model. - -# Arguments -- `k`: a scalar value representing the k-coordinate in k-space -- `α`: a scalar value representing the coupling constant -- `ω`: a scalar value representing the phonon frequency -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. - -# Returns -The coupling strength for the Holstein model. - -# Example -```julia -result = holstein_coupling(2.0, 0.5, 1.0; dims = 3) -println(result) -``` -Expected Output: -6.0 -""" -function holstein_coupling(k, α, ω; dims = 3) - 2 * α * ω * dims -end - - -""" - frohlich_coupling(k, α, ω) - -Calculate the coupling strength for the Frohlich continuum polaron model. - -# Arguments -- `k`: a scalar value representing the k-coordinate in k-space -- `α`: a scalar value representing the coupling constant -- `ω`: a scalar value representing the phonon frequency - -# Returns -The coupling strength for the Frohlich continuum polaron model. - -# Example -```julia -result = frohlich_coupling(2.0, 0.5, 1.0) -println(result) -``` -Expected Output: -`6.0` -""" -function frohlich_coupling(k, α, ω; mb = 1, dims = 3) - r_p = 1 / sqrt(2) - ω^2 * α * r_p * gamma((dims - 1) / 2) * (2√π / k)^(dims - 1) -end - -""" - holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the integrand for the Holstein interaction energy in k-space at finite temperature. - -# Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. - -# Returns -A scalar value representing the integrand for the Holstein interaction energy in k-space at finite temperature. -""" -function holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; dims = 3) - momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) - coupling(k) = holstein_coupling(k, α, ω; dims = dims) - propagator = polaron_propagator(τ, v, w, β) - phonon_propagator(τ, ω, β) * spherical_k_integral(coupling, propagator; dims = dims, limits = [0, momentum_cutoff]) -end - - -""" - holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the integrand for the Holstein interaction energy in k-space at zero temperature. - -# Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. - -# Returns -A scalar value representing the integrand for the Holstein interaction energy in k-space at zero temperature. -""" -function holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω; dims = 3) - momentum_cutoff = gamma(dims / 2 + 1)^(1/ dims) * (2√π) - coupling(k) = holstein_coupling(k, α, ω; dims = dims) - propagator = polaron_propagator(τ, v, w) - phonon_propagator(τ, ω) * spherical_k_integral(coupling, propagator; dims = dims, limits = [0, momentum_cutoff]) -end - - -""" - frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω, β) - -Calculate the integrand for the Frohlich interaction energy in k-space at finite temperaure. - -## Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [0, Inf]. - -## Returns -A scalar value representing the integrand for the Frohlich interaction energy in k-space at finite temperature. -""" -function frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; limits = [0, Inf], dims = 3) - coupling(k) = frohlich_coupling(k, α, ω; dims = dims) - propagator = polaron_propagator(τ, v, w, β) - phonon_propagator(τ, ω, β) * spherical_k_integral(coupling, propagator; limits = limits, dims = dims) -end - -""" - frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω; rₚ = 1, limits = [0, Inf]) - -Calculate the integrand for the Frohlich interaction energy in k-space at finite zero temperaure. - -## Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [0, Inf]. - -## Returns -A scalar value representing the integrand for the Frohlich interaction energy in k-space at zero temperature. -""" -function frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω; limits = [0, Inf], dims = 3) - coupling(k) = frohlich_coupling(k, α, ω; dims = dims) - propagator = polaron_propagator(τ, v, w) - phonon_propagator(τ, ω) * spherical_k_integral(coupling, propagator; limits = limits, dims = dims) -end - - -""" - holstein_interaction_energy_k_space(v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the Holstein polaron interaction energy in k-space at finite temperaure. - -## Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. - -## Returns -A scalar value representing the Holstein polaron interaction energy in k-space at finite temperature. -""" -function holstein_interaction_energy_k_space(v, w, α, ω, β; dims = 3) - interaction_energy, _ = quadgk(τ -> holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; dims = dims), 0, β/2) - return interaction_energy -end - - -""" - holstein_interaction_energy_k_space(v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the Holstein polaron interaction energy in k-space at zero temperaure. - -## Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. - -## Returns -A scalar value representing the Holstein polaron interaction energy in k-space at zero temperature. -""" -function holstein_interaction_energy_k_space(v, w, α, ω; dims = 3) - interaction_energy, _ = quadgk(τ -> holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω; dims = dims), 0, Inf) - return interaction_energy -end - -""" - frohlich_interaction_energy_k_space(v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) - -Calculate the Frohlich polaron interaction energy in k-space at finite temperaure. - -## Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. - -## Returns -A scalar value representing the Frohlich polaron interaction energy in k-space at finite temperature. -""" -function frohlich_interaction_energy_k_space(v, w, α, ω, β; limits = [0, Inf], dims = 3) - interaction_energy, _ = quadgk(τ -> frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; limits = limits, dims = dims), 0, β/2) - return interaction_energy -end - -""" - frohlich_interaction_energy_k_space(v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) - -Calculate the Frohlich polaron interaction energy in k-space at zero temperaure. - -## Arguments -- `τ`: a scalar value representing the imaginary time. -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an optional array of two scalar values representing the lower and upper limits of the integration range in k-space. Default value is [-π, π]. - -## Returns -A scalar value representing the Frohlich polaron interaction energy in k-space at zero temperature. -""" -function frohlich_interaction_energy_k_space(v, w, α, ω; limits = [0, Inf], dims = 3) - interaction_energy, _ = quadgk(τ -> frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω; limits = limits, dims = dims), 0, Inf) - return interaction_energy -end - - -""" - electron_energy(v, w, ω, β) - -Calculate the free electron energy at finite temperature. - -# Arguments -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. - -# Returns -A scalar value representing the calculated free electron energy at finite temperature. - -# Example -```julia -v = 0.5 -w = 1.0 -ω = 2.0 -β = 0.2 -result = electron_energy(v, w, ω, β; dims = 3) -println(result) -``` -Expected Output: -A scalar value representing the calculated free electron energy at finite temperature. -""" -function electron_energy(v, w, ω, β; dims = 3) - if β == Inf - return electron_energy(v, w, ω; dims = dims) - end - -(A(v, w, ω, β) + C(v, w, ω, β)) * dims / 3 -end - - -""" - electron_energy(v, w, ω) - -Calculate the free electron energy at zero temperature. - -# Arguments -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `ω`: a scalar value representing the phonon frequency. - -# Returns -A scalar value representing the calculated free electron energy at zero temperature. - -# Example -```julia -v = 0.5 -w = 1.0 -ω = 2.0 -β = 0.2 -result = electron_energy(v, w, ω) -println(result) -``` -Expected Output: -A scalar value representing the calculated free electron energy at finite temperature. -""" -function electron_energy(v, w, ω; dims = 3) - (v - w)^2 / (4 * v) * ω * dims -end - - -""" - holstein_energy_k_space(v, w, α, ωβ...; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the total energy, kinetic energy, and interaction energy of the Holstein lattice polaron. - -## Arguments -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: (optional) a scalar value representing the inverse temperature. -- `dims`: The number of dimensions of the system (default is 3). -- `rₚ`: The characteristic polaron radius (default is 1). -- `a`: The lattice constant (default is 1). -- `limits`: The limits of integration for the interaction energy calculation (default is [-π, π]). - -## Returns -- `total_energy`: The calculated total polaron energy. -- `kinetic_energy`: The calculated polaron kinetic energy. -- `interaction_energy`: The calculated polaron interaction energy. -""" -function holstein_energy_k_space(v, w, α, ωβ...; dims = 3) - kinetic_energy = electron_energy(v, w, ωβ...; dims = dims) - interaction_energy = -holstein_interaction_energy_k_space(v, w, α, ωβ...; dims = dims) - return kinetic_energy + interaction_energy, kinetic_energy, interaction_energy -end - - -""" - frohlich_energy_k_space(v, w, α, ωβ...; rₚ = 1, limits = [0, Inf]) - -Calculate the total energy, kinetic energy, and interaction energy of the Frohlich lattice polaron. - -## Arguments -- `v`: a scalar value representing a variational paramater. -- `w`: a scalar value representing a variational paramater. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: (optional) a scalar value representing the inverse temperature. -- `rₚ`: The characteristic polaron radius (default is 1). -- `limits`: The limits of integration for the interaction energy calculation (default is [0, Inf]). - -## Returns -- `total_energy`: The calculated total polaron energy. -- `kinetic_energy`: The calculated polaron kinetic energy. -- `interaction_energy`: The calculated polaron interaction energy. -""" -function frohlich_energy_k_space(v, w, α, ωβ...; limits = [0, Inf], dims = 3) - kinetic_energy = electron_energy(v, w, ωβ...; dims = dims) - interaction_energy = -frohlich_interaction_energy_k_space(v, w, α, ωβ...; limits = limits, dims = dims) - return kinetic_energy + interaction_energy, kinetic_energy, interaction_energy -end - - -""" - holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the structure factor in k-space for the Holstein lattice polaron model at finite temperature. - -## Arguments -- `t`: a scalar value representing the real time. -- `v`: a scalar value representing a variational parameter. -- `w`: a scalar value representing a variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. - -## Returns -A scalar value representing the calculated structure factor in k-space for the Holstein lattice polaron model at finite temperature. -""" -function holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = 3) - - momentum_cutoff = gamma(dims / 2 + 1)^(1/ dims) * (2√π) - - coupling(k) = holstein_coupling(k, α, ω; dims = dims) * k^2 - - propagator = polaron_propagator(im * t, v, w, β) - - integral = spherical_k_integral(coupling_one, propagator; limits = [0, momentum_cutoff]) - - prefactor = 2 / 3 * phonon_propagator(im * t, ω, β) - - prefactor * integral -end - - -""" - holstein_structure_factor_k_space(t, v, w, α, ω; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the structure factor in k-space for the Holstein lattice polaron model at zero temperature. - -## Arguments -- `t`: a scalar value representing the real time. -- `v`: a scalar value representing a variational parameter. -- `w`: a scalar value representing a variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. - -## Returns -A scalar value representing the calculated structure factor in k-space for the Holstein lattice polaron model at zero temperature. -""" -function holstein_structure_factor_k_space(t, v, w, α, ω; dims = 3) - - momentum_cutoff = gamma(dims / 2 + 1)^(1/ dims) * (2√π) - - coupling(k) = holstein_coupling(k, α, ω; dims = dims) * k^2 - - propagator = polaron_propagator(im * t, v, w) - - integral = spherical_k_integral(coupling, propagator; limits = [0, momentum_cutoff]) - - prefactor = 2 / 3 * phonon_propagator(im * t, ω) - - prefactor * integral -end - - -""" - frohlich_structure_factor_k_space(t, v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) - -Calculate the structure factor in k-space for the Frohlich continuum polaron model at finite temperature. - -## Arguments -- `t`: a scalar value representing the real time. -- `v`: a scalar value representing a variational parameter. -- `w`: a scalar value representing a variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. - -## Returns -A scalar value representing the calculated structure factor in k-space for the Frohlich continuum polaron model at finite temperature. -""" -function frohlich_structure_factor_k_space(t, v, w, α, ω, β; limits = [0, Inf], dims = 3) - - coupling(k) = frohlich_coupling(k, α, ω; dims = dims) * k^2 - - propagator = polaron_propagator(im * t, v, w, β) - - integral = spherical_k_integral(coupling, propagator; limits = limits, dims = dims) - - prefactor = 2 / 3 * phonon_propagator(im * t, ω, β) - - prefactor * integral -end - - -""" - frohlich_structure_factor_k_space(t, v, w, α, ω; rₚ = 1, limits = [0, Inf]) - -Calculate the structure factor in k-space for the Frohlich continuum polaron model at zero temperature. - -## Arguments -- `t`: a scalar value representing the real time. -- `v`: a scalar value representing the optimal variational parameter. -- `w`: a scalar value representing the optimal variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. - -## Returns -A scalar value representing the calculated structure factor in k-space for the Frohlich continuum polaron model at zero temperature. -""" -function frohlich_structure_factor_k_space(t, v, w, α, ω; limits = [0, Inf], dims = 3) - - coupling(k) = frohlich_coupling(k, α, ω; dims = dims) * k^2 - - propagator = polaron_propagator(im * t, v, w) - - integral = spherical_k_integral(coupling, propagator; limits = limits, dims = dims) - - prefactor = 2 / 3 * phonon_propagator(im * t, ω) - - prefactor * integral -end - - -""" - holstein_memory_function_k_space(Ω, v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the memory function for the Holstein model in k-space at finite temperature and frequency. - -## Arguments -- `Ω`: a scalar value representing the frequency at which the memory function is evaluated. -- `v`: a scalar value representing the optimal variational parameter. -- `w`: a scalar value representing the optimal variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. - -## Returns -A scalar value representing the memory function of the Holstein model in k-space at finite temperature evaluated at the frequency `Ω`. -""" -function holstein_memory_function_k_space(Ω, v, w, α, ω, β; dims = 3) - structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = dims) - return general_memory_function(Ω, structure_factor) -end - - -""" - holstein_memory_function_k_space(Ω, v, w, α, ω; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the memory function for the Holstein model in k-space at zero temperature and finite frequency. - -## Arguments -- `Ω`: a scalar value representing the frequency at which the memory function is evaluated. -- `v`: a scalar value representing the optimal variational parameter. -- `w`: a scalar value representing the optimal variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. - -## Returns -A scalar value representing the memory function of the Holstein model in k-space at zero temperature evaluated at the frequency `Ω`. -""" -function holstein_memory_function_k_space(Ω, v, w, α, ω; dims = 3) - structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ω; dims = dims) - return general_memory_function(Ω, structure_factor, limits = [0, 1e4]) -end - - -""" - holstein_mobility_k_space(v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - -Calculate the DC mobility in k-space for a Holstein polaron system at finite temperature. - -## Arguments -- `v`: a scalar value representing the optimal variational parameter. -- `w`: a scalar value representing the optimal variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `dims`: an optional scalar value representing the dimensionality of the system. Default value is 3. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `a`: an optional scalar value representing the lattice constant. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is 1D Brillouin zone. - - -## Returns -The DC mobility in k-space for the Holstein polaron system at finite temperature. -""" -function holstein_mobility_k_space(v, w, α, ω, β; dims = 3) - structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = dims) - 1 / imag(general_memory_function(structure_factor)) -end - - -""" - frohlich_memory_function_k_space(Ω, v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) - -Calculate the memory function for the Frohlich model in k-space at finite temperature and frequency. - -## Arguments -- `Ω`: a scalar value representing the frequency at which the memory function is evaluated. -- `v`: a scalar value representing the optimal variational parameter. -- `w`: a scalar value representing the optimal variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. - -## Returns -A scalar value representing the memory function of the Frohlich model in k-space at finite temperature evaluated at the frequency `Ω`. -""" -function frohlich_memory_function_k_space(Ω, v, w, α, ω, β; limits = [0, Inf]) - structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ω, β; limits = limits) - return general_memory_function(Ω, structure_factor; limits = [0, Inf]) -end - - -""" - frohlich_memory_function_k_space(Ω, v, w, α, ω; rₚ = 1, limits = [0, Inf]) - -Calculate the memory function for the Frohlich model in k-space at zero temperature and finite frequency. - -## Arguments -- `Ω`: a scalar value representing the frequency at which the memory function is evaluated. -- `v`: a scalar value representing the optimal variational parameter. -- `w`: a scalar value representing the optimal variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. - -## Returns -A scalar value representing the memory function of the Frohlich model in k-space at finite temperature evaluated at the frequency `Ω`. -""" -function frohlich_memory_function_k_space(Ω, v, w, α, ω; limits = [0, Inf]) - structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ω; limits = limits) - return general_memory_function(Ω, structure_factor; limits = [0, 1e4]) -end - - -""" - frohlich_mobility_k_space(v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) - -Calculate the DC mobility in k-space for a Frohlich polaron system at finite temperature. - -## Arguments -- `v`: a scalar value representing the optimal variational parameter. -- `w`: a scalar value representing the optimal variational parameter. -- `α`: a scalar value representing the coupling constant. -- `ω`: a scalar value representing the phonon frequency. -- `β`: a scalar value representing the inverse temperature. -- `rₚ`: an optional scalar value representing the characteristic polaron radius. Default value is 1. -- `limits`: an array of two scalar values representing the lower and upper limits of the integration range in k-space. Default is an infinite sphere. - - -## Returns -The DC mobility in k-space for the Frohlich polaron system at finite temperature. -""" -function frohlich_mobility_k_space(v, w, α, ω, β; limits = [0, Inf]) - structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ω, β; limits = limits) - 1 / imag(general_memory_function(structure_factor; limits = [0, 1e5])) -end - diff --git a/src/LegacyFunctions.jl b/src/LegacyFunctions.jl new file mode 100644 index 0000000..aeed91c --- /dev/null +++ b/src/LegacyFunctions.jl @@ -0,0 +1,91 @@ +# LegacyFunctions.jl +# Some old functions that are still potentially useful and have use in some tests. + +""" + feynmanvw(v, w, αωβ...; upper_limit = Inf64) + +Minimises the multiple phonon mode free energy function for a set of vₚ and wₚ variational parameters. The variational parameters follow the inequality: v₁ > w₁ > v₂ > w₂ > ... > vₙ > wₙ. Generalises `feynmanvw` to multiple variational parameters. + +# Arguments +- `v::Vector{Float64}`: vector of initial v parameters. +- `w::Vector{Float64}`: vector of initial w parameters. +- `α::Union{Float64, Vector{Float64}}`: is the partial dielectric electron-phonon coupling parameter for one or many phonon modes. +- `ω::Union{Float64, Vector{Float64}}`: phonon mode frequencies (2π THz) for one or many phonon modes. +- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) for one or many phonon modes. + +See also [`F`](@ref). +""" +function feynmanvw(v::Vector, w::Vector, αωβ...; upper_limit=1e6) + + if length(v) != length(w) + return error("The number of variational parameters v & w must be equal.") + end + + N_params = length(v) + + Δv = v .- w + initial = vcat(Δv .+ repeat([eps(Float64)], N_params), w) + + # Limits of the optimisation. + lower = fill(0.0, 2 * N_params) + upper = fill(upper_limit, 2 * N_params) + + # The multiple phonon mode free energy function to minimise. + f(x) = frohlich_energy([x[2*n-1] for n in 1:N_params] .+ [x[2*n] for n in 1:N_params], [x[2*n] for n in 1:N_params], αωβ...)[1] + + # Use Optim to optimise the free energy function w.r.t the set of v and w parameters. + solution = Optim.optimize( + Optim.OnceDifferentiable(f, initial; autodiff=:forward), + lower, + upper, + initial, + Fminbox(BFGS()) + ) + + # Extract the v and w parameters that minimised the free energy. + var_params = Optim.minimizer(solution) + + # Separate the v and w parameters into one-dimensional arrays (vectors). + Δv = [var_params[2*n-1] for n in 1:N_params] + w = [var_params[2*n] for n in 1:N_params] + E, A, B, C = frohlich_energy(Δv .+ w, w, αωβ...) + + # if Optim.converged(solution) == false + # @warn "Failed to converge. v = $(Δv .+ w), w = $w, E = $E" + # end + + # Return the variational parameters that minimised the free energy. + return Δv .+ w, w, E, A, B, C +end + +function feynmanvw(v::Real, w::Real, αωβ...; upper = [Inf, Inf], lower = [0, 0]) + + Δv = v .- w + initial = [Δv + eps(Float64), w] + + # The multiple phonon mode free energy function to minimise. + f(x) = frohlich_energy(x[1] .+ x[2], x[2], αωβ...)[1] + + # Use Optim to optimise the free energy function w.r.t the set of v and w parameters. + solution = Optim.optimize( + Optim.OnceDifferentiable(f, initial; autodiff=:forward), + lower, + upper, + initial, + Fminbox(BFGS()) + ) + + # Extract the v and w parameters that minimised the free energy. + Δv, w = Optim.minimizer(solution) + E, A, B, C = frohlich_energy(Δv .+ w, w, αωβ...) + + if Optim.converged(solution) == false + @warn "Failed to converge. v = $(Δv .+ w), w = $w, E = $E" + end + + # Return the variational parameters that minimised the free energy. + return Δv .+ w, w, E, A, B, C +end + +feynmanvw(αωβ...) = feynmanvw(3.4, 2.6, αωβ...) + diff --git a/src/Material.jl b/src/Material.jl index 39f982f..3781af9 100644 --- a/src/Material.jl +++ b/src/Material.jl @@ -1,4 +1,5 @@ # Material.jl +# Calculates inpit variables for the polaron theories from material data and saves the material as a Material type. """ mutable struct Material diff --git a/src/MemoryFunction.jl b/src/MemoryFunction.jl deleted file mode 100644 index 3d5835a..0000000 --- a/src/MemoryFunction.jl +++ /dev/null @@ -1,39 +0,0 @@ -# MemoryFunction.jl -# References: -# [1] R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Mobility of slow electrons in a polar crystal, PhysicalReview127, 1004 (1962) https://doi.org/10.1103/PhysRev.127.1004. -# [2] Devreese, J. De Sitter, and M. Goovaerts, “Optical absorption of polarons in thefeynman-hellwarth-iddings-platzman approximation,”Phys. Rev. B, vol. 5, pp. 2367–2381, Mar 1972. -# [3] F. Peeters and J. Devreese, “Theory of polaron mobility,” in Solid State Physics,pp. 81–133, Elsevier, 1984. - -function frohlich_structure_factor(t, v, w, α, ω, β; dims = 3) - if β == Inf - return frohlich_structure_factor(t, v, w, α, ω; dims = dims) - end - coupling = frohlich_coupling(1, α, ω; dims = dims) - propagator = polaron_propagator(im * t, v, w, β * ω) / 2 - integral = ball_surface(dims) / (2π)^dims * √π / 4 / propagator^(3/2) - return coupling * integral * phonon_propagator(im * t / ω, ω, β) * 2 / dims / ω -end - -function frohlich_structure_factor(t, v, w, α, ω; dims = 3) - coupling = frohlich_coupling(1, α, ω; dims = dims) - propagator = polaron_propagator(im * t, v, w) / 2 - integral = ball_surface(dims) / (2π)^dims * √π / 4 / propagator^(3/2) - return coupling * integral * phonon_propagator(im * t / ω, ω) * 2 / dims / ω -end - -function frohlich_memory_function(Ω, v, w, α, ω, β; dims = 3) - if β == Inf - return frohlich_memory_function(Ω, v, w, α, ω; dims = dims) - end - structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω, β; dims = dims) - return general_memory_function(Ω / ω, structure_factor) -end - -function frohlich_memory_function(Ω, v, w, α, ω; dims = 3) - structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω; dims = dims) - return general_memory_function(Ω / ω, structure_factor, limits = [0, 1e6]) -end - -frohlich_memory_function(Ω, v::Vector, w::Vector, α, ωβ...) = sum(frohlich_memory_function.(Ω, v, w, α, ωβ...)) - -frohlich_memory_function(Ω, v, w, α::Vector, ωβ...) = sum(frohlich_memory_function.(Ω, v, w, α, ωβ...)) diff --git a/src/PolaronFunctions.jl b/src/PolaronFunctions.jl new file mode 100644 index 0000000..cadbb02 --- /dev/null +++ b/src/PolaronFunctions.jl @@ -0,0 +1,223 @@ +# PolaronFunctions.jl + +""" + phonon_propagator(τ, ω, β) + +Calculate the value of the phonon propagator at a given imaginary time (`τ`), phonon frequency (`ω`), and inverse temperature (`β`). + +# Arguments +- `τ`: A scalar representing the imaginary time. +- `ω`: A scalar representing the phonon frequency. +- `β`: A scalar representing the inverse temperature. + +# Returns +The value of the phonon propagator. + +# Example +```julia +τ = 0.5 +ω = 0.2 +β = 1.0 +result = phonon_propagator(τ, ω, β) +println(result) +``` +This example calculates the value of the phonon propagator for `τ = 0.5`, `ω = 0.2`, and `β = 1.0`. The result is then printed. +""" +function phonon_propagator(τ, ω, β) + n = 1 / (exp(β * ω) - 1) + result = n * exp(τ) + (1 + n) * exp(-τ) + if isnan(result) + return phonon_propagator(τ, ω) + else + return result / ω + end +end + +""" + phonon_propagator(τ, ω) + +Calculate the value of the phonon propagator based on the given inputs. + +# Arguments +- `τ`: A scalar representing the imaginary time. +- `ω`: A scalar representing the adiabaticity. + +# Example +```julia +τ = 0.5 +ω = 0.2 +result = phonon_propagator(τ, ω) +println(result) +``` +This example calculates the value of the phonon propagator for the given values of τ and ω. The result is then printed. + +# Returns +The value of the phonon propagator. +""" +function phonon_propagator(τ, ω) + exp(-τ) / ω +end + +""" + vw_variation(energy, initial_v, initial_w; lower_bounds = [0, 0], upper_bounds = [Inf, Inf]) + +This function optimizes the values of `v` and `w` to minimize the energy function `energy(x[1] + x[2], x[2])[1]`. It uses the `Optim` package to perform the optimization and returns the optimized values of `v` and `w`, as well as the minimized energy, kinetic energy, and potential energy. + +## Arguments +- `energy`: A function that takes two arguments `x` and `y` and returns an array of energy components. +- `initial_v`: The initial value of `v`. +- `initial_w`: The initial value of `w`. +- `lower_bounds`: An optional array representing the lower bounds for `v` and `w` optimization. Default is `[0, 0]`. +- `upper_bounds`: An optional array representing the upper bounds for `v` and `w` optimization. Default is `[Inf, Inf]`. + +## Returns +- `Δv + w`: The optimized value of `v`. +- `w`: The optimized value of `w`. +- `energy_minimized`: The minimized energy. +- `kinetic_energy`: The kinetic energy corresponding to the minimized energy. +- `potential_energy`: The potential energy corresponding to the minimized energy. + +## Example +```julia +energy(x, y) = [x^2 + y^2, x^2, y^2] +initial_v = 0.5 +initial_w = 0.2 +lower_bounds = [0, 0] +upper_bounds = [Inf, Inf] +result = vw_variation(energy, initial_v, initial_w; lower_bounds, upper_bounds) +println(result) +``` +This example demonstrates how to use the `vw_variation` function. It defines an energy function `energy(x, y)` that returns an array of energy components. It then calls `vw_variation` with the initial values of `v` and `w`, as well as lower and upper bounds for the optimization. The function optimizes `v` and `w` to minimize the energy and returns the optimized values, as well as the minimized energy, kinetic energy, and potential energy. The result is then printed. +""" +function vw_variation(energy, initial_v, initial_w; lower = [0, 0], upper = [1e4, 1e4], warn = false) + + Δv = initial_v - initial_w # defines a constraint, so that v>w + initial = [Δv + eps(Float64), initial_w] + + # Feynman 1955 athermal action + f(x) = energy(x[1] + x[2], x[2])[1] + + # Use a more efficient optimization algorithm or library to optimise v and w to minimise enthalpy. + solution = Optim.optimize( + Optim.OnceDifferentiable(f, initial; autodiff=:forward), + lower, + upper, + initial, + Fminbox(BFGS()), + # Optim.Options(show_trace = false, g_tol = 1e-12) + ) + + # Get v and w values that minimise the free energy. + Δv, w = Optim.minimizer(solution) # v=Δv+w + energy_minimized = energy(Δv + w, w) + + if !Optim.converged(solution) && warn + @warn "Failed to converge. v = $(Δv .+ w), w = $w, E = $energy_minimized" + end + + # Return variational parameters that minimise the free energy. + return Δv + w, w, energy_minimized... +end + +vw_variation(energy) = vw_variation(energy, 5, 3; lower = [0, 0], upper = [1e4, 1e4]) + +""" + polaron_memory_function(Ω, structure_factor; limits = [0, Inf]) + +This function calculates the integral of a given structure factor with respect to time using the `quadgk` function from the `QuadGK` package. + +## Arguments +- `Ω`: A scalar representing the frequency. +- `structure_factor`: A function that returns the value of the structure factor at a given time. +- `limits` (optional): A 2-element array specifying the lower and upper limits of integration. Default is `[0, Inf]`. + +## Returns +The integral of the structure factor with respect to time. + +## Example +```julia +# Define a structure factor function +function structure_factor(t) + # Implementation of the structure factor + # ... +end + +# Calculate the memory function for a given frequency and structure factor +Ω = 0.5 +limits = [0, 10] +result = general_memory_function(Ω, structure_factor; limits = limits) +println(result) +``` +This example demonstrates how to use the `general_memory_function` to calculate the memory function for a given frequency `Ω` and structure factor function `structure_factor`. The `limits` argument is optional and specifies the lower and upper limits of integration. The result is then printed. +""" +function polaron_memory_function(Ω, structure_factor; limits = [0, Inf]) + integral, _ = quadgk(t -> (1 - exp(im * Ω * t)) / Ω * imag(structure_factor(t)), limits[1], limits[2]) + return integral +end + +""" + general_memory_function(structure_factor; limits = [0, Inf]) + +Calculate the integral of a given function `structure_factor` using the `quadgk` function in Julia. + +## Arguments +- `structure_factor`: A function that takes a single argument `t` and returns a value. +- `limits`: An optional array specifying the lower and upper limits of integration. Default is `[0, Inf]`. + +## Returns +- `integral`: The result of the numerical integration of the function `structure_factor` over the specified limits. + +## Example +```julia +# Define the structure factor function +function structure_factor(t) + return t^2 + 2t + 1 +end + +# Call the general_memory_function with the structure_factor function +result = general_memory_function(structure_factor; limits = [0, 10]) + +println(result) # Output: 383.3333333333333 +``` +""" +function polaron_memory_function(structure_factor; limits = [0, Inf]) + integral, _ = quadgk(t -> -im * t * imag(structure_factor(t)), limits[1], limits[2]) + return integral +end + +# Utility / general use functions + +# QOL function for removing singleton dimensions. +reduce_array(a) = length(a) == 1 ? only(a) : Array(dropdims(a, dims=tuple(findall(size(a) .== 1)...))) + +# Regularized lower gamma function for half dimensions n/2. +function P(n, z) + if n == 1 + erf(sqrt(z)) + elseif n == 2 + 1 - exp(-z) + elseif n > 2 + P_plus_one(n, z) + end +end + +# Recursion relations for higher half dimensions n/2 + m where m ∈ Integer. +function P_plus_one(n, z) + if n<3 + return P(n, z) + else + return P(n-2, z) - exp(-z) * z^(n/2-1) / gamma(n/2) + end +end + +# Surface area of an n-dimensional sphere. +function ball_surface(n) + 2 * π^(n/2) / gamma(n/2) +end + +# n-dimensional spherical integral for polaron theory. +function spherical_k_integral(coupling, propagator; dims = 3, limits = [0, Inf]) + integrand(k) = k^(dims-1) * coupling(k) * exp(-k^2 * propagator / 2) + integral, _ = quadgk(k -> integrand(k), limits[1], limits[2]) + return integral * ball_surface(dims) / (2π)^dims +end \ No newline at end of file diff --git a/src/PolaronMobility.jl b/src/PolaronMobility.jl index 6219b83..60b17de 100644 --- a/src/PolaronMobility.jl +++ b/src/PolaronMobility.jl @@ -8,13 +8,13 @@ PolaronMobility.jl - https://github.com/jarvist/PolaronMobility.jl module PolaronMobility # Type to hold the Frohlich polaron data. -export Polaron, polaron +export FrohlichPolaron, frohlichpolaron # Type to hold material specific data. export Material, material # Type to hold Holstein polaron data. -export Holstein, holstein +export HolsteinPolaron, holsteinpolaron # Frohlich alpha, energy and minimisation functions. export frohlichalpha, ϵ_ionic_mode, frohlich_energy, A, B, C, feynmanvw @@ -32,21 +32,21 @@ export frohlich_complex_impedence, frohlich_complex_conductivity, holstein_compl export frohlich_mobility, Hellwarth_mobility, Kadanoff_mobility_lowT, FHIP_mobility_lowT # Functions to save and load Polaron types to/from .jld format -export save_polaron, load_polaron, save_holstein_polaron, load_holstein_polaron +export save_frohlich_polaron, load_frohlich_polaron, save_holstein_polaron, load_holstein_polaron # Polaron effective mass export polaron_effective_mass # Generalised polaron functions -export general_memory_function, vw_variation, polaron_propagator, phonon_propagator, electron_energy +export polaron_memory_function, vw_variation, polaron_propagator, phonon_propagator, trial_energy # Holstein functions -export holstein_interaction_energy_integrand, holstein_interaction_energy, holstein_energy, holstein_structure_factor, holstein_memory_function, holstein_mobility +export holstein_interaction_energy, holstein_energy, holstein_structure_factor, holstein_memory_function, holstein_mobility export frohlich_structure_factor # K-space functions -export cartesian_k_integrand, spherical_k_integral, holstein_coupling, frohlich_coupling, holstein_interaction_energy_integrand_k_space, frohlich_interaction_energy_integrand_k_space, holstein_interaction_energy_k_space, frohlich_interaction_energy_k_space, holstein_energy_k_space, frohlich_energy_k_space, holstein_structure_factor_k_space, frohlich_structure_factor_k_space, holstein_memory_function_k_space, holstein_mobility_k_space, frohlich_memory_function_k_space, frohlich_mobility_k_space +export spherical_k_integral, holstein_coupling, frohlich_coupling,holstein_interaction_energy_k_space, frohlich_interaction_energy_k_space, holstein_energy_k_space, frohlich_energy_k_space, holstein_structure_factor_k_space, frohlich_structure_factor_k_space, holstein_memory_function_k_space, holstein_mobility_k_space, frohlich_memory_function_k_space, frohlich_mobility_k_space ##### load in library routines... ##### # stdlib @@ -68,50 +68,23 @@ using Unitful using Unitful: @unit, Dimension, Dimensions, NoDims, NoUnits, Units, dimension, uconvert, ustrip # Unitful functions -export puconvert, punit, pustrip, m0_pu, e_pu, ħ_pu, k_pu, ω0_pu, a0_pu, E0_pu, β0_pu, T0_pu, μ0_pu, t0_pu, addunits!, addholsteinunits! - -include("FeynmanTheory.jl") # Frohlich actions + variational functions. -include("HolsteinTheory.jl") # Holstein actions + variational functions. -include("KSpaceTheory.jl") # K-space dependent integral code. -include("HellwarthTheory.jl") # multimode -> equivalent mode. -include("MemoryFunction.jl") # Memory function X calculation. -include("EffectiveMass.jl") # Effective mass (Feynman's ansatz). -include("ResponseFunctions.jl") # Linear reponse functions for polaron. +export puconvert, punit, pustrip, m0_pu, e_pu, ħ_pu, k_pu, ω0_pu, a0_pu, E0_pu, β0_pu, T0_pu, μ0_pu, t0_pu, addunits! + +include("PolaronFunctions.jl") # General use functions for polaron model. +include("LegacyFunctions.jl") # Outdated functions that still have some use. +include("TrialPolaron.jl") # Feynman's Trial Spring-Mass polaron model. +include("FrohlichTheory.jl") # Frohlich polaron theory. +include("HolsteinTheory.jl") # Holstein poalron theory. +include("HellwarthTheory.jl") # Multimode -> equivalent effective mode. include("Material.jl") # Material type and constructors. -include("Polaron.jl") # Polaron type and constructors. -include("PolaronUnits.jl") # Implements internal Feynman 'polaron units', and SI conversions +include("FrohlichPolaron.jl") # Polaron type and constructors. include("HolsteinPolaron.jl") # Holstein type and constructors. +include("PolaronUnits.jl") # Implements internal Feynman 'polaron units', and SI conversions # Register newly defined units with Unitful Unitful.register(PolaronMobility) __init__() = Unitful.register(PolaronMobility) -# QOL function for removing singleton dimensions. -reduce_array(a) = length(a) == 1 ? only(a) : Array(dropdims(a, dims=tuple(findall(size(a) .== 1)...))) - -# Regularized lower gamma function. -function P(n, z) - if n == 1 - erf(sqrt(z)) - elseif n == 2 - 1 - exp(-z) - elseif n > 2 - P_plus_one(n, z) - end -end - -function P_plus_one(n, z) - if n<3 - return P(n, z) - else - return P(n-2, z) - exp(-z) * z^(n/2-1) / gamma(n/2) - end -end - -function ball_surface(n) - 2 * π^(n/2) / gamma(n/2) -end - # Utility functions export reduce_array, P, P_plus_one, ball_surface diff --git a/src/PolaronUnits.jl b/src/PolaronUnits.jl index 6645d75..9aeda1b 100644 --- a/src/PolaronUnits.jl +++ b/src/PolaronUnits.jl @@ -62,7 +62,7 @@ See also: `Unitful.k`, `Unitful.J`, `Unitful.K`. @unit k_pu "k" BoltzmannConstant Unitful.k false """ - PolaronUnits.ω_pu + PolaronUnits.ω0_pu A unit equal to the reduced Planck constant ħ = h / 2π ≈ 1.054,571,8176 × 10^-34 J × s. Printed as "ħ". `Unitful.ħ` is a quantity (with units `J × s`) whereas `PolaronUnits.ħ_pu` is a unit equal to @@ -72,6 +72,18 @@ See also: `Unitful.ħ`, `Unitful.J`, `Unitful.s`. """ @unit ω0_pu "ω₀" PolaronAngularFrequency 1Unitful.THz2π false +""" + PolaronUnits.J0_pu +Transfer integral energy unit +`Unitful.ħ` is a quantity (with units `J × s`) whereas `PolaronUnits.ħ_pu` is a unit equal to +`Unitful.ħ`. +Dimension: 𝐋^2 𝐌 𝐓^-1. +See also: `Unitful.ħ`, `Unitful.J`, `Unitful.s`. +""" +@unit J0_pu "J₀" PolaronTransferIntegral 1Unitful.meV false +@unit ωh0_pu "ωₕ₀" PolaronAdiabticity 1J0_pu/1ħ_pu false +@unit mh0_pu "mₕ₀" PolaronHolsteinMass 1ħ_pu^2/1J0_pu/1Unitful.Å^2 false + # Polaron radius is derived from the base polaron units """ PolaronUnits.a0_pu @@ -84,6 +96,7 @@ Dimension: 𝐋. See also: `Unitful.ε0`, `Unitful.ħ`, `Unitful.me`, `Unitful.q`, `Unitful.m`. """ @unit a0_pu "a₀" PolaronRadius √(1ħ_pu/1m0_pu/1ω0_pu) false +@unit ah0_pu "aₕ₀" PolaronHolsteinRadius √(1ħ_pu/1mh0_pu/1ωh0_pu) false # Polaron energy is derived from the base polaron units """ @@ -97,6 +110,7 @@ Dimension: 𝐋^2 𝐌 𝐓^-2. See also: `Unitful.me`, `Unitful.q`, `Unitful.ε0`, `Unitful.ħ`, `Unitful.J`, `Unitful.eV`, [`UnitfulAtomic.Ry`](@ref). """ @unit E0_pu "E₀" PolaronEnergy 1ħ_pu*1ω0_pu false +@unit Eh0_pu "Eₕ₀" PolaronHolsteinEnergy 1J0_pu false # Polaron thermodynamic temperature is derived from the base polaron units """ @@ -111,8 +125,10 @@ See also: `Unitful.me`, `Unitful.q`, `Unitful.ε0`, `Unitful.ħ`, `Unitful.J`, ` """ @unit T0_pu "T₀" PolaronTemperature 1Unitful.K false @unit β0_pu "β₀" PolaronBeta 1ħ_pu/1k_pu false +@unit βh0_pu "βₕ₀" PolaronHolsteinBeta 1J0_pu/1k_pu false @unit μ0_pu "μ₀" PolaronMobility 1e_pu/1m0_pu/1ω0_pu false +@unit μh0_pu "μₕ₀" PolaronHolsteinMobility 1e_pu/1mh0_pu/1ωh0_pu false @unit t0_pu "t₀" PolaronTime 1Unitful.ns false """ @@ -173,7 +189,7 @@ pustrip(x) = ustrip(puconvert(x)) """ addunits!(polaron::Polaron) """ -function addunits!(polaron::Polaron) +function addunits!(polaron::FrohlichPolaron) polaron.ω = polaron.ω .* ω0_pu polaron.Fs = polaron.Fs .* E0_pu polaron.Fl = polaron.Fl .* E0_pu @@ -215,28 +231,30 @@ function addunits!(polaron::Polaron) polaron.σ = polaron.σ .* punit(u"S") end -function addholsteinunits!(polaron::Holstein) - polaron.ω = polaron.ω .* ω0_pu - polaron.F0 = polaron.F0 .* E0_pu - polaron.K0 = polaron.K0 .* E0_pu - polaron.P0 = polaron.P0 .* E0_pu - polaron.κ0 = polaron.κ0 .* m0_pu * ω0_pu^2 - polaron.M0 = polaron.M0 .* m0_pu - polaron.M0a = polaron.M0a .* m0_pu - polaron.M0r = polaron.M0r .* m0_pu - polaron.R0 = polaron.R0 .* a0_pu +function addunits!(polaron::HolsteinPolaron) + polaron.ω = polaron.ω .* ωh0_pu + polaron.F0 = polaron.F0 .* Eh0_pu + polaron.A0 = polaron.A0 .* Eh0_pu + polaron.B0 = polaron.B0 .* Eh0_pu + polaron.C0 = polaron.C0 .* Eh0_pu + polaron.κ0 = polaron.κ0 .* mh0_pu * ωh0_pu^2 + polaron.M0 = polaron.M0 .* mh0_pu + polaron.M0a = polaron.M0a .* mh0_pu + polaron.M0r = polaron.M0r .* mh0_pu + polaron.R0 = polaron.R0 .* ah0_pu polaron.T = polaron.T .* T0_pu - polaron.β = polaron.β .* β0_pu - polaron.F = polaron.F .* E0_pu - polaron.K = polaron.K .* E0_pu - polaron.P = polaron.P .* E0_pu - polaron.κ = polaron.κ .* m0_pu * ω0_pu^2 - polaron.M = polaron.M .* m0_pu - polaron.Ma = polaron.Ma .* m0_pu - polaron.Mr = polaron.Mr .* m0_pu - polaron.R = polaron.R .* a0_pu - polaron.μ = polaron.μ .* μ0_pu - polaron.χ = polaron.χ .* ω0_pu + polaron.β = polaron.β .* βh0_pu + polaron.F = polaron.F .* Eh0_pu + polaron.A = polaron.A .* Eh0_pu + polaron.B = polaron.B .* Eh0_pu + polaron.C = polaron.C .* Eh0_pu + polaron.κ = polaron.κ .* mh0_pu * ωh0_pu^2 + polaron.M = polaron.M .* mh0_pu + polaron.Ma = polaron.Ma .* mh0_pu + polaron.Mr = polaron.Mr .* mh0_pu + polaron.R = polaron.R .* ah0_pu + polaron.μ = polaron.μ .* μh0_pu + polaron.χ = polaron.χ .* ωh0_pu polaron.z = polaron.z .* punit(u"Ω") polaron.σ = polaron.σ .* punit(u"S") end diff --git a/src/ResponseFunctions.jl b/src/ResponseFunctions.jl deleted file mode 100644 index 2ad5493..0000000 --- a/src/ResponseFunctions.jl +++ /dev/null @@ -1,281 +0,0 @@ -# ReponseFunctions.jl -# - Polaron optical absorption, complex impedence and complex conductivity -# References: -# [1] R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M. Platzman, Mobility of slow electrons in a polar crystal, PhysicalReview127, 1004 (1962) -# [2] Devreese, J. De Sitter, and M. Goovaerts, “Optical absorption of polarons in thefeynman-hellwarth-iddings-platzman approximation,”Phys. Rev. B, vol. 5, pp. 2367–2381, Mar 1972. -# [3] F. Peeters and J. Devreese, “Theory of polaron mobility,” inSolid State Physics,pp. 81–133, Elsevier, 1984. - -""" ----------------------------------------------------------------------- -Polaron absorption coefficient Γ(Ω). ----------------------------------------------------------------------- -""" - -""" - optical_absorption(Ω, β, α, v, w; rtol = 1e-3) - -Calculate the absorption coefficient Γ(Ω) for the polaron at at finite temperatures (Eqn. (11a) in DSG 1972) for a given frequency Ω. - -# Arguments -- `Ω::Float64`: is the frequency (2π THz) of applied electric field. -- `β::Float64`: is the reduced thermodynamic betas. -- `α::Float64`: is the Frohlich alpha coupling parameter. -- `v::Float64`: is the 'v' variational parameter. -- `w::Float64`: is the 'w' variational parameter. -- `rtol`: relative tolerance for the accuracy of any quadrature integrations. - -See DSG 1972: https://doi.org/10.1103/PhysRevB.5.2367. -""" -function optical_absorption(Ω, β, α, v, w) - real(complex_conductivity(Ω, β, α, v, w)) -end - -""" ----------------------------------------------------------------------- -The Complex Impedence and Conductivity of the Polaron. ----------------------------------------------------------------------- -""" - -""" - frohlich_complex_impedence(Ω, β, α, v, w; rtol = 1e-3, T = nothing, verbose = false) - -Calculate the complex impedence Z(Ω) of the polaron at finite temperatures for a given frequency Ω (Eqn. (41) in FHIP 1962). - -# Arguments -- `Ω::Float64`: is the frequency (2π THz) of applied electric field. -- `β::Float64`: is the reduced thermodynamic betas. -- `α::Float64`: is the Frohlich alpha coupling parameter. -- `v::Float64`: is the 'v' variational parameter. -- `w::Float64`: is the 'w' variational parameter. -- `rtol`: relative tolerance for the accuracy of any quadrature integrations. -- `T`: is a token used by `make_polaron()` to keep track of the temperature for printing during a calculation. Do not alter. -- `verbose`: is used by `make_polaron()` to specify whether or not to print. Ignore. - -See FHIP 1962: https://doi.org/10.1103/PhysRev.127.1004. - -""" -function frohlich_complex_impedence(Ω, v, w, α, ωβ...; dims = 3) - -im * (Ω + frohlich_memory_function(Ω, v, w, α, ωβ...; dims = dims)) -end - -""" - frohlich_complex_conductivity(Ω, β, α, v, w; rtol = 1e-3) - -Calculate the complex conductivity σ(Ω) of the polaron at finite temperatures for a given frequency Ω (equal to 1 / Z(Ω) with Z the complex impedence). - -# Arguments -- `Ω::Float64`: is the frequency (2π THz) of applied electric field. -- `β::Float64`: is the reduced thermodynamic betas. -- `α::Float64`: is the Frohlich alpha coupling parameter. -- `v::Float64`: is the 'v' variational parameter. -- `w::Float64`: is the 'w' variational parameter. -- `rtol`: relative tolerance for the accuracy of any quadrature integrations. - -See also [`polaron_complex_impedence`](@ref) -""" -function frohlich_complex_conductivity(Ω, v, w, α, ωβ...; dims = 3) - return 1 / frohlich_complex_impedence(Ω, v, w, α, ωβ...; dims = dims) -end - -""" - inverse_polaron_mobility(v, w, α, ω, β) - -Calculate the inverse of the dc mobility μ of the polaron at finite temperatues (Eqn. (11.5) in F. Peeters and J. Devreese 1984) for a given frequency Ω. - -# Arguments -- `v::Float64`: is the 'v' variational parameter. -- `w::Float64`: is the 'w' variational parameter. -- `α::Float64`: is the Frohlich alpha coupling parameter. -- `ω::Float64`: is the angular phonon frequency. -- `β::Float64`: is the reduced thermodynamic beta. -- `verbose`: is used by `make_polaron()` to specify whether or not to print. Ignore. - -See F. Peeters and J. Devreese 1984: https://doi.org/10.1016/S0081-1947(08)60312-4. - -See also [`polaron_mobility`](@ref), [`polaron_complex_conductivity`](@ref) -""" -function inverse_frohlich_mobility(v, w, α, ω, β; dims = 3) - if β == Inf - return zero(β) - end - structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω, β; dims = dims) - return abs(imag(general_memory_function(structure_factor))) -end - -""" - inverse_polaron_mobility(v, w, α::Vector, ω::Vector, β::Vector) - -inverse of the polaron mobility, but for multiple phonon modes. -""" -inverse_frohlich_mobility(v, w, α::Vector, ω::Vector, β) = sum(inverse_frohlich_mobility.(v, w, α, ω, β)) - -""" - polaron_mobility(v, w, α, ω, β) - -The polaron mobility. - -See also [`inverse_polaron_mobility`](@ref) -""" -frohlich_mobility(v, w, α, ω, β; dims = 3) = 1 ./ inverse_frohlich_mobility(v, w, α, ω, β; dims = dims) - -""" - inverse_FHIP_mobility_lowT(v, w, α, ω, β) - -FHIP low-temperature mobility, final result of Feynman1962. -[1.60] in Devreese2016 page 36; 6th Edition of Frohlich polaron notes (ArXiv). - -See also [`FHIP_mobility_lowT`](@ref) -""" -function inverse_FHIP_mobility_lowT(v, w, α, ω, β) - if β == Inf - return 0 - else - μ = (w / v)^3 * 3 / (4 * ω^2 * α * β) * exp(ω * β) * exp((v^2 - w^2) / (w^2 * v)) - return 1 / μ - end -end - -""" - inverse_FHIP_mobility_lowT(v, w, α::Vector, ω::Vector, β::Vector) - -Inverse FHIP mobility for multiple phonon modes. - -See also [`FHIP_mobility_lowT`](@ref) -""" -inverse_FHIP_mobility_lowT(v, w, α::Vector, ω::Vector, β) = sum(inverse_FHIP_mobility_lowT.(v, w, α, ω, β)) - -""" - FHIP_mobility_lowT(v, w, α, ω, β) - -FHIP mobility in the low-temperature approximation. - -See also [`inverse_FHIP_mobility_lowT`](@ref) -""" -FHIP_mobility_lowT(v, w, α, ω, β) = 1 ./ inverse_FHIP_mobility_lowT(v, w, α, ω, β) - -""" - inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) - -Kadanoff low-temperaure mobility, constructed around Boltzmann equation. -Adds factor of 3 / (2 * β) c.f. FHIP, correcting phonon emission behaviour. - -See also [`Kadanoff_mobility_lowT`](@ref) -""" -function inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) - - if β == Inf - return 0, 0, 0 - else - - # [1.61] in Devreese2016 - Kadanoff's Boltzmann eqn derived mobility. - μ_Devreese2016 = (w / v)^3 / (2 * ω * α) * exp(ω * β) * exp((v^2 - w^2) / (w^2 * v)) - - # From 1963 Kadanoff (particularly on right-hand-side of page 1367), Eqn. (9), - # we define equilibrium number of phonons (just from temperature T and phonon ω): - # N̄ = (exp(β) - 1)^-1. - # But! We find that: - N̄ = exp(-β * ω) - - # is only way to get Kadanoff1963 to be self-consistent with - # FHIP, and later statements (Devreese) of the Kadanoff mobility. - # It suggests that Kadanoff used the wrong identy for Nbar in Eqn. (23b) for - # the Γ₀ function, and should have used a version with the -1 to - # account for Bose / phonon statistics! - - # Between Eqns. (23) and (24) in Kadanoff1963, for small momenta skip intergration - # and sing the fictitious mass M: - M = (v^2 - w^2) / w^2 - - # we get for Γ₀: - Γ₀ = 2 * α * N̄ * (M + 1)^(1 / 2) * exp(-M / v) * ω - - # NB: Kadanoff1963 uses ħ = ω = mb = 1 units. - # Factor of omega to get it as a rate relative to phonon frequency - # Factor of omega*hbar to get it as a rate per energy window. - # Hence, for the mobility from Eqn. (25) in Kadanoff1963 (SI effective mass q / mb): - μ_Kadanoff1963 = 1 / ((M + 1) * Γ₀) - - # # Energy loss is simply energy * rate: - # energy_loss = ω * Γ₀ / 2π - - # # Lifetime is: - # τ = 1 / Γ₀ - - return 1 / μ_Devreese2016, 1 / μ_Kadanoff1963, Γ₀ - end -end - -""" - inverse_Kadanoff_mobility_lowT(v, w, α::Vector, ω::Vector, β::Vector) - -Inverse Kadanoff mobility for multiple phonon modes. - -See also [`Kadanoff_mobility_lowT`](@ref) -""" -inverse_Kadanoff_mobility_lowT(v, w, α::Vector, ω::Vector, β) = map(+, map(x -> 1 ./ x, inverse_Kadanoff_mobility_lowT.(v, w, α, ω, β))...) - -""" - Kadanoff_mobility_lowT(v, w, α, ω, β) - -Kadanoff mobility in the low-temperature approximation. - -See also [`inverse_Kadanoff_mobility_lowT`](@ref) -""" -Kadanoff_mobility_lowT(v, w, α, ω, β) = 1 ./ inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) - -""" - inverse_Hellwarth_mobility(v, w, α, ω, β) - -Calculates the DC mobility using Hellwarth et al. 1999 Eqn. (2). -Directly performs contour integration in Feynman1962, for finite temperature DC mobility. -Eqns. (2) and (1) are going back to the general (pre low-T limit) formulas in Feynman1962. -To evaluate these, you need to do the explicit contour integration to get the polaron self-energy. - -See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. - -See also [`Hellwarth_mobility`](@ref) -""" -function inverse_Hellwarth_mobility(v, w, α, ω, β) - - if β == Inf - return 0, 0 - else - - R = (v^2 - w^2) / (w^2 * v) # inline, page 300 just after Eqn (2) - b = R * ω * β / sinh(ω * β * v / 2) # Feynman1962 version; page 1010, Eqn (47b) - a = sqrt((ω * β / 2)^2 + R * ω * β * coth(ω * β * v / 2)) - k(u, a, b, v) = (u^2 + a^2 - b * cos(v * u) + eps(Float64))^(-3 / 2) * cos(u) # integrand in (2) - K = quadgk(u -> k(u, a, b, v), 0, 1e5)[1] # numerical quadrature integration of (2) - - # Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997 - RHS = α / (3 * sqrt(π)) * (ω * β)^(5 / 2) / sinh(ω * β / 2) * (v^3 / w^3) * K - μ = RHS - - # Hellwarth1999/Biaggio1997, b=0 version... 'Setting b=0 makes less than 0.1% error' - # So let's test this: - K_0 = quadgk(u -> k(u, a, 0, v), 0, 1e5)[1] # Inserted b=0 into k(u, a, b, v). - RHS_0 = α / (3 * sqrt(π)) * (ω * β)^(5 / 2) / sinh(ω * β / 2) * (v^3 / w^3) * K_0 - μ_0 = RHS_0 - - return μ, μ_0 - end -end - -""" - inverse_Hellwarth_mobility(v, w, α::Vector, ω::Vector, β::Vector) - -Inverse Hellwarth mobility for multiple phonon modes. - -See also [`Hellwarth_mobility`](@ref) -""" -inverse_Hellwarth_mobility(v, w, α::Vector, ω::Vector, β) = map(+, map(x -> 1 ./ x, inverse_Hellwarth_mobility.(v, w, α, ω, β))...) - -""" - Hellwarth_mobility(v, w, α, ω, β) - -The Hellwarth polaron mobility. - -See also [`inverse_Hellwarth_mobility`](@ref) -""" -Hellwarth_mobility(v, w, α, ω, β) = 1 ./ inverse_Hellwarth_mobility(v, w, α, ω, β) - diff --git a/src/TrialPolaron.jl b/src/TrialPolaron.jl new file mode 100644 index 0000000..341f694 --- /dev/null +++ b/src/TrialPolaron.jl @@ -0,0 +1,310 @@ +# TrialPolaron.jl +# This file contains all functions associated with calculating properties of the trial polaron system. + +""" + polaron_propagator(τ, v, w, β) + +Calculate the imaginary time polaron Green's function with temperature dependence. + +# Arguments +- `τ`: A scalar representing the imaginary time. +- `v`: A scalar representing a variational parameter. +- `w`: A scalar representing a variational parameter. +- `β`: A scalar representing the inverse temperature. + +# Returns +The value of the polaron propagator. + +# Example +```julia +τ = 0.5 +v = 0.2 +w = 0.1 +β = 1.0 +result = polaron_propagator(τ, v, w, β) +println(result) +``` +This example calculates the polaron propagator for given values of τ, v, w, and β. The result is then printed. +""" +function polaron_propagator(τ, v, w, ω, β) + (v^2 - w^2) / v^3 * (1 - exp(-v * τ)) * (1 - exp(-v * (β * ω - τ))) / (1 - exp(-v * β * ω)) + w^2 / v^2 * τ * (1 - τ / β / ω) + eps(Float64) +end + +""" + polaron_propagator(τ, v, w) + +Calculate the value of the polaron propagator based on the given inputs. + +# Arguments +- `τ::Number`: A scalar representing the imaginary time. +- `v::Number`: A scalar representing a variational parameter. +- `w::Number`: A scalar representing a variational parameter. + +# Returns +The value of the polaron propagator. + +# Example +```julia +τ = 0.5 +v = 0.2 +w = 0.1 +result = polaron_propagator(τ, v, w) +println(result) +``` +This example calculates the polaron propagator for the given values of τ, v, and w. The result is then printed. +""" +function polaron_propagator(τ, v, w, ω) + w^2 * τ / v^2 + (v^2 - w^2) / v^3 * (1 - exp(-v * τ)) + eps(Float64) +end + +""" + D(τ, v, w, β) + +Calculates the recoil function (a generalisation of D(u) in Eqn. (35c) in FHIP 1962) that approximates the exact influence (recoil effects) of the phonon bath on the electron with the influence of the fictitious masses attached by springs to the electron. It appears in the exponent of the intermediate scattering function. + +# Arguments +- `τ::Float64`: is the imaginary time variable. +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. + +See FHIP 1962: https://doi.org/10.1103/PhysRev.127.1004. +""" +function D(τ, v::Vector, w::Vector, ω, β) + return τ * (1 - τ / ω / β) + sum((h_i(i, v, w) / v[i]^2) * ((1 + exp(-v[i] * ω * β) - exp(-v[i] * τ) - exp(v[i] * (τ - ω * β))) / (v[i] * (1 - exp(-v[i] * ω * β))) - τ * (1 - τ / ω / β)) for i in eachindex(v)) +end + +""" + D(τ, v, w) + +Calculates the recoil function at zero-temperature. + +# Arguments +- `τ::Float64`: is the imaginary time variable. +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +""" +function D(τ, v::Vector, w::Vector) + return τ + sum((h_i(i, v, w) / v[i]^2) * ((1 - exp(-v[i] * τ)) / v[i] - τ) for i in eachindex(v)) +end + +# Hellwarth et al. 1999 PRB - Part IV; T-dep of the Feynman variation parameter + +# In Julia we have 'Multiple dispatch', so let's just construct the free +# energies (temperature-dependent) with the same name as before, but with the thermodynamic beta where required. + +# Define Osaka's free-energies (Hellwarth1999 version) as Julia functions. + +""" + A(v, w, ω, β) + +Hellwarth's A expression from Eqn. (62b) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression. + +See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. +""" +A(v, w, ω, β) = β == Inf ? A(v, w, ω) : 3 / β * (log(v / w) - 1 / 2 * log(2π * β * ω +) - (v - w) * ω * β / 2 - log(1 - exp(-v * β * ω)) + log(1 - exp(-w * β * ω))) + +A(v, w, ω::Vector, β) = sum(A.(v, w, ω, β)) / length(ω) + +A(v, w, ω) = -3 * (v - w) / 2 * ω + +A(v, w, ω::Vector) = sum(A.(v, w, ω)) / length(ω) + +""" + C(v, w, ω, β) + +Hellwarth's C expression from Eqn. (62e) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression. + +See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. +""" +C(v, w, ω, β) = 3 / 4 * (v^2 - w^2) * ω / v * (coth(v * β * ω / 2) - 2 / (v * β * ω)) + +C(v, w, ω::Vector, β) = sum(C.(v, w, ω, β)) / length(ω) + +C(v, w, ω) = (3 / (4 * v)) * (v^2 - w^2) * ω + +C(v, w, ω::Vector) = sum(C.(v, w, ω)) / length(ω) + +# Extending the Feynman theory to multiple phonon branches + +# Extending the Feynman theory to multiple variational parameters +# Multiple Parameter Polaron Free Energy +# Calculate the polaron free energy, generalised from Osaka's expression to the case where multiple phonon modes are present in the material. + +""" + κ_i(i, v, w) + +Calculates the spring-constant coupling the electron to the 'ith' fictitious mass that approximates the exact electron-phonon interaction with a harmonic coupling to a massive fictitious particle. + +Required for calculating the polaron free energy. + +Note: Not to be confused with the number of physical phonon branches; many phonon branches could be approximated with one or two etc. fictitious masses for example. The number of fictitious mass does not necessarily need to match the number of phonon branches. + +# Arguments +- `i::Integer`: enumerates the current fictitious mass. +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +""" +function κ_i(i, v::Vector, w::Vector) + κ = v[i]^2 - w[i]^2 + κ *= prod(j != i ? (v[j]^2 - w[i]^2) / (w[j]^2 - w[i]^2) : 1 for j in eachindex(v)) + return κ +end + +""" + h_i(i, v, w) + +Calculates the normal-mode (the eigenmodes) frequency of the coupling between the electron and the `ith' fictitious mass that approximates the exact electron-phonon interaction with a harmonic coupling to a massive fictitious particle. + +Required for calculating the polaron free energy. + +Note: Not to be confused with the number of physical phonon branches; many phonon branches could be approximated with one or two etc. fictitious masses for example. The number of fictitious mass does not necessarily need to match the number of phonon branches. + +# Arguments +- `i::Integer`: enumerates the current fictitious mass. +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +""" +function h_i(i, v::Vector, w::Vector) + h = v[i]^2 - w[i]^2 + h *= prod(j != i ? (w[j]^2 - v[i]^2) / (v[j]^2 - v[i]^2) : 1 for j in eachindex(v)) + return h +end + +""" + C_ij(i, j, v, w) + +Calculates the element to the coupling matrix C_ij (a generalisation of Feynman's `C` coupling variational parameter in Feynman 1955) between the electron and the `ith' and `jth' fictitious masses that approximates the exact electron-phonon interaction with a harmonic coupling to a massive fictitious particle. + +Required for calculating the polaron free energy. + +Note: Not to be confused with the number of physical phonon branches; many phonon branches could be approximated with one or two etc. fictitious masses for example. The number of fictitious mass does not necessarily need to match the number of phonon branches. + +# Arguments +- `i::Integer, j::Integer`: enumerate the current fictitious masses under focus (also the index of the element in the coupling matrix C) +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. + +See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. +""" +function C_ij(i, j, v::Vector, w::Vector) + C = w[i] * κ_i(i, v, w) * h_i(j, v, w) / (4 * (v[j]^2 - w[i]^2)) + return C +end + +""" + C(v, w, β) + +Generalisation of the C function from Eqn. (62e) in Hellwarth et al. 1999. This is the expected value of the trial action taken w.r.t trial action. + +Required for calculating the polaron free energy. + +# Arguments +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. + + +See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. +""" +function C(v::Vector, w::Vector, ω, β) + # Sum over the contributions from each fictitious mass. + s = sum(C_ij(i, j, v, w) / (v[j] * w[i]) * (coth(ω * β * v[j] / 2) - 2 / (ω * β * v[j])) for i in eachindex(v), j in eachindex(w)) + + # Divide by the number of phonon modes to give an average contribution per phonon mode. + return 3 * s * ω +end + +C(v::Vector, w::Vector, ω::Vector, β) = sum(C(v, w, ω[i], β) for i in eachindex(ω)) / length(ω) + +""" + C(v, w) + +Calculates `C(v, w, β)` but at zero-temperature, `β = Inf`. + +# Arguments +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +""" +function C(v::Vector, w::Vector, ω) + # Sum over the contributions from each fictitious mass. + s = sum(C_ij(i, j, v, w) / (v[j] * w[i]) for i in eachindex(v), j in eachindex(w)) + # Divide by the number of phonon modes to give an average contribution per phonon mode. + return 3 * s * ω +end + +C(v::Vector, w::Vector, ω::Vector) = sum(C(v, w, ω[i]) for i in eachindex(ω)) / length(ω) + +""" + A(v, w, β) + +Generalisation of the A function from Eqn. (62b) in Hellwarth et al. 1999. This is the Helmholtz free energy of the trial model. + +Required for calculating the polaron free energy. + +# Arguments +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +- `β::Union{Float64, Vector{Float64}}`: is the reduced thermodynamic temperature ħωⱼ/(kT) associated with the 'jth' phonon mode. + +See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. +""" +function A(v::Vector, w::Vector, ω, β) + # Sum over the contributions from each fictitious mass. + s = -log(2π * ω * β) / 2 + sum(v[i] == w[i] ? 0 : + log(v[i]) - log(w[i]) - ω * β / 2 * (v[i] - w[i]) - log(1 - exp(-v[i] * ω * β)) + log(1 - exp(-w[i] * ω * β)) + for i in eachindex(v)) + # Divide by the number of phonon modes to give an average contribution per phonon mode. + 3 / β * s * ω +end + +A(v::Vector, w::Vector, ω::Vector, β) = sum(A(v, w, ω[i], β) for i in eachindex(ω)) / length(ω) + +""" + A(v, w, n) + +Calculates `A(v, w, β)` but at zero-temperature, `β = Inf`. + +# Arguments +- `v::Vector{Float64}`: is a vector of the v variational parameters. +- `w::Vector{Float64}`: is a vector of the w variational parameters. +""" +function A(v::Vector, w::Vector, ω) + s = sum(v .- w) + return -3 * s / 2 * ω +end + +A(v::Vector, w::Vector, ω::Vector) = sum(A(v, w, ω[i]) for i in eachindex(ω)) / length(ω) + +""" + trial_energy(v, w, ω, β) + +Calculate the free electron energy at finite temperature. + +# Arguments +- `v`: a scalar value representing a variational paramater. +- `w`: a scalar value representing a variational paramater. +- `ω`: a scalar value representing the phonon frequency. +- `β`: a scalar value representing the inverse temperature. + +# Returns +A scalar value representing the calculated free electron energy at finite temperature. + +# Example +```julia +v = 0.5 +w = 1.0 +ω = 2.0 +β = 0.2 +result = electron_energy(v, w, ω, β; dims = 3) +println(result) +``` +Expected Output: +A scalar value representing the calculated free electron energy at finite temperature. +""" +function trial_energy(v, w, ωβ...; dims = 3) + Ar = A(v, w, ωβ...) * dims / 3 + Cr = C(v, w, ωβ...) * dims / 3 + return Ar, Cr +end \ No newline at end of file