git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12992 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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]);
+  seed = force->inumeric(FLERR,arg[5]); 
-  atom_swap_jtype = force->inumeric(FLERR,arg[6]);
+  double temperature = force->numeric(FLERR,arg[6]);
-  seed = force->inumeric(FLERR,arg[7]); 
-  double temperature = force->numeric(FLERR,arg[8]);
   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_atom_list = NULL;
-  local_swap_jatom_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) ||
+  if (nswaptypes < 2)
-      (atom_swap_jtype <= 0 || atom_swap_jtype > atom->ntypes))
+    error->all(FLERR,"Must specify at least 2 types in fix atom/swap command");
-    error->all(FLERR,"Invalid atom type 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) {
+  memory->create(qtype,nswaptypes,"atom/swap:qtype");
-    bool firsti = true;
+  for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++) {
-    bool firstj = true;
+    if (atom->q_flag) {
-    for (int i = 0; i < atom->nlocal; i++) {
+      bool first = true;
-      if (type[i] == atom_swap_itype) {
+      for (int i = 0; i < atom->nlocal; i++) {
-        if (firsti) qitype = atom->q[i];
+        if (type[i] == type_list[iswaptype]) {
-        firsti = false;
+          if (first) qtype[iswaptype] = atom->q[i];
-        if (qitype != atom->q[i])
+          first = false;
-          error->all(FLERR,"All atoms of a swapped type must have the same charge.");
+          if (qtype[iswaptype] != atom->q[i])
-      } else if (type[i] == atom_swap_jtype) {
+            error->all(FLERR,"All atoms of a swapped type must have the same charge.");
-        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.");
       }
     }
   }
   
-  sqrt_mass_ratio_ij = sqrt(atom->mass[atom_swap_itype]/atom->mass[atom_swap_jtype]);
+  memory->create(sqrt_mass_ratio,atom->ntypes+1,atom->ntypes+1,"atom/swap:sqrt_mass_ratio");
-  sqrt_mass_ratio_ji = sqrt(atom->mass[atom_swap_jtype]/atom->mass[atom_swap_itype]);
+  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++)
+  for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++)
-    if (cutsq[atom_swap_itype][ktype] != cutsq[atom_swap_jtype][ktype])
+    for (int jswaptype = 0; jswaptype < nswaptypes; jswaptype++)
-      unequal_cutoffs = true;
+      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<int> (nswaptypes*random_equal->uniform());
+    jtype = type_list[jswaptype];
+    itype = atom->type[i];
+    while (itype == jtype) {
+      jswaptype = static_cast<int> (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;
+  if (i >= 0) {
-  int j = -1;
+    atom->type[i] = jtype;
-  
+    if (atom->q_flag) atom->q[i] = qtype[1];
-  if (semi_grand_flag) {
+  }
-  
+  if (j >= 0) {
-    if (random_equal->uniform() < 0.5) {
+    atom->type[j] = itype;
-      if (niswap == 0) return 0;
+    if (atom->q_flag) atom->q[j] = qtype[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 (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][0] *= sqrt_mass_ratio[itype][jtype];
-        atom->v[i][1] *= sqrt_mass_ratio_ij;
+        atom->v[i][1] *= sqrt_mass_ratio[itype][jtype];
-        atom->v[i][2] *= sqrt_mass_ratio_ij;
+        atom->v[i][2] *= sqrt_mass_ratio[itype][jtype];
       }
       if (j >= 0) {
-        atom->v[j][0] *= sqrt_mass_ratio_ji;
+        atom->v[j][0] *= sqrt_mass_ratio[jtype][itype];
-        atom->v[j][1] *= sqrt_mass_ratio_ji;
+        atom->v[j][1] *= sqrt_mass_ratio[jtype][itype];
-        atom->v[j][2] *= sqrt_mass_ratio_ji;
+        atom->v[j][2] *= sqrt_mass_ratio[jtype][itype];
       }
     }
     return 1;
   } else {
     if (i >= 0) {
-      atom->type[i] = atom_swap_itype;
+      atom->type[i] =  type_list[0];
-      if (atom->q_flag) atom->q[i] = qitype;
+      if (atom->q_flag) atom->q[i] = qtype[0];
     }
     if (j >= 0) {
-      atom->type[j] = atom_swap_jtype;
+      atom->type[j] =  type_list[1];
-      if (atom->q_flag) atom->q[j] = qjtype;
+      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<int> (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++;
         }

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;