From f980e8babf44f3a6fef948c4e03922a41a3e4dcb Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Thu, 19 Dec 2024 09:52:05 -0700 Subject: [PATCH 01/22] Integrated MD-KMC code to current development branch --- src/Depend.sh | 1 + src/MC/Install.sh | 2 + src/MC/fix_neighbor_swap.cpp | 1012 ++++++++++++++++++++++++++++++++++ src/MC/fix_neighbor_swap.h | 146 +++++ 4 files changed, 1161 insertions(+) create mode 100755 src/MC/fix_neighbor_swap.cpp create mode 100755 src/MC/fix_neighbor_swap.h diff --git a/src/Depend.sh b/src/Depend.sh index 9ddb29450d..ba55deb62c 100755 --- a/src/Depend.sh +++ b/src/Depend.sh @@ -146,6 +146,7 @@ fi if (test $1 = "MC") then depend MISC + depend VORONOI fi if (test $1 = "MEAM") then diff --git a/src/MC/Install.sh b/src/MC/Install.sh index ccf6767c4d..9a0095fd34 100755 --- a/src/MC/Install.sh +++ b/src/MC/Install.sh @@ -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 +action fix_neighbor_swap.h action fix_mol_swap.cpp action fix_mol_swap.h action fix_sgcmc.cpp pair_eam.cpp diff --git a/src/MC/fix_neighbor_swap.cpp b/src/MC/fix_neighbor_swap.cpp new file mode 100755 index 0000000000..52730ee6f5 --- /dev/null +++ b/src/MC/fix_neighbor_swap.cpp @@ -0,0 +1,1012 @@ +// clang-format off +/* ---------------------------------------------------------------------- + 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 "dihedral.h" +#include "domain.h" +#include "error.h" +#include "fix.h" +#include "force.h" +#include "group.h" +#include "improper.h" +#include "kspace.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" +#include "pair.h" +#include "random_park.h" +#include "region.h" +#include "update.h" +#include "compute_voronoi_atom.h" + +#include +#include +#include +#include + +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), random_unequal(nullptr), c_pe(nullptr) +{ + if (narg < 10) error->all(FLERR,"Illegal fix neighbor/swap command"); + + // utils::logmesg(lmp,"Starting neighbor/swap\n"); + + 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); + + if (nevery <= 0) error->all(FLERR,"Illegal fix neighbor/swap command"); + if (ncycles < 0) error->all(FLERR,"Illegal fix neighbor/swap command"); + // if (nevery < ncycles) error->all(FLERR,"Illegal fix neighbor/swap command"); + if (seed <= 0) error->all(FLERR,"Illegal fix neighbor/swap command"); + if (temperature <= 0.0) error->all(FLERR,"Illegal fix neighbor/swap command"); + + 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-7,&arg[7]); + + if (voro_flag != 1) error->all(FLERR,"Voro Compute Required for fix neighbor/swap command"); + + // random number generator, same for all procs + + random_equal = new RanPark(lmp,seed); + + // random number generator, not the same for all procs + + random_unequal = 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; + delete random_unequal; +} + +/* ---------------------------------------------------------------------- + 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; + voro_flag = 0; + nswaptypes = 0; + + int iarg = 0; + while (iarg < narg) { + // utils::logmesg(lmp,"Parsing Argument {}\n", arg[iarg]); + 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],"voro") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix neighbor/swap command"); + + int icompute = modify->find_compute(utils::strdup(arg[iarg+1])); + + 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"); + + voro_flag = 1; + iarg += 2; + } 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"); + // utils::logmesg(lmp,"Reading rates\n"); + 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() +{ + // utils::logmesg(lmp,"Initializing neighbor/swap\n"); + 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"); + + // utils::logmesg(lmp,"Done Init\n"); + } +} + +/* ---------------------------------------------------------------------- + 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 + + // utils::logmesg(lmp,"Attempting Swap\n"); + + int nsuccess = 0; + update_iswap_atoms_list(); + // utils::logmesg(lmp,"Updated Atom i List\n"); + for (int i = 0; i < ncycles; i++) nsuccess += attempt_swap(); + + // utils::logmesg(lmp,"\n"); + + // 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); + // utils::logmesg(lmp,"Built neighbor list with {} atoms\n",njswap); + if ( njswap <= 0 ) { + // utils::logmesg(lmp,"No valid neighbors found\n"); + return 0; + } + + // pick a neighbor atom j based on i neighbor list + jtype_selected = -1; + int j = pick_j_swap_neighbor(i); + // utils::logmesg(lmp,"Selected swap neighbor j\n"); + + int itype = type_list[0]; + int jtype = jtype_selected; + // utils::logmesg(lmp,"Selected atom j type {}\n",jtype); + + // utils::logmesg(lmp,"Swapping local atoms {} and ",i); + // utils::logmesg(lmp,"{}\n",j); + + // if ( i >= 0 || j >= 0) { + // utils::logmesg(lmp,"Selected atom j type {}\n",jtype); + // utils::logmesg(lmp,"Swapping atoms {} and ",id[i]); + // utils::logmesg(lmp,"{}\n",id[j]); + // } + + // utils::logmesg(lmp,"Atom type {} and ",itype); + // utils::logmesg(lmp,"{}\n",jtype); + + // Accept swap if types are equal, no change to system + if ( itype == jtype ){ + // utils::logmesg(lmp,"Atoms have same type, no processing needed\n"); + 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; + + // utils::logmesg(lmp,"Swap accepted\n"); + // Swap atom groups if successful swap + + // int groupi = atom->mask[i]; + // int groupj = atom->mask[j]; + + // atom->mask[i] = group->bitmask[groupj]; + // atom->mask[j] = group->bitmask[groupi]; + + 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 + + // utils::logmesg(lmp,"Swap not accepted\n"); + + 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 (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 (global_probability*random_equal->uniform()); + + // utils::logmesg(lmp,"Local Probability {}\n",local_probability); + + // Find which local swap atom corresponds to probability + if ((selected_prob >= prev_probability) && + (selected_prob < prev_probability + local_probability)){ + // utils::logmesg(lmp,"Selected Probability {}\n",selected_prob); + 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]; + // utils::logmesg(lmp,"Selecting Atom j type {}\n",jtype_selected_local); + 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"); + } + + // utils::logmesg(lmp,"Swap atom not found on processor\n"); + 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_x, r_y, r_z; + + // r_x = i[0] - j[0]; + // r_y = i[1] - j[1]; + // r_z = i[2] - j[2]; + + // // Domain::minimum_image(r_x, r_y, r_z); + // double r = sqrt(pow(r_x, 2.) + + // pow(r_y, 2.) + + // pow(r_z, 2.)); + + double r = sqrt(pow((i[0] - j[0]), 2.) + + pow((i[1] - j[1]), 2.) + + pow((i[2] - j[2]), 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; + + // utils::logmesg(lmp,"Searching for atom {}\n",id_center); + + 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; + // utils::logmesg(lmp,"Found neighbor {}\n",id[j]); + break; + } + } + + // If temp_j not on this processor, skip + if ( temp_j < 0 ) continue; + + // // If atom is already in local list, skip + // bool inlist = false; + // for ( int j = 0; j < njswap_local; j++ ){ + // if ( temp_j == local_swap_neighbor_list[j] ){ + // utils::logmesg(lmp,"Skipping atom {}, pair already counted\n",id[temp_j]); + // inlist = true; + // break; + // } + // } + // if (inlist) 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 centr 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(-pow((r/3.), 2.)); + } else { + local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + } + local_probability += local_swap_probability[njswap_local]; + + // utils::logmesg(lmp,"Adding Atom j type {}\n",type[temp_j]); + 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(-pow((r/3.), 2.)); + } else { + local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + } + 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 + // if ( i_center < 0 ){ + // utils::logmesg(lmp,"Initial distance {}\n", r); + 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]); + // utils::logmesg(lmp,"Updated distance {}\n", r); + } + } + // } + + // utils::logmesg(lmp,"Final distance {}\n", r); + + if (rates_flag) { + local_swap_probability[njswap_local] = rate_list[type[temp_j] - 1]*exp(-pow((r/3.), 2.)); + } else{ + local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + } + 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(-pow((r/3.), 2.)); + } else { + local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + } + 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 (buf[m++]); + q[i] = buf[m++]; + } + } else { + for (i = first; i < last; i++) + type[i] = static_cast (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++] = random_unequal->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 (list[n++]); + random_equal->reset(seed); + + seed = static_cast (list[n++]); + random_unequal->reset(seed); + + next_reneighbor = (bigint) ubuf(list[n++]).i; + + nswap_attempts = static_cast(list[n++]); + nswap_successes = static_cast(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"); +} diff --git a/src/MC/fix_neighbor_swap.h b/src/MC/fix_neighbor_swap.h new file mode 100755 index 0000000000..1413de2ef8 --- /dev/null +++ b/src/MC/fix_neighbor_swap.h @@ -0,0 +1,146 @@ +/* -*- 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 voro_flag; // yes = use given voronoi calculation, no = use internal voronoi calculation + 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; + 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 RanPark *random_unequal; + + 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 + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Region ID for fix neighbor/swap does not exist + +Self-explanatory. + +E: Must specify at least 2 types in fix neighbor/swap command + +Self-explanatory. + +E: Need nswaptypes mu values in fix neighbor/swap command + +Self-explanatory. + +E: Invalid atom type in fix neighbor/swap command + +The atom type specified in the neighbor/swap command does not exist. + +E: All atoms of a swapped type must have the same charge. + +Self-explanatory. + +E: At least one atom of each swapped type must be present to define charges. + +Self-explanatory. + +E: All atoms of a swapped type must have same charge. + +Self-explanatory. + +E: Cannot do neighbor/swap on atoms in atom_modify first group + +This is a restriction due to the way atoms are organized in a list to +enable the atom_modify first command. + +*/ From 82569f444826aaad99ff7bd4cdd8d66d1a389f32 Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Thu, 19 Dec 2024 11:29:19 -0700 Subject: [PATCH 02/22] Added documentation files --- doc/src/Commands_fix.rst | 1 + doc/src/fix.rst | 1 + doc/src/fix_neighbor_swap.rst | 169 ++++++++++++++++++++++++++++++++++ 3 files changed, 171 insertions(+) create mode 100644 doc/src/fix_neighbor_swap.rst diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 04d1a9969a..cbadc09880 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -110,6 +110,7 @@ OPT. * :doc:`mvv/tdpd ` * :doc:`neb ` * :doc:`neb/spin ` + * :doc:`neighbor/swap ` * :doc:`nonaffine/displacement ` * :doc:`nph (ko) ` * :doc:`nph/asphere (o) ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 9af607601b..a4a9ac4a86 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -289,6 +289,7 @@ accelerated styles exist. * :doc:`mvv/tdpd ` - constant temperature DPD using the modified velocity-Verlet algorithm * :doc:`neb ` - nudged elastic band (NEB) spring forces * :doc:`neb/spin ` - nudged elastic band (NEB) spring forces for spins +* :doc:`neighbor/swap ` - kinetic Monte Calo atom swapping * :doc:`nonaffine/displacement ` - calculate nonaffine displacement of atoms * :doc:`nph ` - constant NPH time integration via Nose/Hoover * :doc:`nph/asphere ` - NPH for aspherical particles diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst new file mode 100644 index 0000000000..a3eb1c6d82 --- /dev/null +++ b/doc/src/fix_neighbor_swap.rst @@ -0,0 +1,169 @@ +.. index:: fix neighbor/swap + +fix neighbor/swap command +===================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + fix ID group-ID neighbor/swap N X seed T R keyword values ... + +* ID, group-ID are documented in :doc:`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) +* R = scaling swap probability of the MC swaps (distance units) +* one or more keyword/value pairs may be appended to args +* keyword = *types* or *mu* or *ke* or *semi-grand* or *region* + + .. parsed-literal:: + + *types* values = two or more atom types (1-Ntypes or type label) + *ke* value = *no* or *yes* + *no* = no conservation of kinetic energy after atom swaps + *yes* = kinetic energy is conserved after atom swaps + *region* value = region-ID + region-ID = ID of region to use as an exchange/move volume + *diff* values = one atom type + *voro* values = valid voronoi compute id (compute voronoi/atom) + *rates* values = Ntype values to conduct variable diffusion for different atom types (unitless) + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix mc all neighbor/swap 10 160 15238 1000.0 diff 2 voro voroN + fix myFix all neighbor/swap 100 1 12345 298.0 region my_swap_region types 5 6 voro voroN + fix kmc all neighbor/swap 1 100 345 1.0 diff 3 rates 3 1 6 voro voroN + +Description +""""""""""" + +Computes MC evaluations to enable kinetic Monte Carlo (kMC) behavior +during MD simulation through only allowing neighboring atom swaps. +Neighboring atoms are selected using a Voronoi tesselation approach. This +implementation is as described in :ref:`(Tavenner) `. + +The fix is called every *N* timesteps and attempts *X* swaps. The system +is initialized with a random seed, using a temperature *T* for evaluating +the MC energy swaps. The distance-based probability is weighted according +to *R* which sets the radius :math:`r_0` for the weighting + +.. math:: + + p_{ij} = e^{(\frac{r_{ij}}{r_0})^2} + +where :math:`p_{ij}` is the probability of selecting atoms :math:`i` and :math:`j` for an +evaluated swap. + +The keyword *types* is submitted with two or more atom types as the value. +Atoms of the first atom type are swapped with valid neighbors of all the remaining +atom types. + +The keyword *diff* is used for implementation of simulated diffusion of a given atom type +as given by *diff type*. This command selects all atom types as acceptable swap types to a +centrally selected atom of type *type*. This includes the atom type specified by the diff +keyword to account for self-diffusion hops of an atom type with itself. + +Keyword *voro* is currently required, and is implemented as + +.. code-block:: LAMMPS + voro compute-ID + +where *compute-ID* is the ID of a corresponding Voronoi computation with neighbor list, i.e. + +.. code-block:: LAMMPS + + compute compute-ID group-ID voronoi/atom neighbors yes + +The group selected for computing *voro* should correspond to all the potential atoms to +be swapped at the initial step, i.e. + +.. code-block:: LAMMPS + group group-ID type 2 + +for using *fix neighbor/swap* with *diff 2*. + +The keyword *rates* can modify the swap rate for each swapped type by values +where the adjusted rates values are given in order of increasing atom type. +The number of rates provided must equal the number of atom types in the simulaton. +In the third provided example above, a simulation is conducted with three atom types +where the third atom type is the one sampled for attempted swaps. All three atom +types are considered valid swaps, but atoms of type 1 will be selected three times +as often as atoms of type 2. Conversely, atoms of type 3 are six times more likely to +be selected than atoms of type two and twice as likely as atoms of type 1. + +Finally, the *region* keyword is implemented as in other atomic fixes, where +the *region region-ID* command indicates that atom swaps only be considered in the area +given by *region-ID*. If only atoms of certain groups are expected to be in this region, +the corresponding compute voronoi command can be adjusted accordingly. + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This fix writes the state of the fix to :doc:`binary restart files +`. This includes information about the random number generator +seed, the next timestep for MC exchanges, the number of exchange +attempts and successes, etc. See the :doc:`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 ` options are relevant to this +fix. + +This fix computes a global vector of length 2, which can be accessed +by various :doc:`output commands `. 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 ` command. This fix is not invoked during +:doc:`energy minimization `. + +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 ` +doc page for more info. + +Voronoi compute must be enabled on build. See :doc:`compute voronoi/atom `. +A vaild voronoi command which returns neighboring atoms must be used +and referenced with the *voro* keyword. + +When this fix is used with a :doc:`hybrid pair style ` +system, only swaps between atom types of the same sub-style (or +combination of sub-styles) are permitted. + +If this fix is used with systems that do not have per-type masses +(e.g. atom style sphere), the ke flag must be set to off since the implemented +algorithm will not be able to re-scale velocity properly. + +Related commands +"""""""""""""""" + +:doc:`fix nvt `, :doc:`compute voronoi/atom ` +:doc:`delete_atoms `, :doc:`fix gcmc `, +:doc:`fix atom/swap `, :doc:`fix mol/swap `, :doc:`fix sgcmc ` + +Default +""""""" + +The option defaults are *ke* = yes, *diff* = no, *rates* = 1 for +all atom types. + +---------- + +.. Tavenner: + +**(Tavenner)** J Tavenner, M Mendelev, J Lawson, Computational Materials Science, 218, 111929 (2023). From 4004d263048f914dc93f3332cdbd2106ecbe9163 Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Thu, 19 Dec 2024 12:57:17 -0700 Subject: [PATCH 03/22] Spelling check fix --- src/MC/fix_neighbor_swap.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/MC/fix_neighbor_swap.cpp b/src/MC/fix_neighbor_swap.cpp index 52730ee6f5..9bbb777d92 100755 --- a/src/MC/fix_neighbor_swap.cpp +++ b/src/MC/fix_neighbor_swap.cpp @@ -103,7 +103,7 @@ FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) : options(narg-7,&arg[7]); - if (voro_flag != 1) error->all(FLERR,"Voro Compute Required for fix neighbor/swap command"); + if (voro_flag != 1) error->all(FLERR,"Voronoi compute required for fix neighbor/swap command"); // random number generator, same for all procs @@ -194,7 +194,7 @@ void FixNeighborSwap::options(int narg, char **arg) iarg++; } } else if (strcmp(arg[iarg],"voro") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix neighbor/swap command"); + if (iarg + 2 > narg) error->all(FLERR,"Illegal fix neighbor/swap command"); int icompute = modify->find_compute(utils::strdup(arg[iarg+1])); @@ -208,7 +208,7 @@ void FixNeighborSwap::options(int narg, char **arg) voro_flag = 1; iarg += 2; } else if (strcmp(arg[iarg],"diff") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix neighbor/swap command"); + 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; From c847ac1fd45259e4bec3618d2ee2153a213a77fe Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Thu, 19 Dec 2024 13:37:30 -0700 Subject: [PATCH 04/22] Removed trailing whitespace --- doc/src/fix_neighbor_swap.rst | 12 ++++++------ src/MC/fix_neighbor_swap.cpp | 30 +++++++++++++++--------------- src/MC/fix_neighbor_swap.h | 4 ++-- 3 files changed, 23 insertions(+), 23 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index a3eb1c6d82..e97b8bcbf2 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -70,11 +70,11 @@ as given by *diff type*. This command selects all atom types as acceptable swap centrally selected atom of type *type*. This includes the atom type specified by the diff keyword to account for self-diffusion hops of an atom type with itself. -Keyword *voro* is currently required, and is implemented as +Keyword *voro* is currently required, and is implemented as .. code-block:: LAMMPS voro compute-ID - + where *compute-ID* is the ID of a corresponding Voronoi computation with neighbor list, i.e. .. code-block:: LAMMPS @@ -89,17 +89,17 @@ be swapped at the initial step, i.e. for using *fix neighbor/swap* with *diff 2*. -The keyword *rates* can modify the swap rate for each swapped type by values -where the adjusted rates values are given in order of increasing atom type. +The keyword *rates* can modify the swap rate for each swapped type by values +where the adjusted rates values are given in order of increasing atom type. The number of rates provided must equal the number of atom types in the simulaton. In the third provided example above, a simulation is conducted with three atom types where the third atom type is the one sampled for attempted swaps. All three atom types are considered valid swaps, but atoms of type 1 will be selected three times as often as atoms of type 2. Conversely, atoms of type 3 are six times more likely to -be selected than atoms of type two and twice as likely as atoms of type 1. +be selected than atoms of type two and twice as likely as atoms of type 1. Finally, the *region* keyword is implemented as in other atomic fixes, where -the *region region-ID* command indicates that atom swaps only be considered in the area +the *region region-ID* command indicates that atom swaps only be considered in the area given by *region-ID*. If only atoms of certain groups are expected to be in this region, the corresponding compute voronoi command can be adjusted accordingly. diff --git a/src/MC/fix_neighbor_swap.cpp b/src/MC/fix_neighbor_swap.cpp index 9bbb777d92..1bd9cbbc02 100755 --- a/src/MC/fix_neighbor_swap.cpp +++ b/src/MC/fix_neighbor_swap.cpp @@ -195,7 +195,7 @@ void FixNeighborSwap::options(int narg, char **arg) } } else if (strcmp(arg[iarg],"voro") == 0) { if (iarg + 2 > narg) error->all(FLERR,"Illegal fix neighbor/swap command"); - + int icompute = modify->find_compute(utils::strdup(arg[iarg+1])); if (icompute < 0) error->all(FLERR,"Could not find neighbor compute ID"); @@ -481,7 +481,7 @@ int FixNeighborSwap::attempt_swap() } energy_stored = energy_after; - // utils::logmesg(lmp,"Swap accepted\n"); + // utils::logmesg(lmp,"Swap accepted\n"); // Swap atom groups if successful swap // int groupi = atom->mask[i]; @@ -498,7 +498,7 @@ int FixNeighborSwap::attempt_swap() // 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 - // utils::logmesg(lmp,"Swap not accepted\n"); + // utils::logmesg(lmp,"Swap not accepted\n"); if ( i >= 0 ) { atom->type[i] = itype; @@ -613,7 +613,7 @@ double FixNeighborSwap::get_distance(double* i, double* j) // r_x = i[0] - j[0]; // r_y = i[1] - j[1]; // r_z = i[2] - j[2]; - + // // Domain::minimum_image(r_x, r_y, r_z); // double r = sqrt(pow(r_x, 2.) + // pow(r_y, 2.) + @@ -648,7 +648,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) 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"); @@ -656,7 +656,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // Compute voronoi and access neighbor list c_voro->compute_local(); - + voro_neighbor_list = c_voro->array_local; njswap_local = 0; local_probability = 0.0; @@ -664,7 +664,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // utils::logmesg(lmp,"Searching for atom {}\n",id_center); for (int n = 0; n < c_voro->size_local_rows; n++){ - + int temp_j_id = -1; int temp_j = -1; @@ -708,13 +708,13 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) if (atom->mask[temp_j] & groupbit) { if (diff_flag) { // Calculate distance from i to each j, adjust probability of selection - + // Get distance if own centr 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) && @@ -743,7 +743,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) 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) && @@ -776,7 +776,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) if ( i_center >= 0 ){ r = get_distance(x[temp_j], x[i_center]); } - + // Get local id of ghost center atoms // if ( i_center < 0 ){ // utils::logmesg(lmp,"Initial distance {}\n", r); @@ -788,7 +788,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) } } // } - + // utils::logmesg(lmp,"Final distance {}\n", r); if (rates_flag) { @@ -797,7 +797,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); } 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++; @@ -810,7 +810,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) 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) && @@ -823,7 +823,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) local_swap_probability[njswap_local] = rate_list[type[temp_j] - 1]*exp(-pow((r/3.), 2.)); } else { local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); - } + } local_probability += local_swap_probability[njswap_local]; local_swap_type_list[njswap_local] = jswaptype; diff --git a/src/MC/fix_neighbor_swap.h b/src/MC/fix_neighbor_swap.h index 1413de2ef8..e53eb9c1c3 100755 --- a/src/MC/fix_neighbor_swap.h +++ b/src/MC/fix_neighbor_swap.h @@ -63,7 +63,7 @@ class FixNeighborSwap : public Fix { double nswap_attempts; double nswap_successes; - + bool unequal_cutoffs; int atom_swap_nmax; @@ -79,7 +79,7 @@ class FixNeighborSwap : public Fix { 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 RanPark *random_unequal; From 49b2b978e67cd8df42dfddfa7ad5faaa57443fdc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2024 19:20:51 -0500 Subject: [PATCH 05/22] build system updates to skip fix neighbor/swap if VORONOI package is not installed --- cmake/Modules/Packages/MC.cmake | 10 ++++++++++ src/.gitignore | 2 ++ src/MC/Install.sh | 4 ++-- 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/cmake/Modules/Packages/MC.cmake b/cmake/Modules/Packages/MC.cmake index f162254558..2a72a895cf 100644 --- a/cmake/Modules/Packages/MC.cmake +++ b/cmake/Modules/Packages/MC.cmake @@ -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() diff --git a/src/.gitignore b/src/.gitignore index f0554e3bfe..0f0e3759c9 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -912,6 +912,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 diff --git a/src/MC/Install.sh b/src/MC/Install.sh index 9a0095fd34..efe6b7c07b 100755 --- a/src/MC/Install.sh +++ b/src/MC/Install.sh @@ -51,8 +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 -action fix_neighbor_swap.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 From b05172fe9687ba25b65efbb6cd6c47c980d360b1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2024 19:21:35 -0500 Subject: [PATCH 06/22] spelling fixes and documentation formatting corrections --- doc/src/fix.rst | 2 +- doc/src/fix_neighbor_swap.rst | 90 ++++++++++++--------- doc/utils/sphinx-config/false_positives.txt | 1 + 3 files changed, 52 insertions(+), 41 deletions(-) diff --git a/doc/src/fix.rst b/doc/src/fix.rst index a4a9ac4a86..ae02118bec 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -289,7 +289,7 @@ accelerated styles exist. * :doc:`mvv/tdpd ` - constant temperature DPD using the modified velocity-Verlet algorithm * :doc:`neb ` - nudged elastic band (NEB) spring forces * :doc:`neb/spin ` - nudged elastic band (NEB) spring forces for spins -* :doc:`neighbor/swap ` - kinetic Monte Calo atom swapping +* :doc:`neighbor/swap ` - kinetic Monte Carlo (kMC) atom swapping * :doc:`nonaffine/displacement ` - calculate nonaffine displacement of atoms * :doc:`nph ` - constant NPH time integration via Nose/Hoover * :doc:`nph/asphere ` - NPH for aspherical particles diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index e97b8bcbf2..898957c56c 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -1,7 +1,7 @@ .. index:: fix neighbor/swap fix neighbor/swap command -===================== +========================= Syntax """""" @@ -50,61 +50,69 @@ Neighboring atoms are selected using a Voronoi tesselation approach. This implementation is as described in :ref:`(Tavenner) `. The fix is called every *N* timesteps and attempts *X* swaps. The system -is initialized with a random seed, using a temperature *T* for evaluating -the MC energy swaps. The distance-based probability is weighted according -to *R* which sets the radius :math:`r_0` for the weighting +is initialized with a random seed, using a temperature *T* for +evaluating the MC energy swaps. The distance-based probability is +weighted according to *R* which sets the radius :math:`r_0` for the +weighting .. math:: p_{ij} = e^{(\frac{r_{ij}}{r_0})^2} -where :math:`p_{ij}` is the probability of selecting atoms :math:`i` and :math:`j` for an -evaluated swap. +where :math:`p_{ij}` is the probability of selecting atoms :math:`i` and +:math:`j` for an evaluated swap. -The keyword *types* is submitted with two or more atom types as the value. -Atoms of the first atom type are swapped with valid neighbors of all the remaining -atom types. +The keyword *types* is submitted with two or more atom types as the +value. Atoms of the first atom type are swapped with valid neighbors of +all the remaining atom types. -The keyword *diff* is used for implementation of simulated diffusion of a given atom type -as given by *diff type*. This command selects all atom types as acceptable swap types to a -centrally selected atom of type *type*. This includes the atom type specified by the diff -keyword to account for self-diffusion hops of an atom type with itself. +The keyword *diff* is used for implementation of simulated diffusion of +a given atom type as given by *diff type*. This command selects all atom +types as acceptable swap types to a centrally selected atom of type +*type*. This includes the atom type specified by the diff keyword to +account for self-diffusion hops of an atom type with itself. Keyword *voro* is currently required, and is implemented as .. code-block:: LAMMPS - voro compute-ID -where *compute-ID* is the ID of a corresponding Voronoi computation with neighbor list, i.e. + voro compute-ID + +where *compute-ID* is the ID of a corresponding Voronoi computation with +neighbor list, i.e. .. code-block:: LAMMPS compute compute-ID group-ID voronoi/atom neighbors yes -The group selected for computing *voro* should correspond to all the potential atoms to -be swapped at the initial step, i.e. +The group selected for computing *voro* should correspond to all the +potential atoms to be swapped at the initial step, i.e. .. code-block:: LAMMPS - group group-ID type 2 + + group group-ID type 2 for using *fix neighbor/swap* with *diff 2*. -The keyword *rates* can modify the swap rate for each swapped type by values -where the adjusted rates values are given in order of increasing atom type. -The number of rates provided must equal the number of atom types in the simulaton. -In the third provided example above, a simulation is conducted with three atom types -where the third atom type is the one sampled for attempted swaps. All three atom -types are considered valid swaps, but atoms of type 1 will be selected three times -as often as atoms of type 2. Conversely, atoms of type 3 are six times more likely to -be selected than atoms of type two and twice as likely as atoms of type 1. +The keyword *rates* can modify the swap rate for each swapped type by +values where the adjusted rates values are given in order of increasing +atom type. The number of rates provided must equal the number of atom +types in the simulation. In the third provided example above, a +simulation is conducted with three atom types where the third atom type +is the one sampled for attempted swaps. All three atom types are +considered valid swaps, but atoms of type 1 will be selected three times +as often as atoms of type 2. Conversely, atoms of type 3 are six times +more likely to be selected than atoms of type two and twice as likely as +atoms of type 1. -Finally, the *region* keyword is implemented as in other atomic fixes, where -the *region region-ID* command indicates that atom swaps only be considered in the area -given by *region-ID*. If only atoms of certain groups are expected to be in this region, -the corresponding compute voronoi command can be adjusted accordingly. +Finally, the *region* keyword is implemented as in other atomic fixes, +where the *region region-ID* command indicates that atom swaps only be +considered in the area given by *region-ID*. If only atoms of certain +groups are expected to be in this region, the corresponding compute +voronoi command can be adjusted accordingly. Restart, fix_modify, output, run start/stop, minimize info -""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" +"""""""""""""""""""""""""""""""""""""""""""""""""""""""""" This fix writes the state of the fix to :doc:`binary restart files `. This includes information about the random number generator @@ -135,10 +143,11 @@ 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 ` -doc page for more info. +doc page for more info. Also this fix requires that the +:ref:`VORONOI package ` is installed, otherwise the fix +will not be compiled. -Voronoi compute must be enabled on build. See :doc:`compute voronoi/atom `. -A vaild voronoi command which returns neighboring atoms must be used +A valid voronoi command which returns neighboring atoms must be used and referenced with the *voro* keyword. When this fix is used with a :doc:`hybrid pair style ` @@ -146,24 +155,25 @@ system, only swaps between atom types of the same sub-style (or combination of sub-styles) are permitted. If this fix is used with systems that do not have per-type masses -(e.g. atom style sphere), the ke flag must be set to off since the implemented -algorithm will not be able to re-scale velocity properly. +(e.g. atom style sphere), the ke flag must be set to off since the +implemented algorithm will not be able to re-scale velocity properly. Related commands """""""""""""""" :doc:`fix nvt `, :doc:`compute voronoi/atom ` :doc:`delete_atoms `, :doc:`fix gcmc `, -:doc:`fix atom/swap `, :doc:`fix mol/swap `, :doc:`fix sgcmc ` +:doc:`fix atom/swap `, :doc:`fix mol/swap `, +:doc:`fix sgcmc ` Default """"""" -The option defaults are *ke* = yes, *diff* = no, *rates* = 1 for -all atom types. +The option defaults are *ke* = yes, *diff* = no, *rates* = 1 for all +atom types. ---------- -.. Tavenner: +.. _Tavenner: **(Tavenner)** J Tavenner, M Mendelev, J Lawson, Computational Materials Science, 218, 111929 (2023). diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 67db18a17d..4d2b6f5d68 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1812,6 +1812,7 @@ Kloss Kloza kmax Kmax +kMC KMP kmu Knizhnik From 8a89d2fcf6f19b932abc20c440ea2016ac43fe79 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2024 19:22:02 -0500 Subject: [PATCH 07/22] no more error lists in headers --- src/MC/fix_neighbor_swap.h | 43 -------------------------------------- 1 file changed, 43 deletions(-) diff --git a/src/MC/fix_neighbor_swap.h b/src/MC/fix_neighbor_swap.h index e53eb9c1c3..df75505a6d 100755 --- a/src/MC/fix_neighbor_swap.h +++ b/src/MC/fix_neighbor_swap.h @@ -101,46 +101,3 @@ class FixNeighborSwap : public Fix { #endif #endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Region ID for fix neighbor/swap does not exist - -Self-explanatory. - -E: Must specify at least 2 types in fix neighbor/swap command - -Self-explanatory. - -E: Need nswaptypes mu values in fix neighbor/swap command - -Self-explanatory. - -E: Invalid atom type in fix neighbor/swap command - -The atom type specified in the neighbor/swap command does not exist. - -E: All atoms of a swapped type must have the same charge. - -Self-explanatory. - -E: At least one atom of each swapped type must be present to define charges. - -Self-explanatory. - -E: All atoms of a swapped type must have same charge. - -Self-explanatory. - -E: Cannot do neighbor/swap on atoms in atom_modify first group - -This is a restriction due to the way atoms are organized in a list to -enable the atom_modify first command. - -*/ From 3c0d4c8e145140fac912619da6e15bebadcddaff Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2024 19:36:36 -0500 Subject: [PATCH 08/22] add versionadded tag --- doc/src/fix_neighbor_swap.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 898957c56c..3c97400d3f 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -44,6 +44,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Computes MC evaluations to enable kinetic Monte Carlo (kMC) behavior during MD simulation through only allowing neighboring atom swaps. Neighboring atoms are selected using a Voronoi tesselation approach. This From 42b6308e2669ec7d8bcfdbc6c56ccf4505abf334 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2024 19:37:22 -0500 Subject: [PATCH 09/22] remove commented out debug code, replace pow(x,2.0) with MathSpecial::square() --- src/MC/fix_neighbor_swap.cpp | 174 +++++++++-------------------------- 1 file changed, 44 insertions(+), 130 deletions(-) diff --git a/src/MC/fix_neighbor_swap.cpp b/src/MC/fix_neighbor_swap.cpp index 1bd9cbbc02..703207783a 100755 --- a/src/MC/fix_neighbor_swap.cpp +++ b/src/MC/fix_neighbor_swap.cpp @@ -31,6 +31,7 @@ #include "group.h" #include "improper.h" #include "kspace.h" +#include "math_special.h" #include "memory.h" #include "modify.h" #include "neighbor.h" @@ -68,9 +69,7 @@ FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) : sqrt_mass_ratio(nullptr), local_swap_iatom_list(nullptr), random_equal(nullptr), random_unequal(nullptr), c_pe(nullptr) { - if (narg < 10) error->all(FLERR,"Illegal fix neighbor/swap command"); - - // utils::logmesg(lmp,"Starting neighbor/swap\n"); + if (narg < 10) utils::missing_cmd_args(FLERR,"fix neighbor/swap", error); dynamic_group_allow = 1; @@ -88,11 +87,10 @@ FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) : seed = utils::inumeric(FLERR,arg[5],false,lmp); double temperature = utils::numeric(FLERR,arg[6],false,lmp); - if (nevery <= 0) error->all(FLERR,"Illegal fix neighbor/swap command"); - if (ncycles < 0) error->all(FLERR,"Illegal fix neighbor/swap command"); - // if (nevery < ncycles) error->all(FLERR,"Illegal fix neighbor/swap command"); - if (seed <= 0) error->all(FLERR,"Illegal fix neighbor/swap command"); - if (temperature <= 0.0) error->all(FLERR,"Illegal fix neighbor/swap command"); + 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); @@ -170,7 +168,6 @@ void FixNeighborSwap::options(int narg, char **arg) int iarg = 0; while (iarg < narg) { - // utils::logmesg(lmp,"Parsing Argument {}\n", arg[iarg]); 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]); @@ -216,7 +213,6 @@ void FixNeighborSwap::options(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"rates") == 0) { if (iarg+atom->ntypes >= narg) error->all(FLERR,"Illegal fix neighbor/swap command"); - // utils::logmesg(lmp,"Reading rates\n"); iarg++; int i = 0; while (iarg < narg) { @@ -246,7 +242,6 @@ int FixNeighborSwap::setmask() void FixNeighborSwap::init() { - // utils::logmesg(lmp,"Initializing neighbor/swap\n"); c_pe = modify->get_compute_by_id("thermo_pe"); int *type = atom->type; @@ -321,10 +316,7 @@ void FixNeighborSwap::init() 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"); - - // utils::logmesg(lmp,"Done Init\n"); + if (flagall) error->all(FLERR,"Cannot do neighbor/swap on atoms in atom_modify first group"); } } @@ -355,15 +347,10 @@ void FixNeighborSwap::pre_exchange() // attempt Ncycle atom swaps - // utils::logmesg(lmp,"Attempting Swap\n"); - int nsuccess = 0; update_iswap_atoms_list(); - // utils::logmesg(lmp,"Updated Atom i List\n"); for (int i = 0; i < ncycles; i++) nsuccess += attempt_swap(); - // utils::logmesg(lmp,"\n"); - // udpate MC stats nswap_attempts += ncycles; @@ -398,36 +385,17 @@ int FixNeighborSwap::attempt_swap() // build nearest-neighbor list based on atom i build_i_neighbor_list(i); - // utils::logmesg(lmp,"Built neighbor list with {} atoms\n",njswap); - if ( njswap <= 0 ) { - // utils::logmesg(lmp,"No valid neighbors found\n"); - return 0; - } + 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); - // utils::logmesg(lmp,"Selected swap neighbor j\n"); int itype = type_list[0]; int jtype = jtype_selected; - // utils::logmesg(lmp,"Selected atom j type {}\n",jtype); - - // utils::logmesg(lmp,"Swapping local atoms {} and ",i); - // utils::logmesg(lmp,"{}\n",j); - - // if ( i >= 0 || j >= 0) { - // utils::logmesg(lmp,"Selected atom j type {}\n",jtype); - // utils::logmesg(lmp,"Swapping atoms {} and ",id[i]); - // utils::logmesg(lmp,"{}\n",id[j]); - // } - - // utils::logmesg(lmp,"Atom type {} and ",itype); - // utils::logmesg(lmp,"{}\n",jtype); // Accept swap if types are equal, no change to system - if ( itype == jtype ){ - // utils::logmesg(lmp,"Atoms have same type, no processing needed\n"); + if ( itype == jtype ) { return 1; } @@ -480,16 +448,6 @@ int FixNeighborSwap::attempt_swap() } } energy_stored = energy_after; - - // utils::logmesg(lmp,"Swap accepted\n"); - // Swap atom groups if successful swap - - // int groupi = atom->mask[i]; - // int groupj = atom->mask[j]; - - // atom->mask[i] = group->bitmask[groupj]; - // atom->mask[j] = group->bitmask[groupi]; - return 1; } @@ -498,8 +456,6 @@ int FixNeighborSwap::attempt_swap() // 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 - // utils::logmesg(lmp,"Swap not accepted\n"); - if ( i >= 0 ) { atom->type[i] = itype; if (atom->q_flag) atom->q[i] = qtype[0]; @@ -529,7 +485,7 @@ double FixNeighborSwap::energy_full() 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->improper) force->improper->compute(eflag,vflag); } if (force->kspace) force->kspace->compute(eflag,vflag); @@ -576,21 +532,17 @@ int FixNeighborSwap::pick_j_swap_neighbor(int i) // Generate random double from 0 to maximum global probability double selected_prob = static_cast (global_probability*random_equal->uniform()); - // utils::logmesg(lmp,"Local Probability {}\n",local_probability); - // Find which local swap atom corresponds to probability if ((selected_prob >= prev_probability) && - (selected_prob < prev_probability + local_probability)){ - // utils::logmesg(lmp,"Selected Probability {}\n",selected_prob); + (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]){ + 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]; - // utils::logmesg(lmp,"Selecting Atom j type {}\n",jtype_selected_local); MPI_Allreduce(&jtype_selected_local,&jtype_selected,1,MPI_INT,MPI_MAX,world); return j; } @@ -598,7 +550,6 @@ int FixNeighborSwap::pick_j_swap_neighbor(int i) error->all(FLERR,"Did not select local neighbor swap atom"); } - // utils::logmesg(lmp,"Swap atom not found on processor\n"); MPI_Allreduce(&jtype_selected_local,&jtype_selected,1,MPI_INT,MPI_MAX,world); return j; } @@ -608,22 +559,9 @@ int FixNeighborSwap::pick_j_swap_neighbor(int i) double FixNeighborSwap::get_distance(double* i, double* j) { - // double r_x, r_y, r_z; - - // r_x = i[0] - j[0]; - // r_y = i[1] - j[1]; - // r_z = i[2] - j[2]; - - // // Domain::minimum_image(r_x, r_y, r_z); - // double r = sqrt(pow(r_x, 2.) + - // pow(r_y, 2.) + - // pow(r_z, 2.)); - - double r = sqrt(pow((i[0] - j[0]), 2.) + - pow((i[1] - j[1]), 2.) + - pow((i[2] - j[2]), 2.)); - - + double r = sqrt(MathSpecial::square((i[0] - j[0])) + + MathSpecial::square((i[1] - j[1])) + + MathSpecial::square((i[2] - j[2]))); return r; } @@ -661,19 +599,17 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) njswap_local = 0; local_probability = 0.0; - // utils::logmesg(lmp,"Searching for atom {}\n",id_center); - - for (int n = 0; n < c_voro->size_local_rows; n++){ + 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 ){ + 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 ) ){ + ( i_center < 0 ) ) { temp_j_id = voro_neighbor_list[n][0]; temp_j = -1; } else { @@ -681,10 +617,9 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) } // Find which local atom corresponds to neighbor - for (int j = 0; j < nlocal; j++){ - if ( temp_j_id == id[j] ){ + for (int j = 0; j < nlocal; j++) { + if ( temp_j_id == id[j] ) { temp_j = j; - // utils::logmesg(lmp,"Found neighbor {}\n",id[j]); break; } } @@ -692,17 +627,6 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // If temp_j not on this processor, skip if ( temp_j < 0 ) continue; - // // If atom is already in local list, skip - // bool inlist = false; - // for ( int j = 0; j < njswap_local; j++ ){ - // if ( temp_j == local_swap_neighbor_list[j] ){ - // utils::logmesg(lmp,"Skipping atom {}, pair already counted\n",id[temp_j]); - // inlist = true; - // break; - // } - // } - // if (inlist) 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) { @@ -711,51 +635,49 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // Get distance if own centr atom double r = INFINITY; - if ( i_center >= 0 ){ + 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++){ + for (int i=nlocal; i < nlocal+nghost; i++) { if ( (id[i] == id_center) && - (get_distance(x[temp_j], x[i]) < r) ){ + (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(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = rate_list[type[temp_j] - 1]*exp(-MathSpecial::square(r/3.0)); } else { - local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); } local_probability += local_swap_probability[njswap_local]; - - // utils::logmesg(lmp,"Adding Atom j type {}\n",type[temp_j]); 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++){ + 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 ){ + 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++){ + for (int i=nlocal; i < nlocal+nghost; i++) { if ( (id[i] == id_center) && - (get_distance(x[temp_j], x[i]) < r) ){ + (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(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = rate_list[type[temp_j] - 1]*exp(-MathSpecial::square(r/3.0)); } else { - local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); } local_probability += local_swap_probability[njswap_local]; @@ -773,28 +695,20 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // Calculate distance from i to each j, adjust probability of selection // Get distance if own center atom double r = INFINITY; - if ( i_center >= 0 ){ + if ( i_center >= 0 ) { r = get_distance(x[temp_j], x[i_center]); } // Get local id of ghost center atoms - // if ( i_center < 0 ){ - // utils::logmesg(lmp,"Initial distance {}\n", r); - 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]); - // utils::logmesg(lmp,"Updated distance {}\n", r); - } + 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]); } - // } - - // utils::logmesg(lmp,"Final distance {}\n", r); if (rates_flag) { - local_swap_probability[njswap_local] = rate_list[type[temp_j] - 1]*exp(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = rate_list[type[temp_j] - 1]*exp(-MathSpecial::square(r/3.0)); } else{ - local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); } local_probability += local_swap_probability[njswap_local]; @@ -802,27 +716,27 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) local_swap_neighbor_list[njswap_local] = temp_j; njswap_local++; } else { - for (int jswaptype = 1; jswaptype < nswaptypes; jswaptype++){ + 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 ){ + 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++){ + for (int i=nlocal; i < nlocal+nghost; i++) { if ( (id[i] == id_center) && - (get_distance(x[temp_j], x[i]) < r) ){ + (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(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = rate_list[type[temp_j] - 1]*exp(-MathSpecial::square(r/3.0)); } else { - local_swap_probability[njswap_local] = exp(-pow((r/3.), 2.)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); } local_probability += local_swap_probability[njswap_local]; From 60b0ef68a60e1dccb5a7341f7c748985bd27f104 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2024 19:39:19 -0500 Subject: [PATCH 10/22] enable and apply clang-format --- src/MC/fix_neighbor_swap.cpp | 360 +++++++++++++++++------------------ src/MC/fix_neighbor_swap.h | 21 +- 2 files changed, 185 insertions(+), 196 deletions(-) diff --git a/src/MC/fix_neighbor_swap.cpp b/src/MC/fix_neighbor_swap.cpp index 703207783a..95e74bee22 100755 --- a/src/MC/fix_neighbor_swap.cpp +++ b/src/MC/fix_neighbor_swap.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -23,6 +22,7 @@ #include "bond.h" #include "comm.h" #include "compute.h" +#include "compute_voronoi_atom.h" #include "dihedral.h" #include "domain.h" #include "error.h" @@ -39,37 +39,36 @@ #include "random_park.h" #include "region.h" #include "update.h" -#include "compute_voronoi_atom.h" -#include #include #include +#include #include 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"; + "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), random_unequal(nullptr), c_pe(nullptr) + 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), random_unequal(nullptr), c_pe(nullptr) { - if (narg < 10) utils::missing_cmd_args(FLERR,"fix neighbor/swap", error); + if (narg < 10) utils::missing_cmd_args(FLERR, "fix neighbor/swap", error); dynamic_group_allow = 1; @@ -82,34 +81,34 @@ FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) : // 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); + 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); - 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"); + 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"); + 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-7,&arg[7]); + options(narg - 7, &arg[7]); - if (voro_flag != 1) error->all(FLERR,"Voronoi compute required for fix neighbor/swap command"); + if (voro_flag != 1) error->all(FLERR, "Voronoi compute required for fix neighbor/swap command"); // random number generator, same for all procs - random_equal = new RanPark(lmp,seed); + random_equal = new RanPark(lmp, seed); // random number generator, not the same for all procs - random_unequal = new RanPark(lmp,seed); + random_unequal = new RanPark(lmp, seed); // set up reneighboring @@ -130,9 +129,10 @@ FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) : // set comm size needed by this Fix - if (atom->q_flag) comm_forward = 2; - else comm_forward = 1; - + if (atom->q_flag) + comm_forward = 2; + else + comm_forward = 1; } /* ---------------------------------------------------------------------- */ @@ -158,7 +158,7 @@ FixNeighborSwap::~FixNeighborSwap() void FixNeighborSwap::options(int narg, char **arg) { - if (narg < 0) error->all(FLERR,"Illegal fix neighbor/swap command\n"); + if (narg < 0) error->all(FLERR, "Illegal fix neighbor/swap command\n"); ke_flag = 1; diff_flag = 0; @@ -168,64 +168,64 @@ void FixNeighborSwap::options(int narg, char **arg) 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"); + 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); + } 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"); + } 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); + 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],"voro") == 0) { - if (iarg + 2 > narg) error->all(FLERR,"Illegal fix neighbor/swap command"); + } else if (strcmp(arg[iarg], "voro") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix neighbor/swap command"); - int icompute = modify->find_compute(utils::strdup(arg[iarg+1])); + int icompute = modify->find_compute(utils::strdup(arg[iarg + 1])); - if (icompute < 0) error->all(FLERR,"Could not find neighbor compute ID"); + 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"); + 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"); + error->all(FLERR, "Neighbor compute does not give i, j, size as expected"); voro_flag = 1; iarg += 2; - } 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); + } 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"); + } 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); + 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"); + if (i != atom->ntypes) error->all(FLERR, "Illegal fix neighbor/swap command"); + } else + error->all(FLERR, "Illegal fix neighbor/swap command"); } } @@ -246,7 +246,8 @@ void FixNeighborSwap::init() int *type = atom->type; - if (nswaptypes < 2 && !diff_flag) error->all(FLERR,"Must specify at least 2 types in fix neighbor/swap command"); + 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 @@ -257,12 +258,12 @@ void FixNeighborSwap::init() 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"); + 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"); + 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++) { @@ -272,24 +273,27 @@ void FixNeighborSwap::init() 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."); + 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."); + 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); + 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."); + 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"); + 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]); + 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 @@ -314,9 +318,9 @@ void FixNeighborSwap::init() if ((mask[i] == groupbit) && (mask[i] && firstgroupbit)) flag = 1; int flagall; - MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); + 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"); + if (flagall) error->all(FLERR, "Cannot do neighbor/swap on atoms in atom_modify first group"); } } @@ -336,7 +340,7 @@ void FixNeighborSwap::pre_exchange() domain->pbc(); comm->exchange(); comm->borders(); - if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(1); @@ -369,7 +373,7 @@ int FixNeighborSwap::attempt_swap() // int nlocal = atom->nlocal; tagint *id = atom->tag; - if ( niswap == 0 ) return 0; + if (niswap == 0) return 0; // pre-swap energy @@ -385,7 +389,7 @@ int FixNeighborSwap::attempt_swap() // build nearest-neighbor list based on atom i build_i_neighbor_list(i); - if ( njswap <= 0 ) return 0; + if (njswap <= 0) return 0; // pick a neighbor atom j based on i neighbor list jtype_selected = -1; @@ -395,16 +399,14 @@ int FixNeighborSwap::attempt_swap() int jtype = jtype_selected; // Accept swap if types are equal, no change to system - if ( itype == jtype ) { - return 1; - } + if (itype == jtype) { return 1; } // swap their properties - if ( i >= 0 ) { + if (i >= 0) { atom->type[i] = jtype; if (atom->q_flag) atom->q[i] = qtype[jtype_selected]; } - if ( j >= 0 ) { + if (j >= 0) { atom->type[j] = itype; if (atom->q_flag) atom->q[j] = qtype[0]; } @@ -418,7 +420,7 @@ int FixNeighborSwap::attempt_swap() domain->pbc(); comm->exchange(); comm->borders(); - if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(1); } else { @@ -432,16 +434,15 @@ int FixNeighborSwap::attempt_swap() // if swap accepted, return 1 // if ke_flag, rescale atom velocities - if (random_equal->uniform() < - exp(beta*(energy_before - energy_after))) { + if (random_equal->uniform() < exp(beta * (energy_before - energy_after))) { update_iswap_atoms_list(); if (ke_flag) { - if ( i >= 0 ) { + 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 ) { + 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]; @@ -456,11 +457,11 @@ int FixNeighborSwap::attempt_swap() // 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 ) { + if (i >= 0) { atom->type[i] = itype; if (atom->q_flag) atom->q[i] = qtype[0]; } - if ( j >= 0 ) { + if (j >= 0) { atom->type[j] = jtype; if (atom->q_flag) atom->q[j] = qtype[jtype_selected]; } @@ -479,16 +480,16 @@ double FixNeighborSwap::energy_full() if (modify->n_pre_force) modify->pre_force(vflag); - if (force->pair) force->pair->compute(eflag,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->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 (force->kspace) force->kspace->compute(eflag, vflag); if (modify->n_post_force_any) modify->post_force(vflag); @@ -507,15 +508,14 @@ int FixNeighborSwap::pick_i_swap_atom() int id_center_local = -1; int i = -1; - int iwhichglobal = static_cast (niswap*random_equal->uniform()); - if ((iwhichglobal >= niswap_before) && - (iwhichglobal < niswap_before + niswap_local)) { + int iwhichglobal = static_cast(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); + 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); + MPI_Allreduce(&id[i], &id_center, 1, MPI_INT, MPI_MAX, world); } return i; @@ -530,7 +530,7 @@ int FixNeighborSwap::pick_j_swap_neighbor(int i) int jtype_selected_local = -1; // Generate random double from 0 to maximum global probability - double selected_prob = static_cast (global_probability*random_equal->uniform()); + double selected_prob = static_cast(global_probability * random_equal->uniform()); // Find which local swap atom corresponds to probability if ((selected_prob >= prev_probability) && @@ -539,28 +539,26 @@ int FixNeighborSwap::pick_j_swap_neighbor(int i) for (int n = 0; n < njswap_local; n++) { if (search_prob > local_swap_probability[n]) { search_prob -= local_swap_probability[n]; - } - else{ + } 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); + 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"); + error->all(FLERR, "Did not select local neighbor swap atom"); } - MPI_Allreduce(&jtype_selected_local,&jtype_selected,1,MPI_INT,MPI_MAX,world); + MPI_Allreduce(&jtype_selected_local, &jtype_selected, 1, MPI_INT, MPI_MAX, world); return j; } /* ---------------------------------------------------------------------- ------------------------------------------------------------------------- */ -double FixNeighborSwap::get_distance(double* i, double* j) +double FixNeighborSwap::get_distance(double *i, double *j) { - double r = sqrt(MathSpecial::square((i[0] - j[0])) + - MathSpecial::square((i[1] - j[1])) + + double r = sqrt(MathSpecial::square((i[0] - j[0])) + MathSpecial::square((i[1] - j[1])) + MathSpecial::square((i[2] - j[2]))); return r; } @@ -580,16 +578,16 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) 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"); + 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"); + 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"); + local_swap_type_list = + (int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_type_list"); // Compute voronoi and access neighbor list @@ -605,11 +603,10 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) int temp_j = -1; // Find local voronoi entry with selected central atom - if ( (int)voro_neighbor_list[n][0] == id_center ) { + 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 ) ) { + } 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 { @@ -618,39 +615,37 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // Find which local atom corresponds to neighbor for (int j = 0; j < nlocal; j++) { - if ( temp_j_id == id[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 (temp_j < 0) continue; if (region) { - if (region->match(x[temp_j][0],x[temp_j][1],x[temp_j][2]) == 1) { + 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 centr atom double r = INFINITY; - if ( i_center >= 0 ) { - double r = get_distance(x[temp_j], x[i_center]); - } + 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]); + 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/3.0)); + local_swap_probability[njswap_local] = + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); } else { - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); } local_probability += local_swap_probability[njswap_local]; local_swap_type_list[njswap_local] = type[temp_j]; @@ -662,22 +657,20 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // 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]); - } + 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]); + 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/3.0)); + local_swap_probability[njswap_local] = + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); } else { - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); } local_probability += local_swap_probability[njswap_local]; @@ -695,20 +688,19 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // 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]); - } + 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) ) + 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/3.0)); - } else{ - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); + local_swap_probability[njswap_local] = + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); + } else { + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); } local_probability += local_swap_probability[njswap_local]; @@ -721,22 +713,20 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) // 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]); - } + 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]); + 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/3.0)); + local_swap_probability[njswap_local] = + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); } else { - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r/3.0)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); } local_probability += local_swap_probability[njswap_local]; @@ -750,12 +740,12 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) } } - MPI_Allreduce(&njswap_local,&njswap,1,MPI_INT,MPI_SUM,world); - MPI_Scan(&njswap_local,&njswap_before,1,MPI_INT,MPI_SUM,world); + 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); + 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; } @@ -772,8 +762,8 @@ void FixNeighborSwap::update_iswap_atoms_list() 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"); + local_swap_iatom_list = + (int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_iatom_list"); } niswap_local = 0; @@ -781,9 +771,9 @@ void FixNeighborSwap::update_iswap_atoms_list() if (region) { for (int i = 0; i < nlocal; i++) { - if (region->match(x[i][0],x[i][1],x[i][2]) == 1) { + if (region->match(x[i][0], x[i][1], x[i][2]) == 1) { if (atom->mask[i] & groupbit) { - if (type[i] == type_list[0]) { + if (type[i] == type_list[0]) { local_swap_iatom_list[niswap_local] = i; niswap_local++; } @@ -794,7 +784,7 @@ void FixNeighborSwap::update_iswap_atoms_list() } else { for (int i = 0; i < nlocal; i++) { if (atom->mask[i] & groupbit) { - if (type[i] == type_list[0]) { + if (type[i] == type_list[0]) { local_swap_iatom_list[niswap_local] = i; niswap_local++; } @@ -802,16 +792,17 @@ void FixNeighborSwap::update_iswap_atoms_list() } } - MPI_Allreduce(&niswap_local,&niswap,1,MPI_INT,MPI_SUM,world); - MPI_Scan(&niswap_local,&niswap_before,1,MPI_INT,MPI_SUM,world); + 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 FixNeighborSwap::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { - int i,j,m; + int i, j, m; int *type = atom->type; double *q = atom->q; @@ -838,7 +829,7 @@ int FixNeighborSwap::pack_forward_comm(int n, int *list, double *buf, int /*pbc_ void FixNeighborSwap::unpack_forward_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; int *type = atom->type; double *q = atom->q; @@ -848,12 +839,11 @@ void FixNeighborSwap::unpack_forward_comm(int n, int first, double *buf) if (atom->q_flag) { for (i = first; i < last; i++) { - type[i] = static_cast (buf[m++]); + type[i] = static_cast(buf[m++]); q[i] = buf[m++]; } } else { - for (i = first; i < last; i++) - type[i] = static_cast (buf[m++]); + for (i = first; i < last; i++) type[i] = static_cast(buf[m++]); } } @@ -874,7 +864,7 @@ double FixNeighborSwap::compute_vector(int n) double FixNeighborSwap::memory_usage() { - double bytes = (double)atom_swap_nmax * sizeof(int); + double bytes = (double) atom_swap_nmax * sizeof(int); return bytes; } @@ -895,8 +885,8 @@ void FixNeighborSwap::write_restart(FILE *fp) if (comm->me == 0) { int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); + fwrite(&size, sizeof(int), 1, fp); + fwrite(list, sizeof(double), n, fp); } } @@ -909,10 +899,10 @@ void FixNeighborSwap::restart(char *buf) int n = 0; double *list = (double *) buf; - seed = static_cast (list[n++]); + seed = static_cast(list[n++]); random_equal->reset(seed); - seed = static_cast (list[n++]); + seed = static_cast(list[n++]); random_unequal->reset(seed); next_reneighbor = (bigint) ubuf(list[n++]).i; @@ -922,5 +912,5 @@ void FixNeighborSwap::restart(char *buf) 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"); + error->all(FLERR, "Must not reset timestep when restarting fix neighbor/swap"); } diff --git a/src/MC/fix_neighbor_swap.h b/src/MC/fix_neighbor_swap.h index df75505a6d..826ab96a65 100755 --- a/src/MC/fix_neighbor_swap.h +++ b/src/MC/fix_neighbor_swap.h @@ -40,17 +40,17 @@ class FixNeighborSwap : public Fix { 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 voro_flag; // yes = use given voronoi calculation, no = use internal voronoi calculation + 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 voro_flag; // yes = use given voronoi calculation, no = use internal voronoi calculation 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 + class Region *region; // swap region + char *idregion; // swap region id int nswaptypes; int jtype_selected; @@ -68,19 +68,18 @@ class FixNeighborSwap : public Fix { int atom_swap_nmax; double beta; - 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 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 + int *local_swap_type_list; // Type list index of atoms stored on this proc double *local_swap_probability; - class RanPark *random_equal; class RanPark *random_unequal; From dbc930c756696a87463f976d46f34fe19c7a1df6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2024 19:44:28 -0500 Subject: [PATCH 11/22] correct permissions --- src/MC/fix_neighbor_swap.cpp | 0 src/MC/fix_neighbor_swap.h | 0 2 files changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 src/MC/fix_neighbor_swap.cpp mode change 100755 => 100644 src/MC/fix_neighbor_swap.h diff --git a/src/MC/fix_neighbor_swap.cpp b/src/MC/fix_neighbor_swap.cpp old mode 100755 new mode 100644 diff --git a/src/MC/fix_neighbor_swap.h b/src/MC/fix_neighbor_swap.h old mode 100755 new mode 100644 From e6986cbc06208ae22457ad70178da60373c7942a Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Tue, 28 Jan 2025 18:51:07 -0700 Subject: [PATCH 12/22] Removed unused local RNG and restructured reading of command options --- src/MC/fix_neighbor_swap.cpp | 59 ++++++++++++++---------------------- src/MC/fix_neighbor_swap.h | 4 +-- 2 files changed, 23 insertions(+), 40 deletions(-) diff --git a/src/MC/fix_neighbor_swap.cpp b/src/MC/fix_neighbor_swap.cpp index 95e74bee22..d44181a981 100644 --- a/src/MC/fix_neighbor_swap.cpp +++ b/src/MC/fix_neighbor_swap.cpp @@ -66,7 +66,7 @@ static const char cite_fix_neighbor_swap_c[] = 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), random_unequal(nullptr), c_pe(nullptr) + local_swap_iatom_list(nullptr), random_equal(nullptr), c_pe(nullptr) { if (narg < 10) utils::missing_cmd_args(FLERR, "fix neighbor/swap", error); @@ -85,6 +85,17 @@ FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) : 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"); @@ -98,18 +109,12 @@ FixNeighborSwap::FixNeighborSwap(LAMMPS *lmp, int narg, char **arg) : // read options from end of input line - options(narg - 7, &arg[7]); - - if (voro_flag != 1) error->all(FLERR, "Voronoi compute required for fix neighbor/swap command"); + options(narg - 8, &arg[8]); // random number generator, same for all procs random_equal = new RanPark(lmp, seed); - // random number generator, not the same for all procs - - random_unequal = new RanPark(lmp, seed); - // set up reneighboring force_reneighbor = 1; @@ -149,7 +154,6 @@ FixNeighborSwap::~FixNeighborSwap() memory->destroy(local_swap_type_list); delete[] idregion; delete random_equal; - delete random_unequal; } /* ---------------------------------------------------------------------- @@ -163,7 +167,6 @@ void FixNeighborSwap::options(int narg, char **arg) ke_flag = 1; diff_flag = 0; rates_flag = 0; - voro_flag = 0; nswaptypes = 0; int iarg = 0; @@ -190,20 +193,6 @@ void FixNeighborSwap::options(int narg, char **arg) nswaptypes++; iarg++; } - } else if (strcmp(arg[iarg], "voro") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Illegal fix neighbor/swap command"); - - int icompute = modify->find_compute(utils::strdup(arg[iarg + 1])); - - 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"); - - voro_flag = 1; - iarg += 2; } 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"); @@ -630,7 +619,7 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) if (diff_flag) { // Calculate distance from i to each j, adjust probability of selection - // Get distance if own centr atom + // Get distance if own center atom double r = INFINITY; if (i_center >= 0) { double r = get_distance(x[temp_j], x[i_center]); } @@ -643,9 +632,9 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) if (rates_flag) { local_swap_probability[njswap_local] = - rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0)); } else { - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); + 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]; @@ -668,9 +657,9 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) if (rates_flag) { local_swap_probability[njswap_local] = - rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0)); } else { - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / r_0)); } local_probability += local_swap_probability[njswap_local]; @@ -698,9 +687,9 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) if (rates_flag) { local_swap_probability[njswap_local] = - rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0)); } else { - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / r_0)); } local_probability += local_swap_probability[njswap_local]; @@ -724,9 +713,9 @@ void FixNeighborSwap::build_i_neighbor_list(int i_center) if (rates_flag) { local_swap_probability[njswap_local] = - rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / 3.0)); + rate_list[type[temp_j] - 1] * exp(-MathSpecial::square(r / r_0)); } else { - local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / 3.0)); + local_swap_probability[njswap_local] = exp(-MathSpecial::square(r / r_0)); } local_probability += local_swap_probability[njswap_local]; @@ -877,7 +866,6 @@ void FixNeighborSwap::write_restart(FILE *fp) int n = 0; double list[6]; list[n++] = random_equal->state(); - list[n++] = random_unequal->state(); list[n++] = ubuf(next_reneighbor).d; list[n++] = nswap_attempts; list[n++] = nswap_successes; @@ -902,9 +890,6 @@ void FixNeighborSwap::restart(char *buf) seed = static_cast(list[n++]); random_equal->reset(seed); - seed = static_cast(list[n++]); - random_unequal->reset(seed); - next_reneighbor = (bigint) ubuf(list[n++]).i; nswap_attempts = static_cast(list[n++]); diff --git a/src/MC/fix_neighbor_swap.h b/src/MC/fix_neighbor_swap.h index 826ab96a65..821eda1bdc 100644 --- a/src/MC/fix_neighbor_swap.h +++ b/src/MC/fix_neighbor_swap.h @@ -43,7 +43,6 @@ class FixNeighborSwap : public Fix { 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 voro_flag; // yes = use given voronoi calculation, no = use internal voronoi calculation int ncycles; int niswap, njswap; // # of i,j swap atoms on all procs int niswap_local, njswap_local; // # of swap atoms on this proc @@ -67,7 +66,7 @@ class FixNeighborSwap : public Fix { bool unequal_cutoffs; int atom_swap_nmax; - double beta; + 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 @@ -81,7 +80,6 @@ class FixNeighborSwap : public Fix { double *local_swap_probability; class RanPark *random_equal; - class RanPark *random_unequal; class Compute *c_voro; class Compute *c_pe; From a678a3b4742a5fe132babb2ddeea6a14ceabe03e Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Tue, 28 Jan 2025 18:51:52 -0700 Subject: [PATCH 13/22] Initial update of doc file --- doc/src/fix_neighbor_swap.rst | 63 +++++++++++++++++++++++++---------- 1 file changed, 46 insertions(+), 17 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 3c97400d3f..913f479ebe 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -8,7 +8,7 @@ Syntax .. code-block:: LAMMPS - fix ID group-ID neighbor/swap N X seed T R keyword values ... + fix ID group-ID neighbor/swap N X seed T R0 voro keyword values ... * ID, group-ID are documented in :doc:`fix ` command * neighbor/swap = style name of this fix command @@ -16,9 +16,11 @@ Syntax * X = number of swaps to attempt every N steps * seed = random # seed (positive integer) * T = scaling temperature of the MC swaps (temperature units) -* R = scaling swap probability of the MC swaps (distance units) +* R0 = scaling swap probability of the MC swaps (distance units) +* voro = valid voronoi compute id (compute voronoi/atom) * one or more keyword/value pairs may be appended to args -* keyword = *types* or *mu* or *ke* or *semi-grand* or *region* +* keywords *types* and *diff* are mutually exclusive, but one must be specified +* keyword = *types* or *ke* or *region* or *diff* or *rates* .. parsed-literal:: @@ -29,32 +31,35 @@ Syntax *region* value = region-ID region-ID = ID of region to use as an exchange/move volume *diff* values = one atom type - *voro* values = valid voronoi compute id (compute voronoi/atom) - *rates* values = Ntype values to conduct variable diffusion for different atom types (unitless) + *rates* values = V1 V2 . . . Vntypes values to conduct variable diffusion for different atom types (unitless) Examples """""""" .. code-block:: LAMMPS - fix mc all neighbor/swap 10 160 15238 1000.0 diff 2 voro voroN - fix myFix all neighbor/swap 100 1 12345 298.0 region my_swap_region types 5 6 voro voroN - fix kmc all neighbor/swap 1 100 345 1.0 diff 3 rates 3 1 6 voro voroN + fix mc all neighbor/swap 10 160 15238 1000.0 3.0 diff 2 voro voroN + fix myFix all neighbor/swap 100 1 12345 298.0 3.0 region my_swap_region types 5 6 voro voroN + fix kmc all neighbor/swap 1 100 345 1.0 3.0 diff 3 rates 3 1 6 voro voroN Description """"""""""" .. versionadded:: TBD -Computes MC evaluations to enable kinetic Monte Carlo (kMC) behavior -during MD simulation through only allowing neighboring atom swaps. +Computes MC evaluations to enable kinetic Monte Carlo (kMC)-type behavior +during MD simulation through only allowing neighboring atom swaps. This +creates a hybrid type simulation of MDkMC simulation where atoms are only +swapped with their neighbors, but the swapping acceptance is perfomed by +evaluating the change in system energy using the Metropolis Criterion. Neighboring atoms are selected using a Voronoi tesselation approach. This -implementation is as described in :ref:`(Tavenner) `. +implementation is as described in :ref:`(Tavenner 2023) <_TavennerMDkMC>` +as originally intended for simulating accelerated diffusion in an MD context. The fix is called every *N* timesteps and attempts *X* swaps. The system is initialized with a random seed, using a temperature *T* for evaluating the MC energy swaps. The distance-based probability is -weighted according to *R* which sets the radius :math:`r_0` for the +weighted according to *R0* which sets the radius :math:`r_0` for the weighting .. math:: @@ -64,6 +69,10 @@ weighting where :math:`p_{ij}` is the probability of selecting atoms :math:`i` and :math:`j` for an evaluated swap. +Typically, a value around the average nearest-neighbor spacing is appropriate +for *R0*. Since this is simply a proability weighting, behavior is not +particularly sensitive to the exact value of *R0*. + The keyword *types* is submitted with two or more atom types as the value. Atoms of the first atom type are swapped with valid neighbors of all the remaining atom types. @@ -96,6 +105,9 @@ potential atoms to be swapped at the initial step, i.e. for using *fix neighbor/swap* with *diff 2*. +If atoms in the specified group are not in the voro calculated group +they will not be considered for swapping. + The keyword *rates* can modify the swap rate for each swapped type by values where the adjusted rates values are given in order of increasing atom type. The number of rates provided must equal the number of atom @@ -113,6 +125,24 @@ considered in the area given by *region-ID*. If only atoms of certain groups are expected to be in this region, the corresponding compute voronoi command can be adjusted accordingly. +Either the *types* or *diff* keyword must be specified to select atom +types for swapping + +Keywords +"""""""""""""""""""""""""""""""""""""""""""""""""""""""""" +---------------------------------------------------------- + +types = Select random atom matching first type as type I, remaining +atom types are valid for selecting atom J. +diff = Select random atom of this type as atom I, all atoms are valid +for type J. +ke = re-scale velocities when atoms are swapped based on difference in +mass +region = select only atoms I and J from region +rates = pre-factor modification to the J atom selection probability +based on atom type. + + Restart, fix_modify, output, run start/stop, minimize info """""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -149,8 +179,7 @@ doc page for more info. Also this fix requires that the :ref:`VORONOI package ` is installed, otherwise the fix will not be compiled. -A valid voronoi command which returns neighboring atoms must be used -and referenced with the *voro* keyword. +The voronoi command specified by *voro* must return neighboring atoms. When this fix is used with a :doc:`hybrid pair style ` system, only swaps between atom types of the same sub-style (or @@ -171,11 +200,11 @@ Related commands Default """"""" -The option defaults are *ke* = yes, *diff* = no, *rates* = 1 for all +The option defaults are *ke* = yes, *rates* = 1 for all atom types. ---------- -.. _Tavenner: +.. _TavennerMDkMC: -**(Tavenner)** J Tavenner, M Mendelev, J Lawson, Computational Materials Science, 218, 111929 (2023). +**(Tavenner 2023)** J Tavenner, M Mendelev, J Lawson, Computational Materials Science, 218, 111929 (2023). From 94885186b8838d7c0fc091474f7198ac17482f26 Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Thu, 27 Feb 2025 16:12:06 -0700 Subject: [PATCH 14/22] Updated doc with description of kMC algorithm --- doc/src/fix_neighbor_swap.rst | 62 ++++++++++++++++++++++++++++------- 1 file changed, 51 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 913f479ebe..9c87547321 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -20,7 +20,7 @@ Syntax * voro = 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 *ke* or *region* or *diff* or *rates* +* keyword = *types* or *diff* or *rates* or *ke* or *region* .. parsed-literal:: @@ -47,14 +47,55 @@ Description .. versionadded:: TBD -Computes MC evaluations to enable kinetic Monte Carlo (kMC)-type behavior -during MD simulation through only allowing neighboring atom swaps. This -creates a hybrid type simulation of MDkMC simulation where atoms are only -swapped with their neighbors, but the swapping acceptance is perfomed by -evaluating the change in system energy using the Metropolis Criterion. -Neighboring atoms are selected using a Voronoi tesselation approach. This -implementation is as described in :ref:`(Tavenner 2023) <_TavennerMDkMC>` -as originally intended for simulating accelerated diffusion in an MD context. +This fix computes Monte-Carlo (MC) evaluations to enable kinetic +Monte Carlo (kMC)-type behavior during MD simulation through only allowing +neighboring atom swaps. This creates a hybrid type simulation of MDkMC simulation +where atoms are only swapped with their neighbors, but the swapping acceptance is +perfomed by evaluating the change in system energy using the Metropolis Criterion. +Neighboring atoms are selected using a Voronoi tesselation approach. A detailed +explination of the original implementation of this procedure can be found in +:ref:`(Tavenner 2023) <_TavennerMDkMC>` as originally intended for simulating +accelerated diffusion in an MD context. + +Simulating inherently kineticly driven behaviors which rely on rare events +(such as atomic diffusion) is challenging for traditional Molecular Dynamics +approaches since simulations are restricted in their time-scale of events. +Since thermal vibration motion occurs on a timescale much shorter than the movement +of vacancies, such behaviors are challenging to model simultaneously. To address +this challenge, an approach from kMC simulations is adpoted where rare events can +be sampled at selected rates. By selecting such swap behaviors, the process +of atomic diffusion can be approximated during an MD simulation, effectively +decoupling the MD atomic vibrational time and the timescale of atomic hopping. + +To achieve such simulations, this algorithm takes the following approach. First, +the MD simulation is stopped after a given number of steps to perform atom swaps. +Given this instantaneous configuration from the MD simulation, Voronoi neighbors +are computed for all valid swap atoms. From the list of valid swap atoms, one atom +I is selected at random across the entire simulation. One if its Voronoi neighbors +that is a valid atom to swap is then selected. The atom ID is communicated to all +processors, such that if the neighbors are on different processors the swap still +occurs. The two atom types are swapped, and the change in system energy from before +the swap is compared using the Metropolis Criterion. This evaluation of the energy +change is a global calculation, such that it has a computational cost similar to +that of an MD timestep. If the swap is accepted from the Metropolis Criterion, the +atoms remain swapped. Else, the atoms are returned to their original types. This +process of MC evaluation is repeated for a given number of iterations until the +original MD simulation is resumed from the new state, where any successfully +swapped atoms have changed type, though the global system balance is preserved. + +A few key notes regarding this implementation are as follows. The parallel +efficiency of this algorithm is similar to that of other MC approaches. I.e, +due to the additional energy calculations for the MC steps, efficiency is +improved with a smaller number of atoms per processor than standalone MD simulation +since there is more weighting on the calculation of a given atomic domain and +minor additonal communication load. Communication of the atom ids to be swapped +between processors is negligible. Efficiency will additionally be much worse for +pair styles with different per-atom cutoffs, since the neighbor list will need to +be rebuilt between swap events. Limitations are imposed on the Voronoi neighbors +to restrict swapping of atoms which are outside of a reasonable cutoff. + +Input Parameters Usage +""""""""""" The fix is called every *N* timesteps and attempts *X* swaps. The system is initialized with a random seed, using a temperature *T* for @@ -128,9 +169,8 @@ voronoi command can be adjusted accordingly. Either the *types* or *diff* keyword must be specified to select atom types for swapping -Keywords +Keyword Summary """""""""""""""""""""""""""""""""""""""""""""""""""""""""" ----------------------------------------------------------- types = Select random atom matching first type as type I, remaining atom types are valid for selecting atom J. From 6dacf5d52c3fe107ef96e4beba9be682f288c7dc Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Sat, 15 Mar 2025 11:36:20 -0600 Subject: [PATCH 15/22] Update fix_neighbor_swap.rst I fixed some typos and shortened the text a bit. --- doc/src/fix_neighbor_swap.rst | 75 ++++++++++++++--------------------- 1 file changed, 30 insertions(+), 45 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 9c87547321..226926b0fb 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -24,7 +24,7 @@ Syntax .. parsed-literal:: - *types* values = two or more atom types (1-Ntypes or type label) + *types* values = two or more atom types (Integers in range [1,Ntypes] or type labels) *ke* value = *no* or *yes* *no* = no conservation of kinetic energy after atom swaps *yes* = kinetic energy is conserved after atom swaps @@ -47,52 +47,50 @@ Description .. versionadded:: TBD -This fix computes Monte-Carlo (MC) evaluations to enable kinetic -Monte Carlo (kMC)-type behavior during MD simulation through only allowing -neighboring atom swaps. This creates a hybrid type simulation of MDkMC simulation -where atoms are only swapped with their neighbors, but the swapping acceptance is -perfomed by evaluating the change in system energy using the Metropolis Criterion. +This fix computes Monte Carlo (MC) evaluations to enable kinetic +Monte Carlo (kMC) behavior during an MD simulation by allowing +neighboring atoms to swap their positions. This creates a hybrid MD/kMC +simulation +where atoms are swapped only with their neighbors, and the swapping acceptance is +perfomed by evaluating the change in system energy using the Metropolis criterion. Neighboring atoms are selected using a Voronoi tesselation approach. A detailed explination of the original implementation of this procedure can be found in :ref:`(Tavenner 2023) <_TavennerMDkMC>` as originally intended for simulating accelerated diffusion in an MD context. -Simulating inherently kineticly driven behaviors which rely on rare events +Simulating inherently kinetically-limited behaviors which rely on rare events (such as atomic diffusion) is challenging for traditional Molecular Dynamics -approaches since simulations are restricted in their time-scale of events. -Since thermal vibration motion occurs on a timescale much shorter than the movement -of vacancies, such behaviors are challenging to model simultaneously. To address -this challenge, an approach from kMC simulations is adpoted where rare events can +approaches, since thermal vibration motion occurs on a timescale much +shorter than local structural changes such as vacancy hopping. +To address +this challenge, an approach from kMC simulations is adopted where rare events can be sampled at selected rates. By selecting such swap behaviors, the process of atomic diffusion can be approximated during an MD simulation, effectively decoupling the MD atomic vibrational time and the timescale of atomic hopping. To achieve such simulations, this algorithm takes the following approach. First, -the MD simulation is stopped after a given number of steps to perform atom swaps. -Given this instantaneous configuration from the MD simulation, Voronoi neighbors -are computed for all valid swap atoms. From the list of valid swap atoms, one atom -I is selected at random across the entire simulation. One if its Voronoi neighbors -that is a valid atom to swap is then selected. The atom ID is communicated to all -processors, such that if the neighbors are on different processors the swap still -occurs. The two atom types are swapped, and the change in system energy from before -the swap is compared using the Metropolis Criterion. This evaluation of the energy +the MD simulation is stopped *N* steps to attempt *X* atom swaps. +Voronoi neighbors are computed for all valid swap atoms. +For each swap attempt, one atom +I is selected at random from the global list of valid atoms. One if its Voronoi neighbors +that is a valid atom to swap is then selected. The neighbor atom ID is communicated to all +processors, as it may be owned by a different processor. +The two atom types are then swapped, and the change in system energy from before +the swap is used to evaluate the Metropolis acceptance criterion. This evaluation of the energy change is a global calculation, such that it has a computational cost similar to -that of an MD timestep. If the swap is accepted from the Metropolis Criterion, the -atoms remain swapped. Else, the atoms are returned to their original types. This +that of an MD timestep. If the swap is accepted from the Metropolis criterion, the +atoms remain swapped. If the swap is rejected, the atoms are reverted to their original types. This process of MC evaluation is repeated for a given number of iterations until the original MD simulation is resumed from the new state, where any successfully swapped atoms have changed type, though the global system balance is preserved. A few key notes regarding this implementation are as follows. The parallel -efficiency of this algorithm is similar to that of other MC approaches. I.e, -due to the additional energy calculations for the MC steps, efficiency is -improved with a smaller number of atoms per processor than standalone MD simulation -since there is more weighting on the calculation of a given atomic domain and -minor additonal communication load. Communication of the atom ids to be swapped -between processors is negligible. Efficiency will additionally be much worse for -pair styles with different per-atom cutoffs, since the neighbor list will need to -be rebuilt between swap events. Limitations are imposed on the Voronoi neighbors -to restrict swapping of atoms which are outside of a reasonable cutoff. +efficiency of this algorithm is similar to that of other MC approaches, i.e, +the global potential energy must be calculated after each attempted swap. +Efficiency is sensitive to the maximum cutoff distance for the pair style, +since the neighbor list will need to be rebuilt between swap events. +Limitations are imposed on the Voronoi neighbors +to restrict swapping of atoms that are outside of a reasonable cutoff. Input Parameters Usage """"""""""" @@ -111,7 +109,7 @@ where :math:`p_{ij}` is the probability of selecting atoms :math:`i` and :math:`j` for an evaluated swap. Typically, a value around the average nearest-neighbor spacing is appropriate -for *R0*. Since this is simply a proability weighting, behavior is not +for *R0*. Since this is simply a probability weighting, behavior is not particularly sensitive to the exact value of *R0*. The keyword *types* is submitted with two or more atom types as the @@ -121,7 +119,7 @@ all the remaining atom types. The keyword *diff* is used for implementation of simulated diffusion of a given atom type as given by *diff type*. This command selects all atom types as acceptable swap types to a centrally selected atom of type -*type*. This includes the atom type specified by the diff keyword to +*type*. This includes the atom type specified by the *diff* keyword to account for self-diffusion hops of an atom type with itself. Keyword *voro* is currently required, and is implemented as @@ -169,19 +167,6 @@ voronoi command can be adjusted accordingly. Either the *types* or *diff* keyword must be specified to select atom types for swapping -Keyword Summary -"""""""""""""""""""""""""""""""""""""""""""""""""""""""""" - -types = Select random atom matching first type as type I, remaining -atom types are valid for selecting atom J. -diff = Select random atom of this type as atom I, all atoms are valid -for type J. -ke = re-scale velocities when atoms are swapped based on difference in -mass -region = select only atoms I and J from region -rates = pre-factor modification to the J atom selection probability -based on atom type. - Restart, fix_modify, output, run start/stop, minimize info """""""""""""""""""""""""""""""""""""""""""""""""""""""""" From 3a6ea0808daa42d5f681ec76cd53e4a8760989c8 Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Mon, 14 Apr 2025 15:36:33 -0600 Subject: [PATCH 16/22] Merged changes to documentation for further user clarity --- doc/src/fix_neighbor_swap.rst | 268 +++++++++++++++++----------------- 1 file changed, 136 insertions(+), 132 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 9c87547321..da1e89bed5 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -18,19 +18,19 @@ Syntax * T = scaling temperature of the MC swaps (temperature units) * R0 = scaling swap probability of the MC swaps (distance units) * voro = valid voronoi compute id (compute voronoi/atom) -* one or more keyword/value pairs may be appended to args +* two 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 *rates* or *ke* or *region* +* keyword = *types* or *diff* or *ke* or *region* or *rates* .. parsed-literal:: - *types* values = two or more atom types (1-Ntypes or type label) - *ke* value = *no* or *yes* - *no* = no conservation of kinetic energy after atom swaps + *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 - *diff* values = one atom type *rates* values = V1 V2 . . . Vntypes values to conduct variable diffusion for different atom types (unitless) Examples @@ -38,6 +38,7 @@ Examples .. code-block:: LAMMPS + compute voroN all voronoi/atom neighbors yes fix mc all neighbor/swap 10 160 15238 1000.0 3.0 diff 2 voro voroN fix myFix all neighbor/swap 100 1 12345 298.0 3.0 region my_swap_region types 5 6 voro voroN fix kmc all neighbor/swap 1 100 345 1.0 3.0 diff 3 rates 3 1 6 voro voroN @@ -47,140 +48,144 @@ Description .. versionadded:: TBD -This fix computes Monte-Carlo (MC) evaluations to enable kinetic -Monte Carlo (kMC)-type behavior during MD simulation through only allowing -neighboring atom swaps. This creates a hybrid type simulation of MDkMC simulation -where atoms are only swapped with their neighbors, but the swapping acceptance is -perfomed by evaluating the change in system energy using the Metropolis Criterion. -Neighboring atoms are selected using a Voronoi tesselation approach. A detailed -explination of the original implementation of this procedure can be found in -:ref:`(Tavenner 2023) <_TavennerMDkMC>` as originally intended for simulating -accelerated diffusion in an MD context. +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 constrast to the :doc:`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. -Simulating inherently kineticly driven behaviors which rely on rare events -(such as atomic diffusion) is challenging for traditional Molecular Dynamics -approaches since simulations are restricted in their time-scale of events. -Since thermal vibration motion occurs on a timescale much shorter than the movement -of vacancies, such behaviors are challenging to model simultaneously. To address -this challenge, an approach from kMC simulations is adpoted where rare events can -be sampled at selected rates. By selecting such swap behaviors, the process -of atomic diffusion can be approximated during an MD simulation, effectively -decoupling the MD atomic vibrational time and the timescale of atomic hopping. +Neighboring atoms are defined by using a Voronoi tesselation performed +by the :doc:`compute voronoi/atom ` command. +Two atoms are neighbors if their Voronoi cells share a common face +(3d) or edge (2d). -To achieve such simulations, this algorithm takes the following approach. First, -the MD simulation is stopped after a given number of steps to perform atom swaps. -Given this instantaneous configuration from the MD simulation, Voronoi neighbors -are computed for all valid swap atoms. From the list of valid swap atoms, one atom -I is selected at random across the entire simulation. One if its Voronoi neighbors -that is a valid atom to swap is then selected. The atom ID is communicated to all -processors, such that if the neighbors are on different processors the swap still -occurs. The two atom types are swapped, and the change in system energy from before -the swap is compared using the Metropolis Criterion. This evaluation of the energy -change is a global calculation, such that it has a computational cost similar to -that of an MD timestep. If the swap is accepted from the Metropolis Criterion, the -atoms remain swapped. Else, the atoms are returned to their original types. This -process of MC evaluation is repeated for a given number of iterations until the -original MD simulation is resumed from the new state, where any successfully -swapped atoms have changed type, though the global system balance is preserved. +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 few key notes regarding this implementation are as follows. The parallel -efficiency of this algorithm is similar to that of other MC approaches. I.e, -due to the additional energy calculations for the MC steps, efficiency is -improved with a smaller number of atoms per processor than standalone MD simulation -since there is more weighting on the calculation of a given atomic domain and -minor additonal communication load. Communication of the atom ids to be swapped -between processors is negligible. Efficiency will additionally be much worse for -pair styles with different per-atom cutoffs, since the neighbor list will need to -be rebuilt between swap events. Limitations are imposed on the Voronoi neighbors -to restrict swapping of atoms which are outside of a reasonable cutoff. +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. -Input Parameters Usage -""""""""""" +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 fix is called every *N* timesteps and attempts *X* swaps. The system -is initialized with a random seed, using a temperature *T* for -evaluating the MC energy swaps. The distance-based probability is -weighted according to *R0* which sets the radius :math:`r_0` for the -weighting +The algorithm implemented by this fix is as follows. The MD +simulation is paused every *N* stepsA 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) Simliar to other MC algorithms in LAMMPS, an optimized 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 probabiltiy 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 atoms :math:`i` and -:math:`j` for an evaluated swap. +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 proability weighting, the swapping behavior is not very sensitive to +the exact value of *R0*. -Typically, a value around the average nearest-neighbor spacing is appropriate -for *R0*. Since this is simply a proability weighting, behavior is not -particularly sensitive to the exact value of *R0*. +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 *types* is submitted with two or more atom types as the -value. Atoms of the first atom type are swapped with valid neighbors of -all the remaining atom types. +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. -The keyword *diff* is used for implementation of simulated diffusion of -a given atom type as given by *diff type*. This command selects all atom -types as acceptable swap types to a centrally selected atom of type -*type*. This includes the atom type specified by the diff keyword to -account for self-diffusion hops of an atom type with itself. - -Keyword *voro* is currently required, and is implemented as - -.. code-block:: LAMMPS - - voro compute-ID - -where *compute-ID* is the ID of a corresponding Voronoi computation with -neighbor list, i.e. +The keyword *voro* is required. Its *vID* value is the compute-ID of +a :doc:`compute voronoi/atom ` command like +this: .. code-block:: LAMMPS compute compute-ID group-ID voronoi/atom neighbors yes -The group selected for computing *voro* should correspond to all the -potential atoms to be swapped at the initial step, i.e. +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 +potentialy 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. -.. code-block:: LAMMPS +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. - group group-ID type 2 +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 ` command +for details. This means the group-ID for the :doc:`compute +voronoi/atom ` command also need only contain +atoms within the region. -for using *fix neighbor/swap* with *diff 2*. - -If atoms in the specified group are not in the voro calculated group -they will not be considered for swapping. - -The keyword *rates* can modify the swap rate for each swapped type by -values where the adjusted rates values are given in order of increasing -atom type. The number of rates provided must equal the number of atom -types in the simulation. In the third provided example above, a -simulation is conducted with three atom types where the third atom type -is the one sampled for attempted swaps. All three atom types are -considered valid swaps, but atoms of type 1 will be selected three times -as often as atoms of type 2. Conversely, atoms of type 3 are six times -more likely to be selected than atoms of type two and twice as likely as -atoms of type 1. - -Finally, the *region* keyword is implemented as in other atomic fixes, -where the *region region-ID* command indicates that atom swaps only be -considered in the area given by *region-ID*. If only atoms of certain -groups are expected to be in this region, the corresponding compute -voronoi command can be adjusted accordingly. - -Either the *types* or *diff* keyword must be specified to select atom -types for swapping - -Keyword Summary -"""""""""""""""""""""""""""""""""""""""""""""""""""""""""" - -types = Select random atom matching first type as type I, remaining -atom types are valid for selecting atom J. -diff = Select random atom of this type as atom I, all atoms are valid -for type J. -ke = re-scale velocities when atoms are swapped based on difference in -mass -region = select only atoms I and J from region -rates = pre-factor modification to the J atom selection probability -based on atom type. +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, Restart, fix_modify, output, run start/stop, minimize info @@ -215,19 +220,18 @@ 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 ` -doc page for more info. Also this fix requires that the -:ref:`VORONOI package ` is installed, otherwise the fix -will not be compiled. +doc page for more info. Also this fix requires that the :ref:`VORONOI +package ` be installed, otherwise the fix will not be +compiled. -The voronoi command specified by *voro* must return neighboring atoms. - -When this fix is used with a :doc:`hybrid pair style ` -system, only swaps between atom types of the same sub-style (or -combination of sub-styles) are permitted. +The :doc:`compute voronoi/atom ` command +referenced by keyword *voro* 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 flag must be set to off since the -implemented algorithm will not be able to re-scale velocity properly. +(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 """""""""""""""" @@ -240,11 +244,11 @@ Related commands Default """"""" -The option defaults are *ke* = yes, *rates* = 1 for all -atom types. +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). +**(Tavenner 2023)** J Tavenner, M Mendelev, J Lawson, Computational + Materials Science, 218, 111929 (2023). From 51de62ce057f144ed59f5b4737734315e92bf3b4 Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Mon, 14 Apr 2025 15:58:19 -0600 Subject: [PATCH 17/22] Reconcile version changes to documentation and polish --- doc/src/fix_neighbor_swap.rst | 102 ---------------------------------- 1 file changed, 102 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 5f0373d13b..da1e89bed5 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -25,13 +25,8 @@ Syntax .. parsed-literal:: *types* values = two or more atom types (Integers in range [1,Ntypes] or type labels) -<<<<<<< HEAD *diff* values = one atom type *ke* value = *yes* or *no* -======= - *ke* value = *no* or *yes* - *no* = no conservation of kinetic energy after atom swaps ->>>>>>> a3bc1a6c0b09ee7a84167ab8aa626e4a93443846 *yes* = kinetic energy is conserved after atom swaps *no* = no conservation of kinetic energy after atom swaps *region* value = region-ID @@ -53,7 +48,6 @@ Description .. versionadded:: TBD -<<<<<<< HEAD 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 constrast to the :doc:`fix @@ -76,52 +70,6 @@ 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. -======= -This fix computes Monte Carlo (MC) evaluations to enable kinetic -Monte Carlo (kMC) behavior during an MD simulation by allowing -neighboring atoms to swap their positions. This creates a hybrid MD/kMC -simulation -where atoms are swapped only with their neighbors, and the swapping acceptance is -perfomed by evaluating the change in system energy using the Metropolis criterion. -Neighboring atoms are selected using a Voronoi tesselation approach. A detailed -explination of the original implementation of this procedure can be found in -:ref:`(Tavenner 2023) <_TavennerMDkMC>` as originally intended for simulating -accelerated diffusion in an MD context. - -Simulating inherently kinetically-limited behaviors which rely on rare events -(such as atomic diffusion) is challenging for traditional Molecular Dynamics -approaches, since thermal vibration motion occurs on a timescale much -shorter than local structural changes such as vacancy hopping. -To address -this challenge, an approach from kMC simulations is adopted where rare events can -be sampled at selected rates. By selecting such swap behaviors, the process -of atomic diffusion can be approximated during an MD simulation, effectively -decoupling the MD atomic vibrational time and the timescale of atomic hopping. - -To achieve such simulations, this algorithm takes the following approach. First, -the MD simulation is stopped *N* steps to attempt *X* atom swaps. -Voronoi neighbors are computed for all valid swap atoms. -For each swap attempt, one atom -I is selected at random from the global list of valid atoms. One if its Voronoi neighbors -that is a valid atom to swap is then selected. The neighbor atom ID is communicated to all -processors, as it may be owned by a different processor. -The two atom types are then swapped, and the change in system energy from before -the swap is used to evaluate the Metropolis acceptance criterion. This evaluation of the energy -change is a global calculation, such that it has a computational cost similar to -that of an MD timestep. If the swap is accepted from the Metropolis criterion, the -atoms remain swapped. If the swap is rejected, the atoms are reverted to their original types. This -process of MC evaluation is repeated for a given number of iterations until the -original MD simulation is resumed from the new state, where any successfully -swapped atoms have changed type, though the global system balance is preserved. - -A few key notes regarding this implementation are as follows. The parallel -efficiency of this algorithm is similar to that of other MC approaches, i.e, -the global potential energy must be calculated after each attempted swap. -Efficiency is sensitive to the maximum cutoff distance for the pair style, -since the neighbor list will need to be rebuilt between swap events. -Limitations are imposed on the Voronoi neighbors -to restrict swapping of atoms that are outside of a reasonable cutoff. ->>>>>>> a3bc1a6c0b09ee7a84167ab8aa626e4a93443846 Simulating inherently kinetically-limited behaviors which rely on rare events (such as atomic diffusion in a solid) is challenging for @@ -189,15 +137,9 @@ average nearest-neighbor spacing is appropriate. Since this is simply a proability weighting, the swapping behavior is not very sensitive to the exact value of *R0*. -<<<<<<< HEAD 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. -======= -Typically, a value around the average nearest-neighbor spacing is appropriate -for *R0*. Since this is simply a probability weighting, behavior is not -particularly sensitive to the exact value of *R0*. ->>>>>>> a3bc1a6c0b09ee7a84167ab8aa626e4a93443846 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 @@ -205,26 +147,9 @@ 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. -<<<<<<< HEAD The keyword *voro* is required. Its *vID* value is the compute-ID of a :doc:`compute voronoi/atom ` command like this: -======= -The keyword *diff* is used for implementation of simulated diffusion of -a given atom type as given by *diff type*. This command selects all atom -types as acceptable swap types to a centrally selected atom of type -*type*. This includes the atom type specified by the *diff* keyword to -account for self-diffusion hops of an atom type with itself. - -Keyword *voro* is currently required, and is implemented as - -.. code-block:: LAMMPS - - voro compute-ID - -where *compute-ID* is the ID of a corresponding Voronoi computation with -neighbor list, i.e. ->>>>>>> a3bc1a6c0b09ee7a84167ab8aa626e4a93443846 .. code-block:: LAMMPS @@ -251,7 +176,6 @@ for details. This means the group-ID for the :doc:`compute voronoi/atom ` command also need only contain atoms within the region. -<<<<<<< HEAD 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 @@ -262,32 +186,6 @@ 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, -======= -for using *fix neighbor/swap* with *diff 2*. - -If atoms in the specified group are not in the voro calculated group -they will not be considered for swapping. - -The keyword *rates* can modify the swap rate for each swapped type by -values where the adjusted rates values are given in order of increasing -atom type. The number of rates provided must equal the number of atom -types in the simulation. In the third provided example above, a -simulation is conducted with three atom types where the third atom type -is the one sampled for attempted swaps. All three atom types are -considered valid swaps, but atoms of type 1 will be selected three times -as often as atoms of type 2. Conversely, atoms of type 3 are six times -more likely to be selected than atoms of type two and twice as likely as -atoms of type 1. - -Finally, the *region* keyword is implemented as in other atomic fixes, -where the *region region-ID* command indicates that atom swaps only be -considered in the area given by *region-ID*. If only atoms of certain -groups are expected to be in this region, the corresponding compute -voronoi command can be adjusted accordingly. - -Either the *types* or *diff* keyword must be specified to select atom -types for swapping ->>>>>>> a3bc1a6c0b09ee7a84167ab8aa626e4a93443846 Restart, fix_modify, output, run start/stop, minimize info From 4bfdd3eb3481e15aa1fa7e9be746a952ec454b8d Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Mon, 14 Apr 2025 16:22:11 -0600 Subject: [PATCH 18/22] Remove trailing whitespace --- doc/src/fix_neighbor_swap.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index da1e89bed5..274c4ea939 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -48,11 +48,11 @@ Description .. versionadded:: TBD -This fix performs Monte-Carlo (MC) evaluations to enable kinetic -Monte Carlo (kMC)-type behavior during MD simulation by allowing +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 constrast to the :doc:`fix atom/swap ` command which swaps pairs of atoms anywhere -in the simulation domain, the restriction of the MC swapping to +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 @@ -115,7 +115,7 @@ 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 probabiltiy scaling. +use of a distance-weighted probabiltiy scaling. ---------- From 5d9a7d154d02a95ad3c4e32a83105f8529669f49 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 15 Apr 2025 22:48:25 -0400 Subject: [PATCH 19/22] address spelling and formatting issues --- doc/src/fix_neighbor_swap.rst | 71 +++++++++++---------- doc/utils/sphinx-config/false_positives.txt | 1 + 2 files changed, 39 insertions(+), 33 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 274c4ea939..fbe793eef9 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -50,7 +50,7 @@ Description 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 constrast to the :doc:`fix +neighboring atoms to swap their positions. In contrast to the :doc:`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. @@ -68,7 +68,7 @@ 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>` +algorithm can be found in :ref:`(Tavenner 2023) ` where it was used to simulated accelerated diffusion in an MD context. Simulating inherently kinetically-limited behaviors which rely on rare @@ -81,41 +81,46 @@ 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* stepsA 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*. +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. + 1. The cost of a global energy evaluation is similar to that of an MD + timestep. -(2) Simliar to other MC algorithms in LAMMPS, an optimized 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. + 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. + 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 probabiltiy scaling. +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. ---------- @@ -133,8 +138,8 @@ sets the radius :math:`r_0` in this formula 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 proability weighting, the swapping behavior is not very sensitive to +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 keyword *types* takes two or more atom types as its values. Only @@ -157,7 +162,7 @@ this: 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 -potentialy select. I.e. the group-ID used in the voronoi compute should +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. @@ -221,7 +226,7 @@ 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 ` doc page for more info. Also this fix requires that the :ref:`VORONOI -package ` be installed, otherwise the fix will not be +package ` is installed, otherwise the fix will not be compiled. The :doc:`compute voronoi/atom ` command diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index afb0fe3e95..115acc575b 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -4109,6 +4109,7 @@ volpress volumetric von Voro +voro Vorobyov voronoi Voronoi From 6436cc87b7e4412b3748731b04be0cc8a5a36150 Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Fri, 9 May 2025 14:23:03 -0600 Subject: [PATCH 20/22] Corrected inconsistent voro-ID compute references and examples --- doc/src/fix_neighbor_swap.rst | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 274c4ea939..4095f23eef 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -17,7 +17,7 @@ Syntax * 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 = valid voronoi compute id (compute voronoi/atom) +* voro-ID = valid voronoi compute id (compute voronoi/atom) * two 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* @@ -39,9 +39,9 @@ Examples .. code-block:: LAMMPS compute voroN all voronoi/atom neighbors yes - fix mc all neighbor/swap 10 160 15238 1000.0 3.0 diff 2 voro voroN - fix myFix all neighbor/swap 100 1 12345 298.0 3.0 region my_swap_region types 5 6 voro voroN - fix kmc all neighbor/swap 1 100 345 1.0 3.0 diff 3 rates 3 1 6 voro voroN + 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 """"""""""" @@ -137,6 +137,17 @@ average nearest-neighbor spacing is appropriate. Since this is simply a proability 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 ` 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 ` 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. @@ -147,14 +158,6 @@ 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. -The keyword *voro* is required. Its *vID* value is the compute-ID of -a :doc:`compute voronoi/atom ` command like -this: - -.. code-block:: LAMMPS - - compute compute-ID group-ID voronoi/atom neighbors yes - 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 potentialy select. I.e. the group-ID used in the voronoi compute should @@ -225,7 +228,7 @@ package ` be installed, otherwise the fix will not be compiled. The :doc:`compute voronoi/atom ` command -referenced by keyword *voro* must return neighboring atoms as +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 From 17da04f07b335d4e1edd2841a1e746998f1f347f Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Fri, 9 May 2025 14:25:08 -0600 Subject: [PATCH 21/22] Cleaned up language errors --- doc/src/fix_neighbor_swap.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 4095f23eef..6b26f99dff 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -189,6 +189,8 @@ 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 @@ -196,8 +198,8 @@ Restart, fix_modify, output, run start/stop, minimize info This fix writes the state of the fix to :doc:`binary restart files `. This includes information about the random number generator -seed, the next timestep for MC exchanges, the number of exchange -attempts and successes, etc. See the :doc:`read_restart ` +seed, the next timestep for MC exchanges, and the number of exchange +attempts and successes. See the :doc:`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. From ef89edc4c634f28f8f4e16e4518bbc4e3db8b928 Mon Sep 17 00:00:00 2001 From: Jacob Tavenner Date: Fri, 9 May 2025 14:53:27 -0600 Subject: [PATCH 22/22] Additional edits --- doc/src/fix_neighbor_swap.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_neighbor_swap.rst b/doc/src/fix_neighbor_swap.rst index 0f771e8ce9..dffbc93217 100644 --- a/doc/src/fix_neighbor_swap.rst +++ b/doc/src/fix_neighbor_swap.rst @@ -8,7 +8,7 @@ Syntax .. code-block:: LAMMPS - fix ID group-ID neighbor/swap N X seed T R0 voro keyword values ... + fix ID group-ID neighbor/swap N X seed T R0 voro-ID keyword values ... * ID, group-ID are documented in :doc:`fix ` command * neighbor/swap = style name of this fix command @@ -18,7 +18,7 @@ Syntax * 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) -* two or more keyword/value pairs may be appended to args +* 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*