<commit_before>### A Pluto.jl notebook ###
# v0.19.4

using Markdown
using InteractiveUtils

# ╔═╡ 66d2d0e0-ee56-11ec-21f9-5d6ae25cff48
begin
	using Optim
	using QuadGK
	using SpecialFunctions
	using Plots
end

# ╔═╡ 9f6474e5-41f1-4b96-aeee-d299404a3718
md"""
## Feynman's variational theory for the (continious) Frohlich model.
"""

# ╔═╡ 9429ae6b-55e1-4fb2-9a88-ae7636c35bb5
begin
"""
    fF(τ, v, w)
Integrand of Eqn. (31) in Feynman 1955. Part of the overall ground-state energy expression.
See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660.
"""
fF(τ, v, w) = (abs(w^2 * τ + (v^2 - w^2) / v * (1 - exp(-v*τ))))^(-0.5) * exp(-τ)


"""
    AF(v, w, α)
Integral of Eqn. (31) in Feynman 1955. Part of the overall ground-state energy expression.
See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660.
"""
AF(v, w, α) = π^(-0.5) * α * v * quadgk(τ -> fF(τ, v, w), 0.0, Inf)[1]

"""
    F(τ, v, w)
Ground state energy expression. Eqn. (33) in Feynman 1955.
See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660.
"""
F(v, w, α) = (3 / (4 * v)) * (v - w)^2 - AF(v, w, α)

# Let's wrap the Feynman athermal variation approximation in a simple function

"""
    feynmanvw(α; v = 0.0, w = 0.0)
Calculate v and w variational polaron parameters, for the supplied α Frohlich coupling. Returns (v, w).
This version uses the original athermal action (Feynman 1955).
"""
function feynmanvw_old(α; v = 3.0, w = 3.0) # v, w defaults

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = F(x[1] + x[2], x[2], α)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end
end

# ╔═╡ 1cf1e5cc-301c-4779-86cb-5b6436ecebf0
v, w = feynmanvw_old(2.39)

# ╔═╡ 2a609c87-4579-490e-99e7-0d0c25e11333
F(v, w, 2.39)

# ╔═╡ 7bf82f19-e6fc-4c3f-b7e8-a02bebf0d0dd
md"""
# Explicit k-space integral version of Frohlich model.
"""

# ╔═╡ 6cd58deb-fb04-44fc-a854-dce5b5cb8565
function D(t, v, w)
	abs(w^2 / v^2 * t + (v^2 - w^2) / v^3 * (1 - exp(-v * t)))
end

# ╔═╡ f7a07dc5-ba45-4870-b94c-44ae5d85436d
md"""
### Test first that it reproduces the Feynman analytic result. (R = 0).
"""

# ╔═╡ cc237871-08fd-4658-8fa5-8658d9394b3b
md"""
### Now try Frohlich model with k-space integral.
"""

# ╔═╡ 907b2475-0749-4d85-b09e-6eed75d35e68
md"""
### Now the k-space integral for the Holstein can be done and is expressed in terms of an erf function. So lets do it as it will be quicker!
"""

# ╔═╡ e7a4f46e-8fcc-48ef-9c5c-5790b42d7648
"""
Let's now plot the α dependence of Frohlich vs Holstein!
"""

# ╔═╡ 78936a9d-9567-4def-bf89-85be9e62e12d
α = 0.1:0.1:12

# ╔═╡ 8d726e23-5a40-4a10-8405-7c9775a0c9d0
frohlich_params = feynmanvw_old.(α)

# ╔═╡ 2cdfdcbc-a6e7-444c-a9a0-47529ddb091a
begin
	vf = [i[1] for i in frohlich_params]
	wf = [i[2] for i in frohlich_params]
end

# ╔═╡ 35ce9514-521c-4985-a85e-f236f67417c9
frohlich_energies = [F(vf[i], wf[i], α[i]) for i in 1:length(α)]

# ╔═╡ 3f5069df-89e2-4ac7-96e4-f3c79936fdd9
md"""
## Thermal (free energy) version of the theory...
"""

# ╔═╡ 9c1284f7-f503-4c4d-be6a-e4e04043246a
"""
    A(v, w, β)

Hellwarth's A expression from Eqn. (62b) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
A(v, w, β) = 3 / β * (log(v / w) - 1/2 * log(2π * β) - log(sinh(v * β / 2) / sinh(w * β / 2)))

# ╔═╡ b36f7997-f29c-4a92-bb00-9e82ae7d2a05
"""
    Y(x, v, β)

Hellwarth's Y expression from Eqn. (62d) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression. Contained in denominator of the integrand of Eqn. (62c).

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
Y(x, v, β) = 1 / (1 - exp(-v * β)) * (1 + exp(-v * β) - exp(-v * x) - exp(v * (x - β)))

# ╔═╡ d4e1ac9e-84d7-4809-ac1c-c24a9263db84
"""
    f(x, v, w, β)

Integrand of Hellwarth's B expression from Eqn. (62c) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
f_Fro(x, v, w, β) = (exp(β - x) + exp(x)) / sqrt(abs(w^2 * x * (1 - x / β) + Y(x, v, β) * (v^2 - w^2) / v))


# ╔═╡ ca1d9705-b49e-43f8-bd17-ab4a7e28b578
D(x, v, w, β) = w^2 / v^2 * x * (1 - x / β) + Y(x, v, β) * (v^2 - w^2) / v^3

# ╔═╡ 00fc8b3d-c485-43e2-92f0-451539ecf1e6
"""
 k_integral

R is cutoff for Holstei n k-space integral (in k-space spherical volume, R)

if R==0, Frohlich model (power should be 0, as the K^2 cancels with the 1/K^2 dependency in the el-ph matrix element)

Power is the K^Power
"""
function k_integral(t, v, w; power = 0, R = 0)
	if R == 0 # Frohlich
		return quadgk(k -> k^power * exp(-k^2 * D(t, v, w) / 2), -Inf, Inf)[1] / sqrt(2pi)
	else # Holstein, note limits on integral, R=k-space volume
		return quadgk(k -> k^power * exp(-k^2 * D(t, v, w) / 2), 0, R)[1] / sqrt(2π)
	end
end

# ╔═╡ dcf7f24a-6ef4-4613-9778-e2637f5e6c16
function Frohlich_el_integral(α, v, w; power = 0, R = 0)
	α * quadgk(t -> exp(-t) * k_integral(t, v, w; power = power, R = R), 0.0, Inf, rtol = 1e-4)[1] / sqrt(pi)
end

# ╔═╡ 6714b04e-9b26-42d0-89eb-f6c510cd5f56
function Frohlich_energy(α, v, w; power = 0, R = 0)
	3 * (v - w)^2 / (4 * v) - Frohlich_el_integral(α, v, w; power = power, R = R)
end

# ╔═╡ 429d0098-9d78-47e5-915b-cf16c7fce660
function var_params_k(α; v = 3.0, w = 3.0, power = 0, R = 0)
	 # Limits of the optimisation.
    lower = [0.0, 0.3]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = Frohlich_energy(α, x[1] + x[2], x[2]; power = power, R = R)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 71118d82-37fb-4504-b944-c3461de91233
v_k, w_k = var_params_k(2.39; v = v, w = w)

# ╔═╡ 80cd6d37-92d3-4757-a62e-490cffb4b63c
vh_k, wh_k = var_params_k(2.39; v = v, w = w, power = 2, R = (6 * π^2)^(1/3))

# ╔═╡ 2be860db-52b0-4514-872f-db79dc454323
Frohlich_energy(2.39, v_k, w_k)

# ╔═╡ 410d234c-25fe-4752-926e-f6c5d263ef05
Frohlich_energy(2.39, vh_k, wh_k; power = 2, R = (6 * π^2)^(1/3))

# ╔═╡ 079d7328-c2e1-4681-a4e3-e19a53841d25
function holstein_elph_integral(α, v, w; R = (6 * π^2)^(1/3))
	d(t) = D(t, v, w)
	integrand(t) = erf(sqrt(d(t) / 2) * R) / (2 * (d(t))^(3/2)) - R * exp(-R^2 * d(t) / (2)) / (sqrt(2π) * d(t))
	integral = α * quadgk(t -> exp(-t) * integrand(t), 0, Inf, rtol = 1e-4)[1] / sqrt(π)
	return integral
end

# ╔═╡ bc20477e-2ce3-4377-9af0-5d479256175a
function holstein_energy(α, v, w)
	3 * (v - w)^2 / (4 * v) - holstein_elph_integral(α, v, w)
end

# ╔═╡ 380c5c1b-01e1-455b-a757-46c87fc45285
function holstein_variational_solution(α; v = 3.0, w = 3.0)
	 # Limits of the optimisation.
    lower = [0.0, 0.3]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = holstein_energy(α, x[1] + x[2], x[2])

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 650e36ab-7ac8-4d0b-ad74-f5693c2c1586
holstein_energies = [holstein_energy(α[i], vf[i], wf[i]) for i in 1:length(α)]

# ╔═╡ b2454e06-146c-4fc1-8928-70c7982922ec
begin
	pe = plot(α, frohlich_energies .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Fro", linewidth = 2, xlabel = "alpha, α", ylabel = "Energy (meV)")
	plot!(α, holstein_energies .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Hol", linewidth = 2, linestyle = :dash)
end

# ╔═╡ 67154f5a-e61e-4a5e-abf8-215bf8c8646e
savefig(pe, "energy_comparison.pdf")

# ╔═╡ 4e99725c-0895-4bc3-97ba-70eb675c7f47
f_Hol(x, v, w, β) = (exp(β - x) + exp(x)) * (erf(sqrt(D(x, v, w, β) / 2) * (6 * π^2)^(1/3)) / (2 * D(x, v, w, β)^(3/2)) - (6 * π^2)^(1/3) * exp(-D(x, v, w, β) * (6 * π^2)^(2/3) / 2) / (sqrt(2π) * D(x, v, w, β)))

# ╔═╡ e8b8e4de-e548-447e-82c5-979124ce6cf3
"""
    B(v, w, β, α)

