Merge branch 'lammps:develop' into mliappy_unified

This commit is contained in:
Steven Anaya
2022-09-22 02:14:54 -06:00
committed by GitHub
66 changed files with 6682 additions and 4795 deletions

View File

@ -1,7 +1,7 @@
.TH LAMMPS "1" "3 August 2022" "2022-8-3"
.TH LAMMPS "1" "15 September 2022" "2022-9-15"
.SH NAME
.B LAMMPS
\- Molecular Dynamics Simulator. Version 3 August 2022
\- Molecular Dynamics Simulator. Version 15 September 2022
.SH SYNOPSIS
.B lmp

View File

@ -15,7 +15,9 @@
General commands
================
An alphabetic list of general LAMMPS commands.
An alphabetic list of general LAMMPS commands. Note that style
commands with many variants, can be more easily accessed via the small
table above.
.. table_from_list::
:columns: 5

View File

@ -165,6 +165,7 @@ OPT.
* :doc:`orient/fcc <fix_orient>`
* :doc:`orient/eco <fix_orient_eco>`
* :doc:`pafi <fix_pafi>`
* :doc:`pair <fix_pair>`
* :doc:`phonon <fix_phonon>`
* :doc:`pimd <fix_pimd>`
* :doc:`planeforce <fix_planeforce>`

View File

@ -328,7 +328,7 @@ removed so this update is **required** to avoid compilation failure.
Split of fix STORE into fix STORE/GLOBAL and fix STORE/PERATOM
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. versionchanged:: TBD
.. versionchanged:: 15Sep2022
This change splits the GLOBAL and PERATOM modes of fix STORE into two
separate fixes STORE/GLOBAL and STORE/PERATOM. There was very little
@ -387,7 +387,7 @@ This change is **required** or else the code will not compile.
Use Output::get_dump_by_id() instead of Output::find_dump()
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. versionchanged:: TBD
.. versionchanged:: 15Sep2022
The accessor function to individual dump style instances has been changed
from ``Output::find_dump()`` returning the index of the dump instance in

View File

@ -33,46 +33,6 @@ reference state of a bond. Bonds that are created midway into a run,
such as those created by pouring grains using :doc:`fix pour
<fix_pour>`, are initialized on that timestep.
As bonds can be broken between neighbor list builds, the
:doc:`special_bonds <special_bonds>` command works differently for BPM
bond styles. There are two possible settings which determine how pair
interactions work between bonded particles. First, one can turn off
all pair interactions between bonded particles. Unlike :doc:`bond
quartic <bond_quartic>`, this is not done by subtracting pair forces
during the bond computation but rather by dynamically updating the
special bond list. This is the default behavior of BPM bond styles and
is done by updating the 1-2 special bond list as bonds break. To do
this, LAMMPS requires :doc:`newton <newton>` bond off such that all
processors containing an atom know when a bond breaks. Additionally,
one must do either (A) or (B).
A) Use the following special bond settings
.. code-block:: LAMMPS
special_bonds lj 0 1 1 coul 1 1 1
These settings accomplish two goals. First, they turn off 1-3 and 1-4
special bond lists, which are not currently supported for BPMs. As
BPMs often have dense bond networks, generating 1-3 and 1-4 special
bond lists is expensive. By setting the lj weight for 1-2 bonds to
zero, this turns off pairwise interactions. Even though there are no
charges in BPM models, setting a nonzero coul weight for 1-2 bonds
ensures all bonded neighbors are still included in the neighbor list
in case bonds break between neighbor list builds.
B) Alternatively, one can simply overlay pair interactions such that all
bonded particles also feel pair interactions. This can be
accomplished by using the *overlay/pair* keyword present in all bpm
bond styles and by using the following special bond settings
.. code-block:: LAMMPS
special_bonds lj/coul 1 1 1
See the :doc:`Howto <Howto_broken_bonds>` page on broken bonds for
more information.
----------
Currently there are two types of bonds included in the BPM
@ -91,12 +51,6 @@ This also requires a unique integrator :doc:`fix nve/bpm/sphere
<fix_nve_bpm_sphere>` which numerically integrates orientation similar
to :doc:`fix nve/asphere <fix_nve_asphere>`.
To monitor the fracture of bonds in the system, all BPM bond styles
have the ability to record instances of bond breakage to output using
the :doc:`dump local <dump>` command. Additionally, one can use
:doc:`compute nbond/atom <compute_nbond_atom>` to tally the current
number of bonds per atom.
In addition to bond styles, a new pair style :doc:`pair bpm/spring
<pair_bpm_spring>` was added to accompany the bpm/spring bond
style. This pair style is simply a hookean repulsion with similar
@ -104,6 +58,73 @@ velocity damping as its sister bond style.
----------
Bond data can be output using a combination of standard LAMMPS commands.
A list of IDs for bonded atoms can be generated using the
:doc:`compute property/local <compute_property_local>` command.
Various properties of bonds can be computed using the
:doc:`compute bond/local <compute_bond_local>` command. This
command allows one to access data saved to the bond's history
such as the reference length of the bond. More information on
bond history data can be found on the documentation pages for the specific
BPM bond styles. Finally, this data can be output using a :doc:`dump local <dump>`
command. As one may output many columns from the same compute, the
:doc:`dump modify <dump_modify>` *colname* option may be used to provide
more helpful column names. An example of this procedure is found in
/examples/bpm/pour/. External software, such as OVITO, can read these dump
files to render bond data.
----------
As bonds can be broken between neighbor list builds, the
:doc:`special_bonds <special_bonds>` command works differently for BPM
bond styles. There are two possible settings which determine how pair
interactions work between bonded particles. First, one can overlay
pair forces with bond forces such that all bonded particles also
feel pair interactions. This can be accomplished by using the *overlay/pair*
keyword present in all bpm bond styles and by using the following special
bond settings
.. code-block:: LAMMPS
special_bonds lj/coul 1 1 1
Alternatively, one can turn off all pair interactions between bonded
particles. Unlike :doc:`bond quartic <bond_quartic>`, this is not done
by subtracting pair forces during the bond computation but rather by
dynamically updating the special bond list. This is the default behavior
of BPM bond styles and is done by updating the 1-2 special bond list as
bonds break. To do this, LAMMPS requires :doc:`newton <newton>` bond off
such that all processors containing an atom know when a bond breaks.
Additionally, one must use the following special bond settings
.. code-block:: LAMMPS
special_bonds lj 0 1 1 coul 1 1 1
These settings accomplish two goals. First, they turn off 1-3 and 1-4
special bond lists, which are not currently supported for BPMs. As
BPMs often have dense bond networks, generating 1-3 and 1-4 special
bond lists is expensive. By setting the lj weight for 1-2 bonds to
zero, this turns off pairwise interactions. Even though there are no
charges in BPM models, setting a nonzero coul weight for 1-2 bonds
ensures all bonded neighbors are still included in the neighbor list
in case bonds break between neighbor list builds.
To monitor the fracture of bonds in the system, all BPM bond styles
have the ability to record instances of bond breakage to output using
the :doc:`dump local <dump>` command. Since one may frequently output
a list of broken bonds and the time they broke, the
:doc:`dump modify <dump_modify>` option *header no* may be useful to
avoid repeatedly printing the header of the dump file. An example of
this procedure is found in /examples/bpm/impact/. Additionally,
one can use :doc:`compute nbond/atom <compute_nbond_atom>` to tally the
current number of bonds per atom.
See the :doc:`Howto <Howto_broken_bonds>` page on broken bonds for
more information.
----------
While LAMMPS has many utilities to create and delete bonds, *only*
the following are currently compatible with BPM bond styles:

View File

@ -1,7 +1,7 @@
Type labels
===========
.. versionadded:: TBD
.. versionadded:: 15Sep2022
Each atom in LAMMPS has an associated numeric atom type. Similarly,
each bond, angle, dihedral, and improper is assigned a bond type,

View File

@ -24,7 +24,7 @@ Examples
Description
"""""""""""
.. versionadded:: TBD
.. versionadded:: 15Sep2022
The *mesocnt* angle style uses the potential

View File

@ -138,15 +138,14 @@ the *overlay/pair* keyword. These settings require specific
restrictions. Further details can be found in the `:doc: how to
<Howto_BPM>` page on BPMs.
If the *store/local* keyword is used, this fix will track bonds that
If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed
and transferred to an internal fix labeled *fix_ID*. This allows the
local data to be accessed by other LAMMPS commands.
Following any optional keyword/value arguments, a list of one or more
attributes is specified. These include the IDs of the two atoms in
the bond. The other attributes for the two atoms include the timestep
during which the bond broke and the current/initial center of mass
position of the two atoms.
local data to be accessed by other LAMMPS commands. Following this optional
keyword, a list of one or more attributes is specified. These include the
IDs of the two atoms in the bond. The other attributes for the two atoms
include the timestep during which the bond broke and the current/initial
center of mass position of the two atoms.
Data is continuously accumulated over intervals of *N*
timesteps. At the end of each interval, all of the saved accumulated
@ -177,29 +176,38 @@ Restart and other info
This bond style writes the reference state of each bond to
:doc:`binary restart files <restart>`. Loading a restart file will
properly resume bonds.
properly resume bonds. However, the reference state is NOT
written to data files. Therefore reading a data file will not
restore bonds and will cause their reference states to be redefined.
The single() function of these pair styles returns 0.0 for the energy
of a pairwise interaction, since energy is not conserved in these
dissipative potentials. It also returns only the normal component of
the pairwise interaction force.
The accumulated data is not written to restart files and should be
output before a restart file is written to avoid missing data.
The internal fix calculates a local vector or local array depending on the
number of input values. The length of the vector or number of rows in
the array is the number of recorded, lost interactions. If a single
input is specified, a local vector is produced. If two or more inputs
are specified, a local array is produced where the number of columns =
the number of inputs. The vector or array can be accessed by any
command that uses local values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
If the *store/local* option is used, an internal fix will calculate
a local vector or local array depending on the number of input values.
The length of the vector or number of rows in the array is the number
of recorded, broken bonds. If a single input is specified, a local
vector is produced. If two or more inputs are specified, a local array
is produced where the number of columns = the number of inputs. The
vector or array can be accessed by any command that uses local values
from a compute as input. See the :doc:`Howto output <Howto_output>` page
for an overview of LAMMPS output options.
The vector or array will be floating point values that correspond to
the specified attribute.
The single() function of this bond style returns 0.0 for the energy
of a bonded interaction, since energy is not conserved in these
dissipative potentials. It also returns only the normal component of
the bonded interaction force. However, the single() function also
calculates 7 extra bond quantities. The first 4 are data from the
reference state of the bond including the initial distance between particles
:math:`r_0` followed by the :math:`x`, :math:`y`, and :math:`z` components
of the initial unit vector pointing to particle I from particle J. The next 3
quantities (5-7) are the :math:`x`, :math:`y`, and :math:`z` components
of the total force, including normal and tangential contributions, acting
on particle I.
These extra quantities can be accessed by the :doc:`compute bond/local <compute_bond_local>`
command, as *b1*, *b2*, ..., *b7*\ .
Restrictions
""""""""""""

View File

@ -103,15 +103,14 @@ the *overlay/pair* keyword. These settings require specific
restrictions. Further details can be found in the `:doc: how to
<Howto_BPM>` page on BPMs.
If the *store/local* keyword is used, this fix will track bonds that
If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed
and transferred to an internal fix labeled *fix_ID*. This allows the
local data to be accessed by other LAMMPS commands.
Following any optional keyword/value arguments, a list of one or more
attributes is specified. These include the IDs of the two atoms in
the bond. The other attributes for the two atoms include the timestep
during which the bond broke and the current/initial center of mass
position of the two atoms.
local data to be accessed by other LAMMPS commands. Following this optional
keyword, a list of one or more attributes is specified. These include the
IDs of the two atoms in the bond. The other attributes for the two atoms
include the timestep during which the bond broke and the current/initial
center of mass position of the two atoms.
Data is continuously accumulated over intervals of *N*
timesteps. At the end of each interval, all of the saved accumulated
@ -141,28 +140,30 @@ Restart and other info
This bond style writes the reference state of each bond to
:doc:`binary restart files <restart>`. Loading a restart
file will properly resume bonds.
file will properly restore bonds. However, the reference state is NOT
written to data files. Therefore reading a data file will not
restore bonds and will cause their reference states to be redefined.
The single() function of these pair styles returns 0.0 for the energy
of a pairwise interaction, since energy is not conserved in these
dissipative potentials.
The accumulated data is not written to restart files and should be
output before a restart file is written to avoid missing data.
The internal fix calculates a local vector or local array depending on the
number of input values. The length of the vector or number of rows in
the array is the number of recorded, lost interactions. If a single
input is specified, a local vector is produced. If two or more inputs
are specified, a local array is produced where the number of columns =
the number of inputs. The vector or array can be accessed by any
command that uses local values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
output options.
If the *store/local* option is used, an internal fix will calculate
a local vector or local array depending on the number of input values.
The length of the vector or number of rows in the array is the number
of recorded, broken bonds. If a single input is specified, a local
vector is produced. If two or more inputs are specified, a local array
is produced where the number of columns = the number of inputs. The
vector or array can be accessed by any command that uses local values
from a compute as input. See the :doc:`Howto output <Howto_output>` page
for an overview of LAMMPS output options.
The vector or array will be floating point values that correspond to
the specified attribute.
The single() function of this bond style returns 0.0 for the energy
of a bonded interaction, since energy is not conserved in these
dissipative potentials. The single() function also calculates an
extra bond quantity, the initial distance :math:`r_0`. This
extra quantity can be accessed by the
:doc:`compute bond/local <compute_bond_local>` command as *b1*\ .
Restrictions
""""""""""""

View File

@ -22,7 +22,7 @@ Examples
Description
"""""""""""
.. versionadded:: TBD
.. versionadded:: 15Sep2022
The *mesocnt* bond style is a wrapper for the :doc:`harmonic
<bond_harmonic>` style, and uses the potential

View File

@ -13,7 +13,7 @@ Syntax
* ID, group-ID are documented in :doc:`compute <compute>` command
* bond/local = style name of this compute command
* one or more values may be appended
* value = *dist* or *dx* or *dy* or *dz* or *engpot* or *force* or *fx* or *fy* or *fz* or *engvib* or *engrot* or *engtrans* or *omega* or *velvib* or *v_name*
* value = *dist* or *dx* or *dy* or *dz* or *engpot* or *force* or *fx* or *fy* or *fz* or *engvib* or *engrot* or *engtrans* or *omega* or *velvib* or *v_name* or *bN*
.. parsed-literal::
@ -29,6 +29,7 @@ Syntax
*omega* = magnitude of bond angular velocity
*velvib* = vibrational velocity along the bond length
*v_name* = equal-style variable with name (see below)
*bN* = bond style specific quantities for allowed N values
* zero or more keyword/args pairs may be appended
* keyword = *set*
@ -47,7 +48,7 @@ Examples
compute 1 all bond/local engpot
compute 1 all bond/local dist engpot force
compute 1 all bond/local dist fx fy fz
compute 1 all bond/local dist fx fy fz b1 b2
compute 1 all bond/local dist v_distsq set dist d
@ -147,6 +148,19 @@ those quantities via the :doc:`compute reduce <compute_reduce>` command
with thermo output, and the :doc:`fix ave/histo <fix_ave_histo>`
command will histogram the length\ :math:`^2` values and write them to a file.
A bond style may define additional bond quantities which can be
accessed as *b1* to *bN*, where N is defined by the bond style. Most
bond styles do not define any additional quantities, so N = 0. An
example of ones that do are the :doc:`BPM bond styles <Howto_bpm>`
which store the reference state between two particles. See
individual bond styles for details.
When using *bN* with bond style *hybrid*, the output will be the Nth
quantity from the sub-style that computes the bonded interaction
(based on bond type). If that sub-style does not define a *bN*,
the output will be 0.0. The maximum allowed N is the maximum number
of quantities provided by any sub-style.
----------
The local data stored by this command is generated by looping over all

View File

@ -93,7 +93,9 @@ Restrictions
Related commands
""""""""""""""""
:doc:`compute pe <compute_pe>`, :doc:`compute bond <compute_bond>`
:doc:`compute pe <compute_pe>`, :doc:`compute bond <compute_bond>`,
:doc:`fix pair <fix_pair>`
Default
"""""""

View File

@ -167,7 +167,7 @@ triangular particles and define the corner points of each triangle.
In addition, the various per-atom quantities listed above for specific
packages are only accessible by this command.
.. versionchanged:: TBD
.. versionchanged:: 15Sep2022
The *espin* property was previously called *spin*.

View File

@ -17,7 +17,7 @@ Syntax
* one or more keyword/value pairs may be appended
* these keywords apply to various dump styles
* keyword = *append* or *at* or *balance* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap*
* keyword = *append* or *at* or *balance* or *buffer* or *colname* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap*
.. parsed-literal::
@ -65,6 +65,8 @@ Syntax
*refresh* arg = c_ID = compute ID that supports a refresh operation
*scale* arg = *yes* or *no*
*sfactor* arg = coordinate scaling factor (> 0.0)
*skip* arg = v_name
v_name = variable with name which evaluates to non-zero (skip) or 0
*sort* arg = *off* or *id* or N or -N
off = no sorting of per-atom lines within a snapshot
id = sort per-atom lines by atom ID
@ -179,8 +181,8 @@ extra buffering.
.. versionadded:: 4May2022
The *colname* keyword can be used to change the default header keyword
for dump styles: *atom*, *custom*, and *cfg* and their compressed, ADIOS,
and MPIIO variants. The setting for *ID string* replaces the default
for dump styles: *atom*, *custom*, *cfg*, and *local* and their compressed,
ADIOS, and MPIIO variants. The setting for *ID string* replaces the default
text with the provided string. *ID* can be a positive integer when it
represents the column number counting from the left, a negative integer
when it represents the column number from the right (i.e. -1 is the last
@ -694,14 +696,33 @@ most effective when the typical magnitude of position data is between
----------
.. versionadded:: 15Sep2022
The *skip* keyword can be used with all dump styles. It allows a dump
snapshot to be skipped (not written to the dump file), if a condition
is met. The condition is computed by an :doc:`equal-style variable
<variable>`, which should be specified as v_name, where name is the
variable name. If the variable evaluation returns a non-zero value,
then the dump snapshot is skipped. If it returns zero, the dump
proceeds as usual. Note that :doc:`equal-style variable <variable>`
can contain Boolean operators which effectively evaluate as a true
(non-zero) or false (zero) result.
The *skip* keyword can be useful for debugging purposes, e.g. to dump
only on a particular timestep. Or to limit output to conditions of
interest, e.g. only when the force on some atom exceeds a threshold
value.
----------
The *sort* keyword determines whether lines of per-atom output in a
snapshot are sorted or not. A sort value of *off* means they will
typically be written in indeterminate order, either in serial or
parallel. This is the case even in serial if the
:doc:`atom_modify sort <atom_modify>` option is turned on, which it is by
default, to improve performance. A sort value of *id* means sort the output by
atom ID. A sort value of :math:`N` or :math:`-N` means sort the output by the
value in the :math:`N`\ th column of per-atom info in either ascending or
parallel. This is the case even in serial if the :doc:`atom_modify sort
<atom_modify>` option is turned on, which it is by default, to improve
performance. A sort value of *id* means sort the output by atom ID. A
sort value of :math:`N` or :math:`-N` means sort the output by the value
in the :math:`N`\ th column of per-atom info in either ascending or
descending order.
The dump *local* style cannot be sorted by atom ID, since there are
@ -745,8 +766,8 @@ attributes that can be tested for are the same as those that can be
specified in the :doc:`dump custom <dump>` command, with the exception
of the *element* attribute, since it is not a numeric value. Note
that a different attributes can be used than those output by the
:doc:`dump custom <dump>` command. For example, you can output the coordinates
and stress of atoms whose energy is above some threshold.
:doc:`dump custom <dump>` command. For example, you can output the
coordinates and stress of atoms whose energy is above some threshold.
If an atom-style variable is used as the attribute, then it can
produce continuous numeric values or effective Boolean 0/1 values,

View File

@ -312,6 +312,7 @@ accelerated styles exist.
* :doc:`orient/fcc <fix_orient>` - add grain boundary migration force for FCC
* :doc:`orient/eco <fix_orient_eco>` - add generalized grain boundary migration force
* :doc:`pafi <fix_pafi>` - constrained force averages on hyper-planes to compute free energies (PAFI)
* :doc:`pair <fix_pair>` - access per-atom info from pair styles
* :doc:`phonon <fix_phonon>` - calculate dynamical matrix from MD simulations
* :doc:`pimd <fix_pimd>` - Feynman path integral molecular dynamics
* :doc:`planeforce <fix_planeforce>` - constrain atoms to move in a plane

View File

@ -71,6 +71,15 @@ the polymer chains all become interconnected. For this use case, if
angles are defined they should not include bonds between sticker
sites.
.. note::
For the sticker site model, you should set the newton flag for
bonds to "off", via the :doc:`newton on off<newton>` command ("on"
is the default for the 2nd argument). This is to ensure
appropriate randomness in bond selection because the I,J bond will
be stored by both atom I and atom J. LAMMPS cannot check for this,
because it is not aware that a sticker site model is being used.
----------
The bond swapping operation is invoked once every *Nevery* timesteps.
@ -114,11 +123,11 @@ Boltzmann constant, and T is the current temperature of the system.
.. note::
IMPORTANT: Whether the swap is accepted or rejected, no other swaps
are attempted by this processor on this timestep. No other
eligible 4-tuples of atoms are considered. This means that each
processor will perform either a single swap or none on timesteps
this fix is invoked.
Whether the swap is accepted or rejected, no other swaps are
attempted by this processor on this timestep. No other eligible
4-tuples of atoms are considered. This means that each processor
will perform either a single swap or none on timesteps this fix is
invoked.
----------

110
doc/src/fix_pair.rst Normal file
View File

@ -0,0 +1,110 @@
.. index:: fix pair
fix pair command
=======================
Syntax
""""""
.. parsed-literal::
fix ID group-ID pair N pstyle name flag ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* pair = style name of this fix command
* N = invoke this fix once every N timesteps
* pstyle = name of pair style to extract info from (e.g. eam)
* one or more name/flag pairs can be listed
* name = name of quantity the pair style allows extraction of
* flag = 1 if pair style needs to be triggered to produce data for name, 0 if not
Examples
""""""""
.. code-block:: LAMMPS
fix request all pair 100 eam rho 0
fix request all pair 100 amoeba uind 0 uinp 0
Description
"""""""""""
.. versionadded:: 15Sep2022
Extract per-atom quantities from a pair style and store them in this
fix so they can be accessed by other LAMMPS commands, e.g. by a
:doc:`dump <dump>` command or by another :doc:`fix <fix>`,
:doc:`compute <compute>`, or :doc:`variable <variable>` command.
These are example use cases:
* extract per-atom density from :doc:`pair_style eam <pair_eam>` to a dump file
* extract induced dipoles from :doc:`pair_style amoeba <pair_amoeba>` to a dump file
* extract accuracy metrics from a machine-learned potential to trigger output when
a condition is met (see the :doc:`dump_modify skip <dump_modify>` command)
The *N* argument determines how often the fix is invoked.
The *pstyle* argument is the name of the pair style. It can be a
sub-style used in a :doc:`pair_style hybrid <pair_hybrid>` command.
One or more *name/flag* pairs of arguments follow. Each *name* is a
per-atom quantity which the pair style must recognize as an extraction
request. See the doc pages for individual :doc:`pair_styles
<pair_style>` to see what fix pair requests (if any) they support.
The *flag* setting determines whether this fix will also trigger the
pair style to compute the named quantity so it can be extracted. If the
quantity is always computed by the pair style, no trigger is needed;
specify *flag* = 0. If the quantity is not always computed
(e.g. because it is expensive to calculate), then specify *flag* = 1.
This will trigger the quantity to be calculated only on timesteps it is
needed. Again, see the doc pages for individual :doc:`pair_styles
<pair_style>` to determine which fix pair requests (if any) need to be
triggered with a *flag* = 1 setting.
The per-atom data extracted from the pair style is stored by this fix
as either a per-atom vector or array. If there is only one *name*
argument specified and the pair style computes a single value for each
atom, then this fix stores it as a per-atom vector. Otherwise a
per-atom array is created, with its data in the order of the *name*
arguments.
For example, :doc:`pair_style amoeba <pair_amoeba>` allows extraction of
two named quantities: "uind" and "uinp", both of which are 3-vectors for
each atom, i.e. dipole moments. In the example below a 6-column per-atom
array will be created. Columns 1-3 will store the "uind" values;
columns 4-6 will store the "uinp" values.
.. code-block:: LAMMPS
pair_style amoeba
fix ex all pair amoeba 10 uind 0 uinp 0
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 are
relevant to this fix.
As explained above, this fix produces a per-atom vector or array which
can be accessed by various :doc:`output commands <Howto_output>`. If
an array is produced, the number of columns is the sum of the number
of per-atom quantities produced by each *name* argument requested from
the pair style.
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`compute pair <compute_pair>`
Default
"""""""
none

View File

@ -33,7 +33,7 @@ Examples
Description
"""""""""""
.. versionadded:: TBD
.. versionadded:: 15Sep2022
Define alphanumeric type labels to associate with one or more numeric
atom, bond, angle, dihedral or improper types. A collection of type

View File

@ -156,6 +156,20 @@ settings.
----------
.. versionadded:: TBD
The *amoeba* and *hippo* pair styles support extraction of two per-atom
quantities by the :doc:`fix pair <fix_pair>` command. This allows the
quantities to be output to files by the :doc:`dump <dump>` or otherwise
processed by other LAMMPS commands.
The names of the two quantities are "uind" and "uinp" for the induced
dipole moments for each atom. Neither quantity needs to be triggered by
the :doc:`fix pair <fix_pair>` command in order for these pair styles to
calculate it.
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

View File

@ -59,7 +59,7 @@ include this pair interaction and overlay the pair force over the bond
force or to exclude this pair interaction such that the two particles
only interact via the bond force. See discussion of the *overlay/pair*
option for BPM bond styles and the :doc:`special_bonds <special_bonds>`
command in the `:doc: how to <Howto_BPM>` page on BPMs for more details.
command in the :doc:`how to <Howto_bpm>` page on BPMs for more details.
The following coefficients must be defined for each pair of atom types
via the :doc:`pair_coeff <pair_coeff>` command as in the examples

View File

@ -56,8 +56,10 @@ field. This pairwise thermostat can be used in conjunction with any
:doc:`pair style <pair_style>`, and in leiu of per-particle thermostats
like :doc:`fix langevin <fix_langevin>` or ensemble thermostats like
Nose Hoover as implemented by :doc:`fix nvt <fix_nh>`. To use
*dpd/tstat* as a thermostat for another pair style, use the :doc:`pair_style hybrid/overlay <pair_hybrid>` command to compute both the desired
pair interaction and the thermostat for each pair of particles.
*dpd/tstat* as a thermostat for another pair style, use the
:doc:`pair_style hybrid/overlay <pair_hybrid>` command to compute both
the desired pair interaction and the thermostat for each pair of
particles.
For style *dpd*, the force on atom I due to atom J is given as a sum
of 3 terms
@ -68,29 +70,30 @@ of 3 terms
F^C = & A w(r) \\
F^D = & - \gamma w^2(r) (\hat{r_{ij}} \bullet \vec{v_{ij}}) \\
F^R = & \sigma w(r) \alpha (\Delta t)^{-1/2} \\
w(r) = & 1 - r/r_c
w(r) = & 1 - \frac{r}{r_c}
where :math:`F^C` is a conservative force, :math:`F^D` is a dissipative
force, and :math:`F^R` is a random force. :math:`r_{ij}` is a unit
vector in the direction :math:`r_i - r_j`, :math:`v_{ij}` is the vector
difference in velocities of the two atoms :math:`= \vec{v}_i -
\vec{v}_j`, :math:`\alpha` is a Gaussian random number with zero mean and
unit variance, dt is the timestep size, and w(r) is a weighting factor
that varies between 0 and 1. :math:`r_c` is the cutoff. :math:`\sigma`
is set equal to :math:`\sqrt{2 k_B T \gamma}`, where :math:`k_B` is the
Boltzmann constant and T is the temperature parameter in the pair_style
command.
force, and :math:`F^R` is a random force. :math:`\hat{r_{ij}}` is a
unit vector in the direction :math:`r_i - r_j`, :math:`\vec{v_{ij}}` is
the vector difference in velocities of the two atoms :math:`\vec{v}_i -
\vec{v}_j`, :math:`\alpha` is a Gaussian random number with zero mean
and unit variance, *dt* is the timestep size, and :math:`w(r)` is a
weighting factor that varies between 0 and 1. :math:`r_c` is the
pairwise cutoff. :math:`\sigma` is set equal to :math:`\sqrt{2 k_B T
\gamma}`, where :math:`k_B` is the Boltzmann constant and *T* is the
temperature parameter in the pair_style command.
For style *dpd/tstat*, the force on atom I due to atom J is the same
as the above equation, except that the conservative Fc term is
dropped. Also, during the run, T is set each timestep to a ramped
value from Tstart to Tstop.
For style *dpd/tstat*, the force on atom I due to atom J is the same as
the above equation, except that the conservative :math:`F^C` term is
dropped. Also, during the run, *T* is set each timestep to a ramped
value from *Tstart* to *Tstop*.
For style *dpd*, the pairwise energy associated with style *dpd* is
only due to the conservative force term Fc, and is shifted to be zero
at the cutoff distance Rc. The pairwise virial is calculated using
all 3 terms. For style *dpd/tstat* there is no pairwise energy, but
the last two terms of the formula make a contribution to the virial.
For style *dpd*, the pairwise energy associated with style *dpd* is only
due to the conservative force term :math:`F^C`, and is shifted to be
zero at the cutoff distance :math:`r_c`. The pairwise virial is
calculated using all 3 terms. For style *dpd/tstat* there is no
pairwise energy, but the last two terms of the formula make a
contribution to the virial.
For style *dpd*, the following coefficients must be defined for each
pair of atoms types via the :doc:`pair_coeff <pair_coeff>` command as in
@ -146,8 +149,8 @@ I,J pairs must be specified explicitly.
These pair styles do not support the :doc:`pair_modify <pair_modify>`
shift option for the energy of the pair interaction. Note that as
discussed above, the energy due to the conservative Fc term is already
shifted to be 0.0 at the cutoff distance Rc.
discussed above, the energy due to the conservative :math:`F^C` term is already
shifted to be 0.0 at the cutoff distance :math:`r_c`.
The :doc:`pair_modify <pair_modify>` table option is not relevant
for these pair styles.

View File

@ -58,32 +58,27 @@ given as a sum of 3 terms
F^C = & A w(r) \\
F^D = & - \gamma w^2(r) (\hat{r_{ij}} \bullet \vec{v_{ij}}) \\
F^R = & \sigma w(r) \alpha (\Delta t)^{-1/2} \\
w(r) = & 1 - r/r_c
w(r) = & 1 - \frac{r}{r_c}
where :math:`F^C` is a conservative force, :math:`F^D` is a dissipative
force, and :math:`F^R` is a random force. :math:`r_{ij}` is a unit
vector in the direction :math:`r_i - r_j`, :math:`V_{ij} is the vector
difference in velocities of the two atoms :math:`= \vec{v}_i -
\vec{v}_j, :math:`\alpha` is a Gaussian random number with zero mean and
unit variance, dt is the timestep size, and w(r) is a weighting factor
that varies between 0 and 1. Rc is the cutoff. The weighting factor,
:math:`\omega_{ij}`, varies between 0 and 1, and is chosen to have the
following functional form:
force, and :math:`F^R` is a random force. :math:`\hat{r_{ij}}` is a
unit vector in the direction :math:`r_i - r_j`, :math:`\vec{v_{ij}}` is
the vector difference in velocities of the two atoms, :math:`\vec{v}_i -
\vec{v}_j`, :math:`\alpha` is a Gaussian random number with zero mean
and unit variance, *dt* is the timestep size, and :math:`w(r)` is a
weighting factor that varies between 0 and 1, :math:`r_c` is the
pairwise cutoff. Note that alternative definitions of the weighting
function exist, but would have to be implemented as a separate pair
style command.
.. math::
\omega_{ij} = 1 - \frac{r_{ij}}{r_{c}}
Note that alternative definitions of the weighting function exist, but
would have to be implemented as a separate pair style command.
For style *dpd/fdt*, the fluctuation-dissipation theorem defines :math:`\gamma`
to be set equal to :math:`\sigma^2/(2 T)`, where T is the set point
temperature specified as a pair style parameter in the above examples.
The following coefficients must be defined for each pair of atoms types
via the :doc:`pair_coeff <pair_coeff>` command as in the examples above,
or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>` commands:
For style *dpd/fdt*, the fluctuation-dissipation theorem defines
:math:`\gamma` to be set equal to :math:`\sigma^2/(2 T)`, where *T* is the
set point temperature specified as a pair style parameter in the above
examples. The following coefficients must be defined for each pair of
atoms types via the :doc:`pair_coeff <pair_coeff>` command as in the
examples above, or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands:
* A (force units)
* :math:`\sigma` (force\*time\^(1/2) units)
@ -94,9 +89,9 @@ cutoff is used.
Style *dpd/fdt/energy* is used to perform DPD simulations under
isoenergetic and isoenthalpic conditions. The fluctuation-dissipation
theorem defines :math:`\gamma` to be set equal to :math:`sigma^2/(2
\theta)`, where :math:theta` is the average internal temperature for the
pair. The particle internal temperature is related to the particle
theorem defines :math:`\gamma` to be set equal to :math:`\sigma^2/(2
\theta)`, where :math:`\theta` is the average internal temperature for
the pair. The particle internal temperature is related to the particle
internal energy through a mesoparticle equation of state (see :doc:`fix
eos <fix>`). The differential internal conductive and mechanical
energies are computed within style *dpd/fdt/energy* as:
@ -116,15 +111,15 @@ where
\sigma^{2}_{ij} = & 2\gamma_{ij}k_{B}\Theta_{ij} \\
\Theta_{ij}^{-1} = & \frac{1}{2}(\frac{1}{\theta_{i}}+\frac{1}{\theta_{j}})
:math:`\zeta_ij^q` is a second Gaussian random number with zero mean and unit
variance that is used to compute the internal conductive energy. The
fluctuation-dissipation theorem defines :math:`alpha^2` to be set
equal to :math:2k_B\kappa`, where :math:`\kappa` is the mesoparticle thermal
conductivity parameter. The following coefficients must be defined for
each pair of atoms types via the :doc:`pair_coeff <pair_coeff>`
command as in the examples above, or in the data file or restart files
read by the :doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands:
:math:`\zeta_ij^q` is a second Gaussian random number with zero mean and
unit variance that is used to compute the internal conductive
energy. The fluctuation-dissipation theorem defines :math:`alpha^2` to
be set equal to :math:`2k_B\kappa`, where :math:`\kappa` is the
mesoparticle thermal conductivity parameter. The following coefficients
must be defined for each pair of atoms types via the :doc:`pair_coeff
<pair_coeff>` command as in the examples above, or in the data file or
restart files read by the :doc:`read_data <read_data>` or
:doc:`read_restart <read_restart>` commands:
* A (force units)
* :math:`\sigma` (force\*time\^(1/2) units)
@ -135,23 +130,23 @@ The last coefficient is optional. If not specified, the global DPD
cutoff is used.
The pairwise energy associated with styles *dpd/fdt* and
*dpd/fdt/energy* is only due to the conservative force term Fc, and is
shifted to be zero at the cutoff distance Rc. The pairwise virial is
calculated using only the conservative term.
*dpd/fdt/energy* is only due to the conservative force term :math:`F^C`,
and is shifted to be zero at the cutoff distance :math:`r_c`. The
pairwise virial is calculated using only the conservative term.
The forces computed through the *dpd/fdt* and *dpd/fdt/energy* styles
can be integrated with the velocity-Verlet integration scheme or the
Shardlow splitting integration scheme described by :ref:`(Lisal) <Lisal3>`.
In the cases when these pair styles are combined with the
Shardlow splitting integration scheme described by :ref:`(Lisal)
<Lisal3>`. In the cases when these pair styles are combined with the
:doc:`fix shardlow <fix_shardlow>`, these pair styles differ from the
other dpd styles in that the dissipative and random forces are split
from the force calculation and are not computed within the pair style.
Thus, only the conservative force is computed by the pair style,
while the stochastic integration of the dissipative and random forces
are handled through the Shardlow splitting algorithm approach. The
Shardlow splitting algorithm is advantageous, especially when
performing DPD under isoenergetic conditions, as it allows
significantly larger timesteps to be taken.
Thus, only the conservative force is computed by the pair style, while
the stochastic integration of the dissipative and random forces are
handled through the Shardlow splitting algorithm approach. The Shardlow
splitting algorithm is advantageous, especially when performing DPD
under isoenergetic conditions, as it allows significantly larger
timesteps to be taken.
----------
@ -162,8 +157,9 @@ significantly larger timesteps to be taken.
Restrictions
""""""""""""
These commands are part of the DPD-REACT package. They are only
enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
These commands are part of the DPD-REACT package. They are only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Pair styles *dpd/fdt* and *dpd/fdt/energy* require use of the
:doc:`comm_modify vel yes <comm_modify>` option so that velocities are

View File

@ -444,6 +444,20 @@ identical to the FS EAM files (see above).
----------
.. versionadded:: TBD
The *eam*, *eam/alloy*, *eam/fs*, and *eam/he* pair styles support
extraction of two per-atom quantities by the :doc:`fix pair <fix_pair>`
command. This allows the quantities to be output to files by the
:doc:`dump <dump>` or otherwise processed by other LAMMPS commands.
The names of the two quantities are "rho" and "fp" for the density and
derivative of the embedding energy for each atom. Neither quantity
needs to be triggered by the :doc:`fix pair <fix_pair>` command in order
for these pair styles to calculate it.
----------
.. include:: accel_styles.rst
----------
@ -459,21 +473,26 @@ a pair_coeff command with I != J arguments for the eam styles.
This pair style does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options.
The eam pair styles do not write their information to :doc:`binary restart files <restart>`, since it is stored in tabulated potential files.
Thus, you need to re-specify the pair_style and pair_coeff commands in
an input script that reads a restart file.
The eam pair styles do not write their information to :doc:`binary
restart files <restart>`, since it is stored in tabulated potential
files. Thus, you need to re-specify the pair_style and pair_coeff
commands in an input script that reads a restart file.
The eam pair styles can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. They do not support the
*inner*, *middle*, *outer* keywords.
----------
Restrictions
""""""""""""
All of these styles are part of the MANYBODY package. They are only
enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
enabled if LAMMPS was built with that package. See the :doc:`Build
package <Build_package>` page for more info.
Related commands
""""""""""""""""

View File

@ -48,7 +48,7 @@ in LAMMPS. For the exact functional form of the potential and
implementation details, the reader is referred to the original papers
:ref:`(Volkov1) <Volkov1>` and :ref:`(Volkov2) <Volkov2>`.
.. versionchanged:: TBD
.. versionchanged:: 15Sep2022
The potential supports two modes, *segment* and *chain*. By default,
*chain* mode is enabled. In *segment* mode, interactions are
@ -84,7 +84,7 @@ is faster and is enabled by default.
IDs. If this is not possible (e.g. in simulations of CNT rings),
*topology* mode needs to be enabled in the pair_style command.
.. versionadded:: TBD
.. versionadded:: 15Sep2022
In addition to the LJ interactions described above, style
*mesocnt/viscous* explicitly models friction between neighboring

View File

@ -282,25 +282,25 @@ the orientation of a particular atom is the same, regardless of how
many processors are being used. This keyword does not allow use of an
atom-style variable.
.. versionchanged:: TBD
.. versionchanged:: 15Sep2022
Keyword *spin/atom* uses the specified g value to set the magnitude of the
magnetic spin vectors, and the x,y,z values as components of a vector
to set as the orientation of the magnetic spin vectors of the selected
atoms. This keyword was previously called *spin*.
.. versionchanged:: TBD
.. versionchanged:: 15Sep2022
Keyword *spin/atom/random* randomizes the orientation of the magnetic spin
vectors for the selected atoms and sets the magnitude of each to the
specified *Dlen* value. This keyword was previously called *spin/random*.
.. versionadded:: TBD
.. versionadded:: 15Sep2022
Keyword *radius/electron* uses the specified value to set the radius of
electrons or fixed cores.
.. versionadded:: TBD
.. versionadded:: 15Sep2022
Keyword *spin/electron* sets the spin of an electron (+/- 1) or indicates
nuclei (=0), fixed-cores (=2), or pseudo-cores (= 3).

View File

@ -3581,7 +3581,9 @@ UEF
ufm
Uhlenbeck
Ui
uind
uInfParallel
uinp
uk
ul
ulb

File diff suppressed because it is too large Load Diff

View File

@ -1,219 +0,0 @@
LAMMPS (17 Feb 2022)
units lj
dimension 3
boundary f f f
atom_style bpm/sphere
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
newton on off
comm_modify vel yes cutoff 2.6
lattice fcc 1.0
Lattice spacing in x,y,z = 1.5874011 1.5874011 1.5874011
region box block -25 15 -22 22 -22 22
create_box 1 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50
Created orthogonal box = (-39.685026 -34.922823 -34.922823) to (23.811016 34.922823 34.922823)
1 by 2 by 2 MPI processor grid
region disk cylinder x 0.0 0.0 20.0 -0.5 0.5
create_atoms 1 region disk
Created 7529 atoms
using lattice units in orthogonal box = (-39.685026 -34.922823 -34.922823) to (23.811016 34.922823 34.922823)
create_atoms CPU = 0.002 seconds
group plate region disk
7529 atoms in group plate
region ball sphere 8.0 0.0 0.0 6.0
create_atoms 1 region ball
Created 3589 atoms
using lattice units in orthogonal box = (-39.685026 -34.922823 -34.922823) to (23.811016 34.922823 34.922823)
create_atoms CPU = 0.001 seconds
group projectile region ball
3589 atoms in group projectile
displace_atoms all random 0.1 0.1 0.1 134598738
Displacing atoms ...
neighbor 1.0 bin
pair_style gran/hooke/history 1.0 NULL 0.5 NULL 0.1 1
pair_coeff 1 1
fix 1 all nve/bpm/sphere
create_bonds many plate plate 1 0.0 1.5
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2.6
binsize = 1, bins = 64 70 70
2 neighbor lists, perpetual/occasional/extra = 1 1 0
(1) command create_bonds, occasional
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
(2) pair gran/hooke/history, perpetual
attributes: half, newton on, size, history
pair build: half/size/bin/newton
stencil: half/bin/3d
bin: standard
Added 38559 bonds, new total = 38559
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
15 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.001 seconds
create_bonds many projectile projectile 2 0.0 1.5
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
WARNING: Bonds are defined but no bond style is set (../force.cpp:192)
WARNING: Likewise 1-2 special neighbor interactions != 1.0 (../force.cpp:194)
Added 21869 bonds, new total = 60428
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
16 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.001 seconds
neighbor 0.3 bin
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
bond_style bpm/rotational store/local brkbond 100 time id1 id2
bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.2 0.002 0.002
bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.2 0.002 0.002
velocity projectile set -0.05 0.0 0.0
compute nbond all nbond/atom
compute tbond all reduce sum c_nbond
timestep 0.05
thermo_style custom step ke pe pxx pyy pzz c_tbond
thermo 100
thermo_modify lost ignore lost/bond ignore
#dump 1 all custom 100 atomDump id radius x y z c_nbond
dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
dump_modify 2 header no
run 7500
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 1.3
ghost atom cutoff = 2.6
binsize = 0.65, bins = 98 108 108
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair gran/hooke/history, perpetual
attributes: half, newton on, size, history
pair build: half/size/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 33.34 | 33.34 | 33.35 Mbytes
Step KinEng PotEng Pxx Pyy Pzz c_tbond
0 0.00053238861 0 3.8217307e-05 0 0 10.8703
100 0.00053238861 0 3.8217307e-05 1.2166373e-20 1.2308212e-20 10.8703
200 0.00053238861 0 3.8217307e-05 1.2454467e-20 1.2479589e-20 10.8703
300 0.00053238861 0 3.8217307e-05 1.2256702e-20 1.2621091e-20 10.8703
400 0.00053238861 0 3.8217307e-05 1.2170534e-20 1.2751666e-20 10.8703
500 0.00053093549 0 3.8484194e-05 7.645531e-08 1.4861825e-07 10.8703
600 0.00051485902 0 4.0340751e-05 3.8615876e-07 4.8663463e-07 10.869221
700 0.00049942978 0 3.8684008e-05 4.5363318e-07 4.4560229e-07 10.85501
800 0.00049465262 0 3.6604612e-05 1.4755468e-07 2.0062093e-07 10.820651
900 0.00048784775 0 3.5333139e-05 3.5118328e-07 6.6697625e-07 10.769563
1000 0.00048345699 0 3.4702137e-05 7.0312998e-07 4.0218318e-07 10.730347
1100 0.00047945073 0 3.5065961e-05 6.2813891e-07 7.4640359e-07 10.703184
1200 0.00047512604 0 3.4833144e-05 8.5208604e-07 8.7277583e-07 10.686634
1300 0.00047401428 0 3.4236869e-05 1.0321827e-06 7.4545242e-07 10.678
1400 0.00047619121 0 3.4416549e-05 8.7518021e-07 7.3979503e-07 10.671704
1500 0.0004668728 0 3.4495751e-05 1.4077823e-06 1.517373e-06 10.666127
1600 0.00045088371 0 3.3264301e-05 1.8499661e-06 1.9842415e-06 10.66073
1700 0.00044275099 0 3.2471064e-05 1.9028747e-06 2.2248947e-06 10.6458
1800 0.0004424362 0 3.1846336e-05 1.6208492e-06 1.9291602e-06 10.620615
1900 0.00043678957 0 3.1260936e-05 1.4673956e-06 1.6120523e-06 10.603166
2000 0.00042747562 0 3.0652107e-05 1.6455486e-06 1.53127e-06 10.576003
2100 0.0004214344 0 3.0240727e-05 1.8873967e-06 1.5258622e-06 10.539845
2200 0.00041712779 0 3.0329566e-05 1.8846152e-06 1.4971471e-06 10.49937
2300 0.00041095769 0 3.0000572e-05 2.3585924e-06 1.6773177e-06 10.471668
2400 0.00040883568 0 2.9625158e-05 1.9105554e-06 1.8720763e-06 10.45116
2500 0.00040762685 0 2.9441541e-05 1.6848938e-06 1.8877532e-06 10.437309
2600 0.00040579873 0 2.9255988e-05 1.7523874e-06 1.636423e-06 10.422378
2700 0.00040340975 0 2.9035693e-05 1.673158e-06 1.9038932e-06 10.410505
2800 0.00040170914 0 2.8829361e-05 1.6711978e-06 1.9776001e-06 10.400792
2900 0.00040015113 0 2.8614186e-05 1.5982427e-06 1.7994733e-06 10.393416
3000 0.00040029253 0 2.8470718e-05 1.5589166e-06 1.6682302e-06 10.385321
3100 0.00040037329 0 2.8483376e-05 1.2831526e-06 1.4788005e-06 10.378485
3200 0.00040142612 0 2.8481287e-05 1.1577988e-06 1.3495778e-06 10.373988
3300 0.00040105092 0 2.8547009e-05 1.2155138e-06 1.2633439e-06 10.370031
3400 0.00039950673 0 2.8340939e-05 1.1182251e-06 1.1624668e-06 10.364274
3500 0.00039715236 0 2.824813e-05 1.3086462e-06 1.2029185e-06 10.360496
3600 0.00039446552 0 2.8112283e-05 1.1232321e-06 1.0077217e-06 10.353121
3700 0.00039263296 0 2.7927975e-05 1.1083636e-06 1.2091857e-06 10.346645
3800 0.00039061341 0 2.7819957e-05 1.1836841e-06 1.3566272e-06 10.341069
3900 0.00038985051 0 2.7681947e-05 1.3588359e-06 1.4099727e-06 10.329196
4000 0.00038815347 0 2.7492102e-05 1.1111719e-06 1.1700718e-06 10.318043
4100 0.00038651302 0 2.7444105e-05 9.9563429e-07 1.4085969e-06 10.311027
4200 0.00038565809 0 2.7177341e-05 9.5736307e-07 1.0404482e-06 10.299155
4300 0.0003847255 0 2.7029216e-05 9.6204756e-07 1.140804e-06 10.292319
4400 0.0003844421 0 2.6841047e-05 9.6570404e-07 1.2319818e-06 10.286203
4500 0.0003842788 0 2.6633558e-05 9.6452478e-07 1.1954945e-06 10.278287
4600 0.00038365139 0 2.6514403e-05 9.6185846e-07 1.2002452e-06 10.270732
4700 0.00038271503 0 2.6374349e-05 9.4061833e-07 1.1774211e-06 10.264796
4800 0.00038233688 0 2.638398e-05 1.1644119e-06 1.3746239e-06 10.25742
4900 0.00038223496 0 2.6279821e-05 1.1345508e-06 1.4709213e-06 10.246987
5000 0.00038219402 0 2.6188871e-05 1.0115151e-06 1.2024203e-06 10.240511
5100 0.00038195153 0 2.6137945e-05 1.009856e-06 1.1961088e-06 10.236014
5200 0.00038170903 0 2.6103563e-05 1.0046761e-06 1.1881008e-06 10.232236
5300 0.00038194303 0 2.6111938e-05 1.0533375e-06 1.2621634e-06 10.230617
5400 0.00038147407 0 2.6078641e-05 1.082228e-06 1.2915223e-06 10.230098
5500 0.00038156894 0 2.6084488e-05 1.1395485e-06 1.3592644e-06 10.227759
5600 0.00038169434 0 2.6085704e-05 1.1173618e-06 1.3003599e-06 10.2256
5700 0.00038219734 0 2.6095279e-05 1.1026614e-06 1.280455e-06 10.223621
5800 0.00038268758 0 2.6113437e-05 1.1096198e-06 1.2565503e-06 10.222902
5900 0.00038300658 0 2.6131709e-05 1.1123595e-06 1.235992e-06 10.222182
6000 0.00038250316 0 2.606995e-05 1.1590744e-06 1.2888416e-06 10.221123
6100 0.0003821526 0 2.6025605e-05 1.1434025e-06 1.3141861e-06 10.219503
6200 0.00038185711 0 2.5991255e-05 1.1471391e-06 1.3427373e-06 10.219503
6300 0.00038197679 0 2.5996965e-05 1.1338082e-06 1.3315258e-06 10.218604
6400 0.00038232311 0 2.6035805e-05 1.1353407e-06 1.3306683e-06 10.217884
6500 0.00038255543 0 2.6091572e-05 1.1768703e-06 1.3629611e-06 10.217704
6600 0.00038251887 0 2.6068968e-05 1.1808094e-06 1.3969697e-06 10.217344
6700 0.00038177389 0 2.6004288e-05 1.1659866e-06 1.423638e-06 10.218084
6800 0.00038096291 0 2.5969494e-05 1.1377343e-06 1.4348787e-06 10.218103
6900 0.00038090601 0 2.5951546e-05 1.1327767e-06 1.4311663e-06 10.217024
7000 0.00038088094 0 2.5946255e-05 1.1652568e-06 1.4567559e-06 10.215944
7100 0.00038094624 0 2.5972593e-05 1.1558871e-06 1.4692935e-06 10.214684
7200 0.00038168738 0 2.6002e-05 1.1562707e-06 1.4881081e-06 10.212705
7300 0.00038200854 0 2.6038768e-05 1.1339903e-06 1.4808455e-06 10.212345
7400 0.00038187543 0 2.6044759e-05 1.101743e-06 1.4758679e-06 10.213084
7500 0.00038165297 0 2.6004536e-05 1.0892731e-06 1.4872621e-06 10.214742
Loop time of 28.804 on 4 procs for 7500 steps with 11111 atoms
Performance: 1124843.305 tau/day, 260.380 timesteps/s
97.5% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.26977 | 0.28058 | 0.2866 | 1.3 | 0.97
Bond | 22.742 | 23.598 | 24.671 | 16.6 | 81.92
Neigh | 0.54555 | 0.5728 | 0.60272 | 3.2 | 1.99
Comm | 1.4024 | 2.5619 | 3.5079 | 54.8 | 8.89
Output | 0.025307 | 0.025833 | 0.027022 | 0.4 | 0.09
Modify | 1.592 | 1.6506 | 1.7059 | 4.0 | 5.73
Other | | 0.1147 | | | 0.40
Nlocal: 2777.75 ave 2887 max 2682 min
Histogram: 1 0 0 0 2 0 0 0 0 1
Nghost: 1152.5 ave 1189 max 1125 min
Histogram: 1 0 1 0 0 1 0 0 0 1
Neighs: 11515.5 ave 12520 max 10831 min
Histogram: 1 1 0 0 1 0 0 0 0 1
Total # of neighbors = 46062
Ave neighs/atom = 4.1456215
Ave special neighs/atom = 10.214742
Neighbor list builds = 408
Dangerous builds = 0
Total wall time: 0:00:28

View File

@ -0,0 +1,219 @@
LAMMPS (4 May 2022)
units lj
dimension 3
boundary f f f
atom_style bpm/sphere
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
newton on off
comm_modify vel yes cutoff 2.6
lattice fcc 1.0
Lattice spacing in x,y,z = 1.5874011 1.5874011 1.5874011
region box block -25 15 -22 22 -22 22
create_box 1 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50
Created orthogonal box = (-39.685026 -34.922823 -34.922823) to (23.811016 34.922823 34.922823)
1 by 2 by 2 MPI processor grid
region disk cylinder x 0.0 0.0 20.0 -0.5 0.5
create_atoms 1 region disk
Created 7529 atoms
using lattice units in orthogonal box = (-39.685026 -34.922823 -34.922823) to (23.811016 34.922823 34.922823)
create_atoms CPU = 0.002 seconds
group plate region disk
7529 atoms in group plate
region ball sphere 8.0 0.0 0.0 6.0
create_atoms 1 region ball
Created 3589 atoms
using lattice units in orthogonal box = (-39.685026 -34.922823 -34.922823) to (23.811016 34.922823 34.922823)
create_atoms CPU = 0.001 seconds
group projectile region ball
3589 atoms in group projectile
displace_atoms all random 0.1 0.1 0.1 134598738
Displacing atoms ...
neighbor 1.0 bin
pair_style gran/hooke/history 1.0 NULL 0.5 NULL 0.1 1
pair_coeff 1 1
fix 1 all nve/bpm/sphere
create_bonds many plate plate 1 0.0 1.5
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2
ghost atom cutoff = 2.6
binsize = 1, bins = 64 70 70
2 neighbor lists, perpetual/occasional/extra = 1 1 0
(1) command create_bonds, occasional
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
(2) pair gran/hooke/history, perpetual
attributes: half, newton on, size, history
pair build: half/size/bin/newton
stencil: half/bin/3d
bin: standard
Added 38559 bonds, new total = 38559
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
15 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.001 seconds
create_bonds many projectile projectile 2 0.0 1.5
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
WARNING: Bonds are defined but no bond style is set (../force.cpp:192)
WARNING: Likewise 1-2 special neighbor interactions != 1.0 (../force.cpp:194)
Added 21869 bonds, new total = 60428
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
16 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.002 seconds
neighbor 0.3 bin
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
bond_style bpm/rotational store/local brkbond 100 time id1 id2
bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.2 0.002 0.002
bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.2 0.002 0.002
velocity projectile set -0.05 0.0 0.0
compute nbond all nbond/atom
compute tbond all reduce sum c_nbond
timestep 0.05
thermo_style custom step ke pe pxx pyy pzz c_tbond
thermo 100
thermo_modify lost ignore lost/bond ignore
#dump 1 all custom 100 atomDump id radius x y z c_nbond
dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
dump_modify 2 header no
run 7500
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 1.3
ghost atom cutoff = 2.6
binsize = 0.65, bins = 98 108 108
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair gran/hooke/history, perpetual
attributes: half, newton on, size, history
pair build: half/size/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 33.22 | 33.22 | 33.22 Mbytes
Step KinEng PotEng Pxx Pyy Pzz c_tbond
0 0.00053238861 0 3.8217307e-05 -7.6534847e-14 1.9906102e-13 10.8703
100 0.00053238861 0 3.8217307e-05 -4.2831948e-13 5.7428612e-13 10.8703
200 0.00053238861 0 3.8217307e-05 -1.4067913e-12 1.4975493e-12 10.8703
300 0.00053238861 0 3.8217307e-05 -8.77782e-13 3.8245851e-13 10.8703
400 0.00053238861 0 3.8217307e-05 -8.5835238e-13 6.5448014e-13 10.8703
500 0.00053093549 0 4.0340893e-05 4.501394e-07 5.3690512e-07 10.8703
600 0.00051485902 0 6.6761034e-05 4.6258948e-06 5.2285428e-06 10.869221
700 0.00049942978 0 9.5499005e-05 9.3031413e-07 4.5389354e-06 10.85501
800 0.00049465262 0 5.6810167e-05 -5.5619903e-06 -1.6010295e-06 10.820651
900 0.00048784775 0 1.5747004e-05 2.0522042e-05 2.5481542e-05 10.769563
1000 0.00048345699 0 2.1159666e-05 1.2232747e-05 1.4767665e-05 10.730347
1100 0.00047945073 0 5.2779833e-05 -2.6136504e-05 -2.7689007e-05 10.703184
1200 0.00047512604 0 6.3234312e-05 -1.7082136e-05 -2.9178979e-05 10.686634
1300 0.00047401428 0 2.5474717e-05 -7.4782876e-06 -1.9808965e-05 10.678
1400 0.00047619121 0 2.5189461e-05 2.9476227e-06 -4.4149511e-06 10.671704
1500 0.0004668728 0 5.8798861e-05 -2.6192143e-06 1.0805172e-06 10.666127
1600 0.00045088371 0 5.8377694e-05 2.2911269e-06 4.339717e-06 10.66073
1700 0.00044275099 0 4.0766018e-05 2.7748031e-05 2.8202527e-05 10.6458
1800 0.0004424362 0 3.0460351e-05 2.9690484e-05 3.3464889e-05 10.620615
1900 0.00043678957 0 3.809598e-05 2.4448755e-06 1.2651201e-05 10.603166
2000 0.00042747562 0 4.2315209e-05 -7.6783024e-06 -3.3357359e-06 10.576003
2100 0.0004214344 0 2.6171505e-05 5.5424035e-06 -4.1153085e-06 10.539845
2200 0.00041712779 0 3.0796803e-05 1.1362383e-05 7.7103875e-07 10.49937
2300 0.00041095769 0 4.141148e-05 1.4142023e-06 -1.0982633e-05 10.471668
2400 0.00040883568 0 3.5835323e-05 -1.0216635e-05 -2.3669763e-05 10.45116
2500 0.00040762685 0 2.879385e-05 -1.4074395e-05 -1.9314875e-05 10.437309
2600 0.00040579873 0 2.771644e-05 -2.316844e-05 -2.2961097e-05 10.422378
2700 0.00040340975 0 3.4482945e-05 -3.075929e-05 -2.3321344e-05 10.410505
2800 0.00040170914 0 3.1147943e-05 -2.4329639e-05 -1.1787678e-05 10.400792
2900 0.00040015113 0 2.3214992e-05 -1.8068374e-05 3.8127871e-06 10.393416
3000 0.00040029253 0 2.7275906e-05 -7.0762689e-06 1.3782424e-05 10.385321
3100 0.00040037329 0 2.962113e-05 -1.3795312e-06 5.3267315e-06 10.378485
3200 0.00040142612 0 2.86096e-05 4.4289601e-06 -3.3950404e-06 10.373988
3300 0.00040105092 0 2.5423615e-05 9.5689359e-06 -2.5296344e-06 10.370031
3400 0.00039950673 0 2.6405258e-05 9.5776239e-06 -7.3789982e-07 10.364274
3500 0.00039715236 0 3.115696e-05 1.0986224e-05 6.6898207e-06 10.360496
3600 0.00039446552 0 2.8426837e-05 9.6296098e-06 1.5057875e-05 10.353121
3700 0.00039263296 0 2.567768e-05 6.347291e-06 2.4842157e-05 10.346645
3800 0.00039061341 0 2.7635193e-05 5.0165135e-06 2.5989901e-05 10.341069
3900 0.00038985051 0 2.8639932e-05 1.056365e-05 2.4463803e-05 10.329196
4000 0.00038815347 0 2.6613146e-05 2.0316396e-05 2.1434368e-05 10.318043
4100 0.00038651302 0 2.4759493e-05 1.9632349e-05 1.3524912e-05 10.311027
4200 0.00038565809 0 2.5290845e-05 1.9908233e-05 8.3895516e-06 10.299155
4300 0.0003847255 0 2.6461182e-05 1.9385332e-05 5.664874e-06 10.292319
4400 0.0003844421 0 2.4464435e-05 1.4675545e-05 6.8223863e-06 10.286203
4500 0.0003842788 0 2.3080348e-05 7.1469159e-06 6.8395849e-06 10.278287
4600 0.00038365139 0 2.4937615e-05 2.3854793e-06 4.6100509e-06 10.270732
4700 0.00038271503 0 2.431281e-05 -3.127707e-06 3.8840673e-07 10.264796
4800 0.00038233688 0 2.3727372e-05 -3.9575833e-06 1.5658614e-06 10.25742
4900 0.00038223496 0 2.3842519e-05 2.6005447e-06 4.5031882e-06 10.246987
5000 0.00038219402 0 2.4705111e-05 8.2018492e-06 4.0121467e-06 10.240511
5100 0.00038195153 0 2.5112089e-05 9.1687723e-06 3.3825795e-06 10.236014
5200 0.00038170903 0 2.4166312e-05 4.6680528e-06 3.0359762e-06 10.232236
5300 0.00038194303 0 2.4293657e-05 3.067111e-06 2.1353614e-06 10.230617
5400 0.00038147407 0 2.472142e-05 5.2915485e-06 -1.7472423e-06 10.230098
5500 0.00038156894 0 2.4839317e-05 6.6216457e-06 -2.7544087e-06 10.227759
5600 0.00038169434 0 2.4600407e-05 3.8679618e-06 1.2622251e-07 10.2256
5700 0.00038219734 0 2.4459589e-05 2.0025919e-07 -1.1917672e-08 10.223621
5800 0.00038268758 0 2.4768428e-05 8.7913744e-07 -1.4000329e-06 10.222902
5900 0.00038300658 0 2.4827866e-05 3.6621944e-06 -2.2178729e-07 10.222182
6000 0.00038250316 0 2.4309432e-05 4.3755483e-06 2.2736019e-06 10.221123
6100 0.0003821526 0 2.4396115e-05 3.3855804e-06 4.4742551e-06 10.219503
6200 0.00038185711 0 2.4795348e-05 1.7569948e-06 4.3229904e-06 10.219503
6300 0.00038197679 0 2.4817115e-05 1.237731e-07 3.9285574e-06 10.218604
6400 0.00038232311 0 2.4723664e-05 -8.7946112e-07 2.6215619e-06 10.217884
6500 0.00038255543 0 2.4821878e-05 -4.8948257e-07 3.9392146e-07 10.217704
6600 0.00038251887 0 2.4923895e-05 -1.1107041e-07 -8.3900047e-07 10.217344
6700 0.00038177389 0 2.4664007e-05 -6.4474733e-07 -2.1004826e-06 10.218084
6800 0.00038096291 0 2.4262293e-05 -1.7159941e-06 -2.8149643e-06 10.218103
6900 0.00038090601 0 2.4179172e-05 -2.2622259e-06 -2.1268271e-06 10.217024
7000 0.00038088094 0 2.4363443e-05 -2.4652531e-06 -8.209416e-07 10.215944
7100 0.00038094624 0 2.5119358e-05 -1.5432706e-06 1.1465135e-06 10.214684
7200 0.00038168738 0 2.5565338e-05 -1.4052753e-07 3.3146851e-06 10.212705
7300 0.00038200854 0 2.5436565e-05 4.353807e-07 3.3504276e-06 10.212345
7400 0.00038187543 0 2.4764819e-05 6.7297502e-07 1.5923471e-06 10.213084
7500 0.00038165297 0 2.4015684e-05 7.8657712e-07 1.0138435e-06 10.214742
Loop time of 32.2292 on 4 procs for 7500 steps with 11111 atoms
Performance: 1005300.744 tau/day, 232.709 timesteps/s
96.2% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.29763 | 0.30464 | 0.31393 | 1.1 | 0.95
Bond | 25.936 | 26.533 | 27.431 | 11.7 | 82.33
Neigh | 0.60911 | 0.63557 | 0.66475 | 2.8 | 1.97
Comm | 1.7247 | 2.7138 | 3.3823 | 40.5 | 8.42
Output | 0.020408 | 0.020935 | 0.02227 | 0.5 | 0.06
Modify | 1.8299 | 1.8724 | 1.9325 | 2.9 | 5.81
Other | | 0.1491 | | | 0.46
Nlocal: 2777.75 ave 2887 max 2682 min
Histogram: 1 0 0 0 2 0 0 0 0 1
Nghost: 1152.5 ave 1189 max 1125 min
Histogram: 1 0 1 0 0 1 0 0 0 1
Neighs: 11515.5 ave 12520 max 10831 min
Histogram: 1 1 0 0 1 0 0 0 0 1
Total # of neighbors = 46062
Ave neighs/atom = 4.1456215
Ave special neighs/atom = 10.214742
Neighbor list builds = 408
Dangerous builds = 0
Total wall time: 0:00:32

View File

@ -19,6 +19,9 @@ bond_coeff 1 1.0 0.2 0.01 0.01 2.0 0.4 0.02 0.02 0.2 0.04 0.002 0.002
compute nbond all nbond/atom
compute tbond all reduce sum c_nbond
compute bond_ids all property/local batom1 batom2
compute bond_properties all bond/local dist b1
compute_modify thermo_temp dynamic/dof yes
fix 1 all wall/gran hertz/history 1.0 NULL 0.5 NULL 0.1 1 zplane 0.0 NULL
@ -31,5 +34,7 @@ timestep 0.05
thermo_style custom step ke pe pxx pyy pzz c_tbond
thermo 100
#dump 1 all custom 500 atomDump id radius x y z c_nbond mol
#dump 2 all local 500 bondDump c_bond_ids[*] c_bond_properties[*]
#dump_modify 2 colname 1 "id1" colname 2 "id2" colname 3 "r" colname 4 "r0"
run 100000

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -2094,6 +2094,28 @@ void *PairAmoeba::extract(const char *str, int &dim)
return nullptr;
}
/* ----------------------------------------------------------------------
peratom requests from FixPair
return ptr to requested data
also return ncol = # of quantites per atom
0 = per-atom vector
1 or more = # of columns in per-atom array
return NULL if str is not recognized
---------------------------------------------------------------------- */
void *PairAmoeba::extract_peratom(const char *str, int &ncol)
{
if (strcmp(str,"uind") == 0) {
ncol = 3;
return (void *) uind;
} else if (strcmp(str,"uinp") == 0) {
ncol = 3;
return (void *) uinp;
}
return nullptr;
}
/* ----------------------------------------------------------------------
grow local vectors and arrays if necessary
keep them all atom->nmax in length even if ghost storage not needed

View File

@ -50,6 +50,7 @@ class PairAmoeba : public Pair {
void unpack_reverse_grid(int, void *, int, int *) override;
void *extract(const char *, int &) override;
void *extract_peratom(const char *, int &) override;
double memory_usage() override;
protected:

View File

@ -14,6 +14,7 @@
#include "bond_bpm.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix_bond_history.h"
@ -35,7 +36,6 @@ BondBPM::BondBPM(LAMMPS *_lmp) :
id_fix_bond_history(nullptr), id_fix_store_local(nullptr), id_fix_prop_atom(nullptr),
fix_store_local(nullptr), fix_bond_history(nullptr), fix_update_special_bonds(nullptr),
pack_choice(nullptr), output_data(nullptr)
{
overlay_flag = 0;
prop_atom_flag = 0;
@ -256,24 +256,42 @@ double BondBPM::equilibrium_distance(int /*i*/)
{
// Ghost atoms may not yet be communicated, this may only be an estimate
if (r0_max_estimate == 0) {
int type, j;
double delx, dely, delz, r;
double **x = atom->x;
for (int i = 0; i < atom->nlocal; i++) {
for (int m = 0; m < atom->num_bond[i]; m++) {
type = atom->bond_type[i][m];
if (type == 0) continue;
if (!fix_bond_history->restart_reset) {
int type, j;
double delx, dely, delz, r;
double **x = atom->x;
for (int i = 0; i < atom->nlocal; i++) {
for (int m = 0; m < atom->num_bond[i]; m++) {
type = atom->bond_type[i][m];
if (type == 0) continue;
j = atom->map(atom->bond_atom[i][m]);
if (j == -1) continue;
j = atom->map(atom->bond_atom[i][m]);
if (j == -1) continue;
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
domain->minimum_image(delx, dely, delz);
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
domain->minimum_image(delx, dely, delz);
r = sqrt(delx * delx + dely * dely + delz * delz);
if (r > r0_max_estimate) r0_max_estimate = r;
r = sqrt(delx * delx + dely * dely + delz * delz);
if (r > r0_max_estimate) r0_max_estimate = r;
}
}
} else {
int type, j;
double r;
for (int i = 0; i < atom->nlocal; i++) {
for (int m = 0; m < atom->num_bond[i]; m++) {
type = atom->bond_type[i][m];
if (type == 0) continue;
j = atom->map(atom->bond_atom[i][m]);
if (j == -1) continue;
// First value must always be reference length
r = fix_bond_history->get_atom_value(i, m, 0);
if (r > r0_max_estimate) r0_max_estimate = r;
}
}
}
@ -286,6 +304,26 @@ double BondBPM::equilibrium_distance(int /*i*/)
return max_stretch * r0_max_estimate / 1.5;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void BondBPM::write_restart(FILE *fp)
{
fwrite(&overlay_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void BondBPM::read_restart(FILE *fp)
{
if (comm->me == 0)
utils::sfread(FLERR, &overlay_flag, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&overlay_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */
void BondBPM::process_broken(int i, int j)

View File

@ -29,9 +29,8 @@ class BondBPM : public Bond {
void init_style() override;
void settings(int, char **) override;
double equilibrium_distance(int) override;
void write_restart(FILE *) override{};
void read_restart(FILE *) override{};
void write_data(FILE *) override{};
void write_restart(FILE *) override;
void read_restart(FILE *) override;
double single(int, double, int, int, double &) override = 0;
protected:

View File

@ -35,26 +35,19 @@ using namespace LAMMPS_NS;
BondBPMRotational::BondBPMRotational(LAMMPS *_lmp) : BondBPM(_lmp)
{
Kr = nullptr;
Ks = nullptr;
Kt = nullptr;
Kb = nullptr;
Fcr = nullptr;
Fcs = nullptr;
Tct = nullptr;
Tcb = nullptr;
gnorm = nullptr;
gslide = nullptr;
groll = nullptr;
gtwist = nullptr;
partial_flag = 1;
smooth_flag = 1;
single_extra = 7;
svector = new double[7];
}
/* ---------------------------------------------------------------------- */
BondBPMRotational::~BondBPMRotational()
{
delete[] svector;
if (allocated) {
memory->destroy(setflag);
memory->destroy(Kr);
@ -190,7 +183,7 @@ void BondBPMRotational::store_data()
2) P. Mora & Y. Wang Advances in Geomcomputing 2009
---------------------------------------------------------------------- */
double BondBPMRotational::elastic_forces(int i1, int i2, int type, double &Fr, double r_mag,
double BondBPMRotational::elastic_forces(int i1, int i2, int type, double r_mag,
double r0_mag, double r_mag_inv, double * /*rhat*/,
double *r, double *r0, double *force1on2,
double *torque1on2, double *torque2on1)
@ -203,7 +196,7 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double &Fr, d
double q1[4], q2[4];
double q2inv[4], mq[4], mqinv[4], qp21[4], q21[4], qtmp[4];
double rb[3], rb_x_r0[3], s[3], t[3];
double Fs[3], Fsp[3], F_rot[3], Ftmp[3];
double Fr, Fs[3], Fsp[3], F_rot[3], Ftmp[3];
double Ts[3], Tb[3], Tt[3], Tbp[3], Ttp[3], Tsp[3], T_rot[3], Ttmp[3];
double **quat = atom->quat;
@ -372,7 +365,7 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double &Fr, d
Note: n points towards 1 vs pointing towards 2
---------------------------------------------------------------------- */
void BondBPMRotational::damping_forces(int i1, int i2, int type, double &Fr, double *rhat,
void BondBPMRotational::damping_forces(int i1, int i2, int type, double *rhat,
double *r, double *force1on2, double *torque1on2,
double *torque2on1)
{
@ -393,7 +386,6 @@ void BondBPMRotational::damping_forces(int i1, int i2, int type, double &Fr, dou
MathExtra::sub3(vn1, vn2, tmp);
MathExtra::scale3(gnorm[type], tmp);
Fr = MathExtra::lensq3(tmp);
MathExtra::add3(force1on2, tmp, force1on2);
// Damp tangential objective velocities
@ -459,7 +451,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
int i1, i2, itmp, n, type;
double r[3], r0[3], rhat[3];
double rsq, r0_mag, r_mag, r_mag_inv;
double Fr, breaking, smooth;
double breaking, smooth;
double force1on2[3], torque1on2[3], torque2on1[3];
ev_init(eflag, vflag);
@ -515,7 +507,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
// Calculate forces, check if bond breaks
// ------------------------------------------------------//
breaking = elastic_forces(i1, i2, type, Fr, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
breaking = elastic_forces(i1, i2, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
torque1on2, torque2on1);
if (breaking >= 1.0) {
@ -524,7 +516,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
continue;
}
damping_forces(i1, i2, type, Fr, rhat, r, force1on2, torque1on2, torque2on1);
damping_forces(i1, i2, type, rhat, r, force1on2, torque1on2, torque2on1);
if (smooth_flag) {
smooth = breaking * breaking;
@ -557,7 +549,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
torque[i2][2] += torque1on2[2] * smooth;
}
if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, Fr * smooth, r[0], r[1], r[2]);
if (evflag) ev_tally_xyz(i1, i2, nlocal, newton_bond, 0.0, -force1on2[0] * smooth, -force1on2[1] * smooth, -force1on2[2] * smooth, r[0], r[1], r[2]);
}
}
@ -668,11 +660,11 @@ void BondBPMRotational::settings(int narg, char **arg)
for (std::size_t i = 0; i < leftover_iarg.size(); i++) {
iarg = leftover_iarg[i];
if (strcmp(arg[iarg], "smooth") == 0) {
if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command");
if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for smooth");
smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
i += 1;
} else {
error->all(FLERR, "Illegal bond_style command");
error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]);
}
}
}
@ -683,6 +675,9 @@ void BondBPMRotational::settings(int narg, char **arg)
void BondBPMRotational::write_restart(FILE *fp)
{
BondBPM::write_restart(fp);
write_restart_settings(fp);
fwrite(&Kr[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Ks[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Kt[1], sizeof(double), atom->nbondtypes, fp);
@ -703,6 +698,8 @@ void BondBPMRotational::write_restart(FILE *fp)
void BondBPMRotational::read_restart(FILE *fp)
{
BondBPM::read_restart(fp);
read_restart_settings(fp);
allocate();
if (comm->me == 0) {
@ -736,14 +733,23 @@ void BondBPMRotational::read_restart(FILE *fp)
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
proc 0 writes to restart file
------------------------------------------------------------------------- */
void BondBPMRotational::write_data(FILE *fp)
void BondBPMRotational::write_restart_settings(FILE *fp)
{
for (int i = 1; i <= atom->nbondtypes; i++)
fprintf(fp, "%d %g %g %g %g %g %g %g %g %g %g %g %g\n", i, Kr[i], Ks[i], Kt[i], Kb[i], Fcr[i],
Fcs[i], Tct[i], Tcb[i], gnorm[i], gslide[i], groll[i], gtwist[i]);
fwrite(&smooth_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void BondBPMRotational::read_restart_settings(FILE *fp)
{
if (comm->me == 0)
utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&smooth_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */
@ -753,11 +759,12 @@ double BondBPMRotational::single(int type, double rsq, int i, int j, double &ffo
// Not yet enabled
if (type <= 0) return 0.0;
int itmp;
int flipped = 0;
if (atom->tag[j] < atom->tag[i]) {
itmp = i;
int itmp = i;
i = j;
j = itmp;
flipped = 1;
}
double r0_mag, r_mag, r_mag_inv;
@ -779,18 +786,39 @@ double BondBPMRotational::single(int type, double rsq, int i, int j, double &ffo
r_mag_inv = 1.0 / r_mag;
MathExtra::scale3(r_mag_inv, r, rhat);
double breaking, smooth, Fr;
double force1on2[3], torque1on2[3], torque2on1[3];
breaking = elastic_forces(i, j, type, Fr, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
double breaking = elastic_forces(i, j, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
torque1on2, torque2on1);
fforce = Fr;
damping_forces(i, j, type, Fr, rhat, r, force1on2, torque1on2, torque2on1);
fforce += Fr;
damping_forces(i, j, type, rhat, r, force1on2, torque1on2, torque2on1);
fforce = MathExtra::dot3(force1on2, r);
fforce *= -1;
double smooth = 1.0;
if (smooth_flag) {
smooth = breaking * breaking;
smooth = 1.0 - smooth * smooth;
fforce *= smooth;
}
// set single_extra quantities
svector[0] = r0_mag;
if (flipped) {
svector[1] = -r0[0];
svector[2] = -r0[1];
svector[3] = -r0[2];
svector[4] = force1on2[0] * smooth;
svector[5] = force1on2[1] * smooth;
svector[6] = force1on2[2] * smooth;
} else {
svector[1] = r0[0];
svector[2] = r0[1];
svector[3] = r0[2];
svector[4] = -force1on2[0] * smooth;
svector[5] = -force1on2[1] * smooth;
svector[6] = -force1on2[2] * smooth;
}
return 0.0;
}

View File

@ -34,7 +34,8 @@ class BondBPMRotational : public BondBPM {
void settings(int, char **) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
double single(int, double, int, int, double &) override;
protected:
@ -44,9 +45,9 @@ class BondBPMRotational : public BondBPM {
double acos_limit(double);
double elastic_forces(int, int, int, double &, double, double, double, double *, double *,
double elastic_forces(int, int, int, double, double, double, double *, double *,
double *, double *, double *, double *);
void damping_forces(int, int, int, double &, double *, double *, double *, double *, double *);
void damping_forces(int, int, int, double *, double *, double *, double *, double *);
void allocate();
void store_data();

View File

@ -34,12 +34,17 @@ BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) :
{
partial_flag = 1;
smooth_flag = 1;
single_extra = 1;
svector = new double[1];
}
/* ---------------------------------------------------------------------- */
BondBPMSpring::~BondBPMSpring()
{
delete[] svector;
if (allocated) {
memory->destroy(setflag);
memory->destroy(k);
@ -291,11 +296,11 @@ void BondBPMSpring::settings(int narg, char **arg)
for (std::size_t i = 0; i < leftover_iarg.size(); i++) {
iarg = leftover_iarg[i];
if (strcmp(arg[iarg], "smooth") == 0) {
if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command");
if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for smooth");
smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
i += 1;
} else {
error->all(FLERR, "Illegal bond_style command");
error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]);
}
}
}
@ -306,6 +311,9 @@ void BondBPMSpring::settings(int narg, char **arg)
void BondBPMSpring::write_restart(FILE *fp)
{
BondBPM::write_restart(fp);
write_restart_settings(fp);
fwrite(&k[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&ecrit[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&gamma[1], sizeof(double), atom->nbondtypes, fp);
@ -317,6 +325,8 @@ void BondBPMSpring::write_restart(FILE *fp)
void BondBPMSpring::read_restart(FILE *fp)
{
BondBPM::read_restart(fp);
read_restart_settings(fp);
allocate();
if (comm->me == 0) {
@ -332,13 +342,23 @@ void BondBPMSpring::read_restart(FILE *fp)
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
proc 0 writes to restart file
------------------------------------------------------------------------- */
void BondBPMSpring::write_data(FILE *fp)
void BondBPMSpring::write_restart_settings(FILE *fp)
{
for (int i = 1; i <= atom->nbondtypes; i++)
fprintf(fp, "%d %g %g %g\n", i, k[i], ecrit[i], gamma[i]);
fwrite(&smooth_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void BondBPMSpring::read_restart_settings(FILE *fp)
{
if (comm->me == 0)
utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&smooth_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */
@ -377,5 +397,9 @@ double BondBPMSpring::single(int type, double rsq, int i, int j, double &fforce)
fforce *= smooth;
}
// set single_extra quantities
svector[0] = r0;
return 0.0;
}

View File

@ -34,7 +34,8 @@ class BondBPMSpring : public BondBPM {
void settings(int, char **) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
double single(int, double, int, int, double &) override;
protected:

View File

@ -33,8 +33,6 @@ ComputeNBondAtom::ComputeNBondAtom(LAMMPS *_lmp, int narg, char **arg) :
peratom_flag = 1;
size_peratom_cols = 0;
peatomflag = 1;
timeflag = 1;
comm_reverse = 1;
nmax = 0;
@ -51,11 +49,6 @@ ComputeNBondAtom::~ComputeNBondAtom()
void ComputeNBondAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
if (update->eflag_atom != invoked_peratom)
error->all(FLERR, "Per-atom nbond was not tallied on needed timestep");
// grow local nbond array if necessary
// needs to be atom->nmax in length

View File

@ -933,5 +933,28 @@ void *PairEAM::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"scale") == 0) return (void *) scale;
return nullptr;
}
/* ----------------------------------------------------------------------
peratom requests from FixPair
return ptr to requested data
also return ncol = # of quantites per atom
0 = per-atom vector
1 or more = # of columns in per-atom array
return NULL if str is not recognized
---------------------------------------------------------------------- */
void *PairEAM::extract_peratom(const char *str, int &ncol)
{
if (strcmp(str,"rho") == 0) {
ncol = 0;
return (void *) rho;
} else if (strcmp(str,"fp") == 0) {
ncol = 0;
return (void *) fp;
}
return nullptr;
}

View File

@ -53,6 +53,7 @@ class PairEAM : public Pair {
double init_one(int, int) override;
double single(int, int, int, int, double, double, double, double &) override;
void *extract(const char *, int &) override;
void *extract_peratom(const char *, int &) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;

View File

@ -50,6 +50,8 @@ static const char cite_fix_bond_swap[] =
" pages = {12718--12728}\n"
"}\n\n";
#define DELTA_PERMUTE 100
/* ---------------------------------------------------------------------- */
FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
@ -92,11 +94,14 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
modify->add_compute(fmt::format("{} all temp",id_temp));
tflag = 1;
// initialize atom list
// initialize two permutation lists
nmax = 0;
alist = nullptr;
maxpermute = 0;
permute = nullptr;
naccept = foursome = 0;
}
@ -109,9 +114,10 @@ FixBondSwap::~FixBondSwap()
// delete temperature if fix created it
if (tflag) modify->delete_compute(id_temp);
delete [] id_temp;
delete[] id_temp;
memory->destroy(alist);
delete[] permute;
}
/* ---------------------------------------------------------------------- */
@ -238,6 +244,15 @@ void FixBondSwap::post_integrate()
memory->create(alist,nmax,"bondswap:alist");
}
// use randomized permutation of both I and J atoms in double loop below
// this is to avoid any bias in accepted MC swaps based on
// ordering LAMMPS creates on a processor for atoms or their neighbors
// create a random permutation of list of Neligible atoms
// uses one-pass Fisher-Yates shuffle on an initial identity permutation
// output: randomized alist[] vector, used in outer loop to select an I atom
// similar randomized permutation is created for neighbors of each I atom
int neligible = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
@ -246,8 +261,8 @@ void FixBondSwap::post_integrate()
}
int tmp;
for (i = 0; i < neligible; i++) {
j = static_cast<int> (random->uniform() * neligible);
for (i = 0; i < neligible-1; i++) {
j = i + static_cast<int> (random->uniform() * (neligible-i));
tmp = alist[i];
alist[i] = alist[j];
alist[j] = tmp;
@ -256,6 +271,7 @@ void FixBondSwap::post_integrate()
// examine ntest of my eligible atoms for potential swaps
// atom I is randomly selected via atom list
// look at all J neighbors of atom I
// done in a randomized permutation, via neighbor_permutation()
// J must be on-processor (J < nlocal)
// I,J must be in fix group
// I,J must have same molecule IDs
@ -275,8 +291,10 @@ void FixBondSwap::post_integrate()
jlist = firstneigh[i];
jnum = numneigh[i];
neighbor_permutation(jnum);
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j = jlist[permute[jj]];
j &= NEIGHMASK;
if (j >= nlocal) continue;
if ((mask[j] & groupbit) == 0) continue;
@ -431,6 +449,7 @@ void FixBondSwap::post_integrate()
factor = exp(-delta/force->boltz/t_current);
if (random->uniform() < factor) accept = 1;
}
goto done;
}
}
@ -674,7 +693,7 @@ int FixBondSwap::modify_param(int narg, char **arg)
modify->delete_compute(id_temp);
tflag = 0;
}
delete [] id_temp;
delete[] id_temp;
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(id_temp);
@ -737,6 +756,34 @@ double FixBondSwap::angle_eng(int atype, int i, int j, int k)
return force->angle->single(atype,i,j,k);
}
/* ----------------------------------------------------------------------
create a random permutation of one atom's N neighbor list atoms
uses one-pass Fisher-Yates shuffle on an initial identity permutation
output: randomized permute[] vector, used to index neighbors
------------------------------------------------------------------------- */
void FixBondSwap::neighbor_permutation(int n)
{
int i,j,tmp;
if (n > maxpermute) {
delete[] permute;
maxpermute = n + DELTA_PERMUTE;
permute = new int[maxpermute];
}
// Fisher-Yates shuffle
for (i = 0; i < n; i++) permute[i] = i;
for (i = 0; i < n-1; i++) {
j = i + static_cast<int> (random->uniform() * (n-i));
tmp = permute[i];
permute[i] = permute[j];
permute[j] = tmp;
}
}
/* ----------------------------------------------------------------------
return bond swapping stats
n = 1 is # of swaps

View File

@ -46,6 +46,9 @@ class FixBondSwap : public Fix {
int *type;
double **x;
int maxpermute;
int *permute;
class NeighList *list;
class Compute *temperature;
class RanMars *random;
@ -54,6 +57,8 @@ class FixBondSwap : public Fix {
double pair_eng(int, int);
double bond_eng(int, int, int);
double angle_eng(int, int, int, int);
void neighbor_permutation(int);
};
} // namespace LAMMPS_NS

View File

@ -35,7 +35,7 @@ class FixGCMC : public Fix {
double memory_usage() override;
void write_restart(FILE *) override;
void restart(char *) override;
void *extract(const char *, int &);
void *extract(const char *, int &) override;
private:
int molecule_group, molecule_group_bit;

View File

@ -52,6 +52,9 @@ Bond::Bond(LAMMPS *_lmp) : Pointers(_lmp)
born_matrix_enable = 0;
partial_flag = 0;
single_extra = 0;
svector = nullptr;
maxeatom = maxvatom = 0;
eatom = nullptr;
vatom = nullptr;
@ -247,6 +250,91 @@ void Bond::ev_tally(int i, int j, int nlocal, int newton_bond, double ebond, dou
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global or per-atom accumulators
for virial, have delx,dely,delz and fx,fy,fz
------------------------------------------------------------------------- */
void Bond::ev_tally_xyz(int i, int j, int nlocal, int newton_bond,
double ebond, double fx, double fy, double fz,
double delx, double dely, double delz)
{
double ebondhalf,v[6];
if (eflag_either) {
if (eflag_global) {
if (newton_bond) {
energy += ebond;
} else {
ebondhalf = 0.5 * ebond;
if (i < nlocal) energy += ebondhalf;
if (j < nlocal) energy += ebondhalf;
}
}
if (eflag_atom) {
ebondhalf = 0.5 * ebond;
if (newton_bond || i < nlocal) eatom[i] += ebondhalf;
if (newton_bond || j < nlocal) eatom[j] += ebondhalf;
}
}
if (vflag_either) {
v[0] = delx * fx;
v[1] = dely * fy;
v[2] = delz * fz;
v[3] = delx * fy;
v[4] = delx * fz;
v[5] = dely * fz;
if (vflag_global) {
if (newton_bond) {
virial[0] += v[0];
virial[1] += v[1];
virial[2] += v[2];
virial[3] += v[3];
virial[4] += v[4];
virial[5] += v[5];
} else {
if (i < nlocal) {
virial[0] += 0.5 * v[0];
virial[1] += 0.5 * v[1];
virial[2] += 0.5 * v[2];
virial[3] += 0.5 * v[3];
virial[4] += 0.5 * v[4];
virial[5] += 0.5 * v[5];
}
if (j < nlocal) {
virial[0] += 0.5 * v[0];
virial[1] += 0.5 * v[1];
virial[2] += 0.5 * v[2];
virial[3] += 0.5 * v[3];
virial[4] += 0.5 * v[4];
virial[5] += 0.5 * v[5];
}
}
}
if (vflag_atom) {
if (newton_bond || i < nlocal) {
vatom[i][0] += 0.5 * v[0];
vatom[i][1] += 0.5 * v[1];
vatom[i][2] += 0.5 * v[2];
vatom[i][3] += 0.5 * v[3];
vatom[i][4] += 0.5 * v[4];
vatom[i][5] += 0.5 * v[5];
}
if (newton_bond || j < nlocal) {
vatom[j][0] += 0.5 * v[0];
vatom[j][1] += 0.5 * v[1];
vatom[j][2] += 0.5 * v[2];
vatom[j][3] += 0.5 * v[3];
vatom[j][4] += 0.5 * v[4];
vatom[j][5] += 0.5 * v[5];
}
}
}
}
/* ----------------------------------------------------------------------
write a table of bond potential energy/force vs distance to a file
------------------------------------------------------------------------- */

View File

@ -42,6 +42,9 @@ class Bond : protected Pointers {
int reinitflag; // 0 if not compatible with fix adapt
// extract() method may still need to be added
int single_extra; // number of extra single values calculated
double *svector; // vector of extra single quantities
// KOKKOS host/device flag and data masks
ExecutionSpace execution_space;
@ -99,6 +102,7 @@ class Bond : protected Pointers {
}
void ev_setup(int, int, int alloc = 1);
void ev_tally(int, int, int, int, double, double, double, double, double);
void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);
};
} // namespace LAMMPS_NS

View File

@ -50,6 +50,8 @@ BondHybrid::~BondHybrid()
delete[] keywords;
}
delete[] svector;
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
@ -238,6 +240,49 @@ void BondHybrid::settings(int narg, char **arg)
i = jarg;
nstyles++;
}
// set bond flags from sub-style flags
flags();
}
/* ----------------------------------------------------------------------
set top-level bond flags from sub-style flags
------------------------------------------------------------------------- */
void BondHybrid::flags()
{
int m;
// set comm_forward, comm_reverse, comm_reverse_off to max of any sub-style
for (m = 0; m < nstyles; m++) {
if (styles[m]) comm_forward = MAX(comm_forward,styles[m]->comm_forward);
if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse);
if (styles[m]) comm_reverse_off = MAX(comm_reverse_off,
styles[m]->comm_reverse_off);
}
init_svector();
}
/* ----------------------------------------------------------------------
initialize Bond::svector array
------------------------------------------------------------------------- */
void BondHybrid::init_svector()
{
// single_extra = list all sub-style single_extra
// allocate svector
single_extra = 0;
for (int m = 0; m < nstyles; m++)
single_extra = MAX(single_extra,styles[m]->single_extra);
if (single_extra) {
delete[] svector;
svector = new double[single_extra];
}
}
/* ----------------------------------------------------------------------
@ -359,9 +404,27 @@ double BondHybrid::single(int type, double rsq, int i, int j, double &fforce)
{
if (map[type] < 0) error->one(FLERR, "Invoked bond single on bond style none");
if (single_extra) copy_svector(type);
return styles[map[type]]->single(type, rsq, i, j, fforce);
}
/* ----------------------------------------------------------------------
copy Bond::svector data
------------------------------------------------------------------------- */
void BondHybrid::copy_svector(int type)
{
memset(svector,0,single_extra*sizeof(double));
// there is only one style in bond style hybrid for a bond type
Bond *this_style = styles[map[type]];
for (int l = 0; this_style->single_extra; ++l) {
svector[l] = this_style->svector[l];
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */

View File

@ -52,6 +52,10 @@ class BondHybrid : public Bond {
int ***bondlist; // bondlist for each sub-style
void allocate();
void flags();
virtual void init_svector();
virtual void copy_svector(int);
};
} // namespace LAMMPS_NS

View File

@ -34,7 +34,7 @@ using namespace LAMMPS_NS;
#define DELTA 10000
#define EPSILON 1.0e-12
enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE};
enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE,BN};
/* ---------------------------------------------------------------------- */
@ -54,6 +54,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = narg - 3;
bstyle = new int[nvalues];
bindex = new int[nvalues];
vstr = new char*[nvalues];
vvar = new int[nvalues];
@ -80,6 +81,11 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
bstyle[nvalues++] = VARIABLE;
vstr[nvar] = utils::strdup(&arg[iarg][2]);
nvar++;
} else if (arg[iarg][0] == 'b') {
int n = atoi(&arg[iarg][1]);
if (n <= 0) error->all(FLERR, "Invalid keyword {} in compute bond/local command", arg[iarg]);
bstyle[nvalues] = BN;
bindex[nvalues++] = n - 1;
} else break;
}
@ -131,7 +137,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
velflag = 0;
for (int i = 0; i < nvalues; i++) {
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX ||
bstyle[i] == FY || bstyle[i] == FZ) singleflag = 1;
bstyle[i] == FY || bstyle[i] == FZ || bstyle[i] == BN) singleflag = 1;
if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS ||
bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1;
}
@ -151,6 +157,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeBondLocal::~ComputeBondLocal()
{
delete [] bstyle;
delete [] bindex;
for (int i = 0; i < nvar; i++) delete [] vstr[i];
delete [] vstr;
delete [] vvar;
@ -168,6 +175,10 @@ void ComputeBondLocal::init()
if (force->bond == nullptr)
error->all(FLERR,"No bond style is defined for compute bond/local");
for (int i = 0; i < nvalues; i++)
if (bstyle[i] == BN && bindex[i] >= force->bond->single_extra)
error->all(FLERR, "Bond style does not have extra field requested by compute bond/local");
if (nvar) {
for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]);
@ -438,6 +449,9 @@ int ComputeBondLocal::compute_bonds(int flag)
ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
case BN:
ptr[n] = bond->svector[bindex[n]];
break;
}
}
}

View File

@ -39,7 +39,7 @@ class ComputeBondLocal : public Compute {
int singleflag, velflag, ghostvelflag, initflag;
int dvar;
int *bstyle, *vvar;
int *bstyle, *bindex, *vvar;
char *dstr;
char **vstr;

View File

@ -20,11 +20,13 @@
#include "error.h"
#include "fix.h"
#include "group.h"
#include "input.h"
#include "irregular.h"
#include "memory.h"
#include "modify.h"
#include "output.h"
#include "update.h"
#include "variable.h"
#include <cstring>
@ -87,6 +89,9 @@ Dump::Dump(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp)
delay_flag = 0;
write_header_flag = 1;
skipflag = 0;
skipvar = nullptr;
maxfiles = -1;
numfiles = 0;
fileidx = 0;
@ -164,6 +169,7 @@ Dump::~Dump()
delete[] format_bigint_user;
delete[] refresh;
delete[] skipvar;
// format_column_user is deallocated by child classes that use it
@ -298,6 +304,15 @@ void Dump::init()
else error->all(FLERR,"Dump could not find refresh compute ID");
}
// if skipflag, check skip variable
if (skipflag) {
skipindex = input->variable->find(skipvar);
if (skipindex < 0) error->all(FLERR,"Dump skip variable not found");
if (!input->variable->equalstyle(skipindex))
error->all(FLERR,"Variable for dump skip is invalid style");
}
// preallocation for PBC copies if requested
if (pbcflag && atom->nlocal > maxpbc) pbc_allocate();
@ -325,15 +340,6 @@ void Dump::write()
imageint *imagehold;
double **xhold,**vhold;
// if timestep < delaystep, just return
if (delay_flag && update->ntimestep < delaystep) return;
// if file per timestep, open new file
if (multifile) openfile();
if (fp) clearerr(fp);
// simulation box bounds
if (domain->triclinic == 0) {
@ -359,6 +365,24 @@ void Dump::write()
nme = count();
// if timestep < delaystep, just return
// if skip condition is defined and met, just return
// must do both these tests after count() b/c it invokes computes,
// this enables caller to trigger future invocation of needed computes
if (delay_flag && update->ntimestep < delaystep) return;
if (skipflag) {
double value = input->variable->compute_equal(skipindex);
if (value != 0.0) return;
}
// if file per timestep, open new file
// do this after skip check, so no file is opened if skip occurs
if (multifile) openfile();
if (fp) clearerr(fp);
// ntotal = total # of dump lines in snapshot
// nmax = max # of dump lines on any proc
@ -1264,6 +1288,15 @@ void Dump::modify_params(int narg, char **arg)
pbcflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"skip") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
skipflag = 1;
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
delete[] skipvar;
skipvar = utils::strdup(&arg[iarg+1][2]);
} else error->all(FLERR,"Illegal dump_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "dump_modify sort", error);
if (strcmp(arg[iarg+1],"off") == 0) sort_flag = 0;

View File

@ -97,6 +97,10 @@ class Dump : protected Pointers {
char *refresh; // compute ID to invoke refresh() on
int irefresh; // index of compute
int skipflag; // 1 if skip condition defined
char *skipvar; // name of variable to check for skip condition
int skipindex; // index of skip variable
char boundstr[9]; // encoding of boundary flags
char *format; // format string for the file write

View File

@ -82,10 +82,12 @@ int FixBondHistory::setmask()
void FixBondHistory::post_constructor()
{
// Store saved bond quantities for each atom using fix property atom
// Don't copy history to data files because this fix is typically
// not yet instantiated - history is only preserved across restarts
id_fix = utils::strdup(id + std::string("_FIX_PROP_ATOM"));
id_array = utils::strdup(std::string("d2_") + id);
modify->add_fix(fmt::format("{} {} property/atom {} {}", id_fix, group->names[igroup], id_array,
modify->add_fix(fmt::format("{} {} property/atom {} {} writedata no", id_fix, group->names[igroup], id_array,
nbond * ndata));
int tmp1, tmp2;
index = atom->find_custom(&id_array[3], tmp1, tmp2);

343
src/fix_pair.cpp Normal file
View File

@ -0,0 +1,343 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "fix_pair.h"
#include "atom.h"
#include "dump.h"
#include "error.h"
#include "force.h"
#include "fix.h"
#include "input.h"
#include "memory.h"
#include "pair.h"
#include "output.h"
#include "variable.h"
#include "update.h"
#include "variable.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) utils::missing_cmd_args(FLERR, "fix pair", error);
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery < 1) error->all(FLERR,"Illegal fix pair every value: {}", nevery);
pairname = utils::strdup(arg[4]);
pstyle = force->pair_match(pairname,1,0);
if (pstyle == nullptr) error->all(FLERR,"Pair style {} for fix pair not found", pairname);
nfield = (narg-5) / 2;
fieldname = new char*[nfield];
trigger = new int[nfield];
nfield = 0;
int iarg = 5;
while (iarg < narg) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, fmt::format("fix pair {}", arg[iarg]), error);
fieldname[nfield] = utils::strdup(arg[iarg]);
int flag = utils::inumeric(FLERR,arg[iarg+1],true,lmp);
if (flag == 0) trigger[nfield] = 0;
else if (flag == 1) trigger[nfield] = 1;
else error->all(FLERR,"Illegal fix pair {} command flag: {}", arg[iarg], arg[iarg+1]);
nfield++;
iarg += 2;
}
// set trigger names = fieldname + "_flag"
triggername = new char*[nfield];
for (int ifield = 0; ifield < nfield; ifield++) {
if (trigger[ifield] == 0) triggername[ifield] = nullptr;
else triggername[ifield] = utils::strdup(fmt::format("{}_flag", fieldname[ifield]));
}
// extract all fields just to get number of per-atom values
// returned data ptr may be NULL, if pair style has not allocated field yet
// check for recognized field cannot be done until post_force()
// also check if triggername can be extracted as a scalar value
triggerptr = new int*[nfield];
ncols = 0;
for (int ifield = 0; ifield < nfield; ifield++) {
int columns = 0; // set in case fieldname not recognized by pstyle
pstyle->extract_peratom(fieldname[ifield],columns);
if (columns) ncols += columns;
else ncols++;
if (trigger[ifield]) {
int dim;
triggerptr[ifield] = (int *) pstyle->extract(triggername[ifield],dim);
if (!triggerptr[ifield])
error->all(FLERR,"Fix pair pair style cannot extract {}", triggername[ifield]);
if (dim)
error->all(FLERR,"Fix pair pair style {} trigger {} is not a scalar",
pairname, triggername[ifield]);
}
}
// if set peratom_freq = Nevery, then cannot access the per-atom
// values as part of thermo output during minimiziation
// at different frequency or on last step of minimization
// instead set peratom_freq = 1
// ok, since vector/array always have values
// but requires the vector/array be persisted between Nevery steps
// since it may be accessed
peratom_flag = 1;
if (ncols == 1) size_peratom_cols = 0;
else size_peratom_cols = ncols;
peratom_freq = 1;
// perform initial allocation of atom-based array
// register with Atom class
vector = nullptr;
array = nullptr;
grow_arrays(atom->nmax);
atom->add_callback(Atom::GROW);
// zero the vector/array since dump may access it on timestep 0
// zero the vector/array since a variable may access it before first run
int nlocal = atom->nlocal;
if (ncols == 1) {
for (int i = 0; i < nlocal; i++)
vector[i] = 0.0;
} else {
for (int i = 0; i < nlocal; i++)
for (int m = 0; m < ncols; m++)
array[i][m] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
FixPair::~FixPair()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,Atom::GROW);
delete[] pairname;
for (int ifield = 0; ifield < nfield; ifield++) {
delete[] fieldname[ifield];
delete[] triggername[ifield];
}
delete[] fieldname;
delete[] trigger;
delete[] triggername;
delete[] triggerptr;
if (ncols == 1) memory->destroy(vector);
else memory->destroy(array);
}
/* ---------------------------------------------------------------------- */
int FixPair::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
mask |= MIN_PRE_FORCE;
mask |= POST_FORCE;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixPair::init()
{
// insure pair style still exists
pstyle = force->pair_match(pairname,1,0);
if (pstyle == nullptr) error->all(FLERR,"Pair style {} for fix pair not found", pairname);
}
/* ---------------------------------------------------------------------- */
void FixPair::setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixPair::setup_pre_force(int vflag)
{
pre_force(vflag);
}
/* ----------------------------------------------------------------------
trigger pair style computation on steps which are multiples of Nevery
------------------------------------------------------------------------- */
void FixPair::pre_force(int /*vflag*/)
{
if (update->ntimestep % nevery) return;
// set pair style triggers
for (int ifield = 0; ifield < nfield; ifield++)
if (trigger[ifield]) *(triggerptr[ifield]) = 1;
}
/* ---------------------------------------------------------------------- */
void FixPair::min_pre_force(int vflag)
{
pre_force(vflag);
}
/* ----------------------------------------------------------------------
extract results from pair style
------------------------------------------------------------------------- */
void FixPair::post_force(int /*vflag*/)
{
if (update->ntimestep % nevery) return;
// extract pair style fields one by one
// store their values in this fix
int nlocal = atom->nlocal;
int icol = 0;
int columns;
for (int ifield = 0; ifield < nfield; ifield++) {
void *pvoid = pstyle->extract_peratom(fieldname[ifield],columns);
if (pvoid == nullptr)
error->all(FLERR,"Fix pair pair style cannot extract {}",fieldname[ifield]);
if (columns == 0) {
double *pvector = (double *) pvoid;
if (ncols == 1) {
for (int i = 0; i < nlocal; i++)
vector[i] = pvector[i];
} else {
for (int i = 0; i < nlocal; i++)
array[i][icol] = pvector[i];
}
icol++;
} else {
double **parray = (double **) pvoid;
for (int i = 0; i < nlocal; i++)
for (int m = 0; m < columns; m++) {
array[i][icol] = parray[i][m];
icol++;
}
}
}
// unset pair style triggers
for (int ifield = 0; ifield < nfield; ifield++)
if (trigger[ifield]) *(triggerptr[ifield]) = 0;
}
/* ---------------------------------------------------------------------- */
void FixPair::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
allocate atom-based vector or array
------------------------------------------------------------------------- */
void FixPair::grow_arrays(int nmax)
{
if (ncols == 1) {
memory->grow(vector,nmax,"store/state:vector");
vector_atom = vector;
} else {
memory->grow(array,nmax,ncols,"store/state:array");
array_atom = array;
}
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixPair::copy_arrays(int i, int j, int /*delflag*/)
{
if (ncols == 1) {
vector[j] = vector[i];
} else {
for (int m = 0; m < ncols; m++)
array[j][m] = array[i][m];
}
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixPair::pack_exchange(int i, double *buf)
{
if (ncols == 1) {
buf[0] = vector[i];
} else {
for (int m = 0; m < ncols; m++)
buf[m] = array[i][m];
}
return ncols;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixPair::unpack_exchange(int nlocal, double *buf)
{
if (ncols == 1) {
vector[nlocal] = buf[0];
} else {
for (int m = 0; m < ncols; m++)
array[nlocal][m] = buf[m];
}
return ncols;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based vector or array
------------------------------------------------------------------------- */
double FixPair::memory_usage()
{
double bytes = 0.0;
if (ncols == 1) bytes += (double)atom->nmax * sizeof(double);
else bytes += (double)atom->nmax*ncols * sizeof(double);
return bytes;
}

62
src/fix_pair.h Normal file
View File

@ -0,0 +1,62 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(pair,FixPair);
// clang-format on
#else
#ifndef LMP_FIX_PAIR_H
#define LMP_FIX_PAIR_H
#include "fix.h"
namespace LAMMPS_NS {
class FixPair : public Fix {
public:
FixPair(class LAMMPS *, int, char **);
~FixPair() override;
int setmask() override;
void init() override;
void setup(int) override;
void setup_pre_force(int) override;
void pre_force(int) override;
void min_pre_force(int) override;
void post_force(int) override;
void min_post_force(int) override;
void grow_arrays(int) override;
void copy_arrays(int, int, int) override;
int pack_exchange(int, double *) override;
int unpack_exchange(int, double *) override;
double memory_usage() override;
private:
int nevery,nfield,ncols;
char *pairname;
char **fieldname,**triggername;
int *trigger;
int **triggerptr;
class Pair *pstyle;
double *vector;
double **array;
};
} // namespace LAMMPS_NS
#endif
#endif

View File

@ -35,7 +35,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (narg < 4) error->all(FLERR, "Illegal fix property/atom command");
restart_peratom = 1;
wd_section = 1;
wd_section = 1; // can be overwitten using optional arguments
int iarg = 3;
nvalue = narg - iarg;
@ -153,6 +153,10 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix property/atom command");
border = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "writedata") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix property/atom command");
wd_section = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else
error->all(FLERR, "Illegal fix property/atom command");
}

View File

@ -17,6 +17,7 @@
#include "atom_vec.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "pair.h"
@ -26,6 +27,8 @@
using namespace LAMMPS_NS;
using namespace FixConst;
#define DELTA 10000
/* ---------------------------------------------------------------------- */
FixUpdateSpecialBonds::FixUpdateSpecialBonds(LAMMPS *lmp, int narg, char **arg) :
@ -48,6 +51,11 @@ int FixUpdateSpecialBonds::setmask()
void FixUpdateSpecialBonds::setup(int /*vflag*/)
{
// error if more than one fix update/special/bonds
if (modify->get_fix_by_style("UPDATE_SPECIAL_BONDS").size() > 1)
error->all(FLERR,"More than one fix update/special/bonds");
// Require atoms know about all of their bonds and if they break
if (force->newton_bond) error->all(FLERR, "Fix update/special/bonds requires Newton bond off");
@ -61,18 +69,22 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/)
if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 ||
force->special_coul[3] != 1.0)
error->all(FLERR, "Fix update/special/bonds requires special Coulomb weights = 1,1,1");
// Implies neighbor->special_flag = [X, 2, 1, 1]
new_broken_pairs.clear();
broken_pairs.clear();
new_created_pairs.clear();
created_pairs.clear();
}
/* ----------------------------------------------------------------------
Update special bond list and atom bond arrays, empty broken bond list
Update special bond list and atom bond arrays, empty broken/created lists
------------------------------------------------------------------------- */
void FixUpdateSpecialBonds::pre_exchange()
{
int i, j, m, n1, n3;
int i, j, m, n1;
tagint tagi, tagj;
int nlocal = atom->nlocal;
@ -83,21 +95,19 @@ void FixUpdateSpecialBonds::pre_exchange()
for (auto const &it : broken_pairs) {
tagi = it.first;
tagj = it.second;
i = atom->map(tagi);
j = atom->map(tagj);
// remove i from special bond list for atom j and vice versa
// ignore n2, n3 since 1-3, 1-4 special factors required to be 1.0
if (i < nlocal) {
slist = special[i];
n1 = nspecial[i][0];
for (m = 0; m < n1; m++)
if (slist[m] == tagj) break;
n3 = nspecial[i][2];
for (; m < n3 - 1; m++) slist[m] = slist[m + 1];
for (; m < n1 - 1; m++) slist[m] = slist[m + 1];
nspecial[i][0]--;
nspecial[i][1]--;
nspecial[i][2]--;
nspecial[i][1] = nspecial[i][2] = nspecial[i][0];
}
if (j < nlocal) {
@ -105,19 +115,43 @@ void FixUpdateSpecialBonds::pre_exchange()
n1 = nspecial[j][0];
for (m = 0; m < n1; m++)
if (slist[m] == tagi) break;
n3 = nspecial[j][2];
for (; m < n3 - 1; m++) slist[m] = slist[m + 1];
for (; m < n1 - 1; m++) slist[m] = slist[m + 1];
nspecial[j][0]--;
nspecial[j][1]--;
nspecial[j][2]--;
nspecial[j][1] = nspecial[j][2] = nspecial[j][0];
}
}
for (auto const &it : created_pairs) {
tagi = it.first;
tagj = it.second;
i = atom->map(tagi);
j = atom->map(tagj);
// add i to special bond list for atom j and vice versa
// ignore n2, n3 since 1-3, 1-4 special factors required to be 1.0
n1 = nspecial[i][0];
if (n1 >= atom->maxspecial)
error->one(FLERR,"Special list size exceeded in fix update/special/bond");
special[i][n1] = tagj;
nspecial[i][0] += 1;
nspecial[i][1] = nspecial[i][2] = nspecial[i][0];
n1 = nspecial[j][0];
if (n1 >= atom->maxspecial)
error->one(FLERR,"Special list size exceeded in fix update/special/bond");
special[j][n1] = tagi;
nspecial[j][0] += 1;
nspecial[j][1] = nspecial[j][2] = nspecial[j][0];
}
broken_pairs.clear();
created_pairs.clear();
}
/* ----------------------------------------------------------------------
Loop neighbor list and update special bond lists for recently broken bonds
Update special lists for recently broken/created bonds
Assumes appropriate atom/bond arrays were updated, e.g. had called
neighbor->add_temporary_bond(i1, i2, btype);
------------------------------------------------------------------------- */
void FixUpdateSpecialBonds::pre_force(int /*vflag*/)
@ -129,7 +163,7 @@ void FixUpdateSpecialBonds::pre_force(int /*vflag*/)
int nlocal = atom->nlocal;
tagint *tag = atom->tag;
NeighList *list = force->pair->list; // may need to be generalized to work with pair hybrid*
NeighList *list = force->pair->list; // may need to be generalized for pair hybrid*
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -163,7 +197,37 @@ void FixUpdateSpecialBonds::pre_force(int /*vflag*/)
}
}
}
for (auto const &it : new_created_pairs) {
tag1 = it.first;
tag2 = it.second;
i1 = atom->map(tag1);
i2 = atom->map(tag2);
// Loop through atoms of owned atoms i j and update SB bits
if (i1 < nlocal) {
jlist = firstneigh[i1];
jnum = numneigh[i1];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (((j >> SBBITS) & 3) != 0) continue; // Skip bonded pairs
if (tag[j] == tag2) jlist[jj] = j ^ (1 << SBBITS); // Add 1-2 special bond bits
}
}
if (i2 < nlocal) {
jlist = firstneigh[i2];
jnum = numneigh[i2];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (((j >> SBBITS) & 3) != 0) continue; // Skip bonded pairs
if (tag[j] == tag1) jlist[jj] = j ^ (1 << SBBITS); // Add 1-2 special bond bits
}
}
}
new_broken_pairs.clear();
new_created_pairs.clear();
}
/* ---------------------------------------------------------------------- */
@ -174,3 +238,12 @@ void FixUpdateSpecialBonds::add_broken_bond(int i, int j)
new_broken_pairs.push_back(tag_pair);
broken_pairs.push_back(tag_pair);
}
/* ---------------------------------------------------------------------- */
void FixUpdateSpecialBonds::add_created_bond(int i, int j)
{
auto tag_pair = std::make_pair(atom->tag[i], atom->tag[j]);
new_created_pairs.push_back(tag_pair);
created_pairs.push_back(tag_pair);
}

View File

@ -35,12 +35,18 @@ class FixUpdateSpecialBonds : public Fix {
void pre_exchange() override;
void pre_force(int) override;
void add_broken_bond(int, int);
void add_created_bond(int, int);
protected:
// Create two arrays to store bonds broken this timestep (new)
// and since the last neighbor list build
std::vector<std::pair<tagint, tagint>> new_broken_pairs;
std::vector<std::pair<tagint, tagint>> broken_pairs;
// Create two arrays to store newly created this timestep (new)
// and since the last neighbor list build
std::vector<std::pair<tagint, tagint>> new_created_pairs;
std::vector<std::pair<tagint, tagint>> created_pairs;
};
} // namespace LAMMPS_NS

View File

@ -2920,6 +2920,16 @@ bigint Neighbor::get_nneigh_half()
return nneighhalf;
}
/* ----------------------------------------------------------------------
add pair of atoms to bondlist array
will only persist until the next neighbor build
------------------------------------------------------------------------- */
void Neighbor::add_temporary_bond(int i1, int i2, int btype)
{
neigh_bond->add_temporary_bond(i1, i2, btype);
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */

View File

@ -166,6 +166,7 @@ class Neighbor : protected Pointers {
bigint get_nneigh_full(); // return number of neighbors in a regular full neighbor list
bigint get_nneigh_half(); // return number of neighbors in a regular half neighbor list
void add_temporary_bond(int, int, int); // add temporary bond to bondlist array
double memory_usage();
bigint last_setup_bins; // step of last neighbor::setup_bins() call

View File

@ -24,6 +24,7 @@
using namespace LAMMPS_NS;
#define LB_FACTOR 1.5
#define DELTA 10000
/* ---------------------------------------------------------------------- */
@ -207,6 +208,21 @@ void NTopo::dihedral_check(int nlist, int **list)
/* ---------------------------------------------------------------------- */
void NTopo::add_temporary_bond(int i1, int i2, int btype)
{
if (nbondlist == maxbond) {
maxbond += DELTA;
memory->grow(bondlist, maxbond, 3, "neigh_topo:bondlist");
}
bondlist[nbondlist][0] = i1;
bondlist[nbondlist][1] = i2;
bondlist[nbondlist][2] = btype;
nbondlist++;
}
/* ---------------------------------------------------------------------- */
double NTopo::memory_usage()
{
double bytes = 0;

View File

@ -28,6 +28,7 @@ class NTopo : protected Pointers {
virtual void build() = 0;
void add_temporary_bond(int, int, int);
double memory_usage();
protected:

View File

@ -215,6 +215,7 @@ class Pair : protected Pointers {
// specific child-class methods for certain Pair styles
virtual void *extract(const char *, int &) { return nullptr; }
virtual void *extract_peratom(const char *, int &) { return nullptr; }
virtual void swap_eam(double *, double **) {}
virtual void reset_dt() {}
virtual void min_xf_pointers(int, double **, double **) {}

View File

@ -1 +1 @@
#define LAMMPS_VERSION "3 Aug 2022"
#define LAMMPS_VERSION "15 Sep 2022"