support for >ELEMENTS MDI command

examples/mdi/Run.sh

@@ -3,7 +3,7 @@
 # -------------------------------------------------
 # -------------------------------------------------
 
-# Example 1
+# Example 1 = run ab initio MD (AIMD)
 
 # ---
 
@@ -56,7 +56,7 @@ mv log.aimd.engine.plugin log.aimd.engine.plugin.3
 # -------------------------------------------------
 # -------------------------------------------------
 
-# Example 2
+# Example 2 = run LAMMPS, compute QM forces on snapshots from a long run
 
 # ---
 
@@ -116,7 +116,7 @@ mv dump.snapshot.driver.plugin dump.snapshot.driver.plugin.3
 # -------------------------------------------------
 # -------------------------------------------------
 
-# Example 3
+# Example 3 = run LAMMPS, compute QM forces on series of independent systems
 
 # ---
 
@@ -176,7 +176,7 @@ mv dump.series.driver.plugin dump.series.driver.plugin.3
 # -------------------------------------------------
 # -------------------------------------------------
 
-# Example 4
+# Example 4 = Python driver runs a sequence of unrelated LAMMPS calculations
 
 # ---
 
@@ -221,7 +221,7 @@ mpirun -np 3 python3 sequence_driver.py -plugin lammps -mdi "-role DRIVER -name
 # -------------------------------------------------
 # -------------------------------------------------
 
-# Example 5
+# Example 5 = run AIMD with Python driver code and 2 LAMMPS instances as engines
 
 # ---
 

src/MDI/fix_mdi_qm.cpp

@@ -25,6 +25,8 @@ using namespace FixConst;
 
 enum { NATIVE, REAL, METAL };    // LAMMPS units which MDI supports
 
+#define MAXELEMENT 103           // used elsewhere in MDI package
+
 /* ---------------------------------------------------------------------- */
 
 FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
@@ -49,6 +51,7 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
   addflag = 1;
   every = 1;
   connectflag = 1;
+  elements = nullptr;
 
   int iarg = 3;
   while (iarg < narg) {
@@ -75,6 +78,17 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
       else if (strcmp(arg[iarg+1],"no") == 0) connectflag = 0;
       else error->all(FLERR,"Illegal fix mdi/qm command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"elements") == 0) {
+      int ntypes = atom->ntypes;
+      if (iarg+ntypes+1 > narg) error->all(FLERR,"Illegal fix mdi/qm command");
+      delete [] elements;
+      elements = new int[ntypes+1];
+      for (int i = 1; i <= ntypes; i++) {
+        elements[i] = utils::inumeric(FLERR,arg[iarg+i],false,lmp);
+        if (elements[i] < 1 || elements[i] > MAXELEMENT)
+          error->all(FLERR,"Illegal fix mdi/qm command");
+      }
+      iarg += ntypes+1;
     } else error->all(FLERR,"Illegal fix mdi/qm command");
   }
 
@@ -154,6 +168,8 @@ FixMDIQM::~FixMDIQM()
 
   // clean up
 
+  delete[] elements;
+
   memory->destroy(fqm);
 
   memory->destroy(ibuf1);
@@ -220,7 +236,7 @@ void FixMDIQM::init()
     }
   }
 
-  // send natoms, atom types, and simulation box to engine
+  // send natoms, atom types or elements, and simulation box to engine
   // this will trigger setup of a new system
   // subsequent calls in post_force() will be for same system until new init()
 
@@ -232,7 +248,8 @@ void FixMDIQM::init()
   ierr = MDI_Send(&n, 1, MDI_INT, mdicomm);
   if (ierr) error->all(FLERR, "MDI: >NATOMS data");
 
-  send_types();
+  if (elements) send_elements();
+  else send_types();
   send_box();
 }
 
@@ -420,7 +437,7 @@ void FixMDIQM::reallocate()
 }
 
 /* ----------------------------------------------------------------------
-   send numeric atom types to MDI engine
+   send LAMMPS atom types to MDI engine
 ------------------------------------------------------------------------- */
 
 void FixMDIQM::send_types()
@@ -448,6 +465,34 @@ void FixMDIQM::send_types()
   if (ierr) error->all(FLERR, "MDI: >TYPES data");
 }
 
