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

updates to fix_mdi_qm for multiple sims

Browse Source
This commit is contained in:
Steve Plimpton
2022-06-14 10:14:36 -06:00
parent 407e015c80
commit 34863c6c97
5 changed files with 297 additions and 171 deletions
Download Patch File Download Diff File Expand all files Collapse all files

68
doc/src/fix_mdi_qm.rst
View File

@ -13,10 +13,13 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* mdi/qm = style name of this fix command
* zero or more keyword/value pairs may be appended
* keyword = *add* or *every*
* keyword = *virial* or *add* or *every*
.. parsed-literal::
*virial* args = *yes* or *no*
yes = request virial tensor from server code
no = do not request virial tensor from server code
*add* args = *yes* or *no*
yes = add returned value from server code to LAMMPS quantities
no = do not add returned values to LAMMPS quantities
@ -29,6 +32,7 @@ Examples
.. code-block:: LAMMPS
fix 1 all mdi/qm
fix 1 all mdi/qm virial yes
fix 1 all mdi/qm add no every 100
Description
@ -79,13 +83,17 @@ explains how to launch the two codes in either mode.
----------
The *virial* keyword setting of yes or no determines whether
LAMMPS will request the QM code to also compute and return
a 6-element symmetric virial tensor for the system.
The *add* keyword setting of *yes* or *no* determines whether the
energy and forces returned by the QM code will be added to the LAMMPS
internal energy and forces or not. If the setting is *no* then the
default :doc:`fix_modify energy <fix_modify>` and :doc:`fix_modify
virial <fix_modify>` settings are also set to *no* and your input
scripts should not set them to yes. See more details on these
fix_modify settings below.
energy and forces and virial returned by the QM code will be added to
the LAMMPS internal energy and forces and virial or not. If the
setting is *no* then the default :doc:`fix_modify energy <fix_modify>`
and :doc:`fix_modify virial <fix_modify>` settings are also set to
*no* and your input scripts should not set them to yes. See more
details on these fix_modify settings below.
Whatever the setting for the *add* keyword, the QM energy, forces, and
virial will be stored by the fix, so they can be accessed by other
@ -103,30 +111,36 @@ configuration of atoms. The QM code will be called once every
(1) To run an ab initio MD (AIMD) dynamics simulation, or an energy
minimization with QM forces, or a multi-replica NEB calculation, use
*add yes* and *every 1* (the defaults). This is so that every time
LAMMPS needs energy and forces, the QM code will be invoked.
LAMMPS needs energy and forces, the QM code will be invoked.
Both LAMMPS and the QM code should define the same system (simulation
box, atoms and their types) in their respective input scripts. Note
that on this scenario, it may not be necessary for LAMMPS to define a
pair style or use a neighbor list.
LAMMPS will then perform the timestepping for the simulation. At the
point in each timestep when LAMMPS needs the force on each atom, it
communicates with the engine code. It sends the current simulation
box size and shape (if they change dynamically, e.g. during an NPT
simulation), and the current atom coordinates. The engine code
computes quantum forces on each atom and returns them to LAMMPS. If
LAMMPS also needs the system energy and/or virial, it requests those
values from the engine code as well.
LAMMPS will then perform the timestepping or minimization iterations
for the simulation. At the point in each timestep or iteration when
LAMMPS needs the force on each atom, it communicates with the engine
code. It sends the current simulation box size and shape (if they
change dynamically, e.g. during an NPT simulation), and the current
atom coordinates. The engine code computes quantum forces on each
atom and the total energy of the system and returns them to LAMMPS.
Note that if the AIMD simulation is an NPT or NPH model, or the energy
minimization includesf :doc:`fix box relax <fix_box_relaxq>` to
equilibrate the box size/shape, then LAMMPS computes a pressure. This
means the *virial* keyword should be set to *yes* so that the QM
contribution to the pressure can be included.
(2) To run dynamics with a LAMMPS interatomic potential, and evaluate
the QM energy and forces once every 1000 steps, use *add no* and
*every 1000*. This could be useful for using an MD run to generate
randomized configurations which then generate QM data for training a
machine learning potential. A :doc:`dump custom <dump>` command could
be invoked every 1000 steps to dump the atom coordinates and QM forces
to a file. Likewise the QM energy could be output with the
:doc:`thermo_style custom <thermo_style>` command.
randomized configurations which are then passed to the QM code to
produce training data for a machine learning potential. A :doc:`dump
custom <dump>` command could be invoked every 1000 steps to dump the
atom coordinates and QM forces to a file. Likewise the QM energy and
virial could be output with the :doc:`thermo_style custom
<thermo_style>` command.
(3) To do a QM evaulation of energy and forces for a series of *N*
independent systems (simulation box and atoms), use *add no* and
@ -138,10 +152,10 @@ could be initialized by reading them from data files with
:doc:`lattice <lattice>` , :doc:`create_atoms <create_atoms>`,
:doc:`delete_atoms <deletea_atoms>`, and/or :doc:`displace_atoms
random <displace_atoms>` commands to generate a series of different
systems. At the end of the loop perform :doc:`run 0 <run>`
and :doc:`write_dump <write_dump>` commands to invoke the QM code
and output the QM energy and forces. As in (2) this be useful
to generate QM data for training a machine learning potential.
systems. At the end of the loop perform :doc:`run 0 <run>` and
:doc:`write_dump <write_dump>` commands to invoke the QM code and
output the QM energy and forces. As in (2) this be useful to produce
QM data for training a machine learning potential.
----------
@ -228,5 +242,5 @@ Related commands
Default
"""""""
The default for the optional keywords is add = yes, every = 1.
The default for the optional keywords are virial = no, add = yes,
every = 1.

341
src/MDI/fix_mdi_qm.cpp
View File

@ -29,13 +29,6 @@ enum { NATIVE, REAL, METAL }; // LAMMPS units which MDI supports
FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
energy_global_flag = 1;
virial_global_flag = 1;
thermo_energy = thermo_virial = 1;
// check requirements for LAMMPS to work with MDI as an engine
if (atom->tag_enable == 0)
@ -52,12 +45,19 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// optional args
virialflag = 0;
addflag = 1;
every = 1;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"add") == 0) {
if (strcmp(arg[iarg],"virial") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix mdi/qm command");
if (strcmp(arg[iarg],"yes") == 0) virialflag = 1;
else if (strcmp(arg[iarg],"no") == 0) virialflag = 0;
else error->all(FLERR,"Illegal fix mdi/qm command");
iarg += 2;
} else if (strcmp(arg[iarg],"add") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix mdi/qm command");
if (strcmp(arg[iarg],"yes") == 0) addflag = 1;
else if (strcmp(arg[iarg],"no") == 0) addflag = 0;
@ -66,18 +66,40 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
} else if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix mdi/qm command");
every = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
if (every < 0) error->all(FLERR,"Illegal fix mdi/qm command");
if (every <= 0) error->all(FLERR,"Illegal fix mdi/qm command");
iarg += 2;
} else error->all(FLERR,"Illegal fix mdi/qm command");
}
// fix output settings are based on optional keywords
scalar_flag = 1;
global_freq = every;
extscalar = 1;
if (virialflag) {
vector_flag = 1;
size_vector = 6;
extvector = 1;
}
if (addflag) {
energy_global_flag = 1;
virial_global_flag = 1;
thermo_energy = thermo_virial = 1;
}
// mdicomm will be one-time initialized in init()
// cannot be done here for a plugin library, b/c mdi plugin command is later
mdicomm = MDI_COMM_NULL;
// storage for all atoms
// peratom storage, both for nlocal and global natoms
fqm = nullptr;
maxlocal = 0;
ibuf1 = ibuf1all = nullptr;
buf3 = buf3all = nullptr;
maxbuf = 0;
@ -93,6 +115,17 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
unit_conversions();
nprocs = comm->nprocs;
// initialize outputs
qm_energy = 0.0;
if (virialflag) {
for (int i = 0; i < 6; i++) {
qm_virial[i] = 0.0;
virial[i] = 0.0;
}
sumflag = 0;
}
}
/* ---------------------------------------------------------------------- */
@ -109,6 +142,10 @@ FixMDIQM::~FixMDIQM()
// clean up
memory->destroy(fqm);
memory->destroy(ibuf1);
memory->destroy(ibuf1all);
memory->destroy(buf3);
memory->destroy(buf3all);
}
@ -128,26 +165,38 @@ int FixMDIQM::setmask()
void FixMDIQM::init()
{
if (mdicomm != MDI_COMM_NULL) return;
// one-time auto-detect whether engine is stand-alone code or plugin library
// also initializes mdicomm
// plugin = 0/1 for engine = stand-alone code vs plugin library
MDI_Get_communicator(&mdicomm, 0);
// set plugin = 0/1 for engine = stand-alone code vs plugin library
if (mdicomm == MDI_COMM_NULL) {
plugin = 0;
MDI_Accept_communicator(&mdicomm);
if (mdicomm == MDI_COMM_NULL)
error->all(FLERR, "MDI unable to connect to stand-alone engine");
} else {
plugin = 1;
int method;
MDI_Get_method(&method, mdicomm);
if (method != MDI_PLUGIN)
error->all(FLERR, "MDI internal error for plugin engine");
MDI_Get_communicator(&mdicomm, 0);
if (mdicomm == MDI_COMM_NULL) {
plugin = 0;
MDI_Accept_communicator(&mdicomm);
if (mdicomm == MDI_COMM_NULL)
error->all(FLERR, "MDI unable to connect to stand-alone engine");
} else {
plugin = 1;
int method;
MDI_Get_method(&method, mdicomm);
if (method != MDI_PLUGIN)
error->all(FLERR, "MDI internal error for plugin engine");
}
}
// send natoms, atom types, and simulation box to engine
// this will trigger setup of a new system
// subsequent calls in post_force() will be for same system until new init()
int ierr = MDI_Send_command(">NATOMS", mdicomm);
if (ierr) error->all(FLERR, "MDI: >NATOMS command");
int n = static_cast<int> (atom->natoms);
ierr = MDI_Send(&n, 1, MDI_INT, mdicomm);
if (ierr) error->all(FLERR, "MDI: >NATOMS data");
send_types();
send_box();
}
/* ---------------------------------------------------------------------- */
@ -159,59 +208,25 @@ void FixMDIQM::setup(int vflag)
/* ---------------------------------------------------------------------- */
void FixMDIQM::setup_pre_reverse(int eflag, int vflag)
{
pre_reverse(eflag, vflag);
}
/* ----------------------------------------------------------------------
store eflag, so can use it in post_force to request energy
------------------------------------------------------------------------- */
void FixMDIQM::pre_reverse(int eflag, int /*vflag*/)
{
eflag_caller = eflag;
}
/* ---------------------------------------------------------------------- */
void FixMDIQM::post_force(int vflag)
{
int ilocal, ierr;
double cell[9];
int index, ierr;
int eflag = eflag_caller;
ev_init(eflag, vflag);
// skip if timestep is not a multiple of every
if (update->ntimestep % every) return;
// reallocate peratom storage if necessary, both natoms and nlocal
reallocate();
// if simulation box dynamically changes, send current box to MDI engine
if (domain->box_change_size || domain->box_change_shape) {
ierr = MDI_Send_command(">CELL_DISPL", mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command");
cell[0] = domain->boxlo[0] * lmp2mdi_length;
cell[1] = domain->boxlo[1] * lmp2mdi_length;
cell[2] = domain->boxlo[2] * lmp2mdi_length;
ierr = MDI_Send(cell, 3, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL_DISPL data");
ierr = MDI_Send_command(">CELL", mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL command");
cell[0] = domain->boxhi[0] - domain->boxlo[0];
cell[1] = 0.0;
cell[2] = 0.0;
cell[3] = domain->xy;
cell[4] = domain->boxhi[1] - domain->boxlo[1];
cell[5] = 0.0;
cell[6] = domain->xz;
cell[7] = domain->yz;
cell[8] = domain->boxhi[2] - domain->boxlo[2];
ierr = MDI_Send(cell, 9, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL data");
}
if (domain->box_change_size || domain->box_change_shape)
send_box();
// gather all coords, ordered by atomID
reallocate();
memset(buf3, 0, 3 * atom->natoms * sizeof(double));
double **x = atom->x;
@ -219,13 +234,14 @@ void FixMDIQM::post_force(int vflag)
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
ilocal = static_cast<int>(tag[i]) - 1;
buf3[3 * ilocal + 0] = x[i][0] * lmp2mdi_length;
buf3[3 * ilocal + 1] = x[i][1] * lmp2mdi_length;
buf3[3 * ilocal + 2] = x[i][2] * lmp2mdi_length;
index = static_cast<int>(tag[i]) - 1;
buf3[3 * index + 0] = x[i][0] * lmp2mdi_length;
buf3[3 * index + 1] = x[i][1] * lmp2mdi_length;
buf3[3 * index + 2] = x[i][2] * lmp2mdi_length;
}
MPI_Reduce(buf3, buf3all, 3 * atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world);
int n = static_cast<int> (atom->natoms);
MPI_Reduce(buf3, buf3all, 3 * n, MPI_DOUBLE, MPI_SUM, 0, world);
// send current coords to MDI engine
@ -234,42 +250,54 @@ void FixMDIQM::post_force(int vflag)
ierr = MDI_Send(buf3all, 3 * atom->natoms, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: >COORDS data");
// request potential energy from MDI engine
// this triggers engine to perform QM calculation
// qm_energy = fix output for global QM energy
ierr = MDI_Send_command("<PE", mdicomm);
if (ierr) error->all(FLERR, "MDI: <PE command");
ierr = MDI_Recv(&qm_energy, 1, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: <PE data");
MPI_Bcast(&qm_energy, 1, MPI_DOUBLE, 0, world);
qm_energy *= mdi2lmp_energy;
// request forces from MDI engine
// this triggers engine to evaluate forces,energy,stress for current system
ierr = MDI_Send_command("<FORCES", mdicomm);
if (ierr) error->all(FLERR, "MDI: <FORCES command");
ierr = MDI_Recv(buf3, 3 * atom->natoms, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: <FORCES data");
MPI_Bcast(buf3, 3 * atom->natoms, MPI_DOUBLE, 0, world);
MPI_Bcast(buf3, 3 * n, MPI_DOUBLE, 0, world);
// add forces to owned atoms
// use atomID to index into ordered buf3
double **f = atom->f;
// fqm = fix output for peratom QM forces
// use atomID of local atoms to index into ordered buf3
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;
index = static_cast<int>(tag[i]) - 1;
fqm[i][0] = buf3[3 * index + 0] * mdi2lmp_force;
fqm[i][1] = buf3[3 * index + 1] * mdi2lmp_force;
fqm[i][2] = buf3[3 * index + 2] * mdi2lmp_force;
}
// optionally request potential energy from MDI engine
// optionally add forces to owned atoms
// use atomID of local atoms to index into ordered buf3
if (eflag_global) {
ierr = MDI_Send_command("<PE", mdicomm);
if (ierr) error->all(FLERR, "MDI: <PE command");
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;
if (addflag) {
double **f = atom->f;
for (int i = 0; i < nlocal; i++) {
index = static_cast<int>(tag[i]) - 1;
f[i][0] += buf3[3 * index + 0] * mdi2lmp_force;
f[i][1] += buf3[3 * index + 1] * mdi2lmp_force;
f[i][2] += buf3[3 * index + 2] * mdi2lmp_force;
}
}
// optionally request pressure tensor from MDI engine, convert to virial
// divide by nprocs so each proc stores a portion
// MPI_Allreduce is performed in compute_vector()
// qm_virial = fix output for global QM virial
if (vflag_global) {
if (virialflag) {
double ptensor[6];
ierr = MDI_Send_command("<STRESS", mdicomm);
if (ierr) error->all(FLERR, "MDI: <STRESS command");
@ -280,8 +308,16 @@ void FixMDIQM::post_force(int vflag)
double volume = domain->xprd * domain->yprd * domain->zprd;
for (int i = 0; i < 6; i++) {
ptensor[i] *= mdi2lmp_pressure;
virial[i] = ptensor[i] * volume / force->nktv2p / nprocs;
qm_virial[i] = ptensor[i] * volume / force->nktv2p / nprocs;
}
sumflag = 0;
}
// optionally set fix->virial
if (virialflag && addflag) {
for (int i = 0; i < 6; i++)
virial[i] = qm_virial[i];
}
}
@ -298,27 +334,111 @@ void FixMDIQM::min_post_force(int vflag)
double FixMDIQM::compute_scalar()
{
return engine_energy;
return qm_energy;
}
/* ----------------------------------------------------------------------
reallocate storage for all atoms if necessary
virial from MDI engine
------------------------------------------------------------------------- */
double FixMDIQM::compute_vector(int n)
{
// only sum across procs one time
if (sumflag == 0) {
MPI_Allreduce(qm_virial, qm_virial_all, 6, MPI_DOUBLE, MPI_SUM, world);
sumflag = 1;
}
return qm_virial_all[n];
}
/* ----------------------------------------------------------------------
reallocate storage for local and global and atoms if needed
------------------------------------------------------------------------- */
void FixMDIQM::reallocate()
{
if (atom->natoms <= maxbuf) return;
if (atom->nlocal > maxlocal) {
maxlocal = atom->nmax;
memory->destroy(fqm);
memory->create(fqm, maxlocal, 3, "mdi:fqm");
array_atom = fqm;
}
if (3 * atom->natoms > MAXSMALLINT)
error->all(FLERR, "Natoms too large to use with fix mdi/qm");
if (atom->natoms > maxbuf) {
bigint nsize = atom->natoms * 3;
if (nsize > MAXSMALLINT)
error->all(FLERR, "Natoms too large to use with fix mdi/qm");
maxbuf = static_cast<int> (atom->natoms);
memory->destroy(ibuf1);
memory->destroy(buf3);
memory->destroy(buf3all);
memory->create(ibuf1, maxbuf, "mdi:ibuf1");
memory->create(ibuf1all, maxbuf, "mdi:ibuf1all");
memory->create(buf3, 3 * maxbuf, "mdi:buf3");
memory->create(buf3all, 3 * maxbuf, "mdi:buf3all");
}
}
maxbuf = atom->natoms;
/* ----------------------------------------------------------------------
send numeric atom types to MDI engine
------------------------------------------------------------------------- */
memory->destroy(buf3);
memory->destroy(buf3all);
void FixMDIQM::send_types()
{
memset(ibuf1, 0, atom->natoms * sizeof(int));
memory->create(buf3, 3 * maxbuf, "mdi:buf3");
memory->create(buf3all, 3 * maxbuf, "mdi:buf3all");
// use local atomID to index into ordered ibuf1
tagint *tag = atom->tag;
int *type = atom->type;
int nlocal = atom->nlocal;
int index;
for (int i = 0; i < nlocal; i++) {
index = static_cast<int>(tag[i]) - 1;
ibuf1[index] = type[i];
}
int n = static_cast<int> (atom->natoms);
MPI_Reduce(ibuf1, ibuf1all, n, MPI_INT, MPI_SUM, 0, world);
int ierr = MDI_Send(ibuf1all, n, MDI_INT, mdicomm);
if (ierr) error->all(FLERR, "MDI: >TYPES data");
}
/* ----------------------------------------------------------------------
send simulation box size and shape to MDI engine
------------------------------------------------------------------------- */
void FixMDIQM::send_box()
{
double cell[9];
int ierr = MDI_Send_command(">CELL_DISPL", mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command");
cell[0] = domain->boxlo[0] * lmp2mdi_length;
cell[1] = domain->boxlo[1] * lmp2mdi_length;
cell[2] = domain->boxlo[2] * lmp2mdi_length;
ierr = MDI_Send(cell, 3, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL_DISPL data");
ierr = MDI_Send_command(">CELL", mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL command");
cell[0] = domain->boxhi[0] - domain->boxlo[0];
cell[1] = 0.0;
cell[2] = 0.0;
cell[3] = domain->xy;
cell[4] = domain->boxhi[1] - domain->boxlo[1];
cell[5] = 0.0;
cell[6] = domain->xz;
cell[7] = domain->yz;
cell[8] = domain->boxhi[2] - domain->boxlo[2];
ierr = MDI_Send(cell, 9, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: >CELL data");
}
/* ----------------------------------------------------------------------
@ -384,17 +504,4 @@ void FixMDIQM::unit_conversions()
ev_to_hartree / (angstrom_to_bohr * angstrom_to_bohr * angstrom_to_bohr) / force->nktv2p;
mdi2lmp_pressure = 1.0 / lmp2mdi_pressure;
}
// velocity units = distance/time
mdi2lmp_velocity = 1.0;
lmp2mdi_velocity = 1.0;
if (lmpunits == REAL) {
lmp2mdi_velocity = angstrom_to_bohr / (1.0e-15 * second_to_aut);
mdi2lmp_velocity = 1.0 / lmp2mdi_velocity;
} else if (lmpunits == METAL) {
lmp2mdi_velocity = angstrom_to_bohr / (1.0e-12 * second_to_aut);
mdi2lmp_velocity = 1.0 / lmp2mdi_velocity;
}
}

18
src/MDI/fix_mdi_qm.h
View File

@ -33,22 +33,25 @@ class FixMDIQM : public Fix {
void init();
void setup(int);
void setup_pre_reverse(int, int);
void pre_reverse(int, int);
void post_force(int);
void min_post_force(int);
double compute_scalar();
double compute_vector(int);
private:
int nprocs;
int virialflag,addflag,every;
int plugin;
int maxlocal;
int sumflag;
double qm_energy;
int lmpunits;
double qm_virial[6],qm_virial_all[6];
double **fqm;
MDI_Comm mdicomm;
int eflag_caller;
double engine_energy;
int lmpunits;
// unit conversion factors
double lmp2mdi_length, mdi2lmp_length;
@ -60,11 +63,14 @@ class FixMDIQM : public Fix {
// buffers for MDI comm
int maxbuf;
int *ibuf1, *ibuf1all;
double *buf3, *buf3all;
// methods
void reallocate();
void send_types();
void send_box();
void unit_conversions();
};

40
src/MDI/mdi_engine.cpp
View File

@ -135,7 +135,7 @@ MDIEngine::MDIEngine(LAMMPS *_lmp, int narg, char ** /*arg*/) : Pointers(_lmp)
ibuf1 = ibuf1all = nullptr;
maxatom = 0;
sys_natoms = atom->natoms;
sys_natoms = static_cast<int> (atom->natoms);
reallocate();
nsteps = 0;
@ -1239,7 +1239,7 @@ void MDIEngine::receive_velocities()
void MDIEngine::receive_double3(int which)
{
int n = 3 * atom->natoms;
int n = 3 * sys_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);
@ -1352,10 +1352,10 @@ void MDIEngine::send_total_energy()
void MDIEngine::send_labels()
{
auto labels = new char[atom->natoms * MDI_LABEL_LENGTH];
memset(labels, ' ', atom->natoms * MDI_LABEL_LENGTH);
auto labels = new char[sys_natoms * MDI_LABEL_LENGTH];
memset(labels, ' ', sys_natoms * MDI_LABEL_LENGTH);
memset(ibuf1, 0, atom->natoms * sizeof(int));
memset(ibuf1, 0, sys_natoms * sizeof(int));
// use atomID to index into ordered ibuf1
@ -1370,17 +1370,17 @@ void MDIEngine::send_labels()
ibuf1[ilocal] = type[i];
}
MPI_Reduce(ibuf1, ibuf1all, atom->natoms, MPI_INT, MPI_SUM, 0, world);
MPI_Reduce(ibuf1, ibuf1all, sys_natoms, MPI_INT, MPI_SUM, 0, world);
if (comm->me == 0) {
for (int iatom = 0; iatom < atom->natoms; iatom++) {
for (int iatom = 0; iatom < sys_natoms; iatom++) {
std::string label = std::to_string(ibuf1all[iatom]);
int label_len = std::min(int(label.length()), MDI_LABEL_LENGTH);
strncpy(&labels[iatom * MDI_LABEL_LENGTH], label.c_str(), label_len);
}
}
int ierr = MDI_Send(labels, atom->natoms * MDI_LABEL_LENGTH, MDI_CHAR, mdicomm);
int ierr = MDI_Send(labels, sys_natoms * MDI_LABEL_LENGTH, MDI_CHAR, mdicomm);
if (ierr) error->all(FLERR, "MDI: <LABELS data");
delete[] labels;
@ -1393,8 +1393,7 @@ void MDIEngine::send_labels()
void MDIEngine::send_natoms()
{
int natoms = static_cast<int>(atom->natoms);
int ierr = MDI_Send(&natoms, 1, MDI_INT, mdicomm);
int ierr = MDI_Send(&sys_natoms, 1, MDI_INT, mdicomm);
if (ierr != 0) error->all(FLERR, "MDI: <NATOMS data");
}
@ -1435,7 +1434,7 @@ void MDIEngine::send_stress()
void MDIEngine::send_double1(int which)
{
memset(buf1, 0, atom->natoms * sizeof(double));
memset(buf1, 0, sys_natoms * sizeof(double));
// use atomID to index into ordered buf1
@ -1467,9 +1466,9 @@ void MDIEngine::send_double1(int which)
}
}
MPI_Reduce(buf1, buf1all, atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world);
MPI_Reduce(buf1, buf1all, sys_natoms, MPI_DOUBLE, MPI_SUM, 0, world);
int ierr = MDI_Send(buf1all, atom->natoms, MDI_DOUBLE, mdicomm);
int ierr = MDI_Send(buf1all, sys_natoms, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: <double1 data");
}
@ -1481,7 +1480,7 @@ void MDIEngine::send_double1(int which)
void MDIEngine::send_int1(int which)
{
memset(ibuf1, 0, atom->natoms * sizeof(int));
memset(ibuf1, 0, sys_natoms * sizeof(int));
// use atomID to index into ordered ibuf1
@ -1498,9 +1497,9 @@ void MDIEngine::send_int1(int which)
}
}
MPI_Reduce(ibuf1, ibuf1all, atom->natoms, MPI_INT, MPI_SUM, 0, world);
MPI_Reduce(ibuf1, ibuf1all, sys_natoms, MPI_INT, MPI_SUM, 0, world);
int ierr = MDI_Send(ibuf1all, atom->natoms, MDI_INT, mdicomm);
int ierr = MDI_Send(ibuf1all, sys_natoms, MDI_INT, mdicomm);
if (ierr) error->all(FLERR, "MDI: <int1 data");
}
@ -1512,7 +1511,7 @@ void MDIEngine::send_int1(int which)
void MDIEngine::send_double3(int which)
{
memset(buf3, 0, 3 * atom->natoms * sizeof(double));
memset(buf3, 0, 3 * sys_natoms * sizeof(double));
// use atomID to index into ordered buf3
@ -1547,9 +1546,9 @@ void MDIEngine::send_double3(int which)
}
}
MPI_Reduce(buf3, buf3all, 3 * atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world);
MPI_Reduce(buf3, buf3all, 3 * sys_natoms, MPI_DOUBLE, MPI_SUM, 0, world);
int ierr = MDI_Send(buf3all, 3 * atom->natoms, MDI_DOUBLE, mdicomm);
int ierr = MDI_Send(buf3all, 3 * sys_natoms, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: <double3 data");
}
@ -1663,7 +1662,8 @@ void MDIEngine::reallocate()
{
if (sys_natoms <= maxatom) return;
if (3 * sys_natoms > MAXSMALLINT) error->all(FLERR, "Natoms too large to use with mdi engine");
bigint nsize = (bigint) sys_natoms * 3;
if (nsize > MAXSMALLINT) error->all(FLERR, "Natoms too large to use with mdi engine");
maxatom = sys_natoms;

1
src/MDI/mdi_plugin.cpp
View File

@ -19,7 +19,6 @@
#include "mdi_plugin.h"
#include "error.h"
#include "fix_mdi_aimd.h"
#include "input.h"
#include "modify.h"