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

Merge remote-tracking branch 'github/develop' into collected-small-fixes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-08-26 11:37:56 -04:00
parent d9a2fd9f36 121195d19d
commit 4e590507f0
26 changed files with 1170 additions and 1102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
doc/src/compute_pressure.rst
View File

@ -59,11 +59,12 @@ may also contribute to the virial term.
A symmetric pressure tensor, stored as a 6-element vector, is also
calculated by this compute. The six components of the vector are
ordered :math:`xx,` :math:`yy,` :math:`zz,` :math:`xy,` :math:`xz,` :math:`yz.`
The equation for the :math:`(I,J)` components (where :math:`I` and :math:`J`
are :math:`x`, :math:`y`, or :math:`z`) is similar to the above formula,
except that the first term uses components of the kinetic energy tensor and the
second term uses components of the virial tensor:
ordered :math:`xx,` :math:`yy,` :math:`zz,` :math:`xy,` :math:`xz,`
:math:`yz.` The equation for the :math:`(I,J)` components (where
:math:`I` and :math:`J` are :math:`x`, :math:`y`, or :math:`z`) is
similar to the above formula, except that the first term uses
components related to the kinetic energy tensor and the second term
uses components of the virial tensor:
.. math::
@ -75,8 +76,8 @@ calculated. This includes a kinetic energy (temperature) term and the
virial as the sum of pair, bond, angle, dihedral, improper, kspace
(long-range), and fix contributions to the force on each atom. If any
extra keywords are listed, then only those components are summed to
compute temperature or ke and/or the virial. The *virial* keyword
means include all terms except the kinetic energy *ke*\ .
compute temperature or ke and/or the virial. The *virial* keyword means
include all terms except the kinetic energy *ke*\ .
The *pair/hybrid* keyword means to only include contribution
from a sub-style in a *hybrid* or *hybrid/overlay* pair style.
@ -86,26 +87,31 @@ system, including for many-body potentials and accounting for the
effects of periodic boundary conditions are discussed in
:ref:`(Thompson) <Thompson1>`.
The temperature and kinetic energy tensor is not calculated by this
The temperature and kinetic energy tensor are not calculated by this
compute, but rather by the temperature compute specified with the
command. If the kinetic energy is not included in the pressure, than
the temperature compute is not used and can be specified as NULL.
Normally the temperature compute used by compute pressure should
calculate the temperature of all atoms for consistency with the virial
term, but any compute style that calculates temperature can be used
(e.g., one that excludes frozen atoms or other degrees of freedom).
command. See the doc pages for individual compute temp variants for an
explanation of how they calculate temperature and a symmetric tensor
(6-element vector) whose components are twice that of the traditional KE
tensor. That tensor is what appears in the pressure tensor formula
above.
If the kinetic energy is not included in the pressure, than the
temperature compute is not used and can be specified as NULL. Normally
the temperature compute used by compute pressure should calculate the
temperature of all atoms for consistency with the virial term, but any
compute style that calculates temperature can be used (e.g., one that
excludes frozen atoms or other degrees of freedom).
Note that if desired the specified temperature compute can be one that
subtracts off a bias to calculate a temperature using only the thermal
velocity of the atoms (e.g., by subtracting a background streaming
velocity).
See the doc pages for individual :doc:`compute commands <compute>` to determine
which ones include a bias.
velocity). See the doc pages for individual :doc:`compute commands
<compute>` to determine which ones include a bias.
Also note that the :math:`N` in the first formula above is really
degrees-of-freedom divided by :math:`d` = dimensionality, where the DOF value
is calculated by the temperature compute.
See the various :doc:`compute temperature <compute>` styles for details.
degrees-of-freedom divided by :math:`d` = dimensionality, where the
DOF value is calculated by the temperature compute. See the various
:doc:`compute temperature <compute>` styles for details.
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:
@ -136,9 +142,8 @@ The ordering of values in the symmetric pressure tensor is as follows:
:math:`p_{xx},` :math:`p_{yy},` :math:`p_{zz},` :math:`p_{xy},`
:math:`p_{xz},` :math:`p_{yz}.`
The scalar and vector values calculated by this compute are
"intensive". The scalar and vector values will be in pressure
:doc:`units <units>`.
The scalar and vector values calculated by this compute are "intensive".
The scalar and vector values will be in pressure :doc:`units <units>`.
Restrictions
""""""""""""

33
doc/src/compute_temp.rst
View File

@ -48,14 +48,18 @@ the group, :math:`N_\mathrm{fix DOFs}` is the number of degrees of
freedom removed by fix commands (see below), :math:`k_B` is the
Boltzmann constant, and :math:`T` is the resulting computed temperature.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute for use in the computation of a pressure
tensor. The formula for the components of the tensor is the same as the
above expression for :math:`E_\mathrm{kin}`, except that :math:`v_i^2` is
replaced by :math:`v_{i,x} v_{i,y}` for the :math:`xy` component, and so on.
The six components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the *dynamic* option of the
:doc:`compute_modify <compute_modify>` command if this is not the case.
@ -94,16 +98,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length six (KE tensor), which can be accessed by indices
1--6. These values can be used by any command that uses global scalar
or vector values from a compute as input. See the :doc:`Howto output
<Howto_output>` page for an overview of LAMMPS output options.
vector of length six (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

63
doc/src/compute_temp_asphere.rst
View File

@ -41,8 +41,8 @@ translational and rotational kinetic energy. This differs from the
usual :doc:`compute temp <compute_temp>` command, which assumes point
particles with only translational kinetic energy.
Only finite-size particles (aspherical or spherical) can be included
in the group. For 3d finite-size particles, each has six degrees of
Only finite-size particles (aspherical or spherical) can be included in
the group. For 3d finite-size particles, each has six degrees of
freedom (three translational, three rotational). For 2d finite-size
particles, each has three degrees of freedom (two translational, one
rotational).
@ -70,25 +70,39 @@ axis. It will also be the case for biaxial ellipsoids when exactly two
of the semiaxes have the same length and the corresponding relative well
depths are equal.
The translational kinetic energy is computed the same as is described
by the :doc:`compute temp <compute_temp>` command. The rotational
kinetic energy is computed as :math:`\frac12 I \omega^2`, where :math:`I` is
the inertia tensor for the aspherical particle and :math:`\omega` is its
The translational kinetic energy is computed the same as is described by
the :doc:`compute temp <compute_temp>` command. The rotational kinetic
energy is computed as :math:`\frac12 I \omega^2`, where :math:`I` is the
inertia tensor for the aspherical particle and :math:`\omega` is its
angular velocity, which is computed from its angular momentum.
.. note::
For :doc:`2d models <dimension>`, particles are treated as
ellipsoids, not ellipses, meaning their moments of inertia will be the
same as in 3d.
ellipsoids, not ellipses, meaning their moments of inertia will be
the same as in 3d.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute. The formula for the components of the
tensor is the same as the above formula, except that :math:`v^2` and
:math:`\omega^2` are replaced by :math:`v_x v_y` and :math:`\omega_x \omega_y`
for the :math:`xy` component, and the appropriate elements of the moment of
inertia tensor are used. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
:math:`\omega^2` are replaced by :math:`v_x v_y` and :math:`\omega_x
\omega_y` for the :math:`xy` component, and the appropriate elements of
the moment of inertia tensor are used. The six components of the vector
are ordered :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for the
components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` and :math:`\omega^2` are replaced by :math:`v_x v_y`
and :math:`\omega_x \omega_y` for the :math:`xy` component, and so on.
And the appropriate elements of the moment of inertia tensor are used.
Note that because it lacks the 1/2 factor, these tensor components are
twice those of the traditional kinetic energy tensor. The six
components of the vector are ordered :math:`xx`, :math:`yy`, :math:`zz`,
:math:`xy`, :math:`xz`, :math:`yz`.
The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the *dynamic/dof* option of
@ -131,27 +145,26 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input.
See the :doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
vector of length 6 (symmetric tensor), which can be accessed by indices
1--6. These values can be used by any command that uses global scalar
or vector values from a compute as input. See the :doc:`Howto output
<Howto_output>` page for an overview of LAMMPS output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value calculated by this compute is "intensive". The vector
values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""
This compute is part of the ASPHERE package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This compute requires that atoms store angular momentum and a
quaternion as defined by the :doc:`atom_style ellipsoid <atom_style>`
command.
This compute requires that atoms store angular momentum and a quaternion
as defined by the :doc:`atom_style ellipsoid <atom_style>` command.
All particles in the group must be finite-size. They cannot be point
particles, but they can be aspherical or spherical as defined by their

29
doc/src/compute_temp_body.rst
View File

@ -62,12 +62,17 @@ kinetic energy is computed as :math:`\frac12 I \omega^2`, where :math:`I`
is the moment of inertia tensor for the aspherical particle and :math:`\omega`
is its angular velocity, which is computed from its angular momentum.
A kinetic energy tensor, stored as a 6-element vector, is also calculated by
this compute. The formula for the components of the tensor is the same as the
above formula, except that :math:`v^2` and :math:`\omega^2` are
replaced by :math:`v_x v_y` and :math:`\omega_x \omega_y` for the
math:`xy` component, and the appropriate elements of the inertia tensor are
used. The six components of the vector are ordered :math:`xx`, :math:`yy`,
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` and :math:`\omega^2` are replaced by :math:`v_x v_y`
and :math:`\omega_x \omega_y` for the :math:`xy` component, and so on.
And the appropriate elements of the moment of inertia tensor are used.
Note that because it lacks the 1/2 factor, these tensor components are
twice those of the traditional kinetic energy tensor. The six
components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
The number of atoms contributing to the temperature is assumed to be
@ -111,17 +116,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input.
See the :doc:`Howto output <Howto_output>` page for an overview of LAMMPS
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`.
The vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

28
doc/src/compute_temp_chunk.rst
View File

@ -85,12 +85,14 @@ By default, *adof* = 2 or 3 = dimensionality of system, as set via the
:doc:`dimension <dimension>` command, and *cdof* = 0.0.
This gives the usual formula for temperature.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute for use in the computation of a pressure
tensor. The formula for the components of the tensor is the same as
the above formula, except that :math:`v^2` is replaced by
:math:`v_x v_y` for the :math:`xy` component, and so on.
The six components of the vector are ordered :math:`xx`, :math:`yy`,
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute. The formula for the components of the tensor is the
same as the above expression for :math:`E_\mathrm{kin}`, except that
the 1/2 factor is NOT included and the :math:`v_i^2` is replaced by
:math:`v_{i,x} v_{i,y}` for the :math:`xy` component, and so on. Note
that because it lacks the 1/2 factor, these tensor components are
twice those of the traditional kinetic energy tensor. The six
components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
Note that the number of atoms contributing to the temperature is
@ -227,10 +229,10 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input.
See the :doc:`Howto output <Howto_output>` page for an overview of LAMMPS
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
This compute also optionally calculates a global array, if one or more
@ -245,9 +247,9 @@ page for an overview of LAMMPS output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive". The array values are "intensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`. The array values
will be in temperature :doc:`units <units>` for the *temp* value, and in
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`. The array values will be
in temperature :doc:`units <units>` for the *temp* value, and in
energy :doc:`units <units>` for the *kecom* and *internal* values.
Restrictions

31
doc/src/compute_temp_com.rst
View File

@ -44,12 +44,17 @@ where KE is the total kinetic energy of the group of atoms (sum of
simulation, :math:`N` is number of atoms in the group, :math:`k_B` is
the Boltzmann constant, and :math:`T` is the absolute temperature.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute for use in the computation of a pressure
tensor. The formula for the components of the tensor is the same as
the above formula, except that :math:`v^2` is replaced by :math:`v_x v_y`
for the :math:`xy` component, and so on. The six components of the vector are
ordered :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the *dynamic* option of the
@ -81,17 +86,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`.
The vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values is in energy :doc:`units <units>`.
Restrictions
""""""""""""

108
doc/src/compute_temp_cs.rst
View File

@ -31,27 +31,27 @@ on the center-of-mass velocity of atom pairs that are bonded to each
other. This compute is designed to be used with the adiabatic
core/shell model of :ref:`(Mitchell and Fincham) <MitchellFincham1>`.
See the :doc:`Howto coreshell <Howto_coreshell>` page for an overview of
the model as implemented in LAMMPS. Specifically, this compute
enables correct temperature calculation and thermostatting of
core/shell pairs where it is desirable for the internal degrees of
freedom of the core/shell pairs to not be influenced by a thermostat.
A compute of this style can be used by any command that computes a
temperature via :doc:`fix_modify <fix_modify>`
(e.g., :doc:`fix temp/rescale <fix_temp_rescale>`, :doc:`fix npt <fix_nh>`).
the model as implemented in LAMMPS. Specifically, this compute enables
correct temperature calculation and thermostatting of core/shell pairs
where it is desirable for the internal degrees of freedom of the
core/shell pairs to not be influenced by a thermostat. A compute of
this style can be used by any command that computes a temperature via
:doc:`fix_modify <fix_modify>` (e.g., :doc:`fix temp/rescale
<fix_temp_rescale>`, :doc:`fix npt <fix_nh>`).
Note that this compute does not require all ions to be polarized,
hence defined as core/shell pairs. One can mix core/shell pairs and
ions without a satellite particle if desired. The compute will
consider the non-polarized ions according to the physical system.
Note that this compute does not require all ions to be polarized, hence
defined as core/shell pairs. One can mix core/shell pairs and ions
without a satellite particle if desired. The compute will consider the
non-polarized ions according to the physical system.
For this compute, core and shell particles are specified by two
respective group IDs, which can be defined using the
:doc:`group <group>` command. The number of atoms in the two groups
must be the same and there should be one bond defined between a pair
of atoms in the two groups. Non-polarized ions which might also be
included in the treated system should not be included into either of
these groups, they are taken into account by the *group-ID* (second
argument) of the compute.
respective group IDs, which can be defined using the :doc:`group
<group>` command. The number of atoms in the two groups must be the
same and there should be one bond defined between a pair of atoms in the
two groups. Non-polarized ions which might also be included in the
treated system should not be included into either of these groups, they
are taken into account by the *group-ID* (second argument) of the
compute.
The temperature is calculated by the formula
@ -60,52 +60,56 @@ The temperature is calculated by the formula
\text{KE} = \frac{\text{dim}}{2} N k_B T,
where KE is the total kinetic energy of the group of atoms (sum of
:math:`\frac12 m v^2`), dim = 2 or 3 is the dimensionality of the simulation,
:math:`N` is the number of atoms in the group, :math:`k_B` is the Boltzmann
constant, and :math:`T` is the absolute temperature. Note that
the velocity of each core or shell atom used in the KE calculation is
the velocity of the center-of-mass (COM) of the core/shell pair the
atom is part of.
:math:`\frac12 m v^2`), dim = 2 or 3 is the dimensionality of the
simulation, :math:`N` is the number of atoms in the group, :math:`k_B`
is the Boltzmann constant, and :math:`T` is the absolute temperature.
Note that the velocity of each core or shell atom used in the KE
calculation is the velocity of the center-of-mass (COM) of the
core/shell pair the atom is part of.
A kinetic energy tensor, stored as a six-element vector, is also calculated by
this compute for use in the computation of a pressure tensor. The formula for
the components of the tensor is the same as the above formula, except that
:math:`v^2` is replaced by :math:`v_x v_y` for the :math:`xy` component, and so
on. The six components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`. In contrast to the temperature,
the velocity of each core or shell atom is taken individually.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for the
components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
The change this fix makes to core/shell atom velocities is essentially
computing the temperature after a "bias" has been removed from the velocity of
the atoms. This "bias" is the velocity of the atom relative to the
center-of-mass velocity of the core/shell pair. If this compute is used with a
fix command that performs thermostatting then this bias will be subtracted from
each atom, thermostatting of the remaining center-of-mass velocity will be
performed, and the bias will be added back in. This means the thermostatting
will effectively be performed on the core/shell pairs, instead of on the
individual core and shell atoms. Thermostatting fixes that work in this way
include :doc:`fix nvt <fix_nh>`, :doc:`fix temp/rescale <fix_temp_rescale>`,
:doc:`fix temp/berendsen <fix_temp_berendsen>`, and
:doc:`fix langevin <fix_langevin>`.
computing the temperature after a "bias" has been removed from the
velocity of the atoms. This "bias" is the velocity of the atom relative
to the center-of-mass velocity of the core/shell pair. If this compute
is used with a fix command that performs thermostatting then this bias
will be subtracted from each atom, thermostatting of the remaining
center-of-mass velocity will be performed, and the bias will be added
back in. This means the thermostatting will effectively be performed on
the core/shell pairs, instead of on the individual core and shell atoms.
Thermostatting fixes that work in this way include :doc:`fix nvt
<fix_nh>`, :doc:`fix temp/rescale <fix_temp_rescale>`, :doc:`fix
temp/berendsen <fix_temp_berendsen>`, and :doc:`fix langevin
<fix_langevin>`.
The internal energy of core/shell pairs can be calculated by the
:doc:`compute temp/chunk <compute_temp_chunk>` command, if chunks are defined
as core/shell pairs. See the :doc:`Howto coreshell <Howto_coreshell>` doc
page for more discussion on how to do this.
:doc:`compute temp/chunk <compute_temp_chunk>` command, if chunks are
defined as core/shell pairs. See the :doc:`Howto coreshell
<Howto_coreshell>` doc page for more discussion on how to do this.
Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input.
vector of length 6 (symmetric tensor), which can be accessed by indices
1--6. These values can be used by any command that uses global scalar
or vector values from a compute as input.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value calculated by this compute is "intensive". The vector
values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

30
doc/src/compute_temp_deform.rst
View File

@ -73,12 +73,16 @@ simulation, :math:`N` is the number of atoms in the group, :math:`k_B`
is the Boltzmann constant, and :math:`T` is the temperature. Note that
:math:`v` in the kinetic energy formula is the atom's velocity.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute for use in the computation of a pressure
tensor. The formula for the components of the tensor is the same as
the above formula, except that :math:`v^2` is replaced by :math:`v_x v_y` for
the :math:`xy` component, and so on. The six components of the vector are
ordered :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is assumed to be
@ -128,17 +132,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`.
The vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

33
doc/src/compute_temp_deform_eff.rst
View File

@ -29,17 +29,20 @@ model, after subtracting out a streaming velocity induced by the
simulation box changing size and/or shape, for example in a
non-equilibrium MD (NEMD) simulation. The size/shape change is
induced by use of the :doc:`fix deform <fix_deform>` command. A
compute of this style is created by the
:doc:`fix nvt/sllod/eff <fix_nvt_sllod_eff>` command to compute the thermal
temperature of atoms for thermostatting purposes. A compute of this
style can also be used by any command that computes a temperature
(e.g., :doc:`thermo_modify <thermo_modify>`, :doc:`fix npt/eff <fix_nh_eff>`).
compute of this style is created by the :doc:`fix nvt/sllod/eff
<fix_nvt_sllod_eff>` command to compute the thermal temperature of
atoms for thermostatting purposes. A compute of this style can also
be used by any command that computes a temperature (e.g.,
:doc:`thermo_modify <thermo_modify>`, :doc:`fix npt/eff
<fix_nh_eff>`).
The calculation performed by this compute is exactly like that
described by the :doc:`compute temp/deform <compute_temp_deform>`
command, except that the formula for the temperature includes the
radial electron velocity contributions, as discussed by the :doc:`compute temp/eff <compute_temp_eff>` command. Note that only the
translational degrees of freedom for each nuclei or electron are
command, except that the formulas for the temperature (scalar) and
diagonal components of the symmetric tensor (vector) include the
radial electron velocity contributions, as discussed by the
:doc:`compute temp/eff <compute_temp_eff>` command. Note that only
the translational degrees of freedom for each nuclei or electron are
affected by the streaming velocity adjustment. The radial velocity
component of the electrons is not affected.
@ -47,17 +50,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

28
doc/src/compute_temp_partial.rst
View File

@ -44,12 +44,16 @@ constant, and :math:`T` = temperature. The calculation of KE excludes the
is 0. The dim parameter is adjusted to give the correct number of
degrees of freedom.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute for use in the calculation of a pressure
tensor. The formula for the components of the tensor is the same as
the above formula, except that :math:`v^2` is replaced by :math:`v_x v_y` for
the :math:`xy` component, and so on. The six components of the vector are
ordered :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is assumed to be
@ -88,17 +92,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input.
See the :doc:`Howto output <Howto_output>` page for an overview of LAMMPS
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

64
doc/src/compute_temp_profile.rst
View File

@ -97,21 +97,27 @@ center-of-mass velocity across the group in directions where streaming velocity
is *not* subtracted. This can be altered using the *extra* option of the
:doc:`compute_modify <compute_modify>` command.
If the *out* keyword is used with a *tensor* value, which is the default,
a kinetic energy tensor, stored as a six-element vector, is also calculated by
this compute for use in the computation of a pressure tensor. The formula for
the components of the tensor is the same as the above formula, except that
:math:`v^2` is replaced by :math:`v_x v_y` for the :math:`xy` component, and
so on. The six components of the vector are ordered :math:`xx`, :math:`yy`,
If the *out* keyword is used with a *tensor* value, which is the
default, then a symmetric tensor, stored as a six-element vector, is
also calculated by this compute for use in the computation of a
pressure tensor by the :doc:`compute pressue <compute_pressure>`
command. The formula for the components of the tensor is the same as
the above expression for :math:`E_\mathrm{kin}`, except that the 1/2
factor is NOT included and the :math:`v_i^2` is replaced by
:math:`v_{i,x} v_{i,y}` for the :math:`xy` component, and so on. Note
that because it lacks the 1/2 factor, these tensor components are
twice those of the traditional kinetic energy tensor. The six
components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
If the *out* keyword is used with a *bin* value, the count of atoms and
computed temperature for each bin are stored for output, as an array of values,
as described below. The temperature of each bin is calculated as described
above, where the bias velocity is subtracted and only the remaining thermal
velocity of atoms in the bin contributes to the temperature. See the note
below for how the temperature is normalized by the degrees-of-freedom of atoms
in the bin.
If the *out* keyword is used with a *bin* value, the count of atoms
and computed temperature for each bin are stored for output, as an
array of values, as described below. The temperature of each bin is
calculated as described above, where the bias velocity is subtracted
and only the remaining thermal velocity of atoms in the bin
contributes to the temperature. See the note below for how the
temperature is normalized by the degrees-of-freedom of atoms in the
bin.
The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the *dynamic* option of the
@ -166,16 +172,17 @@ Output info
This compute calculates a global scalar (the temperature). Depending
on the setting of the *out* keyword, it also calculates a global
vector or array. For *out* = *tensor*, it calculates a vector of
length 6 (KE tensor), which can be accessed by indices 1--6. For *out*
= *bin* it calculates a global array which has 2 columns and :math:`N` rows,
where :math:`N` is the number of bins. The first column contains the number
of atoms in that bin. The second contains the temperature of that
bin, calculated as described above. The ordering of rows in the array
is as follows. Bins in :math:`x` vary fastest, then :math:`y`, then
:math:`z`. Thus for a :math:`10\times 10\times 10` 3d array of bins, there
will be 1000 rows. The bin with indices :math:`(i_x,i_y,i_z) = (2,3,4)` would
map to row :math:`M = 10^2(i_z-1) + 10(i_y-1) + i_x = 322`, where the rows are
numbered from 1 to 1000 and the bin indices are numbered from 1 to 10 in each
length 6 (symmetric tensor), which can be accessed by indices 1--6.
For *out* = *bin* it calculates a global array which has 2 columns and
:math:`N` rows, where :math:`N` is the number of bins. The first
column contains the number of atoms in that bin. The second contains
the temperature of that bin, calculated as described above. The
ordering of rows in the array is as follows. Bins in :math:`x` vary
fastest, then :math:`y`, then :math:`z`. Thus for a :math:`10\times
10\times 10` 3d array of bins, there will be 1000 rows. The bin with
indices :math:`(i_x,i_y,i_z) = (2,3,4)` would map to row :math:`M =
10^2(i_z-1) + 10(i_y-1) + i_x = 322`, where the rows are numbered from
1 to 1000 and the bin indices are numbered from 1 to 10 in each
dimension.
These values can be used by any command that uses global scalar or
@ -186,9 +193,9 @@ options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive". The array values are "intensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`. The first column
of array values are counts; the values in the second column will be in
The scalar value us in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`. The first column of array
values are counts; the values in the second column will be in
temperature :doc:`units <units>`.
Restrictions
@ -203,7 +210,10 @@ will be for most thermostats.
Related commands
""""""""""""""""
:doc:`compute temp <compute_temp>`, :doc:`compute temp/ramp <compute_temp_ramp>`, :doc:`compute temp/deform <compute_temp_deform>`, :doc:`compute pressure <compute_pressure>`
:doc:`compute temp <compute_temp>`,
:doc:`compute temp/ramp <compute_temp_ramp>`,
:doc:`compute temp/deform <compute_temp_deform>`,
:doc:`compute pressure <compute_pressure>`
Default
"""""""

36
doc/src/compute_temp_ramp.rst
View File

@ -63,12 +63,17 @@ command (e.g., :math:`\AA` for units = real or metal). A
velocity in lattice spacings per unit time). The :doc:`lattice <lattice>`
command must have been previously used to define the lattice spacing.
A kinetic energy tensor, stored as a six-element vector, is also calculated by
this compute for use in the computation of a pressure tensor. The formula for
the components of the tensor is the same as the above formula, except that
:math:`v^2` is replaced by :math:`v_x v_y` for the :math:`xy` component, and
so on. The six components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is assumed to be constant
for the duration of the run; use the *dynamic* option of the
@ -100,17 +105,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""
@ -119,7 +124,10 @@ Restrictions
Related commands
""""""""""""""""
:doc:`compute temp <compute_temp>`, :doc:`compute temp/profie <compute_temp_profile>`, :doc:`compute temp/deform <compute_temp_deform>`, :doc:`compute pressure <compute_pressure>`
:doc:`compute temp <compute_temp>`,
:doc:`compute temp/profile <compute_temp_profile>`,
:doc:`compute temp/deform <compute_temp_deform>`,
:doc:`compute pressure <compute_pressure>`
Default
"""""""

44
doc/src/compute_temp_region.rst
View File

@ -49,12 +49,17 @@ where KE = is the total kinetic energy of the group of atoms (sum of
:math:`N` is the number of atoms in both the group and region, :math:`k_B` is
the Boltzmann constant, and :math:`T` temperature.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute for use in the computation of a pressure
tensor. The formula for the components of the tensor is the same as
the above formula, except that :math:`v^2` is replaced by :math:`v_x v_y`
for the :math:`xy` component, and so on. The six components of the vector are
ordered :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is calculated each
time the temperature is evaluated since it is assumed atoms can
@ -78,12 +83,13 @@ will operate only on atoms that are currently in the geometric region.
Unlike other compute styles that calculate temperature, this compute
does not subtract out degrees-of-freedom due to fixes that constrain
motion, such as :doc:`fix shake <fix_shake>` and :doc:`fix rigid <fix_rigid>`. This is because those degrees of freedom
(e.g., a constrained bond) could apply to sets of atoms that straddle
the region boundary, and hence the concept is somewhat ill-defined.
If needed the number of subtracted degrees of freedom can be set
explicitly using the *extra* option of the
:doc:`compute_modify <compute_modify>` command.
motion, such as :doc:`fix shake <fix_shake>` and :doc:`fix rigid
<fix_rigid>`. This is because those degrees of freedom (e.g., a
constrained bond) could apply to sets of atoms that straddle the
region boundary, and hence the concept is somewhat ill-defined. If
needed the number of subtracted degrees of freedom can be set
explicitly using the *extra* option of the :doc:`compute_modify
<compute_modify>` command.
See the :doc:`Howto thermostat <Howto_thermostat>` page for a
discussion of different ways to compute temperature and perform
@ -93,17 +99,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`.
The vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

27
doc/src/compute_temp_region_eff.rst
View File

@ -32,32 +32,33 @@ temperature (e.g., :doc:`thermo_modify <thermo_modify>`).
The operation of this compute is exactly like that described by the
:doc:`compute temp/region <compute_temp_region>` command, except that
the formula for the temperature itself includes the radial electron
velocity contributions, as discussed by the
:doc:`compute temp/eff <compute_temp_eff>` command.
the formulas for the temperature (scalar) and diagonal components of
the symmetric tensor (vector) include the radial electron velocity
contributions, as discussed by the :doc:`compute temp/eff
<compute_temp_eff>` command.
Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""
This compute is part of the EFF package. It is only enabled if
LAMMPS was built with that package.
See the :doc:`Build package <Build_package>` page for more info.
This compute is part of the EFF package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Related commands
""""""""""""""""

30
doc/src/compute_temp_rotate.rst
View File

@ -43,12 +43,17 @@ where KE is the total kinetic energy of the group of atoms (sum of
:math:`N` is the number of atoms in the group, :math:`k_B` is the Boltzmann
constant, and :math:`T` is the absolute temperature.
A kinetic energy tensor, stored as a six-element vector, is also calculated by
this compute for use in the computation of a pressure tensor. The formula for
the components of the tensor is the same as the above formula, except that
:math:`v^2` is replaced by :math:`v_x v_y` for the :math:`xy` component, and
so on. The six components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the *dynamic* option of the
@ -80,17 +85,16 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1-6.
These values can be used by any command that uses global scalar or
vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
vector of length 6 (symmetric tensor), which can be accessed by
indices 1-6. These values can be used by any command that uses global
scalar or vector values from a compute as input. See the :doc:`Howto
output <Howto_output>` page for an overview of LAMMPS output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

24
doc/src/compute_temp_sphere.rst
View File

@ -77,6 +77,18 @@ tensor is the same as the above formulas, except that :math:`v^2` and
vector are ordered :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`,
:math:`xz`, :math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` and :math:`\omega^2` are replaced by :math:`v_x v_y`
and :math:`\omega_x \omega_y` for the :math:`xy` component, and so on.
Note that because it lacks the 1/2 factor, these tensor components are
twice those of the traditional kinetic energy tensor. The six
components of the vector are ordered :math:`xx`, :math:`yy`,
:math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the *dynamic* option of the
:doc:`compute_modify <compute_modify>` command if this is not the case.
@ -117,17 +129,17 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 6 (KE tensor), which can be accessed by indices 1--6.
These values can be used by any command that uses global scalar or
vector values from a compute as input.
See the :doc:`Howto output <Howto_output>` page for an overview of LAMMPS
vector of length 6 (symmetric tensor), which can be accessed by
indices 1--6. These values can be used by any command that uses
global scalar or vector values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
The scalar value calculated by this compute is "intensive". The
vector values are "extensive".
The scalar value will be in temperature :doc:`units <units>`. The
vector values will be in energy :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The vector
values are in energy :doc:`units <units>`.
Restrictions
""""""""""""

37
doc/src/compute_viscosity_cos.rst
View File

@ -86,12 +86,17 @@ where KE is the total kinetic energy of the group of atoms (sum of
:math:`N` is the number of atoms in the group, :math:`k_B` is the Boltzmann
constant, and :math:`T` is the absolute temperature.
A kinetic energy tensor, stored as a six-element vector, is also
calculated by this compute for use in the computation of a pressure
tensor. The formula for the components of the tensor is the same as
the above formula, except that :math:`v^2` is replaced by :math:`v_x v_y` for
the :math:`xy` component, and so on. The six components of the vector are
ordered :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the
:math:`xy` component, and so on. Note that because it lacks the 1/2
factor, these tensor components are twice those of the traditional
kinetic energy tensor. The six components of the vector are ordered
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`.
The number of atoms contributing to the temperature is assumed to be
constant for the duration of the run; use the *dynamic* option of the
@ -126,21 +131,21 @@ Output info
"""""""""""
This compute calculates a global scalar (the temperature) and a global
vector of length 7, which can be accessed by indices 1--7.
The first six elements of the vector are the KE tensor,
and the seventh is the cosine-shaped velocity amplitude :math:`V`,
which can be used to calculate the reciprocal viscosity, as shown in the example.
These values can be used by any command that uses global scalar or
vector values from a compute as input.
See the :doc:`Howto output <Howto_output>` page for an overview of LAMMPS output options.
vector of length 7, which can be accessed by indices 1--7. The first
six elements of the vector are those of the symmetric tensor discussed
above. The seventh is the cosine-shaped velocity amplitude :math:`V`,
which can be used to calculate the reciprocal viscosity, as shown in
the example. These values can be used by any command that uses global
scalar or vector values from a compute as input. See the :doc:`Howto
output <Howto_output>` page for an overview of LAMMPS output options.
The scalar value calculated by this compute is "intensive". The
first six elements of vector values are "extensive",
and the seventh element of vector values is "intensive".
The scalar value will be in temperature :doc:`units <units>`.
The first six elements of vector values will be in energy :doc:`units <units>`.
The seventh element of vector value will be in velocity :doc:`units <units>`.
The scalar value is in temperature :doc:`units <units>`. The first
six elements of vector values are in energy :doc:`units <units>`. The
seventh element of vector value us in velocity :doc:`units <units>`.
Restrictions
""""""""""""

4
doc/src/fix_bond_react.rst
View File

@ -785,3 +785,7 @@ reset_mol_ids = yes, custom_charges = no, molecule = off, modify_create = *fit a
.. _Gissinger2020:
**(Gissinger2020)** Gissinger, Jensen and Wise, Macromolecules, 53, 22, 9953-9961 (2020).
.. _Gissinger2024:
**(Gissinger2024)** Gissinger, Jensen and Wise, Computer Physics Communications, 304, 109287 (2024).

1
examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.stabilized
View File

@ -44,6 +44,7 @@ thermo 50
fix myrxns all bond/react stabilization yes statted_grp .03 &
react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map rescale_charges yes
fix 1 statted_grp_REACT nvt temp 300 300 100

2
examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.stabilized_variable_probability
View File

@ -47,7 +47,7 @@ thermo 50
fix myrxns all bond/react stabilization yes statted_grp .03 &
react rxn1 all 1 0.0 5.0 mol1 mol2 rxn1_stp1_map prob v_prob1 1234 &
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map prob v_prob2 1234
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map prob v_prob2 1234 rescale_charges yes
fix 1 statted_grp_REACT nvt temp 300 300 100

2
examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.unstabilized
View File

@ -44,7 +44,7 @@ thermo 50
fix myrxns all bond/react stabilization no &
react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map rescale_charges yes
fix 1 all nve/limit .03

21
examples/PACKAGES/reaction/tiny_nylon/rxn1_stp1_reacted.molecule_template
View File

@ -48,27 +48,6 @@ Types
17 hc
18 hc
Charges
1 -0.300000
2 0.000000
3 0.000000
4 0.000000
5 0.000000
6 0.000000
7 0.000000
8 0.000000
9 0.000000
10 0.300000
11 0.000000
12 0.000000
13 0.000000
14 0.000000
15 0.000000
16 0.000000
17 0.000000
18 0.000000
Molecules
1 1

30
examples/PACKAGES/reaction/tiny_nylon/rxn1_stp2_reacted.molecule_template
View File

@ -44,21 +44,21 @@ Types
Charges
1 -0.300000
2 0.000000
3 0.000000
4 0.410000
5 0.000000
6 0.000000
7 0.000000
8 0.000000
9 0.000000
10 0.300000
11 0.000000
12 -0.820000
13 0.000000
14 0.000000
15 0.410000
1 -0.60533
2 -0.01149
3 -0.76306
4 0.38
5 0.29346
6 0.18360
7 0.15396
8 -0.72636
9 -0.27437
10 0.40603
11 -0.65530
12 -0.76
13 0.21423
14 0.18949
15 0.38
Molecules

1263
examples/PACKAGES/reaction/tiny_nylon/tiny_nylon.data
View File

File diff suppressed because it is too large Load Diff

239
src/REACTION/fix_bond_react.cpp
View File

@ -58,7 +58,7 @@ using namespace MathConst;
static const char cite_fix_bond_react[] =
"fix bond/react: reacter.org doi:10.1016/j.polymer.2017.09.038, "
"doi:10.1021/acs.macromol.0c02012\n\n"
"doi:10.1021/acs.macromol.0c02012, doi:10.1016/j.cpc.2024.109287\n\n"
"@Article{Gissinger17,\n"
" author = {J. R. Gissinger and B. D. Jensen and K. E. Wise},\n"
" title = {Modeling Chemical Reactions in Classical Molecular Dynamics Simulations},\n"
@ -75,6 +75,14 @@ static const char cite_fix_bond_react[] =
" volume = 53,\n"
" number = 22,\n"
" pages = {9953--9961}\n"
"}\n\n"
"@Article{Gissinger24,\n"
" author = {J. R. Gissinger, B. D. Jensen, K. E. Wise},\n"
" title = {Molecular Modeling of Reactive Systems with REACTER},\n"
" journal = {Computer Physics Communications},\n"
" year = 2024,\n"
" volume = 304,\n"
" number = 109287\n"
"}\n\n";
static constexpr double BIG = 1.0e20;
@ -225,8 +233,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(reacted_mol,nreacts,"bond/react:reacted_mol");
memory->create(fraction,nreacts,"bond/react:fraction");
memory->create(max_rxn,nreacts,"bond/react:max_rxn");
memory->create(nlocalskips,nreacts,"bond/react:nlocalskips");
memory->create(nghostlyskips,nreacts,"bond/react:nghostlyskips");
memory->create(nlocalkeep,nreacts,"bond/react:nlocalkeep");
memory->create(nghostlykeep,nreacts,"bond/react:nghostlykeep");
memory->create(seed,nreacts,"bond/react:seed");
memory->create(limit_duration,nreacts,"bond/react:limit_duration");
memory->create(rate_limit,3,nreacts,"bond/react:rate_limit");
@ -486,10 +494,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
get_molxspecials();
read_map_file(i);
fclose(fp);
if (ncreate == 0 && onemol->natoms != twomol->natoms)
error->all(FLERR,"Fix bond/react: Reaction templates must contain the same number of atoms");
else if (ncreate > 0 && onemol->natoms + ncreate != twomol->natoms)
error->all(FLERR,"Fix bond/react: Incorrect number of created atoms");
iatomtype[i] = onemol->type[ibonding[i]-1];
jatomtype[i] = onemol->type[jbonding[i]-1];
find_landlocked_atoms(i);
@ -644,8 +648,8 @@ FixBondReact::~FixBondReact()
memory->destroy(fraction);
memory->destroy(seed);
memory->destroy(max_rxn);
memory->destroy(nlocalskips);
memory->destroy(nghostlyskips);
memory->destroy(nlocalkeep);
memory->destroy(nghostlykeep);
memory->destroy(limit_duration);
memory->destroy(var_flag);
memory->destroy(var_id);
@ -716,6 +720,7 @@ int FixBondReact::setmask()
int mask = 0;
mask |= POST_INTEGRATE;
mask |= POST_INTEGRATE_RESPA;
mask |= POST_FORCE;
return mask;
}
@ -872,8 +877,8 @@ void FixBondReact::post_integrate()
reaction_count[i] = 0;
local_rxn_count[i] = 0;
ghostly_rxn_count[i] = 0;
nlocalskips[i] = 0;
nghostlyskips[i] = 0;
nlocalkeep[i] = INT_MAX;
nghostlykeep[i] = INT_MAX;
// update reaction probability
if (var_flag[PROB][i])
fraction[i] = input->variable->compute_equal(var_id[PROB][i]);
@ -1424,10 +1429,13 @@ void FixBondReact::superimpose_algorithm()
MPI_Allreduce(&local_rxn_count[0],&reaction_count[0],nreacts,MPI_INT,MPI_SUM,world);
int rxnflag = 0;
int *delta_rxn;
memory->create(delta_rxn,nreacts,"bond/react:delta_rxn");
if (comm->me == 0)
for (int i = 0; i < nreacts; i++) {
reaction_count_total[i] += reaction_count[i] + ghostly_rxn_count[i];
rxnflag += reaction_count[i] + ghostly_rxn_count[i];
delta_rxn[i] = reaction_count[i] + ghostly_rxn_count[i];
reaction_count_total[i] += delta_rxn[i];
rxnflag += delta_rxn[i];
}
MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world);
@ -1460,42 +1468,43 @@ void FixBondReact::superimpose_algorithm()
if (overstep > 0) {
// let's randomly choose rxns to skip, unbiasedly from local and ghostly
int *local_rxncounts;
int *all_localskips;
int *all_localkeep;
memory->create(local_rxncounts,nprocs,"bond/react:local_rxncounts");
memory->create(all_localskips,nprocs,"bond/react:all_localskips");
memory->create(all_localkeep,nprocs,"bond/react:all_localkeep");
MPI_Gather(&local_rxn_count[i],1,MPI_INT,local_rxncounts,1,MPI_INT,0,world);
if (comm->me == 0) {
int delta_rxn = reaction_count[i] + ghostly_rxn_count[i];
// when using variable input for rate_limit, rate_limit_overstep could be > delta_rxn (below)
// we need to limit overstep to the number of reactions on this timestep
// essentially skipping all reactions, would be more efficient to use a skip_all flag
if (overstep > delta_rxn) overstep = delta_rxn;
if (overstep > delta_rxn[i]) overstep = delta_rxn[i];
int nkeep = delta_rxn[i] - overstep;
int *rxn_by_proc;
memory->create(rxn_by_proc,delta_rxn,"bond/react:rxn_by_proc");
for (int j = 0; j < delta_rxn; j++)
memory->create(rxn_by_proc,delta_rxn[i],"bond/react:rxn_by_proc");
for (int j = 0; j < delta_rxn[i]; j++)
rxn_by_proc[j] = -1; // corresponds to ghostly
int itemp = 0;
for (int j = 0; j < nprocs; j++)
for (int k = 0; k < local_rxncounts[j]; k++)
rxn_by_proc[itemp++] = j;
std::shuffle(&rxn_by_proc[0],&rxn_by_proc[delta_rxn], park_rng);
std::shuffle(&rxn_by_proc[0],&rxn_by_proc[delta_rxn[i]], park_rng);
for (int j = 0; j < nprocs; j++)
all_localskips[j] = 0;
nghostlyskips[i] = 0;
for (int j = 0; j < overstep; j++) {
if (rxn_by_proc[j] == -1) nghostlyskips[i]++;
else all_localskips[rxn_by_proc[j]]++;
all_localkeep[j] = 0;
nghostlykeep[i] = 0;
for (int j = 0; j < nkeep; j++) {
if (rxn_by_proc[j] == -1) nghostlykeep[i]++;
else all_localkeep[rxn_by_proc[j]]++;
}
memory->destroy(rxn_by_proc);
reaction_count_total[i] -= overstep;
}
MPI_Scatter(&all_localskips[0],1,MPI_INT,&nlocalskips[i],1,MPI_INT,0,world);
MPI_Bcast(&nghostlyskips[i],1,MPI_INT,0,world);
MPI_Scatter(&all_localkeep[0],1,MPI_INT,&nlocalkeep[i],1,MPI_INT,0,world);
MPI_Bcast(&nghostlykeep[i],1,MPI_INT,0,world);
memory->destroy(local_rxncounts);
memory->destroy(all_localskips);
memory->destroy(all_localkeep);
}
}
MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world);
memory->destroy(delta_rxn);
// this updates topology next step
next_reneighbor = update->ntimestep;
@ -2965,6 +2974,8 @@ void FixBondReact::update_everything()
int *type = atom->type;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
tagint *tag = atom->tag;
AtomVec *avec = atom->avec;
int **bond_type = atom->bond_type;
tagint **bond_atom = atom->bond_atom;
@ -2977,13 +2988,16 @@ void FixBondReact::update_everything()
memory->create(mark,nmark,"bond/react:mark");
for (int i = 0; i < nmark; i++) mark[i] = 0;
// used when creating atoms
addatomtag = 0;
for (int i = 0; i < nlocal; i++) addatomtag = MAX(addatomtag,tag[i]);
MPI_Allreduce(MPI_IN_PLACE,&addatomtag,1,MPI_LMP_TAGINT,MPI_MAX,world);
addatoms.clear();
// flag used to delete special interactions
int *delflag;
memory->create(delflag,atom->maxspecial,"bond/react:delflag");
tagint *tag = atom->tag;
AtomVec *avec = atom->avec;
// used when creating atoms
int inserted_atoms_flag = 0;
@ -3026,13 +3040,14 @@ void FixBondReact::update_everything()
for (int pass = 0; pass < 2; pass++) {
update_num_mega = 0;
int *iskip = new int[nreacts];
for (int i = 0; i < nreacts; i++) iskip[i] = 0;
int *noccur = new int[nreacts];
for (int i = 0; i < nreacts; i++) noccur[i] = 0;
if (pass == 0) {
for (int i = 0; i < local_num_mega; i++) {
rxnID = (int) local_mega_glove[0][i];
// reactions already shuffled from dedup procedure, so can skip first N
if (iskip[rxnID]++ < nlocalskips[rxnID]) continue;
// wait, this check needs to be after add atoms, because they can also be 'skipped' due to overlap
if (noccur[rxnID] >= nlocalkeep[rxnID]) continue;
// this will be overwritten if reaction skipped by create_atoms below
update_mega_glove[0][update_num_mega] = (tagint) local_mega_glove[0][i];
@ -3043,13 +3058,14 @@ void FixBondReact::update_everything()
if (create_atoms_flag[rxnID] == 1) {
onemol = atom->molecules[unreacted_mol[rxnID]];
twomol = atom->molecules[reacted_mol[rxnID]];
if (insert_atoms(update_mega_glove,update_num_mega)) {
if (insert_atoms_setup(update_mega_glove,update_num_mega)) {
inserted_atoms_flag = 1;
} else { // create aborted
reaction_count_total[rxnID]--;
continue;
}
}
noccur[rxnID]++;
if (rescale_charges_flag[rxnID]) sim_total_charges[update_num_mega] = local_mega_glove[1][i];
update_num_mega++;
@ -3058,7 +3074,7 @@ void FixBondReact::update_everything()
for (int i = 0; i < global_megasize; i++) {
rxnID = (int) global_mega_glove[0][i];
// reactions already shuffled from dedup procedure, so can skip first N
if (iskip[rxnID]++ < nghostlyskips[rxnID]) continue;
if (noccur[rxnID] >= nghostlykeep[rxnID]) continue;
// this will be overwritten if reaction skipped by create_atoms below
update_mega_glove[0][update_num_mega] = (tagint) global_mega_glove[0][i];
@ -3071,29 +3087,48 @@ void FixBondReact::update_everything()
if (create_atoms_flag[rxnID] == 1) {
onemol = atom->molecules[unreacted_mol[rxnID]];
twomol = atom->molecules[reacted_mol[rxnID]];
if (insert_atoms(update_mega_glove,update_num_mega)) {
if (insert_atoms_setup(update_mega_glove,update_num_mega)) {
inserted_atoms_flag = 1;
} else { // create aborted
reaction_count_total[rxnID]--;
continue;
}
}
noccur[rxnID]++;
if (rescale_charges_flag[rxnID]) sim_total_charges[update_num_mega] = global_mega_glove[1][i];
update_num_mega++;
}
}
delete [] iskip;
delete [] noccur;
if (update_num_mega == 0) continue;
// if inserted atoms and global map exists, reset map now instead
// of waiting for comm since other pre-exchange fixes may use it
// invoke map_init() b/c atom count has grown
// do this once after all atom insertions
if (inserted_atoms_flag == 1 && atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
// insert all atoms for all rxns here
if (inserted_atoms_flag == 1) {
// clear to-be-overwritten ghost info
atom->nghost = 0;
atom->avec->clear_bonus();
for (auto & myaddatom : addatoms) {
atom->avec->create_atom(myaddatom.type,myaddatom.x);
int n = atom->nlocal - 1;
atom->tag[n] = myaddatom.tag;
atom->molecule[n] = myaddatom.molecule;
atom->mask[n] = myaddatom.mask;
atom->image[n] = myaddatom.image;
atom->v[n][0] = myaddatom.v[0];
atom->v[n][1] = myaddatom.v[1];
atom->v[n][2] = myaddatom.v[2];
if (atom->rmass) atom->rmass[n]= myaddatom.rmass;
modify->create_attribute(n);
}
// reset atom->map
if (atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
}
}
// mark to-delete atoms
@ -3620,10 +3655,6 @@ void FixBondReact::update_everything()
atom->natoms -= ndel;
// done deleting atoms
// reset mol ids
if (reset_mol_ids_flag) reset_mol_ids->reset();
// something to think about: this could done much more concisely if
// all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG
@ -3651,10 +3682,11 @@ void FixBondReact::update_everything()
}
/* ----------------------------------------------------------------------
insert created atoms
setup for inserting created atoms
atoms for all rxns are actually created all at once in update_everything
------------------------------------------------------------------------- */
int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
int FixBondReact::insert_atoms_setup(tagint **my_update_mega_glove, int iupdate)
{
// inserting atoms based off fix_deposit->pre_exchange
int flag;
@ -3682,14 +3714,9 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
tagint maxtag_all,maxmol_all;
tagint max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
MPI_Allreduce(&max,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]);
MPI_Allreduce(&max,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
tagint maxmol_all;
for (int i = 0; i < nlocal; i++) maxmol_all = MAX(maxmol_all,molecule[i]);
MPI_Allreduce(MPI_IN_PLACE,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
int dimension = domain->dimension;
@ -3786,6 +3813,26 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
if (abortflag) break;
}
}
// also check against previous to-be-added atoms
if (!abortflag) {
for (auto & myaddatom : addatoms) {
for (int m = 0; m < twomol->natoms; m++) {
if (create_atoms[m][rxnID] == 1) {
delx = coords[m][0] - myaddatom.x[0];
dely = coords[m][1] - myaddatom.x[1];
delz = coords[m][2] - myaddatom.x[2];
domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < overlapsq[rxnID]) {
abortflag = 1;
break;
}
}
}
if (abortflag) break;
}
}
MPI_Allreduce(MPI_IN_PLACE,&abortflag,1,MPI_INT,MPI_MAX,world);
if (abortflag) {
memory->destroy(coords);
@ -3794,12 +3841,6 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
}
}
// clear ghost count and any ghost bonus data internal to AtomVec
// same logic as beginning of Comm::exchange()
// do it now b/c inserting atoms will overwrite ghost atoms
atom->nghost = 0;
atom->avec->clear_bonus();
// check if new atoms are in my sub-box or above it if I am highest proc
// if so, add atom to my list via create_atom()
// initialize additional info about the atoms
@ -3842,40 +3883,46 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
}
int root = 0;
addatomtag++;
if (flag) {
struct AddAtom myaddatom;
root = comm->me;
atom->avec->create_atom(twomol->type[m],coords[m]);
int n = atom->nlocal - 1;
atom->tag[n] = maxtag_all + add_count;
myaddatom.type = twomol->type[m];
myaddatom.x[0] = coords[m][0];
myaddatom.x[1] = coords[m][1];
myaddatom.x[2] = coords[m][2];
myaddatom.tag = addatomtag;
// locally update mega_glove
my_update_mega_glove[preID][iupdate] = atom->tag[n];
my_update_mega_glove[preID][iupdate] = myaddatom.tag;
// !! could do better job choosing mol ID for added atoms
if (atom->molecule_flag) {
if (twomol->moleculeflag) {
atom->molecule[n] = maxmol_all + twomol->molecule[m];
myaddatom.molecule = maxmol_all + twomol->molecule[m];
} else {
atom->molecule[n] = maxmol_all + 1;
myaddatom.molecule = maxmol_all + 1;
}
}
atom->mask[n] = 1 | groupbit;
atom->image[n] = imageflags[m];
myaddatom.mask = 1 | groupbit;
myaddatom.image = imageflags[m];
// guess a somewhat reasonable initial velocity based on reaction site
// further control is possible using bond_react_MASTER_group
// compute |velocity| corresponding to a given temperature t, using specific atom's mass
double mymass = atom->rmass ? atom->rmass[n] : atom->mass[twomol->type[m]];
double vtnorm = sqrt(t / (force->mvv2e / (dimension * force->boltz)) / mymass);
v[n][0] = random[rxnID]->uniform();
v[n][1] = random[rxnID]->uniform();
v[n][2] = random[rxnID]->uniform();
double vnorm = sqrt(v[n][0]*v[n][0] + v[n][1]*v[n][1] + v[n][2]*v[n][2]);
v[n][0] = v[n][0]/vnorm*vtnorm;
v[n][1] = v[n][1]/vnorm*vtnorm;
v[n][2] = v[n][2]/vnorm*vtnorm;
modify->create_attribute(n);
myaddatom.rmass = atom->rmass ? twomol->rmass[m] : atom->mass[twomol->type[m]];
double vtnorm = sqrt(t / (force->mvv2e / (dimension * force->boltz)) / myaddatom.rmass);
double myv[3];
myv[0] = random[rxnID]->uniform();
myv[1] = random[rxnID]->uniform();
myv[2] = random[rxnID]->uniform();
double vnorm = sqrt(myv[0]*myv[0] + myv[1]*myv[1] + myv[2]*myv[2]);
myaddatom.v[0] = myv[0]/vnorm*vtnorm;
myaddatom.v[1] = myv[1]/vnorm*vtnorm;
myaddatom.v[2] = myv[2]/vnorm*vtnorm;
addatoms.push_back(myaddatom);
}
// globally update mega_glove and equivalences
MPI_Allreduce(MPI_IN_PLACE,&root,1,MPI_INT,MPI_SUM,world);
@ -3888,12 +3935,11 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
}
// reset global natoms here
// reset atom map elsewhere, after all calls to 'insert_atoms'
// reset atom map elsewhere, after all calls to 'insert_atoms_setup'
atom->natoms += add_count;
if (atom->natoms < 0)
error->all(FLERR,"Too many total atoms");
maxtag_all += add_count;
if (maxtag_all >= MAXTAGINT)
if (addatomtag >= MAXTAGINT)
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
// atom creation successful
memory->destroy(coords);
@ -3970,6 +4016,11 @@ void FixBondReact::read_map_file(int myrxn)
} else break;
}
if (ncreate == 0 && onemol->natoms != twomol->natoms)
error->all(FLERR,"Fix bond/react: Reaction templates must contain the same number of atoms");
else if (ncreate > 0 && onemol->natoms + ncreate != twomol->natoms)
error->all(FLERR,"Fix bond/react: Incorrect number of created atoms");
// grab keyword and skip next line
parse_keyword(0,line,keyword);
@ -4012,6 +4063,13 @@ void FixBondReact::read_map_file(int myrxn)
}
// error check
for (int i = 0; i < onemol->natoms; i++) {
int my_equiv = reverse_equiv[i][1][myrxn];
if (create_atoms[my_equiv-1][myrxn] == 1)
error->all(FLERR,"Fix bond/react: Created atoms cannot also be listed in Equivalences section\n");
}
// error check
if (bondflag == 0 || equivflag == 0)
error->all(FLERR,"Fix bond/react: Map file missing InitiatorIDs or Equivalences section\n");
@ -4071,6 +4129,8 @@ void FixBondReact::CreateAtoms(char *line, int myrxn)
readline(line);
rv = sscanf(line,"%d",&tmp);
if (rv != 1) error->one(FLERR, "CreateIDs section is incorrectly formatted");
if (tmp > twomol->natoms)
error->one(FLERR,"Fix bond/react: Invalid atom ID in CreateIDs section of map file");
create_atoms[tmp-1][myrxn] = 1;
}
if (twomol->xflag == 0)
@ -4331,6 +4391,13 @@ void FixBondReact::post_integrate_respa(int ilevel, int /*iloop*/)
/* ---------------------------------------------------------------------- */
void FixBondReact::post_force(int /*vflag*/)
{
if (reset_mol_ids_flag) reset_mol_ids->reset();
}
/* ---------------------------------------------------------------------- */
int FixBondReact::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{

14
src/REACTION/fix_bond_react.h
View File

@ -46,6 +46,7 @@ class FixBondReact : public Fix {
void init_list(int, class NeighList *) override;
void post_integrate() override;
void post_integrate_respa(int, int) override;
void post_force(int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
@ -62,7 +63,7 @@ class FixBondReact : public Fix {
int *iatomtype, *jatomtype;
int *seed;
double **cutsq, *fraction;
int *max_rxn, *nlocalskips, *nghostlyskips;
int *max_rxn, *nlocalkeep, *nghostlykeep;
tagint lastcheck;
int stabilization_flag;
int reset_mol_ids_flag;
@ -215,7 +216,7 @@ class FixBondReact : public Fix {
void glove_ghostcheck();
void ghost_glovecast();
void update_everything();
int insert_atoms(tagint **, int);
int insert_atoms_setup(tagint **, int);
void unlimit_bond(); // removes atoms from stabilization, and other post-reaction every-step operations
void dedup_mega_gloves(int); //dedup global mega_glove
void write_restart(FILE *) override;
@ -245,6 +246,15 @@ class FixBondReact : public Fix {
std::map<std::set<tagint>, int> atoms2bond; // maps atom pair to index of local bond array
std::vector<std::vector<Constraint>> constraints;
tagint addatomtag;
struct AddAtom {
tagint tag, molecule;
int type, mask;
imageint image;
double rmass, x[3], v[3];
};
std::vector<AddAtom> addatoms;
// DEBUG
void print_bb();