Difference between revisions of "CompSciHW9"
From Predictive Chemistry
| Line 2: | Line 2: | ||
<pre> |
<pre> |
||
| − | 1) Write a complete code to simulate the Lennard-Jones gas in a 2D periodic box |
||
| + | 1) Use the van der Pol oscillator class to simulate 60 second (6000 steps) trajectories with 3 initial conditions, |
||
| − | with box length L = 14 and n=100 particles. Start them off on a 10x10 grid |
||
| + | x, v = [4, 3], [4, 4], and [4, 5] |
||
| − | with Gaussian distributed velocities. Ignore units and assume beta = m = 1. |
||
| + | Here, v = dot x = dx/dt |
||
| − | 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 u_ij^2 - u_ij ) |
||
| − | 2) Run the simulation for 100 steps, and create a plot showing the locations of the atoms |
||
| + | 2) Add a method to calculate the energy of a harmonic oscillator, |
||
| − | every 10 steps. |
||
| + | E(x, v) = (x^2 + v^2)/2 |
||
| − | 3) For every timestep, calculate the kinetic and potential energies. What do you observe |
||
| + | 3) Plot E(t) for all 3 simulations on the same axes. Remember to scale t = arange(6000)*dt. |
||
| − | 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 |
||
| + | 4) Is there a pattern that would let you know when the oscillator is near the steady-state? |
||
| − | for several different values of the numerical timestep, dt. |
||
</pre> |
</pre> |
||
| − | Hints: |
||
| + | Remember, we built the van der Pol integrator on [[CompSciWeek13]]. The complete code is below, |
||
| + | |||
<source lang="python"> |
<source lang="python"> |
||
| − | # Wrap all coordinates in an array to the range [0,L) |
||
| + | #!/usr/bin/env python |
||
| − | print x - L*floor(x/L) |
||
| + | |||
| + | import sys |
||
| + | from numpy import * |
||
| + | from numpy.random import standard_normal |
||
| + | norm = standard_normal |
||
| + | |||
| + | |||
| + | class State: |
||
| + | def __init__(self, x, v, mu, dt): |
||
| + | assert x.shape == v.shape |
||
| + | self.mu = mu |
||
| + | self.dt = dt |
||
| + | self.x = x |
||
| + | self.v = v |
||
| + | # |
||
| + | def step(self, N=1): |
||
| + | for i in xrange(N): |
||
| + | self.v += (self.mu*(1.0-self.x**2)*self.v - self.x)*self.dt |
||
| + | self.x += self.v*self.dt |
||
| + | def __str__(self): |
||
| + | s = "" |
||
| + | for x,v in zip(self.x, self.v): |
||
| + | s += " %f %f"%(x, v) |
||
| + | return s |
||
| − | # 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 |
||
| + | def main(argv): |
||
| − | r = closest_distance(z, x) # matrix of closest distances (n x 3) |
||
| + | assert len(argv) == 2, "Usage: %s <mu>"%argv[0] |
||
| − | r2 = sum(r*r, 1) # vector of squared distances (n) |
||
| + | mu = float(argv[1]) |
||
| − | u = r2**-3 |
||
| + | print mu |
||
| − | print 6 * sum(r * ((2*u*u - u)/sqrt(r2))[:,newaxis], 0) # sum over other atoms |
||
| + | dt = 0.01 |
||
| + | st = State(array([2.0]), array([0.0]), mu, dt) |
||
| + | st.step(10) |
||
| + | print st |
||
| + | |||
| + | if __name__=="__main__": |
||
| + | main(sys.argv) |
||
</source> |
</source> |
||
Latest revision as of 11:29, 13 April 2016
Introduction to Scientific Computing, HW 9. Due Friday, Apr. 15, 2016.
1) Use the van der Pol oscillator class to simulate 60 second (6000 steps) trajectories with 3 initial conditions, x, v = [4, 3], [4, 4], and [4, 5] Here, v = dot x = dx/dt 2) Add a method to calculate the energy of a harmonic oscillator, E(x, v) = (x^2 + v^2)/2 3) Plot E(t) for all 3 simulations on the same axes. Remember to scale t = arange(6000)*dt. 4) Is there a pattern that would let you know when the oscillator is near the steady-state?
Remember, we built the van der Pol integrator on CompSciWeek13. The complete code is below,
<source lang="python">
- !/usr/bin/env python
import sys from numpy import * from numpy.random import standard_normal norm = standard_normal
class State:
def __init__(self, x, v, mu, dt):
assert x.shape == v.shape
self.mu = mu
self.dt = dt
self.x = x
self.v = v
def step(self, N=1):
for i in xrange(N):
self.v += (self.mu*(1.0-self.x**2)*self.v - self.x)*self.dt
self.x += self.v*self.dt
def __str__(self):
s = ""
for x,v in zip(self.x, self.v):
s += " %f %f"%(x, v)
return s
def main(argv):
assert len(argv) == 2, "Usage: %s <mu>"%argv[0] mu = float(argv[1]) print mu dt = 0.01 st = State(array([2.0]), array([0.0]), mu, dt) st.step(10) print st
if __name__=="__main__":
main(sys.argv)
</source>
