Difference between revisions of "CompSciHW10"
From Predictive Chemistry
m |
|||
| Line 34: | Line 34: | ||
u = r2**-3 |
u = r2**-3 |
||
print 6 * sum(r * ((2*u*u - u)/r2)[:,newaxis], 0) # sum over other atoms |
print 6 * sum(r * ((2*u*u - u)/r2)[:,newaxis], 0) # sum over other atoms |
||
| + | </source> |
||
| + | |||
| + | A code that allows you to visualize your (steps x atoms x 2) trajectory is shown below: |
||
| + | |||
| + | <source lang="python"> |
||
| + | #!/usr/bin/env python |
||
| + | from numpy import * |
||
| + | import matplotlib.animation as anim |
||
| + | import pylab as plt |
||
| + | |||
| + | dt = 0.01 |
||
| + | |||
| + | def show_trj(x, dt): |
||
| + | fig = plt.figure() |
||
| + | ax = fig.add_subplot(111, autoscale_on=False, |
||
| + | xlim=(0.0, 10.0), ylim=(0.0, 10.0)) |
||
| + | ax.set_aspect('equal') |
||
| + | ax.grid() |
||
| + | line, = ax.plot([], [], 'o') |
||
| + | template = "time %.1f" |
||
| + | time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes) |
||
| + | |||
| + | def animate(i): |
||
| + | line.set_data(x[i,:,0], x[i,:,1]) |
||
| + | time_text.set_text(template % (i*dt)) |
||
| + | return line, time_text |
||
| + | |||
| + | def init(): |
||
| + | line.set_data([], []) |
||
| + | time_text.set_text('') |
||
| + | return line, time_text |
||
| + | |||
| + | ani = anim.FuncAnimation(fig, animate, \ |
||
| + | arange(len(x)), interval=25, blit=True, \ |
||
| + | init_func=init) |
||
| + | plt.show() |
||
| + | |||
| + | # Replace this with your actual dynamics code! |
||
| + | trj = random.random((100,10,2))*10.0 |
||
| + | |||
| + | show_trj(trj, dt) |
||
</source> |
</source> |
||
Latest revision as of 12:23, 18 April 2016
Intro. Scientific Computing, HW10 - Due Friday, April 22.
1) Write a complete code to simulate the Lennard-Jones gas in a 2D periodic box
with box length L = 14 and n=100 particles. Start them off on a 10x10 grid
with Gaussian distributed velocities. Ignore units and assume beta = m = 1.
The Hamiltonian is given by H = sum_j m v_j^2/2 + 1/2 sum_{i != j} u_ij^2 - u_ij
where u_ij = |x_i - x_j|**-6
The force on each particle, i, is therefore
F_i = sum_{j != i} 6 (x_i - x_j) / |x_i - x_j|^2 ( 2 u_ij^2 - u_ij )
2) Run the simulation for 100 steps, and create a plot showing the locations of the atoms
every 10 steps.
3) For every timestep, calculate the kinetic and potential energies. What do you observe
about the behavior of the potential energy?
4) Make a plot of the total energy vs. time for your 100 step simulation. Overlay these plots
for several different values of the numerical timestep, dt.
Hints: <source lang="python">
- Wrap all coordinates in an array to the range [0,L)
print x - L*floor(x/L)
- Find the closest distance between two points, r_ij
y = x[i] - x[j] print y - L*floor(y/L + 0.5)
- Calculate the LJ force on an atom at point z
r = closest_distance(z, x) # matrix of closest distances (n x 3) r2 = sum(r*r, 1) # vector of squared distances (n) u = r2**-3 print 6 * sum(r * ((2*u*u - u)/r2)[:,newaxis], 0) # sum over other atoms </source>
A code that allows you to visualize your (steps x atoms x 2) trajectory is shown below:
<source lang="python">
- !/usr/bin/env python
from numpy import * import matplotlib.animation as anim import pylab as plt
dt = 0.01
def show_trj(x, dt):
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False,
xlim=(0.0, 10.0), ylim=(0.0, 10.0))
ax.set_aspect('equal')
ax.grid()
line, = ax.plot([], [], 'o')
template = "time %.1f"
time_text = ax.text(0.05, 0.9, , transform=ax.transAxes)
def animate(i):
line.set_data(x[i,:,0], x[i,:,1])
time_text.set_text(template % (i*dt))
return line, time_text
def init():
line.set_data([], [])
time_text.set_text()
return line, time_text
ani = anim.FuncAnimation(fig, animate, \
arange(len(x)), interval=25, blit=True, \
init_func=init)
plt.show()
- Replace this with your actual dynamics code!
trj = random.random((100,10,2))*10.0
show_trj(trj, dt) </source>
