Computational Methods Junior Honours : Checkpoint 5

This checkpoint investigates interacting 2 and 3 particle systems

Simulations with interacting particles are more complex than the previous two particle case in a constant field
in that the force on each particle now depends on the location of all the other particles.

This (expensive) force computation has been behind the design of many supercomputers, including the MD-GRAPE
line, and the IBM BlueGene systems. It is relevant in many areas of science from protein folding through galactic dynamics.

In this course you are limbering up for a simple case of molecular dynamics simulation. The first step will be to
calculate inter-particle forces correctly.

Before we begin, lets review some of the subtleties of syntax you have seen so far in the course:

Static & instance data members

Think about a class representing electrons, you would want the mass and charge to be static, but the position and velocity to be instance data members .

Static data members are shared among all instances of a class.

o class Electron { static double mass ; }     should be read as:  All electrons share a mass

Instance (normal) data members are private to each instance of a class.

o class Electron { Vector3d position; }      should be read as: Each electron has a position

Static & instance methods

Instance (normal) methods are functions that must be called in association with a specific instance of the class.
The reason for this is that the method is allowed to access the data members of the class (e.g. position) just like a local variable.
So, when the method is called, you have to tell it which instance it is to operate on.

An example of an instance method call is below. This is called through a specific instance of the class, namely A:

· Electron A;
A.OperateOnPosition();

When an instance method is called, the instance reference "this" is set up so that both "position" and "this.position" refer to A's position.

A static method does not operate on the instance members of this instance of a class.
It can still access the data members of an instance that is passed as a parameter to the method,
but it can't just access, for example, position or  this.position.

A static method can, however, access static data members, such as Electron.mass as these are properties of the class, not a specific instance.

An example of a static method call is below. This is called through the class, not through an instance.

· Electron.OperateOnMass();

In your arithmetic binary operations for your Vector3d class you used static methods which were passed specific instances
as parameters. Such as Vector3d.add(Vector3d a,Vector3d b). Add could access instance members a.x, and b.x because they
were passed, but it did not directly access just "x" or "this.x" because it is a static method.

Particle Force calculation:

We will write a class that will calculate the gravitational, electro-static, and Lennard-Jones forces between two particles.
This will have independent constants of proportionality, E, G, a, and U to allow either force to be switched off by setting these to zero.

See the background material for more detail.

The ParticeForce class represents fundamental forces that are constant for a simulation. Because of this, it doesn't really make sense to have
two seperate instance of the class that can have two seperate values for the constants. Thus you should make all methods and all
data members of this class static.
This is an example of  a purely static class. This is useful to hold/parcel up constants and code, but doesn't make much conceptual sense to use as a multiple instance object.

•    Write a class called ParticleForce in ParticleForce.java. It should contain
• static double members G, E, a, U + getters and setters
• static Vector3d eForce(Particle3d a,Particle3d b) computes the electrostatic force on a due to b
• static Vector3d gForce(Particle3d a,Particle3d b) computes the gravitational force on a due to b
• static Vector3d ljForce(Particle3d a,Particle3d b) computes the Lennard-Jones force on a due to b
• static Vector3d Force(Particle3d a,Particle3d b) computes the sum of all forces on a due to b

Two Body Systems:

• Write a main programme Particle3dTwoBodyForce.java to read from file the initial state for two particles, integration limits and constants for the three potentials
• This programme should call the ParticleForce.Force method so that it will work for ALL of the force cases by setting them to zero/one

· Advanced : Command line arguments

o

o Make the files opened be specified via command line arguments. It is wise to do this as it will make it easier to manage your datafiles later.

§ Command line arguments are passed in String [] argv to main

§ If you run your programme as 

§ java Particle3dTwoBodyForce file1 file2 file3

then argv[0] will be "file1", etc...

Think about how most conveniently to group the input information into different files.

One of the file name arguments should be the output filename. This will make it easy to keep different output datafiless around.

Check the number of arguments is right with argv.length and an if statement

Two body systems tasks:

· Simulate and save trajectory plots for each particle in a gravitational system

