Merge main into z2 branch

README.md

@@ -1,3 +1,4 @@
+[![DOI](https://zenodo.org/badge/338274763.svg)](https://zenodo.org/doi/10.5281/zenodo.10048249)
 # PFFRGSolver.jl <img src=https://github.com/dominikkiese/PFFRGSolver.jl/blob/main/README/logo.png align="right" height="175" width="250">
 
 **P**seudo-**F**ermion **F**unctional **R**enormalization **G**roup **Solver** <br>

src/Flow/bubbles.jl

@@ -118,7 +118,7 @@ function compute_corrs_kat!(
             s_bound_plus  = s_propagator( s_val)
 
             if abs(s_bound_minus) > 1e-8
-                corrs[1, 1, i] = quadgk(s_propagator, Inf, -s_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / s_bound_minus
+                corrs[1, 1, i] = quadgk(s_propagator, -Inf, -s_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / s_bound_minus
             end 
 
             if abs(s_bound_plus) > 1e-8

src/Lattice/model_lib/model_breathing.jl

@@ -82,64 +82,98 @@ end
         ) :: Nothing
 
 Init Heisenberg model on a breathing pyrochlore lattice with broken C3 symmetry by overwriting the respective bonds.
-Here, J[1] = [J1, J2, δ] specifies the breathing nearest-neighbor couplings J1 (for up tetrahedra), 
-J2 (for down tetrahedra) and δ, which quantifies the C3 breaking perturbation.
+Here, J[1] = [J1, J2, δ1, δ2] specifies the breathing nearest-neighbor couplings J1 (for up tetrahedra), 
+J2 (for down tetrahedra) and δ1 (δ2), which quantifies the C3 breaking perturbation on up (down) tetrahedra.
+J[n] (for n >= 2) specifies additional Heisenberg couplings, these are
+* uniformly initialized if J[n] is a single value
+* initialized in ascending bond distance from the origin, if J[n] is an array of length m
 """
 function init_model_pyrochlore_breathing_c3!(
     J :: Vector{Vector{Float64}},
     l :: Lattice
     ) :: Nothing
 
-    @assert l.name in ["pyrochlore"] "Model requires Pyrochlore lattice."
-    @assert length(J[1]) == 3 "Model requires two nearest neighbor couplings and C3 breaking perturbation."
-    @assert length(J) == 1 "Model supports only nearest neighbor couplings."
+    @assert l.name == "pyrochlore" "Model requires Pyrochlore lattice."
+    @assert length(J[1]) == 4 "Model requires two nearest neighbor couplings and two C3 breaking perturbations."
 
     # iterate over sites and add Heisenberg couplings to lattice bonds
     for i in eachindex(l.sites)
-        # find nearest neighbors
-        nbs = get_nbs(1, l.sites[i], l.sites)
-
-        # treat nearest neighbors according to breathing anisotropy
-        for j in nbs
-            # get basis indices 
-            idxs = (l.sites[j].int[4], l.sites[i].int[4])
-
-            # set coupling to J1 (± δ) for up tetrahedra
-            if (l.sites[j].int - l.sites[i].int)[1 : 3] == [0, 0, 0]
-                # add δ on bonds connecting basis sites 1 and 2
-                if idxs == (1, 2) || idxs == (2, 1)
-                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 1, 1)
-                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 2, 2)
-                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 3, 3)
-                # add δ on bonds connecting basis sites 3 and 4
-                elseif idxs == (3, 4) || idxs == (4, 3)
-                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 1, 1)
-                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 2, 2)
-                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 3, 3)
-                # subtract δ for all other bonds
-                else 
-                    add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 1, 1)
-                    add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 2, 2)
-                    add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 3, 3)
-                end 
-            # set coupling to J2 (± δ) for down tetrahedra
+        for n in eachindex(J)
+            # find nearest neighbors
+            nbs = get_nbs(n, l.sites[i], l.sites)
+
+            # treat nearest neighbors according to breathing anisotropy
+            if n == 1
+                for j in nbs
+                    # get basis indices 
+                    idxs = (l.sites[j].int[4], l.sites[i].int[4])
+
+                    # set coupling to J1 (± δ1) for up tetrahedra
+                    if (l.sites[j].int - l.sites[i].int)[1 : 3] == [0, 0, 0]
+                        # add δ on bonds connecting basis sites 1 and 2
+                        if idxs == (1, 2) || idxs == (2, 1)
+                            add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 1, 1)
+                            add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 2, 2)
+                            add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 3, 3)
+                        # add δ on bonds connecting basis sites 3 and 4
+                        elseif idxs == (3, 4) || idxs == (4, 3)
+                            add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 1, 1)
+                            add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 2, 2)
+                            add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 3, 3)
+                        # subtract δ for all other bonds
+                        else 
+                            add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 1, 1)
+                            add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 2, 2)
+                            add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 3, 3)
+                        end 
+                    # set coupling to J2 (± δ2) for down tetrahedra
+                    else
+                        # add δ on bonds connecting basis sites 1 and 2
+                        if idxs == (1, 2) || idxs == (2, 1)
+                            add_bond!(J[1][2] + J[1][4], l.bonds[i, j], 1, 1)
+                            add_bond!(J[1][2] + J[1][4], l.bonds[i, j], 2, 2)
+                            add_bond!(J[1][2] + J[1][4], l.bonds[i, j], 3, 3)
+                        # add δ on bonds connecting basis sites 3 and 4
+                        elseif idxs == (3, 4) || idxs == (4, 3)
+                            add_bond!(J[1][2] + J[1][4], l.bonds[i, j], 1, 1)
+                            add_bond!(J[1][2] + J[1][4], l.bonds[i, j], 2, 2)
+                            add_bond!(J[1][2] + J[1][4], l.bonds[i, j], 3, 3)
+                        # ignore δ for all other bonds
+                        else 
+                            add_bond!(J[1][2] - J[1][4], l.bonds[i, j], 1, 1)
+                            add_bond!(J[1][2] - J[1][4], l.bonds[i, j], 2, 2)
+                            add_bond!(J[1][2] - J[1][4], l.bonds[i, j], 3, 3)
+                        end 
+                    end
+                end
+            # uniform initialization for n-th nearest neighbor, if no further couplings provided
+            elseif length(J[n]) == 1
+                for j in nbs
+                    add_bond!(J[n][1], l.bonds[i, j], 1, 1)
+                    add_bond!(J[n][1], l.bonds[i, j], 2, 2)
+                    add_bond!(J[n][1], l.bonds[i, j], 3, 3)
+                end
+            # initialize symmetry non-equivalent bonds interactions in ascending bond-length order
             else
-                # add δ on bonds connecting basis sites 1 and 2
-                if idxs == (1, 2) || idxs == (2, 1)
-                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 1, 1)
-                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 2, 2)
-                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 3, 3)
-                # add δ on bonds connecting basis sites 3 and 4
-                elseif idxs == (3, 4) || idxs == (4, 3)
-                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 1, 1)
-                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 2, 2)
-                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 3, 3)
-                # ignore δ for all other bonds
-                else 
-                    add_bond!(J[1][2] - J[1][3], l.bonds[i, j], 1, 1)
-                    add_bond!(J[1][2] - J[1][3], l.bonds[i, j], 2, 2)
-                    add_bond!(J[1][2] - J[1][3], l.bonds[i, j], 3, 3)
-                end 
+                # get bond distances of neighbors
+                dist = Int64[get_metric(l.sites[j], l.sites[i], l.uc) for j in nbs]
+
+                # filter out classes of bond distances
+                nbkinds = sort(unique(dist))
+
+                # sanity check
+                @assert length(J[n]) == length(nbkinds) "$(l.name) has $(length(nbkinds)) inequivalent $(n)-th nearest neighbors, but $(length(J[n])) couplings were supplied."
+
+                for nk in eachindex(nbkinds)
+                    # filter out neighbors with dist == nbkinds[nk]
+                    nknbs = nbs[findall(x -> x == nbkinds[nk], dist)]
+
+                    for j in nknbs
+                        add_bond!(J[n][nk], l.bonds[i, j], 1, 1)
+                        add_bond!(J[n][nk], l.bonds[i, j], 2, 2)
+                        add_bond!(J[n][nk], l.bonds[i, j], 3, 3)
+                    end
+                end
             end
         end
     end

src/Launcher/launcher_1l.jl

@@ -120,7 +120,7 @@ function launch_1l!(
         println("   Performing sanity checks and measurements ...")
 
         # terminate if integration becomes unstable
-        if err >= 50.0
+        if err >= 10.0
             println("   Integration has become unstable, terminating solver ...")
             break
         end
@@ -143,7 +143,7 @@ function launch_1l!(
             dΛ = min(dΛ, Λ - Λ_cp)
 
             # terminate if vertex diverges
-            if get_abs_max(a_inter) > max(min(50.0 / Λ, 1000), 10.0)
+            if get_abs_max(a_inter) >= 1000.0
                 println("   Vertex has diverged, terminating solver ...")
                 t, monotone = measure(symmetry, obs_file, cp_file, Λ, dΛ, χ, χ_tol, t, t0, r, m, a_inter, wt, 0.0)
                 break

src/Launcher/launcher_2l.jl

@@ -123,7 +123,7 @@ function launch_2l!(
         println("   Performing sanity checks and measurements ...")
 
         # terminate if integration becomes unstable
-        if err >= 50.0
+        if err >= 10.0
             println("   Integration has become unstable, terminating solver ...")
             break
         end
@@ -146,7 +146,7 @@ function launch_2l!(
             dΛ = min(dΛ, Λ - Λ_cp)
 
             # terminate if vertex diverges
-            if get_abs_max(a_inter) > max(min(50.0 / Λ, 1000), 10.0)
+            if get_abs_max(a_inter) >= 1000.0
                 println("   Vertex has diverged, terminating solver ...")
                 t, monotone = measure(symmetry, obs_file, cp_file, Λ, dΛ, χ, χ_tol, t, t0, r, m, a_inter, wt, 0.0)
                 break

src/Launcher/launcher_ml.jl

@@ -132,7 +132,7 @@ function launch_ml!(
         println("   Performing sanity checks and measurements ...")
 
         # terminate if integration becomes unstable
-        if err >= 50.0
+        if err >= 10.0
             println("   Integration has become unstable, terminating solver ...")
             break
         end
@@ -155,7 +155,7 @@ function launch_ml!(
             dΛ = min(dΛ, Λ - Λ_cp)
 
             # terminate if vertex diverges
-            if get_abs_max(a_inter) > max(min(50.0 / Λ, 1000), 10.0)
+            if get_abs_max(a_inter) >= 1000.0
                 println("   Vertex has diverged, terminating solver ...")
                 t, monotone = measure(symmetry, obs_file, cp_file, Λ, dΛ, χ, χ_tol, t, t0, r, m, a_inter, wt, 0.0)
                 break

src/Observable/fluctuation_lib/fluctuation_su2.jl

@@ -1,4 +1,27 @@
-# auxiliary function to compute equal time correlator for su2 symmetry
+# auxiliary function to compute equal time correlator for su2 symmetry from correlations
+function compute_eq_corr_su2(
+    m :: Mesh,
+    χ :: Vector{Matrix{Float64}}
+    ) :: Float64
+
+    # compute equal-time on-site correlation 
+    eq_corr = 0.0
+    grid    = m.χ
+
+    @turbo for i in 1 : length(grid) - 1 
+        eq_corr += (grid[i + 1] - grid[i]) * (χ[1][1, i + 1] + χ[1][1, i])
+    end 
+
+    # calculate full correlator
+    eq_corr *= 3.0
+
+    # normalization from Fourier transformation
+    eq_corr /= 2.0 * pi
+
+    return eq_corr
+end
+
+# auxiliary function to compute equal time correlator for su2 symmetry from file
 function compute_eq_corr_su2(
     file    :: HDF5.File,
     Λ       :: Float64
@@ -12,7 +35,7 @@ function compute_eq_corr_su2(
     # compute equal-time on-site correlation 
     eq_corr = 0.0
 
-    for i in 1 : length(m) - 1 
+    @turbo for i in 1 : length(m) - 1 
         eq_corr += (m[i + 1] - m[i]) * (χ[1, i + 1] + χ[1, i])
     end 
 

src/Observable/fluctuation_lib/fluctuation_u1_dm.jl

@@ -1,4 +1,25 @@
-# auxiliary function to compute equal time correlator for u1-dm symmetry
+# auxiliary function to compute equal time correlator for u1-dm symmetry from correlations
+function compute_eq_corr_u1_dm(
+    m :: Mesh,
+    χ :: Vector{Matrix{Float64}}
+    ) :: Float64
+
+    # compute equal-time on-site correlation 
+    eq_corr = 0.0
+    grid    = m.χ
+
+    @turbo for i in 1 : length(grid) - 1 
+        eq_corr += 2.0 * (grid[i + 1] - grid[i]) * (χ[1][1, i + 1] + χ[1][1, i])
+        eq_corr += 1.0 * (grid[i + 1] - grid[i]) * (χ[2][1, i + 1] + χ[2][1, i])
+    end 
+
+    # normalization from Fourier transformation
+    eq_corr /= 2.0 * pi 
+
+    return eq_corr
+end
+
+# auxiliary function to compute equal time correlator for u1-dm symmetry from file
 function compute_eq_corr_u1_dm(
     file    :: HDF5.File,
     Λ       :: Float64
@@ -13,7 +34,7 @@ function compute_eq_corr_u1_dm(
     # compute equal-time on-site correlation 
     eq_corr = 0.0
 
-    for i in 1 : length(m) - 1 
+    @turbo for i in 1 : length(m) - 1 
         eq_corr += 2.0 * (m[i + 1] - m[i]) * (χxx[1, i + 1] + χxx[1, i])
         eq_corr += 1.0 * (m[i + 1] - m[i]) * (χzz[1, i + 1] + χzz[1, i])
     end 

src/Observable/fluctuations.jl

@@ -2,7 +2,29 @@
 include("fluctuation_lib/fluctuation_su2.jl")
 include("fluctuation_lib/fluctuation_u1_dm.jl")
 
-# auxiliary function to compute occupation number fluctuations
+# auxiliary function to compute occupation number fluctuations from correlations
+function compute_fluctuations(
+    symmetry :: String,
+    m        :: Mesh,
+    χ        :: Vector{Matrix{Float64}}
+    )        :: Float64
+
+    # init number fluctuations
+    eq_corr = 0.0
+
+    if symmetry == "su2"
+        eq_corr = compute_eq_corr_su2(m, χ)
+    elseif symmetry == "u1-dm"
+        eq_corr = compute_eq_corr_u1_dm(m, χ)
+    end 
+
+    # compute variance of occupation number operator
+    var = 1.0 - 4.0 * eq_corr / 3.0
+
+    return var
+end
+
+# auxiliary function to compute occupation number fluctuations from file
 function compute_fluctuations(
     file    :: HDF5.File,
     Λ       :: Float64
@@ -28,7 +50,7 @@ function compute_fluctuations(
     return var
 end
 
-# auxiliary function to compute flow of occupation number fluctuations
+# auxiliary function to compute flow of occupation number fluctuations from file
 function compute_fluctuations_flow(
     file :: HDF5.File,
     )    :: NTuple{2, Vector{Float64}}