diff --git a/src/USER-MDI/fix_mdi_engine.cpp b/src/USER-MDI/fix_mdi_engine.cpp index bec2272cd9..3987e84f5b 100644 --- a/src/USER-MDI/fix_mdi_engine.cpp +++ b/src/USER-MDI/fix_mdi_engine.cpp @@ -34,38 +34,31 @@ #include #include +enum{NONE,REAL,METAL}; // LAMMPS units which MDI supports + using namespace LAMMPS_NS; using namespace FixConst; -/****************************************************************************/ +/* ---------------------------------------------------------------------- */ - -/*************************************************************** - * create class and parse arguments in LAMMPS script. Syntax: - * fix ID group-ID mdi_engine [couple ] - ***************************************************************/ FixMDIEngine::FixMDIEngine(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), id_pe(NULL), pe(NULL), id_ke(NULL), ke(NULL) { + if (narg != 3) error->all(FLERR,"Illegal fix mdi command"); - if (narg > 3) - error->all(FLERR,"Illegal fix mdi command"); + // NOTE: real & metal the 2 atomic-scale units LAMMPS has + // I suggest LAMMPS for MDI support both + // real: coords = Ang, eng = Kcal/mole, force = Kcal/mole/Ang + // metal: coords = Ang, eng = eV, force = eV/Ang - // allocate arrays - memory->create(add_force,3*atom->natoms,"mdi:add_force"); - for (int i=0; i< 3*atom->natoms; i++) { - add_force[i] = 0.0; - } + lmpunits = NONE; + if (strcmp(update->unit_style,"real") == 0) lmpunits = REAL; + if (strcmp(update->unit_style,"metal") == 0) lmpunits = METAL; + if (lmpunits == NONE) error->all(FLERR,"MDI requires real or metal units"); - // create a new compute pe style - // id = fix-ID + pe, compute group = all - - master = (comm->me==0) ? 1 : 0; - - // create instance of the Irregular class - irregular = new Irregular(lmp); + // MDI setup most_recent_init = 0; exit_flag = false; @@ -77,50 +70,55 @@ FixMDIEngine::FixMDIEngine(LAMMPS *lmp, int narg, char **arg) : strncpy(target_node, "\0", MDI_COMMAND_LENGTH); strncpy(current_node, "@DEFAULT", MDI_COMMAND_LENGTH); - // create and add a compute for PE calculation + // accept a communicator to the driver + // master = 1 for proc 0, otherwise 0 + + master = (comm->me==0) ? 1 : 0; + + if (master) { + MDI_Accept_Communicator(&driver_socket); + if (driver_socket <= 0) error->all(FLERR,"Unable to connect to driver"); + } else driver_socket = 0; + + // create computes for KE and PE + id_pe = utils::strdup(std::string(id) + "_pe"); modify->add_compute(fmt::format("{} all pe",id_pe)); - // create and add a compute for KE calculation id_ke = utils::strdup(std::string(id) + "_ke"); modify->add_compute(fmt::format("{} all ke",id_ke)); - // accept a communicator to the driver - int ierr; - if (master) { - MDI_Accept_Communicator(&driver_socket); - if (driver_socket <= 0) - error->all(FLERR,"Unable to connect to driver"); - } else driver_socket=0; + // irregular class and data structs used by MDI -} - -/********************************* - * Clean up on deleting the fix. * - *********************************/ -FixMDIEngine::~FixMDIEngine() -{ - modify->delete_compute(id_pe); - modify->delete_compute(id_ke); - delete irregular; - delete [] id_pe; - delete [] id_ke; - delete [] target_command; - delete [] command; + irregular = new Irregular(lmp); + add_force = NULL; } /* ---------------------------------------------------------------------- */ + +FixMDIEngine::~FixMDIEngine() +{ + delete [] target_command; + delete [] command; + delete [] current_node; + delete [] target_node; + + modify->delete_compute(id_pe); + modify->delete_compute(id_ke); + delete irregular; + memory->destroy(add_force); +} + +/* ---------------------------------------------------------------------- */ + int FixMDIEngine::setmask() { int mask = 0; - // MD masks mask |= POST_INTEGRATE; mask |= PRE_REVERSE; mask |= POST_FORCE; - - // Minimizer masks - mask |= MIN_PRE_FORCE; + mask |= MIN_PRE_FORCE; // NOTE: whack this? mask |= MIN_POST_FORCE; return mask; @@ -134,67 +132,76 @@ void FixMDIEngine::exchange_forces() const int * const mask = atom->mask; const int nlocal = atom->nlocal; - // add the forces from the driver - for (int i=0; i < nlocal; ++i) { + // add forces from the driver + + for (int i = 0; i < nlocal; ++i) { if (mask[i] & groupbit) { f[i][0] += add_force[3*(atom->tag[i]-1)+0]; f[i][1] += add_force[3*(atom->tag[i]-1)+1]; f[i][2] += add_force[3*(atom->tag[i]-1)+2]; } } - } /* ---------------------------------------------------------------------- */ void FixMDIEngine::init() { - // Confirm that the required computes are available + // confirm that two required computes are still available + int icompute_pe = modify->find_compute(id_pe); if (icompute_pe < 0) - error->all(FLERR,"Potential energy ID for fix mdi does not exist"); + error->all(FLERR,"Potential energy ID for fix mdi/engine does not exist"); int icompute_ke = modify->find_compute(id_ke); if (icompute_pe < 0) - error->all(FLERR,"Kinetic energy ID for fix mdi does not exist"); + error->all(FLERR,"Kinetic energy ID for fix mdi/engine does not exist"); pe = modify->compute[icompute_pe]; ke = modify->compute[icompute_ke]; - return; + // one-time allocation of add_force array + // NOTE: moved this here b/c natoms may not be defined when Fix is constructed + // NOTE: check that 3*natoms does not overflow a 32-bit int + if (!add_force) { + memory->create(add_force,3*atom->natoms,"mdi/engine:add_force"); + for (int i = 0; i < 3*atom->natoms; i++) add_force[i] = 0.0; + } } /* ---------------------------------------------------------------------- */ void FixMDIEngine::setup(int vflag) { - //compute the potential energy + // NOTE: this seems an odd place to compute these + // I think it would be better to compute these on-demand + // in response to a driver request? + potential_energy = pe->compute_scalar(); kinetic_energy = ke->compute_scalar(); // trigger potential energy computation on next timestep - pe->addstep(update->ntimestep+1); - ke->addstep(update->ntimestep+1); + // NOTE: there is no triggering needed for KE - if ( most_recent_init == 1 ) { // md - // @PRE-FORCES - engine_mode("@PRE-FORCES"); - } + pe->addstep(update->ntimestep+1); + + if (most_recent_init == 1) engine_mode("@PRE-FORCES"); } /* ---------------------------------------------------------------------- */ void FixMDIEngine::min_setup(int vflag) { + // NOTE: this seems an odd place to compute these + potential_energy = pe->compute_scalar(); kinetic_energy = ke->compute_scalar(); - // trigger potential energy computation on next timestep - pe->addstep(update->ntimestep+1); - ke->addstep(update->ntimestep+1); - - // @FORCES engine_mode("@FORCES"); + + // trigger potential energy computation on next timestep + + pe->addstep(update->ntimestep+1); } /* ---------------------------------------------------------------------- */ @@ -208,24 +215,22 @@ void FixMDIEngine::post_integrate() void FixMDIEngine::pre_reverse(int eflag, int vflag) { - // calculate the energy + // NOTE: this seems an odd place to compute these + potential_energy = pe->compute_scalar(); kinetic_energy = ke->compute_scalar(); - // @PRE-FORCES engine_mode("@PRE-FORCES"); // trigger potential energy computation on next timestep - pe->addstep(update->ntimestep+1); - ke->addstep(update->ntimestep+1); -} + pe->addstep(update->ntimestep+1); +} /* ---------------------------------------------------------------------- */ void FixMDIEngine::min_pre_force(int vflag) { - // @COORDS engine_mode("@COORDS"); } @@ -233,35 +238,46 @@ void FixMDIEngine::min_pre_force(int vflag) void FixMDIEngine::min_post_force(int vflag) { - // calculate the energy + // NOTE: this seems an odd place to compute these + potential_energy = pe->compute_scalar(); kinetic_energy = ke->compute_scalar(); - // @FORCES engine_mode("@FORCES"); // trigger potential energy computation on next timestep + pe->addstep(update->ntimestep+1); - ke->addstep(update->ntimestep+1); } /* ---------------------------------------------------------------------- */ void FixMDIEngine::post_force(int vflag) { - if ( most_recent_init == 1 ) { // md - // @FORCES - engine_mode("@FORCES"); - } - else if ( most_recent_init == 2 ) { // optg - // @FORCES - engine_mode("@FORCES"); - } + if (most_recent_init == 1) engine_mode("@FORCES"); + else if (most_recent_init == 2) engine_mode("@FORCES"); + + // NOTE: should this also be done in this method? + // trigger potential energy computation on next timestep + // NOTE: in general, forcing the pair styles to compute PE every step + // is inefficient, would be better to think of another way to do this, + // e.g. at very end of a step, + // possibly by mdi_engine when knows it is needed + + pe->addstep(update->ntimestep+1); } -/* ---------------------------------------------------------------------- */ +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- +// rest of file processes and responds to MDI driver commands +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- -int FixMDIEngine::execute_command(const char *command, MDI_Comm driver_socket) +/* ---------------------------------------------------------------------- + process a single command from driver +---------------------------------------------------------------------- */ + +int FixMDIEngine::execute_command(const char *command, MDI_Comm mdicomm) { /* if (screen) @@ -270,144 +286,139 @@ int FixMDIEngine::execute_command(const char *command, MDI_Comm driver_socket) fprintf(logfile,"MDI command: %s:\n",command); */ - // confirm that this command is supported at this node + // confirm this command is supported at this node + int command_exists = 0; - ierr = MDI_Check_Command_Exists(current_node, command, MDI_COMM_NULL, &command_exists); + ierr = MDI_Check_Command_Exists(current_node, command, + MDI_COMM_NULL, &command_exists); if (ierr != 0) - error->all(FLERR,"Unable to check whether the current command is supported"); - if ( command_exists != 1 ) - error->all(FLERR,"Received a command that is not supported at the current node"); + error->all(FLERR,"MDI: Unable to check whether current command is supported"); + if (command_exists != 1) + error->all(FLERR,"MDI: Received an unsupported at current node"); + + // respond to any driver command + + // send calculation status to the driver + // NOTE: why does STATUS have extra spaces? if (strcmp(command,"STATUS ") == 0 ) { - // send the calculation status to the driver if (master) { - ierr = MDI_Send_Command("READY", driver_socket); + ierr = MDI_Send_Command("READY", mdicomm); if (ierr != 0) - error->all(FLERR,"Unable to return status to driver"); + error->all(FLERR,"MDI: Unable to return status to driver"); } - } - else if (strcmp(command,">NATOMS") == 0 ) { - // receive the number of atoms from the driver + + } else if (strcmp(command,">NATOMS") == 0 ) { if (master) { - ierr = MDI_Recv((char*) &atom->natoms, 1, MDI_INT, driver_socket); + ierr = MDI_Recv((char*) &atom->natoms, 1, MDI_INT, mdicomm); if (ierr != 0) - error->all(FLERR,"Unable to receive number of atoms from driver"); + error->all(FLERR,"MDI: Unable to receive number of atoms from driver"); } MPI_Bcast(&atom->natoms,1,MPI_INT,0,world); - } - else if (strcmp(command,"natoms; - ierr = MDI_Send((char*) &mdi_natoms, 1, MDI_INT64_T, driver_socket); + ierr = MDI_Send((char*) &mdi_natoms, 1, MDI_INT64_T, mdicomm); if (ierr != 0) - error->all(FLERR,"Unable to send number of atoms to driver"); - } - } - else if (strcmp(command,"ntypes, 1, MDI_INT, driver_socket); - if (ierr != 0) - error->all(FLERR,"Unable to send number of atom types to driver"); - } - } - else if (strcmp(command,"COORDS") == 0 ) { - // receive the coordinate information - receive_coordinates(error); - } - else if (strcmp(command,"FORCES") == 0 ) { - // receive the forces from the driver - receive_forces(error); - } - else if (strcmp(command,"+FORCES") == 0 ) { - // receive additional forces from the driver - // these are added prior to SHAKE or other post-processing - add_forces(error); - } - else if (strcmp(command,"@INIT_MD") == 0 ) { - if ( most_recent_init != 0 ) { - error->all(FLERR,"MDI is already performing a simulation"); + error->all(FLERR,"MDI: Unable to send number of atoms to driver"); } - // initialize a new MD simulation + } else if (strcmp(command,"ntypes, 1, MDI_INT, mdicomm); + if (ierr != 0) + error->all(FLERR,"MDI: Unable to send number of atom types to driver"); + } + + } else if (strcmp(command,"CELL command + + } else if (strcmp(command,"COORDS") == 0 ) { + receive_coordinates(error); + + } else if (strcmp(command,"FORCES") == 0 ) { + receive_forces(error); + + // receive additional forces from the driver + // these can be added prior to SHAKE or other post-processing + // NOTE: maybe this is now not necessary? + + } else if (strcmp(command,"+FORCES") == 0 ) { + add_forces(error); + + // initialize new MD simulation or minimization + // return control to return to mdi_engine + + } else if (strcmp(command,"@INIT_MD") == 0 ) { + if (most_recent_init != 0) + error->all(FLERR,"MDI: MDI is already performing a simulation"); most_recent_init = 1; local_exit_flag = true; - } - else if (strcmp(command,"@INIT_OPTG") == 0 ) { - if ( most_recent_init != 0 ) { - error->all(FLERR,"MDI is already performing a simulation"); - } - // initialize a new geometry optimization + // initialize new energy minimization + // return control to return to mdi_engine + + } else if (strcmp(command,"@INIT_OPTG") == 0 ) { + if ( most_recent_init != 0 ) + error->all(FLERR,"MDI: MDI is already performing a simulation"); most_recent_init = 2; local_exit_flag = true; - //optg_init(error); - } - else if (strcmp(command,"@") == 0 ) { + + } else if (strcmp(command,"@") == 0 ) { strncpy(target_node, "\0", MDI_COMMAND_LENGTH); local_exit_flag = true; - } - else if (strcmp(command,"<@") == 0 ) { + + } else if (strcmp(command,"<@") == 0 ) { if (master) { - ierr = MDI_Send(current_node, MDI_NAME_LENGTH, MDI_CHAR, driver_socket); - if (ierr != 0) - error->all(FLERR,"Unable to send node to driver"); + ierr = MDI_Send(current_node, MDI_NAME_LENGTH, MDI_CHAR, mdicomm); + if (ierr != 0) error->all(FLERR,"MDI: Unable to send node to driver"); } - } - else if (strcmp(command,"max_eval = 0; } - } - else if (strcmp(command,"EXIT") == 0 ) { + + } else if (strcmp(command,"EXIT") == 0 ) { // exit the driver code exit_flag = true; @@ -435,10 +446,9 @@ int FixMDIEngine::execute_command(const char *command, MDI_Comm driver_socket) // set the maximum number of force evaluations to 0 update->max_eval = 0; } - } - else { - // the command is not supported - error->all(FLERR,"Unknown command from driver"); + + } else { + error->all(FLERR,"MDI: Unknown command from driver"); } return 0; @@ -455,70 +465,81 @@ char *FixMDIEngine::engine_mode(const char *node) fprintf(logfile,"MDI ENGINE MODE: %i\n",node); */ - // flag to indicate whether the engine should continue listening for commands at this node - strncpy(current_node, node, MDI_COMMAND_LENGTH); - if ( strcmp(target_node,"\0") != 0 and strcmp(target_node, current_node) != 0 ) { + // do not process commands if engine and driver are not at same node + // target_node = node that driver has set via a @ command + // current_node = node that engine (LAMMPS) has set + + strncpy(current_node,node,MDI_COMMAND_LENGTH); + if (strcmp(target_node,"\0") != 0 && strcmp(target_node,current_node) != 0) local_exit_flag = true; - } // register the execute_command function with MDI + // NOTE: does this need to be done multiple times ?? + MDI_Set_execute_command_func(lammps_execute_mdi_command, this); - /* ----------------------------------------------------------------- */ - // Answer commands from the driver - /* ----------------------------------------------------------------- */ + // respond to commands from the driver while (not exit_flag and not local_exit_flag) { // read the next command from the driver + // NOTE: all procs call this, but only proc 0 receives command? + ierr = MDI_Recv_Command(command, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to receive command from driver"); + error->all(FLERR,"MDI: Unable to receive command from driver"); + + // broadcast command to the other MPI tasks - // broadcast the command to the other tasks MPI_Bcast(command,MDI_COMMAND_LENGTH,MPI_CHAR,0,world); // execute the command + this->execute_command(command, driver_socket); // check if the target node is something other than the current node - if ( strcmp(target_node,"\0") != 0 and strcmp(target_node, current_node) != 0 ) { - local_exit_flag = true; - } + if (strcmp(target_node,"\0") != 0 && strcmp(target_node,current_node) != 0 ) + local_exit_flag = true; } - // a local exit has completed, so turn off the local exit flag + // local exit occured so turn off local exit flag + local_exit_flag = false; return command; - } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::receive_coordinates(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The >COORDS MDI command does not support lj units"); - double posconv; - double angstrom_to_bohr; - MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); - posconv=force->angstrom/angstrom_to_bohr; - // create a buffer to hold the coordinates + // NOTE: logic like this everywhere else + + if (lmpunits == REAL) { + double angstrom_to_bohr; + MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); + posconv=force->angstrom/angstrom_to_bohr; + } else if (lmpunits == METAL) { + // ?? + } + + // create buffer to hold all coords + double *buffer; buffer = new double[3*atom->natoms]; if (master) { ierr = MDI_Recv((char*) buffer, 3*atom->natoms, MDI_DOUBLE, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to receive coordinates from driver"); + error->all(FLERR,"MDI: Unable to receive coordinates from driver"); } MPI_Bcast(buffer,3*atom->natoms,MPI_DOUBLE,0,world); // pick local atoms from the buffer + double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -531,13 +552,15 @@ void FixMDIEngine::receive_coordinates(Error* error) // ensure atoms are in current box & update box via shrink-wrap // has to be be done before invoking Irregular::migrate_atoms() // since it requires atoms be inside simulation box + if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); if (domain->triclinic) domain->lamda2x(atom->nlocal); - // move atoms to new processors via irregular() - // only needed if migrate_check() says an atom moves to far + // move atoms to new processors via irregular() only needed if + // migrate_check() says an atom moves too far + if (domain->triclinic) domain->x2lamda(atom->nlocal); if (irregular->migrate_check()) irregular->migrate_atoms(); if (domain->triclinic) domain->lamda2x(atom->nlocal); @@ -545,13 +568,10 @@ void FixMDIEngine::receive_coordinates(Error* error) delete [] buffer; } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_coordinates(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The create(coords,3*atom->natoms,"mdi:coords"); + coords = new double[3*atom->natoms]; coords_reduced = new double[3*atom->natoms]; - // zero the coordinates array - for (int icoord = 0; icoord < 3*atom->natoms; icoord++) { - coords[icoord] = 0.0; - } + // zero coords + + for (int icoord = 0; icoord < 3*atom->natoms; icoord++) coords[icoord] = 0.0; + + // copy local atoms into buffer at correct locations - // pick local atoms from the buffer double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -578,18 +602,20 @@ void FixMDIEngine::send_coordinates(Error* error) coords[3*(atom->tag[i]-1)+2] = x[i][2]/posconv; } - MPI_Reduce(coords, coords_reduced, 3*atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world); + MPI_Reduce(coords,coords_reduced,3*atom->natoms,MPI_DOUBLE,MPI_SUM,0,world); if (master) { - ierr = MDI_Send((char*) coords_reduced, 3*atom->natoms, MDI_DOUBLE, driver_socket); + ierr = MDI_Send((char*) coords_reduced,3*atom->natoms,MDI_DOUBLE, + driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to send coordinates to driver"); + error->all(FLERR,"MDI: Unable to send coordinates to driver"); } delete [] coords; delete [] coords_reduced; } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_charges(Error* error) { @@ -600,11 +626,11 @@ void FixMDIEngine::send_charges(Error* error) charges_reduced = new double[atom->natoms]; // zero the charges array - for (int icharge = 0; icharge < atom->natoms; icharge++) { - charges[icharge] = 0.0; - } + + for (int icharge = 0; icharge < atom->natoms; icharge++) charges[icharge] = 0.0; // pick local atoms from the buffer + double *charge = atom->q; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -617,7 +643,7 @@ void FixMDIEngine::send_charges(Error* error) if (master) { ierr = MDI_Send((char*) charges_reduced, atom->natoms, MDI_DOUBLE, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to send charges to driver"); + error->all(FLERR,"MDI: Unable to send charges to driver"); } delete [] charges; @@ -625,15 +651,15 @@ void FixMDIEngine::send_charges(Error* error) } +/* ---------------------------------------------------------------------- */ + void FixMDIEngine::send_energy(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The all(FLERR,"Unable to send potential energy to driver"); + error->all(FLERR,"MDI: Unable to send potential energy to driver"); } } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_pe(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The all(FLERR,"Unable to send potential energy to driver"); + error->all(FLERR,"MDI: Unable to send potential energy to driver"); } } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_ke(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The all(FLERR,"Unable to send potential energy to driver"); + error->all(FLERR,"MDI: Unable to send potential energy to driver"); } } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_types(Error* error) { int * const type = atom->type; + // NOTE: why is this not supported? + // maybe MDI labels = LAMMPS types? + if (master) { ierr = MDI_Send((char*) type, atom->natoms, MDI_INT, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to send atom types to driver"); + error->all(FLERR,"MDI: Unable to send atom types to driver"); } } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_labels(Error* error) { @@ -724,12 +753,13 @@ void FixMDIEngine::send_labels(Error* error) if (master) { ierr = MDI_Send( labels, atom->natoms * MDI_LABEL_LENGTH, MDI_CHAR, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to send atom types to driver"); + error->all(FLERR,"MDI: Unable to send atom types to driver"); } delete [] labels; } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_masses(Error* error) { @@ -745,6 +775,7 @@ void FixMDIEngine::send_masses(Error* error) } // determine the atomic masses + if (rmass) { for (int iatom=0; iatom < nlocal; iatom++) { mass_by_atom[ atom->tag[iatom] - 1 ] = rmass[iatom]; @@ -759,23 +790,21 @@ void FixMDIEngine::send_masses(Error* error) MPI_Reduce(mass_by_atom, mass_by_atom_reduced, atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world); // send the atomic masses to the driver + if (master) { ierr = MDI_Send((char*) mass_by_atom_reduced, atom->natoms, MDI_DOUBLE, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to send atom masses to driver"); + error->all(FLERR,"MDI: Unable to send atom masses to driver"); } delete [] mass_by_atom; delete [] mass_by_atom_reduced; } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_forces(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The natoms; iforce++) { - forces[iforce] = 0.0; - } + + for (int iforce = 0; iforce < 3*atom->natoms; iforce++) forces[iforce] = 0.0; // if not at a node, calculate the forces + if ( strcmp(current_node, "@DEFAULT") == 0 ) { // certain fixes, such as shake, move the coordinates // to ensure that the coordinates do not change, store a copy @@ -818,6 +847,18 @@ void FixMDIEngine::send_forces(Error* error) update->nsteps = 1; lmp->init(); update->integrate->setup_minimal(1); + + // NOTE: can this be done here instead of below? + + if ( strcmp(current_node, "@DEFAULT") == 0 ) { + // restore the original set of coordinates + double **x_new = atom->x; + for (int i = 0; i < nlocal; i++) { + x_new[i][0] = x_buf[3*i+0]; + x_new[i][1] = x_buf[3*i+1]; + x_new[i][2] = x_buf[3*i+2]; + } + } } // pick local atoms from the buffer @@ -835,17 +876,7 @@ void FixMDIEngine::send_forces(Error* error) if (master) { ierr = MDI_Send((char*) forces_reduced, 3*atom->natoms, MDI_DOUBLE, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to send atom forces to driver"); - } - - if ( strcmp(current_node, "@DEFAULT") == 0 ) { - // restore the original set of coordinates - double **x_new = atom->x; - for (int i = 0; i < nlocal; i++) { - x_new[i][0] = x_buf[3*i+0]; - x_new[i][1] = x_buf[3*i+1]; - x_new[i][2] = x_buf[3*i+2]; - } + error->all(FLERR,"MDI: Unable to send atom forces to driver"); } delete [] forces; @@ -854,13 +885,10 @@ void FixMDIEngine::send_forces(Error* error) } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::receive_forces(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The >FORCES MDI command does not support lj units"); - double angstrom_to_bohr; MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); double kelvin_to_hartree; @@ -877,7 +905,7 @@ void FixMDIEngine::receive_forces(Error* error) if (master) { ierr = MDI_Recv((char*) forces, 3*atom->natoms, MDI_DOUBLE, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to receive atom forces to driver"); + error->all(FLERR,"MDI: Unable to receive atom forces to driver"); } MPI_Bcast(forces,3*atom->natoms,MPI_DOUBLE,0,world); @@ -895,13 +923,15 @@ void FixMDIEngine::receive_forces(Error* error) delete [] forces; } +/* ---------------------------------------------------------------------- */ + +// NOTE: if keeping add_forces (see NOTE above) +// then could make one replace_add_forces method with an extra "mode" arg +// for replace or add +// since these 2 methods are nearly identical void FixMDIEngine::add_forces(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The >+FORCES MDI command does not support lj units"); - double angstrom_to_bohr; MDI_Conversion_Factor("angstrom", "bohr", &angstrom_to_bohr); double kelvin_to_hartree; @@ -918,7 +948,7 @@ void FixMDIEngine::add_forces(Error* error) if (master) { ierr = MDI_Recv((char*) forces, 3*atom->natoms, MDI_DOUBLE, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to receive atom forces to driver"); + error->all(FLERR,"MDI: Unable to receive atom forces to driver"); } MPI_Bcast(forces,3*atom->natoms,MPI_DOUBLE,0,world); @@ -936,13 +966,10 @@ void FixMDIEngine::add_forces(Error* error) delete [] forces; } +/* ---------------------------------------------------------------------- */ void FixMDIEngine::send_cell(Error* error) { - // unable to convert to atomic units if LAMMPS is using lj units - if (strcmp(update->unit_style,"lj") == 0) - error->all(FLERR,"The angstrom * angstrom_to_bohr; for (int icell=0; icell < 9; icell++) { celldata[icell] *= unit_conv; @@ -972,6 +1000,6 @@ void FixMDIEngine::send_cell(Error* error) if (master) { ierr = MDI_Send((char*) celldata, 9, MDI_DOUBLE, driver_socket); if (ierr != 0) - error->all(FLERR,"Unable to send cell dimensions to driver"); + error->all(FLERR,"MDI: Unable to send cell dimensions to driver"); } } diff --git a/src/USER-MDI/fix_mdi_engine.h b/src/USER-MDI/fix_mdi_engine.h index d6509cf29d..11f5c0e900 100644 --- a/src/USER-MDI/fix_mdi_engine.h +++ b/src/USER-MDI/fix_mdi_engine.h @@ -37,6 +37,7 @@ class FixMDIEngine : public Fix { char *engine_mode(const char *node); // receive and update forces + void setup(int); void min_setup(int); void post_integrate(); @@ -50,12 +51,14 @@ class FixMDIEngine : public Fix { double kinetic_energy; // stores kinetic energy // current command + char *command; protected: void exchange_forces(); private: + int lmpunits; // REAL or METAL int master, ierr; int driver_socket; int most_recent_init; // which MDI init command was most recently received? @@ -74,6 +77,7 @@ class FixMDIEngine : public Fix { // -2 - after MD_INIT command followed by @PRE-FORCES (actually @INIT_OPTG?) // command to be executed at the target node + char *target_command; char *id_pe; @@ -82,6 +86,7 @@ class FixMDIEngine : public Fix { class Minimize *minimizer; class Compute *pe; class Compute *ke; + void send_types(Error *); void send_labels(Error *); void send_masses(Error *); @@ -95,7 +100,6 @@ class FixMDIEngine : public Fix { void add_forces(Error *); void receive_forces(Error *); void send_cell(Error *); - }; } diff --git a/src/USER-MDI/mdi_engine.cpp b/src/USER-MDI/mdi_engine.cpp index 1fdf52f950..32e255a525 100644 --- a/src/USER-MDI/mdi_engine.cpp +++ b/src/USER-MDI/mdi_engine.cpp @@ -36,22 +36,19 @@ using namespace LAMMPS_NS; -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + trigger LAMMPS to start acting as an MDI engine + endlessly loop over receiving commands from driver and responding + much of the logic for this is in FixMDIEngine + when EXIT command is received, mdi_engine command exits +---------------------------------------------------------------------- */ -CommandMDIEngine::CommandMDIEngine(LAMMPS *lmp) : Command(lmp) { - return; -} - -CommandMDIEngine::~CommandMDIEngine() { - return; -} - -/* ---------------------------------------------------------------------- */ - -void CommandMDIEngine::command(int narg, char **arg) +void MDIEngine::command(int narg, char **arg) { + // list of nodes and commands that a MDI-compliant MD code should support + + // default node and its commands - // register the default node MDI_Register_Node("@DEFAULT"); MDI_Register_Command("@DEFAULT", "<@"); MDI_Register_Command("@DEFAULT", "FORCES"); MDI_Register_Command("@FORCES", "<@"); @@ -153,7 +155,8 @@ void CommandMDIEngine::command(int narg, char **arg) MDI_Register_Command("@FORCES", "@PRE-FORCES"); MDI_Register_Command("@FORCES", "EXIT"); - // register the coordinates node + // node at POST_INTEGRATE location in timestep + MDI_Register_Node("@COORDS"); MDI_Register_Command("@COORDS", "<@"); MDI_Register_Command("@COORDS", "find_fix_by_style("mdi/engine"); bool added_mdi_engine_fix = false; if (ifix < 0) { - modify->add_fix("_mdiengine_ all mdi/engine"); + modify->add_fix("MDI_ENGINE_INTERNAL all mdi/engine"); added_mdi_engine_fix = true; } // identify the mdi_engine fix + ifix = modify->find_fix_by_style("mdi/engine"); - if (ifix < 0) error->all(FLERR,"The mdi_engine command requires the mdi/engine fix"); mdi_fix = static_cast(modify->fix[ifix]); - /* format for MDI Engine command: - * mdi_engine - */ - if (narg > 0) error->all(FLERR,"Illegal MDI command"); + // check that LAMMPS is setup as a compatible MDI engine + + if (narg > 0) error->all(FLERR,"Illegal mdi_engine command"); if (atom->tag_enable == 0) - error->all(FLERR,"Cannot use MDI command without atom IDs"); + error->all(FLERR,"Cannot use mdi_engine without atom IDs"); if (atom->tag_consecutive() == 0) - error->all(FLERR,"MDI command requires consecutive atom IDs"); + error->all(FLERR,"mdi_engine requires consecutive atom IDs"); + + // endless engine loop, responding to driver commands + + char *command; + + while (1) { + + // mdi_engine command only recognizes three nodes + // DEFAULT, INIT_MD, INIT_OPTG - // begin engine_mode - char *command = NULL; - while ( true ) { - // listen for MDI commands at the default command - // the response to most MDI commands is handled here command = mdi_fix->engine_mode("@DEFAULT"); + + // MDI commands for dynamics or minimization - // MDI commands that involve large-scale program flow are handled here if (strcmp(command,"@INIT_MD") == 0 ) { - // enter MD control loop - int received_exit = mdi_md(); - if ( received_exit == 1 ) { - return; - } - } - if (strcmp(command,"@INIT_OPTG") == 0 ) { - // enter minimizer control loop - int received_exit = mdi_optg(); - if ( received_exit == 1 ) { - return; - } - } - else if (strcmp(command,"EXIT") == 0 ) { - return; - } - else { - error->all(FLERR,fmt::format("MDI node exited with invalid command: {}",command)); - } + command = mdi_md(); + if (strcmp(command,"EXIT")) break; + + } else if (strcmp(command,"@INIT_OPTG") == 0 ) { + command = mdi_optg(); + if (strcmp(command,"EXIT")) break; + + } else if (strcmp(command,"EXIT") == 0) { + break; + + } else + error->all(FLERR, + fmt::format("MDI node exited with " + "invalid command: {}",command)); } - // remove the mdi/engine fix - if (added_mdi_engine_fix) modify->delete_fix("_mdiengine_"); + // remove mdi/engine fix that mdi_engine instantiated - return; + if (added_mdi_engine_fix) modify->delete_fix("MDI_ENGINE_INTERNAL"); } +/* ---------------------------------------------------------------------- + run an MD simulation under control of driver +---------------------------------------------------------------------- */ - -int CommandMDIEngine::mdi_md() +char *MDIEngine::mdi_md() { // initialize an MD simulation - update->whichflag = 1; // 1 for dynamics + + update->whichflag = 1; timer->init_timeout(); update->nsteps = 1; update->ntimestep = 0; @@ -248,119 +253,109 @@ int CommandMDIEngine::mdi_md() update->laststep = update->ntimestep + update->nsteps; update->beginstep = update->firststep; update->endstep = update->laststep; + lmp->init(); - // the MD simulation is now at the @INIT_MD node + // engine is now at @INIT_MD node + char *command = NULL; command = mdi_fix->engine_mode("@INIT_MD"); - if (strcmp(command,"@DEFAULT") == 0 ) { - // return, and flag for @DEFAULT node - return 0; - } - else if (strcmp(command,"EXIT") == 0 ) { - // return, and flag for global exit - return 1; - } + if (strcmp(command,"@DEFAULT") == 0 || strcmp(command,"EXIT") == 0) + return command; + + // setup the MD simulation - // continue the MD simulation update->integrate->setup(1); - // the MD simulation is now at the @FORCES node command = mdi_fix->engine_mode("@FORCES"); - if (strcmp(command,"@DEFAULT") == 0 ) { - // return, and flag for @DEFAULT node - return 0; - } - else if (strcmp(command,"EXIT") == 0 ) { - // return, and flag for global exit - return 1; - } + if (strcmp(command,"@DEFAULT") == 0 || strcmp(command,"EXIT") == 0) + return command; - // do MD iterations until told to exit - while ( true ) { + // run MD one step at a time - // run an MD timestep - update->whichflag = 1; // 1 for dynamics + while (1) { + update->whichflag = 1; timer->init_timeout(); update->nsteps += 1; update->laststep += 1; update->endstep = update->laststep; output->next = update->ntimestep + 1; + + // single MD timestep + update->integrate->run(1); - // get the most recent command the MDI engine received + // done with MD if driver sends @DEFAULT or EXIT + command = mdi_fix->command; - if (strcmp(command,"@DEFAULT") == 0 ) { - // return, and flag for @DEFAULT node - return 0; - } - else if (strcmp(command,"EXIT") == 0 ) { - // return, and flag for global exit - return 1; - } - + if (strcmp(command,"@DEFAULT") == 0 || strcmp(command,"EXIT") == 0) + return command; } + return NULL; } +/* ---------------------------------------------------------------------- + perform minimization under control of driver +---------------------------------------------------------------------- */ - -int CommandMDIEngine::mdi_optg() +char *MDIEngine::mdi_optg() { - char *command = NULL; + // initialize an energy minization - // create instance of the Minimizer class Minimize *minimizer = new Minimize(lmp); - // initialize the minimizer in a way that ensures optimization will continue until MDI exits + // setup the minimizer in a way that ensures optimization + // will continue until MDI driver exits + update->etol = std::numeric_limits::min(); update->ftol = std::numeric_limits::min(); update->nsteps = std::numeric_limits::max(); update->max_eval = std::numeric_limits::max(); - update->whichflag = 2; // 2 for minimization + update->whichflag = 2; update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + update->nsteps; + lmp->init(); + // engine is now at @INIT_OPTG node + + char *command = NULL; command = mdi_fix->engine_mode("@INIT_OPTG"); - if (strcmp(command,"@DEFAULT") == 0 ) { - // return, and flag for @DEFAULT node - return 0; - } - else if (strcmp(command,"EXIT") == 0 ) { - // return, and flag for global exit - return 1; - } + if (strcmp(command,"@DEFAULT") == 0 || strcmp(command,"EXIT") == 0) + return command; + + // setup the minimization update->minimize->setup(); + + // get new command + command = mdi_fix->command; - if (strcmp(command,"@DEFAULT") == 0 ) { - // return, and flag for @DEFAULT node - return 0; - } - else if (strcmp(command,"EXIT") == 0 ) { - // return, and flag for global exit - return 1; - } + if (strcmp(command,"@DEFAULT") == 0 || strcmp(command,"EXIT") == 0) + return command; + + // NOTE: perform minimization forever, for huge # of steps ?? + // if you are expecting driver to terminate the minimization, + // not sure how control will return to this function ?? update->minimize->iterate(update->nsteps); + + // return if driver sends @DEFAULT or EXIT + command = mdi_fix->command; - if (strcmp(command,"@DEFAULT") == 0 ) { - // return, and flag for @DEFAULT node - return 0; - } - else if (strcmp(command,"EXIT") == 0 ) { - // return, and flag for global exit - return 1; - } + if (strcmp(command,"@DEFAULT") == 0 || strcmp(command,"EXIT") == 0) + return command; - error->all(FLERR,fmt::format("MDI reached end of OPTG simulation with invalid command: {}",command)); - return 0; + error->all(FLERR, + fmt::format("MDI reached end of OPTG simulation " + "with invalid command: {}",command)); + return NULL; } diff --git a/src/USER-MDI/mdi_engine.h b/src/USER-MDI/mdi_engine.h index 8e74e05e6a..fb0216b667 100644 --- a/src/USER-MDI/mdi_engine.h +++ b/src/USER-MDI/mdi_engine.h @@ -13,27 +13,28 @@ #ifdef COMMAND_CLASS -CommandStyle(mdi_engine,CommandMDIEngine) +CommandStyle(mdi_engine,MDIEngine) #else -#ifndef LMP_COMMAND_MDI_ENGINE_H -#define LMP_COMMAND_MDI_ENGINE_H +#ifndef LMP_MDI_ENGINE_H +#define LMP_MDI_ENGINE_H #include "command.h" namespace LAMMPS_NS { - class CommandMDIEngine : public Command { +class MDIEngine : public Command { public: - CommandMDIEngine(class LAMMPS *); - virtual ~CommandMDIEngine(); + MDIEngine(LAMMPS *lmp) : Command(lmp) {} + virtual ~MDIEngine() {} void command(int, char **); - int mdi_md(); - int mdi_optg(); private: class FixMDIEngine *mdi_fix; + + char *mdi_md(); + char *mdi_optg(); }; }