add DensityGHQField (#155)

src/MonteCarlo.jl

@@ -118,7 +118,7 @@ export Model, IsingModel
 export HubbardModel, HubbardModelAttractive, HubbardModelRepulsive
 export MonteCarloFlavor, MC, DQMC
 export greens, greens!, lattice, model, parameters
-export DensityHirschField, MagneticHirschField, MagneticGHQField
+export DensityHirschField, MagneticHirschField, DensityGHQField, MagneticGHQField
 
 # For extending
 export AbstractMeasurement, Model

src/flavors/DQMC/fields.jl

@@ -396,13 +396,22 @@ interaction_eltype(::AbstractGHQField{T}) where {T} = T
 interaction_matrix_type(::AbstractGHQField{Float64}, ::Model) = Diagonal{Float64, FVec64}
 interaction_matrix_type(::AbstractGHQField{ComplexF64}, ::Model) = Diagonal{ComplexF64, CVec64}
 
+@inline @bm function accept_local!(
+        mc, f::AbstractGHQField, i, slice, detratio, ΔE_boson, x_new
+    )
+    update_greens!(mc.stack.field_cache, mc.stack.greens, i, size(f.conf, 1))
+    @inbounds f.conf[i, slice] = x_new
+    return nothing
+end
+
+
 """
     MagneticGHQField
 
 This represents a field generated by using (4 node) Gauss Hermite quadrature to 
 solve the integral form of the Hubbard Stratonovich transformed interaction. It 
 is a more general solution than the Hirsch transform, allowing any interaction 
-that can eb written as V ~ ² though the current implementation only allows
+that can be written as V ~ ² though the current implementation only allows
 Hubbard interactions.
 
 In general:
@@ -465,15 +474,70 @@ end
     return detratio * f.γ[x_new] / f.γ[x_old], 0.0, x_new
 end
 
-@inline @bm function accept_local!(
-        mc, f::MagneticGHQField, i, slice, detratio, ΔE_boson, x_new
+
+"""
+    DensityGHQField
+
+This represents a field generated by using (4 node) Gauss Hermite quadrature to 
+solve the integral form of the Hubbard Stratonovich transformed interaction. It 
+is a more general solution than the Hirsch transform, allowing any interaction 
+that can be written as V ~ ² though the current implementation only allows
+Hubbard interactions.
+
+In general:
+exp(Δτ U ²) ~ ∑ₓ γ(x) exp(α η(x) Â)
+with α = sqrt(Δτ U) and x ∈ {-2, -1, 1, 2}
+"""
+struct DensityGHQField{T} <: AbstractGHQField{T}
+    α::T
+    γ::FVec64
+    η::FVec64
+    choices::Matrix{Int8}
+
+    temp_conf::Matrix{Int8}
+    conf::Matrix{Int8}
+end
+
+function DensityGHQField(param::DQMCParameters, model::Model, U::Number = model.U)
+    α = maybe_to_float(sqrt(0.5 * param.delta_tau * ComplexF64(U)))
+    s6 = sqrt(6)
+    gammas = Float64[1 - s6/3, 1 + s6/3, 1 + s6/3, 1 - s6/3]
+    etas = Float64[-sqrt(6 + 2s6), -sqrt(6 - 2s6), sqrt(6 - 2s6), sqrt(6 + 2s6)]
+    choices = Int8[2 3 4; 1 3 4; 1 2 4; 1 2 3]
+
+    DensityGHQField(
+        α, gammas, etas, choices,
+        Matrix{Int8}(undef, length(lattice(model)), param.slices),
+        Matrix{Int8}(undef, length(lattice(model)), param.slices)
     )
-    update_greens!(mc.stack.field_cache, mc.stack.greens, i, size(f.conf, 1))
-    @inbounds f.conf[i, slice] = x_new
+end
+
+nflavors(::DensityGHQField) = 1
+energy_boson(f::DensityGHQField, conf = f.conf) = f.α * sum(conf)
+
+
+# TODO: Maybe worth adding a complex method?
+@inline @bm function interaction_matrix_exp!(f::DensityGHQField, result::Diagonal, slice, power)
+    @inbounds for i in eachindex(result.diag)
+        result.diag[i] = exp(+power * f.α * f.η[f.conf[i, slice]])
+    end
     return nothing
 end
 
 
+@inline @bm function propose_local(mc, f::DensityGHQField, i, slice)
+    x_old = f.conf[i, slice]
+    x_new = @inbounds f.choices[x_old, rand(1:3)]
+
+    ΔE_boson = f.α * (f.η[x_new] - f.η[x_old])
+    exp_ratio = exp(ΔE_boson)
+    mc.stack.field_cache.Δ = exp_ratio - 1.0
+    detratio = calculate_detratio!(mc.stack.field_cache, mc.stack.greens, i)
+
+    return detratio * f.γ[x_new] / f.γ[x_old], ΔE_boson, x_new
+end
+
+
 ################################################################################
 ### Util
 ################################################################################

test/ED/ED_tests.jl

@@ -114,7 +114,7 @@ end
                 thermalization = 1, sweeps = 2, measure_rate=1#, field = field
             )
             print(
-                "  Running DQMC ($(nameof(typeof(model))) $(nameof(field))) " * 
+                "  Running DQMC ($(nameof(typeof(model)))) " * 
                 "β=$(dqmc.parameters.beta)\n    "
             )
 
@@ -391,6 +391,7 @@ end
 
         end
     end
+    println()
 end
 
 
@@ -418,7 +419,7 @@ end
                 )
             )
             print(
-                "  Running DQMC ($(nameof(typeof(model))) $(nameof(field))) " * 
+                "  Running DQMC ($(nameof(typeof(model)))) " * 
                 "β=$(dqmc.parameters.beta), 5k + 5k sweeps\n    "
             )
 
@@ -692,6 +693,7 @@ end
             end
         end
     end
+    println()
 end
 
 
@@ -704,9 +706,9 @@ end
         HubbardModel(2, 2, U = -1.0, t = 1.0),
         HubbardModel(2, 2, U = 1.0, mu = 1.0, t = 1.0)
     )
-    fields = (MagneticHirschField, DensityHirschField, MagneticGHQField)
+    fields = (MagneticHirschField, DensityHirschField, MagneticGHQField, DensityGHQField)
 
-    println("Finite U ED Comparison")
+    println("Greens and Energy vs ED with various fields")
     for model in models, field in fields
         if field == MonteCarlo.choose_field(model)
             continue
@@ -723,7 +725,7 @@ end
             )
             str = model.U >= 0 ? "attractive" : "repulsive"
             print(
-                "  Running DQMC ($mod $(nameof(field))) " * 
+                "  Running DQMC ($str $(nameof(field))) " * 
                 "β=$(dqmc.parameters.beta), 5k + 5k sweeps\n    "
             )
 
@@ -737,8 +739,9 @@ end
             rtol = 2dqmc.parameters.delta_tau^2
             N = length(lattice(model))
 
-            print("    Running ED and checking (tolerance: $atol, $(100rtol)%)\n    ")
-            @time begin
+            # print("    Running ED and checking (tolerance: $atol, $(100rtol)%)\n    ")
+            # @time 
+            begin
                 H = HamiltonMatrix(model)
 
                 @testset "(total) energy" begin