Merge pull request #2611 from MolSSI-MDI/mdi

Add MolSSI Driver Interface Support
This commit is contained in:
Axel Kohlmeyer
2021-05-14 13:18:52 -04:00
committed by GitHub
41 changed files with 2807 additions and 50 deletions

View File

@ -121,7 +121,7 @@ set(STANDARD_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE
PLUGIN QEQ REPLICA RIGID SHOCK SPIN SNAP SRD KIM PYTHON MSCG MPIIO VORONOI
USER-ADIOS USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-MESODPD USER-CGSDK
USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD
USER-LB USER-MANIFOLD USER-MEAMC USER-MESONT USER-MGPT USER-MISC USER-MOFFF
USER-LB USER-MANIFOLD USER-MDI USER-MEAMC USER-MESONT USER-MGPT USER-MISC USER-MOFFF
USER-MOLFILE USER-NETCDF USER-PHONON USER-PLUMED USER-PTM USER-QTB
USER-REACTION USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD USER-SMTBQ USER-SPH
USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM USER-YAFF USER-PACE USER-BROWNIAN)
@ -324,8 +324,8 @@ else()
set(CUDA_REQUEST_PIC)
endif()
foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MOLFILE USER-NETCDF USER-PLUMED USER-QMMM
USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS USER-PACE)
foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MDI USER-MOLFILE USER-NETCDF USER-PLUMED
USER-QMMM USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS USER-PACE)
if(PKG_${PKG_WITH_INCL})
include(Packages/${PKG_WITH_INCL})
endif()

View File

@ -0,0 +1,66 @@
find_package(mdi QUIET)
if(${mdi_FOUND})
set(DOWNLOAD_MDI_DEFAULT OFF)
else()
set(DOWNLOAD_MDI_DEFAULT ON)
endif()
option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an already installed one" ${DOWNLOAD_MDI_DEFAULT})
if(DOWNLOAD_MDI)
message(STATUS "MDI download requested - we will build our own")
set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.2.9.tar.gz" CACHE STRING "URL for MDI tarball")
set(MDI_MD5 "ddfa46d6ee15b4e59cfd527ec7212184" CACHE STRING "MD5 checksum for MDI tarball")
mark_as_advanced(MDI_URL)
mark_as_advanced(MDI_MD5)
set(LAMMPS_LIB_MDI_BIN_DIR ${LAMMPS_LIB_BINARY_DIR}/mdi)
include(ExternalProject)
message(STATUS "Building mdi.")
ExternalProject_Add(mdi_external
URL ${MDI_URL}
URL_MD5 ${MDI_MD5}
UPDATE_COMMAND ""
CMAKE_ARGS ${CMAKE_REQUEST_PIC}
-DCMAKE_INSTALL_PREFIX=${LAMMPS_LIB_MDI_BIN_DIR}
-DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
-DCMAKE_INSTALL_LIBDIR=${CMAKE_INSTALL_LIBDIR}
-DCMAKE_INSTALL_INCLUDEDIR=${CMAKE_INSTALL_INCLUDEDIR}
-DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS}
-Dlanguage=C
CMAKE_CACHE_ARGS -DCMAKE_C_FLAGS:STRING=${CMAKE_C_FLAGS}
-DCMAKE_CXX_FLAGS:STRING=${CMAKE_CXX_FLAGS}
-DTargetOpenMP_FIND_COMPONENTS:STRING=C;CXX)
# Link the lammps library against MDI
target_include_directories(lammps PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/${CMAKE_INSTALL_INCLUDEDIR}/mdi)
target_link_directories(lammps PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/${CMAKE_INSTALL_LIBDIR}/mdi)
target_link_libraries(lammps PRIVATE mdi)
add_dependencies(lammps mdi_external)
# Link the lammps executable against MDI
target_include_directories(lmp PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/${CMAKE_INSTALL_INCLUDEDIR}/mdi)
target_link_directories(lmp PRIVATE ${LAMMPS_LIB_MDI_BIN_DIR}/${CMAKE_INSTALL_LIBDIR}/mdi)
target_link_libraries(lmp PRIVATE mdi)
add_dependencies(lmp mdi_external)
else()
find_package(mdi)
if(NOT mdi_FOUND)
message(FATAL_ERROR "MDI library not found. Help CMake to find it "
"by setting mdi_LIBRARY and mdi_INCLUDE_DIR, or set DOWNLOAD_MDI=ON "
"to download and compile it")
endif()
# Link the lammps library against MDI
target_include_directories(lammps PRIVATE ${mdi_INCLUDE_DIR})
target_link_libraries(lammps PRIVATE ${mdi_LIBRARY})
# Link the lammps executable against MDI
target_include_directories(lmp PRIVATE ${mdi_INCLUDE_DIR})
target_link_libraries(lmp PRIVATE ${mdi_LIBRARY})
endif()
target_compile_definitions(lammps PRIVATE -DLMP_USER_MDI)
target_compile_definitions(lmp PRIVATE -DLMP_USER_MDI)

View File

@ -49,6 +49,7 @@ This is the list of packages that may require additional steps.
* :ref:`USER-COLVARS <user-colvars>`
* :ref:`USER-H5MD <user-h5md>`
* :ref:`USER-INTEL <user-intel>`
* :ref:`USER-MDI <user-mdi>`
* :ref:`USER-MESONT <user-mesont>`
* :ref:`USER-MOLFILE <user-molfile>`
* :ref:`USER-NETCDF <user-netcdf>`
@ -1539,6 +1540,35 @@ TBB and MKL.
----------
.. _user-mdi:
USER-MDI package
-----------------------------
.. tabs::
.. tab:: CMake build
.. code-block:: bash
-D DOWNLOAD_MDI=value # download MDI Library for build, value = no (default) or yes
.. tab:: Traditional make
Before building LAMMPS, you must build the MDI Library in
``lib/mdi``\ . You can do this by executing a command like one
of the following from the ``lib/mdi`` directory:
.. code-block:: bash
$ python Install.py -m gcc # build using gcc compiler
$ python Install.py -m icc # build using icc compiler
The build should produce two files: ``lib/mdi/includelink/mdi.h``
and ``lib/mdi/liblink/libmdi.so``\ .
----------
.. _user-mesont:
USER-MESONT package

View File

@ -67,6 +67,7 @@ An alphabetic list of all general LAMMPS commands.
* :doc:`lattice <lattice>`
* :doc:`log <log>`
* :doc:`mass <mass>`
* :doc:`mdi/engine <mdi_engine>`
* :doc:`message <message>`
* :doc:`minimize <minimize>`
* :doc:`min_modify <min_modify>`

View File

@ -101,6 +101,7 @@ OPT.
* :doc:`lb/viscous <fix_lb_viscous>`
* :doc:`lineforce <fix_lineforce>`
* :doc:`manifoldforce <fix_manifoldforce>`
* :doc:`mdi/engine <fix_mdi_engine>`
* :doc:`meso/move <fix_meso_move>`
* :doc:`momentum (k) <fix_momentum>`
* :doc:`momentum/chunk <fix_momentum>`

View File

@ -23,6 +23,7 @@ General howto
Howto_library
Howto_couple
Howto_client_server
Howto_mdi
Settings howto
==============

132
doc/src/Howto_mdi.rst Normal file
View File

@ -0,0 +1,132 @@
Using LAMMPS with the MDI library for code coupling
===================================================
..note::
This Howto doc page will eventually replace the
:doc:`Howto client/server <Howto_client_server>` doc page.
Client/server coupling of two codes is where one code is the "client"
and sends request messages (data) to a "server" code. The server
responds to each request with a reply message. This enables the two
codes to work in tandem to perform a simulation. LAMMPS can act as
either a client or server code; it does this by using the `MolSSI
Driver Interface (MDI) library
<https://molssi-mdi.github.io/MDI_Library/html/index.html>`_,
developed by the `Molecular Sciences Software Institute (MolSSI)
<https://molssi.org>`_.
Alternate methods for code coupling with LAMMPS are described on the
:doc:`Howto couple <Howto_couple>` doc page.
Some advantages of client/server coupling are that the two codes can run
as stand-alone executables; they need not be linked together. Thus
neither code needs to have a library interface. This also makes it easy
to run the two codes on different numbers of processors. If a message
protocol (format and content) is defined for a particular kind of
simulation, then in principle any code which implements the client-side
protocol can be used in tandem with any code which implements the
server-side protocol. Neither code needs to know what specific other
code it is working with.
In MDI nomenclature, a client code is the "driver", and a server code is
an "engine". One driver code can communicate with one or more instances
of one or more engine codes. Driver and engine codes can be written in
any language: C, C++, Fortran, Python, etc.
In addition to allowing driver and engine(s) running to run as
stand-alone executables, MDI also enables a server code to be a
"plugin" to the client code. In this scenario, server code(s) are
compiled as shared libraries, and one (or more) instances of the
server are instantiated by the driver code. If the driver code runs
in parallel, it can split its MPI communicator into multiple
sub-communicators, and launch each plugin engine instance on a
sub-communicator. Driver processors in that sub-communicator exchange
messages with that engine instance, and can also send MPI messages to
other processors in the driver. The driver code can also destroy
engine instances and re-instantiate them.
The way that a driver communicates with an engine is by making
MDI_Send() and MDI_Recv() calls, which are conceptually similar to
MPI_Send() and MPI_Recv() calls. Each send or receive has a string
which identifies the command name, and optionally some data, which can
be a single value or vector of values of any data type. Inside the
MDI library, data is exchanged between the driver and engine via MPI
calls or sockets. This a run-time choice by the user.
-------------
As an example, LAMMPS and the ``pw.x`` command from Quantum Espresso (a
suite of quantum DFT codes), can work together via the MDI library to
perform an ab initio MD (AIMD) simulation, where LAMMPS runs an MD
simulation and sends a message each timestep to ``pw.x`` asking it to
compute quantum forces on the current configuration of atoms. Here is
how the 2 codes are launched to communicate by MPI:
.. code-block:: bash
% mpirun -np 2 lmp_mpi -mdi "-role DRIVER -name d -method MPI" \
-in in.aimd : -np 16 pw.x -in qe.in -mdi "-role ENGINE -name e -method MPI"
In this case LAMMPS runs on 2 processors (MPI tasks), ``pw.x`` runs on 16
processors.
Here is how the 2 codes are launched to communicate by sockets:
.. code-block:: bash
% mpirun -np 2 lmp_mpi -mdi "-role DRIVER -name d -method TCP -port 8021" -in in.aimd
% mpirun -np 16 pw.x -in qe.in -mdi "-role ENGINE -name e -method TCP -port 8021 -hostname localhost"
These commands could be issued in different windows on a desktop
machine. Or in the same window, if the first command is ended with
"&" so as to run in the background. If "localhost" is replaced by an
IP address, ``pw.x`` could be run on another machine on the same network, or
even on another machine across the country.
After both codes initialize themselves to model the same system, this is
what occurs each timestep:
* LAMMPS send a ">COORDS" message to ``pw.x`` with a 3*N vector of current atom coords
* ``pw.x`` receives the message/coords and computes quantum forces on all the atoms
* LAMMPS send a "<FORCES" message to ``pw.x`` and waits for the result
* ``pw.x`` receives the message (after its computation finishes) and sends a 3*N vector of forces
* LAMMPS receives the forces and time integrates to complete a single timestep
-------------
Examples scripts for using LAMMPS as an MDI engine are in the
examples/mdi directory. See the README file in that directory for
instructions on how to run the examples.
..note::
Work is underway to add commands that allow LAMMPS to be used as an
MDI driver, e.g. for the AIMD example discussed above. Example
scripts for this usage mode will be added the same directory when
available.
If LAMMPS is used as a stand-alone engine it should set up the system
it will be modeling in its input script, then invoke the
:doc:`mdi/engine <mdi_engine>` command. This will put LAMMPS into
engine mode where it waits for messages and data from the driver.
When the driver sends an "EXIT" command, LAMMPS will exit engine mode
and the input script will continue.
If LAMMPS is used as a plugin engine it operates the same way, except
that the driver will pass LAMMPS an input script to initialize itself.
Upon receiving the "EXIT" command, LAMMPS will exit engine mode and the
input script will continue. After finishing execution of the input
script, the instance of LAMMPS will be destroyed.
LAMMPS supports the full set of MD-appropriate engine commands defined
by the MDI library. See the :doc:`mdi/engine <mdi_engine>` doc page for
a list of these.
If those commands are not sufficient for a user-developed driver to use
LAMMPS as an engine, then new commands can be easily added. See these
two files which implement the definition of MDI commands and the logic
for responding to them:
* src/MDI/mdi_engine.cpp
* src/MDI/fix_mdi_engine.cpp

