diff --git a/Hamiltonian.jl b/Hamiltonian.jl
index bfaed69..8697dd2 100644
--- a/Hamiltonian.jl
+++ b/Hamiltonian.jl
@@ -10,26 +10,27 @@ struct Hamiltonian{T}
     N::Int
     L::T
     μ::T
-    ∂1 # Matrix{Complex{T}} or Nothing
-    K_diag # CuTensor{Complex{T}} or Nothing
-    K_mixed # CuTensor{Complex{T}} or Nothing
-    Vs # Array{Complex{T}} or CuArray{Complex{T}}
+    K_partial::Matrix{Complex{T}}
+    K_diag::Union{CuTensor{Complex{T}},Nothing}
+    K_mixed::Union{CuTensor{Complex{T}},Nothing}
+    Vs::Union{Array{Complex{T}},CuArray{Complex{T}}}
     hermitian::Bool
     mode::Hamiltonian_backend
+
     function Hamiltonian{T}(V_twobody::Function, d::Int, n::Int, N::Int, L::T, ϕ::T, μ::T, n_image::Int, mode::Hamiltonian_backend) where {T<:Float}
         @assert mode != gpu_cutensor || CUDA.functional() && CUDA.has_cuda() && CUDA.has_cuda_gpu() "CUDA not available"
         k = -N÷2:N÷2-1
         Vs = calculate_Vs(V_twobody, d, n, N, L, ϕ, n_image)
         hermitian = ϕ == 0.0
-        if mode == cpu_tensor
-            ∂1 = exp(-im * ϕ) .* ∂_1DOF.(L, N, k, k')
-            return new{T}(d, n, N, L, μ, ∂1, nothing, nothing, Vs, hermitian, mode)                
-        elseif mode == gpu_cutensor
-            K_partial = (exp(-im * ϕ) * im / sqrt(2 * μ)) .* ∂_1DOF.(L, N, k, k')
+        K_partial = (exp(-im * ϕ) * im / sqrt(2 * μ)) .* ∂_1DOF.(L, N, k, k')
+        K_diag = nothing
+        K_mixed = nothing
+        if mode == gpu_cutensor
             K_diag = CuTensor(CuArray(K_partial * K_partial), ['a', 'A'])
             K_mixed = CuTensor(CuArray(K_partial), ['a', 'A']) * CuTensor(CuArray(K_partial), ['b', 'B'])
-            return new{T}(d, n, N, L, μ, nothing, K_diag, K_mixed, CuArray(Vs), hermitian, mode)
+            Vs = CuArray(Vs)
         end
+        return new{T}(d, n, N, L, μ, K_partial, K_diag, K_mixed, Vs, hermitian, mode)
     end
 end
 
@@ -45,7 +46,6 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
     # apply V operator
     @. out = H.Vs * v
     # apply K opereator
-    coeff = -1 / (2 * H.μ)
     coords = H.n - 1
     nconList_v_template = -collect(1:H.d*(coords))
     for dim = 1:H.d
@@ -62,8 +62,8 @@ function LinearAlgebra.mul!(out::Array{Complex{T}}, H::Hamiltonian{T}, v::Array{
                     nconList_v[i1] = 1
                 end
                 nconList_v[i2] = 2
-                v_new = @ncon((H.∂1, H.∂1, v), (nconList_1, nconList_2, nconList_v))
-                out = axpy!(coeff, v_new, out)
+                v_new = @ncon((H.K_partial, H.K_partial, v), (nconList_1, nconList_2, nconList_v))
+                out = axpy!(1, v_new, out)
             end
         end
     end