Update tests to match new code changes (e.g. separating v and w from …

…phonon freq. normalisation). Corrected factor of 2 in Holstein el-ph matrix. Removed some questionable `@fastmath` macros causing some errors.
Frost-group · Jun 14, 2024 · 680bc1f · 680bc1f
1 parent b8ecc52
commit 680bc1f
Show file tree

Hide file tree

Showing 7 changed files with 54 additions and 53 deletions.
diff --git a/src/HolsteinPolaron.jl b/src/HolsteinPolaron.jl
@@ -155,14 +155,18 @@ function holsteinpolaron(αrange, Trange, Ωrange; ω=1, ωeff=1, γ=1, mb=1, J=
                 dprocess += 1
             end
 
+        # Update the guesses to keep them close-ish to the true solutions during loops over alphas.
+        if j > 1
+            v_guesses = reduce_array(p["v0"][d, j-1, :])
+            w_guesses = reduce_array(p["w0"][d, j-1, :])
+        end
+
         # Extract the ground-state, athermal polaron properties (energy (enthalpy) and variational parameters v and w).
         # w is also the frequency of oscillation of the SHM trial system composed of the bare particle and fictitous mass.
         # A, B, C are components of the total energy: A is the bare electron energy, B the electron-phonon interaction energy, C is the energy of the harmonic trial system.
         athermal_energy(v, w) = !kspace ? holstein_energy(v, w, α, J, ω; a = a, dims = dims[d]) : holstein_energy_k_space(v, w, α, J, ω; a = a, dims = dims[d])
         v_gs, w_gs, F_gs, A_gs, B_gs, C_gs = vw_variation(athermal_energy, v_guesses, w_guesses)
 
-        # Update the guesses to keep them close-ish to the true solutions during loops over alphas.
-        v_guesses, w_guesses = v_gs, w_gs
 
         # Store the athermal data.
         p["v0"][d, j, :] .= v_gs 

diff --git a/src/HolsteinTheory.jl b/src/HolsteinTheory.jl
@@ -24,7 +24,7 @@ Expected Output:
 6.0
 """
 function holstein_coupling(k, α, J, ω; dims = 3)
-    α * ω * J * dims
+    2 * α * ω * J * dims
 end
 
 holstein_coupling(k, α::Vector, ω::Vector; dims = 3) = sum(holstein_coupling(k, α[j], ω[j]; dims = dims) for j in eachindex(α))
@@ -59,13 +59,12 @@ This example calculates the electron-phonon interaction energy for given values
 """
 function holstein_interaction_energy(v, w, α, J, ωβ...; a = 1, dims = 3)
     coupling = holstein_coupling(1, α, J, ωβ[1]; dims = dims)
-    polaron_radius = a * sqrt(J / ωβ[1])
-	propagator(τ) = length(ωβ) == 1 ? polaron_propagator(τ, v, w) * ωβ[1] * polaron_radius^2 : polaron_propagator(τ, v, w, ωβ[2]) * ωβ[1] * polaron_radius^2
-    q_integral(τ) = (a / 2π) * sqrt(π / propagator(τ)) * erf(π / a * sqrt(propagator(τ)))
+	propagator(τ) = length(ωβ) == 1 ? polaron_propagator(τ, v, w) : polaron_propagator(τ, v, w, ωβ[2])
+    q_integral(τ) = (a / 2π) * sqrt(π / propagator(τ)) * erf(π * sqrt(propagator(τ)))
 	integrand(τ) = phonon_propagator(τ, ωβ...) * q_integral(τ)^dims
     upper_limit = length(ωβ) == 1 ? Inf : ωβ[2] / 2
     integral, _ = quadgk(τ -> integrand(τ), 0, upper_limit)
-    return integral * coupling
+    return integral * coupling / 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(α))
@@ -109,7 +108,7 @@ holstein_interaction_energy_k_space(v, w, α::Vector, ω::Vector; dims = 3) = su
 function holstein_energy(v, w, α, J, ωβ...; a = 1, dims = 3)
 	A, C = length(ωβ) == 1 ? trial_energy(v, w; dims = dims) : trial_energy(v, w, ωβ[2]; dims = dims)
 	B = holstein_interaction_energy(v, w, α, J, ωβ...; a = a, dims = dims)
