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

Merge branch 'develop' into cmake-adjustments

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-05-19 21:03:07 -04:00
parent fd8fb74a88 b3402eaae8
commit ad63eb8bf7
10 changed files with 1284 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
cmake/Modules/Packages/MC.cmake
View File

@ -7,3 +7,13 @@ if(NOT PKG_MANYBODY)
list(REMOVE_ITEM LAMMPS_SOURCES ${LAMMPS_SOURCE_DIR}/MC/fix_sgcmc.cpp)
set_property(TARGET lammps PROPERTY SOURCES "${LAMMPS_SOURCES}")
endif()
# fix neighbor/swap may only be installed if also the VORONOI package is installed
if(NOT PKG_VORONOI)
get_property(LAMMPS_FIX_HEADERS GLOBAL PROPERTY FIX)
list(REMOVE_ITEM LAMMPS_FIX_HEADERS ${LAMMPS_SOURCE_DIR}/MC/fix_neighbor_swap.h)
set_property(GLOBAL PROPERTY FIX "${LAMMPS_FIX_HEADERS}")
get_target_property(LAMMPS_SOURCES lammps SOURCES)
list(REMOVE_ITEM LAMMPS_SOURCES ${LAMMPS_SOURCE_DIR}/MC/fix_neighbor_swap.cpp)
set_property(TARGET lammps PROPERTY SOURCES "${LAMMPS_SOURCES}")
endif()

1
doc/src/Commands_fix.rst
View File

