Commit 85d5440b authored by Valerio's avatar Valerio
Browse files

Solution ex 11 and text last exercise ;)

parent 0583f9b8
%% Cell type:code id: tags:
``` python
from tenpy.models.lattice import Lattice
from tenpy.models.model import CouplingMPOModel
from tenpy.networks.site import SpinSite
from tenpy.tools.params import get_parameter
from tenpy.algorithms import dmrg
from tenpy.networks.mps import MPS
import numpy as np
import matplotlib.pyplot as plt
import time
import random
```
%% Cell type:markdown id: tags:
# Model Definition
%% Cell type:code id: tags:
``` python
#Create a ComplingMPOModel to describe the system.
#The parameters are passed via model_params
class Heisenberg(CouplingMPOModel):
def __init__(self,model_params):
CouplingMPOModel.__init__(self,model_params)
def init_lattice(self, model_params):
#Here initialize the type of lattice considered
lattice = get_parameter(model_params, 'lattice', self.name, False)
return lattice
def init_terms(self, model_params):
D= get_parameter(model_params, 'D', 0., self.name, True)
#Try to get also the J and lambda paramter in the same way
#J= ...
#lam= ...
#Here implement the couplings of the chain
#D term
self.add_onsite(D, 0, 'Sz Sz')
#Heisenberg interaction
self.add_coupling(J, 0, 'Sx', 0, 'Sx', 1,)
#Implement the remaining couplings for the Heisenberg exchange
#Have a look at the documentation of the function add_coupling()
#
#
#Quadratic Heisenberg term, the one with lambda
#Implement this term yourself
```
%% Cell type:markdown id: tags:
# Phase diagram for $\lambda=0$ and varying $D$
%% Cell type:markdown id: tags:
## Define the model
%% Cell type:code id: tags:
``` python
L=2 #Chain length (2 for the iMPS)
J=1 #Anitferromagnetic interaction
lam=0 #Quadratic term
Dmin=-1
Dmax=2
points=30
d=np.linspace(Dmin,Dmax,points,endpoint=True) #Values considered for D
#Definition of the lattice model
#Define the local Hilbert space, have a look at the documentation of SpinSite
#site=SpinSite(...)
#Define the lattice by looking at the documentation of Lattice
#Inifinite MPS requires periodic boundaries for the lattice
#lat=Lattice(...) #We choose an infinite MPS
```
%% Cell type:markdown id: tags:
## Perform iDMRG
%% Cell type:code id: tags:
``` python
entEnt=[]
entSp=[]
stagMag=[]
for D in d:
#Define the paramters of the model
model_params={
'verbose':0,
'J':J,
'D':D,
'lambda':lam,
'lattice':lat
}
#Create the model
chain=Heisenberg(model_params)
#Define the paramters for the DMRG
dmrg_paramsGs = {
'trunc_params': {
'chi_max': 20 #Maximum bond dimension
},
'verbose': 1
}
#Initialize the MPS with an arbitrary product state
product_state=['up',0]
#Have a look at the documentation MPS.from_product_state()
#psiGs=MPS.from_product_state(...)
#Perform iDMRG
#Check documentation of dmrg.run
#info=dmrg.run(...)
#Check documentation for expectation_value(). Which operator do you want to measure for magnetization?
#mag=psiGs.expectation_value(...) #Measure magnetization
#Implement function to compute staggered magnetization from mag
#staggered=...
stagMag.append(staggered) #Compute staggered magnetization
#Check documentation of entanglement_spectrum() and compute it
#spectrum=...
entSp.append(spectrum[0])
#Check documentation of entanglement_entropy() and compute it
#entropy=...
entEnt.append(entropy[0])
```
%% Cell type:markdown id: tags:
## Study staggered magnetization
%% Cell type:code id: tags:
``` python
plt.figure(figsize=(14,5))
plt.plot(d,stagMag,'o',ls='-',c='k')
plt.xlabel('D');
plt.ylabel('Staggered magnetization');
plt.axvspan(-1, -0.33, alpha=0.2, color='red')
plt.axvspan(-0.33, 1, alpha=0.2, color='b')
plt.axvspan(1, 2, alpha=0.2, color='g')
plt.xlim(-1,2)
```
%% Cell type:markdown id: tags:
## Study of the entanglement entropy
%% Cell type:code id: tags:
``` python
plt.figure(figsize=(14,5))
plt.tight_layout()
plt.plot(d,entEnt,'o',ls='-',c='k')
plt.xlabel('D');
plt.ylabel('Entanglement Entropy');
plt.axvspan(-1, -0.33, alpha=0.2, color='red')
plt.axvspan(-0.33, 1, alpha=0.2, color='b')
plt.axvspan(1, 2, alpha=0.2, color='g')
plt.xlim(-1,2);
```
%% Cell type:markdown id: tags:
## Study of the entanglement spectrum
%% Cell type:code id: tags:
``` python
plt.figure(figsize=(14,5))
plt.tight_layout()
for i in range(len(d)):
for j in range(int(len(entSp[i])/2)):
if abs(entSp[i][2*j]-entSp[i][2*j+1])<1e-4:
plt.scatter(d[i]-0.01,entSp[i][2*j],c='k',s=50,marker='.')
plt.scatter(d[i]+0.01,entSp[i][2*j+1],c='k',s=50,marker='.')
plt.plot([d[i]-0.03,d[i]+0.03],[entSp[i][2*j],entSp[i][2*j]],c='k')
else:
plt.scatter(d[i],entSp[i][2*j],c='k',s=50,marker='.')
plt.scatter(d[i],entSp[i][2*j+1],c='k',s=50,marker='.')
plt.plot([d[i]-0.03,d[i]+0.03],[entSp[i][2*j],entSp[i][2*j]],c='k')
plt.plot([d[i]-0.03,d[i]+0.03],[entSp[i][2*j+1],entSp[i][2*j+1]],c='k')
plt.xlabel('D');
plt.ylabel('Entanglement Spectrum');
plt.axvspan(-1.1, -0.33, alpha=0.2, color='red')
plt.axvspan(-0.33, 1, alpha=0.2, color='b')
plt.axvspan(1, 2.1, alpha=0.2, color='g')
plt.xlim(-1.1,2.1);
```
%% Cell type:markdown id: tags:
# Phase diagram for $D=0$ and varying $\lambda$
%% Cell type:code id: tags:
``` python
#Do the same studies for this new situation.
```
%% Cell type:code id: tags:
``` python
%matplotlib inline
from scipy import integrate
from scipy.linalg import expm
from pylab import *
import numpy as np
from svd_robust import svd
from ed_ising import ising_gs
import scipy.sparse.linalg.eigen.arpack as arp
```
%% Cell type:code id: tags:
``` python
""" Conventions:
B[i,a,b] has axes (physical, left virtual, right virtual),
W[a,b,i,j] has axes (virtual left, virtual right, physical out, physical in)
S[i] are schmidt values between sites (i, i+1),
H_bond[i] is the bond hamiltonian between (i,i+1) with (only physical)
axes (out left, out right, in left, in right)"""
def init_fm_mps(L):
""" Return FM Ising MPS"""
d = 2
B = []
s = []
for i in range(L):
B.append(np.zeros([2,1,1])); B[-1][0,0,0]=1
s.append(np.ones([1]))
s.append(np.ones([1]))
return B,s
def init_ising_H_mpo(g,J,L):
""" Returns hamiltonian in MPO form"""
s0 = np.eye(2)
sx = np.array([[0.,1.],[1.,0.]])
sy = np.array([[0.,-1j],[1j,0.]])
sz = np.array([[1.,0.],[0.,-1.]])
d = 2
w_list = []
for i in range(L):
w = np.zeros((3,3,d,d),dtype=np.float)
w[0,:2] = [s0,sz]
w[0:,2] = [g*sx, -J*sz, s0]
w_list.append(np.real(w))
return w_list
def init_ising_H_bond(g,J,L):
""" Returns bond hamiltonian """
sx = np.array([[0.,1.],[1.,0.]])
sy = np.array([[0.,-1j],[1j,0.]])
sz = np.array([[1.,0.],[0.,-1.]])
d = 2
H_bond = []
for i in range(L-2):
H_bond.append(np.reshape(-J*np.kron(sz,sz) + g*np.kron(sx,np.eye(2)),(d,d,d,d)))
H_bond.append(np.reshape(-J*np.kron(sz,sz)+g*np.kron(sx,np.eye(2))+g*np.kron(np.eye(2),sx),(d,d,d,d)))
return H_bond
```
%% Cell type:code id: tags:
``` python
def bond_expectation(B,s,H_bond):
" Expectation value for a bond operator "
E=[]
L = len(B)
for i_bond in range(L-1):
BB = np.tensordot(B[i_bond],B[i_bond+1],axes=(2,1))
sBB = np.tensordot(np.diag(s[i_bond]),BB,axes=(1,1))
C = np.tensordot(sBB,H_bond[i_bond],axes=([1,2],[2,3]))
sBB=np.conj(sBB)
E.append(np.squeeze(np.tensordot(sBB,C,axes=([0,3,1,2],[0,1,2,3]))).item())
return E
def entanglement_entropy(s):
" Returns the half chain entanglement entropy "
S=[]
L = len(s)
for i_bond in range(L):
x=s[i_bond][s[i_bond]>10**(-20)]**2
S.append(-np.inner(np.log(x),x))
return(S)
```
%% Cell type:code id: tags:
``` python
class H_mixed(object):
" Class of the mixed DMRG hamiltonian "
def __init__(self,Lp,Rp,M1,M2,dtype=float):
self.Lp = Lp
self.Rp = Rp
self.M1 = M1
self.M2 = M2
self.d = M1.shape[3]
self.chi1 = Lp.shape[0]
self.chi2 = Rp.shape[0]
self.shape = np.array([self.d**2*self.chi1*self.chi2,self.d**2*self.chi1*self.chi2])
self.dtype = dtype
def matvec(self,x):
x=np.reshape(x,(self.d,self.chi1,self.d,self.chi2)) # i a j b
x=np.tensordot(self.Lp,x,axes=(0,1)) # ap m i j b
x=np.tensordot(x,self.M1,axes=([1,2],[0,2])) # ap j b mp ip
x=np.tensordot(x,self.M2,axes=([3,1],[0,2])) # ap b ip m jp
x=np.tensordot(x,self.Rp,axes=([1,3],[0,2])) # ap ip jp bp
x=np.transpose(x,(1,0,2,3))
x=np.reshape(x,((self.d*self.d)*(self.chi1*self.chi2)))
if(self.dtype==float):
return np.real(x)
else:
return(x)
def update(theta,B,s,chi,H_mpo,Lp_list,Rp_list,chi1,chi3,i1,i2,i3,direction):
"""Updates the environment and the MPS"""
X, Y, Z = svd(theta); Z = Z.T
d = B[0].shape[0]
# Truncate and get the new Environment and MPS
chi2 = np.min([np.sum(Y>10.**(-12)), chi])
X=np.reshape(X[:d*chi1,:chi2],(d,chi1,chi2))
Z=np.transpose(np.reshape(Z[:d*chi3,:chi2],(d,chi3,chi2)),(0,2,1))
if direction == 'R':
Lp = np.tensordot(Lp_list[-1], H_mpo[i1], axes=(2,0))
Lp = np.tensordot(Lp, X, axes=([0,3],[1,0]))
Lp = np.tensordot(Lp, np.conj(X), axes=([0,2],[1,0]))
Lp = np.transpose(Lp,(1,2,0))
Lp_list.append(Lp)
if direction == 'L':
Rp = np.tensordot(H_mpo[i2], Rp_list[-1], axes=(1,2))
Rp = np.tensordot(np.conj(Z),Rp, axes=([0,2],[2,4]))
Rp = np.tensordot(Z,Rp, axes=([0,2],[2,3]))
Rp_list.append(Rp)
s[i2]=Y[:chi2]/np.sqrt(sum(Y[:chi2]**2))
B[i1]=np.transpose(np.tensordot(np.diag(s[i1]**(-1)),X,axes=(1,1)),(1,0,2))
B[i1]=np.tensordot(B[i1], np.diag(s[i2]),axes=(2,1))
B[i2] = Z
def diag(B,s,H,i1,i2,chi1,chi3):
""" Diagonalizes the mixed hamiltonian """
# Get a guess for the ground state based on the old MPS
d = B[0].shape[0]
theta0 = np.tensordot(np.diag(s[i1]),np.tensordot(B[i1],B[i2],axes=(2,1)),axes=(1,1))
theta0 = np.reshape(np.transpose(theta0,(1,0,2,3)),((chi1*chi3)*(d**2)))
# Diagonalize the sparse matrix
e,v = arp.eigsh(H,k=1,which='SA',return_eigenvectors=True,v0=theta0)
v = v.squeeze()
return np.reshape(v,(d*chi1,d*chi3))
def sweep(B,s,chi,H_mpo,Rp_list = None):
""" Performs a Left-Right-Left sweep """
d = B[0].shape[0]
D = H_mpo[0].shape[0]
L = len(B)
dtype = B[0].dtype
if Rp_list == None:
Rp = np.zeros([1,1,D],dtype=dtype)
Rp[0,0,D-1] = 1.
Rp_list = [Rp]
for i in np.arange(L-1,0,-1):
Rp = np.tensordot(B[i], Rp, axes=(2,0))
Rp = np.tensordot(H_mpo[i], Rp, axes=([1,2],[3,0]))
Rp = np.tensordot(np.conj(B[i]), Rp, axes=([0,2],[1,3]))
Rp = np.transpose(Rp,(2,0,1))
Rp_list.append(Rp)
# Right Sweep
Lp = np.zeros([1,1,D],dtype=dtype)
Lp[0,0,0] = 1.
Lp_list = [Lp]
for i_bond in range(L-1):
i1=i_bond
i2=i_bond+1
i3=i_bond+2
chi1 = Lp_list[-1].shape[1]
chi3 = Rp_list[L-2-i1].shape[1]
H = H_mixed(Lp_list[-1],Rp_list[L-2-i1],H_mpo[i1],H_mpo[i2],dtype=float)
theta = diag(B,s,H,i1,i2,chi1,chi3)
update(theta,B,s,chi,H_mpo,Lp_list,Rp_list,chi1,chi3,i1,i2,i3,'R')
# Left Sweep
Rp = np.zeros([1,1,D],dtype=dtype)
Rp[0,0,D-1] = 1.
Rp_list = [Rp]
for i_bond in range(L-2,-1,-1):
i1=i_bond
i2=i_bond+1
i3=i_bond+2
chi1 = Lp_list[i1].shape[1]
chi3 = Rp_list[-1].shape[1]
H = H_mixed(Lp_list[i1],Rp_list[-1],H_mpo[i1],H_mpo[i2],dtype=float)
theta = diag(B,s,H,i1,i2,chi1,chi3)
update(theta,B,s,chi,H_mpo,Lp_list,Rp_list,chi1,chi3,i1,i2,i3,'L')
return Rp_list
```
%% Cell type:code id: tags:
``` python
######## Define the simulation parameter ######################
chi = 10
N = 10
######## Define the Spin operators ###########################
J = 1.
g = 1.
L = 8
```
%% Cell type:code id: tags:
``` python
B,s = init_fm_mps(L)
H_mpo = init_ising_H_mpo(g,J,L)
H_bond = init_ising_H_bond(g,J,L)
Rp_list = None
for step in range(1, N):
Rp_list = sweep(B,s,chi,H_mpo,Rp_list = Rp_list)
print(np.sum(bond_expectation(B,s,H_bond)),"(DMRG)")
print(ising_gs(g,J,L),"(ED)")
```
%% Cell type:code id: tags:
``` python
```
This diff is collapsed.
%% Cell type:code id: tags:
``` python
from tenpy.models.lattice import Lattice
from tenpy.models.model import CouplingMPOModel
from tenpy.networks.site import SpinSite
from tenpy.tools.params import get_parameter
from tenpy.algorithms import dmrg
from tenpy.networks.mps import MPS
import numpy as np
import matplotlib.pyplot as plt
import time
import random
```
%% Cell type:markdown id: tags:
# Model Definition
%% Cell type:code id: tags:
``` python
#Create a ComplingMPOModel to describe the system.
#The parameters are passed via model_params
class Heisenberg(CouplingMPOModel):
def __init__(self,model_params):
CouplingMPOModel.__init__(self,model_params)
def init_lattice(self, model_params):
#Here initialize the type of lattice considered
lattice = get_parameter(model_params, 'lattice', self.name, False)
return lattice
def init_terms(self, model_params):
D= get_parameter(model_params, 'D', 0., self.name, True)
#Try to get also the J and lambda paramter in the same way
#J= ...
#lam= ...
#Here implement the couplings of the chain
#D term
self.add_onsite(D, 0, 'Sz Sz')
#Heisenberg interaction
self.add_coupling(J, 0, 'Sx', 0, 'Sx', 1,)
#Implement the remaining couplings for the Heisenberg exchange
#Have a look at the documentation of the function add_coupling()
#
#
#Quadratic Heisenberg term, the one with lambda
#Implement this term yourself
```
%% Cell type:markdown id: tags:
# Phase diagram for $\lambda=0$ and varying $D$
%% Cell type:markdown id: tags:
## Define the model
%% Cell type:code id: tags:
``` python
L=2 #Chain length (2 for the iMPS)
J=1 #Anitferromagnetic interaction
lam=0 #Quadratic term
Dmin=-1
Dmax=2
points=30
d=np.linspace(Dmin,Dmax,points,endpoint=True) #Values considered for D
#Definition of the lattice model
#Define the local Hilbert space, have a look at the documentation of SpinSite
#site=SpinSite(...)
#Define the lattice by looking at the documentation of Lattice
#Inifinite MPS requires periodic boundaries for the lattice
#lat=Lattice(...) #We choose an infinite MPS
```
%% Cell type:markdown id: tags:
## Perform iDMRG
%% Cell type:code id: tags:
``` python
entEnt=[]
entSp=[]
stagMag=[]
for D in d:
#Define the paramters of the model
model_params={
'verbose':0,
'J':J,