+/* ----------------------------------------------------------------------
+   send elements to MDI engine = atomic numbers for each type
+------------------------------------------------------------------------- */
+
+void FixMDIQM::send_elements()
+{
+  int n = static_cast<int> (atom->natoms);
+  memset(ibuf1, 0, n * sizeof(int));
+
+  // use local atomID to index into ordered ibuf1
+
+  tagint *tag = atom->tag; 
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  
+  int index;
+  for (int i = 0; i < nlocal; i++) {
+    index = static_cast<int>(tag[i]) - 1;
+    ibuf1[index] = elements[type[i]];
+  }
+
+  MPI_Reduce(ibuf1, ibuf1all, n, MPI_INT, MPI_SUM, 0, world);
+
+  int ierr = MDI_Send_command(">ELEMENTS", mdicomm);
+  if (ierr) error->all(FLERR, "MDI: >ELEMENTS command");
+  ierr = MDI_Send(ibuf1all, n, MDI_INT, mdicomm);
+  if (ierr) error->all(FLERR, "MDI: >ELEMETNS data");
+}
 
 /* ----------------------------------------------------------------------
    send simulation box size and shape to MDI engine

src/MDI/fix_mdi_qm.h

@@ -44,6 +44,7 @@ class FixMDIQM : public Fix {
   int plugin;
   int maxlocal;
   int sumflag;
+  int *elements;
 
   double qm_energy;
   int lmpunits;
@@ -70,6 +71,7 @@ class FixMDIQM : public Fix {
 
   void reallocate();
   void send_types();
+  void send_elements();
   void send_box();
   void unit_conversions();
 };

src/MDI/mdi_engine.cpp

@@ -54,6 +54,8 @@ enum { DEFAULT, MD, OPT };       // top-level MDI engine modes
 
 enum { TYPE, CHARGE, MASS, COORD, VELOCITY, FORCE, ADDFORCE };
 
+#define MAXELEMENT 103           // used elsewhere in MDI package
+
 /* ----------------------------------------------------------------------
    trigger LAMMPS to start acting as an MDI engine
    either in standalone mode or plugin mode
@@ -63,17 +65,47 @@ enum { TYPE, CHARGE, MASS, COORD, VELOCITY, FORCE, ADDFORCE };
    when EXIT command is received, mdi engine command exits
 ---------------------------------------------------------------------- */
 
-MDIEngine::MDIEngine(LAMMPS *_lmp, int narg, char ** /*arg*/) : Pointers(_lmp)
+MDIEngine::MDIEngine(LAMMPS *_lmp, int narg, char ** arg) : Pointers(_lmp)
 {
-  if (narg) error->all(FLERR, "Illegal mdi engine command");
-
   // check requirements for LAMMPS to work with MDI as an engine
 
-  if (atom->tag_enable == 0) error->all(FLERR, "Cannot use MDI engine without atom IDs");
+  if (atom->tag_enable == 0) error->all(FLERR, "MDI engine requires atom IDs");
 
   if (atom->natoms && atom->tag_consecutive() == 0)
     error->all(FLERR, "MDI engine requires consecutive atom IDs");
 
+  // optional args
+
+  elements = nullptr;
+
+  int iarg = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"elements") == 0) {
+      int ntypes = atom->ntypes;
+      delete [] elements;
+      elements = new int[ntypes+1];
+      if (iarg+ntypes+1 > narg) error->all(FLERR,"Illegal mdi engine command");
+      for (int i = 1; i <= ntypes; i++) {
+        elements[i] = utils::inumeric(FLERR,arg[iarg+i],false,lmp);
+        if (elements[i] < 0 || elements[i] > MAXELEMENT)
+          error->all(FLERR,"Illegal mdi engine command");
+      }
+      iarg += ntypes+1;
+    } else error->all(FLERR,"Illegal mdi engine command");
+  }
+
+  // error check an MDI element does not map to multiple atom types
+
+  if (elements) {
+    int ntypes = atom->ntypes;
+    for (int i = 1; i < ntypes; i++)
+      for (int j = i+1; j <= ntypes; j++) {
+        if (elements[i] == 0 || elements[j] == 0) continue;
+        if (elements[i] == elements[j]) 
+          error->all(FLERR,"MDI engine element cannot map to multiple types");
+      }
+  }
+
   // confirm LAMMPS is being run as an engine
 
   int role;
