Exercises based on the toy codes
Uncomment and run the cells below (removing the #) when running this notebook on https://colab.research.google.com.
Alternatively, you can run this notebook locally with jupyter, provided that you have the toycodes subfolder from https://github.com/tenpy/tenpy_toycodes in the same folder as your notebook.
[ ]:
#!pip install git+https://github.com/tenpy/tenpy_toycodes.git
# use `pip uninstall tenpy-toycodes` to remove it a gain.
This tutorial focuses on a set of toy codes (using only python with numpy + scipy) that provide a simple implementation of the various MPS algorithms.
You can add your code below by inserting additional cells as neccessary and running them (press Shift+Enter).
DISCLAIMER: the toy codes used are not optimized, and we only use very small bond dimensions here. For state-of-the-art MPS calculations (especially for cylinders towards 2D), chi should be significantly larger, often on the order of several 1000s.
[ ]:
[ ]:
import numpy as np
import scipy
import matplotlib.pyplot as plt
np.set_printoptions(precision=5, suppress=True, linewidth=100)
plt.rcParams['figure.dpi'] = 150
toy codes: a_mps.py
The file a_mps.py defines a SimpleMPS class, that provides methods for expectation values and the entanglement entropy. You can initialize an inital product state MPS with the provided functions init_FM_MPS or init_Neel_MPS:
[ ]:
from tenpy_toycodes.a_mps import init_FM_MPS, init_Neel_MPS, SimpleMPS
[ ]:
psi_FM = init_FM_MPS(L=10, d=2, bc='finite')
print(psi_FM)
SigmaZ = np.diag([1., -1.])
print(psi_FM.site_expectation_value(SigmaZ))
Exercise
Initialize a Neel state MPS. Print the SigmaZ expectation values.
Print the entanglement entropy. What do you expect? Why do you get so many numbers, and not just one? Tipp: read the code ;-)
Extract the half-chain entanglement entropy, i.e., the entropy when cutting the chain into two equal-length halves.
[ ]:
Exercise
Write a function to initialize an MPS representing a product of singlets on neighboring sites. Check the expectation values of \(\sigma^z\), entropies and the correlation functions \(\sigma^z_0 \sigma^z_1\) and \(\sigma^z_0 \sigma^z_2\) for correctness.
[ ]:
toy codes: b_model.py
TFIModel class representing the transverse field Ising model\[H = - J \sum_{i} \sigma^z_i \sigma^z_{i+1} - g \sum_{i} \sigma^x_i\]
that provides both in the form of bond-terms H_bonds (as required for TEBD) and in the form of an MPO H_mpo. You can use H_bonds with SimpleMPS.bond_expectation_values to evalue the energy.
[ ]:
from tenpy_toycodes.b_model import TFIModel
model = TFIModel(L=10, J=1., g=1.2, bc='finite')
print("<H_bonds> = ", psi_FM.bond_expectation_value(model.H_bonds))
print("energy:", np.sum(psi_FM.bond_expectation_value(model.H_bonds)))
# (make sure the model and state have the same length and boundary conditions!)
Exercise
Check the energies for the Neel state and make sure it matches what you expect.
[ ]:
toy codes: c_tebd.py
The file c_tebd.py implements the TEBD algorithm.
You can run TEBD with imaginary time evolutinon to find the ground state like this:
[ ]:
from tenpy_toycodes.c_tebd import calc_U_bonds, run_TEBD
[ ]:
# assuming you defined a `model`
chi_max = 15
psi = init_FM_MPS(model.L, model.d, model.bc)
for dt in [0.1, 0.01, 0.001, 1.e-4, 1.e-5]:
U_bonds = calc_U_bonds(model.H_bonds, dt)
run_TEBD(psi, U_bonds, N_steps=100, chi_max=chi_max, eps=1.e-10)
E = np.sum(psi.bond_expectation_value(model.H_bonds))
print("dt = {dt:.5f}: E = {E:.13f}".format(dt=dt, E=E))
# the `run_TEBD` modified `psi`, so we can now calculate the expectation values from it
print("final bond dimensions: ", psi.get_chi())
mag_x = np.sum(psi.site_expectation_value(model.sigmax))
mag_z = np.sum(psi.site_expectation_value(model.sigmaz))
print("magnetization in X = {mag_x:.5f}".format(mag_x=mag_x))
print("magnetization in Z = {mag_z:.5f}".format(mag_z=mag_z))
Exercise
Read the code of c_tebd.py and try to understand the general structure of how it works.
[ ]:
[ ]:
[ ]:
toy codes: d_dmrg.py
The file d_dmrg.py implements the DMRG algorithm. It can be called like this:
[ ]:
from tenpy_toycodes.d_dmrg import SimpleDMRGEngine, SimpleHeff2
[ ]:
chi_max = 15
psi = init_FM_MPS(model.L, model.d, model.bc)
eng = SimpleDMRGEngine(psi, model, chi_max=chi_max, eps=1.e-10)
for i in range(10):
E_dmrg = eng.sweep()
E = np.sum(psi.bond_expectation_value(model.H_bonds))
print("sweep {i:2d}: E = {E:.13f}".format(i=i + 1, E=E))
print("final bond dimensions: ", psi.get_chi())
mag_x = np.mean(psi.site_expectation_value(model.sigmax))
mag_z = np.mean(psi.site_expectation_value(model.sigmaz))
print("magnetization in X = {mag_x:.5f}".format(mag_x=mag_x))
print("magnetization in Z = {mag_z:.5f}".format(mag_z=mag_z))
Exercise
Read the code of d_dmrg.py and try to understand the general structure of how it works.
Exercise
Compare running the code cells for DMRG and imaginary time evolution TEBD abovce for various parameters of L, J, g, and bc. Which one is faster? Do they always agree?
Note: The transverse field Ising model has a quantum phase transition at \(g/J = 1.\), from a ferro-magnetic phase (g < J) to a paramagnetic phase (g > J). We will map out the full phase diagram below.
[ ]:
[ ]:
Infinite DMRG
The SimpleDMRG code also allows to run infinite DMRG, simply by replacing the bc='finite' for both the model and the MPS. Look at the implementation of d_dmrg.py to see where the differences are.
The L parameter now just indices the number of tensors insite the unit cell of the infinite MPS. It has to be at least 2, since we optimize 2 tensors at once in our DMRG code. Further, note that we now use the mean to calculate densities of observables instead of extensive quantities (like L).
[ ]:
model = TFIModel(L=2, J=1., g=0.2, bc='infinite')
[ ]:
chi_max = 10
psi = init_FM_MPS(model.L, model.d, model.bc)
eng = SimpleDMRGEngine(psi, model, chi_max=chi_max, eps=1.e-7)
for i in range(10):
E_dmrg = eng.sweep()
E = np.mean(psi.bond_expectation_value(model.H_bonds))
print("sweep {i:2d}: E/L = {E:.13f}".format(i=i + 1, E=E))
print("final bond dimensions: ", psi.get_chi())
mag_x = np.mean(psi.site_expectation_value(model.sigmax))
mag_z = np.mean(psi.site_expectation_value(model.sigmaz))
print("magnetization in X = {mag_x:.5f}".format(mag_x=mag_x))
print("magnetization in Z = {mag_z:.5f}".format(mag_z=mag_z))
Exercise
Extend the following function run_DMRG to initialize an MPS compatible with the model, run DMRG on that MPS and return the corresponding state.
[ ]:
def run_DMRG(model, chi_max=50):
print(f"runnning DMRG for L={model.L:d}, g={model.g:.2f}, bc={model.bc}, chi_max={chi_max:d}")
raise NotImplementedError("TODO: this is an Exercise!")
return psi
Exercise
Use that function to check that the average energy density and expectation values for infinite DMRG are independent of the unit cell length L.
[ ]:
Exercise
Use the run_DMRG function above to plot expectation values as a function of g for fixed J=1 with the following code.
Modify it to also plot the average expectation values of sigmaz and sigmax and the sigmaz correlation function between spin 0 and 20.
[ ]:
#gs = [0.1, 0.5, 1.0, 1.5, 2.0]
gs = np.linspace(0., 2., 21)
entropies = []
for g in gs:
model = TFIModel(L=2, J=1., g=g, bc='infinite')
psi = run_DMRG(model, chi_max=50)
entropies.append(np.max(psi.entanglement_entropy()))
[ ]:
plt.plot(gs, entropies, marker='o', label='max entropy $S$')
# plot expecation values of sigmax and sigmaz as well
plt.xlabel('$g/J$')
plt.ylabel('expectation values')
plt.legend(loc='best')
Defining a new model: the XX chain
For the time evolution, we want to consider the XX chain in a staggered field, given by the Hamiltonian
for the usual Pauli matrices \(\sigma^x, \sigma^y, \sigma^z\).
A Jordan-Wigner transformation maps the XX Chain to free fermions, which we can diagonalize exactly with a few lines of python codes that are given in free_fermions_exact.py
[ ]:
from tenpy_toycodes.free_fermions_exact import XX_model_ground_state_energy
print("E_exact = ", XX_model_ground_state_energy(L=10, h_staggered=0.))
Exercise
The following code implements the model for the XX Chain hamiltonian, but the MPO lacks some terms. Fill them in!
Tip: In Python, (-1)**i represents \((-1)^i\).
Tip: For the Hamiltonian
a possible MPO matrix W looks like
Which parts do we need here?
Compare the energies of DMRG and (imaginary) TEBD to check that you got it correct.
Note that you have to use the Neel state as initial state to be in the right charge sector: the |up up up ... up > state already is an eigenstate of this Hamiltonian, so neither DMRG nor TEBD will modify it!
[ ]:
class XXChain:
"""Simple class generating the Hamiltonian of the
The Hamiltonian reads
.. math ::
H = - J \\sum_{i} \\sigma^x_i \\sigma^x_{i+1} - g \\sum_{i} \\sigma^z_i
"""
def __init__(self, L, hs, bc='finite'):
assert bc in ['finite', 'infinite']
self.L, self.d, self.bc = L, 2, bc
self.hs = hs
self.sigmax = np.array([[0., 1.], [1., 0.]]) # Pauli X
self.sigmay = np.array([[0., -1j], [1j, 0.]]) # Pauli Y
self.sigmaz = np.array([[1., 0.], [0., -1.]]) # Pauli Z
self.id = np.eye(2)
self.init_H_bonds()
self.init_H_mpo()
def init_H_bonds(self):
"""Initialize `H_bonds` hamiltonian."""
sx, sy, sz, id = self.sigmax, self.sigmay, self.sigmaz, self.id
d = self.d
nbonds = self.L - 1 if self.bc == 'finite' else self.L
H_list = []
for i in range(nbonds):
hL = hR = 0.5 * self.hs
if self.bc == 'finite':
if i == 0:
hL = self.hs
if i + 1 == self.L - 1:
hR = self.hs
H_bond = np.kron(sx, sx) + np.kron(sy, sy)
H_bond = H_bond - hL * (-1)**i * np.kron(sz, id) - hR * (-1)**(i+1) * np.kron(id, sz)
# H_bond has legs ``i, j, i*, j*``
H_list.append(np.reshape(H_bond, [d, d, d, d]))
self.H_bonds = H_list
# (note: not required for TEBD)
def init_H_mpo(self):
"""Initialize `H_mpo` Hamiltonian."""
w_list = []
for i in range(self.L):
w = np.zeros((4, 4, self.d, self.d), dtype=complex)
w[0, 0] = w[3, 3] = self.id
raise NotImplementedError("add further entries here")
w_list.append(w)
self.H_mpo = w_list
model = XXChain(9, 4., bc='finite')
[ ]:
print("E_exact = ", XX_model_ground_state_energy(model.L, model.hs))
[ ]:
chi_max = 100
psi = init_Neel_MPS(model.L, model.d, model.bc) # important: Neel
for dt in [0.1, 0.01, 0.001, 1.e-4, 1.e-5]:
U_bonds = calc_U_bonds(model.H_bonds, dt)
run_TEBD(psi, U_bonds, N_steps=100, chi_max=chi_max, eps=1.e-10)
E = np.sum(psi.bond_expectation_value(model.H_bonds))
print("dt = {dt:.5f}: E = {E:.13f}".format(dt=dt, E=E))
[ ]:
# DMRG again
chi_max = 100
psi = init_Neel_MPS(model.L, model.d, model.bc) # important: Neel!
eng = SimpleDMRGEngine(psi, model, chi_max=chi_max, eps=1.e-10)
for i in range(10):
eng.sweep()
E = np.sum(psi.bond_expectation_value(model.H_bonds))
#print(psi.get_chi())
print("sweep {i:2d}: E = {E:.13f}".format(i=i + 1, E=E))
[ ]:
Time evolution of the Neel state under the XX chain
We will now consider the Neel state as an initial state for (real-)time evolution, i.e., calculate
under the Hamiltonian of the XX chain (with staggered field).
Again, we can use the mapping to free fermions to calculate some quantities exactly. We will focus on the half chain entanglement entropy.
[ ]:
from tenpy_toycodes.free_fermions_exact import XX_model_time_evolved_entropies
dt = 0.1
tmax = 4.
L = 50
hs = 0.
model = XXChain(L, hs)
times_exact = np.arange(0., tmax, dt)
S_exact = XX_model_time_evolved_entropies(L=L, h_staggered=hs, time_list=times_exact)
Running TEBD works like this:
[ ]:
chi_max = 30
psi = init_Neel_MPS(L, 2, bc='finite')
U_bonds = calc_U_bonds(model.H_bonds, 1.j * dt) # here, imaginary dt = real time evolution
t = 0.
times_tebd = []
S_tebd = []
while t < tmax:
times_tebd.append(t)
S_tebd.append(psi.entanglement_entropy()[(L-1)//2])
print(f"t={t:.2f}")
t += dt
run_TEBD(psi, U_bonds, N_steps=1, chi_max=chi_max, eps=1.e-7)
For reference, one can also run TDVP with the e_tdvp.py as follows. If you know how TDVP works, try to take a look and spot the structure in the code. Otherwise, you can also consider it a black box for now.
The toycode implements the TDVP time evolution for finite MPS with the SimpleTDVPEngine. It implements both the true TDVP with one-site updates in sweep_one_site, and sweep_two_site that allows to grow the bond dimension.
[ ]:
chi_max = 30
t_switch = 1.
psi = init_Neel_MPS(L, 2, bc='finite')
eng = SimpleTDVPEngine(psi, model, chi_max=chi_max, eps=1.e-7)
t = 0.
times_tdvp = []
S_tdvp = []
while t < tmax:
times_tdvp.append(t)
S_tdvp.append(psi.entanglement_entropy()[(L-1)//2])
print(f"t={t:.2f}, chi_max = {max(psi.get_chi()):d}")
t += dt
if t < t_switch:
eng.sweep_two_site(dt)
else:
eng.sweep_one_site(dt)
[ ]:
[ ]:
plt.plot(times_exact, S_exact, label="exact")
plt.plot(times_tebd, S_tebd, label="TEBD")
plt.plot(times_tdvp, S_tdvp, label="TDVP")
plt.xlabel('time $t$')
plt.ylabel('half-chain entropy $S$')
plt.legend(loc='best')
Exercise
Add curves for smaller (or larger) chi to the plot above. (It might take a little while to run for large chi…).
Try switching earlier or later from two-site to one-site TDVP.
What happens if you switch on the staggered field?
NOTE: You can avoid having to rerun the time evolution if you make copies of the results S_tdvp and similar lists.
[ ]:
[ ]:
[ ]:
Advanced exercises - if you’re an expert and have time left ;-)
Obtain the ground state of the transverse field ising model at the critical point with DMRG. Try to plot the corrlation function as a function of
j-i. What form does it have? Is an MPS a good ansatz for that?Calling
SimpleMPS.correlation_functionin a loop overjfor fixediis inefficient for largej-i. Try to rewrite the correlation function to obtain the values for allj>iup to some cutoff in a single loop.Look at the light-cone after a local quench of the XX chain, similarly as is done in
e_tdvp.example_TDVP_tf_ising_lightcone. How does the speed of the light cone change with the staggered field?
[ ]:
[ ]:
[ ]:
[ ]: