diff --git a/doc/src/fix_mdi_qm.rst b/doc/src/fix_mdi_qm.rst index 667dea710a..46671ceee6 100644 --- a/doc/src/fix_mdi_qm.rst +++ b/doc/src/fix_mdi_qm.rst @@ -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 ` 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 `, :doc:`mdi engine ` +:doc:`mdi plugin `, :doc:`mdi engine `, :doc:`fix mdi/qmmm + ` Default """"""" diff --git a/doc/src/fix_mdi_qmmm.rst b/doc/src/fix_mdi_qmmm.rst new file mode 100644 index 0000000000..46a20362fc --- /dev/null +++ b/doc/src/fix_mdi_qmmm.rst @@ -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 ` 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 +`_ as +explained below. + +The code coupling performed by this command is done via the `MDI +Library `_. +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 ` 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 +` 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 ` 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 +` 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 ` 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 ` 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 ` 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 +` 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 ` 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 ` commands. Or, for example, by using the +:doc:`lattice ` , :doc:`create_atoms `, +:doc:`delete_atoms `, and/or :doc:`displace_atoms +random ` commands to generate a series of different +systems. At the end of the loop perform :doc:`run 0 ` and +:doc:`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 +`. + +The :doc:`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 `. The default setting for this fix is +:doc:`fix_modify energy yes `, unless the *add* keyword is +set to *no*, in which case the default setting is *no*. + +The :doc:`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 +`. The default setting for this fix is :doc:`fix_modify +virial yes `, 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 `. 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 `. + +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 `. + +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 `. + +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`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 ` 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 ` *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 +` 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 +` or :doc:`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 ` 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 `, :doc:`mdi engine `, :doc:`fix mdi/qm + ` + +Default +""""""" + +The default for the optional keywords are virial = no and connect = +yes. diff --git a/src/MDI/fix_mdi_qmmm.cpp b/src/MDI/fix_mdi_qmmm.cpp index 537ecfcebd..431991429e 100644 --- a/src/MDI/fix_mdi_qmmm.cpp +++ b/src/MDI/fix_mdi_qmmm.cpp @@ -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("all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: all(FLERR, "MDI: 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) { diff --git a/src/MDI/fix_mdi_qmmm.h b/src/MDI/fix_mdi_qmmm.h index f1c447a9d0..f3902f121b 100644 --- a/src/MDI/fix_mdi_qmmm.h +++ b/src/MDI/fix_mdi_qmmm.h @@ -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); };