tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

finished MDI engine and test script debugging

Browse Source
This commit is contained in:
Steve Plimpton
2022-03-25 15:07:38 -06:00
parent 6e959b6f43
commit 055fefc542
8 changed files with 564 additions and 288 deletions
Download Patch File Download Diff File Expand all files Collapse all files

82
examples/mdi/README
View File

@ -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

255
examples/mdi/aimd_driver.py Normal file
View File

@ -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("<NATOMS",mm_comm)
natoms = mdi.MDI_Recv(1,mdi.MDI_INT,mm_comm)
natoms = world.bcast(natoms,root=0)
# allocate arrays for coordinates and forces
coords = np.zeros(3*natoms,dtype=np.float64)
forces = np.zeros(3*natoms,dtype=np.float64)
# MM engine initializes a new MD simulation
mdi.MDI_Send_command("@INIT_MD",mm_comm)
# -----------------
# compute initial forces for Verlet timestepping
# and initial energy for output on step 0
# -----------------
# MM engine proceeds to @FORCES node in setup()
mdi.MDI_Send_command("@FORCES",mm_comm)
# get coords from MM engine
mdi.MDI_Send_command("<COORDS",mm_comm)
mdi.MDI_Recv(3*natoms,mdi.MDI_DOUBLE,mm_comm,buf=coords)
world.Bcast(coords,root=0)
# send coords to QM 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("<PE",qm_comm)
qm_pe = mdi.MDI_Recv(1,mdi.MDI_DOUBLE,qm_comm)
qm_pe = world.bcast(qm_pe,root=0)
# get forces from QM engine
mdi.MDI_Send_command("<FORCES",qm_comm)
mdi.MDI_Recv(3*natoms,mdi.MDI_DOUBLE,qm_comm,buf=forces)
world.Bcast(forces,root=0)
# send forces to MM engine
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("<KE",mm_comm)
mm_ke = mdi.MDI_Recv(1,mdi.MDI_DOUBLE,mm_comm)
mm_ke = world.bcast(mm_ke,root=0)
# output by driver
# normalize energies by atom count
if me == 0:
print("Step %d: MM energy %g, QM energy %g, Total energy %g" % \
(0,mm_ke/natoms,qm_pe/natoms,(mm_ke+qm_pe)/natoms))
# -----------------
# timestepping loop
# -----------------
for istep in range(nsteps):
# MM engine proceeds to @FORCES node
mdi.MDI_Send_command("@FORCES",mm_comm)
# get coords from MM engine
mdi.MDI_Send_command("<COORDS",mm_comm)
mdi.MDI_Recv(3*natoms,mdi.MDI_DOUBLE,mm_comm,buf=coords)
world.Bcast(coords,root=0)
# send coords to QM 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("<PE",qm_comm)
qm_pe = mdi.MDI_Recv(1,mdi.MDI_DOUBLE,qm_comm)
qm_pe = world.bcast(qm_pe,root=0)
# get forces from QM engine
mdi.MDI_Send_command("<FORCES",qm_comm)
mdi.MDI_Recv(3*natoms,mdi.MDI_DOUBLE,qm_comm,buf=forces)
world.Bcast(forces,root=0)
# send forces to MM engine
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("<KE",mm_comm)
mm_ke = mdi.MDI_Recv(1,mdi.MDI_DOUBLE,mm_comm)
mm_ke = world.bcast(mm_ke,root=0)
# output by driver
# normalize energies by atom count
if me == 0:
print("Step %d: MM energy %g, QM energy %g, Total energy %g" % \
(istep+1,mm_ke/natoms,qm_pe/natoms,(mm_ke+qm_pe)/natoms))
# send EXIT to each engine
mdi.MDI_Send_command("EXIT",mm_comm)
mdi.MDI_Send_command("EXIT",qm_comm)
# ------------------------
# main program
# ------------------------
args = sys.argv[1:]
narg = len(args)
# defaults for command-line args
mdiarg = ""
plugin = ""
plugin_args = ""
nsteps = 5
# parse command-line args
iarg = 0
while iarg < narg:
if args[iarg] == "-mdi":
if iarg+2 > 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("<NAME",mdicomm1)
name1 = mdi.MDI_Recv(mdi.MDI_NAME_LENGTH,mdi.MDI_CHAR,mdicomm1)
name1 = world.bcast(name1,root=0)
mdi.MDI_Send_command("<NAME",mdicomm2)
name2 = mdi.MDI_Recv(mdi.MDI_NAME_LENGTH,mdi.MDI_CHAR,mdicomm2)
name2 = world.bcast(name2,root=0)
if name1 == "MM" and name2 == "QM":
mm_comm = mdicomm1
qm_comm = mdicomm2
elif name1 == "QM" and name2 == "MM":
mm_comm = mdicomm2
qm_comm = mdicomm1
else: error("Two engines have invalid names")
perform_aimd(world,mm_comm,qm_comm)
# LAMMPS engines are plugin libraries
# launch plugins
# NOTE: need to run driver on 2 or more procs
# partition driver into 2 MPI comms
# launch one plugin on each partition
# each with their own callback function
if plugin:
error("Cannot yet run in plugin mode")
mdi.MDI_Init(mdiarg)
world = MPI.COMM_WORLD
plugin_args += " -mdi \"-role ENGINE -name lammps -method LINK\""
mdi.MDI_Launch_plugin(plugin,plugin_args,world,perform_tasks,None)

23
examples/mdi/in.aimd.mm Normal file
View File

@ -0,0 +1,23 @@
# 3d Lennard-Jones melt - LAMMPS as MM engine for aimd_driver.py
variable x index 5
variable y index 5
variable z index 5
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 $x 0 $y 0 $z
create_box 1 box
create_atoms 1 box
mass 1 1.0
velocity all create 1.44 87287 loop geom
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes
fix 1 all nve
mdi engine

22
examples/mdi/in.aimd.qm Normal file
View File

@ -0,0 +1,22 @@
# 3d Lennard-Jones melt - LAMMPS as surrogate QM engine for aimd_driver.py
variable x index 5
variable y index 5
variable z index 5
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 $x 0 $y 0 $z
create_box 1 box
create_atoms 1 box
mass 1 1.0
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes
mdi engine

19
examples/mdi/sequence_driver.py
View File

@ -21,7 +21,7 @@
# style of calculations: single snapshot evals, dynamics, minimization
# default = eval
# -size Nx Ny Nz
# cubic lattice, default = 2 2 2
# cubic lattice, default = 4 4 4
# -rho 0.75 0.1
# reduced density and random variation thereof, default = 0.75 0.1
# -delta 0.1
@ -52,8 +52,6 @@ def error(txt=None):
def perform_tasks(world,mdicomm,dummy):
print("PT start",world,mdicomm,dummy)
me = world.Get_rank()
nprocs = world.Get_size()
@ -88,10 +86,8 @@ def perform_tasks(world,mdicomm,dummy):
# send simulation box to engine
vec = [xhi-xlo,0.0,0.0] + [0.0,yhi-ylo,0.0] + [0.0,0.0,zhi-zlo]
print("PRE-CELL",mdicomm)
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("<PE",mdicomm)
print("POST-PE")
pe = mdi.MDI_Recv(1,mdi.MDI_DOUBLE,mdicomm)
pe = world.bcast(pe,root=0)
@ -199,7 +192,7 @@ def perform_tasks(world,mdicomm,dummy):
mdi.MDI_Send_command("EXIT",mdicomm)
# return needed for plugin callback mode
# return needed for plugin mode
return 0
@ -218,7 +211,7 @@ plugin_args = ""
ncalc = 1
mode = "eval"
nx = ny = nz = 2
nx = ny = nz = 4
rho = 0.75
rhodelta = 0.1
delta = 0.1
@ -312,11 +305,9 @@ if not plugin:
# launch plugin
# MDI will call back to perform_tasks()
print("PRE PLUGIN");
if plugin:
error("Cannot yet run in plugin mode")
mdi.MDI_Init(mdiarg)
world = MPI.COMM_WORLD
plugin_args += " -mdi \"-role ENGINE -name lammps -method LINK\""
print("PRE LAUNCH",plugin_args)
mdi.MDI_Launch_plugin(plugin,plugin_args,world,perform_tasks,None)

