#-------------------------------# # # # Week5: 1D Hubbard Model # # (nicht klausurrelevant) # #-------------------------------# Wir schauen uns zunaechst die Sachen aus dem week5.pdf an. Wir wollen also diese Matrix H anlegen. Das H ist eine |B| x |B| - Matrix, wenn B die Basis bezeichnet, die (1) in week5.pdf angegeben ist. Wir muessen uns zunaechst ueber- legen, wie wir diese Basis in R generieren koennen und durch welches Objekt wir sie dann konkret darstellen koennen: Es sei also N die Anzahl der Teilchen und L die Anzahl der Gitterpunkte. Wir muessen L natuerliche Zahlen n1,...,nL finden, die jeweils bei 0 anfangen duerfen, so dass gilt: n1 + n2 + ... + nL = N (1) Wir addieren L Einsen auf der linken Seite von (1) und ein L auf der rechten Seite von (1) und erhalten: (n1+1) + (n2+1) + ... + (nL+1) = N + L (2) oder m1 + m2 + ... + mL = N + L (3) wobei die m1,...,mL jetzt natuerliche Zahlen sind, die bei 1 anfangen muessen, nicht mehr bei 0. Die natuerliche Zahl N + L auf der rechten Seite von (3) koennen wir uns vorstellen als eine Summe von N+L Einsen, N+L = 1 + 1 + 1 + ... + 1 + 1 (4) N+L - Stueck Die Summe auf der rechten Seite von (4) muessen wir zerlegen in L Summanden. Das koennen wir dadurch erreichen, dass wir sozusagen L-1 Trenn- striche auf die insgesamt N+L-1 Plus-Zeichen verteilen, die Anzahl der Einsen zwischen zwei solchen Trennstrichen definiert dann den Summan- den mj. Wenn man aus eine Gesamtheit von n Sachen k Sachen auswaehlen soll, das ist wie Ziehung der Lottozahlen, da gibt es n ueber k Moeglichkeiten. Also haben wir hier: Anzahl n1,...,nL die (1) erfuellen = Anzahl m1,...,mL die (3) erfuellen = ( N + L - 1 ueber L - 1 ) (5) wobei in der letzten Zeile von (5) also ein Binomialkoeffizient steht. Fuer Binomialkoeffizienten gibt es einige Funk- tionalitaeten in R, zunaechst mal der Binomial- koeffizient selber: n = 10 k = 2 choose(n,k) # = 45 = 10*9/(1*2) = ( 10 ueber 2 ) Daneben gibt es aber noch die combn()-Funktion, und das ist jetzt genau das, was wir brauchen: combn(n,k) t(combn(n,k)) Das tut also aus den Zahlen von 1 bis 10 genau 2 Stueck auswaehlen. Damit koennen wir jetzt die Basis B aus (1) von week5.pdf folgendermassen generieren: #-------------------------------# # set up many body basis # #-------------------------------# L = 4 N = 5 nB = choose(N+L-1,L-1) # = |B| nB # etwa 56 fuer L=4, N=5 B0 = t(combn(N+L-1,L-1)) # t = transponieren, dann nB Zeilen B0 B = matrix(0,nrow=nB,ncol=L) # soll werden Basis B = {(n1,...,nL)| sum nj = N, nj>=0} B for(j in 1:nB) { b = B0[j,] # das b sind jetzt 2 Trennstriche, die bei b[1] und b[2] bb = c( 0 , b , N+L ) # auf die Plus-Zeichen in (4) weiter oben gelegt werden m = diff( bb ) # das sind jetzt die mj aus (3) weiter oben n = m - rep(1,L) # und das dann also die nj aus (1) B[j,] = n } B # ist jetzt die gewuenschte Basis B mit dim(B) # B = {(n1,...,nL)| sum nj = N, nj>=0} #-------------------------------# # END set up many body basis # #-------------------------------# Wir haben die Basis B aus pdf-Gleichung (1) jetzt also durch die Matrix B codiert, jede Zeile von B definiert einen Basisvektor |n1,...,nL> = B[j,] und der Zeilenindex j laeuft dann also von 1 bis |B| = nB. Jetzt koennen wir die Matrix H upsetten. H ist ja eine nB x nB - Matrix, schauen wir uns mal ein bischen an, wie gross die Zahlen nB so sind: wir haben ja nB = choose( N + L - 1 , L - 1 ) und damit bekommt man etwa folgende Werte: L N nB nB^2 3 10 66 4'356 3 20 231 53'361 6 10 3003 9'018'009 6 12 6188 38'291'344 6 15 15504 240'374'016 6 20 53130 2'822'796'900 ..oder etwa auch L N nB nB^2 100 100 4.53*10^58 2.05*10^117 Insbesondere an den letzten Zahlen kann man leicht erkennen: wie effizient man das auch immer irgendwie mit duennbesetzten Matrizen upsetten mag, der Zugang ueber die Matrix- Darstellung des Hubbard-Modells stoesst sehr schnell an Grenzen. Was kann man anstatt dessen machen? -> very active, very hot topic of current research Fixieren wir fuer das folgende vielleicht die folgenden Parameterwerte: L N nB nB^2 6 10 3003 9'018'009 Also die Matrix H, die wir anlegen wollen, ist im wesentlichen ein 3000 X 3000 Matrix mit 9 Mio Elementen. Davon sind aber nur sehr wenige ungleich null, wir werden finden, dass nur 23024 Elemente ungleich null sind. Anstatt jetzt also jedes der 9 Mio Elemente zu speichern, die ja fast alle 0 sind, tun wir also nur diejenigen Elemente speichern, die wirklich ungleich null sind. Der "Hauptindex" von dem neuen Objekt, was wir dann also abspeichern wollen, ist dann nicht mehr ein Doppelindex [i,j] mit 1<=i,j<=3000, sondern der neue Hauptindex nz zaehlt oder indi- ziert dann nur die nonzero elements, 'nz' also fuer 'nonzero'. Die Matrix H ( = h = h0 + hint weiter unten) wird dann durch folgende Objekte codiert: H = vector[nz] Zeile = vector[nz] Spalte = vector[nz] und etwa das nz = achte Element von H, was ungleich null ist, hat Zeilen- und Spaltenindex [i,j] mit i = Zeile[nz] j = Spalte[nz] Ok, codieren wir das jetzt: #-------------------------------# # set up Hamiltonian # #-------------------------------# # first: choose L and N and # set up many body basis -> code above eps = 1 g = 0.25 u = g/N # allocate some storage place for the result: nzmax = 10^7 # non zero matrix elements h0 = rep( 0 , nzmax ) # kinetic energy hint = rep( 0 , nzmax ) # diagonal part rows = rep( 0 , nzmax ) # location in full nB x nB sparse matrix cols = rep( 0 , nzmax ) # location in full nB x nB sparse matrix id = rep( 0 , nzmax ) # soll die Einheitsmatrix werden, # aber in dieser sparse-representation t0 = Sys.time() # monitor calc time nz = 1 # counter of non-zero matrix elements for(i in 1:nB) { m = B[i,] # Basisvektor m for(j in 1:nB) { n = B[j,] # Basisvektor n , wir berechnen # < m | H | n > # diagonal part: if(i==j) { res = 0 for(k in 1:L) { res = res + n[k]*(n[k]-1) # das ist (5) aus pdf-file } h0[nz] = 0 hint[nz] = res id[nz] = 1 rows[nz] = i cols[nz] = j nz = nz + 1 } # off-diagonal part: d = m - n k = which( d != 0 ) # nuetzlicher Befehl.. if( length(k)==2 ) { k1 = k[1] k2 = k[2] if( d[k1]*d[k2] == -1 & k2-k1==1 ) # we have a non zero contribution { if( m[k1] == n[k1]+1 ) { res = sqrt(n[k1]+1)*sqrt(n[k2]) # das ist die erste Zeile von (4).pdf h0[nz] = res hint[nz] = 0 id[nz] = 0 rows[nz] = i cols[nz] = j nz = nz + 1 } if( m[k1] == n[k1]-1 ) { res = sqrt(n[k1])*sqrt(n[k2]+1) # das ist die zweite Zeile von (4).pdf h0[nz] = res hint[nz] = 0 id[nz] = 0 rows[nz] = i cols[nz] = j nz = nz + 1 } } # END we have a non zero contribution } # END if( length(k)==2 ) } # NEXT j, basis element n = B[j,] # nur zur Info: if(i%%100==0) { print( c(i,nB) ) print( Sys.time() - t0 ) } } # NEXT i, basis element m = B[i,] # remove unused elements: if( nz < nzmax ) { h0 = h0[-(nz:nzmax)] hint = hint[-(nz:nzmax)] id = id[-(nz:nzmax)] rows = rows[-(nz:nzmax)] cols = cols[-(nz:nzmax)] } h = eps*h0 + u*hint nz # number of nonzero elements # etwa 23024 for L=6, N=10 Sys.time() - t0 # Finished setting up h #-------------------------------# # END set up Hamiltonian # #-------------------------------# Wir haben jetzt alle Matrix-Elemente berechnet, aber sie sind als Vektor gespeichert. Ok, anschauen koennen wir sie ja schonmal: # info matrix elements: summary(h) plot(h0) plot(hint) plot(h) points(eps*h0,col="green") points(u*hint,col="red") plot( sort(h) ) points(sort(eps*h0),col="green") points(sort(u*hint),col="red") plot(cols,rows,main="Struktur von H") # end info matrix elements #-----------------------------------# # store many-body hamiltonian # # as a sparse matrix # #-----------------------------------# # wenn man eine Matrix als duennbesetzte Matrix # speichern moechte, muss man das Matrix-package # laden: require(Matrix) sH = sparseMatrix(rows,cols,x=h,dims=c(nB,nB)) sH image(sH) H = as.matrix(sH) # info: sH H mode(sH) mode(H) class(sH) class(H) str(sH) str(H) # Rechenzeiten fuer Matrix-Produkte: calctime = Sys.time() res1 = sH %*% sH %*% sH %*% sH %*% sH Sys.time() - calctime calctime = Sys.time() res2 = H %*% H %*% H %*% H %*% H Sys.time() - calctime res1[1:10,1:10] res2[1:10,1:10] # ok, sind identisch.. ############# E N D E ############## Falls jemand noch sehen moechte, wie die Bilder in MatrixDarstellung-1DHubbardModel.pdf generiert worden sind, hier der genaue Code dazu: #----------------------------------------# # # # Die Funktion < a_ell^+ a_ell >(t) # # # #----------------------------------------# L = 2 N = 20 # Jetzt den Code # set up many body basis # # von oben laufen lassen eps = 1 g = 0.5 # Jetzt den Code # set up Hamiltonian # # von oben laufen lassen T = 500 dtplot = 0.02 dt = 0.01 # required dt depends on N and g M = as.integer(T/dt) Mplot = as.integer(T/dtplot) Tplot = (0:Mplot)*dtplot q = as.integer(dtplot/dt) M Mplot Enj = matrix(0, nrow = L, ncol = Mplot+1 ) Enj[,1] = c( N , rep(0,L-1) ) # Anfangszustand, allgemein ein Vektor # c(n1,n2,...,nL) mit Summe nj = N sH = sparseMatrix(rows,cols,x=h,dims=c(nB,nB)) sId = sparseMatrix(rows,cols,x=id,dims=c(nB,nB)) ReU = as.matrix(sId) ImU = as.matrix(0*sId) sH2 = sH %*% sH sH3 = sH2 %*% sH sH4 = sH2 %*% sH2 Rexpidth = as.matrix( sId - 1/2*dt^2*sH2 + 1/24*dt^4*sH4 ) Iexpidth = as.matrix( dt*sH - 1/6*dt^3*sH3 ) t0 = Sys.time() # monitor calctime for(k in 1:M) { # update U: reU = Rexpidth %*% ReU - Iexpidth %*% ImU imU = Rexpidth %*% ImU + Iexpidth %*% ReU ReU = reU ImU = imU # calc En1(t) if( k%%q == 0) { res = rep(0,L) for(i in 1:nB) { # Initial State := |N,0,...,0> = B[nB,] # wir haben diesen Fall # Initial State := |0,N,0> = B[N+1,] j0 = nB #j0 = N/2+1 m = B[i,] abs2 = ReU[i,j0]^2 + ImU[i,j0]^2 res = res + m * abs2 } # res finished kplot = k/q Enj[,1+kplot] = res if( kplot%%(Mplot/10)==0 || kplot==Mplot/100 ) # info { print( c(kplot,Mplot) ) print( Sys.time() - t0 ) } } } plot(Tplot,Enj[1,],ylim=c(0,N)) En = Enj[1,] for(j in 2:L) { points(Tplot,Enj[j,],col=j) En = En + Enj[j,] readline("press enter to continue..") } points(Tplot,En) # das muss = N sein, wenn das nicht konstant ist, # ist das ein Signal dafuer, dass das dt nicht klein # genug gewaehlt wuerde..