diff --git a/Project.toml b/Project.toml
index 8558b400..2747e1c5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -48,13 +48,13 @@ CUDA = "4, 5"
 DelimitedFiles = "1.9"
 JLD2 = "0.4"
 KernelAbstractions = "0.9"
-LoopVectorization = "0.12"
+LoopVectorization = "0.12.170"
 MPI = "0.20"
-Polyester = "0.7.13"
+Polyester = "0.7.14"
 Requires = "1.3"
 StaticArrays = "1.6"
 TimerOutputs = "0.5"
-julia = "1.9"
+julia = "1.9, 1.10"
 
 [extras]
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
diff --git a/src/bias/Bias.jl b/src/bias/Bias.jl
index f2542375..4777a79c 100644
--- a/src/bias/Bias.jl
+++ b/src/bias/Bias.jl
@@ -30,8 +30,7 @@ Container that holds general parameters of bias enhanced sampling, like the kind
 its smearing and filenames/-pointers relevant to the bias. Also holds the specific kind
 of bias (`Metadynamics`, `OPES` or `Parametric` for now). \\
 The `instance` keyword is used in case of PT-MetaD and multiple walkers to assign the
-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.
+correct `usebias` to each stream.
 """
 struct Bias{TCV,TS,TB,TW,T}
     kind_of_cv::TCV
@@ -44,13 +43,14 @@ struct Bias{TCV,TS,TB,TW,T}
     fp::T
 end
 
-function Bias(p::ParameterSet, U; instance=1)
+function Bias(p::ParameterSet, U; use_mpi=false, instance=1)
     @level1("┌ Setting Bias instance $(instance)...")
     kind_of_bias = Unicode.normalize(p.kind_of_bias; casefold=true)
     TCV = get_cvtype_from_parameters(p)
     smearing = StoutSmearing(U, p.numsmears_for_cv, p.rhostout_for_cv)
     is_static = instance == 0 ? true : p.is_static[instance]
     sstr = (is_static || kind_of_bias == "parametric") ? "static" : "dynamic"
+    inum = use_mpi ? myrank+1 : instance
     @level1("|  Type: $(sstr) $(kind_of_bias)")
 
     if kind_of_bias ∈ ["metad", "metadynamics"]
@@ -65,11 +65,11 @@ function Bias(p::ParameterSet, U; instance=1)
 
     @level1("|  CV: $TCV with $(p.numsmears_for_cv) x $(p.rhostout_for_cv) Stout")
 
-    if myrank == 0 && !(bias isa Parametric)
+    if !(bias isa Parametric)
         is_opes = bias isa OPES
         ext = is_opes ? "opes" : "metad"
-        biasfile = p.biasdir * "/stream_$(instance).$(ext)"
-        fp = open(p.measuredir * "/bias_data_$instance.txt", "w")
+        biasfile = myrank==0 ? p.biasdir * "/stream_$(inum).$(ext)" : ""
+        fp = open(p.measuredir * "/bias_data_$inum.txt", "w")
         kinds_of_weights = is_opes ? ["opes"] : p.kinds_of_weights
         str = @sprintf("%-9s\t%-22s", "itrj", "cv")
         print(fp, str)
@@ -80,8 +80,8 @@ function Bias(p::ParameterSet, U; instance=1)
         end
 
         println(fp)
-    elseif myrank == 0 && (bias isa Parametric)
-        fp = open(p.measuredir * "/bias_data_$instance.txt", "w")
+    elseif bias isa Parametric
+        fp = open(p.measuredir * "/bias_data_$inum.txt", "w")
         kinds_of_weights = ["branduardi"]
         str = @sprintf("%-9s\t%-22s\t%-22s\n", "itrj", "cv", "weight_branduardi")
         println(fp, str)
diff --git a/src/bias/weights.jl b/src/bias/weights.jl
index a4e14fdc..40dbea2a 100644
--- a/src/bias/weights.jl
+++ b/src/bias/weights.jl
@@ -14,6 +14,7 @@ function calc_weights(b::Vector{<:Bias}, cv, itrj)
 end
 
 function calc_weights(b::Bias{TCV,TS,TB,T}, cv, itrj) where {TCV,TS,TB,T}
+    b.kinds_of_weights === nothing && return nothing
     str = @sprintf("%-9i\t%+-22.15E", itrj, cv)
     for method in b.kinds_of_weights
         w = calc_weight(b.bias, cv, method)
diff --git a/src/gaugefields/staples.jl b/src/gaugefields/staples.jl
index b5d99733..dfcb8155 100644
--- a/src/gaugefields/staples.jl
+++ b/src/gaugefields/staples.jl
@@ -4,33 +4,33 @@ function staple(::WilsonGaugeAction, U, μ, site)
     return staple_plaq(U, μ, site)
 end
 
-function staple(::SymanzikTreeGaugeAction, U::AbstractArray{SU{3,9,T}}, μ, site) where {T}
-    c₁ = -1//12
+function staple(::SymanzikTreeGaugeAction, U, μ, site)
+    c₁ = float_type(U)(-1/12)
     c₁′ = c₁
     staple_p = staple_plaq(U, μ, site)
     staple_r = staple_rect(U, μ, site)
     return (1 - 8c₁) * staple_p + c₁′ * staple_r
 end
 
-function staple(::SymanzikTadGaugeAction, U::AbstractArray{SU{3,9,T}}, μ, site) where {T}
+function staple(::SymanzikTadGaugeAction, U, μ, site)
     # u0sq = sqrt(plaquette_trace_sum(U))
-    c₁ = -1//12
+    c₁ = float_type(U)(-1/12)
     c₁′ = c₁
     staple_p = staple_plaq(U, μ, site)
     staple_r = staple_rect(U, μ, site)
     return (1 - 8c₁) * staple_p + c₁′ * staple_r
 end
 
-function staple(::IwasakiGaugeAction, U::AbstractArray{SU{3,9,T}}, μ, site) where {T}
-    c₁ = T(-0.331)
+function staple(::IwasakiGaugeAction, U, μ, site)
+    c₁ = float_type(U)(-0.331)
     c₁′ = c₁
     staple_p = staple_plaq(U, μ, site)
     staple_r = staple_rect(U, μ, site)
     return (1 - 8c₁) * staple_p + c₁′ * staple_r
 end
 
-function staple(::DBW2GaugeAction, U::AbstractArray{SU{3,9,T}}, μ, site) where {T}
-    c₁ = T(-1.409)
+function staple(::DBW2GaugeAction, U, μ, site)
+    c₁ = float_type(U)(-1.409)
     c₁′ = c₁
     staple_p = staple_plaq(U, μ, site)
     staple_r = staple_rect(U, μ, site)
diff --git a/src/main/Universe.jl b/src/main/Universe.jl
index 20c300eb..c58b018a 100644
--- a/src/main/Universe.jl
+++ b/src/main/Universe.jl
@@ -69,7 +69,7 @@ function Univ(parameters::ParameterSet; use_mpi=false)
             numinstances = 1
             U = Gaugefield(parameters)
             fermion_actions = init_fermion_actions(parameters, U)
-            bias = Bias(parameters, U)
+            bias = Bias(parameters, U; use_mpi=use_mpi)
         end
     else
         @assert parameters.tempering_enabled == false "tempering can only be enabled with bias"
diff --git a/src/main/runbuild.jl b/src/main/runbuild.jl
index e4c97254..2950d15a 100644
--- a/src/main/runbuild.jl
+++ b/src/main/runbuild.jl
@@ -99,7 +99,7 @@ function build!(parameters, univ, updatemethod, gflow, measurements, measurement
                     friction=0,
                 )
             end
-            @level1("|  Elapsed time:\t$(updatetime) [s]\n-")
+            @level1("|  Elapsed time:\t$(updatetime) [s] @ $(current_time())\n-")
         end
     end
 
@@ -125,7 +125,7 @@ function build!(parameters, univ, updatemethod, gflow, measurements, measurement
                 numaccepts += accepted
             end
 
-            @level1("|  Elapsed time:\t$(updatetime) [s]\n")
+            @level1("|  Elapsed time:\t$(updatetime) [s] @ $(current_time())\n")
             # all procs send their CVs to all other procs and update their copy of the bias
             CVs = MPI.Allgather(U.CV::Float64, comm)
             accepted == true && update_bias!(bias, CVs, itrj, myrank == 0)
diff --git a/src/main/runsim.jl b/src/main/runsim.jl
index aba7036c..54158532 100644
--- a/src/main/runsim.jl
+++ b/src/main/runsim.jl
@@ -187,7 +187,7 @@ function metaqcd!(
                     friction=0,
                 )
             end
-            @level1("|  Elapsed time:\t$(updatetime) [s]\n-")
+            @level1("|  Elapsed time:\t$(updatetime) [s] @ $(current_time())\n-")
         end
     end
 
@@ -214,7 +214,7 @@ function metaqcd!(
             end
 
             print_acceptance_rates(numaccepts, itrj)
-            @level1("|  Elapsed time:\t$(updatetime) [s]")
+            @level1("|  Elapsed time:\t$(updatetime) [s] @ $(current_time())")
 
             save_gaugefield(save_configs, U, itrj)
 
@@ -274,7 +274,7 @@ function metaqcd_PT!(
                     update!(updatemethod, U[i]; bias=NoBias(), metro_test=false, friction=0)
                 end
             end
-            @level1("|  Elapsed time:\t$(updatetime) [s]")
+            @level1("|  Elapsed time:\t$(updatetime) [s] @ $(current_time())")
         end
     end
 
@@ -304,7 +304,7 @@ function metaqcd_PT!(
             end
 
             print_acceptance_rates(numaccepts, itrj)
-            @level1("|  Elapsed time:\t$(updatetime) [s]")
+            @level1("|  Elapsed time:\t$(updatetime) [s] @ $(current_time())")
 
             temper!(U, bias, numaccepts_temper, swap_every, itrj; recalc=!uses_hmc)
 
diff --git a/src/rhmc/AlgRemez.jl b/src/rhmc/AlgRemez.jl
index a65ddabe..7ab226d8 100644
--- a/src/rhmc/AlgRemez.jl
+++ b/src/rhmc/AlgRemez.jl
@@ -7,6 +7,12 @@
 module AlgRemez
 
 using AlgRemez_jll
+using MPI
+
+MPI.Initialized() || MPI.Init()
+const comm = MPI.COMM_WORLD
+const myrank = MPI.Comm_rank(comm)
+const comm_size = MPI.Comm_size(comm)
 
 struct AlgRemezCoeffs{N}
     α0::Float64
@@ -46,13 +52,14 @@ The parameters `y` and `z` must be positive, the approximation to f(x) = x^(-y/z
 the inverse of the approximation to f(x) = x^(y/z).
 """
 function calc_coefficients(y::Int, z::Int, n::Int, lambda_low, lambda_high; precision=42)
-    @assert y > 0 && z > 0 "Inputs y and z need to be positive"
+    myrank==0 && (@assert y > 0 && z > 0 "Inputs y and z need to be positive")
     # Save and reuse already calcuatated approximations
     dirname = pwd() * "/src/rhmc/"
     filename = "approx_$(y)_$(z)_$(n)_$(lambda_low)_$(lambda_high)_$precision.dat"
     filename_err = "error_$(y)_$(z)_$(n)_$(lambda_low)_$(lambda_high)_$precision.dat"
     pathname = dirname * filename
-    if !(filename in readdir(dirname))
+    # If we are using MPI, we only want rank 0 to calculate the coefficients
+    if !(filename in readdir(dirname)) && myrank==0
         try
             run(`$(algremez()) $y $z $n $n $lambda_low $lambda_high $precision`)
         catch _
@@ -64,6 +71,8 @@ function calc_coefficients(y::Int, z::Int, n::Int, lambda_low, lambda_high; prec
         mv("approx.dat", pathname)
         mv("error.dat", dirname * filename_err)
     end
+    MPI.Barrier(comm)
+
     datas = readlines(pathname)
     icount = 3
     αplus0 = parse(Float64, split(datas[icount], "=")[2])