Fix temperature conventions

scripts/01_renormalization_factors_for_sum_rule.jl

@@ -3,7 +3,7 @@ using DrWatson
 
 using Sunny, GLMakie
 include(srcdir("model.jl"))
-include(srcdir("structure_factor_utils.jl"))
+include(srcdir("sum_rule_utils.jl"))
 
 
 # Function to estimate κ at a given temperature 
@@ -65,7 +65,7 @@ sim_params = (;
 global_bounds = (1.0, 2.0)  # Smallest and largest κs (search space)
 thresh = 0.05               # Allowable deviation in estimated sum, relative to reference 
 ref = 16/3                  # Quadratic Casimir of SU(3) for chosen normalization convention
-kTs = 10 .^ range(log10(0.1), log10(20.0), 10)  # 10 temperatures between 0.1 and 100.0, in meV
+kTs = 10 .^ range(log10(0.1), log10(30.0), 15)  # 15 temperatures between 0.1 and 10.0 in K
 
 
 # Estimate κs for chosen temperatures
@@ -85,4 +85,4 @@ wsave(datadir("kappas", "kappas.jld2"), data)
 # Note: For publication quality results, a larger system should be used and many more
 # samples should be collected for each estimate of κ. Additionally, the same
 # sampled initial conditions should be used for each kT to avoid the possibility
-# of locking the binary search due to stochastic effects.
+# of locking the binary search due to stochastic effects.

scripts/02_classical_ordering_temperature.jl

@@ -10,10 +10,12 @@ dur_therm = 15.0    # Conservative thermalization time for all temperatures
 Δt = 0.004          # Step size sufficient to ensure stability
 λ = 0.1             # Phenomenological coupling to thermal bath
 
+
 # Set up a System
 dims = (12, 12, 4)
 gs = 1 seed = 101
 
+
 # Estimate the Bragg intensities at a range of temperatures
 kTs = 10 .^ range(log10(0.1), log10(10), 20) # Temperature range in Kelvin. 
                                              # Manuscript used 88 temperatures between 0.1 and 77 Kelvin.
@@ -56,16 +58,23 @@ for kT in kTs
     push!(σs, std(bragg_samples))
 end
 
+
 # Calculate numerical derivative of the Bragg intensities as a function of temperature.
 Δμs = μs[2:end] .- μs[1:end-1]
 ΔkTs = kTs[2:end] .- kTs[1:end-1]
 Δμ_ΔkT = Δμs ./ ΔkTs
 kTs_centered = (kTs[2:end] .+ kTs[1:end-1]) ./ 2
 
+
 # Plot the results 
 fig = Figure()
 xticks = [1, 2, 4, 8]
 ax1 = Axis(fig[1,1]; xscale=log10, xticks, xlabel="kT (K)", ylabel="I_Bragg/I_Bragg-max")
 ax2 = Axis(fig[1,2]; xscale=log10, xticks, xlabel="kT (K)", ylabel="ΔI_Bragg/ΔkT")
 scatter!(ax1, kTs, μs ./ maximum(μs))
-scatter!(ax2, kTs_centered, Δμ_ΔkT ./ maximum(μs))
+scatter!(ax2, kTs_centered, Δμ_ΔkT ./ maximum(μs))
+
+
+# Inspection of the numerical derivative should show a discontinuity at about 
+# T = 3 K. A more careful calculation on a sufficiently large system will show a
+# result of T = 3.05 K, which we will use in subsequent scripts.

scripts/03_dssfs.jl

@@ -0,0 +1,15 @@
+using DrWatson
+@quickactivate "FiniteTemperatureFeI2"
+
+using Sunny
+include(srcdir("model.jl"))
+include(srcdir("decorrelation_utils.jl"))
+
+# After executing the first script, you should have some estimated κ values
+# for a range of temperatures. In this script we will show how to estimate
+# the dynamical spin structure factor at one of the temperatures for which
+# a κ is available. We begin be loading the earlier results.
+
+@unpack κs, kTs = load(datadir("kappas", "kappas.jld2"))
+
+# In the second script, we found that T_N is approximately 3 K. We'll pick

src/sum_rule_utils.jl

@@ -62,7 +62,7 @@ function estimate_sum(sys::System{N}, κ, kT, sim_params; observables=nothing) w
     ntherm = round(Int, dur_therm/Δt_therm)
 
     # Thermalize the system
-    langevin = Langevin(Δt_therm; λ, kT)
+    langevin = Langevin(Δt_therm; λ, kT = kT*Sunny.meV_per_K)
     for _ in 1:ntherm
         step!(sys, langevin)
     end
@@ -88,5 +88,5 @@ function estimate_sum(sys::System{N}, κ, kT, sim_params; observables=nothing) w
     end
 
     # Return spectral weight per site.
-    return total_spectral_weight(sc; kT) / length(Sunny.eachsite(sys))
+    return total_spectral_weight(sc; kT = kT*Sunny.meV_per_K) / length(Sunny.eachsite(sys))
 end