School of Physics

Modelling and Visualisation in Physics - FAQ

Which order do the steps in the Verlet algorithm come?

The short answer is "it doesn't matter". One of the attractions of the Verlet algorithm for molecular dynamics is that you only have to store three values for each particle (position velocity and acceleration). At each step (or half step) one of these is updated using the current values of the others.

If you write a Verlet method, it is best to start with the new force calculation: then the acceleration can be passed into the method, and the method can be used for all force models (no need to call LennardJonesPotential from inside VerletAlgorithm).

What are the units of the Lennard Jones potential?

Lennard Jones particles have two parameters in the potential, plus a mass. This is only enough to define the dimensions of length, mass and energy, so one simulation of Lennard-Jones is the same as any other - it can't be "fitted" to the behaviours of a particular material. You may find it useful to use "real" units, setting the minimum of the potential to be about 5eV and 2.5Angstroms, with the mass appropriate for an atom (a few tens of proton masses). If you use such units, and SI units for the Boltzmann constant, then the temperature scale will come out to be something sensible. However, you should be aware of numerical over/underflow errors if expressing 1 angstrom in metres and raising to the fourteenth power.

What should be changable in a molecular dynamics without recompiling?

Usually one would want to change the thermodynamic quantities (pressure, temperature, volume, number of particles) without recompiling. These things are typically determined as initial conditions, so do not affect code execution time (much). Changing the interatomic potential and changing the boundary conditions may be difficult to code without additional compute-time cost, so usually it is unnecessary.

How do I start the atoms?

Note that if you have symmetry in the initial positions, all the forces will have the same symmetry, and whatever vibrations occur, the symmetry will not break. Consequently, it is best to start without perfect symmetry.

On the other hand, arranging the particles on a lattice is a convenient way to avoid them overlapping too closely, thus a good trick is to put the particles onto a lattice initially, then assign each a random initial velocity and (small) displacement. The velocity you pick will determine the energy.

Also, with periodic boundaries, its a good idea to start with zero initial velocity, otherwise all your atoms will drift off in one direction. A good way to do this is to calculate the total momentum P after you've assigned random velocities, then subtract P/Nm from each of the N atoms of mass m.

In the NVE ensemble, how do I set the energy and temperature?

With difficulty.

Whatever your initial energy, it will be conserved throughout the run. But it will depend on the initial (random) positions. You can pick a target energy, calculate the potential energy from your initial state, then rescale all the random initial velocities to give the energy that you want.

Setting the temperature is even harder. Equipartition means that the energy is equally split between potential and kinetic, but there is no way to ensure this split at set up time.

There are other ensembles which allow one to control temperature. They all work by rescaling the velocity in some way or another.

How far apart should the atoms be?

The minimum of the potential is a little beyond the length parameter a. The equilibium 3D crystal structure is fcc, and the nearest neighbour separation is a little less than the minimum of the potential, because the more dintant neighbours of an atoms are also attracted to it.

Why isn't the energy conserved?

If total energy isn't conserved at each timestep, the calculation isn't right. Aside from coding errors, possible reasons are: timestep too large, neighbours not being found, boundary conditions not well implemented.

Do hard walls conserve energy?

They should, but need to be carefully implemented. The potential energy depends on position - if you move the particle back into the box, you will change the energy. The kinetic energy depends on the magnitude of the velocity - if you reverse the sign of one component that should be OK. But there is a caveat about timestep - the displacement at each step in MD should be "small relative to change of the potential". A hard wall is an infinite potential... Hard walls certainly do not conserve momentum. Indeed shock wave simulations are done by applying a hard wall boundary at one end of the simulation, and a free surface at the other. This is called a "momentum mirror".

How do I do the Fourier transform?

To get the fourier transform of the velocity autocorrelation function...

First, evaluate the ensemble averaged velocity autocorrelation function C(t)= < v(0).v(t) >. This is a single-valued function of t. You should store this in an array as you'll need it for each frequency.

Second, pick a frequency w and sum c(t) exp(iwt) for each timestep. from the end of the equilibration period to the end of the simulation

This gives a single complex number. The modulus of the complex number tell you the "amount" of frequency w in the signal c(t):

