DMRG Tutorial
The density matrix renormalization group (DMRG) is an algorithm for computing eigenstates of Hamiltonians (or extremal eigenvectors of large, Hermitian matrices). It computes these eigenstates in the matrix product state (MPS) format.
Let's see how to set up and run a DMRG calculation using the ITensor library. We will be interested in finding the ground state of the quantum Hamiltonian $H$ given by:
\[H = \sum_{j=1}^{N-1} \mathbf{S}_{j} \cdot \mathbf{S}_{j+1} = \sum_{j=1}^{N-1} S^z_{j} S^z_{j+1} + \frac{1}{2} S^+_{j} S^-_{j+1} + \frac{1}{2} S^-_{j} S^+_{j+1}\]
This Hamiltonian is known as the one-dimensional Heisenberg model and we will take the spins to be $S=1$ spins (spin-one spins). We will consider the case of $N=100$ and plan to do five sweeps of DMRG (five passes over the system).
ITensor DMRG Code
Let's look at an entire, working ITensor code that will do this calculation then discuss the main steps. If you need help running the code below, see the getting started page on Running ITensor and Julia Codes.
using ITensors, ITensorMPS
let
N = 100
sites = siteinds("S=1",N)
os = OpSum()
for j=1:N-1
os += "Sz",j,"Sz",j+1
os += 1/2,"S+",j,"S-",j+1
os += 1/2,"S-",j,"S+",j+1
end
H = MPO(os,sites)
psi0 = random_mps(sites;linkdims=10)
nsweeps = 5
maxdim = [10,20,100,100,200]
cutoff = [1E-10]
energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff)
return
end
Steps of The Code
The first two lines
using ITensors, ITensorMPS
N = 100
sites = siteinds("S=1",N)
tells the function siteinds
to make an array of ITensor Index objects which have the properties of $S=1$ spins. This means their dimension will be 3 and they will carry the "S=1"
tag, which will enable the next part of the code to know how to make appropriate operators for them.
Try printing out some of these indices to verify their properties:
@show sites[1]
(dim=3|id=482|"S=1,Site,n=1")
The next part of the code builds the Hamiltonian:
os = OpSum()
for j=1:N-1
os += "Sz",j,"Sz",j+1
os += 1/2,"S+",j,"S-",j+1
os += 1/2,"S-",j,"S+",j+1
end
H = MPO(os,sites)
An OpSum
is an object which accumulates Hamiltonian terms such as "Sz",1,"Sz",2
so that they can be summed afterward into a matrix product operator (MPO) tensor network. The line of code H = MPO(os,sites)
constructs the Hamiltonian in the MPO format, with physical indices given by the array sites
.
The line
psi0 = random_mps(sites;linkdims=10)
constructs an MPS psi0
which has the physical indices sites
and a bond dimension of 10. It is made by a random quantum circuit that is reshaped into an MPS, so that it will have as generic and unbiased properties as an MPS of that size can have. This choice can help prevent the DMRG calculation from getting stuck in a local minimum.
The lines
nsweeps = 5
maxdim = [10,20,100,100,200]
cutoff = [1E-10]
define the number of DMRG sweeps (five) we will instruct the code to do, as well as the parameters that will control the speed and accuracy of the DMRG algorithm within each sweep. The array maxdim
limits the maximum MPS bond dimension allowed during each sweep and cutoff
defines the truncation error goal of each sweep (if fewer values are specified than sweeps, the last value is used for all remaining sweeps).
Finally the call
energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff)
runs the DMRG algorithm included in ITensor, using psi0
as an initial guess for the ground state wavefunction. The optimized MPS psi
and its eigenvalue energy
are returned.
After the dmrg
function returns, you can take the returned MPS psi
and do further calculations with it, such as measuring local operators or computing entanglement entropy.