#-----------------------------------#
#                                   #
#     Exact Diagonalization of      #                          
#    Grand Canonical Hamiltonian    #
#                                   #
#-----------------------------------#



##  Start Input Variables  ##

L = 3
M = 2                                    # L x M - Gitter

u = 4
eps = -1.0

beta = 4
dt = 0.01

##   End Input Variables   ##



tk = seq(0,beta,by=dt)
nt = length(tk)

LM = L*M
LMS = 2*L*M                              # Anzahl Gitterpunkte mit Spin, 2*|Gamma|

site = matrix(0,nrow=LMS,ncol=2)         # ncol = d = 2 Raum-Dims
spin = rep("",LMS)
k = 1
for(jy in 1:M)
{
  for(jx in 1:L)
  {
    for(s in c("up","down"))
    {
      site[k,] = c(jx,jy)
      spin[k] = s
      k = k+1
    }
  }
}


H0   = vector(mode = "list", length = LMS)
Hint = vector(mode = "list", length = LMS)
Id   = vector(mode = "list", length = LMS)



#-------------------------------------#
#                                     #
#       LOOP GRAND CANONICAL H        #
#                                     #
#-------------------------------------#
t0 = Sys.time()                          # monitor calc time
for( N in 1:LMS )                        # START Loop Grand Canonical
{


#-------------------------------#
#   set up many body basis      #
#-------------------------------#

nB = choose(LMS,N)                       # = |B|
nB                                       # etwa 12870 fuer L=4, M=2, N=8
                                         # etwa  1820 fuer L=4, M=2, N=4
B0 = t(combn(LMS,N))                
B = matrix(0,nrow=nB,ncol=LMS) 
for(n in 1:nB)
{
  b0 =  B0[n,]                       
  b  =  B[n,] 
  b[b0] = 1
  B[n,] = b
}
B                            


#-------------------------------#
#    set up Hamiltonian H_N     #
#-------------------------------#

# allocate some space 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 )                           

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:LM)
      {
        res = res + n[2*k-1]*n[2*k]
      }
      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 )                          
    if( length(k)==2 )
    {
      k1 = k[1]
      k2 = k[2] 
      Site1 = site[k1,]
      Site2 = site[k2,]
      abstand = sum( (Site1-Site2)^2 )

      if( d[k1]*d[k2] == -1 & abstand==1 & spin[k1]==spin[k2] )       # we have a non zero contribution
      {
        if( m[k1] == n[k1]+1 )
        {
          # ERST: wir vernichten ein Teilchen bei k2
          # DANN: wir erzeugen ein Teilchen bei k1:
          n_less = cumsum(n) - n
          sign_k2 = (-1)^(n_less[k2])
          nn = n 
          nn[k2] = n[k2] - 1
          nn_less  = cumsum(nn) - nn
          sign_k1 = (-1)^(nn_less[k1])
          res = sign_k1 * sign_k2 #sqrt(n[k1]+1)*sqrt(n[k2])

          h0[nz] = res
          hint[nz] = 0 
          id[nz]   = 0
          rows[nz] = i
          cols[nz] = j
          nz = nz + 1
        }
        if( m[k1] == n[k1]-1 )
        {
          # ERST: wir vernichten ein Teilchen bei k1
          # DANN: wir erzeugen ein Teilchen bei k2:
          n_less = cumsum(n) - n
          sign_k1 = (-1)^(n_less[k1])
          nn = n 
          nn[k1] = n[k1] - 1
          nn_less  = cumsum(nn) - nn
          sign_k2 = (-1)^(nn_less[k2])
          res = sign_k1 * sign_k2 #sqrt(n[k1])*sqrt(n[k2]+1)

          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,]

} # 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)]
}


#-----------------------------------#
# use simple matrix representation  #
#    no sparse functionalities      #
#-----------------------------------#
hh0   = matrix(0,nB,nB)
hhint = matrix(0,nB,nB)

for(nz in 1:length(h0))
{
  hh0[rows[nz],cols[nz]]   = h0[nz]
  hhint[rows[nz],cols[nz]] = hint[nz]
}

H0[[N]]   = hh0
Hint[[N]] = hhint
Id[[N]]   = diag(rep(1,nB))

print( c( N , nz , Sys.time() - t0 ) )    

} ## NEXT N in 1:LMS, End Loop Grand Canonical ##

Sys.time() - t0
#-------------------------------------#
#                                     #
#     END LOOP GRAND CANONICAL H      #
#                                     #
#-------------------------------------#



# Grand Canonical:  15.5 min for L=4, M=2
# Grand Canonical:    10 sec for L=3, M=2




#-------------------------------#
#       Grand Canonical         #
#     Expectation Value H       #
#-------------------------------#


TrexpH = rep(1,nt)
TrHexpH = rep(0,nt)


for( N in 1:LMS )
{
  # equation (4.26) of section 4.2:
  HN = eps*H0[[N]] + u*Hint[[N]] - u/2 * N * Id[[N]] + u/4 * LM * Id[[N]]

  res = eigen(HN, symmetric=TRUE, only.values=TRUE)
  evalsN = res[[1]]

  for(k in 1:nt)
  {
    TrexpHN = sum( exp( -tk[k]*evalsN ) )
    TrHNexpHN = sum( evalsN * exp( -tk[k]*evalsN ) )

    TrexpH[k] = TrexpH[k] + TrexpHN
    TrHexpH[k] = TrHexpH[k] + TrHNexpHN
  }
}

EH = TrHexpH / TrexpH /LM 
plot(tk,EH)