diff --git a/prova.jl b/prova.jl
index 559743f..0ec198b 100644
--- a/prova.jl
+++ b/prova.jl
@@ -1 +1,2 @@
+using Revise
 using Scattensor
\ No newline at end of file
diff --git a/src/Scattensor.jl b/src/Scattensor.jl
index bd6d6d0..586e900 100644
--- a/src/Scattensor.jl
+++ b/src/Scattensor.jl
@@ -2,9 +2,12 @@ module Scattensor
 
 # Write your package code here.
 
-include("prova1.jl")
-include("prova2.jl")
-
-export myfunc1, myfunc2
+include("models.jl")
+include("notation.jl")
+include("bloch.jl")
+include("interpolation.jl")
+include("operators.jl")
+include("tensor_networks.jl")
+include("system.jl")
 
 end
\ No newline at end of file
diff --git a/src/bloch.jl b/src/bloch.jl
new file mode 100644
index 0000000..cc615da
--- /dev/null
+++ b/src/bloch.jl
@@ -0,0 +1,293 @@
+using Scattensor
+using LinearAlgebra
+using Plots
+using Optim
+
+"""
+Returns the complete dispersion relation of a quantum system.
+
+Inputs:
+- `𝐇` is the Hamiltonian of the system;
+- `𝐓` is the translation operator of the system.
+
+Outputs:
+- `k` is the `Vector` of momenta associated to the eigenstates;
+- `ℰ` is the `Vector` of energies associated to the eigenstates;
+- `𝛙` is the `Vector` of eigenvectors (`Vector`s) associated to the eigenstates.
+
+Assumptions:
+- `𝐇` must be Hermitian, i.e. 𝐇 = 𝐇†;
+- `𝐇` must be translationally invariant, i.e. [𝐇,𝐓] = 0.
+"""
+function blochStates(
+    𝐇::Matrix{ComplexF64},
+    𝐓::Matrix{ComplexF64}
+    )::Tuple{Vector{Float64}, Vector{Float64}, Vector{Vector{ComplexF64}}}
+
+    @debug "Computing the bloch states from exact diagonalization..."
+
+    dim = size(𝐇)[1] # The dimension of the Hilbert space
+
+    # We compute the groundstate energy E₀ and the matrix product 𝐇𝐓,
+    # where 𝐇 is shifted by E0 - 1
+    E₀ = eigen(𝐇, permute=true).values[1]
+    𝐇′𝐓 = (𝐇 - E₀ * I + I) * 𝐓
+
+    # We compute all the eigenvectors and eigenvalues of HT
+    (𝜆, 𝛙) = eigen(𝐇′𝐓)
+    Nᵥ = length(𝜆) # The number of eigenvalues
+
+    # Extract the phases and moduli from the eigenvalues 𝜆
+    (k, ℰ) = (zeros(Nᵥ), zeros(Nᵥ))
+    for i in eachindex(𝜆)
+        (k[i], ℰ[i]) = (angle(𝜆[i]), abs(𝜆[i]) + E₀ - 1)
+    end
+
+    𝛙 = [𝛙[:,i] for i in 1:dim]
+
+    return k, ℰ, 𝛙
+end
+
+
+"""
+Plots the dispersion relation of the system.
+
+Inputs:
+- `k` is the `Vector` of momenta of the eigenstates;
+- `ℰ` is the list of energies of the eigenstates;
+- `filename` is the name of the file (with extension) where the plot will be saved;
+- `aspectRatio` is the ratio between the width and the height of the plot (default 1:1).
+"""
+function plotDispersionRelation(
+    k::Vector{Float64},
+    ℰ::Vector{Float64},
+    filename::String = "plot.png", 
+    aspectRatio::Float64 = 1.0)
+
+    @debug "Plotting the dispersion relation..."
+
+    # Plotting the dispersion relation
+    plot(k, ℰ, seriestype=:scatter, markersize=4, legend=false, xlabel="k", ylabel="ℰ")
+    plot!(size=(170,170 * aspectRatio), dpi=300) # Setting the size of the plot
+    savefig(filename) # Saving the plot in a file
+end
+
+
+"""
+Finds the bloch states of the first band.
+
+Inputs:
+- `k` is the `Vector` of momenta of the eigenstates;
+- `ℰ` is the list of energies of the eigenstates;
+- `𝛙` is the `Matrix` of eigenvectors of the eigenstates.
+
+Assumptions:
+- The groundstate must be non-degenerate.
+"""
+function firstBandStates(
+    k::Vector{Float64},
+    ℰ::Vector{Float64},
+    𝛙::Vector{Vector{ComplexF64}},
+    L::Integer
+    )::Tuple{Vector{Float64}, Vector{Float64}, Vector{Vector{ComplexF64}}}
+
+    @debug "Computing the first band states..."
+
+    # Constructing the vector of tuples (k, ℰ, 𝛙)
+    𝓣 = [(k[i], ℰ[i], 𝛙[i]) for i in eachindex(k)]
+
+    # Selecting the groundstate and deleting it
+    j = argmin([t[2] for t in 𝓣]) # the groundstate index for 𝓣
+    deleteat!(𝓣, j) # Deleting the element of states at that index
+
+    # Selecting the first band states
+    # The list of the tuple states of the first band
+    𝓣′ = []
+    while length(𝓣) > 0
+        # The index and momentum of the state with the minimum energy
+        j′ = argmin([t[2] for t in 𝓣])
+        k′ = 𝓣[j′][1]
+        # Filling 𝓣′ with the state with the minimum energy
+        push!(𝓣′, 𝓣[j′])
+        # deleting all the elements with same (similar) momentum
+        𝓣 = [t for t in 𝓣 if abs(k′ - t[1]) > π/L]
+        # we repeat the process until 𝓣 is empty
+    end
+
+    # Transforming the a vector of tuples 𝓣′ into a tuple of vectors 𝓑′
+    # the vector of momenta of the first band
+    k = [t[1] for t in 𝓣′]
+    # the vector of energies of the first band
+    ℰ = [t[2] for t in 𝓣′]
+    # the matrix of eigenvectors of the first band
+    𝛙 = [t[3] for t in 𝓣′]
+
+    # Returning the tuple of bloch states of the first band
+    return k, ℰ, 𝛙
+end
+
+
+""" Plots the energy density profile of the state `𝛙`.
+
+Inputs:
+- `𝛙` is a generic qualtum state;
+- `𝐡₁` is the local Hamiltonian of the first site;
+- `𝐓` is the translation operator of the system;
+- `L` is the number of sites of the chain.
+
+Optional inputs:
+- `ℰ₀` is the groundstate energy of the system (default 0.0).
+
+Outputs:
+- The `Vector{Real}` of the energy density profile of the state `𝛙`.
+"""
+function energyDensityArray(
+    𝛙::Vector{ComplexF64}, 
+    𝐡₁::Matrix{ComplexF64}, 
+    𝐓::Matrix{ComplexF64}, 
+    L::Integer;
+    ℰ₀::Float64 = 0.0,
+    )::Vector{Float64}
+
+    @debug "Computing the energy density profile..."
+
+    # We initialize the array
+    arr::Vector{Float64} = [real(𝛙' * 𝐡₁ * 𝛙 - ℰ₀/L)]
+
+    # The loop to compute all the others energies
+    for _ in 2:L
+        𝛙 = 𝐓' * 𝛙
+        push!(arr, real(𝛙' * 𝐡₁ * 𝛙 - ℰ₀/L))
+    end
+
+    return arr
+end
+
+"""
+Computes the maximally localized from a generic set of states.
+
+Inputs:
+- `𝛙` is the set of states (a `Vector` of `Vector`s);
+- `𝐓` is the translation operator of the system;
+- `𝐡` is the local Hamiltonian of the first site;
+- `jᶜ` is the localization position of the maximally localized state;
+- `ℰ₀` is the groundstate energy of the system;
+- `L` is the number of sites of the chain.
+
+Outputs:
+- The maximally localized state.
+"""
+function find_wannier(
+    𝛙::Vector{Vector{ComplexF64}},
+    𝐓::Matrix{ComplexF64},
+    𝐡₁::Matrix{ComplexF64},
+    jᶜ::Int64,
+    ℰ₀::Float64,
+    L::Int64
+    )::Vector{ComplexF64}
+    
+    @debug "Computing the localized Wannier state..."
+
+    # We define the number of states
+    N = length(𝛙)
+
+    # For each vector in 𝛙 we compute the translated vector by 1 unit
+    𝚿 = [deepcopy(𝛙)]
+    𝛙′ = deepcopy(𝛙)
+    for _ in 1:(L-1)
+        𝛙′ = [𝐓' * 𝜓 for 𝜓 in 𝛙′]
+        push!(𝚿, deepcopy(𝛙′))
+    end
+
+    # For each state in 𝚿 we compute the multiplied by the local Hamiltonian
+    𝐡𝚿 = [deepcopy([𝐡₁ * 𝜙 for 𝜙 in 𝛟]) for 𝛟 in 𝚿]
+
+    # Computing the overlap matrix Λ between States and StatesH
+    Λ = [𝚿[j][α]' * 𝐡𝚿[j][β] for α in 1:N, β in 1:N, j in 1:L]
+
+    # We define the distance function
+    Id(j,j′) = 1
+    χ(j,j′) = ((Int64(j - j′ - (L-1)/2) ↻ Int64(L)) - (L-1)/2 - 1)
+    χ²(j,j′) = χ(j,j′)^2
+
+    # We define the function A₀, A₁ and A₂
+    function A(f::Function, θ::Vector{Float64}, j′::Int64)::ComplexF64
+        S1 = sum(exp(im * (θ[β] - θ[α])) * Λ[α,β,j] * f(j,j′) for j in 1:L, α in 1:N, β in 1:N)
+        S2 = sum(ℰ₀ * f(j,j′) for j in 1:L)
+        return (S1 - S2) / L
+    end
+
+    # We define auxiliary functions
+    A₀(θ, j′) = A(Id, θ, j′)
+    A₁(θ, j′) = A(χ, θ, j′)
+    A₂(θ, j′) = A(χ², θ, j′)
+
+    # We define the functional to minimize
+    function 𝑓(θ::Vector{Float64}, j′::Int64)::Float64
+        return real((A₂(θ, j′) / A₀(θ, j′)) - (A₁(θ, j′) / A₀(θ, j′))^2)
+    end
+
+    # We define the restricted functional, defined putting the last θ argument to 0
+    function 𝑓ᵣ(θr::Vector{Float64})::Float64
+        θ = deepcopy(θr)
+        push!(θ, 0.0)
+        return 𝑓(θ, jᶜ)
+    end
+
+    # We minimize the functional
+    result = optimize(𝑓ᵣ, ones(N-1), NelderMead())
+    θ = result.minimizer
+    push!(θ, 0.0)
+
+    println("The minimum is: ", 𝑓(θ, jᶜ))
+
+    # We compute the maximally localized state
+    𝓌 = 1/√L * sum(exp(im * θ[j]) * 𝛙[j] for j in 1:N)
+
+    # We compute the maximally localized state
+    return 𝓌
+end
+
+
+# function find_localized(
+#     𝛙::Vector{Vector{ComplexF64}},
+#     𝐓::Matrix{ComplexF64},
+#     𝐡₁::Matrix{ComplexF64},
+#     jᶜ::Int64,
+#     L::Int64
+#     )::Vector{ComplexF64}
+    
+#     @debug "Computing the localized Wannier state..."
+
+#     # We define the number of states
+#     N = length(𝛙)
+
+#     # For each vector in 𝛙 we compute the translated vector by 1 unit
+#     𝚿 = [deepcopy(𝛙)]
+#     𝛙′ = deepcopy(𝛙)
+#     for _ in 1:(L-1)
+#         𝛙′ = [𝐓 * 𝜓 for 𝜓 in 𝛙′]
+#         push!(𝚿, deepcopy(𝛙′))
+#     end
+
+#     # For each state in 𝚿 we compute the multiplied by the local Hamiltonian
+#     𝐡𝚿 = [deepcopy([𝐡₁ * 𝜙 for 𝜙 in 𝛟]) for 𝛟 in 𝚿]
+
+#     # We define the distance function
+#     χ(j) = ((Int64(j - jᶜ - (L-1)/2) ↻ Int64(L)) - (L-1)/2 - 1)
+#     χ²(j) = χ(j)^2
+
+#     # Computing the overlap matrix Λ between States and StatesH
+#     Λ = [sum(χ(j)^2 * 𝚿[j][α]' * 𝐡𝚿[j][β] for j in 1:L) for α in 1:N, β in 1:N]
+
+#     # We select the first (the smallest) eigenvector z of Λ
+#     eig = eigen(Λ)
+#     z = eig.vectors[:,1]
+
+#     # We compute the maximally localized state
+#     lf = sum(z[α] * 𝛙[α] for α in 1:N)
+#     lf = lf / norm(lf)
+
+#     # We compute the maximally localized state
+#     return lf
+# end
\ No newline at end of file
diff --git a/src/interpolation.jl b/src/interpolation.jl
new file mode 100644
index 0000000..9595f4d
--- /dev/null
+++ b/src/interpolation.jl
@@ -0,0 +1,95 @@
+using LinearSolve
+
+"""
+Interpolates a state 𝛙₁ from a state 𝛙₀, using a linear combinations of compositions
+of local operators 𝐋ⱼ.
+"""
+function interpolate(
+    𝛙₁::Vector{ComplexF64}, # The vector to interpolate
+    𝛙₀::Vector{ComplexF64}, # The vector from which 𝛙₁ is interpolated
+    𝐓::Matrix{ComplexF64}, # The translation operator of the system
+    j₁::Int64, # The position of the left border of the interpolation range
+    j₂::Int64, # The position of the right border of the interpolation range
+    𝐋::Vector{Matrix{ComplexF64}}, # A list of d local operators
+    d::Int64, # The dimension of the local space
+    L::Int64, # The number of sites of the chain
+)
+    # We define the identity
+    Ide = Matrix{ComplexF64}(I, d, d)
+
+    function generate_𝐋₁(n)
+        𝐀 = 𝐋[n]
+        for _ in 2:L
+            𝐀 = 𝐀 ⊗ Ide
+        end
+        return 𝐀
+    end
+
+    𝐋₁ = [generate_𝐋₁(n) for n in 1:d]
+
+    # Apply 𝐋ⱼ to the state 𝛙 in the most efficient way (translating the state instead 
+    # of the operator)
+    function apply(n, j, 𝛙)
+        𝛟 = deepcopy(𝛙)
+        for _ in 1:j
+            𝛟 = 𝐓' * 𝛟
+        end
+        𝛟 = 𝐋₁[n+1] * 𝛟
+        for _ in 1:j
+            𝛟 = 𝐓 * 𝛟
+        end
+        return 𝛟
+    end
+
+    function apply_to_array(n, j, arr)
+        arr′ = deepcopy(arr)
+        arr′[j] = n
+        return arr′
+    end
+
+    𝚿 = [𝛙₀]
+    A = [zeros(L)]
+
+    # Find all the possible combinations of the operators starting from 𝛙₀
+    # We also find an array representing these applications
+    for j in j₁:j₂
+        𝚿 = [apply(n,j,𝛙) for n in 0:(d-1) for 𝛙 in 𝚿]
+        A = [apply_to_array(n,j,a) for n in 0:(d-1) for a in A]
+    end
+
+    # We compute the Hessian matrix
+    Hs = [𝚿[α]' * 𝚿[β] for α in eachindex(𝚿), β in eachindex(𝚿)]
+    b = [𝚿[α]' * 𝛙₁ for α in eachindex(𝚿)]
+    println("Dimension of the Hessian matrix: ", size(Hs))
+    println("Rank of the Hessian matrix: ", rank(Hs))
+
+    # We compute the coefficients, calculating the pseudo-inverse of 𝚿
+    println("Computing the inverse of the Hessian matrix...")
+    c = pinv(Hs) * b
+
+    # We compute the infidelity
+    iFi = 1 + sum(c[α]' * c[β] * Hs[α,β] for α in eachindex(c), β in eachindex(c))
+                 - sum((c[α]' * b[α] + c[α] * b[α]') for α in eachindex(c))
+    println("Infidelity: ", iFi)
+
+    # We create a list of the tuples (coefficients, A)
+    listt = [(c[α], A[α]) for α in eachindex(c)]
+
+    # We sort the list by the coefficients
+    sort!(listt, by = x->abs(x[1]), rev = true)
+
+    # We print the first 20 elements of listt in a nice way
+    # println()
+    # println("The first 20 elements are:")
+    # println()
+    # for α in 1:20
+        # println(listt[α][2], " ", abs(listt[α][1]))
+    # end
+
+    # we create the array of the abs of the coefficients sorted
+    xax = [i for i in 1:length(listt)]
+    cp = [abs(x[1]) for x in listt]
+
+    plot(xax, cp, xaxis=:log, yaxis=:log)
+    savefig("data/coefficients.png")
+end
\ No newline at end of file
diff --git a/src/models.jl b/src/models.jl
new file mode 100644
index 0000000..51b19de
--- /dev/null
+++ b/src/models.jl
@@ -0,0 +1,19 @@
+"""Abstract type for a general model."""
+abstract type Model end
+
+"""Mutable struct for the Ising Model."""
+mutable struct IsingModel <: Model
+    spin_interaction    ::Real
+    transverse_field    ::Real
+    longitudinal_field  ::Real
+end
+
+"""Mutable struct for the Heisenberg Model."""
+mutable struct HeisenbergModel <: Model
+    x_spin_interaction  ::Real
+    y_spin_interaction  ::Real
+    z_spin_interaction  ::Real
+    x_field             ::Real
+    y_field             ::Real
+    z_field             ::Real
+end
\ No newline at end of file
diff --git a/src/notation.jl b/src/notation.jl
new file mode 100644
index 0000000..3403e24
--- /dev/null
+++ b/src/notation.jl
@@ -0,0 +1,12 @@
+using Scattensor
+using LinearAlgebra
+
+"""A shortcut binary notation for the Krnoecher product."""
+function ⊗(A, B)
+    return kron(A,B)
+end
+
+"""A shortcut binary notation for the periodic modulus."""
+function ↻(n::Integer, m::Integer)::Integer
+    n > 0 ? (n-1)%m + 1 : m + n%m
+end
\ No newline at end of file
diff --git a/src/operators.jl b/src/operators.jl
new file mode 100644
index 0000000..c468b91
--- /dev/null
+++ b/src/operators.jl
@@ -0,0 +1,164 @@
+using Scattensor
+
+module Operators
+
+using Scattensor
+using LinearAlgebra
+
+# Shortcuts for the module
+
+𝛔ˣ::Matrix{ComplexF64} = [0 1; 1 0]
+𝛔ᶻ::Matrix{ComplexF64} = [1 0; 0 -1]
+
+𝒪(args::Pair{Matrix{ComplexF64}, Int}...) = local_operator(𝒮, args...)
+
+"""
+Generates the translation operator T for a chain of L sites with local dimension d.
+The system is assumed to be uniform, i.e. the local dimension is the same for all sites.
+The system, in order to perform a translation, must be in periodic boundary conditions.
+
+Inputs:
+- `L` is the number of sites of the chain.
+
+Outputs:
+- The translation operator `T`.
+"""
+function translation_operator(
+    L::Integer, # The length of the system
+    d::Integer # The local dimension
+    )::Matrix{ComplexF64}
+
+    N::Integer = d^L
+    𝐓::Matrix{ComplexF64} = zeros(N,N)
+
+    Lst = []
+    c = 0
+    for i ∈ 1:d
+        lst = []
+        for j in 1:(N/d)
+            c = c + 1
+            push!(lst, c)
+        end
+        push!(Lst, lst)
+    end
+
+    print(Lst)
+
+    for indL in eachindex(Lst)
+        lst = Lst[indL]
+        for ind in eachindex(lst)
+            j = lst[ind]
+            𝐓[j, ((d*(j-1)+1)%N) + indL - 1] = 1
+        end
+    end
+
+    return 𝐓
+end
+
+"""
+Generates the whole matrix corresponding to the action of a product of local operators.
+
+Inputs:
+- `𝒮` is the quantum system;
+- args is a list of pairs (𝐚, j) where 𝐚 is the local operator written in the local 
+space (small matrix) and j is the position of the local operator 𝐚.
+"""
+function local_operator(𝒮::ExactDiagSystem, args::Pair{Matrix{ComplexF64}, Int}...)::Matrix{ComplexF64}
+    L = 𝒮.system_size
+    d = 𝒮.local_dimension
+    N = d^L
+    # We assert that all the matrices have the same dimension of the local space
+    for arg in args
+        @assert (size(arg.first) == (d, d)) "The local operator must have the same dimension of the local space."
+    end
+    # For each arg.second we take the mod L using ↻L
+    for arg in args
+        arg.second = arg.second ↻ L
+    end
+    # We consider the dxd identity matrix
+    𝟙 = Matrix{ComplexF64}(I, d, d)
+    # We create a Vector 𝒱 of Vector of Matrix{ComplexF64, where 𝒱[i] is the vector 
+    # which contains all the operators of the i-th site
+    𝒱 = []
+    for i in 1:L
+        𝒱i = []
+        for arg in args
+            if arg.second == i
+                push!(𝒱i, arg.first)
+            end
+        end
+        # If 𝒱i has no elements, we add the identity
+        if length(𝒱i) == 0
+            push!(𝒱i, 𝟙)
+        end
+        push!(𝒱, 𝒱i)
+    end
+    # We substitute each 𝒱[i] with a Matrix{ComplexF64} which corresponds to the product
+    # of all the operators in 𝒱[i]
+    for i in 1:L
+        𝒱[i] = reduce(*, 𝒱[i])
+    end
+    @assert length(𝒱) == L "Something went wrong in function local_operator."
+    # We generate a tensor product of L dxd matrices: each element of this tensor from 𝒱[i]
+    # is a factor of the tensor product (kronecker product ⊗)
+    𝐀 = reduce(⊗, 𝒱)
+    @assert size(𝐀) == (N, N) "Something went wrong in function local_operator."
+    # We return the matrix 𝐀
+    return 𝐀
+end
+
+
+"""
+Generates the Hamiltonian matrix of the Ising model with transverse and longitudinal fields.
+
+Inputs:
+- `J` is the coupling constant of the spins;
+- `hˣ` is the transverse field;
+- `hᶻ` is the longitudinal field;
+- `L` is the number of sites of the chain.
+"""
+function ising_hamiltonian(
+    𝒮::ExactDiagSystem, # The quantum system
+    J::Real, # The coupling constant of the spin interaction
+    hˣ::Real, # The transverse field
+    hᶻ::Real # The longitudinal field
+    )::Matrix{ComplexF64}
+
+    d = 𝒮.system_size
+    L = 𝒮.local_dimension
+
+    @assert (d == 2) "The local dimension must be 2 in the Ising Model."
+
+    return sum(𝒪(-J/2 * 𝛔ᶻ,i,𝛔ᶻ,i+1) - 𝒪(hˣ*𝛔ˣ + hᶻ*𝛔ᶻ,i) for i in 1:L)
+end
+
+"""
+Generates the matrix corresponding to the local Hamiltonian of the Ising model with 
+transverse and longitudinal fields.
+
+Inputs:
+- `j` is the position of the local Hamiltonian;
+- `J` is the coupling constant of the spins;
+- `hˣ` is the transverse field;
+- `hᶻ` is the longitudinal field;
+- `L` is the number of sites of the chain.
+"""
+function local_hamiltonian(
+    𝒮::ExactDiagSystem, # The quantum system
+    ℳ::IsingModel, # An instance of IsingModel
+    j::Integer, # The position of the local Hamiltonian
+    )::Matrix{ComplexF64}
+
+    d = 𝒮.local_dimension
+    L = 𝒮.system_size
+    N = d^L
+    J = ℳ.spin_interaction
+    hˣ = ℳ.transverse_field
+    hᶻ = ℳ.longitudinal_field
+
+    @assert (d == 2) "The local dimension must be 2 in the Ising Model."
+
+    return 𝒪(-J/4 * 𝛔ᶻ,j-1,𝛔ᶻ,j) - 𝒪(J/4 * 𝛔ᶻ,j,𝛔ᶻ,j+1) - 𝒪(hˣ * 𝛔ˣ + hᶻ * 𝛔ᶻ,j)
+end
+
+end # module Operators
\ No newline at end of file
diff --git a/src/prova1.jl b/src/prova1.jl
deleted file mode 100644
index 67c7aac..0000000
--- a/src/prova1.jl
+++ /dev/null
@@ -1,4 +0,0 @@
-"A nice function 1"
-function myfunc1(x::Int)
-    return 1 * x
-end
\ No newline at end of file
diff --git a/src/prova2.jl b/src/prova2.jl
deleted file mode 100644
index 31ebbc8..0000000
--- a/src/prova2.jl
+++ /dev/null
@@ -1,4 +0,0 @@
-"A nice function 2"
-function myfunc2(x::Int)
-    return 2 * x
-end
\ No newline at end of file
diff --git a/src/system.jl b/src/system.jl
new file mode 100644
index 0000000..5e55335
--- /dev/null
+++ b/src/system.jl
@@ -0,0 +1,52 @@
+"Struct for a quantum many body system on a uniform, on a finite one-dimensional 
+lattice with periodic boundary conditions. This system uses matrices, hence
+it is not suitable for large systems. It is implemented to perform exact
+diagonalization algorithms"
+
+mutable struct ExactDiagSystem
+    # L, the number of sites of the chain
+    system_size                 :: Int64
+    # d, The local dimension
+    local_dimension             :: Int64
+    # 𝐇, the Hamiltonian operator
+    hamiltonian_operator        :: Union{Matrix{ComplexF64}, Missing}
+    # 𝐓, the translation operator
+    translation_operator        :: Matrix{ComplexF64}
+    # Whether [H,T] = 0
+    translational_invariance    :: Union{Bool, Missing} 
+    # 𝛙[i], the eigenstates of the system
+    eigenstates                 :: Union{Vector{Vector{ComplexF64}}, Missing}
+    # ℰ[i], the eigenenergies of the system
+    eigenstates_energies        :: Union{Vector{Float64}, Missing}
+    # k[i], the eigenmomenta of the system
+    eigenstates_momenta         :: Union{Vector{Float64}, Missing}
+    # ℰ₀, the groundstate of the system
+    groundstate                 :: Union{Vector{ComplexF64}, Missing}
+    # The groundstate energy of the system
+    groundstate_energy          :: Union{Float64, Missing}
+    # The Bloch states of the first band
+    first_band_states           :: Union{Vector{Vector{ComplexF64}}, Missing}
+end
+
+function ExactDiagSystem(L::Integer, d::Integer)::ExactDiagSystem
+    𝐓::Matrix{ComplexF64} = translation_operator(L, d)
+    𝒮 = ExactDiagSystem(L,d,missing,𝐓,missing,missing,missing,missing,missing,missing,missing)
+    return 𝒮
+end
+
+presence_𝐇(𝒮::ExactDiagSystem) = !ismissing(𝒮.hamiltonian_operator)
+
+presence_eigenstates(𝒮::ExactDiagSystem) = !ismissing(𝒮.eigenstates)
+
+check_translational_invariance(𝒮::ExactDiagSystem)
+    @assert presence_𝐇(𝒮) "𝐇 not already inserted."
+    𝐇::Matrix{ComplexF64} = 𝒮.hamiltonian_operator
+    𝐓::Matrix{ComplexF64} = 𝒮.translation_operator
+    if (𝐇 * 𝐓 == 𝐓 * 𝐇)
+        @warn "Found hamiltonian not translational invariant."
+        𝒮.translational_invariance = false
+    else
+        @debug "Translational invariance of 𝐇 checked."
+        𝒮.translational_invariance = true 
+    end
+end
\ No newline at end of file
diff --git a/src/tensor_networks.jl b/src/tensor_networks.jl
new file mode 100644
index 0000000..a96888e
--- /dev/null
+++ b/src/tensor_networks.jl
@@ -0,0 +1,196 @@
+using Scattensor
+using ITensors
+using ITensorTDVP
+using ITensorGLMakie
+using Plots
+
+# Electric Hamiltonian Opsum
+function OsHE()::OpSum
+    H = OpSum()
+    for j in 1:L
+        H += -1,"Sz",j
+    end
+    for j in 1:(L-1)
+        H += -2,"Sz",j,"Sz",j+1
+    end
+    return H
+end
+
+# Magnetic Hamiltonian Opsum
+function OsHB()::OpSum
+    H = OpSum()
+    for j in 1:L
+        H += -√8,"Sx",j
+    end
+    return H
+end
+
+# Electric impurity Opsum
+function OsHEimp(p::Integer)::OpSum
+    H = OpSum()
+    H += -1,"Sz",p
+    H += -1,"Sz",p+1
+    H += 5/4,"Sz",p,"Sz",p+1
+    return H
+end
+
+# Magnetic impurity Opsum
+function OsHBimp(p::Integer)::OpSum
+    H = OpSum()
+    c1 = √6 + 2 - √8
+    c2 = √6 + 2
+    H += -c1,"Sx",p
+    H += -c1,"Sx",p+1
+    H += -c2,"Sz",p,"Sx",p+1
+    H += -c2,"Sx",p,"Sz",p+1
+    return H
+end
+
+# Local electric Hamiltonian Opsum
+function OsHE(j::Integer)::OpSum
+    H = OpSum()
+    H += -1,"Sz",j
+    H += -1,"Sz",j-1,"Sz",j
+    H += -1,"Sz",j,"Sz",j+1
+    return H
+end
+
+# Local magnetic Hamiltonian Opsum
+function OsHB(j::Integer)::OpSum
+    H = OpSum()
+    H += -√8,"Sx",j
+    return H
+end
+
+# Naive Wannier creation operator OpSum
+function OsNW(j₀::Integer)::OpSum
+    H = OpSum()
+    H += "Sz",j₀
+    H += "Sx",j₀
+    return H
+end
+
+# Naive wavepacket creation operator OpSum
+function OsNWp(j₀::Integer, k::Float64, σ::Float64, L::Integer, ϵ::Float64)::OpSum
+    H = OpSum()
+    f(j) = exp(im * k * j) * exp(-(j - j₀)^2 / (4 * σ^2))
+    for j in 1:L
+        if abs(f(j)) > ϵ
+            H += f(j),"Sz",j
+            H += f(j),"Sx",j
+        end
+    end
+    return H
+end
+
+# Naive wavepacket creation operator OpSum
+
+println("TNs stuff.")
+
+# System parameters
+L = 100 # Number of sites
+λ = 2 # Coupling constant
+jⁱᵐᵖ = Integer(round(L/2)) # Impurity position
+
+# Wavepacket parameters
+jᵂᴾ = Integer(round(L/4)) # Wavevector initial average position
+kᵂᴾ = π/2 # Wavevector initial average momentum
+σᵂᴾ = 2.0 # Wavevector initial standard deviation
+ϵᵂᴾ = 10^(-10) # Wavevector cutoff for the creation operator
+
+# DMRG parameters
+Nᴰᴹᴿᴳ = 10 # Number of sweeps
+χᴰᴹᴿᴳ = 100 # Maximum bond dimension
+ϵᴰᴹᴿᴳ = 1e-10 # Cutoff
+
+# Simulation parameters
+χᵀᴰⱽᴾ = 100 # Maximum bond dimension for the time evolution
+ϵᵀᴰⱽᴾ = 1e-10 # Cutoff for the time evolution
+dt = 0.5 # Time step
+Δt = 130.0 # Total time of simulation
+
+# We create the sites and the Hamiltonian
+os = λ * OsHE() + (1 - λ) * OsHB() + (λ * OsHEimp(jⁱᵐᵖ) + (1 - λ) * OsHBimp(jⁱᵐᵖ))
+sites = siteinds("S=1/2", L)
+H = MPO(os, sites)
+
+# We find the groundstate
+ℰ₀,ψ₀ = dmrg(H, randomMPS(sites); 
+            nsweeps = Nᴰᴹᴿᴳ, 
+            maxdim = χᴰᴹᴿᴳ, 
+            cutoff = ϵᵀᴰⱽᴾ)
+
+# We generate the creation operator and we apply it to the groundstate
+Loc = MPO(OsNWp(jᵂᴾ, kᵂᴾ, σᵂᴾ, L, ϵᵂᴾ), sites)
+ψ = noprime(Loc * ψ₀)
+normalize!(ψ)
+
+# We time-evolve the system
+function tdvp_time_evolution(
+    H::MPO,
+    ψ::MPS,
+    ψ₀::MPS,
+    dt::Float64,
+    Δt::Float64
+    )
+
+    N = Integer(Δt / dt) # Number of steps
+
+    # We create the list of times 
+    times = [n * dt for n in 1:N]
+    positions = [x for x in 2:L-1]
+
+    # We create a 2D array E of dimensions (N, L-2)
+    En = zeros(N, L-2)
+    EnLog = zeros(N, L-2)
+    Linkdims = zeros(N, L-2)
+    Maxlinkdim = zeros(N)
+
+    for n in 1:N
+        ψ = tdvp(H, ψ, -im * dt, maxlinkdim = 100, cutoff = 1e-10)
+        normalize!(ψ)
+
+        timeElapsed = @elapsed begin
+        for j in 2:(L-1)
+            locEn = λ * OsHE(j) + (1 - λ) * OsHB(j)
+            A = MPO(locEn, sites)
+            ComputedLocEn = real(inner(ψ', A, ψ)) - real(inner(ψ₀', A, ψ₀))
+            En[n,j-1] = ComputedLocEn
+            EnLog[n,j-1] = log10(abs(ComputedLocEn))
+            Linkdims[n,j-1] = linkdim(ψ,j)
+            Maxlinkdim[n] = maxlinkdim(ψ)
+        end
+        end
+
+        # Print information
+        println("Step ", n, " of ", N, " completed.")
+        println("Maximum link dimension of MPS: ", maxlinkdim(ψ))
+        println("Relative Energy: ", real(inner(ψ', H, ψ) - inner(ψ₀', H, ψ₀)))
+        println("Estimated time remained: ", timeElapsed * (N - n), " seconds.")
+    end
+
+    # heatmap(positions, times, En)
+    # savefig("simulation_output/En.png")
+
+    # heatmap(positions, times, EnLog)
+    # savefig("simulation_output/EnLog.png")
+
+    # heatmap(positions, times, Linkdims)
+    # savefig("simulation_output/Linkdims.png")
+
+    titleString = "L = $L, λ = $λ, jⁱᵐᵖ = $jⁱᵐᵖ, jᵂᴾ = $jᵂᴾ, kᵂᴾ = $kᵂᴾ,
+                    σᵂᴾ = $σᵂᴾ, ϵᵂᴾ = $ϵᵂᴾ, Nᴰᴹᴿᴳ = $Nᴰᴹᴿᴳ, χᴰᴹᴿᴳ ≤ $χᴰᴹᴿᴳ,
+                    ϵᴰᴹᴿᴳ = $ϵᴰᴹᴿᴳ, χᵀᴰⱽᴾ ≤ $χᵀᴰⱽᴾ, ϵᵀᴰⱽᴾ = $ϵᵀᴰⱽᴾ, dt = $dt, Δt = $Δt"
+
+    l = @layout [grid(2,2); b{0.2h}]
+
+    title = plot(title = titleString, grid = false, showaxis = false)
+    p1 = heatmap(positions, times, En, title = "Energy density")
+    p2 = heatmap(positions, times, EnLog, title = "Log10 of Energy dens.")
+    p3 = heatmap(positions, times, Linkdims, title = "Link dims")
+    p4 = plot(times, Maxlinkdim, title = "Max link dims")
+    plot(p1, p2, p3, p4, title, layout = l, dpi = 300, size=(1000, 1000))
+    savefig("simulation_output/All.png")
+end
+
+tdvp_time_evolution(H, ψ, ψ₀, dt, Δt)
\ No newline at end of file