# 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

The file b_model.py defines a 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

$H = \sum_{i=0}^{N-2} (\sigma^x_i \sigma^x_{i+1} + \sigma^y_i \sigma^y_{i+1}) - h_s \sum_{i=0}^{N-1} (-1)^i \sigma^z_i = 2 \sum_{i=0}^{N-2} (\sigma^+_i \sigma^-_{i+1} + \sigma^+_i \sigma^-_{i+1}) - h_s \sum_{i=0}^{N-1} (-1)^i \sigma^z_i$

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

$H = \sum_{i=0}^{N-2} (\sigma^x_i \sigma^x_{i+1} + \sigma^y_i \sigma^y_{i+1} + \sigma^z_i \sigma^z_{i+1}) - h \sum_{i=0}^{N-1} \sigma^z_i,$

a possible MPO matrix W looks like

$\begin{split} W = \begin{pmatrix} 1 & \sigma^x & \sigma^y & \sigma^z & -h \sigma^z \\ 0 & 0 & 0 & 0 & \sigma^x \\ 0 & 0 & 0 & 0 & \sigma^y \\ 0 & 0 & 0 & 0 & \sigma^z \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} .\end{split}$

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
.. 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

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

$\vert \psi(t)\rangle =e^{-i H t} \vert\uparrow \downarrow \uparrow \downarrow \dots\rangle$

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_function in a loop over j for fixed i is inefficient for large j-i. Try to rewrite the correlation function to obtain the values for all j>i up 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?

[ ]:



[ ]:



[ ]:



[ ]: