Skip to content
Snippets Groups Projects
Commit 47819331 authored by Donjan Rodic's avatar Donjan Rodic
Browse files

Merge branch 'master' of gitlab.phys.ethz.ch:compphys_lectures/cqp_fs16

parents 6023b384 d0d39cc4
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 @ 47819331
%% 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 @ 47819331
%% 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 @ 47819331
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