#-----------------------------# # # # SDEs from Theorem 4.2 # # with correl functions # # w1 = 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 # 1.1 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 minint1W has to be precalculated with # 2D-Monte-Carlo-Main-Theorem-w1=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-w1=1.txt SolveSDEperCorePlainMC = function(npath, nt, LM, dt, epsh, u, minint1W ) { sqrtdt = sqrt(dt) absu = abs(u) epsu = sign(u) Id = diag(rep(1,LM)) Null = diag(rep(0,LM)) 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: rhouu -> G, Fuu -> F G0 = Null F0 = Null nNA = rep(0,nt) for( i in 1:npath) { G = G0 F = F0 intW = 0 Zk = 1 # solve SDE for( k in 2:nt ) { phi = rnorm(LM) dx = diag( sqrtdt*phi ) # 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 # (chi_on - epsu*chi_off) applied below 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 - minint1W[k]) ) } # next k } # next path result = list(SumWepsZk, SumWuZk, SumZk, SumCspinZk, SumCpairZk, sum(nNA) ) } # END SolveSDEperCorePlainMC # 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% SolveSDEperCorePlainMC(npercore,nt,LM,dt,epsh,u,minint1W) Sys.time() - calctime SumWepsZk = res[[1]] SumWuZk = res[[2]] SumZk = res[[3]] SumCspinZk = res[[4]] SumCpairZk = res[[5]] nNA = res[[6]] nNA E1Weps = SumWepsZk / SumZk /LM E1Wu = SumWuZk / SumZk / LM E1Cspin = matrix(0,LM,nt) E1Cpair = matrix(0,LM,nt) for(j in 2:LM) { E1Cspin[j,] = SumCspinZk[j,] / SumZk E1Cpair[j,] = SumCpairZk[j,] / SumZk * (chiOn[1,j] - epsu*chiOff[1,j]) } E1Cspin[1,] = rep(0,nt) E1Cpair[1,] = rep(0,nt) E1W = E1Weps + E1Wu # 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, E1W ) points(tk, E1Weps, col="green") points(tk, E1Wu, col="red") points(tk, -absu/4 *tanh(tk*absu/4), cex=0.5, col="turquoise" ) # atomistic limit #points(tk, c(0,minint1W[-1]/LM/tk[-1]), col="yellow") ## spin - spin ## info = "all 15 spin-spin correlations on a 4x4 lattice" plot(tk, E1Cspin[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, E1Cspin[j,], col=mycolR[iR], type="l" ) iR = iR+1 } if( sum(site[j,])%%2 == 1 ) { points(tk, E1Cspin[j,], col=mycolB[iB], type="l" ) iB = iB+1 } } ## pair - pair ## info = "all 15 pair-pair correlations on a 4x4 lattice" plot(tk, E1Cpair[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, E1Cpair[j,], col=mycolR[iR], type="l" ) iR = iR+1 } if( sum(site[j,])%%2 == 1 ) { points(tk, E1Cpair[j,], col=mycolB[iB], type="l" ) iB = iB+1 } } stopCluster(mycl)