• G=1
• Sun, mass 1, at origin, velocity adjusted to give zero total momentum
• Planet, mass 0.01, pos = (1,0,0), vel = (0,1,0)
• Verify Keplers third law holds (square of T proportional to cube of R).

· Simulate and save plots of a classical Hydrogen atom

• E=1, Proton mass = 1000, Electron mass = 1, q's=1,-1
• Proton at origin, adjust velocity for momentum zero
• Electron at (1,0,0) velocity (0,1,0).
• Zoom in on the proton orbit

· Advanced : visualization

o Using a parameter in your input file to select one of two possible output file formats

§ Standard xmgrace compatible output

§ VMD xyz output format

§ VMD (visualise molecular dynamics) is a package for displaying three dimensional systems of atoms and molecules.

§ We can use it to visualize our particles.

§ With two particles, the output file should contain

2

BLANK LINE

atom1 x y z

atom2 x y z

2

BLANK LINE

atom1 x y z

atom2 x y z

…. And so on. Three, or many, particles is the obvious extension of this.

· Viewing with VMD

o (assuming your file is called output.dat) View the file using the command vmd –xyz output.dat

Three Body Systems:

· New Java constructs:

o Typing all the n (n-1) = 6 ordered particle pairings to compute the force is clearly going to get tedious. Ponder the cost of this for an entire protein.

o Instead we will to add a new (general) method to ParticleForce. but it is good practice to maintain the old version of the code unaltered too.

o Java allows you enhance old classes with new functionality without either breaking or duplicating the old code. For example:

§ public class Newer extends Older { public static void newMethod(double argument) { does_something_new();} };

• This is a feature called inheritance. Why is this useful?
• Polymorphism allows an instance of class Newer to pretend to be an instance of class Older when interacting with (possibly old) code that treats Older
• Why might polymorphism also be useful if you have, say, four different derived classes Newer1, Newer2, Newer3, Newer4.
• In this case Older would be called the base class for all four derived classes

· Write a new file ParticlesForce.java which has a class ParticlesForce that extends the ParticleForce class implementing the following

• public class ParticlesForce extends ParticleForce
• public static void leapVelocity(Particle3d [],int nparticle, double dt);
• public static void leapPosition(Particle3d [],int nparticle, double dt);
• NO CODE FROM ParticleForce SHOULD BE CUT/PASTED IN - THAT IS WHAT THE EXTENDS IS FOR!!!
• You should use for loops over the n particles to accumulate the force.
• You can use a loop within a loop
•  should avoid attempting to calculate the force of a particle on itself using an if statement.

· Write a new main programme Particle3dThreeBody.java to integrate a three particle system.

o It will be simpler to use an array of Particles declared as

§ Particle3d ps[] = new Particle3d[3];

· For each of the systems below play with the starting conditions to investigate stability - does one

particle escape at large time. Save plots of each. Try changing the velocities and masses.

• Simulate a gravitational system

• G=1
• Sun, mass 1, at origin, velocity adjusted to give zero total momentum
• Planet1, mass 0.01, pos = (1,0,0), vel = (0,1,0)
• Planet2, mass 0.001, pos = (2,0,0), vel = (0,v,0)
• Try v= 0.5, 1, 1.5. What value gives a circular orbit (this will be a different velocity - does this agree with the equations of motion)?

· Simulate a classical Helium atom

• E=1, Nucleus mass = 4000, Electron mass = 1, q's=2,-1,-1
• Proton at origin, adjust velocity for momentum zero
• Electron1 at (1,0,0) velocity (0,sqrt(2),0).
• Electron2 at (-1,0,0) velocity (0,-sqrt(2),0).
• This is stable. Now try to find another stable pair.

· Simulate a three molecule Lennard-Jones system.

· Take U = 4 and a=1, note that r is the vector between the particles.

· Make sure you understand why the powers for the "12-6" potential are 14 and 8, and ensure that the force acts in the correct direction.

· Note that the first term is repulsive, and acting alone pushes the atoms apart, while the second term is cohesive and pulls the atoms together.

• atom 1, pos = (0.7,0,0), vel = 0
• atom 2, pos = (-0.7,0,0), vel = 0
• atom 3, pos = (0,1,0), vel = 0
• What can you say about the position of atom3