Hellwarth's B expression from Eqn. (62c) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
B_Fro(v, w, β, α) = α * v / (sqrt(π) * (exp(β) - 1)) * quadgk(x -> f_Fro(x, v, w, β), 0, β/2, rtol = 1e-5)[1]


# ╔═╡ f87703c9-fa55-420a-b3cc-c6b947ef456a
B_Hol(v, w, β, α) = α / (sqrt(π) * (exp(β) - 1)) * quadgk(x -> f_Hol(x, v, w, β), 0, β/2)[1]

# ╔═╡ c1358b12-4f68-4485-8ee7-73ef672f3ed8
"""
    C(v, w, β)

Hellwarth's C expression from Eqn. (62e) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
C(v, w, β) = 3/4 * (v^2 - w^2) / v * (coth(v * β / 2) - 2 / (v * β))

# ╔═╡ fc699405-e1d6-48dd-aff3-a7ae7fb0114b
"""
    F(v, w, β, α)

Hellwarth's total free energy expression from Eqn. (62a) in Hellwarth et al. 1999 PRB.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
F_Fro(v, w, β, α) = -(A(v, w, β) + B_Fro(v, w, β, α) + C(v, w, β))

# ╔═╡ 9ed1bad6-9ac5-444f-8d02-974bdd955b09
F_Hol(v, w, β, α) = -(A(v, w, β) + B_Hol(v, w, β, α) + C(v, w, β))

# ╔═╡ 5d5bf742-89c2-45e6-81f6-3b22e90b7bd0
"""
Check that it roughly gives the Holstein ground-state energy (high beta)
"""

# ╔═╡ 2455abad-abf2-4b0a-bb3f-30d97de9a9ec
function feynmanvw(α, β; v = 3.0, w = 3.0) # v, w defaults

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = F_Fro(x[1] + x[2], x[2], β, α)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 4d5e8c8e-2225-43fd-a934-0c8dd9f21585
function holstein_variational_solution(α, β; v = 3.0, w = 3.0) # v, w defaults

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = F_Hol(x[1] + x[2], x[2], β, α)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 466bc0e9-4cd2-42f5-8bbc-1317b5c1c7b6
vh, wh = holstein_variational_solution(2.39; v = v, w = w)

# ╔═╡ a1d89e74-7ab6-4c8f-8748-63ede869decb
holstein_energy(2.39, vh, wh)

# ╔═╡ 60254d3a-97b0-4c7e-a6c3-ad9ea0d93c97
F_Hol(vh, wh, 330, 2.39)

# ╔═╡ 5f5785c4-b9be-4246-820d-3e3ccf841cc6
holstein_params = holstein_variational_solution.(α)

# ╔═╡ 2ce6cc8d-0582-4f08-9e2f-5317f2f5b2e0
begin
	vh_new = [i[1] for i in holstein_params]
	wh_new = [i[2] for i in holstein_params]
end

# ╔═╡ f5e2a904-5e2f-4dcd-8691-b7620d24a7f0
begin
	pv = plot(α, vf, label = "Fro v", legend = :topleft, linewidth = 2, xlabel = "alpha, α", ylabel = "v and w", minorgrid = true)
	plot!(α, wf, label = "Fro w", linewidth = 2)
	plot!(α, vh_new, label = "Hol v", linewidth = 2, linestyle = :dash)
	plot!(α, wh_new, label = "Hol w", linewidth = 2, linestyle = :dash)
end

# ╔═╡ a8852b06-3f60-4bac-a610-8b0a1b8cf560
savefig(pv, "vw_comparison.pdf")

# ╔═╡ af11d3d0-b93d-4916-baad-5c81cec45ab4
"""
Temperature dependence free-energy and v/w comparisons
"""

# ╔═╡ 28f3a13d-cee1-4415-be15-4cefc474207d
T_range = 1:400

# ╔═╡ 7cd14da3-799d-4a53-a446-10a6842e5214
hol_params = [holstein_variational_solution(2.39, 2π * 2.25 * 1e12 * 1.05e-34 / T / 1.38e-23) for T in T_range]

# ╔═╡ 2baddbea-169a-4558-9f47-1c43e1269f1f
fro_params = [feynmanvw(2.39, 2π * 2.25 * 1e12 * 1.05e-34 / T / 1.38e-23) for T in T_range]

# ╔═╡ 1d45213c-3921-4a67-a76c-f37db47a14ad
begin
	vh_β = [i[1] for i in hol_params]
	wh_β = [i[2] for i in hol_params]
