lammps/doc/fix_viscosity.txt

"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c

:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)

:line

fix viscosity command :h3

[Syntax:]

fix ID group-ID viscosity N vdim pdim Nbin keyword value ... :pre
  
ID, group-ID are documented in "fix"_fix.html command :ulb,l
viscosity = style name of this fix command :l
N = perform momentum exchange every N steps :l
vdim = {x} or {y} or {z} = which momentum component to exchange :l
pdim = {x} or {y} or {z} = direction of momentum transfer :l
Nbin = # of layers in pdim direction (must be even number) :l

zero or more keyword/value pairs may be appended :l
keyword = {swap} or {target} :l
  {swap} value = Nswap = number of swaps to perform every N steps
  {vtarget} value = V or INF = target velocity of swap partners (velocity units) :pre
:ule

[Examples:]

fix 1 all viscosity 100 x z 20
fix 1 all viscosity 50 x z 20 swap 2 vtarget 1.5 :pre

[Description:]

Use the Muller-Plathe algorithm described in "this
paper"_#Muller-Plathe to exchange momenta between two particles in
different regions of the simulation box every N steps.  This induces a
shear velocity profile in the system.  As described below this enables
a viscosity of the fluid to be calculated.  This algorithm is
sometimes called a reverse non-equilibrium MD (reverse NEMD) approach
to computing viscosity.  This is because the usual NEMD approach is to
impose a shear velocity profile on the system and measure the response
via an off-diagonal component of the stress tensor, which is
proportional to the momentum flux.  In the Muller-Plathe method, the
momentum flux is imposed, and the shear velocity profile is the
system's response.

The simulation box is divided into {Nbin} layers in the {pdim}
direction, where the layer 1 is at the low end of that dimension and
the layer {Nbin} is at the high end.  Every N steps, Nswap pairs of
atoms are chosen in the following manner.  Only atoms in the fix group
are considered.  Nswap atoms in layer 1 with positive velocity
components in the {vdim} direction closest to the target value {V} are
selected.  Similarly, Nswap atoms in the "middle" layer (see below) with
negative velocity components in the {vdim} direction closest to the
negative of the target value {V} are selected.  The two sets of Nswap
atoms are paired up and their {vdim} momenta components are swapped
within each pair.  This resets their velocities, typically in opposite
directions.  Over time, this induces a shear velocity profile in the
system which can be measured using commands such as the following,
which writes the profile to the file tmp.profile:

fix f1 all ave/spatial 100 10 1000 z lower 0.05 vx &
    file tmp.profile units reduced :pre

Note that by default, Nswap = 1 and vtarget = INF, though this can be
changed by the optional {swap} and {vtarget} keywords.  When vtarget =
INF, one or more atoms with the most positive and negative velocity
components are selected.  Setting these parameters appropriately, in
conjunction with the swap rate N, allows the momentum flux rate to be
adjusted across a wide range of values, and the momenta to be
exchanged in large chunks or more smoothly.

The "middle" layer for momenta swapping is defined as the {Nbin}/2 + 1
layer.  Thus if {Nbin} = 20, the two swapping layers are 1 and 11.
This should lead to a symmetric velocity profile since the two layers
are separated by the same distance in both directions in a periodic
sense.  This is why {Nbin} is restricted to being an even number.

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"_#Muller-Plathe for
details.

IMPORTANT NOTE: After equilibration, if the velocity profile you
observe is not linear, then you are likely swapping momentum too
frequently and are not in a regime of linear response.  In this case
you cannot accurately infer a viscosity and should try increasing
the Nevery parameter.

An alternative method for calculating a viscosity is to run a NEMD
simulation, as described in "this section"_Section_howto.html#4_13 of
the manual.  NEMD simulations deform the simmulation box via the "fix
deform"_fix_deform.html command.  Thus they cannot be run on a charged
system using a "PPPM solver"_kspace_style.html since PPPM does not
currently support non-orthogonal boxes.  Using fix viscosity keeps the
box orthogonal; thus it does not suffer from this limitation.

[Restart, fix_modify, output, run start/stop, minimize info:]

No information about this fix is written to "binary restart
files"_restart.html.  None of the "fix_modify"_fix_modify.html options
are relevant to this fix.

This fix computes a global scalar which can be accessed by various
"output commands"_Section_howto.html#4_15.  The scalar is 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.  The scalar value calculated by
this fix is "intensive".

No parameter of this fix can be used with the {start/stop} keywords of
the "run"_run.html command.  This fix is not invoked during "energy
minimization"_minimize.html.

[Restrictions:]

If the masses of all exchange partners are the same, then swaps
conserve both momentum and kinetic energy.  Thus you should not need
to thermostat the system.  If you do use a thermostat, you may want to
apply it only to the non-swapped dimensions (other than {vdim}).

LAMMPS does not check, but you should not use this fix to swap
velocities of atoms that are in constrained molecules, e.g. via "fix
shake"_fix_shake.html or "fix rigid"_fix_rigid.html.  This is because
application of the constraints will alter the amount of transferred
momentum.  You should, however, be able to use flexible molecules.
See the "Maginn paper"_#Maginn for an example of using this algorithm
in a computation of alcohol molecule properties.

When running a simulation with large, massive particles or molecules
in a background solvent, you may want to only exchange momenta bewteen
solvent particles.

[Related commands:]

"fix ave/spatial"_fix_ave_spatial.html, "fix
thermal/conductivity"_fix_thermal_conductivity.html

[Default:]

The option defaults are swap = 1 and vtarget = INF.

:line

:link(Muller-Plathe)
[(Muller-Plathe)] Muller-Plathe, Phys Rev E, 59, 4894-4898 (1999).

:link(Maginn)
[(Maginn)] Kelkar, Rafferty, Maginn, Siepmann, Fluid Phase Equilibria,
260, 218-231 (2007).