@ -112,6 +112,7 @@ 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>`

1
doc/src/fix.rst
View File

@ -291,6 +291,7 @@ 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 Normal file
View File

@ -0,0 +1,264 @@
.. 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).

2
doc/utils/sphinx-config/false_positives.txt
View File

@ -1859,6 +1859,7 @@ Kloss
Kloza
kmax
Kmax
kMC
KMP
kmu
Knizhnik
@ -4131,6 +4132,7 @@ volpress
volumetric
von
Voro
voro
Vorobyov
voronoi
Voronoi

2
src/.gitignore vendored
View File

@ -926,6 +926,8 @@
/fix_msst.h
/fix_neb.cpp
/fix_neb.h
/fix_neighbor_swap.cpp
/fix_neighbor_swap.h
/fix_nh_asphere.cpp
/fix_nh_asphere.h
/fix_nph_asphere.cpp

1
src/Depend.sh
View File

@ -146,6 +146,7 @@ fi
if (test $1 = "MC") then
depend MISC
depend VORONOI
fi
if (test $1 = "MEAM") then

2
src/MC/Install.sh
View File

@ -51,6 +51,8 @@ action fix_charge_regulation.cpp
action fix_charge_regulation.h
action fix_gcmc.cpp
action fix_gcmc.h
action fix_neighbor_swap.cpp compute_voronoi_atom.cpp
action fix_neighbor_swap.h compute_voronoi_atom.h
action fix_mol_swap.cpp
action fix_mol_swap.h
action fix_sgcmc.cpp pair_eam.cpp

901
src/MC/fix_neighbor_swap.cpp Normal file
View File

@ -0,0 +1,901 @@
/* ----------------------------------------------------------------------
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 "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_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 <cctype>
#include <cfloat>
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
static const char cite_fix_neighbor_swap_c[] =
"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://www.sciencedirect.com/science/article/pii/S0927025622006401}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), region(nullptr), idregion(nullptr), type_list(nullptr), qtype(nullptr),
c_voro(nullptr), voro_neighbor_list(nullptr), sqrt_mass_ratio(nullptr),
local_swap_iatom_list(nullptr), random_equal(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;
// 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);
r_0 = utils::inumeric(FLERR, arg[7], false, lmp);
// Voro compute check
int icompute = modify->find_compute(utils::strdup(arg[8]));
if (icompute < 0) error->all(FLERR, "Could not find neighbor compute ID");
c_voro = modify->compute[icompute];
if (c_voro->local_flag == 0)
error->all(FLERR, "Neighbor compute does not compute local info");
if (c_voro->size_local_cols != 3)
error->all(FLERR, "Neighbor compute does not give i, j, size as expected");
if (nevery <= 0) error->all(FLERR, "Illegal fix neighbor/swap command nevery value");
if (ncycles < 0) error->all(FLERR, "Illegal fix neighbor/swap command ncycles value");
if (seed <= 0) error->all(FLERR, "Illegal fix neighbor/swap command seed value");
if (temperature <= 0.0) error->all(FLERR, "Illegal fix neighbor/swap command temperature value");
beta = 1.0 / (force->boltz * temperature);
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 - 8, &arg[8]);
// 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;
voro_neighbor_list = nullptr;
local_swap_iatom_list = nullptr;
local_swap_neighbor_list = nullptr;
local_swap_probability = nullptr;
local_swap_type_list = nullptr;
// 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 random_equal;
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
void FixNeighborSwap::options(int narg, char **arg)
{
if (narg < 0) error->all(FLERR, "Illegal fix neighbor/swap command\n");
ke_flag = 1;
diff_flag = 0;
rates_flag = 0;
nswaptypes = 0;
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg], "region") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix neighbor/swap command");
region = domain->get_region_by_id(arg[iarg + 1]);
if (!region) error->all(FLERR, "Region ID for fix neighbor/swap does not exist");
idregion = utils::strdup(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg], "ke") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix neighbor/swap command");
ke_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "types") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal fix neighbor/swap command");
if (diff_flag != 0) error->all(FLERR, "Illegal fix neighbor/swap command");
iarg++;
nswaptypes = 0;
while (iarg < narg) {
if (isalpha(arg[iarg][0])) break;
if (nswaptypes >= atom->ntypes) error->all(FLERR, "Illegal fix neighbor/swap command");
type_list[nswaptypes] = utils::numeric(FLERR, arg[iarg], false, lmp);
nswaptypes++;
iarg++;
}
} else if (strcmp(arg[iarg], "diff") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix neighbor/swap command");
if (nswaptypes != 0) error->all(FLERR, "Illegal fix neighbor/swap command");
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) error->all(FLERR, "Illegal fix neighbor/swap command");
iarg++;
int i = 0;
while (iarg < narg) {
if (isalpha(arg[iarg][0])) break;
if (i >= atom->ntypes) error->all(FLERR, "Illegal fix neighbor/swap command");
rate_list[i] = utils::numeric(FLERR, arg[iarg], false, lmp);
i++;
iarg++;
}
rates_flag = 1;
if (i != atom->ntypes) error->all(FLERR, "Illegal fix neighbor/swap command");
} else
error->all(FLERR, "Illegal fix neighbor/swap command");
}
}
/* ---------------------------------------------------------------------- */
int FixNeighborSwap::setmask()
{
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNeighborSwap::init()
{
c_pe = modify->get_compute_by_id("thermo_pe");
int *type = atom->type;
if (nswaptypes < 2 && !diff_flag)
error->all(FLERR, "Must specify at least 2 types in fix neighbor/swap command");
// set index and check validity of region
if (idregion) {
region = domain->get_region_by_id(idregion);
if (!region) error->all(FLERR, "Region {} for fix setforce does not exist", idregion);
}
for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++)
if (type_list[iswaptype] <= 0 || type_list[iswaptype] > atom->ntypes)
error->all(FLERR, "Invalid atom type in fix neighbor/swap command");
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) {
qtype[iswaptype] = atom->q[i];
first = 0;
} else if (qtype[iswaptype] != atom->q[i])
error->one(FLERR, "All atoms of a swapped type must have the same charge.");
}
}
}
MPI_Allreduce(&first, &firstall, 1, MPI_INT, MPI_MIN, world);
if (firstall)
error->all(FLERR,
"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, "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, "Cannot do 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()
{
// int nlocal = atom->nlocal;
tagint *id = atom->tag;
if (niswap == 0) return 0;
// pre-swap energy
double energy_before = energy_stored;
// pick a random atom i
int i = pick_i_swap_atom();
// get global id and position of atom i
// get_global_i(i);
// 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(i);
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 id_center_local = -1;
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];
id_center_local = id[i];
MPI_Allreduce(&id[i], &id_center, 1, MPI_INT, MPI_MAX, world);
} else {
MPI_Allreduce(&id[i], &id_center, 1, MPI_INT, MPI_MAX, world);
}
return i;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int FixNeighborSwap::pick_j_swap_neighbor(int i)
{
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, "Did not select local neighbor swap atom");
}
MPI_Allreduce(&jtype_selected_local, &jtype_selected, 1, MPI_INT, MPI_MAX, world);
return j;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
double FixNeighborSwap::get_distance(double *i, double *j)
{
double r = sqrt(MathSpecial::square((i[0] - j[0])) + MathSpecial::square((i[1] - j[1])) +
MathSpecial::square((i[2] - j[2])));
return r;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
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;
if (i_center >= 0) { double r = get_distance(x[temp_j], x[i_center]); }
// Get local id of ghost center atom when ghost
for (int i = nlocal; i < nlocal + nghost; i++) {
if ((id[i] == id_center) && (get_distance(x[temp_j], x[i]) < r)) {
r = get_distance(x[temp_j], x[i]);
}
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0));
} else {
local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 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;
if (i_center >= 0) { double r = get_distance(x[temp_j], x[i_center]); }
// Get local id of ghost center atom when ghost
for (int i = nlocal; i < nlocal + nghost; i++) {
if ((id[i] == id_center) && (get_distance(x[temp_j], x[i]) < r)) {
r = get_distance(x[temp_j], x[i]);
}
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0));
} else {
local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 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;
if (i_center >= 0) { r = get_distance(x[temp_j], x[i_center]); }
// Get local id of ghost center atoms
for (int i = nlocal; i < nlocal + nghost; i++) {
if ((id[i] == id_center) && (get_distance(x[temp_j], x[i]) < r))
r = get_distance(x[temp_j], x[i]);
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0));
} else {
local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 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;
if (i_center >= 0) { double r = get_distance(x[temp_j], x[i_center]); }
// Get local id of ghost center atom when ghost
for (int i = nlocal; i < nlocal + nghost; i++) {
if ((id[i] == id_center) && (get_distance(x[temp_j], x[i]) < r)) {
r = get_distance(x[temp_j], x[i]);
}
}
if (rates_flag) {
local_swap_probability[njswap_local] =
rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0));
} else {
local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 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, "Must not reset timestep when restarting fix neighbor/swap");
}

100
src/MC/fix_neighbor_swap.h Normal file
View File

@ -0,0 +1,100 @@
/* -*- 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
// int global_i_ID; // global id of selected i atom
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, 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;
class Compute *c_voro;
class Compute *c_pe;
void options(int, char **);
int attempt_swap();
double energy_full();
int pick_i_swap_atom();
int pick_j_swap_neighbor(int);
double get_distance(double[3], double[3]);
void build_i_neighbor_list(int);
void update_iswap_atoms_list();
};
} // namespace LAMMPS_NS
#endif
#endif