NEMD


By: Denis J. Evans
Research School of Chemistry
Australian National University
GPO Box 414, Canberra Ph: 61-6-249 3767
A.C.T., 2601 Fax: 61-6-249 0750
Australia Email: Evans@rsc.anu.edu.au
CONTENTS
 1. Introduction.
2. Driving the Programme
2.1. Menus
2.2. INPUT FILE TWODIN
2.3. OUTPUT FILE TWODPROC
2.4. Exporting raster files to NIH Image
3. Technical Details
3.1. Equations of Motion
3.2. The Nosé-Hoover k=1,k=2,k=3,k=4 Thermostats
3.3. Diffusion Under Shear

1. Introduction.


NonEquilibrium Molecular Dynamics

is a computer programme that simulates a two-dimensional system of interacting atoms. The interactions between the atoms is purely repulsive. Depending upon the chosen temperature and density, the system can simulate fluids or solids (at a temperature of 1.0, the system freezes at a density ~0.9 - see display on next page for coexisting crystal and fluid). One can also apply a shear to the system to study the behaviour of a nonequilibrium system under the influence of a uniform applied shear rate. With a thermostat applied, you can the detailed behaviour and structure of a nonequilibrium steady state.

You can use the mouse to change either the conditions under which the fluid is simulated (temperature, density or shear rate), or to change the way in which you view the fluid. One can look at the instantaneous atomic positions, a trace of the past and present atomic positions, the distribution of atoms about any given atom and the distribution of x and y particle velocities. Other display modes show you the history of the pressure and the shear stress and the diffusion coefficient as a function of time.

For more detailed numerical work one can choose the precise state points for study, by editing an input ASCII file, TWODIN , and obtain outputs from a second ASCII file, TWODPROC . Time records of pressure and shear stress are available from tab delimited ASCII files which can then be read into graphing programmes. Similarly the velocity and pair distribution functions can be written to binary files which can be imported into image processing programmes for more detailed analysis.




 Sample graphics screen using the trace display in the fluid - solid coexistence region.


2. Driving the Programme


2.1. menus


Apple
Allows you access to standard Macintosh Desk Accessories. It plays no direct role in the execution of NonEquilibrium Molecular Dynamics . You can however use this menus to find out the name and address of the programme's author.

File

Write out distribution
Writes out the near or far field pair distribution functions, the peculiar velocity distribution or the peculiar kinetic energy histogram depending on which of these distributions/histograms is being displayed in the current NEMD Graphics window. The pair and velocity distributions are written in single byte binary form and can be read by NIH Image or NCSA Image . In order for these programmes to read the image files produced by NEMD you will have to take note of imagsize given in the current values of the NEMD Data window. The peculiar kinetic energy histogram is exported as tab-delimitted ASCII text which can easily be read by Cricket Graph, KaleidaGraph or spreadsheeets for graphing. This permits an easy demonstration of the equilibrium Maxwell-Boltzmann distribution.

Write out time record
Writes out a time record containing:
time Pxy p Pxx-Pyy <Dx(t)2> <Dx(t)2> <Dx(t).Dy(t)>
where Dx(t)=x(t)-x(0) and <...> denotes an average taken over all the particles in the system. The time origin for the displacements is set when the Einstein display is called up from the Display menu. When the shear is nonzero the coordinate displacements are not laboratory displacements but are instead integrals of the SLLOD momenta appearing in the equations of motion for shear flow, Dx(t) = 0

tds px(s)/m, where dx/dt = px/m +gy. The time record is written out as tab delimited text and can be read by word processors and graphics programmes.

Quit
Enables you to exit from the programme.


 Edit
Is only included to allow the Apple menu to function. It serves no purpose in the execution of NonEquilibrium Molecular Dynamics .


 Run Parameters
Most of the items under this menu are self explanatory.

The timestep is the time increment used to solve the differential equations of motion for the system. A typical value for this variable is 0.001-0.003. If the shear rate is increased to values larger than unity the timestep will have to be decreased otherwise the programme will not conserve energy or maintain a constant temperature correctly. If this problem becomes severe the programme can bomb with a series of messages stating that there are too many particles in a cell.

Kplot is the number of timesteps between updates of the graphics and text windows. The programme runs faster for larger values of kplot . The number of timesteps per hour is indicated under Current Values as ts / hr = in the NEMD Data window.

