Ordnungsgemäß kinetisches Monte -Carlo auf einem 1D -Gitter ausführenPython

Python-Programme
Anonymous
 Ordnungsgemäß kinetisches Monte -Carlo auf einem 1D -Gitter ausführen

Post by Anonymous »

Ich versuche, eine Monte-Carlo-Simulation auf einem eindimensionalen Gitter auszuführen. Auf dem Gitter gibt es Partikel, die mit bestimmten Raten nach rechts oder links hüpfen können, die von den Populationen an der entsprechenden Gitterstelle und seinem Nachbarn abhängen. Auch links und sehr rechte Seiten gibt es ein Reservoir mit einer festen Partikelnummer. System der ODEs und erhalten Sie die erwarteten stationären Populationen an jedem Gitterstandort. Ich versuche jedoch auch, mit der dynamischen/kinetischen Monte-Carlo-Methode (Gillespie-Algorithmus) ähnliche stationäre Werte zu erzielen. Ich stecke gerade fest und verstehe nicht, wo ich mit dem Algorithmus falsch lieg. Wenn jemand eine Erfahrung hat, könnten Sie möglicherweise gemeinsame Fallstricke bestimmen? Ich teile meinen Code unten: Der erste Teil ist der ODES -Solver zum Benchmarking der Ergebnisse, und der zweite Teil ist die Monte -Carlo -Simulation, die eine falsche Lösung gibt. Schematische Darstellung des Problems, das aus Scipost Phys entnommen wurde. 16, 029 (2024)

Code: Select all

import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.pyplot import cm
from scipy.integrate import solve_ivp
matplotlib.rcParams['figure.figsize'] = (10, 7)

###--------DEFINING PARAMETERS-----------###########
lattice_length=10
cmap = cm.terrain(np.linspace(0, 0.8, lattice_length))

lattice=np.linspace(1, lattice_length, lattice_length)
initial_condition=np.zeros(lattice_length)

GammaL=1 #hopping rate to the left
GammaR=0 #hopping rate to the right
kappaL=kappaR=1 #hopping rates between the very left and right lattice sites and their corresponding left and right reservoirs

nRight=1 #right reservoir population
nLeft=0 #left reservoir population

######-----------SOLVING RATE EQUATIONS WITH ODE---------############
time_ODE=np.linspace(0,100,10000)

def sol_fun_eq():

def dndt(t,V):

dndt=np.zeros(lattice_length)

J_01=kappaL*V[0]*(1+nLeft)-kappaL*nLeft*(1+V[0])
J_12=GammaL*V[1]*(1+V[0])-GammaR*V[0]*(1+V[1])

J_Last=kappaR*nRight*(1+V[-1])-kappaR*V[-1]*(1+nRight)
J_LastButOne=GammaL*V[-1]*(1+V[-2])-GammaR*V[-2]*(1+V[-1])

dndt[0]=J_12-J_01
dndt[-1]=J_Last-J_LastButOne

for p in range(1,lattice_length-1):
J_in  = GammaL*V[p+1]*(1+V[p])-GammaR*V[p]*(1+V[p+1])
J_out = GammaL*V[p]*(1+V[p-1])-GammaR*V[p-1]*(1+V[p])
dndt[p]= J_in - J_out
return [dndt]

sol = solve_ivp(dndt, [time_ODE[0], time_ODE[-1]], initial_condition, method='LSODA', t_eval=time_ODE)
return sol

solution_eq=sol_fun_eq()

steady_state=np.zeros(lattice_length)
for p in range(lattice_length):
steady_state[p]=solution_eq.y[p][-1]

# Plot the results in time
plt.figure(figsize=(12, 8))
for p in range(lattice_length):
plt.plot(solution_eq.t, solution_eq.y[p], label=f'$n_{{{p}}}$', color=cmap[p])
plt.title("Time traces")
plt.xlabel("Time")
plt.ylabel("Populations")
plt.show()

# Plot the steady state
plt.figure()
plt.plot(lattice, steady_state)
plt.scatter(lattice, steady_state)
plt.xlabel("Lattice sites")
plt.ylabel("Populations")
plt.title("Steady state populations")
plt.show()

###--------------------MONTE CARLO SIMULATION----------##########
timeSize=500
time_MC=np.linspace(0, timeSize, timeSize)
SampleRange=2000 #number of sampling

def MC_sim(dummy):

n=np.zeros(lattice_length)
nInTime=np.zeros((timeSize, lattice_length))

for tIdx in range(timeSize):

rates=[]
event_list=[]
#calculate boundary rates on the very left and right sites
rate_01=kappaL*n[0]*(1+nLeft) #particle jump rate from site #1 to left reservoir
rate_10=kappaL*nLeft*(1+n[0]) #particle jump rate from left reservoir to site #1

rate_LastR=kappaR*nRight*(1+n[-1]) #particle jump rate from right reservoir to last lattice site
rate_RLast=kappaR*n[-1]*(1+nRight) #particle jump rate from last lattice site to right reservoir

rates.extend([rate_01, rate_10, rate_LastR, rate_RLast])
event_list.extend(["rate_01", "rate_10", "rate_LastR", "rate_RLast"])

for i in range(lattice_length-1):

rate_L_p=GammaL*n[i+1]*(1+n[i]) #particle jump rate from i+1 to i site
rate_R_p=GammaR*n[i]*(1+n[i+1]) #particle jump rate from i to i+1 site

rates.extend([rate_L_p, rate_R_p])
event_list.extend([f"L{i}", f"R{i}"])

sumRates = np.sum(rates)
probabilities=rates/sumRates
events=np.random.choice(len(rates),  p=probabilities)

#update the lattice site
if events==0:
n[0]-=1
elif events==1:
n[0]+=1
elif events==2:
n[-1]+=1
elif events==3:
n[-1]-=1
elif events>3:
event_idx=event_list[events]
transition=event_idx[0]
lattice_idx=int(event_idx[1])
if transition=="L":
n[lattice_idx]+=1
n[lattice_idx+1]-=1
elif transition=="R":
n[lattice_idx]-=1
n[lattice_idx+1]+=1
else:
raise Exception("Coding error")
for i in range(lattice_length):
nInTime[tIdx, i]=n[i] #save lattice populations at each timestep
return n, nInTime

nSummed=np.zeros(lattice_length)
nInTimeSummed=np.zeros((timeSize, lattice_length))

for s in range(SampleRange):

n, nInTime = MC_sim(0)

for i in range(lattice_length):

nSummed[i]+=n[i]
nInTimeSummed[:,i]+=nInTime[:,i]

n_avg=nSummed/SampleRange
nInTimeAvg=nInTimeSummed/SampleRange

plt.figure(figsize=(12,7))
for i in range(lattice_length):
plt.plot(time_MC, nInTimeAvg[:, i], label=f"# {i}")
plt.legend()
plt.title("Averaged time traces")
plt.xlabel("Time (arbitrary)")
plt.ylabel("Populations")
plt.show()

plt.figure()
plt.plot(lattice, n_avg)
plt.xlabel("Lattice sites")
plt.title("Averaged steady state populations")
plt.ylabel("Populations")
plt.show()

Quick Reply

Change Text Case: 
   
  • Similar Topics
    Replies
    Views
    Last post