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

Merge pull request #2992 from lammps/molswap

Browse Source 
Add a new fix mol/swap command
This commit is contained in:
Axel Kohlmeyer
2021-10-13 21:33:52 -04:00
committed by GitHub
parent c322064ff3 c5617dc006
commit a1fb6902d5
31 changed files with 13084 additions and 113 deletions
Download Patch File Download Diff File Expand all files Collapse all files

123
src/MC/fix_atom_swap.cpp
View File

@ -74,11 +74,13 @@ FixAtomSwap::FixAtomSwap(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);
beta = 1.0/(force->boltz*temperature);
if (nevery <= 0) error->all(FLERR,"Illegal fix atom/swap command");
if (ncycles < 0) error->all(FLERR,"Illegal fix atom/swap command");
if (seed <= 0) error->all(FLERR,"Illegal fix atom/swap command");
if (temperature <= 0.0) error->all(FLERR,"Illegal fix atom/swap command");
beta = 1.0/(force->boltz*temperature);
memory->create(type_list,atom->ntypes,"atom/swap:type_list");
memory->create(mu,atom->ntypes+1,"atom/swap:mu");
@ -142,7 +144,7 @@ void FixAtomSwap::options(int narg, char **arg)
if (narg < 0) error->all(FLERR,"Illegal fix atom/swap command");
regionflag = 0;
conserve_ke_flag = 1;
ke_flag = 1;
semi_grand_flag = 0;
nswaptypes = 0;
nmutypes = 0;
@ -160,7 +162,7 @@ void FixAtomSwap::options(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg],"ke") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix atom/swap command");
conserve_ke_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
ke_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"semi-grand") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix atom/swap command");
@ -301,6 +303,8 @@ void FixAtomSwap::pre_exchange()
if (next_reneighbor != update->ntimestep) return;
// insure current system is ready to compute energy
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
comm->exchange();
@ -309,8 +313,13 @@ void FixAtomSwap::pre_exchange()
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build(1);
// energy_stored = energy of current state
// will be updated after accepted swaps
energy_stored = energy_full();
// attempt Ncycle atom swaps
int nsuccess = 0;
if (semi_grand_flag) {
update_semi_grand_atoms_list();
@ -320,23 +329,30 @@ void FixAtomSwap::pre_exchange()
for (int i = 0; i < ncycles; i++) nsuccess += attempt_swap();
}
// udpate MC stats
nswap_attempts += ncycles;
nswap_successes += nsuccess;
energy_full();
next_reneighbor = update->ntimestep + nevery;
}
/* ----------------------------------------------------------------------
Note: atom charges are assumed equal and so are not updated
attempt a semd-grand swap of a single atom
compare before/after energy and accept/reject the swap
NOTE: atom charges are assumed equal and so are not updated
------------------------------------------------------------------------- */
int FixAtomSwap::attempt_semi_grand()
{
if (nswap == 0) return 0;
// pre-swap energy
double energy_before = energy_stored;
// pick a random atom and perform swap
int itype,jtype,jswaptype;
int i = pick_semi_grand_atom();
if (i >= 0) {
@ -350,9 +366,12 @@ int FixAtomSwap::attempt_semi_grand()
atom->type[i] = jtype;
}
// 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);
@ -362,6 +381,8 @@ int FixAtomSwap::attempt_semi_grand()
comm->forward_comm_fix(this);
}
// post-swap energy
if (force->kspace) force->kspace->qsum_qsq();
double energy_after = energy_full();
@ -374,10 +395,12 @@ int FixAtomSwap::attempt_semi_grand()
int success_all = 0;
MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
// swap accepted, return 1
if (success_all) {
update_semi_grand_atoms_list();
energy_stored = energy_after;
if (conserve_ke_flag) {
if (ke_flag) {
if (i >= 0) {
atom->v[i][0] *= sqrt_mass_ratio[itype][jtype];
atom->v[i][1] *= sqrt_mass_ratio[itype][jtype];
@ -385,38 +408,35 @@ int FixAtomSwap::attempt_semi_grand()
}
}
return 1;
} else {
if (i >= 0) {
atom->type[i] = itype;
}
if (force->kspace) force->kspace->qsum_qsq();
energy_stored = energy_before;
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_fix(this);
}
}
// swap not accepted, return 0
// restore the swapped atom
// do not need to re-call comm->borders() and rebuild neighbor list
// since will be done on next cycle or in Verlet when this fix finishes
if (i >= 0) atom->type[i] = itype;
if (force->kspace) force->kspace->qsum_qsq();
return 0;
}
/* ----------------------------------------------------------------------
attempt a swap of a pair of atoms
compare before/after energy and accept/reject the swap
------------------------------------------------------------------------- */
int FixAtomSwap::attempt_swap()
{
if ((niswap == 0) || (njswap == 0)) return 0;
// pre-swap energy
double energy_before = energy_stored;
// pick a random pair of atoms
// swap their properties
int i = pick_i_swap_atom();
int j = pick_j_swap_atom();
int itype = type_list[0];
@ -431,6 +451,10 @@ int FixAtomSwap::attempt_swap()
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();
@ -443,13 +467,17 @@ int FixAtomSwap::attempt_swap()
comm->forward_comm_fix(this);
}
// post-swap energy
double energy_after = energy_full();
// swap accepted, return 1
// if ke_flag, rescale atom velocities
if (random_equal->uniform() <
exp(beta*(energy_before - energy_after))) {
update_swap_atoms_list();
energy_stored = energy_after;
if (conserve_ke_flag) {
if (ke_flag) {
if (i >= 0) {
atom->v[i][0] *= sqrt_mass_ratio[itype][jtype];
atom->v[i][1] *= sqrt_mass_ratio[itype][jtype];
@ -461,30 +489,24 @@ int FixAtomSwap::attempt_swap()
atom->v[j][2] *= sqrt_mass_ratio[jtype][itype];
}
}
energy_stored = energy_after;
return 1;
} else {
if (i >= 0) {
atom->type[i] = type_list[0];
if (atom->q_flag) atom->q[i] = qtype[0];
}
if (j >= 0) {
atom->type[j] = type_list[1];
if (atom->q_flag) atom->q[j] = qtype[1];
}
energy_stored = energy_before;
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_fix(this);
}
}
// swap not accepted, return 0
// restore the swapped itype & jtype atoms
// do not need to re-call comm->borders() and rebuild neighbor list
// since will be done on next cycle or in Verlet when this fix finishes
if (i >= 0) {
atom->type[i] = type_list[0];
if (atom->q_flag) atom->q[i] = qtype[0];
}
if (j >= 0) {
atom->type[j] = type_list[1];
if (atom->q_flag) atom->q[j] = qtype[1];
}
return 0;
}
@ -497,7 +519,6 @@ double FixAtomSwap::energy_full()
int eflag = 1;
int vflag = 0;
if (modify->n_pre_neighbor) modify->pre_neighbor();
if (modify->n_pre_force) modify->pre_force(vflag);
if (force->pair) force->pair->compute(eflag,vflag);

