<commit_before>---
layout: post
title: Molecular Dynamics Simulation
---

I recently wrote a discrete molecular dynamics simulation as a project for my simulations course. The Lennard-Jones potential is commonly used to model the interaction of between molecules.

<center><iframe src="//giphy.com/embed/KJes8CXWKg2JO" width="480" height="270" frameBorder="0" class="giphy-embed" allowFullScreen></iframe></center>
<center><i> Example simulation </i></center>

## Prerequisites
* Have Python 2.7 or greater installed
* A basic understanding of programming in Python
* Have [Pygame](http://pygame.org/hifi.html) installed: I used this particular [package](http://anaconda.org/krisvanneste/pygame)

## Lennard-Jones Potential
The Lennard-Jones Potential is defined as:

$$V_{LJ} = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right) ^{12} - \left( \frac{\sigma}{r} \right) ^{6} \right]  $$

where ε is the depth of the potential well, σ is the finite distance at which the inter-particle potential is zero, and r is the distance between the particles ([source](https://en.wikipedia.org/wiki/Lennard-Jones_potential)). 

<center><img src="https://upload.wikimedia.org/wikipedia/commons/5/51/12-6-Lennard-Jones-Potential.svg" width="50%" height="50%" alt="Lennard Jones Potential"></center>
<center><i>By <NAME>, Rainald62 [GFDL (http://www.gnu.org/copyleft/fdl.html), CC-BY-SA-3.0 (http://creativecommons.org/licenses/by-sa/3.0/) or CC BY-SA 2.5-2.0-1.0 (http://creativecommons.org/licenses/by-sa/2.5-2.0-1.0)], via Wikimedia Commons</i></center>

When we differentiate the potential energy with respect to r, we can find the Lennard-Jones Force:

$$F_{LJ} = - 24 \epsilon \left[ 2 \left( \frac{\sigma^{12}}{r^{13}} \right) - 6 \left( \frac{\sigma^{6}}{r^{7}}\right)\right]$$

This form will be more useful when calculating the inter-particle forces.

## Basic Strategy
I simulate n particles by storing the position, speed, and angle of direction for each particle. For every particle, I calculate the LJ force it experiences from every other particle. I sum all a particle's force vectors to determine the overall force and use that overall force vector calculate the new position, angle, and speed of each particle. Once that calculation is complete for all particles, the program renders the screen with the updated position of the particles. Since I jump forward in discrete steps, this is a discrete event simulation. Since I render each step, the simulation looks smooth and continuous as long as there aren't too many particles on screen for my laptop to handle. The upper limit on my laptop seems to be around 200 particles.

Each particle should have n force vectors. n-1 Lennard-Jones force vectors from the n-1 other particles and a vector representing it's current speed and direction.

## An Introduction to Pygame
Pygame is a FOSS Python module for writing video games. While I'm not exactly writing a game, Pygame is going to be really useful for rendering and visualizing our simulation. <NAME> wrote an excellent tutorial on simulating the trajectory of a cannonball in Pygame which can be found [here](http://www.petercollingridge.co.uk/pygame-physics-simulation). Rather than repeating information, I recommend going through parts 1 through 5 in his tutorial to create a basic Pygame window and a Particle class. 
 
## Helper functions
I'm going to need a couple of helper functions. First, I define a function to add two vectors (code from the beginning of part 6 of Peter Collingridge's tutorial). Next, I define a function which returns the angle between two particles. This is essential for figuring out which direction Lennard-Jones force points in. Of course, I'm also going to need a function which calculates the LJ force. A function which returns the scalar distance between two particles will be helpful too.

{% highlight python %}
def add_vector((angle1, length1), (angle2, length2)):
    """Return the resultant sum of two vectors
    Parameters:
        (angle1, length1): Angle and length of first vector
        (angle2, length2): Angle and length of second vector
    """
    x_component = math.sin(angle1) * length1 + math.sin(angle2) * length2
    y_component = math.cos(angle1) * length1 + math.cos(angle2) * length2
    length = math.hypot(x_component, y_component)
    angle = 0.5 * math.pi - math.atan2(y_component, x_component)
    return angle, length

def find_angle(a, b):
    """Find angle between two particles<commit_msg>Grammar fix<commit_after>---
layout: post
title: Molecular Dynamics Simulation
---

I recently wrote a discrete molecular dynamics simulation as a project for my simulations course. The Lennard-Jones potential is commonly used to model the interaction of between molecules.

<center><iframe src="//giphy.com/embed/KJes8CXWKg2JO" width="480" height="270" frameBorder="0" class="giphy-embed" allowFullScreen></iframe></center>
<center><i> Example simulation </i></center>

## Prerequisites
* Have Python 2.7 or greater installed
* A basic understanding of programming in Python
* Have [Pygame](http://pygame.org/hifi.html) installed: I used this particular [package](http://anaconda.org/krisvanneste/pygame)

## Lennard-Jones Potential
The Lennard-Jones Potential is defined as:

$$V_{LJ} = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right) ^{12} - \left( \frac{\sigma}{r} \right) ^{6} \right]  $$

where ε is the depth of the potential well, σ is the finite distance at which the inter-particle potential is zero, and r is the distance between the particles ([source](https://en.wikipedia.org/wiki/Lennard-Jones_potential)). 

<center><img src="https://upload.wikimedia.org/wikipedia/commons/5/51/12-6-Lennard-Jones-Potential.svg" width="50%" height="50%" alt="Lennard Jones Potential"></center>
<center><i>By <NAME>, Rainald62 [GFDL (http://www.gnu.org/copyleft/fdl.html), CC-BY-SA-3.0 (http://creativecommons.org/licenses/by-sa/3.0/) or CC BY-SA 2.5-2.0-1.0 (http://creativecommons.org/licenses/by-sa/2.5-2.0-1.0)], via Wikimedia Commons</i></center>

When we differentiate the potential energy with respect to r, we can find the Lennard-Jones Force:

$$F_{LJ} = - 24 \epsilon \left[ 2 \left( \frac{\sigma^{12}}{r^{13}} \right) - 6 \left( \frac{\sigma^{6}}{r^{7}}\right)\right]$$

This form will be more useful when calculating the inter-particle forces.

## Basic Strategy
I simulate n particles by storing the position, speed, and angle of direction for each particle. For every particle, I calculate the LJ force it experiences from every other particle. I sum all a particle's force vectors to determine the overall force and use that overall force vector calculate the new position, angle, and speed of each particle. Once that calculation is complete for all particles, the program renders the screen with the updated position of the particles. Since I jump forward in discrete steps, this is a discrete event simulation. Since I render each step, the simulation looks smooth and continuous as long as there aren't too many particles on screen for my laptop to handle. The upper limit on my laptop seems to be around 200 particles.

Each particle should have n force vectors. n-1 Lennard-Jones force vectors from the n-1 other particles and a vector representing its current speed and direction.

## An Introduction to Pygame
Pygame is a FOSS Python module for writing video games. While I'm not exactly writing a game, Pygame is going to be really useful for rendering and visualizing our simulation. <NAME> wrote an excellent tutorial on simulating the trajectory of a cannonball in Pygame which can be found [here](http://www.petercollingridge.co.uk/pygame-physics-simulation). Rather than repeating information, I recommend going through parts 1 through 5 in his tutorial to create a basic Pygame window and a Particle class. 
 
## Helper functions
I'm going to need a couple of helper functions. First, I define a function to add two vectors (code from the beginning of part 6 of <NAME>'s tutorial). Next, I define a function which returns the angle between two particles. This is essential for figuring out which direction Lennard-Jones force points in. Of course, I'm also going to need a function which calculates the LJ force. A function which returns the scalar distance between two particles will be helpful too.

{% highlight python %}
def add_vector((angle1, length1), (angle2, length2)):
    """Return the resultant sum of two vectors
    Parameters:
        (angle1, length1): Angle and length of first vector
        (angle2, length2): Angle and length of second vector
    """
    x_component = math.sin(angle1) * length1 + math.sin(angle2) * length2
    y_component = math.cos(angle1) * length1 + math.cos(angle2) * length2
    length = math.hypot(x_component, y_component)
    angle = 0.5 * math.pi - math.atan2(y_component, x_component)
    return angle, length

def find_angle(a, b):
    """Find angle between two particles
{'fraction_non_alphanumeric': 0.09700571941235842, 'fraction_numerical': 0.01794325445777728, 'mean_word_length': 4.427875836883749, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '365620', 'n_tokens_mistral': 2590, 'n_tokens_neox': 2413, 'n_words': 1111}