diff --git a/src/FrohlichPolaron.jl b/src/FrohlichPolaron.jl index a3f1903..9e38dc7 100644 --- a/src/FrohlichPolaron.jl +++ b/src/FrohlichPolaron.jl @@ -9,6 +9,7 @@ Type for storing the polaron information, `x...`. mutable struct FrohlichPolaron α # Fröhlich coupling. αeff # Effective Fröhlich coupling summed for multiple modes. + mb # Band mass T # Temperature. ω # Phonon mode frequencies. ωeff # Effective phonon mode frequency. @@ -70,7 +71,7 @@ Outer constructor for the Polaron type. This function evaluates model data for t julia> polaron(6, 300, 3, 1.0, 3.6, 2.8) ``` """ -function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.1, w_guesses=2.9, dims=3, kspace=false, verbose=false) +function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesses=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. @@ -83,13 +84,6 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v num_d = length(dims) ω = reduce_array(ω) - - Trange = pustrip.(Trange .* 1u"K") - Ωrange = pustrip.(Ωrange .* 1u"THz2π") - ω = pustrip.(ω .* 1u"THz2π") - ωeff = pustrip.(ωeff .* 1u"THz2π") - - a0 = pustrip.(sqrt(Unitful.ħ / (2 * Unitful.me * mb * ω * u"THz2π"))) # For multiple variational modes, ensure that the number of v and w parameters is the same. @assert length(v_guesses) == length(w_guesses) "v and w guesses must be the same length." @@ -99,6 +93,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v p = Dict( "α" => αrange, # Alphas. "αeff" => sum(αrange, dims=2), # Alphas sums. + "mb" => mb, # Band mass. "T" => Trange, # Temperatures. "ω" => ω, # Phonon frequencies. "ωeff" => ωeff, # Effective Phonon frequency. @@ -111,8 +106,8 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v "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. + "Fs" => Vector{Float64}(undef, num_α), # Small alpha (α→0) approximate energy. + "Fl" => Vector{Float64}(undef, 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. @@ -185,7 +180,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Small alpha (α → 0) approximate energy. - F_small = (-αeff - αeff^2 / 81) * ω + F_small = (-αeff - αeff^2 / 81) p["Fs"][j] = F_small # Print small alpha energy. @@ -194,7 +189,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Large alpha (α → ∞) approximate energy. - F_large = (-αeff^2 / 3π - 3 * log(2) - 3 / 4) * ω + F_large = (-αeff^2 / 3π - 3 * log(2) - 3 / 4) p["Fl"][j] = F_large # Print large alpha energy. @@ -203,7 +198,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Small coupling (α → 0) polaron mass approximation. Eqn. (46) in Feynman1955. - M_small = (αeff / 6 + 0.025 * αeff^2) * mb + M_small = (αeff / 6 + 0.025 * αeff^2) p["Ms"][j] = M_small # Print small alpha fictitious mass. @@ -212,7 +207,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Large coupling (α → ∞) polaron mass approximation. Eqn. (47) In Feynman1955. - M_large = 16 * αeff^4 / (81 * π^4) * mb + M_large = 16 * αeff^4 / (81 * π^4) p["Ml"][j] = M_large # Print large alpha fictitious mass. @@ -221,7 +216,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Small coupling (α → 0) polaron radii approximiation. Eqn. (2.5a) in Schultz1959. - R_small = sqrt(3 / (4 / 9 * αeff)) * a0 + R_small = sqrt(3 / (4 / 9 * αeff)) p["Rs"][j] = R_small # Print small alpha polaron radius. @@ -230,7 +225,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Large coupling (α → ∞) polaron radii approximiation. Eqn. (2.5b) in Schultz1959. - R_large = 3 * √(π / 2) * αeff * a0 + R_large = 3 * √(π / 2) * αeff p["Rl"][j] = R_large # Print large alpha polaron radius. @@ -239,7 +234,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Franck-Condon (FC) frequency in large coupling (α → ∞) limit. RHS of pg. 2371 in Devreese1972. - Ω_FC = 4 / 9π * αeff ^2 * ω + Ω_FC = 4 / 9π * αeff ^2 p["ΩFC"][j] = Ω_FC # Print large alpha Franck-Condon peak frequency. @@ -257,7 +252,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v # w is also the frequency of oscillation of the SHM trial system composed of the bare particle and fictitous mass. # A, B, C are components of the total energy: A is the bare electron energy, B the electron-phonon interaction energy, C is the energy of the harmonic trial system. - athermal_energy(v, w) = !kspace ? frohlich_energy(v, w, α, ω; dims = dims[d]) : frohlich_energy_k_space(v, w, α, ω; dims = dims[d]) + athermal_energy(v, w) = !kspace ? frohlich_energy(v, w, α, ω; dims = dims[d], mb = mb) : frohlich_energy_k_space(v, w, α, ω; dims = dims[d]) v_gs, w_gs, F_gs, A_gs, B_gs, C_gs = vw_variation(athermal_energy, v_guesses, w_guesses) @@ -265,7 +260,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v v_guesses, w_guesses = v_gs, w_gs # Store the athermal data. - p["v0"][d, j, :] .= v_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 @@ -287,7 +282,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal - κ_gs = (v_gs .^ 2 .- w_gs .^ 2) .* mb + κ_gs = (v_gs .^ 2 .- w_gs .^ 2) p["κ0"][d, j, :] .= κ_gs # Print athermal fictitious spring constant. @@ -305,7 +300,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Athermal - M_asymp_gs = v_gs ./ w_gs .* mb + M_asymp_gs = v_gs ./ w_gs p["M0a"][d, j, :] .= M_asymp_gs # Print athermal asymptotic fictitious mass. @@ -314,7 +309,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # 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 .* mb + M_reduced_gs = (v_gs .^2 .- w_gs .^2) ./ v_gs .^ 2 p["M0r"][d, j, :] .= M_reduced_gs # Print athermal reduced mass. @@ -323,8 +318,8 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # 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) .* a0 - p["R0"][d, j, :] .= R_gs + R_gs = sqrt.(3 ./ 2 .* M_reduced_gs .* v_gs) + p["R0"][d, j, :] .= R_gs ./ sqrt(ωeff) # Print athermal polaron radius. if verbose @@ -357,7 +352,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v # 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 ? frohlich_energy(v, w, α, ω, β; dims = dims[d]) : frohlich_energy_k_space(v, w, α, ω, β; dims = dims[d]) + thermal_energy(v, w) = !kspace ? frohlich_energy(v, w, α, ω, β; dims = dims[d], mb = mb) : frohlich_energy_k_space(v, w, α, ω, β; dims = dims[d]) v, w, F, A, B, C = vw_variation(thermal_energy, v_guesses, w_guesses) # Update the guesses to keep them close-ish to the true solutions during loops over temperatures. @@ -382,7 +377,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal - κ = (v .^ 2 .- w .^ 2) .* mb + κ = (v .^ 2 .- w .^ 2) p["κ"][i, d, j, :] .= κ # Print spring constants. @@ -400,7 +395,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Thermal - M_asymp = v ./ w .* mb + M_asymp = v ./ w p["Ma"][i, d, j, :] .= M_asymp # Print asymptotic masses. @@ -409,7 +404,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. - M_reduced = (v .^2 .- w .^2) ./ v .^ 2 .* mb + M_reduced = (v .^2 .- w .^2) ./ v .^ 2 p["Mr"][i, d, j, :] .= M_reduced # Print redcued masses. @@ -418,8 +413,8 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v end # Calculate and store polaron radii. - R = sqrt.(3 .* v ./ (v .^ 2 .- w .^ 2) .^ 2) .* a0 - p["R"][i, d, j, :] .= R + R = sqrt.(3 ./ 2 .* M_reduced .* v) + p["R"][i, d, j, :] .= R ./ sqrt(ωeff) # Print polaron radius. if verbose @@ -433,7 +428,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v 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 + μ = !kspace ? frohlich_mobility(v, w, α, ω, β; dims = dims[d]) : frohlich_mobility_k_space(v, w, α, ω, β; dims = dims[d]) p["μ"][i, d, j] = μ # Print DC mobilities. @@ -443,7 +438,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v # 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[1], w[1], α, ω, β) / mb + μ_FHIP = FHIP_mobility_lowT(v[1] ./ ωeff, w[1] ./ ωeff, α, ω, β) p["μFHIP"][i, j] = μ_FHIP # Print low-temperature FHIP mobility. @@ -458,7 +453,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v # 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[1], w[1], α, ω, β) ./ mb + μ_Kadanoff_Devreese, μ_Kadanoff, rel_time = Kadanoff_mobility_lowT(v[1] ./ ωeff, w[1] ./ ωeff, α, ω, β) p["μD"][i, j] = μ_Kadanoff_Devreese p["μK"][i, j] = μ_Kadanoff p["τ"][i, j] = rel_time @@ -474,7 +469,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v # 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[1], w[1], α, ω, β) ./ mb + μ_Hellwarth, μ_Hellwarth_b0 = Hellwarth_mobility(v[1], w[1], α, ω, β) p["μH"][i, j] = μ_Hellwarth p["μH0"][i, j] = μ_Hellwarth_b0 @@ -490,7 +485,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v v, w, β = v_gs, w_gs, Inf p["β"][i] = β p["v"][i, d, j, :] .= v_gs - p["w"][i, d, j, :] .= w_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 @@ -517,6 +512,8 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v for k in eachindex(Ωrange) # E-field frequencies loop. Ω = Ωrange[k] + if !iszero(T) || !iszero(Ω) + # Print E-field frequency. if verbose println(io, "\e[K-----------------------------------------------------------------------") @@ -536,7 +533,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v # Calculate and store polaron complex impedances. - z = -(im * Ω + im * χ) .* mb + z = -(im * Ω + im * χ) p["z"][k, i, d, j] = z # Print complex impedances. @@ -553,7 +550,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v println(io, "\e[KComplex conductivity | σ = ", σ) end - if iszero(T) && iszero(Ω) + else # If zero frequency. p["χ"][k, i, d, j] = Inf + 0 * im @@ -570,7 +567,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v print(io, "\e[2F") end - if verbose print(io, "\e[7F") end + if verbose && (!iszero(T) || !iszero(Ω)) print(io, "\e[7F") end end if verbose && !iszero(T) print(io, "\e[26F") end # Move up 26 lines and erase. end @@ -585,6 +582,7 @@ function frohlichpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v polaron_data = [ p["α"], # Alphas. p["αeff"], # Alphas sums. + p["mb"], # Band mass. p["T"], # Temperatures. p["ω"], # Phonon frequencies. p["ωeff"], # Hellwarth eff phonon frequency. @@ -639,22 +637,22 @@ end """ Single alpha parameter. polaron() expects alpha parameters to be in a Vector. """ -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) +frohlichpolaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No frequency input. """ -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) +frohlichpolaron(αrange, Trange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(αrange, Trange, 0; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No temperature input => 300 K. """ -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) +frohlichpolaron(αrange; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(αrange, 0, 0; ω=ω, ωeff=ωeff, mb=mb,v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No input => α = 1 """ -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) +frohlichpolaron(; ω=1, ωeff=1, mb=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = frohlichpolaron(1, 0, 0; ω=ω, ωeff=ωeff, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ polaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, verbose=false) @@ -674,6 +672,8 @@ function frohlichpolaron(material::Material, TΩrange...; v_guesses=3.11, w_gues phonon_eff_freq = material.feff mb = material.mb + TΩrange = length(TΩrange) == 1 ? TΩrange .* pustrip(1u"K") : TΩrange .* (pustrip(1u"K"),1) + # Generate polaron data from the arbitrary model constructor. p = frohlichpolaron(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) @@ -773,6 +773,7 @@ function save_frohlich_polaron(polaron::FrohlichPolaron, prefix) JLD.save("$prefix.jld", "alpha", polaron.α, "alpha eff", polaron.αeff, + "band mass", polaron.mb, "temperature", polaron.T, "phonon freq", polaron.ω, "phonon freq eff", polaron.ωeff, @@ -837,6 +838,7 @@ function load_frohlich_polaron(polaron_file_path) polaron = FrohlichPolaron( data["alpha"], data["alpha eff"], + data["band mass"], data["temperature"], data["phonon freq"], data["phonon freq eff"], diff --git a/src/FrohlichTheory.jl b/src/FrohlichTheory.jl index f60e204..ff4a4af 100644 --- a/src/FrohlichTheory.jl +++ b/src/FrohlichTheory.jl @@ -125,8 +125,8 @@ println(result) Expected Output: `6.0` """ -function frohlich_coupling(k, α, ω; mb = 1, dims = 3) - r_p = 1 / sqrt(2) +function frohlich_coupling(k, α, ω; dims = 3) + r_p = sqrt(1 / 2) ω^2 * α * r_p * gamma((dims - 1) / 2) * (2√π / k)^(dims - 1) end @@ -142,7 +142,8 @@ 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, ωβ...)) + propagator(τ) = length(ωβ) == 1 ? polaron_propagator(τ, v, w) * ωβ[1] : polaron_propagator(τ, v, w, ωβ[2]) * ωβ[1] + integrand(τ) = phonon_propagator(τ, ωβ...) / sqrt(propagator(τ)) upper_limit = length(ωβ) == 1 ? Inf : ωβ[2] / 2 integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) return coupling * ball_surface(dims) / (2π)^dims * sqrt(π / 2) * integral @@ -173,10 +174,10 @@ A scalar value representing the Frohlich polaron interaction energy in k-space a """ 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, ωβ...) + propagator(τ) = length(ωβ) == 1 ? polaron_propagator(τ, v, w) * ωβ[1] : polaron_propagator(τ, v, w, ωβ[2]) * ωβ[1] 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) + upper_limit = length(ωβ) == 1 ? Inf : ωβ[2] / 2 + integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) return integral end @@ -196,10 +197,10 @@ This generalises the Osaka 1959 (below Eqn. (22)) and Hellwarth. et al 1999 (Eqn 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 +function frohlich_energy(v, w, α, ωβ...; dims = 3, mb = 1) + A, C = length(ωβ) == 1 ? trial_energy(v, w; dims = dims) : trial_energy(v, w, ωβ[2]; dims = dims) + B = frohlich_interaction_energy(v, w, α, ωβ...; dims = dims) + return -(A + B + C), A, B, C end """ @@ -222,7 +223,7 @@ Calculate the total energy, kinetic energy, and interaction energy of the Frohli - `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, ωβ...) + A, C = length(ωβ) == 1 ? trial_energy(v, w; dims = dims) : trial_energy(v, w, ωβ[2]; dims = dims) B = frohlich_interaction_energy_k_space(v, w, α, ωβ...; limits = limits, dims = dims) return -(A + B + C), A, B, C end @@ -240,10 +241,10 @@ The Complex Impedence and Conductivity of the Polaron. # [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 + coupling = frohlich_coupling(1, α, ωβ[1]; dims = dims) * ωβ[1] + propagator = length(ωβ) == 1 ? polaron_propagator(im * t, v, w) * ωβ[1] / 2 : polaron_propagator(im * t, v, w, ωβ[2]) * ωβ[1] / 2 integral = ball_surface(dims) / (2π)^dims * √π / 4 / propagator^(3/2) - return 2 / dims * ωβ[1] * coupling * integral * phonon_propagator(im * t, ωβ...) + return 2 / dims * coupling * integral * phonon_propagator(im * t, ωβ...) end """ @@ -265,10 +266,10 @@ Calculate the structure factor in k-space for the Frohlich continuum polaron mod 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, ωβ...) + coupling(k) = frohlich_coupling(k, α, ωβ[1]; dims = dims) * k^2 * ωβ[1] + propagator = length(ωβ) == 1 ? polaron_propagator(im * t, v, w) * ωβ[1] : polaron_propagator(im * t, v, w, ωβ[2]) * ωβ[1] integral = spherical_k_integral(coupling, propagator; limits = limits, dims = dims) - return 2 / dims * phonon_propagator(im * t, ω, β) * integral + return 2 / dims * phonon_propagator(im * t, ωβ...) * integral end frohlich_structure_factor_k_space(t, v, w, α::Vector, ω::Vector; limits = [0, Inf], dims = 3) = sum(frohlich_structure_factor_k_space(t, v, w, α[j], ω[j]; limits = limits, dims = dims) for j in eachindex(α)) @@ -277,7 +278,7 @@ frohlich_structure_factor_k_space(t, v, w, α::Vector, ω::Vector, β; limits = function frohlich_memory_function(Ω, v, w, α, ωβ...; dims = 3) structure_factor(t) = frohlich_structure_factor(t, v, w, α, ωβ...; dims = dims) - return polaron_memory_function(Ω, structure_factor; limits = [0, 1e4]) + return polaron_memory_function(Ω, structure_factor) end frohlich_memory_function(Ω, v, w, α::Vector, ω::Vector; dims = 3) = sum(frohlich_memory_function(Ω, v, w, α[j], ω[j]; dims = dims) for j in eachindex(α)) @@ -290,7 +291,7 @@ frohlich_memory_function(Ω, v, w, α::Vector, ω::Vector, β; dims = 3) = sum(f 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. +- `Ω`: a scalar value representing the frequency at which the memorpy 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. @@ -367,11 +368,8 @@ See F. Peeters and J. Devreese 1984: https://doi.org/10.1016/S0081-1947(08)60312 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, 1e4]))) + return abs(imag(polaron_memory_function(structure_factor))) end """ @@ -411,9 +409,6 @@ Calculate the DC mobility in k-space for a Frohlich polaron system at finite tem 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 @@ -430,7 +425,7 @@ 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)) + μ = (w / v)^3 * 3 / (4 * ω^2 * α * β) * exp(ω * β) * exp((v^2 - w^2) * ω / (w^2 * v)) return 1 / μ end end @@ -468,7 +463,7 @@ function inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) 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)) + μ_Devreese2016 = (w / v)^3 / (2 * ω^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 ω): @@ -487,7 +482,7 @@ function inverse_Kadanoff_mobility_lowT(v, w, α, ω, β) M = (v^2 - w^2) / w^2 # we get for Γ₀: - Γ₀ = 2 * α * N̄ * (M + 1)^(1 / 2) * exp(-M / v) * ω + Γ₀ = 2 * α * N̄ * (M + 1)^(1 / 2) * exp(-M * ω / v) * ω^2 # NB: Kadanoff1963 uses ħ = ω = mb = 1 units. # Factor of omega to get it as a rate relative to phonon frequency @@ -542,20 +537,20 @@ function inverse_Hellwarth_mobility(v, w, α, ω, β) 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) + 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, Inf, rtol=1e-3)[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 + RHS = α / (3 * sqrt(π)) * (β)^(5 / 2) / sinh(ω * β / 2) * (v^3 / w^3) * K + μ = RHS * ω^(3/2) # 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 + K_0 = quadgk(u -> k(u, a, 0, v), 0, Inf, rtol=1e-3)[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 * ω^(3/2) return μ, μ_0 end diff --git a/src/HolsteinPolaron.jl b/src/HolsteinPolaron.jl index 5010e23..b1a9736 100644 --- a/src/HolsteinPolaron.jl +++ b/src/HolsteinPolaron.jl @@ -9,9 +9,13 @@ Type for storing the Holstein polaron information, `x...`. mutable struct HolsteinPolaron α # Holstein coupling. αeff # Effective Holstein coupling summed for multiple modes. + mb # Band mass. T # Temperature. ω # Phonon mode frequencies. ωeff # Effective phonon mode frequency. + γ # Adiabaticity + J # Transfer integral. + a # Lattice constant. β # Reduced thermodynamic beta ħω₀/kBT. Ω # Electric field frequency. d # Number of spatial dimensions. @@ -47,7 +51,7 @@ mutable struct HolsteinPolaron end end -function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesses=4.0, w_guesses=3.9, dims=3, J=1, a=1, kspace=false, verbose=false) +function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=1, a=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. @@ -60,14 +64,6 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse num_d = length(dims) ω = reduce_array(ω) - - Trange = phustrip.(Trange .* 1u"K") - Ωrange = phustrip.(Ωrange .* 1u"THz") - ω = phustrip.(ω .* J .* u"meV" / Unitful.ħ) - ωeff = phustrip.(ωeff .* J .* u"meV" / Unitful.ħ) - - a0 = phustrip.(sqrt(Unitful.ħ / (2 * Unitful.me * mb * ωeff * 1u"THz2π"))) - # For multiple variational modes, ensure that the number of v and w parameters is the same. @assert length(v_guesses) == length(w_guesses) "v and w guesses must be the same length." @@ -77,9 +73,13 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse p = Dict( "α" => αrange, # Alphas. "αeff" => sum(αrange, dims=2), # Alphas sums. + "mb" => mb, # Band mass "T" => Trange, # Temperatures. "ω" => ω, # Phonon frequencies. "ωeff" => ωeff, # Effective Phonon frequency. + "γ" => γ, # Adiabaticity + "J" => J, # Transfer Integral. + "a" => a, # Lattice constant. "β" => Vector{Float64}(undef, num_T), # Betas. "Ω" => Ωrange, # Photon frequencies. "d" => dims, # Number of spatial dimensions. @@ -165,8 +165,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse v_guesses, w_guesses = v_gs, w_gs # Store the athermal data. - p["v0"][d, j, :] .= v_gs - p["w0"][d, j, :] .= w_gs + p["v0"][d, j, :] .= v_gs ./ ωeff + p["w0"][d, j, :] .= w_gs ./ ωeff p["F0"][d, j] = F_gs p["A0"][d, j] = A_gs p["B0"][d, j] = B_gs @@ -177,8 +177,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse println(io, "\e[K-----------------------------------------------------------------------") println(io, "\e[K Zero Temperature Information: ") println(io, "\e[K-----------------------------------------------------------------------") - println(io, "\e[KVariational parameter | v0 = ", v_gs) - println(io, "\e[KVariational parameter | w0 = ", w_gs) + println(io, "\e[KVariational parameter | v0 = ", v_gs ./ ωeff) + println(io, "\e[KVariational parameter | w0 = ", w_gs ./ ωeff) println(io, "\e[KEnergy | F0 = ", F_gs) println(io, "\e[KElectron energy | A0 = ", A_gs) println(io, "\e[KInteraction energy | B0 = ", B_gs) @@ -187,8 +187,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal - κ_gs = (v_gs .^ 2 .- w_gs .^ 2) .* mb - p["κ0"][d, j, :] .= κ_gs + κ_gs = (v_gs .^ 2 .- w_gs .^ 2) + p["κ0"][d, j, :] .= κ_gs ./ ωeff .^2 # Print athermal fictitious spring constant. if verbose @@ -205,7 +205,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Athermal - M_asymp_gs = v_gs ./ w_gs .* mb + M_asymp_gs = v_gs ./ w_gs p["M0a"][d, j, :] .= M_asymp_gs # Print athermal asymptotic fictitious mass. @@ -214,7 +214,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # 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 .* mb + M_reduced_gs = (v_gs .^2 .- w_gs .^2) / v_gs .^ 2 p["M0r"][d, j, :] .= M_reduced_gs # Print athermal reduced mass. @@ -223,8 +223,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # 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) .* a0 - p["R0"][d, j, :] .= R_gs + R_gs = sqrt.(3 ./ 2 .* M_reduced_gs .* v_gs) + p["R0"][d, j, :] .= R_gs ./ ωeff .^(1/2) # Print athermal polaron radius. if verbose @@ -263,8 +263,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse v_guesses, w_guesses = v, w # Store thermal data. - p["v"][i, d, j, :] .= v - p["w"][i, d, j, :] .= w + p["v"][i, d, j, :] .= v ./ ωeff + p["w"][i, d, j, :] .= w ./ ωeff p["F"][i, d, j] = F p["A"][i, d, j] = A p["B"][i, d, j] = B @@ -272,8 +272,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse # Print thermal data. if verbose - println(io, "\e[KVariational parameter | v = ", v) - println(io, "\e[KVariational parameter | w = ", w) + println(io, "\e[KVariational parameter | v = ", v ./ ωeff) + println(io, "\e[KVariational parameter | w = ", w ./ ωeff) println(io, "\e[KFree energy | F = ", F) println(io, "\e[KElectron energy | A = ", A) println(io, "\e[KInteraction energy | B = ", B) @@ -281,8 +281,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # Calculate and store fictitious spring constants. See after Eqn. (18), pg. 1007 of Feynman1962. Thermal - κ = (v .^ 2 .- w .^ 2) .* mb - p["κ"][i, d, j, :] .= κ + κ = (v .^ 2 .- w .^ 2) + p["κ"][i, d, j, :] .= κ ./ ωeff .^2 # Print spring constants. if verbose @@ -299,7 +299,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # Approximate large coupling asymptotic limit for the polaron mass. Feynman1962. Thermal - M_asymp = v ./ w .* mb + M_asymp = v ./ w p["Ma"][i, d, j, :] .= M_asymp # Print asymptotic masses. @@ -308,7 +308,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # Reduced mass of particle and fictitious mass system. Before eqn. (2.4) in Schultz1959. Athermal. - M_reduced = (v .^2 .- w .^2) / v .^ 2 .* mb + M_reduced = (v .^2 .- w .^2) / v .^ 2 p["Mr"][i, d, j, :] .= M_reduced # Print redcued masses. @@ -317,8 +317,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # Calculate and store polaron radii. - R = sqrt.(3 .* v ./ (v .^ 2 .- w .^ 2) .^ 2) .* a0 - p["R"][i, d, j, :] .= R + R = sqrt.(3 ./ 2 .* v .* M_reduced) + p["R"][i, d, j, :] .= R ./ ωeff .^(1/2) # Print polaron radius. if verbose @@ -332,7 +332,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse end # Calculate and store the DC mobiliies. - μ = !kspace ? holstein_mobility(v, w, α, ω, β; dims = dims[d]) * a^2 * mb : holstein_mobility_k_space(v, w, α, ω, β; dims = dims[d]) * a^2 * mb + μ = !kspace ? holstein_mobility(v, w, α, ω, β; dims = dims[d]) : holstein_mobility_k_space(v, w, α, ω, β; dims = dims[d]) p["μ"][i, d, j] = μ # Print DC mobilities. @@ -344,8 +344,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse # 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["v"][i, d, j, :] .= v_gs ./ ωeff + p["w"][i, d, j, :] .= w_gs ./ ωeff p["F"][i, d, j] = F_gs p["A"][i, d, j] = A_gs p["B"][i, d, j] = B_gs @@ -366,6 +366,8 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse for k in eachindex(Ωrange) # E-field frequencies loop. Ω = Ωrange[k] + if !iszero(T) || !iszero(Ω) + # Print E-field frequency. if verbose println(io, "\e[K-----------------------------------------------------------------------") @@ -385,7 +387,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse # Calculate and store polaron complex impedances. - z = -(im * Ω + im * χ) .* mb + z = -(im * Ω + im * χ) p["z"][k, i, d, j] = z # Print complex impedances. @@ -402,7 +404,7 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse println(io, "\e[KComplex conductivity | σ = ", σ) end - if iszero(T) && iszero(Ω) + else # If zero frequency. p["χ"][k, i, d, j] = Inf + 0 * im @@ -419,14 +421,14 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse print(io, "\e[2F") end - if verbose print(io, "\e[7F") end + if verbose && (!iszero(T) || !iszero(Ω)) print(io, "\e[7F") end end if verbose && !iszero(T) print(io, "\e[20F") end # Move up 20 lines and erase. end if verbose print(io, "\e[15F") end # Move up 20 lines and erase. end # End dimensions loop. - if verbose print(io, "\e[5F") end + if verbose print(io, "\e[5F") end end if verbose print("\e[?25h") end # Show cursor again. @@ -434,10 +436,14 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, v_guesse polaron_data = [ p["α"], # Alphas. p["αeff"], # Alphas sums. + p["mb"], # Band mass. p["T"], # Temperatures. p["ω"], # Phonon frequencies. p["ωeff"], # Hellwarth eff phonon frequency. - p["β"], # Betas. + p["γ"], # Adiabaticity + p["J"], # Transfer integral. + p["a"], # Lattice constant. + p["β"], # Inverse temperatures. p["Ω"], # Photon frequencies. p["d"], # Number of spatial dimensions. p["v0"], # Athermal v params. @@ -475,22 +481,22 @@ end """ Single alpha parameter. holstein() expects alpha parameters to be in a Vector. """ -holsteinpolaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, J=J, a=a,v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) +holsteinpolaron(α::Real, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, γ=γ, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No frequency input => Phonon peak. """ -holsteinpolaron(αrange, Trange; ω=1, ωeff=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(αrange, Trange, 0; ω=ω, ωeff=ωeff, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) +holsteinpolaron(αrange, Trange; ω=1, ωeff=1, γ=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(αrange, Trange, 0; ω=ω, ωeff=ωeff, γ=γ, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No temperature input => 300 K. """ -holsteinpolaron(αrange; ω=1, ωeff=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(αrange, 0, 0; ω=ω, ωeff=ωeff, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) +holsteinpolaron(αrange; ω=1, ωeff=1, γ=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(αrange, 0, 0; ω=ω, ωeff=ωeff, γ=γ, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) """ No input => α = 1 """ -holsteinpolaron(; ω=1, ωeff=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(1, 0, 0; ω=ω, ωeff=ωeff, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) +holsteinpolaron(; ω=1, ωeff=1, γ=1, mb=1, J=1, a=1, v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) = holsteinpolaron(1, 0, 0; ω=ω, ωeff=ωeff, γ=γ, mb=mb, J=J, a=a, 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) @@ -498,7 +504,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 holsteinpolaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, kspace=false, verbose=false) +function holsteinpolaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, kspace=false, verbose=false) # Show material data. if verbose @@ -506,12 +512,18 @@ function holsteinpolaron(material::Material, TΩrange...; v_guesses=3.11, w_gues end # Extract material data from Material type. - phonon_freqs = material.f - phonon_eff_freq = material.feff + ω = material.f + ωeff = material.feff mb = material.mb + J = material.J + γ = material.γ + a = material.a + dims = material.z ./ 2 + + TΩrange = length(TΩrange) == 1 ? TΩrange .* phustrip(1u"K") : TΩrange .* (phustrip(1u"K"),1) # Generate Holstein polaron data from the arbitrary model constructor. - p = holsteinpolaron(material.α', TΩrange...; ω=phonon_freqs, ωeff=phonon_eff_freq, mb=mb, β0=ħ/kB*1e12, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) + p = holsteinpolaron(material.α', TΩrange...; ω=ω, ωeff=ωeff, γ=γ, mb=mb, J=J, a=a, v_guesses=v_guesses, w_guesses=w_guesses, dims=dims, kspace=kspace, verbose=verbose) # Return material-specific, unitful Holstein type. return p @@ -528,7 +540,10 @@ function Base.show(io::IO, ::MIME"text/plain", x::HolsteinPolaron) println(io_lim, "\e[KPhonon frequencies | ωeff = ", x.ωeff, " | ω = ", x.ω) println(io_lim, "\e[KHolstein coupling | αeff = ", x.αeff, " | α = ", x.α) + println(io_lim, "\e[KAdiabaticity | γ = ", x.γ) println(io_lim, "\e[KNumber of spatial dimensions | d = ", x.d) + println(io_lim, "\e[KTransfer integral | J = ", x.J) + println(io_lim, "\e[KLattice constant | a = ", x.a) println("\e[K-----------------------------------------------------------------------") println("\e[K Zero Temperature Information: ") @@ -596,9 +611,13 @@ function save_holstein_polaron(polaron::HolsteinPolaron, prefix) JLD.save("$prefix.jld", "alpha", polaron.α, "alpha eff", polaron.αeff, + "band mass", polaron.mb, "temperature", polaron.T, "phonon freq", polaron.ω, "phonon freq eff", polaron.ωeff, + "adiabaticity", polaron.γ, + "transfer integral", polaron.J, + "lattice constant", polaron.a, "beta", polaron.β, "E-field freq", polaron.Ω, "dims", polaron.d, @@ -647,9 +666,13 @@ function load_holstein_polaron(polaron_file_path) holstein_polaron = HolsteinPolaron( data["alpha"], data["alpha eff"], + data["band mass"], data["temperature"], data["phonon freq"], data["phonon freq eff"], + data["adiabaticity"], + data["transfer integral"], + data["lattice constant"], data["beta"], data["E-field freq"], data["dims"], diff --git a/src/HolsteinTheory.jl b/src/HolsteinTheory.jl index 448008a..738ab68 100644 --- a/src/HolsteinTheory.jl +++ b/src/HolsteinTheory.jl @@ -24,7 +24,7 @@ Expected Output: 6.0 """ function holstein_coupling(k, α, ω; dims = 3) - 2 * α * dims / ω + 2 * α * dims * ω end holstein_coupling(k, α::Vector, ω::Vector; dims = 3) = sum(holstein_coupling(k, α[j], ω[j]; dims = dims) for j in eachindex(α)) @@ -58,13 +58,14 @@ 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) + polaron_radius = 1 / sqrt(2) momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) - coupling = holstein_coupling(1, α, ωβ[1]; dims = dims) - propagator(τ) = polaron_propagator(τ, v, w, ωβ...) * 2 / ωβ[1] - integrand(τ) = phonon_propagator(τ, ωβ...) * P(dims, momentum_cutoff^2 * propagator(τ)) / (propagator(τ))^(dims/2) + coupling = holstein_coupling(1, α, ωβ[1]; dims = dims) * ωβ[1] + propagator(τ) = length(ωβ) == 1 ? polaron_propagator(τ, v, w) * ωβ[1] : polaron_propagator(τ, v, w, ωβ[2]) * ωβ[1] + integrand(τ) = phonon_propagator(τ, ωβ...) * P(dims, momentum_cutoff^2 * polaron_radius^2 * propagator(τ)) / (propagator(τ))^(dims/2) upper_limit = length(ωβ) == 1 ? Inf : ωβ[2] / 2 integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) - return integral * coupling / (4π)^(dims / 2) * dims + return integral * coupling / (4π * polaron_radius^2)^(dims / 2) end holstein_interaction_energy(v, w, α::Vector, ω::Vector, β; dims = 3) = sum(holstein_interaction_energy(v, w, α[j], ω[j], β; dims = dims) for j in eachindex(α)) @@ -83,21 +84,18 @@ Calculate the Holstein polaron interaction energy in k-space at finite temperaur - `ω`: 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) momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) - coupling(k) = holstein_coupling(k, α, ωβ[1]; dims = dims) - propagator(τ) = polaron_propagator(τ, v, w, ωβ...) * 4 + coupling(k) = holstein_coupling(k, α, ωβ[1]; dims = dims) * ωβ[1] + propagator(τ) = length(ωβ) == 1 ? polaron_propagator(τ, v, w) * ωβ[1] : polaron_propagator(τ, v, w, ωβ[2]) * ωβ[1] integrand(τ) = phonon_propagator(τ, ωβ...) * spherical_k_integral(coupling, propagator(τ); dims = dims, limits = [0, momentum_cutoff]) upper_limit = length(ωβ) == 1 ? Inf : ωβ[2] / 2 - integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) * dims - return integral + integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit) + return integral end holstein_interaction_energy_k_space(v, w, α::Vector, ω::Vector, β; dims = 3) = sum(holstein_interaction_energy_k_space(v, w, α[j], ω[j], β; dims = dims) for j in eachindex(α)) @@ -107,9 +105,9 @@ holstein_interaction_energy_k_space(v, w, α::Vector, ω::Vector; dims = 3) = su Total free energy for the Holstein model. Here the k-space integral is evaluated analytically. """ function holstein_energy(v, w, α, ωβ...; dims = 3) - A, C = trial_energy(v, w, ωβ...; dims = dims) + A, C = length(ωβ) == 1 ? trial_energy(v, w; dims = dims) : trial_energy(v, w, ωβ[2]; dims = dims) B = holstein_interaction_energy(v, w, α, ωβ...; dims = dims) - return -(A + B + C) / dims - 2 * dims, A / dims, B / dims, C / dims + return -(A + B + C), A, B, C end """ @@ -134,9 +132,9 @@ Calculate the total energy, kinetic energy, and interaction energy of the Holste - `interaction_energy`: The calculated polaron interaction energy. """ function holstein_energy_k_space(v, w, α, ωβ...; dims = 3) - A, C = trial_energy(v, w, ωβ...; dims = dims) + A, C = length(ωβ) == 1 ? trial_energy(v, w; dims = dims) : trial_energy(v, w, ωβ[2]; dims = dims) B = holstein_interaction_energy_k_space(v, w, α, ωβ...; dims = dims) - return -(A + B + C) / dims - 2 * dims, A, B, C + return -(A + B + C), A, B, C end """ @@ -171,10 +169,10 @@ This example demonstrates how to use the `holstein_structure_factor` function to """ function holstein_structure_factor(t, v, w, α, ωβ...; dims = 3) momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) - coupling = holstein_coupling(1, α, ωβ[1]; dims = dims) - propagator = polaron_propagator(im * t, v, w, ωβ...) * 2 / ωβ[1] - integral = dims / 2 * ball_surface(dims) / (2π)^dims * P_plus_one(dims, propagator * momentum_cutoff^2) / propagator^(dims/2 + 1) - return 2 / dims * ωβ[1] * phonon_propagator(im * t, ωβ...) * coupling * integral + coupling = holstein_coupling(1, α, ωβ[1]; dims = dims) * ωβ[1] + propagator = length(ωβ) == 1 ? polaron_propagator(im * t, v, w) * ωβ[1] : polaron_propagator(im * t, v, w, ωβ[2]) * ωβ[1] + integral = P(dims + 2, propagator * momentum_cutoff^2 / 2) * (propagator)^(-dims/2 - 1) * 2^(1 - dims / 2) * π^(-dims / 2) * dims / 2 + return 2 / dims * coupling * integral * phonon_propagator(im * t, ωβ...) end """ @@ -199,10 +197,10 @@ A scalar value representing the calculated structure factor in k-space for the H """ 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, ωβ...) * 4 / ωβ[1] - integral = spherical_k_integral(coupling_one, propagator; dims = dims, limits = [0, momentum_cutoff]) - return 2 * phonon_propagator(im * t, ωβ...) * integral + coupling(k) = holstein_coupling(1, α, ωβ[1]; dims = dims) * k^2 * ωβ[1] + propagator = length(ωβ) == 1 ? polaron_propagator(im * t, v, w) * ωβ[1] : polaron_propagator(im * t, v, w, ωβ[2]) * ωβ[1] + integral = spherical_k_integral(coupling, propagator; dims = dims, limits = [0, momentum_cutoff]) + return 2 / dims * integral * phonon_propagator(im * t, ωβ...) end """ @@ -226,7 +224,7 @@ In this example, the `holstein_memory_function` is called with the input paramet """ function holstein_memory_function(Ω, v, w, α, ωβ...; dims = 3) structure_factor(t) = holstein_structure_factor(t, v, w, α, ωβ...; dims = dims) - return polaron_memory_function(Ω, structure_factor; limits = [0, 1e5]) + return polaron_memory_function(Ω, structure_factor) end holstein_memory_function(Ω, v, w, α::Vector, ω::Vector, β; dims = 3) = sum(holstein_memory_function(Ω, v, w, α[j], ω[j], β; dims = dims) for j in eachindex(α)) @@ -254,7 +252,7 @@ A scalar value representing the memory function of the Holstein model in k-space """ function holstein_memory_function_k_space(Ω, v, w, α, ωβ...; dims = 3) structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ωβ...; dims = dims) - return polaron_memory_function(Ω, structure_factor; limits = [0, 1e5]) + return polaron_memory_function(Ω, structure_factor) end holstein_memory_function_k_space(Ω, v, w, α::Vector, ω::Vector, β; dims = 3) = sum(holstein_memory_function_k_space(Ω, v, w, α[j], ω[j], β; dims = dims) for j in eachindex(α)) @@ -291,11 +289,8 @@ println(result) This code calculates the mobility using the given parameters and prints the result. """ function inverse_holstein_mobility(v, w, α, ω, β; dims = 3) - if β == Inf - return Inf - end structure_factor(t) = holstein_structure_factor(t, v, w, α, ω, β; dims = dims) - return abs(imag(polaron_memory_function(structure_factor; limits = [0, 1e4]))) + return abs(imag(polaron_memory_function(structure_factor))) end inverse_holstein_mobility(v, w, α::Vector, ω::Vector, β; dims = 3) = sum(inverse_holstein_mobility(v, w, α[j], ω[j], β; dims = dims) for j in eachindex(α)) @@ -332,15 +327,15 @@ Calculate the DC mobility in k-space for a Holstein polaron system at finite tem The DC mobility in k-space for the Holstein polaron system at finite temperature. """ function inverse_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) - return abs(imag(polaron_memory_function(structure_factor; limits = [0, 1e4]))) + return abs(imag(polaron_memory_function(structure_factor))) end inverse_holstein_mobility_k_space(v, w, α::Vector, ω::Vector, β; dims = 3) = sum(inverse_holstein_mobility_k_space(v, w, α[j], ω[j], β; dims = dims) for j in eachindex(α)) +holstein_mobility_k_space(v, w, α, ω, β; dims = 3) = 1 / inverse_holstein_mobility_k_space(v, w, α, ω, β; dims = dims) + +holstein_mobility_k_space(v, w, α, ω; dims = 3) = reduce_array(repeat([Inf], length(α))) function holstein_complex_impedence(Ω, v, w, α, ωβ...; dims = 3) -im * (Ω - holstein_memory_function(v, w, α, ωβ...; dims = dims)) diff --git a/src/LegacyFunctions.jl b/src/LegacyFunctions.jl index aeed91c..3771b90 100644 --- a/src/LegacyFunctions.jl +++ b/src/LegacyFunctions.jl @@ -15,77 +15,82 @@ Minimises the multiple phonon mode free energy function for a set of vₚ and w 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 +feynmanvw(αωβ...) = vw_variation((v,w)->frohlich_energy(v, w, αωβ...)) - N_params = length(v) +feynmanvw(α) = feynmanvw(α, 1) - Δv = v .- w - initial = vcat(Δv .+ repeat([eps(Float64)], N_params), w) +# function feynmanvw(v::Vector, w::Vector, αωβ...; upper_limit=1e6) - # Limits of the optimisation. - lower = fill(0.0, 2 * N_params) - upper = fill(upper_limit, 2 * N_params) +# if length(v) != length(w) +# return error("The number of variational parameters v & w must be equal.") +# end - # 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] +# N_params = length(v) - # 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()) - ) +# Δv = v .- w +# initial = vcat(Δv .+ repeat([eps(Float64)], N_params), w) - # Extract the v and w parameters that minimised the free energy. - var_params = Optim.minimizer(solution) +# # Limits of the optimisation. +# lower = fill(0.0, 2 * N_params) +# upper = fill(upper_limit, 2 * N_params) - # 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, αωβ...) +# # 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] - # if Optim.converged(solution) == false - # @warn "Failed to converge. v = $(Δv .+ w), w = $w, E = $E" - # end +# # 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()) +# ) - # Return the variational parameters that minimised the free energy. - return Δv .+ w, w, E, A, B, C -end +# # Extract the v and w parameters that minimised the free energy. +# var_params = Optim.minimizer(solution) -function feynmanvw(v::Real, w::Real, αωβ...; upper = [Inf, Inf], lower = [0, 0]) +# # 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, αωβ...) +# Δ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 3781af9..d83b244 100644 --- a/src/Material.jl +++ b/src/Material.jl @@ -8,22 +8,35 @@ """ mutable struct Material - ϵo # Optical dielectric constant - ϵs # Static dielectric constant - ϵi # Ionic dielectric contributions - mb # Effective band mass - α # Fröhlich coupling - f # Phonon frequencies - feff # Effective frequency - ir # Infrared activities - V # Unit cell volumes + # General polaron material properties + α # Unitless electron-phonon coupling (Unitless) + mb # Effective band mass (me) + f # Phonon frequencies (THz2π) + feff # Effective frequency (THz2π) + + # Frohlich specific properties for calculating electron-phonon coupling α + ϵo # Optical dielectric constant + ϵs # Static dielectric constant + ϵi # Ionic dielectric contributions + + # Frohlich specific properties for calculating multiple phonon electron-phonon coupling αⱼ + ir # Infrared activities + V # Unit cell volumes + + # Holstein specific properties + γ # Adiabaticity (Unitless) + J # Transfer integral (meV) + a # Lattice constant (Angstroms) + g # Electron-phonon coupling energy (Unitless) + z # Lattice dimensionality (Unitless) + function Material(x...) new(reduce_array.(x)...) end end """ - material(ϵ_optical, ϵ_static, m_eff, phonon_freq) +material(ϵ_optical, ϵ_static, m_eff, phonon_freq) Construct a 'struct Material' object from traditional Frohlich Hamiltonian parameters. @@ -32,13 +45,13 @@ end function material(ϵ_optical, ϵ_static, m_eff, phonon_freq) ϵ_ionic = ϵ_static - ϵ_optical α = frohlichalpha(ϵ_optical, ϵ_static, phonon_freq, m_eff) - return Material(ϵ_optical, ϵ_static, ϵ_ionic, m_eff, α, phonon_freq, phonon_freq, 1, 1) + return Material(α, m_eff, phonon_freq, phonon_freq, ϵ_optical, ϵ_static, ϵ_ionic, 1, 1, 1, 1, 1, 1, 1) end """ - material(ϵ_optic, ϵ_static, m_eff, phonon_freqs, ir_activity, volume) +material(ϵ_optic, ϵ_static, m_eff, phonon_freqs, ir_activity, volume) - Construct a 'struct Material' object from a set of infrared activities and phonon + Construct a 'struct Material' object, for a multiple phonon mode Frohlich polaron, from a set of infrared activities and phonon frequencies, for use with the 'multiple phonon branches' extension of the code. effective_freq calculated with the Hellwarth 'B' scheme; ϵ_ionic and α calculated from @@ -48,7 +61,22 @@ function material(ϵ_optic, ϵ_static, m_eff, phonon_freqs, ir_activity, volume) effective_freq = HellwarthBScheme(hcat(phonon_freqs, ir_activity)) ϵ_ionic = ϵ_ionic_mode.(phonon_freqs, ir_activity, volume) α = frohlichalpha.(ϵ_optic, ϵ_ionic, sum(ϵ_ionic), phonon_freqs, m_eff) - return Material(ϵ_optic, ϵ_static, ϵ_ionic, m_eff, α, phonon_freqs, effective_freq, ir_activity, volume) + return Material(α, m_eff, phonon_freqs, effective_freq, ϵ_optic, ϵ_static, ϵ_ionic, ir_activity, volume, 1, 1, 1, 1, 1) +end + +function material(transfer_integral, coupling_energy, phonon_freq, lattice_constant, lattice_dimensionality) + + transfer_integral *= eV / 1000 + coupling_energy *= eV / 1000 + lattice_constant *= 1e-10 + + phonon_energy = ħ * phonon_freq * 1e12 * 2π + dimensionless_coupling = coupling_energy / phonon_energy + adiabaticity = phonon_energy / transfer_integral + α = dimensionless_coupling^2 * adiabaticity ./ lattice_dimensionality + band_mass = ħ^2 / transfer_integral / lattice_constant^2 / 2 / me + + return Material(α, band_mass, phonon_freq, phonon_freq, 1, 1, 1, 1, 1, adiabaticity, transfer_integral / eV * 1000, lattice_constant * 1e10, dimensionless_coupling, lattice_dimensionality) end function Base.show(io::IO, ::MIME"text/plain", x::Material) @@ -56,15 +84,21 @@ function Base.show(io::IO, ::MIME"text/plain", x::Material) println("\e[K------------------------------------------") println("\e[K Material Information ") println("\e[K------------------------------------------") - println(io_limit, "\e[KOptic dielectric | ϵo = ", x.ϵo) - println(io_limit, "\e[KStatic dielectric | ϵs = ", x.ϵs) - println(io_limit, "\e[KIonic dielectric | ϵi = ", x.ϵi) + println(io_limit, "\e[KUnitless Coupling | α = ", x.α) println(io_limit, "\e[KBand mass | mb = ", x.mb) - println(io_limit, "\e[KFröhlich coupling | α = ", x.α) println(io_limit, "\e[KPhonon frequencies | f = ", x.f) println(io_limit, "\e[KEff Phonon freq | feff = ", x.feff) + println(io_limit, "\e[KOptic dielectric | ϵo = ", x.ϵo) + println(io_limit, "\e[KStatic dielectric | ϵs = ", x.ϵs) + println(io_limit, "\e[KIonic dielectric | ϵi = ", x.ϵi) println(io_limit, "\e[KIR activities | ir = ", x.ir) println(io_limit, "\e[KUnit cell volume | V = ", x.V) + println(io_limit, "\e[KAdiabaticity | γ = ", x.γ) + println(io_limit, "\e[KTransfer integral | J = ", x.J) + println(io_limit, "\e[KLattice constant | a = ", x.a) + println(io_limit, "\e[KEl-ph coupling | g = ", x.g) + println(io_limit, "\e[KLattice Dimensions | z = ", x.z) + println("\e[K-------------------------------------------") end diff --git a/src/PolaronFunctions.jl b/src/PolaronFunctions.jl index dbf202d..d4398ae 100644 --- a/src/PolaronFunctions.jl +++ b/src/PolaronFunctions.jl @@ -89,7 +89,7 @@ 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 = [1e3, 1e3], warn = false) +function vw_variation(energy, initial_v, initial_w; lower = [0, 0], upper = [Inf, Inf], warn = false) Δv = initial_v - initial_w # defines a constraint, so that v>w initial = [Δv + eps(Float64), initial_w] @@ -119,9 +119,9 @@ function vw_variation(energy, initial_v, initial_w; lower = [0, 0], upper = [1e3 return Δv + w, w, energy_minimized... end -vw_variation(energy) = vw_variation(energy, 5, 3; lower = [0, 0], upper = [1e3, 1e3]) +vw_variation(energy) = vw_variation(energy, 3.1, 3; lower = [0, 0], upper = [Inf, Inf]) -function vw_variation(energy, initial_v::Vector, initial_w::Vector; lower = [0, 0], upper = [1e4, 1e4], warn = false) +function vw_variation(energy, initial_v::Vector, initial_w::Vector; lower = [0, 0], upper = [Inf, Inf], warn = false) if length(initial_v) != length(initial_w) return error("The number of variational parameters v & w must be equal.") @@ -198,7 +198,7 @@ function polaron_memory_function(Ω, structure_factor; limits = [0, Inf]) if iszero(Ω) return polaron_memory_function(structure_factor; limits = limits) end - integral, _ = quadgk(t -> (1 - exp(im * Ω * t)) / Ω * imag(structure_factor(t)), limits[1], limits[2]) + integral, _ = quadgk(t -> (1 - exp(im * Ω * t)) / Ω * imag(structure_factor(t)), 0, Inf) return integral end @@ -225,11 +225,11 @@ 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 +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]) + integral, _ = quadgk(t -> -im * t * imag(structure_factor(t)), eps(Float64), Inf) return integral end @@ -264,8 +264,8 @@ function ball_surface(n) 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) +function spherical_k_integral(coupling, propagator; dims = 3, radius = 1/sqrt(2), limits = [0, Inf]) + integrand(k) = k^(dims-1) * coupling(k) * exp(-k^2 * radius^2 * propagator) 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/PolaronUnits.jl b/src/PolaronUnits.jl index a2e41ed..e9b3d37 100644 --- a/src/PolaronUnits.jl +++ b/src/PolaronUnits.jl @@ -228,50 +228,50 @@ pustrip(x) = ustrip(puconvert(x)) addunits!(polaron::Polaron) """ function addunits!(polaron::FrohlichPolaron; unit="pu") - polaron.ω = pustrip.(polaron.ω) .* punit(1Unitful.THz2π) - polaron.ωeff = pustrip.(polaron.ωeff) .* punit(1Unitful.THz2π) - polaron.Fs = pustrip.(polaron.Fs) .* punit(1Unitful.meV) - polaron.Fl = pustrip.(polaron.Fl) .* punit(1Unitful.meV) - polaron.Ms = pustrip.(polaron.Ms) .* punit(1Unitful.me) - polaron.Ml = pustrip.(polaron.Ml) .* punit(1Unitful.me) - polaron.Rs = pustrip.(polaron.Rs) .* punit(1Unitful.Å) - polaron.Rl = pustrip.(polaron.Rl) .* punit(1Unitful.Å) - polaron.ΩFC = pustrip.(polaron.ΩFC) .* punit(1Unitful.THz2π) + polaron.Fs = pustrip.(polaron.Fs .* polaron.ωeff) .* punit(1Unitful.meV) + polaron.Fl = pustrip.(polaron.Fl .* polaron.ωeff) .* punit(1Unitful.meV) + polaron.Ms = pustrip.(polaron.Ms .* polaron.mb) .* punit(1Unitful.me) + polaron.Ml = pustrip.(polaron.Ml .* polaron.mb) .* punit(1Unitful.me) + polaron.Rs = pustrip.(polaron.Rs ./ sqrt.(2 .* polaron.mb .* polaron.ωeff.* 2π)) .* punit(1Unitful.Å) + polaron.Rl = pustrip.(polaron.Rl ./ sqrt.(2 .* polaron.mb .* polaron.ωeff .* 2π)) .* punit(1Unitful.Å) + polaron.ΩFC = pustrip.(polaron.ΩFC) .* punit(1Unitful.THz) polaron.F0 = pustrip.(polaron.F0) .* punit(1Unitful.meV) polaron.A0 = pustrip.(polaron.A0) .* punit(1Unitful.meV) polaron.B0 = pustrip.(polaron.B0) .* punit(1Unitful.meV) polaron.C0 = pustrip.(polaron.C0) .* punit(1Unitful.meV) - polaron.v0 = pustrip.(polaron.v0) .* punit(1Unitful.THz2π) - polaron.w0 = pustrip.(polaron.w0) .* punit(1Unitful.THz2π) - polaron.κ0 = pustrip.(polaron.κ0) .* punit(1u"μN/m") - polaron.M0 = pustrip.(polaron.M0) .* punit(1Unitful.me) - polaron.M0a = pustrip.(polaron.M0a) .* punit(1Unitful.me) - polaron.M0r = pustrip.(polaron.M0r) .* punit(1Unitful.me) - polaron.R0 = pustrip.(polaron.R0) .* punit(1Unitful.Å) - polaron.T = pustrip.(polaron.T) .* punit(1Unitful.K) - polaron.β = pustrip.(polaron.β) .* punit(1 / 1Unitful.K) + polaron.v0 = pustrip.(polaron.v0 ./ 2π) .* punit(1Unitful.THz) + polaron.w0 = pustrip.(polaron.w0 ./ 2π) .* punit(1Unitful.THz) + polaron.κ0 = pustrip.(polaron.κ0 .* polaron.mb ./ 4π^2) .* punit(1u"μN/m") + polaron.M0 = pustrip.(polaron.M0 .* polaron.mb) .* punit(1Unitful.me) + polaron.M0a = pustrip.(polaron.M0a .* polaron.mb) .* punit(1Unitful.me) + polaron.M0r = pustrip.(polaron.M0r .* polaron.mb) .* punit(1Unitful.me) + polaron.R0 = pustrip.(polaron.R0 ./ sqrt.(2 .* polaron.mb .* 2π)) .* punit(1Unitful.Å) polaron.F = pustrip.(polaron.F) .* punit(1Unitful.meV) polaron.A = pustrip.(polaron.A) .* punit(1Unitful.meV) polaron.B = pustrip.(polaron.B) .* punit(1Unitful.meV) polaron.C = pustrip.(polaron.C) .* punit(1Unitful.meV) - polaron.v = pustrip.(polaron.v) .* punit(1Unitful.THz2π) - polaron.w = pustrip.(polaron.w) .* punit(1Unitful.THz2π) - polaron.κ = pustrip.(polaron.κ) .* punit(1u"μN/m") - polaron.M = pustrip.(polaron.M) .* punit(1Unitful.me) - polaron.Ma = pustrip.(polaron.Ma) .* punit(1Unitful.me) - polaron.Mr = pustrip.(polaron.Mr) .* punit(1Unitful.me) - polaron.R = pustrip.(polaron.R) .* punit(1Unitful.Å) - polaron.μ = pustrip.(polaron.μ) .* punit(1u"cm^2/V/s") - polaron.μFHIP = pustrip.(polaron.μFHIP) .* punit(1u"cm^2/V/s") - polaron.μD = pustrip.(polaron.μD) .* punit(1u"cm^2/V/s") - polaron.μK = pustrip.(polaron.μK) .* punit(1u"cm^2/V/s") - polaron.μH = pustrip.(polaron.μH) .* punit(1u"cm^2/V/s") - polaron.μH0 = pustrip.(polaron.μH0) .* punit(1u"cm^2/V/s") + polaron.v = pustrip.(polaron.v ./ 2π) .* punit(1Unitful.THz) + polaron.w = pustrip.(polaron.w ./ 2π) .* punit(1Unitful.THz) + polaron.κ = pustrip.(polaron.κ .* polaron.mb) .* punit(1u"μN/m") + polaron.M = pustrip.(polaron.M .* polaron.mb) .* punit(1Unitful.me) + polaron.Ma = pustrip.(polaron.Ma .* polaron.mb) .* punit(1Unitful.me) + polaron.Mr = pustrip.(polaron.Mr .* polaron.mb) .* punit(1Unitful.me) + polaron.R = pustrip.(polaron.R ./ sqrt.(2 .* polaron.mb .* 2π)) .* punit(1Unitful.Å) + polaron.μ = pustrip.(polaron.μ ./ polaron.mb) .* punit(1u"cm^2/V/s") + polaron.μFHIP = pustrip.(polaron.μFHIP ./ polaron.mb) .* punit(1u"cm^2/V/s") + polaron.μD = pustrip.(polaron.μD ./ polaron.mb) .* punit(1u"cm^2/V/s") + polaron.μK = pustrip.(polaron.μK ./ polaron.mb) .* punit(1u"cm^2/V/s") + polaron.μH = pustrip.(polaron.μH ./ polaron.mb) .* punit(1u"cm^2/V/s") + polaron.μH0 = pustrip.(polaron.μH0 ./ polaron.mb) .* punit(1u"cm^2/V/s") polaron.τ = pustrip.(polaron.τ) .* punit(1u"ns") - polaron.Ω = pustrip.(polaron.Ω) .* punit(1Unitful.THz) - polaron.χ = pustrip.(polaron.χ) .* punit(1Unitful.THz2π) - polaron.z = pustrip.(polaron.z) .* punit(u"Ω") - polaron.σ = pustrip.(polaron.σ) .* punit(u"S") + polaron.χ = pustrip.(polaron.χ .* polaron.mb ./ 2π) .* punit(u"Ω") .* punit(1Unitful.THz) + polaron.z = pustrip.(polaron.z .* polaron.mb) .* punit(u"Ω") + polaron.σ = pustrip.(polaron.σ ./ polaron.mb) .* punit(u"S") + polaron.T = pustrip.(polaron.T) .* punit(1Unitful.K) + polaron.β = pustrip.(polaron.β ./ Unitful.ħ) .* punit(1 / 1Unitful.meV) + polaron.Ω = pustrip.(polaron.Ω ./ 2π) .* punit(1Unitful.THz) + polaron.ω = pustrip.(polaron.ω ./ 2π) .* punit(1Unitful.THz) + polaron.ωeff = pustrip.(polaron.ωeff ./ 2π) .* punit(1Unitful.THz) if unit == "au" auconvert!(polaron) elseif unit == "su" @@ -280,37 +280,37 @@ function addunits!(polaron::FrohlichPolaron; unit="pu") end function addunits!(polaron::HolsteinPolaron; unit="pu") - polaron.ω = phustrip.(polaron.ω) .* phunit(1Unitful.THz2π) - polaron.ωeff = phustrip.(polaron.ωeff) .* phunit(1Unitful.THz2π) - polaron.F0 = phustrip.(polaron.F0) .* phunit(1Unitful.meV) - polaron.A0 = phustrip.(polaron.A0) .* phunit(1Unitful.meV) - polaron.B0 = phustrip.(polaron.B0) .* phunit(1Unitful.meV) - polaron.C0 = phustrip.(polaron.C0) .* phunit(1Unitful.meV) - polaron.v0 = pustrip.(polaron.v0) .* punit(1Unitful.meV / 1Unitful.ħ) - polaron.w0 = pustrip.(polaron.w0) .* punit(1Unitful.meV / 1Unitful.ħ) - polaron.κ0 = phustrip.(polaron.κ0) .* phunit(1u"μN/m") - polaron.M0 = phustrip.(polaron.M0) .* phunit(1Unitful.me) - polaron.M0a = phustrip.(polaron.M0a) .* phunit(1Unitful.me) - polaron.M0r = phustrip.(polaron.M0r) .* phunit(1Unitful.me) - polaron.R0 = phustrip.(polaron.R0) .* phunit(1Unitful.Å) - polaron.T = phustrip.(polaron.T) .* phunit(1Unitful.K) - polaron.β = phustrip.(polaron.β) .* phunit(1 / 1Unitful.K) - polaron.F = phustrip.(polaron.F) .* phunit(1Unitful.meV) - polaron.A = phustrip.(polaron.A) .* phunit(1Unitful.meV) - polaron.B = phustrip.(polaron.B) .* phunit(1Unitful.meV) - polaron.C = phustrip.(polaron.C) .* phunit(1Unitful.meV) - polaron.v = pustrip.(polaron.v) .* punit(1Unitful.meV / 1Unitful.ħ) - polaron.w = pustrip.(polaron.w) .* punit(1Unitful.meV / 1Unitful.ħ) - polaron.κ = phustrip.(polaron.κ) .* phunit(1u"μN/m") - polaron.M = phustrip.(polaron.M) .* phunit(1Unitful.me) - polaron.Ma = phustrip.(polaron.Ma) .* phunit(1Unitful.me) - polaron.Mr = phustrip.(polaron.Mr) .* phunit(1Unitful.me) - polaron.R = phustrip.(polaron.R) .* phunit(1Unitful.Å) - polaron.μ = phustrip.(polaron.μ) .* phunit(1u"cm^2/V/s") + polaron.F0 = phustrip.(polaron.F0 .* polaron.J) .* phunit(1Unitful.meV) + polaron.A0 = phustrip.(polaron.A0 .* polaron.J) .* phunit(1Unitful.meV) + polaron.B0 = phustrip.(polaron.B0 .* polaron.J) .* phunit(1Unitful.meV) + polaron.C0 = phustrip.(polaron.C0 .* polaron.J) .* phunit(1Unitful.meV) + polaron.v0 = pustrip.(polaron.v0 ./ 2π) .* phunit(1Unitful.THz) + polaron.w0 = pustrip.(polaron.w0 ./ 2π) .* phunit(1Unitful.THz) + polaron.κ0 = phustrip.(polaron.κ0 .* polaron.mb) .* phunit(1u"μN/m") + polaron.M0 = phustrip.(polaron.M0 .* polaron.mb) .* phunit(1Unitful.me) + polaron.M0a = phustrip.(polaron.M0a .* polaron.mb) .* phunit(1Unitful.me) + polaron.M0r = phustrip.(polaron.M0r .* polaron.mb) .* phunit(1Unitful.me) + polaron.R0 = phustrip.(polaron.R0 .* polaron.a ./ sqrt(2 * polaron.γ)) .* phunit(1Unitful.Å) + polaron.F = phustrip.(polaron.F .* polaron.J) .* phunit(1Unitful.meV) + polaron.A = phustrip.(polaron.A .* polaron.J) .* phunit(1Unitful.meV) + polaron.B = phustrip.(polaron.B .* polaron.J) .* phunit(1Unitful.meV) + polaron.C = phustrip.(polaron.C .* polaron.J) .* phunit(1Unitful.meV) + polaron.v = pustrip.(polaron.v ./ 2π) .* phunit(1Unitful.THz) + polaron.w = pustrip.(polaron.w ./ 2π) .* phunit(1Unitful.THz) + polaron.κ = phustrip.(polaron.κ .* polaron.mb) .* phunit(1u"μN/m") + polaron.M = phustrip.(polaron.M .* polaron.mb) .* phunit(1Unitful.me) + polaron.Ma = phustrip.(polaron.Ma .* polaron.mb) .* phunit(1Unitful.me) + polaron.Mr = phustrip.(polaron.Mr .* polaron.mb) .* phunit(1Unitful.me) + polaron.R = phustrip.(polaron.R .* polaron.a ./ sqrt(2 * polaron.γ)) .* phunit(1Unitful.Å) + polaron.μ = phustrip.(polaron.μ .* polaron.a^2 ./ polaron.γ) .* phunit(1u"cm^2/V/s") polaron.Ω = phustrip.(polaron.Ω) .* phunit(1Unitful.THz) - polaron.χ = phustrip.(polaron.χ) .* phunit(1Unitful.THz2π) - polaron.z = phustrip.(polaron.z) .* phunit(u"Ω") - polaron.σ = phustrip.(polaron.σ) .* phunit(u"S") + polaron.χ = phustrip.(polaron.χ .* polaron.mb ./ 2π) .* phunit(u"Ω") .* phunit(1Unitful.THz) + polaron.z = phustrip.(polaron.z .* polaron.mb) .* phunit(u"Ω") + polaron.σ = phustrip.(polaron.σ ./ polaron.mb) .* phunit(u"S") + polaron.T = phustrip.(polaron.T) .* phunit(1Unitful.K) + polaron.β = phustrip.(polaron.β ./ Unitful.ħ) .* phunit(1 / 1Unitful.meV) + polaron.ω = phustrip.(polaron.ω .* polaron.J ./ Unitful.ħ ./ 2π) .* phunit(1Unitful.THz) + polaron.ωeff = phustrip.(polaron.ωeff .* polaron.J ./ Unitful.ħ ./ 2π) .* phunit(1Unitful.THz) if unit == "au" auconvert!(polaron) elseif unit == "su" @@ -404,34 +404,34 @@ function auconvert!(polaron::HolsteinPolaron) end function suconvert!(polaron::FrohlichPolaron) - polaron.ω = polaron.ω .|> Unitful.THz2π - polaron.ωeff = polaron.ωeff .|> Unitful.THz2π + polaron.ω = polaron.ω .|> Unitful.THz + polaron.ωeff = polaron.ωeff .|> Unitful.THz polaron.Fs = polaron.Fs .|> Unitful.meV polaron.Fl = polaron.Fl .|> Unitful.meV polaron.Ms = auconvert.(polaron.Ms) polaron.Ml = auconvert.(polaron.Ml) polaron.Rs = polaron.Rs .|> Unitful.Å polaron.Rl = polaron.Rl .|> Unitful.Å - polaron.ΩFC = polaron.ΩFC .|> Unitful.THz2π + polaron.ΩFC = polaron.ΩFC .|> Unitful.THz polaron.F0 = polaron.F0 .|> Unitful.meV polaron.A0 = polaron.A0 .|> Unitful.meV polaron.B0 = polaron.B0 .|> Unitful.meV polaron.C0 = polaron.C0 .|> Unitful.meV - polaron.v0 = polaron.v0 .|> Unitful.THz2π - polaron.w0 = polaron.w0 .|> Unitful.THz2π + polaron.v0 = polaron.v0 .|> Unitful.THz + polaron.w0 = polaron.w0 .|> Unitful.THz polaron.κ0 = polaron.κ0 .|> u"μN/m" polaron.M0 = auconvert.(polaron.M0) polaron.M0a = auconvert.(polaron.M0a) polaron.M0r = auconvert.(polaron.M0r) polaron.R0 = polaron.R0 .|> Unitful.Å polaron.T = polaron.T .|> Unitful.K - polaron.β = polaron.β .|> u"K^-1" + polaron.β = polaron.β .|> u"meV^-1" polaron.F = polaron.F .|> Unitful.meV polaron.A = polaron.A .|> Unitful.meV polaron.B = polaron.B .|> Unitful.meV polaron.C = polaron.C .|> Unitful.meV - polaron.v = polaron.v .|> Unitful.THz2π - polaron.w = polaron.w .|> Unitful.THz2π + polaron.v = polaron.v .|> Unitful.THz + polaron.w = polaron.w .|> Unitful.THz polaron.κ = polaron.κ .|> u"μN/m" polaron.M = auconvert.(polaron.M) polaron.Ma = auconvert.(polaron.Ma) @@ -445,41 +445,41 @@ function suconvert!(polaron::FrohlichPolaron) polaron.μH0 = polaron.μH0 .|> u"cm^2/V/s" polaron.τ = polaron.τ .|> Unitful.ns polaron.Ω = polaron.Ω .|> Unitful.THz - polaron.χ = polaron.χ .|> Unitful.THz2π + polaron.χ = polaron.χ .|> u"kΩ" .* Unitful.THz polaron.z = polaron.z .|> Unitful.kΩ polaron.σ = polaron.σ .|> Unitful.mS end function suconvert!(polaron::HolsteinPolaron) - polaron.ω = polaron.ω .|> Unitful.THz2π - polaron.ωeff = polaron.ωeff .|> Unitful.THz2π + polaron.ω = polaron.ω .|> Unitful.THz + polaron.ωeff = polaron.ωeff .|> Unitful.THz polaron.F0 = polaron.F0 .|> Unitful.meV polaron.A0 = polaron.A0 .|> Unitful.meV polaron.B0 = polaron.B0 .|> Unitful.meV polaron.C0 = polaron.C0 .|> Unitful.meV - polaron.v0 = polaron.v0 .|> Unitful.THz2π - polaron.w0 = polaron.w0 .|> Unitful.THz2π + polaron.v0 = polaron.v0 .|> Unitful.THz + polaron.w0 = polaron.w0 .|> Unitful.THz polaron.κ0 = polaron.κ0 .|> u"μN/m" - polaron.M0 = auconvert.(polaron.M0) - polaron.M0a = auconvert.(polaron.M0a) - polaron.M0r = auconvert.(polaron.M0r) + # polaron.M0 = auconvert.(polaron.M0) + # polaron.M0a = auconvert.(polaron.M0a) + # polaron.M0r = auconvert.(polaron.M0r) polaron.R0 = polaron.R0 .|> Unitful.Å polaron.T = polaron.T .|> Unitful.K - polaron.β = polaron.β .|> u"K^-1" + polaron.β = polaron.β .|> u"meV^-1" polaron.F = polaron.F .|> Unitful.meV polaron.A = polaron.A .|> Unitful.meV polaron.B = polaron.B .|> Unitful.meV polaron.C = polaron.C .|> Unitful.meV - polaron.v = polaron.v .|> Unitful.THz2π - polaron.w = polaron.w .|> Unitful.THz2π + polaron.v = polaron.v .|> Unitful.THz + polaron.w = polaron.w .|> Unitful.THz polaron.κ = polaron.κ .|> u"μN/m" - polaron.M = auconvert.(polaron.M) - polaron.Ma = auconvert.(polaron.Ma) - polaron.Mr = auconvert.(polaron.Mr) + # polaron.M = auconvert.(polaron.M) + # polaron.Ma = auconvert.(polaron.Ma) + # polaron.Mr = auconvert.(polaron.Mr) polaron.R = polaron.R .|> Unitful.Å polaron.μ = polaron.μ .|> u"cm^2/V/s" polaron.Ω = polaron.Ω .|> Unitful.THz - polaron.χ = polaron.χ .|> Unitful.THz2π + polaron.χ = polaron.χ .|> u"kΩ" .* Unitful.THz polaron.z = polaron.z .|> Unitful.kΩ polaron.σ = polaron.σ .|> Unitful.mS end \ No newline at end of file diff --git a/src/TrialPolaron.jl b/src/TrialPolaron.jl index f3841ab..77069ed 100644 --- a/src/TrialPolaron.jl +++ b/src/TrialPolaron.jl @@ -26,8 +26,18 @@ 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 - τ / β)) * ω + 1e-15 +function polaron_propagator(τ, v, w, β) + c = (v^2 - w^2) / v^3 + result = c * (phonon_propagator(0, v, β) - phonon_propagator(τ, v, β)) + (1 - c * v) * τ * (1 - τ / β) + if real(result) <= 0 && iszero(imag(result)) + return eps(Float64) + else + if isnan(result) + return phonon_propagator(τ, ω) + else + return result + end + end end """ @@ -53,8 +63,14 @@ 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 * τ))) * ω + 1e-15 +function polaron_propagator(τ, v, w) + c = (v^2 - w^2) / v^3 + result = (1 - v * c) * τ + c * (phonon_propagator(0, v) - phonon_propagator(τ, v)) + if real(result) <= 0 && iszero(imag(result)) + return eps(Float64) + else + return result + end end """ @@ -70,8 +86,8 @@ Calculates the recoil function (a generalisation of D(u) in Eqn. (35c) in FHIP 1 See FHIP 1962: https://doi.org/10.1103/PhysRev.127.1004. """ -function polaron_propagator(τ, 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)) * ω + 1e-15 +function polaron_propagator(τ, 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 """ @@ -84,8 +100,8 @@ Calculates the recoil function at zero-temperature. - `v::Vector{Float64}`: is a vector of the v variational parameters. - `w::Vector{Float64}`: is a vector of the w variational parameters. """ -function polaron_propagator(τ, 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)) * ω + 1e-15 +function polaron_propagator(τ, 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 @@ -102,13 +118,10 @@ Hellwarth's A expression from Eqn. (62b) in Hellwarth et al. 1999 PRB. Part of t 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, β) = β == Inf ? A(v, w) : 3 / β * (log(v / w) - (v - w) * β / 2 - log(1 - exp(-v * β)) + log(1 - exp(-w * β))) -A(v, w, ω::Vector, β) = sum(A.(v, w, ω, β)) +A(v, w) = -3 * (v - w) / 2 -A(v, w, ω) = -3 * (v - w) / 2 - -A(v, w, ω::Vector) = sum(A.(v, w, ω)) """ C(v, w, ω, β) @@ -117,13 +130,9 @@ Hellwarth's C expression from Eqn. (62e) in Hellwarth et al. 1999 PRB. Part of t 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^2 - w^2) / v * (coth(v * β / 2) - 2 / (v * β)) -C(v, w, ω) = (3 / (4 * v)) * (v^2 - w^2) - -C(v, w, ω::Vector) = sum(C.(v, w, ω)) / length(ω) +C(v, w) = (3 / (4 * v)) * (v^2 - w^2) # Extending the Feynman theory to multiple phonon branches @@ -188,8 +197,7 @@ Note: Not to be confused with the number of physical phonon branches; many phono 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 + return w[i] * κ_i(i, v, w) * h_i(j, v, w) / (4 * (v[j]^2 - w[i]^2)) end """ @@ -207,16 +215,14 @@ Required for calculating the polaron free energy. See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. """ -function C(v::Vector, w::Vector, ω, β) +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)) + 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 * ω + return 3 * s end -C(v::Vector, w::Vector, ω::Vector, β) = sum(C(v, w, ω[i], β) for i in eachindex(ω)) / length(ω) - """ C(v, w) @@ -226,15 +232,13 @@ Calculates `C(v, w, β)` but at zero-temperature, `β = Inf`. - `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, ω) +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 * ω + return 3 * s end -C(v::Vector, w::Vector, ω::Vector) = sum(C(v, w, ω[i]) for i in eachindex(ω)) / length(ω) - """ A(v, w, β) @@ -249,17 +253,15 @@ Required for calculating the polaron free energy. See Hellwarth, R. W., Biaggio, I. (1999): https://doi.org/10.1103/PhysRevB.60.299. """ -function A(v::Vector, w::Vector, ω, β) +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] * ω * β)) + 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 * ω + 3 / β * s end -A(v::Vector, w::Vector, ω::Vector, β) = sum(A(v, w, ω[i], β) for i in eachindex(ω)) / length(ω) - """ A(v, w, n) @@ -269,13 +271,11 @@ Calculates `A(v, w, β)` but at zero-temperature, `β = Inf`. - `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, ω) +function A(v::Vector, w::Vector) s = sum(v .- w) - return -3 * s / 2 * ω + 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, ω, β) @@ -302,8 +302,8 @@ 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 +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