This is a Julia implementation of the pure diffusion Monte Carlo method. The code is inefficient and messy, so please use it only as a reference.
# 概要：拡散モンテカルロ法
# 環境：Windows 10 Pro x64 (i7-4650U, 8.00GB)
# 言語：Julia (Version 1.5.0)
# 作成：2020-05-19 (yyyy-mm-dd)
# 更新：2020-08-05 (yyyy-mm-dd)
# 更新：2020-12-10 (yyyy-mm-dd)
# 文献：https://doi.org/10.1119/1.18168
# 文献：https://doi.org/10.1063/1.4822960
# 共通関数：
function Reset()
println(Name," Run")
global BIN_WIDTH = 0.001
global histogram = [[1.4,1.4+BIN_WIDTH,0,0]] # 左端, 右端, num, num^2
global S .= 1
global S[N_ini+1:end] .= 0
global E_ref = Averaged_Potential_Energy()
global N_rep = count(S.==1)
global Info_File = open( PROGRAM_FILE[1:end-3] *"_"* Name * "_Info.csv", "w" )
global Stat_File = open( PROGRAM_FILE[1:end-3] *"_"* Name * "_Stat.csv", "w" )
write( Info_File, "#Parameter,Value" )
write( Info_File, "\nNumber of Time Steps," * string(h_max) )
write( Info_File, "\nNumber of Dimension," * string(i_max) )
write( Info_File, "\nNumber of Replicas (Maximum)," * string(j_max) )
write( Info_File, "\nNumber of Particles," * string(k_max) )
write( Info_File, "\nTime Step Width," * string(dt) )
write( Info_File, "\nInitial Number of Replicas," * string(N_rep) )
write( Info_File, "\nInitial Reference Energy," * string(E_ref) )
write( Info_File, "\nBin Width of Histgram," * string(bin_w) )
write( Info_File, "\nIndex for Exports," * replace(string(Out), ","=>" ") )
write( Info_File, "\nStandard Deviation of Rand," * replace(string(D), ","=>" ") )
write( Info_File, "\nCharge of Particles," * replace(string(Q), ","=>" ") )
write( Info_File, "\nInitial Status," * replace(string(S), ","=>" ") )
write( Info_File, "\nInitial Displacement," * replace(string(R), ","=>" ") )
close( Info_File )
Label_Line = "#i,t,N_rep,E_ref"
for k1 in 1:k_max-1 # n(n-1)/2通り
for k2 in k1+1:k_max
Label_Line *= ",r" * string(k1) * string(k2)
end
end
write( Stat_File, Label_Line )
end
function List_Export(data_number)
List_File = open( string( PROGRAM_FILE[1:end-3], "_", Name, "_List_", lpad(data_number,5,"0"), ".csv" ) , "w" )
Label_Line = "#i"
for k in 1:k_max
Label_Line *= string(",x", k, ",y", k, ",z", k)
end
write( List_File, Label_Line)
for j in findall(S.==1)
Data_Line = string("\n", j)
for k in 1:k_max
Data_Line *= string(",", R[1,j,k], ",", R[2,j,k], ",", R[3,j,k])
end
write( List_File, Data_Line )
end
close( List_File )
end
function Hist_Update()
Alive = findall(S.==1)
Dist = Distance(1,2)
for x in Dist[Alive]
while !(histogram[1][1] <= x && x < histogram[end][2])
if histogram[1][1] > x
pushfirst!(histogram, [histogram[1][1]-BIN_WIDTH, histogram[1][1], 0, 0])
end
if histogram[end][2] <= x
push!(histogram, [histogram[end][2], histogram[end][2]+BIN_WIDTH, 0, 0])
end
end
# histogram[Int(floor((x - histogram[1][1])/BIN_WIDTH))+1][3] += 1
# histogram[Int(floor((x - histogram[1][1])/BIN_WIDTH))+1][4] += histogram[Int(floor((x - histogram[1][1])/BIN_WIDTH))+1][3]^2
for interval in histogram
if interval[1] <= x && x < interval[2]
interval[3] += 1
interval[4] = interval[3]^2
end
end
end
end
function Hist_Export()
open(string( PROGRAM_FILE[1:end-3], "_", Name, "_Hist.csv" ), "w") do file
println(file, "#left,right,num,num^2")
for i in 1:length(histogram)
println(file, histogram[i][1], ",", histogram[i][2], ",", histogram[i][3], ",", histogram[i][4])
end
end
X = [(interval[1]+interval[2])/2 for interval in histogram]
P = [interval[4] for interval in histogram]
println("DMC : <R> = ", sum(X.*P/sum(P)))
end
function Stat_Export(countdown)
Data_Line = string("\n", countdown, ",", t, ",", N_rep, ",", E_ref)
for k1 in 1:k_max-1 # n(n-1)/2通り
for k2 in k1+1:k_max
Data_Line *= "," * string( sum(Distance(k1,k2)[findall(S.==1)]) / count(S.==1) )
end
end
write( Stat_File, Data_Line )
end
# function GaussianRand(n,m) # 平均0分散1の正規分布に従う乱数をn,m個並べた配列を返す
# return sqrt.(-2.0*log.(1 .-rand(n,m))).*sin.(2.0*pi*rand(n,m))
# end
function Diffusion_Process()
global t += dt
for k in 1:k_max
global R[:,:,k] += D[k] * randn(Float64,i_max,j_max)
# global R[:,:,k] += D[k] * GaussianRand(i_max,j_max)
end
end
#本当に合ってるのか？
function Distance(num1,num2) # 粒子間距離の配列
dX = R[1,:,num1] - R[1,:,num2]
dY = R[2,:,num1] - R[2,:,num2]
dZ = R[3,:,num1] - R[3,:,num2]
return sqrt.(dX .* dX + dY .* dY + dZ .* dZ)
end
function Potential_Energy() # 死んだレプリカも含まれるので注意
U = zeros(Float64,j_max)
for k1 in 1:k_max-1 #findall([k for k in 1:k_max] .!== k1)
for k2 in k1+1:k_max # n(n-1)/2通りの相互作用
U += Q[k1]*Q[k2] ./ Distance(k1,k2)
end
end
return U
end
function Averaged_Potential_Energy()
return sum( Potential_Energy()[findall(S.==1)] ) / count(S.==1)
end
function Birth_Replica(parent_replica)
first_dead = findfirst(S.==0)
# println(parent_replica, "->",first_dead)
global S[first_dead] = 2
for k in 1:k_max
global R[:,first_dead,k] = R[:,parent_replica,k]
end
end
function Branch_Process()
Weight = exp.(-(Potential_Energy().-E_ref)*dt)
mn = floor.(Weight+rand(j_max))
Alive = findall(S.==1)
for j in intersect( findall(mn.==0) , Alive )
global S[j] = 0 # death
end
for j in intersect( findall(mn.==2) , Alive )
Birth_Replica(j)
end
for j in intersect( findall(mn.==3) , Alive )
Birth_Replica(j)
Birth_Replica(j)
end
for j in findall(S.==2) # birth mark
global S[j] = 1
end
end
function Energy_Calculation()
local N_old = N_rep
global N_rep = count(S.==1)
global E_ref += 1.0/dt*(1.0-(N_rep+1.0)/(N_old+1.0)*(N_old+10.0*j_max)/(N_ini+10.0*j_max))
# global E_ave = (1.0-F_ave)*E_ave + F_ave*E_ref
# 参照エネルギーの時間平均が系のエネルギーを表すので移動平均を取れば良い.
end
function Finish()
close( Stat_File )
println(Name," Finish")
end
function Run()
Reset()
Stat_Export(0)
for h in 1:h_max
Diffusion_Process()
Branch_Process()
Energy_Calculation()
Stat_Export(h)
if h > 5000
Hist_Update()
end
end
List_Export(h_max)
Hist_Export()
Finish()
end
# 宣言：
# Out = union([0],[10^i for i in 0:7],5*[10^i for i in 0:6],[100*i for i in 1:10],[1000*i for i in 1:10],[10000*i for i in 1:10]) # Index for Exports e.g. [1 5 50 100 500 1000 5000 10000]
Out = [1 10 100 1000 5000 15000 25000 35000 45000]
const h_max = 25000 # Number of Time Steps
const i_max = 3 # Number of Dimension
const j_max = 8800 # Number of Replicas (Maximum)
const N_ini = 8000 # Initial Number of Replicas
const dt = 0.1 # Time Step Width
const bin_w = 0.2 # Bin Width of Histgram
# 原子：
Names=["Ps","Mu","H","D","T"]
Diffs=[sqrt(dt/1.0000000),sqrt(dt/206.7682830),sqrt(dt/1836.15267343),sqrt(dt/3670.48296788),sqrt(dt/5496.92153573)]
for system in 1:5
global Name = Names[system] # Name of System
global k_max = 2 # Number of Particles (Nuclear & Electron)
global E_ref = 0.0 # Reference Energy
global t = 0.0 # Time
global N_rep = N_ini # Current Number of Replicas
global Q = ones(Float64,k_max) # Charge of Particles
global D = ones(Float64,k_max) # Diffusion Coefficient of Particles
global S = ones(Int64,j_max) # Status of Replicas / 0-not alive, 1-alive, 2-birth mark
global R = (rand(Float64,i_max,j_max,k_max).-0.5)*2 # Displacement of Replicas
global Q[1] = +1.0 # Nuclear
global Q[2] = -1.0 # Electron
global D[1] = Diffs[system] # Positron(1) Muon(206.7682830), Proton(1836.15267343), Deutron(3670.48296788), Triton(5496.92153573)
global D[2] = sqrt(dt) # Electron(1)
global R[:,:,1] .= 0 # Initial Displacement
@time Run()
end
This is a Julia implementation of the pure diffusion Monte Carlo method. The code is inefficient and messy, so please use it only as a reference.