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
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 thecorrelation_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
a possible MPO matrix W looks like
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
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')
[ ]:
[ ]: