ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

use k_B consistently, fix some more math. convert some transcriptions to math

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-02-29 11:37:31 -05:00
parent 07fc624509
commit 6f47f110f1
11 changed files with 274 additions and 258 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
doc/src/compute_ke.rst
View File

@ -28,8 +28,8 @@ Description
Define a computation that calculates the translational kinetic energy Define a computation that calculates the translational kinetic energy
of a group of particles. of a group of particles.
The kinetic energy of each particle is computed as 1/2 m v\^2, where m The kinetic energy of each particle is computed as :math:`\frac{1}{2} m
and v are the mass and velocity of the particle. v^2`, where *m* and *v* are the mass and velocity of the particle.
There is a subtle difference between the quantity calculated by this There is a subtle difference between the quantity calculated by this
compute and the kinetic energy calculated by the *ke* or *etotal* compute and the kinetic energy calculated by the *ke* or *etotal*
@ -37,12 +37,13 @@ keyword used in thermodynamic output, as specified by the
:doc:`thermo_style <thermo_style>` command. For this compute, kinetic :doc:`thermo_style <thermo_style>` command. For this compute, kinetic
energy is "translational" kinetic energy, calculated by the simple energy is "translational" kinetic energy, calculated by the simple
formula above. For thermodynamic output, the *ke* keyword infers formula above. For thermodynamic output, the *ke* keyword infers
kinetic energy from the temperature of the system with 1/2 Kb T of kinetic energy from the temperature of the system with
energy for each degree of freedom. For the default temperature :math:`\frac{1}{2} k_B T` of energy for each degree of freedom. For the
computation via the :doc:`compute temp <compute_temp>` command, these default temperature computation via the :doc:`compute temp
are the same. But different computes that calculate temperature can <compute_temp>` command, these are the same. But different computes
subtract out different non-thermal components of velocity and/or that calculate temperature can subtract out different non-thermal
include different degrees of freedom (translational, rotational, etc). components of velocity and/or include different degrees of freedom
(translational, rotational, etc).
**Output info:** **Output info:**

19
doc/src/compute_ke_atom_eff.rst
View File

@ -18,7 +18,7 @@ Examples
"""""""" """"""""
.. parsed-literal:: .. code-block:: LAMMPS
compute 1 all ke/atom/eff compute 1 all ke/atom/eff
@ -30,11 +30,12 @@ Define a computation that calculates the per-atom translational
group. The particles are assumed to be nuclei and electrons modeled group. The particles are assumed to be nuclei and electrons modeled
with the :doc:`electronic force field <pair_eff>`. with the :doc:`electronic force field <pair_eff>`.
The kinetic energy for each nucleus is computed as 1/2 m v\^2, where m The kinetic energy for each nucleus is computed as :math:`\frac{1}{2} m
corresponds to the corresponding nuclear mass, and the kinetic energy v^2`, where *m* corresponds to the corresponding nuclear mass, and the
for each electron is computed as 1/2 (me v\^2 + 3/4 me s\^2), where me kinetic energy for each electron is computed as :math:`\frac{1}{2} (m_e
and v correspond to the mass and translational velocity of each v^2 + \frac{3}{4} m_e s^2)`, where :math:`m_e` and *v* correspond to the mass
electron, and s to its radial velocity, respectively. and translational velocity of each electron, and *s* to its radial
velocity, respectively.
There is a subtle difference between the quantity calculated by this There is a subtle difference between the quantity calculated by this
compute and the kinetic energy calculated by the *ke* or *etotal* compute and the kinetic energy calculated by the *ke* or *etotal*
@ -43,8 +44,8 @@ keyword used in thermodynamic output, as specified by the
energy is "translational" plus electronic "radial" kinetic energy, energy is "translational" plus electronic "radial" kinetic energy,
calculated by the simple formula above. For thermodynamic output, the calculated by the simple formula above. For thermodynamic output, the
*ke* keyword infers kinetic energy from the temperature of the system *ke* keyword infers kinetic energy from the temperature of the system
with 1/2 Kb T of energy for each (nuclear-only) degree of freedom in with :math:`\frac{1}{2} k_B T` of energy for each (nuclear-only) degree
eFF. of freedom in eFF.
.. note:: .. note::
@ -53,7 +54,7 @@ eFF.
command, as shown in the following example: command, as shown in the following example:
.. parsed-literal:: .. code-block:: LAMMPS
compute effTemp all temp/eff compute effTemp all temp/eff
thermo_style custom step etotal pe ke temp press thermo_style custom step etotal pe ke temp press

33
doc/src/compute_ke_eff.rst
View File

@ -29,11 +29,12 @@ Define a computation that calculates the kinetic energy of motion of a
group of eFF particles (nuclei and electrons), as modeled with the group of eFF particles (nuclei and electrons), as modeled with the
:doc:`electronic force field <pair_eff>`. :doc:`electronic force field <pair_eff>`.
The kinetic energy for each nucleus is computed as 1/2 m v\^2 and the The kinetic energy for each nucleus is computed as :math:`\frac{1}{2} m
kinetic energy for each electron is computed as 1/2(me v\^2 + 3/4 me v^2` and the kinetic energy for each electron is computed as
s\^2), where m corresponds to the nuclear mass, me to the electron :math:`\frac{1}{2}(m_e v^2 + \frac{3}{4} m_e s^2)`, where *m*
mass, v to the translational velocity of each particle, and s to the corresponds to the nuclear mass, :math:`m_e` to the electron mass, *v*
radial velocity of the electron, respectively. to the translational velocity of each particle, and *s* to the radial
velocity of the electron, respectively.
There is a subtle difference between the quantity calculated by this There is a subtle difference between the quantity calculated by this
compute and the kinetic energy calculated by the *ke* or *etotal* compute and the kinetic energy calculated by the *ke* or *etotal*
@ -42,18 +43,20 @@ keyword used in thermodynamic output, as specified by the
energy is "translational" and "radial" (only for electrons) kinetic energy is "translational" and "radial" (only for electrons) kinetic
energy, calculated by the simple formula above. For thermodynamic energy, calculated by the simple formula above. For thermodynamic
output, the *ke* keyword infers kinetic energy from the temperature of output, the *ke* keyword infers kinetic energy from the temperature of
the system with 1/2 Kb T of energy for each degree of freedom. For the system with :math:`\frac{1}{2} k_B T` of energy for each degree of
the eFF temperature computation via the :doc:`compute temp\_eff <compute_temp_eff>` command, these are the same. But freedom. For the eFF temperature computation via the :doc:`compute
different computes that calculate temperature can subtract out temp\_eff <compute_temp_eff>` command, these are the same. But
different non-thermal components of velocity and/or include other different computes that calculate temperature can subtract out different
degrees of freedom. non-thermal components of velocity and/or include other degrees of
freedom.
IMPRORTANT NOTE: The temperature in eFF models should be monitored via .. warning::
the :doc:`compute temp/eff <compute_temp_eff>` command, which can be
printed with thermodynamic output by using the
:doc:`thermo_modify <thermo_modify>` command, as shown in the following
example:
The temperature in eFF models should be monitored via
the :doc:`compute temp/eff <compute_temp_eff>` command, which can be
printed with thermodynamic output by using the
:doc:`thermo_modify <thermo_modify>` command, as shown in the following
example:
.. parsed-literal:: .. parsed-literal::

30
doc/src/compute_pressure.rst
View File

@ -21,7 +21,7 @@ Examples
"""""""" """"""""
.. parsed-literal:: .. code-block:: LAMMPS
compute 1 all pressure thermo_temp compute 1 all pressure thermo_temp
compute 1 all pressure NULL pair bond compute 1 all pressure NULL pair bond
@ -42,20 +42,20 @@ The pressure is computed by the formula
P = \frac{N k_B T}{V} + \frac{\sum_{i}^{N'} r_i \bullet f_i}{dV} P = \frac{N k_B T}{V} + \frac{\sum_{i}^{N'} r_i \bullet f_i}{dV}
where N is the number of atoms in the system (see discussion of DOF where *N* is the number of atoms in the system (see discussion of DOF
below), Kb is the Boltzmann constant, T is the temperature, d is the below), :math:`k_B` is the Boltzmann constant, *T* is the temperature, d
dimensionality of the system (2 or 3 for 2d/3d), and V is the system is the dimensionality of the system (2 or 3 for 2d/3d), and *V* is the
volume (or area in 2d). The second term is the virial, equal to system volume (or area in 2d). The second term is the virial, equal to
-dU/dV, computed for all pairwise as well as 2-body, 3-body, 4-body, -dU/dV, computed for all pairwise as well as 2-body, 3-body, 4-body,
many-body, and long-range interactions, where r\_i and f\_i are the many-body, and long-range interactions, where :math:`r_i` and
position and force vector of atom i, and the black dot indicates a dot :math:`f_i` are the position and force vector of atom *i*, and the black
product. When periodic boundary conditions are used, N' necessarily dot indicates a dot product. When periodic boundary conditions are
includes periodic image (ghost) atoms outside the central box, and the used, N' necessarily includes periodic image (ghost) atoms outside the
position and force vectors of ghost atoms are thus included in the central box, and the position and force vectors of ghost atoms are thus
summation. When periodic boundary conditions are not used, N' = N = included in the summation. When periodic boundary conditions are not
the number of atoms in the system. :doc:`Fixes <fix>` that impose used, N' = N = the number of atoms in the system. :doc:`Fixes <fix>`
constraints (e.g. the :doc:`fix shake <fix_shake>` command) also that impose constraints (e.g. the :doc:`fix shake <fix_shake>` command)
contribute to the virial term. also contribute to the virial term.
A symmetric pressure tensor, stored as a 6-element vector, is also A symmetric pressure tensor, stored as a 6-element vector, is also
calculated by this compute. The 6 components of the vector are calculated by this compute. The 6 components of the vector are
@ -108,7 +108,7 @@ A compute of this style with the ID of "thermo\_press" is created when
LAMMPS starts up, as if this command were in the input script: LAMMPS starts up, as if this command were in the input script:
.. parsed-literal:: .. code-block:: LAMMPS
compute thermo_press all pressure thermo_temp compute thermo_press all pressure thermo_temp

133
doc/src/fix_langevin.rst
View File

@ -51,7 +51,7 @@ Examples
"""""""" """"""""
.. parsed-literal:: .. code-block:: LAMMPS
fix 3 boundary langevin 1.0 1.0 1000.0 699483 fix 3 boundary langevin 1.0 1.0 1000.0 699483
fix 1 all langevin 1.0 1.1 100.0 48279 scale 3 1.5 fix 1 all langevin 1.0 1.1 100.0 48279 scale 3 1.5
@ -67,35 +67,35 @@ performs Brownian dynamics (BD), since the total force on each atom
will have the form: will have the form:
.. parsed-literal:: .. math::
F = Fc + Ff + Fr F = & F_c + F_f + F_r \\
Ff = - (m / damp) v F_f = & - \frac{m}{\mathrm{damp}} v \\
Fr is proportional to sqrt(Kb T m / (dt damp)) F_r \propto & \sqrt{\frac{k_B T m}{dt~\mathrm{damp}}}
Fc is the conservative force computed via the usual inter-particle :math:`F_c` is the conservative force computed via the usual
interactions (:doc:`pair_style <pair_style>`, inter-particle interactions (:doc:`pair_style <pair_style>`,
:doc:`bond_style <bond_style>`, etc). :doc:`bond_style <bond_style>`, etc). The :math:`F_f` and :math:`F_r`
terms are added by this fix on a per-particle basis. See the
:doc:`pair_style dpd/tstat <pair_dpd>` command for a thermostatting
option that adds similar terms on a pairwise basis to pairs of
interacting particles.
The Ff and Fr terms are added by this fix on a per-particle basis. :math:`F_f` is a frictional drag or viscous damping term proportional to
See the :doc:`pair_style dpd/tstat <pair_dpd>` command for a the particle's velocity. The proportionality constant for each atom is
thermostatting option that adds similar terms on a pairwise basis to computed as :math:`\frac{m}{\mathrm{damp}}`, where *m* is the mass of the
pairs of interacting particles. particle and damp is the damping factor specified by the user.
Ff is a frictional drag or viscous damping term proportional to the :math:`F_r` is a force due to solvent atoms at a temperature *T*
particle's velocity. The proportionality constant for each atom is randomly bumping into the particle. As derived from the
computed as m/damp, where m is the mass of the particle and damp is fluctuation/dissipation theorem, its magnitude as shown above is
the damping factor specified by the user. proportional to :math:`\sqrt{\frac{k_B T m}{dt~\mathrm{damp}}}`, where
:math:`k_B` is the Boltzmann constant, *T* is the desired temperature,
Fr is a force due to solvent atoms at a temperature T randomly bumping *m* is the mass of the particle, *dt* is the timestep size, and damp is
into the particle. As derived from the fluctuation/dissipation the damping factor. Random numbers are used to randomize the direction
theorem, its magnitude as shown above is proportional to sqrt(Kb T m / and magnitude of this force as described in :ref:`(Dunweg) <Dunweg1>`,
dt damp), where Kb is the Boltzmann constant, T is the desired where a uniform random number is used (instead of a Gaussian random
temperature, m is the mass of the particle, dt is the timestep size, number) for speed.
and damp is the damping factor. Random numbers are used to randomize
the direction and magnitude of this force as described in
:ref:`(Dunweg) <Dunweg1>`, where a uniform random number is used (instead of
a Gaussian random number) for speed.
Note that unless you use the *omega* or *angmom* keywords, the Note that unless you use the *omega* or *angmom* keywords, the
thermostat effect of this fix is applied to only the translational thermostat effect of this fix is applied to only the translational
@ -107,14 +107,15 @@ thermostatting takes place; see the description below.
.. note:: .. note::
Unlike the :doc:`fix nvt <fix_nh>` command which performs Unlike the :doc:`fix nvt <fix_nh>` command which performs Nose/Hoover
Nose/Hoover thermostatting AND time integration, this fix does NOT thermostatting AND time integration, this fix does NOT perform time
perform time integration. It only modifies forces to effect integration. It only modifies forces to effect thermostatting. Thus
thermostatting. Thus you must use a separate time integration fix, you must use a separate time integration fix, like :doc:`fix nve
like :doc:`fix nve <fix_nve>` to actually update the velocities and <fix_nve>` to actually update the velocities and positions of atoms
positions of atoms using the modified forces. Likewise, this fix using the modified forces. Likewise, this fix should not normally be
should not normally be used on atoms that also have their temperature used on atoms that also have their temperature controlled by another
controlled by another fix - e.g. by :doc:`fix nvt <fix_nh>` or :doc:`fix temp/rescale <fix_temp_rescale>` commands. fix - e.g. by :doc:`fix nvt <fix_nh>` or :doc:`fix temp/rescale
<fix_temp_rescale>` commands.
See the :doc:`Howto thermostat <Howto_thermostat>` doc page for See the :doc:`Howto thermostat <Howto_thermostat>` doc page for
a discussion of different ways to compute temperature and perform a discussion of different ways to compute temperature and perform
@ -154,13 +155,14 @@ the remaining thermal degrees of freedom, and the bias is added back
in. in.
The *damp* parameter is specified in time units and determines how The *damp* parameter is specified in time units and determines how
rapidly the temperature is relaxed. For example, a value of 100.0 rapidly the temperature is relaxed. For example, a value of 100.0 means
means to relax the temperature in a timespan of (roughly) 100 time to relax the temperature in a timespan of (roughly) 100 time units
units (tau or fmsec or psec - see the :doc:`units <units>` command). (:math:`\tau` or fs or ps - see the :doc:`units <units>` command). The
The damp factor can be thought of as inversely related to the damp factor can be thought of as inversely related to the viscosity of
viscosity of the solvent. I.e. a small relaxation time implies a the solvent. I.e. a small relaxation time implies a high-viscosity
hi-viscosity solvent and vice versa. See the discussion about gamma solvent and vice versa. See the discussion about :math:`\gamma` and
and viscosity in the documentation for the :doc:`fix viscous <fix_viscous>` command for more details. viscosity in the documentation for the :doc:`fix viscous <fix_viscous>`
command for more details.
The random # *seed* must be a positive integer. A Marsaglia random The random # *seed* must be a positive integer. A Marsaglia random
number generator is used. Each processor uses the input seed to number generator is used. Each processor uses the input seed to
@ -191,39 +193,40 @@ The rotational temperature of the particles can be monitored by the
:doc:`compute temp/sphere <compute_temp_sphere>` and :doc:`compute temp/asphere <compute_temp_asphere>` commands with their rotate :doc:`compute temp/sphere <compute_temp_sphere>` and :doc:`compute temp/asphere <compute_temp_asphere>` commands with their rotate
options. options.
For the *omega* keyword there is also a scale factor of 10.0/3.0 that For the *omega* keyword there is also a scale factor of
is applied as a multiplier on the Ff (damping) term in the equation :math:`\frac{10.0}{3.0}` that is applied as a multiplier on the
above and of sqrt(10.0/3.0) as a multiplier on the Fr term. This does :math:`F_f` (damping) term in the equation above and of
not affect the thermostatting behavior of the Langevin formalism but :math:`\sqrt{\frac{10.0}{3.0}}` as a multiplier on the :math:`F_r` term.
insures that the randomized rotational diffusivity of spherical This does not affect the thermostatting behavior of the Langevin
particles is correct. formalism but insures that the randomized rotational diffusivity of
spherical particles is correct.
For the *angmom* keyword a similar scale factor is needed which is For the *angmom* keyword a similar scale factor is needed which is
10.0/3.0 for spherical particles, but is anisotropic for aspherical :math:`\frac{10.0}{3.0}` for spherical particles, but is anisotropic for
particles (e.g. ellipsoids). Currently LAMMPS only applies an aspherical particles (e.g. ellipsoids). Currently LAMMPS only applies
isotropic scale factor, and you can choose its magnitude as the an isotropic scale factor, and you can choose its magnitude as the
specified value of the *angmom* keyword. If your aspherical particles specified value of the *angmom* keyword. If your aspherical particles
are (nearly) spherical than a value of 10.0/3.0 = 3.333 is a good are (nearly) spherical than a value of :math:`\frac{10.0}{3.0} =
choice. If they are highly aspherical, a value of 1.0 is as good a 3.\overline{3}` is a good choice. If they are highly aspherical, a
choice as any, since the effects on rotational diffusivity of the value of 1.0 is as good a choice as any, since the effects on rotational
particles will be incorrect regardless. Note that for any reasonable diffusivity of the particles will be incorrect regardless. Note that
scale factor, the thermostatting effect of the *angmom* keyword on the for any reasonable scale factor, the thermostatting effect of the
rotational temperature of the aspherical particles should still be *angmom* keyword on the rotational temperature of the aspherical
valid. particles should still be valid.
The keyword *scale* allows the damp factor to be scaled up or down by The keyword *scale* allows the damp factor to be scaled up or down by
the specified factor for atoms of that type. This can be useful when the specified factor for atoms of that type. This can be useful when
different atom types have different sizes or masses. It can be used different atom types have different sizes or masses. It can be used
multiple times to adjust damp for several atom types. Note that multiple times to adjust damp for several atom types. Note that
specifying a ratio of 2 increases the relaxation time which is specifying a ratio of 2 increases the relaxation time which is
equivalent to the solvent's viscosity acting on particles with 1/2 the equivalent to the solvent's viscosity acting on particles with
diameter. This is the opposite effect of scale factors used by the :math:`\frac{1}{2}` the diameter. This is the opposite effect of scale
:doc:`fix viscous <fix_viscous>` command, since the damp factor in fix factors used by the :doc:`fix viscous <fix_viscous>` command, since the
*langevin* is inversely related to the gamma factor in fix *viscous*\ . damp factor in fix *langevin* is inversely related to the :math:`\gamma`
Also note that the damping factor in fix *langevin* includes the factor in fix *viscous*\ . Also note that the damping factor in fix
particle mass in Ff, unlike fix *viscous*\ . Thus the mass and size of *langevin* includes the particle mass in Ff, unlike fix *viscous*\ .
different atom types should be accounted for in the choice of ratio Thus the mass and size of different atom types should be accounted for
values. in the choice of ratio values.
The keyword *tally* enables the calculation of the cumulative energy The keyword *tally* enables the calculation of the cumulative energy
added/subtracted to the atoms as they are thermostatted. Effectively added/subtracted to the atoms as they are thermostatted. Effectively

35
doc/src/fix_langevin_eff.rst
View File

@ -41,7 +41,7 @@ Examples
"""""""" """"""""
.. parsed-literal:: .. code-block:: LAMMPS
fix 3 boundary langevin/eff 1.0 1.0 10.0 699483 fix 3 boundary langevin/eff 1.0 1.0 10.0 699483
fix 1 all langevin/eff 1.0 1.1 10.0 48279 scale 3 1.5 fix 1 all langevin/eff 1.0 1.1 10.0 48279 scale 3 1.5
@ -55,28 +55,31 @@ this command performs Brownian dynamics (BD), since the total force on
each atom will have the form: each atom will have the form:
.. parsed-literal:: .. math::
F = Fc + Ff + Fr F = & F_c + F_f + F_r \\
Ff = - (m / damp) v F_f = & - \frac{m}{\mathrm{damp}} v \\
Fr is proportional to sqrt(Kb T m / (dt damp)) F_r \propto & \sqrt{\frac{k_B T m}{dt~\mathrm{damp}}}
Fc is the conservative force computed via the usual inter-particle
interactions (:doc:`pair_style <pair_style>`).
The Ff and Fr terms are added by this fix on a per-particle basis. :math:`F_c` is the conservative force computed via the usual
inter-particle interactions (:doc:`pair_style <pair_style>`).
The :math:`F_f` and :math:`F_r` terms are added by this fix on a
per-particle basis.
The operation of this fix is exactly like that described by the :doc:`fix langevin <fix_langevin>` command, except that the thermostatting The operation of this fix is exactly like that described by the
is also applied to the radial electron velocity for electron :doc:`fix langevin <fix_langevin>` command, except that the
particles. thermostatting is also applied to the radial electron velocity for
electron particles.
**Restart, fix\_modify, output, run start/stop, minimize info:** **Restart, fix\_modify, output, run start/stop, minimize info:**
No information about this fix is written to :doc:`binary restart files <restart>`. Because the state of the random number generator No information about this fix is written to :doc:`binary restart files
is not saved in restart files, this means you cannot do "exact" <restart>`. Because the state of the random number generator is not
restarts with this fix, where the simulation continues on the same as saved in restart files, this means you cannot do "exact" restarts with
if no restart had taken place. However, in a statistical sense, a this fix, where the simulation continues on the same as if no restart
restarted simulation should produce the same behavior. had taken place. However, in a statistical sense, a restarted
simulation should produce the same behavior.
The :doc:`fix_modify <fix_modify>` *temp* option is supported by this The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
fix. You can use it to assign a temperature :doc:`compute <compute>` fix. You can use it to assign a temperature :doc:`compute <compute>`

56
doc/src/fix_nve_dotc_langevin.rst
View File

@ -29,7 +29,7 @@ Examples
"""""""" """"""""
.. parsed-literal:: .. code-block:: LAMMPS
fix 1 all nve/dotc/langevin 1.0 1.0 0.03 457145 angmom 10 fix 1 all nve/dotc/langevin 1.0 1.0 0.03 457145 angmom 10
fix 1 all nve/dotc/langevin 0.1 0.1 78.9375 457145 angmom 10 fix 1 all nve/dotc/langevin 0.1 0.1 78.9375 457145 angmom 10
@ -61,33 +61,33 @@ over the standard integrator, permitting slightly larger timestep sizes.
The total force on each atom will have the form: The total force on each atom will have the form:
.. parsed-literal:: .. math::
F = Fc + Ff + Fr F = & F_c + F_f + F_r \\
Ff = - (m / damp) v F_f = & - \frac{m}{\mathrm{damp}} v \\
Fr is proportional to sqrt(Kb T m / (dt damp)) F_r \propto & \sqrt{\frac{k_B T m}{dt~\mathrm{damp}}}
Fc is the conservative force computed via the usual inter-particle :math:`F_c` is the conservative force computed via the usual
interactions (:doc:`pair_style <pair_style>`, inter-particle interactions (:doc:`pair_style <pair_style>`,
:doc:`bond_style <bond_style>`, etc). :doc:`bond_style <bond_style>`, etc). The :math:`F_f` and :math:`F_r`
terms are implicitly taken into account by this fix on a per-particle
basis.
The Ff and Fr terms are implicitly taken into account by this fix :math:`F_f` is a frictional drag or viscous damping term proportional to
on a per-particle basis. the particle's velocity. The proportionality constant for each atom is
computed as :math:`\frac{m}{\mathrm{damp}}`, where *m* is the mass of
the particle and damp is the damping factor specified by the user.
Ff is a frictional drag or viscous damping term proportional to the :math:`F_r` is a force due to solvent atoms at a temperature *T*
particle's velocity. The proportionality constant for each atom is randomly bumping into the particle. As derived from the
computed as m/damp, where m is the mass of the particle and damp is fluctuation/dissipation theorem, its magnitude as shown above is
the damping factor specified by the user. proportional to :math:`\sqrt{\frac{k_B T m}{dt~\mathrm{damp}}}`, where
:math:`k_B` is the Boltzmann constant, *T* is the desired temperature,
Fr is a force due to solvent atoms at a temperature T randomly bumping *m* is the mass of the particle, *dt* is the timestep size, and damp is
into the particle. As derived from the fluctuation/dissipation the damping factor. Random numbers are used to randomize the direction
theorem, its magnitude as shown above is proportional to sqrt(Kb T m / and magnitude of this force as described in :ref:`(Dunweg) <Dunweg5>`,
dt damp), where Kb is the Boltzmann constant, T is the desired where a uniform random number is used (instead of a Gaussian random
temperature, m is the mass of the particle, dt is the timestep size, number) for speed.
and damp is the damping factor. Random numbers are used to randomize
the direction and magnitude of this force as described in
:ref:`(Dunweg) <Dunweg5>`, where a uniform random number is used (instead of
a Gaussian random number) for speed.
---------- ----------
@ -104,7 +104,7 @@ means to relax the temperature in a timespan of (roughly) 0.03 time
units tau (see the :doc:`units <units>` command). units tau (see the :doc:`units <units>` command).
The damp factor can be thought of as inversely related to the The damp factor can be thought of as inversely related to the
viscosity of the solvent, i.e. a small relaxation time implies a viscosity of the solvent, i.e. a small relaxation time implies a
hi-viscosity solvent and vice versa. See the discussion about gamma high-viscosity solvent and vice versa. See the discussion about gamma
and viscosity in the documentation for the :doc:`fix viscous <fix_viscous>` command for more details. and viscosity in the documentation for the :doc:`fix viscous <fix_viscous>` command for more details.
Note that the value 78.9375 in the second example above corresponds Note that the value 78.9375 in the second example above corresponds
to a diffusion constant, which is about an order of magnitude larger to a diffusion constant, which is about an order of magnitude larger
@ -124,10 +124,10 @@ particles are always considered to have a finite size.
The keyword *angmom* enables thermostatting of the rotational degrees of The keyword *angmom* enables thermostatting of the rotational degrees of
freedom in addition to the usual translational degrees of freedom. freedom in addition to the usual translational degrees of freedom.
The scale factor after the *angmom* keyword gives the ratio of the rotational to The scale factor after the *angmom* keyword gives the ratio of the
the translational friction coefficient. rotational to the translational friction coefficient.
An example input file can be found in /examples/USER/cgdna/examples/duplex2/. An example input file can be found in examples/USER/cgdna/examples/duplex2/.
Further details of the implementation and stability of the integrators are contained in :ref:`(Henrich) <Henrich5>`. Further details of the implementation and stability of the integrators are contained in :ref:`(Henrich) <Henrich5>`.
The preprint version of the article can be found `here <PDF/USER-CGDNA.pdf>`_. The preprint version of the article can be found `here <PDF/USER-CGDNA.pdf>`_.

119
doc/src/fix_srd.rst
View File

@ -84,36 +84,37 @@ the implementation and usage of mixture systems (solute particles in
an SRD fluid). See the examples/srd directory for sample input an SRD fluid). See the examples/srd directory for sample input
scripts using SRD particles in both settings. scripts using SRD particles in both settings.
This fix does 2 things: This fix does two things:
(1) It advects the SRD particles, performing collisions between SRD 1. It advects the SRD particles, performing collisions between SRD
and big particles or walls every timestep, imparting force and torque and big particles or walls every timestep, imparting force and torque
to the big particles. Collisions also change the position and to the big particles. Collisions also change the position and
velocity of SRD particles. velocity of SRD particles.
(2) It resets the velocity distribution of SRD particles via random 2. It resets the velocity distribution of SRD particles via random
rotations every N timesteps. rotations every N timesteps.
SRD particles have a mass, temperature, characteristic timestep SRD particles have a mass, temperature, characteristic timestep
dt\_SRD, and mean free path between collisions (lamda). The :math:`dt_{SRD}`, and mean free path between collisions
fundamental equation relating these 4 quantities is (:math:`\lambda`). The fundamental equation relating these 4 quantities
is
.. parsed-literal:: .. math::
lamda = dt_SRD \* sqrt(Kboltz \* Tsrd / mass) \lambda = dt_{SRD} \sqrt{\frac{k_B T_{SRD}}{m}}
The mass of SRD particles is set by the :doc:`mass <mass>` command The mass *m* of SRD particles is set by the :doc:`mass <mass>` command
elsewhere in the input script. The SRD timestep dt\_SRD is N times the elsewhere in the input script. The SRD timestep :math:`dt_{SRD}` is N
step dt defined by the :doc:`timestep <timestep>` command. Big times the step *dt* defined by the :doc:`timestep <timestep>` command.
particles move in the normal way via a time integration :doc:`fix <fix>` Big particles move in the normal way via a time integration :doc:`fix
with a short timestep dt. SRD particles advect with a large timestep <fix>` with a short timestep dt. SRD particles advect with a large
dt\_SRD >= dt. timestep :math:`dt_{SRD} \ge dt`.
If the *lamda* keyword is not specified, the SRD temperature If the *lamda* keyword is not specified, the SRD temperature
*Tsrd* is used in the above formula to compute lamda. If the *lamda* :math:`T_{SRD}` is used in the above formula to compute :math:`\lambda`.
keyword is specified, then the *Tsrd* setting is ignored and the above If the *lamda* keyword is specified, then the *Tsrd* setting is ignored
equation is used to compute the SRD temperature. and the above equation is used to compute the SRD temperature.
The characteristic length scale for the SRD fluid is set by *hgrid* The characteristic length scale for the SRD fluid is set by *hgrid*
which is used to bin SRD particles for purposes of resetting their which is used to bin SRD particles for purposes of resetting their
@ -121,13 +122,14 @@ velocities. Normally hgrid is set to be 1/4 of the big particle
diameter or smaller, to adequately resolve fluid properties around the diameter or smaller, to adequately resolve fluid properties around the
big particles. big particles.
Lamda cannot be smaller than 0.6 \* hgrid, else an error is generated :math:`\lambda` cannot be smaller than 0.6 \* hgrid, else an error is
(unless the *shift* keyword is used, see below). The velocities of generated (unless the *shift* keyword is used, see below). The
SRD particles are bounded by Vmax, which is set so that an SRD velocities of SRD particles are bounded by Vmax, which is set so that an
particle will not advect further than Dmax = 4\*lamda in dt\_SRD. This SRD particle will not advect further than :math:`D_{max} = 4 \lambda` in
means that roughly speaking, Dmax should not be larger than a big :math:`dt_{SRD}`. This means that roughly speaking, :math:`D_{max}`
particle diameter, else SRDs may pass through big particles without should not be larger than a big particle diameter, else SRDs may pass
colliding. A warning is generated if this is the case. through big particles without colliding. A warning is generated if this
is the case.
Collisions between SRD particles and big particles or walls are Collisions between SRD particles and big particles or walls are
modeled as a lightweight SRD point particle hitting a heavy big modeled as a lightweight SRD point particle hitting a heavy big
@ -160,11 +162,12 @@ the big particles appropriately should be used.
The *overlap* keyword should be set to *yes* if two (or more) big The *overlap* keyword should be set to *yes* if two (or more) big
particles can ever overlap. This depends on the pair potential particles can ever overlap. This depends on the pair potential
interaction used for big-big interactions, or could be the case if interaction used for big-big interactions, or could be the case if
multiple big particles are held together as rigid bodies via the :doc:`fix rigid <fix_rigid>` command. If the *overlap* keyword is *no* and multiple big particles are held together as rigid bodies via the
big particles do in fact overlap, then SRD/big collisions can generate :doc:`fix rigid <fix_rigid>` command. If the *overlap* keyword is *no*
an error if an SRD ends up inside two (or more) big particles at once. and big particles do in fact overlap, then SRD/big collisions can
How this error is treated is determined by the *inside* keyword. generate an error if an SRD ends up inside two (or more) big particles
Running with *overlap* set to *no* allows for faster collision at once. How this error is treated is determined by the *inside*
keyword. Running with *overlap* set to *no* allows for faster collision
checking, so it should only be set to *yes* if needed. checking, so it should only be set to *yes* if needed.
The *inside* keyword determines how a collision is treated if the The *inside* keyword determines how a collision is treated if the
@ -258,30 +261,30 @@ warning is generated.
reneighboring. Note that changing the SRD bin size may alter the reneighboring. Note that changing the SRD bin size may alter the
properties of the SRD fluid, such as its viscosity. properties of the SRD fluid, such as its viscosity.
The *shift* keyword determines whether the coordinates of SRD The *shift* keyword determines whether the coordinates of SRD particles
particles are randomly shifted when binned for purposes of rotating are randomly shifted when binned for purposes of rotating their
their velocities. When no shifting is performed, SRD particles are velocities. When no shifting is performed, SRD particles are binned and
binned and the velocity distribution of the set of SRD particles in the velocity distribution of the set of SRD particles in each bin is
each bin is adjusted via a rotation operator. This is a statistically adjusted via a rotation operator. This is a statistically valid
valid operation if SRD particles move sufficiently far between operation if SRD particles move sufficiently far between successive
successive rotations. This is determined by their mean-free path rotations. This is determined by their mean-free path :math:`\lambda`.
lamda. If lamda is less than 0.6 of the SRD bin size, then shifting If :math:`\lambda` is less than 0.6 of the SRD bin size, then shifting
is required. A shift means that all of the SRD particles are shifted is required. A shift means that all of the SRD particles are shifted by
by a vector whose coordinates are chosen randomly in the range [-1/2 a vector whose coordinates are chosen randomly in the range [-1/2 bin
bin size, 1/2 bin size]. Note that all particles are shifted by the size, 1/2 bin size]. Note that all particles are shifted by the same
same vector. The specified random number *shiftseed* is used to vector. The specified random number *shiftseed* is used to generate
generate these vectors. This operation sufficiently randomizes which these vectors. This operation sufficiently randomizes which SRD
SRD particles are in the same bin, even if lamda is small. particles are in the same bin, even if :math:`lambda` is small.
If the *shift* flag is set to *no*\ , then no shifting is performed, but If the *shift* flag is set to *no*\ , then no shifting is performed, but
bin data will be communicated if bins overlap processor boundaries. bin data will be communicated if bins overlap processor boundaries. An
An error will be generated if lamda < 0.6 of the SRD bin size. If the error will be generated if :math:`\lambda < 0.6` of the SRD bin size.
*shift* flag is set to *possible*\ , then shifting is performed only if If the *shift* flag is set to *possible*\ , then shifting is performed
lamda < 0.6 of the SRD bin size. A warning is generated to let you only if :math:`\lambda < 0.6` of the SRD bin size. A warning is
know this is occurring. If the *shift* flag is set to *yes* then generated to let you know this is occurring. If the *shift* flag is set
shifting is performed regardless of the magnitude of lamda. Note that to *yes* then shifting is performed regardless of the magnitude of
the *shiftseed* is not used if the *shift* flag is set to *no*\ , but :math:`\lambda`. Note that the *shiftseed* is not used if the *shift*
must still be specified. flag is set to *no*\ , but must still be specified.
Note that shifting of SRD coordinates requires extra communication, Note that shifting of SRD coordinates requires extra communication,
hence it should not normally be enabled unless required. hence it should not normally be enabled unless required.
@ -398,10 +401,10 @@ Related commands
Default Default
""""""" """""""
The option defaults are lamda inferred from Tsrd, collision = noslip, The option defaults are: *lamda* (:math:`\lambda`) is inferred from *Tsrd*,
overlap = no, inside = error, exact = yes, radius = 1.0, bounce = 0, collision = noslip, overlap = no, inside = error, exact = yes, radius =
search = hgrid, cubic = error 0.01, shift = no, tstat = no, and 1.0, bounce = 0, search = hgrid, cubic = error 0.01, shift = no, tstat =
rescale = yes. no, and rescale = yes.
---------- ----------

56
doc/src/fix_viscous.rst
View File

@ -47,11 +47,11 @@ energy to the system), it has the effect of slowly (or rapidly)
freezing the system; hence it can also be used as a simple energy freezing the system; hence it can also be used as a simple energy
minimization technique. minimization technique.
The damping force F is given by F = - gamma \* velocity. The larger The damping force :math:`F_i` is given by :math:`F_i = - \gamma v_i`.
the coefficient, the faster the kinetic energy is reduced. If the The larger the coefficient, the faster the kinetic energy is reduced.
optional keyword *scale* is used, gamma can scaled up or down by the If the optional keyword *scale* is used, :math:`\gamma` can scaled up or
specified factor for atoms of that type. It can be used multiple down by the specified factor for atoms of that type. It can be used
times to adjust gamma for several atom types. multiple times to adjust :math:`\gamma` for several atom types.
.. note:: .. note::
@ -60,38 +60,40 @@ times to adjust gamma for several atom types.
:doc:`units <units>` options like "real" or "metal" that are not :doc:`units <units>` options like "real" or "metal" that are not
self-consistent. self-consistent.
In a Brownian dynamics context, gamma = Kb T / D, where Kb = In a Brownian dynamics context, :math:`\gamma = \frac{k_B T}{D}`, where
Boltzmann's constant, T = temperature, and D = particle diffusion :math:`k_B =` Boltzmann's constant, *T* = temperature, and *D* = particle
coefficient. D can be written as Kb T / (3 pi eta d), where eta = diffusion coefficient. *D* can be written as :math:`\frac{k_B T}{3 \pi
dynamic viscosity of the frictional fluid and d = diameter of \eta d}`, where :math:`\eta =` dynamic viscosity of the frictional fluid
particle. This means gamma = 3 pi eta d, and thus is proportional to and d = diameter of particle. This means :math:`\gamma = 3 \pi \eta d`,
the viscosity of the fluid and the particle diameter. and thus is proportional to the viscosity of the fluid and the particle
diameter.
In the current implementation, rather than have the user specify a In the current implementation, rather than have the user specify a
viscosity, gamma is specified directly in force/velocity units. If viscosity, :math:`\gamma` is specified directly in force/velocity units.
needed, gamma can be adjusted for atoms of different sizes If needed, :math:`\gamma` can be adjusted for atoms of different sizes
(i.e. sigma) by using the *scale* keyword. (i.e. :math:`\sigma`) by using the *scale* keyword.
Note that Brownian dynamics models also typically include a randomized Note that Brownian dynamics models also typically include a randomized
force term to thermostat the system at a chosen temperature. The :doc:`fix langevin <fix_langevin>` command does this. It has the same force term to thermostat the system at a chosen temperature. The
:doc:`fix langevin <fix_langevin>` command does this. It has the same
viscous damping term as fix viscous and adds a random force to each viscous damping term as fix viscous and adds a random force to each
atom. The random force term is proportional to the sqrt of the chosen atom. The random force term is proportional to the square root of the
thermostatting temperature. Thus if you use fix langevin with a chosen thermostatting temperature. Thus if you use fix langevin with a
target T = 0, its random force term is zero, and you are essentially target :math:`T = 0`, its random force term is zero, and you are
performing the same operation as fix viscous. Also note that the essentially performing the same operation as fix viscous. Also note
gamma of fix viscous is related to the damping parameter of :doc:`fix langevin <fix_langevin>`, however the former is specified in units that the gamma of fix viscous is related to the damping parameter of
of force/velocity and the latter in units of time, so that it can more :doc:`fix langevin <fix_langevin>`, however the former is specified in
easily be used as a thermostat. units of force/velocity and the latter in units of time, so that it can
more easily be used as a thermostat.
---------- ----------
**Restart, fix\_modify, output, run start/stop, minimize info:** **Restart, fix\_modify, output, run start/stop, minimize info:**
No information about this fix is written to :doc:`binary restart files <restart>`. None of the :doc:`fix_modify <fix_modify>` options No information about this fix is written to :doc:`binary restart files
are relevant to this fix. No global or per-atom quantities are stored <restart>`. None of the :doc:`fix_modify <fix_modify>` options are
by this fix for access by various :doc:`output commands <Howto_output>`. relevant to this fix. No global or per-atom quantities are stored by
this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. the :doc:`run <run>` command.

20
doc/src/thermo_style.rst
View File

@ -50,7 +50,7 @@ Syntax
ke = kinetic energy ke = kinetic energy
etotal = total energy (pe + ke) etotal = total energy (pe + ke)
enthalpy = enthalpy (etotal + press\*vol) enthalpy = enthalpy (etotal + press\*vol)
evdwl = VanderWaal pairwise energy (includes etail) evdwl = van der Waals pairwise energy (includes etail)
ecoul = Coulombic pairwise energy ecoul = Coulombic pairwise energy
epair = pairwise energy (evdwl + ecoul + elong) epair = pairwise energy (evdwl + ecoul + elong)
ebond = bond energy ebond = bond energy
@ -59,7 +59,7 @@ Syntax
eimp = improper energy eimp = improper energy
emol = molecular energy (ebond + eangle + edihed + eimp) emol = molecular energy (ebond + eangle + edihed + eimp)
elong = long-range kspace energy elong = long-range kspace energy
etail = VanderWaal energy long-range tail correction etail = van der Waals energy long-range tail correction
vol = volume vol = volume
density = mass density of system density = mass density of system
lx,ly,lz = box lengths in x,y,z lx,ly,lz = box lengths in x,y,z
@ -205,13 +205,13 @@ change the attributes of this potential energy via the
The kinetic energy of the system *ke* is inferred from the temperature The kinetic energy of the system *ke* is inferred from the temperature
of the system with 1/2 Kb T of energy for each degree of freedom. of the system with :math:`\frac{1}{2} k_B T` of energy for each degree
Thus, using different :doc:`compute commands <compute>` for calculating of freedom. Thus, using different :doc:`compute commands <compute>` for
temperature, via the :doc:`thermo_modify temp <thermo_modify>` command, calculating temperature, via the :doc:`thermo_modify temp
may yield different kinetic energies, since different computes that <thermo_modify>` command, may yield different kinetic energies, since
calculate temperature can subtract out different non-thermal different computes that calculate temperature can subtract out different
components of velocity and/or include different degrees of freedom non-thermal components of velocity and/or include different degrees of
(translational, rotational, etc). freedom (translational, rotational, etc).
The potential energy of the system *pe* will include contributions The potential energy of the system *pe* will include contributions
from fixes if the :doc:`fix_modify thermo <fix_modify>` option is set from fixes if the :doc:`fix_modify thermo <fix_modify>` option is set
@ -219,7 +219,7 @@ for a fix that calculates such a contribution. For example, the :doc:`fix wall/
interacting with the wall. See the doc pages for "individual fixes" interacting with the wall. See the doc pages for "individual fixes"
to see which ones contribute. to see which ones contribute.
A long-range tail correction *etail* for the VanderWaal pairwise A long-range tail correction *etail* for the van der Waals pairwise
energy will be non-zero only if the :doc:`pair_modify tail <pair_modify>` option is turned on. The *etail* contribution energy will be non-zero only if the :doc:`pair_modify tail <pair_modify>` option is turned on. The *etail* contribution
is included in *evdwl*\ , *epair*\ , *pe*\ , and *etotal*\ , and the is included in *evdwl*\ , *epair*\ , *pe*\ , and *etotal*\ , and the
corresponding tail correction to the pressure is included in *press* corresponding tail correction to the pressure is included in *press*

14
doc/src/units.rst
View File

@ -57,13 +57,13 @@ is often not simple to do.
For style *lj*\ , all quantities are unitless. Without loss of For style *lj*\ , all quantities are unitless. Without loss of
generality, LAMMPS sets the fundamental quantities mass, :math:`\sigma`, generality, LAMMPS sets the fundamental quantities mass, :math:`\sigma`,
:math:`\epsilon`, and the Boltzmann constant :math:`k_B = 1`. :math:`\epsilon`, and the Boltzmann constant :math:`k_B = 1`. The
The masses, distances, masses, distances, energies you specify are multiples of these
energies you specify are multiples of these fundamental values. The fundamental values. The formulas relating the reduced or unitless
formulas relating the reduced or unitless quantity (with an asterisk) quantity (with an asterisk) to the same quantity with units is also
to the same quantity with units is also given. Thus you can use the given. Thus you can use the mass & :math:`\sigma` & :math:`\epsilon`
mass & sigma & epsilon values for a specific material and convert the values for a specific material and convert the results from a unitless
results from a unitless LJ simulation into physical quantities. LJ simulation into physical quantities.
* mass = mass or *m* * mass = mass or *m*
* distance = :math:`\sigma`, where :math:`x^* = \frac{x}{\sigma}` * distance = :math:`\sigma`, where :math:`x^* = \frac{x}{\sigma}`