From 0b9f8005d6dcde18142b2980d707b6981d2035ac Mon Sep 17 00:00:00 2001 From: Bradley Martin Date: Fri, 8 Dec 2023 11:13:56 +0000 Subject: [PATCH] Add multiple dimensions to the Frohlich model. Fix k-space explicit code. --- src/FeynmanTheory.jl | 37 ++++----- src/HolsteinTheory.jl | 86 +++++++++++---------- src/KSpaceTheory.jl | 162 ++++++++++++++++++++------------------- src/MemoryFunction.jl | 44 +++++------ src/Polaron.jl | 44 ++++++----- src/ResponseFunctions.jl | 19 +++-- 6 files changed, 201 insertions(+), 191 deletions(-) mode change 100755 => 100644 src/HolsteinTheory.jl diff --git a/src/FeynmanTheory.jl b/src/FeynmanTheory.jl index 081267c..295c0ef 100644 --- a/src/FeynmanTheory.jl +++ b/src/FeynmanTheory.jl @@ -34,9 +34,9 @@ Integral of Eqn. (31) in Feynman 1955. Part of the overall ground-state energy e See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. """ - B(v, w, α, ω) = α * ω / √π * quadgk(τ -> phonon_propagator(τ / ω, ω) / sqrt(polaron_propagator(τ, v, w)), 0, Inf)[1] +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) = sum(B.(v, w, α, ω)) +B(v, w, α::Vector, ω::Vector; dims = 3) = sum(B.(v, w, α, ω; dims = dims)) A(v, w, ω) = -3 * (v - w) / 2 * ω @@ -60,7 +60,8 @@ 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, ω, β) = 3 / β / ω * (log(v / w) - 1 / 2 * log(2π * ω * β) - (v - w) * ω * β / 2 - log(1 - exp(-v * ω * β)) + log(1 - exp(-w * ω * β))) * ω +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(ω) @@ -71,7 +72,7 @@ Hellwarth's B expression from Eqn. (62c) in Hellwarth et al. 1999 PRB. Part of t See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299. """ -B(v, w, α, ω, β) = α * ω / √π * quadgk(τ -> phonon_propagator(τ / ω, ω, β) / sqrt(polaron_propagator(τ, v, w, β * ω)), 0, β * ω / 2)[1] +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, α, ω, β)) @@ -82,7 +83,7 @@ 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, ω, β) = 3 / 4 * (v^2 - w^2) * ω / v * (coth(v * β * ω / 2) - 2 / (v * β * ω)) C(v, w, ω::Vector, β) = sum(C.(v, w, ω, β)) / length(ω) @@ -444,15 +445,15 @@ 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, α, ω, β) +function frohlich_energy(v, w, α, ω, β; dims = 3) # Insurance to avoid breaking the integrals with Infinite beta. if β == Inf - return frohlich_energy(v, w, α, ω) + return frohlich_energy(v, w, α, ω; dims = dims) end - Ar = A(v, w, ω, β) + Ar = A(v, w, ω, β) * dims / 3 Br = B(v, w, α, ω, β) - Cr = C(v, w, ω, β) + Cr = C(v, w, ω, β) * dims / 3 # Free energy in units of meV return -(Ar + Br + Cr), Ar, Br, Cr end @@ -470,10 +471,10 @@ Calculates the zero-temperature ground-state energy of the polaron for a materia See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660. """ -function frohlich_energy(v, w, α, ω) - Ar = A(v, w, ω) - Br = B(v, w, α, ω) - Cr = C(v, w, ω) +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 @@ -535,14 +536,10 @@ function feynmanvw(v::Vector, w::Vector, αωβ...; upper_limit=1e6) return Δv .+ w, w, E, A, B, C end -function feynmanvw(v::Real, w::Real, αωβ...; upper_limit=1e6) +function feynmanvw(v::Real, w::Real, αωβ...; upper = [Inf, Inf], lower = [0, 0]) Δv = v .- w - initial = [Δv + eps(Float64), w] - - # Limits of the optimisation. - lower = [0, 0] - upper = [upper_limit, upper_limit] + 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] @@ -568,5 +565,5 @@ function feynmanvw(v::Real, w::Real, αωβ...; upper_limit=1e6) return Δv .+ w, w, E, A, B, C end -feynmanvw(αωβ...; upper_limit=1e6) = feynmanvw(3.4, 2.6, αωβ...; upper_limit=upper_limit) +feynmanvw(αωβ...) = feynmanvw(3.4, 2.6, αωβ...) diff --git a/src/HolsteinTheory.jl b/src/HolsteinTheory.jl old mode 100755 new mode 100644 index f576f63..660e043 --- a/src/HolsteinTheory.jl +++ b/src/HolsteinTheory.jl @@ -145,9 +145,13 @@ 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 * (erf(π * sqrt(propagator)) / 2 / sqrt(π * propagator))^dims + propagator = polaron_propagator(τ, v, w, β * ω) / 2 + phonon_propagator(τ / ω, ω, β) * coupling * P(dims, momentum_cutoff^2 * propagator) / (4π * propagator)^(dims/2) end """ @@ -178,12 +182,12 @@ 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) - coupling = holstein_coupling(1, α, 1; dims = 1) - propagator = polaron_propagator(τ / ω, v * ω, w * ω) / 2 - phonon_propagator(τ / ω, ω) * coupling * (erf(2π * sqrt(propagator * 2 / dims)) / 2 * sqrt(π / propagator) / π√2)^dims - return coupling * phonon_propagator(τ / ω, ω) - + 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) end """ @@ -215,7 +219,10 @@ 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) - integral, _ = quadgk(τ -> holstein_interaction_energy_integrand(τ, v, w, α, ω, β; dims = dims), 0, β/2) + 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 end @@ -254,9 +261,9 @@ 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 = -2 * dims + electron_energy(v, w, ωβ...) + 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 + return (kinetic_energy + potential_energy), kinetic_energy, potential_energy end @@ -291,7 +298,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_bounds = [0, 0], upper_bounds = [Inf, Inf]) +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] @@ -302,26 +309,26 @@ function vw_variation(energy, initial_v, initial_w; lower_bounds = [0, 0], upper # 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_bounds, - upper_bounds, + lower, + upper, initial, Fminbox(BFGS()), - Optim.Options(show_trace = false, g_tol = 1e-12) + # 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, kinetic_energy, potential_energy = energy(Δv + w, 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, kinetic_energy, potential_energy + return Δv + w, w, energy_minimized... end -vw_variation(energy) = vw_variation(energy, 5, 3; lower_bounds = [0, 0], upper_bounds = [Inf, Inf]) +vw_variation(energy) = vw_variation(energy, 5, 3; lower = [0, 0], upper = [1e4, 1e4]) """ general_memory_function(Ω, structure_factor; limits = [0, Inf]) @@ -418,17 +425,14 @@ 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 + momentum_cutoff = gamma(dims / 2 + 1)^(1 / dims) * (2√π) coupling = holstein_coupling(1, α, ω; dims = dims) - - propagator = polaron_propagator(im * t * ω, v, w, β * ω) * 2 / ω - - first_integral = erf(π * sqrt(propagator)) / 4 / √π / propagator^(3/2) - exp(-π^2 * propagator) / propagator / 2 - second_integral = erf(π * sqrt(propagator)) / sqrt(π * propagator) / 2 - - prefactor = 2 * phonon_propagator(im * t, ω, β) - - prefactor * coupling * first_integral * second_integral^(dims - 1) + 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 / ω end """ @@ -460,17 +464,11 @@ println(result) This example calculates the structure factor for the given values of `t`, `v`, `w`, `α`, and `ω`. The result is then printed. """ 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 / ω - - first_integral = erf(π * sqrt(propagator)) / 4 / √π / propagator^(3/2) - exp(-π^2 * propagator) / propagator / 2 - second_integral = erf(π * sqrt(propagator)) / sqrt(π * propagator) / 2 - - prefactor = 2 * phonon_propagator(im * t, ω) - - prefactor * coupling * first_integral * second_integral^(dims - 1) + 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 end @@ -494,8 +492,11 @@ 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) - structure_factor(t) = holstein_structure_factor(t, v, w, α, ω, β; dims = dims) - return general_memory_function(Ω, structure_factor) + 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) end """ @@ -519,7 +520,7 @@ In this example, the `holstein_memory_function` is called with the parameters ` """ 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, 1e5]) + return general_memory_function(Ω / ω, structure_factor, limits = [0, 1e6]) end """ @@ -553,8 +554,11 @@ println(result) This code calculates the mobility using the given parameters and prints the result. """ function holstein_mobility(v, w, α, ω, β; dims = 3) + if β == Inf + return Inf + end structure_factor(t) = holstein_structure_factor(t, v, w, α, ω, β; dims = dims) - abs(1 / imag(general_memory_function(structure_factor, limits = [0, 1e5]))) + abs(1 / imag(general_memory_function(structure_factor))) end function holstein_complex_impedence(Ω, v, w, α, ωβ...; dims = 3) diff --git a/src/KSpaceTheory.jl b/src/KSpaceTheory.jl index 1f55a6d..8a1fbcc 100644 --- a/src/KSpaceTheory.jl +++ b/src/KSpaceTheory.jl @@ -21,8 +21,8 @@ println(result) # print the result Expected Output: 0.5 """ -function cartesian_k_integrand(k, coupling, propagator; rₚ = 1) - coupling(k) * exp(-k^2 * rₚ^2 * propagator / 2) +function cartesian_k_integrand(k, coupling, propagator) + coupling(k) * exp(-k^2 * propagator / 2) end @@ -50,8 +50,8 @@ println(result) # print the result Expected Output: `0.5` """ -function spherical_k_integrand(k, coupling, propagator; rₚ = 1) - 4π * k^2 * coupling(k) * exp(-k^2 * rₚ^2 * propagator / 2) +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 @@ -80,8 +80,8 @@ 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; rₚ = 1, a = 1, limits = [-π, π]) - integral, _ = quadgk(k -> cartesian_k_integrand(k, coupling, propagator; rₚ = rₚ), limits[1] / a, limits[2] / a) +function cartesian_k_integral(coupling, propagator; limits = [-π, π]) + integral, _ = quadgk(k -> cartesian_k_integrand(k, coupling, propagator), limits[1], limits[2]) integral * a / 2π end @@ -110,9 +110,9 @@ 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; rₚ = 1, limits = [0, Inf]) - integral, _ = quadgk(k -> spherical_k_integrand(k, coupling, propagator; rₚ = rₚ), limits[1], limits[2]) - integral / 8π^3 +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 @@ -164,8 +164,9 @@ println(result) Expected Output: `6.0` """ -function frohlich_coupling(k, α, ω) - ω^(3/2) * 2√2 * π * α / k^2 +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 """ @@ -188,10 +189,11 @@ Calculate the integrand for the Holstein interaction energy in k-space at finite # 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, rₚ = 1, a = 1, limits = [-π, π]) +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(τ, ω, β) * cartesian_k_integral(coupling, propagator; rₚ = rₚ, a = a, limits = limits)^dims + propagator = polaron_propagator(τ, v, w, β) + phonon_propagator(τ, ω, β) * spherical_k_integral(coupling, propagator; dims = dims, limits = [0, momentum_cutoff]) end @@ -214,10 +216,11 @@ Calculate the integrand for the Holstein interaction energy in k-space at zero t # 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, rₚ = 1, a = 1, limits = [-π, π]) - coupling(k) = holstein_coupling(k, α, ω; dims = dims) - propagator = polaron_propagator(τ, v, w) - phonon_propagator(τ, ω) * carteisan_k_integral(coupling, propagator; rₚ = rₚ, a = a, limits = limits)^dims +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 @@ -239,10 +242,10 @@ Calculate the integrand for the Frohlich interaction energy in k-space at finite ## 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, α, ω, β; rₚ = 1, limits = [0, Inf]) - coupling(k) = frohlich_coupling(k, α, ω) +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; rₚ = rₚ, limits = limits) + phonon_propagator(τ, ω, β) * spherical_k_integral(coupling, propagator; limits = limits, dims = dims) end """ @@ -262,10 +265,10 @@ Calculate the integrand for the Frohlich interaction energy in k-space at finite ## 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, α, ω; rₚ = 1, limits = [0, Inf]) - coupling(k) = frohlich_coupling(k, α, ω) - propagator = polaron_propagator(τ, v, w, ω) - phonon_propagator(τ, ω) * spherical_k_integral(coupling, propagator; rₚ = rₚ, limits = limits) +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 @@ -289,9 +292,9 @@ Calculate the Holstein polaron interaction energy in k-space at finite temperaur ## 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, rₚ = 1, a = 1, limits = [-π, π]) - interation_energy, _ = quadgk(τ -> holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; dims = dims, rₚ = rₚ, a = a, limits = limits), 0, β/2) - return interation_energy +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 @@ -314,9 +317,9 @@ Calculate the Holstein polaron interaction energy in k-space at zero temperaure. ## 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, rₚ = 1, a = 1, limits = [-π, π]) - interation_energy, _ = quadgk(τ -> holstein_interaction_energy_integrand_k_space(τ, v, w, α, ω; dims = dims, rₚ = rₚ, a = a, limits = limits), 0, Inf) - return interation_energy +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 """ @@ -338,9 +341,9 @@ Calculate the Frohlich polaron interaction energy in k-space at finite temperaur ## Returns A scalar value representing the Frohlich polaron interaction energy in k-space at finite temperature. """ -function frohlich_interaction_energy_k_space(v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) - interation_energy, _ = quadgk(τ -> frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω, β; rₚ = rₚ, limits = limits), 0, β/2) - return interation_energy +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 """ @@ -361,9 +364,9 @@ Calculate the Frohlich polaron interaction energy in k-space at zero temperaure. ## Returns A scalar value representing the Frohlich polaron interaction energy in k-space at zero temperature. """ -function frohlich_interaction_energy_k_space(v, w, α, ω; rₚ = 1, limits = [0, Inf]) - interation_energy, _ = quadgk(τ -> frohlich_interaction_energy_integrand_k_space(τ, v, w, α, ω; rₚ = rₚ, limits = limits), 0, Inf) - return interation_energy +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 @@ -393,8 +396,11 @@ println(result) Expected Output: A scalar value representing the calculated free electron energy at finite temperature. """ -function electron_energy(v, w, ω, β) - -(A(v, w, ω, β) + C(v, w, ω, β)) / 3 +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 @@ -423,8 +429,8 @@ println(result) Expected Output: A scalar value representing the calculated free electron energy at finite temperature. """ -function electron_energy(v, w, ω) - (v - w)^2 / (4 * v) * ω +function electron_energy(v, w, ω; dims = 3) + (v - w)^2 / (4 * v) * ω * dims end @@ -449,9 +455,9 @@ Calculate the total energy, kinetic energy, and interaction energy of the Holste - `kinetic_energy`: The calculated polaron kinetic energy. - `interaction_energy`: The calculated polaron interaction energy. """ -function holstein_energy_k_space(v, w, α, ωβ...; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - kinetic_energy = -2 * dims + electron_energy(v, w, ωβ...) - interaction_energy = -holstein_interaction_energy_k_space(v, w, α, ωβ...; dims = dims, rₚ = rₚ, a = a, limits = limits) +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 @@ -475,9 +481,9 @@ Calculate the total energy, kinetic energy, and interaction energy of the Frohli - `kinetic_energy`: The calculated polaron kinetic energy. - `interaction_energy`: The calculated polaron interaction energy. """ -function frohlich_energy_k_space(v, w, α, ωβ...; rₚ = 1, limits = [0, Inf]) - kinetic_energy = electron_energy(v, w, ωβ...) - interaction_energy = -frohlich_interaction_energy_k_space(v, w, α, ωβ...; rₚ = rₚ, limits = limits) +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 @@ -502,19 +508,19 @@ Calculate the structure factor in k-space for the Holstein lattice polaron model ## 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, rₚ = 1, a = 1, limits = [-π, π]) +function holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = 3) - coupling_one(k) = holstein_coupling(k, α, ω; dims = dims) * k^2 - coupling_two(k) = holstein_coupling(k, α, ω; dims = dims) + 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, β) - first_integral = cartesian_k_integral(coupling_one, propagator; rₚ = rₚ, a = a, limits = limits) - second_integral = cartesian_k_integral(coupling_two, propagator; rₚ = rₚ, a = a, limits = limits) + integral = spherical_k_integral(coupling_one, propagator; limits = [0, momentum_cutoff]) prefactor = 2 / 3 * phonon_propagator(im * t, ω, β) - prefactor * dims * first_integral * second_integral^(dims - 1) + prefactor * integral end @@ -537,19 +543,19 @@ Calculate the structure factor in k-space for the Holstein lattice polaron model ## 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, rₚ = 1, a = 1, limits = [-π, π]) +function holstein_structure_factor_k_space(t, v, w, α, ω; dims = 3) + + momentum_cutoff = gamma(dims / 2 + 1)^(1/ dims) * (2√π) - coupling_one(k) = holstein_coupling(k, α, ω; dims = dims) * k^2 - coupling_two(k) = holstein_coupling(k, α, ω; dims = dims) + coupling(k) = holstein_coupling(k, α, ω; dims = dims) * k^2 propagator = polaron_propagator(im * t, v, w) - first_integral = cartesian_k_integral(coupling_one, propagator; rₚ = rₚ, a = a, limits = limits) - second_integral = cartesian_k_integral(coupling_two, propagator; rₚ = rₚ, a = a, limits = limits) + integral = spherical_k_integral(coupling, propagator; limits = [0, momentum_cutoff]) prefactor = 2 / 3 * phonon_propagator(im * t, ω) - prefactor * dims * first_integral * second_integral^(dims - 1) + prefactor * integral end @@ -571,13 +577,13 @@ Calculate the structure factor in k-space for the Frohlich continuum polaron mod ## 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, α, ω, β; rₚ = 1, limits = [0, Inf]) +function frohlich_structure_factor_k_space(t, v, w, α, ω, β; limits = [0, Inf], dims = 3) - coupling(k) = frohlich_coupling(k, α, ω) * k^2 + coupling(k) = frohlich_coupling(k, α, ω; dims = dims) * k^2 propagator = polaron_propagator(im * t, v, w, β) - integral = spherical_k_integral(coupling, propagator; rₚ = rₚ, limits = limits) + integral = spherical_k_integral(coupling, propagator; limits = limits, dims = dims) prefactor = 2 / 3 * phonon_propagator(im * t, ω, β) @@ -602,13 +608,13 @@ Calculate the structure factor in k-space for the Frohlich continuum polaron mod ## 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, α, ω; rₚ = 1, limits = [0, Inf]) +function frohlich_structure_factor_k_space(t, v, w, α, ω; limits = [0, Inf], dims = 3) - coupling(k) = frohlich_coupling(k, α, ω) * k^2 + coupling(k) = frohlich_coupling(k, α, ω; dims = dims) * k^2 propagator = polaron_propagator(im * t, v, w) - integral = spherical_k_integral(coupling, propagator; rₚ = rₚ, limits = limits) + integral = spherical_k_integral(coupling, propagator; limits = limits, dims = dims) prefactor = 2 / 3 * phonon_propagator(im * t, ω) @@ -636,8 +642,8 @@ Calculate the memory function for the Holstein model in k-space at finite temper ## 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, rₚ = 1, a = 1, limits = [-π, π]) - structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = dims, rₚ = rₚ, a = a, limits = limits) +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 @@ -661,8 +667,8 @@ Calculate the memory function for the Holstein model in k-space at zero temperat ## 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, rₚ = 1, a = 1, limits = [-π, π]) - structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ω; dims = dims, rₚ = rₚ, a = a, limits = limits) +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 @@ -687,8 +693,8 @@ Calculate the DC mobility in k-space for a Holstein polaron system at finite tem ## Returns The DC mobility in k-space for the Holstein polaron system at finite temperature. """ -function holstein_mobility_k_space(v, w, α, ω, β; dims = 3, rₚ = 1, a = 1, limits = [-π, π]) - structure_factor(t) = holstein_structure_factor_k_space(t, v, w, α, ω, β; dims = dims, rₚ = rₚ, a = a, limits = limits) +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 @@ -711,8 +717,8 @@ Calculate the memory function for the Frohlich model in k-space at finite temper ## 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, α, ω, β; rₚ = 1, limits = [0, Inf]) - structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ω, β; rₚ = rₚ, limits = limits) +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 @@ -734,8 +740,8 @@ Calculate the memory function for the Frohlich model in k-space at zero temperat ## 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, α, ω; rₚ = 1, limits = [0, Inf]) - structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ω; rₚ = rₚ, limits = limits) +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 @@ -758,8 +764,8 @@ Calculate the DC mobility in k-space for a Frohlich polaron system at finite tem ## Returns The DC mobility in k-space for the Frohlich polaron system at finite temperature. """ -function frohlich_mobility_k_space(v, w, α, ω, β; rₚ = 1, limits = [0, Inf]) - structure_factor(t) = frohlich_structure_factor_k_space(t, v, w, α, ω, β; rₚ = rₚ, limits = limits) - 1 / imag(general_memory_function(structure_factor; limits = [0, Inf])) +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/MemoryFunction.jl b/src/MemoryFunction.jl index 4efddbc..3d5835a 100644 --- a/src/MemoryFunction.jl +++ b/src/MemoryFunction.jl @@ -4,36 +4,34 @@ # [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, α, ω, β) - - coupling = 2 * α / 3 / √π * ω^2 - - propagator = polaron_propagator(im * t, v, w, β * ω) - - integral = phonon_propagator(im * t / ω, ω, β) / propagator^(3/2) - - coupling * integral +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, α, ω) - - coupling = 2 * α / 3 / √π * ω^2 - - propagator = polaron_propagator(im * t, v, w) - - integral = phonon_propagator(im * t / ω, ω) / propagator^(3/2) - - coupling * integral +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, α, ω, β) - structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω, β) +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, α, ω) - structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω) - return general_memory_function(Ω / ω, structure_factor) +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, α, ωβ...)) diff --git a/src/Polaron.jl b/src/Polaron.jl index 27765a5..6020576 100644 --- a/src/Polaron.jl +++ b/src/Polaron.jl @@ -13,6 +13,7 @@ mutable struct Polaron ωeff # Effective phonon mode frequency. β # Reduced thermodynamic beta ħω₀/kBT. Ω # Electric field frequency. + d # Number of spatial dimensions. v0 # Athermal variational parameter v. w0 # Athermal variational parameter w. F0 # Polaron athermal energy (enthalpy). @@ -68,7 +69,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.11, w_guesses=2.87, verbose=false) +function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.01, w_guesses=2.99, dims=3,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. @@ -94,6 +95,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses "ωeff" => ωeff, # Effective Phonon frequency. "β" => 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. @@ -162,18 +164,12 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses # 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, " | ω = ", ω) - else - println(io, "\e[KPhonon frequencies | ωeff = ", ωeff, " | ω = ", join(round.(first(ω, 2), digits=2), ", ")..., " ... ", join(round.(last(ω, 2), digits=2), ", ")...) - end - end - - # Print α coupling parameters. - if verbose - if num_ω == 1 println(io, "\e[KFröhlich coupling | αeff = ", αeff, " | α = ", join(α, ", ")...) else + 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) end # Small alpha (α → 0) approximate energy. @@ -242,7 +238,8 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses # 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. - v_gs, w_gs, F_gs, A_gs, B_gs, C_gs = feynmanvw(v_guesses, w_guesses, α, ω) + athermal_energy(v, w) = frohlich_energy(v, w, α, ω; dims = dims) + 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 @@ -338,7 +335,8 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses # 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. - v, w, F, A, B, C = feynmanvw(v_guesses, w_guesses, α, ω, β) + thermal_energy(v, w) = frohlich_energy(v, w, α, ω, β; dims = dims) + v, w, F, A, B, C = vw_variation(thermal_energy, v_guesses, w_guesses) # Update the guesses to keep them close-ish to the true solutions during loops over temperatures. v_guesses, w_guesses = v, w @@ -413,7 +411,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses end # Calculate and store the DC mobiliies. - μ = frohlich_mobility(v, w, α, ω, β) / mb + μ = frohlich_mobility(v, w, α, ω, β; dims = dims) / mb p["μ"][i, j] = μ # Print DC mobilities. @@ -479,7 +477,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses end # Calculate and store polaron memory functions (akin to self energy). - χ = frohlich_memory_function(Ω, v, w, α, ω, β) + χ = frohlich_memory_function(Ω, v, w, α, ω, β; dims = dims) p["χ"][k, i, j] = χ # Print memory function. @@ -511,9 +509,9 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses process += 1 # Increment total iterator. end end - if verbose print(io, "\e[26F")end # Move up 26 lines and erase. + if verbose print(io, "\e[26F") end # Move up 26 lines and erase. end - if verbose print(io, "\e[26F") end # Move up 26 lines and erase. + if verbose print(io, "\e[27F") end # Move up 26 lines and erase. end if verbose print("\e[?25h") end # Show cursor again. @@ -526,6 +524,7 @@ function polaron(αrange, Trange, Ωrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses p["ωeff"], # Hellwarth eff phonon frequency. p["β"], # Betas. p["Ω"], # Photon frequencies. + p["d"], # Number of spatial dimensions. p["v0"], # Athermal v params. p["w0"], # Athermal w params. p["F0"], # Athermal energies. @@ -574,22 +573,22 @@ 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, verbose=false) = polaron([α], Trange, Ωrange; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, verbose=verbose) +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) """ No frequency input. """ -polaron(αrange, Trange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, verbose=false) = polaron(αrange, Trange, 1; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, verbose=verbose) +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) """ No temperature input => 300 K. """ -polaron(αrange; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, verbose=false) = polaron(αrange, 300, 1; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, verbose=verbose) +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) """ No input => α = 1 """ -polaron(; ω=1, ωeff=1, mb=1, β0=1, v_guesses=3.11, w_guesses=2.87, verbose=false) = polaron(1, 300, 1; ω=ω, ωeff=ωeff, mb=mb, β0=β0, v_guesses=v_guesses, w_guesses=w_guesses, verbose=verbose) +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) """ polaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, verbose=false) @@ -597,7 +596,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, verbose=false) +function polaron(material::Material, TΩrange...; v_guesses=3.11, w_guesses=2.87, dims=3, verbose=false) # Show material data. if verbose @@ -610,7 +609,7 @@ 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, verbose=verbose) + 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) # Return material-specific, unitful Polaron type. return p @@ -627,6 +626,7 @@ function Base.show(io::IO, ::MIME"text/plain", x::Polaron) println(io_lim, "\e[KPhonon frequencies | ωeff = ", x.ωeff, " | ω = ", x.ω) println(io_lim, "\e[KFröhlich coupling | αeff = ", x.αeff, " | α = ", x.α) + println(io_lim, "\e[KNumber of spatial dimensions | d = ", x.d) println(io_lim, "\e[KSmall α→0 energy | Fs = ", x.Fs) println(io_lim, "\e[KLarge α→∞ energy | Fl = ", x.Fl) println(io_lim, "\e[KSmall α→0 fictitious mass | Ms = ", x.Ms) @@ -712,6 +712,7 @@ function save_polaron(polaron::Polaron, prefix) "phonon freq eff", polaron.ωeff, "beta", polaron.β, "E-field freq", polaron.Ω, + "dims", polaron.d, "v0", polaron.v0, "w0", polaron.w0, "athermal energy", polaron.F0, @@ -777,6 +778,7 @@ function load_polaron(polaron_file_path) data["phonon freq eff"], data["beta"], data["E-field freq"], + data["dims"], data["v0"], data["w0"], data["athermal energy"], diff --git a/src/ResponseFunctions.jl b/src/ResponseFunctions.jl index 86eb47c..2ad5493 100644 --- a/src/ResponseFunctions.jl +++ b/src/ResponseFunctions.jl @@ -54,8 +54,8 @@ Calculate the complex impedence Z(Ω) of the polaron at finite temperatures for See FHIP 1962: https://doi.org/10.1103/PhysRev.127.1004. """ -function frohlich_complex_impedence(Ω, v, w, α, ωβ...) - -im * (Ω + frohlich_memory_function(Ω, v, w, α, ωβ...)) +function frohlich_complex_impedence(Ω, v, w, α, ωβ...; dims = 3) + -im * (Ω + frohlich_memory_function(Ω, v, w, α, ωβ...; dims = dims)) end """ @@ -73,8 +73,8 @@ Calculate the complex conductivity σ(Ω) of the polaron at finite temperatures See also [`polaron_complex_impedence`](@ref) """ -function frohlich_complex_conductivity(Ω, v, w, α, ωβ...) - return 1 / frohlich_complex_impedence(Ω, v, w, α, ωβ...) +function frohlich_complex_conductivity(Ω, v, w, α, ωβ...; dims = 3) + return 1 / frohlich_complex_impedence(Ω, v, w, α, ωβ...; dims = dims) end """ @@ -94,9 +94,12 @@ 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, α, ω, β) - structure_factor(t) = frohlich_structure_factor(t, v, w, α, ω, β) - return abs(imag(general_memory_function(structure_factor))) / ω +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 """ @@ -113,7 +116,7 @@ The polaron mobility. See also [`inverse_polaron_mobility`](@ref) """ -frohlich_mobility(v, w, α, ω, β) = 1 ./ inverse_frohlich_mobility(v, w, α, ω, β) +frohlich_mobility(v, w, α, ω, β; dims = 3) = 1 ./ inverse_frohlich_mobility(v, w, α, ω, β; dims = dims) """ inverse_FHIP_mobility_lowT(v, w, α, ω, β)