From a8a97962d2e207d4a5d1faa8a018fa320d2670fe Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 3 Sep 2021 15:50:17 -0600 Subject: [PATCH] more changes for AIMD testing --- src/MDI/fix_mdi_aimd.cpp | 233 +++++++++++++++++++-- src/MDI/fix_mdi_aimd.h | 25 ++- src/MDI/mdi_engine2.cpp | 436 ++++++++++++++++++++++++--------------- src/MDI/mdi_engine2.h | 6 +- 4 files changed, 507 insertions(+), 193 deletions(-) diff --git a/src/MDI/fix_mdi_aimd.cpp b/src/MDI/fix_mdi_aimd.cpp index f8ebcaac85..c2acc44c94 100644 --- a/src/MDI/fix_mdi_aimd.cpp +++ b/src/MDI/fix_mdi_aimd.cpp @@ -17,10 +17,14 @@ #include "comm.h" #include "domain.h" #include "force.h" +#include "memory.h" +#include "update.h" using namespace LAMMPS_NS; using namespace FixConst; +enum{NATIVE,REAL,METAL}; // LAMMPS units which MDI supports + /* ---------------------------------------------------------------------- */ FixMDIAimd::FixMDIAimd(LAMMPS *lmp, int narg, char **arg) : @@ -35,9 +39,40 @@ FixMDIAimd::FixMDIAimd(LAMMPS *lmp, int narg, char **arg) : virial_global_flag = 1; thermo_energy = thermo_virial = 1; + // check requirements for LAMMPS to work with MDI as an engine + + if (atom->tag_enable == 0) + error->all(FLERR, "Cannot use mdi engined without atom IDs"); + + if (atom->natoms && atom->tag_consecutive() == 0) + error->all(FLERR, "mdi engine requires consecutive atom IDs"); + + // confirm LAMMPS is being run as a driver + + int role; + MDI_Get_role(&role); + if (role != MDI_DRIVER) + error->all(FLERR,"Must invoke LAMMPS as an MDI driver to use fix mdi/aimd"); + + // storage for all atoms + + buf3 = buf3all = nullptr; + maxbuf = 0; + + // set unit conversion factors + + if (strcmp(update->unit_style, "real") == 0) lmpunits = REAL; + else if (strcmp(update->unit_style, "metal") == 0) lmpunits = METAL; + else lmpunits = NATIVE; + + unit_conversions(); + // connect to MDI engine - MDI_Accept_communicator(&engine); + MDI_Accept_communicator(&mdicomm); + if (mdicomm <= 0) error->all(FLERR, "Unable to connect to MDI engine"); + + nprocs = comm->nprocs; } /* ---------------------------------------------------------------------- */ @@ -46,7 +81,13 @@ FixMDIAimd::~FixMDIAimd() { // send exit command to engine - MDI_Send_command("EXIT",engine); + int ierr = MDI_Send_command("EXIT",mdicomm); + if (ierr) error->all(FLERR,"MDI: EXIT command"); + + // clean up + + memory->destroy(buf3); + memory->destroy(buf3all); } /* ---------------------------------------------------------------------- */ @@ -87,43 +128,117 @@ void FixMDIAimd::pre_reverse(int eflag, int /*vflag*/) void FixMDIAimd::post_force(int vflag) { + int ilocal,ierr; + double cell[9]; + int eflag = eflag_caller; ev_init(eflag,vflag); - // send current coords to MDI engine - // NOTE: need to gather them to proc 0 + // if simulation box dynamically changes, send current box to MDI engine - MDI_Send_command(">COORDS",engine); - MDI_Send(&atom->x[0][0],3*atom->nlocal,MDI_DOUBLE,engine); + if (domain->box_change_size || domain->box_change_shape) { + ierr = MDI_Send_command(">CELL_DISPL",mdicomm); + if (ierr) error->all(FLERR,"MDI: >CELL_DISPL command"); + cell[0] = domain->boxlo[0] * lmp2mdi_length; + 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]; + 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]; + ierr = MDI_Send(cell,9,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: >CELL data"); + } + + // gather all coords, ordered by atomID + + reallocate(); + memset(buf3,0,3*atom->natoms*sizeof(double)); + + double **x = atom->x; + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + ilocal = static_cast (tag[i]) - 1; + buf3[3*ilocal+0] = x[i][0] * lmp2mdi_length; + buf3[3*ilocal+1] = x[i][1] * lmp2mdi_length; + buf3[3*ilocal+2] = x[i][2] * lmp2mdi_length; + } + + MPI_Reduce(buf3,buf3all,3*atom->natoms,MPI_DOUBLE,MPI_SUM,0,world); + + // send current coords to MDI engine + + ierr = MDI_Send_command(">COORDS",mdicomm); + if (ierr) error->all(FLERR,"MDI: >COORDS command"); + ierr = MDI_Send(buf3all,3*atom->natoms,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: >COORDS data"); // trigger engine to evaluate forces,energy,pressure for current system - MDI_Send_command("EVAL",engine); + ierr = MDI_Send_command("EVAL",mdicomm); + if (ierr) error->all(FLERR,"MDI: EVAL command"); - // request energy from MDI engine + // request forces from MDI engine - if (eflag_global) { - MDI_Send_command("all(FLERR,"MDI: natoms,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: natoms,MPI_DOUBLE,0,world); + + // add forces to owned atoms + // use atomID to index into ordered buf3 + + double **f = atom->f; + + for (int i = 0; i < nlocal; i++) { + ilocal = static_cast (tag[i]) - 1; + f[i][0] += buf3[3*ilocal+0] * mdi2lmp_force; + f[i][1] += buf3[3*ilocal+1] * mdi2lmp_force; + f[i][2] += buf3[3*ilocal+2] * mdi2lmp_force; } - // request pressure tensor from MDI engine, convert to virial + // optionally request energy from MDI engine + // divide by nprocs so each proc stores a portion + + if (eflag_global) { + ierr = MDI_Send_command("all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: xprd * domain->yprd * domain->zprd; - for (int i = 0; i < 6; i++) - virial[i] = ptensor[i] * volume / force->nktv2p; + for (int i = 0; i < 6; i++) { + ptensor[i] *= mdi2lmp_pressure; + virial[i] = ptensor[i] * volume / force->nktv2p / nprocs; + } } - - // request forces from MDI engine - // NOTE: need to scatter to procs, then add to my local forces - - MDI_Send_command("f[0][0],3*atom->nlocal,MDI_DOUBLE,engine); } /* ---------------------------------------------------------------------- */ @@ -141,3 +256,77 @@ double FixMDIAimd::compute_scalar() { return engine_energy; } + +/* ---------------------------------------------------------------------- + reallocate storage for all atoms if necessary +------------------------------------------------------------------------- */ + +void FixMDIAimd::reallocate() +{ + if (atom->natoms <= maxbuf) return; + + if (3*atom->natoms > MAXSMALLINT) + error->all(FLERR,"Natoms too large to use with fix mdi/aimd"); + + maxbuf = atom->natoms; + + memory->destroy(buf3); + memory->destroy(buf3all); + + memory->create(buf3,3*maxbuf,"mdi:buf3"); + memory->create(buf3all,3*maxbuf,"mdi:buf3all"); +} + +/* ---------------------------------------------------------------------- + MDI to/from LAMMPS conversion factors +------------------------------------------------------------------------- */ + +void FixMDIAimd::unit_conversions() +{ + double angstrom_to_bohr,kelvin_to_hartree,ev_to_hartree; + + MDI_Conversion_factor("angstrom","bohr",&angstrom_to_bohr); + MDI_Conversion_factor("kelvin_energy","hartree",&kelvin_to_hartree); + MDI_Conversion_factor("electron_volt","hartree",&ev_to_hartree); + + // length units + + mdi2lmp_length = 1.0; + lmp2mdi_length = 1.0; + + if (lmpunits == REAL || lmpunits == METAL) { + lmp2mdi_length = angstrom_to_bohr; + mdi2lmp_length = 1.0 / angstrom_to_bohr; + } + + // energy units + + mdi2lmp_energy = 1.0; + lmp2mdi_energy = 1.0; + + if (lmpunits == REAL) { + lmp2mdi_energy = kelvin_to_hartree / force->boltz; + mdi2lmp_energy = force->boltz / kelvin_to_hartree; + } else if (lmpunits == METAL) { + lmp2mdi_energy = ev_to_hartree; + mdi2lmp_energy = 1.0 / ev_to_hartree; + } + + // force units + + mdi2lmp_force = 1.0; + lmp2mdi_force = 1.0; + + if (lmpunits == REAL) { + lmp2mdi_force = (kelvin_to_hartree / force->boltz) / angstrom_to_bohr; + mdi2lmp_force = 1.0 / lmp2mdi_force; + } else if (lmpunits == METAL) { + lmp2mdi_force = ev_to_hartree / angstrom_to_bohr; + mdi2lmp_force = angstrom_to_bohr / ev_to_hartree; + } + + // pressure units + + mdi2lmp_pressure = 1.0; + lmp2mdi_pressure = 1.0; +} diff --git a/src/MDI/fix_mdi_aimd.h b/src/MDI/fix_mdi_aimd.h index 8ea62e1f1e..63db7266cf 100644 --- a/src/MDI/fix_mdi_aimd.h +++ b/src/MDI/fix_mdi_aimd.h @@ -39,10 +39,33 @@ class FixMDIAimd : public Fix { double compute_scalar(); private: + int nprocs; int eflag_caller; double engine_energy; + int lmpunits; + + // MDI communicator + + MDI_Comm mdicomm; + + // unit conversion factors + + double lmp2mdi_length,mdi2lmp_length; + double lmp2mdi_energy,mdi2lmp_energy; + double lmp2mdi_force,mdi2lmp_force; + double lmp2mdi_pressure,mdi2lmp_pressure; + double lmp2mdi_virial,mdi2lmp_virial; + + // buffers for MDI comm + + int maxbuf; + double *buf3,*buf3all; + + // methods + + void reallocate(); + void unit_conversions(); - MDI_Comm engine; }; } diff --git a/src/MDI/mdi_engine2.cpp b/src/MDI/mdi_engine2.cpp index df5aba7fa3..e55a7b44a0 100644 --- a/src/MDI/mdi_engine2.cpp +++ b/src/MDI/mdi_engine2.cpp @@ -29,7 +29,7 @@ #include "force.h" #include "group.h" #include "input.h" -#include "irregular.h" +//#include "irregular.h" #include "library.h" #include "memory.h" #include "min.h" @@ -45,20 +45,24 @@ using namespace LAMMPS_NS; -enum{NATIVE,REAL,METAL}; // LAMMPS units which MDI supports -enum{DEFAULT,MD,OPT}; // top-level MDI engine mode -enum{TYPE,CHARGE,MASS,COORD,VELOCITY,FORCE}; +enum{NATIVE,REAL,METAL}; // LAMMPS units which MDI supports +enum{DEFAULT,MD,OPT}; // top-level MDI engine mode + +// per-atom data which engine commands access + +enum{TYPE,CHARGE,MASS,COORD,VELOCITY,FORCE}; /* ---------------------------------------------------------------------- - mdi command: quit or engine + mdi command: engine + NOTE: may later have other MDI command variants? ---------------------------------------------------------------------- */ void MDIEngine2::command(int narg, char **arg) { - if (narg < 1) error->all(FLERR, "Illegal mdi command"); + if (narg < 1) error->all(FLERR,"Illegal mdi command"); if (strcmp(arg[0],"engine") == 0) mdi_engine(narg-1,&arg[1]); - else error->all(FLERR, "Illegal mdi command"); + else error->all(FLERR,"Illegal mdi command"); } /* ---------------------------------------------------------------------- @@ -72,34 +76,27 @@ void MDIEngine2::command(int narg, char **arg) void MDIEngine2::mdi_engine(int narg, char **arg) { - // args to choose what MDI nodes LAMMPS will interact with - // also whether a fix mdi/engine is needed - // NOTE: add args - - if (narg > 0) error->all(FLERR, "Illegal mdi engine command"); + if (narg > 0) 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 engined without atom IDs"); + error->all(FLERR,"Cannot use MDI engine without atom IDs"); if (atom->natoms && atom->tag_consecutive() == 0) - error->all(FLERR, "mdi engine requires consecutive atom IDs"); - - if (strcmp(update->unit_style, "real") == 0) lmpunits = REAL; - else if (strcmp(update->unit_style, "metal") == 0) lmpunits = METAL; - else lmpunits = NATIVE; + error->all(FLERR,"MDI engine requires consecutive atom IDs"); // confirm LAMMPS is being run as an engine int role; MDI_Get_role(&role); if (role != MDI_ENGINE) - error->all(FLERR, - "Must invoke LAMMPS as an MDI engine to use mdi/engine command"); + error->all(FLERR,"Must invoke LAMMPS as an MDI engine to use mdi engine"); - // if the mdi/engine fix is not already present, add it now - // NOTE: make this optional above + // NOTE: create this fix only when @INIT_MD is triggered? + // if not needed, I dont think it should be defined, + // otherwise it will be be invoked in other use cases + // and switch engine out of DEFAULT node ? //int ifix = modify->find_fix_by_style("mdi/engine"); //bool added_mdi_engine_fix = false; @@ -141,12 +138,6 @@ void MDIEngine2::mdi_engine(int narg, char **arg) id_press = utils::strdup(std::string("MDI_ENGINE") + "_press"); modify->add_compute(fmt::format("{} all pressure thermo_temp", id_press)); - // set unit conversion factors - - unit_conversions(); - - // store pointers to the new computes - int icompute_ke = modify->find_compute(id_ke); int icompute_pe = modify->find_compute(id_pe); int icompute_press = modify->find_compute(id_press); @@ -155,13 +146,22 @@ void MDIEngine2::mdi_engine(int narg, char **arg) pe = modify->compute[icompute_pe]; press = modify->compute[icompute_press]; + // set unit conversion factors + + if (strcmp(update->unit_style, "real") == 0) lmpunits = REAL; + else if (strcmp(update->unit_style, "metal") == 0) lmpunits = METAL; + else lmpunits = NATIVE; + + unit_conversions(); + // irregular class and data structs used by MDI + // NOTE: not clear if irregular comm is ever needed - irregular = new Irregular(lmp); + //irregular = new Irregular(lmp); - buf1 = nullptr; - buf3 = nullptr; - ibuf1 = nullptr; + buf1 = buf1all = nullptr; + buf3 = buf3all = nullptr; + ibuf1 = ibuf1all = nullptr; maxbuf = 0; // define MDI commands that LAMMPS engine recognizes @@ -171,7 +171,7 @@ void MDIEngine2::mdi_engine(int narg, char **arg) // one-time operation to establish a connection with the driver MDI_Accept_communicator(&mdicomm); - if (mdicomm <= 0) error->all(FLERR, "Unable to connect to MDI driver"); + if (mdicomm <= 0) error->all(FLERR,"Unable to connect to MDI driver"); // endless engine loop, responding to driver commands @@ -196,8 +196,8 @@ void MDIEngine2::mdi_engine(int narg, char **arg) break; } else - error->all(FLERR,fmt::format("MDI node exited with invalid command: {}" - ,mdicmd)); + error->all(FLERR, + fmt::format("MDI node exited with invalid command: {}",mdicmd)); } // clean up @@ -206,22 +206,35 @@ void MDIEngine2::mdi_engine(int narg, char **arg) delete [] node_engine; delete [] node_driver; - modify->delete_compute(id_pe); modify->delete_compute(id_ke); + modify->delete_compute(id_pe); modify->delete_compute(id_press); - delete irregular; - memory->destroy(buf3); - memory->destroy(buf1); + delete [] id_ke; + delete [] id_pe; + delete [] id_press; + + //delete irregular; + memory->destroy(ibuf1); + memory->destroy(buf1); + memory->destroy(buf3); - // remove mdi/engine fix that mdi/engine instantiated - // NOTE: make this optional + memory->destroy(ibuf1all); + memory->destroy(buf1all); + memory->destroy(buf3all); + + // remove mdi/engine fix that mdi engine instantiated + // NOTE: decide whether to make this optional, see above //if (added_mdi_engine_fix) modify->delete_fix("MDI_ENGINE_INTERNAL"); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + engine is now at this MDI node + loop over received commands so long as driver also at this node + return when not the case or EXIT command received +---------------------------------------------------------------------- */ void MDIEngine2::engine_node(const char *node) { @@ -251,7 +264,7 @@ void MDIEngine2::engine_node(const char *node) execute_command(mdicmd,mdicomm); - // check if driver node is now something other than engine node + // check if driver node is now somewhere other than engine node if (strcmp(node_driver,"\0") != 0 && strcmp(node_driver,node_engine) != 0) node_match = false; @@ -274,6 +287,7 @@ int MDIEngine2::execute_command(const char *command, MDI_Comm mdicomm) int ierr; // confirm this command is supported at this node + // otherwise is an error int command_exists = 1; if (root) { @@ -286,38 +300,44 @@ int MDIEngine2::execute_command(const char *command, MDI_Comm mdicomm) MPI_Bcast(&command_exists,1,MPI_INT,0,world); if (!command_exists) - error->all(FLERR, "MDI: Received a command unsupported by engine node"); + error->all(FLERR,"MDI: Received a command unsupported by engine node"); + // --------------------------------------- // respond to each possible driver command + // --------------------------------------- - if (strcmp(command, ">NATOMS") == 0) { - // NOTE: needs to be 64-bit value or copied into 64-bit value - ierr = MDI_Recv((char *) &atom->natoms,1,MDI_INT,mdicomm); - if (ierr) - error->all(FLERR, "MDI: Unable to receive number of atoms from driver"); - MPI_Bcast(&atom->natoms, 1, MPI_INT, 0, world); + if (strcmp(command,">NATOMS") == 0) { - } else if (strcmp(command, "natoms; - ierr = MDI_Send((char *) &mdi_natoms, 1, MDI_INT64_T, mdicomm); - if (ierr != 0) - error->all(FLERR, "MDI: Unable to send number of atoms to driver"); + // natoms cannot exceed 32-bit int for use with MDI + + int natoms; + ierr = MDI_Recv(&natoms,1,MDI_INT,mdicomm); + if (ierr) error->all(FLERR,"MDI: >NATOMS data"); + MPI_Bcast(&natoms,1,MPI_INT,0,world); + if (natoms < 0) error->all(FLERR,"MDI received natoms < 0"); + atom->natoms = natoms; + + } else if (strcmp(command," (atom->natoms); + ierr = MDI_Send(&natoms,1,MDI_INT,mdicomm); + if (ierr != 0) error->all(FLERR,"MDI: ntypes, 1, MDI_INT, mdicomm); - if (ierr != 0) - error->all(FLERR, "MDI: Unable to send number of atom types to driver"); + ierr = MDI_Send(&atom->ntypes,1,MDI_INT,mdicomm); + if (ierr != 0) error->all(FLERR, "MDI: +FORCES") == 0) { receive_double3(FORCE,1); - } else if (strcmp(command, "all(FLERR, "MDI: MDI engine is already performing a simulation"); + error->all(FLERR,"MDI: MDI engine is already performing a simulation"); mode = MD; node_match = false; } else if (strcmp(command,"@INIT_OPTG") == 0) { if (mode != DEFAULT) - error->all(FLERR, "MDI: MDI engine is already performing a simulation"); + error->all(FLERR,"MDI: MDI engine is already performing a simulation"); mode = OPT; node_match = false; @@ -377,7 +397,7 @@ int MDIEngine2::execute_command(const char *command, MDI_Comm mdicomm) } else if (strcmp(command,"<@") == 0) { ierr = MDI_Send(node_engine,MDI_NAME_LENGTH,MDI_CHAR,mdicomm); - if (ierr) error->all(FLERR,"MDI: Unable to send node to driver"); + if (ierr) error->all(FLERR,"MDI: <@ data"); } else if (strcmp(command,"@DEFAULT") == 0) { mode = DEFAULT; @@ -446,7 +466,7 @@ int MDIEngine2::execute_command(const char *command, MDI_Comm mdicomm) // ------------------------------------------------------- } else { - error->all(FLERR, "MDI: Unknown command from driver"); + error->all(FLERR,"MDI: Unknown command received from driver"); } return 0; @@ -461,6 +481,7 @@ void MDIEngine2::mdi_commands() { // ------------------------------------ // commands and nodes that an MDI-compliant MD code supports + // NOTE: is all of this correct? // ------------------------------------ // default node and its commands @@ -593,7 +614,7 @@ void MDIEngine2::mdi_commands() MDI_Register_command("@COORDS", "EXIT"); // ------------------------------------ - // custom commands and nodes which LAMMPS adds support for + // custom commands and nodes which LAMMPS supports // max length for a command is current 11 chars in MDI // ------------------------------------ @@ -617,7 +638,7 @@ void MDIEngine2::mdi_commands() void MDIEngine2::mdi_md() { - // initialize an MD simulation + // initialize a new MD simulation update->whichflag = 1; timer->init_timeout(); @@ -638,6 +659,7 @@ void MDIEngine2::mdi_md() // setup the MD simulation update->integrate->setup(1); + engine_node("@FORCES"); if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) return; @@ -720,6 +742,7 @@ void MDIEngine2::mdi_optg() // ---------------------------------------------------------------------- /* ---------------------------------------------------------------------- + >CHARGES command receive vector of 1 double for all atoms atoms are ordered by atomID, 1 to Natoms assumes all atoms already exist @@ -729,8 +752,8 @@ void MDIEngine2::receive_double1(int which) { reallocate(); - int ierr = MDI_Recv((char *) buf1,atom->natoms,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: Unable to receive double1 data from driver"); + int ierr = MDI_Recv(buf1,atom->natoms,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: >double1 data"); MPI_Bcast(buf1,atom->natoms,MPI_DOUBLE,0,world); // extract onwed atom value @@ -751,6 +774,7 @@ void MDIEngine2::receive_double1(int which) } /* ---------------------------------------------------------------------- + >TYPES command receive vector of 1 int for all atoms atoms are ordered by atomID, 1 to Natoms assumes all atoms already exist @@ -760,8 +784,8 @@ void MDIEngine2::receive_int1(int which) { reallocate(); - int ierr = MDI_Recv((char *) ibuf1,atom->natoms,MDI_INT,mdicomm); - if (ierr) error->all(FLERR,"MDI: Unable to receive double1 data from driver"); + int ierr = MDI_Recv(ibuf1,atom->natoms,MDI_INT,mdicomm); + if (ierr) error->all(FLERR,"MDI: >int1 data"); MPI_Bcast(ibuf1,atom->natoms,MPI_INT,0,world); // extract onwed atom value @@ -782,6 +806,7 @@ void MDIEngine2::receive_int1(int which) } /* ---------------------------------------------------------------------- + >COORDS, >FORCES commands receive vector of 3 doubles for all atoms atoms are ordered by atomID, 1 to Natoms assumes all atoms already exist @@ -792,8 +817,8 @@ void MDIEngine2::receive_double3(int which, int addflag) { reallocate(); - int ierr = MDI_Recv((char *) buf3,3*atom->natoms,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: Unable to receive double3 data from driver"); + int ierr = MDI_Recv(buf3,3*atom->natoms,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: >double3 data"); MPI_Bcast(buf3,3*atom->natoms,MPI_DOUBLE,0,world); // extract owned atom values @@ -830,16 +855,11 @@ void MDIEngine2::receive_double3(int which, int addflag) f[i][2] += buf3[3*ilocal+2] * mdi2lmp_force; } } - } else if (which == VELOCITY) { - double **v = atom->v; - for (int i = 0; i < nlocal; i++) { - ilocal = static_cast (tag[i]) - 1; - v[i][0] = buf3[3*ilocal+0] * mdi2lmp_velocity; - v[i][1] = buf3[3*ilocal+1] * mdi2lmp_velocity; - v[i][2] = buf3[3*ilocal+2] * mdi2lmp_velocity; - } } + // NOTE: these operations cannot be done in the middle + // of an arbitrary timestep, only when reneighboring is done + /* // ensure atoms are in current box & update box via shrink-wrap // has to be be done before invoking Irregular::migrate_atoms() @@ -860,6 +880,7 @@ void MDIEngine2::receive_double3(int which, int addflag) } /* ---------------------------------------------------------------------- + natoms,MPI_DOUBLE,MPI_SUM,0,world); - int ierr = MDI_Send((char *) buf1all,atom->natoms,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR, "MDI: Unable to send double1 data to driver"); + int ierr = MDI_Send(buf1all,atom->natoms,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: natoms,MPI_INT,MPI_SUM,0,world); - int ierr = MDI_Send((char *) ibuf1all,atom->natoms,MDI_INT,mdicomm); - if (ierr) error->all(FLERR, "MDI: Unable to send int1 data to driver"); + int ierr = MDI_Send(ibuf1all,atom->natoms,MDI_INT,mdicomm); + if (ierr) error->all(FLERR,"MDI: f; - for (int i = 0; i < nlocal; i++) { - ilocal = static_cast (tag[i]) - 1; - buf3[3*ilocal+0] = f[i][0] * lmp2mdi_velocity; - buf3[3*ilocal+1] = f[i][1] * lmp2mdi_velocity; - buf3[3*ilocal+2] = f[i][2] * lmp2mdi_velocity; - } } MPI_Reduce(buf3,buf3all,3*atom->natoms,MPI_DOUBLE,MPI_SUM,0,world); - int ierr = MDI_Send((char *) buf3all,3*atom->natoms,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR, "MDI: Unable to send double3 data to driver"); + int ierr = MDI_Send(buf3all,3*atom->natoms,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: natoms * MDI_LABEL_LENGTH]; - memset(labels, ' ', atom->natoms * MDI_LABEL_LENGTH); + memset(labels,' ',atom->natoms * MDI_LABEL_LENGTH); + + // NOTE: this loop will not work in parallel for (int iatom = 0; iatom < atom->natoms; iatom++) { std::string label = std::to_string(atom->type[iatom]); @@ -998,15 +1019,18 @@ void MDIEngine2::send_labels() strncpy(&labels[iatom * MDI_LABEL_LENGTH], label.c_str(), label_len); } - int ierr = MDI_Send(labels, atom->natoms * MDI_LABEL_LENGTH, MDI_CHAR, mdicomm); - if (ierr) error->all(FLERR, "MDI: Unable to send atom labels to driver"); + int ierr = MDI_Send(labels,atom->natoms*MDI_LABEL_LENGTH,MDI_CHAR,mdicomm); + if (ierr) error->all(FLERR,"MDI: all(FLERR, "MDI: Unable to send total energy to driver"); + int ierr = MDI_Send(&total_energy,1,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: compute_scalar(); potential_energy *= lmp2mdi_energy; - int ierr = MDI_Send((char *) &potential_energy,1,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR, "MDI: Unable to send potential energy to driver"); + int ierr = MDI_Send(&potential_energy,1,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: compute_scalar(); kinetic_energy *= lmp2mdi_energy; - int ierr = MDI_Send((char *) &kinetic_energy,1,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: Unable to send kinetic energy to driver"); + int ierr = MDI_Send(&kinetic_energy,1,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: all(FLERR, "MDI: Unable to send cell dimensions to driver"); + int ierr = MDI_Send(celldata,9,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: CELL command + reset the simulation box edge vectors + in conjunction with >CELL_DISPL this can adjust box arbitrarily +---------------------------------------------------------------------- */ void MDIEngine2::receive_cell() { double celldata[9]; - int ierr = MDI_Recv((char *) celldata,9,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR, "MDI: Unable to send cell dimensions to driver"); - MPI_Bcast(&celldata[0],9,MPI_DOUBLE,0,world); + int ierr = MDI_Recv(celldata,9,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR, "MDI: >CELL data"); + MPI_Bcast(celldata,9,MPI_DOUBLE,0,world); for (int icell = 0; icell < 9; icell++) - celldata[icell] /= mdi2lmp_length; + celldata[icell] *= mdi2lmp_length; - // require the new cell vector be orthogonal (for now) + // error check that edge vectors match LAMMPS triclinic requirement - double small = std::numeric_limits::min(); - if (abs(celldata[1]) > small or abs(celldata[2]) > small or - abs(celldata[3]) > small or - abs(celldata[5]) > small or abs(celldata[6]) > small or - abs(celldata[7]) > small) { - error->all(FLERR, - "MDI: LAMMPS currently only supports the >CELL command " - "for orthogonal cell vectors"); - } + if (celldata[1] != 0.0 || celldata[2] != 0.0 || celldata[5] != 0.0) + error->all(FLERR,"MDI: Received cell edges are not LAMMPS compatible"); + + // convert atoms to lamda coords before changing box + + domain->x2lamda(atom->nlocal); + + // convert celldata to new boxlo, boxhi, and tilt factors domain->boxhi[0] = celldata[0] + domain->boxlo[0]; domain->boxhi[1] = celldata[4] + domain->boxlo[1]; domain->boxhi[2] = celldata[8] + domain->boxlo[2]; - domain->xy = 0.0; - domain->xz = 0.0; - domain->yz = 0.0; + + domain->xy = celldata[3]; + domain->xz = celldata[6]; + domain->yz = celldata[7]; + + // reset all Domain variables that depend on box size/shape + // convert atoms coords back to new box coords + + domain->set_global_box(); + domain->set_local_box(); + domain->lamda2x(atom->nlocal); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + all(FLERR, "MDI: Unable to send cell displacement to driver"); + int ierr = MDI_Send(celldata,3,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: CELL_DISPL command + reset simulation box origin = lower-left corner +---------------------------------------------------------------------- */ void MDIEngine2::receive_celldispl() { double celldata[3]; - int ierr = MDI_Recv((char *) celldata,3,MDI_DOUBLE,mdicomm); + int ierr = MDI_Recv(celldata,3,MDI_DOUBLE,mdicomm); if (ierr) - error->all(FLERR, "MDI: Unable to receive cell displacement from driver"); - MPI_Bcast(&celldata[0],3,MPI_DOUBLE,0,world); + error->all(FLERR,"MDI: >CELL_DISPLS data"); + MPI_Bcast(celldata,3,MPI_DOUBLE,0,world); for (int icell = 0; icell < 3; icell++) celldata[icell] *= mdi2lmp_length; + // convert atoms to lamda coords before changing box + + domain->x2lamda(atom->nlocal); + + // convert celldata to new boxlo and boxhi + double old_boxlo[3]; old_boxlo[0] = domain->boxlo[0]; old_boxlo[1] = domain->boxlo[1]; old_boxlo[2] = domain->boxlo[2]; - // adjust the values of boxlo and boxhi for the new cell displacement vector - domain->boxlo[0] = celldata[0]; domain->boxlo[1] = celldata[1]; domain->boxlo[2] = celldata[2]; + domain->boxhi[0] += domain->boxlo[0] - old_boxlo[0]; domain->boxhi[1] += domain->boxlo[1] - old_boxlo[1]; domain->boxhi[2] += domain->boxlo[2] - old_boxlo[2]; + + // reset all Domain variables that depend on box origin + // convert atoms coords back to new box coords + + domain->set_global_box(); + domain->set_local_box(); + domain->lamda2x(atom->nlocal); } // ---------------------------------------------------------------------- @@ -1151,13 +1213,13 @@ void MDIEngine2::receive_celldispl() /* ---------------------------------------------------------------------- LENGTH command store received value in length_param - for use by other commands, e.g. that send strings + for use by a subsequent command, e.g. ones that send strings ---------------------------------------------------------------------- */ void MDIEngine2::length_command() { int ierr = MDI_Recv(&length_param,1,MDI_INT,mdicomm); - if (ierr) error->all(FLERR, "MDI LENGTH error"); + if (ierr) error->all(FLERR,"MDI: LENGTH data"); MPI_Bcast(&length_param,1,MPI_INT,0,world); } @@ -1171,7 +1233,7 @@ void MDIEngine2::single_command() { char *cmd = new char[length_param+1]; int ierr = MDI_Recv(cmd,length_param+1,MDI_CHAR,mdicomm); - if (ierr) error->all(FLERR,"MDI COMMAND error"); + if (ierr) error->all(FLERR,"MDI: COMMAND data"); MPI_Bcast(cmd,length_param+1,MPI_CHAR,0,world); cmd[length_param+1] = '\0'; @@ -1190,7 +1252,7 @@ void MDIEngine2::many_commands() { char *cmds = new char[length_param+1]; int ierr = MDI_Recv(cmds, length_param+1, MDI_CHAR, mdicomm); - if (ierr) error->all(FLERR, "MDI COMMANDS error"); + if (ierr) error->all(FLERR,"MDI: COMMANDS data"); MPI_Bcast(cmds,length_param+1,MPI_CHAR,0,world); cmds[length_param+1] = '\0'; @@ -1209,7 +1271,7 @@ void MDIEngine2::infile() { char *infile = new char[length_param+1]; int ierr = MDI_Recv(infile,length_param+1,MDI_CHAR,mdicomm); - if (ierr) error->all(FLERR, "MDI INFILE error"); + if (ierr) error->all(FLERR,"MDI: INFILE data"); MPI_Bcast(infile,length_param+1,MPI_CHAR,0,world); infile[length_param+1] = '\0'; @@ -1218,7 +1280,15 @@ void MDIEngine2::infile() delete [] infile; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + EVAL command + compute forces, energy, pressure of current system + can be called multiple times by driver + for a system that is continuously evolving + distinguishes between first-time call vs + system needs reneighboring vs system does not need reneighboring + does not increment timestep +---------------------------------------------------------------------- */ void MDIEngine2::evaluate() { @@ -1228,6 +1298,9 @@ void MDIEngine2::evaluate() update->integrate->setup(0); } else { + + // insure potential energy and virial are tallied on this step + pe->addstep(update->ntimestep); press->addstep(update->ntimestep); @@ -1241,7 +1314,14 @@ void MDIEngine2::evaluate() } } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + RESET_BOX command + wrapper on library reset_box() method + requires no atoms exist + allows caller to define a new simulation box + NOTE: this will not work in plugin mode, b/c of 3 MDI_Recv() calls + how to effectively do this? +---------------------------------------------------------------------- */ void MDIEngine2::reset_box() { @@ -1258,7 +1338,16 @@ void MDIEngine2::reset_box() lammps_reset_box(lmp,boxlo,boxhi,tilts[0],tilts[1],tilts[2]); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + CREATE_ATOM command + wrapper on library create_atoms() method + requires simulation box be defined + allows caller to define a new set of atoms + with their IDs, types, coords, velocities, image flags + NOTE: this will not work in plugin mode, b/c of multiple MDI_Recv() calls + how to effectively do this? + NOTE: also the memory here is not yet allocated correctly +---------------------------------------------------------------------- */ void MDIEngine2::create_atoms() { @@ -1305,15 +1394,18 @@ void MDIEngine2::create_atoms() } if (!x || !type) - error->all(FLERR,"MDI create_atoms: did not receive atom coords or types"); + error->all(FLERR,"MDI: CREATE_ATOM did not receive atom coords or types"); int ncreate = lammps_create_atoms(lmp,natoms,id,type,x,v,image,1); if (ncreate != natoms) - error->all(FLERR, "MDI create_atoms: created atoms != sent atoms"); + error->all(FLERR, "MDI: CREATE ATOM created atoms != sent atoms"); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + all(FLERR, "MDI: Unable to send pressure to driver"); + if (ierr) error->all(FLERR,"MDI: compute_vector(); for (int i = 0; i < 6; i++) - pvector[i] = press->vector[i] * lmp2mdi_pressure; + ptensor[i] = press->vector[i] * lmp2mdi_pressure; - int ierr = MDI_Send(pvector,6,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR, "MDI: Unable to send ptensor to driver"); + int ierr = MDI_Send(ptensor,6,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: atom->natoms) return; + if (atom->natoms <= maxbuf) return; if (3*atom->natoms > MAXSMALLINT) - error->all(FLERR,"Natoms is too large to use with mdi engine"); + error->all(FLERR,"Natoms too large to use with mdi engine"); maxbuf = atom->natoms; @@ -1369,6 +1468,9 @@ void MDIEngine2::reallocate() memory->create(buf3all,3*maxbuf,"mdi:buf3all"); } +/* ---------------------------------------------------------------------- + MDI to/from LAMMPS conversion factors +------------------------------------------------------------------------- */ void MDIEngine2::unit_conversions() { diff --git a/src/MDI/mdi_engine2.h b/src/MDI/mdi_engine2.h index aeda3f823c..4ef487ad06 100644 --- a/src/MDI/mdi_engine2.h +++ b/src/MDI/mdi_engine2.h @@ -57,7 +57,7 @@ class MDIEngine2 : public Command { int length_param; // LENGTH command value used by other commands - // unit conversion factors; + // unit conversion factors double lmp2mdi_length,mdi2lmp_length; double lmp2mdi_energy,mdi2lmp_energy; @@ -87,12 +87,12 @@ class MDIEngine2 : public Command { void send_double1(int); void send_int1(int); void send_double3(int); - void send_labels(); - void send_energy(); + void send_total_energy(); void send_pe(); void send_ke(); + void send_cell(); void receive_cell(); void send_celldispl();