# Solutions to MPS and model basics

In this notebook, we introduce the SimpleMPS class from tenpy_toycodes/a_mps.py and the model class from tenpy_toycodes/b_model.py.

[1]:
# standard imports and cosmetics

import numpy as np

import matplotlib.pyplot as plt

np.set_printoptions(precision=5, suppress=True, linewidth=100, threshold=50)
plt.rcParams['figure.dpi'] = 150

## SimpleMPS class from tenpy_toycodes/a_mps.py

The file tenpy_toycodes/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 to initialize the state $$\lvert \uparrow \uparrow \cdots \uparrow \uparrow \rangle$$, and - init_Neel_MPS to initialize the Neel state $$\lvert \uparrow \downarrow \cdots \uparrow \downarrow \rangle$$

[2]:
import tenpy_toycodes.a_mps
from tenpy_toycodes.a_mps import SimpleMPS, init_FM_MPS, init_Neel_MPS
[3]:
psi_FM = init_FM_MPS(L=10)
print(psi_FM)
<tenpy_toycodes.a_mps.SimpleMPS object at 0x7fd23a0df6d0>
[4]:
Z = np.diag([1., -1.])

print(psi_FM.site_expectation_value(Z))
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

### Exercise: expectation values and entropy

• Initialize a Neel state MPS. Print the expectation values of Z

[5]:
psi_Neel = init_Neel_MPS(L=10)
print(psi_Neel.site_expectation_value(Z))
[ 1. -1.  1. -1.  1. -1.  1. -1.  1. -1.]
• Print the entanglement entropy. What do you expect? Why do you get so many numbers, and not just one? Tip: read the code ;-)

[6]:
print(psi_Neel.entanglement_entropy())
# one value for cutting at each bond!
[-0. -0. -0. -0. -0. -0. -0. -0. -0.]
• Extract the half-chain entanglement entropy, i.e., the entropy when cutting the chain into two equal-length halves.

[7]:
print("half-chain entropy: ", psi_Neel.entanglement_entropy()[(psi_Neel.L - 1)//2])
# note: (L-1)//2, since we get L-1 values returned by entanglement_entropy()
# a product state has no entanglement!
half-chain entropy:  -0.0
• Read the code of a_mps.py to find out how to get the correlation $$\langle \psi| Z_1 Z_6 |\psi \rangle$$. Try it out!

[8]:
print("<psi|Z_1 Z_6|psi> = ", psi_Neel.correlation_function(Z, 1, Z, 6))
<psi|Z_1 Z_6|psi> =  -1.0
[ ]:

### Exercise: init_PM_MPS()

Write a function init_PM_MPS to initialize the state $$\lvert \rightarrow \rightarrow \cdots \rightarrow \rightarrow \rangle$$, where $$\lvert \rightarrow \rangle = \frac{1}{\sqrt{2}} \big( \lvert\uparrow \rangle + \lvert\downarrow\rangle \big)$$ is the spin-1/2 state pointing in plus x direction.

Tip: the code should be similar to init_FM_MPS.

[9]:
def init_PM_MPS(L, bc='finite'):
"""Return a paramagnetic MPS (= product state with all spins pointing in +x direction)"""
d = 2
B = np.zeros([1, d, 1], dtype=float)
B[0, 0, 0] = 1./np.sqrt(2)
B[0, 1, 0] = 1./np.sqrt(2)
S = np.ones([1], dtype=float)
Bs = [B.copy() for i in range(L)]
Ss = [S.copy() for i in range(L)]
return SimpleMPS(Bs, Ss, bc=bc)
[10]:
psi_PM = init_PM_MPS(L=10)

## Model class from tenpy_toycodes/b_model.py

The file tenpy_toycodes/b_model.py defines a TFIModel class representing the transverse field Ising model

$H = - J \sum_{i} Z_i Z_{i+1} - g \sum_{i} X_i$

It provides the Hamiltonian both in the form of bond-terms H_bonds (as required for TEBD) and in the form of an MPO H_mpo (as required for DMRG). You can use H_bonds with SimpleMPS.bond_expectation_values to evalue the energy:

[11]:
from tenpy_toycodes.b_model import TFIModel

L = 10
J = 1.
g = 1.2
model = TFIModel(L=L, J=J, g=g, 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!)

print("should be", (L-1)* (-J) * (1. * 1.) + L * (-g) * (0.) )
<H_bonds> =  [-1. -1. -1. -1. -1. -1. -1. -1. -1.]
energy: -9.0
should be -9.0

### Exercise

• Find the code where the MPO W for this model is defined to be $$W = \begin{pmatrix} \mathbb{1} & \sigma^Z & -g \sigma^X \\ & & -J \sigma^Z \\ & & \mathbb{1} \end{pmatrix}$$

[12]:
# the code is in the method init_H_mpo() of the TFIModel class
• Check the energies for the other initial states and make sure it matches what you expect.

[13]:
psi_Neel = init_Neel_MPS(L=L)
print("energy:", np.sum(psi_Neel.bond_expectation_value(model.H_bonds)))
print("should be", (L-1)* (-J) * (1. * -1.) + L * (-g) * (0.) )
energy: 9.0
should be 9.0
[14]:
psi_PM = init_PM_MPS(L=L)
print("energy:", np.sum(psi_PM.bond_expectation_value(model.H_bonds)))
print("should be", (L-1)* (-J) * (0. * 0.) + L * (-g) * (1.) )
energy: -11.999999999999993
should be -12.0
[ ]:

### Exercises (optional, if time left - for the experts, and those who want to become them)

• Write an optimized function correlation_function_all_j(psi, op_i, i, op_j, max_j) that returns the same values as the following, naive snippet:

results = []
for j in range(i+1, max_j):
results.append(psi.correlation_function(op_i, i, op_j, j)
return results

This snippet is $$\mathcal{O}(L^2)$$: for each j it calls the correlation_function, which internally also has an $$\mathcal{O}(L)$$ loop for the contractions. You can get this down to a single $$\mathcal{O}(L)$$ loop, if you identify and reuse the parts that are the same in the diagrams for the contractions.

[ ]:

• For benchmarks, it’s useful 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 tenpy_toycodes/free_fermions_exact.py.

[15]:
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.))
E_exact =  -12.053348366664542

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 different states to check that you got it correct. To really test it well, you should also check some states with non-trivial entanglement, e.g. ground states as obtained by DMRG (obtained as discussed in the next exercise notebook), for which you can directly compare to the XX_model_ground_state_energy.

[16]:
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