support two modes of QMMM coupling
This commit is contained in:
Steve Plimpton
@ -8,7 +8,7 @@ Syntax
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID mdi/qm keyword
|
||||
fix ID group-ID mdi/qm keyword value(s) keyword value(s) ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* mdi/qm = style name of this fix command
|
||||
@ -270,7 +270,8 @@ supports, e.g. *lj*, but then no unit conversion is performed.
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`mdi plugin <mdi>`, :doc:`mdi engine <mdi>`
|
||||
:doc:`mdi plugin <mdi>`, :doc:`mdi engine <mdi>`, :doc:`fix mdi/qmmm
|
||||
<fix_mdi_qmmm>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
252
doc/src/fix_mdi_qmmm.rst
Normal file
252
doc/src/fix_mdi_qmmm.rst
Normal file
@ -0,0 +1,252 @@
|
||||
.. index:: fix mdi/qm
|
||||
|
||||
fix mdi/qm command
|
||||
======================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID mdi/qmmm mode keyword value(s) keyword value(s) ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* mdi/qm = style name of this fix command
|
||||
* mode = *direct* or *potential*
|
||||
* zero or more keyword/value pairs may be appended
|
||||
* keyword = *virial* or *add* or *every* or *connect* or *elements*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*virial* args = *yes* or *no*
|
||||
yes = request virial tensor from server code
|
||||
no = do not request virial tensor from server code
|
||||
*connect* args = *yes* or *no*
|
||||
yes = perform a one-time connection to the MDI engine code
|
||||
no = do not perform the connection operation
|
||||
*elements* args = N_1 N_2 ... N_ntypes
|
||||
N_1,N_2,...N_ntypes = atomic number for each of ntypes LAMMPS atom types
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 all mdi/qmmm direct
|
||||
fix 1 all mdi/qmmm potential virial yes
|
||||
fix 1 all mdi/qmmm potential virial yes elements 13 29
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
This command enables LAMMPS to act as a client with another server
|
||||
code to perform a couple QMMM (quantum-mechanics molecular-mechanics)
|
||||
simulation. LAMMPS will perform classical MD (molecular-mechnanics)
|
||||
for the (typically larger) MM portion of the system. A quantum
|
||||
mechanics code will calculate quantum effects for the QM portion of
|
||||
the system. The QM server code must support use of the `MDI Library
|
||||
<https://molssi-mdi.github.io/MDI_Library/html/index.html>`_ as
|
||||
explained below.
|
||||
|
||||
The code coupling performed by this command is done via the `MDI
|
||||
Library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
|
||||
LAMMPS runs as an MDI driver (client), and sends MDI commands to an
|
||||
external MDI engine code (server), in this case a QM code which has
|
||||
support for MDI. See the :doc:`Howto mdi <Howto_mdi>` page for more
|
||||
information about how LAMMPS can operate as either an MDI driver or
|
||||
engine.
|
||||
|
||||
|
||||
|
||||
Q: provide an example where LAMMPS is also the QM code ?
|
||||
similar to fix mdi/qm ?
|
||||
|
||||
The examples/mdi directory contains input scripts using this fix in
|
||||
the various use cases discussed below. In each case, two instances of
|
||||
LAMMPS are used, once as an MDI driver, once as an MDI engine
|
||||
(surrogate for a QM code). The examples/mdi/README file explains how
|
||||
to launch two codes so that they communicate via the MDI library using
|
||||
either MPI or sockets. Any QM code that supports MDI could be used in
|
||||
place of LAMMPS acting as a QM surrogate. See the :doc:`Howto mdi
|
||||
<Howto_mdi>` page for a current list (March 2022) of such QM codes.
|
||||
|
||||
Note that an engine code can support MDI in either or both of two
|
||||
modes. It can be used as a stand-alone code, launched at the same
|
||||
time as LAMMPS. Or it can be used as a plugin library, which LAMMPS
|
||||
loads. See the :doc:`mdi plugin <mdi>` command for how to trigger
|
||||
LAMMPS to load a plugin library. The examples/mdi/README file
|
||||
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 *connect* keyword determines whether this fix performs a one-time
|
||||
connection to the QM code. The default is *yes*. The only time a
|
||||
*no* is needed is if this command is used multiple times in an input
|
||||
script. E.g. if it used inside a loop which also uses the :doc:`clear
|
||||
<clear>` command to destroy the system (including any defined fixes).
|
||||
See the examples/mdi/in.series.driver script as an example of this,
|
||||
where LAMMPS is using the QM code to compute energy and forces for a
|
||||
series of system configurations. In this use case *connect no*
|
||||
is used along with the :doc:`mdi connect and exit <mdi>` command
|
||||
to one-time initiate/terminate the connection outside the loop.
|
||||
|
||||
The *elements* keyword allows specification of what element each
|
||||
LAMMPS atom type corresponds to. This is specified by the atomic
|
||||
number of the element, e.g. 13 for Al. An atomic number must be
|
||||
specified for each of the ntypes LAMMPS atom types. Ntypes is
|
||||
typically specified via the create_box command or in the data file
|
||||
read by the read_data command.
|
||||
|
||||
If this keyword is specified, then this fix will send the MDI
|
||||
">ELEMENTS" command to the engine, to insure the two codes are
|
||||
consistent in their definition of atomic species. If this keyword is
|
||||
not specified, then this fix will send the MDI >TYPES command to the
|
||||
engine. This is fine if both the LAMMPS driver and the MDI engine are
|
||||
initialized so that the atom type values are consistent in both codes.
|
||||
|
||||
----------
|
||||
|
||||
The following 3 example use cases are illustrated in the examples/mdi
|
||||
directory. See its README file for more details.
|
||||
|
||||
(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.
|
||||
|
||||
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 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 includes :doc:`fix box relax <fix_box_relax>` 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 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 evaluation of energy and forces for a series of *N*
|
||||
independent systems (simulation box and atoms), use *add no* and
|
||||
*every 1*. Write a LAMMPS input script which loops over the *N*
|
||||
systems. See the :doc:`Howto multiple <Howto_multiple>` doc page for
|
||||
details on looping and removing old systems. The series of systems
|
||||
could be initialized by reading them from data files with
|
||||
:doc:`read_data <read_data>` commands. Or, for example, by using the
|
||||
:doc:`lattice <lattice>` , :doc:`create_atoms <create_atoms>`,
|
||||
:doc:`delete_atoms <delete_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 produce
|
||||
QM data for training a machine learning potential.
|
||||
|
||||
----------
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
No information about this fix is written to :doc:`binary restart files
|
||||
<restart>`.
|
||||
|
||||
The :doc:`fix_modify <fix_modify>` *energy* option is supported by
|
||||
this fix to add the potential energy computed by the QM code to the
|
||||
global potential energy of the system as part of :doc:`thermodynamic
|
||||
output <thermo_style>`. The default setting for this fix is
|
||||
:doc:`fix_modify energy yes <fix_modify>`, unless the *add* keyword is
|
||||
set to *no*, in which case the default setting is *no*.
|
||||
|
||||
The :doc:`fix_modify <fix_modify>` *virial* option is supported by
|
||||
this fix to add the contribution computed by the QM code to the global
|
||||
pressure of the system as part of :doc:`thermodynamic output
|
||||
<thermo_style>`. The default setting for this fix is :doc:`fix_modify
|
||||
virial yes <fix_modify>`, unless the *add* keyword is set to *no*, in
|
||||
which case the default setting is *no*.
|
||||
|
||||
This fix computes a global scalar which can be accessed by various
|
||||
:doc:`output commands <Howto_output>`. The scalar is the energy
|
||||
returned by the QM code. The scalar value calculated by this fix is
|
||||
"extensive".
|
||||
|
||||
This fix also computes a global vector with of length 6 which contains
|
||||
the symmetric virial tensor values returned by the QM code. It can
|
||||
likewise be accessed by various :doc:`output commands <Howto_output>`.
|
||||
|
||||
The ordering of values in the symmetric virial tensor is as follows:
|
||||
vxx, vyy, vzz, vxy, vxz, vyz. The values will be in pressure
|
||||
:doc:`units <units>`.
|
||||
|
||||
This fix also computes a peratom array with 3 columns which contains
|
||||
the peratom forces returned by the QM code. It can likewise be
|
||||
accessed by various :doc:`output commands <Howto_output>`.
|
||||
|
||||
No parameter of this fix can be used with the *start/stop* keywords of
|
||||
the :doc:`run <run>` command.
|
||||
|
||||
Assuming the *add* keyword is set to *yes* (the default), the forces
|
||||
computed by the QM code are used during an energy minimization,
|
||||
invoked by the :doc:`minimize <minimize>` command.
|
||||
|
||||
.. note::
|
||||
|
||||
If you want the potential energy associated with the QM forces to
|
||||
be included in the total potential energy of the system (the
|
||||
quantity being minimized), you MUST not disable the
|
||||
:doc:`fix_modify <fix_modify>` *energy* option for this fix, which
|
||||
means the *add* keyword should also be set to *yes* (the default).
|
||||
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This command is part of the MDI package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package
|
||||
<Build_package>` page for more info.
|
||||
|
||||
The QM code does not currently compute and return per-atom energy or
|
||||
per-atom virial contributions. So they will not show up as part of
|
||||
the calculations performed by the :doc:`compute pe/atom
|
||||
<compute_pe_atom>` or :doc:`compute stress/atom <compute_stress_atom>`
|
||||
commands.
|
||||
|
||||
To use LAMMPS as an MDI driver in conjunction with other MDI-enabled
|
||||
codes (MD or QM codes), the :doc:`units <units>` command should be
|
||||
used to specify *real* or *metal* units. This will ensure the correct
|
||||
unit conversions between LAMMPS and MDI units. The other code will
|
||||
also perform similar unit conversions into its preferred units.
|
||||
|
||||
LAMMPS can also be used as an MDI driver in other unit choices it
|
||||
supports, e.g. *lj*, but then no unit conversion is performed.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`mdi plugin <mdi>`, :doc:`mdi engine <mdi>`, :doc:`fix mdi/qm
|
||||
<fix_mdi_qm>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
The default for the optional keywords are virial = no and connect =
|
||||
yes.
|
||||
@ -29,6 +29,7 @@ using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
enum { NATIVE, REAL, METAL }; // LAMMPS units which MDI supports
|
||||
enum { DIRECT, POTENTIAL }; // mode of QMMM coupling
|
||||
|
||||
#define MAXELEMENT 103 // used elsewhere in MDI package
|
||||
|
||||
@ -61,13 +62,19 @@ FixMDIQMMM::FixMDIQMMM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
|
||||
if (role != MDI_DRIVER)
|
||||
error->all(FLERR, "Must invoke LAMMPS as an MDI driver to use fix mdi/qmmm");
|
||||
|
||||
// mode arg
|
||||
|
||||
if (strcmp(arg[3]),"direct" == 0) mode = DIRECT;
|
||||
else if (strcmp(arg[3],"potential") == 0) mode = POTENTIAL;
|
||||
error->all(FLERR,"Illegal fix mdi/qmmm command");
|
||||
|
||||
// optional args
|
||||
|
||||
virialflag = 0;
|
||||
connectflag = 1;
|
||||
elements = nullptr;
|
||||
|
||||
int iarg = 3;
|
||||
int iarg = 4;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg], "virial") == 0) {
|
||||
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix mdi/qmmm command");
|
||||
@ -290,8 +297,10 @@ int FixMDIQMMM::setmask()
|
||||
int mask = 0;
|
||||
mask |= POST_NEIGHBOR;
|
||||
mask |= MIN_POST_NEIGHBOR;
|
||||
mask |= PRE_FORCE;
|
||||
mask |= MIN_PRE_FORCE;
|
||||
if (mode == POTENTIAL) {
|
||||
mask |= PRE_FORCE;
|
||||
mask |= MIN_PRE_FORCE;
|
||||
}
|
||||
mask |= POST_FORCE;
|
||||
mask |= MIN_POST_FORCE;
|
||||
return mask;
|
||||
@ -389,13 +398,14 @@ void FixMDIQMMM::init()
|
||||
platform::walltime() - tstart);
|
||||
}
|
||||
|
||||
|
||||
// initial one-time MDI communication with engine
|
||||
|
||||
// send natoms, atom types or elements, and simulation box to engine
|
||||
// confirm engine count of NATOMS is correct
|
||||
// this will trigger setup of a new system
|
||||
// subsequent calls in post_force() will be for same system until new init()
|
||||
|
||||
// NOTE: why is this done here
|
||||
reallocate();
|
||||
|
||||
int natoms_exists;
|
||||
@ -422,6 +432,13 @@ void FixMDIQMMM::init()
|
||||
error->all(FLERR, "MDI: Engine has wrong atom count and does not support >NATOMS command");
|
||||
}
|
||||
|
||||
if (mode == DIRECT) {
|
||||
ierr = MDI_Send_command(">NLATTICE", mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >NLATTICE command");
|
||||
ierr = MDI_Send(&nmm, 1, MDI_INT, mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >NLATTICE data");
|
||||
}
|
||||
|
||||
int elements_exists;
|
||||
int types_exists;
|
||||
ierr = MDI_Check_command_exists("@DEFAULT", ">ELEMENTS", mdicomm, &elements_exists);
|
||||
@ -467,11 +484,15 @@ void FixMDIQMMM::post_neighbor()
|
||||
set_qm2owned();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* ----------------------------------------------------------------------
|
||||
only invoked in POTENTIAL mode
|
||||
calculates Coulomb potential for each QM atom
|
||||
invokes the QM code
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void FixMDIQMMM::pre_force(int vflag)
|
||||
{
|
||||
int ilocal,jlocal;
|
||||
int ilocal,jlocal;
|
||||
double rsq;
|
||||
double delta[3];
|
||||
|
||||
@ -656,9 +677,79 @@ void FixMDIQMMM::pre_force(int vflag)
|
||||
update->minimize->force_clear();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* ----------------------------------------------------------------------
|
||||
different methods invoked for DIRECT vs POTENTIAL mode
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void FixMDIQMMM::post_force(int vflag)
|
||||
{
|
||||
if (mode == DIRECT) post_force_direct(vflag);
|
||||
else if (mode = POTENTIAL post_force_potential(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixMDIQMMM::post_force_direct(int vflag)
|
||||
{
|
||||
if (comm->me == 0) utils::logmesg(lmp, "Calling QM code ...\n");
|
||||
|
||||
MPI_Barrier(world);
|
||||
double tstart = platform::walltime();
|
||||
|
||||
// MDI calls
|
||||
|
||||
// send current coords of QM atoms to MDI engine
|
||||
|
||||
int ierr = MDI_Send_command(">COORDS", mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >COORDS command");
|
||||
ierr = MDI_Send(&xqm[0][0], 3 * nqm, MDI_DOUBLE, mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >COORDS data");
|
||||
|
||||
// send current coords of MM atoms to MDI engine
|
||||
|
||||
int ierr = MDI_Send_command(">CLATTICE", mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >CLATTICE command");
|
||||
ierr = MDI_Send(&xmm[0][0], 3 * nmm, MDI_DOUBLE, mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >CLATTICE data");
|
||||
|
||||
// send charges on MM atoms to MDI engine
|
||||
|
||||
int ierr = MDI_Send_command(">LATTICE", mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >LATTICE command");
|
||||
ierr = MDI_Send(&qmm[0][0], nmm, MDI_DOUBLE, mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: >LATTICE data");
|
||||
|
||||
// request QM 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);
|
||||
|
||||
// request forces on QM atoms from MDI engine
|
||||
// NOTE: will this be forces on all atoms in DIRECT mode ?
|
||||
|
||||
ierr = MDI_Send_command("<FORCES", mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: <FORCES command");
|
||||
ierr = MDI_Recv(&fqm[0][0], 3 * nqm, MDI_DOUBLE, mdicomm);
|
||||
if (ierr) error->all(FLERR, "MDI: <FORCES data");
|
||||
MPI_Bcast(&fqm[0][0], 3 * nqm, MPI_DOUBLE, 0, world);
|
||||
|
||||
// end of MDI calls
|
||||
|
||||
MPI_Barrier(world);
|
||||
if (comm->me == 0)
|
||||
utils::logmesg(lmp, " time = {:.3f} seconds\n",
|
||||
platform::walltime() - tstart);
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixMDIQMMM::post_force_potential(int vflag)
|
||||
{
|
||||
// int ilocal,jlocal;
|
||||
// double rsq,r2inv,rinv,fpair;
|
||||
@ -757,7 +848,10 @@ void FixMDIQMMM::post_force(int vflag)
|
||||
//c_pe->addstep(update->ntimestep+1);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* ----------------------------------------------------------------------
|
||||
NOTE: remove this method, do it via fix mdi/qm ?
|
||||
or is it here for debugging ?
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void FixMDIQMMM::post_force_aimd(int vflag)
|
||||
{
|
||||
|
||||
@ -114,6 +114,8 @@ class FixMDIQMMM : public Fix {
|
||||
void send_box();
|
||||
void unit_conversions();
|
||||
|
||||
void post_force_direct(int);
|
||||
void post_force_potential(int);
|
||||
void post_force_aimd(int);
|
||||
};
|
||||
|
||||
|
||||
Reference in New Issue
Block a user