Restart simulation
Enables the user to start the programme from a crystal. This same result can be achieved by editing the ASCII file TWODIN and changing the value of NTYPE from 3 to 1. If the programme bombs this is the most efficient way of restarting the programme. If this is the case remember to decrease the shear rate (SHEAR ), temperature (TR ), density (DR ) or timestep (DELTA ) that initially caused the programme to bomb .

 Reverse direction of time
Does what it says. This is useful in teaching the reversibility of the equations of motion with or without a thermostat. For example if one starts a zero shear simulation at a temperature and density where the fluid is the stable equilibrium phase, one can observe that the initial crystalline phase melts. If this occurs in a short enough time one can return to the thermodynamically unstable crystal phase by reversing the direction of time. Likewise in a shearing fluid one can defy the Second Law of Thermodynamics for a short while generating a negative viscosity phase where heat is converted into work using reverse direction of time.


Display
Atoms
Shows the atomic positions as a (slow!) movie. The positions are updated at a rate that you choose by setting the value of kplot in the Run Parameters menu. The atoms are coloured according to their individual peculiar kinetic energies, those with higher kinetic energy (relative to the local average velocity) are shown in red). The peculiar velocity of a particle is the velocity measured with respect top the local streaming velocity . One can change the definition of the streaming velocity used for this purpose from within the Thermostat menu. The colour map for the Atoms and the Trace display is linear in the peculiar kinetic energy of the particles. There are 13 bins each kBT/6 wide covering the range of kinetic energies from 0 to 2.166kBT.

 Trace
Plots the atomic positions a pixels which are coloured according to their individual peculiar velocities.. These pixels are not erased as time progresses so that one can trace the past positions of individual atoms. This is very useful for showing the difference between the diffusive behaviour of atoms in the fluid phase and the trapped oscillatory behaviour of atoms in solids.

Kinetic Energy Histogram
Displays the distribution of the peculiar kinetic energies of the atoms. As time evolves the displayed distribution denotes the cumulative time average. At equilibrium the average distribution should be the exponential Boltzmann distribution exp(kinetic energy/temperature). The colours used in the histogram use a different colour map from that used in the atoms and trace displays. For the peculiar kinetic energy histogram there are 62 bins which covering the range of kinetic energies from 0 to 6.2kBT. These bins are coloured with 14 different colours.




Sample Nemd Data window. These parameters were used to generate the pair distribution shown below.
Sample Pair distribution wide field display.
Pair distribution near field
Shows the relative positions of particles in the system. If any given particle is located at the centre of the screen the other particles pack around it generating the so-called pair distribution. Neighbours are excluded from coming too close to the central particle by the repulsive interatomic forces. The plot is coloured according to the probability density - red denotes the highest probability while blue denotes the lowest. In a fluid phase the pair distribution function is isotropic while in the solid a sixfold pattern of spots is revealed. This distribution can be written to disk for importing into image programmes using Write out Distribution from within the file menu. Near field refers to the fact that we only show the first coordination shell.

Pair distribution wide field
Same as above, but but showing the first three coordination shells. The increased range in the pair distribution slows the programme down considerably. This distribution can be written to disk for importing into image programmes using Write out Distribution from within the file menu. Using NIH IMAGE the pair distribution function can be Fourier transformed to give the static structure factor which can be measured experimentally.

Velocity Distribution
Shows the probability distribution of peculiar velocity. The white lines show the zeros of the x and y velocities.The range of peculiar velocities displayed in this window is 4 (kBT/m). This distribution can be written to disk for importing into image programmes using Write out Distribution from within the file menu.



Sample output in the velocity distribution window for a shearing fluid.

Time Record
Shows the past history of the indicated elements of the pressure tensor as a function of time. This is very useful in studying the transient behaviour of systems to a suddenly applied shear rate. One can easily study stress overshoot by this means. The time record can be written out as a tab delimited text file which can be read by word processors and graphics programmes - use write out time record in the file menu.


Einstein Diffusion
Enables you to use Einstein's method of calculating the self diffusion coefficient for a fluid. You can verify that the diffusion coefficient is zero for a solid. In an equilibrium fluid the diffusion process is isotropic and Dxx=Dyy while Dxy=0. In a fluid under shear neither of these conditions is satisfied. In high density fluids under shear the diffusion tensor is anisotropic but very nearly diagonal. At lower densities, Dxy can become quite large. When this display is active the time record which is written out to disk under the File menu includes the mean square displacements of the particles as a function of time.

Thermostat
Allows you to change the thermostats used in the simulation. The thermostat which is currently operating is shown in the Nemd Data window.
None
Is self explanatory. In the absence of an applied shear rate one can observe that Newton's equations of motion conserve energy. This can be used to check whether the current timestep is sufficiently small for energy to be conserved satisfactorily - typically < 0.2% drift in say 2,000 timesteps. If no thermostat is used for a system under shear, one can observe that the work done by the shear is converted into heat, raising the temperature of the system. At equilibrium, time averaging along an unthermostatted trajectory is equivalent to a microcanonical ensemble average.
Gauss
Runs at constant peculiar kinetic energy. In this case the streaming velocity at a point r = (x,y), is assumed to be i gy where i is the x unit vector, g is the value of the strain rate. At equilibrium, time averaging along a Gaussian isokinetic trajectory is equivalent to a microcanonical ensemble average for the kinetic part of phase functions and a canonical ensemble average for the configurational part of phase functions.
Nosé-Hoover
Uses an integral feedback mechanism to constrain the average peculiar kinetic energy. In this case the streaming velocity at a point r = (x,y), is assumed to be i gy where i is the x unit vector, g is the value of the strain rate. At equilibrium, time averaging along a Nosé-Hoover thermostatted trajectory is equivalent to a canonical ensemble average.

Nosé-Hoover k=1, k=2, k=3, k=4
Employs the Nosé-Hoover integral feedback thermostatting mechanism but uses a Fourier series technique to compute the instantaneous streaming velocity. Within Lees-Edwards shearing periodic boundary conditions if v is the laboratory velocity of any particle then v -i gy is a periodic function and can therefore be represented by a Fourier series. k=1,2,3,4 refers to how many Fourier harmonics are used to represent this periodic function. The peculiar velocity of any particle is then defined as the instantaneous difference of v -i gy from the harmonic representation of the periodic velocity field. In this case the streaming velocity is associated with the instantaneous long wavelength velocity modes while the peculiar velocity is assumed to be given by the short wavelength modes. Choosing different values of k allows the user to choose the breakpoint between the 'long' and the short' wavelength modes. At equilibrium the time averages of phase functions and time correlation functions are independent of k. The value of k chosen from this menu is shown under the Current Values in the NEMD Data window.

Palette
Enables you to changes the colour palette for the programme. These changes are not saved to disk.

2.2. INPUT FILE TWODIN

The programme can be driven by editing the input ASCII file TWODIN . A typical input looks like:

	  1.00000  .900000  3.500000E-03
	  1.12246  1  .000000  100.000  0
	  1093076  3  0  1
	  50  20  20  1.00000  .300000



	TR,DR,DELTA
	RCUT,NTIME,SHEAR,TAU,Nklim
	KBMAX,NTYPE,NON,NGAUS
	KSCAL,KPROC,KPLOT,scal,RDEL
Thus precise values of the set temperature density and shear rate may be chosen by using an ASCII editor. A dictionary of variable names used in this file follows:
TR desired temperature
DR desired density
DELTA timestep for integration of equations of motion
RCUT cutoff distance for interaction potential. If this is near 21/6 the programme employs a WCA potential. If it is set at a larger value say, 2.5, the programme will employ a Lennard-Jones potential with the specified cutoff distance. Clearly a Lennard-Jones potential will run more slowly than a WCA potential.
NTIME unused
SHEAR the strain rate, ux/ y
TAU the Nosé-Hoover time constant t, used in the equation of motion for the thermostatting multiplier, a
da/dt = [S pi2/2m]/NkBT - 1)/t
NKLIM Number of harmonics, Nk, used to define the streaming velocity. The number of degrees of freedom in the system is 2N-2-(2Nk+1)2. N is the number of particles moving in 2 dimensions. For example the temperature should be calculated from the equation,
kBT = S pi2/{m[2N-2-(2Nk+1)2]}
KBMAX The max number of timesteps before the programmes stops.
NTYPE Control parameter: 1. means start from crystal This is very important for restarts when the programme is bombing.
 NON control parameter: leave set to zero.
NGAUS control parameter: 0. means turn off thermostats
KSCAL not used
KPROC frequency of writing to disk. Note quitting from the programme automatically invokes a full write to disk of all restart and averaging files.
KPLOT frequency of drawing the screen, in timesteps
SCAL not used
RDEL thickness of neighbour list shell. 0.3 works reasonably well but programme speed may be improved by tuning this value and monitoring programme speed in timesteps/hr in the NEMD DATA window.

2.3. OUTPUT FILE TWODPROC

This file can be opened with a word processor or ASCII file editor to show the accumulated thermodynamic averages. a typical output might be:

          SIMULATION PARAMETERS

     Set Temperature           =     1.0000000
     Density                   =      .9000000
     Current timestep          =      .0035000
     Potential Cutoff distance =     1.1224620
     Number of Particles       =           224
     Simulation Cell Length    =    15.7762120
     Simulation Cell Area      =   248.8889000

     NTYPE  =     3
     NON    =     0          NGAUS  =     1
     No of wavevectors defining streaming velocity  =     0


     KB =    12980     No of steps  =        12980     TIME =        
45.4300

          Simulation Averages

     Temp using Nk-vectors for streaming velocity =   1.00000
     Temp using 0k-vectors for streaming velocity =   1.00000
     Thermostat Mutiplier, Alpha   =   -.00009
     Shear Rate                    =    .00000
     Pxx - Pyy                     =    .06965
     Pressure, 0.5*(Pxx +Pyy)      =   9.93919
     Internal  Energy per Particle =   1.78405
     Potential Energy per Particle =    .78403
     xy-element of pressure tensor =   -.05467

          Diffusion Coefficients
     D(t)  =<(delr(t))**2>/4t     = none
     Dxx(t)=<(delx(t))**2>/4t     = none
     Dyy(t)=<(dely(t))**2>/4t     = none
     Dxy(t)=<(delx(t))*dely(t)>/4t= none

This is all fairly self explanatory. To obtain estimates of the various elements of the diffusion tensor, Einstein Diffusion under the Display menu must have been selected.


2.4. Exporting raster files to NIH Image
 The default size of the NEMD Graphics windows is 256*256. This is convenient for using Fast Fourier Transform routines and computing the static structure factor since 256 is a power of 2. The size of these images can however be changed by the usual method of 'stretching' windows. For exporting image files to NIH Image or NCSA Image, you will need to take note of the value of imagsize as reported in the NEMD Data window.
We will now give the instructions for exporting an image to NIH Image for computing the structure factor. After the image generation is complete, choose Write out Distribution from the File menu. At the prompt for choosing a name for the raster file we suggest that you use the value of imagsize in the file name so that you will not forget it. Then quit from NEMD . Open NIH Image and select import from its file menu. At the prompt select special , 8bit and fill in the value of imagesize in the dialog box. Click OK and the image will be displayed.
If imagsize is a power of 2 you can compute the Fourier Transform.

3. Technical Details


3.1. Equations of Motion
 The 224 particles interact via the WCA pair potential f(r)=4e[(s/r)12-(s/r)6], r<21/6; f(r)=0, r>21/6. The readouts are given in reduced units for which the potential parameters e, s are unity as is the atomic mass, m. The system is subject to a shear rate g = ux/ y and the equations of motion are given by the thermostatted SLLOD equations1.



(1)



