diff --git a/src/MC/fix_atom_swap.cpp b/src/MC/fix_atom_swap.cpp index 8ac5673f6a..fb2c43815e 100644 --- a/src/MC/fix_atom_swap.cpp +++ b/src/MC/fix_atom_swap.cpp @@ -56,7 +56,7 @@ using namespace MathConst; FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (narg < 9) error->all(FLERR,"Illegal fix atom/swap command"); + if (narg < 10) error->all(FLERR,"Illegal fix atom/swap command"); dynamic_group_allow = 1; @@ -71,18 +71,19 @@ FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) : nevery = force->inumeric(FLERR,arg[3]); ncycles = force->inumeric(FLERR,arg[4]); - atom_swap_itype = force->inumeric(FLERR,arg[5]); - atom_swap_jtype = force->inumeric(FLERR,arg[6]); - seed = force->inumeric(FLERR,arg[7]); - double temperature = force->numeric(FLERR,arg[8]); + seed = force->inumeric(FLERR,arg[5]); + double temperature = force->numeric(FLERR,arg[6]); beta = 1.0/(force->boltz*temperature); if (ncycles < 0) error->all(FLERR,"Illegal fix atom/swap command"); if (seed <= 0) error->all(FLERR,"Illegal fix atom/swap command"); + memory->create(delta_mu,atom->ntypes+1,"atom/swap:delta_mu"); + for (int i = 1; i <= atom->ntypes; i++) delta_mu[i] = 0.0; + // read options from end of input line - options(narg-9,&arg[9]); + options(narg-7,&arg[7]); // random number generator, same for all procs @@ -99,8 +100,7 @@ FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) : nswap_successes = 0.0; atom_swap_nmax = 0; - local_swap_iatom_list = NULL; - local_swap_jatom_list = NULL; + local_swap_atom_list = NULL; // set comm size needed by this Fix @@ -146,6 +146,26 @@ void FixAtomSwap::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"yes") == 0) semi_grand_flag = 1; else error->all(FLERR,"Illegal fix atom/swap command"); iarg += 2; + } else if (strcmp(arg[iarg],"types") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix atom/swap command"); + iarg++; + memory->create(type_list,atom->ntypes,"atom/swap:type_list"); + nswaptypes = 0; + while (iarg < narg) { + if (isalpha(arg[iarg][0])) break; + type_list[nswaptypes++] = force->numeric(FLERR,arg[iarg]); + iarg++; + } + } else if (strcmp(arg[iarg],"delta_mu") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix atom/swap command"); + iarg++; + ndeltamutypes = 0; + while (iarg < narg) { + if (isalpha(arg[iarg][0])) break; + delta_mu[ndeltamutypes+2] = force->numeric(FLERR,arg[iarg]); + ndeltamutypes++; + iarg++; + } } else error->all(FLERR,"Illegal fix atom/swap command"); } } @@ -154,6 +174,10 @@ void FixAtomSwap::options(int narg, char **arg) FixAtomSwap::~FixAtomSwap() { + memory->destroy(type_list); + memory->destroy(delta_mu); + memory->destroy(qtype); + memory->destroy(sqrt_mass_ratio); if (regionflag) delete [] idregion; delete random_equal; } @@ -177,39 +201,53 @@ void FixAtomSwap::init() int *type = atom->type; - if ((atom_swap_itype <= 0 || atom_swap_itype > atom->ntypes) || - (atom_swap_jtype <= 0 || atom_swap_jtype > atom->ntypes)) - error->all(FLERR,"Invalid atom type in fix atom/swap command"); + if (nswaptypes < 2) + error->all(FLERR,"Must specify at least 2 types in fix atom/swap command"); + + if (semi_grand_flag) { + if (atom->ntypes-1 != ndeltamutypes) + error->all(FLERR,"Need ntypes-1 delta_mu values in fix atom/swap command"); + } else { + if (nswaptypes != 2) + error->all(FLERR,"Only 2 types allowed when not using semi-grand in fix atom/swap command"); + if (ndeltamutypes != 0) + error->all(FLERR,"Delta_mu not allowed when not using semi-grand in fix atom/swap command"); + } + + 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 atom/swap command"); - if (atom->q_flag) { - bool firsti = true; - bool firstj = true; - for (int i = 0; i < atom->nlocal; i++) { - if (type[i] == atom_swap_itype) { - if (firsti) qitype = atom->q[i]; - firsti = false; - if (qitype != atom->q[i]) - error->all(FLERR,"All atoms of a swapped type must have the same charge."); - } else if (type[i] == atom_swap_jtype) { - if (firstj) qjtype = atom->q[i]; - firstj = false; - if (qjtype != atom->q[i]) - error->all(FLERR,"All atoms of a swapped type must have the same charge."); + memory->create(qtype,nswaptypes,"atom/swap:qtype"); + for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++) { + if (atom->q_flag) { + bool first = true; + for (int i = 0; i < atom->nlocal; i++) { + if (type[i] == type_list[iswaptype]) { + if (first) qtype[iswaptype] = atom->q[i]; + first = false; + if (qtype[iswaptype] != atom->q[i]) + error->all(FLERR,"All atoms of a swapped type must have the same charge."); + } } } } - sqrt_mass_ratio_ij = sqrt(atom->mass[atom_swap_itype]/atom->mass[atom_swap_jtype]); - sqrt_mass_ratio_ji = sqrt(atom->mass[atom_swap_jtype]/atom->mass[atom_swap_itype]); + memory->create(sqrt_mass_ratio,atom->ntypes+1,atom->ntypes+1,"atom/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 ktype = 1; ktype <= atom->ntypes; ktype++) - if (cutsq[atom_swap_itype][ktype] != cutsq[atom_swap_jtype][ktype]) - unequal_cutoffs = true; + 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 @@ -249,10 +287,15 @@ void FixAtomSwap::pre_exchange() neighbor->build(); energy_stored = energy_full(); - update_swap_atoms_list(); - + int nsuccess = 0; - for (int i = 0; i < ncycles; i++) nsuccess += attempt_swap(); + if (semi_grand_flag) { + update_semi_grand_atoms_list(); + for (int i = 0; i < ncycles; i++) nsuccess += attempt_semi_grand(); + } else { + update_swap_atoms_list(); + for (int i = 0; i < ncycles; i++) nsuccess += attempt_swap(); + } nswap_attempts += ncycles; nswap_successes += nsuccess; @@ -260,49 +303,93 @@ void FixAtomSwap::pre_exchange() next_reneighbor = update->ntimestep + nevery; } +/* ---------------------------------------------------------------------- +------------------------------------------------------------------------- */ + +int FixAtomSwap::attempt_semi_grand() +{ + if (nswap == 0) return 0; + + double energy_before = energy_stored; + + int itype,jtype,jswaptype; + double qtmp; + + int i = pick_semi_grand_atom(); + if (i >= 0) { + jswaptype = static_cast (nswaptypes*random_equal->uniform()); + jtype = type_list[jswaptype]; + itype = atom->type[i]; + while (itype == jtype) { + jswaptype = static_cast (nswaptypes*random_equal->uniform()); + jtype = type_list[jswaptype]; + } + atom->type[i] = jtype; + if (atom->q_flag) { + qtmp = atom->q[i]; + atom->q[i] = qtype[jswaptype]; + } + } + + 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(); + } else { + comm->forward_comm_fix(this); + } + + double energy_after = energy_full(); + + if (random_equal->uniform() < + exp(-beta*(energy_after - energy_before + + delta_mu[jtype] - delta_mu[itype]))) { + update_semi_grand_atoms_list(); + energy_stored = energy_after; + if (conserve_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]; + } + } + return 1; + } else { + if (i >= 0) { + atom->type[i] = itype; + if (atom->q_flag) atom->q[i] = qtmp; + } + energy_stored = energy_before; + } + return 0; +} + + /* ---------------------------------------------------------------------- ------------------------------------------------------------------------- */ int FixAtomSwap::attempt_swap() { + if ((niswap == 0) || (njswap == 0)) return 0; + double energy_before = energy_stored; + + int i = pick_i_swap_atom(); + int j = pick_j_swap_atom(); + int itype = type_list[0]; + int jtype = type_list[1]; - int i = -1; - int j = -1; - - if (semi_grand_flag) { - - if (random_equal->uniform() < 0.5) { - if (niswap == 0) return 0; - i = pick_i_swap_atom(); - if (i >= 0) { - atom->type[i] = atom_swap_jtype; - if (atom->q_flag) atom->q[i] = qjtype; - } - } else { - if (njswap == 0) return 0; - j = pick_j_swap_atom(); - if (j >= 0) { - atom->type[j] = atom_swap_itype; - if (atom->q_flag) atom->q[j] = qitype; - } - } - - } else { - - if ((niswap == 0) || (njswap == 0)) return 0; - i = pick_i_swap_atom(); - j = pick_j_swap_atom(); - - if (i >= 0) { - atom->type[i] = atom_swap_jtype; - if (atom->q_flag) atom->q[i] = qjtype; - } - if (j >= 0) { - atom->type[j] = atom_swap_itype; - if (atom->q_flag) atom->q[j] = qitype; - } - + if (i >= 0) { + atom->type[i] = jtype; + if (atom->q_flag) atom->q[i] = qtype[1]; + } + if (j >= 0) { + atom->type[j] = itype; + if (atom->q_flag) atom->q[j] = qtype[0]; } if (unequal_cutoffs) { @@ -325,28 +412,28 @@ int FixAtomSwap::attempt_swap() energy_stored = energy_after; if (conserve_ke_flag) { if (i >= 0) { - atom->v[i][0] *= sqrt_mass_ratio_ij; - atom->v[i][1] *= sqrt_mass_ratio_ij; - atom->v[i][2] *= sqrt_mass_ratio_ij; + 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_ji; - atom->v[j][1] *= sqrt_mass_ratio_ji; - atom->v[j][2] *= sqrt_mass_ratio_ji; + 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]; } } return 1; } else { if (i >= 0) { - atom->type[i] = atom_swap_itype; - if (atom->q_flag) atom->q[i] = qitype; + atom->type[i] = type_list[0]; + if (atom->q_flag) atom->q[i] = qtype[0]; } if (j >= 0) { - atom->type[j] = atom_swap_jtype; - if (atom->q_flag) atom->q[j] = qjtype; + atom->type[j] = type_list[1]; + if (atom->q_flag) atom->q[j] = qtype[1]; } energy_stored = energy_before; - } + } return 0; } @@ -380,6 +467,22 @@ double FixAtomSwap::energy_full() /* ---------------------------------------------------------------------- ------------------------------------------------------------------------- */ +int FixAtomSwap::pick_semi_grand_atom() +{ + int i = -1; + int iwhichglobal = static_cast (nswap*random_equal->uniform()); + if ((iwhichglobal >= nswap_before) && + (iwhichglobal < nswap_before + nswap_local)) { + int iwhichlocal = iwhichglobal - nswap_before; + i = local_swap_atom_list[iwhichlocal]; + } + + return i; +} + +/* ---------------------------------------------------------------------- +------------------------------------------------------------------------- */ + int FixAtomSwap::pick_i_swap_atom() { int i = -1; @@ -409,6 +512,55 @@ int FixAtomSwap::pick_j_swap_atom() return j; } +/* ---------------------------------------------------------------------- + update the list of gas atoms +------------------------------------------------------------------------- */ + +void FixAtomSwap::update_semi_grand_atoms_list() +{ + int nlocal = atom->nlocal; + int *type = atom->type; + double **x = atom->x; + + if (nlocal > atom_swap_nmax) { + memory->sfree(local_swap_atom_list); + atom_swap_nmax = atom->nmax; + local_swap_atom_list = (int *) memory->smalloc(atom_swap_nmax*sizeof(int), + "MCSWAP:local_swap_atom_list"); + } + + nswap_local = 0; + + if (regionflag) { + + for (int i = 0; i < nlocal; i++) { + if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) { + if (atom->mask[i] & groupbit) { + if (type[i] == type_list[0]) { + local_swap_atom_list[nswap_local] = i; + nswap_local++; + } + } + } + } + + } else { + for (int i = 0; i < nlocal; i++) { + if (atom->mask[i] & groupbit) { + if (type[i] == type_list[0]) { + local_swap_atom_list[nswap_local] = i; + nswap_local++; + } + } + } + } + + MPI_Allreduce(&nswap_local,&nswap,1,MPI_INT,MPI_SUM,world); + MPI_Scan(&nswap_local,&nswap_before,1,MPI_INT,MPI_SUM,world); + nswap_before -= nswap_local; +} + + /* ---------------------------------------------------------------------- update the list of gas atoms ------------------------------------------------------------------------- */ @@ -437,10 +589,10 @@ void FixAtomSwap::update_swap_atoms_list() for (int i = 0; i < nlocal; i++) { if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) { if (atom->mask[i] & groupbit) { - if (type[i] == atom_swap_itype) { + if (type[i] == type_list[0]) { local_swap_iatom_list[niswap_local] = i; niswap_local++; - } else if (type[i] == atom_swap_jtype) { + } else if (type[i] == type_list[1]) { local_swap_jatom_list[njswap_local] = i; njswap_local++; } @@ -451,10 +603,10 @@ void FixAtomSwap::update_swap_atoms_list() } else { for (int i = 0; i < nlocal; i++) { if (atom->mask[i] & groupbit) { - if (type[i] == atom_swap_itype) { + if (type[i] == type_list[0]) { local_swap_iatom_list[niswap_local] = i; niswap_local++; - } else if (type[i] == atom_swap_jtype) { + } else if (type[i] == type_list[1]) { local_swap_jatom_list[njswap_local] = i; njswap_local++; } diff --git a/src/MC/fix_atom_swap.h b/src/MC/fix_atom_swap.h index 50eaa63f71..c94f57ac86 100644 --- a/src/MC/fix_atom_swap.h +++ b/src/MC/fix_atom_swap.h @@ -32,10 +32,13 @@ 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 *); @@ -45,17 +48,24 @@ class FixAtomSwap : public Fix { void restart(char *); private: - int atom_swap_itype,atom_swap_jtype,nevery,seed; + int nevery,seed; int conserve_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 swap atoms on all procs + 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 nswap; // # of swap atoms on all procs + int nswap_local; // # of swap atoms on this proc + int nswap_before; // # of swap atoms on procs < this proc int regionflag; // 0 = anywhere in box, 1 = specific region int iregion; // swap region char *idregion; // swap region id + int nswaptypes,ndeltamutypes; + int *type_list; + double *delta_mu; + double nswap_attempts; double nswap_successes; @@ -63,11 +73,12 @@ class FixAtomSwap : public Fix { int atom_swap_nmax; double beta; - double qitype,qjtype; + double *qtype; double energy_stored; - double sqrt_mass_ratio_ij,sqrt_mass_ratio_ji; + double **sqrt_mass_ratio; int *local_swap_iatom_list; int *local_swap_jatom_list; + int *local_swap_atom_list; class RanPark *random_equal;