another attempt at superfluid density/stiffness

docs/src/DQMC/scheduler.md

@@ -1,6 +1,6 @@
 # Update Scheduler
 
-The update scheduler keeps track of and iterates through various Monte Carlo updates. Currently there are two schedulers, `SimpleScheduler` and `AdaptiveScheduler`, and five (full) updates.
+The update scheduler keeps track of and iterates through various Monte Carlo updates. Currently there are two schedulers, `SimpleScheduler` and `AdaptiveScheduler`, and multiple (full) updates.
 
 !!! note
 

src/configurations.jl

@@ -228,8 +228,15 @@ function _load(data, ::Val{:BufferedConfigRecorder})
     if !(data["VERSION"] <= 2)
         throw(ErrorException("Failed to load BufferedConfigRecorder version $(data["VERSION"])"))
     end
+    filename = data["filename"]
+    if filename isa FilePath && filename.is_relative && !isfile(filename.absolute_path)
+        filename = FilePath(
+            true, filename.relative_path, 
+            joinpath(splitdir(data.path)[1], filename.relative_path)
+        )
+    end
     BufferedConfigRecorder(
-        data["filename"], data["buffer"], data["rate"], -1, -1, 
+        filename, data["buffer"], data["rate"], -1, -1, 
         data["total_length"], data["save_idx"]
     )
 end

src/flavors/DQMC/measurements/measurements.jl