@@ -194,6 +226,8 @@ MDIEngine::MDIEngine(LAMMPS *_lmp, int narg, char ** /*arg*/) : Pointers(_lmp)
 
   // clean up
 
+  delete[] elements;
+
   delete[] mdicmd;
   delete[] node_engine;
   delete[] node_driver;
@@ -299,6 +333,11 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm)
   } else if (strcmp(command, ">COORDS") == 0) {
     receive_coords();
 
+  } else if (strcmp(command, ">ELEMENTS") == 0) {
+    if (!elements)
+      error->all(FLERR,"MDI engine command did not define element list");
+    receive_elements();
+
   } else if (strcmp(command, ">FORCES") == 0) {
     receive_double3(FORCE);
 
@@ -479,6 +518,7 @@ void MDIEngine::mdi_commands()
   MDI_Register_command("@DEFAULT", ">CELL_DISPL");
   MDI_Register_command("@DEFAULT", ">CHARGES");
   MDI_Register_command("@DEFAULT", ">COORDS");
+  MDI_Register_command("@DEFAULT", ">ELEMENTS");
   MDI_Register_command("@DEFAULT", ">NATOMS");
   MDI_Register_command("@DEFAULT", ">NSTEPS");
   MDI_Register_command("@DEFAULT", ">TOLERANCE");
@@ -914,7 +954,7 @@ void MDIEngine::evaluate()
 
 /* ----------------------------------------------------------------------
    create a new system
-   >CELL, >NATOMS, >TYPES, >COORDS commands are required
+   >CELL, >NATOMS, >TYPES or >ELEMENTS, >COORDS commands are required
    >CELL_DISPL, >CHARGES, >VELOCITIES commands are optional
 ---------------------------------------------------------------------- */
 
@@ -924,8 +964,8 @@ void MDIEngine::create_system()
 
   if (flag_cell == 0 || flag_natoms == 0 || flag_types == 0 || flag_coords == 0)
     error->all(FLERR,
-               "MDI create_system requires >CELL, >NATOMS, >TYPES, >COORDS "
+               "MDI create_system requires >CELL, >NATOMS, "
-               "MDI commands");
+               ">TYPES or >ELEMENTS, >COORDS MDI commands");
 
   // remove all existing atoms via delete_atoms command
 
@@ -1160,6 +1200,38 @@ void MDIEngine::receive_coords()
   for (int i = 0; i < n; i++) sys_coords[i] *= mdi2lmp_length;
 }
 
+/* ----------------------------------------------------------------------
+   >ELEMENTS command
+   receive elements for each atom = atomic numbers
+   convert to LAMMPS atom types and store in sys_types
+---------------------------------------------------------------------- */
+
+void MDIEngine::receive_elements()
+{
+  actionflag = 0;
+  flag_types = 1;
+  int ierr = MDI_Recv(sys_types, sys_natoms, MDI_INT, mdicomm);
+  if (ierr) error->all(FLERR, "MDI: >ELEMENTS data");
+  MPI_Bcast(sys_types, sys_natoms, MPI_INT, 0, world);
+
+  // convert from element atomic numbers to LAMMPS atom types
+  // use maping provided by mdi engine command
+
+  int ntypes = atom->ntypes;
+  int itype;
+
+  for (int i = 0; i < sys_natoms; i++) {
+    for (itype = 1; itype <= ntypes; itype++) {
+      if (sys_types[i] == elements[itype]) {
+        sys_types[i] = itype;
+        break;
+      }
+    }
+    if (itype > ntypes) 
+      error->all(FLERR,"MDI element not found in element list");
+  }
+}
+
 /* ----------------------------------------------------------------------
    >NATOMS command
    natoms cannot exceed 32-bit int for use with MDI

src/MDI/mdi_engine.h

@@ -70,6 +70,8 @@ class MDIEngine : protected Pointers {
 
   int actionflag;    // 1 if MD or OPTG just completed, else 0
 
+  int *elements;
+
   // buffers for MDI comm
 
   int maxatom;
@@ -106,6 +108,7 @@ class MDIEngine : protected Pointers {
   void receive_cell_displ();
   void receive_charges();
   void receive_coords();
+  void receive_elements();
   void receive_natoms();
   void receive_nsteps();
   void receive_tolerance();