We now present results of a direct numerical test the Kawasaki representation of the nonlinear isothermal response. We show that phase averages calculated using the explicitly normalised Kawasaki distribution function agree with those calculated directly from computer simulation.

The system we consider is the thermostatted NEMD simulation of planar Couette flow using the isothermal SLLOD algorithm (§6.3). As remarked ealier, the primary difficulty in using the Kawasaki expression in numerical calculations arises because it involves calculating an extensive exponential. For a 100-particle Lennard-Jones triple point system we would have to average quantities of the order of, e²⁰⁰, to determine the viscosity. Larger system sizes would involve proportionately larger exponents! The simulations presented here attempt to reduce these difficulties by using two strategies: they use a comparatively small number of particles, N=18 in two dimensions, and they were carried out at a low density, ρ^*=0.1, where the viscosity is ~50 times smaller than its triple point value. For small systems it is necessary to take into consideration terms of order, 1∕N, in the definition of the temperature,T = (Σ _i p_i ²∕m )∕(dN-d-1), and the shear stress, P_xyV=.(dN-d)∕(dN-d-1)Σ _i p_xi p_yi∕m )-(1/2)Σ _ijy_ijF_xij.

The first order equations of motion were solved using the 4th order Runge-Kutta method with a reduced timestep of 0.005. The reduced shear rate γ ^* = 1, and the reduced temperature was also set to unity.

The simulation consisted of a single equilibrium trajectory. At regular intervals (every 399 timesteps) the current configuration was used to construct four different configurations using the trajectory mappings described in §7.4. Each of these configurations was used as an initial starting point for a non-equilibrium simulation of 400 timesteps, with a negative timestep and reduced shear rate γ = 1. Time dependent averages were calculated, with the time being measured since the last equilibrium starting state. The aim was to produce the Kawasaki averages by exactly programming the dynamics in the Kawasaki distribution function (equation 7.24).

The phase space integral of the bare Kawasaki distribution function f(t), equation (7.24), is

(7.72)

Z(0) is the phase integral of the equilibrium distribution function which is equal to unity since f(0) is the normalised equilibrium distribution function. It is interesting to consider the rate of change of Z(t) after the external field is switched on. Using manipulations based on the reversible Liouville equation we can show that,

(7.73)

The last equality is a consequence of the Schrödinger-Heisenberg equivalence (§3.3). This implies that the bare Kawasaki distribution function is normalised for all times t. This is a direct result of the reversibility of the classical equations of motion. In Figure 7.10 we present the numerical results for Z(t). Figure 7.10 shows that equation (7.73) is clearly false. The normalisation is unity only for a time of the order of the Lyapunov time for the system. After this time the normalisation decreases rapidly. The explanation of this apparent paradox is that the analysis used to establish (7.73) is based on the reversible Liouville equation. The simulation used to generate the results shown in Figure 7.10 is, however, not time reversible. Trajectories which imply a long periods (compared to the Lyapunov time) of entropy decrease are mechanically unstable both in nature and in computer simulations. Because it is impossible to integrate the equations of motion exactly, these entropy decreasing trajectories are not observed for times longer than the Lyapunov time which characterises the irreversible instability of the equations of motion.

The form of the function, Z(t), shown in Figure 7.10, is determined by the accuracy with which the calculations are carried out. In principle by using ever more powerful computers one could, by increasing the word length and by decreasing the integration time step, ensure that the computed Z(t) stayed close to unity for longer and longer times. The exact result is that Z(t)=1. For a hard sphere system, the time over which the trajectory accuracy is better than a set tolerance only grows like, -ln(ε^1/λ) where λ is the largest Lyapunov exponent for the system and ε is the magnitude of the least significant digit representable on the computer. However our ability to numerically extend the times over which Z(t)~1, is much worse than this analysis implies. As we compute (7.72) for longer times, the variance in <exp[-βF_{e 0}∫^t ds J(-s)]> grows exponentially in time, regardless of the accuracy with which the trajectories are computed!

We have discussed the Kawasaki normalization in terms of numerical procedures. However exactly the same arguments apply to the experimental observation of the normalization. In nature, the problems in observing Z(t)~1 for long times result from uncontrolled external perturbations on the system rather than from numerical errors. However numerical error can be regarded as a particular form of external perturbation (ε above, would then be a measure of the background noise level). Of course the act of observation itself is a source of ‘external’ noise.

The results in Figure 7.10, show that the computed bare Kawasaki distribution function is not be properly normalised. Thus we should not surprised to see that the bare Kawasaki expression for the average shear stress is inconsistent with the results of direct calculation as is shown in Figure 7.11.

The obvious way around this problem is to explicitly normalise the distribution function (Morriss and Evans, 1987). The explicitly normalised form is

(7.74)

The renormalized average of the shear stress is then

(7.75)

We used computer simulation to compare the direct NEMD averages, and the bare and renormalized Kawasaki expressions for the time dependent average shear stress in a fluid. The results shown in Figure 7.11 are very encouraging. The renormalized Kawasaki result (denoted 'Kawasaki') agrees with that calculated directly and with the TTCF result. This is despite the fact that the normalisation has decreased by nearly two orders of magnitude at t^* = 2.0. The results show that the bare Kawasaki result is incorrect. It is two orders of magnitude smaller than the correct results.

Incidentally Figure 7.11 shows extraoridinarily close agreement (~0.2% for 0<t^*<2) between the TTCF prediction and direct NEMD. The agreement between direct NEMD and TTCF results for both the hydrostatic pressure and the normal stress difference is of a similar order. This indicates that one does not need to take the thermodynamic limit for the TTCF or GK formulae to be valid. Provided correct expressions are used for the temperature and the various thermodynamic fluxes, 18 particles seems sufficient.

Clearly no one should plan to use the renormalized Kawasaki formalism as a routine means of computing transport coefficients. It is probably the least efficient known method for computing nonequilibrium averages. The Kawasaki formalism is however, a very important theoretical tool. It was of crucial importance to the development of nonlinear response theory and it provides an extremely useful basis for subsequent theoretical derivations. As we will see in Chapter 9, the renormalized Kawasaki formalism, in contrast to the TTCF formalism, is very useful in providing an understanding of fluctuations in nonequilibrium steady states.