#-----------------------------# # # # SDEs from Theorem 4.2 # # with correl functions # # w2 = 1 # #-----------------------------# require(doParallel) require(doRNG) ncores = detectCores() # 20 on my pc ncores = 16 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)} ) set.seed(123456) ### Start Input Variables ### npath = 40000 # 2.8 minutes for 40K on my pc L = 4 M = 4 # L x M Gitter u = 4 vep = -1.0 # varepsilon, hopping parameter T = 2 dt = 0.01 bc = "pbc" # boundary conditions #bc = "none" # In Addition, the quantity minint2W has to be precalculated with # 2D-Monte-Carlo-Main-Theorem-w2=1.txt with, say, 10K MC is sufficient ### End Input Variables ### tk = seq(0,T,by=dt) nt = length(tk) 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 } } 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 epsu = sign(u) absu = abs(u) 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 } #---------------------# # SDE # #---------------------# # the quantity minintW is used only for numerical reasons: # # sum_i[ G_i * exp( -intW_i ) ] / sum_i[ exp( -intW_i ) ] = # sum_i[ G_i * exp( -(intW_i - minintW) ) ] / sum_i[ exp( -(intW_i - minintW) ) ] # # and has to be precalculated with 2D-Monte-Carlo-Main-Theorem-w2=1.txt SolveSDEperCorePlainMC2 = function(npath, nt, LM, dt, epsh, u, minint2W ) { LMS = 2*LM sqrtdt = sqrt(dt) absu = abs(u) epsu = sign(u) Id = diag(rep(1,LM)) Null = diag(rep(0,LM)) Id2 = diag(rep(1,LMS)) Null2 = diag(rep(0,LMS)) epsepsh = rbind( cbind( epsh , Null ) , cbind( Null , epsh ) ) SumWepsZk = rep(0,nt) SumWuZk = rep(0,nt) SumZk = rep(0,nt) SumCspinZk = matrix(0, nrow=LM, ncol=nt) SumCpairZk = matrix(0, nrow=LM, ncol=nt) # Notation: rho -> G, F -> F G0 = Null2 F0 = Null2 nNA = rep(0,nt) for( i in 1:npath) { G = G0 F = F0 intW = 0 Zk = 1 # solve SDE for( k in 2:nt ) { xi = rnorm(LM) dy = diag( sqrtdt*xi ) 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 MC-Sums # G = Gnew F = Fnew W = Weps + Wu SumWepsZk[k] = SumWepsZk[k] + Weps*Zk SumWuZk[k] = SumWuZk[k] + Wu*Zk SumZk[k] = SumZk[k] + Zk SumCspinZk[,k] = SumCspinZk[,k] + Cspin*Zk SumCpairZk[,k] = SumCpairZk[,k] + Cpair*Zk intW = intW + W*dt Zk = exp( -(intW - minint2W[k]) ) } # next k } # next path result = list(SumWepsZk, SumWuZk, SumZk, SumCspinZk, SumCpairZk, sum(nNA) ) } # END SolveSDEperCorePlainMC2 # AddResMC = function(list1,list2) { res1 = list1[[1]] + list2[[1]] res2 = list1[[2]] + list2[[2]] res3 = list1[[3]] + list2[[3]] res4 = list1[[4]] + list2[[4]] res5 = list1[[5]] + list2[[5]] res6 = list1[[6]] + list2[[6]] res = list(res1,res2,res3,res4,res5,res6) return(res) } npercore = as.integer(npath/ncores) calctime = Sys.time() res = foreach(i=1:ncores, .combine='AddResMC') %dorng% SolveSDEperCorePlainMC2(npercore,nt,LM,dt,epsh,u,minint2W) Sys.time() - calctime SumWepsZk = res[[1]] SumWuZk = res[[2]] SumZk = res[[3]] SumCspinZk = res[[4]] SumCpairZk = res[[5]] nNA = res[[6]] nNA E2Weps = SumWepsZk / SumZk /LM E2Wu = SumWuZk / SumZk / LM E2Cspin = matrix(0,LM,nt) E2Cpair = matrix(0,LM,nt) for(j in 2:LM) { E2Cspin[j,] = SumCspinZk[j,] / SumZk E2Cpair[j,] = SumCpairZk[j,] / SumZk } E2Cspin[1,] = rep(0,nt) E2Cpair[1,] = rep(0,nt) E2W = E2Weps + E2Wu # farbpaletten: ramp = colorRamp(c("blue","cyan")) mycolB = rgb( ramp(seq(0,1,length = 8)), max = 255) ramp = colorRamp(c("darkred","magenta")) mycolR = rgb( ramp(seq(0,1,length = 8)), max = 255) ## Energies Weps, Wu, W ## plot(tk, E2W ) points(tk, E2Weps, col="green") points(tk, E2Wu, col="red") points(tk, -absu/4 *tanh(tk*absu/4), cex=0.5, col="turquoise" ) # atomistic limit #points(tk, c(0,minint2W[-1]/LM/tk[-1]), col="yellow") ## spin - spin ## info = "all 15 spin-spin correlations on a 4x4 lattice" plot(tk, E2Cspin[1,], type="l", ylim=c(-0.12,0.12), main=info, font.main=10, xlab="beta", ylab="Cspin( i , j )" ) iB=1 iR=1 for(j in 2:LM) { if( sum(site[j,])%%2 == 0 ) { points(tk, E2Cspin[j,], col=mycolR[iR], type="l" ) iR = iR+1 } if( sum(site[j,])%%2 == 1 ) { points(tk, E2Cspin[j,], col=mycolB[iB], type="l" ) iB = iB+1 } } ## pair - pair ## info = "all 15 pair-pair correlations on a 4x4 lattice" plot(tk, E2Cpair[1,], type="l", ylim=c(-0.06,0.06), main=info, font.main=10, xlab="beta", ylab="Cpair( i , j )" ) iB=1 iR=1 for(j in 2:LM) { if( sum(site[j,])%%2 == 0 ) { points(tk, E2Cpair[j,], col=mycolR[iR], type="l" ) iR = iR+1 } if( sum(site[j,])%%2 == 1 ) { points(tk, E2Cpair[j,], col=mycolB[iB], type="l" ) iB = iB+1 } } stopCluster(mycl)