We recently did some simulations to calculate the viscosity using the Muller-Plathe method. Using units=real, when we divide the momentum flux (j) over the velocity gradient, I thought to obtain the viscosity directly in Poise units (without any further conversion). However, this way I get crazy values for the viscosity.
Is there some conversion factor that I am missing?
I guess the output of the fix, namely the total momentum, is in units of
(grams/mole) * (A/fs). Then one has to divide it by time (in fs) and
cross-sectional area (in A^2) to get the momentum flux. Finally the latter has
to be divided by the shear rate (in fs^-1). The result should therefore be in
units of (grams/mole)/(A*fs). To get the viscosity in Pa s = kg/(m s), one has
to apply the following conversion factor:
1 gram/mole/A/fs = 10^-3 / 6.02e23 kg * 10^10 m^-1 * 10^15 s^-1
I tried to use the Muller-Plathe method recently and to get results consistent
with other methods, I had to divide the results by 2, namely really apply the
definition (eq 3) of Muller-Plathe article. Is it right? In this case, maybe
this should be precised in the documentation...
Laurent is correct about the units - you need to convert mass*velocity
(grams/mole) * (A/fs) to Poise.
As the doc page for fix viscosity states:
As described below, the total momentum transferred by these velocity swaps is computed by the fix and >can be output. Dividing this quantity by time and the cross-sectional area of the simulation box yields a >momentum flux. The ratio of momentum flux to the slope of the shear velocity profile is the viscosity of >the fluid, in appopriate units. See the Muller-Plathe paper for details.
The cummulative momentum transferred between the bottom and middle of the simulation box (in the >pdim direction) is stored as a scalar quantity by this fix. This quantity is zeroed when the fix is defined and >accumlates thereafter, once every N steps. The units of the quantity are momentum = mass*velocity.
You also need to include the factor of 2 that the MP paper includes in
its formula, since the
system is relaxing in 2 directions due to periodicity. (I'll note
this in the doc page).
LAMMPS does not do these things for you, b/c you might use this fix in
You might apply it to a subset of atoms, or use the flux quantity for
some other calculation.
If you want a viscosity, you need to use the correct formulas and
do the correct unit conversion.