diff --git a/bench_mm.jl b/bench_mm.jl
index d5c00e0..cf21ba5 100755
--- a/bench_mm.jl
+++ b/bench_mm.jl
@@ -4,7 +4,7 @@ function time_mm()
     L = rand(100, 100)
     R = rand(100, 100)
     @time begin
-        @matmul M1[x, σ, α] := sum(β) L[x, β] * M[β, σ, α] 
+        @matmul M1[x, σ, α] := sum(β) L[x, β] * M[β, σ, α]
         @matmul MM[x, σ, y] := sum(α) M1[x, σ, α] * R[α, y]
     end
 end
diff --git a/benchmark/args.jl b/benchmark/args.jl
index b161ed0..516cd3d 100644
--- a/benchmark/args.jl
+++ b/benchmark/args.jl
@@ -9,4 +9,4 @@ T = Float64
 n = 2
 A = rand(T, 2, 2)
 
-my_svd(A, full=true)
+my_svd(A, full = true)
diff --git a/benchmark/cuda_matrix_mul.jl b/benchmark/cuda_matrix_mul.jl
index d11f25c..165cf5d 100644
--- a/benchmark/cuda_matrix_mul.jl
+++ b/benchmark/cuda_matrix_mul.jl
@@ -12,18 +12,18 @@ A = CUDA.CUSPARSE.CuSparseMatrixCSR(Ptr, Ind, Val, (100, 100))
 B = CUDA.rand(Float64, 100, 100)
 C = CUDA.CUSPARSE.CuSparseMatrixCSC(Ptr, Ind, Val, (100, 100))
 
-A*B # no scalar indexing
-CUDA.@allowscalar B*A # scalar indexing
+A * B # no scalar indexing
+CUDA.@allowscalar B * A # scalar indexing
 
-C*B # no scalar indexing
-CUDA.@allowscalar B*C # scalar indexing
+C * B # no scalar indexing
+CUDA.@allowscalar B * C # scalar indexing
 
-A'*B # no scalar indexing
-CUDA.@allowscalar B*A' # scalar indexing
+A' * B # no scalar indexing
+CUDA.@allowscalar B * A' # scalar indexing
 
-transpose(A)*B # no scalar indexing
-CUDA.@allowscalar B*transpose(A) # scalar indexing
+transpose(A) * B # no scalar indexing
+CUDA.@allowscalar B * transpose(A) # scalar indexing
 # problem is when we multiply dense x sparse
 
- D = rand(Float64, (100, 100))
- CUDA.@allowscalar D*A # scalar indexing
+D = rand(Float64, (100, 100))
+CUDA.@allowscalar D * A # scalar indexing
diff --git a/benchmark/memoization_cusparse.jl b/benchmark/memoization_cusparse.jl
index ea2798e..01f51b3 100644
--- a/benchmark/memoization_cusparse.jl
+++ b/benchmark/memoization_cusparse.jl
@@ -5,12 +5,15 @@ using BenchmarkTools
 
 # Functions from constactions_cuda/sparse.jl which are not exported
 
-@memoize Dict function aux_cusparse(::Type{R}, n::Int64) where R <: Real
+@memoize Dict function aux_cusparse(::Type{R}, n::Int64) where {R<:Real}
     println("entering aux_cusparse function")
     CuArray(1:n+1), CUDA.ones(R, n)
 end
 
-@memoize Dict function CUDA.CUSPARSE.CuSparseMatrixCSC(::Type{R}, p::Vector{Int}) where R <: Real
+@memoize Dict function CUDA.CUSPARSE.CuSparseMatrixCSC(
+    ::Type{R},
+    p::Vector{Int},
+) where {R<:Real}
     println("entering cusparse")
     n = length(p)
     cn, co = aux_cusparse(R, n)
@@ -18,7 +21,7 @@ end
 end
 
 
-function CuSparseMatrixCSC_no_memo(::Type{R}, p::Vector{Int}) where R <: Real
+function CuSparseMatrixCSC_no_memo(::Type{R}, p::Vector{Int}) where {R<:Real}
     println("entering no memo")
     n = length(p)
     cn, co = aux_cusparse(R, n)
@@ -39,5 +42,5 @@ p2 = sort(rand(1:5000, 10000000))
 @time F = CUDA.CUSPARSE.CuSparseMatrixCSC(Float64, p2)
 CUDA.memory_status()
 Memoization.empty_all_caches!()
-CUDA.memory_status() 
+CUDA.memory_status()
 # clearing memoization caches doeas not free GPU memory
diff --git a/benchmark/memoization_test.jl b/benchmark/memoization_test.jl
index f08efa6..750d6fb 100644
--- a/benchmark/memoization_test.jl
+++ b/benchmark/memoization_test.jl
@@ -3,22 +3,25 @@ using Memoization
 using CUDA
 
 
-@memoize Dict function example_cuda_array(::Type{R}, size::Int64) where R <: Real
+@memoize Dict function example_cuda_array(::Type{R}, size::Int64) where {R<:Real}
     CUDA.rand(R, (size, size))
 end
 
 
-@memoize Dict function example_array(::Type{R}, size::Int64) where R <: Real
+@memoize Dict function example_array(::Type{R}, size::Int64) where {R<:Real}
     rand(R, size, size)
 end
 
 
-@memoize Dict function aux_cusparse(::Type{R}, n::Int64) where R <: Real
+@memoize Dict function aux_cusparse(::Type{R}, n::Int64) where {R<:Real}
     CuArray(1:n+1), CUDA.ones(R, n)
 end
 
 
-@memoize Dict function CUDA.CUSPARSE.CuSparseMatrixCSC(::Type{R}, p::Vector{Int}) where R <: Real
+@memoize Dict function CUDA.CUSPARSE.CuSparseMatrixCSC(
+    ::Type{R},
+    p::Vector{Int},
+) where {R<:Real}
     n = length(p)
     cn, co = aux_cusparse(R, n)
     CUDA.CUSPARSE.CuSparseMatrixCSC(cn, CuArray(p), co, (maximum(p), n))
diff --git a/benchmark/mulit_gpu.jl b/benchmark/mulit_gpu.jl
index 3fe6ae0..16748c7 100644
--- a/benchmark/mulit_gpu.jl
+++ b/benchmark/mulit_gpu.jl
@@ -1,9 +1,11 @@
 using CUDA
 
-function move_to_CUDA(a::Array{T, N}) where {T, N}
+function move_to_CUDA(a::Array{T,N}) where {T,N}
     buf_a = Mem.alloc(Mem.Unified, sizeof(a))
-    d_a = unsafe_wrap(CuArray{T, N}, convert(CuPtr{T}, buf_a), size(a))
-    finalizer(d_a) do _ Mem.free(buf_a) end
+    d_a = unsafe_wrap(CuArray{T,N}, convert(CuPtr{T}, buf_a), size(a))
+    finalizer(d_a) do _
+        Mem.free(buf_a)
+    end
     copyto!(d_a, a)
     d_a
 end
diff --git a/benchmark/mulit_gpu2.jl b/benchmark/mulit_gpu2.jl
index 5824e4e..bbbcf8c 100644
--- a/benchmark/mulit_gpu2.jl
+++ b/benchmark/mulit_gpu2.jl
@@ -5,6 +5,6 @@ n = 100
 gpus = Int(length(devices()))
 
 a = rand(T, n, n, gpus)
-a_d = cu(a, unified=true)
+a_d = cu(a, unified = true)
 
-a_d 
+a_d
diff --git a/benchmark/psvd.jl b/benchmark/psvd.jl
index ff5de62..2f30108 100644
--- a/benchmark/psvd.jl
+++ b/benchmark/psvd.jl
@@ -11,17 +11,17 @@ cut = 8
 
 mat = rand(100, 100);
 U, S, V = svd(mat);
-S = exp.(collect(0:N-1) * log(4/5));
+S = exp.(collect(0:N-1) * log(4 / 5));
 
 mat = U * Diagonal(S) * V';
 U, S, V = svd(mat);
 
-U, S, V  = U[:, 1:cut], S[1:cut], V[:, 1:cut]
+U, S, V = U[:, 1:cut], S[1:cut], V[:, 1:cut]
 mat1 = U * Diagonal(S) * V'
 println(S[1:cut])
 println(norm(mat - mat1))
 
-Up, Sp, Vp = psvd(mat, rank=2*cut)
+Up, Sp, Vp = psvd(mat, rank = 2 * cut)
 
 mat2 = Up[:, 1:cut] * Diagonal(Sp[1:cut]) * Vp[:, 1:cut]'
 
@@ -29,36 +29,36 @@ println(Sp[1:cut])
 println(Sp[1:cut] - S[1:cut])
 println(norm(mat - mat2))
 
-   # Vp = V
+# Vp = V
 
