#-----------------------------#
#                             #
#   Temperature Zero Limit:   #
#       ODE,  w1 = 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   ###            

Phi = seq(0,+1,by=0.01)                        # 1.4 min with these phi's 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   ###                       




NPhi = length(Phi)
Phi = matrix(Phi, nrow=ncores, ncol=ceiling(NPhi/ncores), byrow=TRUE )
NPhi
Phi                                            # 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,])
}

chiOn = matrix(0,LM,LM)
for(k in 1:LM)
{
  i = site[k,]
  for(ell in 1:LM)
  {
    j = site[ell,]
    
    if( sum(i) %% 2 == 0 & sum(j) %% 2 == 0 )
    {
      chiOn[k,ell] = 1
    }
    if( sum(i) %% 2 == 1 & sum(j) %% 2 == 1 )
    {
      chiOn[k,ell] = 1
    }
  }
}
chiOff = matrix(1,LM,LM) - chiOn

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          #
#---------------------#

SolveODEperCore = function(Phi, ncore, nt, LM, dt, epsh, u )
{

phi = Phi[ncore,]
nphi = length(phi)
result = vector("list", length = nphi)

absu = abs(u)
epsu = sign(u)
Id   = diag(rep(1,LM))
Null  = diag(rep(0,LM))

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: rhouu -> G,  Fuu -> F

G0 = Null                                    
F0 = Null

nNA = rep(0,nt)

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

  # solve ODE
  for( k in 2:nt )
  {
    pphi = phi[i] * SignDiagLM
    dx = dt * diag( pphi )

    # Calculate DeltaG and DeltaF #

    DG = diag(diag(G))
    dh = -epsh*dt + absu*dt * DG/2 - sqrt(absu) * dx  
    dhG = dh %*% G
    Mat = (Id-F) %*% dhG 
    DeltaG = ( (Id-F) %*% dh %*% (Id+F) - G %*% dhG )/2 
    DeltaF = ( Mat - t(Mat) )/2

    # END Calculate DeltaG and DeltaF #

    Gnew = G + DeltaG 
    Fnew = F + DeltaF

    # Calculate energies, correls and check path #

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

    diagG = diag(G)
    cspin = ( outer(diagG,diagG)*(1+epsu)^2 + 2*F^2 - 2*G^2 )/4      # LM x LM  Matrix
    cpair = ( F + G )*( F - epsu*G )/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(phi[i], Wt, Wepst, Wut, Cspint, Cpairt, sum(nNA) )

} # next i in 1:nphi

return(result)

} ##  END SolveSDEperCore  ##


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


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



phi  = rep(0,NPhi)
V    = rep(0,NPhi)
Wbar = rep(0,NPhi)
WT   = rep(0,NPhi)
sumNA = 0

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

plot(phi,V)

imin   = which.min(V)
phimin = phi[imin]
Wmin   = WT[imin]

phimin
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
for(j in 2:LM)
{
  CpairT[j] = CpairT[j] * (chiOn[1,j] - epsu*chiOff[1,j])
}
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)