diff --git a/cmake/Modules/Packages/ML-IAP.cmake b/cmake/Modules/Packages/ML-IAP.cmake index 63f91ba8d3..0e711cde8f 100644 --- a/cmake/Modules/Packages/ML-IAP.cmake +++ b/cmake/Modules/Packages/ML-IAP.cmake @@ -25,16 +25,18 @@ if(MLIAP_ENABLE_PYTHON) endif() set(MLIAP_BINARY_DIR ${CMAKE_BINARY_DIR}/cython) - set(MLIAP_CYTHON_SRC ${LAMMPS_SOURCE_DIR}/ML-IAP/mliap_model_python_couple.pyx) - get_filename_component(MLIAP_CYTHON_BASE ${MLIAP_CYTHON_SRC} NAME_WE) + file(GLOB MLIAP_CYTHON_SRC ${LAMMPS_SOURCE_DIR}/ML-IAP/*.pyx) file(MAKE_DIRECTORY ${MLIAP_BINARY_DIR}) - add_custom_command(OUTPUT ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.h - COMMAND ${CMAKE_COMMAND} -E copy_if_different ${MLIAP_CYTHON_SRC} ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx - COMMAND ${Cythonize_EXECUTABLE} -3 ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx - WORKING_DIRECTORY ${MLIAP_BINARY_DIR} - MAIN_DEPENDENCY ${MLIAP_CYTHON_SRC} - COMMENT "Generating C++ sources with cythonize...") + foreach(MLIAP_CYTHON_FILE ${MLIAP_CYTHON_SRC}) + get_filename_component(MLIAP_CYTHON_BASE ${MLIAP_CYTHON_FILE} NAME_WE) + add_custom_command(OUTPUT ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.h + COMMAND ${CMAKE_COMMAND} -E copy_if_different ${MLIAP_CYTHON_FILE} ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx + COMMAND ${Cythonize_EXECUTABLE} -3 ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx + WORKING_DIRECTORY ${MLIAP_BINARY_DIR} + MAIN_DEPENDENCY ${MLIAP_CYTHON_FILE} + COMMENT "Generating C++ sources with cythonize...") + target_sources(lammps PRIVATE ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp) + endforeach() target_compile_definitions(lammps PRIVATE -DMLIAP_PYTHON) - target_sources(lammps PRIVATE ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp) target_include_directories(lammps PRIVATE ${MLIAP_BINARY_DIR}) endif() diff --git a/doc/src/Fortran.rst b/doc/src/Fortran.rst index dd848a812e..77ab447c7c 100644 --- a/doc/src/Fortran.rst +++ b/doc/src/Fortran.rst @@ -3,7 +3,9 @@ The ``LIBLAMMPS`` Fortran Module The ``LIBLAMMPS`` module provides an interface to call LAMMPS from a Fortran code. It is based on the LAMMPS C-library interface and -requires a Fortran 2003 compatible compiler to be compiled. +requires a Fortran 2003 compatible compiler to be compiled. It is +designed to be self-contained and not require any support functions +written in C, C++, or Fortran. While C libraries have a defined binary interface (ABI) and can thus be used from multiple compiler versions from different vendors for as long @@ -19,12 +21,20 @@ for a simple program using the Fortran interface would be: mpifort -o testlib.x lammps.f90 testlib.f90 -L. -llammps Please note, that the MPI compiler wrapper is only required when the -calling the library from an MPI parallel code. Please also note the -order of the source files: the ``lammps.f90`` file needs to be compiled -first, since it provides the ``LIBLAMMPS`` module that is imported by -the Fortran code using the interface. A working example code can be -found together with equivalent examples in C and C++ in the -``examples/COUPLE/simple`` folder of the LAMMPS distribution. +calling the library from an MPI parallel code. Otherwise, using the +fortran compiler (gfortran, ifort, flang, etc.) will suffice. It may be +necessary to link to additional libraries depending on how LAMMPS was +configured and whether the LAMMPS library :doc:`was compiled as a static +or shared library `. + +If the LAMMPS library itself has been compiled with MPI support, the +resulting executable will still be able to run LAMMPS in parallel with +``mpirun`` or equivalent. Please also note that the order of the source +files matters: the ``lammps.f90`` file needs to be compiled first, since +it provides the ``LIBLAMMPS`` module that is imported by the Fortran +code using the interface. A working example code can be found together +with equivalent examples in C and C++ in the ``examples/COUPLE/simple`` +folder of the LAMMPS distribution. .. versionadded:: 9Oct2020 @@ -49,17 +59,18 @@ found together with equivalent examples in C and C++ in the Creating or deleting a LAMMPS object ************************************ -With the Fortran interface the creation of a :cpp:class:`LAMMPS +With the Fortran interface, the creation of a :cpp:class:`LAMMPS ` instance is included in the constructor for creating the :f:func:`lammps` derived type. To import the definition of -that type and its type bound procedures you need to add a ``USE +that type and its type bound procedures, you need to add a ``USE LIBLAMMPS`` statement. Internally it will call either :cpp:func:`lammps_open_fortran` or :cpp:func:`lammps_open_no_mpi` from the C library API to create the class instance. All arguments are optional and :cpp:func:`lammps_mpi_init` will be called automatically, -if it is needed. Similarly, a possible call to :cpp:func:`lammps_finalize` -is integrated into the :f:func:`close` function and triggered with -the optional logical argument set to ``.true.``. Here is a simple example: +if it is needed. Similarly, a possible call to +:cpp:func:`lammps_mpi_finalize` is integrated into the :f:func:`close` +function and triggered with the optional logical argument set to +``.true.``. Here is a simple example: .. code-block:: fortran @@ -80,11 +91,11 @@ the optional logical argument set to ``.true.``. Here is a simple example: END PROGRAM testlib It is also possible to pass command line flags from Fortran to C/C++ and -thus make the resulting executable behave similar to the standalone -executable (it will ignore the `-in/-i` flag, though). This allows to -use the command line to configure accelerator and suffix settings, +thus make the resulting executable behave similarly to the standalone +executable (it will ignore the `-in/-i` flag, though). This allows +using the command line to configure accelerator and suffix settings, configure screen and logfile output, or to set index style variables -from the command line and more. Here is a correspondingly adapted +from the command line and more. Here is a correspondingly adapted version of the previous example: .. code-block:: fortran @@ -117,13 +128,13 @@ version of the previous example: -------------------- Executing LAMMPS commands -========================= +************************* Once a LAMMPS instance is created, it is possible to "drive" the LAMMPS -simulation by telling LAMMPS to read commands from a file, or pass +simulation by telling LAMMPS to read commands from a file or to pass individual or multiple commands from strings or lists of strings. This -is done similar to how it is implemented in the `C-library -` interface. Before handing off the calls to the +is done similarly to how it is implemented in the :doc:`C-library +interface `. Before handing off the calls to the C-library interface, the corresponding Fortran versions of the calls (:f:func:`file`, :f:func:`command`, :f:func:`commands_list`, and :f:func:`commands_string`) have to make a copy of the strings passed as @@ -165,6 +176,57 @@ Below is a small demonstration of the uses of the different functions: --------------- +Accessing system properties +*************************** + +The C-library interface allows the :doc:`extraction of different kinds +of information ` about the active simulation +instance and also - in some cases - to apply modifications to it. In +some cases, the C-library interface makes pointers to internal data +structures accessible, thus when accessing them from Fortran, special +care is needed to avoid data corruption and crashes. Thus please see +the documentation of the individual type bound procedures for details. + +Below is an example demonstrating some of the possible uses. + +.. code-block:: fortran + + PROGRAM testprop + USE LIBLAMMPS + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int64_t + TYPE(lammps) :: lmp + INTEGER(kind=8) :: natoms + REAL(c_double), POINTER :: dt + INTEGER(c_int64_t), POINTER :: ntimestep + REAL(kind=8) :: pe, ke + + lmp = lammps() + CALL lmp%file('in.sysinit') + natoms = INT(lmp%get_natoms(),8) + WRITE(6,'(A,I8,A)') 'Running a simulation with', natoms, ' atoms' + WRITE(6,'(I8,A,I8,A,I3,A)') lmp%extract_setting('nlocal'), ' local and', & + lmp%extract_setting('nghost'), ' ghost atom. ', & + lmp%extract_setting('ntypes'), ' atom types' + + CALL lmp%command('run 2 post no') + dt = lmp%extract_global('dt') + ntimestep = lmp%extract_global('ntimestep') + WRITE(6,'(A,I4,A,F4.1,A)') 'At step:', ntimestep, ' Changing timestep from', dt, ' to 0.5' + dt = 0.5 + CALL lmp%command('run 2 post no') + + WRITE(6,'(A,I4)') 'At step:', ntimestep + pe = lmp%get_thermo('pe') + ke = lmp%get_thermo('ke') + PRINT*, 'PE = ', pe + PRINT*, 'KE = ', ke + + CALL lmp%close(.TRUE.) + + END PROGRAM testprop + +--------------- + The ``LIBLAMMPS`` module API **************************** @@ -178,14 +240,24 @@ of the contents of the ``LIBLAMMPS`` Fortran interface to LAMMPS. class instance that any of the included calls are forwarded to. :f c_ptr handle: reference to the LAMMPS class - :f close: :f:func:`close` - :f version: :f:func:`version` - :f file: :f:func:`file` - :f command: :f:func:`command` - :f commands_list: :f:func:`commands_list` - :f commands_string: :f:func:`commands_string` + :f subroutine close: :f:func:`close` + :f subroutine error: :f:func:`error` + :f function version: :f:func:`version` + :f subroutine file: :f:func:`file` + :f subroutine command: :f:func:`command` + :f subroutine commands_list: :f:func:`commands_list` + :f subroutine commands_string: :f:func:`commands_string` + :f function get_natoms: :f:func:`get_natoms` + :f function get_thermo: :f:func:`get_thermo` + :f subroutine extract_box: :f:func:`extract_box` + :f subroutine reset_box: :f:func:`reset_box` + :f subroutine memory_usage: :f:func:`memory_usage` + :f function extract_setting: :f:func:`extract_setting` + :f function extract_global: :f:func:`extract_global` -.. f:function:: lammps(args[,comm]) +-------- + +.. f:function:: lammps([args][,comm]) This is the constructor for the Fortran class and will forward the arguments to a call to either :cpp:func:`lammps_open_fortran` @@ -198,10 +270,31 @@ of the contents of the ``LIBLAMMPS`` Fortran interface to LAMMPS. If *comm* is not provided, ``MPI_COMM_WORLD`` is assumed. For more details please see the documentation of :cpp:func:`lammps_open`. - :p character(len=*) args(*) [optional]: arguments as list of strings + :o character(len=\*) args(\*) [optional]: arguments as list of strings :o integer comm [optional]: MPI communicator :r lammps: an instance of the :f:type:`lammps` derived type + .. note:: + + The ``MPI_F08`` module, which defines Fortran 2008 bindings for MPI, + is not directly supported by this interface due to the complexities of + supporting both the ``MPI_F08`` and ``MPI`` modules at the same time. + However, you should be able to use the ``MPI_VAL`` member of the + ``MPI_comm`` derived type to access the integer value of the + communicator, such as in + + .. code-block:: Fortran + + PROGRAM testmpi + USE LIBLAMMPS + USE MPI_F08 + TYPE(lammps) :: lmp + lmp = lammps(MPI_COMM_SELF%MPI_VAL) + END PROGRAM testmpi + +Procedures Bound to the lammps Derived Type +=========================================== + .. f:subroutine:: close([finalize]) This method will close down the LAMMPS instance through calling @@ -211,6 +304,20 @@ of the contents of the ``LIBLAMMPS`` Fortran interface to LAMMPS. :o logical finalize [optional]: shut down the MPI environment of the LAMMPS library if true. +-------- + +.. f:subroutine:: error(error_type, error_text) + + This method is a wrapper around the :cpp:func:`lammps_error` function and will dispatch + an error through the LAMMPS Error class. + + .. versionadded:: TBD + + :p integer error_type: constant to select which Error class function to call + :p character(len=\*) error_text: error message + +-------- + .. f:function:: version() This method returns the numeric LAMMPS version like :cpp:func:`lammps_version` @@ -224,25 +331,243 @@ of the contents of the ``LIBLAMMPS`` Fortran interface to LAMMPS. This method will call :cpp:func:`lammps_file` to have LAMMPS read and process commands from a file. - :p character(len=*) filename: name of file with LAMMPS commands + :p character(len=\*) filename: name of file with LAMMPS commands + +-------- .. f:subroutine:: command(cmd) This method will call :cpp:func:`lammps_command` to have LAMMPS execute a single command. - :p character(len=*) cmd: single LAMMPS command + :p character(len=\*) cmd: single LAMMPS command + +-------- .. f:subroutine:: commands_list(cmds) This method will call :cpp:func:`lammps_commands_list` to have LAMMPS execute a list of input lines. - :p character(len=*) cmd(:): list of LAMMPS input lines + :p character(len=\*) cmd(:): list of LAMMPS input lines + +-------- .. f:subroutine:: commands_string(str) This method will call :cpp:func:`lammps_commands_string` to have LAMMPS execute a block of commands from a string. - :p character(len=*) str: LAMMPS input in string + :p character(len=\*) str: LAMMPS input in string + +-------- + +.. f:function:: get_natoms() + + This function will call :cpp:func:`lammps_get_natoms` and return the number + of atoms in the system. + + :r real(c_double): number of atoms + +-------- + +.. f:function:: get_thermo(name) + + This function will call :cpp:func:`lammps_get_thermo` and return the value + of the corresponding thermodynamic keyword. + + .. versionadded:: TBD + + :p character(len=\*) name: string with the name of the thermo keyword + :r real(c_double): value of the requested thermo property or `0.0_c_double` + +-------- + +.. f:subroutine:: extract_box([boxlo][, boxhi][, xy][, yz][, xz][, pflags][, boxflag]) + + This subroutine will call :cpp:func:`lammps_extract_box`. All + parameters are optional, though obviously at least one should be + present. The parameters *pflags* and *boxflag* are stored in LAMMPS + as integers, but should be declared as ``LOGICAL`` variables when + calling from Fortran. + + .. versionadded:: TBD + + :o real(c_double) boxlo [dimension(3),optional]: vector in which to store + lower-bounds of simulation box + :o real(c_double) boxhi [dimension(3),optional]: vector in which to store + upper-bounds of simulation box + :o real(c_double) xy [optional]: variable in which to store *xy* tilt factor + :o real(c_double) yz [optional]: variable in which to store *yz* tilt factor + :o real(c_double) xz [optional]: variable in which to store *xz* tilt factor + :o logical pflags [dimension(3),optional]: vector in which to store + periodicity flags (``.TRUE.`` means periodic in that dimension) + :o logical boxflag [optional]: variable in which to store boolean denoting + whether the box will change during a simulation + (``.TRUE.`` means box will change) + +.. note:: + + Note that a frequent use case of this function is to extract only one or + more of the options rather than all seven. For example, assuming "lmp" + represents a properly-initialized LAMMPS instance, the following code will + extract the periodic box settings into the variable "periodic": + + .. code-block:: Fortran + + ! code to start up + logical :: periodic(3) + ! code to initialize LAMMPS / run things / etc. + call lmp%extract_box(pflags = periodic) + +-------- + +.. f:subroutine:: reset_box(boxlo, boxhi, xy, yz, xz) + + This subroutine will call :cpp:func:`lammps_reset_box`. All parameters + are required. + + .. versionadded:: TBD + + :p real(c_double) boxlo [dimension(3)]: vector of three doubles containing + the lower box boundary + :p real(c_double) boxhi [dimension(3)]: vector of three doubles containing + the upper box boundary + :p real(c_double) xy: *x--y* tilt factor + :p real(c_double) yz: *y--z* tilt factor + :p real(c_double) xz: *x--z* tilt factor + +-------- + +.. f:subroutine:: memory_usage(meminfo) + + This subroutine will call :cpp:func:`lammps_memory_usage` and store the + result in the three-element array *meminfo*. + + .. versionadded:: TBD + + :p real(c_double) meminfo [dimension(3)]: vector of three doubles in which + to store memory usage data + +-------- + +.. f:function:: get_mpi_comm() + + This function returns a Fortran representation of the LAMMPS "world" + communicator. + + .. versionadded:: TBD + + :r integer: Fortran integer equivalent to the MPI communicator LAMMPS is + using + + .. note:: + + The C library interface currently returns type ``int`` instead of + type ``MPI_Fint``, which is the C type corresponding to Fortran + ``INTEGER`` types of the default kind. On most compilers, these + are the same anyway, but this interface exchanges values this way + to avoid warning messages. + + .. note:: + + The `MPI_F08` module, which defines Fortran 2008 bindings for MPI, + is not directly supported by this function. However, you should be + able to convert between the two using the `MPI_VAL` member of the + communicator. For example, + + .. code-block:: fortran + + USE MPI_F08 + USE LIBLAMMPS + TYPE (LAMMPS) :: lmp + TYPE (MPI_Comm) :: comm + ! ... [commands to set up LAMMPS/etc.] + comm%MPI_VAL = lmp%get_mpi_comm() + + should assign an `MPI_F08` communicator properly. + +-------- + +.. f:function:: extract_setting(keyword) + + Query LAMMPS about global settings. See the documentation for the + :cpp:func:`lammps_extract_setting` function from the C library. + + .. versionadded:: TBD + + :p character(len=\*) keyword: string containing the name of the thermo keyword + :r integer(c_int): value of the queried setting or :math:`-1` if unknown + +-------- + +.. f:function:: extract_global(name) + + This function calls :cpp:func:`lammps_extract_global` and returns + either a string or a pointer to internal global LAMMPS data, + depending on the data requested through *name*. + + .. versionadded:: TBD + + Note that this function actually does not return a value, but rather + associates the pointer on the left side of the assignment to point to + internal LAMMPS data (with the exception of string data, which are + copied and returned as ordinary Fortran strings). Pointers must be of + the correct data type to point to said data (typically + ``INTEGER(c_int)``, ``INTEGER(c_int64_t)``, or ``REAL(c_double)``) + and have compatible kind and rank. The pointer being associated with + LAMMPS data is type-, kind-, and rank-checked at run-time via an + overloaded assignment operator. The pointers returned by this + function are generally persistent; therefore it is not necessary to + call the function again, unless a :doc:`clear` command has been + issued, which wipes out and recreates the contents of the + :cpp:class:`LAMMPS ` class. + + For example, + + .. code-block:: fortran + + PROGRAM demo + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int64_t + USE LIBLAMMPS + TYPE(lammps) :: lmp + INTEGER(c_int), POINTER :: nlocal + INTEGER(c_int64_t), POINTER :: ntimestep + CHARACTER(LEN=10) :: units + REAL(c_double), POINTER :: dt + lmp = lammps() + ! other commands + nlocal = lmp%extract_global('nlocal') + ntimestep = lmp%extract_global('ntimestep') + dt = lmp%extract_global('dt') + units = lmp%extract_global('units') + ! more commands + lmp.close(.TRUE.) + END PROGRAM demo + + would extract the number of atoms on this processor, the current time step, + the size of the current time step, and the units being used into the + variables *nlocal*, *ntimestep*, *dt*, and *units*, respectively. + + .. note:: + + if this function returns a string, the string must have + length greater than or equal to the length of the string (not including the + terminal NULL character) that LAMMPS returns. If the variable's length is + too short, the string will be truncated. As usual in Fortran, strings + are padded with spaces at the end. + + :p character(len=\*) name: string with the name of the extracted property + :r polymorphic: pointer to LAMMPS data. The left-hand side of the assignment + should be either a string (if expecting string data) or a C-compatible + pointer (e.g., ``INTEGER (c_int), POINTER :: nlocal``) to the extracted + property. If expecting vector data, the pointer should have dimension ":". + +.. warning:: + + Modifying the data in the location pointed to by the returned pointer + may lead to inconsistent internal data and thus may cause failures or + crashes or bogus simulations. In general it is thus usually better + to use a LAMMPS input command that sets or changes these parameters. + Those will take care of all side effects and necessary updates of + settings derived from such settings. diff --git a/doc/src/Howto_spins.rst b/doc/src/Howto_spins.rst index 87e30a5f6f..8ca220f596 100644 --- a/doc/src/Howto_spins.rst +++ b/doc/src/Howto_spins.rst @@ -30,9 +30,11 @@ can be coupled to another Langevin thermostat applied to the atoms using :doc:`fix langevin ` in order to simulate thermostatted spin-lattice systems. -The magnetic Gilbert damping can also be applied using :doc:`fix langevin/spin `. It allows to either dissipate -the thermal energy of the Langevin thermostat, or to perform a -relaxation of the magnetic configuration toward an equilibrium state. +The magnetic damping can also be applied +using :doc:`fix langevin/spin `. +It allows to either dissipate the thermal energy of the Langevin +thermostat, or to perform a relaxation of the magnetic configuration +toward an equilibrium state. The command :doc:`fix setforce/spin ` allows to set the components of the magnetic precession vectors (while erasing and @@ -52,9 +54,11 @@ All the computed magnetic properties can be output by two main commands. The first one is :doc:`compute spin `, that enables to evaluate magnetic averaged quantities, such as the total magnetization of the system along x, y, or z, the spin temperature, or -the magnetic energy. The second command is :doc:`compute property/atom `. It enables to output all the -per atom magnetic quantities. Typically, the orientation of a given -magnetic spin, or the magnetic force acting on this spin. +the magnetic energy. The second command +is :doc:`compute property/atom `. +It enables to output all the per atom magnetic quantities. Typically, +the orientation of a given magnetic spin, or the magnetic force +acting on this spin. ---------- diff --git a/doc/src/Library_create.rst b/doc/src/Library_create.rst index 8043819891..8ccc2e80ce 100644 --- a/doc/src/Library_create.rst +++ b/doc/src/Library_create.rst @@ -11,6 +11,7 @@ This section documents the following functions: - :cpp:func:`lammps_mpi_finalize` - :cpp:func:`lammps_kokkos_finalize` - :cpp:func:`lammps_python_finalize` +- :cpp:func:`lammps_error` -------------------- @@ -115,3 +116,8 @@ calling program. .. doxygenfunction:: lammps_python_finalize :project: progguide + +----------------------- + +.. doxygenfunction:: lammps_error + :project: progguide diff --git a/doc/src/Library_properties.rst b/doc/src/Library_properties.rst index e023c78185..a5c9c79c64 100644 --- a/doc/src/Library_properties.rst +++ b/doc/src/Library_properties.rst @@ -15,21 +15,21 @@ This section documents the following functions: -------------------- -The library interface allows extraction of different kinds of -information about the active simulation instance and also -modifications to it. This enables combining of a LAMMPS simulation -with other processing and simulation methods computed by the calling -code, or by another code that is coupled to LAMMPS via the library -interface. In some cases the data returned is direct reference to the -original data inside LAMMPS, cast to a void pointer. In that case the -data needs to be cast to a suitable pointer for the calling program to -access it, and you may need to know the correct dimensions and -lengths. This also means you can directly change those value(s) from -the calling program, e.g. to modify atom positions. Of course, this -should be done with care. When accessing per-atom data, please note -that this data is the per-processor **local** data and is indexed -accordingly. Per-atom data can change sizes and ordering at every -neighbor list rebuild or atom sort event as atoms migrate between +The library interface allows the extraction of different kinds of +information about the active simulation instance and also - in some +cases - to apply modifications to it. This enables combining of a +LAMMPS simulation with other processing and simulation methods computed +by the calling code, or by another code that is coupled to LAMMPS via +the library interface. In some cases the data returned is direct +reference to the original data inside LAMMPS, cast to a void pointer. +In that case the data needs to be cast to a suitable pointer for the +calling program to access it, and you may need to know the correct +dimensions and lengths. This also means you can directly change those +value(s) from the calling program, e.g. to modify atom positions. Of +course, this should be done with care. When accessing per-atom data, +please note that this data is the per-processor **local** data and is +indexed accordingly. Per-atom data can change sizes and ordering at +every neighbor list rebuild or atom sort event as atoms migrate between sub-domains and processors. .. code-block:: C diff --git a/doc/src/Library_utility.rst b/doc/src/Library_utility.rst index da64e3b8f0..9d16daae0b 100644 --- a/doc/src/Library_utility.rst +++ b/doc/src/Library_utility.rst @@ -19,6 +19,7 @@ functions. They do not directly call the LAMMPS library. - :cpp:func:`lammps_force_timeout` - :cpp:func:`lammps_has_error` - :cpp:func:`lammps_get_last_error_message` +- :cpp:func:`lammps_python_api_version` The :cpp:func:`lammps_free` function is a clean-up function to free memory that the library had allocated previously via other function @@ -100,3 +101,9 @@ where such memory buffers were allocated that require the use of .. doxygenfunction:: lammps_get_last_error_message :project: progguide + +----------------------- + +.. doxygenfunction:: lammps_python_api_version + :project: progguide + diff --git a/doc/src/fix_langevin_spin.rst b/doc/src/fix_langevin_spin.rst index 0926881eef..0ca291a732 100644 --- a/doc/src/fix_langevin_spin.rst +++ b/doc/src/fix_langevin_spin.rst @@ -40,7 +40,7 @@ the following stochastic differential equation: \times\left( \vec{\omega}_{i} \times\vec{s}_{i} \right) \right) with :math:`\lambda` the transverse damping, and :math:`\eta` a random vector. -This equation is referred to as the stochastic Landau-Lifshitz-Gilbert (sLLG) +This equation is referred to as the stochastic Landau-Lifshitz (sLL) equation. The components of :math:`\eta` are drawn from a Gaussian probability @@ -49,7 +49,7 @@ the external thermostat T (in K in metal units). More details about this implementation are reported in :ref:`(Tranchida) `. -Note: due to the form of the sLLG equation, this fix has to be defined just +Note: due to the form of the sLL equation, this fix has to be defined just before the nve/spin fix (and after all other magnetic fixes). As an example: diff --git a/doc/src/min_modify.rst b/doc/src/min_modify.rst index 6b416a9b6a..1b57a39064 100644 --- a/doc/src/min_modify.rst +++ b/doc/src/min_modify.rst @@ -24,7 +24,7 @@ Syntax inf = max force component across all 3-vectors max = max force norm across all 3-vectors *alpha_damp* value = damping - damping = fictitious Gilbert damping for spin minimization (adim) + damping = fictitious magnetic damping for spin minimization (adim) *discrete_factor* value = factor factor = discretization factor for adaptive spin timestep (adim) *integrator* value = *eulerimplicit* or *verlet* @@ -109,9 +109,9 @@ norm is replaced by the spin-torque norm. Keywords *alpha_damp* and *discrete_factor* only make sense when a :doc:`min_spin ` command is declared. -Keyword *alpha_damp* defines an analog of a magnetic Gilbert -damping. It defines a relaxation rate toward an equilibrium for -a given magnetic system. +Keyword *alpha_damp* defines an analog of a magnetic damping. +It defines a relaxation rate toward an equilibrium for a given +magnetic system. Keyword *discrete_factor* defines a discretization factor for the adaptive timestep used in the *spin* minimization. See :doc:`min_spin ` for more information about those diff --git a/doc/src/min_spin.rst b/doc/src/min_spin.rst index d9572f4463..9b6841ae8c 100644 --- a/doc/src/min_spin.rst +++ b/doc/src/min_spin.rst @@ -39,7 +39,7 @@ timestep, according to: \frac{d \vec{s}_{i}}{dt} = \lambda\, \vec{s}_{i} \times\left( \vec{\omega}_{i} \times\vec{s}_{i} \right) -with :math:`\lambda` a damping coefficient (similar to a Gilbert +with :math:`\lambda` a damping coefficient (similar to a magnetic damping). :math:`\lambda` can be defined by setting the *alpha_damp* keyword with the :doc:`min_modify ` command. diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst index 32737eaf50..3993a4701e 100644 --- a/doc/src/pair_mliap.rst +++ b/doc/src/pair_mliap.rst @@ -10,8 +10,8 @@ Syntax pair_style mliap ... keyword values ... -* two keyword/value pairs must be appended -* keyword = *model* or *descriptor* +* one or two keyword/value pairs must be appended +* keyword = *model* or *descriptor* or *unified* .. parsed-literal:: @@ -21,6 +21,9 @@ Syntax *descriptor* values = style filename style = *sna* or *so3* filename = name of file containing descriptor definitions + *unified* values = filename ghostneigh_flag + filename = name of file containing serialized unified Python object + ghostneigh_flag = 0/1 to turn off/on inclusion of ghost neighbors in neighbors list Examples """""""" @@ -30,20 +33,24 @@ Examples pair_style mliap model linear InP.mliap.model descriptor sna InP.mliap.descriptor pair_style mliap model quadratic W.mliap.model descriptor sna W.mliap.descriptor pair_style mliap model nn Si.nn.mliap.model descriptor so3 Si.nn.mliap.descriptor + pair_style mliap unified mliap_unified_lj_Ar.pkl 0 pair_coeff * * In P Description """"""""""" Pair style *mliap* provides a general interface to families of -machine-learning interatomic potentials. It allows separate definitions +machine-learning interatomic potentials. It allows separate definitions of the interatomic potential functional form (*model*) and the geometric -quantities that characterize the atomic positions (*descriptor*). By -defining *model* and *descriptor* separately, it is possible to use many -different models with a given descriptor, or many different descriptors -with a given model. The pair style currently supports only *sna* and *so3* -descriptor styles, but it is is straightforward to add new descriptor -styles. +quantities that characterize the atomic positions (*descriptor*). + +By defining *model* and *descriptor* separately, it is possible to use +many different models with a given descriptor, or many different +descriptors with a given model. The pair style currently supports only +*sna* and *so3* descriptor styles, but it is straightforward to add new +descriptor styles. By using the *unified* keyword, it is possible to +define a Python model that combines functionalities of both *model* and +*descriptor*. The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap `, including the linear, quadratic, and @@ -63,11 +70,11 @@ basis :ref:`(Bartok) ` and :ref:`(Zagaceta) `. The available models are *linear* and *nn*. The pair_style *mliap* command must be followed by two keywords *model* -and *descriptor* in either order. A single *pair_coeff* command is also -required. The first 2 arguments must be \* \* so as to span all LAMMPS -atom types. This is followed by a list of N arguments that specify the -mapping of MLIAP element names to LAMMPS atom types, where N is the -number of LAMMPS atom types. +and *descriptor* in either order, or the one keyword *unified*. A +single *pair_coeff* command is also required. The first 2 arguments +must be \* \* so as to span all LAMMPS atom types. This is followed by +a list of N arguments that specify the mapping of MLIAP element names to +LAMMPS atom types, where N is the number of LAMMPS atom types. The *model* keyword is followed by the model style. This is followed by a single argument specifying the model filename containing the @@ -115,14 +122,15 @@ The detail of *nn* module implementation can be found at :ref:`(Yanxon) ` +cutoff manually, such as in the following example. + + +.. code-block:: LAMMPS + + variable ninteractions equal 2 + variable cutdist equal 7.5 + variable skin equal 1.0 + variable commcut equal (${ninteractions}*${cutdist})+${skin} + neighbor ${skin} bin + comm_modify cutoff ${commcut} + + +.. note:: + + To load a model from memory + (i.e. an existing python object), call `lammps.mliap.load_unified(unified)` + from python, and then specify the filename as "EXISTS". When using LAMMPS via + the library mode, you will need to call `lammps.mliappy.activate_mliappy(lmp)` + on the active LAMMPS object before the pair style is defined. This call locates + and loads the mliap-specific python module that is built into LAMMPS. + Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index ab2272635f..f18ed74ff4 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1080,6 +1080,7 @@ filesystem filesystems Fily Fincham +Fint fingerprintconstants fingerprintsperelement Finnis @@ -1092,6 +1093,7 @@ flagHI flaglog flagN flagVF +flang fld floralwhite Florez @@ -1203,6 +1205,7 @@ gewald Gezelter Gflop gfortran +ghostneigh ghostwhite Giacomo gif diff --git a/examples/SPIN/test_problems/README b/examples/SPIN/test_problems/README index 17e0935a35..66de9fa88e 100644 --- a/examples/SPIN/test_problems/README +++ b/examples/SPIN/test_problems/README @@ -44,7 +44,9 @@ directory. energy versus temperature. Comparison to the Langevin function results (computed by the python script). Note: This example is a reworked version of a test problem - provided by Martin Kroger (ETHZ). + provided by Martin Kroger (ETHZ). + Note 2: Two versions of this test are deployed, one at low + damping (0.01) and one at large damping (1.0). - validation_nve: simulates a small assembly of magnetic atoms (54). The atoms are diff --git a/examples/SPIN/test_problems/run_all.sh b/examples/SPIN/test_problems/run_all.sh index fb682e3ef5..51f11dc6d6 100755 --- a/examples/SPIN/test_problems/run_all.sh +++ b/examples/SPIN/test_problems/run_all.sh @@ -12,13 +12,19 @@ cd validation_damped_precession/ rm res_lammps.dat res_llg.dat cd .. -# test 3: langevin, damping and Zeeman -cd validation_langevin_precession/ +# test 3: langevin, damping and Zeeman, low damping (0.01) +cd validation_langevin_precession_d0.01/ ./run-test-prec.sh rm average_spin test-prec-spin.in res_lammps.dat res_langevin.dat cd .. -# test 4: NVE run, test Etot preservation +# test 4: langevin, damping and Zeeman, large damping (1.0) +cd validation_langevin_precession_d1.0/ +./run-test-prec.sh +rm average_spin test-prec-spin.in res_lammps.dat res_langevin.dat +cd .. + +# test 5: NVE run, test Etot preservation cd validation_nve/ ./run-test-nve.sh rm nve_spin_lattice.pdf res_lammps.dat diff --git a/examples/SPIN/test_problems/validation_langevin_precession/langevin.py b/examples/SPIN/test_problems/validation_langevin_precession_d0.01/langevin.py similarity index 100% rename from examples/SPIN/test_problems/validation_langevin_precession/langevin.py rename to examples/SPIN/test_problems/validation_langevin_precession_d0.01/langevin.py diff --git a/examples/SPIN/test_problems/validation_langevin_precession/plot_precession.py b/examples/SPIN/test_problems/validation_langevin_precession_d0.01/plot_precession.py similarity index 100% rename from examples/SPIN/test_problems/validation_langevin_precession/plot_precession.py rename to examples/SPIN/test_problems/validation_langevin_precession_d0.01/plot_precession.py diff --git a/examples/SPIN/test_problems/validation_langevin_precession/run-test-prec.sh b/examples/SPIN/test_problems/validation_langevin_precession_d0.01/run-test-prec.sh similarity index 100% rename from examples/SPIN/test_problems/validation_langevin_precession/run-test-prec.sh rename to examples/SPIN/test_problems/validation_langevin_precession_d0.01/run-test-prec.sh diff --git a/examples/SPIN/test_problems/validation_langevin_precession/test-prec-spin.template b/examples/SPIN/test_problems/validation_langevin_precession_d0.01/test-prec-spin.template similarity index 100% rename from examples/SPIN/test_problems/validation_langevin_precession/test-prec-spin.template rename to examples/SPIN/test_problems/validation_langevin_precession_d0.01/test-prec-spin.template diff --git a/examples/SPIN/test_problems/validation_langevin_precession_d1.0/langevin.py b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/langevin.py new file mode 100755 index 0000000000..0e9b634eeb --- /dev/null +++ b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/langevin.py @@ -0,0 +1,22 @@ +#!/usr/bin/env python3 + +import numpy as np, pylab, tkinter +import matplotlib.pyplot as plt +import mpmath as mp + +mub=5.78901e-5 # Bohr magneton (eV/T) +kb=8.617333262145e-5 # Boltzman constant (eV/K) +g=2.0 # Lande factor (adim) + +Hz=10.0 # mag. field (T) + +#Definition of the Langevin function +def func(t): + return mp.coth(g*mub*Hz/(kb*t))-1.0/(g*mub*Hz/(kb*t)) + +npoints=200 +ti=0.01 +tf=20.0 +for i in range (0,npoints): + temp=ti+i*(tf-ti)/npoints + print('%lf %lf %lf' % (temp,func(temp),-g*mub*Hz*func(temp))) diff --git a/examples/SPIN/test_problems/validation_langevin_precession_d1.0/plot_precession.py b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/plot_precession.py new file mode 100755 index 0000000000..0141af56a0 --- /dev/null +++ b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/plot_precession.py @@ -0,0 +1,37 @@ +#!/usr/bin/env python3 + +#Program fitting the exchange interaction +#Model curve: Bethe-Slater function +import numpy as np, pylab, tkinter +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit +from decimal import * +import sys, string, os + + +argv = sys.argv +if len(argv) != 3: + print("Syntax: ./plot_precession.py res_lammps.dat res_langevin.dat") + sys.exit() + +lammps_file = sys.argv[1] +langevin_file = sys.argv[2] + +T_lmp,S_lmp,E_lmp = np.loadtxt(lammps_file, skiprows=0, usecols=(0,2,3),unpack=True) +T_lan,S_lan,E_lan = np.loadtxt(langevin_file, skiprows=0, usecols=(0,1,2),unpack=True) + +plt.figure() +plt.subplot(211) +plt.ylabel('') +plt.plot(T_lmp, S_lmp, 'b-', label='LAMMPS') +plt.plot(T_lan, S_lan, 'r--', label='Langevin') + +plt.subplot(212) +plt.ylabel('E (in eV)') +plt.plot(T_lmp, E_lmp, 'b-', label='LAMMPS') +plt.plot(T_lan, E_lan, 'r--', label='Langevin') + +plt.xlabel('T (in K)') +pylab.xlim([0,20.0]) +plt.legend() +plt.show() diff --git a/examples/SPIN/test_problems/validation_langevin_precession_d1.0/run-test-prec.sh b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/run-test-prec.sh new file mode 100755 index 0000000000..ecdbb2c156 --- /dev/null +++ b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/run-test-prec.sh @@ -0,0 +1,33 @@ +#!/bin/bash + +tempi=0.0 +tempf=20.0 + +rm res_*.dat + +# compute Lammps +N=20 +for (( i=0; i<$N; i++ )) +do + temp="$(echo "$tempi+$i*($tempf-$tempi)/$N" | bc -l)" + sed s/temperature/${temp}/g test-prec-spin.template > \ + test-prec-spin.in + + # test standard Lammps + ./../../../../src/lmp_serial -in test-prec-spin.in + + # test spin/kk with Kokkos Lammps + # mpirun -np 1 ../../../../src/lmp_kokkos_mpi_only \ + # -k on -sf kk -in test-prec-spin.in + + Hz="$(tail -n 1 average_spin | awk -F " " '{print $3}')" + sz="$(tail -n 1 average_spin | awk -F " " '{print $5}')" + en="$(tail -n 1 average_spin | awk -F " " '{print $6}')" + echo $temp $Hz $sz $en >> res_lammps.dat +done + +# compute Langevin +python3 langevin.py > res_langevin.dat + +# plot results +python3 plot_precession.py res_lammps.dat res_langevin.dat diff --git a/examples/SPIN/test_problems/validation_langevin_precession_d1.0/test-prec-spin.template b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/test-prec-spin.template new file mode 100644 index 0000000000..94e7a02fee --- /dev/null +++ b/examples/SPIN/test_problems/validation_langevin_precession_d1.0/test-prec-spin.template @@ -0,0 +1,48 @@ +#LAMMPS in.run + +units metal +atom_style spin +# atom_style spin/kk +atom_modify map array +boundary p p p + +# read_data singlespin.data + +lattice sc 3.0 +region box block 0.0 1.0 0.0 1.0 0.0 1.0 +create_box 1 box +create_atoms 1 box + +mass 1 1.0 +set type 1 spin 1.0 0.0 0.0 1.0 + +# defines a pair/style for neighbor list, but do not use it +pair_style spin/exchange 4.0 +pair_coeff * * exchange 1.0 0.0 0.0 1.0 + +group bead type 1 + +variable H equal 10.0 +variable Kan equal 0.0 +variable Temperature equal temperature +variable RUN equal 1000000 + +fix 1 all nve/spin lattice no +fix 2 all precession/spin zeeman ${H} 0.0 0.0 1.0 anisotropy ${Kan} 0.0 0.0 1.0 +fix_modify 2 energy yes +# fix 3 all langevin/spin ${Temperature} 0.01 12345 +fix 3 all langevin/spin ${Temperature} 1.0 12345 + +compute compute_spin all spin +compute outsp all property/atom spx spy spz sp +compute magsz all reduce ave c_outsp[3] + +thermo 50000 +thermo_style custom step time temp vol pe c_compute_spin[5] etotal + +variable magnetic_energy equal c_compute_spin[5] + +fix avespin all ave/time 1 ${RUN} ${RUN} v_Temperature v_H v_Kan c_magsz v_magnetic_energy file average_spin + +timestep 0.1 +run ${RUN} diff --git a/examples/mliap/README b/examples/mliap/README index c16fa218e0..070ce86bfd 100644 --- a/examples/mliap/README +++ b/examples/mliap/README @@ -1,6 +1,6 @@ -This directory contains multiple examples of -machine-learning potentials defined using the -MLIAP package in LAMMPS. The input files +This directory contains multiple examples of +machine-learning potentials defined using the +ML-IAP package in LAMMPS. The input files are described below. in.mliap.snap.Ta06A @@ -21,14 +21,14 @@ Run EME-SNAP, equivalent to examples/snap/in.snap.InP.JCPA2020 in.mliap.snap.compute --------------------- -Generate the A matrix, the gradients (w.r.t. coefficients) -of total potential energy, forces, and stress tensor for +Generate the A matrix, the gradients (w.r.t. coefficients) +of total potential energy, forces, and stress tensor for linear SNAP, equivalent to in.snap.compute in.mliap.quadratic.compute -------------------------- -Generate the A matrix, the gradients (w.r.t. coefficients) -of total potential energy, forces, and stress tensor for +Generate the A matrix, the gradients (w.r.t. coefficients) +of total potential energy, forces, and stress tensor for for quadratic SNAP, equivalent to in.snap.compute.quadratic in.mliap.pytorch.Ta06A @@ -41,8 +41,8 @@ This example can be run in two different ways: 1: Running a LAMMPS executable: in.mliap.pytorch.Ta06A First run ``python convert_mliap_Ta06A.py``. It creates -a PyTorch energy model that replicates the -SNAP Ta06A potential and saves it in the file +a PyTorch energy model that replicates the +SNAP Ta06A potential and saves it in the file "Ta06A.mliap.pytorch.model.pt". You can then run the example as follows @@ -52,7 +52,7 @@ You can then run the example as follows The resultant log.lammps output should be identical to that generated by in.mliap.snap.Ta06A. -If this fails, see the instructions for building the MLIAP package +If this fails, see the instructions for building the ML-IAP package with Python support enabled. Also, confirm that the LAMMPS Python embedded Python interpreter is working by running ../examples/in.python. @@ -62,7 +62,7 @@ working by running ../examples/in.python. Before testing this, ensure that the previous method (running a LAMMPS executable) works. -You can run the example in serial: +You can run the example in serial: `python mliap_pytorch_Ta06A.py` @@ -80,25 +80,25 @@ the script will exit with an error message. in.mliap.pytorch.relu1hidden ---------------------------- This example demonstrates a simple neural network potential -using PyTorch and SNAP descriptors. +using PyTorch and SNAP descriptors. `lmp -in in.mliap.pytorch.relu1hidden -echo both` -It was trained on just the energy component (no forces) of -the data used in the original SNAP Ta06A potential for -tantalum (Thompson, Swiler, Trott, Foiles, Tucker, +It was trained on just the energy component (no forces) of +the data used in the original SNAP Ta06A potential for +tantalum (Thompson, Swiler, Trott, Foiles, Tucker, J Comp Phys, 285, 316 (2015).). Because of the very small amount -of energy training data, it uses just 1 hidden layer with -a ReLU activation function. It is not expected to be -very accurate for forces. +of energy training data, it uses just 1 hidden layer with +a ReLU activation function. It is not expected to be +very accurate for forces. -NOTE: Unlike the previous example, this example uses -a pre-built PyTorch file `Ta06A.mliap.pytorch.model.pt`. +NOTE: Unlike the previous example, this example uses +a pre-built PyTorch file `Ta06A.mliap.pytorch.model.pt`. It is read using `torch.load`, which implicitly uses the Python `pickle` module. -This is known to be insecure. It is possible to construct malicious -pickle data that will execute arbitrary code during unpickling. Never -load data that could have come from an untrusted source, or that +This is known to be insecure. It is possible to construct malicious +pickle data that will execute arbitrary code during unpickling. Never +load data that could have come from an untrusted source, or that could have been tampered with. Only load data you trust. in.mliap.nn.Ta06A @@ -109,4 +109,54 @@ in.mliap.nn.cu ------------------------- Run a neural network potential for Cu, a combination of SNAP descriptors and the "nn" model style +in.mliap.unified.lj.Ar +----------------------- +This reproduces the output of examples/melt/in.melt, +but using the ``unified`` keyword and +the Python coupling to PyTorch. +This example can be run in two different ways: + +1: Running a LAMMPS executable: in.mliap.unified.lj.Ar + +You can run the example as follows + +`lmp -in in.mliap.unified.lj.Ar -echo both` + +The resultant thermo output should be identical to that generated +by in.melt. + +The input first runs some python code that creates an ML-IAP unified model +that replicates the ``lj/cut`` pair style with parameters as defined in +examples/melt/in/melt and saves it in the file "mliap_unified_lj_Ar.pkl". + +If this fails, see the instructions for building the ML-IAP package +with Python support enabled. Also, confirm that the LAMMPS Python +embedded Python interpreter is working by running ../examples/in.python +and that the LAMMPS python module is installed and working. + +2: Running a Python script: mliap_unified_lj_Ar.py + +Before testing this, ensure that the previous method +(running a LAMMPS executable) works. + +You can run the example in serial: + +`python mliap_unified_lj_Ar.py` + +or in parallel: + +`mpirun -np 4 python mliap_unified_lj_Ar.py` + +The resultant log.lammps output should be identical to that generated +by in.melt and in.mliap.unified.lj.Ar. + +in.mliap.so3.Ni_Mo +------------------ + +Example of linear model with SO3 descriptors for Ni/Mo + +in.mliap.so3.nn.Si +------------------ + +Example of NN model with SO3 descriptors for Si diff --git a/examples/mliap/in.mliap.unified.lj.Ar b/examples/mliap/in.mliap.unified.lj.Ar new file mode 100644 index 0000000000..190ea1288f --- /dev/null +++ b/examples/mliap/in.mliap.unified.lj.Ar @@ -0,0 +1,53 @@ + +# create pickle of unified model + +python create_pickle here """ +import lammps +import lammps.mliap + +from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ + +def create_pickle(): + unified = MLIAPUnifiedLJ(["Ar"]) + unified.pickle('mliap_unified_lj_Ar.pkl') +""" +python create_pickle invoke + +# 3d Lennard-Jones melt + +units lj +atom_style atomic + +lattice fcc 0.8442 +region box block 0 10 0 10 0 10 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 3.0 87287 loop geom + +pair_style mliap unified mliap_unified_lj_Ar.pkl 0 +pair_coeff * * Ar + +neighbor 0.3 bin +neigh_modify every 20 delay 0 check no + +fix 1 all nve + +#dump id all atom 50 dump.melt + +#dump 2 all image 25 image.*.jpg type type & +# axes yes 0.8 0.02 view 60 -30 +#dump_modify 2 pad 3 + +#dump 3 all movie 1 movie.mpg type type & +# axes yes 0.8 0.02 view 60 -30 +#dump_modify 3 pad 3 + +#dump 4 all custom 1 forces.xyz fx fy fz + +thermo 50 +run 250 + +# clean up +shell rm -f mliap_unified_lj_Ar.pkl diff --git a/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.1 b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.1 new file mode 100644 index 0000000000..1059b65538 --- /dev/null +++ b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.1 @@ -0,0 +1,107 @@ +LAMMPS (15 Sep 2022) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98) + using 1 OpenMP thread(s) per MPI task + +# create pickle of unified model + +python create_pickle here """ +import lammps +import lammps.mliap +from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ +def create_pickle(): + unified = MLIAPUnifiedLJ(["Ar"]) + unified.pickle('mliap_unified_lj_Ar.pkl') +""" +python create_pickle invoke + +# 3d Lennard-Jones melt + +units lj +atom_style atomic + +lattice fcc 0.8442 +Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962 +region box block 0 10 0 10 0 10 +create_box 1 box +Created orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962) + 1 by 1 by 1 MPI processor grid +create_atoms 1 box +Created 4000 atoms + using lattice units in orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962) + create_atoms CPU = 0.001 seconds +mass 1 1.0 + +velocity all create 3.0 87287 loop geom + +pair_style mliap unified mliap_unified_lj_Ar.pkl 0 +pair_coeff * * Ar + +neighbor 0.3 bin +neigh_modify every 20 delay 0 check no + +fix 1 all nve + +#dump id all atom 50 dump.melt + +#dump 2 all image 25 image.*.jpg type type # axes yes 0.8 0.02 view 60 -30 +#dump_modify 2 pad 3 + +#dump 3 all movie 1 movie.mpg type type # axes yes 0.8 0.02 view 60 -30 +#dump_modify 3 pad 3 + +#dump 4 all custom 1 forces.xyz fx fy fz + +thermo 50 +run 250 +Neighbor list info ... + update: every = 20 steps, delay = 0 steps, check = no + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 2.8 + ghost atom cutoff = 2.8 + binsize = 1.4, bins = 12 12 12 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair mliap, perpetual + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 12.7 | 12.7 | 12.7 Mbytes + Step Temp E_pair E_mol TotEng Press + 0 3 -6.7733681 0 -2.2744931 -3.7033504 + 50 1.6842865 -4.8082494 0 -2.2824513 5.5666131 + 100 1.6712577 -4.7875609 0 -2.281301 5.6613913 + 150 1.6444751 -4.7471034 0 -2.2810074 5.8614211 + 200 1.6471542 -4.7509053 0 -2.2807916 5.8805431 + 250 1.6645597 -4.7774327 0 -2.2812174 5.7526089 +Loop time of 3.48563 on 1 procs for 250 steps with 4000 atoms + +Performance: 30984.374 tau/day, 71.723 timesteps/s +79.5% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 3.3348 | 3.3348 | 3.3348 | 0.0 | 95.67 +Neigh | 0.13055 | 0.13055 | 0.13055 | 0.0 | 3.75 +Comm | 0.0089845 | 0.0089845 | 0.0089845 | 0.0 | 0.26 +Output | 0.00018194 | 0.00018194 | 0.00018194 | 0.0 | 0.01 +Modify | 0.0090026 | 0.0090026 | 0.0090026 | 0.0 | 0.26 +Other | | 0.002151 | | | 0.06 + +Nlocal: 4000 ave 4000 max 4000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 5506 ave 5506 max 5506 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +FullNghs: 303576 ave 303576 max 303576 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 303576 +Ave neighs/atom = 75.894 +Neighbor list builds = 12 +Dangerous builds not checked + +# clean up +shell rm -f mliap_unified_lj_Ar.pkl +Total wall time: 0:00:03 diff --git a/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.4 b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.4 new file mode 100644 index 0000000000..236e415ae6 --- /dev/null +++ b/examples/mliap/log.23Sep22.mliap.unified.lj.Ar.g++.4 @@ -0,0 +1,107 @@ +LAMMPS (15 Sep 2022) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98) + using 1 OpenMP thread(s) per MPI task + +# create pickle of unified model + +python create_pickle here """ +import lammps +import lammps.mliap +from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ +def create_pickle(): + unified = MLIAPUnifiedLJ(["Ar"]) + unified.pickle('mliap_unified_lj_Ar.pkl') +""" +python create_pickle invoke + +# 3d Lennard-Jones melt + +units lj +atom_style atomic + +lattice fcc 0.8442 +Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962 +region box block 0 10 0 10 0 10 +create_box 1 box +Created orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962) + 1 by 2 by 2 MPI processor grid +create_atoms 1 box +Created 4000 atoms + using lattice units in orthogonal box = (0 0 0) to (16.795962 16.795962 16.795962) + create_atoms CPU = 0.000 seconds +mass 1 1.0 + +velocity all create 3.0 87287 loop geom + +pair_style mliap unified mliap_unified_lj_Ar.pkl 0 +pair_coeff * * Ar + +neighbor 0.3 bin +neigh_modify every 20 delay 0 check no + +fix 1 all nve + +#dump id all atom 50 dump.melt + +#dump 2 all image 25 image.*.jpg type type # axes yes 0.8 0.02 view 60 -30 +#dump_modify 2 pad 3 + +#dump 3 all movie 1 movie.mpg type type # axes yes 0.8 0.02 view 60 -30 +#dump_modify 3 pad 3 + +#dump 4 all custom 1 forces.xyz fx fy fz + +thermo 50 +run 250 +Neighbor list info ... + update: every = 20 steps, delay = 0 steps, check = no + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 2.8 + ghost atom cutoff = 2.8 + binsize = 1.4, bins = 12 12 12 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair mliap, perpetual + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 5.774 | 5.774 | 5.774 Mbytes + Step Temp E_pair E_mol TotEng Press + 0 3 -6.7733681 0 -2.2744931 -3.7033504 + 50 1.6842865 -4.8082494 0 -2.2824513 5.5666131 + 100 1.6712577 -4.7875609 0 -2.281301 5.6613913 + 150 1.6444751 -4.7471034 0 -2.2810074 5.8614211 + 200 1.6471542 -4.7509053 0 -2.2807916 5.8805431 + 250 1.6645597 -4.7774327 0 -2.2812174 5.7526089 +Loop time of 1.18059 on 4 procs for 250 steps with 4000 atoms + +Performance: 91479.359 tau/day, 211.758 timesteps/s +78.3% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 1.0795 | 1.1007 | 1.1197 | 1.7 | 93.23 +Neigh | 0.035618 | 0.03627 | 0.037413 | 0.4 | 3.07 +Comm | 0.019611 | 0.038646 | 0.060389 | 9.4 | 3.27 +Output | 0.00013034 | 0.0001421 | 0.00017083 | 0.0 | 0.01 +Modify | 0.0037233 | 0.0037641 | 0.0038121 | 0.1 | 0.32 +Other | | 0.001086 | | | 0.09 + +Nlocal: 1000 ave 1008 max 987 min +Histogram: 1 0 0 0 0 0 1 0 1 1 +Nghost: 2711.25 ave 2728 max 2693 min +Histogram: 1 0 0 0 0 2 0 0 0 1 +Neighs: 0 ave 0 max 0 min +Histogram: 4 0 0 0 0 0 0 0 0 0 +FullNghs: 75894 ave 77414 max 74214 min +Histogram: 1 0 0 1 0 0 0 1 0 1 + +Total # of neighbors = 303576 +Ave neighs/atom = 75.894 +Neighbor list builds = 12 +Dangerous builds not checked + +# clean up +shell rm -f mliap_unified_lj_Ar.pkl +Total wall time: 0:00:01 diff --git a/examples/mliap/mliap_numpy_Ta06A.py b/examples/mliap/mliap_numpy_Ta06A.py new file mode 100644 index 0000000000..25c928cd89 --- /dev/null +++ b/examples/mliap/mliap_numpy_Ta06A.py @@ -0,0 +1,101 @@ +before_loading =\ +"""# Demonstrate MLIAP/PyTorch interface to linear SNAP potential + +# Initialize simulation + +variable nsteps index 100 +variable nrep equal 4 +variable a equal 3.316 +units metal + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable ny equal ${nrep} +variable nz equal ${nrep} + +boundary p p p + +lattice bcc $a +region box block 0 ${nx} 0 ${ny} 0 ${nz} +create_box 1 box +create_atoms 1 box + +mass 1 180.88 + +# choose potential + +# DATE: 2014-09-05 UNITS: metal CONTRIBUTOR: Aidan Thompson athomps@sandia.gov CITATION: Thompson, Swiler, Trott, Foiles and Tucker, arxiv.org, 1409.3880 (2014) + +# Definition of SNAP potential Ta_Cand06A +# Assumes 1 LAMMPS atom type + +variable zblcutinner equal 4 +variable zblcutouter equal 4.8 +variable zblz equal 73 + +# Specify hybrid with SNAP, ZBL + +pair_style hybrid/overlay & +zbl ${zblcutinner} ${zblcutouter} & +mliap model mliappy LATER & +descriptor sna Ta06A.mliap.descriptor +pair_coeff 1 1 zbl ${zblz} ${zblz} +pair_coeff * * mliap Ta +""" +after_loading =\ +""" + +# Setup output + +compute eatom all pe/atom +compute energy all reduce sum c_eatom + +compute satom all stress/atom NULL +compute str all reduce sum c_satom[1] c_satom[2] c_satom[3] +variable press equal (c_str[1]+c_str[2]+c_str[3])/(3*vol) + +thermo_style custom step temp epair c_energy etotal press v_press +thermo 10 +thermo_modify norm yes + +# Set up NVE run + +timestep 0.5e-3 +neighbor 1.0 bin +neigh_modify once no every 1 delay 0 check yes + +# Run MD + +velocity all create 300.0 4928459 loop geom +fix 1 all nve +run ${nsteps} +""" + +import numpy as np + +class LinearModel(): + def __init__(self,file): + coeffs = np.genfromtxt(file,skip_header=6) + self.bias = coeffs[0] + self.weights = coeffs[1:] + self.n_params = len(coeffs) + self.n_descriptors = len(self.weights) + self.n_elements = 1 + def __call__(self,elems,bispectrum,beta,energy): + energy[:] = bispectrum @ self.weights + self.bias + beta[:] = self.weights + + +mymodel = LinearModel("Ta06A.mliap.model") + +import lammps + +lmp = lammps.lammps(cmdargs=['-echo','both']) + +import lammps.mliap +lammps.mliap.activate_mliappy(lmp) +lmp.commands_string(before_loading) +lammps.mliap.load_model(mymodel) +lmp.commands_string(after_loading) + diff --git a/examples/mliap/mliap_unified_lj_Ar.py b/examples/mliap/mliap_unified_lj_Ar.py new file mode 100644 index 0000000000..2b404bd98e --- /dev/null +++ b/examples/mliap/mliap_unified_lj_Ar.py @@ -0,0 +1,65 @@ +# Demonstrate how to load a unified model from python. +# This is essentially the same as in.mliap.unified.lj.Ar +# except that python is the driving program, and lammps +# is in library mode. + +before_loading =\ +"""# 3d Lennard-Jones melt + +units lj +atom_style atomic + +lattice fcc 0.8442 +region box block 0 10 0 10 0 10 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 3.0 87287 loop geom +""" +after_loading =\ +""" + +pair_style mliap unified EXISTS +pair_coeff * * Ar + +neighbor 0.3 bin +neigh_modify every 20 delay 0 check no + +fix 1 all nve + +thermo 50 +run 250 +""" + +import lammps + +lmp = lammps.lammps(cmdargs=['-echo','both']) + +# Before defining the pair style, you must do the following: +import lammps.mliap +lammps.mliap.activate_mliappy(lmp) +# Otherwise, when running lammps in library mode, +# you will get an error + +# Setup the simulation to just before declaring +# the mliap unified pair style +lmp.commands_string(before_loading) + +# Define the model however you like. In this example +# we simply import the unified L-J example from mliap +from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ +unified = MLIAPUnifiedLJ(["Ar"]) + +# You can also load the model from a pickle file. +# import pickle +# with open('mliap_unified_lj_Ar.pkl', 'rb') as pfile: +# unified = pickle.load(pfile) + +# Connect the L-J model to the mliap unified pair style. +lammps.mliap.load_unified(unified) + + +# Run the simulation with the mliap unified pair style +# Use pre-loaded model by specifying model filename as "EXISTS" +lmp.commands_string(after_loading) diff --git a/fortran/README b/fortran/README index 653b6a966e..57d163e197 100644 --- a/fortran/README +++ b/fortran/README @@ -3,9 +3,9 @@ and allows the LAMMPS library interface to be invoked from Fortran codes. It requires a Fortran compiler that supports the Fortran 2003 standard. This interface is based on and supersedes the previous Fortran interfaces -in the examples/COUPLE/fortran* folders. But is fully supported by the +in the examples/COUPLE/fortran* folders, but is fully supported by the LAMMPS developers and included in the documentation and unit testing. Details on this Fortran interface and how to build programs using it -are in the manual in the doc/html/pg_fortran.html file. +are in the manual in the doc/html/Fortran.html file. diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index 144fd15652..98378c833a 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -19,117 +19,358 @@ ! Karl D. Hammond ! University of Missouri, 2012-2020 ! -! The Fortran module tries to follow the API of the C-library interface -! closely, but like the Python wrapper it employs an object oriented -! approach. To accommodate the object oriented approach, all exported -! subroutine and functions have to be implemented in Fortran to then -! call the interfaced C style functions with adapted calling conventions -! as needed. The C-library interfaced functions retain their names -! starting with "lammps_" while the Fortran versions start with "lmp_". +! The Fortran module tries to follow the API of the C library interface +! closely, but like the Python wrapper it employs an object-oriented +! approach. To accommodate the object-oriented approach, all exported +! subroutines and functions have to be implemented in Fortran and +! call the interfaced C-style functions with adapted calling conventions +! as needed. The C library interface functions retain their names +! starting with "lammps_", while the Fortran versions start with "lmp_". ! MODULE LIBLAMMPS USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_ptr, c_null_ptr, c_loc, & - c_int, c_char, c_null_char, c_double, c_size_t, c_f_pointer + c_int, c_int64_t, c_char, c_null_char, c_double, c_size_t, c_f_pointer + USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY : ERROR_UNIT IMPLICIT NONE PRIVATE - PUBLIC :: lammps + PUBLIC :: lammps, ASSIGNMENT(=) + + ! Data type constants for extracting data from global, atom, compute, and fix + ! + ! Must be kept in sync with the equivalent declarations in + ! src/library.h and python/lammps/constants.py + ! + ! NOT part of the API (the part the user sees) + INTEGER (c_int), PARAMETER :: & + LAMMPS_INT = 0, & ! 32-bit integer (array) + LAMMPS_INT_2D = 1, & ! two-dimensional 32-bit integer array + LAMMPS_DOUBLE = 2, & ! 64-bit double (array) + LAMMPS_DOUBLE_2D = 3, & ! two-dimensional 64-bit double array + LAMMPS_INT64 = 4, & ! 64-bit integer (array) + LAMMPS_INT64_2D = 5, & ! two-dimensional 64-bit integer array + LAMMPS_STRING = 6 ! C-String TYPE lammps TYPE(c_ptr) :: handle CONTAINS PROCEDURE :: close => lmp_close + PROCEDURE :: error => lmp_error PROCEDURE :: file => lmp_file PROCEDURE :: command => lmp_command PROCEDURE :: commands_list => lmp_commands_list PROCEDURE :: commands_string => lmp_commands_string - PROCEDURE :: version => lmp_version PROCEDURE :: get_natoms => lmp_get_natoms + PROCEDURE :: get_thermo => lmp_get_thermo + PROCEDURE :: extract_box => lmp_extract_box + PROCEDURE :: reset_box => lmp_reset_box + PROCEDURE :: memory_usage => lmp_memory_usage + PROCEDURE :: get_mpi_comm => lmp_get_mpi_comm + PROCEDURE :: extract_setting => lmp_extract_setting + PROCEDURE :: extract_global => lmp_extract_global + PROCEDURE :: version => lmp_version + PROCEDURE :: is_running => lmp_is_running END TYPE lammps INTERFACE lammps - MODULE PROCEDURE lmp_open + MODULE PROCEDURE lmp_open END INTERFACE lammps + ! Constants to use in working with lammps_data + ENUM, BIND(C) + ENUMERATOR :: DATA_INT, DATA_INT_1D, DATA_INT_2D + ENUMERATOR :: DATA_INT64, DATA_INT64_1D, DATA_INT64_2D + ENUMERATOR :: DATA_DOUBLE, DATA_DOUBLE_1D, DATA_DOUBLE_2D + ENUMERATOR :: DATA_STRING + END ENUM + + ! Derived type for receiving LAMMPS data (in lieu of the ability to type cast + ! pointers) + TYPE lammps_data + INTEGER(c_int) :: datatype + INTEGER(c_int), POINTER :: i32 + INTEGER(c_int), DIMENSION(:), POINTER :: i32_vec + INTEGER(c_int64_t), POINTER :: i64 + INTEGER(c_int64_t), DIMENSION(:), POINTER :: i64_vec + REAL(c_double), POINTER :: r64 + REAL(c_double), DIMENSION(:), POINTER :: r64_vec + CHARACTER(LEN=:), ALLOCATABLE :: str + END TYPE lammps_data + + ! This overloads the assignment operator (=) so that assignments of the + ! form + ! nlocal = extract_global('nlocal') + ! which are of the form "pointer to double = type(lammps_data)" result in + ! re-associating the pointer on the left with the appropriate piece of + ! LAMMPS data (after checking type-compatibility) + INTERFACE ASSIGNMENT(=) + MODULE PROCEDURE assign_int_to_lammps_data, assign_int64_to_lammps_data, & + assign_intvec_to_lammps_data, & + assign_double_to_lammps_data, assign_doublevec_to_lammps_data, & + assign_string_to_lammps_data + END INTERFACE + ! interface definitions for calling functions in library.cpp INTERFACE - FUNCTION lammps_open(argc, argv, comm) BIND(C, name='lammps_open_fortran') - IMPORT :: c_ptr, c_int - INTEGER(c_int), VALUE, INTENT(in) :: argc, comm - TYPE(c_ptr), DIMENSION(*), INTENT(in) :: argv - TYPE(c_ptr) :: lammps_open - END FUNCTION lammps_open + FUNCTION lammps_open(argc, argv, comm) BIND(C,name='lammps_open_fortran') + IMPORT :: c_ptr, c_int + IMPLICIT NONE + INTEGER(c_int), VALUE, INTENT(IN) :: argc, comm + TYPE(c_ptr), DIMENSION(*), INTENT(IN) :: argv + TYPE(c_ptr) :: lammps_open + END FUNCTION lammps_open - FUNCTION lammps_open_no_mpi(argc, argv, handle) BIND(C, name='lammps_open_no_mpi') - IMPORT :: c_ptr, c_int - INTEGER(c_int), VALUE, INTENT(in) :: argc - TYPE(c_ptr), DIMENSION(*), INTENT(in) :: argv - TYPE(c_ptr), VALUE, INTENT(in) :: handle - TYPE(c_ptr) :: lammps_open_no_mpi - END FUNCTION lammps_open_no_mpi + FUNCTION lammps_open_no_mpi(argc, argv, handle) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + INTEGER(c_int), VALUE, INTENT(IN) :: argc + TYPE(c_ptr), DIMENSION(*), INTENT(IN) :: argv + TYPE(c_ptr), VALUE, INTENT(IN) :: handle + TYPE(c_ptr) :: lammps_open_no_mpi + END FUNCTION lammps_open_no_mpi - SUBROUTINE lammps_close(handle) BIND(C, name='lammps_close') - IMPORT :: c_ptr - TYPE(c_ptr), VALUE :: handle - END SUBROUTINE lammps_close + SUBROUTINE lammps_close(handle) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + END SUBROUTINE lammps_close - SUBROUTINE lammps_mpi_init() BIND(C, name='lammps_mpi_init') - END SUBROUTINE lammps_mpi_init + SUBROUTINE lammps_mpi_init() BIND(C) + END SUBROUTINE lammps_mpi_init - SUBROUTINE lammps_mpi_finalize() BIND(C, name='lammps_mpi_finalize') - END SUBROUTINE lammps_mpi_finalize + SUBROUTINE lammps_mpi_finalize() BIND(C) + END SUBROUTINE lammps_mpi_finalize - SUBROUTINE lammps_kokkos_finalize() BIND(C, name='lammps_kokkos_finalize') - END SUBROUTINE lammps_kokkos_finalize + SUBROUTINE lammps_kokkos_finalize() BIND(C) + END SUBROUTINE lammps_kokkos_finalize - SUBROUTINE lammps_file(handle, filename) BIND(C, name='lammps_file') - IMPORT :: c_ptr - TYPE(c_ptr), VALUE :: handle - TYPE(c_ptr), VALUE :: filename - END SUBROUTINE lammps_file + SUBROUTINE lammps_error(handle, error_type, error_text) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + INTEGER(c_int), VALUE :: error_type + TYPE(c_ptr), VALUE :: error_text + END SUBROUTINE lammps_error - SUBROUTINE lammps_command(handle, cmd) BIND(C, name='lammps_command') - IMPORT :: c_ptr - TYPE(c_ptr), VALUE :: handle - TYPE(c_ptr), VALUE :: cmd - END SUBROUTINE lammps_command + SUBROUTINE lammps_file(handle, filename) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + TYPE(c_ptr), VALUE :: filename + END SUBROUTINE lammps_file - SUBROUTINE lammps_commands_list(handle, ncmd, cmds) BIND(C, name='lammps_commands_list') - IMPORT :: c_ptr, c_int - TYPE(c_ptr), VALUE :: handle - INTEGER(c_int), VALUE, INTENT(in) :: ncmd - TYPE(c_ptr), DIMENSION(*), INTENT(in) :: cmds - END SUBROUTINE lammps_commands_list + SUBROUTINE lammps_command(handle, cmd) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + TYPE(c_ptr), VALUE :: cmd + END SUBROUTINE lammps_command - SUBROUTINE lammps_commands_string(handle, str) BIND(C, name='lammps_commands_string') - IMPORT :: c_ptr - TYPE(c_ptr), VALUE :: handle - TYPE(c_ptr), VALUE :: str - END SUBROUTINE lammps_commands_string + SUBROUTINE lammps_commands_list(handle, ncmd, cmds) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + INTEGER(c_int), VALUE, INTENT(IN) :: ncmd + TYPE(c_ptr), DIMENSION(*), INTENT(IN) :: cmds + END SUBROUTINE lammps_commands_list - FUNCTION lammps_malloc(size) BIND(C, name='malloc') - IMPORT :: c_ptr, c_size_t - INTEGER(c_size_t), value :: size - TYPE(c_ptr) :: lammps_malloc - END FUNCTION lammps_malloc + SUBROUTINE lammps_commands_string(handle, str) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + TYPE(c_ptr), VALUE :: str + END SUBROUTINE lammps_commands_string - SUBROUTINE lammps_free(ptr) BIND(C, name='lammps_free') - IMPORT :: c_ptr - TYPE(c_ptr), VALUE :: ptr - END SUBROUTINE lammps_free + FUNCTION lammps_get_natoms(handle) BIND(C) + IMPORT :: c_ptr, c_double + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + REAL(c_double) :: lammps_get_natoms + END FUNCTION lammps_get_natoms - FUNCTION lammps_version(handle) BIND(C, name='lammps_version') - IMPORT :: c_ptr, c_int - TYPE(c_ptr), VALUE :: handle - INTEGER(c_int) :: lammps_version - END FUNCTION lammps_version + FUNCTION lammps_get_thermo(handle,name) BIND(C) + IMPORT :: c_ptr, c_double + IMPLICIT NONE + REAL(c_double) :: lammps_get_thermo + TYPE(c_ptr), VALUE :: handle + TYPE(c_ptr), VALUE :: name + END FUNCTION lammps_get_thermo + + SUBROUTINE lammps_extract_box(handle,boxlo,boxhi,xy,yz,xz,pflags, & + boxflag) BIND(C) + IMPORT :: c_ptr, c_double, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, boxlo, boxhi, xy, yz, xz, pflags, & + boxflag + END SUBROUTINE lammps_extract_box + + SUBROUTINE lammps_reset_box(handle,boxlo,boxhi,xy,yz,xz) BIND(C) + IMPORT :: c_ptr, c_double + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + REAL(c_double), DIMENSION(3) :: boxlo, boxhi + REAL(c_double), VALUE :: xy, yz, xz + END SUBROUTINE lammps_reset_box + + SUBROUTINE lammps_memory_usage(handle,meminfo) BIND(C) + IMPORT :: c_ptr, c_double + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + REAL(c_double), DIMENSION(*) :: meminfo + END SUBROUTINE lammps_memory_usage + + FUNCTION lammps_get_mpi_comm(handle) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + INTEGER(c_int) :: lammps_get_mpi_comm + END FUNCTION lammps_get_mpi_comm + + FUNCTION lammps_extract_setting(handle,keyword) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, keyword + INTEGER(c_int) :: lammps_extract_setting + END FUNCTION lammps_extract_setting + + FUNCTION lammps_extract_global_datatype(handle,name) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, name + INTEGER(c_int) :: lammps_extract_global_datatype + END FUNCTION lammps_extract_global_datatype + + FUNCTION c_strlen (str) BIND(C,name='strlen') + IMPORT :: c_ptr, c_size_t + IMPLICIT NONE + TYPE(c_ptr), VALUE :: str + INTEGER(c_size_t) :: c_strlen + END FUNCTION c_strlen + + FUNCTION lammps_extract_global(handle, name) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, name + TYPE(c_ptr) :: lammps_extract_global + END FUNCTION lammps_extract_global + + !INTEGER (c_int) FUNCTION lammps_extract_atom_datatype + + !(generic) lammps_extract_atom + + !(generic) lammps_extract_compute + + !(generic) lammps_extract_fix + + !(generic) lammps_extract_variable + + !INTEGER (c_int) lammps_set_variable + + !SUBROUTINE lammps_gather_atoms + + !SUBROUTINE lammps_gather_atoms_concat + + !SUBROUTINE lammps_gather_atoms_subset + + !SUBROUTINE lammps_scatter_atoms + + !SUBROUTINE lammps_scatter_atoms_subset + + !SUBROUTINE lammps_gather_bonds + + !SUBROUTINE lammps_gather + + !SUBROUTINE lammps_gather_concat + + !SUBROUTINE lammps_gather_subset + + !SUBROUTINE lammps_scatter_subset + + !(generic / id, type, and image are special) / requires LAMMPS_BIGBIG + !INTEGER (C_int) FUNCTION lammps_create_atoms + + !INTEGER (C_int) FUNCTION lammps_find_pair_neighlist + + !INTEGER (C_int) FUNCTION lammps_find_fix_neighlist + + !INTEGER (C_int) FUNCTION lammps_find_compute_neighlist + + !INTEGER (C_int) FUNCTION lammps_neighlist_num_elements + + !SUBROUTINE lammps_neighlist_element_neighbors + + FUNCTION lammps_version(handle) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + INTEGER(c_int) :: lammps_version + END FUNCTION lammps_version + + !SUBROUTINE lammps_get_os_info + + !LOGICAL FUNCTION lammps_config_has_mpi_support + !LOGICAL FUNCTION lammps_config_has_gzip_support + !LOGICAL FUNCTION lammps_config_has_png_support + !LOGICAL FUNCTION lammps_config_has_jpeg_support + !LOGICAL FUNCTION lammps_config_has_ffmpeg_support + !LOGICAL FUNCTION lammps_config_has_exceptions + !LOGICAL FUNCTION lammps_config_has_package + !INTEGER (C_int) FUNCTION lammps_config_package_count + !SUBROUTINE lammps_config_package_name + + !LOGICAL FUNCTION lammps_config_accelerator + !LOGICAL FUNCTION lammps_has_gpu_device + !SUBROUTINE lammps_get_gpu_device + + !LOGICAL FUNCTION lammps_has_id + !INTEGER (C_int) FUNCTION lammps_id_count + !SUBROUTINE lammps_id_name + + !INTEGER (C_int) FUNCTION lammps_plugin_count + !SUBROUTINE lammps_plugin_name + + !Both of these use LAMMPS_BIGBIG + !INTEGER (LAMMPS_imageint) FUNCTION lammps_encode_image_flags + !SUBROUTINE lammps_decode_image_flags + + !SUBROUTINE lammps_set_fix_external_callback ! may have trouble.... + !FUNCTION lammps_fix_external_get_force() ! returns real(c_double) (:) + + !SUBROUTINE lammps_fix_external_set_energy_global + !SUBROUTINE lammps_fix_external_set_energy_peratom + !SUBROUTINE lammps_fix_external_set_virial_global + !SUBROUTINE lammps_fix_external_set_virial_peratom + !SUBROUTINE lammps_fix_external_set_vector_length + !SUBROUTINE lammps_fix_external_set_vector + + !SUBROUTINE lammps_flush_buffers + + FUNCTION lammps_malloc(size) BIND(C, name='malloc') + IMPORT :: c_ptr, c_size_t + IMPLICIT NONE + INTEGER(c_size_t), VALUE :: size + TYPE(c_ptr) :: lammps_malloc + END FUNCTION lammps_malloc + + SUBROUTINE lammps_free(ptr) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: ptr + END SUBROUTINE lammps_free + + INTEGER(c_int) FUNCTION lammps_is_running(handle) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + END FUNCTION lammps_is_running + + !SUBROUTINE lammps_force_timeout + + !LOGICAL FUNCTION lammps_has_error + + !INTEGER (c_int) FUNCTION lammps_get_last_error_message - FUNCTION lammps_get_natoms(handle) BIND(C, name='lammps_get_natoms') - IMPORT :: c_ptr, c_double - TYPE(c_ptr), VALUE :: handle - REAL(c_double) :: lammps_get_natoms - END FUNCTION lammps_get_natoms END INTERFACE CONTAINS @@ -138,9 +379,8 @@ CONTAINS ! Constructor for the LAMMPS class. ! Combined wrapper around lammps_open_fortran() and lammps_open_no_mpi() TYPE(lammps) FUNCTION lmp_open(args, comm) - IMPLICIT NONE - INTEGER, INTENT(in), OPTIONAL :: comm - CHARACTER(len=*), INTENT(in), OPTIONAL :: args(:) + INTEGER, INTENT(IN), OPTIONAL :: comm + CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: args(:) TYPE(c_ptr), ALLOCATABLE :: argv(:) INTEGER(c_int) :: i, c_comm, argc @@ -173,9 +413,8 @@ CONTAINS ! Combined Fortran wrapper around lammps_close() and lammps_mpi_finalize() SUBROUTINE lmp_close(self, finalize) - IMPLICIT NONE CLASS(lammps) :: self - LOGICAL, INTENT(in), OPTIONAL :: finalize + LOGICAL, INTENT(IN), OPTIONAL :: finalize CALL lammps_close(self%handle) @@ -187,22 +426,20 @@ CONTAINS END IF END SUBROUTINE lmp_close - INTEGER FUNCTION lmp_version(self) - IMPLICIT NONE + ! equivalent function to lammps_error() + SUBROUTINE lmp_error(self, error_type, error_text) CLASS(lammps) :: self + INTEGER :: error_type + CHARACTER(len=*) :: error_text + TYPE(c_ptr) :: str - lmp_version = lammps_version(self%handle) - END FUNCTION lmp_version - - DOUBLE PRECISION FUNCTION lmp_get_natoms(self) - IMPLICIT NONE - CLASS(lammps) :: self - - lmp_get_natoms = lammps_get_natoms(self%handle) - END FUNCTION lmp_get_natoms + str = f2c_string(error_text) + CALL lammps_error(self%handle, error_type, str) + CALL lammps_free(str) + END SUBROUTINE lmp_error + ! equivalent function to lammps_file() SUBROUTINE lmp_file(self, filename) - IMPLICIT NONE CLASS(lammps) :: self CHARACTER(len=*) :: filename TYPE(c_ptr) :: str @@ -214,7 +451,6 @@ CONTAINS ! equivalent function to lammps_command() SUBROUTINE lmp_command(self, cmd) - IMPLICIT NONE CLASS(lammps) :: self CHARACTER(len=*) :: cmd TYPE(c_ptr) :: str @@ -226,9 +462,8 @@ CONTAINS ! equivalent function to lammps_commands_list() SUBROUTINE lmp_commands_list(self, cmds) - IMPLICIT NONE CLASS(lammps) :: self - CHARACTER(len=*), INTENT(in), OPTIONAL :: cmds(:) + CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: cmds(:) TYPE(c_ptr), ALLOCATABLE :: cmdv(:) INTEGER :: i, ncmd @@ -250,7 +485,6 @@ CONTAINS ! equivalent function to lammps_commands_string() SUBROUTINE lmp_commands_string(self, str) - IMPLICIT NONE CLASS(lammps) :: self CHARACTER(len=*) :: str TYPE(c_ptr) :: tmp @@ -260,13 +494,274 @@ CONTAINS CALL lammps_free(tmp) END SUBROUTINE lmp_commands_string + ! equivalent function to lammps_get_natoms + DOUBLE PRECISION FUNCTION lmp_get_natoms(self) + CLASS(lammps) :: self + + lmp_get_natoms = lammps_get_natoms(self%handle) + END FUNCTION lmp_get_natoms + + ! equivalent function to lammps_get_thermo + REAL (C_double) FUNCTION lmp_get_thermo(self,name) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*) :: name + TYPE(C_ptr) :: Cname + + Cname = f2c_string(name) + lmp_get_thermo = lammps_get_thermo(self%handle, Cname) + CALL lammps_free(Cname) + END FUNCTION lmp_get_thermo + + ! equivalent subroutine to lammps_extract_box + SUBROUTINE lmp_extract_box(self, boxlo, boxhi, xy, yz, xz, pflags, boxflag) + CLASS(lammps), INTENT(IN) :: self + REAL(c_double), INTENT(OUT), TARGET, OPTIONAL :: boxlo(3), boxhi(3) + REAL(c_double), INTENT(OUT), TARGET, OPTIONAL :: xy, yz, xz + LOGICAL, INTENT(OUT), OPTIONAL :: pflags(3), boxflag + INTEGER(c_int), TARGET :: C_pflags(3), C_boxflag + TYPE (c_ptr) :: ptr(7) + + ptr = c_null_ptr + IF ( PRESENT(boxlo) ) ptr(1) = C_LOC(boxlo(1)) + IF ( PRESENT(boxhi) ) ptr(2) = C_LOC(boxhi(1)) + IF ( PRESENT(xy) ) ptr(3) = C_LOC(xy) + IF ( PRESENT(yz) ) ptr(4) = C_LOC(yz) + IF ( PRESENT(xz) ) ptr(5) = C_LOC(xz) + IF ( PRESENT(pflags) ) ptr(6) = C_LOC(C_pflags(1)) + IF ( PRESENT(boxflag) ) ptr(7) = C_LOC(C_boxflag) + CALL lammps_extract_box(self%handle, ptr(1), ptr(2), ptr(3), ptr(4), & + ptr(5), ptr(6), ptr(7)) + IF ( PRESENT(pflags) ) pflags = ( C_pflags /= 0_C_int ) + IF ( PRESENT(boxflag) ) boxflag = ( C_boxflag /= 0_C_int ) + END SUBROUTINE lmp_extract_box + + ! equivalent function to lammps_reset_box + SUBROUTINE lmp_reset_box(self, boxlo, boxhi, xy, yz, xz) + CLASS(lammps), INTENT(IN) :: self + REAL(C_double), INTENT(IN) :: boxlo(3), boxhi(3), xy, yz, xz + + CALL lammps_reset_box(self%handle, boxlo, boxhi, xy, yz, xz) + END SUBROUTINE lmp_reset_box + + ! equivalent function to lammps_memory_usage + SUBROUTINE lmp_memory_usage(self,meminfo) + CLASS(lammps), INTENT(IN) :: self + INTEGER, PARAMETER :: MEMINFO_ELEM = 3 + REAL (c_double), DIMENSION(MEMINFO_ELEM), INTENT(OUT) :: meminfo + + CALL lammps_memory_usage(self%handle,meminfo) + END SUBROUTINE lmp_memory_usage + + ! equivalent function to lammps_get_mpi_comm + INTEGER FUNCTION lmp_get_mpi_comm(self) + CLASS(lammps), INTENT(IN) :: self + + lmp_get_mpi_comm = lammps_get_mpi_comm(self%handle) + END FUNCTION lmp_get_mpi_comm + + ! equivalent function to lammps_extract_setting + INTEGER (c_int) FUNCTION lmp_extract_setting(self, keyword) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: keyword + TYPE(c_ptr) :: Ckeyword + + Ckeyword = f2c_string(keyword) + lmp_extract_setting = lammps_extract_setting(self%handle, Ckeyword) + CALL lammps_free(Ckeyword) + END FUNCTION lmp_extract_setting + + ! equivalent function to lammps_extract_global + ! the assignment is actually overloaded so as to bind the pointers to + ! lammps data based on the information available from LAMMPS + FUNCTION lmp_extract_global(self, name) RESULT (global_data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + TYPE(lammps_data) :: global_data + + INTEGER(c_int) :: datatype + TYPE(c_ptr) :: Cname, Cptr + INTEGER(c_size_t) :: length, i + CHARACTER(KIND=c_char, LEN=1), DIMENSION(:), POINTER :: Fptr + + ! Determine vector length + ! FIXME Is there a way to get the length of the vector from C rather + ! than defining it here AND in the Python API? + SELECT CASE (name) + CASE ('boxlo','boxhi','sublo','subhi','sublo_lambda','subhi_lambda', & + 'periodicity') + length = 3 + CASE DEFAULT + length = 1 + ! string cases are overridden later + END SELECT + + Cname = f2c_string(name) + datatype = lammps_extract_global_datatype(self%handle, Cname) + ! above could be c_null_ptr in place of self%handle...doesn't matter + Cptr = lammps_extract_global(self%handle, Cname) + CALL lammps_free(Cname) + + SELECT CASE (datatype) + CASE (LAMMPS_INT) + IF ( length == 1 ) THEN + global_data%datatype = DATA_INT + CALL C_F_POINTER(Cptr, global_data%i32) + ELSE + global_data%datatype = DATA_INT_1D + CALL C_F_POINTER(Cptr, global_data%i32_vec, [length]) + END IF + CASE (LAMMPS_INT64) + IF ( length == 1 ) THEN + global_data%datatype = DATA_INT64 + CALL C_F_POINTER(Cptr, global_data%i64) + ELSE + global_data%datatype = DATA_INT64_1D + CALL C_F_POINTER(Cptr, global_data%i64_vec, [length]) + END IF + CASE (LAMMPS_DOUBLE) + IF ( length == 1 ) THEN + global_data%datatype = DATA_DOUBLE + CALL C_F_POINTER(Cptr, global_data%r64) + ELSE + global_data%datatype = DATA_DOUBLE_1D + CALL C_F_POINTER(Cptr, global_data%r64_vec, [length]) + END IF + CASE (LAMMPS_STRING) + global_data%datatype = DATA_STRING + length = c_strlen(Cptr) + CALL C_F_POINTER(Cptr, Fptr, [length]) + ALLOCATE ( CHARACTER(LEN=length) :: global_data%str ) + FORALL ( I=1:length ) + global_data%str(i:i) = Fptr(i) + END FORALL + CASE DEFAULT + ! FIXME convert to use symbolic constants later + CALL lmp_error(self, 6, 'Unknown pointer type in extract_global') + END SELECT + END FUNCTION + + ! equivalent function to lammps_version() + INTEGER FUNCTION lmp_version(self) + CLASS(lammps) :: self + + lmp_version = lammps_version(self%handle) + END FUNCTION lmp_version + + ! equivalent function to lammps_is_running + LOGICAL FUNCTION lmp_is_running(self) + CLASS(lammps) :: self + + lmp_is_running = ( lammps_is_running(self%handle) /= 0_C_int ) + END FUNCTION lmp_is_running + + ! ---------------------------------------------------------------------- + ! functions to assign user-space pointers to LAMMPS data + ! ---------------------------------------------------------------------- + SUBROUTINE assign_int_to_lammps_data (lhs, rhs) + INTEGER(c_int), INTENT(OUT), POINTER :: lhs + CLASS(lammps_data), INTENT(IN) :: rhs + + IF ( rhs%datatype == DATA_INT ) THEN + lhs => rhs%i32 + ELSE + CALL assignment_error(rhs%datatype, 'scalar int') + END IF + END SUBROUTINE assign_int_to_lammps_data + + SUBROUTINE assign_int64_to_lammps_data (lhs, rhs) + INTEGER(c_int64_t), INTENT(OUT), POINTER :: lhs + CLASS(lammps_data), INTENT(IN) :: rhs + + IF ( rhs%datatype == DATA_INT64 ) THEN + lhs => rhs%i64 + ELSE + CALL assignment_error(rhs%datatype, 'scalar long int') + END IF + END SUBROUTINE assign_int64_to_lammps_data + + SUBROUTINE assign_intvec_to_lammps_data (lhs, rhs) + INTEGER(c_int), DIMENSION(:), INTENT(OUT), POINTER :: lhs + CLASS(lammps_data), INTENT(IN) :: rhs + + IF ( rhs%datatype == DATA_INT_1D ) THEN + lhs => rhs%i32_vec + ELSE + CALL assignment_error(rhs%datatype, 'vector of ints') + END IF + END SUBROUTINE assign_intvec_to_lammps_data + + SUBROUTINE assign_double_to_lammps_data (lhs, rhs) + REAL(c_double), INTENT(OUT), POINTER :: lhs + CLASS(lammps_data), INTENT(IN) :: rhs + + IF ( rhs%datatype == DATA_DOUBLE ) THEN + lhs => rhs%r64 + ELSE + CALL assignment_error(rhs%datatype, 'scalar double') + END IF + END SUBROUTINE assign_double_to_lammps_data + + SUBROUTINE assign_doublevec_to_lammps_data (lhs, rhs) + REAL(c_double), DIMENSION(:), INTENT(OUT), POINTER :: lhs + CLASS(lammps_data), INTENT(IN) :: rhs + + IF ( rhs%datatype == DATA_DOUBLE_1D ) THEN + lhs => rhs%r64_vec + ELSE + CALL assignment_error(rhs%datatype, 'vector of doubles') + END IF + END SUBROUTINE assign_doublevec_to_lammps_data + + SUBROUTINE assign_string_to_lammps_data (lhs, rhs) + CHARACTER(LEN=*), INTENT(OUT) :: lhs + CLASS(lammps_data), INTENT(IN) :: rhs + + IF ( rhs%datatype == DATA_STRING ) THEN + lhs = rhs%str + ELSE + CALL assignment_error(rhs%datatype, 'string') + END IF + END SUBROUTINE assign_string_to_lammps_data + + SUBROUTINE assignment_error (type1, type2) + INTEGER (c_int) :: type1 + CHARACTER (LEN=*) :: type2 + INTEGER, PARAMETER :: ERROR_CODE = 1 + CHARACTER (LEN=:), ALLOCATABLE :: str1 + + SELECT CASE (type1) + CASE (DATA_INT) + str1 = 'scalar int' + CASE (DATA_INT_1D) + str1 = 'vector of ints' + CASE (DATA_INT_2D) + str1 = 'matrix of ints' + CASE (DATA_INT64) + str1 = 'scalar long int' + CASE (DATA_INT64_1D) + str1 = 'vector of long ints' + CASE (DATA_INT64_2D) + str1 = 'matrix of long ints' + CASE (DATA_DOUBLE) + str1 = 'scalar double' + CASE (DATA_DOUBLE_1D) + str1 = 'vector of doubles' + CASE (DATA_DOUBLE_2D) + str1 = 'matrix of doubles' + CASE DEFAULT + str1 = 'that type' + END SELECT + WRITE (ERROR_UNIT,'(A)') 'Cannot associate ' // str1 // ' with ' // type2 + STOP ERROR_CODE + END SUBROUTINE assignment_error + ! ---------------------------------------------------------------------- ! local helper functions ! copy fortran string to zero terminated c string ! ---------------------------------------------------------------------- FUNCTION f2c_string(f_string) RESULT(ptr) - CHARACTER (len=*), INTENT(in) :: f_string - CHARACTER (len=1, kind=c_char), POINTER :: c_string(:) + CHARACTER(LEN=*), INTENT(IN) :: f_string + CHARACTER(LEN=1, KIND=c_char), POINTER :: c_string(:) TYPE(c_ptr) :: ptr INTEGER(c_size_t) :: i, n @@ -279,3 +774,5 @@ CONTAINS c_string(n+1) = c_null_char END FUNCTION f2c_string END MODULE LIBLAMMPS + +! vim: ts=2 sts=2 sw=2 et diff --git a/python/lammps/constants.py b/python/lammps/constants.py index a50d58b28f..26bb92626a 100644 --- a/python/lammps/constants.py +++ b/python/lammps/constants.py @@ -13,29 +13,35 @@ # various symbolic constants to be used # in certain calls to select data formats -LAMMPS_AUTODETECT = None -LAMMPS_INT = 0 -LAMMPS_INT_2D = 1 -LAMMPS_DOUBLE = 2 -LAMMPS_DOUBLE_2D = 3 -LAMMPS_INT64 = 4 -LAMMPS_INT64_2D = 5 -LAMMPS_STRING = 6 +LAMMPS_AUTODETECT = None +LAMMPS_INT = 0 +LAMMPS_INT_2D = 1 +LAMMPS_DOUBLE = 2 +LAMMPS_DOUBLE_2D = 3 +LAMMPS_INT64 = 4 +LAMMPS_INT64_2D = 5 +LAMMPS_STRING = 6 # these must be kept in sync with the enums in library.h -LMP_STYLE_GLOBAL = 0 -LMP_STYLE_ATOM = 1 -LMP_STYLE_LOCAL = 2 +LMP_STYLE_GLOBAL = 0 +LMP_STYLE_ATOM = 1 +LMP_STYLE_LOCAL = 2 -LMP_TYPE_SCALAR = 0 -LMP_TYPE_VECTOR = 1 -LMP_TYPE_ARRAY = 2 -LMP_SIZE_VECTOR = 3 -LMP_SIZE_ROWS = 4 -LMP_SIZE_COLS = 5 +LMP_TYPE_SCALAR = 0 +LMP_TYPE_VECTOR = 1 +LMP_TYPE_ARRAY = 2 +LMP_SIZE_VECTOR = 3 +LMP_SIZE_ROWS = 4 +LMP_SIZE_COLS = 5 -LMP_VAR_EQUAL = 0 -LMP_VAR_ATOM = 1 +LMP_ERROR_WARNING = 0 +LMP_ERROR_ONE = 1 +LMP_ERROR_ALL = 2 +LMP_ERROR_WORLD = 4 +LMP_ERROR_UNIVERSE = 8 + +LMP_VAR_EQUAL = 0 +LMP_VAR_ATOM = 1 # ------------------------------------------------------------------------- diff --git a/python/lammps/core.py b/python/lammps/core.py index 930a40a4b0..aa4aae13db 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -167,6 +167,8 @@ class lammps(object): self.lib.lammps_flush_buffers.argtypes = [c_void_p] self.lib.lammps_free.argtypes = [c_void_p] + self.lib.lammps_error.argtypes = [c_void_p, c_int, c_char_p] + self.lib.lammps_file.argtypes = [c_void_p, c_char_p] self.lib.lammps_file.restype = None @@ -486,6 +488,26 @@ class lammps(object): # ------------------------------------------------------------------------- + def error(self, error_type, error_text): + """Forward error to the LAMMPS Error class. + + This is a wrapper around the :cpp:func:`lammps_error` function of the C-library interface. + + .. versionadded:: TBD + + :param error_type: + :type error_type: int + :param error_text: + :type error_text: string + """ + if error_text: error_text = error_text.encode() + else: error_text = "(unknown error)".encode() + + with ExceptionCheck(self): + self.lib.lammps_error(self.lmp, error_type, error_text) + + # ------------------------------------------------------------------------- + def version(self): """Return a numerical representation of the LAMMPS version in use. @@ -1622,8 +1644,8 @@ class lammps(object): """Return a string with detailed information about any devices that are usable by the GPU package. - This is a wrapper around the :cpp:func:`lammps_get_gpu_device_info` - function of the C-library interface. + This is a wrapper around the :cpp:func:`lammps_get_gpu_device_info` + function of the C-library interface. :return: GPU device info string :rtype: string diff --git a/python/lammps/mliap/__init__.py b/python/lammps/mliap/__init__.py index cfe3fb6b38..b888eb55ce 100644 --- a/python/lammps/mliap/__init__.py +++ b/python/lammps/mliap/__init__.py @@ -5,6 +5,7 @@ import sysconfig import ctypes import platform +import warnings py_ver = sysconfig.get_config_vars('VERSION')[0] OS_name = platform.system() @@ -25,8 +26,10 @@ except Exception as e: raise OSError("Unable to locate python shared library") from e if not pylib.Py_IsInitialized(): - raise RuntimeError("This interpreter is not compatible with python-based mliap for LAMMPS.") + warnings.warn("This interpreter is not compatible with python-based MLIAP for LAMMPS. " + "Attempting to activate the MLIAP-python coupling from python may result " + "in undefined behavior.") +else: + from .loader import load_model, load_unified, activate_mliappy del sysconfig, ctypes, library, pylib - -from .loader import load_model, activate_mliappy diff --git a/python/lammps/mliap/loader.py b/python/lammps/mliap/loader.py index dff791bfc1..e1e60c4f50 100644 --- a/python/lammps/mliap/loader.py +++ b/python/lammps/mliap/loader.py @@ -17,27 +17,58 @@ import sys +import importlib import importlib.util import importlib.machinery +import importlib.abc + +from ctypes import pythonapi, c_int, c_void_p, py_object + +# This dynamic loader imports a python module embedded in a shared library. +# The default value of api_version is 1013 because it has been stable since 2006. +class DynamicLoader(importlib.abc.Loader): + def __init__(self,module_name,library,api_version=1013): + self.api_version = api_version + + attr = "PyInit_"+module_name + initfunc = getattr(library,attr) + # c_void_p is standin for PyModuleDef * + initfunc.restype = c_void_p + initfunc.argtypes = () + self.module_def = initfunc() + + def create_module(self, spec): + createfunc = pythonapi.PyModule_FromDefAndSpec2 + # c_void_p is standin for PyModuleDef * + createfunc.argtypes = c_void_p, py_object, c_int + createfunc.restype = py_object + module = createfunc(self.module_def, spec, self.api_version) + return module + + def exec_module(self, module): + execfunc = pythonapi.PyModule_ExecDef + # c_void_p is standin for PyModuleDef * + execfunc.argtypes = py_object, c_void_p + execfunc.restype = c_int + result = execfunc(module, self.module_def) + if result<0: + raise ImportError() def activate_mliappy(lmp): try: - # Begin Importlib magic to find the embedded python module - # This is needed because the filename for liblammps does not - # match the spec for normal python modules, wherein - # file names match with PyInit function names. - # Also, python normally doesn't look for extensions besides '.so' - # We fix both of these problems by providing an explict - # path to the extension module 'mliap_model_python_couple' in + library = lmp.lib + module_names = ["mliap_model_python_couple", "mliap_unified_couple"] + api_version = library.lammps_python_api_version() - path = lmp.lib._name - loader = importlib.machinery.ExtensionFileLoader('mliap_model_python_couple', path) - spec = importlib.util.spec_from_loader('mliap_model_python_couple', loader) - module = importlib.util.module_from_spec(spec) - sys.modules['mliap_model_python_couple'] = module - spec.loader.exec_module(module) - # End Importlib magic to find the embedded python module + for module_name in module_names: + # Make Machinery + loader = DynamicLoader(module_name,library,api_version) + spec = importlib.util.spec_from_loader(module_name,loader) + # Do the import + module = importlib.util.module_from_spec(spec) + sys.modules[module_name] = module + spec.loader.exec_module(module) except Exception as ee: raise ImportError("Could not load ML-IAP python coupling module.") from ee @@ -50,3 +81,11 @@ def load_model(model): ) from ie mliap_model_python_couple.load_from_python(model) +def load_unified(model): + try: + import mliap_unified_couple + except ImportError as ie: + raise ImportError("ML-IAP python module must be activated before loading\n" + "the pair style. Call lammps.mliap.activate_mliappy(lmp)." + ) from ie + mliap_unified_couple.load_from_python(model) diff --git a/python/lammps/mliap/mliap_unified_abc.py b/python/lammps/mliap/mliap_unified_abc.py new file mode 100644 index 0000000000..3243df85c4 --- /dev/null +++ b/python/lammps/mliap/mliap_unified_abc.py @@ -0,0 +1,28 @@ +from abc import ABC, abstractmethod +import pickle + +class MLIAPUnified(ABC): + """Abstract base class for MLIAPUnified.""" + + def __init__(self): + self.interface = None + self.element_types = None + self.ndescriptors = None + self.nparams = None + self.rcutfac = None + + @abstractmethod + def compute_gradients(self, data): + """Compute gradients.""" + + @abstractmethod + def compute_descriptors(self, data): + """Compute descriptors.""" + + @abstractmethod + def compute_forces(self, data): + """Compute forces.""" + + def pickle(self, fname): + with open(fname, 'wb') as fp: + pickle.dump(self, fp) diff --git a/python/lammps/mliap/mliap_unified_lj.py b/python/lammps/mliap/mliap_unified_lj.py new file mode 100644 index 0000000000..37e7170d4a --- /dev/null +++ b/python/lammps/mliap/mliap_unified_lj.py @@ -0,0 +1,44 @@ +from .mliap_unified_abc import MLIAPUnified +import numpy as np + + +class MLIAPUnifiedLJ(MLIAPUnified): + """Test implementation for MLIAPUnified.""" + + def __init__(self, element_types, epsilon=1.0, sigma=1.0, rcutfac=1.25): + super().__init__() + self.element_types = element_types + self.ndescriptors = 1 + self.nparams = 3 + # Mimicking the LJ pair-style: + # pair_style lj/cut 2.5 + # pair_coeff * * 1 1 + self.epsilon = epsilon + self.sigma = sigma + self.rcutfac = rcutfac + + def compute_gradients(self, data): + """Test compute_gradients.""" + + def compute_descriptors(self, data): + """Test compute_descriptors.""" + + def compute_forces(self, data): + """Test compute_forces.""" + eij, fij = self.compute_pair_ef(data) + data.update_pair_energy(eij) + data.update_pair_forces(fij) + + def compute_pair_ef(self, data): + rij = data.rij + + r2inv = 1.0 / np.sum(rij ** 2, axis=1) + r6inv = r2inv * r2inv * r2inv + + lj1 = 4.0 * self.epsilon * self.sigma**12 + lj2 = 4.0 * self.epsilon * self.sigma**6 + + eij = r6inv * (lj1 * r6inv - lj2) + fij = r6inv * (3.0 * lj2 - 6.0 * lj2 * r6inv) * r2inv + fij = fij[:, np.newaxis] * rij + return eij, fij diff --git a/python/lammps/mliap/pytorch.py b/python/lammps/mliap/pytorch.py index d699c239b0..51b953fd61 100644 --- a/python/lammps/mliap/pytorch.py +++ b/python/lammps/mliap/pytorch.py @@ -80,10 +80,10 @@ class TorchWrapper(torch.nn.Module): n_params : torch.nn.Module (None) Number of NN model parameters - + device : torch.nn.Module (None) Accelerator device - + dtype : torch.dtype (torch.float64) Dtype to use on device """ @@ -325,6 +325,6 @@ class ElemwiseModels(torch.nn.Module): per_atom_attributes = torch.zeros(elems.size(dim=0), dtype=self.dtype) given_elems, elem_indices = torch.unique(elems, return_inverse=True) for i, elem in enumerate(given_elems): - self.subnets[elem].to(self.dtype) + self.subnets[elem].to(self.dtype) per_atom_attributes[elem_indices == i] = self.subnets[elem](descriptors[elem_indices == i]).flatten() return per_atom_attributes diff --git a/src/ML-IAP/Install.sh b/src/ML-IAP/Install.sh index 3a1a7493aa..27c9b19b01 100755 --- a/src/ML-IAP/Install.sh +++ b/src/ML-IAP/Install.sh @@ -42,6 +42,7 @@ done # Install cython pyx file only if also Python is available action mliap_model_python_couple.pyx python_impl.cpp +action mliap_unified_couple.pyx python_impl.cpp # edit 2 Makefile.package files to include/exclude package info @@ -57,7 +58,7 @@ if (test $1 = 1) then include ..\/..\/lib\/python\/Makefile.mliap_python ' ../Makefile.package.settings fi - cythonize -3 ../mliap_model_python_couple.pyx + cythonize -3 ../mliap_model_python_couple.pyx ../mliap_unified_couple.pyx fi elif (test $1 = 0) then @@ -65,7 +66,8 @@ elif (test $1 = 0) then if (test -e ../Makefile.package) then sed -i -e 's/[^ \t]*-DMLIAP_PYTHON[^ \t]* //g' ../Makefile.package fi - rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h + rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h \ + ../mliap_unified_couple.cpp ../mliap_unified_couple.h sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings elif (test $1 = 2) then @@ -73,7 +75,8 @@ elif (test $1 = 2) then if (test -e ../Makefile.package) then sed -i -e 's/[^ \t]*-DMLIAP_PYTHON[^ \t]* //g' ../Makefile.package fi - rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h + rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h \ + ../mliap_unified_couple.cpp ../mliap_unified_couple.h sed -i -e '/^include.*python.*mliap_python.*$/d' ../Makefile.package.settings if (test -e ../Makefile.package) then sed -i -e 's|^PKG_INC =[ \t]*|&-DMLIAP_PYTHON |' ../Makefile.package @@ -85,8 +88,9 @@ elif (test $1 = 2) then include ..\/..\/lib\/python\/Makefile.mliap_python ' ../Makefile.package.settings fi - cythonize -3 ../mliap_model_python_couple.pyx + cythonize -3 ../mliap_model_python_couple.pyx ../mliap_unified_couple.pyx else - rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h + rm -f ../mliap_model_python_couple.cpp ../mliap_model_python_couple.h \ + ../mliap_unified_couple.cpp ../mliap_unified_couple.h fi fi diff --git a/src/ML-IAP/mliap_data.cpp b/src/ML-IAP/mliap_data.cpp index e96e87ba11..a91e66b211 100644 --- a/src/ML-IAP/mliap_data.cpp +++ b/src/ML-IAP/mliap_data.cpp @@ -31,11 +31,12 @@ MLIAPData::MLIAPData(LAMMPS *lmp, int gradgradflag_in, int *map_in, class MLIAPModel* model_in, class MLIAPDescriptor* descriptor_in, class PairMLIAP* pairmliap_in) : - Pointers(lmp), gradforce(nullptr), betas(nullptr), + Pointers(lmp), f(nullptr), gradforce(nullptr), betas(nullptr), descriptors(nullptr), eatoms(nullptr), gamma(nullptr), gamma_row_index(nullptr), gamma_col_index(nullptr), egradient(nullptr), - numneighs(nullptr), iatoms(nullptr), ielems(nullptr), jatoms(nullptr), jelems(nullptr), - rij(nullptr), graddesc(nullptr), model(nullptr), descriptor(nullptr), list(nullptr) + numneighs(nullptr), iatoms(nullptr), pair_i(nullptr), ielems(nullptr), + jatoms(nullptr), jelems(nullptr), elems(nullptr), rij(nullptr), + graddesc(nullptr), model(nullptr), descriptor(nullptr), list(nullptr) { gradgradflag = gradgradflag_in; map = map_in; @@ -85,10 +86,12 @@ MLIAPData::~MLIAPData() memory->destroy(gradforce); memory->destroy(iatoms); + memory->destroy(pair_i); memory->destroy(ielems); memory->destroy(numneighs); memory->destroy(jatoms); memory->destroy(jelems); + memory->destroy(elems); memory->destroy(rij); memory->destroy(graddesc); } @@ -107,6 +110,7 @@ void MLIAPData::init() void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_in) { list = list_in; + f = atom->f; double **x = atom->x; int *type = atom->type; @@ -114,21 +118,25 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - // grow nmax gradforce array if necessary + // grow nmax gradforce, elems arrays if necessary if (atom->nmax > nmax) { nmax = atom->nmax; memory->grow(gradforce,nmax,size_gradforce, "MLIAPData:gradforce"); + memory->grow(elems,nmax,"MLIAPData:elems"); } - // clear gradforce array + // clear gradforce, elems arrays int nall = atom->nlocal + atom->nghost; - for (int i = 0; i < nall; i++) + ntotal = nall; + for (int i = 0; i < nall; i++) { for (int j = 0; j < size_gradforce; j++) { gradforce[i][j] = 0.0; } + elems[i] = 0; + } // grow arrays if necessary @@ -153,8 +161,9 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i grow_neigharrays(); + npairs = 0; int ij = 0; - for (int ii = 0; ii < nlistatoms; ii++) { + for (int ii = 0; ii < natomneigh; ii++) { const int i = ilist[ii]; const double xtmp = x[i][0]; @@ -178,6 +187,7 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i const int jelem = map[jtype]; if (rsq < descriptor->cutsq[ielem][jelem]) { + pair_i[ij] = i; jatoms[ij] = j; jelems[ij] = jelem; rij[ij][0] = delx; @@ -190,6 +200,12 @@ void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_i iatoms[ii] = i; ielems[ii] = ielem; numneighs[ii] = ninside; + npairs += ninside; + } + + for (int i = 0; i < nall; i++) { + const int itype = type[i]; + elems[i] = map[itype]; } eflag = eflag_in; @@ -204,12 +220,14 @@ void MLIAPData::grow_neigharrays() { // grow neighbor atom arrays if necessary - - if (natomneigh_max < nlistatoms) { - memory->grow(iatoms,nlistatoms,"MLIAPData:iatoms"); - memory->grow(ielems,nlistatoms,"MLIAPData:ielems"); - memory->grow(numneighs,nlistatoms,"MLIAPData:numneighs"); - natomneigh_max = nlistatoms; + natomneigh = list->inum; + if (list->ghost == 1) + natomneigh += list->gnum; + if (natomneigh_max < natomneigh) { + memory->grow(iatoms,natomneigh,"MLIAPData:iatoms"); + memory->grow(ielems,natomneigh,"MLIAPData:ielems"); + memory->grow(numneighs,natomneigh,"MLIAPData:numneighs"); + natomneigh_max = natomneigh; } // grow neighbor arrays if necessary @@ -221,7 +239,7 @@ void MLIAPData::grow_neigharrays() int *type = atom->type; int nneigh = 0; - for (int ii = 0; ii < nlistatoms; ii++) { + for (int ii = 0; ii < natomneigh; ii++) { const int i = ilist[ii]; const double xtmp = x[i][0]; @@ -249,6 +267,7 @@ void MLIAPData::grow_neigharrays() } if (nneigh_max < nneigh) { + memory->grow(pair_i,nneigh,"MLIAPData:pair_i"); memory->grow(jatoms,nneigh,"MLIAPData:jatoms"); memory->grow(jelems,nneigh,"MLIAPData:jelems"); memory->grow(rij,nneigh,3,"MLIAPData:rij"); @@ -264,6 +283,7 @@ double MLIAPData::memory_usage() bytes += (double)nelements*nparams*sizeof(double); // egradient bytes += (double)nmax*size_gradforce*sizeof(double); // gradforce + bytes += (double)nmax*sizeof(int); // elems if (gradgradflag == 1) { bytes += (double)natomgamma_max* @@ -282,6 +302,7 @@ double MLIAPData::memory_usage() bytes += (double)natomneigh_max*sizeof(int); // ielems bytes += (double)natomneigh_max*sizeof(int); // numneighs + bytes += (double)nneigh_max*sizeof(int); // pair_i bytes += (double)nneigh_max*sizeof(int); // jatoms bytes += (double)nneigh_max*sizeof(int); // jelems bytes += (double)nneigh_max*3*sizeof(double); // rij" diff --git a/src/ML-IAP/mliap_data.h b/src/ML-IAP/mliap_data.h index 258261fa1a..562dcd63cb 100644 --- a/src/ML-IAP/mliap_data.h +++ b/src/ML-IAP/mliap_data.h @@ -35,6 +35,7 @@ class MLIAPData : protected Pointers { int size_gradforce; int yoffset, zoffset; int ndims_force, ndims_virial; + double **f; double **gradforce; double **betas; // betas for all atoms in list double **descriptors; // descriptors for all atoms in list @@ -56,15 +57,20 @@ class MLIAPData : protected Pointers { // data structures for mliap neighbor list // only neighbors strictly inside descriptor cutoff - int nlistatoms; // current number of atoms in neighborlist + int ntotal; // total number of owned and ghost atoms on this proc + int nlistatoms; // current number of atoms in local atom lists int nlistatoms_max; // allocated size of descriptor array + int natomneigh; // current number of atoms and ghosts in atom neighbor arrays int natomneigh_max; // allocated size of atom neighbor arrays int *numneighs; // neighbors count for each atom int *iatoms; // index of each atom int *ielems; // element of each atom int nneigh_max; // number of ij neighbors allocated + int npairs; // number of ij neighbor pairs + int *pair_i; // index of each i atom for each ij pair int *jatoms; // index of each neighbor int *jelems; // element of each neighbor + int *elems; // element of each atom in or not in the neighborlist double **rij; // distance vector of each neighbor double ***graddesc; // descriptor gradient w.r.t. each neighbor int eflag; // indicates if energy is needed diff --git a/src/ML-IAP/mliap_descriptor.cpp b/src/ML-IAP/mliap_descriptor.cpp index ff8cfb43a9..ce8abf4e56 100644 --- a/src/ML-IAP/mliap_descriptor.cpp +++ b/src/ML-IAP/mliap_descriptor.cpp @@ -25,7 +25,7 @@ using namespace LAMMPS_NS; MLIAPDescriptor::MLIAPDescriptor(LAMMPS *lmp) : Pointers(lmp), ndescriptors(0), nelements(0), elements(nullptr), cutsq(nullptr), - radelem(nullptr), wjelem(nullptr) + cutghost(nullptr), radelem(nullptr), wjelem(nullptr) { cutmax = 0.0; } @@ -37,6 +37,7 @@ MLIAPDescriptor::~MLIAPDescriptor() for (int i = 0; i < nelements; i++) delete[] elements[i]; delete[] elements; memory->destroy(cutsq); + memory->destroy(cutghost); memory->destroy(radelem); memory->destroy(wjelem); } diff --git a/src/ML-IAP/mliap_descriptor.h b/src/ML-IAP/mliap_descriptor.h index 156f8d1e92..6db872ccaf 100644 --- a/src/ML-IAP/mliap_descriptor.h +++ b/src/ML-IAP/mliap_descriptor.h @@ -29,13 +29,14 @@ class MLIAPDescriptor : protected Pointers { virtual void init() = 0; virtual double memory_usage(); - int ndescriptors; // number of descriptors - int nelements; // # of unique elements - char **elements; // names of unique elements - double **cutsq; // nelem x nelem rcutsq values - double cutmax; // maximum cutoff needed - double *radelem; // element radii - double *wjelem; // elements weights + int ndescriptors; // number of descriptors + int nelements; // # of unique elements + char **elements; // names of unique elements + double **cutsq; // nelem x nelem rcutsq values + double **cutghost; // cutoff for each ghost pair + double cutmax; // maximum cutoff needed + double *radelem; // element radii + double *wjelem; // elements weights protected: }; diff --git a/src/ML-IAP/mliap_model_python.h b/src/ML-IAP/mliap_model_python.h index 4243c67332..915e79e53c 100644 --- a/src/ML-IAP/mliap_model_python.h +++ b/src/ML-IAP/mliap_model_python.h @@ -21,20 +21,20 @@ namespace LAMMPS_NS { class MLIAPModelPython : public MLIAPModel { public: MLIAPModelPython(LAMMPS *, char * = nullptr); - ~MLIAPModelPython(); - virtual int get_nparams(); - virtual int get_gamma_nnz(class MLIAPData *); - virtual void compute_gradients(class MLIAPData *); - virtual void compute_gradgrads(class MLIAPData *); - virtual void compute_force_gradients(class MLIAPData *); - virtual double memory_usage(); + ~MLIAPModelPython() override; + int get_nparams() override; + int get_gamma_nnz(class MLIAPData *) override; + void compute_gradients(class MLIAPData *) override; + void compute_gradgrads(class MLIAPData *) override; + void compute_force_gradients(class MLIAPData *) override; + double memory_usage() override; void connect_param_counts(); // If possible convert this to protected/private and // and figure out how to declare cython fn // load_from_python as a friend. int model_loaded; protected: - virtual void read_coeffs(char *); + void read_coeffs(char *) override; private: }; diff --git a/src/ML-IAP/mliap_unified.cpp b/src/ML-IAP/mliap_unified.cpp new file mode 100644 index 0000000000..6bcdc88289 --- /dev/null +++ b/src/ML-IAP/mliap_unified.cpp @@ -0,0 +1,291 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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: Steven Ray Anaya (LANL) +------------------------------------------------------------------------- */ + +#ifdef MLIAP_PYTHON + +#include "mliap_unified.h" +#include + +#include "error.h" +#include "lmppython.h" +#include "memory.h" +#include "mliap_data.h" +#include "mliap_unified_couple.h" +#include "pair_mliap.h" +#include "python_compat.h" +#include "utils.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +MLIAPDummyDescriptor::MLIAPDummyDescriptor(LAMMPS *_lmp) : MLIAPDescriptor(_lmp) {} + +MLIAPDummyDescriptor::~MLIAPDummyDescriptor() +{ + // manually decrement borrowed reference from Python + Py_DECREF(unified_interface); +} + +/* ---------------------------------------------------------------------- + invoke compute_descriptors from Cython interface + ---------------------------------------------------------------------- */ + +void MLIAPDummyDescriptor::compute_descriptors(class MLIAPData *data) +{ + PyGILState_STATE gstate = PyGILState_Ensure(); + compute_descriptors_python(unified_interface, data); + if (PyErr_Occurred()) { + PyErr_Print(); + PyErr_Clear(); + PyGILState_Release(gstate); + lmp->error->all(FLERR, "Running mliappy unified compute_descriptors failure."); + } + PyGILState_Release(gstate); +} + +/* ---------------------------------------------------------------------- + invoke compute_forces from Cython interface + ---------------------------------------------------------------------- */ + +void MLIAPDummyDescriptor::compute_forces(class MLIAPData *data) +{ + PyGILState_STATE gstate = PyGILState_Ensure(); + compute_forces_python(unified_interface, data); + if (PyErr_Occurred()) { + PyErr_Print(); + PyErr_Clear(); + PyGILState_Release(gstate); + lmp->error->all(FLERR, "Running mliappy unified compute_forces failure."); + } + PyGILState_Release(gstate); +} + +// not implemented +void MLIAPDummyDescriptor::compute_force_gradients(class MLIAPData *) +{ + error->all(FLERR, "compute_force_gradients not implemented"); +} + +// not implemented +void MLIAPDummyDescriptor::compute_descriptor_gradients(class MLIAPData *) +{ + error->all(FLERR, "compute_descriptor_gradients not implemented"); +} + +void MLIAPDummyDescriptor::init() +{ + memory->create(radelem, nelements, "mliap_dummy_descriptor:radelem"); + for (int ielem = 0; ielem < nelements; ielem++) { radelem[ielem] = 1; } + + double cut; + cutmax = 0.0; + memory->create(cutsq, nelements, nelements, "mliap/descriptor/dummy:cutsq"); + memory->create(cutghost, nelements, nelements, "mliap/descriptor/dummy:cutghost"); + for (int ielem = 0; ielem < nelements; ielem++) { + // rcutfac set from python, is global cutoff for all elements + cut = 2.0 * radelem[ielem] * rcutfac; + if (cut > cutmax) cutmax = cut; + cutsq[ielem][ielem] = cut * cut; + cutghost[ielem][ielem] = cut * cut; + for (int jelem = ielem + 1; jelem < nelements; jelem++) { + cut = (radelem[ielem] + radelem[jelem]) * rcutfac; + cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut * cut; + cutghost[ielem][jelem] = cutghost[jelem][ielem] = cut * cut; + } + } +} + +void MLIAPDummyDescriptor::set_elements(char **elems, int nelems) +{ + nelements = nelems; + elements = new char *[nelems]; + for (int i = 0; i < nelems; i++) { elements[i] = utils::strdup(elems[i]); } +} + +/* ---------------------------------------------------------------------- */ + +MLIAPDummyModel::MLIAPDummyModel(LAMMPS *lmp, char *coefffilename) : MLIAPModel(lmp, coefffilename) +{ + nonlinearflag = 1; +} + +MLIAPDummyModel::~MLIAPDummyModel() +{ + // manually decrement borrowed reference from Python + Py_DECREF(unified_interface); +} + +int MLIAPDummyModel::get_nparams() +{ + return nparams; +} + +int MLIAPDummyModel::get_gamma_nnz(class MLIAPData *) +{ + // TODO: get_gamma_nnz + return 0; +} + +/* ---------------------------------------------------------------------- + invoke compute_gradients from Cython interface + ---------------------------------------------------------------------- */ + +void MLIAPDummyModel::compute_gradients(class MLIAPData *data) +{ + PyGILState_STATE gstate = PyGILState_Ensure(); + compute_gradients_python(unified_interface, data); + if (PyErr_Occurred()) { + PyErr_Print(); + PyErr_Clear(); + PyGILState_Release(gstate); + lmp->error->all(FLERR, "Running mliappy unified compute_gradients failure."); + } + PyGILState_Release(gstate); +} + +// not implemented +void MLIAPDummyModel::compute_gradgrads(class MLIAPData *) +{ + error->all(FLERR, "compute_gradgrads not implemented"); +} + +// not implemented +void MLIAPDummyModel::compute_force_gradients(class MLIAPData *) +{ + error->all(FLERR, "compute_force_gradients not implemented"); +} + +/* ---------------------------------------------------------------------- + memory usage unclear due to Cython/Python implementation + ---------------------------------------------------------------------- */ + +double MLIAPDummyModel::memory_usage() +{ + // TODO: implement memory usage in Cython(?) + return 0; +} + +// not implemented +void MLIAPDummyModel::read_coeffs(char *) +{ + error->all(FLERR, "read_coeffs not implemented"); +} + +/* ---------------------------------------------------------------------- + build the unified interface object, connect to dummy model and descriptor + ---------------------------------------------------------------------- */ + +MLIAPBuildUnified_t LAMMPS_NS::build_unified(char *unified_fname, MLIAPData *data, LAMMPS *lmp, + char *coefffilename) +{ + lmp->python->init(); + PyGILState_STATE gstate = PyGILState_Ensure(); + + PyObject *pyMain = PyImport_AddModule("__main__"); + + if (!pyMain) { + PyGILState_Release(gstate); + lmp->error->all(FLERR, "Could not initialize embedded Python"); + } + + PyObject *unified_module = PyImport_ImportModule("mliap_unified_couple"); + + if (PyErr_Occurred()) { + PyErr_Print(); + PyErr_Clear(); + PyGILState_Release(gstate); + lmp->error->all(FLERR, "Loading mliappy unified module failure."); + } + + // Connect dummy model, dummy descriptor, data to Python unified + MLIAPDummyModel *model = new MLIAPDummyModel(lmp, coefffilename); + MLIAPDummyDescriptor *descriptor = new MLIAPDummyDescriptor(lmp); + + PyObject *unified_interface = mliap_unified_connect(unified_fname, model, descriptor); + if (PyErr_Occurred()) { + PyErr_Print(); + PyErr_Clear(); + PyGILState_Release(gstate); + lmp->error->all(FLERR, "Running mliappy unified module failure."); + } + + // Borrowed references must be manually incremented + model->unified_interface = unified_interface; + Py_INCREF(unified_interface); + descriptor->unified_interface = unified_interface; + Py_INCREF(unified_interface); + + PyGILState_Release(gstate); + + MLIAPBuildUnified_t build = {data, descriptor, model}; + return build; +} + +/* ---------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + set energy for ij atom pairs + ---------------------------------------------------------------------- */ + +void LAMMPS_NS::update_pair_energy(MLIAPData *data, double *eij) +{ + double e_total = 0.0; + const auto nlistatoms = data->nlistatoms; + for (int ii = 0; ii < nlistatoms; ii++) data->eatoms[ii] = 0; + + for (int ii = 0; ii < data->npairs; ii++) { + int i = data->pair_i[ii]; + double e = 0.5 * eij[ii]; + + // must not count any contribution where i is not a local atom + if (i < nlistatoms) { + data->eatoms[i] += e; + e_total += e; + } + } + data->energy = e_total; +} + +/* ---------------------------------------------------------------------- + set forces for ij atom pairs + ---------------------------------------------------------------------- */ + +void LAMMPS_NS::update_pair_forces(MLIAPData *data, double *fij) +{ + const auto nlistatoms = data->nlistatoms; + double **f = data->f; + for (int ii = 0; ii < data->npairs; ii++) { + int ii3 = ii * 3; + int i = data->pair_i[ii]; + int j = data->jatoms[ii]; + + // must not count any contribution where i is not a local atom + if (i < nlistatoms) { + f[i][0] += fij[ii3]; + f[i][1] += fij[ii3 + 1]; + f[i][2] += fij[ii3 + 2]; + f[j][0] -= fij[ii3]; + f[j][1] -= fij[ii3 + 1]; + f[j][2] -= fij[ii3 + 2]; + + if (data->vflag) data->pairmliap->v_tally(i, j, &fij[ii3], data->rij[ii]); + } + } +} + +#endif diff --git a/src/ML-IAP/mliap_unified.h b/src/ML-IAP/mliap_unified.h new file mode 100644 index 0000000000..adc900d8c1 --- /dev/null +++ b/src/ML-IAP/mliap_unified.h @@ -0,0 +1,69 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 LMP_MLIAP_UNIFIED_H +#define LMP_MLIAP_UNIFIED_H + +#include "mliap_data.h" +#include "mliap_descriptor.h" +#include "mliap_model.h" + +#include + +namespace LAMMPS_NS { + +class MLIAPDummyDescriptor : public MLIAPDescriptor { + public: + MLIAPDummyDescriptor(LAMMPS *); + ~MLIAPDummyDescriptor() override; + void compute_descriptors(class MLIAPData *) override; + void compute_forces(class MLIAPData *) override; + void compute_force_gradients(class MLIAPData *) override; + void compute_descriptor_gradients(class MLIAPData *) override; + void init() override; + void set_elements(char **, int); + + PyObject *unified_interface; // MLIAPUnifiedInterface + double rcutfac; +}; + +class MLIAPDummyModel : public MLIAPModel { + public: + MLIAPDummyModel(LAMMPS *, char * = nullptr); + ~MLIAPDummyModel() override; + int get_nparams() override; + int get_gamma_nnz(class MLIAPData *) override; + void compute_gradients(class MLIAPData *) override; + void compute_gradgrads(class MLIAPData *) override; + void compute_force_gradients(class MLIAPData *) override; + double memory_usage() override; + + PyObject *unified_interface; + + protected: + void read_coeffs(char *) override; +}; + +struct MLIAPBuildUnified_t { + MLIAPData *data; + MLIAPDummyDescriptor *descriptor; + MLIAPDummyModel *model; +}; + +MLIAPBuildUnified_t build_unified(char *, MLIAPData *, LAMMPS *, char * = NULL); +void update_pair_energy(MLIAPData *, double *); +void update_pair_forces(MLIAPData *, double *); + +} // namespace LAMMPS_NS + +#endif diff --git a/src/ML-IAP/mliap_unified_couple.pyx b/src/ML-IAP/mliap_unified_couple.pyx new file mode 100644 index 0000000000..3fde99a25e --- /dev/null +++ b/src/ML-IAP/mliap_unified_couple.pyx @@ -0,0 +1,393 @@ +# cython: language_level=3 +# distutils: language = c++ + +import pickle +import numpy as np +import lammps.mliap + +cimport cython +from cpython.ref cimport PyObject +from libc.stdlib cimport malloc, free + + +cdef extern from "lammps.h" namespace "LAMMPS_NS": + cdef cppclass LAMMPS: + pass + + +cdef extern from "mliap_data.h" namespace "LAMMPS_NS": + cdef cppclass MLIAPData: + # ----- may not need ----- + int size_array_rows + int size_array_cols + int natoms + int yoffset + int zoffset + int ndims_force + int ndims_virial + # -END- may not need -END- + int size_gradforce + # ----- write only ----- + double ** f + double ** gradforce + double ** betas # betas for all atoms in list + double ** descriptors # descriptors for all atoms in list + double * eatoms # energies for all atoms in list + double energy + # -END- write only -END- + int ndescriptors # number of descriptors + int nparams # number of model parameters per element + int nelements # number of elements + + # data structures for grad-grad list (gamma) + + # ----- ignore for now ----- + int gamma_nnz # number of non-zero entries in gamma + double ** gamma # gamma element + int ** gamma_row_index # row (parameter) index + int ** gamma_col_index # column (descriptor) index + double * egradient # energy gradient w.r.t. parameters + # -END- ignore for now -END- + + # data structures for mliap neighbor list + # only neighbors strictly inside descriptor cutoff + + int ntotal # total number of owned and ghost atoms on this proc + int nlistatoms # current number of atoms in local atom lists + int natomneigh # current number of atoms and ghosts in atom neighbor arrays + int * numneighs # neighbors count for each atom + int * iatoms # index of each atom + int * pair_i # index of each i atom for each ij pair + int * ielems # element of each atom + int nneigh_max # number of ij neighbors allocated + int npairs # number of ij neighbor pairs + int * jatoms # index of each neighbor + int * jelems # element of each neighbor + int * elems # element of each atom in or not in the neighborlist + double ** rij # distance vector of each neighbor + # ----- write only ----- + double *** graddesc # descriptor gradient w.r.t. each neighbor + # -END- write only -END- + int eflag # indicates if energy is needed + int vflag # indicates if virial is needed + + +cdef extern from "mliap_unified.h" namespace "LAMMPS_NS": + cdef cppclass MLIAPDummyDescriptor: + MLIAPDummyDescriptor(PyObject *, LAMMPS *) except + + int ndescriptors # number of descriptors + int nelements # # of unique elements + char **elements # names of unique elements + double cutmax # maximum cutoff needed + double rcutfac + double *radelem # element radii + + void compute_descriptors(MLIAPData *) + void compute_forces(MLIAPData *) + void set_elements(char **, int) + + cdef cppclass MLIAPDummyModel: + MLIAPDummyModel(PyObject *, LAMMPS *, char * = NULL) except + + int ndescriptors # number of descriptors + int nparams # number of parameters per element + int nelements; # # of unique elements + + void compute_gradients(MLIAPData *) + + cdef void update_pair_energy(MLIAPData *, double *) except + + cdef void update_pair_forces(MLIAPData *, double *) except + + + +LOADED_MODEL = None + + +# @property sans getter +def write_only_property(fset): + return property(fget=None, fset=fset) + + +# Cython implementation of MLIAPData +# Automatically converts between C arrays and numpy when needed +cdef class MLIAPDataPy: + cdef MLIAPData * data + + def __cinit__(self): + self.data = NULL + + def update_pair_energy(self, eij): + cdef double[:] eij_arr = eij + update_pair_energy(self.data, &eij_arr[0]) + + def update_pair_forces(self, fij): + cdef double[:, ::1] fij_arr = fij + update_pair_forces(self.data, &fij_arr[0][0]) + + @property + def f(self): + if self.data.f is NULL: + return None + return np.asarray( &self.data.f[0][0]) + + @property + def size_gradforce(self): + return self.data.size_gradforce + + @write_only_property + def gradforce(self, value): + if self.data.gradforce is NULL: + raise ValueError("attempt to set NULL gradforce") + cdef double[:, :] gradforce_view = &self.data.gradforce[0][0] + cdef double[:, :] value_view = value + gradforce_view[:] = value_view + + @write_only_property + def betas(self, value): + if self.data.betas is NULL: + raise ValueError("attempt to set NULL betas") + cdef double[:, :] betas_view = &self.data.betas[0][0] + cdef double[:, :] value_view = value + betas_view[:] = value_view + + @write_only_property + def descriptors(self, value): + if self.data.descriptors is NULL: + raise ValueError("attempt to set NULL descriptors") + cdef double[:, :] descriptors_view = &self.data.descriptors[0][0] + cdef double[:, :] value_view = value + descriptors_view[:] = value_view + + @write_only_property + def eatoms(self, value): + if self.data.eatoms is NULL: + raise ValueError("attempt to set NULL eatoms") + cdef double[:] eatoms_view = &self.data.eatoms[0] + cdef double[:] value_view = value + eatoms_view[:] = value_view + + @write_only_property + def energy(self, value): + self.data.energy = value + + @property + def ndescriptors(self): + return self.data.ndescriptors + + @property + def nparams(self): + return self.data.nparams + + @property + def nelements(self): + return self.data.nelements + + # data structures for grad-grad list (gamma) + + @property + def gamma_nnz(self): + return self.data.gamma_nnz + + @property + def gamma(self): + if self.data.gamma is NULL: + return None + return np.asarray( &self.data.gamma[0][0]) + + @property + def gamma_row_index(self): + if self.data.gamma_row_index is NULL: + return None + return np.asarray( &self.data.gamma_row_index[0][0]) + + @property + def gamma_col_index(self): + if self.data.gamma_col_index is NULL: + return None + return np.asarray( &self.data.gamma_col_index[0][0]) + + @property + def egradient(self): + if self.data.egradient is NULL: + return None + return np.asarray( &self.data.egradient[0]) + + # data structures for mliap neighbor list + # only neighbors strictly inside descriptor cutoff + + @property + def ntotal(self): + return self.data.ntotal + + @property + def elems(self): + if self.data.elems is NULL: + return None + return np.asarray( &self.data.elems[0]) + + @property + def nlistatoms(self): + return self.data.nlistatoms + + @property + def natomneigh(self): + return self.data.natomneigh + + @property + def numneighs(self): + if self.data.numneighs is NULL: + return None + return np.asarray( &self.data.numneighs[0]) + + @property + def iatoms(self): + if self.data.iatoms is NULL: + return None + return np.asarray( &self.data.iatoms[0]) + + @property + def ielems(self): + if self.data.ielems is NULL: + return None + return np.asarray( &self.data.ielems[0]) + + @property + def npairs(self): + return self.data.npairs + + @property + def pair_i(self): + if self.data.pair_i is NULL: + return None + return np.asarray( &self.data.pair_i[0]) + + @property + def pair_j(self): + return self.jatoms + + @property + def jatoms(self): + if self.data.jatoms is NULL: + return None + return np.asarray( &self.data.jatoms[0]) + + @property + def jelems(self): + if self.data.jelems is NULL: + return None + return np.asarray( &self.data.jelems[0]) + + @property + def rij(self): + if self.data.rij is NULL: + return None + return np.asarray( &self.data.rij[0][0]) + + @write_only_property + def graddesc(self, value): + if self.data.graddesc is NULL: + raise ValueError("attempt to set NULL graddesc") + cdef double[:, :, :] graddesc_view = &self.data.graddesc[0][0][0] + cdef double[:, :, :] value_view = value + graddesc_view[:] = value_view + + @property + def eflag(self): + return self.data.eflag + + @property + def vflag(self): + return self.data.vflag + + +# Interface between C and Python compute functions +cdef class MLIAPUnifiedInterface: + cdef MLIAPDummyModel * model + cdef MLIAPDummyDescriptor * descriptor + cdef unified_impl + + def __init__(self, unified_impl): + self.model = NULL + self.descriptor = NULL + self.unified_impl = unified_impl + + def compute_gradients(self, data): + self.unified_impl.compute_gradients(data) + + def compute_descriptors(self, data): + self.unified_impl.compute_descriptors(data) + + def compute_forces(self, data): + self.unified_impl.compute_forces(data) + + +cdef public void compute_gradients_python(unified_int, MLIAPData *data) except * with gil: + pydata = MLIAPDataPy() + pydata.data = data + unified_int.compute_gradients(pydata) + + +cdef public void compute_descriptors_python(unified_int, MLIAPData *data) except * with gil: + pydata = MLIAPDataPy() + pydata.data = data + unified_int.compute_descriptors(pydata) + + +cdef public void compute_forces_python(unified_int, MLIAPData *data) except * with gil: + pydata = MLIAPDataPy() + pydata.data = data + unified_int.compute_forces(pydata) + + +# Create a MLIAPUnifiedInterface and connect it to the dummy model, descriptor +cdef public object mliap_unified_connect(char *fname, MLIAPDummyModel * model, + MLIAPDummyDescriptor * descriptor) with gil: + str_fname = fname.decode('utf-8') + if str_fname == 'EXISTS': + if LOADED_MODEL is None: + raise ValueError("No unified model loaded") + unified = LOADED_MODEL + elif str_fname.endswith(".pt") or str_fname.endswith('.pth'): + import torch + unified = torch.load(str_fname) + else: + with open(str_fname, 'rb') as pfile: + unified = pickle.load(pfile) + + unified_int = MLIAPUnifiedInterface(unified) + unified_int.model = model + unified_int.descriptor = descriptor + + unified.interface = unified_int + + if unified.ndescriptors is None: + raise ValueError("no descriptors set") + + unified_int.descriptor.ndescriptors = unified.ndescriptors + unified_int.descriptor.rcutfac = unified.rcutfac + unified_int.model.ndescriptors = unified.ndescriptors + unified_int.model.nparams = unified.nparams + + if unified.element_types is None: + raise ValueError("no element type set") + + cdef int nelements = len(unified.element_types) + cdef char **elements = malloc(nelements * sizeof(char*)) + + if not elements: + raise MemoryError("failed to allocate memory for element names") + + cdef char *elem_name + for i, elem in enumerate(unified.element_types): + elem_name_bytes = elem.encode('UTF-8') + elem_name = elem_name_bytes + elements[i] = &elem_name[0] + unified_int.descriptor.set_elements(elements, nelements) + unified_int.model.nelements = nelements + + free(elements) + return unified_int + + +# For pre-loading a Python model +def load_from_python(unified): + global LOADED_MODEL + LOADED_MODEL = unified diff --git a/src/ML-IAP/pair_mliap.cpp b/src/ML-IAP/pair_mliap.cpp index 179649caff..af8768568e 100644 --- a/src/ML-IAP/pair_mliap.cpp +++ b/src/ML-IAP/pair_mliap.cpp @@ -25,14 +25,17 @@ #include "mliap_model_nn.h" #include "mliap_model_quadratic.h" #ifdef MLIAP_PYTHON +#include "mliap_unified.h" #include "mliap_model_python.h" #endif #include "atom.h" +#include "comm.h" #include "error.h" #include "force.h" #include "memory.h" #include "neighbor.h" +#include "neigh_request.h" #include #include @@ -64,6 +67,7 @@ PairMLIAP::~PairMLIAP() if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); + memory->destroy(cutghost); memory->destroy(map); } } @@ -78,10 +82,12 @@ void PairMLIAP::compute(int eflag, int vflag) // consistency checks if (data->ndescriptors != model->ndescriptors) - error->all(FLERR, "Incompatible model and descriptor descriptor count"); + error->all(FLERR, "Inconsistent model and descriptor descriptor count: {} vs {}", + model->ndescriptors, data->ndescriptors); if (data->nelements != model->nelements) - error->all(FLERR, "Incompatible model and descriptor element count"); + error->all(FLERR, "Inconsistent model and descriptor element count: {} vs {}", + model->nelements, data->nelements); ev_init(eflag, vflag); data->generate_neighdata(list, eflag, vflag); @@ -93,11 +99,11 @@ void PairMLIAP::compute(int eflag, int vflag) // compute E_i and beta_i = dE_i/dB_i for all i in list model->compute_gradients(data); - e_tally(data); // calculate force contributions beta_i*dB_i/dR_j descriptor->compute_forces(data); + e_tally(data); // calculate stress @@ -115,6 +121,7 @@ void PairMLIAP::allocate() memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutghost,n+1,n+1,"pair:cutghost"); memory->create(map,n+1,"pair:map"); } @@ -124,8 +131,7 @@ void PairMLIAP::allocate() void PairMLIAP::settings(int narg, char ** arg) { - if (narg < 4) - error->all(FLERR,"Illegal pair_style command"); + if (narg < 2) utils::missing_cmd_args(FLERR, "pair_style mliap", error); // set flags for required keywords @@ -140,47 +146,70 @@ void PairMLIAP::settings(int narg, char ** arg) while (iarg < narg) { if (strcmp(arg[iarg],"model") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model", error); + if (modelflag) error->all(FLERR,"Illegal multiple pair_style mliap model definition"); if (strcmp(arg[iarg+1],"linear") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model linear", error); model = new MLIAPModelLinear(lmp,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"quadratic") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model quadratic", error); model = new MLIAPModelQuadratic(lmp,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"nn") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model nn", error); model = new MLIAPModelNN(lmp,arg[iarg+2]); iarg += 3; -#ifdef MLIAP_PYTHON } else if (strcmp(arg[iarg+1],"mliappy") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); +#ifdef MLIAP_PYTHON + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap mliappy", error); model = new MLIAPModelPython(lmp,arg[iarg+2]); iarg += 3; +#else + error->all(FLERR,"Using pair_style mliap model mliappy requires ML-IAP with python support"); #endif - } else error->all(FLERR,"Illegal pair_style mliap command"); + } else error->all(FLERR,"Unknown pair_style mliap model keyword: {}", arg[iarg]); modelflag = 1; } else if (strcmp(arg[iarg],"descriptor") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor", error); + if (descriptorflag) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definition"); if (strcmp(arg[iarg+1],"sna") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor sna", error); descriptor = new MLIAPDescriptorSNAP(lmp,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg+1],"so3") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor so3", error); descriptor = new MLIAPDescriptorSO3(lmp,arg[iarg+2]); iarg += 3; } else error->all(FLERR,"Illegal pair_style mliap command"); descriptorflag = 1; + } else if (strcmp(arg[iarg], "unified") == 0) { +#ifdef MLIAP_PYTHON + if (modelflag) error->all(FLERR,"Illegal multiple pair_style mliap model definitions"); + if (descriptorflag) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definitions"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap unified", error); + MLIAPBuildUnified_t build = build_unified(arg[iarg+1], data, lmp); + if (iarg+3 > narg) { + ghostneigh = 0; + } else { + ghostneigh = utils::logical(FLERR, arg[iarg+2], false, lmp); + } + + iarg += 3; + model = build.model; + descriptor = build.descriptor; + modelflag = 1; + descriptorflag = 1; +#else + error->all(FLERR,"Using pair_style mliap unified requires ML-IAP with python support"); +#endif } else - error->all(FLERR,"Illegal pair_style mliap command"); + error->all(FLERR,"Unknown pair_style mliap keyword: {}", arg[iarg]); } if (modelflag == 0 || descriptorflag == 0) - error->all(FLERR,"Illegal pair_style command"); - + error->all(FLERR,"Incomplete pair_style mliap setup: need model and descriptor, or unified"); } /* ---------------------------------------------------------------------- @@ -196,11 +225,6 @@ void PairMLIAP::coeff(int narg, char **arg) char* type2 = arg[1]; char** elemtypes = &arg[2]; - // insure I,J args are * * - - if (strcmp(type1,"*") != 0 || strcmp(type2,"*") != 0) - error->all(FLERR,"Incorrect args for pair coefficients"); - // read args that map atom types to elements // map[i] = which element the Ith atom type is, -1 if not mapped // map[0] is not used @@ -313,7 +337,11 @@ void PairMLIAP::init_style() // need a full neighbor list - neighbor->add_request(this, NeighConst::REQ_FULL); + if (ghostneigh == 1) { + neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_GHOST); + } else { + neighbor->add_request(this, NeighConst::REQ_FULL); + } } @@ -324,7 +352,9 @@ void PairMLIAP::init_style() double PairMLIAP::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); - return sqrt(descriptor->cutsq[map[i]][map[j]]); + double cutmax = sqrt(descriptor->cutsq[map[i]][map[j]]); + cutghost[i][j] = cutghost[j][i] = 2.0 * cutmax + neighbor->skin; + return cutmax; } /* ---------------------------------------------------------------------- @@ -338,6 +368,7 @@ double PairMLIAP::memory_usage() int n = atom->ntypes+1; bytes += (double)n*n*sizeof(int); // setflag bytes += (double)n*n*sizeof(int); // cutsq + bytes += (double)n*n*sizeof(int); // cutghost bytes += (double)n*sizeof(int); // map bytes += descriptor->memory_usage(); // Descriptor object bytes += model->memory_usage(); // Model object diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index 5a7098814b..aba8e88e10 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -29,10 +29,12 @@ #ifdef MLIAP_PYTHON #include "mliap_model_python.h" +#include "mliap_unified.h" // The above should somehow really be included in the next file. // We could get around this with cython --capi-reexport-cincludes // However, that exposes -too many- headers. #include "mliap_model_python_couple.h" +#include "mliap_unified_couple.h" #endif using namespace LAMMPS_NS; @@ -66,6 +68,9 @@ PythonImpl::PythonImpl(LAMMPS *lmp) : Pointers(lmp) // This -must- happen before python is initialized. int err = PyImport_AppendInittab("mliap_model_python_couple", PyInit_mliap_model_python_couple); if (err) error->all(FLERR, "Could not register MLIAPPY embedded python module."); + + err = PyImport_AppendInittab("mliap_unified_couple", PyInit_mliap_unified_couple); + if (err) error->all(FLERR, "Could not register MLIAPPY unified embedded python module."); #endif Py_Initialize(); diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp index ba553335bb..b674598265 100644 --- a/src/SPIN/compute_spin.cpp +++ b/src/SPIN/compute_spin.cpp @@ -43,7 +43,7 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), pair(nullptr), spin_pairs(nullptr) + Compute(lmp, narg, arg), lockprecessionspin(nullptr), pair(nullptr), spin_pairs(nullptr) { if ((narg != 3) && (narg != 4)) error->all(FLERR,"Illegal compute compute/spin command"); @@ -68,7 +68,8 @@ ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) : ComputeSpin::~ComputeSpin() { memory->destroy(vector); - delete [] spin_pairs; + delete[] spin_pairs; + delete[] lockprecessionspin; } /* ---------------------------------------------------------------------- */ @@ -96,7 +97,7 @@ void ComputeSpin::init() else npairs = hybrid->nstyles; for (int i = 0; ipair_match("^spin",0,i)) { - npairspin ++; + npairspin++; } } } @@ -135,14 +136,18 @@ void ComputeSpin::init() } } - // ptrs FixPrecessionSpin classes + // set ptrs for fix precession/spin styles - int iforce; - for (iforce = 0; iforce < modify->nfix; iforce++) { - if (utils::strmatch(modify->fix[iforce]->style,"^precession/spin")) { - precession_spin_flag = 1; - lockprecessionspin = dynamic_cast(modify->fix[iforce]); - } + auto precfixes = modify->get_fix_by_style("^precession/spin"); + nprecspin = precfixes.size(); + + if (nprecspin > 0) { + lockprecessionspin = new FixPrecessionSpin *[nprecspin]; + precession_spin_flag = 1; + + int i = 0; + for (auto &ifix : precfixes) + lockprecessionspin[i++] = dynamic_cast(ifix); } } @@ -191,7 +196,9 @@ void ComputeSpin::compute_vector() // update magnetic precession energies if (precession_spin_flag) { - magenergy += lockprecessionspin->emag[i]; + for (int k = 0; k < nprecspin; k++) { + magenergy += lockprecessionspin[k]->emag[i]; + } } // update magnetic pair interactions diff --git a/src/SPIN/compute_spin.h b/src/SPIN/compute_spin.h index 1bcebfec67..45aac3c89f 100644 --- a/src/SPIN/compute_spin.h +++ b/src/SPIN/compute_spin.h @@ -40,7 +40,8 @@ class ComputeSpin : public Compute { // pointers to magnetic fixes - class FixPrecessionSpin *lockprecessionspin; + int nprecspin; + class FixPrecessionSpin **lockprecessionspin; // pointers to magnetic pair styles diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp index f8aacc0c5c..397c3c618f 100644 --- a/src/SPIN/fix_langevin_spin.cpp +++ b/src/SPIN/fix_langevin_spin.cpp @@ -109,7 +109,7 @@ void FixLangevinSpin::init() double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) double kb = force->boltz; // eV/K - D = (alpha_t*gil_factor*kb*temp); + D = (alpha_t*(1.0+(alpha_t)*(alpha_t))*kb*temp); D /= (hbar*dts); sigma = sqrt(2.0*D); } diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index e2b0c1e422..f6ed867c63 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -152,7 +152,6 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]); pack_choice[i] = &ComputePropertyAtom::pack_spx; } else if (strcmp(arg[iarg],"spy") == 0) { - error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]); if (!atom->sp_flag) error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]); pack_choice[i] = &ComputePropertyAtom::pack_spy; diff --git a/src/library.cpp b/src/library.cpp index fc3700c786..16381a089d 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -57,6 +57,10 @@ #include "exceptions.h" #endif +#if defined(LMP_PYTHON) +#include +#endif + using namespace LAMMPS_NS; // for printing the non-null pointer argument warning only once @@ -416,6 +420,89 @@ void lammps_python_finalize() Python::finalize(); } + +/* ---------------------------------------------------------------------- */ + +/** Call a LAMMPS Error class function + * +\verbatim embed:rst + +This function is a wrapper around functions in the ``Error`` to print an +error message and then stop LAMMPS. + +The *error_type* parameter selects which function to call. It is a sum +of constants from :cpp:enum:`_LMP_ERROR_CONST`. If the value does not +match any valid combination of constants a warning is printed and the +function returns. + +.. versionadded:: TBD + +\endverbatim + * + * \param handle pointer to a previously created LAMMPS instance + * \param error_type parameter to select function in the Error class + * \param error_text error message */ + +void lammps_error(void *handle, int error_type, const char *error_text) +{ + auto lmp = (LAMMPS *) handle; + + BEGIN_CAPTURE + { + switch (error_type) { + case LMP_ERROR_WARNING: + lmp->error->warning("(library)", 0, error_text); + break; + case LMP_ERROR_ONE: + lmp->error->one("(library)", 0, error_text); + break; + case LMP_ERROR_ALL: + lmp->error->all("(library)", 0, error_text); + break; + case LMP_ERROR_WARNING|LMP_ERROR_WORLD: + lmp->error->warning("(library)", 0, error_text); + break; + case LMP_ERROR_ONE|LMP_ERROR_WORLD: + lmp->error->one("(library)", 0, error_text); + break; + case LMP_ERROR_ALL|LMP_ERROR_WORLD: + lmp->error->all("(library)", 0, error_text); + break; + case LMP_ERROR_WARNING|LMP_ERROR_UNIVERSE: + lmp->error->universe_warn("(library)", 0, error_text); + break; + case LMP_ERROR_ONE|LMP_ERROR_UNIVERSE: + lmp->error->universe_one("(library)", 0, error_text); + break; + case LMP_ERROR_ALL|LMP_ERROR_UNIVERSE: + lmp->error->universe_all("(library)", 0, error_text); + break; + default: + auto mesg = fmt::format("Unknown error type {} for message: {}", error_type, error_text); + lmp->error->warning("(library)", 0, mesg); + } + } + END_CAPTURE + +#if defined(LAMMPS_EXCEPTIONS) + // with enabled exceptions the above code will simply throw an + // exception and record the error message. So we have to explicitly + // stop here like we do in main.cpp + if (lammps_has_error(handle)) { + if (error_type & 1) { + lammps_kokkos_finalize(); + lammps_python_finalize(); + MPI_Abort(lmp->universe->uworld, 1); + } else if (error_type & 2) { + lammps_kokkos_finalize(); + lammps_python_finalize(); + lammps_mpi_finalize(); + exit(1); + } + } +#endif +} + // ---------------------------------------------------------------------- // Library functions to process commands // ---------------------------------------------------------------------- @@ -1156,8 +1243,8 @@ Please also see :cpp:func:`lammps_extract_setting`, may lead to inconsistent internal data and thus may cause failures or crashes or bogus simulations. In general it is thus usually better to use a LAMMPS input command that sets or changes these parameters. - Those will takes care of all side effects and necessary updates of - settings derived from such settings. Where possible a reference to + Those will take care of all side effects and necessary updates of + settings derived from such settings. Where possible, a reference to such a command or a relevant section of the manual is given below. The following tables list the supported names, their data types, length @@ -5571,6 +5658,29 @@ int lammps_get_last_error_message(void *handle, char *buffer, int buf_size) { return 0; } +/* ---------------------------------------------------------------------- */ + +/** Return API version of embedded Python interpreter + +\verbatim embed:rst +This function is used by the ML-IAP python code (mliappy) to verify +the API version of the embedded python interpreter of the PYTHON +package. It returns -1 if the PYTHON package is not enabled. + +.. versionadded:: TBD + +\endverbatim + * + * \return PYTHON_API_VERSION constant of the python interpreter or -1 */ + +int lammps_python_api_version() { +#if defined(LMP_PYTHON) + return PYTHON_API_VERSION; +#else + return -1; +#endif +} + // Local Variables: // fill-column: 72 // End: diff --git a/src/library.h b/src/library.h index 311219e5ba..1eec57898e 100644 --- a/src/library.h +++ b/src/library.h @@ -75,6 +75,18 @@ enum _LMP_TYPE_CONST { LMP_SIZE_COLS = 5 /*!< return number of columns */ }; +/** Error codes to select the suitable function in the Error class + * + * Must be kept in sync with the equivalent constants in lammps/constants.py */ + +enum _LMP_ERROR_CONST { + LMP_ERROR_WARNING = 0, /*!< call Error::warning() */ + LMP_ERROR_ONE = 1, /*!< called from one MPI rank */ + LMP_ERROR_ALL = 2, /*!< called from all MPI ranks */ + LMP_ERROR_WORLD = 4, /*!< error on Comm::world */ + LMP_ERROR_UNIVERSE = 8 /*!< error on Comm::universe */ +}; + /* Ifdefs to allow this file to be included in C and C++ programs */ #ifdef __cplusplus @@ -97,6 +109,8 @@ void lammps_mpi_finalize(); void lammps_kokkos_finalize(); void lammps_python_finalize(); +void lammps_error(void *handle, int error_type, const char *error_text); + /* ---------------------------------------------------------------------- * Library functions to process commands * ---------------------------------------------------------------------- */ @@ -256,6 +270,8 @@ void lammps_force_timeout(void *handle); int lammps_has_error(void *handle); int lammps_get_last_error_message(void *handle, char *buffer, int buf_size); +int lammps_python_api_version(); + #ifdef __cplusplus } #endif diff --git a/unittest/c-library/test_library_open.cpp b/unittest/c-library/test_library_open.cpp index c75fd3738b..1cd65d843d 100644 --- a/unittest/c-library/test_library_open.cpp +++ b/unittest/c-library/test_library_open.cpp @@ -198,3 +198,32 @@ TEST(lammps_open_fortran, no_args) if (verbose) std::cout << output; MPI_Comm_free(&mycomm); } + +TEST(lammps_open_no_mpi, lammps_error) +{ + const char *args[] = {"liblammps", "-log", "none", "-nocite"}; + char **argv = (char **)args; + int argc = sizeof(args) / sizeof(char *); + + ::testing::internal::CaptureStdout(); + void *alt_ptr; + void *handle = lammps_open_no_mpi(argc, argv, &alt_ptr); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_EQ(handle, alt_ptr); + LAMMPS_NS::LAMMPS *lmp = (LAMMPS_NS::LAMMPS *)handle; + + EXPECT_EQ(lmp->world, MPI_COMM_WORLD); + EXPECT_EQ(lmp->infile, stdin); + EXPECT_NE(lmp->screen, nullptr); + EXPECT_EQ(lmp->logfile, nullptr); + EXPECT_EQ(lmp->citeme, nullptr); + EXPECT_EQ(lmp->suffix_enable, 0); + + EXPECT_STREQ(lmp->exename, "liblammps"); + EXPECT_EQ(lmp->num_package, 0); + ::testing::internal::CaptureStdout(); + lammps_error(handle, 0, "test_warning"); + output = ::testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, HasSubstr("WARNING: test_warning")); +} + diff --git a/unittest/force-styles/CMakeLists.txt b/unittest/force-styles/CMakeLists.txt index d2345ceaa7..294ece3a1f 100644 --- a/unittest/force-styles/CMakeLists.txt +++ b/unittest/force-styles/CMakeLists.txt @@ -79,6 +79,7 @@ else() set(FORCE_TEST_ENVIRONMENT PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}) endif() list(APPEND FORCE_TEST_ENVIRONMENT "PYTHONUNBUFFERED=1") +list(APPEND FORCE_TEST_ENVIRONMENT "PYTHONDONTWRITEBYTECODE=1") list(APPEND FORCE_TEST_ENVIRONMENT "OMP_PROC_BIND=false") list(APPEND FORCE_TEST_ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}") get_property(BUILD_IS_MULTI_CONFIG GLOBAL PROPERTY GENERATOR_IS_MULTI_CONFIG) @@ -213,7 +214,7 @@ foreach(TEST ${FIX_TIMESTEP_TESTS}) if(WIN32) set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}\\\;${LAMMPS_PYTHON_DIR}") else() - set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH}") + set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:${LAMMPS_PYTHON_DIR}:$ENV{PYTHONPATH};PYTHONDONTWRITEBYTECODE=1") endif() set_tests_properties(${TNAME} PROPERTIES LABELS "${TEST_TAGS}") endforeach() @@ -231,7 +232,7 @@ foreach(TEST ${DIHEDRAL_TESTS}) continue() endif() add_test(NAME ${TNAME} COMMAND test_dihedral_style ${TEST}) - set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}") + set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH};PYTHONDONTWRITEBYTECODE=1") set_tests_properties(${TNAME} PROPERTIES LABELS "${TEST_TAGS}") endforeach() @@ -248,6 +249,14 @@ foreach(TEST ${IMPROPER_TESTS}) continue() endif() add_test(NAME ${TNAME} COMMAND test_improper_style ${TEST}) - set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH}") + set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR};PYTHONPATH=${TEST_INPUT_FOLDER}:$ENV{PYTHONPATH};PYTHONDONTWRITEBYTECODE=1") set_tests_properties(${TNAME} PROPERTIES LABELS "${TEST_TAGS}") endforeach() + +if(MLIAP_ENABLE_PYTHON AND (NOT WIN32)) + add_executable(test_mliappy_unified test_mliappy_unified.cpp) + target_link_libraries(test_mliappy_unified PRIVATE lammps GTest::GMockMain) + add_test(NAME TestMliapPyUnified COMMAND test_mliappy_unified) + set_tests_properties(TestMliapPyUnified PROPERTIES ENVIRONMENT "PYTHONPATH=${LAMMPS_PYTHON_DIR};PYTHONDONTWRITEBYTECODE=1") +endif() + diff --git a/unittest/force-styles/test_mliappy_unified.cpp b/unittest/force-styles/test_mliappy_unified.cpp new file mode 100644 index 0000000000..1eb800a954 --- /dev/null +++ b/unittest/force-styles/test_mliappy_unified.cpp @@ -0,0 +1,116 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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. +------------------------------------------------------------------------- */ + +#include "library.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +const char pickle[] = "python create_pickle here \"\"\"\n" + "import lammps\n" + "import lammps.mliap\n" + "from lammps.mliap.mliap_unified_lj import MLIAPUnifiedLJ\n" + "def create_pickle():\n" + " unified = MLIAPUnifiedLJ(['Ar'])\n" + " unified.pickle('mliap_unified_lj_Ar.pkl')\n" + "\"\"\"\n"; + +const char first[] = "units lj\n" + "atom_style atomic\n" + "lattice fcc 0.8442\n" + "region box block 0 2 0 2 0 2\n" + "create_box 1 box\n" + "create_atoms 1 box\n" + "mass 1 1.0\n" + "velocity all create 3.0 87287 loop geom\n"; + +const char second[] = "neighbor 0.3 bin\n" + "neigh_modify every 20 delay 0 check no\n" + "fix 1 all nve\n" + "run 2 post no\n"; + +namespace LAMMPS_NS { + +TEST(MliapUnified, VersusLJMelt) +{ + const char *lmpargv[] = {"melt", "-log", "none", "-nocite"}; + int lmpargc = sizeof(lmpargv) / sizeof(const char *); + + void *ljmelt = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr); + void *mliap = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr); + + lammps_commands_string(ljmelt, first); + lammps_command(ljmelt, "pair_style lj/cut 2.5"); + lammps_command(ljmelt, "pair_coeff * * 1.0 1.0"); + lammps_commands_string(ljmelt, second); + + lammps_command(mliap, pickle); + lammps_command(mliap, "python create_pickle invoke"); + + lammps_commands_string(mliap, first); + lammps_command(mliap, "pair_style mliap unified mliap_unified_lj_Ar.pkl 0"); + lammps_command(mliap, "pair_coeff * * Ar"); + lammps_commands_string(mliap, second); + + double lj_pe = lammps_get_thermo(ljmelt, "pe"); + double ml_pe = lammps_get_thermo(mliap, "pe"); + EXPECT_NEAR(lj_pe, ml_pe, 1.0e-14); + double lj_ke = lammps_get_thermo(ljmelt, "ke"); + double ml_ke = lammps_get_thermo(mliap, "ke"); + EXPECT_NEAR(lj_ke, ml_ke, 1.0e-14); + double lj_press = lammps_get_thermo(ljmelt, "press"); + double ml_press = lammps_get_thermo(mliap, "press"); + EXPECT_NEAR(lj_press, ml_press, 1.0e-14); + + lammps_command(mliap, "shell rm mliap_unified_lj_Ar.pkl"); + lammps_close(ljmelt); + lammps_close(mliap); +} + +TEST(MliapUnified, VersusLJMeltGhost) +{ + const char *lmpargv[] = {"melt", "-log", "none", "-nocite"}; + int lmpargc = sizeof(lmpargv) / sizeof(const char *); + + void *ljmelt = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr); + void *mliap = lammps_open_no_mpi(lmpargc, (char **)lmpargv, nullptr); + + lammps_commands_string(ljmelt, first); + lammps_command(ljmelt, "pair_style lj/cut 2.5"); + lammps_command(ljmelt, "pair_coeff * * 1.0 1.0"); + lammps_commands_string(ljmelt, second); + + lammps_command(mliap, pickle); + lammps_command(mliap, "python create_pickle invoke"); + + lammps_commands_string(mliap, first); + lammps_command(mliap, "pair_style mliap unified mliap_unified_lj_Ar.pkl 1"); + lammps_command(mliap, "pair_coeff * * Ar"); + lammps_commands_string(mliap, second); + + double lj_pe = lammps_get_thermo(ljmelt, "pe"); + double ml_pe = lammps_get_thermo(mliap, "pe"); + EXPECT_NEAR(lj_pe, ml_pe, 1.0e-14); + double lj_ke = lammps_get_thermo(ljmelt, "ke"); + double ml_ke = lammps_get_thermo(mliap, "ke"); + EXPECT_NEAR(lj_ke, ml_ke, 1.0e-14); + double lj_press = lammps_get_thermo(ljmelt, "press"); + double ml_press = lammps_get_thermo(mliap, "press"); + EXPECT_NEAR(lj_press, ml_press, 1.0e-14); + + lammps_command(mliap, "shell rm mliap_unified_lj_Ar.pkl"); + lammps_close(ljmelt); + lammps_close(mliap); +} + +} // namespace LAMMPS_NS diff --git a/unittest/fortran/CMakeLists.txt b/unittest/fortran/CMakeLists.txt index 9e46cef23e..d1f632f959 100644 --- a/unittest/fortran/CMakeLists.txt +++ b/unittest/fortran/CMakeLists.txt @@ -40,6 +40,23 @@ if(CMAKE_Fortran_COMPILER) add_executable(test_fortran_commands wrap_commands.cpp test_fortran_commands.f90) target_link_libraries(test_fortran_commands PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain) add_test(NAME FortranCommands COMMAND test_fortran_commands) + + add_executable(test_fortran_get_thermo wrap_get_thermo.cpp test_fortran_get_thermo.f90) + target_link_libraries(test_fortran_get_thermo PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain) + add_test(NAME FortranGetThermo COMMAND test_fortran_get_thermo) + + add_executable(test_fortran_box wrap_box.cpp test_fortran_box.f90) + target_link_libraries(test_fortran_box PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain) + add_test(NAME FortranBox COMMAND test_fortran_box) + + add_executable(test_fortran_properties wrap_properties.cpp test_fortran_properties.f90 test_fortran_commands.f90) + target_link_libraries(test_fortran_properties PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain) + add_test(NAME FortranProperties COMMAND test_fortran_properties) + + add_executable(test_fortran_extract_global wrap_extract_global.cpp test_fortran_extract_global.f90) + target_link_libraries(test_fortran_extract_global PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain) + add_test(NAME FortranExtractGlobal COMMAND test_fortran_extract_global) + else() message(STATUS "Skipping Tests for the LAMMPS Fortran Module: no Fortran compiler") endif() diff --git a/unittest/fortran/test_fortran_box.f90 b/unittest/fortran/test_fortran_box.f90 new file mode 100644 index 0000000000..c21918afd0 --- /dev/null +++ b/unittest/fortran/test_fortran_box.f90 @@ -0,0 +1,142 @@ +MODULE keepbox + USE liblammps + IMPLICIT NONE + TYPE(LAMMPS) :: lmp + CHARACTER(len=40), DIMENSION(3), PARAMETER :: demo_input = & + [ CHARACTER(len=40) :: & + 'region box block 0 $x 0 2 0 2', & + 'create_box 1 box', & + 'create_atoms 1 single 1.0 1.0 ${zpos}' ] + CHARACTER(len=40), DIMENSION(2), PARAMETER :: cont_input = & + [ CHARACTER(len=40) :: & + 'create_atoms 1 single &', & + ' 0.2 0.1 0.1' ] +END MODULE keepbox + +FUNCTION f_lammps_with_args() BIND(C, name="f_lammps_with_args") + USE ISO_C_BINDING, ONLY: c_ptr + USE liblammps + USE keepbox, ONLY: lmp + IMPLICIT NONE + TYPE(c_ptr) :: f_lammps_with_args + + CHARACTER(len=12), DIMENSION(12), PARAMETER :: args = & + [ CHARACTER(len=12) :: 'liblammps', '-log', 'none', & + '-echo','screen','-nocite','-var','zpos','1.5','-var','x','2'] + + lmp = lammps(args) + f_lammps_with_args = lmp%handle +END FUNCTION f_lammps_with_args + +SUBROUTINE f_lammps_close() BIND(C, name="f_lammps_close") + USE ISO_C_BINDING, ONLY: c_null_ptr + USE liblammps + USE keepbox, ONLY: lmp + IMPLICIT NONE + + CALL lmp%close() + lmp%handle = c_null_ptr +END SUBROUTINE f_lammps_close + +SUBROUTINE f_lammps_box_setup () BIND(C) + USE liblammps + USE keepbox, ONLY : lmp, demo_input + IMPLICIT NONE + + CALL lmp%commands_list(demo_input) +END SUBROUTINE f_lammps_box_setup + +SUBROUTINE f_lammps_delete_everything() BIND(C) + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + + CALL lmp%command("delete_atoms group all"); +END SUBROUTINE f_lammps_delete_everything + +FUNCTION f_lammps_extract_box_xlo () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_extract_box_xlo + REAL (c_double) :: boxdim(3) + + CALL lmp%extract_box(boxlo=boxdim) + f_lammps_extract_box_xlo = boxdim(1) +END FUNCTION f_lammps_extract_box_xlo + +FUNCTION f_lammps_extract_box_xhi () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_extract_box_xhi + REAL (c_double) :: boxdim(3) + + CALL lmp%extract_box(boxhi=boxdim) + f_lammps_extract_box_xhi = boxdim(1) +END FUNCTION f_lammps_extract_box_xhi + +FUNCTION f_lammps_extract_box_ylo () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_extract_box_ylo + REAL (c_double) :: boxdim(3) + + CALL lmp%extract_box(boxlo=boxdim) + f_lammps_extract_box_ylo = boxdim(2) +END FUNCTION f_lammps_extract_box_ylo + +FUNCTION f_lammps_extract_box_yhi () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_extract_box_yhi + REAL (c_double) :: boxdim(3) + + CALL lmp%extract_box(boxhi=boxdim) + f_lammps_extract_box_yhi = boxdim(2) +END FUNCTION f_lammps_extract_box_yhi + +FUNCTION f_lammps_extract_box_zlo () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_extract_box_zlo + REAL (c_double) :: boxdim(3) + + CALL lmp%extract_box(boxlo=boxdim) + f_lammps_extract_box_zlo = boxdim(2) +END FUNCTION f_lammps_extract_box_zlo + +FUNCTION f_lammps_extract_box_zhi () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_extract_box_zhi + REAL (c_double) :: boxdim(3) + + CALL lmp%extract_box(boxhi=boxdim) + f_lammps_extract_box_zhi = boxdim(2) +END FUNCTION f_lammps_extract_box_zhi + +SUBROUTINE f_lammps_reset_box_2x () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE liblammps + USE keepbox, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: newlo(3), newhi(3), xy, yz, xz + + xy = 0.0_c_double + yz = 0.0_c_double + xz = 0.0_c_double + newlo = [-1.0_c_double, -1.0_c_double, -1.0_c_double] + newhi = [3.0_c_double, 3.0_c_double, 3.0_c_double] + CALL lmp%reset_box(newlo, newhi, xy, yz, xz) +END SUBROUTINE f_lammps_reset_box_2x diff --git a/unittest/fortran/test_fortran_commands.f90 b/unittest/fortran/test_fortran_commands.f90 index b5cffe698f..d85b0f183e 100644 --- a/unittest/fortran/test_fortran_commands.f90 +++ b/unittest/fortran/test_fortran_commands.f90 @@ -1,5 +1,6 @@ MODULE keepcmds USE liblammps + IMPLICIT NONE TYPE(LAMMPS) :: lmp CHARACTER(len=40), DIMENSION(3), PARAMETER :: demo_input = & [ CHARACTER(len=40) :: & diff --git a/unittest/fortran/test_fortran_extract_global.f90 b/unittest/fortran/test_fortran_extract_global.f90 new file mode 100644 index 0000000000..9a8552d677 --- /dev/null +++ b/unittest/fortran/test_fortran_extract_global.f90 @@ -0,0 +1,491 @@ +MODULE keepglobal + USE liblammps + TYPE(LAMMPS) :: lmp + CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: demo_input = & + [ CHARACTER(len=40) :: & + 'region box block 0 $x 0 3 0 4', & + 'create_box 1 box', & + 'create_atoms 1 single 1.0 1.0 ${zpos}' ] + CHARACTER(LEN=40), DIMENSION(2), PARAMETER :: cont_input = & + [ CHARACTER(len=40) :: & + 'create_atoms 1 single &', & + ' 0.2 0.1 0.1' ] + CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: pair_input = & + [ CHARACTER(LEN=40) :: & + 'pair_style lj/cut 2.5', & + 'pair_coeff 1 1 1.0 1.0', & + 'mass 1 1.0' ] +END MODULE keepglobal + +FUNCTION f_lammps_with_args() BIND(C, name="f_lammps_with_args") + USE ISO_C_BINDING, ONLY: c_ptr + USE liblammps + USE keepglobal, ONLY: lmp + IMPLICIT NONE + TYPE(c_ptr) :: f_lammps_with_args + CHARACTER(len=12), DIMENSION(12), PARAMETER :: args = & + [ CHARACTER(len=12) :: 'liblammps', '-log', 'none', & + '-echo','screen','-nocite','-var','zpos','1.5','-var','x','2'] + + lmp = lammps(args) + f_lammps_with_args = lmp%handle +END FUNCTION f_lammps_with_args + +SUBROUTINE f_lammps_close() BIND(C, name="f_lammps_close") + USE ISO_C_BINDING, ONLY: c_null_ptr + USE liblammps + USE keepglobal, ONLY: lmp + IMPLICIT NONE + + CALL lmp%close() + lmp%handle = c_null_ptr +END SUBROUTINE f_lammps_close + +SUBROUTINE f_lammps_setup_extract_global () BIND(C) + USE LIBLAMMPS + USE keepglobal, ONLY : lmp, demo_input, cont_input, pair_input + IMPLICIT NONE + + CALL lmp%commands_list(demo_input) + CALL lmp%commands_list(cont_input) + CALL lmp%commands_list(pair_input) + CALL lmp%command('run 0') +END SUBROUTINE f_lammps_setup_extract_global + +SUBROUTINE f_lammps_setup_full_extract_global () BIND(C) + USE LIBLAMMPS + USE keepglobal, ONLY : lmp + IMPLICIT NONE + INTERFACE + SUBROUTINE f_lammps_setup_extract_global () BIND(C) + END SUBROUTINE f_lammps_setup_extract_global + END INTERFACE + + CALL lmp%command('atom_style full') + CALL f_lammps_setup_extract_global + CALL lmp%command('bond_style zero') + CALL lmp%command('angle_style zero') + CALL lmp%command('dihedral_style zero') + CALL lmp%command('run 0') +END SUBROUTINE f_lammps_setup_full_extract_global + +FUNCTION f_lammps_extract_global_units () BIND(C) RESULT(success) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE LIBLAMMPS + USE keepglobal, ONLY : lmp + IMPLICIT NONE + INTEGER (C_int) :: success + CHARACTER (LEN=16) :: units + + ! passing strings from Fortran to C is icky, so we do the test here and + ! report the result instead + units = lmp%extract_global('units') + IF ( TRIM(units) == 'lj' ) THEN + success = 1_C_int + ELSE + success = 0_C_int + END IF +END FUNCTION f_lammps_extract_global_units + +FUNCTION f_lammps_extract_global_ntimestep () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: ntimestep + INTEGER (C_int) :: f_lammps_extract_global_ntimestep + + ntimestep = lmp%extract_global("ntimestep") + f_lammps_extract_global_ntimestep = ntimestep +END FUNCTION f_lammps_extract_global_ntimestep +FUNCTION f_lammps_extract_global_ntimestep_big () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int64_t), POINTER :: ntimestep + INTEGER (C_int64_t) :: f_lammps_extract_global_ntimestep_big + + ntimestep = lmp%extract_global("ntimestep") + f_lammps_extract_global_ntimestep_big = ntimestep +END FUNCTION f_lammps_extract_global_ntimestep_big + +FUNCTION f_lammps_extract_global_dt () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double), POINTER :: dt + REAL (C_double) :: f_lammps_extract_global_dt + + dt = lmp%extract_global("dt") + f_lammps_extract_global_dt = dt +END FUNCTION f_lammps_extract_global_dt + +SUBROUTINE f_lammps_extract_global_boxlo (C_boxlo) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double), DIMENSION(3) :: C_boxlo + REAL (C_double), DIMENSION(:), POINTER :: boxlo + + boxlo = lmp%extract_global("boxlo") + C_boxlo = boxlo +END SUBROUTINE f_lammps_extract_global_boxlo + +SUBROUTINE f_lammps_extract_global_boxhi (C_boxhi) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double), DIMENSION(3) :: C_boxhi + REAL (C_double), DIMENSION(:), POINTER :: boxhi + + boxhi = lmp%extract_global("boxhi") + C_boxhi = boxhi +END SUBROUTINE f_lammps_extract_global_boxhi + +FUNCTION f_lammps_extract_global_boxxlo () BIND(C) RESULT(C_boxxlo) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_boxxlo + REAL (C_double), POINTER :: boxxlo + + boxxlo = lmp%extract_global("boxxlo") + C_boxxlo = boxxlo +END FUNCTION f_lammps_extract_global_boxxlo + +FUNCTION f_lammps_extract_global_boxxhi () BIND(C) RESULT(C_boxxhi) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_boxxhi + REAL (C_double), POINTER :: boxxhi + + boxxhi = lmp%extract_global("boxxhi") + C_boxxhi = boxxhi +END FUNCTION f_lammps_extract_global_boxxhi + +FUNCTION f_lammps_extract_global_boxylo () BIND(C) RESULT(C_boxylo) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_boxylo + REAL (C_double), POINTER :: boxylo + + boxylo = lmp%extract_global("boxylo") + C_boxylo = boxylo +END FUNCTION f_lammps_extract_global_boxylo + +FUNCTION f_lammps_extract_global_boxyhi () BIND(C) RESULT(C_boxyhi) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_boxyhi + REAL (C_double), POINTER :: boxyhi + + boxyhi = lmp%extract_global("boxyhi") + C_boxyhi = boxyhi +END FUNCTION f_lammps_extract_global_boxyhi + +FUNCTION f_lammps_extract_global_boxzlo () BIND(C) RESULT(C_boxzlo) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_boxzlo + REAL (C_double), POINTER :: boxzlo + + boxzlo = lmp%extract_global("boxzlo") + C_boxzlo = boxzlo +END FUNCTION f_lammps_extract_global_boxzlo + +FUNCTION f_lammps_extract_global_boxzhi () BIND(C) RESULT(C_boxzhi) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_boxzhi + REAL (C_double), POINTER :: boxzhi + + boxzhi = lmp%extract_global("boxzhi") + C_boxzhi = boxzhi +END FUNCTION f_lammps_extract_global_boxzhi + +SUBROUTINE f_lammps_extract_global_periodicity (C_periodicity) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), DIMENSION(3) :: C_periodicity + INTEGER (C_int), DIMENSION(:), POINTER :: periodicity + + periodicity = lmp%extract_global("periodicity") + C_periodicity = periodicity +END SUBROUTINE f_lammps_extract_global_periodicity + +FUNCTION f_lammps_extract_global_triclinic () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: triclinic + INTEGER (C_int) :: f_lammps_extract_global_triclinic + + triclinic = lmp%extract_global("triclinic") + f_lammps_extract_global_triclinic = triclinic +END FUNCTION f_lammps_extract_global_triclinic + +FUNCTION f_lammps_extract_global_xy () BIND(C) RESULT(C_xy) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_xy + REAL (C_double), POINTER :: xy + + xy = lmp%extract_global("xy") + C_xy = xy +END FUNCTION f_lammps_extract_global_xy + +FUNCTION f_lammps_extract_global_xz () BIND(C) RESULT(C_xz) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_xz + REAL (C_double), POINTER :: xz + + xz = lmp%extract_global("xz") + C_xz = xz +END FUNCTION f_lammps_extract_global_xz + +FUNCTION f_lammps_extract_global_yz () BIND(C) RESULT(C_yz) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_yz + REAL (C_double), POINTER :: yz + + yz = lmp%extract_global("yz") + C_yz = yz +END FUNCTION f_lammps_extract_global_yz + +FUNCTION f_lammps_extract_global_natoms () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: natoms + INTEGER (C_int) :: f_lammps_extract_global_natoms + + natoms = lmp%extract_global("natoms") + f_lammps_extract_global_natoms = natoms +END FUNCTION f_lammps_extract_global_natoms +FUNCTION f_lammps_extract_global_natoms_big () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int64_t), POINTER :: natoms + INTEGER (C_int64_t) :: f_lammps_extract_global_natoms_big + + natoms = lmp%extract_global("natoms") + f_lammps_extract_global_natoms_big = natoms +END FUNCTION f_lammps_extract_global_natoms_big + +FUNCTION f_lammps_extract_global_nbonds () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: nbonds + INTEGER (C_int) :: f_lammps_extract_global_nbonds + + nbonds = lmp%extract_global("nbonds") + f_lammps_extract_global_nbonds = nbonds +END FUNCTION f_lammps_extract_global_nbonds +FUNCTION f_lammps_extract_global_nbonds_big () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int64_t), POINTER :: nbonds + INTEGER (C_int64_t) :: f_lammps_extract_global_nbonds_big + + nbonds = lmp%extract_global("nbonds") + f_lammps_extract_global_nbonds_big = nbonds +END FUNCTION f_lammps_extract_global_nbonds_big + +FUNCTION f_lammps_extract_global_nangles () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: nangles + INTEGER (C_int) :: f_lammps_extract_global_nangles + + nangles = lmp%extract_global("nangles") + f_lammps_extract_global_nangles = nangles +END FUNCTION f_lammps_extract_global_nangles +FUNCTION f_lammps_extract_global_nangles_big () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int64_t), POINTER :: nangles + INTEGER (C_int64_t) :: f_lammps_extract_global_nangles_big + + nangles = lmp%extract_global("nangles") + f_lammps_extract_global_nangles_big = nangles +END FUNCTION f_lammps_extract_global_nangles_big + +FUNCTION f_lammps_extract_global_ndihedrals () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: ndihedrals + INTEGER (C_int) :: f_lammps_extract_global_ndihedrals + + ndihedrals = lmp%extract_global("ndihedrals") + f_lammps_extract_global_ndihedrals = ndihedrals +END FUNCTION f_lammps_extract_global_ndihedrals +FUNCTION f_lammps_extract_global_ndihedrals_big () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int64_t), POINTER :: ndihedrals + INTEGER (C_int64_t) :: f_lammps_extract_global_ndihedrals_big + + ndihedrals = lmp%extract_global("ndihedrals") + f_lammps_extract_global_ndihedrals_big = ndihedrals +END FUNCTION f_lammps_extract_global_ndihedrals_big + +FUNCTION f_lammps_extract_global_nimpropers () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: nimpropers + INTEGER (C_int) :: f_lammps_extract_global_nimpropers + + nimpropers = lmp%extract_global("nimpropers") + f_lammps_extract_global_nimpropers = nimpropers +END FUNCTION f_lammps_extract_global_nimpropers +FUNCTION f_lammps_extract_global_nimpropers_big () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int64_t + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int64_t), POINTER :: nimpropers + INTEGER (C_int64_t) :: f_lammps_extract_global_nimpropers_big + + nimpropers = lmp%extract_global("nimpropers") + f_lammps_extract_global_nimpropers_big = nimpropers +END FUNCTION f_lammps_extract_global_nimpropers_big + + +FUNCTION f_lammps_extract_global_ntypes () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: ntypes + INTEGER (C_int) :: f_lammps_extract_global_ntypes + + ntypes = lmp%extract_global("ntypes") + f_lammps_extract_global_ntypes = ntypes +END FUNCTION f_lammps_extract_global_ntypes + +FUNCTION f_lammps_extract_global_nlocal () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: nlocal + INTEGER (C_int) :: f_lammps_extract_global_nlocal + + nlocal = lmp%extract_global("nlocal") + f_lammps_extract_global_nlocal = nlocal +END FUNCTION f_lammps_extract_global_nlocal + +FUNCTION f_lammps_extract_global_nghost () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: nghost + INTEGER (C_int) :: f_lammps_extract_global_nghost + + nghost = lmp%extract_global("nghost") + f_lammps_extract_global_nghost = nghost +END FUNCTION f_lammps_extract_global_nghost + +FUNCTION f_lammps_extract_global_nmax () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int), POINTER :: nmax + INTEGER (C_int) :: f_lammps_extract_global_nmax + + nmax = lmp%extract_global("nmax") + f_lammps_extract_global_nmax = nmax +END FUNCTION f_lammps_extract_global_nmax + +FUNCTION f_lammps_extract_global_boltz () BIND(C) RESULT(C_k_B) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_k_B + REAL (C_double), POINTER :: k_B + + k_B = lmp%extract_global("boltz") + C_k_B = k_B +END FUNCTION f_lammps_extract_global_boltz + +FUNCTION f_lammps_extract_global_hplanck () BIND(C) RESULT(C_h) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: C_h + REAL (C_double), POINTER :: h + + h = lmp%extract_global("boltz") + C_h = h +END FUNCTION f_lammps_extract_global_hplanck + +FUNCTION f_lammps_extract_global_angstrom () BIND(C) RESULT(Angstrom) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: Angstrom + REAL (C_double), POINTER :: A + + A = lmp%extract_global("angstrom") + Angstrom = A +END FUNCTION f_lammps_extract_global_angstrom + +FUNCTION f_lammps_extract_global_femtosecond () BIND(C) RESULT(fs) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE keepglobal, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + REAL (C_double) :: fs + REAL (C_double), POINTER :: femtosecond + + femtosecond = lmp%extract_global("femtosecond") + fs = femtosecond +END FUNCTION f_lammps_extract_global_femtosecond diff --git a/unittest/fortran/test_fortran_get_thermo.f90 b/unittest/fortran/test_fortran_get_thermo.f90 new file mode 100644 index 0000000000..e96c964e7c --- /dev/null +++ b/unittest/fortran/test_fortran_get_thermo.f90 @@ -0,0 +1,174 @@ +MODULE keepthermo + USE liblammps + IMPLICIT NONE + TYPE(LAMMPS) :: lmp + CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: demo_input = & + [ CHARACTER(len=40) :: & + 'region box block 0 $x 0 3 0 4', & + 'create_box 1 box', & + 'create_atoms 1 single 1.0 1.0 ${zpos}' ] + CHARACTER(LEN=40), DIMENSION(2), PARAMETER :: cont_input = & + [ CHARACTER(len=40) :: & + 'create_atoms 1 single &', & + ' 0.2 0.1 0.1' ] + CHARACTER(LEN=40), DIMENSION(3), PARAMETER :: pair_input = & + [ CHARACTER(LEN=40) :: & + 'pair_style lj/cut 2.5', & + 'pair_coeff 1 1 1.0 1.0', & + 'mass 1 1.0' ] +END MODULE keepthermo + +FUNCTION f_lammps_with_args() BIND(C) + USE ISO_C_BINDING, ONLY: c_ptr + USE liblammps + USE keepthermo, ONLY: lmp + IMPLICIT NONE + TYPE(c_ptr) :: f_lammps_with_args + + CHARACTER(len=12), DIMENSION(12), PARAMETER :: args = & + [ CHARACTER(len=12) :: 'liblammps', '-log', 'none', & + '-echo','screen','-nocite','-var','zpos','1.5','-var','x','2'] + + lmp = lammps(args) + f_lammps_with_args = lmp%handle +END FUNCTION f_lammps_with_args + +SUBROUTINE f_lammps_close() BIND(C) + USE ISO_C_BINDING, ONLY: c_null_ptr + USE liblammps + USE keepthermo, ONLY: lmp + IMPLICIT NONE + + CALL lmp%close() + lmp%handle = c_null_ptr +END SUBROUTINE f_lammps_close + +SUBROUTINE f_lammps_get_thermo_setup () BIND(C) + USE liblammps + USE keepthermo, ONLY : lmp, demo_input, cont_input, pair_input + IMPLICIT NONE + + CALL lmp%commands_list(demo_input) + CALL lmp%commands_list(cont_input) + CALL lmp%commands_list(pair_input) +END SUBROUTINE f_lammps_get_thermo_setup + +FUNCTION f_lammps_get_thermo_natoms () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_natoms + + f_lammps_get_thermo_natoms = lmp%get_thermo('atoms') +END FUNCTION f_lammps_get_thermo_natoms + +FUNCTION f_lammps_get_thermo_dt () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_dt + + f_lammps_get_thermo_dt = lmp%get_thermo('dt') +END FUNCTION f_lammps_get_thermo_dt + +FUNCTION f_lammps_get_thermo_vol () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_vol + + f_lammps_get_thermo_vol = lmp%get_thermo('vol') +END FUNCTION f_lammps_get_thermo_vol + +FUNCTION f_lammps_get_thermo_lx () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_lx + + f_lammps_get_thermo_lx = lmp%get_thermo('lx') +END FUNCTION f_lammps_get_thermo_lx + +FUNCTION f_lammps_get_thermo_ly () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_ly + + f_lammps_get_thermo_ly = lmp%get_thermo('ly') +END FUNCTION f_lammps_get_thermo_ly + +FUNCTION f_lammps_get_thermo_lz () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_lz + + f_lammps_get_thermo_lz = lmp%get_thermo('lz') +END FUNCTION f_lammps_get_thermo_lz + +FUNCTION f_lammps_get_thermo_xlo () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_xlo + + f_lammps_get_thermo_xlo = lmp%get_thermo('xlo') +END FUNCTION f_lammps_get_thermo_xlo + +FUNCTION f_lammps_get_thermo_xhi () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_xhi + + f_lammps_get_thermo_xhi = lmp%get_thermo('xhi') +END FUNCTION f_lammps_get_thermo_xhi + +FUNCTION f_lammps_get_thermo_ylo () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_ylo + + f_lammps_get_thermo_ylo = lmp%get_thermo('ylo') +END FUNCTION f_lammps_get_thermo_ylo + +FUNCTION f_lammps_get_thermo_yhi () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_yhi + + f_lammps_get_thermo_yhi = lmp%get_thermo('yhi') +END FUNCTION f_lammps_get_thermo_yhi + +FUNCTION f_lammps_get_thermo_zlo () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_zlo + + f_lammps_get_thermo_zlo = lmp%get_thermo('zlo') +END FUNCTION f_lammps_get_thermo_zlo + +FUNCTION f_lammps_get_thermo_zhi () BIND (C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_double + USE liblammps + USE keepthermo, ONLY : lmp + IMPLICIT NONE + REAL (c_double) :: f_lammps_get_thermo_zhi + + f_lammps_get_thermo_zhi = lmp%get_thermo('zhi') +END FUNCTION f_lammps_get_thermo_zhi diff --git a/unittest/fortran/test_fortran_properties.f90 b/unittest/fortran/test_fortran_properties.f90 new file mode 100644 index 0000000000..00e083b301 --- /dev/null +++ b/unittest/fortran/test_fortran_properties.f90 @@ -0,0 +1,52 @@ +FUNCTION f_lammps_version () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE liblammps + USE keepcmds, ONLY : lmp + IMPLICIT NONE + INTEGER (C_int) :: f_lammps_version + + f_lammps_version = lmp%version() +END FUNCTION f_lammps_version + +SUBROUTINE f_lammps_memory_usage (meminfo) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_double + USE liblammps + USE keepcmds, ONLY : lmp + IMPLICIT NONE + REAL (C_double), DIMENSION(3), INTENT(OUT) :: meminfo + + CALL lmp%memory_usage(meminfo) +END SUBROUTINE f_lammps_memory_usage + +FUNCTION f_lammps_get_mpi_comm () BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int + USE liblammps + USE keepcmds, ONLY : lmp + IMPLICIT NONE + INTEGER (C_int) :: f_lammps_get_mpi_comm + + f_lammps_get_mpi_comm = lmp%get_mpi_comm() +END FUNCTION f_lammps_get_mpi_comm + +FUNCTION f_lammps_extract_setting (Cstr) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_int, C_char + USE keepcmds, ONLY : lmp + USE LIBLAMMPS + IMPLICIT NONE + INTEGER (C_int) :: f_lammps_extract_setting + CHARACTER (KIND=C_char, LEN=1), DIMENSION(*), INTENT(IN) :: Cstr + INTEGER :: strlen, i + CHARACTER (LEN=:), ALLOCATABLE :: Fstr + + i = 1 + DO WHILE (Cstr(i) /= ACHAR(0)) + i = i + 1 + END DO + strlen = i + allocate ( CHARACTER(LEN=strlen) :: Fstr) + FORALL (i=1:strlen) + Fstr(i:i) = Cstr(i) + END FORALL + f_lammps_extract_setting = lmp%extract_setting(Fstr) + deallocate (Fstr) +END FUNCTION f_lammps_extract_setting diff --git a/unittest/fortran/wrap_box.cpp b/unittest/fortran/wrap_box.cpp new file mode 100644 index 0000000000..8678816658 --- /dev/null +++ b/unittest/fortran/wrap_box.cpp @@ -0,0 +1,64 @@ +// unit tests for extracting box dimensions fom a LAMMPS instance through the Fortran wrapper + +#include "lammps.h" +#include +#include + +#include "gtest/gtest.h" + +// prototypes for fortran reverse wrapper functions +extern "C" { +void *f_lammps_with_args(); +void f_lammps_close(); +void f_lammps_box_setup(); +double f_lammps_extract_box_xlo(); +double f_lammps_extract_box_xhi(); +double f_lammps_extract_box_ylo(); +double f_lammps_extract_box_yhi(); +double f_lammps_extract_box_zlo(); +double f_lammps_extract_box_zhi(); +void f_lammps_delete_everything(); +void f_lammps_reset_box_2x(); +} + +class LAMMPS_commands : public ::testing::Test { +protected: + LAMMPS_NS::LAMMPS *lmp; + LAMMPS_commands() = default; + ~LAMMPS_commands() override = default; + + void SetUp() override + { + ::testing::internal::CaptureStdout(); + lmp = (LAMMPS_NS::LAMMPS *)f_lammps_with_args(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 8).c_str(), "LAMMPS ("); + } + void TearDown() override + { + ::testing::internal::CaptureStdout(); + f_lammps_close(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 16).c_str(), "Total wall time:"); + lmp = nullptr; + } +}; + +TEST_F(LAMMPS_commands, get_thermo) +{ + f_lammps_box_setup(); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_xlo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_xhi(), 2.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_ylo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_yhi(), 2.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_zlo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_zhi(), 2.0); + f_lammps_delete_everything(); + f_lammps_reset_box_2x(); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_xlo(),-1.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_xhi(), 3.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_ylo(),-1.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_yhi(), 3.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_zlo(),-1.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_box_zhi(), 3.0); +}; diff --git a/unittest/fortran/wrap_extract_global.cpp b/unittest/fortran/wrap_extract_global.cpp new file mode 100644 index 0000000000..adf3986073 --- /dev/null +++ b/unittest/fortran/wrap_extract_global.cpp @@ -0,0 +1,177 @@ +// unit tests for extracting global data from a LAMMPS instance through the +// Fortran wrapper + +#include "lammps.h" +#include "library.h" +#include +#include +#include +#include + +#include "gtest/gtest.h" + +// prototypes for Fortran reverse wrapper functions +extern "C" { +void *f_lammps_with_args(); +void f_lammps_close(); +void f_lammps_setup_extract_global(); +void f_lammps_setup_full_extract_global(); +int f_lammps_extract_global_units(); +int f_lammps_extract_global_ntimestep(); +int64_t f_lammps_extract_global_ntimestep_big(); +double f_lammps_extract_global_dt(); +void f_lammps_extract_global_boxlo(double[3]); +void f_lammps_extract_global_boxhi(double[3]); +double f_lammps_extract_global_boxxlo(); +double f_lammps_extract_global_boxylo(); +double f_lammps_extract_global_boxzlo(); +double f_lammps_extract_global_boxxhi(); +double f_lammps_extract_global_boxyhi(); +double f_lammps_extract_global_boxzhi(); +void f_lammps_extract_global_periodicity(int[3]); +int f_lammps_extract_global_triclinic(); +double f_lammps_extract_global_xy(); +double f_lammps_extract_global_yz(); +double f_lammps_extract_global_xz(); +int f_lammps_extract_global_natoms(); +int64_t f_lammps_extract_global_natoms_big(); +int f_lammps_extract_global_nbonds(); +int64_t f_lammps_extract_global_nbonds_big(); +int f_lammps_extract_global_nangles(); +int64_t f_lammps_extract_global_nangles_big(); +int f_lammps_extract_global_ndihedrals(); +int64_t f_lammps_extract_global_ndihedrals_big(); +int f_lammps_extract_global_nimpropers(); +int64_t f_lammps_extract_global_nimpropers_big(); +int f_lammps_extract_global_ntypes(); +int f_lammps_extract_global_nlocal(); +int f_lammps_extract_global_nghost(); +int f_lammps_extract_global_nmax(); +double f_lammps_extract_global_boltz(); +double f_lammps_extract_global_hplanck(); +double f_lammps_extract_global_angstrom(); +double f_lammps_extract_global_femtosecond(); +} + +class LAMMPS_extract_global : public ::testing::Test { +protected: + LAMMPS_NS::LAMMPS *lmp; + LAMMPS_extract_global() = default; + ~LAMMPS_extract_global() override = default; + + void SetUp() override + { + ::testing::internal::CaptureStdout(); + lmp = (LAMMPS_NS::LAMMPS *)f_lammps_with_args(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 8).c_str(), "LAMMPS ("); + } + void TearDown() override + { + ::testing::internal::CaptureStdout(); + f_lammps_close(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 16).c_str(), "Total wall time:"); + lmp = nullptr; + } +}; + +TEST_F(LAMMPS_extract_global, units) +{ + f_lammps_setup_extract_global(); + EXPECT_EQ(f_lammps_extract_global_units(), 1); +}; + +TEST_F(LAMMPS_extract_global, ntimestep) +{ + f_lammps_setup_extract_global(); +#ifdef LAMMPS_SMALLSMALL + EXPECT_EQ(f_lammps_extract_global_ntimestep(), 0); +#else + EXPECT_EQ(f_lammps_extract_global_ntimestep_big(), 0l); +#endif +}; + +TEST_F(LAMMPS_extract_global, dt) +{ + f_lammps_setup_extract_global(); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_dt(), 0.005); +}; + +TEST_F(LAMMPS_extract_global, boxprops) +{ + f_lammps_setup_extract_global(); + double boxlo[3], boxhi[3]; + f_lammps_extract_global_boxlo(boxlo); + EXPECT_DOUBLE_EQ(boxlo[0], 0.0); + EXPECT_DOUBLE_EQ(boxlo[1], 0.0); + EXPECT_DOUBLE_EQ(boxlo[2], 0.0); + f_lammps_extract_global_boxhi(boxhi); + EXPECT_DOUBLE_EQ(boxhi[0], 2.0); + EXPECT_DOUBLE_EQ(boxhi[1], 3.0); + EXPECT_DOUBLE_EQ(boxhi[2], 4.0); + + EXPECT_DOUBLE_EQ(f_lammps_extract_global_boxxlo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_boxxhi(), 2.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_boxylo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_boxyhi(), 3.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_boxzlo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_boxzhi(), 4.0); + + int periodicity[3]; + f_lammps_extract_global_periodicity(periodicity); + EXPECT_EQ(periodicity[0], 1); + EXPECT_EQ(periodicity[1], 1); + EXPECT_EQ(periodicity[2], 1); + + EXPECT_EQ(f_lammps_extract_global_triclinic(), 0); + + EXPECT_DOUBLE_EQ(f_lammps_extract_global_xy(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_yz(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_xz(), 0.0); +}; + +TEST_F(LAMMPS_extract_global, atomprops) +{ + f_lammps_setup_extract_global(); +#ifdef LAMMPS_SMALLSMALL + EXPECT_EQ(f_lammps_extract_global_natoms(), 2); + EXPECT_EQ(f_lammps_extract_global_nbonds(), 0); + EXPECT_EQ(f_lammps_extract_global_nangles(), 0); + EXPECT_EQ(f_lammps_extract_global_ndihedrals(), 0); +#else + EXPECT_EQ(f_lammps_extract_global_natoms_big(), 2l); + EXPECT_EQ(f_lammps_extract_global_nbonds_big(), 0l); + EXPECT_EQ(f_lammps_extract_global_nangles_big(), 0l); + EXPECT_EQ(f_lammps_extract_global_ndihedrals_big(), 0l); +#endif + + EXPECT_EQ(f_lammps_extract_global_ntypes(), 1); + EXPECT_EQ(f_lammps_extract_global_nlocal(), 2); + EXPECT_EQ(f_lammps_extract_global_nghost(), 41); + EXPECT_EQ(f_lammps_extract_global_nmax(), 16384); + + EXPECT_DOUBLE_EQ(f_lammps_extract_global_boltz(), 1.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_hplanck(), 1.0); + EXPECT_DOUBLE_EQ(f_lammps_extract_global_angstrom(), 1.0); + + EXPECT_DOUBLE_EQ(f_lammps_extract_global_femtosecond(), 1.0); +}; + +TEST_F(LAMMPS_extract_global, fullprops) +{ + if (! lammps_has_style(lmp, "atom", "full")) GTEST_SKIP(); + // This is not currently the world's most convincing test.... + f_lammps_setup_full_extract_global(); +#ifdef LAMMPS_SMALLSMALL + EXPECT_EQ(f_lammps_extract_global_natoms(), 2); + EXPECT_EQ(f_lammps_extract_global_nbonds(), 0); + EXPECT_EQ(f_lammps_extract_global_nangles(), 0); + EXPECT_EQ(f_lammps_extract_global_ndihedrals(), 0); +#else + EXPECT_EQ(f_lammps_extract_global_natoms_big(), 2l); + EXPECT_EQ(f_lammps_extract_global_nbonds_big(), 0l); + EXPECT_EQ(f_lammps_extract_global_nangles_big(), 0l); + EXPECT_EQ(f_lammps_extract_global_ndihedrals_big(), 0l); +#endif +} diff --git a/unittest/fortran/wrap_get_thermo.cpp b/unittest/fortran/wrap_get_thermo.cpp new file mode 100644 index 0000000000..71b51609eb --- /dev/null +++ b/unittest/fortran/wrap_get_thermo.cpp @@ -0,0 +1,67 @@ +// unit tests for getting thermodynamic output from a LAMMPS instance through the Fortran wrapper + +#include "lammps.h" +#include +#include + +#include "gtest/gtest.h" + +// prototypes for fortran reverse wrapper functions +extern "C" { +void *f_lammps_with_args(); +void f_lammps_close(); +void f_lammps_get_thermo_setup(); +double f_lammps_get_thermo_natoms(); +double f_lammps_get_thermo_dt(); +double f_lammps_get_thermo_vol(); +double f_lammps_get_thermo_lx(); +double f_lammps_get_thermo_ly(); +double f_lammps_get_thermo_lz(); +double f_lammps_get_thermo_xlo(); +double f_lammps_get_thermo_xhi(); +double f_lammps_get_thermo_ylo(); +double f_lammps_get_thermo_yhi(); +double f_lammps_get_thermo_zlo(); +double f_lammps_get_thermo_zhi(); +} + +class LAMMPS_thermo : public ::testing::Test { +protected: + LAMMPS_NS::LAMMPS *lmp; + LAMMPS_thermo() = default; + ~LAMMPS_thermo() override = default; + + void SetUp() override + { + ::testing::internal::CaptureStdout(); + lmp = (LAMMPS_NS::LAMMPS *)f_lammps_with_args(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 8).c_str(), "LAMMPS ("); + } + void TearDown() override + { + ::testing::internal::CaptureStdout(); + f_lammps_close(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 16).c_str(), "Total wall time:"); + lmp = nullptr; + } +}; + +TEST_F(LAMMPS_thermo, get_thermo) +{ + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_natoms(), 0.0); + f_lammps_get_thermo_setup(); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_natoms(), 2.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_dt(), 0.005); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_vol(), 24.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_lx(), 2.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_ly(), 3.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_lz(), 4.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_xlo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_xhi(), 2.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_ylo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_yhi(), 3.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_zlo(), 0.0); + EXPECT_DOUBLE_EQ(f_lammps_get_thermo_zhi(), 4.0); +}; diff --git a/unittest/fortran/wrap_properties.cpp b/unittest/fortran/wrap_properties.cpp new file mode 100644 index 0000000000..8ecd9346dc --- /dev/null +++ b/unittest/fortran/wrap_properties.cpp @@ -0,0 +1,109 @@ +// unit tests for getting LAMMPS properties through the Fortran wrapper + +#include "lammps.h" +//#include // for stdin, stdout +#include "library.h" +#include +#include + +#include "gtest/gtest.h" + +// prototypes for fortran reverse wrapper functions +extern "C" { +void *f_lammps_with_args(); +void f_lammps_close(); +int f_lammps_version(); +void f_lammps_memory_usage(double*); +int f_lammps_get_mpi_comm(); +int f_lammps_extract_setting(const char*); +} + +class LAMMPS_properties : public ::testing::Test { +protected: + LAMMPS_NS::LAMMPS *lmp; + LAMMPS_properties() = default; + ~LAMMPS_properties() override = default; + + void SetUp() override + { + ::testing::internal::CaptureStdout(); + lmp = (LAMMPS_NS::LAMMPS *)f_lammps_with_args(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 8).c_str(), "LAMMPS ("); + } + void TearDown() override + { + ::testing::internal::CaptureStdout(); + f_lammps_close(); + std::string output = ::testing::internal::GetCapturedStdout(); + EXPECT_STREQ(output.substr(0, 16).c_str(), "Total wall time:"); + lmp = nullptr; + } +}; + +TEST_F(LAMMPS_properties, version) +{ + EXPECT_LT(20200917, f_lammps_version()); +}; + +TEST_F(LAMMPS_properties, memory_usage) +{ +// copied from c-library, with a two-character modification + double meminfo[3]; + f_lammps_memory_usage(meminfo); + EXPECT_GT(meminfo[0], 0.0); +#if defined(__linux__) || defined(_WIN32) + EXPECT_GE(meminfo[1], 0.0); +#endif +#if (defined(__linux__) || defined(__APPLE__) || defined(_WIN32)) && !defined(__INTEL_LLVM_COMPILER) + EXPECT_GT(meminfo[2], 0.0); +#endif +}; + +TEST_F(LAMMPS_properties, get_mpi_comm) +{ + int f_comm = f_lammps_get_mpi_comm(); + if ( lammps_config_has_mpi_support() ) + EXPECT_GE(f_comm, 0); + else + EXPECT_EQ(f_comm, -1); +}; + +TEST_F(LAMMPS_properties, extract_setting) +{ +#if defined(LAMMPS_SMALLSMALL) + EXPECT_EQ(f_lammps_extract_setting("bigint"), 4); +#else + EXPECT_EQ(f_lammps_extract_setting("bigint"), 8); +#endif +#if defined(LAMMPS_BIGBIG) + EXPECT_EQ(f_lammps_extract_setting("tagint"), 8); + EXPECT_EQ(f_lammps_extract_setting("imageint"), 8); +#else + EXPECT_EQ(f_lammps_extract_setting("tagint"), 4); + EXPECT_EQ(f_lammps_extract_setting("imageint"), 4); +#endif + + EXPECT_EQ(f_lammps_extract_setting("box_exist"), 0); + EXPECT_EQ(f_lammps_extract_setting("dimension"), 3); + EXPECT_EQ(f_lammps_extract_setting("world_size"), 1); + EXPECT_EQ(f_lammps_extract_setting("world_rank"), 0); + EXPECT_EQ(f_lammps_extract_setting("universe_size"), 1); + EXPECT_EQ(f_lammps_extract_setting("universe_rank"), 0); + EXPECT_GT(f_lammps_extract_setting("nthreads"), 0); + EXPECT_EQ(f_lammps_extract_setting("newton_pair"), 1); + EXPECT_EQ(f_lammps_extract_setting("newton_bond"), 1); + + EXPECT_EQ(f_lammps_extract_setting("ntypes"), 0); + EXPECT_EQ(f_lammps_extract_setting("nbondtypes"), 0); + EXPECT_EQ(f_lammps_extract_setting("nangletypes"), 0); + EXPECT_EQ(f_lammps_extract_setting("ndihedraltypes"), 0); + EXPECT_EQ(f_lammps_extract_setting("nimpropertypes"), 0); + + EXPECT_EQ(f_lammps_extract_setting("molecule_flag"), 0); + EXPECT_EQ(f_lammps_extract_setting("q_flag"), 0); + EXPECT_EQ(f_lammps_extract_setting("mu_flag"), 0); + EXPECT_EQ(f_lammps_extract_setting("rmass_flag"), 0); + EXPECT_EQ(f_lammps_extract_setting("UNKNOWN"), -1); + +}; diff --git a/unittest/python/CMakeLists.txt b/unittest/python/CMakeLists.txt index 0aebc6fc7c..78b0d18446 100644 --- a/unittest/python/CMakeLists.txt +++ b/unittest/python/CMakeLists.txt @@ -44,6 +44,7 @@ else() endif() list(APPEND PYTHON_TEST_ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}") list(APPEND PYTHON_TEST_ENVIRONMENT "PYTHONUNBUFFERED=1") +list(APPEND PYTHON_TEST_ENVIRONMENT "PYTHONDONTWRITEBYTECODE=1") if(APPLE) list(APPEND PYTHON_TEST_ENVIRONMENT "DYLD_LIBRARY_PATH=${LAMMPS_LIB_PATH}:$ENV{DYLD_LIBRARY_PATH}") elseif(WIN32) diff --git a/unittest/python/python-open.py b/unittest/python/python-open.py index 328745ded0..3c9165b8c8 100644 --- a/unittest/python/python-open.py +++ b/unittest/python/python-open.py @@ -45,11 +45,21 @@ class PythonOpen(unittest.TestCase): def testWithArgs(self): """Create LAMMPS instance with a few arguments""" - lmp=lammps(name=self.machine, - cmdargs=['-nocite','-sf','opt','-log','none']) + lmp=lammps(name=self.machine,cmdargs=['-nocite','-sf','opt','-log','none']) self.assertIsNot(lmp.lmp,None) self.assertEqual(lmp.opened,1) + def testError(self): + """Print warning message through LAMMPS Error class""" + lmp=lammps(name=self.machine,cmdargs=['-nocite','-log','none','-screen','tmp.error.output']) + lmp.error(0,'test_warning') + lmp.close() + with open('tmp.error.output','r') as f: + output = f.read() + self.assertTrue('LAMMPS' in output) + self.assertTrue('Total wall time' in output) + self.assertTrue('WARNING: test_warning' in output) + def testContextManager(self): """Automatically clean up LAMMPS instance""" with lammps(name=self.machine) as lmp: