make new/incremental quantum trigger more robust, including fix GCMC and atom/swap

src/MC/fix_atom_swap.cpp

@@ -103,6 +103,8 @@ FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) :
 
   // zero out counters
 
+  mc_active = 0;
+  
   nswap_attempts = 0.0;
   nswap_successes = 0.0;
 
@@ -304,6 +306,8 @@ void FixAtomSwap::pre_exchange()
 
   if (next_reneighbor != update->ntimestep) return;
 
+  mc_active = 1;
+  
   // ensure current system is ready to compute energy
 
   if (domain->triclinic) domain->x2lamda(atom->nlocal);
@@ -336,6 +340,8 @@ void FixAtomSwap::pre_exchange()
   nswap_successes += nsuccess;
 
   next_reneighbor = update->ntimestep + nevery;
+
+  mc_active = 0;
 }
 
 /* ----------------------------------------------------------------------
@@ -815,3 +821,18 @@ void FixAtomSwap::restart(char *buf)
   if (ntimestep_restart != update->ntimestep)
     error->all(FLERR, "Must not reset timestep when restarting fix atom/swap");
 }
+
+/* ----------------------------------------------------------------------
+   extract variable which stores whether MC is active or not
+     active = MC moves are taking place
+     not active = normal MD is taking place
+------------------------------------------------------------------------- */
+
+void *FixAtomSwap::extract(const char *name, int &dim)
+{
+  if (strcmp(name,"mc_active") == 0) {
+    dim = 0;
+    return (void *) &mc_active;
+  }
+  return nullptr;
+}

src/MC/fix_atom_swap.h

@@ -37,6 +37,7 @@ class FixAtomSwap : public Fix {
   double memory_usage() override;
   void write_restart(FILE *) override;
   void restart(char *) override;
+  void *extract(const char *, int &) override;
 
  private:
   int nevery, seed;

src/MC/fix_gcmc.cpp

@@ -220,6 +220,8 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
 
   // zero out counters
 
+  mc_active = 0;
+
   ntranslation_attempts = 0.0;
   ntranslation_successes = 0.0;
   nrotation_attempts = 0.0;
@@ -718,6 +720,8 @@ void FixGCMC::pre_exchange()
 
   if (next_reneighbor != update->ntimestep) return;
 
+  mc_active = 1;
+  
   xlo = domain->boxlo[0];
   xhi = domain->boxhi[0];
   ylo = domain->boxlo[1];
@@ -795,7 +799,10 @@ void FixGCMC::pre_exchange()
       }
     }
   }
+  
   next_reneighbor = update->ntimestep + nevery;
