diff --git a/examples/mdi/README b/examples/mdi/README index 81999ae9f7..0b67d50217 100644 --- a/examples/mdi/README +++ b/examples/mdi/README @@ -21,8 +21,8 @@ The example run commands below have variants for these options. ------------------------------------------------- ------------------------------------------------- -* Example #1 = AIMD with LAMMPS as both a driver and engine - As an engine, LAMMPS is acting as a surrogate for a quantum code. +* Example #1 = run AIMD with 2 instances of LAMMPS as driver and engine + as an engine, LAMMPS is a surrogate for a quantum code Note that the 2 input scripts in.aimd.alone and in.aimd.driver have an option for running in NVE vs NPT mode. Comment in/out @@ -32,10 +32,13 @@ changed in the 3rd input script in.aimd.engine. --- Run the entire calculation with a single instance of LAMMPS by itself - results should be identical to running in driver/engine mode + results should be identical to running this example with MDI % lmp_mpi < in.aimd.alone +With MDI, the thermo output of the driver should match the thermo +output of the in.aimd.alone script. + --- Run with TCP: 1 proc each @@ -67,14 +70,14 @@ Run with MPI: 3 procs + 4 procs ------------------------------------------------- ------------------------------------------------- -* Example #2 = Use a Python driver code to run a sequence of independent - LAMMPS calculations +* Example #2 = Python driver runs sequence of unrelated LAMMPS calculations + calcs can be single-point, MD runs, or minimizations -Note that the sequence_driver.py code allows for optional switches in -addition to -mdi (required) and the -plugin and -plugin_args swithces -which are used to link to an engine as a plugin library. The example -run commands below just usu the default values fo rhte optional -switches. They are also explained the top of the file; the info is +The sequence_driver.py code allows for optional switches in addition +to -mdi (required) and the -plugin and -plugin_args switches which are +used to link to an engine as a plugin library. The example run +commands below just use the default values of the optional switches. +The switches are also explained the top of the file; the info is copied here: # -n 10 @@ -139,3 +142,62 @@ mpiexec -n 1 python3 plugin_driver.py --plugin_name "lammps" --mdi "-role DRIVER Run in plugin mode: 3 procs % mpirun -np 3 python3 sequence_driver.py -plugin lammps -mdi "-role DRIVER -name sequence -method LINK -plugin_path /home/sjplimp/lammps/git/src" -plugin_args "-log log.sequence -in in.sequence" + +------------------------------------------------- +------------------------------------------------- + +* Example #3 = run AIMD with Python driver code and 2 LAMMPS instances + first LAMMPS instance is MM = performs MD timesteps + second LAMMPS instance is surrogate QM = computes forces + +The aimd_driver.py code allows for an optional switch in addition to +-mdi (required) and the -plugin and -plugin_args swiches which are +used to link to the 2 engines as a plugin libraries. The example run +commands below use the default values of the optional switch. The +switches are also explained the top of the file; the info is copied +here: + +# -nsteps 5 +# number of timesteps in dynamics runs, default = 5 + +--- + +Run the entire calculation with a single instance of LAMMPS by itself + results should be identical to running this example with MDI + +% lmp_mpi < in.aimd.alone + +With MDI, the driver prints the QM and Total energies. These should +match the PotEng and TotEng output of the in.aimd.alone script. + +--- + +Run with TCP: 1 proc each + +% python3 aimd_driver.py -mdi "-role DRIVER -name aimd -method TCP -port 8021" + +% lmp_mpi -mdi "-role ENGINE -name MM -method TCP -port 8021 -hostname localhost" -log log.aimd.mm -in in.aimd.mm + +% lmp_mpi -mdi "-role ENGINE -name QM -method TCP -port 8021 -hostname localhost" -log log.aimd.qm -in in.aimd.qm + +--- + +Run with TCP: 2 procs + 2 procs + 3 procs + +% mpirun -np 2 python3 aimd_driver.py -mdi "-role DRIVER -name aimd -method TCP -port 8021" + +% mpirun -np 2 lmp_mpi -mdi "-role ENGINE -name MM -method TCP -port 8021 -hostname localhost" -log log.aimd.mm -in in.aimd.mm + +% mpirun -np 3 lmp_mpi -mdi "-role ENGINE -name QM -method TCP -port 8021 -hostname localhost" -log log.aimd.qm -in in.aimd.qm + +--- + +Run with MPI: 1 proc each + +% mpirun -np 1 python3 aimd_driver.py -mdi "-role DRIVER -name aimd -method MPI" : -np 1 lmp_mpi -mdi "-role ENGINE -name MM -method MPI" -log log.aimd.mm -in in.aimd.mm : -np 1 lmp_mpi -mdi "-role ENGINE -name QM -method MPI" -log log.aimd.qm -in in.aimd.qm + +--- + +Run with MPI: 2 procs + 2 procs + 3 procs + +% mpirun -np 2 python3 aimd_driver.py -mdi "-role DRIVER -name aimd -method MPI" : -np 2 lmp_mpi -mdi "-role ENGINE -name MM -method MPI" -log log.aimd.mm -in in.aimd.mm : -np 3 lmp_mpi -mdi "-role ENGINE -name QM -method MPI" -log log.aimd.qm -in in.aimd.qm diff --git a/examples/mdi/aimd_driver.py b/examples/mdi/aimd_driver.py new file mode 100644 index 0000000000..2d8fe10c1a --- /dev/null +++ b/examples/mdi/aimd_driver.py @@ -0,0 +1,255 @@ +# MDI driver to perform an AIMD simulation +# using one instance of LAMMPS as the MD timestepper +# using second instance of LAMMPS as a QM surrogate to compute forces + +# NOTE: this script is derived from the MDI_AIMD_Driver.cpp code +# included in the MDI distribution +# it alters the timestepping to match a velocity Verlet algorithm +# forces are computed once before timestepping beings +# both the @COORDS and @FORCES nodes are triggered in the MM code +# as the appropriate places to extract COORDS and provide FORCES + +# Syntax: python3 aimd_driver.py switch arg switch arg ... +# possible switches: +# -mdi "-role DRIVER ..." +# required switch +# example for stand-alone mode: +# -mdi "-role DRIVER -name sequence -method TCP -port 8021" +# example for plugin mode: +# -mdi "-role DRIVER -name sequemce -method LINK +# -plugin_path /home/sjplimp/lammps/src/" +# -plugin name +# name of plugin library, only when using plugin mode +# -plugin_args arglist +# args to add when launching plugin library, only when using plugin mode +# enclose arglist in quotes if multiple words +# -nsteps 5 +# number of timesteps, default = 5 + +import sys,math,random +import mdi +import numpy as np +from mpi4py import MPI + +# error message + +def error(txt=None): + if txt: raise Exception(txt) + raise Exception("Syntax: python3 aimd_driver.py switch arg switch arg ...") + +# run an AIMD simulation + +def perform_aimd(world,mm_comm,qm_comm): + + me = world.Get_rank() + nprocs = world.Get_size() + + # receive number of atoms from the MM engine + + mdi.MDI_Send_command("COORDS",qm_comm) + mdi.MDI_Send(coords,3*natoms,mdi.MDI_DOUBLE,qm_comm) + + # get QM potential energy + + mdi.MDI_Send_command("FORCES",mm_comm) + mdi.MDI_Send(forces,3*natoms,mdi.MDI_DOUBLE,mm_comm) + + # get MM kinetic energy + + mdi.MDI_Send_command("COORDS",qm_comm) + mdi.MDI_Send(coords,3*natoms,mdi.MDI_DOUBLE,qm_comm) + + # get QM potential energy + + mdi.MDI_Send_command("FORCES",mm_comm) + mdi.MDI_Send(forces,3*natoms,mdi.MDI_DOUBLE,mm_comm) + + # MM engine proceeds to @ENDSTEP node + # so that KE will be for fully updated velocity + + mdi.MDI_Send_command("@ENDSTEP",mm_comm) + + # get MM kinetic energy + + mdi.MDI_Send_command(" narg: error() + mdiarg = args[iarg+1] + iarg += 2 + elif args[iarg] == "-plugin": + if iarg+2 > narg: error() + plugin = args[iarg+1] + iarg += 2 + elif args[iarg] == "-plugin_args": + if iarg+2 > narg: error() + plugin_args = args[iarg+1] + iarg += 2 + elif args[iarg] == "-nsteps": + if iarg+2 > narg: error() + nsteps = int(args[iarg+1]) + if nsteps < 0: error() + iarg += 2 + else: error() + +if not mdiarg: error() + +# LAMMPS engines are stand-alone codes +# world = MPI communicator for just this driver +# invoke perform_tasks() directly + +if not plugin: + mdi.MDI_Init(mdiarg) + world = mdi.MDI_MPI_get_world_comm() + + # connect to 2 engines, determine which is MM vs QM + + mdicomm1 = mdi.MDI_Accept_Communicator() + mdicomm2 = mdi.MDI_Accept_Communicator() + + mdi.MDI_Send_command("CELL",mdicomm) mdi.MDI_Send(vec,9,mdi.MDI_DOUBLE,mdicomm) - print("POST-CELL") # create atoms on perfect lattice @@ -147,16 +143,13 @@ def perform_tasks(world,mdicomm,dummy): elif mode == "min": mdi.MDI_Send_command("@INIT_OPTG",mdicomm) mdi.MDI_Send_command(">TOLERANCE",mdicomm) - params = [1.0e-4,1.0e-4,100.0,100.0] + params = [tol,tol,1000.0,1000.0] mdi.MDI_Send(params,4,mdi.MDI_DOUBLE,mdicomm) mdi.MDI_Send_command("@DEFAULT",mdicomm) # request potential energy - print("PRE-PE") - mdi.MDI_Send_command("all(FLERR,"MDI: compute[icompute_pe]; press = modify->compute[icompute_press]; + //pe = modify->get_compute_by_id("thermo_pe"); + //press = modify->get_compute_by_id("thermo_press"); + // irregular class used if >COORDS change dramatically irregular = new Irregular(lmp); @@ -167,8 +169,6 @@ void MDIEngine::mdi_engine(int narg, char **arg) MDI_Accept_communicator(&mdicomm); if (mdicomm <= 0) error->all(FLERR,"Unable to connect to MDI driver"); - printf("ENG post accept MDI comm\n"); - // endless engine loop, responding to driver commands mode = DEFAULT; @@ -178,7 +178,7 @@ void MDIEngine::mdi_engine(int narg, char **arg) while (1) { // top-level mdi engine only recognizes three nodes - // DEFAULT, INIT_MD, INIT_OPTG, INIT_SYS, EVAL + // DEFAULT, INIT_MD, INIT_OPTG engine_node("@DEFAULT"); @@ -207,7 +207,6 @@ void MDIEngine::mdi_engine(int narg, char **arg) delete [] node_engine; delete [] node_driver; - modify->delete_compute(id_ke); modify->delete_compute(id_pe); modify->delete_compute(id_press); @@ -232,8 +231,6 @@ void MDIEngine::engine_node(const char *node) { int ierr; - printf("ENG ENODE %s\n",node); - // do not process commands if engine and driver request are not the same strncpy(node_engine,node,MDI_COMMAND_LENGTH); @@ -248,13 +245,9 @@ void MDIEngine::engine_node(const char *node) // read the next command from the driver // all procs call this, but only proc 0 receives the command - printf("ENG PRE-RECV %d\n",mdicomm); - ierr = MDI_Recv_command(mdicmd,mdicomm); if (ierr) error->all(FLERR,"MDI: Unable to receive command from driver"); - printf("ENG POST-RECV %s\n",mdicmd); - // broadcast command to the other MPI tasks MPI_Bcast(mdicmd,MDI_COMMAND_LENGTH,MPI_CHAR,0,world); @@ -299,7 +292,8 @@ int MDIEngine::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 command {} unsupported by engine node {}", + command,node_engine); // --------------------------------------- // respond to MDI standard commands @@ -318,6 +312,9 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) } else if (strcmp(command,">COORDS") == 0) { receive_coords(); + } else if (strcmp(command,">FORCES") == 0) { + receive_double3(FORCE); + } else if (strcmp(command,">NATOMS") == 0) { receive_natoms(); @@ -346,11 +343,11 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) send_double3(COORD); } else if (strcmp(command,"etol = std::numeric_limits::max(); - update->ftol = std::numeric_limits::max(); - update->max_eval = 0; - } - } else if (strcmp(command,"@COORDS") == 0) { strncpy(node_driver,command,MDI_COMMAND_LENGTH); node_match = false; @@ -426,14 +415,6 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) } else if (strcmp(command,"EXIT") == 0) { exit_command = true; - // if minimization in progress, force it to quit - - if (mode == OPT) { - update->etol = std::numeric_limits::max(); - update->ftol = std::numeric_limits::max(); - update->max_eval = 0; - } - // ------------------------------------------------------- // custom LAMMPS commands // ------------------------------------------------------- @@ -458,7 +439,7 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) // ------------------------------------------------------- } else { - error->all(FLERR,"MDI: Unknown command received from driver"); + error->all(FLERR,"MDI: Unknown command {} received from driver",command); } return 0; @@ -506,22 +487,13 @@ void MDIEngine::mdi_commands() MDI_Register_command("@DEFAULT", "COMMAND"); MDI_Register_command("@DEFAULT", "COMMANDS"); MDI_Register_command("@DEFAULT", "INFILE"); - MDI_Register_command("@DEFAULT", ">TOLERANCE"); MDI_Register_command("@DEFAULT", "CELL"); - MDI_Register_command("@INIT_MD", ">CELL_DISPL"); - MDI_Register_command("@INIT_MD", ">CHARGES"); - MDI_Register_command("@INIT_MD", ">COORDS"); - MDI_Register_command("@INIT_MD", ">NATOMS"); MDI_Register_command("@INIT_MD", ">NITERATE"); - MDI_Register_command("@INIT_MD", ">TYPES"); - MDI_Register_command("@INIT_MD", ">VELOCITIES"); MDI_Register_command("@INIT_MD", "@"); MDI_Register_command("@INIT_MD", "@DEFAULT"); MDI_Register_command("@INIT_MD", "@COORDS"); @@ -533,24 +505,15 @@ void MDIEngine::mdi_commands() MDI_Register_node("@INIT_OPTG"); MDI_Register_command("@INIT_OPTG", "<@"); - MDI_Register_command("@INIT_OPTG", "CELL"); - MDI_Register_command("@INIT_OPTG", ">CELL_DISPL"); - MDI_Register_command("@INIT_OPTG", ">CHARGES"); - MDI_Register_command("@INIT_OPTG", ">COORDS"); - MDI_Register_command("@INIT_OPTG", ">NATOMS"); MDI_Register_command("@INIT_OPTG", ">TOLERANCE"); - MDI_Register_command("@INIT_OPTG", ">TYPES"); - MDI_Register_command("@INIT_OPTG", ">VELOCITIES"); MDI_Register_command("@INIT_OPTG", "@"); MDI_Register_command("@INIT_OPTG", "@DEFAULT"); MDI_Register_command("@INIT_OPTG", "@COORDS"); MDI_Register_command("@INIT_OPTG", "@FORCES"); - MDI_Register_command("@INIT_OPTG", "@ENDSTEP"); MDI_Register_command("@INIT_OPTG", "EXIT"); // node at POST_INTEGRATE location in timestep - // only if fix MDI/ENGINE is instantiated + // only used if fix MDI/ENGINE is instantiated MDI_Register_node("@COORDS"); MDI_Register_command("@COORDS", "<@"); @@ -564,17 +527,19 @@ void MDIEngine::mdi_commands() MDI_Register_command("@COORDS", "EXIT"); // node at POST_FORCE location in timestep - // only if fix MDI/ENGINE is instantiated + // only used if fix MDI/ENGINE is instantiated MDI_Register_node("@FORCES"); - MDI_Register_command("@FORCES", "<@"); - MDI_Register_command("@FORCES", "FORCES"); MDI_Register_callback("@FORCES", ">+FORCES"); + MDI_Register_command("@FORCES", "<@"); + MDI_Register_command("@FORCES", "FORCES"); MDI_Register_command("@FORCES", "@"); MDI_Register_command("@FORCES", "@DEFAULT"); MDI_Register_command("@FORCES", "@COORDS"); @@ -583,14 +548,15 @@ void MDIEngine::mdi_commands() MDI_Register_command("@FORCES", "EXIT"); // node at END_OF_STEP location in timestep - // only if fix MDI/ENGINE is instantiated + // only used if fix MDI/ENGINE is instantiated MDI_Register_node("@ENDSTEP"); MDI_Register_command("@ENDSTEP", "<@"); - MDI_Register_command("@FORCES", "NITERATE command sets # of timesteps + either for NITERATE steps or one step at a time + latter is controlled by driver ---------------------------------------------------------------------- */ void MDIEngine::mdi_md() { - // engine is now at @INIT_MD node - // receive >NITERATE or other commands driver wishes to send - // @RUN_MD or @DEFAULT command from driver will trigger the simulation - - niterate = 0; - - engine_node("@INIT_MD"); - - if (strcmp(mdicmd,"EXIT") == 0) return; - - // create or update system if requested - // assume only incremental changes in atom coords + // create or update system if requested prior to @INIT_MD int flag_create = flag_natoms | flag_types; if (flag_create) create_system(); @@ -628,40 +584,15 @@ void MDIEngine::mdi_md() if (flag_velocities) adjust_velocities(); } - // perform the MD simulation + // engine is now at @INIT_MD node + // receive >NITERATE command if driver sends, else niterate = -1 + // any @ command from driver will start the simulation - update->whichflag = 1; - timer->init_timeout(); + niterate = -1; - update->nsteps = niterate; - update->beginstep = update->firststep = update->ntimestep; - update->endstep = update->laststep = update->ntimestep + update->nsteps; - - lmp->init(); - update->integrate->setup(1); - - timer->init(); - timer->barrier_start(); - update->integrate->run(niterate); - timer->barrier_stop(); - - update->integrate->cleanup(); + engine_node("@INIT_MD"); + if (strcmp(mdicmd,"EXIT") == 0) return; - engine_node("@RUN_MD"); - - // clear flags - - flag_natoms = flag_types = 0; - flag_cell = flag_cell_displ = flag_charges = flag_coords = flag_velocities = 0; -} - -/* ---------------------------------------------------------------------- - run MD simulation under control of driver one step at a time - use of fix MDI/ENGINE allows MDI comm within timesteps ----------------------------------------------------------------------- */ - -void MDIEngine::mdi_md_old() -{ // add an instance of fix MDI/ENGINE // delete the instance before this method returns @@ -670,77 +601,65 @@ void MDIEngine::mdi_md_old() (FixMDIEngine *) modify->get_fix_by_id("MDI_ENGINE_INTERNAL"); mdi_fix->mdi_engine = this; - // initialize a new MD simulation + // initialize LAMMPS and setup() the simulation + // set nsteps to niterate if >= 0, else set to 1 update->whichflag = 1; timer->init_timeout(); - update->nsteps = 1; - update->ntimestep = 0; - update->firststep = update->ntimestep; - update->laststep = update->ntimestep + update->nsteps; - update->beginstep = update->firststep; - update->endstep = update->laststep; + + update->nsteps = (niterate >= 0) ? niterate : 1; + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->ntimestep + update->nsteps; lmp->init(); - - // engine is now at @INIT_MD node - // receive any commands driver wishes to send - - engine_node("@INIT_MD"); - - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) { - modify->delete_fix("MDI_ENGINE_INTERNAL"); - return; - } - - // setup the MD simulation - update->integrate->setup(1); - // run MD one step at a time until driver sends @DEFAULT or EXIT - // driver can communicate with LAMMPS within each timestep - // by sending a node command which matches a method in FixMDIEngine + timer->init(); + timer->barrier_start(); - while (true) { - update->whichflag = 1; - timer->init_timeout(); - update->nsteps += 1; - update->laststep += 1; - update->endstep = update->laststep; - output->next = update->ntimestep + 1; + // if niterate >= 0, run for niterate steps + // else if niterate < 0: + // run one step at a time forever + // driver triggers exit with @ command other than @COORDS,@FORCES,@ENDSTEP - // single MD timestep + if (niterate >= 0) { + update->integrate->run(niterate); - update->integrate->run(1); + } else { - // driver triggers end of MD loop by sending @DEFAULT or EXIT + while (true) { + update->nsteps += 1; + update->laststep += 1; + update->endstep = update->laststep; + output->next = update->ntimestep + 1; - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) { - modify->delete_fix("MDI_ENGINE_INTERNAL"); - return; + update->integrate->run(1); + + if (strcmp(mdicmd,"@COORDS") != 0 && + strcmp(mdicmd,"@FORCES") != 0 && + strcmp(mdicmd,"@ENDSTEP") != 0) break; } } + + timer->barrier_stop(); + update->integrate->cleanup(); + modify->delete_fix("MDI_ENGINE_INTERNAL"); + + // clear flags + + flag_natoms = flag_types = 0; + flag_cell = flag_cell_displ = flag_charges = flag_coords = flag_velocities = 0; } /* ---------------------------------------------------------------------- - perform minimization for at most Niterate steps - >TOLERANCE command sets tolerances + perform minimization + either to convergence using >TOLERANCE settings or one iteration at a time + latter is controlled by driver ---------------------------------------------------------------------- */ void MDIEngine::mdi_optg() { - // engine is now at @INIT_OPTG node - // receive >TOLERANCE or other commands driver wishes to send - // @DEFAULT command from driver will trigger the minimization - - etol = ftol = 1.0e-6; - niterate = max_eval = 1000; - - engine_node("@INIT_OPTG"); - - if (strcmp(mdicmd,"EXIT") == 0) return; - - // create or update system if requested + // create or update system if requested prior to @INIT_OPTG int flag_create = flag_natoms | flag_types; if (flag_create) create_system(); @@ -751,42 +670,16 @@ void MDIEngine::mdi_optg() if (flag_velocities) adjust_velocities(); } - // perform the minmization + // engine is now at @INIT_OPTG node + // receive >TOLERANCE if driver sends - update->etol = etol; - update->ftol = ftol; - update->nsteps = niterate; - update->max_eval = max_eval; + etol = ftol = 1.0e-6; + niterate = -1; + max_eval = std::numeric_limits::max(); - update->whichflag = 2; - timer->init_timeout(); + engine_node("@INIT_OPTG"); + if (strcmp(mdicmd,"EXIT") == 0) return; - update->beginstep = update->firststep = update->ntimestep; - update->endstep = update->laststep = update->firststep + update->nsteps; - - lmp->init(); - update->minimize->setup(); - - timer->init(); - timer->barrier_start(); - update->minimize->run(update->nsteps); - timer->barrier_stop(); - - update->minimize->cleanup(); - - // clear flags - - flag_natoms = flag_types = 0; - flag_cell = flag_cell_displ = flag_charges = flag_coords = flag_velocities = 0; -} - -/* ---------------------------------------------------------------------- - perform minimization under control of driver one iteration at a time - use of fix MDI/ENGINE allows MDI comm within iteration ----------------------------------------------------------------------- */ - -void MDIEngine::mdi_optg_old() -{ // add an instance of fix MDI/ENGINE // delete the instance before this method returns @@ -795,48 +688,60 @@ void MDIEngine::mdi_optg_old() (FixMDIEngine *) modify->get_fix_by_id("MDI_ENGINE_INTERNAL"); mdi_fix->mdi_engine = this; - // set tolerances to epsilon and iteration limits huge - // allows MDI driver to force an exit - - 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(); + // initialize LAMMPS and setup() the simulation + // set nsteps to niterate if >= 0 via >TOLERANCE, else set to huge value update->whichflag = 2; + timer->init_timeout(); + + update->nsteps = (niterate >= 0) ? niterate : max_eval; update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + update->nsteps; + update->etol = etol; + update->ftol = ftol; + update->max_eval = max_eval; + lmp->init(); - - // engine is now at @INIT_OPTG node - // receive any commands driver wishes to send - - engine_node("@INIT_OPTG"); - - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) { - modify->delete_fix("MDI_ENGINE_INTERNAL"); - return; - } - - // setup the minimization - update->minimize->setup(); - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) return; + timer->init(); + timer->barrier_start(); - // start minimization - // when the driver sends @DEFAULT or EXIT minimizer tolerances are - // set to large values to force it to exit + // if niterate >= 0, minimize for at most niterate iterations + // else if niterate < 0: + // run one iteration at a time forever + // driver triggers exit with @ command other than @COORDS,@FORCES + // two issues with running in this mode: + // @COORDS and @FORCES are not triggered per min iteration + // but also for line search evals + // if driver triggers exit on step that is not multiple of thermo output + // then energy/virial not computed, and minimize->iterate(update->nsteps); + if (niterate >= 0) { + update->minimize->run(niterate); - // return if driver sends @DEFAULT or EXIT + } else { + niterate = std::numeric_limits::max(); + update->etol = 0.0; + update->ftol = 0.0; - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) { - modify->delete_fix("MDI_ENGINE_INTERNAL"); - return; + while (1) { + update->minimize->run(1); + + if (strcmp(mdicmd,"@COORDS") != 0 && + strcmp(mdicmd,"@FORCES") != 0) break; + } } + + timer->barrier_stop(); + update->minimize->cleanup(); + modify->delete_fix("MDI_ENGINE_INTERNAL"); + + // clear flags + + flag_natoms = flag_types = 0; + flag_cell = flag_cell_displ = flag_charges = flag_coords = flag_velocities = 0; } /* ---------------------------------------------------------------------- @@ -847,14 +752,14 @@ void MDIEngine::mdi_optg_old() void MDIEngine::evaluate() { - // create or update system if requested - int flag_create = flag_natoms | flag_types; + int flag_other = flag_cell | flag_cell_displ | flag_charges | + flag_coords | flag_velocities; + + // create or update system if requested + if (flag_create) create_system(); - else { - int flag_any = flag_cell | flag_cell_displ | flag_charges | - flag_coords | flag_velocities; - if (!flag_any) return; + else if (flag_other) { if (flag_cell || flag_cell_displ) adjust_box(); if (flag_charges) adjust_charges(); if (flag_coords) adjust_coords(); @@ -875,8 +780,8 @@ void MDIEngine::evaluate() // this can only be done if comm->style is not tiled // also requires atoms be in box and lamda coords (for triclinic) // finally invoke setup_minimal(1) to trigger exchange() & reneigh() - // NOTE: what this logic still lacks is invoking migrate_atoms() - // if necessary for comm->style tiled, not easy to detect + // NOTE: what this logic still lacks for comm->style tiled, + // is when to invoke migrate_atoms() if necessary, not easy to detect if (flag_create || neighbor->ago < 0) { update->whichflag = 1; @@ -884,7 +789,7 @@ void MDIEngine::evaluate() update->integrate->setup(1); update->whichflag = 0; - } else { + } else if (flag_other) { update->ntimestep++; pe->addstep(update->ntimestep); press->addstep(update->ntimestep); @@ -933,7 +838,8 @@ void MDIEngine::create_system() if (flag_cell == 0 || flag_natoms == 0 || flag_types == 0 || flag_coords == 0) error->all(FLERR, - "@INIT_SYS requires >CELL, >NATOMS, >TYPES, >COORDS MDI commands"); + "MDI create_system requires >CELL, >NATOMS, >TYPES, >COORDS " + "MDI commands"); // remove all existing atoms via delete_atoms command @@ -1203,6 +1109,37 @@ void MDIEngine::receive_velocities() sys_coords[i] * mdi2lmp_velocity; } +/* ---------------------------------------------------------------------- + >FORCES command + receive vector of 3 doubles for all atoms + atoms are ordered by atomID, 1 to Natoms +---------------------------------------------------------------------- */ + +void MDIEngine::receive_double3(int which) +{ + int n = 3*atom->natoms; + int ierr = MDI_Recv(buf3,n,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: tag; + int nlocal = atom->nlocal; + + int ilocal; + + if (which == FORCE) { + 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; + } + } +} + // ---------------------------------------------------------------------- // ---------------------------------------------------------------------- // MDI "<" driver commands that request data diff --git a/src/MDI/mdi_engine.h b/src/MDI/mdi_engine.h index 41dfa61bf3..2e0a619d59 100644 --- a/src/MDI/mdi_engine.h +++ b/src/MDI/mdi_engine.h @@ -35,24 +35,20 @@ class MDIEngine : public Command { void engine_node(const char *node); private: - int lmpunits; // REAL or METAL or NATIVE - int root; // 1 for proc 0, otherwise 0 + int lmpunits; // REAL or METAL or NATIVE + int root; // 1 for proc 0, otherwise 0 + + MDI_Comm mdicomm; // MDI communicator // state of MDI engine - int mode; // which mode engine is in (DEFAULT,MD,OPTG,etc) + int mode; // which mode engine is in (DEFAULT,MD,OPTG) char *mdicmd; // current MDI command being processed char *node_engine; // which node engine is at char *node_driver; // which node driver has requested bool node_match; // true if driver and engine node currently match bool exit_command; // true if EXIT command received from driver - MDI_Comm mdicomm; - - char *id_ke,*id_pe,*id_press; - class Compute *ke,*pe,*press; - class Irregular *irregular; - // unit conversion factors double lmp2mdi_length,mdi2lmp_length; @@ -62,7 +58,8 @@ class MDIEngine : public Command { double lmp2mdi_pressure,mdi2lmp_pressure; double lmp2mdi_virial,mdi2lmp_virial; - // system state for MDI + // flags for data received by engine + // not acted on until a request to send