View File

@ -82,6 +82,7 @@ page gives those details.
* :ref:`USER-INTEL <PKG-USER-INTEL>`
* :ref:`USER-LB <PKG-USER-LB>`
* :ref:`USER-MANIFOLD <PKG-USER-MANIFOLD>`
* :ref:`USER-MDI <PKG-USER-MDI>`
* :ref:`USER-MEAMC <PKG-USER-MEAMC>`
* :ref:`USER-MESODPD <PKG-USER-MESODPD>`
* :ref:`USER-MESONT <PKG-USER-MESONT>`
@ -1791,6 +1792,28 @@ Waltham, MA, USA)
----------
.. _PKG-USER-MDI:
USER-MDI package
----------------
**Contents:**
A LAMMPS command and fix to allow client-server coupling of LAMMPS to
other atomic or molecular simulation codes via the `MolSSI Driver Interface
(MDI) library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
**Author:** Taylor Barnes - MolSSI, taylor.a.barnes at gmail.com
**Supporting info:**
* src/USER-MDI/README
* :doc:`mdi/engine <mdi_engine>`
* :doc:`fix mdi/engine <fix_mdi_engine>`
* examples/USER/mdi
----------
.. _PKG-USER-MEAMC:
USER-MEAMC package

View File

@ -65,6 +65,8 @@ package:
+------------------------------------------------+-----------------------------------------------------------------+-------------------------------------------------------------------------------+------------------------------------------------------+---------+
| :ref:`USER-MANIFOLD <PKG-USER-MANIFOLD>` | motion on 2d surfaces | :doc:`fix manifoldforce <fix_manifoldforce>` | USER/manifold | no |
+------------------------------------------------+-----------------------------------------------------------------+-------------------------------------------------------------------------------+------------------------------------------------------+---------+
| :ref:`USER-MDI <PKG-USER-MDI>` | client-server coupling | :doc:`MDI Howto <Howto_mdi>` | USER/mdi | ext |
+------------------------------------------------+-----------------------------------------------------------------+-------------------------------------------------------------------------------+------------------------------------------------------+---------+
| :ref:`USER-MEAMC <PKG-USER-MEAMC>` | modified EAM potential (C++) | :doc:`pair_style meam/c <pair_meamc>` | meamc | no |
+------------------------------------------------+-----------------------------------------------------------------+-------------------------------------------------------------------------------+------------------------------------------------------+---------+
| :ref:`USER-MESODPD <PKG-USER-MESODPD>` | mesoscale DPD models | :doc:`pair_style edpd <pair_mesodpd>` | USER/mesodpd | no |

View File

@ -34,7 +34,7 @@ Examples
comm_modify mode multi reduce/multi
comm_modify mode multi group solvent
comm_modift mode multi cutoff/multi 1 10.0 cutoff/multi 2*4 15.0
comm_modify mode multi cutoff/multi 1 10.0 cutoff/multi 2*4 15.0
comm_modify vel yes
comm_modify mode single cutoff 5.0 vel yes
comm_modify cutoff/multi * 0.0
@ -75,8 +75,8 @@ communicated. in *multi/old*, a similar technique is used but atoms
are grouped by atom type. See the :doc:`neighbor multi <neighbor>` and
:doc:`neighbor multi/old <neighbor>` commands for
neighbor list construction options that may also be beneficial for
simulations of this kind. The *multi* communiction mode is only compatable
with the *multi* neighbor style. The *multi/old* communication mode is compatble
simulations of this kind. The *multi* communication mode is only compatible
with the *multi* neighbor style. The *multi/old* communication mode is comparable
with both the *multi* and *multi/old* neighbor styles.
The *cutoff* keyword allows you to extend the ghost cutoff distance
@ -97,18 +97,18 @@ warning is printed, if this bond based estimate is larger than the
communication cutoff used.
The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to
communication mode *multi* instead. Since the communication
cutoffs are determined per atom collections, a collection specifier is needed and
cutoff for one or multiple collections can be extended. Also ranges of collections
using the usual asterisk notation can be given.
Collections are indexed from 1 to N where N is the total number of collections.
communication mode *multi* instead. Since the communication cutoffs are
determined per atom collections, a collection specifier is needed and
cutoff for one or multiple collections can be extended. Also ranges of
collections using the usual asterisk notation can be given. Collections
are indexed from 1 to N where N is the total number of collections.
Note that the arguments for *cutoff/multi* are parsed right before each
simulation to account for potential changes in the number of collections.
Custom cutoffs are preserved between runs but if collections are redefined,
one may want to respecify communication cutoffs.
For granular pair styles,the default cutoff is set to the sum of the
current maximum atomic radii for each collection.
The *cutoff/multi/old* option is similar to *cutoff/multi* except it
simulation to account for potential changes in the number of
collections. Custom cutoffs are preserved between runs but if
collections are redefined, one may want to re-specify the communication
cutoffs. For granular pair styles,the default cutoff is set to the sum
of the current maximum atomic radii for each collection. The
*cutoff/multi/old* option is similar to *cutoff/multi* except it
operates on atom types as opposed to collections.
The *reduce/multi* option applies to *multi* and sets the communication
@ -147,7 +147,7 @@ ghost cutoff should be set.
In the last scenario, a :doc:`fix <fix>` or :doc:`compute <compute>` or
:doc:`pairwise potential <pair_style>` needs to calculate with ghost
atoms beyond the normal pairwise cutoff for some computation it
performs (e.g. locate neighbors of ghost atoms in a multibody pair
performs (e.g. locate neighbors of ghost atoms in a manybody pair
potential). Setting the ghost cutoff appropriately can insure it will
find the needed atoms.

View File

@ -59,6 +59,7 @@ Commands
lattice
log
mass
mdi_engine
message
min_modify
min_spin

View File

@ -0,0 +1,59 @@
.. index:: fix move
fix mdi/engine command
======================
Syntax
""""""
.. parsed-literal::
fix ID group-ID mdi/engine
* ID, group-ID are documented in :doc:`fix <fix>` command
* mdi/engine = style name of this fix command
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all mdi/engine
Description
"""""""""""
This fix is used along with the :doc:`mdi/engine <mdi_engine>` command
to enable LAMMPS to use the `MDI Library
<https://molssi-mdi.github.io/MDI_Library/html/index.html>`_ to run as
an MDI engine. The fix provides hooks that enable MDI driver codes to
communicate with LAMMPS at various points within a LAMMPS timestep.
It is not generally necessary to add this fix to a LAMMPS input file,
even when using the :doc:`mdi/engine <mdi_engine>` command. If the
:doc:`mdi/engine <mdi_engine>` command is executed and this fix is not
present, it will automatically be added and applied as a new fix for
all atoms for the duration of the command. Thus it is only necessary
to add this fix to an input file when you want to modify the group-ID
or the ordering of this fix relative to other fixes in the input script.
For more information about running LAMMPS as an MDI engine, see the
:doc:`mdi/engine <mdi_engine>` command and the :doc:`Howto mdi
<Howto_mdi>` doc page.
Restrictions
""""""""""""
This command is part of the USER-MDI package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`mdi/engine <mdi_engine>`
Default
"""""""
none

88
doc/src/mdi_engine.rst Normal file
View File

