Thermal conductivity has proven to be one of the most difficult transport coefficients to calculate. Green-Kubo calculations are notoriously difficult to perform. Natural NEMD where one might simulate heat flow between walls maintained at different temperatures (Tenenbaum, Ciccotti & Gallico [1982]) is also fraught with major difficulties. Molecules stack against the walls leading to a major change in the microscopic fluid structure. This means that the results can be quite different from those characteristic of the bulk fluid. In order to measure a statistically significant heat flux, one must use enormously large temperature gradients. These gradients are so large that the absolute temperature of the system may change by 50% in a few tens of Ångstroms. The thermal conductivity that one obtains from such simulations is an average over the wide range of temperatures and densities present in the simulation cell.

We will now describe the most efficient presently known algorithm for calculating the thermal conductivity, (Evans, 1982b). This technique is synthetic, in that a fictitious field replaces the temperature gradient as the force driving the heat flux. Unlike real heat flow, this technique is homogeneous with no temperature or density gradients. We start with the Green-Kubo expression for the thermal conductivity (§4.4),

(6.58)

where J_Qz, is the z component of the heat flux vector. It appears to be impossible to construct a Hamiltonian algorithm for the calculation of thermal conductivity. This is because the equations of motion so obtained are discontinuous when used in conjunction with periodic boundary conditions. We shall instead invent an external field and its coupling to the phase of the N-particle system so that the heat flux generated by this external field is trivially related to the magnitude of the heat flux induced by a real temperature gradient.

Aided by the realisation that the heat flux vector is the diffusive energy flux, computed in a co-moving coordinate frame (see equation 3.151), we proposed the following equations of motion,

(6.59)

(6.60)

where E_i is the energy of particle i and,

(6.61)

the instantaneous average energy per particle.

There is no known Hamiltonian which generates these equations but they do satisfy AIΓ. This means that linear response theory can be applied in a straightforward fashion. The equations of motion are momentum preserving, homogeneous and compatible with the usual periodic boundary conditions. It is clear from the term (E_i - E^bar) F(t) that these equations of motion will drive a heat current. A particle whose energy is greater than the average energy will experience a force in the direction of F, while a particle whose energy is lower than the average will experience a force in the -F direction. Hotter particles are driven with the field; colder particles are driven against the field.

If the total momentum is zero it will be conserved and the dissipation is

(6.62)

Using linear response theory we have

(6.63)

Consider a field F=(0,0,F_z), then taking the limit t→∞ we find that the ratio of the induced heat flux to the product of the absolute temperature and the magnitude of the external field is in fact the thermal conductivity.

(6.64)

In the linear limit the effect the heat field has on the system is identical to that of a logarithmic temperature gradient (F = ∂lnT∕∂z). The theoretical justification for this algorithm is tied to linear response theory. No meaning is known for the finite field susceptibility.

In 1983 Gillan and Dixon introduced a slightly different synthetic method for computing the thermal conductivity (Gillan and Dixon, 1983). Although their algorithm is considerably more complex to apply in computer simulations, their equations of motion look quite similar to those given above. Gillan’s synthetic algorithm is of some theoretical interest since it is the only known algorithm which violates momentum conservation and AIΓ, (MacGowan and Evans, 1986b).

Figure 6.12 shows the thermal conductivity of the triple point Lennard-Jones fluid computed as a function of the strength of the heat field. We also show the experimental data for argon assuming that argon can be modelled by the standard Lennard-Jones model (ε∕k_B=119.8K, σ=3.405Å). The experimental uncertainties are so large that that if we used an accurate potential function, we could calculate the thermal conductivity more accurately than it can be measured.