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

add more features to mol/swap, sync with atom/swap

Browse Source
This commit is contained in:
Steve Plimpton
2021-10-12 09:56:51 -06:00
parent eedd953258
commit 0a98ff3c38
4 changed files with 221 additions and 86 deletions
Download Patch File Download Diff File Expand all files Collapse all files

119
src/MC/fix_atom_swap.cpp
View File

@ -144,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;
@ -162,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");
@ -303,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();
@ -311,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();
@ -322,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) {
@ -352,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);
@ -364,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();
@ -376,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];
@ -387,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];
@ -433,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();
@ -445,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];
@ -463,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;
}
@ -499,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);

2
src/MC/fix_atom_swap.h
View File

@ -40,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

165
src/MC/fix_mol_swap.cpp
View File

@ -35,12 +35,14 @@
using namespace LAMMPS_NS;
using namespace FixConst;
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
FixMolSwap::FixMolSwap(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 9) error->all(FLERR,"Illegal fix mol/swap command");
if (narg < 9) error->all(FLERR,"Illegal fix mol/swap command");
vector_flag = 1;
size_vector = 2;
@ -58,13 +60,32 @@ FixMolSwap::FixMolSwap(LAMMPS *lmp, int narg, char **arg) :
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);
@ -79,12 +100,17 @@ FixMolSwap::FixMolSwap(LAMMPS *lmp, int narg, char **arg) :
// zero out counters
nswap_attempts = 0.0;
nswap_successes = 0.0;
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 (atom->q_flag) comm_forward = 2;
if (qflag) comm_forward = 2;
else comm_forward = 1;
}
@ -141,6 +167,54 @@ void FixMolSwap::init()
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
@ -152,7 +226,7 @@ void FixMolSwap::init()
}
/* ----------------------------------------------------------------------
perform Ncycle Monte Carlo swaps
perform ncycle Monte Carlo swaps
------------------------------------------------------------------------- */
void FixMolSwap::pre_exchange()
@ -178,11 +252,13 @@ void FixMolSwap::pre_exchange()
// attempt Ncycle molecule swaps
int nsuccess = 0;
for (int m = 0; m < ncycles; m++) nsuccess += attempt_swap();
int naccept = 0;
for (int m = 0; m < ncycles; m++) naccept += attempt_swap();
nswap_attempts += ncycles;
nswap_successes += nsuccess;
// udpate MC stats
nswap_attempt += ncycles;
nswap_accept += naccept;
next_reneighbor = update->ntimestep + nevery;
}
@ -199,13 +275,15 @@ int FixMolSwap::attempt_swap()
double energy_before = energy_stored;
// pick a random molecule
// swap all of its eligible itype & jtype atoms
// 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;
@ -213,8 +291,23 @@ int FixMolSwap::attempt_swap()
for (int i = 0; i < nlocal; i++) {
if (molecule[i] != molID) continue;
if (!(mask[i] & groupbit)) continue;
if (type[i] == itype) type[i] = jtype;
else if (type[i] == jtype) type[i] = itype;
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
@ -236,23 +329,37 @@ int FixMolSwap::attempt_swap()
double energy_after = energy_full();
// swap accepted
// swap accepted, return 1
if (random->uniform() < exp(beta*(energy_before - energy_after))) {
energy_stored = energy_after;
return 1;
}
// swap not accepted
// 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 is done
// since will be done on next cycle or in Verlet when this fix finishes
} else {
for (int i = 0; i < nlocal; i++) {
if (molecule[i] != molID) continue;
if (!(mask[i] & groupbit)) continue;
if (type[i] == itype) type[i] = jtype;
else if (type[i] == jtype) type[i] = itype;
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;
}
}
}
@ -300,7 +407,7 @@ int FixMolSwap::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*
m = 0;
if (atom->q_flag) {
if (qflag) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = type[j];
@ -328,7 +435,7 @@ void FixMolSwap::unpack_forward_comm(int n, int first, double *buf)
m = 0;
last = first + n;
if (atom->q_flag) {
if (qflag) {
for (i = first; i < last; i++) {
type[i] = static_cast<int> (buf[m++]);
q[i] = buf[m++];
@ -345,8 +452,8 @@ void FixMolSwap::unpack_forward_comm(int n, int first, double *buf)
double FixMolSwap::compute_vector(int n)
{
if (n == 0) return nswap_attempts;
if (n == 1) return nswap_successes;
if (n == 0) return nswap_attempt;
if (n == 1) return nswap_accept;
return 0.0;
}
@ -360,8 +467,8 @@ void FixMolSwap::write_restart(FILE *fp)
double list[6];
list[n++] = random->state();
list[n++] = ubuf(next_reneighbor).d;
list[n++] = nswap_attempts;
list[n++] = nswap_successes;
list[n++] = nswap_attempt;
list[n++] = nswap_accept;
list[n++] = ubuf(update->ntimestep).d;
if (comm->me == 0) {
@ -385,8 +492,8 @@ void FixMolSwap::restart(char *buf)
next_reneighbor = (bigint) ubuf(list[n++]).i;
nswap_attempts = static_cast<int>(list[n++]);
nswap_successes = static_cast<int>(list[n++]);
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)

21
src/MC/fix_mol_swap.h
View File

@ -42,12 +42,21 @@ class FixMolSwap : public Fix {
int itype, jtype;
double temperature;
bool unequal_cutoffs;
tagint minmol,maxmol;
double nswap_attempts;
double nswap_successes;
double beta;
double energy_stored;
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;