end

# ╔═╡ 3d423a59-b295-4026-8666-e08c38f01747
begin
	vf_β = [i[1] for i in fro_params]
	wf_β = [i[2] for i in fro_params]
end

# ╔═╡ 0e771039-19b8-4cec-bb12-8ddcb4b6975e
hol_energy_β = [F_Hol(vh_β[T], wh_β[T], 2π * 2.25 * 1e12 * 1.05e-34 / T_range[T] / 1.38e-23, 2.39) for T in 1:length(T_range)]

# ╔═╡ 07e267b3-03e8-4b1f-af0a-6c5ea6c57684
fro_energy_β = [F_Fro(vf_β[T], wf_β[T], 2π * 2.25 * 1e12 * 1.05e-34 / T_range[T] / 1.38e-23, 2.39) for T in 1:length(T_range)]

# ╔═╡ 1d08e7de-026b-4b1c-ba29-d12ce9f2b4ef
begin
	pvt = plot(T_range, vf_β, label = "Fro v", linewidth = 2, legend = :topleft, xlabel = "Temperature (K)", ylabel = "v and w params")
	plot!(T_range, wf_β, label = "Fro w", linewidth = 2)
	plot!(T_range, vh_β, label = "Hol v", linestyle = :dash, linewidth = 2)
	plot!(T_range, wh_β, label = "Hol w", linestyle = :dash, linewidth = 2)
end

