from IPython.display import clear_output
import numpy as np
from MD_utils import lennardjones,simple_molecule_vis
from scipy.constants import N_A
epsilon=0.9959
sigma=0.3405
k_B = 0.008311
M = 39.948
plt = simple_molecule_vis(bs=1)
import numpy as np
# bs jest rozmiarem kostki w której znajdują się cząstki.
nx = 8
dx = sigma *1.2
bs = dx*nx
box = np.array([bs,bs,bs])
try:
plt.update_box(bs=bs)
print "no plot object"
except:
pass
Nparticles = nx**3
print Nparticles
U = np.zeros((3,Nparticles))
l = 0
for i in range(nx):
for j in range(nx):
for k in range(nx):
U[:,l] = (dx*i-dx*(nx/2-0.50),dx*j-dx*(nx/2-0.5),dx*k-dx*(nx/2-0.5))
l+=1
V = np.zeros_like(U)
plt.pkts.points_colors = 0xff0000*np.ones(Nparticles)
plt.pkts.points_positions = U.T.copy()
plt.pkts.point_size = 3
no plot object 512
box = np.array([bs,bs,bs])
%time Epot, F, Vir = lennardjones(U, box,sigma = 0.3405, epsilon=0.9959)
CPU times: user 4.55 ms, sys: 112 µs, total: 4.67 ms Wall time: 4.58 ms
V = 0.1*(np.random.randn(3,Nparticles)-0.5)
plt.plot.display()
%%time
# Start simulation
import time
n_steps = 1520
dt = 0.005
N = Nparticles
box = np.array([bs,bs,bs])
plt.update_box(bs=bs)
(epot,F,Vir) = lennardjones(U,box,sigma = 0.3405, epsilon=0.9959)
traj = []
for i in range(n_steps):
U += V * dt + 0.5 * F/M * dt * dt
F0 = F[:]
(epot, F, Vir) = lennardjones(U, box,sigma = 0.3405, epsilon=0.9959)
V += 0.5 * (F + F0)/M * dt
U -= bs*np.rint(U/bs)
traj.append(U[:,233].copy())
if i%10==0:
#time.sleep(0.1)
T = M*np.sum(V**2)/(k_B*(3*N-6))
P = 1/bs**3*( N*k_B*T + 1/3.*Vir )
Ek = 0.5*M*np.sum(V**2)
plt.pkts.points_positions = U.T.copy()
plt.pkts.points_colors = (np.sum((V**2),axis=0)/10.0* 0xffffff).astype(np.int64)
plt.update_box(bs=bs)
clear_output(wait=True)
print i,epot,Ek,epot+Ek,"T=",T,"P=",P
print "Vir:",Vir
print "bs:",bs
traj = np.array(traj)
1510 -2352.57522965 585.612548004 -1766.96268165 T= 92.1076403198 P= -199.355972869 Vir: -22064.7173548 bs: 3.2688 CPU times: user 6.55 s, sys: 270 ms, total: 6.82 s Wall time: 6.65 s
from k3d import K3D
plt_traj = K3D.points(traj.copy(),color=0xff0000, point_size=1.2)
plt.plot += plt_traj
plt_traj.point_size = 3