mirror of https://github.com/lammps/lammps.git
Merge branch 'develop' into fortran-tinkering
This commit is contained in:
Axel Kohlmeyer
commit
0f2a7d3f33
|
|
@ -8,8 +8,8 @@ option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an al
|
|||
|
||||
if(DOWNLOAD_MDI)
|
||||
message(STATUS "MDI download requested - we will build our own")
|
||||
set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.11.tar.gz" CACHE STRING "URL for MDI tarball")
|
||||
set(MDI_MD5 "3791fe5081405c14aac07d4687f1cc58" CACHE STRING "MD5 checksum for MDI tarball")
|
||||
set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.12.tar.gz" CACHE STRING "URL for MDI tarball")
|
||||
set(MDI_MD5 "7a222353ae8e03961d5365e6cd48baee" CACHE STRING "MD5 checksum for MDI tarball")
|
||||
mark_as_advanced(MDI_URL)
|
||||
mark_as_advanced(MDI_MD5)
|
||||
enable_language(C)
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
""""""""""""
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
""""""""""""
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 *skip* 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::
|
||||
|
||||
|
|
@ -181,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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -44,14 +44,15 @@ ellipsoidal and spherical particle via the formulas
|
|||
U_r = & 4 \epsilon ( \varrho^{12} - \varrho^6) \\
|
||||
\varrho = & \frac{\sigma}{ h_{12} + \gamma \sigma}
|
||||
|
||||
where A1 and A2 are the transformation matrices from the simulation box
|
||||
frame to the body frame and :math:`r_{12}` is the center to center
|
||||
vector between the particles. :math:`U_r` controls the shifted distance
|
||||
dependent interaction based on the distance of closest approach of the
|
||||
two particles (:math:`h_{12}`) and the user-specified shift parameter
|
||||
gamma. When both particles are spherical, the formula reduces to the
|
||||
usual Lennard-Jones interaction (see details below for when Gay-Berne
|
||||
treats a particle as "spherical").
|
||||
where :math:`\mathbf{A}_1` and :math:`\mathbf{A}_2` are the
|
||||
transformation matrices from the simulation box frame to the body frame
|
||||
and :math:`r_{12}` is the center to center vector between the particles.
|
||||
:math:`U_r` controls the shifted distance dependent interaction based on
|
||||
the distance of closest approach of the two particles (:math:`h_{12}`)
|
||||
and the user-specified shift parameter :math:`\gamma`. When both
|
||||
particles are spherical, the formula reduces to the usual Lennard-Jones
|
||||
interaction (see details below for when Gay-Berne treats a particle as
|
||||
"spherical").
|
||||
|
||||
For large uniform molecules it has been shown that the energy
|
||||
parameters are approximately representable in terms of local contact
|
||||
|
|
@ -74,8 +75,9 @@ listed below and in `this supplementary document <PDF/pair_gayberne_extra.pdf>`_
|
|||
|
||||
Use of this pair style requires the NVE, NVT, or NPT fixes with the
|
||||
*asphere* extension (e.g. :doc:`fix nve/asphere <fix_nve_asphere>`) in
|
||||
order to integrate particle rotation. Additionally, :doc:`atom_style ellipsoid <atom_style>` should be used since it defines the
|
||||
rotational state and the size and shape of each ellipsoidal particle.
|
||||
order to integrate particle rotation. Additionally, :doc:`atom_style
|
||||
ellipsoid <atom_style>` should be used since it defines the rotational
|
||||
state and the size and shape of each ellipsoidal particle.
|
||||
|
||||
The following coefficients must be defined for each pair of atoms
|
||||
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
|
||||
|
|
@ -129,7 +131,7 @@ that type. e.g. in a "pair_coeff I J" command.
|
|||
|
||||
.. note::
|
||||
|
||||
If the :math:`\epsilon` a = b = c for an atom type, and if the shape
|
||||
If the :math:`\epsilon_{a}` = :math:`\epsilon_{b}` = :math:`\epsilon_{c}` for an atom type, and if the shape
|
||||
of the particle itself is spherical, meaning its 3 shape parameters
|
||||
are all the same, then the particle is treated as an LJ sphere by the
|
||||
Gay-Berne potential. This is significant because if two LJ spheres
|
||||
|
|
@ -137,7 +139,7 @@ that type. e.g. in a "pair_coeff I J" command.
|
|||
their interaction energy/force using the specified epsilon and sigma
|
||||
as the standard LJ parameters. This is much cheaper to compute than
|
||||
the full Gay-Berne formula. To treat the particle as a LJ sphere
|
||||
with sigma = D, you should normally set :math:`\epsilon` a = b = c =
|
||||
with sigma = D, you should normally set :math:`\epsilon_{a}` = :math:`\epsilon_{b}` = :math:`\epsilon_{c}` =
|
||||
1.0, set the pair_coeff :math:`\sigma = D`, and also set the 3 shape
|
||||
parameters for the particle to D. The one exception is that if the 3
|
||||
shape parameters are set to 0.0, which is a valid way in LAMMPS to
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
|
@ -32,12 +32,13 @@ make lib-mdi args="-m mpi" # build MDI lib with same settings as in the mpi Make
|
|||
|
||||
# settings
|
||||
|
||||
version = "1.4.11"
|
||||
version = "1.4.12"
|
||||
url = "https://github.com/MolSSI-MDI/MDI_Library/archive/v%s.tar.gz" % version
|
||||
|
||||
# known checksums for different MDI versions. used to validate the download.
|
||||
checksums = { \
|
||||
'1.4.11' : '3791fe5081405c14aac07d4687f1cc58', \
|
||||
'1.4.12' : '7a222353ae8e03961d5365e6cd48baee', \
|
||||
}
|
||||
|
||||
# print error message or help
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -47,7 +47,7 @@ static int cmptagint(const void *p1, const void *p2)
|
|||
|
||||
void Group2Ndx::command(int narg, char **arg)
|
||||
{
|
||||
FILE *fp;
|
||||
FILE *fp = nullptr;
|
||||
|
||||
if (narg < 1) error->all(FLERR,"Illegal group2ndx command");
|
||||
|
||||
|
|
@ -57,8 +57,7 @@ void Group2Ndx::command(int narg, char **arg)
|
|||
if (comm->me == 0) {
|
||||
fp = fopen(arg[0], "w");
|
||||
if (fp == nullptr)
|
||||
error->one(FLERR,"Cannot open index file for writing: {}",
|
||||
utils::getsyserror());
|
||||
error->one(FLERR,"Cannot open index file for writing: {}", utils::getsyserror());
|
||||
utils::logmesg(lmp,"Writing groups to index file {}:\n",arg[0]);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -51,6 +51,7 @@ FixNVTSllodEff::FixNVTSllodEff(LAMMPS *lmp, int narg, char **arg) :
|
|||
modify->add_compute(fmt::format("{} {} tmp/deform/eff",
|
||||
id_temp,group->names[igroup]));
|
||||
tcomputeflag = 1;
|
||||
nondeformbias = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
|
|||
|
|
@ -398,7 +398,7 @@ void EwaldDisp::reallocate()
|
|||
h[0] = unit[0]*ix;
|
||||
h[1] = unit[5]*ix+unit[1]*iy;
|
||||
h[2] = unit[4]*ix+unit[3]*iy+unit[2]*iz;
|
||||
if ((*(flag++) = h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=gsqmx)) ++nkvec;
|
||||
if ((*(flag++) = (h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=gsqmx))) ++nkvec;
|
||||
}
|
||||
|
||||
if (nkvec>nkvec_max) {
|
||||
|
|
|
|||
|
|
@ -140,6 +140,8 @@ FixChargeRegulation::FixChargeRegulation(LAMMPS *lmp, int narg, char **arg) :
|
|||
nsalt_successes = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixChargeRegulation::~FixChargeRegulation() {
|
||||
|
||||
memory->destroy(ptype_ID);
|
||||
|
|
@ -153,14 +155,23 @@ FixChargeRegulation::~FixChargeRegulation() {
|
|||
int igroupall = group->find("all");
|
||||
neighbor->exclusion_group_group_delete(exclusion_group, igroupall);
|
||||
}
|
||||
|
||||
if (groupstrings) {
|
||||
for (int i = 0; i < ngroups; ++i) delete[] groupstrings[i];
|
||||
memory->destroy(groupstrings);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixChargeRegulation::setmask() {
|
||||
int mask = 0;
|
||||
mask |= PRE_EXCHANGE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::init() {
|
||||
|
||||
triclinic = domain->triclinic;
|
||||
|
|
@ -246,6 +257,8 @@ void FixChargeRegulation::init() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::pre_exchange() {
|
||||
|
||||
if (next_reneighbor != update->ntimestep) return;
|
||||
|
|
@ -372,16 +385,15 @@ void FixChargeRegulation::pre_exchange() {
|
|||
next_reneighbor = update->ntimestep + nevery;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::forward_acid() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
double dummyp[3];
|
||||
double pos[3];
|
||||
pos[0] = 0;
|
||||
pos[1] = 0;
|
||||
pos[2] = 0; // acid/base particle position
|
||||
double pos_all[3];
|
||||
double dummyp[3] = {0.0, 0.0, 0.0};
|
||||
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
|
||||
double pos_all[3] = {0.0, 0.0, 0.0};
|
||||
int m1 = -1, m2 = -1;
|
||||
|
||||
m1 = get_random_particle(acid_type, 0, 0, dummyp);
|
||||
|
|
@ -432,17 +444,16 @@ void FixChargeRegulation::forward_acid() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::backward_acid() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
int mask_tmp;
|
||||
double dummyp[3];
|
||||
double pos[3];
|
||||
pos[0] = 0;
|
||||
pos[1] = 0;
|
||||
pos[2] = 0; // acid/base particle position
|
||||
double pos_all[3];
|
||||
double dummyp[3] = {0.0, 0.0, 0.0};
|
||||
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
|
||||
double pos_all[3] = {0.0, 0.0, 0.0};
|
||||
int m1 = -1, m2 = -1;
|
||||
|
||||
m1 = get_random_particle(acid_type, -1, 0, dummyp);
|
||||
|
|
@ -510,16 +521,15 @@ void FixChargeRegulation::backward_acid() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::forward_base() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
double dummyp[3];
|
||||
double pos[3];
|
||||
pos[0] = 0;
|
||||
pos[1] = 0;
|
||||
pos[2] = 0; // acid/base particle position
|
||||
double pos_all[3];
|
||||
double dummyp[3] = {0.0, 0.0, 0.0};
|
||||
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
|
||||
double pos_all[3] = {0.0, 0.0, 0.0};
|
||||
int m1 = -1, m2 = -1;
|
||||
|
||||
m1 = get_random_particle(base_type, 0, 0, dummyp);
|
||||
|
|
@ -570,17 +580,16 @@ void FixChargeRegulation::forward_base() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::backward_base() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
double dummyp[3];
|
||||
double dummyp[3] = {0.0, 0.0, 0.0};
|
||||
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
|
||||
double pos_all[3] = {0.0, 0.0, 0.0};
|
||||
int mask_tmp;
|
||||
double pos[3];
|
||||
pos[0] = 0;
|
||||
pos[1] = 0;
|
||||
pos[2] = 0; // acid/base particle position
|
||||
double pos_all[3];
|
||||
int m1 = -1, m2 = -1;
|
||||
|
||||
m1 = get_random_particle(base_type, 1, 0, dummyp);
|
||||
|
|
@ -647,11 +656,13 @@ void FixChargeRegulation::backward_base() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::forward_ions() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
double dummyp[3];
|
||||
double dummyp[3] = {0.0, 0.0, 0.0};
|
||||
int m1 = -1, m2 = -1;
|
||||
factor = volume_rx * volume_rx * c10pI_plus * c10pI_minus /
|
||||
((1 + ncation) * (1 + nanion));
|
||||
|
|
@ -682,13 +693,14 @@ void FixChargeRegulation::forward_ions() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::backward_ions() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor;
|
||||
int mask1_tmp = 0, mask2_tmp = 0;
|
||||
double dummyp[3];
|
||||
double dummyp[3] = {0.0, 0.0, 0.0};
|
||||
int m1 = -1, m2 = -1;
|
||||
|
||||
m1 = get_random_particle(cation_type, +1, 0, dummyp);
|
||||
|
|
@ -764,11 +776,13 @@ void FixChargeRegulation::backward_ions() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::forward_ions_multival() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor = 1;
|
||||
double dummyp[3];
|
||||
double dummyp[3] = {0.0, 0.0, 0.0};
|
||||
|
||||
// particle ID array for all ions to be inserted
|
||||
auto mm = std::unique_ptr<int[]>(new int[salt_charge_ratio + 1]);
|
||||
|
|
@ -822,11 +836,13 @@ void FixChargeRegulation::forward_ions_multival() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::backward_ions_multival() {
|
||||
|
||||
double energy_before = energy_stored;
|
||||
double factor = 1;
|
||||
double dummyp[3]; // dummy particle
|
||||
double dummyp[3] = {0.0, 0.0, 0.0}; // dummy particle
|
||||
// particle ID array for all deleted ions
|
||||
auto mm = std::unique_ptr<int[]>(new int[salt_charge_ratio + 1]);
|
||||
// charge array for all deleted ions
|
||||
|
|
@ -943,6 +959,8 @@ void FixChargeRegulation::backward_ions_multival() {
|
|||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixChargeRegulation::insert_particle(int ptype, double charge, double rd, double *target) {
|
||||
|
||||
// insert a particle of type (ptype) with charge (charge) within distance (rd) of (target)
|
||||
|
|
@ -1013,6 +1031,8 @@ int FixChargeRegulation::insert_particle(int ptype, double charge, double rd, do
|
|||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixChargeRegulation::get_random_particle(int ptype, double charge, double rd, double *target) {
|
||||
|
||||
// returns a randomly chosen particle of type (ptype) with charge (charge)
|
||||
|
|
@ -1077,6 +1097,8 @@ int FixChargeRegulation::get_random_particle(int ptype, double charge, double rd
|
|||
return -1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixChargeRegulation::energy_full() {
|
||||
if (triclinic) domain->x2lamda(atom->nlocal);
|
||||
domain->pbc();
|
||||
|
|
@ -1136,6 +1158,8 @@ double FixChargeRegulation::energy_full() {
|
|||
return total_energy;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixChargeRegulation::particle_number_xrd(int ptype, double charge, double rd, double *target) {
|
||||
|
||||
int count = 0;
|
||||
|
|
@ -1165,6 +1189,8 @@ int FixChargeRegulation::particle_number_xrd(int ptype, double charge, double rd
|
|||
return count_sum;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixChargeRegulation::particle_number(int ptype, double charge) {
|
||||
|
||||
int count = 0;
|
||||
|
|
@ -1177,6 +1203,8 @@ int FixChargeRegulation::particle_number(int ptype, double charge) {
|
|||
return count_sum;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixChargeRegulation::compute_vector(int n) {
|
||||
if (n == 0) {
|
||||
return nacid_attempts + nbase_attempts + nsalt_attempts;
|
||||
|
|
@ -1253,6 +1281,8 @@ void FixChargeRegulation::restart(char *buf)
|
|||
error->all(FLERR,"Must not reset timestep when restarting fix gcmc");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::setThermoTemperaturePointer() {
|
||||
int ifix = -1;
|
||||
ifix = modify->find_fix(idftemp);
|
||||
|
|
@ -1266,6 +1296,8 @@ void FixChargeRegulation::setThermoTemperaturePointer() {
|
|||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixChargeRegulation::assign_tags() {
|
||||
// Assign tags to ions with zero tags
|
||||
if (atom->tag_enable) {
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -48,20 +48,18 @@
|
|||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using namespace MathConst;
|
||||
using MathConst::MY_2PI;
|
||||
|
||||
#define MAXENERGYTEST 1.0e50
|
||||
|
||||
enum{EXCHATOM,EXCHMOL}; // exchmode
|
||||
enum { EXCHATOM, EXCHMOL }; // exchmode
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
region(nullptr), idregion(nullptr), full_flag(false), local_gas_list(nullptr),
|
||||
molcoords(nullptr),molq(nullptr), molimage(nullptr), random_equal(nullptr)
|
||||
Fix(lmp, narg, arg), region(nullptr), idregion(nullptr), full_flag(false), molcoords(nullptr),
|
||||
molq(nullptr), molimage(nullptr), random_equal(nullptr), c_pe(nullptr)
|
||||
{
|
||||
if (narg < 8) error->all(FLERR,"Illegal fix widom command");
|
||||
if (narg < 8) utils::missing_cmd_args(FLERR, "fix widom", error);
|
||||
|
||||
if (atom->molecular == Atom::TEMPLATE)
|
||||
error->all(FLERR,"Fix widom does not (yet) work with atom_style template");
|
||||
|
|
@ -83,11 +81,11 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
|
|||
seed = utils::inumeric(FLERR,arg[6],false,lmp);
|
||||
insertion_temperature = utils::numeric(FLERR,arg[7],false,lmp);
|
||||
|
||||
if (nevery <= 0) error->all(FLERR,"Illegal fix widom command");
|
||||
if (ninsertions < 0) error->all(FLERR,"Illegal fix widom command");
|
||||
if (seed <= 0) error->all(FLERR,"Illegal fix widom command");
|
||||
if (nevery <= 0) error->all(FLERR,"Invalid fix widom every argument: {}", nevery);
|
||||
if (ninsertions < 0) error->all(FLERR,"Invalid fix widom insertions argument: {}", ninsertions);
|
||||
if (seed <= 0) error->all(FLERR,"Invalid fix widom seed argument: {}", seed);
|
||||
if (insertion_temperature < 0.0)
|
||||
error->all(FLERR,"Illegal fix widom command");
|
||||
error->all(FLERR,"Invalid fix widom temperature argument: {}", insertion_temperature);
|
||||
|
||||
// read options from end of input line
|
||||
|
||||
|
|
@ -99,13 +97,12 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
// error checks on region and its extent being inside simulation box
|
||||
|
||||
region_xlo = region_xhi = region_ylo = region_yhi =
|
||||
region_zlo = region_zhi = 0.0;
|
||||
region_xlo = region_xhi = region_ylo = region_yhi = region_zlo = region_zhi = 0.0;
|
||||
if (region) {
|
||||
if (region->bboxflag == 0)
|
||||
error->all(FLERR,"Fix widom region does not support a bounding box");
|
||||
error->all(FLERR,"Fix widom region {} does not support a bounding box", region->id);
|
||||
if (region->dynamic_check())
|
||||
error->all(FLERR,"Fix widom region cannot be dynamic");
|
||||
error->all(FLERR,"Fix widom region {} cannot be dynamic", region->id);
|
||||
|
||||
region_xlo = region->extent_xlo;
|
||||
region_xhi = region->extent_xhi;
|
||||
|
|
@ -117,7 +114,7 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (region_xlo < domain->boxlo[0] || region_xhi > domain->boxhi[0] ||
|
||||
region_ylo < domain->boxlo[1] || region_yhi > domain->boxhi[1] ||
|
||||
region_zlo < domain->boxlo[2] || region_zhi > domain->boxhi[2])
|
||||
error->all(FLERR,"Fix widom region extends outside simulation box");
|
||||
error->all(FLERR,"Fix widom region {} extends outside simulation box", region->id);
|
||||
|
||||
// estimate region volume using MC trials
|
||||
|
||||
|
|
@ -132,8 +129,8 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
|
|||
inside++;
|
||||
}
|
||||
|
||||
double max_region_volume = (region_xhi - region_xlo) *
|
||||
(region_yhi - region_ylo) * (region_zhi - region_zlo);
|
||||
double max_region_volume = (region_xhi - region_xlo) * (region_yhi - region_ylo)
|
||||
* (region_zhi - region_zlo);
|
||||
|
||||
region_volume = max_region_volume * static_cast<double>(inside) / static_cast<double>(attempts);
|
||||
}
|
||||
|
|
@ -141,28 +138,25 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
|
|||
// error check and further setup for exchmode = EXCHMOL
|
||||
|
||||
if (exchmode == EXCHMOL) {
|
||||
if (onemols[imol]->xflag == 0)
|
||||
error->all(FLERR,"Fix widom molecule must have coordinates");
|
||||
if (onemols[imol]->typeflag == 0)
|
||||
error->all(FLERR,"Fix widom molecule must have atom types");
|
||||
if (onemol->xflag == 0)
|
||||
error->all(FLERR,"Fix widom molecule {} must have coordinates", onemol->id);
|
||||
if (onemol->typeflag == 0)
|
||||
error->all(FLERR,"Fix widom molecule {} must have atom types", onemol->id);
|
||||
if (nwidom_type != 0)
|
||||
error->all(FLERR,"Atom type must be zero in fix widom mol command");
|
||||
if (onemols[imol]->qflag == 1 && atom->q == nullptr)
|
||||
error->all(FLERR,"Fix widom molecule has charges, but atom style does not");
|
||||
if (onemol->qflag == 1 && atom->q == nullptr)
|
||||
error->all(FLERR,"Fix widom molecule {} has charges, but atom style does not", onemol->id);
|
||||
|
||||
if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
|
||||
error->all(FLERR,"Fix widom molecule template ID must be same "
|
||||
"as atom_style template ID");
|
||||
onemols[imol]->check_attributes();
|
||||
onemol->check_attributes();
|
||||
}
|
||||
|
||||
if (charge_flag && atom->q == nullptr)
|
||||
error->all(FLERR,"Fix Widom atom has charge, but atom style does not");
|
||||
error->all(FLERR,"Fix widom atom has charge, but atom style does not");
|
||||
|
||||
// setup of array of coordinates for molecule insertion
|
||||
|
||||
if (exchmode == EXCHATOM) natoms_per_molecule = 1;
|
||||
else natoms_per_molecule = onemols[imol]->natoms;
|
||||
else natoms_per_molecule = onemol->natoms;
|
||||
nmaxmolatoms = natoms_per_molecule;
|
||||
grow_molecule_arrays(nmaxmolatoms);
|
||||
|
||||
|
|
@ -173,7 +167,7 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
// zero out counters
|
||||
widom_nmax = 0;
|
||||
local_gas_list = nullptr;
|
||||
ave_widom_chemical_potential = 0.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -201,17 +195,16 @@ void FixWidom::options(int narg, char **arg)
|
|||
|
||||
int iarg = 0;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"mol") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix widom command");
|
||||
imol = atom->find_molecule(arg[iarg+1]);
|
||||
if (imol == -1)
|
||||
error->all(FLERR,"Molecule template ID for fix widom does not exist");
|
||||
if (atom->molecules[imol]->nset > 1 && comm->me == 0)
|
||||
error->warning(FLERR,"Molecule template for "
|
||||
"fix widom has multiple molecules");
|
||||
if (strcmp(arg[iarg],"mol") == 0) {
|
||||
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix widom mol", error);
|
||||
auto onemols = atom->get_molecule_by_id(arg[iarg+1]);
|
||||
if (onemols.size() == 0)
|
||||
error->all(FLERR,"Molecule template ID {} for fix widom does not exist", arg[iarg+1]);
|
||||
if (onemols.size() > 1 && comm->me == 0)
|
||||
error->warning(FLERR,"Molecule template {} for fix widom has multiple molecules; "
|
||||
"will use only the first molecule", arg[iarg+1]);
|
||||
exchmode = EXCHMOL;
|
||||
onemols = atom->molecules;
|
||||
nmol = onemols[imol]->nset;
|
||||
onemol = onemols[0];
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"region") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix widom command");
|
||||
|
|
@ -243,7 +236,6 @@ FixWidom::~FixWidom()
|
|||
delete[] idregion;
|
||||
delete random_equal;
|
||||
|
||||
memory->destroy(local_gas_list);
|
||||
memory->destroy(molcoords);
|
||||
memory->destroy(molq);
|
||||
memory->destroy(molimage);
|
||||
|
|
@ -292,7 +284,7 @@ void FixWidom::init()
|
|||
}
|
||||
}
|
||||
|
||||
if (full_flag) c_pe = modify->compute[modify->find_compute("thermo_pe")];
|
||||
if (full_flag) c_pe = modify->get_compute_by_id("thermo_pe");
|
||||
|
||||
if (exchmode == EXCHATOM) {
|
||||
if (nwidom_type <= 0 || nwidom_type > atom->ntypes)
|
||||
|
|
@ -360,22 +352,21 @@ void FixWidom::init()
|
|||
|
||||
if (exchmode == EXCHMOL) {
|
||||
|
||||
onemols[imol]->compute_mass();
|
||||
onemols[imol]->compute_com();
|
||||
gas_mass = onemols[imol]->masstotal;
|
||||
for (int i = 0; i < onemols[imol]->natoms; i++) {
|
||||
onemols[imol]->x[i][0] -= onemols[imol]->com[0];
|
||||
onemols[imol]->x[i][1] -= onemols[imol]->com[1];
|
||||
onemols[imol]->x[i][2] -= onemols[imol]->com[2];
|
||||
onemol->compute_mass();
|
||||
onemol->compute_com();
|
||||
gas_mass = onemol->masstotal;
|
||||
for (int i = 0; i < onemol->natoms; i++) {
|
||||
onemol->x[i][0] -= onemol->com[0];
|
||||
onemol->x[i][1] -= onemol->com[1];
|
||||
onemol->x[i][2] -= onemol->com[2];
|
||||
}
|
||||
onemols[imol]->com[0] = 0;
|
||||
onemols[imol]->com[1] = 0;
|
||||
onemols[imol]->com[2] = 0;
|
||||
onemol->com[0] = 0;
|
||||
onemol->com[1] = 0;
|
||||
onemol->com[2] = 0;
|
||||
|
||||
} else gas_mass = atom->mass[nwidom_type];
|
||||
|
||||
if (gas_mass <= 0.0)
|
||||
error->all(FLERR,"Illegal fix widom gas mass <= 0");
|
||||
if (gas_mass <= 0.0) error->all(FLERR,"Illegal fix widom gas mass <= 0");
|
||||
|
||||
// check that no deletable atoms are in atom->firstgroup
|
||||
// deleting such an atom would not leave firstgroup atoms first
|
||||
|
|
@ -443,7 +434,6 @@ void FixWidom::pre_exchange()
|
|||
atom->nghost = 0;
|
||||
comm->borders();
|
||||
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
|
||||
update_gas_atoms_list();
|
||||
|
||||
if (full_flag) {
|
||||
energy_stored = energy_full();
|
||||
|
|
@ -525,12 +515,12 @@ void FixWidom::attempt_atomic_insertion()
|
|||
if (!domain->inside(coord))
|
||||
error->one(FLERR,"Fix widom put atom outside box");
|
||||
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
|
||||
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
|
||||
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
|
||||
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
|
||||
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
|
||||
} else {
|
||||
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
|
||||
}
|
||||
|
||||
if (proc_flag) {
|
||||
|
|
@ -545,9 +535,7 @@ void FixWidom::attempt_atomic_insertion()
|
|||
double incr_chem_pot = (inst_chem_pot - ave_widom_chemical_potential);
|
||||
ave_widom_chemical_potential += incr_chem_pot / (imove + 1);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -571,13 +559,13 @@ void FixWidom::attempt_molecule_insertion()
|
|||
com_coord[2] = region_zlo + random_equal->uniform() *
|
||||
(region_zhi-region_zlo);
|
||||
while (region->match(com_coord[0],com_coord[1],
|
||||
com_coord[2]) == 0) {
|
||||
com_coord[2]) == 0) {
|
||||
com_coord[0] = region_xlo + random_equal->uniform() *
|
||||
(region_xhi-region_xlo);
|
||||
(region_xhi-region_xlo);
|
||||
com_coord[1] = region_ylo + random_equal->uniform() *
|
||||
(region_yhi-region_ylo);
|
||||
(region_yhi-region_ylo);
|
||||
com_coord[2] = region_zlo + random_equal->uniform() *
|
||||
(region_zhi-region_zlo);
|
||||
(region_zhi-region_zlo);
|
||||
region_attempt++;
|
||||
if (region_attempt >= max_region_attempts) return;
|
||||
}
|
||||
|
|
@ -622,7 +610,7 @@ void FixWidom::attempt_molecule_insertion()
|
|||
auto procflag = new bool[natoms_per_molecule];
|
||||
|
||||
for (int i = 0; i < natoms_per_molecule; i++) {
|
||||
MathExtra::matvec(rotmat,onemols[imol]->x[i],molcoords[i]);
|
||||
MathExtra::matvec(rotmat,onemol->x[i],molcoords[i]);
|
||||
molcoords[i][0] += com_coord[0];
|
||||
molcoords[i][1] += com_coord[1];
|
||||
molcoords[i][2] += com_coord[2];
|
||||
|
|
@ -641,40 +629,37 @@ void FixWidom::attempt_molecule_insertion()
|
|||
procflag[i] = false;
|
||||
if (triclinic == 0) {
|
||||
if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] &&
|
||||
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
|
||||
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) procflag[i] = true;
|
||||
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
|
||||
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) procflag[i] = true;
|
||||
} else {
|
||||
domain->x2lamda(xtmp,lamda);
|
||||
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) procflag[i] = true;
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) procflag[i] = true;
|
||||
}
|
||||
|
||||
if (procflag[i]) {
|
||||
int ii = -1;
|
||||
if (onemols[imol]->qflag == 1) {
|
||||
ii = atom->nlocal + atom->nghost;
|
||||
if (ii >= atom->nmax) atom->avec->grow(0);
|
||||
atom->q[ii] = onemols[imol]->q[i];
|
||||
if (onemol->qflag == 1) {
|
||||
ii = atom->nlocal + atom->nghost;
|
||||
if (ii >= atom->nmax) atom->avec->grow(0);
|
||||
atom->q[ii] = onemol->q[i];
|
||||
}
|
||||
insertion_energy += energy(ii,onemols[imol]->type[i],-1,xtmp);
|
||||
insertion_energy += energy(ii,onemol->type[i],-1,xtmp);
|
||||
}
|
||||
}
|
||||
|
||||
double insertion_energy_sum = 0.0;
|
||||
MPI_Allreduce(&insertion_energy,&insertion_energy_sum,1,
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// the insertion_energy_sum is the variable with the energy of inserting one molecule
|
||||
double inst_chem_pot = exp(-insertion_energy_sum*beta);
|
||||
double incr_chem_pot = (inst_chem_pot - ave_widom_chemical_potential);
|
||||
ave_widom_chemical_potential += incr_chem_pot / (imove + 1);
|
||||
|
||||
|
||||
delete[] procflag;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -727,12 +712,12 @@ void FixWidom::attempt_atomic_insertion_full()
|
|||
if (!domain->inside(coord))
|
||||
error->one(FLERR,"Fix widom put atom outside box");
|
||||
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
|
||||
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
|
||||
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
|
||||
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
|
||||
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
|
||||
} else {
|
||||
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
|
||||
}
|
||||
|
||||
if (proc_flag) {
|
||||
|
|
@ -770,11 +755,7 @@ void FixWidom::attempt_atomic_insertion_full()
|
|||
if (proc_flag) atom->nlocal--;
|
||||
if (force->kspace) force->kspace->qsum_qsq();
|
||||
if (force->pair->tail_flag) force->pair->reinit();
|
||||
|
||||
update_gas_atoms_list();
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -803,20 +784,13 @@ void FixWidom::attempt_molecule_insertion_full()
|
|||
double com_coord[3];
|
||||
if (region) {
|
||||
int region_attempt = 0;
|
||||
com_coord[0] = region_xlo + random_equal->uniform() *
|
||||
(region_xhi-region_xlo);
|
||||
com_coord[1] = region_ylo + random_equal->uniform() *
|
||||
(region_yhi-region_ylo);
|
||||
com_coord[2] = region_zlo + random_equal->uniform() *
|
||||
(region_zhi-region_zlo);
|
||||
while (region->match(com_coord[0],com_coord[1],
|
||||
com_coord[2]) == 0) {
|
||||
com_coord[0] = region_xlo + random_equal->uniform() *
|
||||
(region_xhi-region_xlo);
|
||||
com_coord[1] = region_ylo + random_equal->uniform() *
|
||||
(region_yhi-region_ylo);
|
||||
com_coord[2] = region_zlo + random_equal->uniform() *
|
||||
(region_zhi-region_zlo);
|
||||
com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
|
||||
com_coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
|
||||
com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
|
||||
while (region->match(com_coord[0],com_coord[1], com_coord[2]) == 0) {
|
||||
com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
|
||||
com_coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
|
||||
com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
|
||||
region_attempt++;
|
||||
if (region_attempt >= max_region_attempts) return;
|
||||
}
|
||||
|
|
@ -861,7 +835,7 @@ void FixWidom::attempt_molecule_insertion_full()
|
|||
|
||||
for (int i = 0; i < natoms_per_molecule; i++) {
|
||||
double xtmp[3];
|
||||
MathExtra::matvec(rotmat,onemols[imol]->x[i],xtmp);
|
||||
MathExtra::matvec(rotmat,onemol->x[i],xtmp);
|
||||
xtmp[0] += com_coord[0];
|
||||
xtmp[1] += com_coord[1];
|
||||
xtmp[2] += com_coord[2];
|
||||
|
|
@ -876,40 +850,39 @@ void FixWidom::attempt_molecule_insertion_full()
|
|||
int proc_flag = 0;
|
||||
if (triclinic == 0) {
|
||||
if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] &&
|
||||
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
|
||||
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) proc_flag = 1;
|
||||
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
|
||||
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) proc_flag = 1;
|
||||
} else {
|
||||
domain->x2lamda(xtmp,lamda);
|
||||
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
|
||||
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
|
||||
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
|
||||
}
|
||||
|
||||
if (proc_flag) {
|
||||
atom->avec->create_atom(onemols[imol]->type[i],xtmp);
|
||||
atom->avec->create_atom(onemol->type[i],xtmp);
|
||||
int m = atom->nlocal - 1;
|
||||
|
||||
atom->image[m] = imagetmp;
|
||||
atom->molecule[m] = insertion_molecule;
|
||||
if (maxtag_all+i+1 >= MAXTAGINT)
|
||||
error->all(FLERR,"Fix widom ran out of available atom IDs");
|
||||
error->all(FLERR,"Fix widom ran out of available atom IDs");
|
||||
atom->tag[m] = maxtag_all + i + 1;
|
||||
atom->v[m][0] = 0;
|
||||
atom->v[m][1] = 0;
|
||||
atom->v[m][2] = 0;
|
||||
|
||||
atom->add_molecule_atom(onemols[imol],i,m,maxtag_all);
|
||||
atom->add_molecule_atom(onemol,i,m,maxtag_all);
|
||||
modify->create_attribute(m);
|
||||
}
|
||||
}
|
||||
|
||||
atom->natoms += natoms_per_molecule;
|
||||
if (atom->natoms < 0)
|
||||
error->all(FLERR,"Too many total atoms");
|
||||
atom->nbonds += onemols[imol]->nbonds;
|
||||
atom->nangles += onemols[imol]->nangles;
|
||||
atom->ndihedrals += onemols[imol]->ndihedrals;
|
||||
atom->nimpropers += onemols[imol]->nimpropers;
|
||||
if (atom->natoms < 0) error->all(FLERR,"Too many total atoms");
|
||||
atom->nbonds += onemol->nbonds;
|
||||
atom->nangles += onemol->nangles;
|
||||
atom->ndihedrals += onemol->ndihedrals;
|
||||
atom->nimpropers += onemol->nimpropers;
|
||||
if (atom->map_style != Atom::MAP_NONE) atom->map_init();
|
||||
atom->nghost = 0;
|
||||
if (triclinic) domain->x2lamda(atom->nlocal);
|
||||
|
|
@ -924,10 +897,10 @@ void FixWidom::attempt_molecule_insertion_full()
|
|||
double incr_chem_pot = (inst_chem_pot - ave_widom_chemical_potential);
|
||||
ave_widom_chemical_potential += incr_chem_pot / (imove + 1);
|
||||
|
||||
atom->nbonds -= onemols[imol]->nbonds;
|
||||
atom->nangles -= onemols[imol]->nangles;
|
||||
atom->ndihedrals -= onemols[imol]->ndihedrals;
|
||||
atom->nimpropers -= onemols[imol]->nimpropers;
|
||||
atom->nbonds -= onemol->nbonds;
|
||||
atom->nangles -= onemol->nangles;
|
||||
atom->ndihedrals -= onemol->ndihedrals;
|
||||
atom->nimpropers -= onemol->nimpropers;
|
||||
atom->natoms -= natoms_per_molecule;
|
||||
|
||||
int i = 0;
|
||||
|
|
@ -939,9 +912,6 @@ void FixWidom::attempt_molecule_insertion_full()
|
|||
}
|
||||
if (force->kspace) force->kspace->qsum_qsq();
|
||||
if (force->pair->tail_flag) force->pair->reinit();
|
||||
|
||||
update_gas_atoms_list();
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -1061,92 +1031,6 @@ double FixWidom::energy_full()
|
|||
return total_energy;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
update the list of gas atoms
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixWidom::update_gas_atoms_list()
|
||||
{
|
||||
int nlocal = atom->nlocal;
|
||||
int *mask = atom->mask;
|
||||
tagint *molecule = atom->molecule;
|
||||
double **x = atom->x;
|
||||
|
||||
if (atom->nmax > widom_nmax) {
|
||||
memory->sfree(local_gas_list);
|
||||
widom_nmax = atom->nmax;
|
||||
local_gas_list = (int *) memory->smalloc(widom_nmax*sizeof(int),
|
||||
"Widom:local_gas_list");
|
||||
}
|
||||
|
||||
ngas_local = 0;
|
||||
|
||||
if (region) {
|
||||
|
||||
if (exchmode == EXCHMOL) {
|
||||
|
||||
tagint maxmol = 0;
|
||||
for (int i = 0; i < nlocal; i++) maxmol = MAX(maxmol,molecule[i]);
|
||||
tagint maxmol_all;
|
||||
MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
|
||||
auto comx = new double[maxmol_all];
|
||||
auto comy = new double[maxmol_all];
|
||||
auto comz = new double[maxmol_all];
|
||||
for (int imolecule = 0; imolecule < maxmol_all; imolecule++) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (molecule[i] == imolecule) {
|
||||
mask[i] |= molecule_group_bit;
|
||||
} else {
|
||||
mask[i] &= molecule_group_inversebit;
|
||||
}
|
||||
}
|
||||
double com[3];
|
||||
com[0] = com[1] = com[2] = 0.0;
|
||||
group->xcm(molecule_group,gas_mass,com);
|
||||
|
||||
// remap unwrapped com into periodic box
|
||||
|
||||
domain->remap(com);
|
||||
comx[imolecule] = com[0];
|
||||
comy[imolecule] = com[1];
|
||||
comz[imolecule] = com[2];
|
||||
}
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (region->match(comx[molecule[i]],
|
||||
comy[molecule[i]],comz[molecule[i]]) == 1) {
|
||||
local_gas_list[ngas_local] = i;
|
||||
ngas_local++;
|
||||
}
|
||||
}
|
||||
}
|
||||
delete[] comx;
|
||||
delete[] comy;
|
||||
delete[] comz;
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (region->match(x[i][0],x[i][1],x[i][2]) == 1) {
|
||||
local_gas_list[ngas_local] = i;
|
||||
ngas_local++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
local_gas_list[ngas_local] = i;
|
||||
ngas_local++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(&ngas_local,&ngas,1,MPI_INT,MPI_SUM,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return acceptance ratios
|
||||
------------------------------------------------------------------------- */
|
||||
|
|
|
|||
|
|
@ -40,7 +40,6 @@ class FixWidom : public Fix {
|
|||
double energy(int, int, tagint, double *);
|
||||
double molecule_energy(tagint);
|
||||
double energy_full();
|
||||
void update_gas_atoms_list();
|
||||
double compute_vector(int) override;
|
||||
double memory_usage() override;
|
||||
void write_restart(FILE *) override;
|
||||
|
|
@ -53,10 +52,7 @@ class FixWidom : public Fix {
|
|||
int exclusion_group, exclusion_group_bit;
|
||||
int nwidom_type, nevery, seed;
|
||||
int ninsertions;
|
||||
int ngas; // # of gas atoms on all procs
|
||||
int ngas_local; // # of gas atoms on this proc
|
||||
int exchmode; // exchange ATOM or MOLECULE
|
||||
int movemode; // move ATOM or MOLECULE
|
||||
class Region *region; // widom region
|
||||
char *idregion; // widom region id
|
||||
bool charge_flag; // true if user specified atomic charge
|
||||
|
|
@ -71,14 +67,13 @@ class FixWidom : public Fix {
|
|||
int max_region_attempts;
|
||||
double gas_mass;
|
||||
double insertion_temperature;
|
||||
double beta, sigma, volume;
|
||||
double beta, volume;
|
||||
double charge;
|
||||
double xlo, xhi, ylo, yhi, zlo, zhi;
|
||||
double region_xlo, region_xhi, region_ylo, region_yhi, region_zlo, region_zhi;
|
||||
double region_volume;
|
||||
double energy_stored; // full energy of old/current configuration
|
||||
double *sublo, *subhi;
|
||||
int *local_gas_list;
|
||||
double **cutsq;
|
||||
double **molcoords;
|
||||
double *molq;
|
||||
|
|
@ -93,8 +88,7 @@ class FixWidom : public Fix {
|
|||
|
||||
class Atom *model_atom;
|
||||
|
||||
class Molecule **onemols;
|
||||
int imol, nmol;
|
||||
class Molecule *onemol;
|
||||
int triclinic; // 0 = orthog box, 1 = triclinic
|
||||
|
||||
class Compute *c_pe;
|
||||
|
|
|
|||
|
|
@ -56,9 +56,9 @@ FixNVTSllodOMP::FixNVTSllodOMP(LAMMPS *lmp, int narg, char **arg) :
|
|||
// id = fix-ID + temp
|
||||
|
||||
id_temp = utils::strdup(std::string(id) + "_temp");
|
||||
modify->add_compute(fmt::format("{} {} temp/deform",
|
||||
id_temp,group->names[igroup]));
|
||||
modify->add_compute(fmt::format("{} {} temp/deform",id_temp,group->names[igroup]));
|
||||
tcomputeflag = 1;
|
||||
nondeformbias = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -79,8 +79,7 @@ void FixNVTSllodOMP::init()
|
|||
for (i = 0; i < modify->nfix; i++)
|
||||
if (utils::strmatch(modify->fix[i]->style,"^deform")) {
|
||||
if ((dynamic_cast<FixDeform *>(modify->fix[i]))->remapflag != Domain::V_REMAP)
|
||||
error->all(FLERR,"Using fix nvt/sllod/omp with inconsistent fix "
|
||||
"deform remap option");
|
||||
error->all(FLERR,"Using fix nvt/sllod/omp with inconsistent fix deform remap option");
|
||||
break;
|
||||
}
|
||||
if (i == modify->nfix)
|
||||
|
|
|
|||
|
|
@ -779,9 +779,11 @@ std::string Atom::get_style()
|
|||
std::string retval = atom_style;
|
||||
if (retval == "hybrid") {
|
||||
auto avec_hybrid = dynamic_cast<AtomVecHybrid *>(avec);
|
||||
for (int i = 0; i < avec_hybrid->nstyles; i++) {
|
||||
retval += ' ';
|
||||
retval += avec_hybrid->keywords[i];
|
||||
if (avec_hybrid) {
|
||||
for (int i = 0; i < avec_hybrid->nstyles; i++) {
|
||||
retval += ' ';
|
||||
retval += avec_hybrid->keywords[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
return retval;
|
||||
|
|
|
|||
88
src/bond.cpp
88
src/bond.cpp
|
|
@ -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
|
||||
------------------------------------------------------------------------- */
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
------------------------------------------------------------------------- */
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -51,9 +51,9 @@ FixNVTSllod::FixNVTSllod(LAMMPS *lmp, int narg, char **arg) :
|
|||
// id = fix-ID + temp
|
||||
|
||||
id_temp = utils::strdup(std::string(id) + "_temp");
|
||||
modify->add_compute(fmt::format("{} {} temp/deform",
|
||||
id_temp,group->names[igroup]));
|
||||
modify->add_compute(fmt::format("{} {} temp/deform",id_temp,group->names[igroup]));
|
||||
tcomputeflag = 1;
|
||||
nondeformbias = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -74,8 +74,7 @@ void FixNVTSllod::init()
|
|||
for (i = 0; i < modify->nfix; i++)
|
||||
if (strncmp(modify->fix[i]->style,"deform",6) == 0) {
|
||||
if ((dynamic_cast<FixDeform *>(modify->fix[i]))->remapflag != Domain::V_REMAP)
|
||||
error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform "
|
||||
"remap option");
|
||||
error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform remap option");
|
||||
break;
|
||||
}
|
||||
if (i == modify->nfix)
|
||||
|
|
|
|||
|
|
@ -84,7 +84,7 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
|
|||
|
||||
for (int ifield = 0; ifield < nfield; ifield++) {
|
||||
int columns = 0; // set in case fieldname not recognized by pstyle
|
||||
void *pvoid = pstyle->extract_peratom(fieldname[ifield],columns);
|
||||
pstyle->extract_peratom(fieldname[ifield],columns);
|
||||
if (columns) ncols += columns;
|
||||
else ncols++;
|
||||
if (trigger[ifield]) {
|
||||
|
|
@ -120,7 +120,8 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
|
|||
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
|
||||
// zero the vector/array since a variable may access it before first ru
|
||||
// initialize lasttime so step 0 will trigger/extract
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
|
|
@ -132,6 +133,8 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
|
|||
for (int m = 0; m < ncols; m++)
|
||||
array[i][m] = 0.0;
|
||||
}
|
||||
|
||||
lasttime = -1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -188,6 +191,13 @@ void FixPair::setup(int vflag)
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPair::min_setup(int vflag)
|
||||
{
|
||||
setup(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPair::setup_pre_force(int vflag)
|
||||
{
|
||||
pre_force(vflag);
|
||||
|
|
@ -195,11 +205,13 @@ void FixPair::setup_pre_force(int vflag)
|
|||
|
||||
/* ----------------------------------------------------------------------
|
||||
trigger pair style computation on steps which are multiples of Nevery
|
||||
lasttime prevents mulitiple triggers by min linesearch on same iteration
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixPair::pre_force(int /*vflag*/)
|
||||
{
|
||||
if (update->ntimestep % nevery) return;
|
||||
if (update->ntimestep == lasttime) return;
|
||||
|
||||
// set pair style triggers
|
||||
|
||||
|
|
@ -215,12 +227,15 @@ void FixPair::min_pre_force(int vflag)
|
|||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
extract results from pair style
|
||||
extract results from pair style on steps which are multiples of Nevery
|
||||
lasttime prevents mulitiple extracts by min linesearch on same iteration
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixPair::post_force(int /*vflag*/)
|
||||
{
|
||||
if (update->ntimestep % nevery) return;
|
||||
if (update->ntimestep == lasttime) return;
|
||||
lasttime = update->ntimestep;
|
||||
|
||||
// extract pair style fields one by one
|
||||
// store their values in this fix
|
||||
|
|
|
|||
|
|
@ -31,6 +31,7 @@ class FixPair : public Fix {
|
|||
int setmask() override;
|
||||
void init() override;
|
||||
void setup(int) override;
|
||||
void min_setup(int) override;
|
||||
void setup_pre_force(int) override;
|
||||
void pre_force(int) override;
|
||||
void min_pre_force(int) override;
|
||||
|
|
@ -46,6 +47,7 @@ class FixPair : public Fix {
|
|||
|
||||
private:
|
||||
int nevery,nfield,ncols;
|
||||
bigint lasttime;
|
||||
char *pairname;
|
||||
char **fieldname,**triggername;
|
||||
int *trigger;
|
||||
|
|
|
|||
|
|
@ -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");
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -99,9 +99,8 @@ NeighRequest::NeighRequest(LAMMPS *_lmp) : Pointers(_lmp)
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
NeighRequest::NeighRequest(LAMMPS *_lmp, int idx, void *ptr, int num) : NeighRequest(_lmp)
|
||||
NeighRequest::NeighRequest(LAMMPS *_lmp, void *ptr, int num) : NeighRequest(_lmp)
|
||||
{
|
||||
index = idx;
|
||||
requestor = ptr;
|
||||
requestor_instance = num;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -29,7 +29,6 @@ class NeighRequest : protected Pointers {
|
|||
friend class FixIntel;
|
||||
|
||||
protected:
|
||||
int index; // index of which neigh request this is
|
||||
void *requestor; // class that made request
|
||||
int requestor_instance; // instance of that class (only Fix, Compute, Pair)
|
||||
int id; // ID of request as stored by requestor
|
||||
|
|
@ -122,7 +121,7 @@ class NeighRequest : protected Pointers {
|
|||
// methods
|
||||
public:
|
||||
NeighRequest(class LAMMPS *);
|
||||
NeighRequest(class LAMMPS *, int, void *, int);
|
||||
NeighRequest(class LAMMPS *, void *, int);
|
||||
NeighRequest(NeighRequest *);
|
||||
~NeighRequest() override;
|
||||
|
||||
|
|
|
|||
|
|
@ -2163,13 +2163,13 @@ int Neighbor::request(void *requestor, int instance)
|
|||
memory->srealloc(requests,maxrequest*sizeof(NeighRequest *), "neighbor:requests");
|
||||
}
|
||||
|
||||
requests[nrequest] = new NeighRequest(lmp, nrequest, requestor, instance);
|
||||
requests[nrequest] = new NeighRequest(lmp, requestor, instance);
|
||||
nrequest++;
|
||||
return nrequest-1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
called by other classes to request a pairwise neighbor list
|
||||
add_request() is called by other classes to request a pairwise neighbor list
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
NeighRequest *Neighbor::add_request(Pair *requestor, int flags)
|
||||
|
|
@ -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
|
||||
------------------------------------------------------------------------- */
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -28,6 +28,7 @@ class NTopo : protected Pointers {
|
|||
|
||||
virtual void build() = 0;
|
||||
|
||||
void add_temporary_bond(int, int, int);
|
||||
double memory_usage();
|
||||
|
||||
protected:
|
||||
|
|
|
|||
|
|
@ -33,7 +33,7 @@ if(CMAKE_Fortran_COMPILER)
|
|||
add_library(flammps STATIC ${LAMMPS_FORTRAN_MODULE})
|
||||
|
||||
add_executable(test_fortran_create wrap_create.cpp test_fortran_create.f90)
|
||||
target_link_libraries(test_fortran_create PRIVATE lammps MPI::MPI_Fortran GTest::GTestMain)
|
||||
target_link_libraries(test_fortran_create PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
|
||||
target_include_directories(test_fortran_create PRIVATE "${LAMMPS_SOURCE_DIR}/../fortran")
|
||||
add_test(NAME FortranOpen COMMAND test_fortran_create)
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue