# Frenkel-Kontorova Model
#Please keep in mind to run the model multiple times. Sometimes the restoration force or time steps
#aren't enough to restore the lattice back into the lattice back.
import random
import matplotlib.pyplot as plt
XMAX = 15
YMAX = 15
delta = 0.01
dt = 0.1
tmax = 1500
def actualx(x): return x % XMAX
def actualy(y): return y % YMAX
def metric(x, x2): return min(abs(x - x2), XMAX - abs(x - x2))
class Particle:
def __init__(self, x, y, vx, vy, q, m):
self.x = x
self.y = y
self.vx = vx
self.vy = vy
self.q = q
self.m = m
def potential(self, x, y):
dx = metric(x, self.x)
dy = metric(y, self.y)
r = (dx**2 + dy**2)**0.5
return self.q / r if r != 0 else 0
class Defect:
def __init__(self, x, y, strength):
self.x = x
self.y = y
self.strength = strength
def force_on(self, px, py):
dx = metric(px, self.x)
dy = metric(py, self.y)
r2 = dx**2 + dy**2
if r2 == 0: return 0, 0
return (-0.5*self.strength * dx / r2, -0.5*self.strength * dy / r2)
def chargelistgenerator(m, q, xmax, ymax, vxmax, vymax, numberofparticles):
return [Particle(random.uniform(0, xmax), random.uniform(0, ymax),
random.uniform(-vxmax, vxmax), random.uniform(-vymax, vymax), q, m)
for _ in range(numberofparticles)]
def generate_defects(num_defects, strength):
return [Defect(random.uniform(0, XMAX), random.uniform(0, YMAX), strength)
for _ in range(num_defects)]
def gridPotentialAtPoint(x, y, chargelist):
return sum(charge.potential(x, y) for charge in chargelist)
def gridForceAtPointX(x, y, q, chargelist, Edrive):
return q*Edrive - (q * (gridPotentialAtPoint(actualx(x + delta), actualy(y), chargelist) -
gridPotentialAtPoint(actualx(x - delta), actualy(y), chargelist)) / (2 * delta))
def gridForceAtPointY(x, y, q, chargelist):
return - (q * (gridPotentialAtPoint(actualx(x), actualy(y + delta), chargelist) -
gridPotentialAtPoint(actualx(x), actualy(y - delta), chargelist)) / (2 * delta))
def particleVelocityUpdate(particle, chargelist, defects, Edrive):
fx = gridForceAtPointX(particle.x, particle.y, particle.q, chargelist, Edrive)
fy = gridForceAtPointY(particle.x, particle.y, particle.q, chargelist)
#account for force from defects
for defect in defects:
dfx, dfy = defect.force_on(particle.x, particle.y)
fx += dfx
fy += dfy
particle.vx += fx / particle.m
particle.vy += fy / particle.m
particle.vx *= 0.999
particle.vy *= 0.999
def particlePositionUpdate(particle, deltat):
particle.x = actualx(particle.x + particle.vx * deltat)
particle.y = actualy(particle.y + particle.vy * deltat)
def run_simulation(Edrive):
particles = chargelistgenerator(m=1, q=0.5, xmax=XMAX, ymax=YMAX,
vxmax=0.05, vymax=0.05, numberofparticles=20)
defects = generate_defects(num_defects=11, strength=0.01)
t = 0
while t < tmax:
for p in particles:
particleVelocityUpdate(p, particles, defects, Edrive)
for p in particles:
particlePositionUpdate(p, dt)
t += dt
return particles, defects
particles_before, defects = run_simulation(Edrive=0)
particles_after, _ = run_simulation(Edrive=450)
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
for i, (particles, title) in enumerate(zip([particles_before, particles_after],
['Before Edrive', 'After Edrive'])):
x_vals = [p.x for p in particles]
y_vals = [p.y for p in particles]
dx_vals = [d.x for d in defects]
dy_vals = [d.y for d in defects]
axs[i].scatter(x_vals, y_vals, c='blue', label='Particles')
axs[i].scatter(dx_vals, dy_vals, c='red', marker='x', s=100, label='Defects')
axs[i].set_xlim(0, XMAX)
axs[i].set_ylim(0, YMAX)
axs[i].set_title(title)
axs[i].set_aspect('equal')
axs[i].grid(True)
plt.suptitle("Frenkel-Kontorova: Before and After Applying Driving Force", fontsize=16)
plt.show()
No comments:
Post a Comment