#-----------------------------#
#                             #
#   Temperature Zero Limit:   #
#       ODE,  w2 = 1          #
#                             #
#-----------------------------#


require(doParallel)
require(doRNG)

ncores = detectCores()                         # 20 on my pc           
ncores = 15                             
mycl = makeCluster(ncores)
registerDoParallel(mycl)

# if you are using Microsoft R Open with Intel MathKernelLibrary:
# set the number of threads to 1 for each cluster process:

clusterEvalQ( cl = mycl , {setMKLthreads(1)} )




###   Start Input Variables   ###            

Xi = seq(0,0.5,by=0.02)                        # for these xi's: 2.6 min on my pc
                       
L = 12
M = 12                                         # L x M Gitter

u = 4
vep = -1.0                                     # varepsilon, hopping parameter

T = 20                            
dt = 0.01

bc = "pbc"                                     # boundary conditions
#bc = "none"

###   End Input Variables   ###                       




NXi = length(Xi)
Xi = matrix(Xi, nrow=ncores, ncol=ceiling(NXi/ncores), byrow=TRUE )
NXi
Xi                                           # kann man so machen..

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

absu = abs(u)
epsu = sign(u)                          

LM   = L*M

site = matrix(0,nrow=LM,ncol=2)
k = 1
for(jy in 1:M)
{
  for(jx in 1:L)
  {
    site[k,] = c(jx,jy)
    k = k+1
  }
}

SignDiagLM = rep(0,LM)
for(k in 1:LM)
{
  SignDiagLM[k] = (-1)^sum(site[k,])
}

eps = matrix(0,LM,LM) 
eps0 = matrix(0,LM,LM)                                    

for(k in 1:LM)
{
  for(ell in 1:LM)
  {
    abstand = sum( (site[k,]-site[ell,])^2 )
    if( abstand==1  )
    {
      eps[k,ell] = vep 
      eps0[k,ell] = vep
    }

    # periodic bc's
    d = site[k,]-site[ell,]
    dx = d[1]
    dy = d[2]
    if( abs(dx) == L-1 & dy==0 )
    {
      eps[k,ell] = vep 
    }
    if( dx==0 & abs(dy) == M-1 )
    {
      eps[k,ell] = vep 
    }
    # end periodic bc's
  }
}

if(bc=="none")
{
  epsh = eps0
}
if(bc=="pbc")
{
  epsh = eps
}



#---------------------#
#        ODE          #
#---------------------#

SolveODEperCore2 = function(Xi, ncore, nt, LM, dt, epsh, u )
{

xi = Xi[ncore,]
nxi = length(xi)
result = vector("list", length = nxi)
LMS = 2*LM

Id   = diag(rep(1,LM))
Null  = diag(rep(0,LM))
Id2   = diag(rep(1,LMS))
Null2  = diag(rep(0,LMS))

absu = abs(u)
epsu = sign(u)
epsepsh = rbind( cbind( epsh , Null ) , cbind( Null , epsh ) )

Wepst = rep(0,nt) 
Wut = rep(0,nt)
Wt = rep(0,nt)
Cspint = matrix(0, nrow=LM, ncol=nt)
Cpairt = matrix(0, nrow=LM, ncol=nt)


# Notation: rho -> G,  F -> F

G0 = Null2                                    
F0 = Null2

nNA = rep(0,nt)

for( i in 1:nxi)
{
  G = G0
  F = F0

  # solve ODE
  for( k in 2:nt )
  {
    if( epsu > 0 )
    {  
      xxi = xi[i] * SignDiagLM  
    }
    if( epsu <= 0 )
    {  
      xxi = xi[i] * rep(1,LM)  
    }
    dy = dt * diag( xxi )
    dBy = rbind( cbind(Null,dy) , cbind(epsu*dy,Null) ) 

    # Calculate DeltaG and DeltaF #

    Guu = G[1:LM,1:LM]
    Gud = G[1:LM,(LM+1):LMS]
    Fuu = F[1:LM,1:LM]
    Fud = F[1:LM,(LM+1):LMS]
    Gdu = epsu*Gud
    Fdu = epsu*Fud
    Gdd = Guu
    Fdd = Fuu
    tiF = rbind( cbind( Fdd , Fdu ) , cbind( Fud , Fuu ) )
    DGud = diag(diag(Gud))
    DGdu = diag(diag(Gdu))
    temp1 = cbind( Null , DGud )
    temp2 = cbind( epsu*DGud , Null )
    DG = rbind( temp1 , temp2 )
    dh = -epsepsh*dt + absu*dt * DG/2 - sqrt(absu) * dBy  
    dht = t(dh)
    Gdht = G %*% dht
    Mat = Gdht %*% (Id2+F) 
    DeltaG = ( (Id2-F) %*% dh %*% (Id2+tiF) - Gdht %*% G )/2 
    DeltaF = ( t(Mat) - Mat )/2

    # END Calculate DeltaG and DeltaF #

    Gnew = G + DeltaG 
    Fnew = F + DeltaF

    # Calculate energies, correls and check path #
    Guu = Gnew[1:LM,1:LM]
    Gud = Gnew[1:LM,(LM+1):LMS]
    Fuu = Fnew[1:LM,1:LM]
    Fud = Fnew[1:LM,(LM+1):LMS]

    Weps = ( sum( diag( epsh%*%Guu ) ) )
    Wu   = -absu/4 * sum( diag(Gud)^2 )
    nrgdensity = ( abs(Weps) + abs(Wu) )/LM

    diagGud = diag(Gud)
    cspin = (Fuu^2 - Guu^2 - epsu*Fud^2 + epsu*Gud^2)/2                        # LM x LM matrix
    cpair = (Fuu^2 + Guu^2 - (1+epsu)/2*(Fud^2 + Gud^2) + (1-epsu)/2*outer(diagGud,diagGud) )/4
    Cspin = cspin[1,]                                                  # Cspin(j0,j), j0 = (1,1) 
    Cpair = cpair[1,]                                                  # Cpair(j0,j), j0 = (1,1)

    # KEEP PATH ?
    if( is.na(nrgdensity) || nrgdensity > 10 )
    {
      nNA[k] = nNA[k] + 1
      break              # ignore path, leave k-loop, start new MC-path
    }
    # END Calculate energies, correls and check path # 

    # Update G and F and Wt and Ct #
    G = Gnew
    F = Fnew
    W = Weps + Wu

    Wepst[k] = Weps           
    Wut[k] = Wu
    Wt[k] = W
    Cspint[,k] = Cspin 
    Cpairt[,k] = Cpair

  } # next k

  result[[i]] = list(xi[i], Wt, Wepst, Wut, Cspint, Cpairt, sum(nNA) )

} # next i in 1:nxi

return(result)

} ##  END SolveSDEperCore2  ##


AddResODE = function(list1,list2)
{
  res = c( list1 , list2 )
  return(res)
}


calctime = Sys.time()
result = foreach(i=1:ncores, .combine='AddResODE') %dorng% SolveODEperCore2(Xi,i,nt,LM,dt,epsh,u)
Sys.time() - calctime



xi   = rep(0,NXi)
V    = rep(0,NXi)
Wbar = rep(0,NXi)
WT   = rep(0,NXi)
sumNA = 0

for( i in 1:NXi)
{
  res = result[[i]]
  xi[i] = res[[1]]
  W = res[[2]]
  WT[i] = W[nt]/LM
  Wbar[i] = sum(W)*dt / T / LM
  V[i] = Wbar[i] + xi[i]^2/2
  sumNA = sumNA + res[[7]]
}
sumNA

plot(xi,V)

imin  = which.min(V)
ximin = xi[imin]
Wmin  = WT[imin]

ximin
Wmin




##  spin-spin and pair-pair correlations (L=M)  ##

resMin = result[[imin]]
CspinT = resMin[[5]][,nt]
CpairT = resMin[[6]][,nt]
CspinT[1] = 0
CpairT[1] = 0
cspin = round( matrix(CspinT, nrow=M, ncol=L, byrow=TRUE) , digits=2 ) 
cpair = round( matrix(CpairT, nrow=M, ncol=L, byrow=TRUE) , digits=2 ) 
cspin
cpair






stopCluster(mycl)