diff --git a/examples/mdi/Run.sh b/examples/mdi/Run.sh index 3817841103..3b70a8bf79 100644 --- a/examples/mdi/Run.sh +++ b/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 # --- diff --git a/src/MDI/fix_mdi_qm.cpp b/src/MDI/fix_mdi_qm.cpp index 38b7270651..1c725a0871 100644 --- a/src/MDI/fix_mdi_qm.cpp +++ b/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 (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(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 diff --git a/src/MDI/fix_mdi_qm.h b/src/MDI/fix_mdi_qm.h index 1db6a0cce2..11180925b1 100644 --- a/src/MDI/fix_mdi_qm.h +++ b/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(); }; diff --git a/src/MDI/mdi_engine.cpp b/src/MDI/mdi_engine.cpp index 3683e3f3c5..79bfbe097e 100644 --- a/src/MDI/mdi_engine.cpp +++ b/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); @@ -323,7 +362,7 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) else receive_double3(VELOCITY); - // ----------------------------------------------- + // ----------------------------------------------- } else if (strcmp(command, "<@") == 0) { ierr = MDI_Send(node_engine, MDI_NAME_LENGTH, MDI_CHAR, mdicomm); @@ -372,9 +411,9 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) } else if (strcmp(command, "all(FLERR, "MDI: MDI engine is already performing a simulation"); @@ -419,14 +458,14 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) strncpy(node_driver, command, MDI_COMMAND_LENGTH); node_match = false; - // exit command + // exit command } else if (strcmp(command, "EXIT") == 0) { exit_command = true; - // ------------------------------------------------------- - // custom LAMMPS commands - // ------------------------------------------------------- + // ------------------------------------------------------- + // custom LAMMPS commands + // ------------------------------------------------------- } else if (strcmp(command, "NBYTES") == 0) { nbytes_command(); @@ -439,9 +478,9 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) } else if (strcmp(command, "all(FLERR, "MDI: Unknown command {} received from driver", command); @@ -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 commands"); + "MDI create_system requires >CELL, >NATOMS, " + ">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 diff --git a/src/MDI/mdi_engine.h b/src/MDI/mdi_engine.h index ea840fba02..242755e3b7 100644 --- a/src/MDI/mdi_engine.h +++ b/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();