MatrixProductStates.jl

• MatrixProductStates.jl

This is a package-in-progress in which I am implementing the [[https://en.wikipedia.org/wiki/Density_matrix_renormalization_group][DMRG]] algorithm over matrix product states as explained in Schollwöck's [[https://www.sciencedirect.com/science/article/pii/S0003491610001752][The density-matrix renormalization group in the age of matrix product states]]. A similar project has been undertaken in [[https://github.com/0/LatticeSweeper.jl][LatticeSweeper.jl]].

I'm not longer actively developthing this library, but I think it still has value as an educational resource for those wanting to learn DMRG. Julia made it possible for me to implement this package using syntax which is very close to the math written in the online literature.

To acquire this package, simply open a ~julia~ repl (obtained from https://julialang.org/downloads/) and type #+BEGIN_SRC julia using Pkg; Pkg.add("https://github.com/MasonProtter/MatrixProductStates.jl.git") #+END_SRC

** Example: Transverse Field Ising Model #+HTML: <details><summary>Click me!</summary> #+HTML: <p> Suppose we didn't realize the one dimensional transverse field Ising model was exactly solvable and we wanted to study it with DMRG.

The TFIM Hamiltonian is written #+BEGIN_SRC H = - ∑ᵢ σᶻᵢσᶻᵢ₊₁ - ∑ᵢ g σˣᵢ #+END_SRC which in MPO form can be written as #+BEGIN_SRC H = W¹ W² W³... Wᴸ⁻¹ Wᴸ [ 𝟙 𝟘 𝟘] [ 𝟙 𝟘 𝟘] [ 𝟙 𝟘 𝟘] [ 𝟙 ] = [-gσˣ σᶻ 𝟙] | -σᶻ 𝟘 𝟘| | -σᶻ 𝟘 𝟘| ... | -σᶻ 𝟘 𝟘| |-σᶻ | [-gσˣ σᶻ 𝟙] [-gσˣ σᶻ 𝟙] [-gσˣ σᶻ 𝟙] [-gσˣ] #+END_SRC We can study this Hamiltonian using MatrixProductStates.jl as follows:

First, make a function for generating the Hamiltonian given a coupling strength ~g = h/J~ and a system length ~L~: #+BEGIN_SRC julia using MatrixProductStates

function H_TFIM(g, L) id = [1 0; 0 1] σˣ = [0 1; 1 0] σᶻ = [1 0; 0 -1] W_tnsr = zeros(Complex{Float64}, 3, 3, 2, 2) W_tnsr[1, 1, :, :] = id
W_tnsr[2, 1, :, :] = -σᶻ
W_tnsr[3, 1, :, :] = -g*σˣ W_tnsr[3, 2, :, :] = σᶻ
W_tnsr[3, 3, :, :] = id

return MPO(W_tnsr, L) # MPO will assume that W¹ = W_tnsr[end:end, :, :, :] and Wᴸ = W_tnsr[:, 1:1, :, :]

end #+END_SRC

*** Ground State Suppose we want to know the ground state of this system for ~g=0.8~ and ~L=12~ and we have no idea what the MPS form of the ground state looks like a-priori. #+BEGIN_SRC julia g = 1.1; L = 12;

d = 2; # This is the local Hilbert space dimension for each site Dcut = 100; # This is the maximum bond dimension we'll allow our matrix product state to take

H = H_TFIM(g, L) ψ = randn(MPS{L, Complex{Float64}}, Dcut, d) # Generate a completely randomized matrix product state

ϕ, E₀ = ground_state(ψ, H, quiet=true) #Set quiet to false (the deault) to turn off notifications about the algorithm's progress #+END_SRC We now have the ground state ~ϕ~, and an estimate of it's energy eigenvalue ~E₀~!

Note that 12 sites can be easily studied with far less computational cost as an exact diagonalization, but I didn't want to suggest doing something like ~L=50~ right off the bat since that took ~90 minutes on my machine.

We can make sure that this state's energy matches our estimate: #+BEGIN_SRC julia julia> ϕ' * H * ϕ ≈ E₀ # computing ⟨ϕ|H|ϕ⟩ true #+END_SRC and we can varify that it's approximately an eigenstate: #+BEGIN_SRC julia julia> ϕ' * H * H * ϕ ≈ (ϕ' * H * ϕ)^2 # computing ⟨ϕ| H^2 |ϕ⟩ ≈ (⟨ϕ|H|ϕ⟩)^2 true #+END_SRC

*** Correlators We can take advantage of the ~two_point_correlator~ function to study spin-spin correlations in the TFIM #+BEGIN_SRC julia :exports both using UnicodePlots

σᶻ = [1 0 0 -1]

zz(i, j) = two_point_correlator(i=>σᶻ, j=>σᶻ, 12)

js = 2:12

zzs = [realize(ϕ'*zz(1, j)*ϕ) for j in js] #realize will convert complex numbers with a small imaginary part to real.

lineplot(js, zzs, canvas=DotCanvas, ylim=[0, 1.01], width=80, height=30, ylabel="⟨σᶻ₁σᶻⱼ⟩", xlabel="lattice site j", title="Spin-Spin Correlation for g = $g") #+END_SRC

#+RESULTS: #+BEGIN_EXAMPLE Spin-Spin Correlation for g = 1.1 ┌────────────────────────────────────────────────────────────────────────────────┐ 1.01 │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ ⟨σᶻ₁σᶻⱼ⟩ │ │ │ │ │: │ │ '. │ │ '. │ │ '. │ │ '. │ │ ''. │ │ ''.. │ │ ''... │ │ ''.... │ │ ''''.... │ │ '''''....... │ │ '''''''......... │ │ '''''''''........│ 0 │ │ └────────────────────────────────────────────────────────────────────────────────┘ 2 12 lattice site j #+END_EXAMPLE which shows exponentially decaying correlations in the ground state, as expected for ~g > 1~. We can also redo our calculation in the ordered phase: #+BEGIN_SRC julia :exports both g = 0.8;

H = H_TFIM(g, L)

ϕ, Eₒ = ground_state(ψ, H, quiet=true)

ordered_zzs = [realize(ϕ'*zz(1, j)*ϕ) for j in js]

lineplot(js, realize.(ordered_zzs), canvas=DotCanvas, ylim=[0, 1.01], width=80, height=30, ylabel="⟨σᶻ₁σᶻⱼ⟩", xlabel="lattice site j", title="Spin-Spin Correlation for g = $g") #+END_SRC

#+RESULTS: #+BEGIN_EXAMPLE Spin-Spin Correlation for g = 0.8 ┌────────────────────────────────────────────────────────────────────────────────┐ 1.01 │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │. │ │ ''. │ │ ''.. │ │ '''.... │ ⟨σᶻ₁σᶻⱼ⟩ │ ''''''......... │ │ ''''''''''''........... │ │ '''''''...... │ │ '''.... │