Make multiple walkers work with all the refactoring

Manifest.toml

@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.9.4"
+julia_version = "1.10.4"
 manifest_format = "2.0"
 project_hash = "d50fba1849d6062e02cd5e68cb6158f402e39f01"
 
@@ -127,7 +127,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+0"
+version = "1.1.1+0"
 
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
@@ -200,7 +200,7 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 [[deps.GMP_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
-version = "6.2.1+2"
+version = "6.2.1+6"
 
 [[deps.HostCPUFeatures]]
 deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
@@ -296,9 +296,14 @@ uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
 version = "8.4.0+0"
 
 [[deps.LibGit2]]
-deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[deps.LibGit2_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"]
+uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
+version = "1.6.4+0"
+
 [[deps.LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
@@ -332,7 +337,7 @@ version = "0.12.170"
 [[deps.MPFR_jll]]
 deps = ["Artifacts", "GMP_jll", "Libdl"]
 uuid = "3a97d323-0669-5f0c-9066-3539efd106a3"
-version = "4.1.1+4"
+version = "4.2.0+1"
 
 [[deps.MPI]]
 deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"]
@@ -384,7 +389,7 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.2+0"
+version = "2.28.2+1"
 
 [[deps.MicrosoftMPI_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -397,7 +402,7 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2022.10.11"
+version = "2023.1.10"
 
 [[deps.MuladdMacro]]
 git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab"
@@ -420,7 +425,7 @@ weakdeps = ["Adapt"]
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.21+4"
+version = "0.3.23+4"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
@@ -436,7 +441,7 @@ version = "1.6.3"
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.9.2"
+version = "1.10.0"
 
 [[deps.PkgVersion]]
 deps = ["Pkg"]
@@ -477,7 +482,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[deps.Random]]
-deps = ["SHA", "Serialization"]
+deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.RecipesBase]]
@@ -521,6 +526,7 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+version = "1.10.0"
 
 [[deps.Static]]
 deps = ["IfElse"]
@@ -561,7 +567,7 @@ version = "1.4.2"
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.StrideArraysCore]]
 deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface", "ThreadingUtilities"]
@@ -574,9 +580,9 @@ deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [[deps.SuiteSparse_jll]]
-deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
+deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "5.10.1+6"
+version = "7.2.1+1"
 
 [[deps.TOML]]
 deps = ["Dates"]
@@ -645,12 +651,12 @@ version = "0.21.67"
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.13+0"
+version = "1.2.13+1"
 
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+0"
+version = "5.8.0+1"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -660,4 +666,4 @@ version = "1.52.0+1"
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.4.0+0"
+version = "17.4.0+2"

src/bias/Bias.jl

@@ -33,14 +33,14 @@ The `instance` keyword is used in case of PT-MetaD and multiple walkers to assig
 correct `usebias` to each stream. `has_fp` indicates whether the stream prints to file
 at any point, since only rank 0 should print in case of MPI usage.
 """
-struct Bias{TCV,TS,TB,T}
+struct Bias{TCV,TS,TB,TW,T}
     kind_of_cv::TCV
     smearing::TS
     is_static::Bool
     bias::TB
     biasfile::String
     write_bias_every::Int64
-    kinds_of_weights::Vector{String}
+    kinds_of_weights::TW
     fp::T
 end
 

src/bias/metadynamics.jl

@@ -202,13 +202,13 @@ function metad_from_file(p::ParameterSet, usebias)
     if usebias == ""
         bin_vals = range(p.cvlims[1], p.cvlims[2]; step=p.bin_width)
         values = zero(bin_vals)
-        println("\t>> initialized as zeros")
+        @level1("|  initialized as zeros")
         return bin_vals, values
     else
         values, _ = readdlm(usebias, Float64; header=true)
         bin_vals = range(p.cvlims[1], p.cvlims[2]; step=p.bin_width)
         @assert length(values[:, 2]) == length(bin_vals) "your bias doesn't match parameters"
-        println("\t>> initialized from \"$(usebias)\"")
+        @level1("|  initialized from \"$(usebias)\"")
         return bin_vals, values[:, 2]
     end
 end

src/diracoperators/wilson.jl

@@ -15,6 +15,7 @@ A Wilson Dirac operator with gauge background is created by applying it to a `Ga
 - `TF`: Type of the `Fermionfield` used to store intermediate results when using the 
         Hermitian version of the operator
 - `TG`: Type of the underlying `Gaugefield`
+- `C`: Boolean declaring whether the operator is clover improved or not
 """
 struct WilsonDiracOperator{B,T,C,TF,TG} <: AbstractDiracOperator
     U::TG

src/diracoperators/wilson_eo.jl

@@ -1,4 +1,3 @@
-# TODO
 """
     WilsonEOPreDiracOperator(::Abstractfield, mass; anti_periodic=true)
     WilsonEOPreDiracOperator(
@@ -18,12 +17,15 @@ A Wilson Dirac operator with gauge background is created by applying it to a `Ga
 - `TF`: Type of the `Fermionfield` used to store intermediate results when using the 
         Hermitian version of the operator
 - `TG`: Type of the underlying `Gaugefield`
+- `C`: Boolean declaring whether the operator is clover improved or not
 """
 struct WilsonEOPreDiracOperator{B,T,C,TF,TG,TC,TO} <: AbstractDiracOperator
     U::TG
     temp::TF # temp for storage of intermediate result for DdaggerD operator
     Fμν::TC # Holds the fieldstrength in case we use clover-improved fermions (C==true)
     D_ee_oo::TO
+    D_oo_inv::TO
+    mass::Float64
     κ::Float64
     r::Float64
     csw::Float64
@@ -35,40 +37,44 @@ struct WilsonEOPreDiracOperator{B,T,C,TF,TG,TC,TO} <: AbstractDiracOperator
         κ = 1 / (2mass + 8)
         U = nothing
         C = csw == 0 ? false : true
+        NX, NY, NZ, NT = dims(U)
         temp = even_odd(Fermionfield(f))
-        D_ee_oo = KA.zeros(B(), SMatrix{12,12,Complex{T},144}, dims(U)...)
+        D_ee_oo = KA.zeros(B(), SMatrix{6,6,Complex{T},36}, 2, NX, NY, NZ, NT)
+        D_oo_inv = KA.zeros(B(), SMatrix{6,6,Complex{T},36}, 2, NX, NY, NZ, NT÷2)
 
         TG = Nothing
         TF = typeof(temp)
         TO = typeof(D_ee_oo)
-        return new{B,T,C,TF,TG,TO}(U, temp, D_ee_oo, mass, κ, r, csw, anti_periodic)
+        return new{B,T,C,TF,TG,TO}(U, temp, D_ee_oo, D_oo_inv, mass, κ, r, csw, anti_periodic)
     end
 
     function WilsonEOPreDiracOperator(
         D::WilsonDiracOperator{B,T,C,TF}, U::Gaugefield{B,T}
     ) where {B,T,C,TF}
         check_dims(U, D.temp)
+        NX, NY, NZ, NT = dims(U)
         temp = even_odd(D.temp)
-        D_ee_oo = KA.zeros(B(), SMatrix{12,12,Complex{T},144}, dims(U)...)
+        D_ee_oo = KA.zeros(B(), SMatrix{6,6,Complex{T},36}, 2, NX, NY, NZ, NT)
+        D_oo_inv = KA.zeros(B(), SMatrix{6,6,Complex{T},36}, 2, NX, NY, NZ, NT÷2)
+        calc_ee_oo!(D_ee_oo, D_oo_inv, U, D.mass)
 
         TF_new = typeof(temp)
         TG = typeof(U)
         TO = typeof(D_ee_oo)
         return new{B,T,C,TF_new,TG,TO}(
-            U, temp, D_ee_oo, D.mass, D.κ, D.r, D.csw, D.anti_periodic
+            U, temp, D_ee_oo, D_oo_inv, D.mass, D.κ, D.r, D.csw, D.anti_periodic
         )
     end
 
     function WilsonEOPreDiracOperator(
         D::WilsonEOPreDiracOperator{B,T,C,TF}, U::Gaugefield{B,T}
     ) where {B,T,C,TF}
         check_dims(U, D.temp.parent)
-        D_ee_oo = D.D_ee_oo
-        calc_ee_oo!(D_ee_oo, U)
+        calc_ee_oo!(D.D_ee_oo, D.D_oo_inv, U, D.mass)
         TG = typeof(U)
-        TO = typeof(D_ee_oo)
+        TO = typeof(D.D_ee_oo)
         return new{B,T,C,TF,TG,TO}(
-            U, D.temp, D_ee_oo, D.mass, D.κ, D.r, D.csw, D.anti_periodic
+            U, D.temp, D.D_ee_oo, D.D_oo_inv, D.mass, D.κ, D.r, D.csw, D.anti_periodic
         )
     end
 end
@@ -259,12 +265,10 @@ function LinearAlgebra.mul!(
     U = D.U
     check_dims(ψ_eo, ϕ_eo, U)
     anti = D.anti_periodic
-    D_ee_oo = D.D_ee_oo
+    D_oo_inv = D.D_oo_inv
 
-    # ψₒ = Dₒₑϕₑ
-    mul_oe!(ψ_eo, U, ϕ_eo, D_ee_oo, anti, true, false)
-    # ψₑ = DₑₒDₒₑϕₑ
-    mul_eo!(ψ_eo, U, ψ_eo, anti, false, false)
+    mul_oe!(ψ_eo, U, ϕ_eo, anti, true, Val(false)) # ψₒ = Dₒₑϕₑ
+    mul!(D_oo_inv, ψ_eo) # ψₒ = Dₒₒ⁻¹Dₒₑϕₑ
     axpy!(D_ee_oo, ϕ_eo, ψ_eo) # ψₑ = Dₑₑϕₑ - DₑₒDₒₒ⁻¹Dₒₑϕₑ
     return nothing
 end
@@ -276,12 +280,10 @@ function LinearAlgebra.mul!(
     U = D.parent.U
     check_dims(ψ_eo, ϕ_eo, U)
     anti = D.parent.anti_periodic
-    D_ee_oo = D.D_ee_oo
+    D_oo_inv = D.D_oo_inv
 
-    # ψₒ = Dₒₑϕₑ
-    mul_oe!(ψ_eo, U, ϕ_eo, D_ee_oo, anti, true, true)
-    # ψₑ = DₑₒDₒₑϕₑ
-    mul_eo!(ψ_eo, U, ψ_eo, anti, false, true)
+    mul_oe!(ψ_eo, U, ϕ_eo, anti, true, Val(true)) # ψₒ = Dₒₑϕₑ
+    mul!(D_oo_inv, ψ_eo) # ψₒ = Dₒₒ⁻¹Dₒₑϕₑ
     axpy!(D_ee_oo, ϕ_eo, ψ_eo) # ψₑ = Dₑₑϕₑ - DₑₒDₒₒ⁻¹Dₒₑϕₑ
     return nothing
 end
@@ -331,41 +333,50 @@ function mul_eo!(
     end
 end
 
-function calc_oo_inv!(D_oo_inv, U::Gaugefield{CPU,T}, mass) where {T}
-    @assert size(D_oo_inv) == dims(U)
+function calc_ee_oo!(D_ee_oo, D_oo_inv, U::Gaugefield{CPU,T}, mass) where {T}
+    check_dims(D_ee_oo, U)
     mass_term = Complex{T}(4 + mass)
     fdims = dims(U)
-    NV = volume(U)
+    NV = U.NV
 
     @batch for site in eachindex(U)
-        isodd(site) || continue
         _site = eo_site(site, fdims..., NV)
-
         M = @SMatrix fill(mass_term, 12, 12)
+
+        # TODO: make σμν_spin_mul return two 6x6 matrices instead of one 12x12
         C₁₂ = clover_square(U, 1, 2, site, 1)
         F₁₂ = antihermitian(C₁₂)
-        M += ckron(F₁₂, σμν_spin_mul(ϕ[o_site], Val(1), Val(2)))
+        M += ckron(F₁₂, σμν_spin_mul(ϕ[_site], Val(1), Val(2)))
 
         C₁₃ = clover_square(U, 1, 3, site, 1)
         F₁₃ = antihermitian(C₁₃)
-        M += ckron(F₁₃, σμν_spin_mul(ϕ[o_site], Val(1), Val(3)))
+        M += ckron(F₁₃, σμν_spin_mul(ϕ[_site], Val(1), Val(3)))
 
         C₁₄ = clover_square(U, 1, 4, site, 1)
         F₁₄ = antihermitian(C₁₄)
-        M += ckron(F₁₄, σμν_spin_mul(ϕ[o_site], Val(1), Val(4)))
+        M += ckron(F₁₄, σμν_spin_mul(ϕ[_site], Val(1), Val(4)))
 
         C₂₃ = clover_square(U, 2, 3, site, 1)
         F₂₃ = antihermitian(C₂₃)
-        M += ckron(F₂₃, σμν_spin_mul(ϕ[o_site], Val(2), Val(3)))
+        M += ckron(F₂₃, σμν_spin_mul(ϕ[_site], Val(2), Val(3)))
 
         C₂₄ = clover_square(U, 2, 4, site, 1)
         F₂₄ = antihermitian(C₂₄)
-        M += ckron(F₂₄, σμν_spin_mul(ϕ[o_site], Val(2), Val(4)))
+        M += ckron(F₂₄, σμν_spin_mul(ϕ[_site], Val(2), Val(4)))
 
         C₃₄ = clover_square(U, 3, 4, site, 1)
         F₃₄ = antihermitian(C₃₄)
-        M += ckron(F₃₄, σμν_spin_mul(ϕ[o_site], Val(3), Val(4)))
-        D_oo_inv[o_site] = inv(M)
+        M += ckron(F₃₄, σμν_spin_mul(ϕ[_site], Val(3), Val(4)))
+
+        M1 = SMatrix{6,6,Complex{T},36}(view(M, 1:6, 1:6))
+        M2 = SMatrix{6,6,Complex{T},36}(view(M, 7:12, 7:12))
+        D_ee_oo[1, _site] = M1 
+        D_ee_oo[2, _site] = M2 
+
+        if isodd(site)
+            D_oo_inv[1, _site] = inv(M1)
+            D_oo_inv[2, _site] = inv(M2)
+        end
     end
 end
 

src/gaugefields/fermionfields.jl

@@ -164,7 +164,11 @@ even_odd(f::Fermionfield) = EvenOdd(f)
 
 struct EvenOdd{B,T,A,ND} <: Abstractfield{B,T,A}
     parent::Fermionfield{B,T,A,ND}
-    EvenOdd(f::Fermionfield{B,T,A,ND}) where {B,T,A,ND} = new{B,T,A,ND}(f)
+    function EvenOdd(f::Fermionfield{B,T,A,ND}) where {B,T,A,ND}
+        _, _, _, NT = dims(f)
+        @assert iseven(NT) "Need even time extent for even-odd preconditioning"
+        return new{B,T,A,ND}(f)
+    end
 end
 
 @inline dims(f::EvenOdd) = dims(f.parent)

src/main/Universe.jl

@@ -32,6 +32,8 @@ end
 function Univ(parameters::ParameterSet; use_mpi=false)
     @level1("[ Running MetaQCD.jl version $(PACKAGE_VERSION)\n")
     @level1("┌ Setting Universe...")
+    @level1("|  MPI $(ifelse(use_mpi, "enabled", "disabled"))")
+    @level1("|  Number of instances: $(parameters.numinstances)")
     NX, NY, NZ, NT = parameters.L
     β = parameters.beta
     backend = parameters.backend