Using the reduced units q i stands for q i/s, p i for p i/(me)1/2, t for t (e/ms2)1/2 with (e/ms2)1/2 ( 5 x 1011Hz for argon) and g for g(ms2/e)1/2. i is a unit vector in the x-direction. The p i are peculiar momenta, defined in terms of the peculiar velocities, i.e., the velocities of the particles with respect to the (local) fluid velocity ux(y)=gy. The last term in (1) represents an thermostat which removes heat from the system so as to keep the kinetic temperature, Tk= (1/NkB)Si(pi2/2m) constant. Pxy is the xy element of the pressure tensor1, PxyV = Si[pxipyi/m +Fxiyi]. The equations (1) are equivalent to Newton's equations of motion in the presence of frictional forces -ap i, if a constant shear rate g is imposed on the fluid at t = 0+. This is illustrated in the figure below.
The equations of motion (1) are supplemented by Lees-Edwards periodic boundary conditions.[1]. In their adiabatic form the SLLOD equations of motion give an exact description of adiabatic planar Couette flow arbitrarily far from equilibrium.
These time varying periodic boundary conditions can be represented in a number of apparently different but equivalent ways. One such diagrammatic representation is shown below. The blue hatched square is the primitive cell (i.e. represents the coordinates of the set of N particles that we arbitrarily choose to define the N-particle system configuration. The red hatched square is the Minimum image cell of particle (1,j). It is a cell of the same size as the primitive cell but it is centred on particle (1,j). Within the infinite array of periodic images of particles, it contains those that are closest to (1,j). The particles that can interact with (1,j) are chosen from those within the minimum image cell of (1,j). Thus (1,j) can interact with (9,i) but not with (1,i) even though (1,i) is in the primitive cell.



During the course of time particles diffuse and convect through the system. If a particle moves through a vertical face of the primitive cell it is replaced in the primitive cell by one of its images moving with the same SLLOD and laboratory velocities, with the same y coordinate but with an x-coordinate L different from its own x-coordinate. If a particle moves out of the primitive cell through an horizontal face, matters are more complex. It will be replaced by an image of itself but in this case the image particle differs in y by L. Because however of the time dependent horizontal motion of the cells above and below the primitive cell, the x-coordinate can also be different from that of the particle in the primitive cell that it replaces.
Further, although it enters the primitive cell with the same SLLOD momentum as its image its laboratory momentum will be different. From the figure you can deduce that the SLLOD momentum of any particle is identical to that of each of its images. This is obviously not so for the laboratory velocities.
Lees-Edwards periodic boundary conditions.
 3.2. The Nosé-Hoover k=1,k=2,k=3,k=4 Thermostats
Consider a fluid of N particles of mass, m, interacting via a pair potential, f(r), in volume V (density, n = N/V) in planar Couette flow. The properties of the system can be simulated with the SLLOD algorithm1, supplemented by Lees-Edwards periodic boundary conditions, which are known to give an exact description of adiabatic planar Couette flow. If the algorithm is supplemented with a thermostat to keep the kinetic temperature, T, constant, the SLLOD equations for isothermal flow in a Cartesian coordinate system with a gradient in the y-direction and flow in the x-direction can be written;

(2)

where the term involving a thermostatting multiplier, a, in (2) represents the thermostat whose form is determined from Gauss' principle of least constraint.
In Eq.(2), F i = - S j f(rij)/ r i , i is the unit vector in the flow direction and p i is the SLLOD momentum of particle i. The velocity u S(r i) is the local SLLOD streaming velocity at the position of particle i, defined in terms of the local number density, n(r ) S id(r -r i), and the momentum current, J S(r ) S ip id(r -r i) mn(r )u S(r ). Note that the streaming velocity measured in the laboratory coordinate frame, u L(r ), is the sum of nonzero and the zero wavevector components: viz, u L(r ) = u S(r ) + i gy, where g is formally the zero wavevector component of the strain rate which governs the motion of the infinite array of shearing Lees-Edwards1 unit cells.
At low Reynolds number the stable laboratory streaming velocity is i gy and the average SLLOD streaming velocity is zero. Hence, at low Reynolds number, the peculiar velocity of any particle, i, is given by dr i/dt - i gyi = p i/m, and the kinetic temperature is defined in the usual way, namely; T = S i(p i2)/2mNkB, where kB is Boltzmann's constant.
At higher Reynolds number we cannot make a priori any prediction about the finite wavevector components of the streaming velocity, it is clear that the SLLOD momenta, {p i}, defined in Eq.(2), are periodic functions of displacement along the time dependent Lees-Edwards axes. The Lees-Edwards periodic boundary condition implies that the SLLOD momentum of any particle is identical to the SLLOD momentum of any image of that same particle. Thus the instantaneous finite wavevector components of the streaming velocity can be computed from a Fourier series analysis as the simulation evolves in time.
The number density, n(r ), and SLLOD momentum current, J S(r ), are periodic functions of the primitive lattice vectors describing the instantaneous alignment of the Lees-Edwards periodic boundaries. These lattice vectors are time dependent and non-orthogonal. We compute n(k ) and J S(k ) for a number of k -vectors given by (Nx, Ny) = (2imax+1, 2jmax+1). On transforming from the k -domain back to real space, the SLLOD streaming velocity u S(r ) is evaluated from its defining equation. If imax,jmax are equal and nonzero for each Cartesian component of u S, the algorithm allows the system to take on any possible form for the instantaneous streaming velocity as a function of position as long as the length scale for variations in the streaming velocity is greater than L/(2imax) lmin.
In this programme we set, imax = jmax for both Cartesian components of u S. Thus Nosé-Hoover k=1 refers to the situation where imax = imax = 1, Nosé-Hoover k=2 refers to the situation where imax = imax = 2,... etc.
Once the local streaming velocity is known, the kinetic temperature follows;

 (3)
where Nk = Nx.Ny is the total number of k-vectors used to determine the local streaming velocity. We use the identity at equilibrium of the temperature from (3) with the usual equipartition relation for the temperature, to verify our coding of the PUT algorithm.
The simulations use a Nosé-Hoover thermostat to keep the temperature constant; that is, the multiplier,a, of Eq. (2) fixed the average peculiar kinetic energy, defined by Eq. (3), to the preset value, Ts , from the feedback equation,

(4)

 where the reciprocal of t the feedback time, is set in the ASCII file TWODIN .

3.3. Diffusion Under Shear
The self diffusion coefficient describes the diffusive flux, J , of labelled atoms induced to flow by the presence of a concentration gradient, -- n, in the number density of those labelled atoms. The convective contributions to the flux are removed. Thus the macroscopic equations describing diffusion under shear are the defining constitutive relation,

(5)

and the convective diffusion equation,

(6)

When a single component fluid is subject to steady planar Couette flow, g= ux/ y, the self diffusion coefficient becomes a self diffusion tensor (SDT) with shear rate dependent components. In Cartesian coordinates the steady state diffusion tensor can be written as,

(7)

The xz,yz,zx and zy are identically zero because of symmetry. In previous work2 the following Green-Kubo relations for the various components of this tensor were derived.

(8)

where a,b=x,y,z. The ensemble average is taken over shearing steady states, hence the subscript g. We can operationally define this ensemble by considering some equilibrium ensemble at t = - , being subject for all subsequent time, to thermostatted planar Couette flow (equations (1,2)). At t=0 the system is assumed to have completely relaxed to a steady state. Thus the time correlation function referred to in (8) correlates the Cartesian components of momentum at times 0,t long after the transients leading to the establishment of the steady state, have decayed.
Note that the pa in this equation are the SLLOD momenta of (1). These relations can be integrated to give the Einstein relations for the mean square displacement. The diagonal components become

(9)

where

(10)

is the convected, Lagrangian position of particle i. It is not possible to obtain separate Einstein relations for the two nonzero off diagonal elements of the diffusion tensor3 it is only possible to obtain an Einstein relation for the sum of these two elements,

(11)

In general the self diffusion tensor for a fluid under Couette flow is non symmetric. The x mass current generated by a y-concentration gradient is not necessarily equal to the y-current generated by an x-gradient. However, it is easy to prove that by measuring the time dependence of concentrations alone, it is only possible to determine the symmetric part of the diffusion tensor3.
An alternative route to the Green-Kubo relations for the diagonal components of the Self Diffusion Tensor is to derive the above Mean Square Displacements from the macroscopic convective diffusion equation3,4,5. One then assumes that the macroscopic and microscopic MSDs are the same and GK relations are thus recovered from the MSDs instead of the the other way around. Unfortunately, it is not possible to derive any GK relations for the off diagonal elements by this method.
Final remarks: at equilibrium (g=0), the relations given above reduce to the standard equilibrium Green-Kubo and Einstein relations for the self diffusion tensor. When any of the Nosé-Hoover k=1,k=2,k=3,k=4 Thermostats are selected the expressions given above for the diffusion tensor are evaluated with the p i above being peculiar momenta evaluated relative to the streaming velocity obtained from the Fourier-Series fit of the streaming velocity using k=1,2,3 or 4.
References

[1] D.J. Evans and G.P. Morriss, Statistical Mechanics of Nonequilibrium Liquids, (Academic Press, London, 1990).
[2] D.J. Evans, Phys. Rev. A 44 , 3630 (1991).
[3] S. Sarman, D.J. Evans and P.T. Cummings, J. Chem. Phys. 95 , 8675 (1991).
[4] P.T. Cummings, B.Y. Wang, D.J. Evans and K.J. Fraser, J. Chem. Phys. 94 , 2149 (1991).
[5] T.G.M. van de Ven, Colloidal Hydrodynamics , (Academic Press, London, 1989).

12th January 1996