A good test is to use an example, say c(t)=A sin (at) + B cos (bt). The Fourier transform of this should give you delta functions at w=a and w=b (actually, the delta functions will have a finite width inversely proportional to the length of time the signal runs for.

Now, take different frequencies between the reciprocal of the timestep and the reciprocal of the whole length of the run. Actual vibration periods should be between these values.

For a small system, the only frequencies you get are those associated with phonons having a wavelength which is an integer fraction of the box size. Consequently the FT of the velocity autocorrelation function should be a few sharp peaks. The actual frequency can be approximated from that of a pair of atoms: the second derivative of the LJ function at the minimum. In units where u=4 and a=1, it comes out to be about 500 (check for yourself - being able to estimate orders of magnitude can be a crucial aid to debugging its the sort of thing you might get asked in an exam!). So you should look at freqencies between (say) 100 and 5000.

How do I set up a 2D crystal??

The square lattice is unstable in 2d, you need to set up a hexagonal lattice. If your box edges can be different

How do I profile my code?

If you run the code with "java -profile myCode" instead of "java myCode" then it will produce profiling statistics in a file called java.prof. Sun's documentation could be clearer, and there's a nice summary of optimisation methods here

Other environments provide other profilers. For example java -Xrunhprof produces more info in a file called java.hprof.txt see java -Xrunhprof:help .

Should the Autocorrelation be in a particle method?

It could be, see implementation from CompMeth

Why not use Runge-Kutta?

The Runge-Kutta algorithm gives a very accurate numerical integration of the trajectory. However, unlike Verlet, it is not time-reversible, and since time and energy are conjugate variables, time-reversibility gives excellent energy conservation which is more important in MD than accurate trajectories. Nth order Runge-Kutta also requires storage on more derivatives. Formally, Hamiltonian systems are best integrated using symplectic integrators (look it up!) which, even for finite timestep, give the exact integration of a Hamiltonian (though not the one you think you're integrating.

SIR Model

How do I choose which sites to update?

There are basically three options you might consider, one of which is wrong. The simplest (and the wrong one) is typewriter - loop over all sites sequentially. The problem with this is that it can break left-right symmetry- an infection can spread rapidly from left to right, but not from right to left. Synchronous updating involves a loop over all sites sequentially, updating each site based on the state of the whole system at time t, and copying the new state (t+1) into another temporary array. Once all sites are updated, the new array is pasted into the old one. When implementing this, be careful that your temporary array is instantiated as a separate object - the value of an object in java is its reference, so setting MyObject newObject = oldObject will not create a new object. The disadvantage with synchronous updating is that you need double the memory to store the state. Asynchronous updating involves picking a site at random to update. While simple, it has the additional overhead of a call to Math.random. In equilibrium Monte Carlo (when the update rules are based on energies from a Hamiltonian, such as the Ising model), Synchronous and Asynchronous updating give the same equilibrium state. In non-equilibrium systems (when the rules specify the model, without other constraints)

What are the initial conditions?

There is no clean answer to this.

Infections can be endemic i.e. persist peculiar to a district or particular locality, To study this type of infection you would need to have different environments: e.g. the bottom half of the grid have lower infection rate than the upper half.

epidemic i.e. spreading from an initial source to cover the whole region, and perhaps then dying out . To study this you use one initially infected site, and plot the number of infections against time. Assuming the infection ultimately dies out (e.g. if there is no recovery rate) then an epidemic is one where the number of infections is proportional to the size of the system. The final state of an epidemic has no remaining infected sites.

pandemic i.e. covers the whole region Here the system has a steady state with a finite number of infections. To study this it doesn't matter what starting conditions you use (except for a few extreme special cases, like all reisisant or none infected), the system will reach the same equilibrium.

What are the boundary conditions?

Although we don't attempt analytic solutions in this course, physicists like to do that and it is usually easiest with periodic boundary conditions (such that all sites are equivalent - otherwise the boundary sites obey different rules.). Real epidemiologists study real world geometries with boundaries like the UK coastline and replacing the square lattice with a network of connections. We'll look at networks in the final checkpoint.

How do I visualise a 4-dimensional space

In SIR, you may want to visualise the phase diagram of a system in terms of four variables, I, R, p & q. There are many ways of doing this and you should explore which gives the most easily understandable pictures.

You can break the variables into external and internal groups. In SIR p and q are the things you set, mean values of I and R are what you measure. e.g. you can easily use colour to represent the level of I on a 2D surface where x and y represent p and q, and another graph to represent R. Alternately, there are methods in swing which allow you to plot a 2D surface (there doesn't seem to be a standard yet, but if you dig around on the web there are several). Or you can export the data p q I R into a graphics package such as gnuplot.

What do I measure?

For endemic disease, we are interested in the steady state of the system. Regardless of the initial conditions, for a sufficiently large system after sufficiently long the mean numbers of S, I and R cells becomes steady. It is up to you to determine what "sufficiently large/long" means - it depends on the parameters, the update scheme and initial conditions. Experimentation is the only way to determine it. Note that for a finite system there is always a probability that the infection will die out by chance.

Beware of graphs which jump around a lot - this is often caused by finite size effects causing the infection to fail to get started.

For epidemic (high infectivity, long resistance), we are interested in how many sites get infected before the disease dies out. The length of time is then defined by the end of the epidemic. In a true epidemic, the number of infections is proportional to the size of the system. Conversely, if the disease dies out without becoming epidemic, the number of infections is a constant.

Can I see an example GUI with bells and whistles on?

Yes, here's a Mandelbrot set viewer and the source code There are some applets at the NANIA site. If you root around in "View Source" on the browser you may find the source code, and even trace classes back to previous MVP students! Try finding the name of the .class file, then accessing the equivalent .java file.

Is the java Math.random() any good?

No. Its neither particularly random nor particularly fast. A better routine is the Mersenne Twister. You can find java versions on the web - try downloading one and using the profiler to compare speed with Math.random().

Maxwells Equations

Do I need to do both relaxation and matrix methods?

No, you can do either. Gauss Siedel is easiest.

How big is that A matrix?

If your system is on a grid of length L, the A matrix is an L*L*L squared matrix. You really don't want to store it in your program...

What does the graphics mean?

If you plot field as a function of iterations, you get a nice looking diffusion towards the solution. This has no physical meaning, though it can help to see how the relaxation goes..

How do I invert the lower triangle for Gauss-Siedel?

You dont have to. The updating step is

u(n) = [b - Lu(n) - Uu(n-1)] /D

If you overwrite u(n) with its new values, the Lu(n) term automatically picks up terms of u which have already been updated.

What's wrong with periodic boundaries?

The potential is the integral of the field, so there is a constant of integration. In fact the physics is totally unchanged by adding a constant (or a constant slope) to the potential. Depending on how you do the minimisation, you might get lucky and hit a solution, or the constants might keep increasing for ever.

The only way to resolve this is to fix the value of the potential somewhere. Analytically the (implicit) choice is usually at infinity: this is obviously not an option in a program. Hence fixing the field to zero at the edges of a cube is one possibility (particle in a box) or at all eight corners where it should anyway be equal by symmetry is another (array of charges).

What's all this with the vector potential?

Cast in terms of the vector potential curl.B=curl.curl.A = J. From vector calculus we know that curl.curl.A = grad.div.A - lap.A = - lap.A since div.A is zero by definition. Thus the problem of solving for A is simply the Laplace equation three times over, one for each component of A. Once we have A everywhere, its straightforward to generate B=curl.A . and that's why the vector potential is worth knowing about.

Ising Model

How many neighbours should I use?

The Ising model is almost always run on a square lattice, with four neighbours. This gives the cleanest phase transition and has been solved analytically by Onsager. With four neighbours, the model is essentially identical if we change E to -E (the low T phase is just a check patterns). Conversely with six or eight neighbours, the order-disorder, (anti)ferromagnet-paramagnet transition is very different.

What about Temperature?

The Boltzmann probability of finding a macrostate with energy U is exp(-U/kT). The Ising model energy depends on the microstate via a single parameter, E. Thus the different behaviours of an Ising system depend on the number of like/unlike neighbours times a single parameter (E/kT). In the problem as posed, we have implicitly set kT=1, and will access different macrostates by changing E.

So am I changing energy or temperature?

Yes. (It amounts to the same thing)

Random updates or typewriter-to-copy?

Here it shouldn't make any difference. Since we're sampling an integral, there's no concept of "time" - unlike the SIR model the Ising simulation is strictly at equilibrium, so the answers really should be the same. Equilibrium is tied up with detailed balance (the probability of visiting a microstate at any time is proportional to its Boltzmann factor. You may have seen this idea in Statistical Mechanics.

What is the interesting range for E?

That's the question you have to resolve- but remember that it might be negative.

Can I update many sites at once ?

Yes. Two popular methods for this are Swendson-Wang Physical Review Letters 58 86 (1987) and Wolff (Physical Review Letters 62 361 (1989) In SW we create a bond between like neighbours with probability 1-exp(-2E/kT). A cluster is a set of points connected by bonds. Once the whole system is divided into clusters, we then set all sites in each cluster to a random value +-1 In Wolff we start with a random site, and add similar neighbour with probability 1-exp(-2E/kT). When we've built the cluster, flip it. This is slightly quicker than SW, as it is more likely to pick large clusters for flipping