@ -0,0 +1,88 @@
.. index:: mdi_engine
mdi_engine command
==================
Syntax
""""""
.. parsed-literal::
mdi_engine
Description
"""""""""""
This command is used to have LAMMPS act as a server with another
client code to effectively couple the two codes together in
client/server mode.
More specifically, this command causes LAMMPS to begin using the `MDI
Library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_
to run as an MDI engine (server), responding to commands made by an
external MDI driver code (client). See the :doc:`Howto mdi
<Howto_mdi>` doc page for more information about how LAMMPS can work
as both an MDI driver or engine.
General information about launching codes that communicate using the
MDI Library can be found in the `corresponding page
<https://molssi-mdi.github.io/MDI_Library/html/library_page.html#library_launching_sec>`_
of the MDI Library's documentation.
----------
This command should typically be used in an input script after LAMMPS
has setup the system it is going to model in collaboration with the
driver code. Depending on how the driver code tells the LAMMPS engine
to exit, other commands can be executed after this command, but
typically it should be used at the end of the LAMMPS input script.
To act as a MD-based MDI engine, this is the list of MDI commands from
a driver code which LAMMPS currently recognizes. See more details
about these commands in the `MDI library documentation
<https://molssi-mdi.github.io/MDI_Library/html/mdi_standard.html>`_
.. NOTE: Taylor - is this the best link for this info? Can we flesh this
.. out with the full list of supported commands? Maybe the distinction
.. of what "node" the commands refer to is not needed in this table?
.. list-table::
:widths: 20 80
:header-rows: 1
* - Command name
- Action
* - >NATOMS
- Driver sends the number of atoms in the system
* - <NATOMS
- Driver requests the number of atoms in the system
* - <COORDS
- Driver requests 3*N double-precision atom coordinates
* - >FORCES
- Driver sends 3*N double-precision atom forces
* - <COORDS
- Driver requests 3*N double-precision atom forces
* - EXIT
- Driver tells the engine (LAMMPS) to exit engine mode
If these commands are not sufficient to support what a driver which
you write needs, additional commands can be defined by simply using a
new command name not in this list. Code to support the new command
needs to be added to the USER-MDI package within LAMMPS; see its
src/USER-MDI/mdi_engine.cpp and fix_mdi_engine.cpp files.
Restrictions
""""""""""""
This command is part of the USER-MDI package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`fix mdi/engine <fix_mdi_engine>`
Default
"""""""
None

View File

@ -221,14 +221,15 @@ from 1 to n (inclusive). A trailing asterisk means all types from n to M
(inclusive). A middle asterisk means all types from m to n (inclusive).
Note that all atom types must be included in exactly one of the N collections.
The *collection/interval* option provides a similar capability.
This command allows a user to define collections by specifying a
series of cutoff intervals. LAMMPS will automatically sort atoms into these intervals
based on their type-dependent cutoffs or their finite size.
You must first specify the number of collections N to be
defined followed by N values representing the upper cutoff of each interval.
This command is particularly useful for granular pairstyles where the interaction distance
of particles depends on their radius and may not depend on their atom type.
The *collection/interval* option provides a similar capability. This
command allows a user to define collections by specifying a series of
cutoff intervals. LAMMPS will automatically sort atoms into these
intervals based on their type-dependent cutoffs or their finite size.
You must first specify the number of collections N to be defined
followed by N values representing the upper cutoff of each interval.
This command is particularly useful for granular pair styles where the
interaction distance of particles depends on their radius and may not
depend on their atom type.
Restrictions
""""""""""""

View File

@ -47,6 +47,7 @@ Agnolin
Ai
Aidan
aij
aimd
airebo
Aj
ajs
@ -196,6 +197,7 @@ Ballenegger
Bammann
Banna
Barashev
barnes
barostat
Barostats
barostatted
@ -1207,6 +1209,7 @@ Halver
Hamaker
Hamel
Hammerschmidt
Hanley
haptic
Hara
Harpertown
@ -1269,6 +1272,7 @@ holonomic
Homebrew
hooke
Hookean
hostname
hotpink
Houlle
howto
@ -1728,6 +1732,7 @@ Lmpsdata
lmptype
LMT
ln
localhost
localTemp
localvectors
Loewen
@ -1858,6 +1863,8 @@ mc
McLachlan
md
mdf
MDI
mdi
mdpd
mDPD
meam
@ -3143,6 +3150,7 @@ Tanmoy
Tartakovsky
taskset
taubi
taylor
tb
tchain
Tchain

24
examples/USER/mdi/README Normal file
View File

@ -0,0 +1,24 @@
This dir contains scripts that demonstrate how to use LAMMPS as an
MDI engine. LAMMPS as an engine performs the MD timestepping.
The driver is a simple Python script. Every timestep the driver
sends one or more commands to LAMMPS.
--------------
The Script.sh file has comands to perform some very simple example
runs.
--------------
More complex calculations using LAMMPS as an MDI engine will
typically require the use of an MDI driver. Several MDI drivers
support calculations with LAMMPS, and include:
Ab Initio Molecular Dynamics (AIMD) Driver:
https://github.com/MolSSI-MDI/MDI_AIMD_Driver
Nudged Elastic Band (NEB) Driver:
https://github.com/MolSSI-MDI/MDI_NEB_Driver
Metadynamics Driver:
https://github.com/MolSSI-MDI/MDI_Metadynamics

View File

@ -0,0 +1,16 @@
#!/bin/bash
# sample launch scripts
# TCP, running LAMMPS on one proc
python driver.py -mdi "-name driver -role DRIVER -method TCP -port 8021" &
../../../src/lmp_mdi -mdi "-name LAMMPS -role ENGINE -method TCP -port 8021 -hostname localhost" -in lammps.in > lammps.out &
wait
# TCP, running LAMMPS on two procs
python driver.py -mdi "-name driver -role DRIVER -method TCP -port 8021" &
mpiexec -n 2 ../../../src/lmp_mdi -mdi "-name LAMMPS -role ENGINE -method TCP -port 8021 -hostname localhost" -in lammps.in > lammps.out &
wait

View File

@ -0,0 +1,24 @@
import sys
import mdi
use_mpi4py = False
try:
from mpi4py import MPI
use_mpi4py = True
except:
pass
# Initialize the MDI Library
mdi.MDI_Init(sys.argv[2])
# Connect to the engine
comm = mdi.MDI_Accept_communicator()
# Determine the name of the engine
mdi.MDI_Send_Command("<NAME", comm)
name = mdi.MDI_Recv(mdi.MDI_NAME_LENGTH, mdi.MDI_CHAR, comm)
print("Engine name: " + str(name))
# Send the "EXIT" command to the engine
mdi.MDI_Send_Command("EXIT", comm)

View File

@ -0,0 +1,92 @@
LAMMPS data file for water
24 atoms
16 bonds
8 angles
0 dihedrals
0 impropers
2 atom types
1 bond types
1 angle types
0 dihedral types
0 improper types
0.0 4.9325 xlo xhi
0.0 4.9325 ylo yhi
0.0 4.9325 zlo zhi
Masses
1 15.9994
2 1.008
Pair Coeffs
1 0.102 3.188
2 0.000 0.000
Bond Coeffs
1 450 0.9572
Angle Coeffs
1 55.0 104.52
Atoms
1 0 1 -0.83400 2.17919 0.196156 4.15513
2 0 2 0.41700 2.29785 4.8353 0.126003
3 0 2 0.41700 1.82037 1.07996 4.23498
4 0 1 -0.83400 4.65839 0.120414 0.305758
5 0 2 0.41700 4.67446 -0.0220991 4.29186
6 0 2 0.41700 4.28188 0.994196 0.410515
7 0 1 -0.83400 3.65045 2.40907 0.344349
8 0 2 0.41700 3.52052 2.1838 4.35565
9 0 2 0.41700 4.26579 3.14208 0.327669
10 0 1 -0.83400 1.21327 2.62177 4.15519
11 0 2 0.41700 1.47452 3.53837 4.0667
12 0 2 0.41700 1.20743 2.46396 0.16677
13 0 1 -0.83400 4.45777 4.47325 2.74192
14 0 2 0.41700 4.53396 4.49652 1.78804
15 0 2 0.41700 4.21354 3.56943 2.94119
16 0 1 -0.83400 2.04119 4.41585 1.64725
17 0 2 0.41700 2.26934 4.77582 2.50434
18 0 2 0.41700 1.69079 3.54574 1.83793
19 0 1 -0.83400 3.73384 1.97964 2.81949
20 0 2 0.41700 3.41083 2.22014 1.95113
21 0 2 0.41700 3.91914 1.04272 2.75561
22 0 1 -0.83400 1.20859 2.09853 1.68186
23 0 2 0.41700 1.01865 2.25693 2.60655
24 0 2 0.41700 1.16884 1.14674 1.58832
Bonds
1 1 1 2
2 1 1 3
3 1 4 5
4 1 4 6
5 1 7 8
6 1 7 9
7 1 10 11
8 1 10 12
9 1 13 14
10 1 13 15
11 1 16 17
12 1 16 18
13 1 19 20
14 1 19 21
15 1 22 23
16 1 22 24
Angles
1 1 2 1 3
2 1 5 4 6
3 1 8 7 9
4 1 11 10 12
5 1 14 13 15
6 1 17 16 18
7 1 20 19 21
8 1 23 22 24

View File

@ -0,0 +1,28 @@
units real
neigh_modify delay 0 every 1 check yes
atom_style full
bond_style harmonic
angle_style harmonic
pair_style lj/cut/coul/long 10.0
pair_modify mix arithmetic
kspace_style pppm 1e-4
special_bonds amber
atom_modify sort 0 0
read_data lammps.data
timestep 1.0
dump 1 all custom 1 dump.lammpstrj id element xu yu zu
dump 2 all custom 1 dump.force id element fx fy fz
dump 3 all xyz 1 dump.xyz
dump_modify 1 element O H
dump_modify 2 element O H
thermo_style multi
thermo 1
fix 1 all nvt temp 300.0 300.0 70.0
mdi/engine

69
lib/mdi/.gitignore vendored Normal file
View File

@ -0,0 +1,69 @@
# Prerequisites
*.d
# Object files
*.o
*.ko
*.obj
*.elf
# Linker output
*.ilk
*.map
*.exp
# Precompiled Headers
*.gch
*.pch
# Libraries
*.lib
*.a
*.la
*.lo
# Shared objects (inc. Windows DLLs)
*.dll
*.so
*.so.*
*.dylib
# Executables
*.exe
*.out
*.app
*.i*86
*.x86_64
*.hex
# Debug files
*.dSYM/
*.su
*.idb
*.pdb
# Kernel Module Compile Results
*.mod*
*.cmd
.tmp_versions/
modules.order
Module.symvers
Mkfile.old
dkms.conf
# Emacs temporary files
*~
# MacOS
.DS_Store
# pip
*.egg-info
# MDI files
MDI_Library/
build/
libmdi.a
mdi.h
includelink
liblink

211
lib/mdi/Install.py Normal file
View File

@ -0,0 +1,211 @@
#!/usr/bin/env python
# install.py tool to do a generic build of a library
# soft linked to by many of the lib/Install.py files
# used to automate the steps described in the corresponding lib/README
from __future__ import print_function
import sys,os,subprocess
import glob
sys.path.append('..')
from install_helpers import checkmd5sum
# help message
help = """
Syntax from src dir: make lib-libname args="-m machine -e suffix"
Syntax from lib dir: python Install.py -m machine -e suffix
libname = name of lib dir (e.g. atc, h5md, meam, poems, etc)
specify -m and optionally -e, order does not matter
-m = peform a clean followed by "make -f Makefile.machine"
machine = suffix of a lib/Makefile.* file
-e = set EXTRAMAKE variable in Makefile.machine to Makefile.lammps.suffix
does not alter existing Makefile.machine
Examples:
make lib-poems args="-m serial" # build POEMS lib with same settings as in the serial Makefile in src
make lib-colvars args="-m mpi" # build USER-COLVARS lib with same settings as in the mpi Makefile in src
make lib-meam args="-m ifort" # build MEAM lib with custom Makefile.ifort (using Intel Fortran)
"""
# settings
version = "1.2.9"
url = "https://github.com/MolSSI-MDI/MDI_Library/archive/v%s.tar.gz" % version
# known checksums for different MDI versions. used to validate the download.
checksums = { \
'1.2.7' : '2f3177b30ccdbd6ae28ea3bdd5fed0db', \
'1.2.9' : 'ddfa46d6ee15b4e59cfd527ec7212184', \
}
# print error message or help
def error(str=None):
if not str: print(help)
else: print("ERROR",str)
sys.exit()
# expand to full path name
# process leading '~' or relative path
def fullpath(path):
return os.path.abspath(os.path.expanduser(path))
def which(program):
def is_exe(fpath):
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
fpath, fname = os.path.split(program)
if fpath:
if is_exe(program):
return program
else:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
exe_file = os.path.join(path, program)
if is_exe(exe_file):
return exe_file
return None
def geturl(url,fname):
success = False
if which('curl') != None:
cmd = 'curl -L -o "%s" %s' % (fname,url)
try:
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
success = True
except subprocess.CalledProcessError as e:
print("Calling curl failed with: %s" % e.output.decode('UTF-8'))
if not success and which('wget') != None:
cmd = 'wget -O "%s" %s' % (fname,url)
try:
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
success = True
except subprocess.CalledProcessError as e:
print("Calling wget failed with: %s" % e.output.decode('UTF-8'))
if not success:
error("Failed to download source code with 'curl' or 'wget'")
return
# parse args
args = sys.argv[1:]
nargs = len(args)
if nargs == 0: error()
machine = None
extraflag = 0
iarg = 0
while iarg < nargs:
if args[iarg] == "-m":
if iarg+2 > nargs: error()
machine = args[iarg+1]
iarg += 2
elif args[iarg] == "-e":
if iarg+2 > nargs: error()
extraflag = 1
suffix = args[iarg+1]
iarg += 2
else: error()
# set lib from working dir
cwd = os.getcwd()
lib = os.path.basename(cwd)
# download and unpack MDI_Library tarball
homepath = "."
homedir = "%s/MDI_Library" % homepath
print("Downloading MDI_Library ...")
mditar = "%s/v%s.tar.gz" % (homepath,version)
geturl(url, mditar)
# verify downloaded archive integrity via md5 checksum, if known.
if version in checksums:
if not checkmd5sum(checksums[version], mditar):
sys.exit("Checksum for MDI library does not match")
print("Unpacking MDI_Library tarball ...")
if os.path.exists("%s/v%s" % (homepath,version)):
cmd = 'rm -rf "%s/v%s"' % (homepath,version)
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
cmd = 'cd "%s"; tar -xzvf v%s.tar.gz' % (homepath,version)
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
os.remove("%s/v%s.tar.gz" % (homepath,version))
if os.path.basename(homedir) != version:
if os.path.exists(homedir):
cmd = 'rm -rf "%s"' % homedir
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
os.rename("%s/MDI_Library-%s" % (homepath,version),homedir)
# create Makefile.auto as copy of Makefile.machine
# reset EXTRAMAKE if requested
if not os.path.exists("Makefile.%s" % machine):
error("lib/%s/Makefile.%s does not exist" % (lib,machine))
lines = open("Makefile.%s" % machine,'r').readlines()
fp = open("Makefile.auto",'w')
has_extramake = False
for line in lines:
words = line.split()
if len(words) == 3 and words[0] == "EXTRAMAKE" and words[1] == '=':
has_extramake = True
if extraflag:
line = line.replace(words[2],"Makefile.lammps.%s" % suffix)
fp.write(line)
fp.close()
# make the library via Makefile.auto optionally with parallel make
try:
import multiprocessing
n_cpus = multiprocessing.cpu_count()
except:
n_cpus = 1
print("Building lib%s.so ..." % lib)
cmd = "make -f Makefile.auto clean; make -f Makefile.auto -j%d" % n_cpus
txt = subprocess.check_output(cmd,shell=True,stderr=subprocess.STDOUT)
print(txt.decode('UTF-8'))
# create 2 links in lib/mdi to MDI Library src dir
print("Creating links to MDI Library include and lib files")
if os.path.isfile("includelink") or os.path.islink("includelink"):
os.remove("includelink")
if os.path.isfile("liblink") or os.path.islink("liblink"):
os.remove("liblink")
os.symlink(os.path.join(homedir, 'MDI_Library'), 'includelink')
os.symlink(os.path.join(homepath, 'build', 'MDI_Library'), 'liblink')
# Append the -rpath option to Makefile.lammps
dir_path = os.path.dirname(os.path.realpath(__file__))
rpath_option = "-Wl,-rpath=" + str(dir_path) + "/liblink"
makefile_lammps = open(str(dir_path) + "/Makefile.lammps", "a")
makefile_lammps.write(str(rpath_option) + "\n")
makefile_lammps.close()
shared_files = glob.glob( os.path.join( homepath, "liblink", "lib%s.so*" % lib) )
if len(shared_files) > 0:
print("Build was successful")
else:
error("Build of lib/%s/lib%s.so was NOT successful" % (lib,lib))
if has_extramake and not os.path.exists("Makefile.lammps"):
print("lib/%s/Makefile.lammps was NOT created" % lib)

18
lib/mdi/Makefile.g++ Normal file
View File

@ -0,0 +1,18 @@
SHELL = /bin/sh
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.empty
# ------ MAKE PROCEDURE ------
lib: $(OBJ)
mkdir -p build
cd build; cmake -Dlibtype=SHARED -Dlanguage=CXX -D CMAKE_C_COMPILER=gcc -D CMAKE_CXX_COMPILER=g++ ../MDI_Library; make
@cp $(EXTRAMAKE) Makefile.lammps
# ------ CLEAN ------
clean:
-rm *.o *.h $(LIB)
-rm -r build

18
lib/mdi/Makefile.gcc Normal file
View File

@ -0,0 +1,18 @@
SHELL = /bin/sh
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.empty
# ------ MAKE PROCEDURE ------
lib: $(OBJ)
mkdir -p build
cd build; cmake -Dlibtype=SHARED -Dlanguage=C -D CMAKE_C_COMPILER=gcc ../MDI_Library; make
@cp $(EXTRAMAKE) Makefile.lammps
# ------ CLEAN ------
clean:
-rm *.o *.h $(LIB)
-rm -r build

18
lib/mdi/Makefile.icc Normal file
View File

@ -0,0 +1,18 @@
SHELL = /bin/sh
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.empty
# ------ MAKE PROCEDURE ------
lib: $(OBJ)
mkdir -p build
cd build; cmake -Dlibtype=SHARED -Dlanguage=C -D CMAKE_C_COMPILER=icc -D CMAKE_CXX_COMPILER=icc ../MDI_Library; make
@cp $(EXTRAMAKE) Makefile.lammps
# ------ CLEAN ------
clean:
-rm *.o *.h $(LIB)
-rm -r build

View File

@ -0,0 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
mdi_SYSINC =
mdi_SYSLIB =
mdi_SYSPATH =

18
lib/mdi/Makefile.mpi Normal file
View File

@ -0,0 +1,18 @@
SHELL = /bin/sh
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.empty
# ------ MAKE PROCEDURE ------
lib: $(OBJ)
mkdir -p build
cd build; cmake -Dlibtype=SHARED -Dlanguage=CXX -D CMAKE_C_COMPILER=mpicc -D CMAKE_CXX_COMPILER=mpicxx ../MDI_Library; make
@cp $(EXTRAMAKE) Makefile.lammps
# ------ CLEAN ------
clean:
-rm *.o *.h $(LIB)
-rm -r build

7
src/.gitignore vendored
View File

@ -27,6 +27,13 @@
/*_ssa.h
/*_ssa.cpp
/fix_mdi_engine.cpp
/fix_mdi_engine.h
/library_mdi.cpp
/library_mdi.h
/mdi_engine.cpp
/mdi_engine.h
/fix_brownian*.cpp
/fix_brownian*.h
/fix_propel_self.cpp

View File

@ -54,20 +54,21 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \
PACKUSER = user-adios user-atc user-awpmd user-brownian user-bocs user-cgdna \
user-cgsdk user-colvars user-diffraction user-dpd user-drude \
user-eff user-fep user-h5md user-intel user-lb user-manifold \
user-meamc user-mesodpd user-mesont user-mgpt user-misc \
user-mdi user-meamc user-mesodpd user-mesont user-mgpt user-misc \
user-mofff \user-molfile user-netcdf user-omp user-phonon \
user-pace user-plumed user-ptm user-qmmm user-qtb user-quip \
user-reaction user-reaxc user-scafacos user-smd user-smtbq \
user-sdpd user-sph user-tally user-uef user-vtk user-yaff
PACKLIB = compress gpu kim kokkos latte message mpiio mscg poems python voronoi \
user-adios user-atc user-awpmd user-colvars user-h5md user-lb user-mesont \
user-molfile user-netcdf user-pace user-plumed user-qmmm user-quip \
user-adios user-atc user-awpmd user-colvars user-h5md user-lb user-mdi \
user-mesont user-molfile user-netcdf user-pace user-plumed user-qmmm user-quip \
user-scafacos user-smd user-vtk
PACKSYS = compress mpiio python user-lb
PACKINT = gpu kokkos message poems user-atc user-awpmd user-colvars user-mesont
PACKINT = gpu kokkos message poems user-atc user-awpmd user-colvars user-mesont \
user-mdi
PACKEXT = kim latte mscg voronoi \
user-adios user-h5md user-molfile user-netcdf user-pace user-plumed \

71
src/USER-MDI/Install.sh Executable file
View File

@ -0,0 +1,71 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# enforce using portable C locale
LC_ALL=C
export LC_ALL
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files with no dependencies
for file in *.cpp *.h; do
test -f ${file} && action $file
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*mdi[^ \t]* //g' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/mdi/includelink |' ../Makefile.package
sed -i -e 's/[^ \t]*MDI[^ \t]* //' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-DLMP_USER_MDI |' ../Makefile.package
sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/mdi/liblink |' ../Makefile.package
sed -i -e 's|^PKG_LIB =[ \t]*|&-lmdi |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(mdi_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(mdi_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(mdi_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*mdi.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/mdi\/Makefile.lammps
' ../Makefile.package.settings
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*mdi[^ \t]* //g' ../Makefile.package
sed -i -e 's/[^ \t]*MDI[^ \t]* //' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*mdi.*$/d' ../Makefile.package.settings
fi
fi

24
src/USER-MDI/README Normal file
View File

@ -0,0 +1,24 @@
The USER-MDI package adds an mdi/engine command which enables LAMMPS
to operate as a MolSSI Driver Interface (MDI) engine, responding to
commands from an external MDI driver.
It uses the MDI Library, which is available at
https://github.com/MolSSI-MDI/MDI_Library. The MDI Library is
developed and maintained by Taylor Barnes (tbarnes1@vt.edu) at the
Molecular Sciences Software Institute (MolSSI).
For more information about using MDI with LAMMPS, see the LAMMPS
documentation for the mdi/engine command and the mdi/engine fix.
For general purpose information about MDI, see the MDI Library
documentation:
https://molssi-mdi.github.io/MDI_Library/html/index.html
The MDI Library is required in order to use this package, and can be
built using the Install.py file in lib/mdi. For example:
python Install.py -m gcc
python Install.py -m icc
When using CMake to build LAMMPS the CMake configuration can use
either a preinstalled MDI library or download and compile a
private copy.

View File

@ -0,0 +1,953 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Taylor Barnes (MolSSI)
MolSSI Driver Interface (MDI) support for LAMMPS
------------------------------------------------------------------------- */
#include "fix_mdi_engine.h"
#include "library_mdi.h"
#include "atom.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "irregular.h"
#include "library.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include "verlet.h"
#include <limits>
#include <string.h>
enum { NONE, REAL, METAL }; // LAMMPS units which MDI supports
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixMDIEngine::FixMDIEngine(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), id_pe(nullptr), pe(nullptr), id_ke(nullptr), ke(nullptr)
{
if (narg != 3) error->all(FLERR, "Illegal fix mdi command");
// The 2 atomic-scale units LAMMPS has are:
// real: coords = Ang, eng = Kcal/mole, force = Kcal/mole/Ang
// metal: coords = Ang, eng = eV, force = eV/Ang
lmpunits = NONE;
if (strcmp(update->unit_style, "real") == 0) lmpunits = REAL;
if (strcmp(update->unit_style, "metal") == 0) lmpunits = METAL;
if (lmpunits == NONE) error->all(FLERR, "MDI requires real or metal units");
// MDI setup
most_recent_init = 0;
exit_flag = false;
local_exit_flag = false;
target_command = new char[MDI_COMMAND_LENGTH + 1];
command = new char[MDI_COMMAND_LENGTH + 1];
current_node = new char[MDI_COMMAND_LENGTH];
target_node = new char[MDI_COMMAND_LENGTH];
strncpy(target_node, "\0", MDI_COMMAND_LENGTH);
strncpy(current_node, "@DEFAULT", MDI_COMMAND_LENGTH);
// register the execute_command function with MDI
MDI_Set_execute_command_func(lammps_execute_mdi_command, this);
// accept a communicator to the driver
// master = 1 for proc 0, otherwise 0
master = (comm->me == 0) ? 1 : 0;
MDI_Accept_communicator(&driver_socket);
if (driver_socket <= 0) error->all(FLERR, "Unable to connect to driver");
// create computes for KE and PE
id_pe = utils::strdup(std::string(id) + "_pe");
modify->add_compute(fmt::format("{} all pe", id_pe));
id_ke = utils::strdup(std::string(id) + "_ke");
modify->add_compute(fmt::format("{} all ke", id_ke));
// irregular class and data structs used by MDI
irregular = new Irregular(lmp);
add_force = nullptr;
}
/* ---------------------------------------------------------------------- */
FixMDIEngine::~FixMDIEngine()
{
delete[] target_command;
delete[] command;
delete[] current_node;
delete[] target_node;
modify->delete_compute(id_pe);
modify->delete_compute(id_ke);
delete irregular;
memory->destroy(add_force);
}
/* ---------------------------------------------------------------------- */
int FixMDIEngine::setmask()
{
int mask = 0;
mask |= POST_INTEGRATE;
mask |= POST_FORCE;
mask |= MIN_PRE_FORCE;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::exchange_forces()
{
double **f = atom->f;
const int *const mask = atom->mask;
const int nlocal = atom->nlocal;
// add forces from the driver
for (int i = 0; i < nlocal; ++i) {
if (mask[i] & groupbit) {
f[i][0] += add_force[3 * (atom->tag[i] - 1) + 0];
f[i][1] += add_force[3 * (atom->tag[i] - 1) + 1];
f[i][2] += add_force[3 * (atom->tag[i] - 1) + 2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::init()
{
// confirm that two required computes are still available
int icompute_pe = modify->find_compute(id_pe);
if (icompute_pe < 0) error->all(FLERR, "Potential energy ID for fix mdi/engine does not exist");
int icompute_ke = modify->find_compute(id_ke);
if (icompute_pe < 0) error->all(FLERR, "Kinetic energy ID for fix mdi/engine does not exist");
pe = modify->compute[icompute_pe];
ke = modify->compute[icompute_ke];
// one-time allocation of add_force array
if (!add_force) {
int64_t ncoords = 3 * atom->natoms;
memory->create(add_force, ncoords, "mdi/engine:add_force");
for (int64_t i = 0; i < ncoords; i++) add_force[i] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::min_setup(int vflag)
{
engine_mode("@FORCES");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::post_integrate()
{
engine_mode("@COORDS");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::min_pre_force(int vflag)
{
engine_mode("@COORDS");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::min_post_force(int vflag)
{
engine_mode("@FORCES");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::post_force(int vflag)
{
if (most_recent_init == 1)
engine_mode("@FORCES");
else if (most_recent_init == 2)
engine_mode("@FORCES");
}
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
// rest of file processes and responds to MDI driver commands
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
process a single command from driver
---------------------------------------------------------------------- */
int FixMDIEngine::execute_command(const char *command, MDI_Comm mdicomm)
{
// confirm this command is supported at this node
int command_exists = 1;
if (master) {
ierr = MDI_Check_command_exists(current_node, command, MDI_COMM_NULL, &command_exists);
}
if (ierr != 0) error->all(FLERR, "MDI: Unable to check whether current command is supported");
if (command_exists != 1)
error->all(FLERR, "MDI: Received a command that is unsupported at current node");
// respond to any driver command
if (strcmp(command, ">NATOMS") == 0) {
ierr = MDI_Recv((char *) &atom->natoms, 1, MDI_INT, mdicomm);
if (ierr != 0) error->all(FLERR, "MDI: Unable to receive number of atoms from driver");
MPI_Bcast(&atom->natoms, 1, MPI_INT, 0, world);
} else if (strcmp(command, "<NATOMS") == 0) {
int64_t mdi_natoms = atom->natoms;
ierr = MDI_Send((char *) &mdi_natoms, 1, MDI_INT64_T, mdicomm);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send number of atoms to driver");
} else if (strcmp(command, "<NTYPES") == 0) {
ierr = MDI_Send((char *) &atom->ntypes, 1, MDI_INT, mdicomm);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send number of atom types to driver");
} else if (strcmp(command, "<TYPES") == 0) {
send_types(error);
} else if (strcmp(command, "<LABELS") == 0) {
send_labels(error);
} else if (strcmp(command, "<MASSES") == 0) {
send_masses(error);
} else if (strcmp(command, "<CELL") == 0) {
send_cell(error);
} else if (strcmp(command, ">CELL") == 0) {
receive_cell(error);
} else if (strcmp(command, "<CELL_DISPL") == 0) {
send_celldispl(error);
} else if (strcmp(command, ">CELL_DISPL") == 0) {
receive_celldispl(error);
} else if (strcmp(command, ">COORDS") == 0) {
receive_coordinates(error);
} else if (strcmp(command, "<COORDS") == 0) {
send_coordinates(error);
} else if (strcmp(command, "<CHARGES") == 0) {
send_charges(error);
} else if (strcmp(command, "<ENERGY") == 0) {
send_energy(error);
} else if (strcmp(command, "<FORCES") == 0) {
send_forces(error);
// replace current forces with forces received from the driver
} else if (strcmp(command, ">FORCES") == 0) {
receive_forces(error, 0);
// add forces received from the driver to current forces
} else if (strcmp(command, ">+FORCES") == 0) {
receive_forces(error, 1);
// initialize new MD simulation or minimization
// return control to return to mdi/engine
} else if (strcmp(command, "@INIT_MD") == 0) {
if (most_recent_init != 0) error->all(FLERR, "MDI: MDI is already performing a simulation");
most_recent_init = 1;
local_exit_flag = true;
// initialize new energy minimization
// return control to return to mdi/engine
} else if (strcmp(command, "@INIT_OPTG") == 0) {
if (most_recent_init != 0) error->all(FLERR, "MDI: MDI is already performing a simulation");
most_recent_init = 2;
local_exit_flag = true;
} else if (strcmp(command, "@") == 0) {
strncpy(target_node, "\0", MDI_COMMAND_LENGTH);
local_exit_flag = true;
} else if (strcmp(command, "<@") == 0) {
ierr = MDI_Send(current_node, MDI_NAME_LENGTH, MDI_CHAR, mdicomm);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send node to driver");
} else if (strcmp(command, "<KE") == 0) {
send_ke(error);
} else if (strcmp(command, "<PE") == 0) {
send_pe(error);
} else if (strcmp(command, "@DEFAULT") == 0) {
most_recent_init = 0;
strncpy(target_node, "@DEFAULT", MDI_COMMAND_LENGTH);
local_exit_flag = true;
// are we in the middle of a geometry optimization?
if (most_recent_init == 2) {
// ensure that the energy and force tolerances are met
update->etol = std::numeric_limits<double>::max();
update->ftol = std::numeric_limits<double>::max();
// set the maximum number of force evaluations to 0
update->max_eval = 0;
}
} else if (strcmp(command, "@COORDS") == 0) {
strncpy(target_node, "@COORDS", MDI_COMMAND_LENGTH);
local_exit_flag = true;
} else if (strcmp(command, "@FORCES") == 0) {
strncpy(target_node, "@FORCES", MDI_COMMAND_LENGTH);
local_exit_flag = true;
} else if (strcmp(command, "EXIT") == 0) {
// exit the driver code
exit_flag = true;
// are we in the middle of a geometry optimization?
if (most_recent_init == 2) {
// ensure that the energy and force tolerances are met
update->etol = std::numeric_limits<double>::max();
update->ftol = std::numeric_limits<double>::max();
// set the maximum number of force evaluations to 0
update->max_eval = 0;
}
} else {
error->all(FLERR, "MDI: Unknown command from driver");
}
return 0;
}
/* ---------------------------------------------------------------------- */
char *FixMDIEngine::engine_mode(const char *node)
{
/*
if (screen)
fprintf(screen,"MDI ENGINE MODE: %i\n",node);
if (logfile)
fprintf(logfile,"MDI ENGINE MODE: %i\n",node);
*/
// do not process commands if engine and driver are not at same node
// target_node = node that driver has set via a @ command
// current_node = node that engine (LAMMPS) has set
strncpy(current_node, node, MDI_COMMAND_LENGTH);
if (strcmp(target_node, "\0") != 0 && strcmp(target_node, current_node) != 0)
local_exit_flag = true;
// respond to commands from the driver
while (not exit_flag and not local_exit_flag) {
// read the next command from the driver
// all procs call this, but only proc 0 receives the command
ierr = MDI_Recv_command(command, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to receive command from driver");
// broadcast command to the other MPI tasks
MPI_Bcast(command, MDI_COMMAND_LENGTH, MPI_CHAR, 0, world);
// execute the command
this->execute_command(command, driver_socket);
// check if the target node is something other than the current node
if (strcmp(target_node, "\0") != 0 && strcmp(target_node, current_node) != 0)
local_exit_flag = true;
}
// local exit occured so turn off local exit flag
local_exit_flag = false;
return command;
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::receive_coordinates(Error *error)
{
// get conversion factor to atomic units
double posconv;
// real: coords = Ang, eng = Kcal/mole, force = Kcal/mole/Ang
// metal: coords = Ang, eng = eV, force = eV/Ang
if (lmpunits == REAL) {
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
posconv = force->angstrom / angstrom_to_bohr;
} else if (lmpunits == METAL) {
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
posconv = force->angstrom / angstrom_to_bohr;
}
// create buffer to hold all coords
double *buffer;
buffer = new double[3 * atom->natoms];
ierr = MDI_Recv((char *) buffer, 3 * atom->natoms, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to receive coordinates from driver");
MPI_Bcast(buffer, 3 * atom->natoms, MPI_DOUBLE, 0, world);
// pick local atoms from the buffer
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
x[i][0] = buffer[3 * (atom->tag[i] - 1) + 0] * posconv;
x[i][1] = buffer[3 * (atom->tag[i] - 1) + 1] * posconv;
x[i][2] = buffer[3 * (atom->tag[i] - 1) + 2] * posconv;
}
// ensure atoms are in current box & update box via shrink-wrap
// has to be be done before invoking Irregular::migrate_atoms()
// since it requires atoms be inside simulation box
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
if (domain->triclinic) domain->lamda2x(atom->nlocal);
// move atoms to new processors via irregular() only needed if
// migrate_check() says an atom moves too far
if (domain->triclinic) domain->x2lamda(atom->nlocal);
if (irregular->migrate_check()) irregular->migrate_atoms();
if (domain->triclinic) domain->lamda2x(atom->nlocal);
delete[] buffer;
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_coordinates(Error *error)
{
// get conversion factor to atomic units
double posconv;
if (lmpunits == REAL) {
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
posconv = force->angstrom / angstrom_to_bohr;
} else if (lmpunits == METAL) {
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
posconv = force->angstrom / angstrom_to_bohr;
}
int64_t ncoords = 3 * atom->natoms;
double *coords;
double *coords_reduced;
memory->create(coords, ncoords, "mdi/engine:coords");
memory->create(coords_reduced, ncoords, "mdi/engine:coords_reduced");
// zero coords
for (int64_t icoord = 0; icoord < ncoords; icoord++) coords[icoord] = 0.0;
// copy local atoms into buffer at correct locations
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
coords[3 * (atom->tag[i] - 1) + 0] = x[i][0] / posconv;
coords[3 * (atom->tag[i] - 1) + 1] = x[i][1] / posconv;
coords[3 * (atom->tag[i] - 1) + 2] = x[i][2] / posconv;
}
MPI_Reduce(coords, coords_reduced, 3 * atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world);
ierr = MDI_Send((char *) coords_reduced, 3 * atom->natoms, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send coordinates to driver");
memory->destroy(coords);
memory->destroy(coords_reduced);
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_charges(Error *error)
{
double *charges;
double *charges_reduced;
memory->create(charges, atom->natoms, "mdi/engine:charges");
memory->create(charges_reduced, atom->natoms, "mdi/engine:charges_reduced");
// zero the charges array
for (int icharge = 0; icharge < atom->natoms; icharge++) charges[icharge] = 0.0;
// pick local atoms from the buffer
double *charge = atom->q;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) { charges[atom->tag[i] - 1] = charge[i]; }
MPI_Reduce(charges, charges_reduced, atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world);
ierr = MDI_Send((char *) charges_reduced, atom->natoms, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send charges to driver");
memory->destroy(charges);
memory->destroy(charges_reduced);
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_energy(Error *error)
{
// get conversion factor to atomic units
double energy_conv;
if (lmpunits == REAL) {
double kelvin_to_hartree;
MDI_Conversion_factor("kelvin_energy", "hartree", &kelvin_to_hartree);
energy_conv = kelvin_to_hartree / force->boltz;
} else if (lmpunits == METAL) {
double ev_to_hartree;
MDI_Conversion_factor("electron_volt", "hartree", &ev_to_hartree);
energy_conv = ev_to_hartree;
}
double kelvin_to_hartree;
MDI_Conversion_factor("kelvin_energy", "hartree", &kelvin_to_hartree);
double potential_energy = pe->compute_scalar();
double kinetic_energy = ke->compute_scalar();
double total_energy;
double *send_energy = &total_energy;
// convert the energy to atomic units
potential_energy *= energy_conv;
kinetic_energy *= energy_conv;
total_energy = potential_energy + kinetic_energy;
ierr = MDI_Send((char *) send_energy, 1, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send potential energy to driver");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_pe(Error *error)
{
// get conversion factor to atomic units
double energy_conv;
if (lmpunits == REAL) {
double kelvin_to_hartree;
MDI_Conversion_factor("kelvin_energy", "hartree", &kelvin_to_hartree);
energy_conv = kelvin_to_hartree / force->boltz;
} else if (lmpunits == METAL) {
double ev_to_hartree;
MDI_Conversion_factor("electron_volt", "hartree", &ev_to_hartree);
energy_conv = ev_to_hartree;
}
double potential_energy = pe->compute_scalar();
double *send_energy = &potential_energy;
// convert the energy to atomic units
potential_energy *= energy_conv;
ierr = MDI_Send((char *) send_energy, 1, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send potential energy to driver");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_ke(Error *error)
{
// get conversion factor to atomic units
double energy_conv;
if (lmpunits == REAL) {
double kelvin_to_hartree;
MDI_Conversion_factor("kelvin_energy", "hartree", &kelvin_to_hartree);
energy_conv = kelvin_to_hartree / force->boltz;
} else if (lmpunits == METAL) {
double ev_to_hartree;
MDI_Conversion_factor("electron_volt", "hartree", &ev_to_hartree);
energy_conv = ev_to_hartree;
}
double kinetic_energy = ke->compute_scalar();
double *send_energy = &kinetic_energy;
// convert the energy to atomic units
kinetic_energy *= energy_conv;
ierr = MDI_Send((char *) send_energy, 1, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send potential energy to driver");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_types(Error *error)
{
int *const type = atom->type;
ierr = MDI_Send((char *) type, atom->natoms, MDI_INT, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send atom types to driver");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_labels(Error *error)
{
char *labels = new char[atom->natoms * MDI_LABEL_LENGTH];
memset(labels, ' ', atom->natoms * MDI_LABEL_LENGTH);
for (int iatom = 0; iatom < atom->natoms; iatom++) {
std::string label = std::to_string(atom->type[iatom]);
int label_len = std::min(int(label.length()), MDI_LABEL_LENGTH);
strncpy(&labels[iatom * MDI_LABEL_LENGTH], label.c_str(), label_len);
}
ierr = MDI_Send(labels, atom->natoms * MDI_LABEL_LENGTH, MDI_CHAR, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send atom types to driver");
delete[] labels;
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_masses(Error *error)
{
double *const rmass = atom->rmass;
double *const mass = atom->mass;
int *const type = atom->type;
int nlocal = atom->nlocal;
double *mass_by_atom;
double *mass_by_atom_reduced;
memory->create(mass_by_atom, atom->natoms, "mdi/engine:mass_by_atom");
memory->create(mass_by_atom_reduced, atom->natoms, "mdi/engine:mass_by_atom_reduced");
for (int iatom = 0; iatom < atom->natoms; iatom++) { mass_by_atom[iatom] = 0.0; }
// determine the atomic masses
if (rmass) {
for (int iatom = 0; iatom < nlocal; iatom++) {
mass_by_atom[atom->tag[iatom] - 1] = rmass[iatom];
}
} else {
for (int iatom = 0; iatom < nlocal; iatom++) {
mass_by_atom[atom->tag[iatom] - 1] = mass[type[iatom]];
}
}
MPI_Reduce(mass_by_atom, mass_by_atom_reduced, atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world);
// send the atomic masses to the driver
ierr = MDI_Send((char *) mass_by_atom_reduced, atom->natoms, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send atom masses to driver");
memory->destroy(mass_by_atom);
memory->destroy(mass_by_atom_reduced);
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_forces(Error *error)
{
// get conversion factor to atomic units
double force_conv;
if (lmpunits == REAL) {
double kelvin_to_hartree;
double angstrom_to_bohr;
MDI_Conversion_factor("kelvin_energy", "hartree", &kelvin_to_hartree);
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
force_conv = (kelvin_to_hartree / force->boltz) * (force->angstrom / angstrom_to_bohr);
} else if (lmpunits == METAL) {
double ev_to_hartree;
double angstrom_to_bohr;
MDI_Conversion_factor("electron_volt", "hartree", &ev_to_hartree);
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
force_conv = ev_to_hartree / angstrom_to_bohr;
}
double *forces;
double *forces_reduced;
double *x_buf;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int64_t ncoords = 3 * atom->natoms;
memory->create(forces, ncoords, "mdi/engine:forces");
memory->create(forces_reduced, ncoords, "mdi/engine:forces_reduced");
x_buf = new double[3 * nlocal];
// zero the forces array
for (int iforce = 0; iforce < 3 * atom->natoms; iforce++) forces[iforce] = 0.0;
// if not at a node, calculate the forces
if (strcmp(current_node, "@DEFAULT") == 0) {
// certain fixes, such as shake, move the coordinates
// to ensure that the coordinates do not change, store a copy
double **x = atom->x;
for (int i = 0; i < nlocal; i++) {
x_buf[3 * i + 0] = x[i][0];
x_buf[3 * i + 1] = x[i][1];
x_buf[3 * i + 2] = x[i][2];
}
// calculate the forces
update->whichflag = 1; // 1 for dynamics
update->nsteps = 1;
lmp->init();
update->integrate->setup_minimal(1);
if (strcmp(current_node, "@DEFAULT") == 0) {
// restore the original set of coordinates
double **x_new = atom->x;
for (int i = 0; i < nlocal; i++) {
x_new[i][0] = x_buf[3 * i + 0];
x_new[i][1] = x_buf[3 * i + 1];
x_new[i][2] = x_buf[3 * i + 2];
}
}
}
// pick local atoms from the buffer
double **f = atom->f;
for (int i = 0; i < nlocal; i++) {
forces[3 * (atom->tag[i] - 1) + 0] = f[i][0] * force_conv;
forces[3 * (atom->tag[i] - 1) + 1] = f[i][1] * force_conv;
forces[3 * (atom->tag[i] - 1) + 2] = f[i][2] * force_conv;
}
// reduce the forces onto rank 0
MPI_Reduce(forces, forces_reduced, 3 * atom->natoms, MPI_DOUBLE, MPI_SUM, 0, world);
// send the forces through MDI
ierr = MDI_Send((char *) forces_reduced, 3 * atom->natoms, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send atom forces to driver");
memory->destroy(forces);
memory->destroy(forces_reduced);
delete[] x_buf;
}
/* ---------------------------------------------------------------------- */
// Receive forces from the driver
// mode = 0: replace current forces with forces from driver
// mode = 1: add forces from driver to current forces
void FixMDIEngine::receive_forces(Error *error, int mode)
{
// get conversion factor to atomic units
double force_conv;
if (lmpunits == REAL) {
double kelvin_to_hartree;
double angstrom_to_bohr;
MDI_Conversion_factor("kelvin_energy", "hartree", &kelvin_to_hartree);
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
force_conv = (kelvin_to_hartree / force->boltz) * (force->angstrom / angstrom_to_bohr);
} else if (lmpunits == METAL) {
double ev_to_hartree;
double angstrom_to_bohr;
MDI_Conversion_factor("electron_volt", "hartree", &ev_to_hartree);
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
force_conv = ev_to_hartree / angstrom_to_bohr;
}
int64_t ncoords = 3 * atom->natoms;
double *forces;
memory->create(forces, ncoords, "mdi/engine:forces");
ierr = MDI_Recv((char *) forces, 3 * atom->natoms, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to receive atom forces to driver");
MPI_Bcast(forces, 3 * atom->natoms, MPI_DOUBLE, 0, world);
// pick local atoms from the buffer
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (mode == 0) { // Replace
for (int i = 0; i < nlocal; i++) {
f[i][0] = forces[3 * (atom->tag[i] - 1) + 0] / force_conv;
f[i][1] = forces[3 * (atom->tag[i] - 1) + 1] / force_conv;
f[i][2] = forces[3 * (atom->tag[i] - 1) + 2] / force_conv;
}
} else {
for (int i = 0; i < nlocal; i++) {
f[i][0] += forces[3 * (atom->tag[i] - 1) + 0] / force_conv;
f[i][1] += forces[3 * (atom->tag[i] - 1) + 1] / force_conv;
f[i][2] += forces[3 * (atom->tag[i] - 1) + 2] / force_conv;
}
}
memory->destroy(forces);
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_cell(Error *error)
{
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
double celldata[9];
celldata[0] = domain->boxhi[0] - domain->boxlo[0];
celldata[1] = 0.0;
celldata[2] = 0.0;
celldata[3] = domain->xy;
celldata[4] = domain->boxhi[1] - domain->boxlo[1];
celldata[5] = 0.0;
celldata[6] = domain->xz;
celldata[7] = domain->yz;
celldata[8] = domain->boxhi[2] - domain->boxlo[2];
// convert the units to bohr
double unit_conv = force->angstrom * angstrom_to_bohr;
for (int icell = 0; icell < 9; icell++) { celldata[icell] *= unit_conv; }
ierr = MDI_Send((char *) celldata, 9, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send cell dimensions to driver");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::receive_cell(Error *error)
{
double celldata[9];
// receive the new cell vector from the driver
ierr = MDI_Recv((char *) celldata, 9, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send cell dimensions to driver");
MPI_Bcast(&celldata[0], 9, MPI_DOUBLE, 0, world);
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
double unit_conv = force->angstrom * angstrom_to_bohr;
for (int icell = 0; icell < 9; icell++) { celldata[icell] /= unit_conv; }
// ensure that the new cell vector is orthogonal
double small = std::numeric_limits<double>::min();
if (abs(celldata[1]) > small or abs(celldata[2]) > small or abs(celldata[3]) > small or
abs(celldata[5]) > small or abs(celldata[6]) > small or abs(celldata[7]) > small) {
error->all(FLERR,
"MDI: LAMMPS currently only supports the >CELL command for orthogonal cell vectors");
}
// set the new LAMMPS cell dimensions
// This only works for orthogonal cell vectors.
// Supporting the more general case would be possible,
// but considerably more complex.
domain->boxhi[0] = celldata[0] + domain->boxlo[0];
domain->boxhi[1] = celldata[4] + domain->boxlo[1];
domain->boxhi[2] = celldata[8] + domain->boxlo[2];
domain->xy = 0.0;
domain->xz = 0.0;
domain->yz = 0.0;
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::send_celldispl(Error *error)
{
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
double celldata[3];
celldata[0] = domain->boxlo[0];
celldata[1] = domain->boxlo[1];
celldata[2] = domain->boxlo[2];
// convert the units to bohr
double unit_conv = force->angstrom * angstrom_to_bohr;
for (int icell = 0; icell < 3; icell++) { celldata[icell] *= unit_conv; }
ierr = MDI_Send((char *) celldata, 3, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to send cell displacement to driver");
}
/* ---------------------------------------------------------------------- */
void FixMDIEngine::receive_celldispl(Error *error)
{
// receive the cell displacement from the driver
double celldata[3];
ierr = MDI_Recv((char *) celldata, 3, MDI_DOUBLE, driver_socket);
if (ierr != 0) error->all(FLERR, "MDI: Unable to receive cell displacement from driver");
MPI_Bcast(&celldata[0], 3, MPI_DOUBLE, 0, world);
double angstrom_to_bohr;
MDI_Conversion_factor("angstrom", "bohr", &angstrom_to_bohr);
double unit_conv = force->angstrom * angstrom_to_bohr;
double old_boxlo[3];
old_boxlo[0] = domain->boxlo[0];
old_boxlo[1] = domain->boxlo[1];
old_boxlo[2] = domain->boxlo[2];
// adjust the values of boxlo and boxhi for the new cell displacement vector
domain->boxlo[0] = celldata[0] / unit_conv;
domain->boxlo[1] = celldata[1] / unit_conv;
domain->boxlo[2] = celldata[2] / unit_conv;
domain->boxhi[0] += domain->boxlo[0] - old_boxlo[0];
domain->boxhi[1] += domain->boxlo[1] - old_boxlo[1];
domain->boxhi[2] += domain->boxlo[2] - old_boxlo[2];
}

View File

@ -0,0 +1,143 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(mdi/engine, FixMDIEngine);
// clang-format on
#else
#ifndef LMP_FIX_MDI_ENGINE_H
#define LMP_FIX_MDI_ENGINE_H
#include "fix.h"
#include "mdi.h"
namespace LAMMPS_NS {
class FixMDIEngine : public Fix {
public:
FixMDIEngine(class LAMMPS *, int, char **);
~FixMDIEngine();
int setmask();
void init();
int execute_command(const char *command, MDI_Comm driver_socket);
char *engine_mode(const char *node);
// receive and update forces
//void setup(int);
void min_setup(int);
void post_integrate();
void post_force(int);
void min_pre_force(int); //@COORDS
void min_post_force(int); //@FORCES
double *add_force; // stores forces added using +FORCE command
double potential_energy; // stores potential energy
double kinetic_energy; // stores kinetic energy
// current command
char *command;
protected:
void exchange_forces();
private:
int lmpunits; // REAL or METAL
int master, ierr;
int driver_socket;
int most_recent_init; // which MDI init command was most recently received?
// 0 - none
// 1 - MD
// 2 - OPTG
bool exit_flag;
bool local_exit_flag;
char *current_node;
char *target_node; // is the code supposed to advance to a particular node?
// 0 - none
// 1 - @COORDS (before pre-force calculation)
// 2 - @PRE-FORCES (before final force calculation)
// 3 - @FORCES (before time integration)
// -1 - after MD_INIT command
// -2 - after MD_INIT command followed by @PRE-FORCES (actually @INIT_OPTG?)
// command to be executed at the target node
char *target_command;
char *id_pe;
char *id_ke;
class Irregular *irregular;
class Minimize *minimizer;
class Compute *pe;
class Compute *ke;
void send_types(Error *);
void send_labels(Error *);
void send_masses(Error *);
void receive_coordinates(Error *);
void send_coordinates(Error *);
void send_charges(Error *);
void send_energy(Error *);
void send_forces(Error *);
void send_pe(Error *);
void send_ke(Error *);
void receive_forces(Error *, int);
void send_cell(Error *);
void receive_cell(Error *);
void send_celldispl(Error *);
void receive_celldispl(Error *);
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory.
E: Potential energy ID for fix mdi does not exist
Self-explanatory.
E: Cannot use MDI command without atom IDs
Self-explanatory.
E: MDI command requires consecutive atom IDs
Self-explanatory.
E: Unable to connect to driver
Self-explanatory.
E: Unable to ... driver
Self-explanatory.
E: Unknown command from driver
The driver sent a command that is not supported by the LAMMPS
interface. In some cases this might be because a nonsensical
command was sent (i.e. "SCF"). In other cases, the LAMMPS
interface might benefit from being expanded.
*/

View File

@ -0,0 +1,113 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
// ----------------------------------------------------------------------
// MolSSI Driver Interface functions
// ----------------------------------------------------------------------
#include "library_mdi.h"
// needed to enable MPI support
#define LAMMPS_LIB_MPI 1
#include "library.h"
#include "fix_mdi_engine.h"
#include <cstring>
using namespace LAMMPS_NS;
/** Initialize an instance of LAMMPS as an MDI plugin
*
\verbatim embed:rst
This function is called by the MolSSI Driver Interface library (MDI)
when LAMMPS is run as a plugin, and should not otherwise be used.
The function initializes MDI, then creates and initializes an instance
of LAMMPS. The command-line arguments ``argc`` and ``argv`` used to
initialize LAMMPS are recieved from MDI. The LAMMPS instance runs an
input file, which must include the ``mdi/engine`` command; when LAMMPS
executes this command, it will begin listening for commands from the
driver. The name of the input file is obtained from the ``-in``
command-line argument, which must be provided by the MDI driver.
\endverbatim
* \param command string buffer corresponding to the command to be executed
* \param comm MDI communicator that can be used to communicated with the driver.
* \param class_obj pointer to an instance of an mdi/engine fix cast to ``void *``.
* \return 0 on no error. */
int MDI_Plugin_init_lammps()
{
// initialize MDI
int mdi_argc;
char **mdi_argv;
if (MDI_Plugin_get_argc(&mdi_argc)) MPI_Abort(MPI_COMM_WORLD, 1);
if (MDI_Plugin_get_argv(&mdi_argv)) MPI_Abort(MPI_COMM_WORLD, 1);
if (MDI_Init(&mdi_argc, &mdi_argv)) MPI_Abort(MPI_COMM_WORLD, 1);
// get the MPI intra-communicator for this code
MPI_Comm mpi_world_comm = MPI_COMM_WORLD;
if (MDI_MPI_get_world_comm(&mpi_world_comm)) MPI_Abort(MPI_COMM_WORLD, 1);
// find the -in argument
int iarg = 0;
char *filename;
bool found_filename = false;
while (iarg < mdi_argc && !found_filename) {
if ((strcmp(mdi_argv[iarg], "-in") == 0) || (strcmp(mdi_argv[iarg], "-i") == 0)) {
if (iarg + 2 > mdi_argc) MPI_Abort(MPI_COMM_WORLD, 1);
filename = mdi_argv[iarg + 1];
found_filename = true;
// remove -in argument from the command list
mdi_argc -= 2;
for (int jarg = iarg; jarg < mdi_argc; jarg++) mdi_argv[jarg] = mdi_argv[jarg + 2];
}
iarg++;
}
if (!found_filename) MPI_Abort(MPI_COMM_WORLD, 1);
// create and run a LAMMPS instance
void *lmp = nullptr;
if (lammps_config_has_mpi_support() > 0)
lmp = lammps_open(mdi_argc, mdi_argv, mpi_world_comm, nullptr);
else
lmp = lammps_open_no_mpi(mdi_argc, mdi_argv, nullptr);
lammps_file(lmp, filename);
lammps_close(lmp);
return 0;
}
/** Execute an MDI command
*
\verbatim embed:rst
This function is called by the MolSSI Driver Interface Library (MDI)
when LAMMPS is run as a plugin, and should not otherwise be used.
The function executes a single command from an external MDI driver.
\endverbatim
* \param command string buffer corresponding to the command to be executed
* \param comm MDI communicator that can be used to communicated with the driver.
* \param class_obj pointer to an instance of an mdi/engine fix cast to ``void *``.
* \return 0 on no error, 1 on error. */
int lammps_execute_mdi_command(const char *command, MDI_Comm comm, void *class_obj)
{
FixMDIEngine *mdi_fix = (FixMDIEngine *) class_obj;
return mdi_fix->execute_command(command, comm);
}

View File

@ -0,0 +1,26 @@
/* -*- c -*- ------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LAMMPS_LIBRARY_MDI_H
#define LAMMPS_LIBRARY_MDI_H
/* C style library calls to LAMMPS when a LAMMPS shared library is
* used as a plugin through MolSSI Driver Interface (MDI). */
#include <mdi.h>
extern "C" {
int MDI_Plugin_init_lammps();
int lammps_execute_mdi_command(const char *, MDI_Comm, void *);
}
#endif

354
src/USER-MDI/mdi_engine.cpp Normal file
View File

@ -0,0 +1,354 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Taylor Barnes (MolSSI)
MolSSI Driver Interface (MDI) support for LAMMPS
------------------------------------------------------------------------- */
#include "mdi_engine.h"
#include "atom.h"
#include "error.h"
#include "fix_mdi_engine.h"
#include "force.h"
#include "mdi.h"
#include "min.h"
#include "minimize.h"
#include "modify.h"
#include "output.h"
#include "timer.h"
#include "update.h"
#include "verlet.h"
#include <limits>
#include <string.h>
using namespace LAMMPS_NS;
/* ----------------------------------------------------------------------
trigger LAMMPS to start acting as an MDI engine
endlessly loop over receiving commands from driver and responding
much of the logic for this is in FixMDIEngine
when EXIT command is received, mdi/engine command exits
---------------------------------------------------------------------- */
void MDIEngine::command(int narg, char **arg)
{
// list of nodes and commands that a MDI-compliant MD code should support
// default node and its commands
MDI_Register_node("@DEFAULT");
MDI_Register_command("@DEFAULT", "<@");
MDI_Register_command("@DEFAULT", "<CELL");
MDI_Register_command("@DEFAULT", "<CELL_DISPL");
MDI_Register_command("@DEFAULT", "<CHARGES");
MDI_Register_command("@DEFAULT", "<COORDS");
MDI_Register_command("@DEFAULT", "<LABELS");
MDI_Register_command("@DEFAULT", "<MASSES");
MDI_Register_command("@DEFAULT", "<NATOMS");
MDI_Register_command("@DEFAULT", "<TYPES");
MDI_Register_command("@DEFAULT", ">CELL");
MDI_Register_command("@DEFAULT", ">CELL_DISPL");
MDI_Register_command("@DEFAULT", ">COORDS");
MDI_Register_command("@DEFAULT", "@INIT_MD");
MDI_Register_command("@DEFAULT", "@INIT_OPTG");
MDI_Register_command("@DEFAULT", "EXIT");
// node for setting up and running a dynamics simulation
MDI_Register_node("@INIT_MD");
MDI_Register_command("@INIT_MD", "<@");
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", "<ENERGY");
MDI_Register_command("@INIT_MD", "<FORCES");
MDI_Register_command("@INIT_MD", "<KE");
MDI_Register_command("@INIT_MD", "<LABELS");
MDI_Register_command("@INIT_MD", "<MASSES");
MDI_Register_command("@INIT_MD", "<NATOMS");
MDI_Register_command("@INIT_MD", "<PE");
MDI_Register_command("@INIT_MD", "<TYPES");
MDI_Register_command("@INIT_MD", ">CELL");
MDI_Register_command("@INIT_MD", ">CELL_DISPL");
MDI_Register_command("@INIT_MD", ">COORDS");
MDI_Register_command("@INIT_MD", ">FORCES");
MDI_Register_command("@INIT_MD", ">+FORCES");
MDI_Register_command("@INIT_MD", "@");
MDI_Register_command("@INIT_MD", "@COORDS");
MDI_Register_command("@INIT_MD", "@DEFAULT");
MDI_Register_command("@INIT_MD", "@FORCES");
MDI_Register_command("@INIT_MD", "EXIT");
// node for setting up and running a minimization
MDI_Register_node("@INIT_OPTG");
MDI_Register_command("@INIT_OPTG", "<@");
MDI_Register_command("@INIT_OPTG", "<CELL");
MDI_Register_command("@INIT_OPTG", "<CELL_DISPL");
MDI_Register_command("@INIT_OPTG", "<CHARGES");
MDI_Register_command("@INIT_OPTG", "<COORDS");
MDI_Register_command("@INIT_OPTG", "<ENERGY");
MDI_Register_command("@INIT_OPTG", "<FORCES");
MDI_Register_command("@INIT_OPTG", "<KE");
MDI_Register_command("@INIT_OPTG", "<LABELS");
MDI_Register_command("@INIT_OPTG", "<MASSES");
MDI_Register_command("@INIT_OPTG", "<NATOMS");
MDI_Register_command("@INIT_OPTG", "<PE");
MDI_Register_command("@INIT_OPTG", "<TYPES");
MDI_Register_command("@INIT_OPTG", ">CELL");
MDI_Register_command("@INIT_OPTG", ">CELL_DISPL");
MDI_Register_command("@INIT_OPTG", ">COORDS");
MDI_Register_command("@INIT_OPTG", ">FORCES");
MDI_Register_command("@INIT_OPTG", ">+FORCES");
MDI_Register_command("@INIT_OPTG", "@");
MDI_Register_command("@INIT_OPTG", "@COORDS");
MDI_Register_command("@INIT_OPTG", "@DEFAULT");
MDI_Register_command("@INIT_OPTG", "@FORCES");
MDI_Register_command("@INIT_OPTG", "EXIT");
// node at POST_FORCE location in timestep
MDI_Register_node("@FORCES");
MDI_Register_callback("@FORCES", ">FORCES");
MDI_Register_callback("@FORCES", ">+FORCES");
MDI_Register_command("@FORCES", "<@");
MDI_Register_command("@FORCES", "<CELL");
MDI_Register_command("@FORCES", "<CELL_DISPL");
MDI_Register_command("@FORCES", "<CHARGES");
MDI_Register_command("@FORCES", "<COORDS");
MDI_Register_command("@FORCES", "<ENERGY");
MDI_Register_command("@FORCES", "<FORCES");
MDI_Register_command("@FORCES", "<KE");
MDI_Register_command("@FORCES", "<LABELS");
MDI_Register_command("@FORCES", "<MASSES");
MDI_Register_command("@FORCES", "<NATOMS");
MDI_Register_command("@FORCES", "<PE");
MDI_Register_command("@FORCES", "<TYPES");
MDI_Register_command("@FORCES", ">CELL");
MDI_Register_command("@FORCES", ">CELL_DISPL");
MDI_Register_command("@FORCES", ">COORDS");
MDI_Register_command("@FORCES", ">FORCES");
MDI_Register_command("@FORCES", ">+FORCES");
MDI_Register_command("@FORCES", "@");
MDI_Register_command("@FORCES", "@COORDS");
MDI_Register_command("@FORCES", "@DEFAULT");
MDI_Register_command("@FORCES", "@FORCES");
MDI_Register_command("@FORCES", "EXIT");
// node at POST_INTEGRATE location in timestep
MDI_Register_node("@COORDS");
MDI_Register_command("@COORDS", "<@");
MDI_Register_command("@COORDS", "<CELL");
MDI_Register_command("@COORDS", "<CELL_DISPL");
MDI_Register_command("@COORDS", "<CHARGES");
MDI_Register_command("@COORDS", "<COORDS");
MDI_Register_command("@COORDS", "<ENERGY");
MDI_Register_command("@COORDS", "<FORCES");
MDI_Register_command("@COORDS", "<KE");
MDI_Register_command("@COORDS", "<LABELS");
MDI_Register_command("@COORDS", "<MASSES");
MDI_Register_command("@COORDS", "<NATOMS");
MDI_Register_command("@COORDS", "<PE");
MDI_Register_command("@COORDS", "<TYPES");
MDI_Register_command("@COORDS", ">CELL");
MDI_Register_command("@COORDS", ">CELL_DISPL");
MDI_Register_command("@COORDS", ">COORDS");
MDI_Register_command("@COORDS", ">FORCES");
MDI_Register_command("@COORDS", ">+FORCES");
MDI_Register_command("@COORDS", "@");
MDI_Register_command("@COORDS", "@COORDS");
MDI_Register_command("@COORDS", "@DEFAULT");
MDI_Register_command("@COORDS", "@FORCES");
MDI_Register_command("@COORDS", "EXIT");
// if the mdi/engine fix is not already present, add it now
int ifix = modify->find_fix_by_style("mdi/engine");
bool added_mdi_engine_fix = false;
if (ifix < 0) {
modify->add_fix("MDI_ENGINE_INTERNAL all mdi/engine");
added_mdi_engine_fix = true;
}
// identify the mdi_engine fix
ifix = modify->find_fix_by_style("mdi/engine");
mdi_fix = static_cast<FixMDIEngine *>(modify->fix[ifix]);
// check that LAMMPS is setup as a compatible MDI engine
if (narg > 0) error->all(FLERR, "Illegal mdi/engine command");
if (atom->tag_enable == 0) error->all(FLERR, "Cannot use mdi/engine without atom IDs");
if (atom->tag_consecutive() == 0) error->all(FLERR, "mdi/engine requires consecutive atom IDs");
// endless engine loop, responding to driver commands
char *command;
while (1) {
// mdi/engine command only recognizes three nodes
// DEFAULT, INIT_MD, INIT_OPTG
command = mdi_fix->engine_mode("@DEFAULT");
// MDI commands for dynamics or minimization
if (strcmp(command, "@INIT_MD") == 0) {
command = mdi_md();
if (strcmp(command, "EXIT")) break;
} else if (strcmp(command, "@INIT_OPTG") == 0) {
command = mdi_optg();
if (strcmp(command, "EXIT")) break;
} else if (strcmp(command, "EXIT") == 0) {
break;
} else
error->all(FLERR,
fmt::format("MDI node exited with "
"invalid command: {}",
command));
}
// remove mdi/engine fix that mdi/engine instantiated
if (added_mdi_engine_fix) modify->delete_fix("MDI_ENGINE_INTERNAL");
}
/* ----------------------------------------------------------------------
run an MD simulation under control of driver
---------------------------------------------------------------------- */
char *MDIEngine::mdi_md()
{
// initialize an MD simulation
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;
lmp->init();
// engine is now at @INIT_MD node
char *command = nullptr;
command = mdi_fix->engine_mode("@INIT_MD");
if (strcmp(command, "@DEFAULT") == 0 || strcmp(command, "EXIT") == 0) return command;
// setup the MD simulation
update->integrate->setup(1);
command = mdi_fix->engine_mode("@FORCES");
if (strcmp(command, "@DEFAULT") == 0 || strcmp(command, "EXIT") == 0) return command;
// run MD one step at a time
while (1) {
update->whichflag = 1;
timer->init_timeout();
update->nsteps += 1;
update->laststep += 1;
update->endstep = update->laststep;
output->next = update->ntimestep + 1;
// single MD timestep
update->integrate->run(1);
// done with MD if driver sends @DEFAULT or EXIT
command = mdi_fix->command;
if (strcmp(command, "@DEFAULT") == 0 || strcmp(command, "EXIT") == 0) return command;
}
return nullptr;
}
/* ----------------------------------------------------------------------
perform minimization under control of driver
---------------------------------------------------------------------- */
char *MDIEngine::mdi_optg()
{
// initialize an energy minization
Minimize *minimizer = new Minimize(lmp);
// setup the minimizer in a way that ensures optimization
// will continue until MDI driver exits
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();
update->whichflag = 2;
update->beginstep = update->firststep = update->ntimestep;
update->endstep = update->laststep = update->firststep + update->nsteps;
lmp->init();
// engine is now at @INIT_OPTG node
char *command = nullptr;
command = mdi_fix->engine_mode("@INIT_OPTG");
if (strcmp(command, "@DEFAULT") == 0 || strcmp(command, "EXIT") == 0) return command;
// setup the minimization
update->minimize->setup();
// get new command
command = mdi_fix->command;
if (strcmp(command, "@DEFAULT") == 0 || strcmp(command, "EXIT") == 0) return command;
// Start a minimization, which is configured to run (essentially)
// infinite steps. When the driver sends the EXIT command,
// the minimizer's energy and force tolerances are set to
// extremely large values, causing the minimization to end.
update->minimize->iterate(update->nsteps);
// return if driver sends @DEFAULT or EXIT
command = mdi_fix->command;
if (strcmp(command, "@DEFAULT") == 0 || strcmp(command, "EXIT") == 0) return command;
error->all(FLERR,
fmt::format("MDI reached end of OPTG simulation "
"with invalid command: {}",
command));
return nullptr;
}

53
src/USER-MDI/mdi_engine.h Normal file
View File

@ -0,0 +1,53 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMMAND_CLASS
// clang-format off
CommandStyle(mdi/engine, MDIEngine);
// clang-format on
#else
#ifndef LMP_MDI_ENGINE_H
#define LMP_MDI_ENGINE_H
#include "command.h"
namespace LAMMPS_NS {
class MDIEngine : public Command {
public:
MDIEngine(LAMMPS *lmp) : Command(lmp) {}
virtual ~MDIEngine() {}
void command(int, char **);
private:
class FixMDIEngine *mdi_fix;
char *mdi_md();
char *mdi_optg();
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
*/

View File

@ -673,7 +673,8 @@ void lammps_extract_box(void *handle, double *boxlo, double *boxhi,
This function sets the simulation box dimensions (upper and lower bounds
and tilt factors) from the provided data and then re-initializes the box
information and all derived settings.
information and all derived settings. It may only be called before atoms
are created.
\endverbatim
*
@ -692,12 +693,16 @@ void lammps_reset_box(void *handle, double *boxlo, double *boxhi,
BEGIN_CAPTURE
{
// error if box does not exist
if ((lmp->domain->box_exist == 0)
&& (lmp->comm->me == 0)) {
lmp->error->warning(FLERR,"Calling lammps_reset_box without a box");
if (lmp->atom->natoms > 0)
lmp->error->all(FLERR,"Calling lammps_reset_box not supported when atoms exist");
// warn and do nothing if no box exists
if (lmp->domain->box_exist == 0) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,"Ignoring call to lammps_reset_box without a box");
return;
}
domain->boxlo[0] = boxlo[0];
domain->boxlo[1] = boxlo[1];
domain->boxlo[2] = boxlo[2];

View File

@ -12,10 +12,11 @@
------------------------------------------------------------------------- */
#include "lammps.h"
#include "input.h"
#include <mpi.h>
#include <cstdlib>
#include <mpi.h>
#if defined(LAMMPS_TRAP_FPE) && defined(_GNU_SOURCE)
#include <fenv.h>
@ -25,6 +26,41 @@
#include "exceptions.h"
#endif
// import real or dummy calls to MolSSI Driver Interface library
#if defined(LMP_USER_MDI)
// true interface to MDI
#include <mdi.h>
#else
// dummy interface to MDI
// needed for compiling when MDI is not installed
typedef int MDI_Comm;
static int MDI_Init(int *argc, char ***argv)
{
return 0;
}
static int MDI_Initialized(int *flag)
{
return 0;
}
static int MDI_MPI_get_world_comm(void *world_comm)
{
return 0;
}
static int MDI_Plugin_get_argc(int *argc)
{
return 0;
}
static int MDI_Plugin_get_argv(char ***argv)
{
return 0;
}
#endif
using namespace LAMMPS_NS;
/* ----------------------------------------------------------------------
@ -33,11 +69,24 @@ using namespace LAMMPS_NS;
int main(int argc, char **argv)
{
MPI_Init(&argc,&argv);
MPI_Init(&argc, &argv);
// enable trapping selected floating point exceptions.
// this uses GNU extensions and is only tested on Linux
// therefore we make it depend on -D_GNU_SOURCE, too.
// initialize MDI or MDI dummy interface
int mdi_flag;
if (MDI_Init(&argc, &argv)) MPI_Abort(MPI_COMM_WORLD, 1);
if (MDI_Initialized(&mdi_flag)) MPI_Abort(MPI_COMM_WORLD, 1);
// get the MPI communicator that spans all ranks running LAMMPS
// when using MDI, this may be a subset of MPI_COMM_WORLD
MPI_Comm lammps_comm = MPI_COMM_WORLD;
if (mdi_flag)
if (MDI_MPI_get_world_comm(&lammps_comm)) MPI_Abort(MPI_COMM_WORLD, 1);
// enable trapping selected floating point exceptions.
// this uses GNU extensions and is only tested on Linux
// therefore we make it depend on -D_GNU_SOURCE, too.
#if defined(LAMMPS_TRAP_FPE) && defined(_GNU_SOURCE)
fesetenv(FE_NOMASK_ENV);
@ -49,31 +98,31 @@ int main(int argc, char **argv)
#ifdef LAMMPS_EXCEPTIONS
try {
LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD);
LAMMPS *lammps = new LAMMPS(argc, argv, lammps_comm);
lammps->input->file();
delete lammps;
} catch(LAMMPSAbortException &ae) {
} catch (LAMMPSAbortException &ae) {
MPI_Abort(ae.universe, 1);
} catch(LAMMPSException &e) {
MPI_Barrier(MPI_COMM_WORLD);
} catch (LAMMPSException &e) {
MPI_Barrier(lammps_comm);
MPI_Finalize();
exit(1);
} catch(fmt::format_error &fe) {
fprintf(stderr,"fmt::format_error: %s\n", fe.what());
} catch (fmt::format_error &fe) {
fprintf(stderr, "fmt::format_error: %s\n", fe.what());
MPI_Abort(MPI_COMM_WORLD, 1);
exit(1);
}
#else
try {
LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD);
LAMMPS *lammps = new LAMMPS(argc, argv, lammps_comm);
lammps->input->file();
delete lammps;
} catch(fmt::format_error &fe) {
fprintf(stderr,"fmt::format_error: %s\n", fe.what());
} catch (fmt::format_error &fe) {
fprintf(stderr, "fmt::format_error: %s\n", fe.what());
MPI_Abort(MPI_COMM_WORLD, 1);
exit(1);
}
#endif
MPI_Barrier(MPI_COMM_WORLD);
MPI_Barrier(lammps_comm);
MPI_Finalize();
}

View File

@ -54,6 +54,7 @@ PLUMED_URL="https://github.com/plumed/plumed2/releases/download/v2.7.0/plumed-sr
PACELIB_URL="https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2021.4.9.tar.gz"
LATTE_URL="https://github.com/lanl/LATTE/archive/v1.2.2.tar.gz"
SCAFACOS_URL="https://github.com/scafacos/scafacos/releases/download/v1.0.1/scafacos-1.0.1.tar.gz"
MDI_URL="https://github.com/MolSSI-MDI/MDI_Library/archive/v1.2.9.tar.gz"
GTEST_FILENAME="gtest-1.10.0.tar.gz"
MATHJAX_FILENAME="mathjax-3.1.2.tar.gz"
@ -81,6 +82,7 @@ TARBALLS=(
PACELIB_URL
LATTE_URL
SCAFACOS_URL
MDI_URL
)
###############################################################################

View File

@ -161,6 +161,8 @@ TEST_F(LibraryProperties, box)
boxhi[1] = 7.3;
xy = 0.1;
if (!verbose) ::testing::internal::CaptureStdout();
// lammps_reset_box() may only be called without atoms
lammps_command(lmp, "delete_atoms group all bond yes");
lammps_reset_box(lmp, boxlo, boxhi, xy, yz, xz);
if (!verbose) ::testing::internal::GetCapturedStdout();
lammps_extract_box(lmp, boxlo, boxhi, &xy, &yz, &xz, pflags, &boxflag);