@@ -160,57 +160,78 @@ function current_current_susceptibility(
     li = wrapper === nothing ? lattice_iterator : wrapper{lattice_iterator}
     Measurement(dqmc, model, greens_iterator, li, cc_kernel; kwargs...)
 end
+
 function superfluid_density(
         dqmc::DQMC, model::Model, Ls = size(lattice(model)); 
-        K = 1+length(neighbors(lattice(model), 1)), 
+        K = let
+            dirs = directions(lattice(model))
+            L2 = dot(dirs[2], dirs[2])
+            i = findfirst(v -> dot(v, v) > L2 + 100eps(L2), dirs)
+            i === nothing ? length(lattice) : i - 1
+        end, 
         capacity = _default_capacity(dqmc),
         obs = LogBinner(ComplexF64(0), capacity=capacity),
         kwargs...
     )
     @assert K > 1
-    dirs = directions(lattice(model))
-
-    # Note: this only works with 2d vecs
-    # Note: longs and trans are epsilon-vectors, i.e. they're the smallest 
-    # step we can take in discrete reciprocal space
-    lvecs = lattice_vectors(lattice(model))
-    uc_vecs = lvecs ./ Ls
-    prefactor = 2pi / dot(lvecs...)
-    rvecs = map(lv -> prefactor * cross([lv; 0], [0, 0, 1])[[1, 2]], uc_vecs)
-
-    # Note: this only works for 2d lattices... maybe
-    longs = []
-    trans = []
-    dir_idxs = []
-    for i in 2:K
-        if -uc_vecs[1] ≈ dirs[i] || uc_vecs[1] ≈ dirs[i]
-            push!(longs, rvecs[1])
-            push!(trans, rvecs[2])
-            push!(dir_idxs, i)
-        elseif -uc_vecs[2] ≈ dirs[i] || uc_vecs[2] ≈ dirs[i]
-            push!(longs, rvecs[2])
-            push!(trans, rvecs[1])
-            push!(dir_idxs, i)
-        else
-            # We kinda need the i, j from R = i*a_1 + j*a_2 here so we can
-            # construct the matching reciprocal vectors here
-            @error("Skipping $(dirs[i]) - not a nearest neighbor")
+
+    # Collect nearest neighbors directions
+    NNs = directions(lattice(model))[2:K]
+
+    # cross(v, e_z) = R90 * v
+    R90 = [0.0 1.0; -1.0 0.0]
+
+    # reciprocal lattice vectors
+    # We use two collections to generate r_1, r_2, -r_1, -r_2, r_i + r_j, ..., 
+    # r_i + 2r_j, ... with i != j
+    # r_i + 2r_j is required for Honeycomb
+    rvecs = let 
+        lvecs = lattice_vectors(lattice(model))
+        uc_vecs = lvecs ./ Ls
+        prefactor = 2pi / (uc_vecs[1][1] * uc_vecs[2][2] - uc_vecs[1][2] * uc_vecs[2][1])
+        rvecs = map(lv -> prefactor * (R90 * lv), uc_vecs)
+        rvecs = vcat(-rvecs, rvecs)
+    end
+    rvecs_ext = vcat([[0, 0]], rvecs, 2rvecs)
+
+    # collect longitudinal (parallel to NN dir) and transversal (perpendicular
+    # to NN dir) reciprocal vectors
+    longs = similar(NNs)
+    trans = similar(NNs)
+    ϵ = 100.0 * eps(dot(NNs[1], NNs[1]) * dot(rvecs_ext[end], rvecs_ext[end]))
+    for k in eachindex(NNs)
+        u = normalize(NNs[k])
+        u_perp = R90 * u
+
+        for i in eachindex(rvecs_ext), j in eachindex(rvecs)
+            i == j && continue # skip elongations
+
+            q = rvecs_ext[i] .+ rvecs[j]
+            abs(q[1]) < ϵ && abs(q[2]) < ϵ && continue # skip zero vectors
+
+            perpendicular = dot(q, u_perp)
+            parallel = dot(q, u)
+
+            # These are directed for now. (There should be symmetry here?) 
+            # If q || u, then parallel = ⟨q, u⟩ has to be positive (same 
+            # direction) and perpendicular = ⟨q, u_perp⟩ has to be close to 0.
+            # We also minimize the length of q
+            if parallel > ϵ && abs(perpendicular) < ϵ 
+                if !(isassigned(longs, k) && dot(longs[k], longs[k]) < dot(q, q))
+                    longs[k] = q
+                end
+            elseif perpendicular > ϵ && abs(parallel) < ϵ
+                if !(isassigned(trans, k) && dot(trans[k], trans[k]) < dot(q, q))
+                    trans[k] = q
+                end
+            end
         end
     end
 
-    #=
-    longs = normalize.(dirs[2:K]) * 1/L
-    trans = map(dirs[2:K]) do v
-        n = [normalize(v)..., 0]
-        u = cross([0,0,1], n)
-        u[1:2] / L
-    end
-    longs .*= 2pi
-    trans .*= 2pi
-    =#
-    li = SuperfluidDensity{EachLocalQuadBySyncedDistance{K}}(
-        dir_idxs, longs, trans
-    )
+    @assert all(i -> isassigned(longs, i), eachindex(NNs)) "Reciprocal vectors not fully defined: longs = $longs; NNs = $NNs"
+    @assert all(i -> isassigned(trans, i), eachindex(NNs)) "Reciprocal vectors not fully defined: trans = $trans; NNs = $NNs"
+
+    li = SuperfluidDensity{EachLocalQuadBySyncedDistance{K}}(2:K, longs, trans)
     # Measurement(dqmc, model, CombinedGreensIterator, li, cc_kernel, obs=obs; kwargs...)
     current_current_susceptibility(dqmc, model, lattice_iterator = li)
 end
@@ -542,53 +563,55 @@ function cc_kernel(mc, model, sites::NTuple{4}, packed_greens::NTuple{4})
     for σ1 in (0, N), σ2 in (0, N)
         s1 = src1 + σ1; t1 = trg1 + σ1
         s2 = src2 + σ2; t2 = trg2 + σ2
+        # raw version
+        # -= from i²
+        # output -= 
+        #     + T[t1, s1] * T[t2, s2] * (dagger(Gll)[t1, s1] * dagger(G00)[t2, s2] + dagger(G0l)[t1, s2] * Gl0[s1, t2]) +
+        #     - T[t1, s1] * T[s2, t2] * (dagger(Gll)[t1, s1] * dagger(G00)[s2, t2] + dagger(G0l)[t1, t2] * Gl0[s1, s2]) +
+        #     - T[s1, t1] * T[t2, s2] * (dagger(Gll)[s1, t1] * dagger(G00)[t2, s2] + dagger(G0l)[s1, s2] * Gl0[t1, t2]) +
+        #     + T[s1, t1] * T[s2, t2] * (dagger(Gll)[s1, t1] * dagger(G00)[s2, t2] + dagger(G0l)[s1, t2] * Gl0[t1, s2])
+        # Replacing constant daggers
+        # output -= 
+        #     + T[t1, s1] * T[t2, s2] * ((I[s1, t1] - Gll[s1, t1]) * (I[s2, t2] - G00[s2, t2]) + dagger(G0l)[t1, s2] * Gl0[s1, t2]) +
+        #     - T[t1, s1] * T[s2, t2] * ((I[s1, t1] - Gll[s1, t1]) * (I[t2, s2] - G00[t2, s2]) + dagger(G0l)[t1, t2] * Gl0[s1, s2]) +
+        #     - T[s1, t1] * T[t2, s2] * ((I[t1, s1] - Gll[t1, s1]) * (I[s2, t2] - G00[s2, t2]) + dagger(G0l)[s1, s2] * Gl0[t1, t2]) +
+        #     + T[s1, t1] * T[s2, t2] * ((I[t1, s1] - Gll[t1, s1]) * (I[t2, s2] - G00[t2, s2]) + dagger(G0l)[s1, t2] * Gl0[t1, s2])
+        # Seperating terms
+        # output -= 
+        #     + T[t1, s1] * (I[s1, t1] - Gll[s1, t1]) * T[t2, s2] * (I[s2, t2] - G00[s2, t2]) +
+        #     - T[t1, s1] * (I[s1, t1] - Gll[s1, t1]) * T[s2, t2] * (I[t2, s2] - G00[t2, s2]) +
+        #     - T[s1, t1] * (I[t1, s1] - Gll[t1, s1]) * T[t2, s2] * (I[s2, t2] - G00[s2, t2]) +
+        #     + T[s1, t1] * (I[t1, s1] - Gll[t1, s1]) * T[s2, t2] * (I[t2, s2] - G00[t2, s2]) +
+        #     + T[t1, s1] * T[t2, s2] * dagger(G0l)[t1, s2] * Gl0[s1, t2] +
+        #     - T[t1, s1] * T[s2, t2] * dagger(G0l)[t1, t2] * Gl0[s1, s2] +
+        #     - T[s1, t1] * T[t2, s2] * dagger(G0l)[s1, s2] * Gl0[t1, t2] +
+        #     + T[s1, t1] * T[s2, t2] * dagger(G0l)[s1, t2] * Gl0[t1, s2]
+        # Combining
+        output -= 
+            (T[t1, s1] * (I[s1, t1] - Gll[s1, t1]) - T[s1, t1] * (I[t1, s1] - Gll[t1, s1])) * 
+            (T[t2, s2] * (I[s2, t2] - G00[s2, t2]) - T[s2, t2] * (I[t2, s2] - G00[t2, s2])) +
+            + T[t1, s1] * T[t2, s2] * dagger(G0l)[t1, s2] * Gl0[s1, t2] +
+            - T[t1, s1] * T[s2, t2] * dagger(G0l)[t1, t2] * Gl0[s1, s2] +
+            - T[s1, t1] * T[t2, s2] * dagger(G0l)[s1, s2] * Gl0[t1, t2] +
+            + T[s1, t1] * T[s2, t2] * dagger(G0l)[s1, t2] * Gl0[t1, s2]
+
+        # previous version:
         # Note: if H is real and Hermitian, T can be pulled out and the I's cancel
         # Note: This matches crstnbr/dqmc if H real, Hermitian
-        # Note: I for G0l and Gl0 auto-cancels
-        output -= (
-                T[t2, s2] * (I[s2, t2] - Gll[s2, t2]) - 
-                T[s2, t2] * (I[t2, s2] - Gll[t2, s2])
-            ) * (
-                T[t1, s1] * (I[s1, t1] - G00[s1, t1]) - 
-                T[s1, t1] * (I[t1, s1] - G00[t1, s1])
-            ) +
-            - T[t2, s2] * T[t1, s1] * G0l[s1, t2] * Gl0[s2, t1] +
-            + T[t2, s2] * T[s1, t1] * G0l[t1, t2] * Gl0[s2, s1] +
-            + T[s2, t2] * T[t1, s1] * G0l[s1, s2] * Gl0[t2, t1] +
-            - T[s2, t2] * T[s1, t1] * G0l[t1, s2] * Gl0[t2, s1]
+        # Note: I for G0l and Gl0 auto-cancels <- this isn't correct because l can be 0
+        # output -= (
+        #         T[t2, s2] * (I[s2, t2] - Gll[s2, t2]) - 
+        #         T[s2, t2] * (I[t2, s2] - Gll[t2, s2])
+        #     ) * (
+        #         T[t1, s1] * (I[s1, t1] - G00[s1, t1]) - 
+        #         T[s1, t1] * (I[t1, s1] - G00[t1, s1])
+        #     ) +
+        #     - T[t2, s2] * T[t1, s1] * G0l[s1, t2] * Gl0[s2, t1] +
+        #     + T[t2, s2] * T[s1, t1] * G0l[t1, t2] * Gl0[s2, s1] +
+        #     + T[s2, t2] * T[t1, s1] * G0l[s1, s2] * Gl0[t2, t1] +
+        #     - T[s2, t2] * T[s1, t1] * G0l[t1, s2] * Gl0[t2, s1]
     end
 
-    # OLD PARTIALLY OUTDATED
-
-    # This should compute 
-    # ⟨j_{trg1-src1}(src1, τ) j_{trg2-src2}(src2, 0)⟩
-    # where (trg-src) picks a direction (e.g. NN directions)
-    # and (src1-src2) is the distance vector that the Fourier transform applies to
-    # From dos Santos: Introduction to Quantum Monte Carlo Simulations
-    # j_{trg-src}(src, τ) = it \sum\sigma (c^\dagger(trg,\sigma, \tau) c(src, \sigma, \tau) - c^\dagger(src, \sigma, \tau) c(trg, \sigma \tau))
-    # where i -> src, i+x -> trg as a generalization
-    # and t is assumed to be hopping matrix element, generalizing to
-    # = i \sum\sigma (T[trg, src] c^\dagger(trg,\sigma, \tau) c(src, \sigma, \tau) - T[src, trg] c^\dagger(src, \sigma, \tau) c(trg, \sigma \tau))
-
-            # Why no I? - delta_0l = 0
-            # T[t1, s1] * T[t2, s2] * (I[s2, t1] - G0l[s2, t1]) * Gl0[s1, t2] -
-            # T[s1, t1] * T[t2, s2] * (I[s2, s1] - G0l[s2, s1]) * Gl0[t1, t2] -
-            # T[t1, s1] * T[s2, t2] * (I[t2, t1] - G0l[t2, t1]) * Gl0[s1, s2] +
-            # T[s1, t1] * T[s2, t2] * (I[t2, s1] - G0l[t2, s1]) * Gl0[t1, s2]
-
-        # Uncompressed Wicks expansion
-        # output += T[t1, s1] * T[t2, s2] *
-        #     ((I[s1, t1] - Gll[s1, t1]) * (I[s2, t2] - G00[s2, t2]) +
-        #     (I[s2, t1] - G0l[s2, t1]) * Gl0[s1, t2])
-        # output -= T[s1, t1] * T[t2, s2] *
-        #     ((I[t1, s1] - Gll[t1, s1]) * (I[s2, t2] - G00[s2, t2]) +
-        #     (I[s2, s1] - G0l[s2, s1]) * Gl0[t1, t2])
-        # output -= T[t1, s1] * T[s2, t2] *
-        #     ((I[s1, t1] - Gll[s1, t1]) * (I[t2, s2] - G00[t2, s2]) +
-        #     (I[t2, t1] - G0l[t2, t1]) * Gl0[s1, s2])
-        # output += T[s1, t1] * T[s2, t2] *
-        #     ((I[t1, s1] - Gll[t1, s1]) * (I[t2, s2] - G00[t2, s2]) +
-        #     (I[t2, s1] - G0l[t2, s1]) * Gl0[t1, s2])
     output
 end
 

src/models/HubbardModel/HubbardModelAttractive.jl

@@ -326,19 +326,29 @@ function cc_kernel(mc, ::HubbardModelAttractive, sites::NTuple{4}, pg::NTuple{4}
     # up-up counts, down-down counts, mixed only on 11s or 22s
     s1 = src1; t1 = trg1
     s2 = src2; t2 = trg2
-    output = -(
-        4.0 * (
-            T[t2, s2] * (I[s2, t2] - Gll[s2, t2]) - 
-            T[s2, t2] * (I[t2, s2] - Gll[t2, s2])
-        ) * (
-            T[t1, s1] * (I[s1, t1] - G00[s1, t1]) - 
-            T[s1, t1] * (I[t1, s1] - G00[t1, s1])
-        ) +
-        - 2.0 * T[t2, s2] * T[t1, s1] * G0l[s1, t2] * Gl0[s2, t1] +
-        + 2.0 * T[t2, s2] * T[s1, t1] * G0l[t1, t2] * Gl0[s2, s1] +
-        + 2.0 * T[s2, t2] * T[t1, s1] * G0l[s1, s2] * Gl0[t2, t1] +
-        - 2.0 * T[s2, t2] * T[s1, t1] * G0l[t1, s2] * Gl0[t2, s1]
-    )
+
+    output = -4.0 * 
+        (T[t1, s1] * (I[s1, t1] - Gll[s1, t1]) - T[s1, t1] * (I[t1, s1] - Gll[t1, s1])) * 
+        (T[t2, s2] * (I[s2, t2] - G00[s2, t2]) - T[s2, t2] * (I[t2, s2] - G00[t2, s2])) +
+        - 2.0 * T[t1, s1] * T[t2, s2] * dagger(G0l)[t1, s2] * Gl0[s1, t2] +
+        + 2.0 * T[t1, s1] * T[s2, t2] * dagger(G0l)[t1, t2] * Gl0[s1, s2] +
+        + 2.0 * T[s1, t1] * T[t2, s2] * dagger(G0l)[s1, s2] * Gl0[t1, t2] +
+        - 2.0 * T[s1, t1] * T[s2, t2] * dagger(G0l)[s1, t2] * Gl0[t1, s2]
+
+
+    # output = -(
+    #     4.0 * (
+    #         T[t2, s2] * (I[s2, t2] - Gll[s2, t2]) - 
+    #         T[s2, t2] * (I[t2, s2] - Gll[t2, s2])
+    #     ) * (
+    #         T[t1, s1] * (I[s1, t1] - G00[s1, t1]) - 
+    #         T[s1, t1] * (I[t1, s1] - G00[t1, s1])
+    #     ) +
+    #     - 2.0 * T[t2, s2] * T[t1, s1] * G0l[s1, t2] * Gl0[s2, t1] +
+    #     + 2.0 * T[t2, s2] * T[s1, t1] * G0l[t1, t2] * Gl0[s2, s1] +
+    #     + 2.0 * T[s2, t2] * T[t1, s1] * G0l[s1, s2] * Gl0[t2, t1] +
+    #     - 2.0 * T[s2, t2] * T[s1, t1] * G0l[t1, s2] * Gl0[t2, s1]
+    # )
     # output = 4.0 * 
     #     (T[s1, t1] * Gll[t1, s1] - T[t1, s1] * Gll[s1, t1]) * 
     #     (T[s2, t2] * G00[t2, s2] - T[t2, s2] * G00[s2, t2]) +

test/parallel.jl

@@ -35,7 +35,7 @@ addprocs(2)
     println("Creating dqmc")
     dqmc = DQMC(
         model, 
-        beta=5.0, recorder=Discarder, scheduler = scheduler,
+        beta=5.0, recorder=Discarder(), scheduler = scheduler,
         thermalization=10000, sweeps=10000
     )