We will now give another example of the usefulness of the Norton ensemble. Suppose we wish to calculate the yield stress of a Bingham plastic - a solid with a yield stress. If we use the SLLOD method outlined above the Bingham plastic will always yield simply because the strain rate is an input into the simulation. It would not be easy to determine the yield stress from such a calculation. For simulating yield phenomena one would prefer the shear stress as the input variable. If this were the case simulations could be run for a series of incremented values of the shear stress. If the stress was less than the yield stress, the solid would strain elastically under the stress. Once the yield stress was exceeded, the material would shear.

Here we discuss a simple method for performing NEMD simulations in the stress ensemble. We will use this as an opportunity to illustrate the use the Nosé-Hoover feedback mechanism. We will also derive linear response expressions for the viscosity within the context of the Norton ensemble. The equations of motion for shear flow, thermostatted using the Nosé-Hoover thermostat are

(6.98)

(6.99)

(6.100)

Using the Nosé-Hoover feedback mechanism we relate the rate of change of the strain rate, γ, to the degree to which the instantaneous shear stress, -P_xy(Γ) differs from a specified mean value, -S_xy(t). We therefore determine the strain rate from the differential equation,

(6.101)

If the instantaneous stress is greater (ie more negative) than the specified value, the strain rate will decrease in an attempt to make the two stresses more nearly equal. The relaxation constant Q_γ should be chosen so that the timescale for feedback fluctuations is roughly equal to the natural relaxation time of the system.

From the equations of motion, the time derivative of the internal energy H₀=∑_ip_i²∕2m+Φ, is easily seen to be,

(6.102)

The Nosé constant stress, constant temperature dynamics satisfy a Liouville equation in which phase space behaves as a compressible 6N+2 dimensional fluid. The equilibrium distribution function is a function of the 3N particle coordinates, the 3N particle momenta, the thermostatting multiplier ξ, and strain rate γ, f₀=f₀(Γ,ξ,γ). The Liouville equation for this system is then

(6.103)

Since

and

then

(6.104)

the phase space compression factor Λ(Γ) is easily seen to be -3Nξ. If we consider the time derivative of the extended internal energy H₀+½Q_ξξ²+Q_γγ² we find that

(6.105)

If we consider the situation at equilibrium when the set value of the shear stress, -S_xy(t), is zero and K₀ = 3N∕2β, the Liouville equation becomes

(6.106)

Integrating both sides with respect to time gives the equilibrium distribution function for the constant stress Norton ensemble to be

(6.107)

The equilibrium distribution function is thus a generalised canonical distribution, permitting strain rate fluctuations. Indeed the mean square strain rate is

(6.108)

so the amplitude of the strain rate fluctuations are controlled by the adjustable constant Q_γ.

We wish to calculate the linear response of an equilibrium ensemble of systems (characterised by the distribution f₀, at time t=0), to an externally imposed time dependent shear stress, -S_xy(t). For the Nosé-Hoover feedback mechanism the external field is the mean shear stress, and it appears explicitly in the equations of motion (Hood, Evans and Morriss, 1987). This is in contrast to the more difficult Gaussian case (Brown and Clarke, 1986). For the Gaussian feedback mechanism the numerical value of the constraint variable does not usually appear explicitly in the equations of motion. This is a natural consequence of the differential nature of the Gaussian feedback scheme.

The linear response of an arbitrary phase variable B(Γ) to an applied time dependent external field is given by

(6.109)

where iL₀ is the equilibrium (Nosé-Hoover thermostatted) f-Liouvillean and iΔL(s) = iL(s)-iL₀ where iL(s) is the full field dependent thermostatted f-Liouvillean. It only remains to calculate iΔL(s) f₀. Using the equations of motion and the equilibrium distribution function obtained previously we see that,

(6.110)

Here we make explicit reference to the phase dependence of γ, and the explicit time dependence of the external field S_xy(t). The quantity -S_xy(t)Vγ(Γ) is the adiabatic derivative of the extended internal energy, E=H₀+½Q_γγ².

Combining these results the linear response of the phase variable B is

(6.111)

In order to compute the shear viscosity of the system we need to calculate the time dependence of the thermodynamic force and flux which appear in the defining constitutive relation for shear viscosity. Because of the presence of the Nosé-Hoover relaxation time, controlled by the parameter Q_γ, the actual shear stress in the system -P_xy(Γ), does not match the externally imposed shear stress S_xy(t), instantaneously. To compute the shear viscosity we need to know the precise relation between P_xy and γ, not that between S_xy and the strain rate. The two quantities of interest are easily computed from (6.111).

(6.112)

(6.113)

Fourier-Laplace transforming we obtain the frequency dependent linear response relations

(6.114)

(6.115)

where the Fourier-Laplace transform of χ(t) is defined to be

(6.116)

The linear constitutive relation for the frequency dependent shear viscosity is (§2.4),

(6.117)

so that the frequency dependent viscosity is

(6.118)

This expression shows that the complex frequency dependent shear viscosity is given by ratio of two susceptibilities. However, these two different time correlation functions can be related by using the Nosé-Hoover equation of motion (6.101),

(6.119)

In the frequency domain this relation becomes,

(6.120)

The frequency dependent shear viscosity in the constant stress ensemble can be written as,

(6.121)

In a similar way it is possible to write the frequency dependent viscosity in terms of either the Norton ensemble stress autocorrelation function, or the Norton ensemble stress-strain cross correlation function. Using equation (4.10), the stress autocorrelation function can be related to the strain autocorrelation function using the relation,

(6.122)

In the frequency domain this becomes,

(6.123)

Substituting this equation into (6.121) gives,

(6.124)

In terms of the cross correlation function, the frequency dependent viscosity is

(6.125)

In Figure 6.14 we show the results of a test of the theory given above. Hood, Evans and Morriss (1987) computed the strain rate autocorrelation function in the Norton ensemble and the stress autocorrelation function in the Thévenin ensemble. They then used equation (6.121) to predict the strain rate autocorrelation function on the basis of their Thévenin ensemble data. The system studied was the Lennard-Jones triple point fluid. The smooth curves denote the autocorrelation function computed in the Norton ensemble and the points give the predictions from the Thévenin ensemble data. The two sets of data are in statistical agreement. This analysis shows that in spite of the fact that the damping constant Q_γ, has a profound influence on the time dependent fluctuations in the system, the theory given above correctly relates the Q_γ-dependent fluctuations of strain rate and stress to the Q_γ-independent, frequency dependent viscosity.

Figures 6.15-17 show the various Norton ensemble susceptibilities as a function of frequency. The system is the Lennard-Jones triple point fluid.