Dynamics

So far all we have done is displayed the positions of the balls after minimisation calculations. To see the atoms move in real time subject to the Lennard-Jones potential, we need to use the force calculated to work out the velocity of each atom, then move the atoms accordingly. This is placed in a loop, and so the atoms interact with each other in real time.

To start off, we need to initialise the velocity. For now, we'll initialise it to 0.

<Initialise Velocity>=
veloc = []
for atom in atoms:
    veloc.append((0,0,0))

Used below.

Now we need a routine to move the atoms. We use a variable tstep to represent the amount of time that occurs between each step. Of course, the greater the value of tstep, the less accurate the simulation will be.

<Update Coordinates>=
for i in range(len(atoms)):
    atoms[i] = (atoms[i][0] + tstep * (veloc[i][0] +
                    (tstep * newForce[i][0] / 2.0)),
                atoms[i][1] + tstep * (veloc[i][1] +
                    (tstep * newForce[i][1] / 2.0)),
                atoms[i][2] + tstep * (veloc[i][2] +
                    (tstep * newForce[i][2] / 2.0)))

Used below.

The velocity needs to be updated in each step, this is done using the following routine.

<Update Velocity>=
for i in range(len(veloc)):
    veloc[i] = (veloc[i][0] + tstep *
                    (newForce[i][0] + oldForce[i][0]) / 2.0,
                veloc[i][1] + tstep * 
                    (newForce[i][1] + oldForce[i][1]) / 2.0,
                veloc[i][2] + tstep * 
                    (newForce[i][2] + oldForce[i][2]) / 2.0)

Used below.

Now we'll package the main loop and include our visualisation routines. We use the variable nstep to represent the number of steps the simulation will do before exiting.

<MD Initialisation and Main Loop>=
atoms = createSimpleCube(rside, nside)
<Initialise Velocity>
newForce = oldForce = evalForce(atoms)
<Create Balls>

for istep in range(nstep):
    <Update Coordinates>
    newForce = evalForce(atoms)
    <Update Velocity>
    oldForce = newForce
    <Move Balls>

Used below.

*

Putting it all together

Finally, we need to put it all together, and add a simple command line interface to setting the values of tstep and nstep for the simulator, and rside and nside for both the simulator and the minimiser. For now, to make the command simpler, we'll keep the values in the tuples of rside and nside the same, hence we will only be working with perfect cubes.

<Minimise.py>=
from visual import *
import sys
import math

<Create Simple Cube>
<Evaluate Force>
<Evaluate Potential>
<Next Atoms>
<Golden Section Search>

try:
    r = float(sys.argv[1])
    rside = (r, r, r)
    n = int(sys.argv[2])
    nside = (n, n, n)
except (ValueError, IndexError):
    print "Usage: python %s rside nside" % sys.argv[0]
    sys.exit(1)

scene.autoscale = 0
scene.title = "Molecular dynamics Minimiser"

<Minimise Initialisation and Loop>
print "The minimum potential is %e" % totalPotential

This code is written to a file (or else not used).

*

<SimpleMD.py>=
from visual import *
import sys

<Create Simple Cube>
<Evaluate Force>

try:
    tstep = float(sys.argv[3])
    nstep = int(sys.argv[4])
    r = float(sys.argv[1])
    rside = (r, r, r)
    n = int(sys.argv[2])
    nside = (n, n, n)
except (ValueError, IndexError):
    print "Usage: python %s rside nside tstep nstep" % sys.argv[0]
    sys.exit(1)

scene.autoscale = 0
scene.title = "Molecular dynamics simulator"

<MD Initialisation and Main Loop>

This code is written to a file (or else not used).

• <Create Balls>: U1, D2, U3
• <Create Simple Cube>: D1, U2, U3
• <Evaluate Force>: D1, U2, U3
• <Evaluate Potential>: D1, U2
• <Golden Section Search>: D1, U2
• <Initialise Velocity>: D1, U2
• <MD Initialisation and Main Loop>: D1, U2
• <Minimise Initialisation and Loop>: D1, U2
• <Minimise.py>: D1
• <Move Balls>: U1, D2, U3
• <Next Atoms>: D1, U2
• <SimpleMD.py>: D1
• <Update Coordinates>: D1, U2
• <Update Velocity>: D1, U2

*