Commit 234a0771 authored by Trung Hieu NGUYEN's avatar Trung Hieu NGUYEN

new code

parent aad9dac1
#! /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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment