#--------------------------------------# # # # Week11: Die Matrix-Darstellung # # des 1D Bose-Hubbard Modells # # # #--------------------------------------# Wir schauen uns zunaechst die Sachen aus dem week11.pdf an. Wir wollen die Matrix H aus Gleichung (15) anlegen. Das H ist eine |BN| x |BN| - Matrix, wenn BN die Basis bezeichnet, die (9) im week11 angegeben ist. Wir muessen uns also zunaechst ueberlegen, 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 (9) von week11.pdf folgendermassen generieren: #-------------------------------# # set up many body basis # #-------------------------------# L = 4 N = 6 nB = choose(N+L-1,L-1) # = |B| nB # etwa 84 fuer L=4, N=6 B0 = t(combn(N+L-1,L-1)) # t = transponieren, dann nB Zeilen B0 # N+L-1 = 6+4-1 = 9 fuer L=4, N=6 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 L-1 Trennstriche, die bei b[1],..,b[L-1] 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 BN aus pdf-Gleichung (9) 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 |BN| = 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 (14) 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 (13).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 (13).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 head(sH) image(sH) H = as.matrix(sH) # als normale Matrix # info: sH[1:10,1:10] H[1:10,1:10] mode(sH) mode(H) class(sH) class(H) str(sH) str(H) #---------------------------------# # Eigenwerte und Eigenvektoren: # #---------------------------------# Fuer n x n Matrizen mit n <= 10000 ist die eigen()-Funktion sehr performant, wenn man eine R-Version mit Intel-MathKernelLibrary (MKL) benutzt: (die Rechenzeiten fuer sH oder H sind im wesentlichen identisch) calctime = Sys.time() res = eigen(sH, symmetric=TRUE, only.values=TRUE) Sys.time() - calctime # etwa 0.9 sec mit MKL (10 cores), # etwa 9.2 sec ohne MKL calctime = Sys.time() res = eigen(sH, symmetric=TRUE) Sys.time() - calctime # etwa 2.3 sec mit MKL, # etwa 34.4 sec ohne MKL calctime = Sys.time() res = eigen(sH) Sys.time() - calctime # etwa 2.5 sec mit MKL, # etwa 34.6 sec ohne MKL evals = res[[1]] # res[[2]] ist die Matrix der Eigenvektoren plot(evals) plot(rev(evals)) summary(evals) sH0 = sparseMatrix(rows, cols, x = h0, dims = c(nB,nB) ) res0 = eigen(sH0, symmetric=TRUE, only.values=TRUE) evals0 = rev(res0[[1]]) plot(evals0) summary(evals0) # checken wir eben noch gegen die exakten, analytischen # Eigenwerte von H0: wir hatten die Vielteilchen-Basis BN # durch die Matrix B codiert und jede Zeile von B definiert # einen Basisvektor |n1,...,nL> = B[k,] mit k=1,..,nB; # die Einteilchen-EW's lambda sind durch Aufg1 vom UeBlatt2 # gegeben: lambda = 2*cos((1:L)*pi/(L+1)) evexakt = rep(0,nB) for(k in 1:nB) { n = B[k,] E = 0 # many body energy eigenvalue for(j in 1:L) { E = E + n[j]*lambda[j] } evexakt[k] = E } evexakt = sort(evexakt) plot(evals0) points(evexakt,col="red") # das passt summary(evals0-evexakt) # sind numerische Nullen..