+
+  mc_active = 0;
 }
 
 /* ----------------------------------------------------------------------
@@ -2583,11 +2590,18 @@ void FixGCMC::grow_molecule_arrays(int nmolatoms) {
 
 
 /* ----------------------------------------------------------------------
+   extract variable which stores whether MC is active or not
+     active = MC moves are taking place
+     not active = normal MD is taking place
    extract variable which stores index of exclusion group
 ------------------------------------------------------------------------- */
 
 void *FixGCMC::extract(const char *name, int &dim)
 {
+  if (strcmp(name,"mc_active") == 0) {
+    dim = 0;
+    return (void *) &mc_active;
+  }
   if (strcmp(name,"exclusion_group") == 0) {
     dim = 0;
     return (void *) &exclusion_group;

src/MC/fix_gcmc.h

@@ -74,6 +74,8 @@ class FixGCMC : public Fix {
   double ninsertion_attempts;
   double ninsertion_successes;
 
+  int mc_active;              // 1 during MC trials, otherwise 0
+  
   int gcmc_nmax;
   int max_region_attempts;
   double gas_mass;

src/MDI/fix_mdi_qm.cpp

@@ -58,8 +58,8 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
   addflag = 1;
   every = 1;
   connectflag = 1;
-  exclude = 0;
+  mcflag = 0;
-  id_exclude = nullptr;
+  id_mcfix = nullptr;
   elements = nullptr;
 
   int iarg = 3;
@@ -96,11 +96,11 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
       else
         error->all(FLERR, "Illegal fix mdi/qm command");
       iarg += 2;
-    } else if (strcmp(arg[iarg],"exclude") == 0) {
+    } else if (strcmp(arg[iarg],"mc") == 0) {
       if (iarg+2 > narg) error->all(FLERR, "Illegal fix mdi/qm command");
-      exclude = 1;
+      mcflag = 1;
-      delete[] id_exclude;
+      delete[] id_mcfix;
-      id_exclude = utils::strdup(arg[iarg+1]);
+      id_mcfix = utils::strdup(arg[iarg+1]);
       iarg += 2;
 
     } else if (strcmp(arg[iarg], "elements") == 0) {
@@ -230,7 +230,7 @@ FixMDIQM::~FixMDIQM()
 
   // clean up
 
-  delete[] id_exclude;
+  delete[] id_mcfix;
   delete[] elements;
   
   memory->destroy(qmIDs);
@@ -324,48 +324,90 @@ void FixMDIQM::init()
     MPI_Bcast(&types_exists, 1, MPI_INT, 0, world);
   }
 
-  // extract pointer to exclusion_group variable from id_exclude
+  // extract pointers to MC variables from id_mcfix
-  // exclusion_group = index of a group the Fix defines
+  // mc_active = 1/0 for MC moves on/off (required)
-  // the group flags atoms to be excluded from QM calculation
+  // exclusion_group = index of a group the Fix defines to exclude atoms (optional)
 
-  if (exclude) {
+  if (mcflag) {
-    Fix *f_exclude = modify->get_fix_by_id(id_exclude);
+    Fix *f_mc = modify->get_fix_by_id(id_mcfix);
-    if (!f_exclude) error->all(FLERR,"Fix mdi/qm could not find exclude fix ID {}", id_exclude);
+    if (!f_mc) error->all(FLERR,"Fix mdi/qm could not find Monte Carlo fix ID {}", id_mcfix);
     int dim;
-    exclusion_group_ptr = (int *) f_exclude->extract("exclusion_group", dim);
+    mc_active_ptr = (int *) f_mc->extract("mc_active", dim);
-    if (!exclusion_group_ptr || dim != 0)
+    if (!mc_active_ptr || dim != 0)
-      error->all(FLERR,"Fix mdi/qm could not query exclude_group of fix ID {}", id_exclude);
+      error->all(FLERR,"Fix mdi/qm could not query mc_active from Monte Carlo fix ID {}", id_mcfix);
+    exclusion_group_ptr = (int *) f_mc->extract("exclusion_group", dim);
   }
 
-  // nqm = # of QM atoms
+  // determine whether a new vs incremental QM calc is needed
-  // if hasn't changed, then continuing previous simulation
+  // new when first run or when
-  // else setup new simulation with MDI engine
+  //   atom count or elements/types or box has changed between runs
+  // otherwise incremental = subsequent run of same system
+  
+  int new_system = 0;
 
+  // check if quantum atom count has changed
+  // on first run, old count is 0
+  
+  int nqm_old = nqm;
   nqm = set_nqm();
   
-  if (nqm != nqm_last) {
+  if (nqm != nqm_old) {
     if (nqm > max_nqm) reallocate();
-
-    // create list of QM atom IDs in sorted order
-    // set mapping of my local atoms into global list
-
     create_qm_list();
     set_qm2owned();
-
+    new_system = 1;
-    // initial send of info to MDI engine
+  }
-    // values that often won't change for AIMD simulations
-    // if not sending elements or types, assume engine initialized itself
-    
-    send_box();
-    send_natoms();
   
-    if (elements && elements_exists) {
+  // check if box has changed
+
+  if (new_system) set_box();
+  else {
+    double old_cell[9],old_cell_displ[3];
+    memcpy(old_cell,qm_cell,9*sizeof(double));
+    memcpy(old_cell_displ,qm_cell_displ,3*sizeof(double));
+    set_box();
+    for (int i = 0; i < 9; i++)
+      if (qm_cell[i] != old_cell[i]) new_system = 1;
+    for (int i = 0; i < 3; i++)
+      if (qm_cell_displ[i] != old_cell_displ[i]) new_system = 1;
+  }
+
+  // check if atom elements or types have changed
+
+  if (elements && elements_exists) {
+    if (new_system) set_eqm();
+    else {
+      int *eqm_old;
+      memory->create(eqm_old,nqm,"mdi/qm:eqm_old");
+      memcpy(eqm_old,eqm,nqm*sizeof(int));
       set_eqm();
-      send_elements();
+      for (int i = 0; i < nqm; i++)
-    } else if (types_exists) {
+        if (eqm[i] != eqm_old[i]) new_system = 1;
-      set_tqm();
+      memory->destroy(eqm_old);
-      send_types();
     }
+  } else if (types_exists) {
+    if (new_system) set_tqm();
+    else {
+      int *tqm_old;
+      memory->create(tqm_old,nqm,"mdi/qm:tqm_old");
+      memcpy(tqm_old,tqm,nqm*sizeof(int));
+      set_tqm();
+      for (int i = 0; i < nqm; i++)
+        if (tqm[i] != tqm_old[i]) new_system = 1;
+      memory->destroy(tqm_old);
+    }
+  }
+
+  // new system, so send setup info to MDI engine
+  // values that often won't change for AIMD simulations
+  // if not sending elements or types, assume engine initialized itself
+
+  if (new_system) {
+    send_natoms();
+    send_box();
+    if (elements && elements_exists) send_elements();
+    else if (types_exists) send_types();
+    nqm_last = nqm;
   }
 }
 
@@ -408,17 +450,27 @@ void FixMDIQM::post_force(int vflag)
     memory->create(array_atom,nmax,3,"mdi/qm:array_atom");
   }
 
-  // check count of QM atoms, could be different due to exclusions
+  // determine whether a new vs incremental QM calc is needed
-  // if so, need to reset QM data structures
+  // new when atom count has changed (deposit, evaporate)
+  //   or when MC fix is active (insertion, deletion, large moves)
+  // incremental when a system is slowly evolving (AIMD)
   
-  nqm = set_nqm();
+  int new_system = 0;
+  if (nqm != nqm_last) new_system = 1;
+  else if (mcflag && *mc_active_ptr) new_system = 1;
   
-  if (nqm != nqm_last) {
+  // send new system info to MDI engine: atom count and elements/types
-    nqm_last = nqm;
+  // reset QM data structures if atom count has changed
+  
+  if (new_system) {
+    nqm = set_nqm();
     if (nqm > max_nqm) reallocate();
     create_qm_list();
     set_qm2owned();
+
     send_natoms();
+    set_box();
+    send_box();
     if (elements && elements_exists) {
       set_eqm();
       send_elements();
@@ -426,11 +478,16 @@ void FixMDIQM::post_force(int vflag)
       set_tqm();
       send_types();
     }
-  }
 
+    nqm_last = nqm;
+
+  // incremental system
   // if simulation box dynamically changes, send current box to MDI engine
 
-  if (domain->box_change_size || domain->box_change_shape) send_box();
+  } else if (domain->box_change_size || domain->box_change_shape) {
+    set_box();
+    send_box();
+  }
 
   // send current coords of QM atoms to MDI engine
   
@@ -476,9 +533,9 @@ void FixMDIQM::post_force(int vflag)
   }
 
   // array_atom = fix output for peratom QM forces
-  // if exclude, some atoms may not be QM atoms, so zero array_atom first
+  // if nexclude, some atoms are not QM atoms, zero array_atom first
 
-  if (exclude) {
+  if (nexclude) {
     int nlocal = atom->nlocal;
     for (int i = 0; i < nlocal; i++) {
       array_atom[i][0] = 0.0;
@@ -626,7 +683,7 @@ void FixMDIQM::reallocate()
 
 /* ----------------------------------------------------------------------
    ncount = # of QM atoms
-   can be less than all atoms if exclude flag is set   
+   can be less than all atoms if MC flag is set   
    return ncount to set nqm
 ------------------------------------------------------------------------- */
 
@@ -635,7 +692,7 @@ int FixMDIQM::set_nqm()
   int ncount = atom->natoms;
   nexclude = 0;
 
-  if (exclude) {
+  if (mcflag && exclusion_group_ptr) {
     int nexclude_mine = 0;
 
     int exclusion_group = *exclusion_group_ptr;
@@ -756,6 +813,29 @@ void FixMDIQM::set_qm2owned()
 
 /* ---------------------------------------------------------------------- */
 
+void FixMDIQM::set_box()
+{
+  qm_cell_displ[0] = domain->boxlo[0] * lmp2mdi_length;
+  qm_cell_displ[1] = domain->boxlo[1] * lmp2mdi_length;
+  qm_cell_displ[2] = domain->boxlo[2] * lmp2mdi_length;
+
+  qm_cell[0] = domain->boxhi[0] - domain->boxlo[0];
+  qm_cell[1] = 0.0;
+  qm_cell[2] = 0.0;
+  qm_cell[3] = domain->xy;
+  qm_cell[4] = domain->boxhi[1] - domain->boxlo[1];
+  qm_cell[5] = 0.0;
+  qm_cell[6] = domain->xz;
+  qm_cell[7] = domain->yz;
+  qm_cell[8] = domain->boxhi[2] - domain->boxlo[2];
+
+  // convert cell units to bohr
+
+  for (int icell = 0; icell < 9; icell++) qm_cell[icell] *= lmp2mdi_length;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixMDIQM::set_xqm()
 {
   for (int i = 0; i < nqm; i++) {
@@ -883,35 +963,17 @@ void FixMDIQM::send_elements()
 void FixMDIQM::send_box()
 {
   int ierr;
-  double cell[9];
 
   if (celldispl_exists) {
     ierr = MDI_Send_command(">CELL_DISPL", mdicomm);
     if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command");
-    cell[0] = domain->boxlo[0] * lmp2mdi_length;
+    ierr = MDI_Send(qm_cell_displ, 3, MDI_DOUBLE, mdicomm);
-    cell[1] = domain->boxlo[1] * lmp2mdi_length;
-    cell[2] = domain->boxlo[2] * lmp2mdi_length;
-    ierr = MDI_Send(cell, 3, MDI_DOUBLE, mdicomm);
     if (ierr) error->all(FLERR, "MDI: >CELL_DISPL data");
   }
 
   ierr = MDI_Send_command(">CELL", mdicomm);
   if (ierr) error->all(FLERR, "MDI: >CELL command");
-  cell[0] = domain->boxhi[0] - domain->boxlo[0];
+  ierr = MDI_Send(qm_cell, 9, MDI_DOUBLE, mdicomm);
-  cell[1] = 0.0;
-  cell[2] = 0.0;
-  cell[3] = domain->xy;
-  cell[4] = domain->boxhi[1] - domain->boxlo[1];
-  cell[5] = 0.0;
-  cell[6] = domain->xz;
-  cell[7] = domain->yz;
-  cell[8] = domain->boxhi[2] - domain->boxlo[2];
-
-  // convert cell units to bohr
-  
-  for (int icell = 0; icell < 9; icell++) cell[icell] *= lmp2mdi_length;
-
-  ierr = MDI_Send(cell, 9, MDI_DOUBLE, mdicomm);
   if (ierr) error->all(FLERR, "MDI: >CELL data");
 }
 

src/MDI/fix_mdi_qm.h

@@ -50,11 +50,13 @@ class FixMDIQM : public Fix {
   int nprocs;
   int every, virialflag, addflag, connectflag;
   int plugin;
-  int sumflag,exclude;
+  int sumflag,mcflag;
-  char *id_exclude;
+  char *id_mcfix;
-  int *exclusion_group_ptr;
+  int *mc_active_ptr,*exclusion_group_ptr;
   int *elements;
 
+  double qm_cell[9],qm_cell_displ[3];
+  
   double qm_energy;
   double qm_virial[9], qm_virial_symmetric[6];
 
@@ -92,6 +94,7 @@ class FixMDIQM : public Fix {
   int set_nqm();
   void create_qm_list();
   void set_qm2owned();
+  void set_box();
   void set_xqm();
   void set_eqm();
   void set_tqm();