Latest revision as of 13:51, 6 March 2017

Two-Body Problems

This page explains the numerical simulation of two-body problems suitable for a first numerical simulation.

The coordinates of two point particles can be represented by two vectors, <math>\vec x_1 \in \mathbb R^3</math> and <math>\vec x_2 \in \mathbb R^3</math>. They have corresponding masses and velocities (<math>m_i,\vec v_i</math>).

Because the interaction between the particles is invariant to translations, the energy of interaction can only depend on the difference between coordinates, <math>\vec r = \vec x_2 - \vec x_1</math>. Even better, since the interaction should be invariant to rotation, the potential energy should really only depend on the length of <math>\vec r</math>, which we denote <math>r = |\vec r|</math>. You must memorize the formula for computing the magnitude of a vector,

<math>|\vec r|^2 = \sum_i \vec r_i^2 = r_x^2 + r_y^2 + r_z^2</math>.

Common Energy Functions

Harmonic Oscillator

The most common energy function is undoubtedly the potential energy of the harmonic oscillator,

<math>U(r) = \frac{k}{2} (r - r_0)^2</math>.

It has two parameters, k, and <math>r_0</math>. Because of its commonly usage to represents the energy of a bond between atoms, <math>r_0</math> is sometimes referred to as the bond length, and k as the bond stretch constant.

Lennard-Jones

The second most common pair-energy function is arguably the Lennard-Jones form,

<math>U(r) = 4 \epsilon\left( \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right)</math>.

This form is used to represent a pair of atoms that are not bonded. It also has two parameters. <math>\epsilon</math> is the `well-depth' and <math>\sigma</math> is the distance where the energy crosses from positive to negative. There is an alternate way to write this energy by re-defining <math>R = \sigma 2^{1/6}</math> and re-factoring the 4 inside the parentheses.

Other Forms

There are many other forms that try to represent the actual physical interaction energy between pairs of atoms more realistically than the primitive approximations above. You should look these up on your own

• Wikipedia:Morse potential
• Wikipedia:Buckingham potential
• Wikipedia:Coulomb potential
• Wikipedia:Gravitational potential
• FENE (finitely extensible nonlinear elastic) Potential

For the FENE potential (and the others as well) an excellent reference is the Gromacs User Manual and high-level summaries of these functions are also present in the LAMMPS Documentation.

Derivatives of the Energy

According to Newton's law, we must compute the force as,

<math>\vec F = -\vec \nabla U</math>, where <math>U</math> is the potential energy function.

In component form, <math>\vec \nabla_i = \frac{\partial}{\partial x_i}</math> is the derivative with respect to the <math> i^\mathrm{th}</math> component.

To find the derivative of <math>U(r)</math> with respect to <math>\vec x_i</math>, the chain rule has to be used at least twice. The first one is easy, and can be summarized in vector form,

<math>\vec \nabla_{x_2} = -\vec \nabla_{x_2} = \nabla_r</math>

In English, <math>\vec x_2</math> increases as much as <math>\vec r</math> does, and <math>x_1</math> decreases as much as <math>\vec r</math> increases.

The second part is harder, since it requires going from <math>\vec r</math> dependence to <math>\vec r</math> dependence.

It is written in vector form as,

<math> \vec \nabla_r U(r) = \frac{\vec r}{r} \frac{dU(r)}{dr}</math>

Example Implemenentation

This example shows the steps needed to run a 2-particle simulation in outline form. It is not complete, since the functions are not yet implemented.

<source lang="python"> from numpy import * import pylab as plt

1. Understand system.

def plot_atoms(x):

   # Make a plot of all atoms.
   #plt.plot([x0[0], x1[0]], [x0[1], x1[1]], 'o')
   plt.plot([x[0,0]], [x[0,1]], 'o')
   plt.plot([x[1,0]], [x[1,1]], 'o')
   # cruft
   plt.gca().set_xlim([0,10])
   plt.gca().set_ylim([0,10])

   plt.show()

def mag(r):

   return sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])

1. Use physics.

def calc_energy(x):

   # coords -> distance vectors
   r = x[1] - x[0]

   # distance vectors -> lengths
   print(mag(r))
   print(r)
   # lengths -> energies

def main():

   # Define starting locations
   # x has 2 indices:
   #  index 0 - atom number
   #  index 1 - x,y, or z
   x = array([ [5,5,-1],
               [7,8,0] ])
   plot_atoms(x)
   calc_energy(x)

main() </source>