tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2788 from jtclemm/track-contacts

Browse Source 
Track properties of pair interactions
This commit is contained in:
Axel Kohlmeyer
2021-07-20 19:38:21 -04:00
committed by GitHub
parent 2718596d65 63794aa698
commit 96389213ab
24 changed files with 9625 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_fix.rst
View File

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

1
doc/src/Commands_pair.rst
View File

@ -273,6 +273,7 @@ OPT.
* :doc:`tip4p/cut (o) <pair_coul>`
* :doc:`tip4p/long (o) <pair_coul>`
* :doc:`tip4p/long/soft (o) <pair_fep_soft>`
* :doc:`tracker <pair_tracker>`
* :doc:`tri/lj <pair_tri_lj>`
* :doc:`ufm (got) <pair_ufm>`
* :doc:`vashishta (gko) <pair_vashishta>`

2
doc/src/Examples.rst
View File

@ -150,6 +150,8 @@ Lowercase directories
+-------------+------------------------------------------------------------------+
| threebody | regression test input for a variety of manybody potentials |
+-------------+------------------------------------------------------------------+
| tracker | track interactions in LJ melt |
+-------------+------------------------------------------------------------------+
| vashishta | use of the Vashishta potential |
+-------------+------------------------------------------------------------------+
| voronoi | Voronoi tesselation via compute voronoi/atom command |

12
doc/src/dump_modify.rst
View File

@ -66,7 +66,7 @@ Syntax
*unwrap* arg = *yes* or *no*
* these keywords apply only to the *image* and *movie* :doc:`styles <dump_image>`
* keyword = *acolor* or *adiam* or *amap* or *backcolor* or *bcolor* or *bdiam* or *boxcolor* or *color* or *bitrate* or *framerate*
* keyword = *acolor* or *adiam* or *amap* or *backcolor* or *bcolor* or *bdiam* or *boxcolor* or *color* or *bitrate* or *framerate* or *header*
.. parsed-literal::
@ -113,6 +113,9 @@ Syntax
rate = target bitrate for movie in kbps
*framerate* arg = fps
fps = frames per second for movie
*header* arg = *yes* or *no*
*yes* to write the header
*no* to not write the header
* these keywords apply only to the */gz* and */zstd* dump styles
* keyword = *compression_level*
@ -977,6 +980,13 @@ images less frequently.
----------
The *header* keyword toggles whether the dump file will include a header.
Excluding a header will reduce the size of the dump file for fixes such as
:doc:`fix pair/tracker <fix_pair_tracker>` which do not require the information
typically written to the header.
----------
The COMPRESS package offers both GZ and Zstd compression variants of styles
atom, custom, local, cfg, and xyz. When using these styles the compression
level can be controlled by the :code:`compression_level` parameter. File names

1
doc/src/fix.rst
View File

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

4
doc/src/fix_oneway.rst
View File

@ -50,7 +50,9 @@ the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minim
Restrictions
""""""""""""
none
This fix is part of the MISC package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""

124
doc/src/fix_pair_tracker.rst Normal file
View File

@ -0,0 +1,124 @@
.. index:: fix pair/tracker
fix pair/tracker command
========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID pair/tracker N attribute1 attribute2 ... keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* pair/tracker = style name of this fix command
* N = prepare data for output every this many timesteps
* one or more attributes may be appended
.. parsed-literal::
possible attributes = id1 id2 time/created time/broken time/total
rmin rave x y z
.. parsed-literal::
id1, id2 = IDs of the 2 atoms in each pair interaction
time/created = the time that the 2 atoms began interacting
time/broken = the time that the 2 atoms stopped interacting
time/total = the total time the 2 atoms interacted
r/min = the minimum radial distance between the 2 atoms during the interaction
r/ave = the average radial distance between the 2 atoms during the interaction
x, y, z = the center of mass position of the 2 atoms when they stopped interacting
* zero or more keyword/value pairs may be appended
* keyword = *time/min* or *type/include*
.. parsed-literal::
*time/min* value = T
T = minimum interaction time
*type/include* value = arg1 arg2
arg = separate lists of types (see below)
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all pair/tracker 1000 id1 id2 time/min 100
fix 1 all pair/tracker 1000 time/created time/broken type/include 1 * type/include 2 3,4
Description
"""""""""""
Tracks properties of pairwise interactions between two atoms and records data
whenever the atoms move beyond the interaction cutoff.
Must be used in conjunction with :doc:`pair tracker <pair_tracker>`.
Data is accumulated over a span of *N* timesteps before being deleted.
The number of datums generated, aggregated across all processors, equals
the number of broken interactions. Interactions are only included if both
atoms are included in the specified fix group. Additional filters can be
applied using the *time/min* or *type/include* keywords described below.
.. note::
For extremely long-lived interactions, the calculation of *r/ave* may not be
correct due to double overflow.
The *time/min* keyword sets a minimum amount of time that an interaction must
persist to be included. This setting can be used to censor short-lived interactions.
The *type/include* keyword filters interactions based on the types of the two atoms.
Data is only saved for interactions between atoms with types in the two lists.
Each list consists of a series of type
ranges separated by commas. The range can be specified as a
single numeric value, or a wildcard asterisk can be used to specify a range
of values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For
example, if M = the number of atom types, then an asterisk with no numeric
values means all types from 1 to M. A leading asterisk means all types
from 1 to n (inclusive). A trailing asterisk means all types from n to M
(inclusive). A middle asterisk means all types from m to n (inclusive).
Multiple *type/include* keywords may be added.
----------
Restart, fix_modify, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command.
Output info
"""""""""""
This compute 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>` doc page for an overview of LAMMPS output
options.
The vector or array values will be doubles that correspond to the
specified attribute.
Restrictions
""""""""""""
Must be used in conjunction with :doc:`pair style tracker <pair_tracker>`.
This fix is part of the MISC package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`pair tracker <pair_tracker>`
Default
"""""""
none

6
doc/src/pair_granular.rst
View File

@ -665,12 +665,6 @@ then LAMMPS will use that cutoff for the specified atom type
combination, and automatically set pairwise cutoffs for the remaining
atom types.
If two particles are moving away from each other while in contact, there
is a possibility that the particles could experience an effective attractive
force due to damping. If the *limit_damping* keyword is used, this option
will zero out the normal component of the force if there is an effective
attractive force. This keyword cannot be used with the JKR or DMT models.
----------
.. include:: accel_styles.rst

1
doc/src/pair_style.rst
View File

@ -338,6 +338,7 @@ accelerated styles exist.
* :doc:`tip4p/cut <pair_coul>` - Coulomb for TIP4P water w/out LJ
* :doc:`tip4p/long <pair_coul>` - long-range Coulomb for TIP4P water w/out LJ
* :doc:`tip4p/long/soft <pair_fep_soft>` -
* :doc:`tracker <pair_tracker>` - monitor information about pairwise interactions
* :doc:`tri/lj <pair_tri_lj>` - LJ potential between triangles
* :doc:`ufm <pair_ufm>` -
* :doc:`vashishta <pair_vashishta>` - Vashishta 2-body and 3-body potential

97
doc/src/pair_tracker.rst Normal file
View File

@ -0,0 +1,97 @@
.. index:: pair_style tracker
pair_style tracker command
==========================
Syntax
""""""
.. code-block:: LAMMPS
pair_style tracker keyword
* zero or more keyword/arg pairs may be appended
* keyword = *finite*
.. parsed-literal::
*finite* value = none
pair style uses atomic diameters to identify contacts
Examples
""""""""
.. code-block:: LAMMPS
pair_style hybrid/overlay tracker ...
pair_coeff 1 1 tracker 2.0
pair_style hybrid/overlay tracker finite ...
pair_coeff * * tracker
fix 1 all pair/tracker 1000 time/created time/broken
dump 1 all local 1000 dump.local f_1[1] f_1[2]
dump_modify 1 write_header no
Description
"""""""""""
Style *tracker* monitors information about pairwise interactions.
It does not calculate any forces on atoms.
:doc:`Pair hybrid/overlay <pair_hybrid>` can be used to combine this pair
style with another pair style. Style *tracker* must be used in conjunction
with about :doc:`fix pair_tracker <fix_pair_tracker>` which contains
information on what data can be output.
If the *finite* keyword is not defined, the following coefficients must be
defined for each pair of atom types via the :doc:`pair_coeff <pair_coeff>`
command as in the examples above, or in the data file or restart files
read by the :doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands, or by mixing as described below:
* cutoff (distance units)
If the *finite* keyword is defined, no coefficients may be defined.
Interaction cutoffs are alternatively calculated based on the
diameter of finite particles.
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
For atom type pairs I,J and I != J, the cutoff coefficient and cutoff
distance for this pair style can be mixed. The cutoff is always mixed via a
*geometric* rule. The cutoff is mixed according to the pair_modify
mix value. The default mix value is *geometric*\ . See the
"pair_modify" command for details.
This pair style writes its information to :doc:`binary restart files <restart>`, so
pair_style and pair_coeff commands do not need
to be specified in an input script that reads a restart file.
The :doc:`pair_modify <pair_modify>` shift, table, and tail options
are not relevant for this pair style.
----------
Restrictions
""""""""""""
A corresponding :doc:`fix pair_tracker <fix_pair_tracker>` must be defined
to use this pair style.
This pair style is currently incompatible with granular pair styles that extend
beyond the contact (e.g. JKR and DMT).
This fix is part of the MISC package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`fix pair_tracker <fix_pair_tracker>`
Default
"""""""
none

8122
examples/tracker/contact_history.dat Normal file
View File

File diff suppressed because it is too large Load Diff

34
examples/tracker/in.track Normal file
View File

@ -0,0 +1,34 @@
# 3d Lennard-Jones melt with tracking
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 4 0 4 0 4
create_box 2 box
create_atoms 2 box
mass 1 1.0
mass 2 1.0
velocity all create 2.0 87287 loop geom
pair_style hybrid/overlay lj/cut 2.5 tracker
pair_coeff * * lj/cut 1.0 1.0 2.5
pair_coeff * * tracker 2.5
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
fix 2 all pair/tracker 1000 id1 id2 time/created time/broken time/total time/min 0.5
fix 3 all ave/histo 1000 1 1000 0 6 30 f_2[5] mode vector kind local file lifetime_hist.dat
dump 1 all local 1000 contact_history.dat f_2[1] f_2[2] f_2[3] f_2[4]
dump_modify 1 header no
# compute 1 all property/local patom1 patom2
# dump 2 all local 1 pairs.dat c_1[1] c_1[2]
# dump 3 all atom 1 atoms.dat
thermo 50
run 1000 #0

65
examples/tracker/lifetime_hist.dat Normal file
View File

@ -0,0 +1,65 @@
# Histogrammed data for fix 3
# TimeStep Number-of-bins Total-counts Missing-counts Min-value Max-value
# Bin Coord Count Count/Total
0 30 0 0 1e+20 -1e+20
1 0.1 0 0
2 0.3 0 0
3 0.5 0 0
4 0.7 0 0
5 0.9 0 0
6 1.1 0 0
7 1.3 0 0
8 1.5 0 0
9 1.7 0 0
10 1.9 0 0
11 2.1 0 0
12 2.3 0 0
13 2.5 0 0
14 2.7 0 0
15 2.9 0 0
16 3.1 0 0
17 3.3 0 0
18 3.5 0 0
19 3.7 0 0
20 3.9 0 0
21 4.1 0 0
22 4.3 0 0
23 4.5 0 0
24 4.7 0 0
25 4.9 0 0
26 5.1 0 0
27 5.3 0 0
28 5.5 0 0
29 5.7 0 0
30 5.9 0 0
1000 30 8122 0 0.5 5
1 0.1 0 0
2 0.3 0 0
3 0.5 910 0.112041
4 0.7 1253 0.154272
5 0.9 953 0.117336
6 1.1 747 0.0919724
7 1.3 559 0.0688254
8 1.5 501 0.0616843
9 1.7 421 0.0518345
10 1.9 356 0.0438316
11 2.1 300 0.0369367
12 2.3 281 0.0345974
13 2.5 242 0.0297956
14 2.7 226 0.0278257
15 2.9 175 0.0215464
16 3.1 168 0.0206846
17 3.3 162 0.0199458
18 3.5 129 0.0158828
19 3.7 151 0.0185915
20 3.9 137 0.0168678
21 4.1 98 0.012066
22 4.3 104 0.0128047
23 4.5 83 0.0102192
24 4.7 77 0.00948042
25 4.9 86 0.0105885
26 5.1 3 0.000369367
27 5.3 0 0
28 5.5 0 0
29 5.7 0 0
30 5.9 0 0

107
examples/tracker/log.8Apr21.tracker.g++.1 Normal file
View File

@ -0,0 +1,107 @@
LAMMPS (8 Apr 2021)
# 3d Lennard-Jones melt with tracking
units lj
atom_style atomic
lattice fcc 0.8442
Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962
region box block 0 4 0 4 0 4
create_box 2 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (6.7183848 6.7183848 6.7183848)
1 by 1 by 1 MPI processor grid
create_atoms 2 box
Created 256 atoms
create_atoms CPU = 0.001 seconds
mass 1 1.0
mass 2 1.0
velocity all create 2.0 87287 loop geom
pair_style hybrid/overlay lj/cut 2.5 tracker
pair_coeff * * lj/cut 1.0 1.0 2.5
pair_coeff * * tracker 2.5
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
fix 2 all pair/tracker 1000 id1 id2 time/created time/broken time/total time/min 0.5
fix 3 all ave/histo 1000 1 1000 0 6 30 f_2[5] mode vector kind local file lifetime_hist.dat
dump 1 all local 1000 contact_history.dat f_2[1] f_2[2] f_2[3] f_2[4]
dump_modify 1 header no
# compute 1 all property/local patom1 patom2
# dump 2 all local 1 pairs.dat c_1[1] c_1[2]
# dump 3 all atom 1 atoms.dat
thermo 50
run 1000 #0
Neighbor list info ...
update every 20 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
(2) pair tracker, perpetual
attributes: half, newton on, history
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 10.46 | 10.46 | 10.46 Mbytes
Step Temp E_pair E_mol TotEng Press
0 2 -6.7733681 0 -3.7850868 -4.5535126
50 1.0540428 -5.3670859 0 -3.7921977 2.3386375
100 1.0402713 -5.3493439 0 -3.7950323 2.4553748
150 1.0570745 -5.3738961 0 -3.7944781 2.3767396
200 1.0431846 -5.3518647 0 -3.7932002 2.5010135
250 1.070121 -5.3902744 0 -3.791363 2.4908658
300 1.0667723 -5.3866302 0 -3.7927224 2.3589344
350 1.000601 -5.2859643 0 -3.7909257 2.9065274
400 0.99256113 -5.2738812 0 -3.7908553 2.8595867
450 1.0482542 -5.357452 0 -3.7912128 2.4707397
500 1.0196176 -5.3123538 0 -3.7889017 2.7230338
550 0.98274535 -5.2586303 0 -3.7902706 2.9156947
600 1.0683914 -5.3863229 0 -3.789996 2.3002719
650 1.0130779 -5.303917 0 -3.7902362 2.8726423
700 1.0583333 -5.3737358 0 -3.792437 2.5770307
750 0.98274506 -5.2612464 0 -3.792887 2.9447027
800 1.0294191 -5.332001 0 -3.7939042 2.5293193
850 0.99240027 -5.2735754 0 -3.7907899 2.7672711
900 1.0293488 -5.3306241 0 -3.7926323 2.6054041
950 0.97137182 -5.2424403 0 -3.7910742 3.129989
1000 1.0009431 -5.2864286 0 -3.7908788 2.7536598
Loop time of 0.310363 on 1 procs for 1000 steps with 256 atoms
Performance: 1391918.252 tau/day, 3222.033 timesteps/s
100.7% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.22132 | 0.22132 | 0.22132 | 0.0 | 71.31
Neigh | 0.03458 | 0.03458 | 0.03458 | 0.0 | 11.14
Comm | 0.0087938 | 0.0087938 | 0.0087938 | 0.0 | 2.83
Output | 0.014075 | 0.014075 | 0.014075 | 0.0 | 4.54
Modify | 0.02818 | 0.02818 | 0.02818 | 0.0 | 9.08
Other | | 0.003411 | | | 1.10
Nlocal: 256.000 ave 256 max 256 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1327.00 ave 1327 max 1327 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 9612.00 ave 9612 max 9612 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 9612
Ave neighs/atom = 37.546875
Neighbor list builds = 50
Dangerous builds not checked
Total wall time: 0:00:00

4
src/.gitignore vendored
View File

@ -722,6 +722,8 @@
/fix_orient_eco.h
/fix_orient_fcc.cpp
/fix_orient_fcc.h
/fix_pair_tracker.cpp
/fix_pair_tracker.h
/fix_peri_neigh.cpp
/fix_peri_neigh.h
/fix_phonon.cpp
@ -1171,6 +1173,8 @@
/pair_tip4p_long.h
/pair_tip4p_long_soft.cpp
/pair_tip4p_long_soft.h
/pair_tracker.h
/pair_tracker.cpp
/pair_tri_lj.cpp
/pair_tri_lj.h
/pair_yukawa_colloid.cpp

369
src/MISC/fix_pair_tracker.cpp Normal file
View File

@ -0,0 +1,369 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "fix_pair_tracker.h"
#include "atom.h"
#include "atom_vec.h"
#include "error.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "tokenizer.h"
#include "update.h"
#include <string.h>
using namespace LAMMPS_NS;
using namespace FixConst;
#define DELTA 1000
/* ---------------------------------------------------------------------- */
FixPairTracker::FixPairTracker(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), nvalues(0), vector(NULL), array(NULL), pack_choice(NULL)
{
if (narg < 3) error->all(FLERR, "Illegal fix pair/tracker command");
local_flag = 1;
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
if (nevery <= 0) error->all(FLERR, "Illegal fix pair/tracker command");
local_freq = nevery;
// If optional arguments included, this will be oversized
nvalues = narg - 4;
pack_choice = new FnPtrPack[nvalues];
tmin = -1;
type_filter = nullptr;
int iarg = 4;
nvalues = 0;
while (iarg < narg) {
if (strcmp(arg[iarg], "id1") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_id1;
} else if (strcmp(arg[iarg], "id2") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_id2;
} else if (strcmp(arg[iarg], "time/created") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_time_created;
} else if (strcmp(arg[iarg], "time/broken") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_time_broken;
} else if (strcmp(arg[iarg], "time/total") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_time_total;
} else if (strcmp(arg[iarg], "x") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_x;
} else if (strcmp(arg[iarg], "y") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_y;
} else if (strcmp(arg[iarg], "z") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_z;
} else if (strcmp(arg[iarg], "r/min") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_rmin;
} else if (strcmp(arg[iarg], "r/ave") == 0) {
pack_choice[nvalues++] = &FixPairTracker::pack_rave;
} else if (strcmp(arg[iarg], "time/min") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command");
tmin = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg++;
} else if (strcmp(arg[iarg], "type/include") == 0) {
if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command");
int ntypes = atom->ntypes;
int i, j, itype, jtype, in, jn, infield, jnfield;
int inlo, inhi, jnlo, jnhi;
char *istr, *jstr;
if (!type_filter) {
memory->create(type_filter, ntypes + 1, ntypes + 1, "fix/pair/tracker:type_filter");
for (i = 0; i <= ntypes; i++) {
for (j = 0; j <= ntypes; j++) { type_filter[i][j] = 0; }
}
}
in = strlen(arg[iarg + 1]) + 1;
istr = new char[in];
strcpy(istr, arg[iarg + 1]);
std::vector<std::string> iwords = Tokenizer(istr, ",").as_vector();
infield = iwords.size();
jn = strlen(arg[iarg + 2]) + 1;
jstr = new char[jn];
strcpy(jstr, arg[iarg + 2]);
std::vector<std::string> jwords = Tokenizer(jstr, ",").as_vector();
jnfield = jwords.size();
for (i = 0; i < infield; i++) {
const char *ifield = iwords[i].c_str();
utils::bounds(FLERR, ifield, 1, ntypes, inlo, inhi, error);
for (j = 0; j < jnfield; j++) {
const char *jfield = jwords[j].c_str();
utils::bounds(FLERR, jfield, 1, ntypes, jnlo, jnhi, error);
for (itype = inlo; itype <= inhi; itype++) {
for (jtype = jnlo; jtype <= jnhi; jtype++) {
type_filter[itype][jtype] = 1;
type_filter[jtype][itype] = 1;
}
}
}
}
delete[] istr;
delete[] jstr;
iarg += 2;
} else
error->all(FLERR, "Invalid keyword in fix pair/tracker command");
iarg++;
}
if (nvalues == 1)
size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0;
ncount = 0;
vector = NULL;
array = NULL;
}
/* ---------------------------------------------------------------------- */
FixPairTracker::~FixPairTracker()
{
delete[] pack_choice;
memory->destroy(vector);
memory->destroy(array);
memory->destroy(type_filter);
}
/* ---------------------------------------------------------------------- */
int FixPairTracker::setmask()
{
int mask = 0;
mask |= POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::init()
{
// Set size of array/vector
ncount = 0;
if (ncount > nmax) reallocate(ncount);
size_local_rows = ncount;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::lost_contact(int i, int j, double time_tmp, double nstep_tmp, double rsum_tmp,
double rmin_tmp)
{
double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
if ((time - time_tmp) < tmin) return;
if (type_filter) {
int *type = atom->type;
if (type_filter[type[i]][type[j]] == 0) return;
}
int *mask = atom->mask;
if (!(mask[i] & groupbit)) return;
if (!(mask[j] & groupbit)) return;
if (ncount == nmax) reallocate(ncount);
index_i = i;
index_j = j;
rmin = rmin_tmp;
rsum = rsum_tmp;
time_initial = time_tmp;
nstep_initial = nstep_tmp;
// fill vector or array with local values
if (nvalues == 1) {
(this->*pack_choice[0])(0);
} else {
for (int k = 0; k < nvalues; k++) { (this->*pack_choice[k])(k); }
}
ncount += 1;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::post_force(int /*vflag*/)
{
if (update->ntimestep % nevery == 0) {
size_local_rows = ncount;
ncount = 0;
}
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::reallocate(int n)
{
// grow vector or array
while (nmax <= n) nmax += DELTA;
if (nvalues == 1) {
memory->grow(vector, nmax, "fix_pair_tracker:vector");
vector_local = vector;
} else {
memory->grow(array, nmax, nvalues, "fix_pair_tracker:array");
array_local = array;
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double FixPairTracker::memory_usage()
{
double bytes = nmax * nvalues * sizeof(double);
bytes += nmax * 2 * sizeof(int);
return bytes;
}
/* ----------------------------------------------------------------------
one method for every keyword fix pair/tracker can output
the atom property is packed into a local vector or array
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_time_created(int n)
{
if (nvalues == 1)
vector[ncount] = time_initial;
else
array[ncount][n] = time_initial;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_time_broken(int n)
{
double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
if (nvalues == 1)
vector[ncount] = time;
else
array[ncount][n] = time;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_time_total(int n)
{
double time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
if (nvalues == 1)
vector[ncount] = time - time_initial;
else
array[ncount][n] = time - time_initial;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_id1(int n)
{
tagint *tag = atom->tag;
if (nvalues == 1)
vector[ncount] = tag[index_i];
else
array[ncount][n] = tag[index_i];
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_id2(int n)
{
tagint *tag = atom->tag;
if (nvalues == 1)
vector[ncount] = tag[index_j];
else
array[ncount][n] = tag[index_j];
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_x(int n)
{
double **x = atom->x;
if (nvalues == 1)
vector[ncount] = (x[index_i][0] + x[index_j][0]) / 2;
else
array[ncount][n] = (x[index_i][0] + x[index_j][0]) / 2;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_y(int n)
{
double **x = atom->x;
if (nvalues == 1)
vector[ncount] = (x[index_i][1] + x[index_j][1]) / 2;
else
array[ncount][n] = (x[index_i][1] + x[index_j][1]) / 2;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_z(int n)
{
double **x = atom->x;
if (nvalues == 1)
vector[ncount] = (x[index_i][2] + x[index_j][2]) / 2;
else
array[ncount][n] = (x[index_i][2] + x[index_j][2]) / 2;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_rmin(int n)
{
if (nvalues == 1)
vector[ncount] = rmin;
else
array[ncount][n] = rmin;
}
/* ---------------------------------------------------------------------- */
void FixPairTracker::pack_rave(int n)
{
if (nvalues == 1)
vector[ncount] = rsum / (update->ntimestep - nstep_initial);
else
array[ncount][n] = rsum / (update->ntimestep - nstep_initial);
}

85
src/MISC/fix_pair_tracker.h Normal file
View File

@ -0,0 +1,85 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(pair/tracker,FixPairTracker);
// clang-format on
#else
#ifndef LMP_FIX_PAIR_TRACKING_H
#define LMP_FIX_PAIR_TRACKING_H
#include "fix.h"
namespace LAMMPS_NS {
class FixPairTracker : public Fix {
public:
FixPairTracker(class LAMMPS *, int, char **);
~FixPairTracker();
int setmask();
void init();
void post_force(int);
double memory_usage();
void lost_contact(int, int, double, double, double, double);
private:
int nvalues, nmax;
int index_i, index_j;
double tmin, rmin, rsum, time_initial, nstep_initial;
double *vector;
double **array;
int **type_filter;
int ncount;
void reallocate(int);
typedef void (FixPairTracker::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
void pack_id1(int);
void pack_id2(int);
void pack_time_created(int);
void pack_time_broken(int);
void pack_time_total(int);
void pack_x(int);
void pack_y(int);
void pack_z(int);
void pack_rmin(int);
void pack_rave(int);
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Invalid keyword in fix pair/tracker command
Self-explanatory.
*/

478
src/MISC/pair_tracker.cpp Normal file
View File

@ -0,0 +1,478 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "pair_tracker.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "fix.h"
#include "fix_dummy.h"
#include "fix_neigh_history.h"
#include "fix_pair_tracker.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "update.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 1;
no_virial_fdotr_compute = 1;
neighprev = 0;
history = 1;
size_history = 4;
nondefault_history_transfer = 1;
finitecutflag = 0;
// create dummy fix as placeholder for FixNeighHistory
// this is so final order of Modify:fix will conform to input script
fix_history = nullptr;
modify->add_fix("NEIGH_HISTORY_TRACK_DUMMY all DUMMY");
fix_dummy = (FixDummy *) modify->fix[modify->nfix - 1];
}
/* ---------------------------------------------------------------------- */
PairTracker::~PairTracker()
{
if (!fix_history)
modify->delete_fix("NEIGH_HISTORY_TRACK_DUMMY");
else
modify->delete_fix("NEIGH_HISTORY_TRACK");
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
delete[] onerad_dynamic;
delete[] onerad_frozen;
delete[] maxrad_dynamic;
delete[] maxrad_frozen;
}
}
/* ---------------------------------------------------------------------- */
void PairTracker::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, time;
double radi, radj, radsum, rsq, r;
int *ilist, *jlist, *numneigh, **firstneigh;
int *touch, **firsttouch;
double *data, *alldata, **firstdata;
int updateflag = 1;
if (update->setupflag) updateflag = 0;
ev_init(eflag, vflag);
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
firsttouch = fix_history->firstflag;
firstdata = fix_history->firstvalue;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (finitecutflag) radi = radius[i];
itype = type[i];
touch = firsttouch[i];
alldata = firstdata[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
r = sqrt(rsq);
if (finitecutflag) {
radj = radius[j];
radsum = radi + radj;
if (rsq >= radsum * radsum) {
data = &alldata[size_history * jj];
if (touch[jj] == 1) {
fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]);
}
touch[jj] = 0;
data[0] = 0.0; // initial time
data[1] = 0.0; // initial timestep
data[2] = 0.0; // sum of r, may overflow
data[3] = 0.0; // min of r
} else {
data = &alldata[size_history * jj];
if (touch[jj] == 0) {
time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
data[0] = time;
data[1] = (double) update->ntimestep;
data[2] = r;
data[3] = r;
} else if (updateflag) {
data[2] += r;
if (data[3] > r) data[3] = r;
}
touch[jj] = 1;
}
} else {
jtype = type[j];
if (rsq >= cutsq[itype][jtype]) {
data = &alldata[size_history * jj];
if (touch[jj] == 1) {
fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]);
}
touch[jj] = 0;
data[0] = 0.0; // initial time
data[1] = 0.0; // initial timestep
data[2] = 0.0; // sum of r, may overflow
data[3] = 0.0; // min of r
} else {
data = &alldata[size_history * jj];
if (touch[jj] == 0) {
time = update->atime + (update->ntimestep - update->atimestep) * update->dt;
data[0] = time;
data[1] = (double) update->ntimestep;
data[2] = r;
data[3] = r;
} else if (updateflag) {
data[2] += r;
if (data[3] > r) data[3] = r;
}
touch[jj] = 1;
}
}
}
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairTracker::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag, n + 1, n + 1, "pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++) setflag[i][j] = 0;
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(cut, n + 1, n + 1, "pair:cut");
onerad_dynamic = new double[n + 1];
onerad_frozen = new double[n + 1];
maxrad_dynamic = new double[n + 1];
maxrad_frozen = new double[n + 1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairTracker::settings(int narg, char **arg)
{
if (narg != 0 && narg != 1) error->all(FLERR, "Illegal pair_style command");
if (narg == 1) {
if (strcmp(arg[0], "finite") == 0)
finitecutflag = 1;
else
error->all(FLERR, "Illegal pair_style command");
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairTracker::coeff(int narg, char **arg)
{
if (narg > 2 && finitecutflag) error->all(FLERR, "Incorrect args for pair coefficients");
if (narg != 3 && !finitecutflag) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double cut_one = 0.0;
if (!finitecutflag) cut_one = utils::numeric(FLERR, arg[2], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
setflag[i][j] = 1;
cut[i][j] = cut_one;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairTracker::init_style()
{
int i;
// error and warning checks
if (!atom->radius_flag && finitecutflag)
error->all(FLERR, "Pair tracker requires atom attribute radius for finite cutoffs");
// need a history neigh list
int irequest = neighbor->request(this, instance_me);
if (finitecutflag) {
neighbor->requests[irequest]->size = 1;
neighbor->requests[irequest]->history = 1;
// history flag won't affect results, but match granular pairstyles
// so neighborlist can be copied to reduce overhead
}
// if history is stored and first init, create Fix to store history
// it replaces FixDummy, created in the constructor
// this is so its order in the fix list is preserved
if (fix_history == nullptr) {
char dnumstr[16];
sprintf(dnumstr, "%d", size_history);
char **fixarg = new char *[4];
fixarg[0] = (char *) "NEIGH_HISTORY_TRACK";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "NEIGH_HISTORY";
fixarg[3] = dnumstr;
modify->replace_fix("NEIGH_HISTORY_TRACK_DUMMY", 4, fixarg, 1);
delete[] fixarg;
int ifix = modify->find_fix("NEIGH_HISTORY_TRACK");
fix_history = (FixNeighHistory *) modify->fix[ifix];
fix_history->pair = this;
fix_history->use_bit_flag = 0;
}
if (finitecutflag) {
if (force->pair->beyond_contact)
error->all(FLERR,
"Pair tracker incompatible with granular pairstyles that extend beyond contact");
// check for FixPour and FixDeposit so can extract particle radii
int ipour;
for (ipour = 0; ipour < modify->nfix; ipour++)
if (strcmp(modify->fix[ipour]->style, "pour") == 0) break;
if (ipour == modify->nfix) ipour = -1;
int idep;
for (idep = 0; idep < modify->nfix; idep++)
if (strcmp(modify->fix[idep]->style, "deposit") == 0) break;
if (idep == modify->nfix) idep = -1;
// set maxrad_dynamic and maxrad_frozen for each type
// include future FixPour and FixDeposit particles as dynamic
int itype;
for (i = 1; i <= atom->ntypes; i++) {
onerad_dynamic[i] = onerad_frozen[i] = 0.0;
if (ipour >= 0) {
itype = i;
onerad_dynamic[i] = *((double *) modify->fix[ipour]->extract("radius", itype));
}
if (idep >= 0) {
itype = i;
onerad_dynamic[i] = *((double *) modify->fix[idep]->extract("radius", itype));
}
}
double *radius = atom->radius;
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (mask[i] & freeze_group_bit)
onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]], radius[i]);
else
onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]], radius[i]);
MPI_Allreduce(&onerad_dynamic[1], &maxrad_dynamic[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
MPI_Allreduce(&onerad_frozen[1], &maxrad_frozen[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
}
int ifix = modify->find_fix("NEIGH_HISTORY_TRACK");
if (ifix < 0) error->all(FLERR, "Could not find pair fix neigh history ID");
fix_history = (FixNeighHistory *) modify->fix[ifix];
ifix = modify->find_fix_by_style("pair/tracker");
if (ifix < 0) error->all(FLERR, "Cannot use pair tracker without fix pair/tracker");
fix_pair_tracker = (FixPairTracker *) modify->fix[ifix];
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairTracker::init_one(int i, int j)
{
if (!allocated) allocate();
// always mix prefactors geometrically
if (setflag[i][j] == 0) { cut[i][j] = mix_distance(cut[i][i], cut[j][j]); }
cut[j][i] = cut[i][j];
// if finite, cutoff = sum of max I,J radii for
// dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
double cutoff;
if (finitecutflag) {
cutoff = maxrad_dynamic[i] + maxrad_dynamic[j];
cutoff = MAX(cutoff, maxrad_frozen[i] + maxrad_dynamic[j]);
cutoff = MAX(cutoff, maxrad_dynamic[i] + maxrad_frozen[j]);
} else {
cutoff = cut[i][j];
}
return cutoff;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairTracker::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) { fwrite(&cut[i][j], sizeof(double), 1, fp); }
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairTracker::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) { utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); }
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairTracker::write_restart_settings(FILE *fp)
{
fwrite(&mix_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairTracker::read_restart_settings(FILE *fp)
{
if (comm->me == 0) { utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); }
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */
double PairTracker::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, double /*rsq*/,
double /*factor_coul*/, double /*factor_lj*/, double &/*fforce*/)
{
return 0.0;
}
/* ----------------------------------------------------------------------
transfer history during fix/neigh/history exchange
only needed if any history entries i-j are not just negative of j-i entries
------------------------------------------------------------------------- */
void PairTracker::transfer_history(double *source, double *target)
{
for (int i = 0; i < size_history; i++) target[i] = source[i];
}
/* ----------------------------------------------------------------------
self-interaction range of particle if finite particles
------------------------------------------------------------------------- */
double PairTracker::atom2cut(int i)
{
double cut = atom->radius[i] * 2;
return cut;
}
/* ----------------------------------------------------------------------
maximum interaction range for two finite particles
------------------------------------------------------------------------- */
double PairTracker::radii2cut(double r1, double r2)
{
double cut = r1 + r2;
return cut;
}

91
src/MISC/pair_tracker.h Normal file
View File

@ -0,0 +1,91 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(tracker,PairTracker);
// clang-format on
#else
#ifndef LMP_PAIR_TRACKER_H
#define LMP_PAIR_TRACKER_H
#include "pair.h"
namespace LAMMPS_NS {
class PairTracker : public Pair {
public:
PairTracker(class LAMMPS *);
virtual ~PairTracker();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
double atom2cut(int);
double radii2cut(double, double);
protected:
int sizeflag;
int history;
int size_history;
int neighprev;
double **cut;
double *onerad_dynamic, *onerad_frozen;
double *maxrad_dynamic, *maxrad_frozen;
int freeze_group_bit;
class FixDummy *fix_dummy;
class FixNeighHistory *fix_history;
class FixPairTracker *fix_pair_tracker;
void transfer_history(double *, double *);
void allocate();
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair tracker requires atom attribute radius for finite cutoffs
The atom style defined does not have these attributes.
E: Could not find pair fix neigh history ID
The associated fix neigh/history is missing
E: Cannot use pair tracker without fix pair/tracker
This pairstyle requires one to define a pair/tracker fix
*/

10
src/dump.cpp
View File

@ -84,6 +84,7 @@ Dump::Dump(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp)
unit_flag = 0;
unit_count = 0;
delay_flag = 0;
write_header_flag = 1;
maxfiles = -1;
numfiles = 0;
@ -378,7 +379,7 @@ void Dump::write()
if (multiproc)
MPI_Allreduce(&bnme,&nheader,1,MPI_LMP_BIGINT,MPI_SUM,clustercomm);
if (filewriter) write_header(nheader);
if (filewriter && write_header_flag) write_header(nheader);
// insure buf is sized for packing and communicating
// use nmax to insure filewriter proc can receive info from others
@ -932,6 +933,13 @@ void Dump::modify_params(int narg, char **arg)
else delay_flag = 0;
iarg += 2;
} else if (strcmp(arg[iarg],"header") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) write_header_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) write_header_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
int idump;

2
src/dump.h
View File

@ -79,6 +79,8 @@ class Dump : protected Pointers {
int unit_flag; // 1 if dump should contain unit information
int unit_count; // # of times the unit information was written
int delay_flag; // 1 if delay output until delaystep
int write_header_flag; // 1 if write header, 0 if not
bigint delaystep;
int refreshflag; // 1 if dump_modify refresh specified

11
src/fix_neigh_history.cpp
View File

@ -45,6 +45,7 @@ FixNeighHistory::FixNeighHistory(LAMMPS *lmp, int narg, char **arg) :
create_attribute = 1;
maxexchange_dynamic = 1;
use_bit_flag = 1;
newton_pair = force->newton_pair;
@ -624,13 +625,21 @@ void FixNeighHistory::post_neighbor()
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (use_bit_flag) {
rflag = sbmask(j) | pair->beyond_contact;
j &= NEIGHMASK;
jlist[jj] = j;
} else {
rflag = 1;
j &= NEIGHMASK;
}
// rflag = 1 if r < radsum in npair_size() method
// rflag = 1 if r < radsum in npair_size() method or if pair interactions extend further
// preserve neigh history info if tag[j] is in old-neigh partner list
// this test could be more geometrically precise for two sphere/line/tri
// if use_bit_flag is turned off, always record data since not all npair classes
// apply a mask for history (and they could use the bits for special bonds)
if (rflag) {
jtag = tag[j];

1
src/fix_neigh_history.h
View File

@ -28,6 +28,7 @@ class FixNeighHistory : public Fix {
public:
int nlocal_neigh; // nlocal at last time neigh list was built
int nall_neigh; // ditto for nlocal+nghost
int use_bit_flag; // flag whether this fix uses the extra bit in the nlist
int **firstflag; // ptr to each atom's neighbor flsg
double **firstvalue; // ptr to each atom's values
class Pair *pair; // ptr to pair style that uses neighbor history

4
src/pair_hybrid.cpp
View File

@ -597,6 +597,10 @@ void PairHybrid::init_style()
}
}
// check beyond contact (set during pair coeff) before init style
for (istyle = 0; istyle < nstyles; istyle++)
if (styles[istyle]->beyond_contact) beyond_contact = 1;
// each sub-style makes its neighbor list request(s)
for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style();