-    return -(A + B + C), A, B, C
+    return -2 * dims -(A + B + C), A, B, C
 end 
 
 """

diff --git a/src/PolaronFunctions.jl b/src/PolaronFunctions.jl
@@ -194,7 +194,7 @@ function polaron_memory_function(Ω, structure_factor; limits = [0, Inf])
     if iszero(Ω)
       return polaron_memory_function(structure_factor; limits = limits)
     end
-    @fastmath integral, _ = quadgk(t -> (1 - exp(im * Ω * t)) / Ω * imag(structure_factor(t)), limits..., rtol=1e-4)
+    integral, _ = quadgk(t -> (1 - exp(im * Ω * t)) / Ω * imag(structure_factor(t)), limits..., rtol=1e-4)
     return integral
 end
 
@@ -225,7 +225,7 @@ println(result)  # Output: 383.3333333333333
 ```
 """
 function polaron_memory_function(structure_factor; limits = [0, Inf])
-  @fastmath integral, _ = quadgk(t -> -im * t * imag(structure_factor(t)), limits..., rtol=1e-5)
+  integral, _ = quadgk(t -> -im * t * imag(structure_factor(t)), limits..., rtol=1e-5)
   return integral
 end
 

diff --git a/test/FeynmanAthermal.jl b/test/FeynmanAthermal.jl
@@ -16,13 +16,13 @@ rows(M::Matrix) = map(x->reshape(getindex(M, x, :), :, size(M)[2]), 1:size(M)[1]
 
 for r in rows(Schultz)
     α,vSchultz,wSchultz,ESchultz=r # Unpack row of results
-    v,w,E=feynmanvw(3.1, 3.0, α, 1.0) # performs the optimisation
+    v,w,E=feynmanvw(α, 1.0) # performs the optimisation
 
     println("α=$α v=$v w=$w E=$E  | Schultz: v=$vSchultz w=$wSchultz E=$ESchultz")
 
     @test v ≈ vSchultz atol=0.1 # Strangely these need more tolerance than the Energies
     @test w ≈ wSchultz atol=0.1
-    @test E ≈ ESchultz atol=0.001
+    @test E ≈ ESchultz atol=0.01
 end
 
 end

diff --git a/test/FrostPolaronMobility2017.jl b/test/FrostPolaronMobility2017.jl
@@ -23,13 +23,13 @@
     @test ustrip(MAPIe_polaron.μ |> u"cm^2/V/s") ≈ 136.42 rtol=0.02
 
     # Test variational parameters
-    @test MAPIe_polaron.v ≈ 19.86 rtol=0.02
-    @test MAPIe_polaron.w ≈ 16.96 rtol=0.02
+    @test ustrip(MAPIe_polaron.v) / phonon_freq ≈ 19.86 rtol=0.02
+    @test ustrip(MAPIe_polaron.w) / phonon_freq ≈ 16.96 rtol=0.02
 
     # Same for the MAPI holes @ 300 K
     @test ustrip(MAPIh_polaron.μ |> u"cm^2/V/s") ≈ 94.15 rtol=0.02
-    @test MAPIh_polaron.v ≈ 20.09  rtol=0.02
-    @test MAPIh_polaron.w ≈ 16.81  rtol=0.02
+    @test ustrip(MAPIh_polaron.v) / phonon_freq ≈ 20.09  rtol=0.02
+    @test ustrip(MAPIh_polaron.w) / phonon_freq ≈ 16.81  rtol=0.02
 
 end
 
diff --git a/test/MultipleBranches.jl b/test/MultipleBranches.jl
@@ -169,9 +169,9 @@
     println("alphas: ", alphas)
     println("sum alphas: ", sum(alphas))
 
-    println("Feynman v,w for alphas: ", feynmanvw.(3.1, 3.0, alphas, 1.0))
+    println("Feynman v,w for alphas: ", feynmanvw.(alphas, 1.0))
 
-    mobilityproblem = hcat(alphas, feynmanvw.(3.1, 3.0, alphas, 1.0), MAPI[:, 1])
+    mobilityproblem = hcat(alphas, feynmanvw.(alphas, 1.0), MAPI[:, 1])
     println("Specify mobility problem: ", mobilityproblem)
 
     """
@@ -219,7 +219,7 @@
         alphas = frohlichPartial.(eachrow(f_dielectric), ϵ_o=ϵ_o, ϵ_s=ϵ_o + sum(splat), meff=meff)
         βred = ħ * MAPI[:, 1] * 1E12 * 2π / (kB * T)
 
-        mobilityproblem = hcat(alphas, feynmanvw.(3.1, 3.0, alphas, 1.0, βred), MAPI[:, 1])
+        mobilityproblem = hcat(alphas, feynmanvw.(alphas, 1.0, βred), MAPI[:, 1])
         inverse_μ = Hellwarth1999mobilityRHS.(eachrow(mobilityproblem), meff, T)
         μ = sum(inverse_μ)^-1
 

diff --git a/test/MultiplePhonons.jl b/test/MultiplePhonons.jl
@@ -58,18 +58,20 @@
 
     # Ionic dielectric and decomposed alpha
 
-    ω = 2π .* phonon_freq
+    ω = phonon_freq .* 2π
     ϵ_ionic = [ϵ_ionic_mode(i, j, volume) for (i, j) in zip(phonon_freq, ir_activity)]
 
     α = [frohlichalpha(ϵ_optic, i, sum(ϵ_ionic), j, m_eff) for (i, j) in zip(ϵ_ionic, phonon_freq)]
 
+    Hellwarth_B_freq = HellwarthBScheme(hcat(phonon_freq, ir_activity))
+
     @testset "Ionic dielectric and decomposed alphas" begin
 
-        @test ϵ_ionic ≈ [0.2999680470756664, 0.0247387647244569, 0.2543132184018061, 0.16621617310133838, 2.3083204422506296, 3.0707979601813267, 3.2026782087486407, 0.892674135624958, 0.8861579096771846, 0.7209278375756829, 0.40199759805819046, 0.6183279038315278, 0.7666391525823296, 0.1555444994147378, 1.1627710200840813] rtol = 1e-3
+        @test ϵ_ionic ≈ [0.299961772906185, 0.024738247285335784, 0.2543078991544757, 0.16621269650291148, 2.3082721611293127, 3.07073373098414, 3.2026112211275524, 0.8926554643402222, 0.8861393746866522, 0.7209127585581645, 0.40198918982576065, 0.6183149708071556, 0.7666231174611848, 0.1555412460263828, 1.1627466994188664] rtol = 1e-3
 
-        @test α ≈ [0.03401013445306177, 0.002850969846883158, 0.03075081562607006, 0.02275292381052607, 0.33591418423943553, 0.46526818717696034, 0.5046331098089347, 0.14223560721522646, 0.16083871312929882, 0.1622911897190622, 0.09123913334086006, 0.14070972715961402, 0.18159786148348117, 0.039500396977008016, 0.3487735470060312] rtol = 1e-3
+        @test α ≈ [0.03400925262932455, 0.002850895926184945, 0.030750018310811124, 0.0227523338667158, 0.33590547456785985, 0.4652561235806758, 0.5046200255485525, 0.14223191929288398, 0.1608345428607191, 0.16228698179028636, 0.09123676766854395, 0.1407060788006921, 0.1815931529662271, 0.03949937280029587, 0.34876450391350566] rtol = 1e-3
 
-        @test sum(α) ≈ 2.663366500992453 rtol = 1e-3
+        @test sum(α) ≈ 2.6632974445232787 rtol = 1e-3
 
         println('\n', "Ionic dielectric and alphas:")
         println("-------------------------------")
@@ -81,20 +83,18 @@
 
     # Variations
 
-    athermal_energy(v, w) = frohlich_energy(v, w, α, ω)
-    v_0, w_0, E_0, A_0, B_0, C_0 = vw_variation(athermal_energy, 4, 3) # Athermal
+    v_0, w_0, E_0, A_0, B_0, C_0 = feynmanvw(α, ω) # Athermal
 
     β = ħ / (kB * 300) * 1e12
-    thermal_energy(v, w) = frohlich_energy(v, w, α, ω, β)
-    v, w, E = vw_variation(thermal_energy, v_0, w_0) # Thermal
+    v, w, E, A, B, C = feynmanvw(α, ω, β) # Thermal
 
     @testset "Multiple mode variations" begin
 
-        @test v_0[1] ≈ 3.292283619446986 rtol = 1e-3
-        @test w_0[1] ≈ 2.679188425097246 rtol = 1e-3
+        @test v_0[1] / Hellwarth_B_freq ≈ 16.283245452806543 rtol = 1e-3
+        @test w_0[1] / Hellwarth_B_freq ≈ 12.93471218329606 rtol = 1e-3
 
-        @test v[1] ≈ 35.19211042393129 rtol = 1e-3
-        @test w[1] ≈ 32.454157668863225 rtol = 1e-3
+        @test v[1] / Hellwarth_B_freq ≈ 123.90465526352604 rtol = 1e-3
+        @test w[1] / Hellwarth_B_freq ≈ 106.7711782241155 rtol = 1e-3
 
         println("\nVariational parameters:")
         println("Athermal v = $(v_0[1]) | athermal w = $(w_0[1])")
@@ -109,9 +109,9 @@
 
     @testset "Multiple mode energies" begin
 
-        @test E_0 ≈ -19.51150570496287 rtol = 1e-3
+        @test E_0 ≈ -19.54914965269833 rtol = 1e-3
 
-        @test E ≈ -42.79764110613318 rtol = 1e-3
+        @test E ≈ -63.315841394023664 rtol = 1e-3
 
         println("\nPolaron energies")
         println("0K energy: $E_0 meV")
@@ -124,7 +124,7 @@
 
     @testset "Multiple mode mobility" begin
 
-        @test μ ≈ 160.39270615550004 rtol = 1e-3
+        @test μ ≈ 145.35917679067404  rtol = 1e-3
 
         println("\nMobility at 300K: $μ")
     end
@@ -137,9 +137,9 @@
 
     @testset "Multiple mode impedence" begin
 
-        @test Z_dc ≈ 91.38457766169273 + 0.0im rtol = 1e-3
+        @test Z_dc ≈ 100.83185205569625 + 0.0im  rtol = 1e-3
 
-        @test Z ≈ 76.18909903384355 - 23.457492887508455im rtol = 1e-3
+        @test Z ≈ 82.51429805705051 - 18.759832224634735im rtol = 1e-3
 
         println("\nComplex Impedence:")
         println("DC limit: $Z_dc")
@@ -154,9 +154,9 @@
 
     @testset "Multiple mode conductivity" begin
 
-        @test σ_dc ≈ 0.0109427654598571 - 0.0im rtol = 1e-3
+        @test σ_dc ≈ 0.009917501063529335 - 0.0im rtol = 1e-3
 
-        @test σ ≈ 0.011988781430646566 + 0.003691167879729921im rtol = 1e-3
+        @test σ ≈ 0.011523473106500403 + 0.0026198904579370222im rtol = 1e-3
 
         println("\nComplex Conductivity:")
         println("DC limit: $σ_dc")
@@ -165,8 +165,6 @@
 
     # Single Mode Tests
 
-    Hellwarth_B_freq = HellwarthBScheme(hcat(phonon_freq, ir_activity))
-
     @testset "Hellwarth effective frequency" begin
 
         @test Hellwarth_B_freq ≈ 2.25 rtol = 1e-3
@@ -181,13 +179,13 @@
         display(singlemode_polaron)
 
         @test singlemode_polaron.α ≈ 2.393156008589176 rtol = 1e-3
-        @test singlemode_polaron.v ≈ [3.3086321909253815, 19.847527941504232] rtol = 1e-3
-        @test singlemode_polaron.w ≈ [2.6634018089397014, 16.948211346396874] rtol = 1e-3
-        @test singlemode_polaron.F ≈ [-5.571343967852482, -8.585426348039439] rtol = 1e-3
+        @test singlemode_polaron.v ≈ [7.449619526820557, 44.68801621691182] rtol = 1e-3
+        @test singlemode_polaron.w ≈ [5.996837775164936, 38.160005552975356] rtol = 1e-3
+        @test singlemode_polaron.F ≈ [-5.5702304614638845, -16.23756639407282] rtol = 1e-3
         @test singlemode_polaron.χ ≈ [1.8185958394302872 + 2.8215456392868816im, -2.6431955087396575 + 16.602938172029397im] rtol = 1e-3
-        @test singlemode_polaron.z ≈ [0.3385854767144258 - 0.5782315007316344im, 1.9923525806435276 - 0.0428165389512411im] rtol = 1e-3
-        @test singlemode_polaron.σ ≈ [0.7541017043762221 + 1.287844252674551im, 0.5016874943626112 + 0.010781486345548893im] rtol = 1e-3
-        @test singlemode_polaron.μ ≈ [Inf, 0.4873858565327472] rtol = 1e-3
+        @test singlemode_polaron.z ≈ [2.821550935533812 - 4.818576894764461im, 16.60369379393479 - 0.356819708523505im] rtol = 1e-3
+        @test singlemode_polaron.σ ≈ [0.09049281752137871 + 0.15454145950697004im, 0.06019975971109561 + 0.0012937157827582215im] rtol = 1e-3
+        @test singlemode_polaron.μ ≈ [Inf, 0.058486537513714805] rtol = 1e-3
     end
 
     # Multimode Tests
@@ -202,13 +200,13 @@
 
         @test multimode_polaron.α ≈ [0.03401013640864639, 0.0028509700108140822, 0.030750817394243672, 0.02275292511882046, 0.33591420355451984, 0.46526821392990697, 0.5046331388253664, 0.14223561539378177, 0.16083872237753374, 0.16229119905081466, 0.0912391385871153, 0.14070973525043112, 0.18159787192536828, 0.03950039924828304, 0.3487735670605295] rtol = 1e-3
         @test sum(multimode_polaron.α) ≈ 2.663366654136176 rtol = 1e-3
-        @test multimode_polaron.v ≈ [3.292268504574272, 35.19207894065109] rtol = 1e-3
-        @test multimode_polaron.w ≈ [2.679193078382525, 32.45420279136673] rtol = 1e-3
-        @test multimode_polaron.F ≈ [-4.719036156508013, -10.357755243151143] rtol = 1e-3
-        @test multimode_polaron.χ ≈ [1.0521976989565067 + 2.089861861859462im, -2.1193321531426945 + 14.077489512387219im] rtol = 1e-3
-        @test multimode_polaron.z ≈ [0.2507834234231354 - 0.4862637238747808im, 1.6892987414864662 - 0.10568014162287666im] rtol = 1e-3
-        @test multimode_polaron.σ ≈ [0.8377746271072964 + 1.624427182564005im, 0.5896539523503038 + 0.036887917845742565im] rtol = 1e-3
-        @test multimode_polaron.μ ≈ [Inf, 0.5729621483218189] rtol = 1e-3
+        @test multimode_polaron.v ≈ [5.835086623282841, 44.40143374273715] rtol = 1e-3
+        @test multimode_polaron.w ≈ [4.635142386961883, 38.26163944180995] rtol = 1e-3
+        @test multimode_polaron.F ≈ [-4.728298068562315, -15.318577804140983]  rtol = 1e-3
+        @test multimode_polaron.χ ≈ [1.0083819418968094 + 2.262844921017003im, -2.5135009156247414 + 15.484273072268337im] rtol = 1e-3
+        @test multimode_polaron.z ≈ [2.262844921017003 - 4.008381941896809im, 15.484273072268337 - 0.4864990843752586im] rtol = 1e-3
+        @test multimode_polaron.σ ≈ [0.10680047179649949 + 0.18918533857934283im, 0.0645179673929453 + 0.002027084637162277im] rtol = 1e-3
+        @test multimode_polaron.μ ≈ [Inf, 0.062313430372364004] rtol = 1e-3
     end
 
     print('\n')