3
src/MDI/fix_mdi_aimd.cpp
View File

@ -207,7 +207,6 @@ void FixMDIAimd::post_force(int vflag)
}
// optionally request potential energy from MDI engine
// divide by nprocs so each proc stores a portion
if (eflag_global) {
ierr = MDI_Send_command("<PE",mdicomm);
@ -215,7 +214,7 @@ void FixMDIAimd::post_force(int vflag)
ierr = MDI_Recv(&engine_energy,1,MDI_DOUBLE,mdicomm);
if (ierr) error->all(FLERR,"MDI: <PE data");
MPI_Bcast(&engine_energy,1,MPI_DOUBLE,0,world);
engine_energy *= mdi2lmp_energy / nprocs;
engine_energy *= mdi2lmp_energy;
}
// optionally request pressure tensor from MDI engine, convert to virial

397
src/MDI/mdi_engine.cpp
View File

@ -46,8 +46,7 @@
using namespace LAMMPS_NS;
enum{NATIVE,REAL,METAL}; // LAMMPS units which MDI supports
enum{DEFAULT,MD,OPT,SYS,EVAL}; // top-level MDI engine mode
enum{EVALPOINT,EVALMD,EVALOPT}; // EVAL mode
enum{DEFAULT,MD,OPT}; // top-level MDI engine modes
// per-atom data which engine commands access
@ -126,6 +125,9 @@ void MDIEngine::mdi_engine(int narg, char **arg)
pe = modify->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,"<ENERGY") == 0) {
evaluate();
if (strcmp(node_engine,"@DEFAULT") == 0) evaluate();
send_total_energy();
} else if (strcmp(command,"<FORCES") == 0) {
evaluate();
if (strcmp(node_engine,"@DEFAULT") == 0) evaluate();
send_double3(FORCE);
} else if (strcmp(command,"<LABELS") == 0) {
@ -363,11 +360,11 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm)
send_natoms();
} else if (strcmp(command,"<PE") == 0) {
evaluate();
if (strcmp(node_engine,"@DEFAULT") == 0) evaluate();
send_pe();
} else if (strcmp(command,"<STRESS") == 0) {
evaluate();
if (strcmp(node_engine,"@DEFAULT") == 0) evaluate();
send_stress();
} else if (strcmp(command,"<TYPES") == 0) {
@ -401,14 +398,6 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm)
strncpy(node_driver,command,MDI_COMMAND_LENGTH);
node_match = false;
// if minimization in progress, force it to quit
if (mode == OPT) {
update->etol = std::numeric_limits<double>::max();
update->ftol = std::numeric_limits<double>::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<double>::max();
update->ftol = std::numeric_limits<double>::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", "<KE");
// node for setting up and running a dynamics simulation
MDI_Register_node("@INIT_MD");
MDI_Register_command("@INIT_MD", "<@");
MDI_Register_command("@INIT_MD", "<NATOMS");
MDI_Register_command("@INIT_MD", ">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", "<NATOMS");
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", "<ENERGY");
MDI_Register_command("@FORCES", "<FORCES");
MDI_Register_command("@FORCES", "<PE");
MDI_Register_command("@FORCES", "<STRESS");
MDI_Register_callback("@FORCES", ">FORCES");
MDI_Register_callback("@FORCES", ">+FORCES");
MDI_Register_command("@FORCES", "<@");
MDI_Register_command("@FORCES", "<COORDS");
MDI_Register_command("@FORCES", "<ENERGY");
MDI_Register_command("@FORCES", "<FORCES");
MDI_Register_command("@FORCES", "<KE");
MDI_Register_command("@FORCES", "<PE");
MDI_Register_command("@FORCES", "<STRESS");
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", "<ENERGY");
MDI_Register_command("@FORCES", "<FORCES");
MDI_Register_command("@FORCES", "<PE");
MDI_Register_command("@FORCES", "<STRESS");
MDI_Register_command("@ENDSTEP", "<ENERGY");
MDI_Register_command("@ENDSTEP", "<FORCES");
MDI_Register_command("@ENDSTEP", "<KE");
MDI_Register_command("@ENDSTEP", "<PE");
MDI_Register_command("@ENDSTEP", "<STRESS");
MDI_Register_command("@ENDSTEP", "@");
MDI_Register_command("@ENDSTEP", "@DEFAULT");
MDI_Register_command("@ENDSTEP", "@COORDS");
@ -601,23 +567,13 @@ void MDIEngine::mdi_commands()
/* ----------------------------------------------------------------------
run MD simulation
>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<int>::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<double>::min();
update->ftol = std::numeric_limits<double>::min();
update->nsteps = std::numeric_limits<int>::max();
update->max_eval = std::numeric_limits<int>::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 <PE, <STRESS will fail
update->minimize->iterate(update->nsteps);
if (niterate >= 0) {
update->minimize->run(niterate);
// return if driver sends @DEFAULT or EXIT
} else {
niterate = std::numeric_limits<int>::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: <double3 data");
MPI_Bcast(buf3,n,MPI_DOUBLE,0,world);
// use atomID to index into ordered buf3
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
int ilocal;
if (which == FORCE) {
double **f = atom->f;
for (int i = 0; i < nlocal; i++) {
ilocal = static_cast<int> (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

51
src/MDI/mdi_engine.h
View File

@ -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 <ENERGY,<FORCES,<PE,<STRESS
int flag_natoms,flag_types;
int flag_cell,flag_cell_displ;
@ -73,11 +70,11 @@ class MDIEngine : public Command {
double *sys_charges,*sys_coords,*sys_velocities;
double sys_cell[9],sys_cell_displ[3];
int niterate;
int niterate; // settings for MD and OPTG
int max_eval;
double etol,ftol;
int nbytes; // NBYTES command value used by other commands
int nbytes; // NBYTES command value used for length by other commands
// buffers for MDI comm
@ -86,15 +83,21 @@ class MDIEngine : public Command {
double *buf3,*buf3all;
int *ibuf1,*ibuf1all;
// other classes used by MDI
char *id_ke,*id_pe,*id_press; // computes invoked by MDI
class Compute *ke,*pe,*press;
class Irregular *irregular; // irregular comm if new COORDS
// are highly displaced
// class methods
void mdi_engine(int, char **);
void mdi_commands();
void mdi_md();
void mdi_md_old();
void mdi_optg();
void mdi_optg_old();
void evaluate();
void create_system();
@ -104,36 +107,20 @@ class MDIEngine : public Command {
void adjust_velocities();
void receive_cell();
void receive_cell_default();
void receive_cell_sys();
void receive_cell_displ();
void receive_cell_displ_default();
void receive_cell_displ_sys();
void receive_charges();
void receive_charges_sys();
void receive_coords();
void receive_coords_sys();
void receive_natoms();
void receive_natoms_default();
void receive_natoms_sys();
void receive_types();
void receive_types_sys();
void receive_velocities();
void receive_velocities_sys();
void receive_double3(int);
void send_cell();
void send_cell_displ();
void send_total_energy();
void send_labels();
void send_natoms();
void send_pe();
void send_stress();