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

added link to tDMRG example

parent 6a8162b6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/heisenberg_tDMRG/heisenberg_groundstate.ipynb
+ 7
− 0
View file @ 2fc60942
%% Cell type:markdown id: tags:
Algorithm after http://arxiv.org/pdf/1008.3477v2.pdf chapter 7.3.1
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
from numpy import transpose as tr, conjugate as co from numpy import transpose as tr, conjugate as co
from scipy.linalg import expm, svd from scipy.linalg import expm, svd
import math import math
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
%matplotlib inline %matplotlib inline
plt.rcParams['figure.figsize'] = 16, 9 plt.rcParams['figure.figsize'] = 16, 9
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Some useful functions ### Some useful functions
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def dot(A,B): def dot(A,B):
""" Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions """ """ Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions """
s1 = A.shape s1 = A.shape
s2 = B.shape s2 = B.shape
return np.dot(A,B).reshape((s1[0],s2[1])) return np.dot(A,B).reshape((s1[0],s2[1]))
def truncated_svd(thetas,chi): def truncated_svd(thetas,chi):
""" Does an svd on two-site matrix thetas, and truncates to chi or the last nonzero singular value """ """ 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) U, S, Vh = svd(thetas,full_matrices = False)
# trunkieren # trunkieren
ind = np.where(np.isclose(np.cumsum(S[::-1])[::-1],0))[0] ind = np.where(np.isclose(np.cumsum(S[::-1])[::-1],0))[0]
if len(ind)>0: if len(ind)>0:
chi_tr = min(chi,max(1,ind[0])) chi_tr = min(chi,max(1,ind[0]))
else: else:
chi_tr = chi chi_tr = chi
S=S[:chi_tr] S=S[:chi_tr]
w=1.-np.sum(S**2) w=1.-np.sum(S**2)
S/=math.sqrt(1-w) S/=math.sqrt(1-w)
U=U[:,:chi_tr] U=U[:,:chi_tr]
Vh=Vh[:chi_tr,:] Vh=Vh[:chi_tr,:]
return U, S, Vh return U, S, Vh
def left_canonize(thetas,s,chi, return_S = False): 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 """ """ 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] da, dg = thetas.shape[2], thetas.shape[3]
thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes
U, S, Vh = truncated_svd(thetas,chi) U, S, Vh = truncated_svd(thetas,chi)
db = len(S) db = len(S)
U = U.reshape((da,s,db)).transpose((1,0,2)) U = U.reshape((da,s,db)).transpose((1,0,2))
Vh = Vh.reshape((db,dg,s)).transpose((2,0,1)) Vh = Vh.reshape((db,dg,s)).transpose((2,0,1))
for s2 in range(s): for s2 in range(s):
Vh[s2] = dot(np.diag(S),Vh[s2]) Vh[s2] = dot(np.diag(S),Vh[s2])
if return_S: if return_S:
return U, Vh, S return U, Vh, S
else: else:
return U, Vh return U, Vh
def right_canonize(thetas,s,chi, return_S = False): 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 """ """ 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] da, dg = thetas.shape[2], thetas.shape[3]
thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes
U, S, Vh = truncated_svd(thetas,chi) U, S, Vh = truncated_svd(thetas,chi)
db = len(S) db = len(S)
U = U.reshape((da,s,db)).transpose((1,0,2)) U = U.reshape((da,s,db)).transpose((1,0,2))
Vh = Vh.reshape((db,dg,s)).transpose((2,0,1)) Vh = Vh.reshape((db,dg,s)).transpose((2,0,1))
for s1 in range(s): for s1 in range(s):
U[s1] = dot(U[s1],np.diag(S)) U[s1] = dot(U[s1],np.diag(S))
if return_S: if return_S:
return U, Vh, S return U, Vh, S
else: else:
return U, Vh return U, Vh
def apply_two_site_H(theta,H,s): def apply_two_site_H(theta,H,s):
""" Applies a two-site operator on the two-site matrix theta """ """ Applies a two-site operator on the two-site matrix theta """
thetas = np.zeros_like(theta) thetas = np.zeros_like(theta)
for s1p in range(s): for s1p in range(s):
for s2p in range(s): for s2p in range(s):
for s1 in range(s): for s1 in range(s):
for s2 in range(s): for s2 in range(s):
thetas[s1p,s2p] += H[s1p*s+s2p,s1*s+s2] * theta[s1,s2] thetas[s1p,s2p] += H[s1p*s+s2p,s1*s+s2] * theta[s1,s2]
return thetas return thetas
def combine_two_matrices(M1,M2,s): def combine_two_matrices(M1,M2,s):
""" Combines the two neighbouring one-site matrices M1 and M2 into a two-site matrix theta """ """ 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] da, db, dg=M1.shape[1], M1.shape[2], M2.shape[2]
assert M1.shape[2] == M2.shape[1] assert M1.shape[2] == M2.shape[1]
theta = np.zeros((s,s,da,dg),dtype=complex) theta = np.zeros((s,s,da,dg),dtype=complex)
for s1 in range(s): for s1 in range(s):
for s2 in range(s): for s2 in range(s):
theta[s1,s2]=dot(M1[s1],M2[s2]) theta[s1,s2]=dot(M1[s1],M2[s2])
return theta return theta
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Definition of system and initializiation ### Definition of system and initializiation
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
J = 1. J = 1.
L = 20 # Length of chain L = 20 # Length of chain
s = 2 # Local dimension of Hilbert space s = 2 # Local dimension of Hilbert space
dt = 0.05 # The imaginary timestep (should be 0.01 or smaller for good accuracy) 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 chi = 60 # The maximum matrix dimension, from which on the matrices are truncated
nmax = 10000 # Maximum number of iterations nmax = 10000 # Maximum number of iterations
nskip = 10 # Check energy convergence after nskip steps nskip = 10 # Check energy convergence after nskip steps
# Two-site Hamiltonian # Two-site Hamiltonian
H = np.array([[J/4,0,0,0], H = np.array([[J/4,0,0,0],
[0,-J/4,J/2,0], [0,-J/4,J/2,0],
[0,J/2,-J/4,0], [0,J/2,-J/4,0],
[0,0,0,J/4]]) [0,0,0,J/4]])
# Imaginary time evolution operator # Imaginary time evolution operator
exp_beta_H=expm(-H*dt); exp_beta_H=expm(-H*dt);
# antiferromagnetic starting configuration # antiferromagnetic starting configuration
M = [] M = []
for i in range(L): for i in range(L):
ar = np.zeros((2,1,1),dtype=complex) ar = np.zeros((2,1,1),dtype=complex)
ar[0,0,0] = i%2 ar[0,0,0] = i%2
ar[1,0,0] = (i+1)%2 ar[1,0,0] = (i+1)%2
M.append(ar) M.append(ar)
even = np.array(range(0,L-1,2)) even = np.array(range(0,L-1,2))
odd = np.array(range(L+L%2-2,0,-2)) odd = np.array(range(L+L%2-2,0,-2))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Imaginary time evolution ### Imaginary time evolution
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
enold = 0. # for checking energy convergence enold = 0. # for checking energy convergence
for n in range(nmax): for n in range(nmax):
# ++++ Trotter scheme: first even bonds, than odd bonds # ++++ Trotter scheme: first even bonds, than odd bonds
# === apply exp_beta_H on all even bonds === # === apply exp_beta_H on all even bonds ===
for j in even: for j in even:
# go from left to right # go from left to right
theta = combine_two_matrices(M[j],M[j+1],s) theta = combine_two_matrices(M[j],M[j+1],s)
thetas = apply_two_site_H(theta,exp_beta_H,s) thetas = apply_two_site_H(theta,exp_beta_H,s)
M[j], M[j+1] = left_canonize(thetas,s,chi) M[j], M[j+1] = left_canonize(thetas,s,chi)
# advance left canonization by a further step # advance left canonization by a further step
if j < L-2: if j < L-2:
theta = combine_two_matrices(M[j+1],M[j+2],s) theta = combine_two_matrices(M[j+1],M[j+2],s)
M[j+1], M[j+2] = left_canonize(theta,s,chi) M[j+1], M[j+2] = left_canonize(theta,s,chi)
# === apply exp_beta_H on all even bonds === # === apply exp_beta_H on all even bonds ===
# === renormalize the state on the last site === # === renormalize the state on the last site ===
da = M[L-1].shape[1] da = M[L-1].shape[1]
theta = M[L-1].reshape((s*da,1)) theta = M[L-1].reshape((s*da,1))
U, _, _ = truncated_svd(theta,chi) # throw away norm U, _, _ = truncated_svd(theta,chi) # throw away norm
M[L-1] = U.reshape((s,da,1)) M[L-1] = U.reshape((s,da,1))
# === renormalize the state on the last site === # === renormalize the state on the last site ===
if L%2 == 0: if L%2 == 0:
# the right-most matrix has to be right canonized in this case # the right-most matrix has to be right canonized in this case
theta = combine_two_matrices(M[L-2],M[L-1],s) theta = combine_two_matrices(M[L-2],M[L-1],s)
M[L-2], M[L-1] = right_canonize(theta,s,chi) M[L-2], M[L-1] = right_canonize(theta,s,chi)
# === apply exp_beta_H on all odd bonds === # === apply exp_beta_H on all odd bonds ===
for j in odd: for j in odd:
# go from right to left # go from right to left
theta = combine_two_matrices(M[j-1],M[j],s) theta = combine_two_matrices(M[j-1],M[j],s)
thetas = apply_two_site_H(theta,exp_beta_H,s) thetas = apply_two_site_H(theta,exp_beta_H,s)
M[j-1], M[j] = right_canonize(thetas,s,chi) M[j-1], M[j] = right_canonize(thetas,s,chi)
# advance right canonization by a further step # advance right canonization by a further step
if j>1: if j>1:
theta = combine_two_matrices(M[j-2],M[j-1],s) theta = combine_two_matrices(M[j-2],M[j-1],s)
M[j-2], M[j-1] = right_canonize(theta,s,chi) M[j-2], M[j-1] = right_canonize(theta,s,chi)
# === apply exp_beta_H on all odd bonds === # === apply exp_beta_H on all odd bonds ===
# === renormalize the state on the first site === # === renormalize the state on the first site ===
db = M[0].shape[2] db = M[0].shape[2]
theta = M[0].reshape((1,s*db)) theta = M[0].reshape((1,s*db))
_, _, Vh = truncated_svd(theta,chi) # throw away norm _, _, Vh = truncated_svd(theta,chi) # throw away norm
M[0] = Vh.reshape((s,1,db)) M[0] = Vh.reshape((s,1,db))
# === renormalize the state on the first site === # === renormalize the state on the first site ===
# ++++ Calculate energy and check convergence ++++ # ++++ Calculate energy and check convergence ++++
if n%nskip == 0: if n%nskip == 0:
# Measure Energy via sum of two-site Operators H # Measure Energy via sum of two-site Operators H
en = 0. en = 0.
for j in range(L-1): for j in range(L-1):
theta = combine_two_matrices(M[j],M[j+1],s) theta = combine_two_matrices(M[j],M[j+1],s)
thetas = apply_two_site_H(theta,H,s) thetas = apply_two_site_H(theta,H,s)
for s1 in range(s): for s1 in range(s):
for s2 in range(s): for s2 in range(s):
en += np.trace(dot(thetas[s1,s2],co(tr(theta[s1,s2])))).real 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 M[j], M[j+1] = left_canonize(theta,s,chi) # Proceed with left canonization
# make right canonized MPS for next tDMRG step # make right canonized MPS for next tDMRG step
for j in range(L-1,0,-1): for j in range(L-1,0,-1):
theta = combine_two_matrices(M[j-1],M[j],s) theta = combine_two_matrices(M[j-1],M[j],s)
M[j-1], M[j] = right_canonize(theta,s,chi) M[j-1], M[j] = right_canonize(theta,s,chi)
print(n,'dE = ',abs(en-enold)) print(n,'dE = ',abs(en-enold))
if abs(en-enold) < 1e-8: if abs(en-enold) < 1e-8:
print("Converged after "+str(n)+" steps!") print("Converged after "+str(n)+" steps!")
print("Ground-state energy: " + str(en)) print("Ground-state energy: " + str(en))
break break
enold = en enold = en
# === check convergence === # === check convergence ===
if n == nmax: if n == nmax:
print("Maximum iterations reached! Convergence criterium not satisfied!") print("Maximum iterations reached! Convergence criterium not satisfied!")
``` ```
%% Output %% Output
0 dE = 5.2022859299 0 dE = 5.2022859299
10 dE = 2.41237932615 10 dE = 2.41237932615
20 dE = 0.616048162871 20 dE = 0.616048162871
30 dE = 0.19945192398 30 dE = 0.19945192398
40 dE = 0.0856110936378 40 dE = 0.0856110936378
50 dE = 0.0446638884127 50 dE = 0.0446638884127
60 dE = 0.0266683370958 60 dE = 0.0266683370958
70 dE = 0.0176509913134 70 dE = 0.0176509913134
80 dE = 0.0127047223848 80 dE = 0.0127047223848
90 dE = 0.00977408061919 90 dE = 0.00977408061919
100 dE = 0.00789488605981 100 dE = 0.00789488605981
110 dE = 0.00658572945765 110 dE = 0.00658572945765
120 dE = 0.00559930606822 120 dE = 0.00559930606822
130 dE = 0.00480731763688 130 dE = 0.00480731763688
140 dE = 0.00414298945949 140 dE = 0.00414298945949
150 dE = 0.0035712463933 150 dE = 0.0035712463933
160 dE = 0.00307305023898 160 dE = 0.00307305023898
170 dE = 0.00263722905494 170 dE = 0.00263722905494
180 dE = 0.00225631290084 180 dE = 0.00225631290084
190 dE = 0.00192449985237 190 dE = 0.00192449985237
200 dE = 0.00163673118822 200 dE = 0.00163673118822
210 dE = 0.00138832203542 210 dE = 0.00138832203542
220 dE = 0.0011748521653 220 dE = 0.0011748521653
230 dE = 0.000992164333084 230 dE = 0.000992164333084
240 dE = 0.000836395474847 240 dE = 0.000836395474847
250 dE = 0.000704007231025 250 dE = 0.000704007231025
260 dE = 0.000591803752473 260 dE = 0.000591803752473
270 dE = 0.000496934572579 270 dE = 0.000496934572579
280 dE = 0.000416884722746 280 dE = 0.000416884722746
290 dE = 0.000349455574773 290 dE = 0.000349455574773
300 dE = 0.000292739848431 300 dE = 0.000292739848431
310 dE = 0.000245093656 310 dE = 0.000245093656
320 dE = 0.000205107728403 320 dE = 0.000205107728403
330 dE = 0.000171579325436 330 dE = 0.000171579325436
340 dE = 0.000143485782885 340 dE = 0.000143485782885
350 dE = 0.000119960236445 350 dE = 0.000119960236445
360 dE = 0.000100269773123 360 dE = 0.000100269773123
370 dE = 8.37960556659e-05 370 dE = 8.37960556659e-05
380 dE = 7.00183361975e-05 380 dE = 7.00183361975e-05
390 dE = 5.84987014349e-05 390 dE = 5.84987014349e-05
400 dE = 4.88693502501e-05 400 dE = 4.88693502501e-05
410 dE = 4.08216810968e-05 410 dE = 4.08216810968e-05
420 dE = 3.40969704631e-05 420 dE = 3.40969704631e-05
430 dE = 2.84784768727e-05 430 dE = 2.84784768727e-05
440 dE = 2.37847443216e-05 440 dE = 2.37847443216e-05
450 dE = 1.98639157318e-05 450 dE = 1.98639157318e-05
460 dE = 1.65889561448e-05 460 dE = 1.65889561448e-05
470 dE = 1.38536333161e-05 470 dE = 1.38536333161e-05
480 dE = 1.15691357934e-05 480 dE = 1.15691357934e-05
490 dE = 9.661232097e-06 490 dE = 9.661232097e-06
500 dE = 8.06788848529e-06 500 dE = 8.06788848529e-06
510 dE = 6.73727353018e-06 510 dE = 6.73727353018e-06
520 dE = 5.62608618715e-06 520 dE = 5.62608618715e-06
530 dE = 4.69815473814e-06 530 dE = 4.69815473814e-06
540 dE = 3.92326507992e-06 540 dE = 3.92326507992e-06
550 dE = 3.27618108287e-06 550 dE = 3.27618108287e-06
560 dE = 2.73582595511e-06 560 dE = 2.73582595511e-06
570 dE = 2.28459734686e-06 570 dE = 2.28459734686e-06
580 dE = 1.90779520892e-06 580 dE = 1.90779520892e-06
590 dE = 1.59314354775e-06 590 dE = 1.59314354775e-06
600 dE = 1.33039094941e-06 600 dE = 1.33039094941e-06
610 dE = 1.11097674527e-06 610 dE = 1.11097674527e-06
620 dE = 9.27752324742e-07 620 dE = 9.27752324742e-07
630 dE = 7.747481785e-07 630 dE = 7.747481785e-07
640 dE = 6.46979531282e-07 640 dE = 6.46979531282e-07
650 dE = 5.40283849659e-07 650 dE = 5.40283849659e-07
660 dE = 4.51185265149e-07 660 dE = 4.51185265149e-07
670 dE = 3.76781240874e-07 670 dE = 3.76781240874e-07
680 dE = 3.1464805339e-07 680 dE = 3.1464805339e-07
690 dE = 2.62761810532e-07 690 dE = 2.62761810532e-07
700 dE = 2.19432413573e-07 700 dE = 2.19432413573e-07
710 dE = 1.83248596741e-07 710 dE = 1.83248596741e-07
720 dE = 1.53031841421e-07 720 dE = 1.53031841421e-07
730 dE = 1.27798031713e-07 730 dE = 1.27798031713e-07
740 dE = 1.06725389415e-07 740 dE = 1.06725389415e-07
750 dE = 8.91276528137e-08 750 dE = 8.91276528137e-08
760 dE = 7.44317638635e-08 760 dE = 7.44317638635e-08
770 dE = 6.21591595973e-08 770 dE = 6.21591595973e-08
780 dE = 5.19102449914e-08 780 dE = 5.19102449914e-08
790 dE = 4.33512816755e-08 790 dE = 4.33512816755e-08
800 dE = 3.62035965651e-08 800 dE = 3.62035965651e-08
810 dE = 3.02344727032e-08 810 dE = 3.02344727032e-08
820 dE = 2.52495695463e-08 820 dE = 2.52495695463e-08
830 dE = 2.10865902517e-08 830 dE = 2.10865902517e-08
840 dE = 1.76100023452e-08 840 dE = 1.76100023452e-08
850 dE = 1.47066288037e-08 850 dE = 1.47066288037e-08
860 dE = 1.22819674431e-08 860 dE = 1.22819674431e-08
870 dE = 1.02570609783e-08 870 dE = 1.02570609783e-08
880 dE = 8.56601900523e-09 880 dE = 8.56601900523e-09
Converged after 880 steps! Converged after 880 steps!
Ground-state energy: -8.68010558621 Ground-state energy: -8.68010558621
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Measure and plot the Entanglement Entropy between all sites ### Measure and plot the Entanglement Entropy between all sites
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
ent = np.zeros((L-1)) ent = np.zeros((L-1))
## warning measurment assumes right canonized MPS ## warning measurment assumes right canonized MPS
for j in range(L-1): for j in range(L-1):
theta = combine_two_matrices(M[j],M[j+1],s) theta = combine_two_matrices(M[j],M[j+1],s)
M[j], M[j+1], S = left_canonize(theta,s,chi, return_S = True) M[j], M[j+1], S = left_canonize(theta,s,chi, return_S = True)
ent[j] = -np.sum(S**2*np.log(S**2)) ent[j] = -np.sum(S**2*np.log(S**2))
plt.figure() plt.figure()
plt.plot(ent,'-o') plt.plot(ent,'-o')
plt.xlabel('site i') plt.xlabel('site i')
plt.ylabel(r'$S(\rho_i)$') plt.ylabel(r'$S(\rho_i)$')
plt.show() plt.show()
``` ```
%% Output %% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......
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