Problem-1: A particle undergoes Brownian Motion with displacement rule
∆𝑥~𝑁(0, 2𝐷∆𝑡),
Where the diffusion coefficient is 𝐷 = 1 and the time step is ∆𝑡 = 0.01. Now answer the following
questions.
i) Write a python code to simulate five independent trajectories of Brownian motion from
𝑡 = 0 𝑡𝑜 𝑡 = 2. [Plot all five trajectories on the same graph]
ii) For 𝑀 = 5000 simulated particles, compute the mean squared displacement 〈𝑥 2 (𝑡)〉 at
multiple time points. Plot 〈𝑥 2 (𝑡)〉 vs time t.
Solution:
Given (what the rule means)
You are told the displacement (increment) in one small time step is:
Time setup
You simulate from (t=0) to (t=2).
Step 1: How many steps?
So you will generate 200 random increments.
Step 2: Build time points
You have (N+1) time points (including (t=0)):
Five trajectories
Idea (what we do)
For each trajectory:
Step A1: Generate increments
Generate 200 random values:
So position is the cumulative sum of increments.
Step A3: Plot all 5 on the same graph.
Code (with comments)
import numpy as np
import [Link] as plt
# ---------------------------
# Step 1: Parameters
# ---------------------------
D = 1.0
dt = 0.01
T = 2.0
# ---------------------------
# Step 2: Time grid
# ---------------------------
N = int(T / dt) # number of steps = 200
t = [Link](0, T, N + 1) # time points from 0 to 2
# ---------------------------
# Step 3: Standard deviation of step
# Δx ~ N(0, 2Ddt) => std = sqrt(2Ddt)
# ---------------------------
sigma = [Link](2 * D * dt)
# ---------------------------
# Step 4: Simulate 5 paths
# ---------------------------
num_paths = 5
[Link](figsize=(9, 5))
for i in range(num_paths):
# Step 4.1: Generate random increments (200 of them)
dx = sigma * [Link](N)
# Step 4.2: Convert increments to position using cumulative sum
x = [Link](N + 1) # x(0)=0
x[1:] = [Link](dx) # x(t_k) = sum of dx up to k
# Step 4.3: Plot the trajectory
[Link](t, x, label=f"Trajectory {i+1}")
[Link]("Time t")
[Link]("Position x(t)")
[Link]("Five Brownian Motion Trajectories (D=1, dt=0.01)")
[Link](True)
[Link]()
plt.tight_layout()
[Link]()
iii) Mean Squared Displacement (MSD) for
M = 5000 particles
What is MSD?
At time (t), each particle has a position (x(t)).
MSD is:
Meaning:
square each particle’s position at time (t)
average those squares
Why it’s important
For 1D Brownian motion, the theoretical result is:
So your simulated MSD plot should be close to a straight line with slope 2.
Efficient simulation idea (vectorized)
Instead of simulating particles one-by-one, we simulate all increments in one array:
dx shape: ((M, N)) = (5000, 200)
Then cumulative sum along time axis gives positions.
Code (with comments + theory comparison)
import numpy as np
import [Link] as plt
# ---------------------------
# Step 1: Parameters
# ---------------------------
D = 1.0
dt = 0.01
T = 2.0
N = int(T / dt)
t = [Link](0, T, N + 1)
M = 5000
sigma = [Link](2 * D * dt)
# ---------------------------
# Step 2: Generate all increments at once
# dx has shape (M, N)
# ---------------------------
dx = sigma * [Link](M, N)
# ---------------------------
# Step 3: Convert increments to positions x(t)
# x has shape (M, N+1)
# ---------------------------
x = [Link]((M, N + 1))
x[:, 1:] = [Link](dx, axis=1)
# ---------------------------
# Step 4: Compute MSD at each time point
# msd[k] = mean of x(:,k)^2
# ---------------------------
msd = [Link](x**2, axis=0)
# ---------------------------
# Step 5: Theoretical MSD line = 2Dt
# ---------------------------
msd_theory = 2 * D * t
# ---------------------------
# Step 6: Plot
# ---------------------------
[Link](figsize=(9, 5))
[Link](t, msd, label="Simulated ⟨x²(t)⟩")
[Link](t, msd_theory, linestyle="--", label="Theory: 2Dt")
[Link]("Time t")
[Link]("Mean Squared Displacement ⟨x²(t)⟩")
[Link]("MSD vs Time (M=5000, D=1, dt=0.01)")
[Link](True)
[Link]()
plt.tight_layout()
[Link]()
Problem-2: In a very simple 1D Particle-In-Cell (PIC) simulation, electrons are represented as
particles, and the electric field is defined on grid points.
i) Creates 10 electrons with random positions between 𝑥= 0 and x = 1.
ii) Deposits their charge onto a grid of 20 cells (charge density array).
You have electrons as particles with positions xpx_pxp in [0,1][0,1][0,1].
You have a 1D grid with 20 cells spanning the same domain.
You want to compute a charge density array on the grid by “depositing” each particle’s
charge onto nearby grid points.
To keep it simple (and still PIC-like), we’ll use Cloud-In-Cell (CIC) deposition:
Each particle shares its charge between the two nearest grid points based on distance
(linear weights).
This avoids noisy “all charge into one cell” behavior.
import numpy as np
# ----------------------------
# (i) Create 10 electrons in [0,1]
# ----------------------------
[Link](42) # for reproducibility
Np = 10
L = 1.0
x_p = [Link](Np) * L # uniform in [0,1)
q = -1.0 # electron charge (normalized units)
# ----------------------------
# (ii) Deposit charge on a grid of 20 cells
# ----------------------------
Nc = 20
dx = L / Nc
# grid points (Nc+1 points)
x_grid = [Link](0, L, Nc + 1)
# charge density array on grid points
rho = [Link](Nc + 1)
# Cloud-In-Cell (CIC) deposition
for xp in x_p:
# find left index i such that x_i <= xp < x_{i+1}
i = int([Link](xp / dx))
# special case: if xp == L, clamp to last cell
# (with rand(), xp is < L, but keep this for safety)
if i == Nc:
i = Nc - 1
# left grid point position
x_i = i * dx
# fractional distance inside the cell
f = (xp - x_i) / dx # in [0,1)
# weights to left and right grid points
w_left = 1.0 - f
w_right = f
# deposit charge density contribution
rho[i] += (q * w_left) / dx
rho[i + 1] += (q * w_right) / dx
# ----------------------------
# Display results
# ----------------------------
print("Particle positions (x_p):")
print([Link](x_p, 4))
print("\nCharge density rho on grid points (length Nc+1 = 21):")
print([Link](rho, 4))
# Optional check: total charge should be close to Np*q
# Integral of rho over domain ~ sum(rho)*dx (rough check on grid points)
total_charge_approx = [Link](rho) * dx
print("\nApprox total charge from rho:", total_charge_approx)
print("Expected total charge:", Np * q)
Problem-3:
Consider a 1D Ising chain with energy:
𝐸 = −𝐽 ∑ 𝑆𝑖 𝑆𝑖+1, 𝑆𝑖 = ±1
𝑖
at temperature 𝑇 = 2, 𝐽 = 1.
i) Write Python code to simulate one Metropolis update for a spin.
ii) Simulate 100 spins for 10000 Monte Carlo steps.
Plot magnetization vs steps (explain expected behavior).
Solution:
import numpy as np
def metropolis_update(spins, J, T): // for single spin
"""
Perform one Metropolis update (one attempted spin flip).
"""
N = len(spins)
# Pick a random spin
i = [Link](N)
# Periodic boundary conditions
left = spins[(i - 1) % N]
right = spins[(i + 1) % N]
# Energy change for flipping spin i
delta_E = 2 * J * spins[i] * (left + right)
# Metropolis criterion
if delta_E <= 0:
spins[i] *= -1
else:
if [Link]() < [Link](-delta_E / T):
spins[i] *= -1
// For 100 spins
import numpy as np
import [Link] as plt
# -------------------------
# Parameters
# -------------------------
N = 100
J = 1.0
T = 2.0
MC_steps = 10000
# -------------------------
# Initialize spins randomly
# -------------------------
spins = [Link]([-1, 1], size=N)
magnetization = []
# -------------------------
# Monte Carlo simulation
# -------------------------
for step in range(MC_steps):
# One Monte Carlo step = N Metropolis updates
for _ in range(N):
metropolis_update(spins, J, T)
# Compute magnetization
M = [Link](spins) / N
[Link](M)
# -------------------------
# Plot magnetization
# -------------------------
[Link](figsize=(9, 5))
[Link](magnetization, linewidth=1)
[Link]("Monte Carlo Steps")
[Link]("Magnetization M")
[Link]("Magnetization vs Monte Carlo Steps (1D Ising, T=2)")
[Link](True)
plt.tight_layout()
[Link]()
Problem-4:
In a simple 1D tight-binding model with nearest-neighbor hopping 𝑡 = 1, the energy 7
dispersion relation is:
𝐸(𝑘) = −2𝑡𝑐𝑜𝑠(𝑘)
i) Generate 200 values of k from −𝜋 𝑡𝑜 + 𝜋.
ii) Compute the band energy E(k).
iii) Plot the resulting band structure curve.
-
Solution:
import numpy as np
import [Link] as plt
# Given parameter
t = 1.0
# (1) Generate 200 k values from -pi to +pi
k = [Link](-[Link], [Link], 200)
# (2) Compute E(k) = -2 t cos(k)
E = -2 * t * [Link](k)
# (3) Plot the band structure
[Link](figsize=(8, 5))
[Link](k, E, linewidth=2)
[Link]("k")
[Link]("E(k)")
[Link]("1D Tight-Binding Band Structure (t=1)")
[Link](True)
plt.tight_layout()
[Link]()