Difference between revisions of "CompSciHW9"

From Predictive Chemistry
Jump to: navigation, search
 
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 12: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">

  1. !/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>