From 234a07717dec6dd1c440dedc8758c88e543bc740 Mon Sep 17 00:00:00 2001 From: Trung Hieu NGUYEN Date: Sun, 22 Mar 2020 20:35:00 +0100 Subject: [PATCH] new code --- chuan_harmonic_verlet_vr.py | 202 ++++++++++++++++++++++++++++++++++++ 1 file changed, 202 insertions(+) create mode 100644 chuan_harmonic_verlet_vr.py diff --git a/chuan_harmonic_verlet_vr.py b/chuan_harmonic_verlet_vr.py new file mode 100644 index 0000000..a472e38 --- /dev/null +++ b/chuan_harmonic_verlet_vr.py @@ -0,0 +1,202 @@ +#! /usr/bin/env python3 + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation +from matplotlib.patches import Circle, Rectangle +import math + +###################################### +# System parameters +###################################### +L = 4.0 +N = 22 +m = 1.0 + +###################################### +# Initial conditions: +# arrange particles as regural array +###################################### +# look for nside s.t.: nside*nside > N +x = np.empty(N, float) +y = np.empty(N, float) +R = np.empty(N, float) + +nside = int(math.ceil(math.sqrt(N))) +# grid spacing +dx = L/nside +for i in range(N): + rm = i % nside + rn = (i-rm)/nside + x[i] = (rn+0.5)*dx + y[i] = (rm+0.5)*dx + if i%2==0: + R[i] = 0.5 + else: + R[i] = 0.4 + +#px =np.zeros(N, float) +#py =np.zeros(N, float) +px =np.random.normal(0.0, 1.0, N) +py =np.random.normal(0.0, 1.0, N) +###################################### +# Difference vector between 2 atoms +###################################### +def difference_vector(i,j): + xij = x[j]-x[i] + yij = y[j]-y[i] + nx = math.floor(0.5-xij/L) + ny = math.floor(0.5-yij/L) + return xij + nx * L, yij + ny * L + +###################################### +# Interaction potential +###################################### +# parameters +k_part = 500 + +def pair_potential_energy(i,j): + (xij, yij) = difference_vector(i,j) + rij = math.sqrt(xij**2+yij**2) + D = R[i]+R[j] + if D > rij: + return 0.5*k_part*(rij-D)**2 + else: + return 0.0 + +def total_potential_energy(): + E = 0 + for i in range(N): + for j in range(i): + E += pair_potential_energy(i,j) + return E + +def virial(): + value = 0.0 + for i in range(N): + for j in range(i): + (xij, yij) = difference_vector(i,j) + rij = math.sqrt(xij**2+yij**2) + D = R[i]+R[j] + if D > rij: + value += k_part*(rij-D)*rij + return value + +def compute_forces(): + fx = np.zeros(N) + fy = np.zeros(N) + + for i in range(N): + for j in range(i): + (xij, yij) = difference_vector(i,j) + rij = math.sqrt(xij**2+yij**2) + D = R[i]+R[j] + if D > rij: + fij_x = k_part*(rij-D)*xij/rij + fij_y = k_part*(rij-D)*yij/rij + fx[i] += fij_x + fy[i] += fij_y + fx[j] -= fij_x + fy[j] -= fij_y + return fx, fy + + +###################################### +# Observables +###################################### +rmax = min(L/2,10) +Nbox = 100 +dr = rmax/Nbox +paircounter = np.zeros(Nbox, np.int32) +sample_nb = 0 + +def update_pair_correlation_counter(): + global sample_nb + for i in range(N): + for j in range(i): + (xij, yij) = difference_vector(i,j) + rij = math.sqrt(xij*xij+yij*yij) + n = int(math.floor(rij/dr)) + if n < Nbox: + paircounter[n] += 1 + sample_nb += 1 + +###################################### +# Verlet algorithm +###################################### +kT = 0.5 +dt = 0.01 +t = 0 +def verlet_step(): + global x, y, px, py, t + + x += 0.5*dt*px/m + y += 0.5*dt*py/m + + x = np.mod(x,L) + y = np.mod(y,L) + + (fx, fy) = compute_forces() + + px += dt*fx + py += dt*fy + + # velocity rescaling***** + if t > 100: + K = 0.5*np.sum(px*px+py*py)/m + lam = math.sqrt(N*kT/K) + px *= lam + py *= lam + #************************ + + + x += 0.5*dt*px/m + y += 0.5*dt*py/m + t += dt +''' +###################################### +# Drawing +###################################### +fig = plt.figure(figsize=(7,7)) +margin = 0.7 +ax = fig.add_axes([0, 0, 1, 1], frameon=False) +ax.set_xlim(-margin, L+margin), ax.set_xticks([]) +ax.set_ylim(-margin, L+margin), ax.set_yticks([]) +ax.add_artist(Rectangle((0,0),L,L,fill=False,edgecolor='black')) +circles = [ ax.add_artist(Circle((x[i],y[i]), R[i], edgecolor='none', facecolor='maroon',alpha=0.5)) for i in range(N) ] + +###################################### +# Image particles +###################################### +def update(framenb): + print(t) + verlet_step() + for i in range(N): + circles[i].center = x[i], y[i] +animation = FuncAnimation(fig, update, interval=200) +plt.show() + +''' + +###################################### +# Printout data +###################################### +dataE = open('dataE.dat','w') +dataP = open('dataP.dat','w') +dataK = open('dataK.dat','w') +tmax = 500 +nmax = int(tmax/dt) +for n in range(nmax): + + verlet_step() + + if (n%10)==0: + E = total_potential_energy() + K = 0.5*np.sum(px*px+py*py)/m + P = K/L/L-virial()/L/L/2 + print(n*dt, E, file=dataE) + dataE.flush() + print(n*dt, K/N, file=dataK) + dataK.flush() + print(n*dt, P, file=dataP) + dataP.flush() \ No newline at end of file -- GitLab