22
src/MC/fix_atom_swap.h
View File

@ -17,8 +17,8 @@ FixStyle(atom/swap,FixAtomSwap);
// clang-format on
#else
#ifndef LMP_FIX_MCSWAP_H
#define LMP_FIX_MCSWAP_H
#ifndef LMP_FIX_ATOM_SWAP_H
#define LMP_FIX_ATOM_SWAP_H
#include "fix.h"
@ -31,14 +31,6 @@ class FixAtomSwap : public Fix {
int setmask();
void init();
void pre_exchange();
int attempt_semi_grand();
int attempt_swap();
double energy_full();
int pick_semi_grand_atom();
int pick_i_swap_atom();
int pick_j_swap_atom();
void update_semi_grand_atoms_list();
void update_swap_atoms_list();
int pack_forward_comm(int, int *, double *, int, int *);
void unpack_forward_comm(int, int, double *);
double compute_vector(int);
@ -48,7 +40,7 @@ class FixAtomSwap : public Fix {
private:
int nevery, seed;
int conserve_ke_flag; // yes = conserve ke, no = do not conserve ke
int ke_flag; // yes = conserve ke, no = do not conserve ke
int semi_grand_flag; // yes = semi-grand canonical, no = constant composition
int ncycles;
int niswap, njswap; // # of i,j swap atoms on all procs
@ -85,6 +77,14 @@ class FixAtomSwap : public Fix {
class Compute *c_pe;
void options(int, char **);
int attempt_semi_grand();
int attempt_swap();
double energy_full();
int pick_semi_grand_atom();
int pick_i_swap_atom();
int pick_j_swap_atom();
void update_semi_grand_atoms_list();
void update_swap_atoms_list();
};
} // namespace LAMMPS_NS

501
src/MC/fix_mol_swap.cpp Normal file
View File

@ -0,0 +1,501 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "fix_mol_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 "improper.h"
#include "kspace.h"
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include "random_park.h"
#include "update.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
FixMolSwap::FixMolSwap(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), random(nullptr), c_pe(nullptr)
{
if (narg < 9) error->all(FLERR,"Illegal fix mol/swap command");
vector_flag = 1;
size_vector = 2;
global_freq = 1;
extvector = 0;
restart_global = 1;
time_depend = 1;
// parse args
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
ncycles = utils::inumeric(FLERR,arg[4],false,lmp);
itype = utils::inumeric(FLERR,arg[5],false,lmp);
jtype = utils::inumeric(FLERR,arg[6],false,lmp);
seed = utils::inumeric(FLERR,arg[7],false,lmp);
double temperature = utils::numeric(FLERR,arg[8],false,lmp);
// optional args
ke_flag = 1;
int iarg = 9;
while (iarg < narg) {
if (strcmp(arg[iarg],"ke") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix mol/swap command");
ke_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else error->all(FLERR,"Illegal fix mol/swap command");
}
// error check
if (nevery <= 0) error->all(FLERR,"Illegal fix mol/swap command");
if (ncycles < 0) error->all(FLERR,"Illegal fix mol/swap command");
if (itype == jtype) error->all(FLERR,"Illegal fix mol/swap command");
if (itype <= 0 || itype > atom->ntypes ||
jtype <= 0 || jtype > atom->ntypes)
error->all(FLERR,"Fix mol/swap atom types are invalid");
if (seed <= 0) error->all(FLERR,"Illegal fix mol/swap command");
if (temperature <= 0.0) error->all(FLERR,"Illegal fix mol/swap command");
if (ke_flag && atom->rmass)
error->all(FLERR,"Cannot conserve kinetic energy with fix mol/swap "
"unless per-type masses");
beta = 1.0/(force->boltz*temperature);
// random number generator, same for all procs
random = new RanPark(lmp,seed);
// set up reneighboring
force_reneighbor = 1;
next_reneighbor = update->ntimestep + 1;
// zero out counters
nswap_attempt = 0.0;
nswap_accept = 0.0;
// qflag = 1 if charges are defined
if (atom->q_flag) qflag = 1;
else qflag = 0;
// set comm size needed by this Fix
if (qflag) comm_forward = 2;
else comm_forward = 1;
}
/* ---------------------------------------------------------------------- */
FixMolSwap::~FixMolSwap()
{
delete random;
}
/* ---------------------------------------------------------------------- */
int FixMolSwap::setmask()
{
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixMolSwap::init()
{
// c_pe = compute used to calculate before/after potential energy
char *id_pe = (char *) "thermo_pe";
int ipe = modify->find_compute(id_pe);
c_pe = modify->compute[ipe];
// minmol = smallest molID with atoms of itype or jtype
// maxmol = largest molID with atoms of itype or jtype
// require: molID > 0, atoms in fix group
int *mask = atom->mask;
int *type = atom->type;
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
tagint minmol_me = MAXTAGINT;
tagint maxmol_me = 0;
for (int i = 0; i < nlocal; i++) {
if (molecule[i] == 0) continue;
if (!(mask[i] & groupbit)) continue;
if (molecule[i] < minmol_me) {
if (type[i] == itype || type[i] == jtype) minmol_me = molecule[i];
}
if (molecule[i] > maxmol_me) {
if (type[i] == itype || type[i] == jtype) maxmol_me = molecule[i];
}
}
MPI_Allreduce(&minmol_me,&minmol,1,MPI_LMP_TAGINT,MPI_MIN,world);
MPI_Allreduce(&maxmol_me,&maxmol,1,MPI_LMP_TAGINT,MPI_MAX,world);
// if ke_flag, check if itype/jtype masses are different
// if yes, pre-calcuate velocity scale factors that keep KE constant
// if no, turn ke_flag off
if (ke_flag) {
double *mass = atom->mass;
if (mass[itype] != mass[jtype]) {
i2j_vscale = sqrt(mass[itype]/mass[jtype]);
j2i_vscale = sqrt(mass[jtype]/mass[itype]);
} else ke_flag = 0;
}
// if qflag, check if all charges for itype are identical, ditto for jtype
// if not, cannot swap charges, issue warning
if (qflag) {
double *q = atom->q;
double iqone,jqone;
iqone = jqone = -BIG;
for (int i = 0; i < nlocal; i++) {
if (molecule[i] == 0) continue;
if (!(mask[i] & groupbit)) continue;
if (type[i] == itype) iqone = q[i];
if (type[i] == jtype) jqone = q[i];
}
MPI_Allreduce(&iqone,&iq,1,MPI_DOUBLE,MPI_MAX,world);
MPI_Allreduce(&jqone,&jq,1,MPI_DOUBLE,MPI_MAX,world);
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (molecule[i] == 0) continue;
if (!(mask[i] & groupbit)) continue;
if (type[i] == itype && q[i] != iq) flag = 1;
if (type[i] == jtype && q[i] != jq) flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall) qflag = 0;
if (!qflag && comm->me == 0)
error->warning(FLERR,"Cannot swap charges in fix mol/swap");
}
// check if all cutoffs involving itype and jtype are the same
// if not, reneighboring will be needed after swaps
double **cutsq = force->pair->cutsq;
unequal_cutoffs = false;
for (int ktype = 1; ktype <= atom->ntypes; ktype++)
if (cutsq[itype][ktype] != cutsq[jtype][ktype]) unequal_cutoffs = true;
}
/* ----------------------------------------------------------------------
perform ncycle Monte Carlo swaps
------------------------------------------------------------------------- */
void FixMolSwap::pre_exchange()
{
// just return if should not be called on this timestep
if (next_reneighbor != update->ntimestep) return;
// insure 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 molecule swaps
int naccept = 0;
for (int m = 0; m < ncycles; m++) naccept += attempt_swap();
// udpate MC stats
nswap_attempt += ncycles;
nswap_accept += naccept;
next_reneighbor = update->ntimestep + nevery;
}
/* ----------------------------------------------------------------------
attempt a swap of atom types within a random molecule
compare before/after energy and accept/reject the swap
------------------------------------------------------------------------- */
int FixMolSwap::attempt_swap()
{
// pre-swap energy
double energy_before = energy_stored;
// pick a random molecule
// swap properties of all its eligible itype & jtype atoms
tagint molID =
minmol + static_cast<tagint> (random->uniform() * (maxmol-minmol+1));
if (molID > maxmol) molID = maxmol;
int *mask = atom->mask;
double **v = atom->v;
double *q = atom->q;
int *type = atom->type;
tagint *molecule = atom->molecule;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (molecule[i] != molID) continue;
if (!(mask[i] & groupbit)) continue;
if (type[i] == itype) {
type[i] = jtype;
if (qflag) q[i] = jq;
if (ke_flag) {
v[i][0] *= i2j_vscale;
v[i][1] *= i2j_vscale;
v[i][2] *= i2j_vscale;
}
} else if (type[i] == jtype) {
type[i] = itype;
if (qflag) q[i] = iq;
if (ke_flag) {
v[i][0] *= j2i_vscale;
v[i][1] *= j2i_vscale;
v[i][2] *= j2i_vscale;
}
}
}
// 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);
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_fix(this);
}
// post-swap energy
double energy_after = energy_full();
// swap accepted, return 1
if (random->uniform() < exp(beta*(energy_before - energy_after))) {
energy_stored = energy_after;
return 1;
}
// swap not accepted, return 0
// restore the swapped itype & jtype atoms
// do not need to re-call comm->borders() and rebuild neighbor list
// since will be done on next cycle or in Verlet when this fix finishes
for (int i = 0; i < nlocal; i++) {
if (molecule[i] != molID) continue;
if (!(mask[i] & groupbit)) continue;
if (type[i] == itype) {
type[i] = jtype;
if (qflag) q[i] = jq;
if (ke_flag) {
v[i][0] *= i2j_vscale;
v[i][1] *= i2j_vscale;
v[i][2] *= i2j_vscale;
}
} else if (type[i] == jtype) {
type[i] = itype;
if (qflag) q[i] = iq;
if (ke_flag) {
v[i][0] *= j2i_vscale;
v[i][1] *= j2i_vscale;
v[i][2] *= j2i_vscale;
}
}
}
return 0;
}
/* ----------------------------------------------------------------------
compute system potential energy before or after a swap
------------------------------------------------------------------------- */
double FixMolSwap::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) modify->post_force(vflag);
update->eflag_global = update->ntimestep;
double total_energy = c_pe->compute_scalar();
return total_energy;
}
/* ---------------------------------------------------------------------- */
int FixMolSwap::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 (qflag) {
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 FixMolSwap::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 (qflag) {
for (i = first; i < last; i++) {
type[i] = static_cast<int> (buf[m++]);
q[i] = buf[m++];
}
} else {
for (i = first; i < last; i++)
type[i] = static_cast<int> (buf[m++]);
}
}
/* ----------------------------------------------------------------------
return acceptance ratio
------------------------------------------------------------------------- */
double FixMolSwap::compute_vector(int n)
{
if (n == 0) return nswap_attempt;
if (n == 1) return nswap_accept;
return 0.0;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixMolSwap::write_restart(FILE *fp)
{
int n = 0;
double list[6];
list[n++] = random->state();
list[n++] = ubuf(next_reneighbor).d;
list[n++] = nswap_attempt;
list[n++] = nswap_accept;
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 FixMolSwap::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
seed = static_cast<int> (list[n++]);
random->reset(seed);
next_reneighbor = (bigint) ubuf(list[n++]).i;
nswap_attempt = static_cast<int>(list[n++]);
nswap_accept = static_cast<int>(list[n++]);
bigint ntimestep_restart = (bigint) ubuf(list[n++]).i;
if (ntimestep_restart != update->ntimestep)
error->all(FLERR,"Must not reset timestep when restarting fix mol/swap");
}

122
src/MC/fix_mol_swap.h Normal file
View File

@ -0,0 +1,122 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(mol/swap,FixMolSwap);
// clang-format on
#else
#ifndef LMP_FIX_MOL_SWAP_H
#define LMP_FIX_MOL_SWAP_H
#include "fix.h"
namespace LAMMPS_NS {
class FixMolSwap : public Fix {
public:
FixMolSwap(class LAMMPS *, int, char **);
~FixMolSwap();
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);
void write_restart(FILE *);
void restart(char *);
private:
int nevery, ncycles, seed;
int itype, jtype;
double temperature;
int ke_flag; // 1 if kinetic energy is also swapped
double i2j_vscale; // scale factors for velocity to keep KE constant
double j2i_vscale;
int qflag; // 1 if charge is also swapped
double iq,jq; // charge values for all itype,jtype atoms
bool unequal_cutoffs; // 1 if itype and jtype have any different cutoffs
tagint minmol,maxmol; // range of mol IDs selected for swaps
double nswap_attempt; // cummulative stats on MC attempts and accepts
double nswap_accept;
double beta; // 1/kT
double energy_stored; // energy of current state as swaps are accepted
class RanPark *random;
class Compute *c_pe;
int attempt_swap();
double energy_full();
};
} // 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 atom/swap does not exist
Self-explanatory.
E: Must specify at least 2 types in fix atom/swap command
Self-explanatory.
E: Need nswaptypes mu values in fix atom/swap command
Self-explanatory.
E: Only 2 types allowed when not using semi-grand in fix atom/swap command
Self-explanatory.
E: Mu not allowed when not using semi-grand in fix atom/swap command
Self-explanatory.
E: Invalid atom type in fix atom/swap command
The atom type specified in the atom/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 atom/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.
*/