-   C = mat * Vp
-   println(size(C))
-   Ut, _ = qr(C)
-   Ut = Ut[:, 1:cut]
-   println(size(Ut))
-   C = Ut' * mat
-   Vp, _ = qr(C')
-   Vp = Vp[:, 1:cut]
+C = mat * Vp
+println(size(C))
+Ut, _ = qr(C)
+Ut = Ut[:, 1:cut]
+println(size(Ut))
+C = Ut' * mat
+Vp, _ = qr(C')
+Vp = Vp[:, 1:cut]
 
 
 
-   C = mat * Vp
-   Uf, Sf, Vf = svd(C);
-   Uf, Sf, Vf  = Uf[:, 1:cut], Sf[1:cut], Vf[:, 1:cut]
-   mat3 = Uf * Diagonal(Sf) * Vf' * Vp'
-   println(Sf - S[1:cut])
-   println(norm(mat - mat3))
+C = mat * Vp
+Uf, Sf, Vf = svd(C);
+Uf, Sf, Vf = Uf[:, 1:cut], Sf[1:cut], Vf[:, 1:cut]
+mat3 = Uf * Diagonal(Sf) * Vf' * Vp'
+println(Sf - S[1:cut])
+println(norm(mat - mat3))
 
-   nothing
+nothing
 
 
 iter = 5
 Up, Sp, Vp = [], [], []
-for i in 1:iter
-   Utemp, Stemp, Vtemp = psvd(mat, rank=2*cut)
-   push!(Up, Utemp)
-   push!(Sp, Stemp)
-   push!(Vp, Vtemp)
+for i = 1:iter
+    Utemp, Stemp, Vtemp = psvd(mat, rank = 2 * cut)
+    push!(Up, Utemp)
+    push!(Sp, Stemp)
+    push!(Vp, Vtemp)
 end
 
 Ups = hcat(Up...)
@@ -81,6 +81,6 @@ println(S2)
 mat4 = U2 * Diagonal(S2) * V2'
 
 
-println(norm(mat1-mat2))
-println(norm(mat1-mat3))
-println(norm(mat1-mat4))
\ No newline at end of file
+println(norm(mat1 - mat2))
+println(norm(mat1 - mat3))
+println(norm(mat1 - mat4))
diff --git a/benchmark/sparse_mul.jl b/benchmark/sparse_mul.jl
index ef0489c..dec7c10 100644
--- a/benchmark/sparse_mul.jl
+++ b/benchmark/sparse_mul.jl
@@ -1,7 +1,7 @@
 using CUDA
 using LinearAlgebra
 
-function CUDA.:*(Md::DenseCuMatrix{T}, Mcsc::CUSPARSE.CuSparseMatrixCSC{T}) where T
+function CUDA.:*(Md::DenseCuMatrix{T}, Mcsc::CUSPARSE.CuSparseMatrixCSC{T}) where {T}
     ret = CUDA.zeros(T, size(Mcsc, 1), size(Md, 2))
     CUSPARSE.mm!('N', 'N', one(T), Mcsc, Md, zero(T), ret, 'O')
     ret
diff --git a/benchmark/sparse_mul_bench.jl b/benchmark/sparse_mul_bench.jl
index 81f520a..79deb4b 100644
--- a/benchmark/sparse_mul_bench.jl
+++ b/benchmark/sparse_mul_bench.jl
@@ -2,7 +2,7 @@ using CUDA
 using LinearAlgebra
 using SparseArrays
 
-function dense_x_CSC(Md::DenseCuMatrix{T}, Mcsc::CUSPARSE.CuSparseMatrixCSC{T}) where T
+function dense_x_CSC(Md::DenseCuMatrix{T}, Mcsc::CUSPARSE.CuSparseMatrixCSC{T}) where {T}
     ret = CUDA.zeros(T, size(Md, 1), size(Mcsc, 2))
     CUSPARSE.mm!('N', 'N', one(T), Mcsc, Md, zero(T), ret, 'O')
     ret
diff --git a/benchmark/svd.jl b/benchmark/svd.jl
index de54530..52cc593 100644
--- a/benchmark/svd.jl
+++ b/benchmark/svd.jl
@@ -8,26 +8,26 @@ using FameSVD
 
 
 # C = A * B
-struct MyTensor{T <: Number}
-    A::Array{T, 2}
-    B::Array{T, 2}
+struct MyTensor{T<:Number}
+    A::Array{T,2}
+    B::Array{T,2}
 end
 
 Base.Array(ten::MyTensor) = ten.A * ten.B
 
 # this is for tsvd to work
-Base.eltype(ten::MyTensor{T}) where T = T
+Base.eltype(ten::MyTensor{T}) where {T} = T
 Base.size(ten::MyTensor) = (size(ten.A, 1), size(ten.B, 2))
 Base.size(ten::MyTensor, n::Int) = size(ten)[n]
-Base.adjoint(ten::MyTensor{T}) where T = MyTensor{T}(adjoint(ten.B), adjoint(ten.A))
-Base.:(*)(ten::MyTensor{T}, v::Vector{T}) where T = (ten.A * (ten.B * v))
+Base.adjoint(ten::MyTensor{T}) where {T} = MyTensor{T}(adjoint(ten.B), adjoint(ten.A))
+Base.:(*)(ten::MyTensor{T}, v::Vector{T}) where {T} = (ten.A * (ten.B * v))
 
 # this is for psvd to work
 LinearAlgebra.ishermitian(ten::MyTensor) = ishermitian(ten.A) && ishermitian(ten.B)
 LinearAlgebra.mul!(y, ten::MyTensor, x) = mul!(y, ten.B, ten.A * x)
 
-n = 2 ^ 12
-cut = 2 ^ 6
+n = 2^12
+cut = 2^6
 T = Float64
 
 ten = MyTensor(rand(T, n, n), rand(T, n, n))
@@ -37,7 +37,7 @@ println("tsvd:")
 @time U, Σ, V = tsvd(ten, cut)
 
 println("psvd:")
-@time U, Σ, V = psvd(ten, rank=cut)
+@time U, Σ, V = psvd(ten, rank = cut)
 
 println("Array:")
 println("tsvd:")
@@ -55,7 +55,7 @@ end
 println("psvd:")
 @time begin
     C = Array(ten)
-    U, Σ, V = psvd(C, rank=cut)
+    U, Σ, V = psvd(C, rank = cut)
 end
 
 println("rsvd:")
diff --git a/benchmark/svd2.jl b/benchmark/svd2.jl
index 3ce6cca..84b0e20 100644
--- a/benchmark/svd2.jl
+++ b/benchmark/svd2.jl
@@ -8,31 +8,32 @@ using FameSVD
 
 
 # C = A * B
-struct MyTensor{T <: Number}
-    A::Array{T, 2}
-    B::Array{T, 2}
+struct MyTensor{T<:Number}
+    A::Array{T,2}
+    B::Array{T,2}
 end
 
-struct AMyTensor{T <: Number}
-    A::Array{T, 2}
-    B::Array{T, 2}
+struct AMyTensor{T<:Number}
+    A::Array{T,2}
+    B::Array{T,2}
 end
 
 
 Base.Array(ten::MyTensor) = kron(ten.A, ten.B)
 
 # this is for tsvd to work
-Base.eltype(ten::MyTensor{T}) where T = T
-Base.size(ten::MyTensor) = (size(ten.A, 1) * size(ten.B, 1), size(ten.A, 2) * size(ten.B, 2))
+Base.eltype(ten::MyTensor{T}) where {T} = T
+Base.size(ten::MyTensor) =
+    (size(ten.A, 1) * size(ten.B, 1), size(ten.A, 2) * size(ten.B, 2))
 Base.size(ten::MyTensor, n::Int) = size(ten)[n]
 # Base.adjoint(ten::MyTensor{T}) where T = MyTensor{T}(adjoint(ten.A), adjoint(ten.B))
 
 # Base.:(*)(ten::MyTensor{T}, v::Vector{T}) where T = (kron(ten.A, ten.B) * v)
 
-Base.adjoint(ten::MyTensor{T}) where T = AMyTensor{T}(ten.A, ten.B)
+Base.adjoint(ten::MyTensor{T}) where {T} = AMyTensor{T}(ten.A, ten.B)
 
 
-function Base.:(*)(ten::MyTensor{T}, v::Vector{T}) where T 
+function Base.:(*)(ten::MyTensor{T}, v::Vector{T}) where {T}
     println("M")
     vv = reshape(v, size(ten.A, 2), size(ten.B, 2))
     println(size(vv))
@@ -40,7 +41,7 @@ function Base.:(*)(ten::MyTensor{T}, v::Vector{T}) where T
     reshape(x, size(ten.A, 1) * size(ten.B, 1))
 end
 
-function Base.:(*)(ten::AMyTensor{T}, v::Vector{T}) where T 
+function Base.:(*)(ten::AMyTensor{T}, v::Vector{T}) where {T}
     println("A")
     vv = reshape(v, size(ten.A, 1), size(ten.B, 1))
     println(size(vv))
@@ -54,7 +55,7 @@ end
 LinearAlgebra.ishermitian(ten::MyTensor) = false
 
 
-function LinearAlgebra.mul!(y, ten::MyTensor, v) 
+function LinearAlgebra.mul!(y, ten::MyTensor, v)
     println("K")
     vv = reshape(v, size(ten.A, 2), size(ten.B, 2), :)
     println(size(vv))
@@ -62,7 +63,7 @@ function LinearAlgebra.mul!(y, ten::MyTensor, v)
     y[:, :] = reshape(x, size(ten.A, 1) * size(ten.B, 1), :)
 end
 
-function LinearAlgebra.mul!(y, ten::AMyTensor, v) 
+function LinearAlgebra.mul!(y, ten::AMyTensor, v)
     println("L")
     vv = reshape(v, size(ten.A, 1), size(ten.B, 1), :)
     println(size(vv))
@@ -71,17 +72,17 @@ function LinearAlgebra.mul!(y, ten::AMyTensor, v)
 end
 
 
-n = 2 ^ 2
-cut = 2 ^ 1
+n = 2^2
+cut = 2^1
 T = Float64
 
-ten = MyTensor(rand(T, n+1, n), rand(T, n+2, n-1))
+ten = MyTensor(rand(T, n + 1, n), rand(T, n + 2, n - 1))
 
 println("tsvd:")
 @time U, Σ1, V = tsvd(ten, cut)
 
 println("psvd:")
-@time U, Σ2, V = psvd(ten, rank=cut)
+@time U, Σ2, V = psvd(ten, rank = cut)
 
 println("svd:")
 @time begin
@@ -98,4 +99,4 @@ end
 println(Σ1)
 println(Σ2)
 println(Σ3[1:cut])
-# println(Σ4)
\ No newline at end of file
+# println(Σ4)
diff --git a/docs/make.jl b/docs/make.jl
index cd5c87d..d417f52 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,19 +3,18 @@ using Documenter, SpinGlassTensors
 _pages = [
     "User guide" => "index.md",
     "Matrix Product States and Matrix Product Operations" => "mpo.md",
-    "API Reference" => "api.md"
+    "API Reference" => "api.md",
 ]
 # ============================
 
-format = Documenter.HTML(edit_link = "master",
-                         prettyurls = get(ENV, "CI", nothing) == "true",
-)
+format =
+    Documenter.HTML(edit_link = "master", prettyurls = get(ENV, "CI", nothing) == "true")
 
 # format = Documenter.LaTeX(platform="none")
 
 makedocs(
-    sitename="SpinGlassTensors.jl",
+    sitename = "SpinGlassTensors.jl",
     modules = [SpinGlassTensors],
     pages = _pages,
-    format = format
-    )
\ No newline at end of file
+    format = format,
+)
diff --git a/src/SpinGlassTensors.jl b/src/SpinGlassTensors.jl
index 3895657..f8d54ad 100644
--- a/src/SpinGlassTensors.jl
+++ b/src/SpinGlassTensors.jl
@@ -1,43 +1,43 @@
 module SpinGlassTensors
-    using LinearAlgebra, MKL
-    using TensorOperations, TensorCast
-    using LowRankApprox, TSVD
-    using CUDA, CUDA.CUSPARSE
-    using cuTENSOR
-    using NNlib
-    using Memoization
-    using SparseArrays
-    using DocStringExtensions
-    using Base.Cartesian
+using LinearAlgebra, MKL
+using TensorOperations, TensorCast
+using LowRankApprox, TSVD
+using CUDA, CUDA.CUSPARSE
+using cuTENSOR
+using NNlib
+using Memoization
+using SparseArrays
+using DocStringExtensions
+using Base.Cartesian
 
-    import Base.Prehashed
- #   using SpinGlassNetworks
+import Base.Prehashed
+#   using SpinGlassNetworks
 
-    CUDA.allowscalar(false)
+CUDA.allowscalar(false)
 
-    include("projectors.jl")
-    include("base.jl")
-    include("linear_algebra_ext.jl")
-    include("utils/utils.jl")
-    include("./mps/base.jl")
-    include("./mps/transpose.jl")
-    include("./mps/dot.jl")
-    include("./mps/identity.jl")
-    include("./mps/utils.jl")
-    include("./mps/rand.jl")
-    include("transfer.jl")
-    include("environment.jl")
-    include("utils/memory.jl")
-    include("./mps/canonise.jl")
-    include("variational.jl")
-    include("zipper.jl")
-    include("gauges.jl")
-    include("contractions/sparse.jl")
-    include("contractions/dense.jl")
-    include("contractions/central.jl")
-    include("contractions/diagonal.jl")
-    include("contractions/site.jl")
-    include("contractions/virtual.jl")
+include("projectors.jl")
+include("base.jl")
+include("linear_algebra_ext.jl")
+include("utils/utils.jl")
+include("./mps/base.jl")
+include("./mps/transpose.jl")
+include("./mps/dot.jl")
+include("./mps/identity.jl")
+include("./mps/utils.jl")
+include("./mps/rand.jl")
+include("transfer.jl")
+include("environment.jl")
+include("utils/memory.jl")
+include("./mps/canonise.jl")
+include("variational.jl")
+include("zipper.jl")
+include("gauges.jl")
+include("contractions/sparse.jl")
+include("contractions/dense.jl")
+include("contractions/central.jl")
+include("contractions/diagonal.jl")
+include("contractions/site.jl")
+include("contractions/virtual.jl")
 
 
 end # module
diff --git a/src/base.jl b/src/base.jl
index b2a564f..8a0dc39 100644
--- a/src/base.jl
+++ b/src/base.jl
@@ -1,7 +1,6 @@
 # base.jl: This file defines basic tensor structures to be used with SpinGlassEngine
 
-export
-    Tensor,
+export Tensor,
     SiteTensor,
     VirtualTensor,
     DiagonalTensor,
@@ -9,24 +8,29 @@ export
     CentralOrDiagonal,
     dense_central
 
-abstract type AbstractSparseTensor{T, N} end
+abstract type AbstractSparseTensor{T,N} end
 
-mutable struct SiteTensor{T <: Real, N} <: AbstractSparseTensor{T, N}
+mutable struct SiteTensor{T<:Real,N} <: AbstractSparseTensor{T,N}
     lp::PoolOfProjectors
     loc_exp::AbstractVector{T}
-    projs::NTuple{4, Int} # pl, pt, pr, pb
+    projs::NTuple{4,Int} # pl, pt, pr, pb
     dims::Dims{N}
 
-    function SiteTensor(lp::PoolOfProjectors, loc_exp, projs::NTuple{4, Vector{Int}})
+    function SiteTensor(lp::PoolOfProjectors, loc_exp, projs::NTuple{4,Vector{Int}})
         T = eltype(loc_exp)
         ks = Tuple(add_projector!(lp, p) for p ∈ projs)
         dims = size.(Ref(lp), ks)
-        new{T, 4}(lp, loc_exp, ks, dims)
+        new{T,4}(lp, loc_exp, ks, dims)
     end
 
-    function SiteTensor(lp::PoolOfProjectors, loc_exp, projs::NTuple{4, Int}, dims::NTuple{4, Int})
+    function SiteTensor(
+        lp::PoolOfProjectors,
+        loc_exp,
+        projs::NTuple{4,Int},
+        dims::NTuple{4,Int},
+    )
         T = eltype(loc_exp)
-        new{T, 4}(lp, loc_exp, projs, dims)
+        new{T,4}(lp, loc_exp, projs, dims)
     end
 end
 
@@ -36,7 +40,7 @@ function mpo_transpose(ten::SiteTensor)
     SiteTensor(ten.lp, ten.loc_exp, ten.projs[perm], ten.dims[perm])
 end
 
-mutable struct CentralTensor{T <: Real, N} <: AbstractSparseTensor{T, N}
+mutable struct CentralTensor{T<:Real,N} <: AbstractSparseTensor{T,N}
     e11::AbstractMatrix{T}
     e12::AbstractMatrix{T}
     e21::AbstractMatrix{T}
@@ -48,66 +52,80 @@ mutable struct CentralTensor{T <: Real, N} <: AbstractSparseTensor{T, N}
         @assert s11[1] == s12[1] && s21[1] == s22[1] && s11[2] == s21[2] && s12[2] == s22[2]
         dims = (s11[1] * s21[1], s11[2] * s12[2])
         T = promote_type(eltype.((e11, e12, e21, e22))...)
-        new{T, 2}(e11, e12, e21, e22, dims)
+        new{T,2}(e11, e12, e21, e22, dims)
     end
 end
 
-Base.adjoint(M::CentralTensor{R, 2}) where R <: Real = CentralTensor(M.e11', M.e21', M.e12', M.e22')
+Base.adjoint(M::CentralTensor{R,2}) where {R<:Real} =
+    CentralTensor(M.e11', M.e21', M.e12', M.e22')
 
-mpo_transpose(ten::CentralTensor) = CentralTensor(permutedims.((ten.e11, ten.e21, ten.e12, ten.e22), Ref((2, 1)))...)
+mpo_transpose(ten::CentralTensor) =
+    CentralTensor(permutedims.((ten.e11, ten.e21, ten.e12, ten.e22), Ref((2, 1)))...)
 
-const MatOrCentral{T, N} = Union{AbstractMatrix{T}, CentralTensor{T, N}}
+const MatOrCentral{T,N} = Union{AbstractMatrix{T},CentralTensor{T,N}}
 
 # TODO: to be removed eventually
 function dense_central(ten::CentralTensor)
-    @cast V[(u1, u2), (d1, d2)] := ten.e11[u1, d1] * ten.e21[u2, d1] * ten.e12[u1, d2] * ten.e22[u2, d2]
+    @cast V[(u1, u2), (d1, d2)] :=
+        ten.e11[u1, d1] * ten.e21[u2, d1] * ten.e12[u1, d2] * ten.e22[u2, d2]
     V ./ maximum(V)
 end
 dense_central(ten::AbstractArray) = ten
 
-mutable struct DiagonalTensor{T <: Real, N} <: AbstractSparseTensor{T, N}
-    e1::MatOrCentral{T, N}
-    e2::MatOrCentral{T, N}
+mutable struct DiagonalTensor{T<:Real,N} <: AbstractSparseTensor{T,N}
+    e1::MatOrCentral{T,N}
+    e2::MatOrCentral{T,N}
     dims::Dims{N}
 
     function DiagonalTensor(e1, e2)
         dims = (size(e1, 1) * size(e2, 1), size(e1, 2) * size(e2, 2))
         T = promote_type(eltype.((e1, e2))...)
-        new{T, 2}(e1, e2, dims)
+        new{T,2}(e1, e2, dims)
     end
 end
 
 mpo_transpose(ten::DiagonalTensor) = DiagonalTensor(mpo_transpose.((ten.e2, ten.e1))...)
 
-mutable struct VirtualTensor{T <: Real, N} <: AbstractSparseTensor{T, N}
+mutable struct VirtualTensor{T<:Real,N} <: AbstractSparseTensor{T,N}
     lp::PoolOfProjectors
-    con::MatOrCentral{T, 2}
-    projs::NTuple{6, Int}  # == (p_lb, p_l, p_lt, p_rb, p_r, p_rt)
+    con::MatOrCentral{T,2}
+    projs::NTuple{6,Int}  # == (p_lb, p_l, p_lt, p_rb, p_r, p_rt)
     dims::Dims{N}
 
-    function VirtualTensor(lp::PoolOfProjectors, con, projs::NTuple{6, Vector{Int}})
+    function VirtualTensor(lp::PoolOfProjectors, con, projs::NTuple{6,Vector{Int}})
         T = eltype(con)
         ks = Tuple(add_projector!(lp, p) for p ∈ projs)
-        dims = (length(lp, ks[2]), size(lp, ks[3]) * size(lp, ks[6]),
-                length(lp, ks[5]), size(lp, ks[1]) * size(lp, ks[4]))
-        new{T, 4}(lp, con, ks, dims)
+        dims = (
+            length(lp, ks[2]),
+            size(lp, ks[3]) * size(lp, ks[6]),
+            length(lp, ks[5]),
+            size(lp, ks[1]) * size(lp, ks[4]),
+        )
+        new{T,4}(lp, con, ks, dims)
     end
 
-    function VirtualTensor(lp::PoolOfProjectors, con, projs::NTuple{6, Int}, dims::NTuple{4, Int})
+    function VirtualTensor(
+        lp::PoolOfProjectors,
+        con,
+        projs::NTuple{6,Int},
+        dims::NTuple{4,Int},
+    )
         T = eltype(con)
-        new{T, 4}(lp, con, projs, dims)
+        new{T,4}(lp, con, projs, dims)
     end
 end
 
-mpo_transpose(ten::VirtualTensor) = VirtualTensor(ten.lp, ten.con, ten.projs[[3, 2, 1, 6, 5, 4]], ten.dims[[1, 4, 3, 2]])
-mpo_transpose(ten::AbstractArray{T, 4}) where T = permutedims(ten, (1, 4, 3, 2))
-mpo_transpose(ten::AbstractArray{T, 2}) where T = permutedims(ten, (2, 1))
+mpo_transpose(ten::VirtualTensor) =
+    VirtualTensor(ten.lp, ten.con, ten.projs[[3, 2, 1, 6, 5, 4]], ten.dims[[1, 4, 3, 2]])
+mpo_transpose(ten::AbstractArray{T,4}) where {T} = permutedims(ten, (1, 4, 3, 2))
+mpo_transpose(ten::AbstractArray{T,2}) where {T} = permutedims(ten, (2, 1))
 
-const SparseTensor{T, N} = Union{SiteTensor{T, N}, VirtualTensor{T, N}, CentralTensor{T, N}, DiagonalTensor{T, N}}
-const Tensor{T, N} = Union{AbstractArray{T, N}, SparseTensor{T, N}}
-const CentralOrDiagonal{T, N} = Union{CentralTensor{T, N}, DiagonalTensor{T, N}}
+const SparseTensor{T,N} =
+    Union{SiteTensor{T,N},VirtualTensor{T,N},CentralTensor{T,N},DiagonalTensor{T,N}}
+const Tensor{T,N} = Union{AbstractArray{T,N},SparseTensor{T,N}}
+const CentralOrDiagonal{T,N} = Union{CentralTensor{T,N},DiagonalTensor{T,N}}
 
-Base.eltype(ten::Tensor{T, N}) where {T, N} = T
-Base.ndims(ten::Tensor{T, N}) where {T, N} = N
+Base.eltype(ten::Tensor{T,N}) where {T,N} = T
+Base.ndims(ten::Tensor{T,N}) where {T,N} = N
 Base.size(ten::SparseTensor, n::Int) = ten.dims[n]
 Base.size(ten::SparseTensor) = ten.dims
diff --git a/src/contractions/central.jl b/src/contractions/central.jl
index 411b280..9397449 100644
--- a/src/contractions/central.jl
+++ b/src/contractions/central.jl
@@ -1,22 +1,19 @@
 # contractions with CentralTensor on CPU and CUDA
 
-export
-    contract_tensor3_matrix,
-    contract_matrix_tensor3,
-    update_reduced_env_right
-    # my_batched_mul!
+export contract_tensor3_matrix, contract_matrix_tensor3, update_reduced_env_right
+# my_batched_mul!
 
-function contract_tensor3_matrix(LE::Tensor{R, 3}, M::CentralTensor{R, 2}) where R <: Real
+function contract_tensor3_matrix(LE::Tensor{R,3}, M::CentralTensor{R,2}) where {R<:Real}
     contract_tensor3_central(LE, M.e11, M.e12, M.e21, M.e22)
 end
 
-function contract_matrix_tensor3(M::CentralTensor{R, 2}, RE::Tensor{R, 3}) where R <: Real
+function contract_matrix_tensor3(M::CentralTensor{R,2}, RE::Tensor{R,3}) where {R<:Real}
     contract_tensor3_central(RE, M.e11', M.e21', M.e12', M.e22')
 end
 
-function update_reduced_env_right(RR::Tensor{R, 2}, M::CentralTensor{R, 2}) where R <: Real
+function update_reduced_env_right(RR::Tensor{R,2}, M::CentralTensor{R,2}) where {R<:Real}
     RR = reshape(RR, size(RR, 1), 1, size(RR, 2))
-    dropdims(contract_matrix_tensor3(M, RR), dims=2)
+    dropdims(contract_matrix_tensor3(M, RR), dims = 2)
 end
 
 
@@ -33,44 +30,53 @@ function contract_tensor3_central(LE, e11, e12, e21, e22)
         LE = reshape(reshape(LE, sbt * sl1 * sr1, sl2) * e22, (sbt, sl1, sr1, sr2))  # [tb, l1, r1, r2]
         LE .*= reshape(e11, 1, sl1, sr1, 1)  # [tb, l1, r1, r2] .* [:, l1, r1, :]
         LE .*= reshape(e12, 1, sl1, 1, sr2)  # [tb, l1, r1, r2] .* [:, l1, :, r2]
-        LE = sum(LE, dims=2)
+        LE = sum(LE, dims = 2)
     elseif sr1 <= sr2 && sl2 <= sl1
         LE = permutedims(reshape(LE, (sbt, sl1, sl2)), (1, 3, 2))  # [tb, l2, l1]
         LE = reshape(LE, sbt, sl2, 1, sl1) .* reshape(e11', 1, 1, sr1, sl1)  # [tb, l2, r1, l1]
         LE = reshape(reshape(LE, sbt * sl2 * sr1, sl1) * e12, (sbt, sl2, sr1, sr2))  # [tb, l2, r1, r2]
         LE .*= reshape(e21, 1, sl2, sr1, 1)  # [tb, l2, r1, r2] .* [:, l2, r1, :]
         LE .*= reshape(e22, 1, sl2, 1, sr2)  # [tb, l2, r1, r2] .* [:, l2, :, r2]
-        LE = sum(LE, dims=2)
+        LE = sum(LE, dims = 2)
     elseif sr2 <= sr1 && sl1 <= sl2
         LE = reshape(LE, sbt, sl1, 1, sl2) .* reshape(e22', 1, 1, sr2, sl2)  # [tb, l1, r2, l2]
         LE = reshape(reshape(LE, sbt * sl1 * sr2, sl2) * e21, (sbt, sl1, sr2, sr1))  # [tb, l1, r2, r1]
         LE .*= reshape(e11, 1, sl1, 1, sr1)  # [tb, l1, r2, r1] .* [:, l1, :, r1]
         LE .*= reshape(e12, 1, sl1, sr2, 1)  # [tb, l1, r2, r1] .* [:, l1, r2, :]
-        LE = permutedims(dropdims(sum(LE, dims=2), dims=2), (1, 3, 2))
+        LE = permutedims(dropdims(sum(LE, dims = 2), dims = 2), (1, 3, 2))
     else # sr2 <= sr1 && sl2 <= sl1
         LE = permutedims(reshape(LE, (sbt, sl1, sl2)), (1, 3, 2))  # [tb, l2, l1]
         LE = reshape(LE, sbt, sl2, 1, sl1) .* reshape(e12', 1, 1, sr2, sl1)  # [tb, l2, r2, l1]
         LE = reshape(reshape(LE, sbt * sl2 * sr2, sl1) * e11, (sbt, sl2, sr2, sr1))  # [tb, l2, r2, r1]
         LE .*= reshape(e21, 1, sl2, 1, sr1)  # [tb, l2, r2, r1] .* [:, l2, :, r1]
         LE .*= reshape(e22, 1, sl2, sr2, 1)  # [tb, l2, r2, r1] .* [:, l2, r2, :]
-        LE = permutedims(dropdims(sum(LE, dims=2), dims=2), (1, 3, 2))
+        LE = permutedims(dropdims(sum(LE, dims = 2), dims = 2), (1, 3, 2))
     end
     reshape(LE, sb, st, sr1 * sr2)
 end
 
-function batched_mul!(newLE::Tensor{R, 3}, LE::Tensor{R, 3}, M::AbstractArray{R, 2}) where R <: Real
+function batched_mul!(
+    newLE::Tensor{R,3},
+    LE::Tensor{R,3},
+    M::AbstractArray{R,2},
+) where {R<:Real}
     N1, N2 = size(M)
     new_M = CUDA.CuArray(M)  # TODO: this is a hack to solve problem with types;
     new_M = reshape(new_M, (N1, N2, 1))
     NNlib.batched_mul!(newLE, LE, new_M)
 end
 
-function batched_mul!(newLE::Tensor{R, 3}, LE::Tensor{R, 3}, M::CentralTensor{R, 2}) where R <: Real
+function batched_mul!(
+    newLE::Tensor{R,3},
+    LE::Tensor{R,3},
+    M::CentralTensor{R,2},
+) where {R<:Real}
     sb, _, st = size(LE)
     sl1, sl2, sr1, sr2 = size(M.e11, 1), size(M.e22, 1), size(M.e11, 2), size(M.e22, 2)
     sinter = sb * st * max(sl1 * sl2 * min(sr1, sr2), sr1 * sr2 * min(sl1, sl2))
     if sl1 * sl2 * sr1 * sr2 < sinter
-        @cast E[(l1, l2), (r1, r2)] := M.e11[l1, r1] * M.e21[l2, r1] * M.e12[l1, r2] * M.e22[l2, r2]
+        @cast E[(l1, l2), (r1, r2)] :=
+            M.e11[l1, r1] * M.e21[l2, r1] * M.e12[l1, r2] * M.e22[l2, r2]
         E = reshape(E, (sl1 * sl2, sr1 * sr2, 1))
         NNlib.batched_mul!(newLE, LE, E)
     elseif sr1 <= sr2 && sl1 <= sl2
diff --git a/src/contractions/dense.jl b/src/contractions/dense.jl
index c2aa845..faa6f97 100644
--- a/src/contractions/dense.jl
+++ b/src/contractions/dense.jl
@@ -3,16 +3,19 @@
 #     update_reduced_env_right2
 
 const MatrixOrCuMatrix{R} = Union{
-    CuMatrix{R}, Matrix{R}, Diagonal{R, CuArray{R, 1, Mem.DeviceBuffer}}, Diagonal{R, Vector{R}}
+    CuMatrix{R},
+    Matrix{R},
+    Diagonal{R,CuArray{R,1,Mem.DeviceBuffer}},
+    Diagonal{R,Vector{R}},
 }
 
-function contract_tensor3_matrix(A::Tensor{R, 3}, M::MatrixOrCuMatrix{R}) where R <: Real
+function contract_tensor3_matrix(A::Tensor{R,3}, M::MatrixOrCuMatrix{R}) where {R<:Real}
     sl1, sl2, sl3 = size(A)
     A = reshape(A, sl1 * sl2, sl3)
     reshape(A * M, sl1, sl2, :)
 end
 
-function contract_matrix_tensor3(M::MatrixOrCuMatrix{R}, A::Tensor{R, 3}) where R <: Real
+function contract_matrix_tensor3(M::MatrixOrCuMatrix{R}, A::Tensor{R,3}) where {R<:Real}
     sl1, sl2, sl3 = size(A)
     A = reshape(A, sl1 * sl2, sl3)
     reshape(A * M', sl1, sl2, :)
@@ -25,8 +28,14 @@ end
       |    |
         -- B --
 """
-function update_env_left(LE::S, A::S, M::T, B::S) where {S <: Tensor{R, 3}, T <: Tensor{R, 4}} where R <: Real
-    @tensor order = (ot, α, oc, β, ob) LE[nb, nt, nc] := LE[ob, ot, oc] * A[ot, nt, α] * M[oc, α, nc, β] * B[ob, nb, β] # TODO: split the line
+function update_env_left(
+    LE::S,
+    A::S,
+    M::T,
+    B::S,
+) where {S<:Tensor{R,3},T<:Tensor{R,4}} where {R<:Real}
+    @tensor order = (ot, α, oc, β, ob) LE[nb, nt, nc] :=
+        LE[ob, ot, oc] * A[ot, nt, α] * M[oc, α, nc, β] * B[ob, nb, β] # TODO: split the line
 end
 
 """
@@ -36,7 +45,11 @@ end
       |    |
         -- B --
 """
-function update_env_left(LE::T, A::S, B::S) where {S <: Tensor{R, 3}, T <: Tensor{R, 2}} where R <: Real
+function update_env_left(
+    LE::T,
+    A::S,
+    B::S,
+) where {S<:Tensor{R,3},T<:Tensor{R,2}} where {R<:Real}
     @tensor order = (ot, α, ob) LE[nb, nt] := LE[ob, ot] * A[ot, nt, α] * B[ob, nb, α]
 end
 
@@ -47,7 +60,7 @@ end
       |
 
 """
-function update_env_left(LE::T, A::S) where {S <: Tensor{R, 3}, T <: Tensor{R, 2}} where R <: Real
+function update_env_left(LE::T, A::S) where {S<:Tensor{R,3},T<:Tensor{R,2}} where {R<:Real}
     @tensor A[nb, nt, nc] := LE[nb, ot] * A[ot, nt, nc]
 end
 
@@ -58,8 +71,14 @@ end
          |    |
       -- B --
 """
-function update_env_right(RE::S, A::S, M::T, B::S) where {T <: Tensor{R, 4}, S <: Tensor{R, 3}} where R <: Real
-    @tensor order = (ot, α, oc, β, ob) RE[nb, nt, nc] := RE[ob, ot, oc] * A[nt, ot, α] * M[nc, α, oc, β] * B[nb, ob, β]
+function update_env_right(
+    RE::S,
+    A::S,
+    M::T,
+    B::S,
+) where {T<:Tensor{R,4},S<:Tensor{R,3}} where {R<:Real}
+    @tensor order = (ot, α, oc, β, ob) RE[nb, nt, nc] :=
+        RE[ob, ot, oc] * A[nt, ot, α] * M[nc, α, oc, β] * B[nb, ob, β]
 end
 
 """
@@ -69,7 +88,11 @@ end
          |    |
       -- B --
 """
-function update_env_right(RE::T, A::S, B::S) where {T <: Tensor{R, 2}, S <: Tensor{R, 3}} where R <: Real
+function update_env_right(
+    RE::T,
+    A::S,
+    B::S,
+) where {T<:Tensor{R,2},S<:Tensor{R,3}} where {R<:Real}
     @tensor order = (ot, α, ob) RE[nb, nt] := RE[ob, ot] * A[nt, ot, α] * B[nb, ob, α]
 end
 
@@ -80,7 +103,7 @@ end
               |
 
 """
-function update_env_right(RE::S, C::S) where {S <: Tensor{R, 3}} where R <: Real
+function update_env_right(RE::S, C::S) where {S<:Tensor{R,3}} where {R<:Real}
     @tensor order = (ot, oc) RR[nb, nt] := RE[nb, ot, oc] * C[nt, ot, oc]
 end
 
@@ -90,8 +113,14 @@ end
    |    |    |
      -- B --
 """
-function project_ket_on_bra(LE::S, B::S, M::T, RE::S) where {T <: Tensor{R, 4}, S <: Tensor{R, 3}} where R <: Real
-    @tensor order = (ol, lc, oc, or, rc) A[nl, nr, nc] := LE[ol, nl, lc] * B[ol, or, oc] * M[lc, nc, rc, oc] * RE[or, nr, rc]
+function project_ket_on_bra(
+    LE::S,
+    B::S,
+    M::T,
+    RE::S,
+) where {T<:Tensor{R,4},S<:Tensor{R,3}} where {R<:Real}
+    @tensor order = (ol, lc, oc, or, rc) A[nl, nr, nc] :=
+        LE[ol, nl, lc] * B[ol, or, oc] * M[lc, nc, rc, oc] * RE[or, nr, rc]
 end
 
 """
@@ -99,7 +128,11 @@ end
    |    |    |
      -- B --
 """
-function project_ket_on_bra(LE::T, B::S, RE::T) where {T <: Tensor{R, 2}, S <: Tensor{R, 3}} where R <: Real
+function project_ket_on_bra(
+    LE::T,
+    B::S,
+    RE::T,
+) where {T<:Tensor{R,2},S<:Tensor{R,3}} where {R<:Real}
     @tensor order = (ol, or) A[nl, nr, nc] := LE[ol, nl] * B[ol, or, nc] * RE[or, nr]
 end
 
@@ -107,7 +140,10 @@ end
    |      |
   LE ---- RE --
 """
-function project_ket_on_bra(LE::T, RE::S) where {T <: Tensor{R, 2}, S <: Tensor{R, 3}} where R <: Real
+function project_ket_on_bra(
+    LE::T,
+    RE::S,
+) where {T<:Tensor{R,2},S<:Tensor{R,3}} where {R<:Real}
     @tensor A[nl, nr, nc] := LE[ol, nl] * RE[ol, nr, nc]
 end
 
@@ -118,13 +154,22 @@ end
       |    |
    -- B ---
 """
-function update_reduced_env_right(RE::Tensor{R, 2}, m::Int, M::MpoTensor{R, 4}, B::Tensor{R, 3}) where R <: Real
+function update_reduced_env_right(
+    RE::Tensor{R,2},
+    m::Int,
+    M::MpoTensor{R,4},
+    B::Tensor{R,3},
+) where {R<:Real}
     K = zeros(R, size(M, 2))
     K[m] = one(R)
-    if typeof(RE) <: CuArray K = CuArray(K) end
+    if typeof(RE) <: CuArray
+        K = CuArray(K)
+    end
     K = reshape(K, 1, 1, size(K, 1))
-    for v ∈ M.top K = contract_tensor3_matrix(K, v) end
-    K = dropdims(K, dims=(1, 2))
+    for v ∈ M.top
+        K = contract_tensor3_matrix(K, v)
+    end
+    K = dropdims(K, dims = (1, 2))
 
     for v ∈ reverse(M.bot)
         B = contract_matrix_tensor3(v, B)   # TODO do we ever enter here? in mpo layers that we have now, we don't
@@ -133,19 +178,27 @@ function update_reduced_env_right(RE::Tensor{R, 2}, m::Int, M::MpoTensor{R, 4},
 end
 
 function update_reduced_env_right(
-    K::Tensor{R, 1}, RE::Tensor{R, 2}, M::Tensor{R, 4}, B::Tensor{R, 3}
-) where R <: Real
+    K::Tensor{R,1},
+    RE::Tensor{R,2},
+    M::Tensor{R,4},
+    B::Tensor{R,3},
+) where {R<:Real}
     @tensor order = (d, β, γ, α) RE[x, y] := K[d] * M[y, d, β, γ] * B[x, α, γ] * RE[α, β]
 end
 
-function update_reduced_env_right(RR::S, M0::S) where S <: Tensor{<:Real, 2}
+function update_reduced_env_right(RR::S, M0::S) where {S<:Tensor{<:Real,2}}
     @tensor RR[x, y] := M0[y, z] * RR[x, z]
 end
 
-function contract_tensors43(B::Tensor{R, 4}, A::Tensor{R, 3}) where R <: Real
+function contract_tensors43(B::Tensor{R,4}, A::Tensor{R,3}) where {R<:Real}
     @matmul C[(x, y), (b, a), z] := sum(σ) B[y, z, a, σ] * A[x, b, σ]
 end
 
-function corner_matrix(C::S, M::T, B::S) where {S <: Tensor{R, 3}, T <: Tensor{R, 4}} where R <: Real
-    @tensor order = (rr, mb, mr) Cnew[ll, ml, tt, mt] := M[ml, mt, mr, mb] * B[ll, rr, mb] * C[rr, tt, mr] 
+function corner_matrix(
+    C::S,
+    M::T,
+    B::S,
+) where {S<:Tensor{R,3},T<:Tensor{R,4}} where {R<:Real}
+    @tensor order = (rr, mb, mr) Cnew[ll, ml, tt, mt] :=
+        M[ml, mt, mr, mb] * B[ll, rr, mb] * C[rr, tt, mr]
 end
diff --git a/src/contractions/diagonal.jl b/src/contractions/diagonal.jl
index 280f042..ce41966 100644
--- a/src/contractions/diagonal.jl
+++ b/src/contractions/diagonal.jl
@@ -1,6 +1,6 @@
 # diagonal.jl: contractions with DiagonalTensor on CPU and CUDA
 
-function contract_tensor3_matrix(B::Tensor{R, 3}, C::DiagonalTensor{R}) where R <: Real
+function contract_tensor3_matrix(B::Tensor{R,3}, C::DiagonalTensor{R}) where {R<:Real}
     @cast B[l, (r, s1), s2] := B[l, r, (s1, s2)] (s2 ∈ 1:size(C.e2, 1))
     B = contract_tensor3_matrix(B, C.e2)
     @cast B[l, r, s1, q2] := B[l, (r, s1), q2] (s1 ∈ 1:size(C.e1, 1))
@@ -10,7 +10,7 @@ function contract_tensor3_matrix(B::Tensor{R, 3}, C::DiagonalTensor{R}) where R
     @cast B[l, r, (q2, q1)] := B[l, (r, q2), q1] (q2 ∈ 1:size(C.e2, 2))
 end
 
-function contract_matrix_tensor3(C::DiagonalTensor{R}, B::Tensor{R, 3}) where R <: Real
+function contract_matrix_tensor3(C::DiagonalTensor{R}, B::Tensor{R,3}) where {R<:Real}
     @cast B[l, (r, s2), s1] := B[l, r, (s2, s1)] (s1 ∈ 1:size(C.e1, 2))
     B = contract_matrix_tensor3(C.e1, B)
     @cast B[l, r, s2, q1] := B[l, (r, s2), q1] (s2 ∈ 1:size(C.e2, 2))
diff --git a/src/contractions/site.jl b/src/contractions/site.jl
index 1867c7e..68c495f 100644
--- a/src/contractions/site.jl
+++ b/src/contractions/site.jl
@@ -4,71 +4,126 @@
 # cf. https://discourse.julialang.org/t/correct-implementation-of-cuarrays-slicing-operations/90600
 
 function contract_sparse_with_three(
-    lp, X1::S, X2::S, X3::S, loc_exp::T, k1::Q, k2::Q, k3::Q, kout::Q
-) where {S <: Tensor{R, 3}, T <: Tensor{R, 1}, Q <: Integer} where R <: Real
-s1, s2, _ = size(X1)
-s3, s4, _ = size(X3)
-
-device = typeof(loc_exp) <: CuArray ? :GPU : :CPU
-p1 = get_projector!(lp, k1, device)
-p2 = get_projector!(lp, k2, device)
-p3 = get_projector!(lp, k3, device)
-
-total_memory = 2^32 # TODO add better handling for this; also depending on device
-batch_size = max(Int(floor(total_memory / (8 * (s1 * s2 + s2 * s3 + s3 * s4 + s4 * s1 + min(s1 * s3, s2 * s4))))), 1)
-batch_size = Int(2^floor(log2(batch_size) + 1e-6))
-
-total_size = length(p1)
-batch_size = min(batch_size, total_size)
-
-onGPU = typeof(loc_exp) <: CuArray
-out = onGPU ? CUDA.zeros(R, size(lp, kout), s1, s4) : zeros(R, size(lp, kout), s1, s4)
-
-from = 1
-while from <= total_size
-    to = min(total_size, from + batch_size - 1)
-
-    vp1 = @view p1[from:to]
-    vp2 = @view p2[from:to]
-    vp3 = @view p3[from:to]
-
-    X1p = X1[:, :, vp1]
-    X2p = X2[:, :, vp2]
-    X3p = X3[:, :, vp3]
-
-    if s1 * s3 < s2 * s4
-        Xtmp = batched_mul(X1p, X2p)
-        outp = batched_mul(Xtmp, X3p)
-    else
-        Xtmp = batched_mul(X2p, X3p)
-        outp = batched_mul(X1p, Xtmp)
-    end
+    lp,
+    X1::S,
+    X2::S,
+    X3::S,
+    loc_exp::T,
+    k1::Q,
+    k2::Q,
+    k3::Q,
+    kout::Q,
+) where {S<:Tensor{R,3},T<:Tensor{R,1},Q<:Integer} where {R<:Real}
+    s1, s2, _ = size(X1)
+    s3, s4, _ = size(X3)
+
+    device = typeof(loc_exp) <: CuArray ? :GPU : :CPU
+    p1 = get_projector!(lp, k1, device)
+    p2 = get_projector!(lp, k2, device)
+    p3 = get_projector!(lp, k3, device)
 
-    le = @view loc_exp[from:to]
-    outp .*= reshape(le, 1, 1, :)
-    outpp = reshape(outp, s1 * s4, :)
-    ipr, rf, rt = SparseCSC(R, lp, kout, device; from, to)
-    @inbounds out[rf:rt, :, :] .+= reshape(ipr * outpp', :, s1, s4)
-    from = to + 1
-end
-permutedims(out, (2, 3, 1))
+    total_memory = 2^32 # TODO add better handling for this; also depending on device
+    batch_size = max(
+        Int(
+            floor(
+                total_memory /
+                (8 * (s1 * s2 + s2 * s3 + s3 * s4 + s4 * s1 + min(s1 * s3, s2 * s4))),
+            ),
+        ),
+        1,
+    )
+    batch_size = Int(2^floor(log2(batch_size) + 1e-6))
+
+    total_size = length(p1)
+    batch_size = min(batch_size, total_size)
+
+    onGPU = typeof(loc_exp) <: CuArray
+    out = onGPU ? CUDA.zeros(R, size(lp, kout), s1, s4) : zeros(R, size(lp, kout), s1, s4)
+
+    from = 1
+    while from <= total_size
+        to = min(total_size, from + batch_size - 1)
+
+        vp1 = @view p1[from:to]
+        vp2 = @view p2[from:to]
+        vp3 = @view p3[from:to]
+
+        X1p = X1[:, :, vp1]
+        X2p = X2[:, :, vp2]
+        X3p = X3[:, :, vp3]
+
+        if s1 * s3 < s2 * s4
+            Xtmp = batched_mul(X1p, X2p)
+            outp = batched_mul(Xtmp, X3p)
+        else
+            Xtmp = batched_mul(X2p, X3p)
+            outp = batched_mul(X1p, Xtmp)
+        end
+
+        le = @view loc_exp[from:to]
+        outp .*= reshape(le, 1, 1, :)
+        outpp = reshape(outp, s1 * s4, :)
+        ipr, rf, rt = SparseCSC(R, lp, kout, device; from, to)
+        @inbounds out[rf:rt, :, :] .+= reshape(ipr * outpp', :, s1, s4)
+        from = to + 1
+    end
+    permutedims(out, (2, 3, 1))
 end
 
-function update_env_left(LE::S, A::S, M::T, B::S) where {S <: Tensor{R, 3}, T <: SiteTensor{R, 4}} where R <: Real
-    contract_sparse_with_three(M.lp, permutedims(B, (2, 1, 3)), LE, A, M.loc_exp, M.projs[[4, 1, 2, 3]]...)
+function update_env_left(
+    LE::S,
+    A::S,
+    M::T,
+    B::S,
+) where {S<:Tensor{R,3},T<:SiteTensor{R,4}} where {R<:Real}
+    contract_sparse_with_three(
+        M.lp,
+        permutedims(B, (2, 1, 3)),
+        LE,
+        A,
+        M.loc_exp,
+        M.projs[[4, 1, 2, 3]]...,
+    )
 end
 
-function update_env_right(RE::S, A::S, M::SiteTensor{R, 4}, B::S) where {S <: Tensor{R, 3}} where R <: Real
-    contract_sparse_with_three(M.lp, B, RE, permutedims(A, (2, 1, 3)), M.loc_exp, M.projs[[4, 3, 2, 1]]...)
+function update_env_right(
+    RE::S,
+    A::S,
+    M::SiteTensor{R,4},
+    B::S,
+) where {S<:Tensor{R,3}} where {R<:Real}
+    contract_sparse_with_three(
+        M.lp,
+        B,
+        RE,
+        permutedims(A, (2, 1, 3)),
+        M.loc_exp,
+        M.projs[[4, 3, 2, 1]]...,
+    )
 end
 
-function project_ket_on_bra(LE::S, B::S, M::SiteTensor{R, 4}, RE::S) where {S <: Tensor{R, 3}} where R <: Real
-    contract_sparse_with_three(M.lp, permutedims(LE, (2, 1, 3)), B, RE, M.loc_exp, M.projs[[1, 4, 3, 2]]...)
+function project_ket_on_bra(
+    LE::S,
+    B::S,
+    M::SiteTensor{R,4},
+    RE::S,
+) where {S<:Tensor{R,3}} where {R<:Real}
+    contract_sparse_with_three(
+        M.lp,
+        permutedims(LE, (2, 1, 3)),
+        B,
+        RE,
+        M.loc_exp,
+        M.projs[[1, 4, 3, 2]]...,
+    )
 end
 
 function update_reduced_env_right(
-    K::Tensor{R, 1}, RE::Tensor{R, 2}, M::SiteTensor{R, 4}, B::Tensor{R, 3}
-) where R <: Real
+    K::Tensor{R,1},
+    RE::Tensor{R,2},
+    M::SiteTensor{R,4},
+    B::Tensor{R,3},
+) where {R<:Real}
     device = typeof(M.loc_exp) <: CuArray ? :GPU : :CPU
     s1, s2, _ = size(B)
 
@@ -79,7 +134,9 @@ function update_reduced_env_right(
     batch_size = max(Int(floor(total_memory / (8 * (s1 * s2 + s1 + s2 + 1)))), 1)
     batch_size = Int(2^floor(log2(batch_size) + 1e-6))
 
-    out = typeof(M.loc_exp) <: CuArray ? CUDA.zeros(R, size(M.lp, k1), s1) : zeros(R, size(M.lp, k1), s1)
+    out =
+        typeof(M.loc_exp) <: CuArray ? CUDA.zeros(R, size(M.lp, k1), s1) :
+        zeros(R, size(M.lp, k1), s1)
     RE = reshape(RE, size(RE, 1), 1, size(RE, 2))
 
     from = 1
@@ -95,7 +152,7 @@ function update_reduced_env_right(
         @inbounds Bp = B[:, :, vp4]
         le = @view M.loc_exp[from:to]
 
-        outp = dropdims(Bp ⊠ REp, dims=2)
+        outp = dropdims(Bp ⊠ REp, dims = 2)
         outp .*= reshape(le .* Kp, 1, :)
 
         ipr, rf, rt = SparseCSC(R, M.lp, k1, device; from, to)
@@ -105,7 +162,7 @@ function update_reduced_env_right(
     permutedims(out, (2, 1))
 end
 
-function contract_tensors43(M::SiteTensor{R, 4}, B::Tensor{R, 3}) where R <: Real
+function contract_tensors43(M::SiteTensor{R,4}, B::Tensor{R,3}) where {R<:Real}
     device = typeof(M.loc_exp) <: CuArray ? :GPU : :CPU
     p4 = get_projector!(M.lp, M.projs[4], device)
     sb1, sb2, _ = size(B)
@@ -118,7 +175,11 @@ function contract_tensors43(M::SiteTensor{R, 4}, B::Tensor{R, 3}) where R <: Rea
     reshape(out, sb1 * sm1, sb2 * sm3, sm2)
 end
 
-function corner_matrix(C::S, M::T, B::S) where {S <: Tensor{R, 3}, T <: SiteTensor{R, 4}} where R <: Real
+function corner_matrix(
+    C::S,
+    M::T,
+    B::S,
+) where {S<:Tensor{R,3},T<:SiteTensor{R,4}} where {R<:Real}
     device = typeof(M.loc_exp) <: CuArray ? :GPU : :CPU
     projs = [get_projector!(M.lp, x, device) for x in M.projs]
     @inbounds Bp = B[:, :, projs[4]]
@@ -128,7 +189,7 @@ function corner_matrix(C::S, M::T, B::S) where {S <: Tensor{R, 3}, T <: SiteTens
     @cast outp[(x, y), z] := outp[x, y, z]
     sm1, sm2 = maximum(projs[1]), maximum(projs[2])
     @inbounds p12 = projs[1] .+ (projs[2] .- 1) .* sm1
-    ip12 = SparseCSC(R, p12; mp=sm1 * sm2)
+    ip12 = SparseCSC(R, p12; mp = sm1 * sm2)
     out = reshape(ip12 * outp', sm1, maximum(projs[2]), size(B, 1), size(C, 2))
     permutedims(out, (3, 1, 4, 2))
 end
diff --git a/src/contractions/sparse.jl b/src/contractions/sparse.jl
index 6375ea3..a20af83 100644
--- a/src/contractions/sparse.jl
+++ b/src/contractions/sparse.jl
@@ -10,7 +10,7 @@ end
 =#
 
 # TODO shouldn't we have CSR format instead?
-function SparseCSC(::Type{R}, p::CuArray{Int64, 1}; mp=nothing) where R <: Real
+function SparseCSC(::Type{R}, p::CuArray{Int64,1}; mp = nothing) where {R<:Real}
     n = length(p)
     if mp == nothing
         mp = maximum(p)
@@ -20,7 +20,7 @@ function SparseCSC(::Type{R}, p::CuArray{Int64, 1}; mp=nothing) where R <: Real
     CuSparseMatrixCSC(cn, p, co, (mp, n))
 end
 
-function SparseCSC(::Type{R}, p::Vector{Int64}; mp=nothing) where R <: Real
+function SparseCSC(::Type{R}, p::Vector{Int64}; mp = nothing) where {R<:Real}
     n = length(p)
     if mp == nothing
         mp = maximum(p)
@@ -31,13 +31,13 @@ function SparseCSC(::Type{R}, p::Vector{Int64}; mp=nothing) where R <: Real
 end
 
 @memoize Dict function SparseCSC(
-    ::Type{T}, 
-    lp::PoolOfProjectors, 
-    k1::R, 
-    k2::R, 
-    k3::R, 
-    device::Symbol
-    ) where {T <: Real, R <: Int}
+    ::Type{T},
+    lp::PoolOfProjectors,
+    k1::R,
+    k2::R,
+    k3::R,
+    device::Symbol,
+) where {T<:Real,R<:Int}
     p1 = get_projector!(lp, k1) #, device)
     p2 = get_projector!(lp, k2) #, device)
     p3 = get_projector!(lp, k3) #, device)
@@ -47,17 +47,17 @@ end
     if device == :GPU
         p = CuArray(p)
     end
-    SparseCSC(T, p; mp=s1 * s2 * s3)
+    SparseCSC(T, p; mp = s1 * s2 * s3)
 end
 
 @memoize Dict function SparseCSC(
-    ::Type{R}, 
-    lp::PoolOfProjectors, 
-    k::Int, 
-    device::Symbol; 
-    from::Int=1, 
-    to::Int=length(lp, k)
-    ) where R <: Real
+    ::Type{R},
+    lp::PoolOfProjectors,
+    k::Int,
+    device::Symbol;
+    from::Int = 1,
+    to::Int = length(lp, k),
+) where {R<:Real}
     p = get_projector!(lp, k)
     pp = @view p[from:to]
     rf = minimum(pp)
diff --git a/src/contractions/virtual.jl b/src/contractions/virtual.jl
index b042b49..ebaa4a3 100644
--- a/src/contractions/virtual.jl
+++ b/src/contractions/virtual.jl
@@ -28,8 +28,8 @@ function proj_2step_12(lp, (k1, k2), k3, device)
 
     if device == :CPU
         p12 = CuArray(p12)
-        p1  = CuArray(p1)
-        p2  = CuArray(p2)
+        p1 = CuArray(p1)
+        p2 = CuArray(p2)
     end
 
     pf1 = p12 .+ s12 .* (p3 .- 1)
@@ -53,8 +53,8 @@ function proj_2step_23(lp, k1, (k2, k3), device)
 
     if device == :CPU
         p23 = CuArray(p23)
-        p2  = CuArray(p2)
-        p3  = CuArray(p3)
+        p2 = CuArray(p2)
+        p3 = CuArray(p3)
     end
 
     pf1 = p1 .+ s1 .* (p23 .- 1)
@@ -63,7 +63,7 @@ function proj_2step_23(lp, k1, (k2, k3), device)
     pf1, pf2, s23
 end
 
-function merge_projectors_inter(lp, p1, p2, p3, onGPU; order="1_23")
+function merge_projectors_inter(lp, p1, p2, p3, onGPU; order = "1_23")
     s1 = size(lp, p1)
     s2 = size(lp, p2)
     device = onGPU ? :GPU : :CPU
@@ -84,7 +84,12 @@ function merge_projectors_inter(lp, p1, p2, p3, onGPU; order="1_23")
     p123, p2_3, s23
 end
 
-function update_env_left(LE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <: Tensor{R, 3}} where R <: Real
+function update_env_left(
+    LE::S,
+    A::S,
+    M::VirtualTensor{R,4},
+    B::S,
+) where {S<:Tensor{R,3}} where {R<:Real}
     p_lb, p_lc, p_lt, p_rb, p_rc, p_rt = M.projs
     slb, srb = size(B, 1), size(B, 2)
     slt, srt = size(A, 1), size(A, 2)
@@ -100,8 +105,10 @@ function update_env_left(LE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
     Lout = alloc_zeros(R, onGPU, (srb, srt, src))
 
     if slpb * srpt >= slpt * srpb
-        pl_b_ct, pl_c_t, slpct = merge_projectors_inter(M.lp, p_lb, p_lc, p_lt, onGPU; order="1_23")
-        pr_bc_t, pr_b_c, srpbc = merge_projectors_inter(M.lp, p_rt, p_rb, p_rc, onGPU; order="23_1")
+        pl_b_ct, pl_c_t, slpct =
+            merge_projectors_inter(M.lp, p_lb, p_lc, p_lt, onGPU; order = "1_23")
+        pr_bc_t, pr_b_c, srpbc =
+            merge_projectors_inter(M.lp, p_rt, p_rb, p_rc, onGPU; order = "23_1")
 
         B2 = permutedims(B, (1, 3, 2, 4))  # [lb, lpb, rb, rpb]
         B2 = reshape(B2, (slb * slpb, srb * srpb))  # [(lb, lpb), (rb, rpb)]
@@ -112,7 +119,7 @@ function update_env_left(LE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
         tmp5 = alloc_undef(R, onGPU, (srb * srpb, srpc, slpt))
         tmp8 = alloc_undef(R, onGPU, (srb * srpbc, srpt))
 
-        for ilt ∈ 1 : slt
+        for ilt ∈ 1:slt
             tmp1[:, pl_b_ct] = (@view LE[:, ilt, :])  # [lb, (lpb, lpct)]
             mul!(tmp2, B2', reshape(tmp1, (slb * slpb, slpct)))  # [(rb, rpb), lpct]
             tmp3[:, pl_c_t] = tmp2  # [(rb, rpb), (lpc, lpt)]
@@ -120,15 +127,17 @@ function update_env_left(LE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
             batched_mul!(tmp5, tmp4, M.con)
             tmp6 = reshape(tmp5, (srb, srpb * srpc, slpt))  # [rb, (rpb, rpc), lpt]
             tmp7 = reshape(tmp6[:, pr_b_c, :], (srb * srpbc, slpt))  # [(rb, rpbc), lpt]
-            for irt ∈ 1 : srt
+            for irt ∈ 1:srt
                 mul!(tmp8, tmp7, (@view A[ilt, irt, :, :]))
                 tmp9 = reshape(tmp8, (srb, srpbc * srpt))
                 Lout[:, irt, :] .+= tmp9[:, pr_bc_t]  # [rb, rc]
             end
         end
     else
-        pl_t_cb, pl_c_b, slpcb = merge_projectors_inter(M.lp, p_lt, p_lc, p_lb, onGPU; order="1_23")
-        pr_tc_b, pr_t_c, srptc = merge_projectors_inter(M.lp, p_rb, p_rt, p_rc, onGPU; order="23_1")
+        pl_t_cb, pl_c_b, slpcb =
+            merge_projectors_inter(M.lp, p_lt, p_lc, p_lb, onGPU; order = "1_23")
+        pr_tc_b, pr_t_c, srptc =
+            merge_projectors_inter(M.lp, p_rb, p_rt, p_rc, onGPU; order = "23_1")
 
         A2 = permutedims(A, (1, 3, 2, 4))  # [lt, lpt, rt, rpt]
         A2 = reshape(A2, (slt * slpt, srt * srpt))  # [(lt, lpt), (rt, rpt)]
@@ -139,7 +148,7 @@ function update_env_left(LE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
         tmp5 = alloc_undef(R, onGPU, (srt * srpt, srpc, slpb))
         tmp8 = alloc_zeros(R, onGPU, (srt * srptc, srpb))
 
-        for ilb ∈ 1 : slb
+        for ilb ∈ 1:slb
             tmp1[:, pl_t_cb] = (@view LE[ilb, :, :])  # [lt, (lpt, lpcb)]
             mul!(tmp2, A2', reshape(tmp1, (slt * slpt, slpcb)))  # [(rt, rpt), lpcb]
             tmp3[:, pl_c_b] = tmp2
@@ -147,7 +156,7 @@ function update_env_left(LE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
             batched_mul!(tmp5, tmp4, M.con)  # [(rt, rpt), lpb, rpc]
             tmp6 = reshape(tmp5, (srt, srpt * srpc, slpb))  # [(rt, rpt * rpc), lcb]
             tmp7 = reshape(tmp6[:, pr_t_c, :], (srt * srptc, slpb))  # [(rt, rptc), lpb]
-            for irb ∈ 1 : srb
+            for irb ∈ 1:srb
                 mul!(tmp8, tmp7, (@view B[ilb, irb, :, :]))
                 tmp9 = reshape(tmp8, (srt, srptc * srpb))
                 Lout[irb, :, :] .+= tmp9[:, pr_tc_b]  # [rt, rc]
@@ -158,7 +167,12 @@ function update_env_left(LE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
 end
 
 
-function project_ket_on_bra(LE::S, B::S, M::VirtualTensor{R, 4}, RE::S) where {S <: Tensor{R, 3}} where R <: Real
+function project_ket_on_bra(
+    LE::S,
+    B::S,
+    M::VirtualTensor{R,4},
+    RE::S,
+) where {S<:Tensor{R,3}} where {R<:Real}
     p_lb, p_lc, p_lt, p_rb, p_rc, p_rt = M.projs
     slb, slt = size(LE, 1), size(LE, 2)
     srb, srt = size(RE, 1), size(RE, 2)
@@ -173,15 +187,17 @@ function project_ket_on_bra(LE::S, B::S, M::VirtualTensor{R, 4}, RE::S) where {S
     LR = alloc_zeros(R, onGPU, (slt, srt, slpt, srpt))
 
     if slpb >= srpb
-        pl_b_ct, pl_c_t, slpct = merge_projectors_inter(M.lp, p_lb, p_lc, p_lt, onGPU; order="1_23")
-        pr_bc_t, pr_b_c, srpbc = merge_projectors_inter(M.lp, p_rt, p_rb, p_rc, onGPU; order="23_1")
+        pl_b_ct, pl_c_t, slpct =
+            merge_projectors_inter(M.lp, p_lb, p_lc, p_lt, onGPU; order = "1_23")
+        pr_bc_t, pr_b_c, srpbc =
+            merge_projectors_inter(M.lp, p_rt, p_rb, p_rc, onGPU; order = "23_1")
 
         tmp1 = alloc_zeros(R, onGPU, (slb, slpb * slpct))
         tmp2 = alloc_undef(R, onGPU, (srb * srpb, slpct))
         tmp3 = alloc_zeros(R, onGPU, (srb * srpb, slpc * slpt))
         tmp5 = alloc_undef(R, onGPU, (srb * srpb, srpc, slpt))
         tmp8 = alloc_zeros(R, onGPU, (srb, srpbc * srpt))
-        for ilt ∈ 1 : slt
+        for ilt ∈ 1:slt
             tmp1[:, pl_b_ct] = (@view LE[:, ilt, :])  # [lb, (lpb, lpct)]
             mul!(tmp2, B2', reshape(tmp1, (slb * slpb, slpct)))  # [(rb, rpb), lpct]
             tmp3[:, pl_c_t] = tmp2  # [(rb, rpb), (lpc, lpt)]
@@ -189,21 +205,23 @@ function project_ket_on_bra(LE::S, B::S, M::VirtualTensor{R, 4}, RE::S) where {S
             batched_mul!(tmp5, tmp4, M.con)  # [(rb, rpb), rpc, lpt]
             tmp6 = reshape(tmp5, (srb, srpb * srpc, slpt))  # [rb, (rpb, rpc), lpt]
             tmp7 = reshape(tmp6[:, pr_b_c, :], (srb * srpbc, slpt))  # [(rb, rpbc), lpt]
-            for irt ∈ 1 : srt
+            for irt ∈ 1:srt
                 tmp8[:, pr_bc_t] = (@view RE[:, irt, :])  # [rb, (rpbc, rpt)]
                 LR[ilt, irt, :, :] = tmp7' * reshape(tmp8, (srb * srpbc, srpt))  # [lpt, rpt]
             end
         end
     else
-        pr_b_ct, pr_c_t, srpct = merge_projectors_inter(M.lp, p_rb, p_rc, p_rt, onGPU; order="1_23")
-        pl_bc_t, pl_b_c, slpbc = merge_projectors_inter(M.lp, p_lt, p_lb, p_lc, onGPU; order="23_1")
+        pr_b_ct, pr_c_t, srpct =
+            merge_projectors_inter(M.lp, p_rb, p_rc, p_rt, onGPU; order = "1_23")
+        pl_bc_t, pl_b_c, slpbc =
+            merge_projectors_inter(M.lp, p_lt, p_lb, p_lc, onGPU; order = "23_1")
 
         tmp1 = alloc_zeros(R, onGPU, (srb, srpb * srpct))
         tmp2 = alloc_undef(R, onGPU, (slb * slpb, srpct))
         tmp3 = alloc_zeros(R, onGPU, (slb * slpb, srpc * srpt))
         tmp5 = alloc_undef(R, onGPU, (slb * slpb, slpc, srpt))
         tmp8 = alloc_zeros(R, onGPU, (slb, slpbc * slpt))
-        for irt ∈ 1 : srt
+        for irt ∈ 1:srt
             tmp1[:, pr_b_ct] = (@view RE[:, irt, :])  # [rb, (rpb, rpct)]
             mul!(tmp2, B2, reshape(tmp1, (srb * srpb, srpct)))  # [(lb, lpb), rpct]
             tmp3[:, pr_c_t] = tmp2  # [(lb, lpb), (rpc, rpt)]
@@ -211,7 +229,7 @@ function project_ket_on_bra(LE::S, B::S, M::VirtualTensor{R, 4}, RE::S) where {S
             batched_mul!(tmp5, tmp4, M.con')  # [(lb, lpb), lpc, rpt]
             tmp6 = reshape(tmp5, (slb, slpb * slpc, srpt))  # [lb, (lpb, lpc), rpt]
             tmp7 = reshape(tmp6[:, pl_b_c, :], (slb * slpbc, srpt))  # [(lb, lpbc), rpt]
-            for ilt ∈ 1 : slt
+            for ilt ∈ 1:slt
                 tmp8[:, pl_bc_t] = (@view LE[:, ilt, :])  # [lb, (lpbc, lpt)]
                 LR[ilt, irt, :, :] = reshape(tmp8, (slb * slpbc, slpt))' * tmp7  # [lct, rct]
             end
@@ -221,7 +239,12 @@ function project_ket_on_bra(LE::S, B::S, M::VirtualTensor{R, 4}, RE::S) where {S
 end
 
 
-function update_env_right(RE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <: Tensor{R, 3}} where R <: Real
+function update_env_right(
+    RE::S,
+    A::S,
+    M::VirtualTensor{R,4},
+    B::S,
+) where {S<:Tensor{R,3}} where {R<:Real}
     p_lb, p_lc, p_lt, p_rb, p_rc, p_rt = M.projs
     slb, srb = size(B, 1), size(B, 2)
     slt, srt = size(A, 1), size(A, 2)
@@ -240,8 +263,10 @@ function update_env_right(RE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
         B2 = permutedims(B, (1, 3, 2, 4))  # [lb, lpb, rb, rpb]
         B2 = reshape(B2, (slb * slpb, srb * srpb))  # [(lb, lpb), (rb, rpb)]
 
-        pr_b_ct, pr_c_t, srpct = merge_projectors_inter(M.lp, p_rb, p_rc, p_rt, onGPU; order="1_23")
-        pl_bc_t, pl_b_c, slpbc = merge_projectors_inter(M.lp, p_lt, p_lb, p_lc, onGPU; order="23_1")
+        pr_b_ct, pr_c_t, srpct =
+            merge_projectors_inter(M.lp, p_rb, p_rc, p_rt, onGPU; order = "1_23")
+        pl_bc_t, pl_b_c, slpbc =
+            merge_projectors_inter(M.lp, p_lt, p_lb, p_lc, onGPU; order = "23_1")
 
         tmp1 = alloc_zeros(R, onGPU, (srb, srpb * srpct))
         tmp2 = alloc_undef(R, onGPU, (slb * slpb, srpct))
@@ -249,7 +274,7 @@ function update_env_right(RE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
         tmp5 = alloc_undef(R, onGPU, (slb * slpb, slpc, srpt))
         tmp8 = alloc_undef(R, onGPU, (slb * slpbc, slpt))
 
-        for irt ∈ 1 : srt
+        for irt ∈ 1:srt
             tmp1[:, pr_b_ct] = (@view RE[:, irt, :])  # [rb, (rpb, rpct)]
             mul!(tmp2, B2, reshape(tmp1, (srb * srpb, srpct)))  # [(lb, lpb), rpct]
             tmp3[:, pr_c_t] = tmp2  # [(lb, lpb), (rpc, rpt)]
@@ -257,7 +282,7 @@ function update_env_right(RE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
             batched_mul!(tmp5, tmp4, M.con')
             tmp6 = reshape(tmp5, (slb, slpb * slpc, srpt))  # [lb, (lpb, lpc), rpt]
             tmp7 = reshape(tmp6[:, pl_b_c, :], (slb * slpbc, srpt))  # [(lb, lpbc), rpt]
-            for ilt ∈ 1 : slt
+            for ilt ∈ 1:slt
                 mul!(tmp8, tmp7, (@view A[ilt, irt, :, :])')
                 tmp9 = reshape(tmp8, (slb, slpbc * slpt))
                 Rout[:, ilt, :] .+= tmp9[:, pl_bc_t]
@@ -267,15 +292,17 @@ function update_env_right(RE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
         A2 = permutedims(A, (1, 3, 2, 4))  # [lt, lpt, rt, rpt]
         A2 = reshape(A2, (slt * slpt, srt * srpt))  # [(lt, lpt), (rt, rpt)]
 
-        pr_t_cb, pr_c_b, srpcb = merge_projectors_inter(M.lp, p_rt, p_rc, p_rb, onGPU; order="1_23")
-        pl_tc_b, pl_t_c, slptc = merge_projectors_inter(M.lp, p_lb, p_lt, p_lc, onGPU; order="23_1")
+        pr_t_cb, pr_c_b, srpcb =
+            merge_projectors_inter(M.lp, p_rt, p_rc, p_rb, onGPU; order = "1_23")
+        pl_tc_b, pl_t_c, slptc =
+            merge_projectors_inter(M.lp, p_lb, p_lt, p_lc, onGPU; order = "23_1")
 
         tmp1 = alloc_zeros(R, onGPU, (srt, srpt * srpcb))
         tmp2 = alloc_undef(R, onGPU, (slt * slpt, srpcb))
         tmp3 = alloc_zeros(R, onGPU, (slt * slpt, srpc * srpb))
         tmp5 = alloc_undef(R, onGPU, (slt * slpt, slpc, srpb))
         tmp8 = alloc_undef(R, onGPU, (slt * slptc, slpb))
-        for irb ∈ 1 : srb
+        for irb ∈ 1:srb
             tmp1[:, pr_t_cb] = (@view RE[irb, :, :])  # [rt, (rpt, rpcb)]
             mul!(tmp2, A2, reshape(tmp1, (srt * srpt, srpcb)))  # [(lt, lpt), rpcb]
             tmp3[:, pr_c_b] = tmp2  # [(lt, lpt), (rpc, rpb)]
@@ -283,7 +310,7 @@ function update_env_right(RE::S, A::S, M::VirtualTensor{R, 4}, B::S) where {S <:
             batched_mul!(tmp5, tmp4, M.con')  # [(lt, lpt), lpc, rpb]
             tmp6 = reshape(tmp5, (slt, slpt * slpc, srpb))  # [lt, (lpt, lpc), rpb]
             tmp7 = reshape(tmp6[:, pl_t_c, :], (slt * slptc, srpb))  # [(lb, lptc), rpb]
-            for ilb ∈ 1 : slb
+            for ilb ∈ 1:slb
                 mul!(tmp8, tmp7, (@view B[ilb, irb, :, :])')
                 tmp9 = reshape(tmp8, (slt, slptc * slpb))
                 Rout[ilb, :, :] .+= tmp9[:, pl_tc_b]
@@ -295,8 +322,11 @@ end
 
 
 function update_reduced_env_right(
-    K::Tensor{R, 1}, RE::Tensor{R, 2}, M::VirtualTensor{R, 4}, B::Tensor{R, 3}
-) where R <: Real
+    K::Tensor{R,1},
+    RE::Tensor{R,2},
+    M::VirtualTensor{R,4},
+    B::Tensor{R,3},
+) where {R<:Real}
     p_lb, p_lc, p_lt, p_rb, p_rc, p_rt = M.projs
 
     slb, srb = size(B, 1), size(B, 2)
@@ -311,8 +341,10 @@ function update_reduced_env_right(
     B2 = reshape(B2, (slb * slpb, srb * srpb))  # [(lb, lpb), (rb, rpb)]
 
     if srpb * slpt >= srpt * slpb
-        pr_b_ct, pr_c_t, srpct = merge_projectors_inter(M.lp, p_rb, p_rc, p_rt, onGPU; order="1_23")
-        pl_bc_t, pl_b_c, slpbc = merge_projectors_inter(M.lp, p_lt, p_lb, p_lc, onGPU; order="23_1")
+        pr_b_ct, pr_c_t, srpct =
+            merge_projectors_inter(M.lp, p_rb, p_rc, p_rt, onGPU; order = "1_23")
+        pl_bc_t, pl_b_c, slpbc =
+            merge_projectors_inter(M.lp, p_lt, p_lb, p_lc, onGPU; order = "23_1")
 
         tmp1 = alloc_zeros(R, onGPU, (srb, srpb * srpct))
         tmp4 = alloc_zeros(R, onGPU, (slb * slpb, srpc * srpt))
@@ -330,8 +362,10 @@ function update_reduced_env_right(
         tmp10 = reshape(tmp9, (slb, slpbc * slpt))
         Rtemp = tmp10[:, pl_bc_t]
     else
-        pr_bc_t, pr_b_c, srpbc = merge_projectors_inter(M.lp, p_rt, p_rb, p_rc, onGPU; order="23_1")
-        pl_b_ct, pl_c_t, slpct = merge_projectors_inter(M.lp, p_lb, p_lc, p_lt, onGPU; order="1_23")
+        pr_bc_t, pr_b_c, srpbc =
+            merge_projectors_inter(M.lp, p_rt, p_rb, p_rc, onGPU; order = "23_1")
+        pl_b_ct, pl_c_t, slpct =
+            merge_projectors_inter(M.lp, p_lb, p_lc, p_lt, onGPU; order = "1_23")
 
         tmp1 = alloc_zeros(R, onGPU, (srb, srpbc * srpt))
         tmp4 = alloc_zeros(R, onGPU, (srb, srpb * srpc, slpt))
@@ -353,7 +387,7 @@ function update_reduced_env_right(
 end
 
 
-function contract_tensors43(M::VirtualTensor{R, 4}, B::Tensor{R, 3}) where R <: Real
+function contract_tensors43(M::VirtualTensor{R,4}, B::Tensor{R,3}) where {R<:Real}
     p_lb, p_l, p_lt, p_rb, p_r, p_rt = M.projs
 
     slb, srb = size(B, 1), size(B, 2)
@@ -378,13 +412,18 @@ function contract_tensors43(M::VirtualTensor{R, 4}, B::Tensor{R, 3}) where R <:
     else
         pls = contract_tensor3_matrix(pls, M.con)
     end
-    @tensor order=(lb, c, rb) MB[l, lt, r, rt] := pls[lt, lb, c] * prs[rt, rb, c] * B[l, r, lb, rb]
+    @tensor order = (lb, c, rb) MB[l, lt, r, rt] :=
+        pls[lt, lb, c] * prs[rt, rb, c] * B[l, r, lb, rb]
     MB = reshape(MB, slb, slct, slcp, srb, srct, srcp)
     MB = permutedims(MB, (1, 3, 4, 6, 2, 5))
     reshape(MB, (slb * slcp, srb * srcp, slct * srct))
 end
 
-function corner_matrix(C::S, M::T, B::S) where {S <: Tensor{R, 3}, T <: VirtualTensor{R, 4}} where R <: Real
+function corner_matrix(
+    C::S,
+    M::T,
+    B::S,
+) where {S<:Tensor{R,3},T<:VirtualTensor{R,4}} where {R<:Real}
     slb, srb = size(B, 1), size(B, 2)
     srcc, stc = size(C, 2), size(C, 3)
     V = contract_tensors43(M, B)
@@ -392,4 +431,4 @@ function corner_matrix(C::S, M::T, B::S) where {S <: Tensor{R, 3}, T <: VirtualT
     V = reshape(V, (vl, srb, stc, vt))
     @tensor Cnew[vl, vt, vrr] := V[vl, srb, stc, vt] * C[srb, vrr, stc]
     reshape(Cnew, (slb, :, srcc, vt))
-end
\ No newline at end of file
+end
diff --git a/src/environment.jl b/src/environment.jl
index 99a7189..b18caf5 100644
--- a/src/environment.jl
+++ b/src/environment.jl
@@ -1,26 +1,27 @@
-export
-    Environment,
-    EnvironmentMixed,
-    left_nbrs_site,
-    right_nbrs_site
+export Environment, EnvironmentMixed, left_nbrs_site, right_nbrs_site
 
 abstract type AbstractEnvironment end
 
-mutable struct EnvironmentMixed{T <: Real} <: AbstractEnvironment
+mutable struct EnvironmentMixed{T<:Real} <: AbstractEnvironment
     bra::QMps{T}  # mps that is to be optimized
     mpo::QMpo{T}
     ket::QMps{T}
-    C::Tensor{T, 3}
-    site  # position of C is at: site - epsilon  ::Union(Sites, :central)
+    C::Tensor{T,3}
+    site::Any  # position of C is at: site - epsilon  ::Union(Sites, :central)
     env::Dict
     onGPU::Bool
 
-    function EnvironmentMixed(bra::QMps{T}, C::Tensor{T, 3}, mpo::QMpo{T}, ket::QMps{T};) where T <: Real
+    function EnvironmentMixed(
+        bra::QMps{T},
+        C::Tensor{T,3},
+        mpo::QMpo{T},
+        ket::QMps{T};
+    ) where {T<:Real}
         onGPU = bra.onGPU && mpo.onGPU && ket.onGPU
         @assert bra.sites == ket.sites && issubset(bra.sites, mpo.sites)
         id3 = onGPU ? CUDA.ones(T, 1, 1, 1) : ones(T, 1, 1, 1)
         id2 = onGPU ? CUDA.ones(T, 1, 1) : ones(T, 1, 1)
-        env0 = Dict{Any, Any}((bra.sites[1], :left) => id2, (bra.sites[end], :right) => id3)
+        env0 = Dict{Any,Any}((bra.sites[1], :left) => id2, (bra.sites[end], :right) => id3)
         env = new{T}(bra, mpo, ket, C, last(bra.sites) + 1, env0, onGPU)  #
         update_env_left!.(Ref(env), env.bra.sites)
         env
@@ -46,14 +47,14 @@ function clear_env_containing_site!(env::EnvironmentMixed, site)
     end
 end
 
-mutable struct Environment{T <: Real} <: AbstractEnvironment
+mutable struct Environment{T<:Real} <: AbstractEnvironment
     bra::QMps{T}  # mps that is to be optimized
     mpo::QMpo{T}
     ket::QMps{T}
     env::Dict
     log_norms::Dict
 
-    function Environment(bra::QMps{T}, mpo::QMpo{T}, ket::QMps{T}) where T <: Real
+    function Environment(bra::QMps{T}, mpo::QMpo{T}, ket::QMps{T}) where {T<:Real}
         onGPU = bra.onGPU && mpo.onGPU && ket.onGPU
         @assert bra.sites == ket.sites && issubset(bra.sites, mpo.sites)
         id = onGPU ? CUDA.ones(T, 1, 1, 1) : ones(T, 1, 1, 1)
@@ -99,10 +100,14 @@ function update_env_left(
     LE::S,
     A::S,
     M::T,
-    B::S
-    ) where {S <: AbstractArray{R, 3}, T <: MpoTensor{R, 4}} where R <: Real
-    for v ∈ M.top A = contract_tensor3_matrix(A, v) end
-    for v ∈ reverse(M.bot) B = contract_matrix_tensor3(v, B) end
+    B::S,
+) where {S<:AbstractArray{R,3},T<:MpoTensor{R,4}} where {R<:Real}
+    for v ∈ M.top
+        A = contract_tensor3_matrix(A, v)
+    end
+    for v ∈ reverse(M.bot)
+        B = contract_matrix_tensor3(v, B)
+    end
     update_env_left(LE, A, M.ctr, B)
 end
 
@@ -122,7 +127,7 @@ function update_env_left!(env::Environment, site::Site)
     push!(env.log_norms, (site, :left) => nLL)
 end
 
-function update_env_left!(env::EnvironmentMixed{T}, site) where T    # site::Union(Sites, :central)
+function update_env_left!(env::EnvironmentMixed{T}, site) where {T}    # site::Union(Sites, :central)
     if site == first(env.bra.sites)
         if env.site == first(env.bra.sites)
             LL = env.onGPU ? CUDA.ones(T, 1, 1, 1) : ones(T, 1, 1, 1)
@@ -167,10 +172,17 @@ end
       -- B --
 """
 function update_env_right(
-    RE::S, A::S1, M::T, B::S
-) where {T <: MpoTensor{R, 4}, S <: AbstractArray{R, 3}, S1 <: AbstractArray{R, 3}} where R <: Real
-    for v ∈ M.top  A = contract_tensor3_matrix(A, v) end
-    for v ∈ reverse(M.bot) B = contract_matrix_tensor3(v, B) end
+    RE::S,
+    A::S1,
+    M::T,
+    B::S,
+) where {T<:MpoTensor{R,4},S<:AbstractArray{R,3},S1<:AbstractArray{R,3}} where {R<:Real}
+    for v ∈ M.top
+        A = contract_tensor3_matrix(A, v)
+    end
+    for v ∈ reverse(M.bot)
+        B = contract_matrix_tensor3(v, B)
+    end
     update_env_right(RE, A, M.ctr, B)
 end
 
@@ -190,7 +202,7 @@ function update_env_right!(env::Environment, site::Site)
     push!(env.log_norms, (site, :right) => nRR)
 end
 
-function update_env_right!(env::EnvironmentMixed{T}, site) where T   # site::Union(Sites, :central)
+function update_env_right!(env::EnvironmentMixed{T}, site) where {T}   # site::Union(Sites, :central)
     if site == last(env.bra.sites)
         if env.site > last(env.bra.sites)
             RR = env.onGPU ? CUDA.ones(T, 1, 1) : ones(T, 1, 1)
@@ -231,26 +243,42 @@ function project_ket_on_bra(
     LE::S,
     B::S,
     M::T,
-    RE::S
-    ) where {S <: AbstractArray{R, 3}, T <: MpoTensor{R, 4}} where R <: Real
-    for v ∈ reverse(M.bot) B = contract_matrix_tensor3(v, B) end
+    RE::S,
+) where {S<:AbstractArray{R,3},T<:MpoTensor{R,4}} where {R<:Real}
+    for v ∈ reverse(M.bot)
+        B = contract_matrix_tensor3(v, B)
+    end
     B = project_ket_on_bra(LE, B, M.ctr, RE)
-    for v ∈ reverse(M.top) B = contract_matrix_tensor3(v, B) end
+    for v ∈ reverse(M.top)
+        B = contract_matrix_tensor3(v, B)
+    end
     B
 end
 
 
 project_ket_on_bra(env::Environment, site::Site) = project_ket_on_bra(
-    env.env[(site, :left)], env.ket[site], env.mpo[site], env.env[(site, :right)]
+    env.env[(site, :left)],
+    env.ket[site],
+    env.mpo[site],
+    env.env[(site, :right)],
 )
 
 function project_ket_on_bra(env::EnvironmentMixed, site)
     if site == :central
         B = project_ket_on_bra(env.env[(site, :left)], env.env[(site, :right)])
     elseif site >= env.site
-        B = project_ket_on_bra(env.env[(site, :left)], env.ket[site], env.mpo[site], env.env[(site, :right)])
+        B = project_ket_on_bra(
+            env.env[(site, :left)],
+            env.ket[site],
+            env.mpo[site],
+            env.env[(site, :right)],
+        )
     else
-        B = project_ket_on_bra(env.env[(site, :left)], env.ket[site], env.env[(site, :right)])
+        B = project_ket_on_bra(
+            env.env[(site, :left)],
+            env.ket[site],
+            env.env[(site, :right)],
+        )
     end
     B
 end
diff --git a/src/gauges.jl b/src/gauges.jl
index a847d06..c55cc78 100644
--- a/src/gauges.jl
+++ b/src/gauges.jl
@@ -1,11 +1,9 @@
 
 # gauges.jl: This file provides basic functions to optimize gauges for the PEPS network. CUDA is supported.
 
-export
-    optimize_gauges_for_overlaps!!,
-    overlap_density_matrix
+export optimize_gauges_for_overlaps!!, overlap_density_matrix
 
-function update_rq!(ψ::QMps{T}, AT::Array{T, 3}, i::Site) where T <: Real
+function update_rq!(ψ::QMps{T}, AT::Array{T,3}, i::Site) where {T<:Real}
     @cast ATR[x, (σ, y)] := AT[x, σ, y]
     RT, QT = rq_fact(ATR)
     RT ./= maximum(abs.(RT))
@@ -14,7 +12,7 @@ function update_rq!(ψ::QMps{T}, AT::Array{T, 3}, i::Site) where T <: Real
     RT
 end
 
-function update_rq!(ψ::QMps{T}, AT::CuArray{T, 3}, i::Site) where T <: Real
+function update_rq!(ψ::QMps{T}, AT::CuArray{T,3}, i::Site) where {T<:Real}
     @cast ATR[x, (σ, y)] := AT[x, σ, y]
     RT, QT = rq_fact(ATR)
     RT ./= maximum(abs.(RT))
@@ -23,7 +21,7 @@ function update_rq!(ψ::QMps{T}, AT::CuArray{T, 3}, i::Site) where T <: Real
     RT
 end
 
-function update_qr!(ψ::QMps{T}, AT::Array{T, 3}, i::Site) where T <: Real
+function update_qr!(ψ::QMps{T}, AT::Array{T,3}, i::Site) where {T<:Real}
     @cast ATR[(x, σ), y] := AT[x, σ, y]
     QT, RT = qr_fact(ATR)
     RT ./= maximum(abs.(RT))
@@ -32,7 +30,7 @@ function update_qr!(ψ::QMps{T}, AT::Array{T, 3}, i::Site) where T <: Real
     RT
 end
 
-function update_qr!(ψ::QMps{T}, AT::CuArray{T, 3}, i::Site) where T <: Real
+function update_qr!(ψ::QMps{T}, AT::CuArray{T,3}, i::Site) where {T<:Real}
     @cast ATR[(x, σ), y] := AT[x, σ, y]
     QT, RT = qr_fact(ATR)
     RT ./= maximum(abs.(RT))
@@ -41,7 +39,12 @@ function update_qr!(ψ::QMps{T}, AT::CuArray{T, 3}, i::Site) where T <: Real
     RT
 end
 
-function _gauges_right_sweep!!!(ψ_top::QMps{R}, ψ_bot::QMps{R}, gauges::Dict; tol=1E-12) where R <: Real
+function _gauges_right_sweep!!!(
+    ψ_top::QMps{R},
+    ψ_bot::QMps{R},
+    gauges::Dict;
+    tol = 1E-12,
+) where {R<:Real}
     RT = ψ_top.onGPU && ψ_bot.onGPU ? CUDA.ones(R, 1, 1) : ones(R, 1, 1)
     RB = copy(RT)
     for i ∈ ψ_top.sites
@@ -66,7 +69,12 @@ function _gauges_right_sweep!!!(ψ_top::QMps{R}, ψ_bot::QMps{R}, gauges::Dict;
     end
 end
 
-function _gauges_left_sweep!!!(ψ_top::QMps{R}, ψ_bot::QMps{R}, gauges::Dict; tol=1E-12) where R <: Real
+function _gauges_left_sweep!!!(
+    ψ_top::QMps{R},
+    ψ_bot::QMps{R},
+    gauges::Dict;
+    tol = 1E-12,
+) where {R<:Real}
     RT = ψ_top.onGPU && ψ_bot.onGPU ? CUDA.ones(R, 1, 1) : ones(R, 1, 1)
     RB = copy(RT)
     for i ∈ reverse(ψ_top.sites)
@@ -91,7 +99,12 @@ function _gauges_left_sweep!!!(ψ_top::QMps{R}, ψ_bot::QMps{R}, gauges::Dict; t
     end
 end
 
-function optimize_gauges_for_overlaps!!(ψ_top::QMps{T}, ψ_bot::QMps{T}, tol=1E-8, max_sweeps::Int=4) where T <: Real
+function optimize_gauges_for_overlaps!!(
+    ψ_top::QMps{T},
+    ψ_bot::QMps{T},
+    tol = 1E-8,
+    max_sweeps::Int = 4,
+) where {T<:Real}
     onGPU = ψ_top.onGPU && ψ_bot.onGPU
     canonise!(ψ_top, :right)
     canonise!(ψ_bot, :right)
@@ -104,12 +117,14 @@ function optimize_gauges_for_overlaps!!(ψ_top::QMps{T}, ψ_bot::QMps{T}, tol=1E
         overlap_new = dot(ψ_top, ψ_bot)
         Δ = overlap_new / overlap_old
         overlap_old = overlap_new
-        if abs(Δ - one(T)) < tol break end
+        if abs(Δ - one(T)) < tol
+            break
+        end
     end
     gauges
 end
 
-function overlap_density_matrix(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where T <: Real
+function overlap_density_matrix(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where {T<:Real}
     @assert ψ.sites == ϕ.sites
     C = _overlap_forward(ϕ, ψ, k)
     D = _overlap_backwards(ϕ, ψ, k)
@@ -117,7 +132,7 @@ function overlap_density_matrix(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where T <: Re
     @tensor E[x, y] := C[b, a] * conj(B)[b, x, β] * A[a, y, α] * D[β, α]
 end
 
-function _overlap_forward(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where T <: Real
+function _overlap_forward(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where {T<:Real}
     C = ϕ.onGPU && ψ.onGPU ? CUDA.ones(T, 1, 1) : ones(T, 1, 1)
     i = ψ.sites[1]
     while i < k
@@ -128,7 +143,7 @@ function _overlap_forward(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where T <: Real
     C
 end
 
-function _overlap_backwards(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where T <: Real
+function _overlap_backwards(ϕ::QMps{T}, ψ::QMps{T}, k::Site) where {T<:Real}
     D = ϕ.onGPU && ψ.onGPU ? CUDA.ones(T, 1, 1) : ones(T, 1, 1)
     i = ψ.sites[end]
     while i > k
diff --git a/src/linear_algebra_ext.jl b/src/linear_algebra_ext.jl
index 7b66a0c..baea8ec 100644
--- a/src/linear_algebra_ext.jl
+++ b/src/linear_algebra_ext.jl
@@ -2,30 +2,33 @@
 # linear_algebra_ext.jl: This file provides basic functions to perform custom SVD, and QR.
 #                        Both are calculated on CPU, but can be transferd to GPU if need be.
 
-export
-    rq_fact,
-    qr_fact,
-    svd_fact
+export rq_fact, qr_fact, svd_fact
 
-@inline phase(d::T; atol=eps()) where T <: Real = isapprox(d, zero(T), atol=atol) ? one(T) : d / abs(d)
-@inline phase(d::AbstractArray; atol=eps()) = map(x -> phase(x; atol=atol), d)
+@inline phase(d::T; atol = eps()) where {T<:Real} =
+    isapprox(d, zero(T), atol = atol) ? one(T) : d / abs(d)
+@inline phase(d::AbstractArray; atol = eps()) = map(x -> phase(x; atol = atol), d)
 
-function svd_fact(A::AbstractMatrix{T}, Dcut::Int=typemax(Int), tol=eps(T); kwargs...) where T <: Real
+function svd_fact(
+    A::AbstractMatrix{T},
+    Dcut::Int = typemax(Int),
+    tol = eps(T);
+    kwargs...,
+) where {T<:Real}
     U, Σ, V = svd(A; kwargs...)
     δ = min(Dcut, sum(Σ .> Σ[1] * max(eps(), tol)))
     U, Σ, V = U[:, 1:δ], Σ[1:δ], V[:, 1:δ]
     Σ ./= sqrt(sum(Σ .^ 2))
-    ϕ = reshape(phase(diag(U); atol=tol), 1, :)
+    ϕ = reshape(phase(diag(U); atol = tol), 1, :)
     U .* ϕ, Σ, V .* ϕ
 end
 
 function qr_fact(
-    M::AbstractMatrix{T}, 
-    Dcut::Int=typemax(Int), 
-    tol::T=eps(); 
-    toGPU::Bool=true, 
-    kwargs...
-    ) where T <: Real
+    M::AbstractMatrix{T},
+    Dcut::Int = typemax(Int),
+    tol::T = eps();
+    toGPU::Bool = true,
+    kwargs...,
+) where {T<:Real}
     q, r = qr_fix(qr(Array(M); kwargs...))
     if Dcut >= size(q, 2)
         toGPU && return CuArray.((q, r))
@@ -37,18 +40,18 @@ function qr_fact(
 end
 
 function rq_fact(
-    M::AbstractMatrix{T}, 
-    Dcut::Int=typemax(Int), 
-    tol::T=eps(); 
-    toGPU::Bool=true, 
-    kwargs...
-    ) where T <: Real
-    q, r = qr_fact(M', Dcut, tol; toGPU=toGPU, kwargs...)
+    M::AbstractMatrix{T},
+    Dcut::Int = typemax(Int),
+    tol::T = eps();
+    toGPU::Bool = true,
+    kwargs...,
+) where {T<:Real}
+    q, r = qr_fact(M', Dcut, tol; toGPU = toGPU, kwargs...)
     toGPU && return CuArray.((r', q'))
     r', q'
 end
 
-function qr_fix(QR_fact; tol::T=eps()) where T <: Real
-    ϕ = phase(diag(QR_fact.R); atol=tol)
+function qr_fix(QR_fact; tol::T = eps()) where {T<:Real}
+    ϕ = phase(diag(QR_fact.R); atol = tol)
     QR_fact.Q * Diagonal(ϕ), ϕ .* QR_fact.R
 end
diff --git a/src/mps/base.jl b/src/mps/base.jl
index afb6f85..fe72822 100644
--- a/src/mps/base.jl
+++ b/src/mps/base.jl
@@ -1,17 +1,12 @@
 # ./mps/base.jl: This file provides basic definitions of custom Matrix Product States / Operators.
 
-export
-    Site,
-    Sites,
-    AbstractTensorNetwork,
-    MpoTensor,
-    QMpsOrMpo
+export Site, Sites, AbstractTensorNetwork, MpoTensor, QMpsOrMpo
 
 abstract type AbstractTensorNetwork{T} end
 
-const Site = Union{Int, Rational{Int}}
-const Sites = NTuple{N, Site} where N
-const TensorMap{T} = Dict{Site, Union{Tensor{T, 2}, Tensor{T, 3}, Tensor{T, 4}}}  # 2 and 4 - mpo;  3 - mps
+const Site = Union{Int,Rational{Int}}
+const Sites = NTuple{N,Site} where {N}
+const TensorMap{T} = Dict{Site,Union{Tensor{T,2},Tensor{T,3},Tensor{T,4}}}  # 2 and 4 - mpo;  3 - mps
 
 """
 A mutable struct representing a Matrix Product Operator (MPO) tensor in a tensor network.
@@ -29,10 +24,10 @@ The top and bottom legs are vectors of two-dimensional tensors (`Tensor{T, 2}`).
 The central tensor is of type `Tensor{T, N}` or `Nothing` if not present. 
 The dimensions of the MPO tensor are specified by `dims`. 
 """
-mutable struct MpoTensor{T <: Real, N}
-    top::Vector{Tensor{T, 2}}  # N == 2 top = []
-    ctr::Union{Tensor{T, N}, Nothing}
-    bot::Vector{Tensor{T, 2}}  # N == 2 bot = []
+mutable struct MpoTensor{T<:Real,N}
+    top::Vector{Tensor{T,2}}  # N == 2 top = []
+    ctr::Union{Tensor{T,N},Nothing}
+    bot::Vector{Tensor{T,2}}  # N == 2 bot = []
     dims::Dims{N}
 end
 
@@ -41,10 +36,10 @@ end
 
 #  ## Arguments
 # - `ten::TensorMap{T}`: A dictionary mapping `Site` indices to tensors of type `Tensor{T, 2}`, `Tensor{T, 3}`, or `Tensor{T, 4}`. The key `0` represents the central tensor (`ctr`) of the MPO.
-    
+
 # ## Returns
 # - An instance of `MpoTensor{T, nn}` representing the Matrix Product Operator.
-    
+
 # ## Description
 # The `MpoTensor` function constructs a Matrix Product Operator tensor from a `TensorMap`, 
 # where the keys are `Site` indices and the values are tensors of appropriate dimensions. 
@@ -52,11 +47,11 @@ end
 # and separating tensors into the top, central, and bottom parts. 
 # The central tensor is identified by the key `0`. 
 # The resulting `MpoTensor` encapsulates the tensors along with their dimensions.
-        
+
 # ## Exceptions
 # - Throws a `DomainError` if the central tensor (`ctr`) has dimensions other than 2 or 4.
 # """ 
-function MpoTensor(ten::TensorMap{T}) where T
+function MpoTensor(ten::TensorMap{T}) where {T}
     sk = sort(collect(keys(ten)))
     top = [ten[k] for k ∈ sk if k < 0]
     bot = [ten[k] for k ∈ sk if k > 0]
@@ -73,41 +68,43 @@ function MpoTensor(ten::TensorMap{T}) where T
             dims = size(ctr)
         elseif nn == 4
             dims = (
-                size(ctr, 1), isempty(top) ? size(ctr, 2) : size(top[1], 1),
-                size(ctr, 3), isempty(bot) ? size(ctr, 4) : size(bot[end], 2)
+                size(ctr, 1),
+                isempty(top) ? size(ctr, 2) : size(top[1], 1),
+                size(ctr, 3),
+                isempty(bot) ? size(ctr, 4) : size(bot[end], 2),
             )
         else
             throw(DomainError(ndims(ctr), "MpoTensor should have ndims 2 or 4"))
         end
     end
-    MpoTensor{T, nn}(top, ctr, bot, dims)
+    MpoTensor{T,nn}(top, ctr, bot, dims)
 end
 
-Base.eltype(ten::MpoTensor{T, N}) where {T, N} = T
-Base.ndims(ten::MpoTensor{T, N}) where {T, N} = N
+Base.eltype(ten::MpoTensor{T,N}) where {T,N} = T
+Base.ndims(ten::MpoTensor{T,N}) where {T,N} = N
 Base.size(ten::MpoTensor, n::Int) = ten.dims[n]
 Base.size(ten::MpoTensor) = ten.dims
 
-const MpoTensorMap{T} = Dict{Site, MpoTensor{T}}
+const MpoTensorMap{T} = Dict{Site,MpoTensor{T}}
 
 for (S, M) ∈ ((:QMpo, :MpoTensorMap), (:QMps, :TensorMap))
     @eval begin
         export $S, $M
-        mutable struct $S{F <: Real} <: AbstractTensorNetwork{F}
+        mutable struct $S{F<:Real} <: AbstractTensorNetwork{F}
             tensors::$M{F}
             sites::Vector{Site}
             onGPU::Bool
 
-            function $S(ten::$M{F}; onGPU::Bool=false) where F
+            function $S(ten::$M{F}; onGPU::Bool = false) where {F}
                 new{F}(ten, sort(collect(keys(ten))), onGPU)
             end
         end
     end
 end
 
-const QMpsOrMpo{T} = Union{QMpo{T}, QMps{T}}
+const QMpsOrMpo{T} = Union{QMpo{T},QMps{T}}
 
 @inline Base.getindex(ψ::QMpsOrMpo, i) = getindex(ψ.tensors, i)
 @inline Base.setindex!(ψ::QMpsOrMpo, A, i::Site) = ψ.tensors[i] = A
-@inline Base.eltype(ψ::QMpsOrMpo{T}) where T = T
-@inline Base.copy(ψ::QMps) = QMps(copy(ψ.tensors), onGPU=ψ.onGPU)
\ No newline at end of file
+@inline Base.eltype(ψ::QMpsOrMpo{T}) where {T} = T
+@inline Base.copy(ψ::QMps) = QMps(copy(ψ.tensors), onGPU = ψ.onGPU)
diff --git a/src/mps/canonise.jl b/src/mps/canonise.jl
index e393e45..ab0a381 100644
--- a/src/mps/canonise.jl
+++ b/src/mps/canonise.jl
@@ -1,14 +1,10 @@
 
 # canonise.jl: This file provides basic function to left / right truncate / canonise MPS. CUDA is supported.
 
-export
-    canonise!,
-    truncate!,
-    canonise_truncate!,
-    measure_spectrum
+export canonise!, truncate!, canonise_truncate!, measure_spectrum
 
 
-function measure_spectrum(ψ::QMps{T}) where T <: Real
+function measure_spectrum(ψ::QMps{T}) where {T<:Real}
     # Assume that ψ is left_canonical
     @assert is_left_normalized(ψ)
     R = ones(T, 1, 1)
@@ -17,7 +13,7 @@ function measure_spectrum(ψ::QMps{T}) where T <: Real
         B = permutedims(Array(ψ[i]), (1, 3, 2)) # [x, σ, α]
         @matmul M[x, σ, y] := sum(α) B[x, σ, α] * R[α, y]
         @cast M[x, (σ, y)] := M[x, σ, y]
-        Dcut, tolS = 100000, 0.
+        Dcut, tolS = 100000, 0.0
         U, S, _ = svd_fact(Array(M), Dcut, tolS)
         push!(schmidt, i => S)
         R = U * Diagonal(S)
@@ -27,7 +23,13 @@ end
 
 
 
-function truncate!(ψ::QMps{T}, s::Symbol, Dcut::Int=typemax(Int), tolS::T=eps(); kwargs...) where T <: Real
+function truncate!(
+    ψ::QMps{T},
+    s::Symbol,
+    Dcut::Int = typemax(Int),
+    tolS::T = eps();
+    kwargs...,
+) where {T<:Real}
     @assert s ∈ (:left, :right)
     if s == :right
         _right_sweep!(ψ; kwargs...)
@@ -42,7 +44,13 @@ canonise!(ψ::QMps, s::Symbol) = canonise!(ψ, Val(s))
 canonise!(ψ::QMps, ::Val{:right}) = _left_sweep!(ψ, typemax(Int))
 canonise!(ψ::QMps, ::Val{:left}) = _right_sweep!(ψ, typemax(Int))
 
-function canonise_truncate!(ψ::QMps, type::Symbol, Dcut::Int=typemax(Int), tolS=eps(); kwargs...)
+function canonise_truncate!(
+    ψ::QMps,
+    type::Symbol,
+    Dcut::Int = typemax(Int),
+    tolS = eps();
+    kwargs...,
+)
     if type == :right
         _left_sweep!(ψ, Dcut, tolS; kwargs...)
     elseif type == :left
@@ -52,7 +60,12 @@ function canonise_truncate!(ψ::QMps, type::Symbol, Dcut::Int=typemax(Int), tolS
     end
 end
 
-function _right_sweep!(ψ::QMps{T}, Dcut::Int=typemax(Int), tolS::T=eps(T); kwargs...) where T <: Real
+function _right_sweep!(
+    ψ::QMps{T},
+    Dcut::Int = typemax(Int),
+    tolS::T = eps(T);
+    kwargs...,
+) where {T<:Real}
     R = ψ.onGPU ? CUDA.ones(T, 1, 1) : ones(T, 1, 1)
     for i ∈ ψ.sites
         A = ψ[i]
@@ -66,7 +79,12 @@ function _right_sweep!(ψ::QMps{T}, Dcut::Int=typemax(Int), tolS::T=eps(T); kwar
     end
 end
 
-function _left_sweep!(ψ::QMps{T}, Dcut::Int=typemax(Int), tolS::T=eps(T); kwargs...) where T <: Real
+function _left_sweep!(
+    ψ::QMps{T},
+    Dcut::Int = typemax(Int),
+    tolS::T = eps(T);
+    kwargs...,
+) where {T<:Real}
     R = ψ.onGPU ? CUDA.ones(T, 1, 1) : ones(T, 1, 1)
     for i ∈ reverse(ψ.sites)
         B = permutedims(ψ[i], (1, 3, 2)) # [x, σ, α]
diff --git a/src/mps/dot.jl b/src/mps/dot.jl
index e07ef52..8cefc22 100644
--- a/src/mps/dot.jl
+++ b/src/mps/dot.jl
@@ -6,7 +6,7 @@ LinearAlgebra.norm(ψ::QMps) = sqrt(abs(dot(ψ, ψ)))
 Base.:(*)(ϕ::QMps, ψ::QMps) = dot(ϕ, ψ)
 Base.:(*)(W::QMpo, ψ::QMps) = dot(W, ψ)
 
-function LinearAlgebra.dot(ψ::QMps{T}, ϕ::QMps{T}) where T <: Real
+function LinearAlgebra.dot(ψ::QMps{T}, ϕ::QMps{T}) where {T<:Real}
     @assert ψ.sites == ϕ.sites
     C = ψ.onGPU && ϕ.onGPU ? CUDA.ones(T, 1, 1) : ones(T, 1, 1)
     for i ∈ ϕ.sites
@@ -16,13 +16,17 @@ function LinearAlgebra.dot(ψ::QMps{T}, ϕ::QMps{T}) where T <: Real
     tr(C)
 end
 
-function LinearAlgebra.dot(ψ::QMpo{R}, ϕ::QMps{R}) where R <: Real  
+function LinearAlgebra.dot(ψ::QMpo{R}, ϕ::QMps{R}) where {R<:Real}
     D = TensorMap{R}()
     for i ∈ reverse(ϕ.sites)
         M, B = ψ[i], ϕ[i]
-        for v ∈ reverse(M.bot) B = contract_matrix_tensor3(v, B) end
+        for v ∈ reverse(M.bot)
+            B = contract_matrix_tensor3(v, B)
+        end
         B = contract_tensors43(M.ctr, B)
-        for v ∈ reverse(M.top) B = contract_matrix_tensor3(v, B) end
+        for v ∈ reverse(M.top)
+            B = contract_matrix_tensor3(v, B)
+        end
 
         mps_li = left_nbrs_site(i, ϕ.sites)
         mpo_li = left_nbrs_site(i, ψ.sites)
@@ -42,6 +46,8 @@ function LinearAlgebra.dot(ψ::QMpo{R}, ϕ::QMps{R}) where R <: Real
     QMps(D; onGPU = ψ.onGPU && ϕ.onGPU)
 end
 
-contract_tensor3_matrix(B::AbstractArray{T, 3}, M::MpoTensor{T, 2}) where T <: Real = contract_tensor3_matrix(B, M.ctr)
-contract_matrix_tensor3(M::MpoTensor{T, 2}, B::AbstractArray{T, 3}) where T <: Real = contract_matrix_tensor3(M.ctr, B)
-contract_tensors43(B::Nothing, A::AbstractArray{T, 3}) where T <: Real = A
+contract_tensor3_matrix(B::AbstractArray{T,3}, M::MpoTensor{T,2}) where {T<:Real} =
+    contract_tensor3_matrix(B, M.ctr)
+contract_matrix_tensor3(M::MpoTensor{T,2}, B::AbstractArray{T,3}) where {T<:Real} =
+    contract_matrix_tensor3(M.ctr, B)
+contract_tensors43(B::Nothing, A::AbstractArray{T,3}) where {T<:Real} = A
diff --git a/src/mps/identity.jl b/src/mps/identity.jl
index 37afde2..c07c7bc 100644
--- a/src/mps/identity.jl
+++ b/src/mps/identity.jl
@@ -1,11 +1,14 @@
 # ./mps/identity.jl: This file provides custom MPS Identity. Note, this approach is easier than
 #                    trying to overload the universal identity operator, I, from LinearAlgebra.
 
-export
-    local_dims,
-    IdentityQMps
+export local_dims, IdentityQMps
 
-function IdentityQMps(::Type{T}, loc_dims::Dict, Dmax::Int=1; onGPU=true) where T <: Real
+function IdentityQMps(
+    ::Type{T},
+    loc_dims::Dict,
+    Dmax::Int = 1;
+    onGPU = true,
+) where {T<:Real}
     _zeros = onGPU ? CUDA.zeros : zeros
     id = TensorMap{T}(keys(loc_dims) .=> _zeros.(T, Dmax, Dmax, values(loc_dims)))
 
@@ -21,11 +24,11 @@ function IdentityQMps(::Type{T}, loc_dims::Dict, Dmax::Int=1; onGPU=true) where
     for (site, ld) ∈ loc_dims
         id[site][1, 1, :] .= 1 / sqrt(ld)
     end
-    QMps(id; onGPU=onGPU)
+    QMps(id; onGPU = onGPU)
 end
 
 function local_dims(mpo::QMpo, dir::Symbol)
     @assert dir ∈ (:down, :up)
     dim = dir == :down ? 4 : 2
-    Dict{Site, Int}(k => size(mpo[k], dim) for k ∈ mpo.sites if ndims(mpo[k]) == 4)
+    Dict{Site,Int}(k => size(mpo[k], dim) for k ∈ mpo.sites if ndims(mpo[k]) == 4)
 end
diff --git a/src/mps/rand.jl b/src/mps/rand.jl
index 42f76ee..b96ab9d 100644
--- a/src/mps/rand.jl
+++ b/src/mps/rand.jl
@@ -1,7 +1,12 @@
 
 # ./mps/rand.jl: This file provides methods to generate random MPS / MPO
 
-function Base.rand(::Type{QMps{T}}, loc_dims::Dict, Dmax::Int=1; onGPU=false) where T <: Real
+function Base.rand(
+    ::Type{QMps{T}},
+    loc_dims::Dict,
+    Dmax::Int = 1;
+    onGPU = false,
+) where {T<:Real}
     id = TensorMap{T}(keys(loc_dims) .=> rand.(T, Dmax, Dmax, values(loc_dims)))
     site_min, ld_min = minimum(loc_dims)
     site_max, ld_max = maximum(loc_dims)
@@ -14,23 +19,27 @@ function Base.rand(::Type{QMps{T}}, loc_dims::Dict, Dmax::Int=1; onGPU=false) wh
     onGPU ? move_to_CUDA!(QMps(id)) : QMps(id)
 end
 
-function Base.rand(
-    ::Type{CentralTensor{T}}, s::Vector{Int}
-) where T <: Real
-    CentralTensor(Real.(rand(s[1], s[5])), Real.(rand(s[2], s[6])), Real.(rand(s[3], s[7])), Real.(rand(s[4], s[8])))
+function Base.rand(::Type{CentralTensor{T}}, s::Vector{Int}) where {T<:Real}
+    CentralTensor(
+        Real.(rand(s[1], s[5])),
+        Real.(rand(s[2], s[6])),
+        Real.(rand(s[3], s[7])),
+        Real.(rand(s[4], s[8])),
+    )
 end
 
 function Base.rand(
-    ::Type{SiteTensor{T}}, lp::PoolOfProjectors, l::Int, D::NTuple
-) where T <: Real
+    ::Type{SiteTensor{T}},
+    lp::PoolOfProjectors,
+    l::Int,
+    D::NTuple,
+) where {T<:Real}
     loc_exp = rand(l)
     projs = D
 
     SiteTensor(lp, loc_exp, projs)
 end
 
-function Base.rand(
-    ::Type{QMpo{T}}, loc_dims::Dict; onGPU::Bool=false
-) where T <:Real
+function Base.rand(::Type{QMpo{T}}, loc_dims::Dict; onGPU::Bool = false) where {T<:Real}
     QMpo(MpoTensorMap{T}(loc_dims))
 end
diff --git a/src/mps/transpose.jl b/src/mps/transpose.jl
index 54f8b2d..dbb1ab1 100644
--- a/src/mps/transpose.jl
+++ b/src/mps/transpose.jl
@@ -2,17 +2,20 @@
 # ./mps/transpose.jl: This file defines what it means to transpse MPO. Note, this should not be
 #                     done by overloading Base.transpose for QMpo to avoid overloading (Array)'.
 
-function Base.transpose(ψ::QMpo{T}) where T <: Real
-    QMpo(MpoTensorMap{T}(keys(ψ.tensors) .=> mpo_transpose.(values(ψ.tensors))); onGPU=ψ.onGPU)
+function Base.transpose(ψ::QMpo{T}) where {T<:Real}
+    QMpo(
+        MpoTensorMap{T}(keys(ψ.tensors) .=> mpo_transpose.(values(ψ.tensors)));
+        onGPU = ψ.onGPU,
+    )
 end
 
-mpo_transpose(M::MpoTensor{T, 2}) where T <: Real = M
+mpo_transpose(M::MpoTensor{T,2}) where {T<:Real} = M
 
-function mpo_transpose(M::MpoTensor{T, 4}) where T <: Real
-    MpoTensor{T, 4}(
+function mpo_transpose(M::MpoTensor{T,4}) where {T<:Real}
+    MpoTensor{T,4}(
         mpo_transpose.(reverse(M.bot)),
         mpo_transpose(M.ctr),
         mpo_transpose.(reverse(M.top)),
-        M.dims[[1, 4, 3, 2]]
+        M.dims[[1, 4, 3, 2]],
     )
 end
diff --git a/src/mps/utils.jl b/src/mps/utils.jl
index e40e2c4..5ab360d 100644
--- a/src/mps/utils.jl
+++ b/src/mps/utils.jl
@@ -1,13 +1,7 @@
 # ./mps/aux.jl: This file provides auxiliary functions to verify various MPS properties.
 
-export
-    bond_dimension,
-    bond_dimensions,
-    is_consistent,
-    is_left_normalized,
-    is_right_normalized,
-    length,
-    size
+export bond_dimension,
+    bond_dimensions, is_consistent, is_left_normalized, is_right_normalized, length, size
 
 @inline bond_dimension(ψ::QMpsOrMpo) = maximum(size.(values(ψ.tensors), 1))
 @inline bond_dimensions(ψ::QMpsOrMpo) = [size(ψ.tensors[n]) for n ∈ ψ.sites]
@@ -23,12 +17,16 @@ function is_consistent(ψ::QMps)
         @assert size(ψ.tensors[s1], 2) == size(ψ.tensors[s2], 1) "Incorrect link between $i and $(i+1)."
     end
     dev = which_device(ψ)
-    if  ψ.onGPU @assert :GPU ∈ dev && :CPU ∉ dev end
-    if !ψ.onGPU @assert :GPU ∉ dev && :CPU ∈ dev end
+    if ψ.onGPU
+        @assert :GPU ∈ dev && :CPU ∉ dev
+    end
+    if !ψ.onGPU
+        @assert :GPU ∉ dev && :CPU ∈ dev
+    end
     true
 end
 
-function eye(::Type{T}, dim; toGPU::Bool=false) where T
+function eye(::Type{T}, dim; toGPU::Bool = false) where {T}
     v = ones(T, dim)
     toGPU && return cu(spdiagm(v))
     Diagonal(v)
@@ -36,28 +34,34 @@ end
 
 function is_left_normalized(ψ::QMps, ::Val{false})
     all(
-        eye(eltype(ψ), size(A, 2); toGPU = false) ≈ @tensor(Id[x, y] := A[α, x, σ] * A[α, y, σ];  order = (α, σ)) for A ∈ values(ψ.tensors) # TODO: split the line
-        )
+        eye(eltype(ψ), size(A, 2); toGPU = false) ≈
+        @tensor(Id[x, y] := A[α, x, σ] * A[α, y, σ]; order = (α, σ)) for
+        A ∈ values(ψ.tensors) # TODO: split the line
+    )
 end
 
 function is_left_normalized(ψ::QMps, ::Val{true})
     all(
-        eye(eltype(ψ), size(A, 2); toGPU = true) ≈ @cutensor(Id[x, y] := A[α, x, σ] * A[α, y, σ]) for A ∈ values(ψ.tensors) # TODO: split the line
-        )
+        eye(eltype(ψ), size(A, 2); toGPU = true) ≈
+        @cutensor(Id[x, y] := A[α, x, σ] * A[α, y, σ]) for A ∈ values(ψ.tensors) # TODO: split the line
+    )
 end
 
 is_left_normalized(ψ::QMps) = is_left_normalized(ψ, Val(ψ.onGPU))
 
 function is_right_normalized(ψ::QMps, ::Val{false})
     all(
-        eye(eltype(ψ), size(B, 1); toGPU = false) ≈ @tensor(Id[x, y] := B[x, α, σ] * B[y, α, σ];  order = (α, σ)) for B ∈ values(ψ.tensors) # TODO: split the line
-        )
+        eye(eltype(ψ), size(B, 1); toGPU = false) ≈
+        @tensor(Id[x, y] := B[x, α, σ] * B[y, α, σ]; order = (α, σ)) for
+        B ∈ values(ψ.tensors) # TODO: split the line
+    )
 end
 
 function is_right_normalized(ψ::QMps, ::Val{true})
     all(
-        eye(eltype(ψ), size(B, 1); toGPU = true) ≈ @cutensor(Id[x, y] := B[x, α, σ] * B[y, α, σ]) for B ∈ values(ψ.tensors) # TODO: split the line
-        )
+        eye(eltype(ψ), size(B, 1); toGPU = true) ≈
+        @cutensor(Id[x, y] := B[x, α, σ] * B[y, α, σ]) for B ∈ values(ψ.tensors) # TODO: split the line
+    )
 end
 
 is_right_normalized(ψ::QMps) = is_right_normalized(ψ, Val(ψ.onGPU))
diff --git a/src/projectors.jl b/src/projectors.jl
index 5b6f0cf..fa21437 100644
--- a/src/projectors.jl
+++ b/src/projectors.jl
@@ -1,10 +1,6 @@
-export
-    PoolOfProjectors,
-    get_projector!,
-    add_projector!,
-    empty!
+export PoolOfProjectors, get_projector!, add_projector!, empty!
 
-const Proj{T} = Union{Vector{T}, CuArray{T, 1}}
+const Proj{T} = Union{Vector{T},CuArray{T,1}}
 
 """
 $(TYPEDSIGNATURES)
@@ -24,23 +20,19 @@ The data is provided as a dictionary that maps site indices to projectors stored
 The `sizes` dictionary is automatically populated based on the maximum projector size for each site.
 - `PoolOfProjectors{T}() where T`: Create an empty `PoolOfProjectors` with no projectors initially stored.    
 """
-struct PoolOfProjectors{T <: Integer}
-    data::Dict{Symbol, Dict{Int, Proj{T}}}
+struct PoolOfProjectors{T<:Integer}
+    data::Dict{Symbol,Dict{Int,Proj{T}}}
     default_device::Symbol
-    sizes::Dict{Int, Int}
-
-    PoolOfProjectors(
-        data::Dict{Int, Dict{Int, Vector{T}}}
-        ) where T = new{T}(Dict(:CPU => data), 
-                            :CPU,
-                            Dict{Int, Int}(k => maximum(v) for (k, v) ∈ data))
-    PoolOfProjectors{T}() where T = new{T}(Dict(:CPU => Dict{Int, Proj{T}}()),
-                                           :CPU,
-                                           Dict{Int, Int}())
+    sizes::Dict{Int,Int}
+
+    PoolOfProjectors(data::Dict{Int,Dict{Int,Vector{T}}}) where {T} =
+        new{T}(Dict(:CPU => data), :CPU, Dict{Int,Int}(k => maximum(v) for (k, v) ∈ data))
+    PoolOfProjectors{T}() where {T} =
+        new{T}(Dict(:CPU => Dict{Int,Proj{T}}()), :CPU, Dict{Int,Int}())
 end
 
 
-Base.eltype(lp::PoolOfProjectors{T}) where T = T
+Base.eltype(lp::PoolOfProjectors{T}) where {T} = T
 Base.length(lp::PoolOfProjectors) = length(lp.data[lp.default_device])
 Base.length(lp::PoolOfProjectors, device::Symbol) = length(lp.data[device])
 
@@ -55,7 +47,7 @@ This function removes all projectors stored on the specified computational devic
 - `lp::PoolOfProjectors`: The `PoolOfProjectors` object containing projectors.
 - `device::Symbol`: The computational device for which projectors should be emptied (e.g., `:CPU`, `:GPU`).
 """
-function Base.empty!(lp::PoolOfProjectors, device::Symbol) 
+function Base.empty!(lp::PoolOfProjectors, device::Symbol)
     if device ∈ keys(lp.data)
         empty!(lp.data[device])
     end
@@ -64,7 +56,8 @@ end
 Base.length(lp::PoolOfProjectors, index::Int) = length(lp.data[lp.default_device][index])
 Base.size(lp::PoolOfProjectors, index::Int) = lp.sizes[index]
 
-get_projector!(lp::PoolOfProjectors, index::Int) = get_projector!(lp, index, lp.default_device)
+get_projector!(lp::PoolOfProjectors, index::Int) =
+    get_projector!(lp, index, lp.default_device)
 
 """
 $(TYPEDSIGNATURES)
@@ -84,9 +77,13 @@ If the projector does not exist in the pool, it creates a new one and stores it
 # Returns:
 - `Proj{T}`: The projector of type `T` associated with the specified index and device.
 """
-function get_projector!(lp::PoolOfProjectors{T}, index::Int, device::Symbol) where T <: Integer
+function get_projector!(
+    lp::PoolOfProjectors{T},
+    index::Int,
+    device::Symbol,
+) where {T<:Integer}
     if device ∉ keys(lp.data)
-        push!(lp.data, device => Dict{Int, Proj{T}}())
+        push!(lp.data, device => Dict{Int,Proj{T}}())
     end
 
     if index ∉ keys(lp.data[device])
@@ -118,7 +115,7 @@ The index can be used to retrieve the projector later using `get_projector!`.
 # Returns:
 - `Int`: The unique index associated with the added projector in the pool.
 """
-function add_projector!(lp::PoolOfProjectors{T}, p::Proj) where T <: Integer
+function add_projector!(lp::PoolOfProjectors{T}, p::Proj) where {T<:Integer}
     if lp.default_device == :CPU
         p = Array{T}(p)
     elseif lp.default_device == :GPU
@@ -140,4 +137,4 @@ function add_projector!(lp::PoolOfProjectors{T}, p::Proj) where T <: Integer
         push!(lp.sizes, key => maximum(p))
     end
     key
-end
\ No newline at end of file
+end
diff --git a/src/transfer.jl b/src/transfer.jl
index c759703..7fcac54 100644
--- a/src/transfer.jl
+++ b/src/transfer.jl
@@ -2,15 +2,12 @@
 # transfer.jl: This file provides rules of how to transfer tensors to GPU. Note, NOT all of
 #              tensor's coponents are moved from CPU to GPU and most tensors are generated
 #              on CPU due to the size of clustered Hamiltonian.
-export
-    which_device,
-    move_to_CUDA!,
-    move_to_CPU!
+export which_device, move_to_CUDA!, move_to_CPU!
 
-move_to_CUDA!(ten::Array{T, N}) where {T, N} = CuArray(ten) #cu(ten, unified=true)
+move_to_CUDA!(ten::Array{T,N}) where {T,N} = CuArray(ten) #cu(ten, unified=true)
 
 
-move_to_CUDA!(ten::Union{CuArray{T, N}, Nothing}) where {T, N} = ten
+move_to_CUDA!(ten::Union{CuArray{T,N},Nothing}) where {T,N} = ten
 move_to_CUDA!(ten::Diagonal) = Diagonal(move_to_CUDA!(diag(ten)))
 
 function move_to_CUDA!(ten::CentralTensor)
@@ -40,24 +37,32 @@ function move_to_CUDA!(ten::SiteTensor)
 end
 
 function move_to_CUDA!(ten::MpoTensor)
-    for i ∈ 1:length(ten.top) ten.top[i] = move_to_CUDA!(ten.top[i]) end
-    for i ∈ 1:length(ten.bot) ten.bot[i] = move_to_CUDA!(ten.bot[i]) end
+    for i ∈ 1:length(ten.top)
+        ten.top[i] = move_to_CUDA!(ten.top[i])
+    end
+    for i ∈ 1:length(ten.bot)
+        ten.bot[i] = move_to_CUDA!(ten.bot[i])
+    end
     ten.ctr = move_to_CUDA!(ten.ctr)
     ten
 end
 
-function move_to_CUDA!(ψ::Union{QMpo{T}, QMps{T}}) where T
-    for k ∈ keys(ψ.tensors) ψ[k] = move_to_CUDA!(ψ[k]) end
+function move_to_CUDA!(ψ::Union{QMpo{T},QMps{T}}) where {T}
+    for k ∈ keys(ψ.tensors)
+        ψ[k] = move_to_CUDA!(ψ[k])
+    end
     ψ.onGPU = true
     ψ
 end
 
-move_to_CPU!(ten::CuArray{T, N}) where {T, N} = Array(ten)
-move_to_CPU!(ten::Union{Array{T, N}, Nothing}) where {T, N} = ten
+move_to_CPU!(ten::CuArray{T,N}) where {T,N} = Array(ten)
+move_to_CPU!(ten::Union{Array{T,N},Nothing}) where {T,N} = ten
 move_to_CPU!(ten::Diagonal) = Diagonal(move_to_CPU!(diag(ten)))
 
-function move_to_CPU!(ψ::QMps{T}) where T
-    for k ∈ keys(ψ.tensors) ψ[k] = move_to_CPU!(ψ[k]) end
+function move_to_CPU!(ψ::QMps{T}) where {T}
+    for k ∈ keys(ψ.tensors)
+        ψ[k] = move_to_CPU!(ψ[k])
+    end
     ψ.onGPU = false
     ψ
 end
@@ -65,12 +70,15 @@ end
 
 
 which_device(ten::Nothing) = Set()
-which_device(ψ::Union{QMpo{T}, QMps{T}}) where T = union(which_device.(values(ψ.tensors))...)
-which_device(ten::MpoTensor) = union(which_device(ten.ctr), which_device.(ten.top)..., which_device.(ten.bot)...)
+which_device(ψ::Union{QMpo{T},QMps{T}}) where {T} =
+    union(which_device.(values(ψ.tensors))...)
+which_device(ten::MpoTensor) =
+    union(which_device(ten.ctr), which_device.(ten.top)..., which_device.(ten.bot)...)
 which_device(ten::DiagonalTensor) = union(which_device.((ten.e1, ten.e2))...)
 which_device(ten::VirtualTensor) = union(which_device.((ten.con,))...) # TODO cf. 1)  ten.projs
-which_device(ten::CentralTensor) = union(which_device.((ten.e11, ten.e12, ten.e21, ten.e22))...)
+which_device(ten::CentralTensor) =
+    union(which_device.((ten.e11, ten.e12, ten.e21, ten.e22))...)
 which_device(ten::SiteTensor) = union(which_device.((ten.loc_exp,))...) # TODO cf. 2)  ten.projs
-which_device(ten::Array{T, N}) where {T, N} = Set((:CPU, ))
-which_device(ten::CuArray{T, N}) where {T, N} = Set((:GPU, ))
+which_device(ten::Array{T,N}) where {T,N} = Set((:CPU,))
+which_device(ten::CuArray{T,N}) where {T,N} = Set((:GPU,))
 which_device(ten::Diagonal) = which_device(diag(ten))
diff --git a/src/utils/memory.jl b/src/utils/memory.jl
index cd3e12b..06f33c5 100644
--- a/src/utils/memory.jl
+++ b/src/utils/memory.jl
@@ -1,19 +1,21 @@
-export
-    measure_memory,
-    format_bytes
+export measure_memory, format_bytes
 
 measure_memory(ten::AbstractArray) = [Base.summarysize(ten), 0]  # [CPU_memory, GPU_memory]
 measure_memory(ten::CuArray) = [0, prod(size(ten)) * sizeof(eltype(ten))]
-measure_memory(ten::SparseMatrixCSC) = sum(measure_memory.([ten.colptr, ten.rowval, ten.nzval]))
-measure_memory(ten::CuSparseMatrixCSC) = sum(measure_memory.([ten.colPtr, ten.rowVal, ten.nzVal]))
-measure_memory(ten::CuSparseMatrixCSR) = sum(measure_memory.([ten.rowPtr, ten.colVal, ten.nzVal]))
+measure_memory(ten::SparseMatrixCSC) =
+    sum(measure_memory.([ten.colptr, ten.rowval, ten.nzval]))
+measure_memory(ten::CuSparseMatrixCSC) =
+    sum(measure_memory.([ten.colPtr, ten.rowVal, ten.nzVal]))
+measure_memory(ten::CuSparseMatrixCSR) =
+    sum(measure_memory.([ten.rowPtr, ten.colVal, ten.nzVal]))
 measure_memory(ten::Diagonal) = measure_memory(diag(ten))
-measure_memory(ten::SiteTensor) = sum(measure_memory.([ten.loc_exp, ]))  # ten.projs...]))
-measure_memory(ten::CentralTensor) = sum(measure_memory.([ten.e11, ten.e12, ten.e21, ten.e22]))
+measure_memory(ten::SiteTensor) = sum(measure_memory.([ten.loc_exp]))  # ten.projs...]))
+measure_memory(ten::CentralTensor) =
+    sum(measure_memory.([ten.e11, ten.e12, ten.e21, ten.e22]))
 measure_memory(ten::DiagonalTensor) = sum(measure_memory.([ten.e1, ten.e2]))
-measure_memory(ten::VirtualTensor) = sum(measure_memory.([ten.con, ]))  # ten.projs...]))
+measure_memory(ten::VirtualTensor) = sum(measure_memory.([ten.con]))  # ten.projs...]))
 measure_memory(ten::MpoTensor) = sum(measure_memory.([ten.top..., ten.ctr, ten.bot...]))
-measure_memory(ten::Union{QMps, QMpo}) = sum(measure_memory.(values(ten.tensors)))
+measure_memory(ten::Union{QMps,QMpo}) = sum(measure_memory.(values(ten.tensors)))
 measure_memory(env::Environment) = sum(measure_memory.(values(env.env)))
 measure_memory(env::EnvironmentMixed) = sum(measure_memory.(values(env.env)))
 measure_memory(lp::PoolOfProjectors) = sum([measure_memory(da) for da ∈ values(lp.data)])
@@ -28,15 +30,18 @@ function format_bytes(bytes, decimals::Int = 2, k::Int = 1024)
     dm = decimals < 0 ? 0 : decimals
     sizes = ["Bytes", "KB", "MB", "GB", "TB", "PB", "EB", "ZB", "YB"]
     i = convert(Int, floor(log(bytes) / log(k)))
-    string(round((bytes / ^(k, i)), digits=dm)) * " " * sizes[i+1]
+    string(round((bytes / ^(k, i)), digits = dm)) * " " * sizes[i+1]
 end
 
-function measure_memory(caches::IdDict{Any, Any}, bytes::Bool = true)
-    memoization_memory = bytes ? Dict{Any, Vector{String}}() : Dict{Any, Vector{Int64}}()
+function measure_memory(caches::IdDict{Any,Any}, bytes::Bool = true)
+    memoization_memory = bytes ? Dict{Any,Vector{String}}() : Dict{Any,Vector{Int64}}()
     for key in keys(caches)
         push!(
-            memoization_memory, key => bytes ? format_bytes.(measure_memory(caches[key])) : measure_memory(caches[key])
-            )
+            memoization_memory,
+            key =>
+                bytes ? format_bytes.(measure_memory(caches[key])) :
+                measure_memory(caches[key]),
+        )
     end
     memoization_memory
 end
diff --git a/src/utils/utils.jl b/src/utils/utils.jl
index e3f39dc..f498204 100644
--- a/src/utils/utils.jl
+++ b/src/utils/utils.jl
@@ -22,7 +22,7 @@ It can be either `:PE` for 'Projector Energy)' order (default) or `:EP` for 'Ene
   - `E`: An array of energy values.
   - `P`: A permutation matrix representing projectors.
 """
-function rank_reveal(energy, order=:PE) #TODO: add type
+function rank_reveal(energy, order = :PE) #TODO: add type
     @assert order ∈ (:PE, :EP)
     dim = order == :PE ? 1 : 2
     E, idx = unique_dims(energy, dim)
@@ -37,14 +37,18 @@ end
 
         # Compute hash for each row
         k = 0
-        @nloops $N i A d->(if d == dim; k = i_d; end) begin
+        @nloops $N i A d -> (
+            if d == dim
+                k = i_d
+            end
+        ) begin
             @inbounds hashes[k] = hash(hashes[k], hash((@nref $N A i)))
         end
 
         # Collect index of first row for each hash
         uniquerow = similar(Array{Int}, axes(A, dim))
         firstrow = Dict{Prehashed,Int}()
-        for k = axes(A, dim)
+        for k in axes(A, dim)
             uniquerow[k] = get!(firstrow, Prehashed(hashes[k]), k)
         end
         uniquerows = collect(values(firstrow))
@@ -52,12 +56,14 @@ end
         # Check for collisions
         collided = falses(axes(A, dim))
         @inbounds begin
-            @nloops $N i A d->(if d == dim
-                k = i_d
-                j_d = uniquerow[k]
-            else
-                j_d = i_d
-            end) begin
+            @nloops $N i A d -> (
+                if d == dim
+                    k = i_d
+                    j_d = uniquerow[k]
+                else
+                    j_d = i_d
+                end
+            ) begin
                 if (@nref $N A j) != (@nref $N A i)
                     collided[k] = true
                 end
@@ -69,7 +75,7 @@ end
             while any(collided)
                 # Collect index of first row for each collided hash
                 empty!(firstrow)
-                for j = axes(A, dim)
+                for j in axes(A, dim)
                     collided[j] || continue
                     uniquerow[j] = get!(firstrow, Prehashed(hashes[j]), j)
                 end
@@ -79,7 +85,7 @@ end
 
                 # Check for collisions
                 fill!(nowcollided, false)
-                @nloops $N i A d->begin
+                @nloops $N i A d -> begin
                     if d == dim
                         k = i_d
                         j_d = uniquerow[k]
@@ -96,6 +102,7 @@ end
             end
         end
 
-        (@nref $N A d->d == dim ? sort!(uniquerows) : (axes(A, d))), indexin(uniquerow, uniquerows)
+        (@nref $N A d -> d == dim ? sort!(uniquerows) : (axes(A, d))),
+        indexin(uniquerow, uniquerows)
     end
-end
\ No newline at end of file
+end
diff --git a/src/variational.jl b/src/variational.jl
index 56e93f7..9af479f 100644
--- a/src/variational.jl
+++ b/src/variational.jl
@@ -3,18 +3,16 @@
 # variational.jl: This file provides basic functions to perform variational compression for MPS.
 #                 If the MPS is moved to the GPU, its compression will be performed on the device.
 
-export
-    variational_compress!,
-    variational_sweep!
+export variational_compress!, variational_sweep!
 
 function variational_compress!(
-    bra::QMps{T}, 
-    mpo::QMpo{T}, 
-    ket::QMps{T}, 
-    tol=1E-10, 
-    max_sweeps::Int=4, 
-    kwargs...
-    ) where T <: Real
+    bra::QMps{T},
+    mpo::QMpo{T},
+    ket::QMps{T},
+    tol = 1E-10,
+    max_sweeps::Int = 4,
+    kwargs...,
+) where {T<:Real}
     @assert is_left_normalized(bra)
     env = Environment(bra, mpo, ket)
     overlap = Inf
@@ -78,10 +76,22 @@ function _right_sweep_var_site!(env::Environment, site::Site; kwargs...)
 end
 
 # TODO those 2 functions are to be removed eventually
-function variational_sweep!(bra::QMps{T}, mpo::QMpo{T}, ket::QMps{T}, ::Val{:left}; kwargs...) where T <: Real
+function variational_sweep!(
+    bra::QMps{T},
+    mpo::QMpo{T},
+    ket::QMps{T},
+    ::Val{:left};
+    kwargs...,
+) where {T<:Real}
     _right_sweep_var!(Environment(bra, mpo, ket); kwargs...)
 end
 
-function variational_sweep!(bra::QMps{T}, mpo::QMpo{T}, ket::QMps{T}, ::Val{:right}; kwargs...) where T <: Real
+function variational_sweep!(
+    bra::QMps{T},
+    mpo::QMpo{T},
+    ket::QMps{T},
+    ::Val{:right};
+    kwargs...,
+) where {T<:Real}
     _left_sweep_var!(Environment(bra, mpo, ket); kwargs...)
 end
diff --git a/src/zipper.jl b/src/zipper.jl
index 9274679..658092d 100644
--- a/src/zipper.jl
+++ b/src/zipper.jl
@@ -1,12 +1,9 @@
-export
-    zipper,
-    corner_matrix,
-    CornerTensor
+export zipper, corner_matrix, CornerTensor
 
- struct CornerTensor{T <: Real}
-    C::Tensor{T, 3}
-    M::MpoTensor{T, 4}
-    B::Tensor{T, 3}
+struct CornerTensor{T<:Real}
+    C::Tensor{T,3}
+    M::MpoTensor{T,4}
+    B::Tensor{T,3}
 
     function CornerTensor(C, M, B)
         T = promote_type(eltype.((C, M, B))...)
@@ -14,12 +11,12 @@ export
     end
 end
 
-struct Adjoint{T, S <: CornerTensor}
+struct Adjoint{T,S<:CornerTensor}
     parent::S
 
-    function Adjoint{T}(ten::CornerTensor{S}) where {T, S}
+    function Adjoint{T}(ten::CornerTensor{S}) where {T,S}
         F = promote_type(T, S)
-        new{F, CornerTensor{F}}(ten)
+        new{F,CornerTensor{F}}(ten)
     end
 end
 
@@ -34,8 +31,8 @@ function zipper(
     iters_var = 1,
     Dtemp_multiplier = 2,
     depth::Int = 0,
-    kwargs...
-    ) where R <: Real
+    kwargs...,
+) where {R<:Real}
     onGPU = ψ.onGPU && ϕ.onGPU
     @assert is_left_normalized(ϕ)
 
@@ -60,18 +57,19 @@ function zipper(
             CM = CornerTensor(C, ψ[i], out[i])
 
             Urs, Srs, Vrs = [], [], []
-            for i in 1:iters_rand
-               Utemp, Stemp, Vtemp = svd_corner_matrix(CM, method, Dtemp, tol; toGPU=false, kwargs...)
-               push!(Urs, Utemp)
-               push!(Srs, Stemp)
-               push!(Vrs, Vtemp)
+            for i = 1:iters_rand
+                Utemp, Stemp, Vtemp =
+                    svd_corner_matrix(CM, method, Dtemp, tol; toGPU = false, kwargs...)
+                push!(Urs, Utemp)
+                push!(Srs, Stemp)
+                push!(Vrs, Vtemp)
             end
 
             Ur = hcat(Urs...)
             Vr = hcat(Vrs...)
             Sr = vcat(Srs...) ./ iters_rand
-            QU, RU = qr_fact(Ur, Dtemp * iters_rand, 0.0; toGPU=false, kwargs...)
-            QV, RV = qr_fact(Vr, Dtemp * iters_rand, 0.0; toGPU=false, kwargs...)
+            QU, RU = qr_fact(Ur, Dtemp * iters_rand, 0.0; toGPU = false, kwargs...)
+            QV, RV = qr_fact(Vr, Dtemp * iters_rand, 0.0; toGPU = false, kwargs...)
             Ur, Sr, Vr = svd_fact(RU * Diagonal(Sr) * RV', Dtemp, tol; kwargs...)
             # Ur = QU * Ur
             Vr = QV * Vr
@@ -80,7 +78,7 @@ function zipper(
                 Vr = CuArray(Vr)
             end
 
-            for _ = 1 : iters_svd
+            for _ = 1:iters_svd
                 # CM * Vr
                 x = reshape(Vr, size(CM.C, 2), size(CM.M, 2), :)
                 x = permutedims(x, (3, 1, 2))
@@ -120,7 +118,7 @@ function zipper(
             C = onGPU ? CUDA.ones(R, 1, 1, 1) : ones(R, 1, 1, 1)
         end
 
-        for _ in 1:iters_var
+        for _ = 1:iters_var
             env.site = i
             update_env_right!(env, i)
             env.C = C
@@ -187,28 +185,39 @@ end
 
 function Base.Array(CM::CornerTensor)  # change name, or be happy with "psvd(Array(Array(CM))"
     B, M, C = CM.B, CM.M, CM.C
-    for v ∈ reverse(M.bot) B = contract_matrix_tensor3(v, B) end
+    for v ∈ reverse(M.bot)
+        B = contract_matrix_tensor3(v, B)
+    end
     Cnew = corner_matrix(C, M.ctr, B)
     @cast Cnew[(t1, t2), t3, t4] := Cnew[t1, t2, t3, t4]
-    for v ∈ reverse(M.top) Cnew = contract_matrix_tensor3(v, Cnew) end
+    for v ∈ reverse(M.top)
+        Cnew = contract_matrix_tensor3(v, Cnew)
+    end
     @cast Cnew[t12, (t3, t4)] := Cnew[t12, t3, t4]
 end
 
-Base.Array(CM::Adjoint{T, CornerTensor{T}}) where T = adjoint(Array(CM.ten))
+Base.Array(CM::Adjoint{T,CornerTensor{T}}) where {T} = adjoint(Array(CM.ten))
 
 # TODO rethink this function
-function svd_corner_matrix(CM, method::Symbol, Dcut::Int, tol::Real; toGPU::Bool=true, kwargs...)
+function svd_corner_matrix(
+    CM,
+    method::Symbol,
+    Dcut::Int,
+    tol::Real;
+    toGPU::Bool = true,
+    kwargs...,
+)
     if method == :svd
         U, Σ, V = svd_fact(Array(Array(CM)), Dcut, tol; kwargs...)
     elseif method == :psvd
-        U, Σ, V = psvd(Array(Array(CM)), rank=Dcut)
+        U, Σ, V = psvd(Array(Array(CM)), rank = Dcut)
     elseif method == :psvd_sparse
-        U, Σ, V = psvd(CM, rank=Dcut)
+        U, Σ, V = psvd(CM, rank = Dcut)
     elseif method == :tsvd
         U, Σ, V = tsvd(Array(CM), Dcut; kwargs...)
     elseif method == :tsvd_sparse
         v0 = 2 .* rand(eltype(CM), size(CM, 1)) .- 1
-        U, Σ, V = tsvd(CM, Dcut, initvec=v0; kwargs...)
+        U, Σ, V = tsvd(CM, Dcut, initvec = v0; kwargs...)
     else
         throw(ArgumentError("Wrong method $method"))
     end
@@ -218,13 +227,15 @@ end
 
 # this is for psvd to work
 LinearAlgebra.ishermitian(ten::CornerTensor) = false
-Base.size(ten::CornerTensor) = (size(ten.B, 1) * size(ten.M, 1), size(ten.C, 2) * size(ten.M, 2))
+Base.size(ten::CornerTensor) =
+    (size(ten.B, 1) * size(ten.M, 1), size(ten.C, 2) * size(ten.M, 2))
 Base.size(ten::CornerTensor, n::Int) = size(ten)[n]
-Base.eltype(ten::CornerTensor{T}) where T = T
-Base.adjoint(ten::CornerTensor{T}) where T = Adjoint{T}(ten)
+Base.eltype(ten::CornerTensor{T}) where {T} = T
+Base.adjoint(ten::CornerTensor{T}) where {T} = Adjoint{T}(ten)
 
 CuArrayifneeded(ten::CornerTensor, x) = typeof(ten.B) <: CuArray ? CuArray(x) : x
-CuArrayifneeded(ten::Adjoint{T, CornerTensor{T}}, x) where T = CuArrayifneeded(ten.parent, x)
+CuArrayifneeded(ten::Adjoint{T,CornerTensor{T}}, x) where {T} =
+    CuArrayifneeded(ten.parent, x)
 
 
 function LinearAlgebra.mul!(y, ten::CornerTensor, x)
@@ -232,18 +243,25 @@ function LinearAlgebra.mul!(y, ten::CornerTensor, x)
     x = reshape(x, size(ten.C, 2), size(ten.M, 2), :)
     x = permutedims(x, (3, 1, 2))
     yp = update_env_right(ten.C, x, ten.M, ten.B)
-    y[:, :] .= Array(reshape(permutedims(yp, (1, 3, 2)), size(ten.B, 1) * size(ten.M, 1), :))
+    y[:, :] .=
+        Array(reshape(permutedims(yp, (1, 3, 2)), size(ten.B, 1) * size(ten.M, 1), :))
 end
 
-function LinearAlgebra.mul!(y, ten::Adjoint{T, CornerTensor{T}}, x) where T <: Real
+function LinearAlgebra.mul!(y, ten::Adjoint{T,CornerTensor{T}}, x) where {T<:Real}
     x = CuArrayifneeded(ten, x)  # CuArray(x)  # TODO this an ugly patch
     x = reshape(x, size(ten.parent.B, 1), size(ten.parent.M, 1), :)
     x = permutedims(x, (1, 3, 2))
     yp = project_ket_on_bra(x, ten.parent.B, ten.parent.M, ten.parent.C)
-    y[:, :] .= Array(reshape(permutedims(yp, (2, 3, 1)), size(ten.parent.C, 2) * size(ten.parent.M, 2), :))
+    y[:, :] .= Array(
+        reshape(
+            permutedims(yp, (2, 3, 1)),
+            size(ten.parent.C, 2) * size(ten.parent.M, 2),
+            :,
+        ),
+    )
 end
 
-function Base.:(*)(ten::CornerTensor{T}, x) where T
+function Base.:(*)(ten::CornerTensor{T}, x) where {T}
     x = CuArrayifneeded(ten, x)  # CuArray(x)  # TODO this an ugly patch
     x = reshape(x, 1, size(ten.C, 2), size(ten.M, 2))
     yp = update_env_right(ten.C, x, ten.M, ten.B)
@@ -251,7 +269,7 @@ function Base.:(*)(ten::CornerTensor{T}, x) where T
     Array(out)  # TODO this an ugly patch
 end
 
-function Base.:(*)(ten::Adjoint{T, CornerTensor{T}}, x) where T <: Real
+function Base.:(*)(ten::Adjoint{T,CornerTensor{T}}, x) where {T<:Real}
     x = CuArrayifneeded(ten, x)  # CuArray(x)  # TODO this an ugly patch
     x = reshape(x, size(ten.parent.B, 1), 1, size(ten.parent.M, 1))
     yp = project_ket_on_bra(x, ten.parent.B, ten.parent.M, ten.parent.C)
diff --git a/test/attic/canonise.jl b/test/attic/canonise.jl
index 17eb5e5..24330a8 100644
--- a/test/attic/canonise.jl
+++ b/test/attic/canonise.jl
@@ -1,7 +1,7 @@
 T = Float64
 D = 16
 
-sites = [1, 3//2, 2, 5//2, 3, 7//2, 4]
+sites = [1, 3 // 2, 2, 5 // 2, 3, 7 // 2, 4]
 d = [1, 2, 2, 2, 4, 2, 2]
 
 id = Dict(j => d[i] for (i, j) in enumerate(sites))
diff --git a/test/attic/compressions.jl b/test/attic/compressions.jl
index 615cc31..a73e460 100644
--- a/test/attic/compressions.jl
+++ b/test/attic/compressions.jl
@@ -10,7 +10,7 @@
     tol = 1E-10
 
     ψ = rand(QMps{T}, sites, D, d)
-    W = rand(QMpo{T}, [1,2,3,4], 2, 4)
+    W = rand(QMpo{T}, [1, 2, 3, 4], 2, 4)
 
     bra = ψ
     ket = ψ
diff --git a/test/attic/contractions.jl b/test/attic/contractions.jl
index 61f0d9a..4f5a3ea 100644
--- a/test/attic/contractions.jl
+++ b/test/attic/contractions.jl
@@ -27,20 +27,20 @@
 
             ϕ.tensors[ψ.sites[1]] *= 1 / norm(ϕ)
             @test dot(ϕ, ϕ) ≈ 1
-        end 
-    end 
-    
+        end
+    end
+
     @testset "dot product of MPS with MPO" begin
-        B = randn(Float64, 4,2,3)
-        A = randn(Float64, 2,2)  
-        C = randn(Float64, 2,2,2,2)
+        B = randn(Float64, 4, 2, 3)
+        A = randn(Float64, 2, 2)
+        C = randn(Float64, 2, 2, 2, 2)
         O2 = random_QMpo(sites, D, d)
 
-        @testset "contract_left gives correct sizes" begin 
-            @test size(contract_left(B, A)) == (4,2,3)
+        @testset "contract_left gives correct sizes" begin
+            @test size(contract_left(B, A)) == (4, 2, 3)
         end
 
-        @testset "contract_tensors43 gives correct sizes" begin 
+        @testset "contract_tensors43 gives correct sizes" begin
             @test size(contract_tensors43(C, B)) == (8, 2, 6)
         end
 
@@ -58,4 +58,4 @@
             @test size(F[2]) == size(G[2]) == (4, 1, 2)
         end
     end
-end
\ No newline at end of file
+end
diff --git a/test/attic/environment.jl b/test/attic/environment.jl
index f8b8b56..d8557cd 100644
--- a/test/attic/environment.jl
+++ b/test/attic/environment.jl
@@ -1,10 +1,10 @@
 @testset "Environment" begin
-    sites = [1, 1//2, 2, 3, 7//2, 4, 5, 6]
+    sites = [1, 1 // 2, 2, 3, 7 // 2, 4, 5, 6]
     site = 3
     @testset "left_nbrs_site gives correct left neighbor of a given site" begin
         @test left_nbrs_site(site, sites) == 2
     end
     @testset "left_nbrs_site gives correct right neighbor of a given site" begin
-        @test right_nbrs_site(site, sites) == 7//2
+        @test right_nbrs_site(site, sites) == 7 // 2
     end
 end
diff --git a/test/attic/linear_algebra_ext.jl b/test/attic/linear_algebra_ext.jl
index 2496135..7bbefce 100644
--- a/test/attic/linear_algebra_ext.jl
+++ b/test/attic/linear_algebra_ext.jl
@@ -11,20 +11,20 @@ using LowRankApprox
 
     δ = min(Dcut, size(Σ1)...)
     U1 = U1[:, 1:δ]
-    Σ1 = Σ1[1:δ] 
+    Σ1 = Σ1[1:δ]
     V1 = V1[:, 1:δ]
 
     U2, Σ2, V2 = svd(a)
 
     δ = min(Dcut, size(Σ2)...)
     U2 = U2[:, 1:δ]
-    Σ2 = Σ2[1:δ] 
+    Σ2 = Σ2[1:δ]
     V2 = V2[:, 1:δ]
-    
+
     r1 = U1 * Diagonal(Σ1) * V1'
     r2 = U2 * Diagonal(Σ2) * V2'
 
-    @test norm(r1-r2) < tol
+    @test norm(r1 - r2) < tol
 end
 
 
@@ -36,11 +36,11 @@ end
 
     a = rand(D, D)
 
-    U1, Σ1, V1 = psvd(a, rank=Dcut, atol=1E-16, rtol=1E-16)
-    U2, Σ2, V2 = psvd(a, rank=Dcut, atol=1E-16, rtol=1E-16)
+    U1, Σ1, V1 = psvd(a, rank = Dcut, atol = 1E-16, rtol = 1E-16)
+    U2, Σ2, V2 = psvd(a, rank = Dcut, atol = 1E-16, rtol = 1E-16)
 
     r1 = U1 * Diagonal(Σ1) * V1'
     r2 = U2 * Diagonal(Σ2) * V2'
 
-    @test norm(r1-r2) < tol
-end
\ No newline at end of file
+    @test norm(r1 - r2) < tol
+end
diff --git a/test/attic/memoization.jl b/test/attic/memoization.jl
index b3fc415..13bbdff 100644
--- a/test/attic/memoization.jl
+++ b/test/attic/memoization.jl
@@ -16,7 +16,7 @@ T = Float64
 
     Random.seed!(69)
 
-    ψ = randn(MPS{T}, sites, D, d) 
+    ψ = randn(MPS{T}, sites, D, d)
     ϕ = copy(ψ)
 
     σ = [1, 2, 3]
@@ -30,8 +30,8 @@ T = Float64
     end
 
     @testset "Cached results are equal to the ones computed during first call" begin
-        @test all(env->env==envs[1], envs)
-    end    
+        @test all(env -> env == envs[1], envs)
+    end
 end
 
 
@@ -40,7 +40,7 @@ end
 
     Random.seed!(42)
 
-    ψ = randn(MPS{T}, sites, D, d) 
+    ψ = randn(MPS{T}, sites, D, d)
     ϕ = copy(ψ)
 
     W = randn(MPO{T}, sites, D, d)
@@ -48,17 +48,15 @@ end
 
     σ = [2, 1, 1, 2, 1]
     η = copy(σ)
-    
-    envs = [
-        right_env(mps, mpo, v) for mps in (ψ, ϕ) for mpo ∈ (W, V) for v ∈ (σ, η)
-    ]
-    
+
+    envs = [right_env(mps, mpo, v) for mps in (ψ, ϕ) for mpo ∈ (W, V) for v ∈ (σ, η)]
+
     @testset "Calls made with equal arguments result in a cache hit" begin
         # Should result in 6 calls total (top-level + 5 recursive)
         @test length(memoize_cache(right_env)) == 6
     end
 
     @testset "Cached results are equal to the ones computed during first call" begin
-        @test all(env->env==envs[1], envs)
+        @test all(env -> env == envs[1], envs)
     end
 end
diff --git a/test/attic/mps.jl b/test/attic/mps.jl
index 4c60a6e..f00a778 100644
--- a/test/attic/mps.jl
+++ b/test/attic/mps.jl
@@ -3,7 +3,7 @@
     T = Float64
     D = 16
 
-    sites = [1, 3//2, 2, 5//2, 3, 7//2, 4]
+    sites = [1, 3 // 2, 2, 5 // 2, 3, 7 // 2, 4]
     d = [1, 2, 2, 2, 4, 2, 2]
 
     id = Dict(j => d[i] for (i, j) in enumerate(sites))
@@ -13,7 +13,7 @@
 
         @testset "has correct number of sites" begin
             @test length(ψ) == maximum(sites)
-            @test size(ψ) == (maximum(sites), )
+            @test size(ψ) == (maximum(sites),)
         end
 
         @testset "has correct type" begin
@@ -26,8 +26,15 @@
 
         @testset "has correct bonds" begin
             @test bond_dimension(ψ) ≈ D
-            @test bond_dimensions(ψ) ==
-                [(1, d[1], D), (D, d[2], D), (D, d[3], D), (D, d[4], D), (D, d[5], D), (D, d[6], D), (D, d[7], 1)]
+            @test bond_dimensions(ψ) == [
+                (1, d[1], D),
+                (D, d[2], D),
+                (D, d[3], D),
+                (D, d[4], D),
+                (D, d[5], D),
+                (D, d[6], D),
+                (D, d[7], 1),
+            ]
             @test verify_bonds(ψ) === nothing
         end
 
diff --git a/test/attic/runtests.jl b/test/attic/runtests.jl
index b4c1515..f43ee30 100644
--- a/test/attic/runtests.jl
+++ b/test/attic/runtests.jl
@@ -12,7 +12,7 @@ my_tests = [
     #"mps.jl",
     "canonise.jl",
     #"environment.jl"
-    ]
+]
 
 
 for my_test in my_tests
diff --git a/test/canonise.jl b/test/canonise.jl
index ab29115..89e44d7 100644
--- a/test/canonise.jl
+++ b/test/canonise.jl
@@ -1,6 +1,6 @@
 D = 16
 
-sites = [1, 3//2, 2, 5//2, 3, 7//2, 4]
+sites = [1, 3 // 2, 2, 5 // 2, 3, 7 // 2, 4]
 d = [1, 2, 2, 2, 4, 2, 2]
 
 id = Dict(j => d[i] for (i, j) in enumerate(sites))
@@ -37,12 +37,15 @@ end
 
 @testset "Measure spectrum($T)" for T in (Float32, Float64)
     svd_mps = TensorMap{T}(
-        1 => T[-0.694389933025747 -0.7195989305943268;;; 0.7195989305943268 -0.6943899330257469],
-        2 => T[0.7071067811865477; 0.0;;; -7.850461536237973e-17; 0.7071067811865477]
+        1 => T[
+            -0.694389933025747 -0.7195989305943268;;;
+            0.7195989305943268 -0.6943899330257469
+        ],
+        2 => T[0.7071067811865477; 0.0;;; -7.850461536237973e-17; 0.7071067811865477],
     )
     ψ = QMps(svd_mps)
     @test is_left_normalized(ψ)
     A = measure_spectrum(ψ)
-    @test A[1] ≈ [1.0] 
-    @test A[2] ≈ [0.7071067811865476, 0.7071067811865475] 
-end
\ No newline at end of file
+    @test A[1] ≈ [1.0]
+    @test A[2] ≈ [0.7071067811865476, 0.7071067811865475]
+end
diff --git a/test/projectors.jl b/test/projectors.jl
index bc42434..d255328 100644
--- a/test/projectors.jl
+++ b/test/projectors.jl
@@ -38,15 +38,15 @@
             k = add_projector!(lp, p1)
             k = add_projector!(lp, p2)
 
-            @test typeof(get_projector!(lp, 2, :CPU)) <: Array{T, 1}
-            @test typeof(get_projector!(lp, 1, :GPU)) <: CuArray{T, 1}
+            @test typeof(get_projector!(lp, 2, :CPU)) <: Array{T,1}
+            @test typeof(get_projector!(lp, 1, :GPU)) <: CuArray{T,1}
             @test length(lp, :GPU) == 1
             @test length(lp, :CPU) == 2
 
-            @test typeof(get_projector!(lp, 1, :GPU)) <: CuArray{T, 1}
+            @test typeof(get_projector!(lp, 1, :GPU)) <: CuArray{T,1}
             @test length(lp, :GPU) == 1
 
-            @test typeof(get_projector!(lp, 2, :GPU)) <: CuArray{T, 1}
+            @test typeof(get_projector!(lp, 2, :GPU)) <: CuArray{T,1}
             @test length(lp, :GPU) == 2
         end
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index a460d90..fb5d8cf 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,11 +9,7 @@ disable_logging(LogLevel(1))
 
 using Test
 
-my_tests = [
-    "canonise.jl",
-    "variational.jl",
-    "projectors.jl"
-    ]
+my_tests = ["canonise.jl", "variational.jl", "projectors.jl"]
 
 for my_test in my_tests
     include(my_test)
diff --git a/test/variational.jl b/test/variational.jl
index f06e628..a4e34fc 100644
--- a/test/variational.jl
+++ b/test/variational.jl
@@ -4,14 +4,29 @@ D2 = ([1, 1], [1, 2], [1, 1], [1, 2])
 S = Float64
 rand_central = rand(CentralTensor{S}, [1, 1, 1, 1, 1, 1, 1, 1])
 map1 = MpoTensor(
-    TensorMap{S}(Dict(-1//2 => rand_central, 0 => rand(SiteTensor{S}, PoolOfProjectors{Integer}(), l, D1)))
-    )
+    TensorMap{S}(
+        Dict(
+            -1 // 2 => rand_central,
+            0 => rand(SiteTensor{S}, PoolOfProjectors{Integer}(), l, D1),
+        ),
+    ),
+)
 map2 = MpoTensor(
-    TensorMap{S}(Dict(-1//2 => rand_central, 0 => rand(SiteTensor{S}, PoolOfProjectors{Integer}(), l, D2)))
-    )
+    TensorMap{S}(
+        Dict(
+            -1 // 2 => rand_central,
+            0 => rand(SiteTensor{S}, PoolOfProjectors{Integer}(), l, D2),
+        ),
+    ),
+)
 map3 = MpoTensor(
-    TensorMap{S}(Dict(-1//2 => rand_central, 0 => rand(SiteTensor{S}, PoolOfProjectors{Integer}(), l, D1)))
-    )
+    TensorMap{S}(
+        Dict(
+            -1 // 2 => rand_central,
+            0 => rand(SiteTensor{S}, PoolOfProjectors{Integer}(), l, D1),
+        ),
+    ),
+)
 mpomap = Dict(1 => map1, 2 => map2, 3 => map3)
 
 D = 2