# ╔═╡ a6ff82d8-ac2f-4034-ad64-c3800d67ecfa
begin
	pft = plot(T_range, fro_energy_β .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Fro", xlabel = "Temperature (K)", ylabel = "Free energy (ħω)", linewidth = 2)
	plot!(T_range, hol_energy_β .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Hol", linewidth = 2, linestyle = :dash)
end

# ╔═╡ a0cd088c-0412-4b1d-ba32-441f59e20c3c
savefig(pvt, "vw_comparison_temp.pdf")

# ╔═╡ 35864bca-4813-45ab-b579-ba9b40492189
savefig(pft, "free_energy_comparison.pdf")

# ╔═╡ a835bb67-f4d0-4b42-8cda-de6b44936ca2
"""
Have a look at effective mass
"""

# ╔═╡ 5ac9e949-4b55-42f5-9dd4-1de517051112
function fro_eff_mass(v, w, α)
	integrand(τ) = τ^2 * exp(-τ) / (D(τ, v, w))^(3/2)
	integral = α * quadgk(τ -> integrand(τ), 0.0, Inf)[1] / (3 * sqrt(π))
	return 1 + integral
end

# ╔═╡ 6a848cac-dc22-4419-8833-7279a4940834
function hol_eff_mass(v, w, α)
	R = (6 * π^2)^(1/3)
	integrand(τ) = τ^2 * exp(-τ) * (3 * erf(sqrt(D(τ, v, w) / 2) * R) / (2 * D(τ, v, w)^(5/2)) - R * exp(-D(τ, v, w) * R^2 / 2) * (D(τ, v, w) * R^2 + 3) / (sqrt(2π) * D(τ, v, w)^2))
	integral = α * quadgk(τ -> integrand(τ), 0.0, Inf)[1] / (3 * sqrt(π))
	return 1 + integral
end

# ╔═╡ 2041244d-1a95-4e9f-846b-0f64d5fca357
α[1:30]

# ╔═╡ 4d9a6e4c-0b35-4a55-ac12-4c21e52603b7
fro_m = [fro_eff_mass(vf[i], wf[i], α[i]) for i in 1:length(α)]

# ╔═╡ 4a4d9cd9-b991-40c9-9e09-18ab10f0e7ba
hol_m = [hol_eff_mass(vh_new[i], wh_new[i], α[i]) for i in 1:length(α)]

# ╔═╡ 308b415e-8965-4d31-b374-078add023990
begin
	pm = plot(α, 1 ./ fro_m, label = "Fro", xlabel = "α", ylabel = "m₀ / m*", linewidth = 2)
	plot!(α, 1 ./ hol_m, label = "Hol", linewidth = 2, linestyle = :dash)
end

# ╔═╡ 4b001539-7208-4100-ae75-6c1232f540f8
savefig(pm, "mass_comparison_alpha.pdf")

# ╔═╡ 614fe5ba-f265-420d-b914-0ecadedfbacc
md"""<commit_msg>Update holstein_variational.jl<commit_after>### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# ╔═╡ 66d2d0e0-ee56-11ec-21f9-5d6ae25cff48
begin
	using Optim
	using QuadGK
	using SpecialFunctions
	using Plots
end

# ╔═╡ 9f6474e5-41f1-4b96-aeee-d299404a3718
md"""
## Feynman's variational theory for the (continious) Frohlich model.
"""

# ╔═╡ 9429ae6b-55e1-4fb2-9a88-ae7636c35bb5
begin
"""
    fF(τ, v, w)
Integrand of Eqn. (31) in Feynman 1955. Part of the overall ground-state energy expression.
See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660.
"""
fF(τ, v, w) = (abs(w^2 * τ + (v^2 - w^2) / v * (1 - exp(-v*τ))))^(-0.5) * exp(-τ)


"""
    AF(v, w, α)
Integral of Eqn. (31) in Feynman 1955. Part of the overall ground-state energy expression.
See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660.
"""
AF(v, w, α) = π^(-0.5) * α * v * quadgk(τ -> fF(τ, v, w), 0.0, Inf)[1]

"""
    F(τ, v, w)
Ground state energy expression. Eqn. (33) in Feynman 1955.
See Feynman 1955: http://dx.doi.org/10.1103/PhysRev.97.660.
"""
F(v, w, α) = (3 / (4 * v)) * (v - w)^2 - AF(v, w, α)

# Let's wrap the Feynman athermal variation approximation in a simple function

"""
    feynmanvw(α; v = 0.0, w = 0.0)
Calculate v and w variational polaron parameters, for the supplied α Frohlich coupling. Returns (v, w).
This version uses the original athermal action (Feynman 1955).
"""
function feynmanvw_old(α; v = 3.0, w = 3.0) # v, w defaults

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = F(x[1] + x[2], x[2], α)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end
end

# ╔═╡ 1cf1e5cc-301c-4779-86cb-5b6436ecebf0
v, w = feynmanvw_old(2.39)

# ╔═╡ 2a609c87-4579-490e-99e7-0d0c25e11333
F(v, w, 2.39)

# ╔═╡ 7bf82f19-e6fc-4c3f-b7e8-a02bebf0d0dd
md"""
# Explicit k-space integral version of Frohlich model.
"""

# ╔═╡ 6cd58deb-fb04-44fc-a854-dce5b5cb8565
function D(t, v, w)
	abs(w^2 / v^2 * t + (v^2 - w^2) / v^3 * (1 - exp(-v * t)))
end

# ╔═╡ f7a07dc5-ba45-4870-b94c-44ae5d85436d
md"""
### Test first that it reproduces the Feynman analytic result. (R = 0).
"""

# ╔═╡ cc237871-08fd-4658-8fa5-8658d9394b3b
md"""
### Now try Frohlich model with k-space integral.
"""

# ╔═╡ 907b2475-0749-4d85-b09e-6eed75d35e68
md"""
### Now the k-space integral for the Holstein can be done and is expressed in terms of an erf function. So lets do it as it will be quicker!
"""

# ╔═╡ e7a4f46e-8fcc-48ef-9c5c-5790b42d7648
"""
Let's now plot the α dependence of Frohlich vs Holstein!
"""

# ╔═╡ 78936a9d-9567-4def-bf89-85be9e62e12d
α = 0.1:0.1:12

# ╔═╡ 8d726e23-5a40-4a10-8405-7c9775a0c9d0
frohlich_params = feynmanvw_old.(α)

# ╔═╡ 2cdfdcbc-a6e7-444c-a9a0-47529ddb091a
begin
	vf = [i[1] for i in frohlich_params]
	wf = [i[2] for i in frohlich_params]
end

# ╔═╡ 35ce9514-521c-4985-a85e-f236f67417c9
frohlich_energies = [F(vf[i], wf[i], α[i]) for i in 1:length(α)]

# ╔═╡ 3f5069df-89e2-4ac7-96e4-f3c79936fdd9
md"""
## Thermal (free energy) version of the theory...
"""

# ╔═╡ 9c1284f7-f503-4c4d-be6a-e4e04043246a
"""
    A(v, w, β)

Hellwarth's A expression from Eqn. (62b) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
A(v, w, β) = 3 / β * (log(v / w) - 1/2 * log(2π * β) - log(sinh(v * β / 2) / sinh(w * β / 2)))

# ╔═╡ b36f7997-f29c-4a92-bb00-9e82ae7d2a05
"""
    Y(x, v, β)

Hellwarth's Y expression from Eqn. (62d) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression. Contained in denominator of the integrand of Eqn. (62c).

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
Y(x, v, β) = 1 / (1 - exp(-v * β)) * (1 + exp(-v * β) - exp(-v * x) - exp(v * (x - β)))

# ╔═╡ d4e1ac9e-84d7-4809-ac1c-c24a9263db84
"""
    f(x, v, w, β)

Integrand of Hellwarth's B expression from Eqn. (62c) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
f_Fro(x, v, w, β) = (exp(β - x) + exp(x)) / sqrt(abs(w^2 * x * (1 - x / β) + Y(x, v, β) * (v^2 - w^2) / v))


# ╔═╡ ca1d9705-b49e-43f8-bd17-ab4a7e28b578
D(x, v, w, β) = w^2 / v^2 * x * (1 - x / β) + Y(x, v, β) * (v^2 - w^2) / v^3

# ╔═╡ 00fc8b3d-c485-43e2-92f0-451539ecf1e6
"""
 k_integral

R is cutoff for Holstei n k-space integral (in k-space spherical volume, R)

if R==0, Frohlich model (power should be 0, as the K^2 cancels with the 1/K^2 dependency in the el-ph matrix element)

Power is the K^Power
"""
function k_integral(t, v, w; power = 0, R = 0)
	if R == 0 # Frohlich
		return quadgk(k -> k^power * exp(-k^2 * D(t, v, w) / 2), -Inf, Inf)[1] / sqrt(2pi)
	else # Holstein, note limits on integral, R=k-space volume
		return quadgk(k -> k^power * exp(-k^2 * D(t, v, w) / 2), 0, R)[1] / sqrt(2) / pi
	end
end

# ╔═╡ dcf7f24a-6ef4-4613-9778-e2637f5e6c16
function Frohlich_el_integral(α, v, w; power = 0, R = 0)
	α * quadgk(t -> exp(-t) * k_integral(t, v, w; power = power, R = R), 0.0, Inf, rtol = 1e-4)[1] / sqrt(pi)
end

# ╔═╡ 6714b04e-9b26-42d0-89eb-f6c510cd5f56
function Frohlich_energy(α, v, w; power = 0, R = 0)
	3 * (v - w)^2 / (4 * v) - Frohlich_el_integral(α, v, w; power = power, R = R)
end

# ╔═╡ 429d0098-9d78-47e5-915b-cf16c7fce660
function var_params_k(α; v = 3.0, w = 3.0, power = 0, R = 0)
	 # Limits of the optimisation.
    lower = [0.0, 0.3]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = Frohlich_energy(α, x[1] + x[2], x[2]; power = power, R = R)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 71118d82-37fb-4504-b944-c3461de91233
v_k, w_k = var_params_k(2.39; v = v, w = w)

# ╔═╡ 80cd6d37-92d3-4757-a62e-490cffb4b63c
vh_k, wh_k = var_params_k(2.39; v = v, w = w, power = 2, R = (6 * π^2)^(1/3))

# ╔═╡ 2be860db-52b0-4514-872f-db79dc454323
Frohlich_energy(2.39, v_k, w_k)

# ╔═╡ 410d234c-25fe-4752-926e-f6c5d263ef05
Frohlich_energy(2.39, vh_k, wh_k; power = 2, R = (6 * π^2)^(1/3))

# ╔═╡ 079d7328-c2e1-4681-a4e3-e19a53841d25
function holstein_elph_integral(α, v, w; R = (6 * π^2)^(1/3))
	d(t) = D(t, v, w)
	integrand(t) = erf(sqrt(d(t) / 2) * R) / (2 * (d(t))^(3/2)) - R * exp(-R^2 * d(t) / (2)) / (sqrt(2π) * d(t))
	integral = α * quadgk(t -> exp(-t) * integrand(t), 0, Inf, rtol = 1e-4)[1] / sqrt(π / 2)
	return integral
end

# ╔═╡ bc20477e-2ce3-4377-9af0-5d479256175a
function holstein_energy(α, v, w)
	3 * (v - w)^2 / (4 * v) - holstein_elph_integral(α, v, w)
end

# ╔═╡ 380c5c1b-01e1-455b-a757-46c87fc45285
function holstein_variational_solution(α; v = 3.0, w = 3.0)
	 # Limits of the optimisation.
    lower = [0.0, 0.3]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = holstein_energy(α, x[1] + x[2], x[2])

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 650e36ab-7ac8-4d0b-ad74-f5693c2c1586
holstein_energies = [holstein_energy(α[i], vf[i], wf[i]) for i in 1:length(α)]

# ╔═╡ b2454e06-146c-4fc1-8928-70c7982922ec
begin
	pe = plot(α, frohlich_energies .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Fro", linewidth = 2, xlabel = "alpha, α", ylabel = "Energy (meV)")
	plot!(α, holstein_energies .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Hol", linewidth = 2, linestyle = :dash)
end

# ╔═╡ 67154f5a-e61e-4a5e-abf8-215bf8c8646e
savefig(pe, "energy_comparison.pdf")

# ╔═╡ 4e99725c-0895-4bc3-97ba-70eb675c7f47
f_Hol(x, v, w, β) = (exp(β - x) + exp(x)) * (erf(sqrt(D(x, v, w, β) / 2) * (6 * π^2)^(1/3)) / (2 * D(x, v, w, β)^(3/2)) - (6 * π^2)^(1/3) * exp(-D(x, v, w, β) * (6 * π^2)^(2/3) / 2) / (sqrt(2π) * D(x, v, w, β)))

# ╔═╡ e8b8e4de-e548-447e-82c5-979124ce6cf3
"""
    B(v, w, β, α)

Hellwarth's B expression from Eqn. (62c) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
B_Fro(v, w, β, α) = α * v / (sqrt(π) * (exp(β) - 1)) * quadgk(x -> f_Fro(x, v, w, β), 0, β/2, rtol = 1e-5)[1]


# ╔═╡ f87703c9-fa55-420a-b3cc-c6b947ef456a
B_Hol(v, w, β, α) = α / (sqrt(2) * pi * (exp(β) - 1)) * quadgk(x -> f_Hol(x, v, w, β), 0, β/2)[1]

# ╔═╡ c1358b12-4f68-4485-8ee7-73ef672f3ed8
"""
    C(v, w, β)

Hellwarth's C expression from Eqn. (62e) in Hellwarth et al. 1999 PRB. Part of the overall free energy expression.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
C(v, w, β) = 3/4 * (v^2 - w^2) / v * (coth(v * β / 2) - 2 / (v * β))

# ╔═╡ fc699405-e1d6-48dd-aff3-a7ae7fb0114b
"""
    F(v, w, β, α)

Hellwarth's total free energy expression from Eqn. (62a) in Hellwarth et al. 1999 PRB.

See Hellwarth et a. 1999: https://doi.org/10.1103/PhysRevB.60.299.
"""
F_Fro(v, w, β, α) = -(A(v, w, β) + B_Fro(v, w, β, α) + C(v, w, β))

# ╔═╡ 9ed1bad6-9ac5-444f-8d02-974bdd955b09
F_Hol(v, w, β, α) = -(A(v, w, β) + B_Hol(v, w, β, α) + C(v, w, β))

# ╔═╡ 5d5bf742-89c2-45e6-81f6-3b22e90b7bd0
"""
Check that it roughly gives the Holstein ground-state energy (high beta)
"""

# ╔═╡ 2455abad-abf2-4b0a-bb3f-30d97de9a9ec
function feynmanvw(α, β; v = 3.0, w = 3.0) # v, w defaults

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = F_Fro(x[1] + x[2], x[2], β, α)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 4d5e8c8e-2225-43fd-a934-0c8dd9f21585
function holstein_variational_solution(α, β; v = 3.0, w = 3.0) # v, w defaults

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]
    Δv = v - w # defines a constraint, so that v>w
    initial = [Δv + 0.01, w]

    # Feynman 1955 athermal action 
    f(x) = F_Hol(x[1] + x[2], x[2], β, α)

    # Use Optim to optimise v and w to minimise enthalpy.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    Δv, w = Optim.minimizer(solution) # v=Δv+w

    # Return variational parameters that minimise the free energy.
    return Δv + w, w
end

# ╔═╡ 466bc0e9-4cd2-42f5-8bbc-1317b5c1c7b6
vh, wh = holstein_variational_solution(2.39; v = v, w = w)

# ╔═╡ a1d89e74-7ab6-4c8f-8748-63ede869decb
holstein_energy(2.39, vh, wh)

# ╔═╡ 60254d3a-97b0-4c7e-a6c3-ad9ea0d93c97
F_Hol(vh, wh, 330, 2.39)

# ╔═╡ 5f5785c4-b9be-4246-820d-3e3ccf841cc6
holstein_params = holstein_variational_solution.(α)

# ╔═╡ 2ce6cc8d-0582-4f08-9e2f-5317f2f5b2e0
begin
	vh_new = [i[1] for i in holstein_params]
	wh_new = [i[2] for i in holstein_params]
end

# ╔═╡ f5e2a904-5e2f-4dcd-8691-b7620d24a7f0
begin
	pv = plot(α, vf, label = "Fro v", legend = :topleft, linewidth = 2, xlabel = "alpha, α", ylabel = "v and w", minorgrid = true)
	plot!(α, wf, label = "Fro w", linewidth = 2)
	plot!(α, vh_new, label = "Hol v", linewidth = 2, linestyle = :dash)
	plot!(α, wh_new, label = "Hol w", linewidth = 2, linestyle = :dash)
end

# ╔═╡ a8852b06-3f60-4bac-a610-8b0a1b8cf560
savefig(pv, "vw_comparison.pdf")

# ╔═╡ af11d3d0-b93d-4916-baad-5c81cec45ab4
"""
Temperature dependence free-energy and v/w comparisons
"""

# ╔═╡ 28f3a13d-cee1-4415-be15-4cefc474207d
T_range = 1:400

# ╔═╡ 7cd14da3-799d-4a53-a446-10a6842e5214
hol_params = [holstein_variational_solution(2.39, 2π * 2.25 * 1e12 * 1.05e-34 / T / 1.38e-23) for T in T_range]

# ╔═╡ 2baddbea-169a-4558-9f47-1c43e1269f1f
fro_params = [feynmanvw(2.39, 2π * 2.25 * 1e12 * 1.05e-34 / T / 1.38e-23) for T in T_range]

# ╔═╡ 1d45213c-3921-4a67-a76c-f37db47a14ad
begin
	vh_β = [i[1] for i in hol_params]
	wh_β = [i[2] for i in hol_params]
end

# ╔═╡ 3d423a59-b295-4026-8666-e08c38f01747
begin
	vf_β = [i[1] for i in fro_params]
	wf_β = [i[2] for i in fro_params]
end

# ╔═╡ 0e771039-19b8-4cec-bb12-8ddcb4b6975e
hol_energy_β = [F_Hol(vh_β[T], wh_β[T], 2π * 2.25 * 1e12 * 1.05e-34 / T_range[T] / 1.38e-23, 2.39) for T in 1:length(T_range)]

# ╔═╡ 07e267b3-03e8-4b1f-af0a-6c5ea6c57684
fro_energy_β = [F_Fro(vf_β[T], wf_β[T], 2π * 2.25 * 1e12 * 1.05e-34 / T_range[T] / 1.38e-23, 2.39) for T in 1:length(T_range)]

# ╔═╡ 1d08e7de-026b-4b1c-ba29-d12ce9f2b4ef
begin
	pvt = plot(T_range, vf_β, label = "Fro v", linewidth = 2, legend = :topleft, xlabel = "Temperature (K)", ylabel = "v and w params")
	plot!(T_range, wf_β, label = "Fro w", linewidth = 2)
	plot!(T_range, vh_β, label = "Hol v", linestyle = :dash, linewidth = 2)
	plot!(T_range, wh_β, label = "Hol w", linestyle = :dash, linewidth = 2)
end

# ╔═╡ a6ff82d8-ac2f-4034-ad64-c3800d67ecfa
begin
	pft = plot(T_range, fro_energy_β .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Fro", xlabel = "Temperature (K)", ylabel = "Free energy (ħω)", linewidth = 2)
	plot!(T_range, hol_energy_β .* 1000 * 1.05e-34 / 1.6e-19 * 1e12 * 2.25 * 2π, label = "Hol", linewidth = 2, linestyle = :dash)
end

# ╔═╡ a0cd088c-0412-4b1d-ba32-441f59e20c3c
savefig(pvt, "vw_comparison_temp.pdf")

# ╔═╡ 35864bca-4813-45ab-b579-ba9b40492189
savefig(pft, "free_energy_comparison.pdf")

# ╔═╡ a835bb67-f4d0-4b42-8cda-de6b44936ca2
"""
Have a look at effective mass
"""

# ╔═╡ 5ac9e949-4b55-42f5-9dd4-1de517051112
function fro_eff_mass(v, w, α)
	integrand(τ) = τ^2 * exp(-τ) / (D(τ, v, w))^(3/2)
	integral = α * quadgk(τ -> integrand(τ), 0.0, Inf)[1] / (3 * sqrt(π))
	return 1 + integral
end

# ╔═╡ 6a848cac-dc22-4419-8833-7279a4940834
function hol_eff_mass(v, w, α)
	R = (6 * π^2)^(1/3)
	integrand(τ) = τ^2 * exp(-τ) * (3 * erf(sqrt(D(τ, v, w) / 2) * R) / (2 * D(τ, v, w)^(5/2)) - R * exp(-D(τ, v, w) * R^2 / 2) * (D(τ, v, w) * R^2 + 3) / (sqrt(2π) * D(τ, v, w)^2))
	integral = α * quadgk(τ -> integrand(τ), 0.0, Inf)[1] / (3 * sqrt(2) * pi)
	return 1 + integral
end

# ╔═╡ 2041244d-1a95-4e9f-846b-0f64d5fca357
α[1:30]

# ╔═╡ 4d9a6e4c-0b35-4a55-ac12-4c21e52603b7
fro_m = [fro_eff_mass(vf[i], wf[i], α[i]) for i in 1:length(α)]

# ╔═╡ 4a4d9cd9-b991-40c9-9e09-18ab10f0e7ba
hol_m = [hol_eff_mass(vh_new[i], wh_new[i], α[i]) for i in 1:length(α)]

# ╔═╡ 308b415e-8965-4d31-b374-078add023990
begin
	pm = plot(α, 1 ./ fro_m, label = "Fro", xlabel = "α", ylabel = "m₀ / m*", linewidth = 2)
	plot!(α, 1 ./ hol_m, label = "Hol", linewidth = 2, linestyle = :dash)
end

# ╔═╡ 4b001539-7208-4100-ae75-6c1232f540f8
savefig(pm, "mass_comparison_alpha.pdf")

# ╔═╡ 614fe5ba-f265-420d-b914-0ecadedfbacc
md"""
{'fraction_non_alphanumeric': 0.17783785463472235, 'fraction_numerical': 0.13302880581709706, 'mean_word_length': 3.1865487186158448, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 31, 'https://': 12, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4717009', 'n_tokens_mistral': 18376, 'n_tokens_neox': 15775, 'n_words': 3712}