The first NEMD algorithm for self-diffusion was devised by Holian (Erpenbeck and Wood, 1977). In this elegant scheme the self-diffusion coefficient was evaluated as the limiting value of the mutual diffusion coefficient as the two species become identical. In this limit the two species differ only by a colour label which plays no role in their subsequent dynamics but which is reset in a probabilistic fashion as particles cross a labelling plane. A concentration gradient in coloured species is set up and the mutual diffusion coefficient is calculated from the constitutive relation (colour current/colour gradient). If the labels or colours of the atoms are ignored, the simulation is an ordinary equilibrium molecular dynamics simulation. If one calculates the species density as a function of position, the periodic boundary conditions imply that it is a periodic saw tooth profile. Exactly how sharp the teeth are, is not clear. The technique is inhomogeneous and is not applicable to mutual diffusion of species which are really different molecules. If the species are really distinct, the relabelling process will obviously generate discontinuities in pressure and energy.

The techniques we will describe are homogeneous. They do not create concentration gradients or coupled temperature gradients as does the Holian scheme. The algorithms can be extended to calculate mutual diffusion or thermal diffusion coefficients of actual mixtures (MacGowan and Evans, 1986a and Evans and MacGowan, 1987).

We begin by considering the Green-Kubo relation for the self diffusion coefficient (§4.1):

(6.8)

We design a Hamiltonian so that the susceptibility of the colour current to the magnitude of the perturbing colour field is closely related to the single-particle velocity autocorrelation function (6.8). Consider the colour Hamiltonian (Evans et. al., 1983)

(6.9)

where H₀ is the unperturbed Hamiltonian. The c_i are called colour charges. We call this property colour rather than charge to emphasise that H₀ is independent of the set of colour charges {c_i}. At equilibrium, in the absence of the colour field, the dynamics is colour blind. For simplicity we consider an even number of particles N, with

(6.10)

The response we consider is the colour current density J_x,

(6.11)

Since we are dealing with a Hamiltonian system, AIΓ (§5.3), is automatically satisfied. The dissipation function is

(6.12)

Linear response theory therefore predicts that (§5.1 & §5.3),

(6.13)

where the propagator implicit in J_x(t-s) is the field free equilibrium propagator. (Were we considering electrical rather than colour conductivity, equation (6.13) would give the Kubo expression for the electrical conductivity.) To obtain the diffusion coefficient we need to relate the colour current autocorrelation function to the single particle velocity autocorrelation function. This relation, as we shall see, depends slightly on the choice of the equilibrium ensemble. If we choose the canonical ensemble then

(6.14)

In the thermodynamic limit, for the canonical ensemble, if j≠i, then ⟨v_xi(t)v_xj(0)⟩=0,∀t. This is clear since if c is the sound speed, v_xj(0) can only be correlated with other particles within its sound cone (ie a volume with radius, ct). In the thermodynamic limit there will always be infinitely more particles outside the sound cone than within it. Since the particles outside this cone cannot possibly be correlated with particle i, we find that,

(6.15)

Combining this equation with the Green-Kubo relation for self diffusion gives,

(6.16)

If we are working within the molecular dynamics ensemble in which the total linear momentum of the system is zero, then v_xi is not independent of v_xj . In this case there is an order N^-1 correction to this equation and the self diffusion coefficient becomes (Evans et. al., 1983),

(6.17)

In the absence of a thermostat the order of the limits in (6.17) and (6.10) is important. They cannot be reversed. If a thermostat is applied to the system a trivial application of the results of §5.3 allows the limits to be taken in either order.

As an example of the use of thermostats we will now derive the Gaussian isokinetic version of the colour diffusion algorithm. Intuitively it is easy to see that as the heating effect is nonlinear (that is O(F²)), it does not effect the linear response. The equations of motion we employ are:

(6.18)

and

(6.19)

where the Gaussian multiplier required to thermostat the system is obtained from the constraint equation

(6.20)

In this definition of the temperature we calculate the peculiar particle velocities relative to the streaming velocity of each species. If one imagined that the two species are physically separated, then this definition of the temperature is independent of the bulk velocity of the two species. In the absence of this definition of the peculiar kinetic energy, the thermostat and the colour field would work against each other and the temperature would have an explicit quadratic dependence on the colour current. Combining (6.12 & 13) we identify the thermostatting multiplier as

(6.21)

In the original paper, (Evans, et.al., 1983), the thermostat was only applied to the components of the velocity which were orthogonal to the colour field. It can be shown that the linear response of these two systems is identical, provided the systems are at the same state point (in particular if the systems have the same temperature).

The algorithm is homogeneous since if we translate particle i and its interacting neighbours, the total force on i remains unchanged. The algorithm is also consistent with ordinary periodic boundary conditions (Figure 6.1). There is no change in the colour charge of particles if they enter or leave the primitive cell. It may seem paradoxical that we can measure diffusion coefficients without the presence of concentration gradients, however we have replaced the chemical potential gradient which drives real diffusion with a fictitious colour field. A gradient in chemical potential implies a composition gradient and a coupled temperature gradient. Our colour field acts homogeneously and leads to no temperature or density gradients. Linear response theory, when applied to our fictitious colour field, tells us how the transport properties of our fictitious mechanical system relate to the thermal transport process of diffusion.

By applying a sinusoidal colour field F(t)=F₀e^iωt, we can calculate the entire equilibrium velocity autocorrelation function. Noting the amplitude and the relative phase of the colour current we can calculate the complex frequency dependent susceptibility

(6.22)

An inverse Fourier-Laplace transform gives of χ(t) gives the velocity autocorrelation function.

Figure 6.2 shows the results of computer simulations of the diffusion coefficient for the 108 particle Lennard-Jones fluid at a reduced temperature of 1.08 and a reduced density of 0.85. The open circles were obtained using the algorithm outlined in this section (Evans et. al., 1983) which is based on equation (6.17). We see the colour conductivity (left y-axis) and the diffusion coefficient (right y-axis), plotted as a function of the colour current. The self diffusion coefficient is obtained by extrapolating the current to zero. The arrow denoted ‘EMD’, shows the results of equilibrium molecular dynamics where the diffusion coefficient was obtained (Levesque and Verlet, 1970), by integrating the velocity autocorrelation function (§4.1). The nonequilibrium and nonequilibrium simulations are in statistical agreement with each other.

Also shown in Figure 6.2, are the results of simulations performed at constant colour current, rather than constant applied colour field. We will return to this matter when we describe Norton ensemble methods in §6.6.

The filled in squares are the results of nonequilibrium simulations which were performed at constant colour current rather than constant applied colour field. The constant current methods will be described in more detail in §6.8. Briefly, one treats the colour field F(t) as a Lagrange multiplier whose value is chosen in such a way that the colour current is a constant of the motion. It is clear from the diagram that the constant current and constant colour field simulations are also in statistical agreement with each other.

In terms of computational efficiency, the self diffusion coefficient, being a single particle property, is far more efficiently computed from equilibrium simulations rather than from the algorithm given above. The algorithm we have outlined above is useful for pedagogical reasons. It is the simplest NEMD algorithm. It is also the basis for developing algorithms for the mutual diffusion coefficients of mixtures (Evans and MacGowan, 1987). The mutual diffusion coefficient, being a collective transport property, is difficult to calculate using equilibrium molecular dynamics (Erpenbeck, 1989). If the two coloured species are distinct electrically charged species, the colour conductivity is actually the electrical conductivity and the algorithm given above provides a simple means for its calculation.