import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation %matplotlib inline

n = 25 #total number of particles box_l = 10 #box size

def lj_potential(r, epsilon, sigma): return 4*epsilon*(sigma**12/r**12-sigma**6/r**6) r = np.linspace(1,4,100) #radius vector; 2<x<5, 100 points epsilons = [1,2,5] #range of energy min sigma = 1.5 #distance to zero crossing point force_r = [] for epsilon in epsilons: plt.plot(r, lj_potential(r, epsilon, sigma)) plt.xlabel("r (Å)") plt.ylabel("U(rij)") plt.ylim(-6,6) plt.legend(["epsilon=1", "epsilon=2", "epsilon=5"])

def compute_force(r, epsilon, sigma): return 4*epsilon*((12/r**13)*(sigma**12)-(6/r**7)*(sigma**6))

# TODO modify this function to specify the initial temperature def initialize(beta=1.0, m=1.0): xs = np.zeros([n,2]) for i in range(5): for j in range(5): xs[5*i+j, 0] = i * box_l/5 - box_l/2.0 + box_l/10. xs[5*i+j, 1] = j * box_l/5 - box_l/2.0 + box_l/10. vs = np.random.normal(0,m/(beta),n*2).reshape(n,2) # temperature set here! return xs, vs

xs, vs = initialize()

# helper functions for plotting ppa = 72.*8/20 def plot_conf(xs): fig, ax = plt.subplots(figsize=(8,8)) ax.scatter(xs[:,0], xs[:,1], s=ppa**2) ax.set_xlim([-5,5]) ax.set_ylim([-5,5])

plot_conf(xs)

# distance function with periodic boundary conditions def d(x0, x1): dxs = x1-x0 dxs = dxs - box_l * (np.abs(dxs)>box_l/2.0) * np.sign(dxs) return np.sqrt( (np.sum( (dxs)**2 ) ) ) def lj_energy(xs, epsilon=1.0, sigma=1.0): energy = 0. # TODO! for i in range(n): for j in range(i+1,n): r = d(xs[i], xs[j]) en = lj_potential(r, epsilon, sigma) energy += en return energy def lj_force(xs, epsilon=1.0, sigma=1.0): force = np.zeros([n, n, 2]) # TODO! for i in range(n): for j in range(i+1,n): r = d(xs[i], xs[j]) force[i,j] = compute_force(r, epsilon, sigma) force[j,i] = -force[i,j] return force # can do both in one loop, of course def lj_both(xs, epsilon=1.0, sigma=1.0, r_cutoff=3.0): energy = 0 force = np.zeros([n, n, 2]) #en = [] # TODO for i in range(n): for j in range(i+1,n): r = d(xs[i], xs[j]) if r > r_cutoff: continue force[i,j] = compute_force(r, epsilon, sigma) force[i,j] = -force[i,j] en = lj_potential(r, epsilon, sigma) energy += en return force, energy

m = 1. def integrate_step(xs, vs, forces, dt=5e-4): #TODO velocity half step vs = vs + m**-2*forces*dt*0.5 #TODO position step xs = xs + vs*dt # correct for periodic boundary conditions after updating positions xs = xs - box_l * (np.abs(xs)>(box_l/2.)) * np.sign(xs) force, energy = lj_both(xs) #TODO velocity half step vs = vs + m**-2*forces*dt*0.5 #print(xs) return xs, vs, force, energy

def run_traj(xs, vs, n_steps, dt=1e-4): traj = [] velocities = [] # compute the initial forces force, energy = lj_both(xs) for step in range(n_steps): xs, vs, force, energy = integrate_step(xs, vs, force, dt=dt) traj.append(xs.copy()) velocities.append(vs.copy()) return traj, velocities

xs, vs = initialize() traj = run_traj(xs, vs, 5000, dt=0.1)

# set up an animation fig = plt.figure() ax = plt.gca() ax.set_xlim([-5,5]) ax.set_ylim([-5,5]) scat = ax.scatter(xs[:,0], xs[:,1], s=ppa**2) def init(): scat.set_offsets(np.array(traj[0]).reshape(25*25,2)) return scat, def animate(i): scat.set_offsets((traj[i]).reshape(25*25,2)) return scat, anim = FuncAnimation(fig, animate, init_func=init, blit=True, interval=0.05) anim.save("mygif.gif", writer="pillow")

# TODO compute the average kinetic energy as a function of trajectory time for T = 1kb, 2kb, 5kb x1s, v1s = initialize() x2s, v2s = initialize(beta = 1/2) x5s, v5s = initialize(beta = 1/5) def KE_avg(velocity): KE = [] for step in velocity: KE.append(np.mean(0.5*m*step**2)) return KE def run_traj_new(xs, vs, n_steps, dt=5e-4): traj = [] vels = [] # compute the initial forces force, energy = lj_both(xs) for step in range(n_steps): xs, vs, force, energy = integrate_step(xs, vs, force, dt=dt) traj.append(xs.copy()) vels.append(vs.copy()) return traj, vels x1s, v1s = initialize(beta = 1) traj, vel = run_traj_new(x1s, v1s, 5000, dt = 1) values_1 = KE_avg(vel) x2s, v2s = initialize(beta=1/2) traj, v = run_traj_new(x2s,v2s,5000,dt=1) values_2 = KE_avg(v) x5s, v5s = initialize(beta=1/2) traj, v = run_traj_new(x5s,v5s,5000,dt=1) values_5 = KE_avg(v) plt.scatter(np.linspace(0, 4998, 4999), values_1[1:]) plt.scatter(np.linspace(0, 4998, 4999), values_2[1:]) plt.scatter(np.linspace(0, 4998, 4999), values_5[1:])

# TODO compute g(r) from a simulation of duration 5000 steps