Skip to content
Snippets Groups Projects
Commit d0d39cc4 authored by Georg Wolfgang Winkler's avatar Georg Wolfgang Winkler
Browse files

added ex13

parent ad80a105
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
exercises/ex13_skeleton/heisenberg_dmrg_skeleton.ipynb 0 → 100644
+ 222
− 0
View file @ d0d39cc4
%% Cell type:code id: tags:
``` python
import numpy as np
from numpy import transpose as tr, conjugate as co
from scipy.linalg import expm, svd
from scipy.sparse.linalg import eigsh, LinearOperator
import math
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = 16, 9
```
%% Cell type:markdown id: tags:
### Useful functions
%% Cell type:code id: tags:
``` python
def dot(A,B):
""" Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions """
s1 = A.shape
s2 = B.shape
return np.dot(A,B).reshape((s1[0],s2[1]))
def left_canonize(M1,M2,return_S = False):
""" Left normalizes M1 into A matrix, M2 looses its canonization"""
s, da, db = M1.shape
U, S, Vh = svd(M1.reshape((s*da,db)))
U = U.reshape((s,da,U.shape[1]))[:,:,:db] # is now A matrix
M2 = np.tensordot(dot(np.diag(S[:db]),Vh[:db,:]),M2,axes=((1),(1))) # looses its canonization
if return_S:
return U, M2, S
else:
return U, M2
def right_canonize(M1,M2,return_S = False):
""" Right normalizes M2 into B matrix, M1 looses its canonization """
s, da, db = M2.shape
U, S, Vh = svd(M2.transpose((1,0,2)).reshape((da,s*db)))
Vh = Vh.reshape((Vh.shape[0],s,db)).transpose((1,0,2))[:,:da,:] # is now B matrix
M1 = np.tensordot(M1,dot(U[:,:da],np.diag(S[:da])),axes=((2),(0))) # looses its canonization
if return_S:
return M1, Vh, S
else:
return M1, Vh
def add_site_to_R(Rexp,M,W):
## to be completed
return Rexp
def add_site_to_L(Lexp,M,W):
## to be completed
return Lexp
def Hpsi(L,W,R,M):
"""
Calculates the "matrix-vector product" at arbitrary site from
L ... L-expression, belonging to the "matrix"
W ... W-array at site l, belonging to the "matrix"
R ... R-expression, belonging to the "matrix"
M ... M-array, the "vector"
"""
## to be completed
```
%% Cell type:markdown id: tags:
### Definitions and construction of MPO
%% Cell type:code id: tags:
``` python
J = 1.
L = 40 # Length of chain
s = 2 # Local dimension of Hilbert space
chi = 60 # The maximum matrix dimension, from which on the matrices are truncated
nmax = 10000 # Maximum number of iterations
## Construct the MPO for the Heisenberg model below
# ...
```
%% Cell type:markdown id: tags:
### Random guess as groundstate
%% Cell type:code id: tags:
``` python
## random starting configuration
M = [np.random.standard_normal((s,min(chi,s**i,s**(L-i)),min(chi,s**(i+1),s**(L-1-i)))) for i in range(L)]
## need a right canonized MPS
for i in range(L-1,0,-1):
M[i-1], M[i] = right_canonize(M[i-1],M[i])
## normalization, throw away the left part
_, M[0] = right_canonize(np.ones((1,1,1)),M[0])
```
%% Cell type:markdown id: tags:
### DMRG sweeps
%% Cell type:code id: tags:
``` python
## get initial R expression
Rexp = [np.ones((1,1,1))]
for i in range(1,L):
Rexp.append(add_site_to_R(Rexp[i-1],M[-i],W[-i]))
enold = 0
for n in range(nmax):
## right sweep
Lexp = [np.ones((1,1,1))]
for i in range(L-1):
print('.', end="")
# ----define a linear operator with square dimension s*D^2 via providing a matrix vector product---
dD2 = s*Lexp[i].shape[0]*Rexp[L-1-i].shape[0]
Hop = LinearOperator((dD2,dD2),matvec= lambda M: Hpsi(Lexp[i],W[i],Rexp[L-1-i],M))
# ----Use an iterative eigenvalue solver, M[i] needs to be flattened into a vector---
en,V = eigsh(Hop,k=1,v0=M[i].flatten(),tol=1e-3 if n < 2 else 0)
# ----Reshape M[i] back into its 3D array shape----
M[i] = V.reshape((s,Lexp[i].shape[0],Rexp[L-1-i].shape[0]))
# ----Left normalize M[i]---
M[i], M[i+1] = left_canonize(M[i],M[i+1])
# ----add site to L-expression----
Lexp.append(add_site_to_L(Lexp[i],M[i],W[i]))
print('')
## left sweep
Rexp = [np.ones((1,1,1))]
for i in range(L-1,0,-1):
print('.', end="")
## to be completed
print('')
# === check convergence ===
print(n,'dE = ',abs(en-enold))
if abs(en-enold) < 1e-8:
print("Converged after "+str(n)+" sweeps!")
print("Ground-state energy: " + str(en))
break
enold = en
# === check convergence ===
```
%% Cell type:code id: tags:
``` python
```
exercises/ex13_skeleton/heisenberg_imaginary_time_evolution.ipynb 0 → 100644
+ 364
− 0
View file @ d0d39cc4
%% Cell type:code id: tags:
``` python
import numpy as np
from numpy import transpose as tr, conjugate as co
from scipy.linalg import expm, svd
import math
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = 16, 9
```
%% Cell type:markdown id: tags:
### Some useful functions
%% Cell type:code id: tags:
``` python
def dot(A,B):
""" Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions """
s1 = A.shape
s2 = B.shape
return np.dot(A,B).reshape((s1[0],s2[1]))
def truncated_svd(thetas,chi):
""" Does an svd on two-site matrix thetas, and truncates to chi or the last nonzero singular value """
U, S, Vh = svd(thetas,full_matrices = False)
# trunkieren
ind = np.where(np.isclose(np.cumsum(S[::-1])[::-1],0))[0]
if len(ind)>0:
chi_tr = min(chi,max(1,ind[0]))
else:
chi_tr = chi
S=S[:chi_tr]
w=1.-np.sum(S**2)
S/=math.sqrt(1-w)
U=U[:,:chi_tr]
Vh=Vh[:chi_tr,:]
return U, S, Vh
def left_canonize(thetas,s,chi, return_S = False):
""" Splits up a two-site matrix thetas into two one-site matrices such that the left one is canonized """
da, dg = thetas.shape[2], thetas.shape[3]
thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes
U, S, Vh = truncated_svd(thetas,chi)
db = len(S)
U = U.reshape((da,s,db)).transpose((1,0,2))
Vh = Vh.reshape((db,dg,s)).transpose((2,0,1))
for s2 in range(s):
Vh[s2] = dot(np.diag(S),Vh[s2])
if return_S:
return U, Vh, S
else:
return U, Vh
def right_canonize(thetas,s,chi, return_S = False):
""" Splits up a two-site matrix thetas into two one-site matrices such that the right one is canonized """
da, dg = thetas.shape[2], thetas.shape[3]
thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes
U, S, Vh = truncated_svd(thetas,chi)
db = len(S)
U = U.reshape((da,s,db)).transpose((1,0,2))
Vh = Vh.reshape((db,dg,s)).transpose((2,0,1))
for s1 in range(s):
U[s1] = dot(U[s1],np.diag(S))
if return_S:
return U, Vh, S
else:
return U, Vh
def apply_two_site_H(theta,H,s):
""" Applies a two-site operator on the two-site matrix theta """
thetas = np.zeros_like(theta)
for s1p in range(s):
for s2p in range(s):
for s1 in range(s):
for s2 in range(s):
thetas[s1p,s2p] += H[s1p*s+s2p,s1*s+s2] * theta[s1,s2]
return thetas
def combine_two_matrices(M1,M2,s):
""" Combines the two neighbouring one-site matrices M1 and M2 into a two-site matrix theta """
da, db, dg=M1.shape[1], M1.shape[2], M2.shape[2]
assert M1.shape[2] == M2.shape[1]
theta = np.zeros((s,s,da,dg),dtype=complex)
for s1 in range(s):
for s2 in range(s):
theta[s1,s2]=dot(M1[s1],M2[s2])
return theta
```
%% Cell type:markdown id: tags:
### Definition of system and initializiation
%% Cell type:code id: tags:
``` python
J = 1.
L = 20 # Length of chain
s = 2 # Local dimension of Hilbert space
dt = 0.05 # The imaginary timestep (should be 0.01 or smaller for good accuracy)
chi = 60 # The maximum matrix dimension, from which on the matrices are truncated
nmax = 10000 # Maximum number of iterations
nskip = 10 # Check energy convergence after nskip steps
# Two-site Hamiltonian
H = np.array([[J/4,0,0,0],
[0,-J/4,J/2,0],
[0,J/2,-J/4,0],
[0,0,0,J/4]])
# Imaginary time evolution operator
exp_beta_H=expm(-H*dt);
# antiferromagnetic starting configuration
M = []
for i in range(L):
ar = np.zeros((2,1,1),dtype=complex)
ar[0,0,0] = i%2
ar[1,0,0] = (i+1)%2
M.append(ar)
even = np.array(range(0,L-1,2))
odd = np.array(range(L+L%2-2,0,-2))
```
%% Cell type:markdown id: tags:
### Imaginary time evolution
%% Cell type:code id: tags:
``` python
enold = 0. # for checking energy convergence
for n in range(nmax):
# ++++ Trotter scheme: first even bonds, then odd bonds
# === apply exp_beta_H on all even bonds ===
for j in even:
# go from left to right
theta = combine_two_matrices(M[j],M[j+1],s)
thetas = apply_two_site_H(theta,exp_beta_H,s)
M[j], M[j+1] = left_canonize(thetas,s,chi)
# advance left canonization by a further step
if j < L-2:
theta = combine_two_matrices(M[j+1],M[j+2],s)
M[j+1], M[j+2] = left_canonize(theta,s,chi)
# === apply exp_beta_H on all even bonds ===
# === renormalize the state on the last site ===
da = M[L-1].shape[1]
theta = M[L-1].reshape((s*da,1))
U, _, _ = truncated_svd(theta,chi) # throw away norm
M[L-1] = U.reshape((s,da,1))
# === renormalize the state on the last site ===
if L%2 == 0:
# the right-most matrix has to be right canonized in this case
theta = combine_two_matrices(M[L-2],M[L-1],s)
M[L-2], M[L-1] = right_canonize(theta,s,chi)
# === apply exp_beta_H on all odd bonds ===
for j in odd:
# go from right to left
theta = combine_two_matrices(M[j-1],M[j],s)
thetas = apply_two_site_H(theta,exp_beta_H,s)
M[j-1], M[j] = right_canonize(thetas,s,chi)
# advance right canonization by a further step
if j>1:
theta = combine_two_matrices(M[j-2],M[j-1],s)
M[j-2], M[j-1] = right_canonize(theta,s,chi)
# === apply exp_beta_H on all odd bonds ===
# === renormalize the state on the first site ===
db = M[0].shape[2]
theta = M[0].reshape((1,s*db))
_, _, Vh = truncated_svd(theta,chi) # throw away norm
M[0] = Vh.reshape((s,1,db))
# === renormalize the state on the first site ===
# ++++ Calculate energy and check convergence ++++
if n%nskip == 0:
# Measure Energy via sum of two-site Operators H
en = 0.
for j in range(L-1):
theta = combine_two_matrices(M[j],M[j+1],s)
thetas = apply_two_site_H(theta,H,s)
for s1 in range(s):
for s2 in range(s):
en += np.trace(dot(thetas[s1,s2],co(tr(theta[s1,s2])))).real
M[j], M[j+1] = left_canonize(theta,s,chi) # Proceed with left canonization
# make right canonized MPS for next tDMRG step
for j in range(L-1,0,-1):
theta = combine_two_matrices(M[j-1],M[j],s)
M[j-1], M[j] = right_canonize(theta,s,chi)
print(n,'dE = ',abs(en-enold))
if abs(en-enold) < 1e-8:
print("Converged after "+str(n)+" steps!")
print("Ground-state energy: " + str(en))
break
enold = en
# === check convergence ===
if n == nmax:
print("Maximum iterations reached! Convergence criterium not satisfied!")
```
%% Cell type:markdown id: tags:
### Implement the Heisenberg MPO below
You can build it as a list of four-dimensional arrays with indexing $\hat W^{[i]} = W_{\sigma'_i\, \sigma_i\, b_{l-1}\, b_l}$. It can be represented as "matrices" with indices $(b_{l-1}, b_l)$ which have single-site operators as entries having $(\sigma'_l,\sigma_l)$ indices. For the Heiseberg model
\begin{equation}
H = \sum_{i=1}^{L-1} \left( \frac{J}{2} S_i^+ S_{i+1}^- + \frac{J}{2} S_i^- S_{i+1}^+ + J S_i^z S_{i+1}^z\right) ,
\end{equation}
the matrices are
\begin{equation}
\hat W^{[i]} = \begin{bmatrix}
\hat I & 0 & 0 & 0 & 0 \\
S^+_i & 0 & 0 & 0 & 0 \\
S^-_i & 0 & 0 & 0 & 0 \\
S^z_i & 0 & 0 & 0 & 0 \\
0 & (J/2) S^-_i & (J/2) S^+_i & J S^z_i & \hat I
\end{bmatrix},
\end{equation}
and for the first and last site
\begin{equation}
\hat W^{[1]} = \begin{bmatrix}
0 & (J/2)S^-_0 & (J/2) S^+_0 & J S^z_0 & \hat I
\end{bmatrix}, \ \ \hat W^{[L]} = \begin{bmatrix}
\hat I \\ S^+_L \\ S^-_L \\ S^z_L \\ 0
\end{bmatrix}.
\end{equation}
%% Cell type:code id: tags:
``` python
## Generate MPO
W = []
```
%% Cell type:markdown id: tags:
### Build iteratively the L- or R-expression and evaluate the energy of the MPS using above MPO
You can use a three-dimensional array for $R$ with the indexing $R_{a'_l\, b_l\, a_l}$.
%% Cell type:code id: tags:
``` python
def add_site_to_R(Rexp,M,W):
"""Rexp... 3-d array
M ... M-array of the site to be added
W ... W-array of the site to be added"""
## to be completed
return Rexp
def add_site_to_L(Lexp,M,W):
## to be completed
return Lexp
## Calculate energy via iteratively constructing the R-expression
Rexp = np.ones((1,1,1))
for i in range(1,L+1):
Rexp = add_site_to_R(Rexp,M[-i],W[-i])
print('Energy from R-expression: ',np.squeeze(Rexp)) # should be total energy
## Calculate energy via iteratively constructing the L-expression
# ...
```
%% Cell type:code id: tags:
``` python
```
exercises/exercise13.pdf 0 → 100644
+ 0
− 0
View file @ d0d39cc4
File added
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment