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

Merge pull request #4614 from akohlmey/collected-small-fixes

Browse Source 
Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer
2025-06-12 23:59:06 -04:00
committed by GitHub
parent 66684deb9f 349def92d2
commit 949bac6b4a
22 changed files with 80 additions and 2009 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_fix.rst
View File

@ -113,7 +113,6 @@ OPT.
* :doc:`mvv/tdpd <fix_mvv_dpd>`
* :doc:`neb <fix_neb>`
* :doc:`neb/spin <fix_neb_spin>`
* :doc:`neighbor/swap <fix_neighbor_swap>`
* :doc:`nonaffine/displacement <fix_nonaffine_displacement>`
* :doc:`nph (ko) <fix_nh>`
* :doc:`nph/asphere (o) <fix_nph_asphere>`

4
doc/src/Install_git.rst
View File

@ -21,8 +21,8 @@ You can follow the LAMMPS development on 4 different git branches:
* **develop** : this branch follows the ongoing development and is
updated with every merge commit of a pull request
* **release** : this branch is updated with every "feature release";
updates are always "fast-forward" merges from *develop*
* **release** : this branch is updated with every "feature release"
and updates are always "fast-forward" merges from *develop*
* **maintenance** : this branch collects back-ported bug fixes from the
*develop* branch to the *stable* branch. It is used to update the
*stable* branch for "stable update releases".

1
doc/src/fix.rst
View File

@ -292,7 +292,6 @@ accelerated styles exist.
* :doc:`mvv/tdpd <fix_mvv_dpd>` - constant temperature DPD using the modified velocity-Verlet algorithm
* :doc:`neb <fix_neb>` - nudged elastic band (NEB) spring forces
* :doc:`neb/spin <fix_neb_spin>` - nudged elastic band (NEB) spring forces for spins
* :doc:`neighbor/swap <fix_neighbor_swap>` - kinetic Monte Carlo (kMC) atom swapping
* :doc:`nonaffine/displacement <fix_nonaffine_displacement>` - calculate nonaffine displacement of atoms
* :doc:`nph <fix_nh>` - constant NPH time integration via Nose/Hoover
* :doc:`nph/asphere <fix_nph_asphere>` - NPH for aspherical particles

264
doc/src/fix_neighbor_swap.rst
View File

@ -1,264 +0,0 @@
.. index:: fix neighbor/swap
fix neighbor/swap command
=========================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID neighbor/swap N X seed T R0 voro-ID keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* neighbor/swap = style name of this fix command
* N = invoke this fix every N steps
* X = number of swaps to attempt every N steps
* seed = random # seed (positive integer)
* T = scaling temperature of the MC swaps (temperature units)
* R0 = scaling swap probability of the MC swaps (distance units)
* voro-ID = valid voronoi compute id (compute voronoi/atom)
* one or more keyword/value pairs may be appended to args
* keywords *types* and *diff* are mutually exclusive, but one must be specified
* keyword = *types* or *diff* or *ke* or *region* or *rates*
.. parsed-literal::
*types* values = two or more atom types (Integers in range [1,Ntypes] or type labels)
*diff* values = one atom type
*ke* value = *yes* or *no*
*yes* = kinetic energy is conserved after atom swaps
*no* = no conservation of kinetic energy after atom swaps
*region* value = region-ID
region-ID = ID of region to use as an exchange/move volume
*rates* values = V1 V2 . . . Vntypes values to conduct variable diffusion for different atom types (unitless)
Examples
""""""""
.. code-block:: LAMMPS
compute voroN all voronoi/atom neighbors yes
fix mc all neighbor/swap 10 160 15238 1000.0 3.0 voroN diff 2
fix myFix all neighbor/swap 100 1 12345 298.0 3.0 voroN region my_swap_region types 5 6
fix kmc all neighbor/swap 1 100 345 1.0 3.0 voroN diff 3 rates 3 1 6
Description
"""""""""""
.. versionadded:: TBD
This fix performs Monte-Carlo (MC) evaluations to enable kinetic
Monte Carlo (kMC)-type behavior during MD simulation by allowing
neighboring atoms to swap their positions. In contrast to the :doc:`fix
atom/swap <fix_atom_swap>` command which swaps pairs of atoms anywhere
in the simulation domain, the restriction of the MC swapping to
neighbors enables a hybrid MD/kMC-like simulation.
Neighboring atoms are defined by using a Voronoi tesselation performed
by the :doc:`compute voronoi/atom <compute_voronoi_atom>` command.
Two atoms are neighbors if their Voronoi cells share a common face
(3d) or edge (2d).
The selection of a swap neighbor is made using a distance-based
criterion for weighting the selection probability of each swap, in the
same manner as kMC selects a next event using relative probabilities.
The acceptance or rejection of each swap is determined via the
Metropolis criterion after evaluating the change in system energy due
to the swap.
A detailed explanation of the original implementation of this
algorithm can be found in :ref:`(Tavenner 2023) <TavennerMDkMC>`
where it was used to simulated accelerated diffusion in an MD context.
Simulating inherently kinetically-limited behaviors which rely on rare
events (such as atomic diffusion in a solid) is challenging for
traditional MD since its relatively short timescale will not naturally
sample many events. This fix addresses this challenge by allowing rare
neighbor hopping events to be sampled in a kMC-like fashion at a much
faster rate (set by the specified *N* and *X* parameters). This enables
the processes of atomic diffusion to be approximated during an MD
simulation, effectively decoupling the MD atomic vibrational timescale
and the atomic hopping (kMC event) timescale.
The algorithm implemented by this fix is as follows:
- The MD simulation is paused every *N* steps
- A Voronoi tesselation is performed for the current atom configuration.
- Then *X* atom swaps are attempted, one after the other.
- For each swap, an atom *I* is selected randomly from the list of
atom types specified by either the *types* or *diff* keywords.
- One of *I*'s Voronoi neighbors *J* is selected using the
distance-weighted probability for each neighbor detailed below.
- The *I,J* atom IDs are communicated to all processors so that a
global energy evaluation can be performed for the post-swap state
of the system.
- The swap is accepted or rejected based on the Metropolis criterion
using the energy change of the system and the specified temperature
*T*.
Here are a few comments on the computational cost of the swapping
algorithm.
1. The cost of a global energy evaluation is similar to that of an MD
timestep.
2. Similar to other MC algorithms in LAMMPS, improved parallel
efficiency is achieved with a smaller number of atoms per
processor than would typically be used in an standard MD
simulation. This is because the per-energy evaluation cost
increases relative to the balance of MD/MC steps as indicated by
1., but the communication cost remains relatively constant for a
given number of MD steps.
3. The MC portion of the simulation will run dramatically slower if
the pair style uses different cutoffs for different atom types (or
type pairs). This is because each atom swap then requires a
rebuild of the neighbor list to ensure the post-swap global energy
can be computed correctly.
Limitations are imposed on selection of *I,J* atom pairs to avoid
swapping of atoms which are outside of a reasonable cutoff (e.g. due to
a Voronoi tesselation near free surfaces) though the use of a
distance-weighted probability scaling.
----------
This section gives more details on other arguments and keywords.
The random number generator (RNG) used by all the processors for MC
operations is initialized with the specified *seed*.
The distance-based probability is weighted by the specified *R0* which
sets the radius :math:`r_0` in this formula
.. math::
p_{ij} = e^{(\frac{r_{ij}}{r_0})^2}
where :math:`p_{ij}` is the probability of selecting atom :math:`j` to
swap with atom :math:`i`. Typically, a value for *R0* around the
average nearest-neighbor spacing is appropriate. Since this is simply a
probability weighting, the swapping behavior is not very sensitive to
the exact value of *R0*.
The required *voro-ID* value is the compute-ID of a
:doc:`compute voronoi/atom <compute_voronoi_atom>` command like
this:
.. code-block:: LAMMPS
compute compute-ID group-ID voronoi/atom neighbors yes
It must return per-atom list of valid neighbor IDs as in the
:doc:`compute voronoi/atom <compute_voronoi_atom>` command.
The keyword *types* takes two or more atom types as its values. Only
atoms *I* of the first atom type will be selected. Only atoms *J* of the
remaining atom types will be considered as potential swap partners.
The keyword *diff* take a single atom type as its value. Only atoms
*I* of the that atom type will be selected. Atoms *J* of all
remaining atom types will be considered as potential swap partners.
This includes the atom type specified with the *diff* keyword to
account for self-diffusive hops between two atoms of the same type.
Note that the *neighbors yes* option must be enabled for use with this
fix. The group-ID should include all the atoms which this fix will
potentially select. I.e. the group-ID used in the voronoi compute should
include the same atoms as that indicated by the *types* keyword. If the
*diff* keyword is used, the group-ID should include atoms of all types
in the simulation.
The keyword *ke* takes *yes* (default) or *no* as its value. It two
atoms are swapped with different masses, then a value of *yes* will
rescale their respective velocities to conserve the kinetic energy of
the system. A value of *no* will perform no rescaling, so that
kinetic energy is not conserved. See the restriction on this keyword
below.
The *region* keyword takes a *region-ID* as its value. If specified,
then only atoms *I* and *J* within the geometric region will be
considered as swap partners. See the :doc:`region <region>` command
for details. This means the group-ID for the :doc:`compute
voronoi/atom <compute_voronoi_atom>` command also need only contain
atoms within the region.
The keyword *rates* can modify the swap rate based on the type of atom
*J*. Ntype values must be specified, where Ntype = the number of atom
types in the system. Each value is used to scale the probability
weighting given by the equation above. In the third example command
above, a simulation has 3 atoms types. Atom *I*s of type 1 are
eligible for swapping. Swaps may occur with atom *J*s of all 3 types.
Assuming all *J* atoms are equidistant from an atom *I*, *J* atoms of
type 1 will be 3x more likely to be selected as a swap partner than
atoms of type 2. And *J* atoms of type 3 will be 6.5x more likely to
be selected than atoms of type 2. If the *rates* keyword is not used,
all atom types will be treated with the same probability during selection
of swap attempts.
Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix writes the state of the fix to :doc:`binary restart files
<restart>`. This includes information about the random number generator
seed, the next timestep for MC exchanges, and the number of exchange
attempts and successes. See the :doc:`read_restart <read_restart>`
command for info on how to re-specify a fix in an input script that
reads a restart file, so that the operation of the fix continues in an
uninterrupted fashion.
None of the :doc:`fix_modify <fix_modify>` options are relevant to this
fix.
This fix computes a global vector of length 2, which can be accessed
by various :doc:`output commands <Howto_output>`. The vector values are
the following global cumulative quantities:
#. swap attempts
#. swap accepts
The vector values calculated by this fix are "intensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is part of the MC package. It is only enabled if LAMMPS was
built with that package. See the :doc:`Build package <Build_package>`
doc page for more info. Also this fix requires that the :ref:`VORONOI
package <PKG-VORONOI>` is installed, otherwise the fix will not be
compiled.
The :doc:`compute voronoi/atom <compute_voronoi_atom>` command
referenced by the required voro-ID must return neighboring atoms as
illustrated in the examples above.
If this fix is used with systems that do not have per-type masses
(e.g. atom style sphere), the *ke* keyword must be set to *off* since
the implemented algorithm will not be able to re-scale velocities
properly.
Related commands
""""""""""""""""
:doc:`fix nvt <fix_nh>`, :doc:`compute voronoi/atom <compute_voronoi_atom>`
:doc:`delete_atoms <delete_atoms>`, :doc:`fix gcmc <fix_gcmc>`,
:doc:`fix atom/swap <fix_atom_swap>`, :doc:`fix mol/swap <fix_mol_swap>`,
:doc:`fix sgcmc <fix_sgcmc>`
Default
"""""""
The option defaults are *ke* = yes and *rates* = 1 for all atom types.
----------
.. _TavennerMDkMC:
**(Tavenner 2023)** J Tavenner, M Mendelev, J Lawson, Computational
Materials Science, 218, 111929 (2023).

11
doc/utils/sphinx-config/conf.py.in
View File

@ -370,6 +370,17 @@ latex_elements = {
{%
\hypersetup{pageanchor=false}% avoid duplicate destination warnings
\begin{titlepage}%
\sffamily\Large
The LAMMPS developers are thinking about dropping the PDF format version of
the LAMMPS manual. This would allow us to focus on the HTML version, use
HTML-only features, and skip checking if the documentation source files,
especially the embedded mathematical expressions, are compatible with \LaTeX{} output.
Please let us know how you feel about this change by sending an email to
\texttt{developers@lammps.org} stating whether you agree or disagree with
removing support for the PDF format version of the manual and optionally
provide arguments for your preference.
\clearpage
\sffamily\bfseries
\begingroup % for PDF information dictionary
\def\endgraf{ }\def\and{\& }%

1
examples/PACKAGES/neighbor-swap/MoCoNiVFeAlCr_2nn.meam
View File

@ -1 +0,0 @@
../../../potentials/MoCoNiVFeAlCr_2nn.meam

46
examples/PACKAGES/neighbor-swap/in.KMC_pulse_center
View File

@ -1,46 +0,0 @@
# May 2025
# Test script for MD-KMC accelerated diffusion testing in LAMMPS
# Created by Jacob Tavenner, Baylor University
# Initiation -------------------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
# Atom Definition --------------------------------
lattice fcc 3.762
region whole block 0 1 0 1 0 1
create_box 2 whole
create_atoms 1 region whole
replicate 6 16 6
region puck block INF INF 7 9 INF INF
set region puck type 2
# Force Fields -----------------------------------
pair_style meam
pair_coeff * * library_2nn.meam Mo Co Ni V Fe Al Cr MoCoNiVFeAlCr_2nn.meam Ni Cr
# Settings ---------------------------------------
timestep 0.002
thermo 100
# Computations -----------------------------------
compute voroN all voronoi/atom neighbors yes
run 0
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe
# Execution --------------------------------------
velocity all create 2400 908124 loop geom
fix temp all npt temp 1000 1000 1000 aniso 0 0 1
fix mc all neighbor/swap 50 12 1340723 1000 3 voroN diff 2
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe f_mc[*]
#dump dump2 all custom 5000 dump.edge-3_Ni-Cr.* id type x y z c_eng c_csym
run 1000
#write_data pulse_center.data

47
examples/PACKAGES/neighbor-swap/in.KMC_pulse_edge
View File

@ -1,47 +0,0 @@
# May 2025
# Test script for MD-KMC accelerated diffusion testing in LAMMPS
# Created by Jacob Tavenner, Baylor University
# Initiation -------------------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
# Atom Definition --------------------------------
lattice fcc 3.762
region whole block 0 1 0 1 0 1
create_box 2 whole
create_atoms 1 region whole
replicate 6 16 6
region puck block INF INF INF 2 INF INF
set region puck type 2
# Force Fields -----------------------------------
pair_style meam
pair_coeff * * library_2nn.meam Mo Co Ni V Fe Al Cr MoCoNiVFeAlCr_2nn.meam Ni Cr
# Settings ---------------------------------------
timestep 0.002
thermo 100
# Computations -----------------------------------
compute voroN all voronoi/atom neighbors yes
run 0
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe
# Execution --------------------------------------
velocity all create 2400 908124 loop geom
fix temp all npt temp 1000 1000 1000 aniso 0 0 1
fix mc all neighbor/swap 50 12 1340723 1000 3 voroN diff 2
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe f_mc[*]
#dump dump2 all custom 5000 dump.edge-3_Ni-Cr.* id type x y z c_eng c_csym
run 1000
#write_data pulse_end.data

1
examples/PACKAGES/neighbor-swap/library_2nn.meam
View File

@ -1 +0,0 @@
../../../potentials/library_2nn.meam

154
examples/PACKAGES/neighbor-swap/log.22May22.KMC_pulse_center.g++.1
View File

@ -1,154 +0,0 @@
LAMMPS (2 Apr 2025 - Development - patch_2Apr2025-384-g88bc7dc720-modified)
using 1 OpenMP thread(s) per MPI task
# May 2025
# Test script for MD-KMC accelerated diffusion testing in LAMMPS
# Created by Jacob Tavenner, Baylor University
# Initiation -------------------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
# Atom Definition --------------------------------
lattice fcc 3.762
Lattice spacing in x,y,z = 3.762 3.762 3.762
region whole block 0 1 0 1 0 1
create_box 2 whole
Created orthogonal box = (0 0 0) to (3.762 3.762 3.762)
1 by 1 by 1 MPI processor grid
create_atoms 1 region whole
Created 4 atoms
using lattice units in orthogonal box = (0 0 0) to (3.762 3.762 3.762)
create_atoms CPU = 0.000 seconds
replicate 6 16 6
Replication is creating a 6x16x6 = 576 times larger system...
orthogonal box = (0 0 0) to (22.572 60.192 22.572)
1 by 1 by 1 MPI processor grid
2304 atoms
replicate CPU = 0.000 seconds
region puck block INF INF 7 9 INF INF
set region puck type 2
Setting atom values ...
360 settings made for type
# Force Fields -----------------------------------
pair_style meam
pair_coeff * * library_2nn.meam Mo Co Ni V Fe Al Cr MoCoNiVFeAlCr_2nn.meam Ni Cr
Reading MEAM library file library_2nn.meam with DATE: 2024-08-08
Reading MEAM potential file MoCoNiVFeAlCr_2nn.meam with DATE: 2024-08-08
# Settings ---------------------------------------
timestep 0.002
thermo 100
# Computations -----------------------------------
compute voroN all voronoi/atom neighbors yes
run 0
WARNING: No fixes with time integration, atoms won't move
For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60)
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.8
ghost atom cutoff = 6.8
binsize = 3.4, bins = 7 18 7
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 13.32 | 13.32 | 13.32 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -9674.3728 0 -9674.3728 -212400.94
Loop time of 1.202e-06 on 1 procs for 0 steps with 2304 atoms
0.0% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0 | 0 | 0 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0 | 0 | 0 | 0.0 | 0.00
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 1.202e-06 | | |100.00
Nlocal: 2304 ave 2304 max 2304 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4735 ave 4735 max 4735 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 99072 ave 99072 max 99072 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 198144 ave 198144 max 198144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 198144
Ave neighs/atom = 86
Neighbor list builds = 0
Dangerous builds = 0
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe
# Execution --------------------------------------
velocity all create 2400 908124
fix temp all npt temp 1000 1000 1000 aniso 0 0 1
fix mc all neighbor/swap 50 12 1340723 1000 3 voroN diff 2
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe f_mc[*]
#dump dump2 all custom 5000 dump.edge-3_Ni-Cr.* id type x y z c_eng c_csym
run 1000
Per MPI rank memory allocation (min/avg/max) = 13.32 | 13.32 | 13.32 Mbytes
Step Temp Press Pxx Pyy Pzz Lx Ly Lz Volume PotEng f_mc[1] f_mc[2]
0 2400 -187517.52 -187403.07 -187750.14 -187399.35 22.572 60.192 22.572 30667.534 -9674.3728 0 0
100 1664.9956 14000 14280.682 15095.077 12624.241 21.635315 57.726568 21.64791 27036.778 -9592.8978 24 22
200 1560.0093 -5452.2434 -5749.5816 -2957.4228 -7649.7258 21.734212 58.085959 21.724853 27426.596 -9562.8822 48 45
300 1586.4553 2030.9253 2776.4677 775.50538 2540.803 21.678654 58.101753 21.654423 27275.215 -9571.1308 72 66
400 1603.6896 -223.16773 156.17673 -478.47929 -347.20061 21.701021 58.098904 21.657752 27306.213 -9576.4456 96 90
500 1618.236 -925.51874 -1640.9078 451.6228 -1587.2713 21.718334 58.042685 21.666081 27312.054 -9581.2045 120 110
600 1581.9995 290.10126 1359.1314 1407.5434 -1896.371 21.679813 58.086147 21.692118 27316.815 -9570.4803 144 132
700 1568.3261 1387.3472 938.81523 2159.3686 1063.8577 21.685928 58.075626 21.67273 27295.153 -9566.2914 168 155
800 1607.1531 46.792964 -453.90265 -1533.3908 2127.6723 21.685188 58.202356 21.628338 27297.753 -9577.7848 192 177
900 1573.4747 -84.225488 548.90935 -1356.7479 555.16208 21.69634 58.150052 21.651847 27316.908 -9567.7039 216 196
1000 1609.2136 1215.0833 764.08936 3301.0811 -419.92053 21.683731 58.000401 21.68726 27275.31 -9578.2843 240 219
Loop time of 31.6263 on 1 procs for 1000 steps with 2304 atoms
Performance: 5.464 ns/day, 4.393 hours/ns, 31.619 timesteps/s, 72.851 katom-step/s
99.2% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 28.487 | 28.487 | 28.487 | 0.0 | 90.07
Neigh | 0.22789 | 0.22789 | 0.22789 | 0.0 | 0.72
Comm | 0.010808 | 0.010808 | 0.010808 | 0.0 | 0.03
Output | 0.00033526 | 0.00033526 | 0.00033526 | 0.0 | 0.00
Modify | 2.8963 | 2.8963 | 2.8963 | 0.0 | 9.16
Other | | 0.003905 | | | 0.01
Nlocal: 2304 ave 2304 max 2304 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4750 ave 4750 max 4750 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 130023 ave 130023 max 130023 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 260046 ave 260046 max 260046 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 260046
Ave neighs/atom = 112.86719
Neighbor list builds = 65
Dangerous builds = 0
#write_data pulse_center.data
Total wall time: 0:00:31

154
examples/PACKAGES/neighbor-swap/log.22May22.KMC_pulse_center.g++.4
View File

@ -1,154 +0,0 @@
LAMMPS (2 Apr 2025 - Development - patch_2Apr2025-384-g88bc7dc720-modified)
using 1 OpenMP thread(s) per MPI task
# May 2025
# Test script for MD-KMC accelerated diffusion testing in LAMMPS
# Created by Jacob Tavenner, Baylor University
# Initiation -------------------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
# Atom Definition --------------------------------
lattice fcc 3.762
Lattice spacing in x,y,z = 3.762 3.762 3.762
region whole block 0 1 0 1 0 1
create_box 2 whole
Created orthogonal box = (0 0 0) to (3.762 3.762 3.762)
1 by 2 by 2 MPI processor grid
create_atoms 1 region whole
Created 4 atoms
using lattice units in orthogonal box = (0 0 0) to (3.762 3.762 3.762)
create_atoms CPU = 0.000 seconds
replicate 6 16 6
Replication is creating a 6x16x6 = 576 times larger system...
orthogonal box = (0 0 0) to (22.572 60.192 22.572)
1 by 4 by 1 MPI processor grid
2304 atoms
replicate CPU = 0.000 seconds
region puck block INF INF 7 9 INF INF
set region puck type 2
Setting atom values ...
360 settings made for type
# Force Fields -----------------------------------
pair_style meam
pair_coeff * * library_2nn.meam Mo Co Ni V Fe Al Cr MoCoNiVFeAlCr_2nn.meam Ni Cr
Reading MEAM library file library_2nn.meam with DATE: 2024-08-08
Reading MEAM potential file MoCoNiVFeAlCr_2nn.meam with DATE: 2024-08-08
# Settings ---------------------------------------
timestep 0.002
thermo 100
# Computations -----------------------------------
compute voroN all voronoi/atom neighbors yes
run 0
WARNING: No fixes with time integration, atoms won't move
For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60)
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.8
ghost atom cutoff = 6.8
binsize = 3.4, bins = 7 18 7
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 9.636 | 9.636 | 9.636 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -9674.3728 0 -9674.3728 -212400.94
Loop time of 1.422e-06 on 4 procs for 0 steps with 2304 atoms
35.2% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0 | 0 | 0 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0 | 0 | 0 | 0.0 | 0.00
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 1.422e-06 | | |100.00
Nlocal: 576 ave 576 max 576 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 2131 ave 2131 max 2131 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 24768 ave 24768 max 24768 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 49536 ave 49536 max 49536 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Total # of neighbors = 198144
Ave neighs/atom = 86
Neighbor list builds = 0
Dangerous builds = 0
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe
# Execution --------------------------------------
velocity all create 2400 908124
fix temp all npt temp 1000 1000 1000 aniso 0 0 1
fix mc all neighbor/swap 50 12 1340723 1000 3 voroN diff 2
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe f_mc[*]
#dump dump2 all custom 5000 dump.edge-3_Ni-Cr.* id type x y z c_eng c_csym
run 1000
Per MPI rank memory allocation (min/avg/max) = 9.636 | 9.636 | 9.636 Mbytes
Step Temp Press Pxx Pyy Pzz Lx Ly Lz Volume PotEng f_mc[1] f_mc[2]
0 2400 -187517.52 -187403.09 -187750.05 -187399.42 22.572 60.192 22.572 30667.534 -9674.3728 0 0
100 1668.8754 13300.763 12419.304 15568.772 11914.212 21.636248 57.724775 21.647685 27036.823 -9594.7526 24 23
200 1584.9699 -5686.0414 -4741.8496 -5914.7681 -6401.5064 21.729384 58.060532 21.730736 27415.923 -9571.0639 48 46
300 1582.0473 2806.2983 3413.4122 2716.0124 2289.4702 21.6679 58.033587 21.694744 27280.402 -9570.5549 72 69
400 1582.5825 845.29268 -849.61221 2123.5339 1261.9563 21.676298 58.14253 21.656418 27293.905 -9570.7948 96 93
500 1591.7285 -501.17955 1151.9743 -1719.3712 -936.14174 21.696367 58.157211 21.648308 27315.839 -9573.5089 120 116
600 1610.708 -821.74669 -1002.4957 291.88502 -1754.6294 21.730338 58.008213 21.661226 27304.8 -9579.5573 144 138
700 1598.5176 -590.00633 -1844.42 408.97706 -334.57602 21.712908 57.96131 21.698129 27307.281 -9575.8973 168 162
800 1584.3478 330.16711 666.88818 74.698331 248.91482 21.650908 58.045055 21.719838 27295.933 -9571.9268 192 186
900 1557.9946 1471.1207 2124.6512 1526.9937 761.71731 21.645578 58.156083 21.681637 27293.323 -9564.4385 216 207
1000 1582.5312 379.57005 -602.96446 2696.737 -955.06238 21.655418 58.231248 21.649581 27300.598 -9571.9879 240 227
Loop time of 9.1632 on 4 procs for 1000 steps with 2304 atoms
Performance: 18.858 ns/day, 1.273 hours/ns, 109.132 timesteps/s, 251.440 katom-step/s
98.5% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.867 | 7.9923 | 8.1311 | 4.3 | 87.22
Neigh | 0.054997 | 0.057518 | 0.060145 | 1.0 | 0.63
Comm | 0.017529 | 0.14801 | 0.27408 | 29.5 | 1.62
Output | 0.00015963 | 0.00017216 | 0.00020869 | 0.0 | 0.00
Modify | 0.95227 | 0.96325 | 0.9917 | 1.7 | 10.51
Other | | 0.001983 | | | 0.02
Nlocal: 576 ave 609 max 540 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 2161.5 ave 2173 max 2151 min
Histogram: 1 0 1 0 0 0 1 0 0 1
Neighs: 32450.2 ave 35422 max 29271 min
Histogram: 2 0 0 0 0 0 0 0 0 2
FullNghs: 64900.5 ave 70800 max 58684 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 259602
Ave neighs/atom = 112.67448
Neighbor list builds = 62
Dangerous builds = 0
#write_data pulse_center.data
Total wall time: 0:00:09

155
examples/PACKAGES/neighbor-swap/log.22May22.KMC_pulse_edge.g++.1
View File

@ -1,155 +0,0 @@
LAMMPS (2 Apr 2025 - Development - patch_2Apr2025-384-g88bc7dc720-modified)
using 1 OpenMP thread(s) per MPI task
# May 2025
# Test script for MD-KMC accelerated diffusion testing in LAMMPS
# Created by Jacob Tavenner, Baylor University
# Initiation -------------------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
# Atom Definition --------------------------------
lattice fcc 3.762
Lattice spacing in x,y,z = 3.762 3.762 3.762
region whole block 0 1 0 1 0 1
create_box 2 whole
Created orthogonal box = (0 0 0) to (3.762 3.762 3.762)
1 by 1 by 1 MPI processor grid
create_atoms 1 region whole
Created 4 atoms
using lattice units in orthogonal box = (0 0 0) to (3.762 3.762 3.762)
create_atoms CPU = 0.000 seconds
replicate 6 16 6
Replication is creating a 6x16x6 = 576 times larger system...
orthogonal box = (0 0 0) to (22.572 60.192 22.572)
1 by 1 by 1 MPI processor grid
2304 atoms
replicate CPU = 0.000 seconds
region puck block INF INF INF 2 INF INF
set region puck type 2
Setting atom values ...
360 settings made for type
# Force Fields -----------------------------------
pair_style meam
pair_coeff * * library_2nn.meam Mo Co Ni V Fe Al Cr MoCoNiVFeAlCr_2nn.meam Ni Cr
Reading MEAM library file library_2nn.meam with DATE: 2024-08-08
Reading MEAM potential file MoCoNiVFeAlCr_2nn.meam with DATE: 2024-08-08
# Settings ---------------------------------------
timestep 0.002
thermo 100
# Computations -----------------------------------
compute voroN all voronoi/atom neighbors yes
run 0
WARNING: No fixes with time integration, atoms won't move
For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60)
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.8
ghost atom cutoff = 6.8
binsize = 3.4, bins = 7 18 7
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 13.32 | 13.32 | 13.32 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -9674.3728 0 -9674.3728 -212400.94
Loop time of 1.232e-06 on 1 procs for 0 steps with 2304 atoms
81.2% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0 | 0 | 0 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0 | 0 | 0 | 0.0 | 0.00
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 1.232e-06 | | |100.00
Nlocal: 2304 ave 2304 max 2304 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4735 ave 4735 max 4735 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 99072 ave 99072 max 99072 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 198144 ave 198144 max 198144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 198144
Ave neighs/atom = 86
Neighbor list builds = 0
Dangerous builds = 0
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe
# Execution --------------------------------------
velocity all create 2400 908124 loop geom
fix temp all npt temp 1000 1000 1000 aniso 0 0 1
fix mc all neighbor/swap 50 12 1340723 1000 3 voroN diff 2
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe f_mc[*]
#dump dump2 all custom 5000 dump.edge-3_Ni-Cr.* id type x y z c_eng c_csym
run 1000
Per MPI rank memory allocation (min/avg/max) = 13.32 | 13.32 | 13.32 Mbytes
Step Temp Press Pxx Pyy Pzz Lx Ly Lz Volume PotEng f_mc[1] f_mc[2]
0 2400 -187517.52 -187464.47 -188202.62 -186885.48 22.572 60.192 22.572 30667.534 -9674.3728 0 0
100 1665.6154 14281.316 14426.547 14555.867 13861.534 21.637238 57.719793 21.637281 27022.733 -9594.4303 24 24
200 1603.3309 -7325.7341 -8878.1524 -5333.0485 -7766.0015 21.710246 58.122827 21.725933 27415.106 -9577.4545 48 48
300 1603.2974 207.19165 1983.4565 -1841.9518 480.07024 21.678227 58.079126 21.674033 27288.745 -9577.6391 72 69
400 1600.1515 810.95054 1087.969 802.04946 542.83316 21.683731 58.045848 21.678505 27285.662 -9576.6508 96 92
500 1629.8313 -2808.1005 -3197.9357 310.89931 -5537.265 21.683924 58.090375 21.697076 27330.229 -9585.5435 120 113
600 1598.8232 -67.845623 -1573.0718 -1526.7607 2896.2957 21.70213 58.12191 21.653853 27313.504 -9576.4147 144 137
700 1607.2185 154.66718 -1777.2469 2566.4705 -325.22208 21.712408 57.971553 21.678708 27287.033 -9579.1772 168 158
800 1582.559 -891.23631 -632.46037 -636.88203 -1404.3665 21.671936 58.127004 21.678224 27308.594 -9571.6663 192 180
900 1586.7172 -617.17083 -2495.5378 -2302.8766 2946.9018 21.658489 58.181921 21.668968 27305.771 -9572.9641 216 204
1000 1607.563 -389.8113 810.4908 298.84287 -2278.7676 21.624573 58.076745 21.724272 27283.183 -9579.5034 240 227
Loop time of 31.7733 on 1 procs for 1000 steps with 2304 atoms
Performance: 5.439 ns/day, 4.413 hours/ns, 31.473 timesteps/s, 72.514 katom-step/s
99.2% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 28.604 | 28.604 | 28.604 | 0.0 | 90.02
Neigh | 0.21293 | 0.21293 | 0.21293 | 0.0 | 0.67
Comm | 0.010645 | 0.010645 | 0.010645 | 0.0 | 0.03
Output | 0.00033194 | 0.00033194 | 0.00033194 | 0.0 | 0.00
Modify | 2.9411 | 2.9411 | 2.9411 | 0.0 | 9.26
Other | | 0.00448 | | | 0.01
Nlocal: 2304 ave 2304 max 2304 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4748 ave 4748 max 4748 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 130301 ave 130301 max 130301 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 260602 ave 260602 max 260602 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 260602
Ave neighs/atom = 113.10851
Neighbor list builds = 62
Dangerous builds = 0
#write_data pulse_end.data
Total wall time: 0:00:31

155
examples/PACKAGES/neighbor-swap/log.22May22.KMC_pulse_edge.g++.4
View File

@ -1,155 +0,0 @@
LAMMPS (2 Apr 2025 - Development - patch_2Apr2025-384-g88bc7dc720-modified)
using 1 OpenMP thread(s) per MPI task
# May 2025
# Test script for MD-KMC accelerated diffusion testing in LAMMPS
# Created by Jacob Tavenner, Baylor University
# Initiation -------------------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
# Atom Definition --------------------------------
lattice fcc 3.762
Lattice spacing in x,y,z = 3.762 3.762 3.762
region whole block 0 1 0 1 0 1
create_box 2 whole
Created orthogonal box = (0 0 0) to (3.762 3.762 3.762)
1 by 2 by 2 MPI processor grid
create_atoms 1 region whole
Created 4 atoms
using lattice units in orthogonal box = (0 0 0) to (3.762 3.762 3.762)
create_atoms CPU = 0.000 seconds
replicate 6 16 6
Replication is creating a 6x16x6 = 576 times larger system...
orthogonal box = (0 0 0) to (22.572 60.192 22.572)
1 by 4 by 1 MPI processor grid
2304 atoms
replicate CPU = 0.000 seconds
region puck block INF INF INF 2 INF INF
set region puck type 2
Setting atom values ...
360 settings made for type
# Force Fields -----------------------------------
pair_style meam
pair_coeff * * library_2nn.meam Mo Co Ni V Fe Al Cr MoCoNiVFeAlCr_2nn.meam Ni Cr
Reading MEAM library file library_2nn.meam with DATE: 2024-08-08
Reading MEAM potential file MoCoNiVFeAlCr_2nn.meam with DATE: 2024-08-08
# Settings ---------------------------------------
timestep 0.002
thermo 100
# Computations -----------------------------------
compute voroN all voronoi/atom neighbors yes
run 0
WARNING: No fixes with time integration, atoms won't move
For more information see https://docs.lammps.org/err0028 (src/verlet.cpp:60)
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.8
ghost atom cutoff = 6.8
binsize = 3.4, bins = 7 18 7
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 9.636 | 9.636 | 9.636 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -9674.3728 0 -9674.3728 -212400.94
Loop time of 1.53e-06 on 4 procs for 0 steps with 2304 atoms
65.4% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0 | 0 | 0 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0 | 0 | 0 | 0.0 | 0.00
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 1.53e-06 | | |100.00
Nlocal: 576 ave 576 max 576 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 2131 ave 2131 max 2131 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 24768 ave 24768 max 24768 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 49536 ave 49536 max 49536 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Total # of neighbors = 198144
Ave neighs/atom = 86
Neighbor list builds = 0
Dangerous builds = 0
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe
# Execution --------------------------------------
velocity all create 2400 908124 loop geom
fix temp all npt temp 1000 1000 1000 aniso 0 0 1
fix mc all neighbor/swap 50 12 1340723 1000 3 voroN diff 2
thermo_style custom step temp press pxx pyy pzz lx ly lz vol pe f_mc[*]
#dump dump2 all custom 5000 dump.edge-3_Ni-Cr.* id type x y z c_eng c_csym
run 1000
Per MPI rank memory allocation (min/avg/max) = 9.636 | 9.636 | 9.636 Mbytes
Step Temp Press Pxx Pyy Pzz Lx Ly Lz Volume PotEng f_mc[1] f_mc[2]
0 2400 -187517.52 -187464.47 -188202.62 -186885.48 22.572 60.192 22.572 30667.534 -9674.3728 0 0
100 1665.569 14271.813 14638.855 14316.569 13860.016 21.63675 57.721065 21.637799 27023.366 -9594.291 24 24
200 1598.6479 -6990.8349 -8574.1986 -5033.6147 -7364.6916 21.708963 58.123129 21.724821 27412.223 -9575.7322 48 47
300 1604.388 456.43285 1926.408 -1214.1721 657.0626 21.673369 58.090421 21.671716 27285.018 -9577.698 72 70
400 1601.1591 1303.6721 703.88473 1137.6607 2069.471 21.684004 58.049595 21.671161 27278.522 -9576.4811 96 94
500 1623.6044 -2243.2478 -2084.532 320.87709 -4966.0885 21.686171 58.097101 21.695911 27334.758 -9583.1878 120 118
600 1587.2041 421.60034 190.88741 -328.76599 1402.6796 21.712439 58.086039 21.655927 27312.229 -9572.559 144 141
700 1591.2923 32.327829 -2893.2353 1839.7574 1150.4614 21.719102 57.999862 21.666164 27292.974 -9573.9009 168 165
800 1580.8587 -105.51079 654.26389 -160.04168 -810.75457 21.670225 58.109245 21.684683 27306.229 -9570.6482 192 186
900 1570.7648 1290.088 1252.3689 255.62548 2362.2695 21.68101 58.100507 21.658755 27283.051 -9567.9864 216 209
1000 1598.1483 -125.35291 -3626.5479 3404.789 -154.29983 21.720146 57.952942 21.686111 27297.313 -9576.2975 240 231
Loop time of 9.17241 on 4 procs for 1000 steps with 2304 atoms
Performance: 18.839 ns/day, 1.274 hours/ns, 109.023 timesteps/s, 251.188 katom-step/s
98.1% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.7477 | 8.0143 | 8.1344 | 5.5 | 87.37
Neigh | 0.050543 | 0.056882 | 0.05986 | 1.6 | 0.62
Comm | 0.069784 | 0.16898 | 0.40996 | 34.2 | 1.84
Output | 0.00015612 | 0.0001707 | 0.00021249 | 0.0 | 0.00
Modify | 0.90628 | 0.93003 | 0.96157 | 2.2 | 10.14
Other | | 0.002053 | | | 0.02
Nlocal: 576 ave 614 max 505 min
Histogram: 1 0 0 0 0 0 1 0 0 2
Nghost: 2165.75 ave 2204 max 2132 min
Histogram: 1 0 0 0 2 0 0 0 0 1
Neighs: 32430.8 ave 35552 max 26564 min
Histogram: 1 0 0 0 0 0 1 0 0 2
FullNghs: 64861.5 ave 71111 max 53164 min
Histogram: 1 0 0 0 0 0 1 0 0 2
Total # of neighbors = 259446
Ave neighs/atom = 112.60677
Neighbor list builds = 62
Dangerous builds = 0
#write_data pulse_end.data
Total wall time: 0:00:09

14
src/KOKKOS/pair_mliap_kokkos.cpp
View File

@ -85,14 +85,24 @@ void PairMLIAPKokkos<DeviceType>::compute(int eflag, int vflag)
error->all(FLERR, "Incompatible model and descriptor element count");
ev_init(eflag, vflag, 0);
if (eflag_atom && (int)k_eatom.h_view.extent(0) < maxeatom) {
if (eflag_atom) {
if ((int)k_eatom.h_view.extent(0) < maxeatom) {
memoryKK->destroy_kokkos(k_eatom,eatom);
memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
} else {
Kokkos::deep_copy(k_eatom.d_view,0);
k_eatom.modify<DeviceType>();
}
}
if (vflag_atom && (int)k_vatom.h_view.extent(0) < maxeatom) {
if (vflag_atom) {
if ((int)k_vatom.h_view.extent(0) < maxeatom) {
memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxeatom,6,"pair:eatom");
} else {
Kokkos::deep_copy(k_vatom.d_view,0);
k_vatom.modify<DeviceType>();
}
}
data->generate_neighdata(list, eflag, vflag);

921
src/MC/fix_neighbor_swap.cpp
View File

@ -1,921 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Jacob Tavenner
------------------------------------------------------------------------- */
#include "fix_neighbor_swap.h"
#include "angle.h"
#include "atom.h"
#include "bond.h"
#include "citeme.h"
#include "comm.h"
#include "compute.h"
#include "compute_voronoi_atom.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "group.h"
#include "improper.h"
#include "kspace.h"
#include "math_extra.h"
#include "math_special.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include "random_park.h"
#include "region.h"
#include "update.h"
#include <cfloat>
#include <cmath>
#include <cstring>
#include <unordered_set>
using namespace LAMMPS_NS;
using namespace FixConst;
using MathExtra::distsq3;
using MathSpecial::square;
static const char cite_fix_neighbor_swap[] =
"fix neighbor/swap command: doi:10.1016/j.commatsci.2022.111929\n\n"
"@Article{Tavenner2023111929,\n"
" author = {Jacob P. Tavenner and Mikhail I. Mendelev and John W. Lawson},\n"
" title = {Molecular dynamics based kinetic Monte Carlo simulation for accelerated "
"diffusion},\n"
" journal = {Computational Materials Science},\n"
" year = {2023},\n"
" volume = {218},\n"
" pages = {111929}\n"
" url = {https://dx.doi.org/10.1016/j.commatsci.2022.111929}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), region(nullptr), idregion(nullptr), type_list(nullptr), rate_list(nullptr),
qtype(nullptr), sqrt_mass_ratio(nullptr), voro_neighbor_list(nullptr),
local_swap_iatom_list(nullptr), local_swap_neighbor_list(nullptr),
local_swap_type_list(nullptr), local_swap_probability(nullptr), random_equal(nullptr),
id_voro(nullptr), c_voro(nullptr), c_pe(nullptr)
{
if (narg < 10) utils::missing_cmd_args(FLERR, "fix neighbor/swap", error);
dynamic_group_allow = 1;
vector_flag = 1;
size_vector = 2;
global_freq = 1;
extvector = 0;
restart_global = 1;
time_depend = 1;
ke_flag = 1;
diff_flag = 0;
rates_flag = 0;
nswaptypes = 0;
if (lmp->citeme) lmp->citeme->add(cite_fix_neighbor_swap);
// required args
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
ncycles = utils::inumeric(FLERR, arg[4], false, lmp);
seed = utils::inumeric(FLERR, arg[5], false, lmp);
double temperature = utils::numeric(FLERR, arg[6], false, lmp);
double r_0 = utils::inumeric(FLERR, arg[7], false, lmp);
if (nevery <= 0)
error->all(FLERR, 3, "Illegal fix neighbor/swap command nevery value: {}", nevery);
if (ncycles < 0)
error->all(FLERR, 4, "Illegal fix neighbor/swap command ncycles value: {}", ncycles);
if (seed <= 0) error->all(FLERR, 5, "Illegal fix neighbor/swap command seed value: {}", seed);
if (temperature <= 0.0)
error->all(FLERR, 6, "Illegal fix neighbor/swap command temperature value: {}", temperature);
if (r_0 <= 0.0) error->all(FLERR, 7, "Illegal fix neighbor/swap command R0 value: {}", r_0);
// Voro compute check
id_voro = utils::strdup(arg[8]);
c_voro = modify->get_compute_by_id(id_voro);
if (!c_voro) error->all(FLERR, 8, "Could not find compute voronoi ID {}", id_voro);
if (c_voro->local_flag == 0)
error->all(FLERR, 8, "Voronoi compute {} does not compute local info", id_voro);
if (c_voro->size_local_cols != 3)
error->all(FLERR, 8, "Voronoi compute {} does not compute i, j, sizes as expected", id_voro);
beta = 1.0 / (force->boltz * temperature);
inv_r_0 = 1.0 / r_0;
memory->create(type_list, atom->ntypes, "neighbor/swap:type_list");
memory->create(rate_list, atom->ntypes, "neighbor/swap:rate_list");
// read options from end of input line
options(narg - 9, &arg[9]);
// random number generator, same for all procs
random_equal = new RanPark(lmp, seed);
// set up reneighboring
force_reneighbor = 1;
next_reneighbor = update->ntimestep + 1;
// zero out counters
nswap_attempts = 0.0;
nswap_successes = 0.0;
atom_swap_nmax = 0;
// set comm size needed by this Fix
if (atom->q_flag)
comm_forward = 2;
else
comm_forward = 1;
}
/* ---------------------------------------------------------------------- */
FixNeighborSwap::~FixNeighborSwap()
{
memory->destroy(type_list);
memory->destroy(rate_list);
memory->destroy(qtype);
memory->destroy(sqrt_mass_ratio);
memory->destroy(local_swap_iatom_list);
memory->destroy(local_swap_neighbor_list);
memory->destroy(local_swap_probability);
memory->destroy(local_swap_type_list);
delete[] idregion;
delete[] id_voro;
delete random_equal;
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
static const std::unordered_set<std::string> known_keywords = {"region", "ke", "types", "diff",
"rates"};
static bool is_keyword(const std::string &arg)
{
return known_keywords.find(arg) != known_keywords.end();
}
void FixNeighborSwap::options(int narg, char **arg)
{
if (narg < 0) utils::missing_cmd_args(FLERR, "fix neighbor/swap", error);
int ioffset = 9; // first 9 arguments are fixed and handled in constructor
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg], "region") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix neighbor/swap region", error);
delete[] idregion;
idregion = utils::strdup(arg[iarg + 1]);
region = domain->get_region_by_id(idregion);
if (!region)
error->all(FLERR, iarg + 1 + ioffset, "Region ID {} for fix neighbor/swap does not exist",
idregion);
iarg += 2;
} else if (strcmp(arg[iarg], "ke") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix neighbor/swap ke", error);
ke_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "types") == 0) {
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix neighbor/swap types", error);
if (diff_flag)
error->all(FLERR, iarg + ioffset, "Cannot use 'diff' and 'types' keywords together");
iarg++;
nswaptypes = 0;
while (iarg < narg) {
if (is_keyword(arg[iarg])) break;
if (nswaptypes >= atom->ntypes)
error->all(FLERR, iarg + ioffset, "Too many arguments to fix neighbor/swap types");
type_list[nswaptypes] = utils::expand_type_int(FLERR, arg[iarg], Atom::ATOM, lmp);
nswaptypes++;
iarg++;
}
} else if (strcmp(arg[iarg], "diff") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix neighbor/swap diff", error);
if (diff_flag) error->all(FLERR, iarg + ioffset, "Cannot use 'diff' keyword multiple times");
if (nswaptypes != 0)
error->all(FLERR, iarg + ioffset, "Cannot use 'diff' and 'types' keywords together");
type_list[nswaptypes] = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
diff_flag = 1;
nswaptypes++;
iarg += 2;
} else if (strcmp(arg[iarg], "rates") == 0) {
if (iarg + atom->ntypes >= narg)
utils::missing_cmd_args(FLERR, "fix neighbor/swap rates", error);
iarg++;
int i = 0;
while (iarg < narg) {
if (is_keyword(arg[iarg])) break;
if (i >= atom->ntypes) error->all(FLERR, "Too many values for fix neighbor/swap rates");
rate_list[i] = utils::numeric(FLERR, arg[iarg], false, lmp);
i++;
iarg++;
}
rates_flag = 1;
if (i != atom->ntypes)
error->all(FLERR, "Fix neighbor/swap rates keyword must have exactly {} arguments",
atom->ntypes);
} else {
error->all(FLERR, "Unknown fix neighbor/swap keyword: {}", arg[iarg]);
}
}
// checks
if (!nswaptypes && !diff_flag)
error->all(FLERR, Error::NOLASTLINE,
"Must specify at either 'types' or 'diff' keyword with fix neighbor/swap");
if (nswaptypes < 2 && !diff_flag)
error->all(FLERR, Error::NOLASTLINE,
"Must specify at least 2 atom types in fix neighbor/swap 'types' keyword");
}
/* ---------------------------------------------------------------------- */
int FixNeighborSwap::setmask()
{
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNeighborSwap::init()
{
c_pe = modify->get_compute_by_id("thermo_pe");
if (!c_pe) error->all(FLERR, Error::NOLASTLINE, "Could not find 'thermo_pe' compute");
c_voro = modify->get_compute_by_id(id_voro);
if (!c_voro)
error->all(FLERR, Error::NOLASTLINE, "Could not find compute voronoi ID {}", id_voro);
// set index and check validity of region
if (idregion) {
region = domain->get_region_by_id(idregion);
if (!region)
error->all(FLERR, Error::NOLASTLINE, "Region {} for fix neighbor/swap does not exist",
idregion);
}
for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++)
if (type_list[iswaptype] <= 0 || type_list[iswaptype] > atom->ntypes)
error->all(FLERR, Error::NOLASTLINE, "Invalid atom type in fix neighbor/swap command");
int *type = atom->type;
if (atom->q_flag) {
double qmax, qmin;
int firstall, first;
memory->create(qtype, nswaptypes, "neighbor/swap:qtype");
for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++) {
first = 1;
for (int i = 0; i < atom->nlocal; i++) {
if (atom->mask[i] & groupbit) {
if (type[i] == type_list[iswaptype]) {
if (first > 0) {
qtype[iswaptype] = atom->q[i];
first = 0;
} else if (qtype[iswaptype] != atom->q[i])
first = -1;
}
}
}
MPI_Allreduce(&first, &firstall, 1, MPI_INT, MPI_MIN, world);
if (firstall < 0)
error->all(FLERR, Error::NOLASTLINE,
"All atoms of a swapped type must have the same charge");
if (firstall > 0)
error->all(FLERR, Error::NOLASTLINE,
"At least one atom of each swapped type must be present to define charges");
if (first) qtype[iswaptype] = -DBL_MAX;
MPI_Allreduce(&qtype[iswaptype], &qmax, 1, MPI_DOUBLE, MPI_MAX, world);
if (first) qtype[iswaptype] = DBL_MAX;
MPI_Allreduce(&qtype[iswaptype], &qmin, 1, MPI_DOUBLE, MPI_MIN, world);
if (qmax != qmin)
error->all(FLERR, Error::NOLASTLINE, "All atoms of a swapped type must have same charge.");
}
}
memory->create(sqrt_mass_ratio, atom->ntypes + 1, atom->ntypes + 1,
"neighbor/swap:sqrt_mass_ratio");
for (int itype = 1; itype <= atom->ntypes; itype++)
for (int jtype = 1; jtype <= atom->ntypes; jtype++)
sqrt_mass_ratio[itype][jtype] = sqrt(atom->mass[itype] / atom->mass[jtype]);
// check to see if itype and jtype cutoffs are the same
// if not, reneighboring will be needed between swaps
double **cutsq = force->pair->cutsq;
unequal_cutoffs = false;
for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++)
for (int jswaptype = 0; jswaptype < nswaptypes; jswaptype++)
for (int ktype = 1; ktype <= atom->ntypes; ktype++)
if (cutsq[type_list[iswaptype]][ktype] != cutsq[type_list[jswaptype]][ktype])
unequal_cutoffs = true;
// check that no swappable atoms are in atom->firstgroup
// swapping such an atom might not leave firstgroup atoms first
if (atom->firstgroup >= 0) {
int *mask = atom->mask;
int firstgroupbit = group->bitmask[atom->firstgroup];
int flag = 0;
for (int i = 0; i < atom->nlocal; i++)
if ((mask[i] == groupbit) && (mask[i] && firstgroupbit)) flag = 1;
int flagall;
MPI_Allreduce(&flag, &flagall, 1, MPI_INT, MPI_SUM, world);
if (flagall)
error->all(FLERR, Error::NOLASTLINE,
"Cannot use fix neighbor/swap on atoms in atom_modify first group");
}
}
/* ----------------------------------------------------------------------
attempt Monte Carlo swaps
------------------------------------------------------------------------- */
void FixNeighborSwap::pre_exchange()
{
// just return if should not be called on this timestep
if (next_reneighbor != update->ntimestep) return;
// ensure current system is ready to compute energy
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build(1);
// energy_stored = energy of current state
// will be updated after accepted swaps
energy_stored = energy_full();
// attempt Ncycle atom swaps
int nsuccess = 0;
update_iswap_atoms_list();
for (int i = 0; i < ncycles; i++) nsuccess += attempt_swap();
// udpate MC stats
nswap_attempts += ncycles;
nswap_successes += nsuccess;
next_reneighbor = update->ntimestep + nevery;
}
/* ----------------------------------------------------------------------
attempt a swap of a pair of atoms
compare before/after energy and accept/reject the swap
------------------------------------------------------------------------- */
int FixNeighborSwap::attempt_swap()
{
if (niswap == 0) return 0;
// pre-swap energy
double energy_before = energy_stored;
// pick a random atom i
int i = pick_i_swap_atom();
// build nearest-neighbor list based on atom i
build_i_neighbor_list(i);
if (njswap <= 0) return 0;
// pick a neighbor atom j based on i neighbor list
jtype_selected = -1;
int j = pick_j_swap_neighbor();
int itype = type_list[0];
int jtype = jtype_selected;
// Accept swap if types are equal, no change to system
if (itype == jtype) { return 1; }
// swap their properties
if (i >= 0) {
atom->type[i] = jtype;
if (atom->q_flag) atom->q[i] = qtype[jtype_selected];
}
if (j >= 0) {
atom->type[j] = itype;
if (atom->q_flag) atom->q[j] = qtype[0];
}
// if unequal_cutoffs, call comm->borders() and rebuild neighbor list
// else communicate ghost atoms
// call to comm->exchange() is a no-op but clears ghost atoms
if (unequal_cutoffs) {
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build(1);
} else {
comm->forward_comm(this);
}
// post-swap energy
double energy_after = energy_full();
// if swap accepted, return 1
// if ke_flag, rescale atom velocities
if (random_equal->uniform() < exp(beta * (energy_before - energy_after))) {
update_iswap_atoms_list();
if (ke_flag) {
if (i >= 0) {
atom->v[i][0] *= sqrt_mass_ratio[itype][jtype];
atom->v[i][1] *= sqrt_mass_ratio[itype][jtype];
atom->v[i][2] *= sqrt_mass_ratio[itype][jtype];
}
if (j >= 0) {
atom->v[j][0] *= sqrt_mass_ratio[jtype][itype];
atom->v[j][1] *= sqrt_mass_ratio[jtype][itype];
atom->v[j][2] *= sqrt_mass_ratio[jtype][itype];
}
}
energy_stored = energy_after;
return 1;
}
// swap not accepted, return 0
// restore the swapped itype & jtype atoms
// do not need to re-call comm->borders() and rebuild neighbor list
// since will be done on next cycle or in Verlet when this fix finishes
if (i >= 0) {
atom->type[i] = itype;
if (atom->q_flag) atom->q[i] = qtype[0];
}
if (j >= 0) {
atom->type[j] = jtype;
if (atom->q_flag) atom->q[j] = qtype[jtype_selected];
}
return 0;
}
/* ----------------------------------------------------------------------
compute system potential energy
------------------------------------------------------------------------- */
double FixNeighborSwap::energy_full()
{
int eflag = 1;
int vflag = 0;
if (modify->n_pre_force) modify->pre_force(vflag);
if (force->pair) force->pair->compute(eflag, vflag);
if (atom->molecular != Atom::ATOMIC) {
if (force->bond) force->bond->compute(eflag, vflag);
if (force->angle) force->angle->compute(eflag, vflag);
if (force->dihedral) force->dihedral->compute(eflag, vflag);
if (force->improper) force->improper->compute(eflag, vflag);
}
if (force->kspace) force->kspace->compute(eflag, vflag);
if (modify->n_post_force_any) modify->post_force(vflag);
update->eflag_global = update->ntimestep;
double total_energy = c_pe->compute_scalar();
return total_energy;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int FixNeighborSwap::pick_i_swap_atom()
{
tagint *id = atom->tag;
int i = -1;
int iwhichglobal = static_cast<int>(niswap * random_equal->uniform());
if ((iwhichglobal >= niswap_before) && (iwhichglobal < niswap_before + niswap_local)) {
int iwhichlocal = iwhichglobal - niswap_before;
i = local_swap_iatom_list[iwhichlocal];
MPI_Allreduce(&id[i], &id_center, 1, MPI_INT, MPI_MAX, world);
} else {
id_center = -1;
}
return i;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int FixNeighborSwap::pick_j_swap_neighbor()
{
int j = -1;
int jtype_selected_local = -1;
// Generate random double from 0 to maximum global probability
double selected_prob = static_cast<double>(global_probability * random_equal->uniform());
// Find which local swap atom corresponds to probability
if ((selected_prob >= prev_probability) &&
(selected_prob < prev_probability + local_probability)) {
double search_prob = selected_prob - prev_probability;
for (int n = 0; n < njswap_local; n++) {
if (search_prob > local_swap_probability[n]) {
search_prob -= local_swap_probability[n];
} else {
j = local_swap_neighbor_list[n];
jtype_selected_local = local_swap_type_list[n];
MPI_Allreduce(&jtype_selected_local, &jtype_selected, 1, MPI_INT, MPI_MAX, world);
return j;
}
}
error->all(FLERR, Error::NOLASTLINE, "Did not select local neighbor swap atom");
}
MPI_Allreduce(&jtype_selected_local, &jtype_selected, 1, MPI_INT, MPI_MAX, world);
return j;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixNeighborSwap::build_i_neighbor_list(int i_center)
{
int nghost = atom->nghost;
int nlocal = atom->nlocal;
int *type = atom->type;
double **x = atom->x;
tagint *id = atom->tag;
// Allocate local_swap_neighbor_list size
memory->sfree(local_swap_neighbor_list);
atom_swap_nmax = atom->nmax;
local_swap_neighbor_list =
(int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_neighbor_list");
memory->sfree(local_swap_probability);
local_swap_probability = (double *) memory->smalloc(atom_swap_nmax * sizeof(double),
"MCSWAP:local_swap_probability_list");
memory->sfree(local_swap_type_list);
local_swap_type_list =
(int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_type_list");
// Compute voronoi and access neighbor list
c_voro->compute_local();
voro_neighbor_list = c_voro->array_local;
njswap_local = 0;
local_probability = 0.0;
for (int n = 0; n < c_voro->size_local_rows; n++) {
int temp_j_id = -1;
int temp_j = -1;
// Find local voronoi entry with selected central atom
if ((int) voro_neighbor_list[n][0] == id_center) {
temp_j_id = voro_neighbor_list[n][1];
temp_j = -1;
} else if (((int) voro_neighbor_list[n][1] == id_center) && (i_center < 0)) {
temp_j_id = voro_neighbor_list[n][0];
temp_j = -1;
} else {
continue;
}
// Find which local atom corresponds to neighbor
for (int j = 0; j < nlocal; j++) {
if (temp_j_id == id[j]) {
temp_j = j;
break;
}
}
// If temp_j not on this processor, skip
if (temp_j < 0) continue;
if (region) {
if (region->match(x[temp_j][0], x[temp_j][1], x[temp_j][2]) == 1) {
if (atom->mask[temp_j] & groupbit) {
if (diff_flag) {
// Calculate distance from i to each j, adjust probability of selection
// Get distance if own center atom
double r = INFINITY;
// Get local id of ghost center atom when ghost
for (int i = nlocal; i < nlocal + nghost; i++) {
double rtmp = sqrt(distsq3(x[temp_j], x[i]));
if ((id[i] == id_center) && (rtmp < r)) r = rtmp;
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-square(r * inv_r_0));
} else {
local_swap_probability[njswap_local] = exp(-square(r * inv_r_0));
}
local_probability += local_swap_probability[njswap_local];
local_swap_type_list[njswap_local] = type[temp_j];
local_swap_neighbor_list[njswap_local] = temp_j;
njswap_local++;
} else {
for (int jswaptype = 1; jswaptype < nswaptypes; jswaptype++) {
if (type[temp_j] == type_list[jswaptype]) {
// Calculate distance from i to each j, adjust probability of selection
// Get distance if own center atom
double r = INFINITY;
// Get local id of ghost center atom when ghost
for (int i = nlocal; i < nlocal + nghost; i++) {
double rtmp = sqrt(distsq3(x[temp_j], x[i]));
if ((id[i] == id_center) && (rtmp < r)) r = rtmp;
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-square(r * inv_r_0));
} else {
local_swap_probability[njswap_local] = exp(-square(r * inv_r_0));
}
local_probability += local_swap_probability[njswap_local];
local_swap_type_list[njswap_local] = jswaptype;
local_swap_neighbor_list[njswap_local] = temp_j;
njswap_local++;
}
}
}
}
}
} else {
if (atom->mask[temp_j] & groupbit) {
if (diff_flag) {
// Calculate distance from i to each j, adjust probability of selection
// Get distance if own center atom
double r = INFINITY;
// Get local id of ghost center atoms
for (int i = nlocal; i < nlocal + nghost; i++) {
double rtmp = sqrt(distsq3(x[temp_j], x[i]));
if ((id[i] == id_center) && (rtmp < r)) r = rtmp;
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-square(r * inv_r_0));
} else {
local_swap_probability[njswap_local] = exp(-square(r * inv_r_0));
}
local_probability += local_swap_probability[njswap_local];
local_swap_type_list[njswap_local] = type[temp_j];
local_swap_neighbor_list[njswap_local] = temp_j;
njswap_local++;
} else {
for (int jswaptype = 1; jswaptype < nswaptypes; jswaptype++) {
if (type[temp_j] == type_list[jswaptype]) {
// Calculate distance from i to each j, adjust probability of selection
// Get distance if own center atom
double r = INFINITY;
// Get local id of ghost center atom when ghost
for (int i = nlocal; i < nlocal + nghost; i++) {
double rtmp = sqrt(distsq3(x[temp_j], x[i]));
if ((id[i] == id_center) && (rtmp < r)) r = rtmp;
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-square(r * inv_r_0));
} else {
local_swap_probability[njswap_local] = exp(-square(r * inv_r_0));
}
local_probability += local_swap_probability[njswap_local];
local_swap_type_list[njswap_local] = jswaptype;
local_swap_neighbor_list[njswap_local] = temp_j;
njswap_local++;
}
}
}
}
}
}
MPI_Allreduce(&njswap_local, &njswap, 1, MPI_INT, MPI_SUM, world);
MPI_Scan(&njswap_local, &njswap_before, 1, MPI_INT, MPI_SUM, world);
njswap_before -= njswap_local;
MPI_Allreduce(&local_probability, &global_probability, 1, MPI_DOUBLE, MPI_SUM, world);
MPI_Scan(&local_probability, &prev_probability, 1, MPI_DOUBLE, MPI_SUM, world);
prev_probability -= local_probability;
}
/* ----------------------------------------------------------------------
update the list of swap atoms
------------------------------------------------------------------------- */
void FixNeighborSwap::update_iswap_atoms_list()
{
int nlocal = atom->nlocal;
int *type = atom->type;
double **x = atom->x;
if (atom->nmax > atom_swap_nmax) {
memory->sfree(local_swap_iatom_list);
atom_swap_nmax = atom->nmax;
local_swap_iatom_list =
(int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_iatom_list");
}
niswap_local = 0;
if (region) {
for (int i = 0; i < nlocal; i++) {
if (region->match(x[i][0], x[i][1], x[i][2]) == 1) {
if (atom->mask[i] & groupbit) {
if (type[i] == type_list[0]) {
local_swap_iatom_list[niswap_local] = i;
niswap_local++;
}
}
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (atom->mask[i] & groupbit) {
if (type[i] == type_list[0]) {
local_swap_iatom_list[niswap_local] = i;
niswap_local++;
}
}
}
}
MPI_Allreduce(&niswap_local, &niswap, 1, MPI_INT, MPI_SUM, world);
MPI_Scan(&niswap_local, &niswap_before, 1, MPI_INT, MPI_SUM, world);
niswap_before -= niswap_local;
}
/* ---------------------------------------------------------------------- */
int FixNeighborSwap::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
int * /*pbc*/)
{
int i, j, m;
int *type = atom->type;
double *q = atom->q;
m = 0;
if (atom->q_flag) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = type[j];
buf[m++] = q[j];
}
} else {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = type[j];
}
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixNeighborSwap::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
int *type = atom->type;
double *q = atom->q;
m = 0;
last = first + n;
if (atom->q_flag) {
for (i = first; i < last; i++) {
type[i] = static_cast<int>(buf[m++]);
q[i] = buf[m++];
}
} else {
for (i = first; i < last; i++) type[i] = static_cast<int>(buf[m++]);
}
}
/* ----------------------------------------------------------------------
return acceptance ratio
------------------------------------------------------------------------- */
double FixNeighborSwap::compute_vector(int n)
{
if (n == 0) return nswap_attempts;
if (n == 1) return nswap_successes;
return 0.0;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixNeighborSwap::memory_usage()
{
double bytes = (double) atom_swap_nmax * sizeof(int);
return bytes;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixNeighborSwap::write_restart(FILE *fp)
{
int n = 0;
double list[6];
list[n++] = random_equal->state();
list[n++] = ubuf(next_reneighbor).d;
list[n++] = nswap_attempts;
list[n++] = nswap_successes;
list[n++] = ubuf(update->ntimestep).d;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size, sizeof(int), 1, fp);
fwrite(list, sizeof(double), n, fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixNeighborSwap::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
seed = static_cast<int>(list[n++]);
random_equal->reset(seed);
next_reneighbor = (bigint) ubuf(list[n++]).i;
nswap_attempts = static_cast<int>(list[n++]);
nswap_successes = static_cast<int>(list[n++]);
bigint ntimestep_restart = (bigint) ubuf(list[n++]).i;
if (ntimestep_restart != update->ntimestep)
error->all(FLERR, Error::NOLASTLINE,
"Must not reset timestep when restarting fix neighbor/swap");
}

98
src/MC/fix_neighbor_swap.h
View File

@ -1,98 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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(neighbor/swap,FixNeighborSwap);
// clang-format on
#else
#ifndef LMP_FIX_NEIGH_SWAP_H
#define LMP_FIX_NEIGH_SWAP_H
#include "fix.h"
namespace LAMMPS_NS {
class FixNeighborSwap : public Fix {
public:
FixNeighborSwap(class LAMMPS *, int, char **);
~FixNeighborSwap();
int setmask();
void init();
void pre_exchange();
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
double compute_vector(int);
double memory_usage();
void write_restart(FILE *);
void restart(char *);
private:
int nevery, seed;
int ke_flag; // yes = conserve ke, no = do not conserve ke
int diff_flag; // yes = simulate diffusion of central atom, no = swap only to certain types
int rates_flag; // yes = use modified type rates, no = swap rates are equivilent across types
int ncycles;
int niswap, njswap; // # of i,j swap atoms on all procs
int niswap_local, njswap_local; // # of swap atoms on this proc
int niswap_before, njswap_before; // # of swap atoms on procs < this proc
class Region *region; // swap region
char *idregion; // swap region id
int nswaptypes;
int jtype_selected;
int id_center;
double x_center;
double y_center;
double z_center;
int *type_list;
double *rate_list;
double nswap_attempts;
double nswap_successes;
bool unequal_cutoffs;
int atom_swap_nmax;
double beta, inv_r_0;
double local_probability; // Total swap probability stored on this proc
double global_probability; // Total swap probability across all proc
double prev_probability; // Swap probability on proc < this proc
double *qtype;
double energy_stored;
double **sqrt_mass_ratio;
double **voro_neighbor_list;
int *local_swap_iatom_list;
int *local_swap_neighbor_list;
int *local_swap_type_list; // Type list index of atoms stored on this proc
double *local_swap_probability;
class RanPark *random_equal;
char *id_voro;
class Compute *c_voro, *c_pe;
void options(int, char **);
int attempt_swap();
double energy_full();
int pick_i_swap_atom();
int pick_j_swap_neighbor();
void build_i_neighbor_list(int);
void update_iswap_atoms_list();
};
} // namespace LAMMPS_NS
#endif
#endif

7
src/tokenizer.cpp
View File

@ -19,6 +19,8 @@
#include "fmt/format.h"
#include "utils.h"
#include <cmath>
#include <cstdlib>
#include <utility>
using namespace LAMMPS_NS;
@ -365,6 +367,11 @@ double ValueTokenizer::next_double()
return val;
// rethrow exceptions from std::stod()
} catch (std::out_of_range const &) {
// could be a denormal number. try again with std::strtod().
char *end;
auto val = std::strtod(current.c_str(), &end);
// return value of denormal
if ((val > -HUGE_VAL) && (val < HUGE_VAL)) return val;
throw InvalidFloatException(current);
} catch (std::invalid_argument const &) {
throw InvalidFloatException(current);

15
src/utils.cpp
View File

@ -34,6 +34,8 @@
#include <cctype>
#include <cerrno>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <stdexcept>
@ -272,6 +274,11 @@ void utils::print(FILE *fp, const std::string &mesg)
fputs(mesg.c_str(), fp);
}
void utils::print(const std::string &mesg)
{
fputs(mesg.c_str(), stdout);
}
void utils::fmtargs_print(FILE *fp, fmt::string_view format, fmt::format_args args)
{
print(fp, fmt::vformat(format, args));
@ -529,7 +536,7 @@ double utils::numeric(const char *file, int line, const std::string &str, bool d
lmp->error->all(file, line, msg);
}
double rv = 0;
double rv = 0.0;
auto msg = fmt::format("Floating point number {} in input script or data file is invalid", buf);
try {
std::size_t endpos;
@ -546,6 +553,12 @@ double utils::numeric(const char *file, int line, const std::string &str, bool d
else
lmp->error->all(file, line, msg);
} catch (std::out_of_range const &) {
// could be a denormal number. try again with std::strtod().
char *end;
rv = std::strtod(buf.c_str(), &end);
// return value if denormal
if ((rv > -HUGE_VAL) && (rv < HUGE_VAL)) return rv;
msg = fmt::format("Floating point number {} in input script or data file is out of range", buf);
if (do_abort)
lmp->error->one(file, line, msg);

22
src/utils.h
View File

@ -151,7 +151,7 @@ output are compressed to a single blank by calling :cpp:func:`strcompress()`
\endverbatim
*
* This function implements a version of fprintf() that uses {fmt} formatting
* This function implements a version of (f)printf() that uses {fmt} formatting
*
* \param fp stdio FILE pointer
* \param format format string of message to be printed
@ -163,12 +163,32 @@ output are compressed to a single blank by calling :cpp:func:`strcompress()`
}
/*! \overload
*
* Print to stdout without specifying the FILE pointer.
*
* \param mesg string with message to be printed */
template <typename... Args> void print(const std::string &format, Args &&...args)
{
fmtargs_print(stdout, format, fmt::make_format_args(args...));
}
/*! \overload
*
* Print string message without format
*
* \param fp stdio FILE pointer
* \param mesg string with message to be printed */
void print(FILE *fp, const std::string &mesg);
/*! \overload
*
* Print string message without format to stdout
*
* \param mesg string with message to be printed */
void print(const std::string &mesg);
/*! Return text redirecting the user to a specific paragraph in the manual
*
* The LAMMPS manual contains detailed explanations for errors and

3
tools/lammps-gui/help_index.table
View File

@ -482,6 +482,7 @@ fix_aveforce.html fix aveforce
fix_ave_grid.html fix ave/grid
fix_ave_histo.html fix ave/histo
fix_ave_histo.html fix ave/histo/weight
fix_ave_moments.html fix ave/moments
fix_ave_time.html fix ave/time
fix_balance.html fix balance
fix_bocs.html fix bocs
@ -542,6 +543,7 @@ fix_flow_gauss.html fix flow/gauss
fix_freeze.html fix freeze
fix_freeze.html fix freeze/kk
fix_gcmc.html fix gcmc
fix_gjf.html fix gjf
fix_gld.html fix gld
fix_gle.html fix gle
fix_gravity.html fix gravity
@ -732,6 +734,7 @@ fix_saed_vtk.html fix saed/vtk
fix_setforce.html fix setforce
fix_setforce.html fix setforce/kk
fix_setforce.html fix setforce/spin
fix_set.html fix set
fix_sgcmc.html fix sgcmc
fix_shake.html fix rattle
fix_shake.html fix shake

2
tools/lammps-gui/lammps-gui.appdata.xml
View File

@ -57,7 +57,7 @@
<releases>
<release version="1.6.14" timestamp="1747828753">
<description>
...
Must set en_US.UTF-8 locale on macOS since it lacks support for C.UTF-8
</description>
</release>
<release version="1.6.13" timestamp="1743734509">

5
tools/lammps-gui/main.cpp
View File

@ -36,8 +36,13 @@ int main(int argc, char *argv[])
qRegisterMetaTypeStreamOperators<QList<QString>>("QList<QString>");
#endif
#ifndef Q_OS_MACOS
// enforce using the plain ASCII C locale with UTF-8 encoding within the GUI.
qputenv("LC_ALL", "C.UTF-8");
#else
// macOS does not support "C" locale with UTF-8 encoding, but Qt requires UTF-8
qputenv("LC_ALL", "en_US.UTF-8");
#endif
QApplication app(argc, argv);
QCoreApplication::setOrganizationName("The LAMMPS Developers");