This chapter lists and describes the code employed in the calculation of non-linear Burnett coefficients. The program was initially designed for a canonical ensemble to be produced, using a standard Nosé-Hoover thermostat. Each member of this canonical ensemble was then used as a starting point for a computation at constant current, with both current and temperature fixed by a Gaussian feedback mechanism, as initially suggested by Evans and Lynden-Bell [1988]. Because we required the computation to accurately reflect the initial condition for each of the constant flux evolutions, we needed to use a self-starting differential equation solver, such as the fourth order Runge-Kutta method. Later on, we realised that a single equilibrium trajectory was all that was required.
Instead of modifying an existing molecular dynamics program, I decided to code from scratch a set of modules that can be mixed and matched to obtain the desired algorithm. The idea was to have a subroutine that integrates differential equations, and to pass to it a function that implements the equations of motion. So a typical subroutine call to the integrator will be (for example a Runge-Kutta algorithm)
CALL RK( X, F, N, H, NTSTEP )This subroutine integrates a first order differential equation of the form
. In the above subroutine call, SPM_quotX"
is a one dimensional array, SPM_quotF(I,X)" is a function returning the
ith component of 
, SPM_quotN" is the number of
components of SPM_quotX", SPM_quotH" is the stepsize, and SPM_quotNTSTEPS" is
the number of timesteps to be evolved. The result of the computation is
returned in SPM_quotX".
In molecular dynamics simulations, the SPM_quotX" to be integrated is just
the phase-space position 
. From the perspective of the
integrator, it does not matter in which order the components of 
occur, but practically it works best to organize the
x-components of the positions of all the particles, followed by the
y-components and z-components, then the three momentum components are
similarly arranged, as in fig 
.
Figure: Internal format of SPM_quotGAMMA"
This allows one to declare arrays containing only the x-components of
position for example, by means of SPM_quotEQUIVALENCE" statements, or common
blocks.